Multiprocessing Example¶
In this example we will use Python Multiprocessing to perform several propagations in parallel.
First perform all necessary imports:
%matplotlib inline
from godot.core import tempo, util
from godot.model import common, geometry
from godot import cosmos
import numpy as np
import multiprocessing as mp
from dataclasses import dataclass
# avoid verbose output
util.suppressLogger()
Now we will prepare a function to be called in several processess.
First define a special output dataclass for this function:
@dataclass
class PropagationCase:
id: str
epoch0:tempo.Epoch
state0: np.ndarray
area: float
mass: float
cd: float
epoch1:tempo.Epoch
state1: np.ndarray
Our function takes the universe configuration as input and constructs a new universe for each process.
Then it constructs a Ballistic Propagator and adds an interrupt event in case the spacecraft altitude drops below a certain value.
Note that in the universe file we created several constants sc_mass, sc_area and sc_cd with nominally zero values. These values are overwritten in each local process with values passed as arguments.
def prop(uni_cfg, max_duration, dyn, case):
print(f"Performing Propagation {case.id}")
# ccreate universe
uni = cosmos.Universe(uni_cfg)
# reset values of sc parameters
uni.evaluables.get('sc_area').set(case.area)
uni.evaluables.get('sc_mass').set(case.mass)
uni.evaluables.get('sc_cd').set(case.cd)
# reference date for propagation
ref = tempo.Epoch('2020-01-01T00:00:00 TDB')
# construct propagator
pro = cosmos.BallisticPropagator(
uni,
'sc',
dyn,
ref,
"Earth",
1e-9
)
# add event for atmosphere entry at 6538 km radius (~160 km radius)
pro.addEvent(
common.InterruptEvent(
common.ValueTrigger(
geometry.Range(uni.frames, 'Earth', 'sc'),
common.ConstantScalarTimeEvaluable(6538))))
# perform propagation
pro.compute(
case.state0,
case.epoch0 - ref,
case.epoch0 + max_duration)
# return final epoch and state in icrf and itrf with input values
case.epoch1 = pro.range().end()
case.state1 = uni.frames.vector6('Earth', 'sc', 'ICRF', pro.range().end())
return case
We read the universe file, set the dynamics to use and then define the data for a number of cases to propagate:
# read universe file
config_uni = cosmos.util.load_yaml('universe.yml')
# define dynamics to use
dyn = 'dynamics'
# define some cases to propagate
e0 = tempo.Epoch('2020-01-12T00:00:00.00 TDB')
x0 = np.array([7200.0, 200.0, 35.0, 0.0, 2.0, 7.8])
cd0 = 1.5
mass0 = 10.0
area0 = 1.0
cases = []
for i in range(200):
x = x0
x[:3] += np.random.normal(0.0, 1e-1, (3))
x[3:] += np.random.normal(0.0, 1e-4, (3))
cd = np.maximum(1.0, cd0 + np.random.normal(0.0, 1.0))
mass = np.maximum(1.0, mass0 + np.random.normal(0.0, 5.0))
area = np.maximum(0.1, area0 + np.random.normal(0.0, 0.5))
cases.append(
PropagationCase(
i, e0, x, area, mass, cd, None, None
)
)
# set the maximum duration to 4 days
max_duration = 4*tempo.SecondsInDay
Now we use the Python multiprocessing package to propagate the cases 4 at a time (assuming you have the resouces available):
npool = 4
pool = mp.Pool(processes=npool)
output = [
pool.apply_async(
prop,
args=(config_uni, max_duration, dyn, case)
)
for case in cases
]
Finally we can print out the results
results = [p.get() for p in output]
for res in results:
print(f"{res.id} {res.epoch1} {res.state1[:3]}")
Plot a histogram of orbital lifetime for the cases.
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams.update({'figure.figsize':(7,5), 'figure.dpi':100})
# Plot Histogram on x
x = [(res.epoch1 - res.epoch0)/tempo.SecondsInDay for res in results]
plt.hist(x, bins=50)
plt.xlabel('Days')
plt.gca().set(title='Orbital Lifetime', ylabel='Frequency');
