Jobs/Burst based substepping

[gotmachine]

Replace the managed thread based substepping sim with a Jobs/Burst implementation.

Overview of the substep sim :

In SubStepSim.OnFixedUpdate() (called by Kerbalism.FixedUpdate() before anything else has been done) :

• Wait for your background worker to have finished (if he hasn't already)
• For every simulated vessel, from the collection used by the background workers to store "future steps" for this vessel, remove all steps whose UT < the current UT, and store them in a collection accessible from the VesselData class with managed code, and that can safely be read/write from the main thread.
• Update your "job-safe" versions of all vessels/bodies and their orbits with the current KSP position/frame/parameters
• Compute maxFutureUT, the maximum time in the future that might be attained in the next fixedupdate at max timewarp : currentUT + (TimeWarp.fetch.warpRates[7] * 0.02) + one extra substep just in case
• Start your background worker, computing all steps for all vessels, starting from the last already stored step UT up to maxFutureUT

Later in Kerbalism.FixedUpdate(), VesselData.EnvironmentUpdate will be called.
The elapsedSeconds parameter is checked against the substep duration "constant" (60s for example).

If elapsedSecond < substepDuration * 2.0 :

• Remove all steps whose UT < current UT from the "main thread-safe" step collection that you updated in SubStepSim.OnFixedUpdate()
• Create a new step
• Directly copy the vessel position and bodies position from the current ones as defined by KSP to the "job-compatible" data structures (or, in case that's possible, you can use the ones you already updated in SubStepSim.OnFixedUpdate())
• Run the step evaluation (occlusion, fluxes, etc)
• Wait for the results and copy them back in the corresponding VesselData managed data structures

else

• Get the "main thread-safe" collection of steps for that vessel that you updated in SubStepSim.OnFixedUpdate()
• Iterate on those to compute scalar factors or average values for the results (occlusion, fluxes, etc)
• Copy those results in the corresponding VesselData managed data structures
• Clear the "main thread-safe" step collection

Data : current implementation performance

• Single threaded blocking implementation (main thread will wait for the worker to finish)
• 30s substeps
• 20 vessels in orbit around various stock bodies
• Space center scene, max stock warp (100.000x)

Save : 20 vessels.zip

Profiler results (2C/[email protected], Cinebench R20 MT : 1135 / ST : 453) :

Relevant metrics :

• Sim thread : substepping thread processing time
• Sim lag : time spent blocking the main thread waiting for the sim thread to finish
• Sim sync : time spent synchronizing data between KSP/Kerbalism objects and the thread safe objects
• Sim processing : time spent processing the substeps in the main thread

[gotmachine]

Inclusion of radiation evaluation

Like irradiances, radiation evaluation is also subject to huge inaccuracies at high timewarp rates. Ideally, we should move it to the substep simulation.
Said plainly, that whole portion of code :

Kerbalism4/src/Kerbalism/Database/VesselData/VesselData.cs

Lines 979 to 1002 in 1e390eb

	gammaTransparency = Sim.GammaTransparency(Vessel.mainBody, Vessel.altitude);
	bool new_innerBelt, new_outerBelt, new_magnetosphere;
	radiation = Radiation.Compute(this, position, EnvGammaTransparency, mainStar.sunlightFactor, out blackout, out new_magnetosphere, out new_innerBelt, out new_outerBelt, out interstellar);
	if (new_innerBelt != innerBelt || new_outerBelt != outerBelt || new_magnetosphere != magnetosphere)
	{
	innerBelt = new_innerBelt;
	outerBelt = new_outerBelt;
	magnetosphere = new_magnetosphere;
	if (IsSimulated) API.OnRadiationFieldChanged.Notify(Vessel, innerBelt, outerBelt, magnetosphere);
	}
	thermosphere = Sim.InsideThermosphere(Vessel);
	exosphere = Sim.InsideExosphere(Vessel);
	if (Storm.InProgress(this))
	{
	double sunActivity = Radiation.Info(mainStar.Star.body).SolarActivity(false) / 2.0;
	stormRadiation = PreferencesRadiation.Instance.StormRadiation * mainStar.sunlightFactor * (sunActivity + 0.5);
	}
	else
	{
	stormRadiation = 0.0;
	}

In terms of final output, this translate into getting substepped averages (instead of "current position") for the following VesselData radiation rates (rad/s):

/// <summary> ambiant radiation from the inner belt</summary>
public virtual double EnvRadiationInnerBelt { get; }

/// <summary> ambiant radiation from the outer belt</summary>
public virtual double EnvRadiationOuterBelt { get; }

/// <summary> ambiant radiation from the magnetopause</summary>
public virtual double EnvRadiationMagnetopause { get; }

/// <summary> directional radiation from the surface</summary>
public virtual double EnvRadiationBodies { get; }

/// <summary> directional radiation from the sun(s)</summary>
public virtual double EnvRadiationSolar { get; }

 /// <summary> directional radiation from solar storm</summary>
public virtual double EnvStormRadiation { get; }

/// <summary> proportion of ionizing radiation not blocked by atmosphere</summary>
public virtual double EnvGammaTransparency  { get; }

The current code is a bit messy. The core method to port is the Radiation.Compute() method (which is currently a bit messy and was pending a refactor/cleanup) :

Kerbalism4/src/Kerbalism/Radiation/Radiation.cs

Lines 618 to 770 in 1e390eb

	/// <summary> return the total environent radiation at position specified </summary>
	public static double Compute(VesselData vd, Vector3d position, double gammaTransparency, double sunlight, out bool blackout,
	out bool magnetosphere, out bool innerBelt, out bool outerBelt, out bool interstellar)
	{
	// prepare out parameters
	blackout = false;
	magnetosphere = false;
	innerBelt = false;
	outerBelt = false;
	interstellar = false;
	// no-op when Radiation is disabled
	if (!Features.Radiation) return 0.0;
	// store stuff
	Space gsm;
	Vector3 p;
	double D;
	double r;
	// accumulate radiation
	double radiation = 0.0;
	double mainBodyRadiation = 0.0;
	CelestialBody body = vd.Vessel.mainBody;
	vd.radiationInnerBelt = 0.0;
	vd.radiationOuterBelt = 0.0;
	vd.radiationMagnetopause = 0.0;
	vd.radiationBodies = 0.0;
	vd.radiationSolar = 0.0;
	while (body != null)
	{
	// Compute radiation values from overlapping 3d fields (belts + magnetospheres)
	RadiationBody rb = Info(body);
	RadiationModel mf = rb.model;
	// activity is [-0.15..1.05]
	var activity = rb.SolarActivity(false);
	if (mf.HasField())
	{
	// transform to local space once
	var scaled_position = ScaledSpace.LocalToScaledSpace(position);
	// generate radii-normalized GSM space
	gsm = GsmSpace(rb, true);
	// move the point in GSM space
	p = gsm.TransformIn(scaled_position);
	// accumulate radiation and determine pause/belt flags
	if (mf.has_inner)
	{
	D = mf.InnerFunc(p);
	innerBelt |= D < 0;
	// allow for radiation field to grow/shrink with solar activity
	D -= activity * 0.25 / mf.inner_radius;
	r = RadiationInBelt(D, mf.inner_radius, rb.radiation_inner_gradient);
	vd.radiationInnerBelt += r * rb.radiation_inner * (1 + activity * 0.3);
	}
	if (mf.has_outer)
	{
	D = mf.OuterFunc(p);
	outerBelt |= D < 0;
	// allow for radiation field to grow/shrink with solar activity
	D -= activity * 0.25 / mf.outer_radius;
	r = RadiationInBelt(D, mf.outer_radius, rb.radiation_outer_gradient);
	vd.radiationOuterBelt += r * rb.radiation_outer * (1 + activity * 0.3);
	}
	if (mf.has_pause)
	{
	gsm = GsmSpace(rb, false);
	p = gsm.TransformIn(scaled_position);
	D = mf.PauseFunc(p);
	vd.radiationMagnetopause += Lib.Clamp(D / -0.1332f, 0.0f, 1.0f) * rb.RadiationPause();
	magnetosphere |= D < 0.0f && !Sim.IsStar(rb.body); //< ignore heliopause
	interstellar |= D > 0.0f && Sim.IsStar(rb.body); //< outside heliopause
	}
	}
	if (rb.radiation_surface > 0)
	{
	if (Sim.IsBodyVisible(vd.Vessel, position, body, vd.VisibleBodies, out _, out double distance))
	{
	var r0 = RadiationR0(rb);
	var r1 = DistanceRadiation(r0, distance);
	// clamp to max. surface radiation.
	// when loading on a rescaled system, the vessel can appear to be within the sun for a few ticks
	var rad = Math.Min(r1, rb.radiation_surface);
	if (Sim.IsStar(body))
	{
	vd.radiationSolar += rad;
	}
	else
	{
	if (body == vd.MainBody)
	mainBodyRadiation += rad; // TODO : this should be affected by the body gamma transparency, but the other way around
	else
	vd.radiationBodies += rad;
	}
	}
	}
	// avoid loops in the chain
	body = (body.referenceBody != null && body.referenceBody.referenceBody == body) ? null : body.referenceBody;
	}
	// add all sources and apply gamma transparency if inside atmosphere
	vd.radiationInnerBelt *= gammaTransparency;
	vd.radiationOuterBelt *= gammaTransparency;
	vd.radiationMagnetopause *= gammaTransparency;
	vd.radiationSolar *= gammaTransparency;
	vd.radiationBodies *= gammaTransparency;
	radiation += vd.radiationInnerBelt;
	radiation += vd.radiationOuterBelt;
	radiation += vd.radiationMagnetopause;
	radiation += vd.radiationSolar;
	radiation += vd.radiationBodies;
	// add extern radiation
	radiation += Settings.ExternRadiation * gammaTransparency;
	// main body surface radiation is unaffected by gamma transparency
	radiation += mainBodyRadiation;
	vd.radiationBodies += mainBodyRadiation;
	#if DEBUG_RADIATION
	if (v.loaded) Lib.Log("Radiation " + v + " after gamma: " + Lib.HumanReadableRadiation(radiation) + " transparency: " + gamma_transparency);
	#endif
	//var passiveShielding = PassiveShield.Total(v);
	//shieldedRadiation -= passiveShielding;
	// clamp radiation to positive range
	// note: we avoid radiation going to zero by using a small positive value
	radiation = Math.Max(radiation, 0.0);
	#if DEBUG_RADIATION
	if (v.loaded) Lib.Log("Radiation " + v + " after clamp: " + Lib.HumanReadableRadiation(radiation) + " shielded " + Lib.HumanReadableRadiation(shieldedRadiation));
	#endif
	// return radiation
	return radiation;
	}

You will notice that for each of the above variable, there is a matching boolean. This is because we used to only output the "total" aggregate radiation, and booleans to identify what sources were "active". I quick-and-dirty changed that recently so we have access to the detail rate of each source, but I was also in the process of refactoring the whole thing so EnvRadiationSolar is tracked per sun, and EnvRadiationBodies per body (because they are directional, and I was planning to take that into account, at least for solar radiation). Said otherwise, I intended to move EnvRadiationSolar to the StarFlux structure. But since you already are in the process of extending/refactoring StarFlux to a BodyFlux structure (and processing all bodies internally), we should likely do both from the ground up.

In terms of re-implementation, a quick overview :

You will need to make a "job-safe" copy of the RadiationBody and RadiationModel data structures. All values in those classes are constants, so you only need to make that once. RadiationModel are essentially "templates", so I suggest you aggregate the data in a single "per body" data structure (and in fact, you should use the same data structure for info about stars too, instead of re-getting it every iteration. This is also constant data.)

Belts and magnetopause radiation are computed by evaluating the vessel position against the fields shape, which are represented by math functions. The called methods are a bit all over the place for no reason other than technical debt, but you should be able to rebuild the execution path quite easily by following the execution of Radiation.Compute() (it's entirely self contained, to my knowledge there is nothing modifying the data that this depends on).

Solar/bodies radiation are directional, and affected by other bodies occlusion (and also, they are the same value in RadiationModel). From a physics POV, they are identical to the solar irradiance, so you can likely just merge that into your existing code.

Storm radiation is a bit trickier. As mentioned in #5, I intended to more or less entirely scrap the current per-vessel-random implementation, possibly in favor of generating "cone" storms, where every vessel/body inside the cone is subject to the storm. For now, I would suggest to postpone including storm radiation in the substep sim until this is done.

[darco89]

beloved modders, i just found out about kerbalism 4. are you still working on it? Cheers!