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

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:

  1. Create the propagators and the handler
  2. Create the parallelizer
  3. Propagate for 1 day
  4. Add a firing
  5. 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:

  1. Create the propagators and the handler
  2. Create the parallelizer
  3. Add a firing
  4. Propagate

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

Best regards,

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

import orekit;

from orekit.pyhelpers import 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 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;

# =============================================================================
# =============================================================================
mu = Constants.WGS84_EARTH_MU;
earth = OneAxisEllipsoid(Re,
                          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);

# =============================================================================
# =============================================================================
n, m = 4, 4;
nonSphericalGravity = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
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),
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];

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

integrator0 = DormandPrince853Integrator(minStep, maxStep, absTol, relTol);

integrator1 = DormandPrince853Integrator(minStep, maxStep, absTol, relTol);

propagator0 = NumericalPropagator(integrator0);
for force in forces:

propagator1 = NumericalPropagator(integrator1);
for force in forces:

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

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

parallelizer.propagate(initDate, finalDate);


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

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


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):
    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);
        self.raanDiff = raanDiff;
    def getRaanDiff(self) -> list:
        return self.raanDiff;
    def finish(self, finalStates):

Hi @Emiliano,

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

Hello @MaximeJ,

Yes, unfortunately I still run into this issue

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 ( 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… :frowning:
This one will be difficult to fix because multi-threading is hard to debug.

Hello @MaximeJ,

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

Best regards,

1 Like

Hello @Emiliano,

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

Best regards,