Modifying one of the propagators handled by parallelizer causes code to stall

Emiliano · May 3, 2023, 3:13pm

Hello everyone,
As I mentioned in this thread, I’m running into a weird issue.
I’m trying to propagate a set of satellites with a parallelizer, and while doing so I also want to add maneuvers to the satellites and see how they evolve relative to each other. However there is something weird happening. If I run a first “simple” propagation (without adding any firing to any satellite) the parallelizer works fine and it returns the results I’m expecting. However, if I later add a firing to one of the propagators and run the propagate method at parallelizer level, the propagation stalls. It doesn’t raise any error, nothing. It just doesn’t do anything. The console appears busy, but the propagation doesn’t start (or at least that what I think). In particular, I tried adding a print statement to my MultiSatStepHandler so as to verify whether something is getting done. While something gets printed during the first “simple” propagation, nothing gets printed after calling the second propagation (after adding the firing).
Schematically, the propagation is done this way:

Create the propagators and the handler
Create the parallelizer
Propagate for 1 day
Add a firing
Propagate for 1 day

The last step is the one that doesn’t run.
To troubleshoot a bit (thinking that maybe modifying the propagators after feeding them to the parallelizer wasn’t allowed) I tried to use this scheme:

Create the propagators and the handler
Create the parallelizer
Add a firing
Propagate

This setup works. The propagation runs and returns the expected results.

Best regards,
Emiliano

Emiliano · May 3, 2023, 3:23pm

Here is a simple version of my code that runs into this issue:

import orekit;
orekit.initVM();

from orekit.pyhelpers import setup_orekit_curdir;
setup_orekit_curdir();

from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator;

from org.orekit.attitudes import BodyCenterPointing;
from org.orekit.bodies import CelestialBodyFactory, OneAxisEllipsoid;
from org.orekit.forces.drag import DragForce, IsotropicDrag;
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel, OceanTides, SolidTides, ThirdBodyAttraction;
from org.orekit.forces.gravity.potential import GravityFieldFactory;
from org.orekit.forces.maneuvers import ConstantThrustManeuver;
from org.orekit.forces.radiation import SolarRadiationPressure, IsotropicRadiationSingleCoefficient;
from org.orekit.frames import FramesFactory;
from org.orekit.models.earth.atmosphere import HarrisPriester;
from org.orekit.orbits import CircularOrbit, PositionAngle;
from org.orekit.propagation import PropagatorsParallelizer, SpacecraftState;
from org.orekit.propagation.numerical import NumericalPropagator;
from org.orekit.propagation.semianalytical.dsst.forces import DSSTAtmosphericDrag, DSSTSolarRadiationPressure, DSSTTesseral, DSSTThirdBody, DSSTZonal;
from org.orekit.time import AbsoluteDate, TimeScalesFactory;
from org.orekit.utils import Constants, IERSConventions;

from globalStepHandlers import GlobalCustomStepHandler;

from math import pi;
from numpy import degrees;
from java.util import Arrays;

# %% INITIAL STATES
# =============================================================================
# =============================================================================
Re = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
mu = Constants.WGS84_EARTH_MU;
earth = OneAxisEllipsoid(Re,
                          Constants.WGS84_EARTH_FLATTENING,
                          FramesFactory.getITRF(IERSConventions.IERS_2010, True));
earthCenterAttitudeLaw0 = BodyCenterPointing(FramesFactory.getEME2000(), earth);
earthCenterAttitudeLaw1 = BodyCenterPointing(FramesFactory.getEME2000(), earth);

mass = 150.;
area = 1.;
CD = 2.2;
CR = 1.2;

a0 = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 500e3;
a1 = a0 + 5e3;
ex = 1e-4;
ey = 1e-4;
i = 1.5;
raan = pi/6;
u = pi/18;

orbit0 = CircularOrbit(a0, ex, ey, i, raan, u, PositionAngle.MEAN,
                        FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH, Constants.WGS84_EARTH_MU);
initState0 = SpacecraftState(orbit0, mass);
orbit1 = CircularOrbit(a1, ex, ey, i, raan, u, PositionAngle.MEAN,
                        FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH, Constants.WGS84_EARTH_MU);
initState1 = SpacecraftState(orbit1, mass);

# %% DYNAMICAL MODEL
# =============================================================================
# =============================================================================
n, m = 4, 4;
nonSphericalGravity = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
                                                        GravityFieldFactory.getConstantNormalizedProvider(n,m));        
sunThirdBody = ThirdBodyAttraction(CelestialBodyFactory.getSun());
moonThirdBody = ThirdBodyAttraction(CelestialBodyFactory.getMoon());
dragForce = DragForce(HarrisPriester(CelestialBodyFactory.getSun(), earth), IsotropicDrag(area, CD));
srpForce = SolarRadiationPressure(CelestialBodyFactory.getSun(),
                                  Re, IsotropicRadiationSingleCoefficient(area, CR));
oceanicTides = OceanTides(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
                          Re, mu, n, m, IERSConventions.IERS_2010, 
                          TimeScalesFactory.getUT1(IERSConventions.IERS_2010, True));
solidTides = SolidTides(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
                        Re, mu, GravityFieldFactory.getConstantNormalizedProvider(n, m).getTideSystem(),
                        IERSConventions.IERS_2010, TimeScalesFactory.getUT1(IERSConventions.IERS_2010, True),
                        [CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon()]);
forces = [nonSphericalGravity, sunThirdBody, moonThirdBody, dragForce, srpForce, oceanicTides, solidTides];

zonalDsst = DSSTZonal(GravityFieldFactory.getConstantUnnormalizedProvider(n,m));
tesseralDsst = DSSTTesseral(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
                            Constants.WGS84_EARTH_ANGULAR_VELOCITY,
                            GravityFieldFactory.getConstantUnnormalizedProvider(n,m));
sunThirdBodyDsst = DSSTThirdBody(CelestialBodyFactory.getSun(), mu);
moonThirdBodyDsst = DSSTThirdBody(CelestialBodyFactory.getMoon(), mu);
dragForceDsst = DSSTAtmosphericDrag(HarrisPriester(CelestialBodyFactory.getSun(), earth),
                                    IsotropicDrag(area, CD), mu);
srpForceDsst = DSSTSolarRadiationPressure(CelestialBodyFactory.getSun(), Re,
                                          IsotropicRadiationSingleCoefficient(area, CR), mu);
forcesDsst = [zonalDsst, tesseralDsst, sunThirdBodyDsst, moonThirdBodyDsst, dragForceDsst, srpForceDsst];

# %% PROPAGATORS
# =============================================================================
# =============================================================================
minStep = 1e-2;
maxStep = 1e+2;
absTol = 1e-9;
relTol = 1e-9;

integrator0 = DormandPrince853Integrator(minStep, maxStep, absTol, relTol);
integrator0.setInitialStepSize(60.0);

integrator1 = DormandPrince853Integrator(minStep, maxStep, absTol, relTol);
integrator1.setInitialStepSize(60.0);

propagator0 = NumericalPropagator(integrator0);
propagator0.setAttitudeProvider(earthCenterAttitudeLaw0);
for force in forces:
    propagator0.addForceModel(force);
propagator0.setInitialState(initState0);

propagator1 = NumericalPropagator(integrator1);
propagator1.setAttitudeProvider(earthCenterAttitudeLaw1);
for force in forces:
    propagator1.addForceModel(force);
propagator1.setInitialState(initState1);

myHandler = GlobalCustomStepHandler();
parallelizer = PropagatorsParallelizer(Arrays.asList([propagator0, propagator1]), myHandler);

day = 24*60*60.;
initDate = orbit0.getDate();
finalDate = initDate.shiftedBy(day);

parallelizer.propagate(initDate, finalDate);

print(degrees(myHandler.getRaanDiff()));

dateNow = finalDate;
firingDate = dateNow.shiftedBy(day/3);
firing = ConstantThrustManeuver(firingDate, 1e3, 1e-1, 1e3, Vector3D.PLUS_I);
propagator1.addForceModel(firing);

parallelizer.propagate(dateNow, dateNow.shiftedBy(day));

print(degrees(myHandler.getRaanDiff()));

This is my step handler:

from org.orekit.orbits import CircularOrbit;
from org.orekit.propagation.sampling import PythonMultiSatStepHandler, OrekitStepInterpolator;
from org.orekit.propagation.semianalytical.dsst import DSSTPropagator;

class GlobalCustomStepHandler(PythonMultiSatStepHandler):
    def init(self, states0, t):
        pass;
        
    def handleStep(self, interpolators):
        raanDiff = [];
        interpolators = list(interpolators);

        stateChief = OrekitStepInterpolator.cast_(interpolators[0]).getCurrentState();
        orbitChief = CircularOrbit(stateChief.getOrbit());
        raanChief = orbitChief.getRightAscensionOfAscendingNode();
        
        for ii in range(1, len(interpolators)):
            stateNow = OrekitStepInterpolator.cast_(interpolators[ii]).getCurrentState();
            orbitNow = CircularOrbit(stateNow.getOrbit());
            raanDiff.append(orbitNow.getRightAscensionOfAscendingNode() - raanChief);
        print(raanDiff);
        self.raanDiff = raanDiff;
        
    def getRaanDiff(self) -> list:
        return self.raanDiff;
        
    def finish(self, finalStates):
        pass;

MaximeJ · June 6, 2023, 3:15pm

Hi @Emiliano,

Sorry it looks like we didn’t answer your post…
Are you still experimenting this issue ?

Emiliano · June 7, 2023, 8:10am

Hello @MaximeJ,

Yes, unfortunately I still run into this issue

MaximeJ · June 7, 2023, 12:18pm

Hello @Emiliano,

Again, sorry for the delay.

I’ve reproduced your issue in Java, and it looks like a bug to me.
Could you please open an issue for this on the forge?

The behavior is a bit erratic and sometimes even without the maneuver it stalls.
I’ve used a profiler and here is what it says:

I don’t know exactly what’s happening and why a deadlock occurs…

I’ve made it work though, but it implies re-instantiating both propagators and the parallelizer.
See a Java working example (TestParallelizerAndManeuver.java). I’ve simplified your code but it should be easy to reproduce with Python.
If you comment on the working part this should stall again, even without the maneuver…
This one will be difficult to fix because multi-threading is hard to debug.

Emiliano · June 12, 2023, 7:42am

Hello @MaximeJ,

Thank you for your insights, I opened an issue on the forge as you asked me. Thank you for your help!

Best regards,
Emiliano

MaximeJ · June 12, 2023, 1:21pm

Hello @Emiliano,

Good news is I think I found what the problem was, and how to fix it
See merge request 366.

Best regards,
Maxime