Merge pull request #56 from rahil-makadia/dev

safer parallel programming with ephemeris memory maps

grss/version.txt

@@ -1 +1 @@
-3.8.2
+3.8.3

include/grss.h

@@ -13,6 +13,7 @@ void PropSimulation::integrate() {
     }
 
     // integrate the system
+    this->map_ephemeris();
     this->preprocess();
     if (!this->parallelMode) furnsh_c(this->DEkernelPath.c_str());
     gr15(this);

include/simulation.h

@@ -10,7 +10,7 @@ class PropSimulation;
  * @brief Get the name of body-fixed frame for a given SPICE ID.
  */
 void get_baseBodyFrame(const int &spiceId, const real &tMjdTDB,
-                       ConstSpiceChar *&baseBodyFrame);
+                       std::string &baseBodyFrame);
 
 /**
  * @brief Get the observer state for a given time.

include/spk.h

@@ -3,14 +3,7 @@
 #ifndef SPK_H
 #define SPK_H
 
-#include <fcntl.h>
-#include <sys/mman.h>
-#include <sys/stat.h>
-#include <unistd.h>
-#include <cstring>
-#include <iostream>
-#include <vector>
-#include <stdexcept>
+#include "utilities.h"
 
 /**
  * @brief Structure to hold a single time,pos,vel,accel record from an SPK file.
@@ -157,4 +150,10 @@ void spk_calc(DafInfo *pl, double epoch, int spiceId, double *out_x,
 void get_spk_state(const int &spiceId, const double &t0_mjd, Ephemeris &ephem,
                    double state[9]);
 
+/**
+ * @brief Using cspice to get the rotation matrix from one frame to another.
+ */
+void get_pck_rotMat(const std::string &from, const std::string &to,
+                    const real &et, std::vector<std::vector<real>> &rotMat);
+
 #endif

include/timeconvert.h

@@ -1,7 +1,7 @@
 #ifndef TIME_H
 #define TIME_H
 
-#include "utilities.h"
+#include "spk.h"
 
 /**
  * @brief Convert Julian Date to TDB ephemeris time.

include/utilities.h

@@ -1,6 +1,16 @@
 #ifndef UTILITIES_H
 #define UTILITIES_H
 
+#include <fcntl.h>
+#include <sys/mman.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include "SpiceUsr.h"
+
+#include <cstring>
+#include <iostream>
+#include <vector>
+#include <stdexcept>
 #include <algorithm>
 #include <cmath>
 #include <limits>
@@ -22,9 +32,6 @@ using std::sqrt;
 using std::tan;
 using std::tanh;
 
-#include "SpiceUsr.h"
-#include "spk.h"
-
 /**
  * @brief Real type to be used for floating point calculations.
  * Choose between double and long double.

src/approach.cpp

@@ -81,9 +81,11 @@ void check_ca_or_impact(PropSimulation *propSim, const real &tOld,
                     if (!propSim->parallelMode) impact.get_impact_parameters(propSim);
                     impact.impact = true;
                     propSim->impactParams.push_back(impact);
-                    std::cout << "Impact detected at MJD " << tImp << " TDB. "
-                              << propSim->integBodies[i].name
-                              << " collided with " << bodyj->name << "!" << std::endl;
+                    if (!propSim->parallelMode) {
+                        std::cout << "Impact detected at MJD " << tImp << " TDB. "
+                                << propSim->integBodies[i].name
+                                << " collided with " << bodyj->name << "!" << std::endl;
+                    }
                 }
                 // check close approach
                 real radialVel, radialVelOld;
@@ -632,18 +634,12 @@ void CloseApproachParameters::print_summary(int prec){
  * @param[in] propSim PropSimulation object.
  */
 void ImpactParameters::get_impact_parameters(PropSimulation *propSim){
-    ConstSpiceChar *baseBodyFrame;
+    std::string baseBodyFrame;
     get_baseBodyFrame(this->centralBodySpiceId, this->t, baseBodyFrame);
     real tET;
     mjd_to_et(this->t, tET);
-    SpiceDouble rotMatSpice[6][6];
-    sxform_c("J2000", baseBodyFrame, tET, rotMatSpice);
     std::vector<std::vector<real>> rotMat(6, std::vector<real>(6));
-    for (size_t i = 0; i < 6; i++) {
-        for (size_t j = 0; j < 6; j++) {
-            rotMat[i][j] = (real)rotMatSpice[i][j];
-        }
-    }
+    get_pck_rotMat("J2000", baseBodyFrame, tET, rotMat);
     std::vector<real> impactRelStateInertial(6);
     impactRelStateInertial[0] = this->xRel[0]*propSim->consts.du2m/1.0e3L;
     impactRelStateInertial[1] = this->xRel[1]*propSim->consts.du2m/1.0e3L;

src/parallel.cpp

@@ -12,6 +12,13 @@ std::vector<PropSimulation> propSim_parallel_omp(
     size_t numBodies = allBodies.size();
     std::vector<PropSimulation> allSims(numBodies, refSim);
 
+    // make sure refSim is in parallel mode
+    if (!refSim.parallelMode) {
+        throw std::runtime_error(
+            "ERROR: The reference simulation must be in parallel mode to "
+            "propagate in parallel.");
+    }
+
     // parallel for loop to first create an integBody for each entry in the
     // allBodies vector, then integrate each integBody using the reference
     // simulation
@@ -48,5 +55,13 @@ std::vector<PropSimulation> propSim_parallel_omp(
             allSims[i] = sim;
         }
     }
+    // in serial, compute the body-fixed impact coordinates (lat,lon,alt) for each instance
+    furnsh_c(refSim.DEkernelPath.c_str());
+    for (size_t i = 0; i < allSims.size(); i++) {
+        for (size_t j = 0; j < allSims[i].impactParams.size(); j++) {
+            allSims[i].impactParams[j].get_impact_parameters(&(allSims[i]));
+        }
+    }
+    unload_c(refSim.DEkernelPath.c_str());
     return allSims;
 }

src/simulation.cpp

@@ -6,7 +6,7 @@
  * @param[out] baseBodyFrame Name of the body-fixed frame.
  */
 void get_baseBodyFrame(const int &spiceId, const real &tMjdTDB,
-                       ConstSpiceChar *&baseBodyFrame) {
+                       std::string &baseBodyFrame) {
     switch (spiceId) {
         case 10:
             baseBodyFrame = "IAU_SUN";
@@ -63,7 +63,7 @@ void get_observer_state(const real &tObsMjd,
                         const std::vector<real> &observerInfo,
                         PropSimulation *propSim, const bool tObsInUTC,
                         std::vector<real> &observerState) {
-    SpiceInt baseBody = observerInfo[0];
+    int baseBody = observerInfo[0];
     if ((int)observerInfo[0] == 500) baseBody = 399;
     if (baseBody == 0) {
         observerState[0] = 0.0L;
@@ -94,16 +94,10 @@ void get_observer_state(const real &tObsMjd,
         observerState[5] = (real) baseBodyState[5] + observerInfo[6]/propSim->consts.duptu2mps;
         return;
     }
-    ConstSpiceChar *baseBodyFrame;
+    std::string baseBodyFrame;
     get_baseBodyFrame((int)observerInfo[0], tObsMjdTDB, baseBodyFrame);
-    SpiceDouble rotMatSpice[6][6];
-    sxform_c(baseBodyFrame, "J2000", secPastJ2000, rotMatSpice);
     std::vector<std::vector<real>> rotMat(6, std::vector<real>(6));
-    for (size_t i = 0; i < 6; i++) {
-        for (size_t j = 0; j < 6; j++) {
-            rotMat[i][j] = (real)rotMatSpice[i][j];
-        }
-    }
+    get_pck_rotMat(baseBodyFrame, "J2000", secPastJ2000, rotMat);
     real lon = observerInfo[1];
     real lat = observerInfo[2];
     real rho = observerInfo[3];
@@ -680,14 +674,12 @@ PropSimulation::PropSimulation(std::string name, real t0,
         default:
             throw std::invalid_argument(
                 "The defaultSpiceBodies argument is only defined for no "
-                "default "
-                "SPICE bodies (case 0) or DE431 (case 431) or DE441 (case "
-                "441).");
+                "preloaded SPICE bodies (case 0) or DE431 (case 431) or DE441 "
+                "(case 441).");
             break;
     }
     this->ephem.mbPath = kernel_mb;
     this->ephem.sbPath = kernel_sb;
-    this->map_ephemeris();
 }
 
 /**
@@ -797,6 +789,7 @@ void PropSimulation::prepare_for_evaluation(
         tEval.size(), std::vector<real>(6, 0.0L));
     std::vector<int> radarObserver = std::vector<int>(tEval.size(), 0);
     if (!this->parallelMode) {
+        this->map_ephemeris();
         furnsh_c(this->DEkernelPath.c_str());
         if (tEval.size() != 0) {
             for (size_t i = 0; i < tEval.size(); i++) {
@@ -813,12 +806,10 @@ void PropSimulation::prepare_for_evaluation(
                 }
                 get_observer_state(tEval[i], observerInfo[i], this,
                                 this->tEvalUTC, xObserver[i]);
-                // std::cout << xObserver[i][0] << " " << xObserver[i][1] << " " <<
-                // xObserver[i][2] << " " << xObserver[i][3] << " " <<
-                // xObserver[i][4] << " " << xObserver[i][5] << std::endl;
             }
         }
         unload_c(this->DEkernelPath.c_str());
+        this->unmap_ephemeris();
     }
 
     if (this->tEval.size() == 0) {
@@ -866,8 +857,8 @@ std::vector<real> PropSimulation::get_spiceBody_state(const real t, const std::s
     }
     if (this->ephem.mb == nullptr || this->ephem.sb == nullptr){
         throw std::invalid_argument(
-            "Ephemeris kernels are not loaded. Memory map the ephemeris using "
-            "PropSimulation.map_ephemeris() method first.");
+            "get_spiceBody_state: Ephemeris kernels are not loaded. Memory map "
+            "the ephemeris using PropSimulation.map_ephemeris() method first.");
     }
     double spiceState[9];
     get_spk_state(spiceId, t, this->ephem, spiceState);
@@ -913,18 +904,6 @@ void PropSimulation::add_integ_body(IntegBody body) {
             " TDB which is different from the simulation initial time: MJD " +
             std::to_string(this->integParams.t0) + " TDB.");
     }
-    if (body.isCometary) {
-        double sunState[9];
-        get_spk_state(10, body.t0, this->ephem, sunState);
-        body.pos[0] += sunState[0];
-        body.pos[1] += sunState[1];
-        body.pos[2] += sunState[2];
-        body.vel[0] += sunState[3];
-        body.vel[1] += sunState[4];
-        body.vel[2] += sunState[5];
-    }
-    body.initState = {body.pos[0], body.pos[1], body.pos[2],
-                      body.vel[0], body.vel[1], body.vel[2]};
     body.radius /= this->consts.du2m;
     this->integBodies.push_back(body);
     this->integParams.nInteg++;
@@ -1098,15 +1077,25 @@ void PropSimulation::preprocess() {
     if (!this->isPreprocessed) {
         this->t = this->integParams.t0;
         for (size_t i = 0; i < this->integParams.nInteg; i++) {
+            if (this->integBodies[i].isCometary) {
+                double sunState[9];
+                get_spk_state(10, this->integBodies[i].t0, this->ephem, sunState);
+                this->integBodies[i].pos[0] += sunState[0];
+                this->integBodies[i].pos[1] += sunState[1];
+                this->integBodies[i].pos[2] += sunState[2];
+                this->integBodies[i].vel[0] += sunState[3];
+                this->integBodies[i].vel[1] += sunState[4];
+                this->integBodies[i].vel[2] += sunState[5];
+            }
             for (size_t j = 0; j < 3; j++) {
-                this->xInteg.push_back(integBodies[i].pos[j]);
+                this->xInteg.push_back(this->integBodies[i].pos[j]);
             }
             for (size_t j = 0; j < 3; j++) {
-                this->xInteg.push_back(integBodies[i].vel[j]);
+                this->xInteg.push_back(this->integBodies[i].vel[j]);
             }
-            if (integBodies[i].propStm) {
-                for (size_t j = 0; j < integBodies[i].stm.size(); j++) {
-                    this->xInteg.push_back(integBodies[i].stm[j]);
+            if (this->integBodies[i].propStm) {
+                for (size_t j = 0; j < this->integBodies[i].stm.size(); j++) {
+                    this->xInteg.push_back(this->integBodies[i].stm[j]);
                 }
             }
         }
@@ -1150,9 +1139,6 @@ void PropSimulation::extend(real tf, std::vector<real> tEvalNew,
     this->radarObs.clear();
     this->radarPartials.clear();
 
-    // map the ephemeris again
-    this->map_ephemeris();
-
     // first prepare for integration and then integrate
     this->integParams.t0 = this->t;
     this->interpParams.tStack.push_back(this->t);

src/spk.cpp

@@ -109,7 +109,7 @@ DafInfo* daf_init(const std::string &path, const std::string &type) {
         close(fd);
         return NULL;
     }
-    // std::cout << "record.summary.nsum: " << record.summary.nsum << std::endl;
+    // std::cout << "record.summaries.nsum: " << record.summaries.nsum << std::endl;
     // std::cout << "record.file.nd: " << record.file.nd << std::endl;
     // std::cout << "record.file.ni: " << record.file.ni << std::endl;
     // std::cout << "nc: " << nc << std::endl;
@@ -282,13 +282,11 @@ void spk_calc(DafInfo *pl, double epoch, int spiceId, double *out_x,
     U[2] = 4.0;
     for (p = 2; p < P; p++) {
         T[p] = 2.0 * z * T[p - 1] - T[p - 2];
-        // Not used at the moment:
         S[p] = 2.0 * z * S[p - 1] + 2.0 * T[p - 1] - S[p - 2];
     }
     for (p = 3; p < P; p++) {
         U[p] = 2.0 * z * U[p - 1] + 4.0 * S[p - 1] - U[p - 2];
     }
-    // double c = (0.125 / 2 / 86400.0);
     double t1 = 32.0;
     int niv;
     switch (spiceId) {
@@ -365,6 +363,11 @@ void spk_calc(DafInfo *pl, double epoch, int spiceId, double *out_x,
  */
 void get_spk_state(const int &spiceId, const double &t0_mjd, Ephemeris &ephem,
                    double state[9]) {
+    if (ephem.mb == nullptr || ephem.sb == nullptr){
+        throw std::invalid_argument(
+            "get_spk_state: Ephemeris kernels are not loaded. Memory map "
+            "the ephemeris using PropSimulation.map_ephemeris() method first.");
+    }
     bool smallBody = spiceId > 1000000;
     DafInfo *infoToUse;
     if (smallBody) {
@@ -461,3 +464,20 @@ void get_spk_state(const int &spiceId, const double &t0_mjd, Ephemeris &ephem,
     ephem.cache[ephem.nextIdxToWrite].items[cacheIdx].ay = state[7];
     ephem.cache[ephem.nextIdxToWrite].items[cacheIdx].az = state[8];
 }
+
+/**
+ * @param from Frame to rotate from.
+ * @param to Frame to rotate to.
+ * @param et TDB ephemeris time in seconds past J2000 to get the rotation matrix at.
+ * @param rotMat Rotation matrix from 'from' to 'to'.
+ */
+void get_pck_rotMat(const std::string &from, const std::string &to,
+                    const real &et, std::vector<std::vector<real>> &rotMat) {
+    SpiceDouble rotMatSpice[6][6];
+    sxform_c(from.c_str(), to.c_str(), et, rotMatSpice);
+    for (size_t i = 0; i < 6; i++) {
+        for (size_t j = 0; j < 6; j++) {
+            rotMat[i][j] = (real)rotMatSpice[i][j];
+        }
+    }
+}