Merge branch 'main' of https://github.com/VirtualPlanetaryLaboratory/…

…vplanet
VirtualPlanetaryLaboratory · Feb 6, 2024 · 1fd0b47 · 1fd0b47
2 parents fff87ae + 9c4b5d5
commit 1fd0b47
Show file tree

Hide file tree

Showing 62 changed files with 3,930 additions and 1,987 deletions.
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
@@ -2,9 +2,7 @@ name: docs
 
 on:
   push:
-    branches: [master, main]
-  #pull_request: TODO: Enable me
-  #  branches:[master, main]
+    branches: [main]
 
 jobs:
   tests:

diff --git a/.github/workflows/pip-install.yml b/.github/workflows/pip-install.yml
@@ -1,8 +1,7 @@
 name: pip-install
 on:
   release:
-    types:
-      - published
+    types: [published]
 
 jobs:
   build_wheels:
@@ -11,7 +10,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-20.04, ubuntu-22.04, acos-11, macos-12, macos-13]
-        python: ['3.7', '3.8', '3.9', '3.10', '3.11']
+        python: ['3.6', '3.7', '3.8', '3.9', '3.10', '3.11']
         exclude:
           - os: ubuntu-22.04
             python-version: '3.6'

diff --git a/.github/workflows/tests-linux.yml b/.github/workflows/tests-linux.yml
@@ -1,8 +1,6 @@
 name: tests-linux
 
 on:
-  push:
-    branches: [main]
   pull_request:
     branches: [main]
   workflow_dispatch:

diff --git a/.github/workflows/tests-macos.yml b/.github/workflows/tests-macos.yml
@@ -1,7 +1,5 @@
 name: tests-macos
 on:
-  push:
-    branches: [main]
   pull_request:
     branches: [main]
   workflow_dispatch:
@@ -14,7 +12,7 @@ jobs:
       matrix:
         os: [macos-11, macos-12, macos-13]
         python-version: ['3.6', '3.7', '3.8', '3.9', '3.10', '3.11']
-          
+
       fail-fast: false
     steps:
       - uses: actions/checkout@v4
@@ -41,4 +39,3 @@ jobs:
       - name: Get unique id
         uses: Tiryoh/gha-jobid-action@v1
         id: jobs
-
diff --git a/docs/conf.py b/docs/conf.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
-# planetplanet documentation build configuration file, created by
+# vplanet documentation build configuration file, created by
 # sphinx-quickstart on Wed Aug  9 19:49:05 2017.
 #
 # This file is execfile()d with the current directory set to its

diff --git a/environment.yml b/environment.yml
@@ -15,6 +15,7 @@ dependencies:
       - pytest-env>=0.6.2
       - pytest-cov>=2.11.1
       - sphinx
+      - sphinx-rtd-theme
       - breathe>=4.30.0
       - sphinxcontrib-bibtex>=2.2.0
       - coverage>=5.5

diff --git a/src/atmesc.c b/src/atmesc.c
diff --git a/src/eqtide.c b/src/eqtide.c
@@ -1490,7 +1490,8 @@ void VerifyCPL(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
       /* Tidal Q was also set */
       if (control->Io.iVerbose >= VERBINPUT) {
         fprintf(stderr,
-                "WARNING: Phase lag model selected, but both %s and %s set in file %s. "
+                "WARNING: Phase lag model selected, but both %s and %s set in "
+                "file %s. "
                 "Using %s = %lf and ignoring %s.\n",
                 options[OPT_TIDALTAU].cName, options[OPT_TIDALQ].cName,
                 options[OPT_TIDALTAU].cFile[iBody + 1],
@@ -1591,10 +1592,9 @@ void VerifyCPL(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
 void VerifyPerturbersEqtide(BODY *body, FILES *files, OPTIONS *options,
                             UPDATE *update, int iNumBodies, int iBody) {
   int iPert, iBodyPert, iVar, ok;
-  int *bFound = malloc(iNumBodies*sizeof(int));
+  int *bFound = malloc(iNumBodies * sizeof(int));
 
   for (iBody = 0; iBody < iNumBodies; iBody++) {
-    fprintf(stderr,"Body: %s\n",body[iBody].cName);
     if (body[iBody].bEqtide) {
       if (body[iBody].iTidePerts > 0) {
         for (iPert = 0; iPert < body[iBody].iTidePerts; iPert++) {
@@ -2665,8 +2665,8 @@ void WriteEqRotPer(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
                    UNITS *units, UPDATE *update, int iBody, double *dTmp,
                    char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
-  *dTmp = fdFreqToPer(fdEqRotRate(
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
+  *dTmp        = fdFreqToPer(fdEqRotRate(
         body, iBody, body[iOrbiter].dMeanMotion, body[iOrbiter].dEccSq,
         control->Evolve.iEqtideModel, control->Evolve.bDiscreteRot));
 
@@ -2683,7 +2683,7 @@ void WriteEqRotPerCont(BODY *body, CONTROL *control, OUTPUT *output,
                        SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                        double *dTmp, char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
 
   // To CPL, or to CTL? That is the question
   if (control->Evolve.iEqtideModel == CPL) {
@@ -2708,7 +2708,7 @@ void WriteEqRotPerDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
                            SYSTEM *system, UNITS *units, UPDATE *update,
                            int iBody, double *dTmp, char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdFreqToPer(fdCPLEqRotRateDiscrete(body[iOrbiter].dMeanMotion,
                                                body[iOrbiter].dEccSq));
@@ -2729,10 +2729,10 @@ void WriteEqRotRate(BODY *body, CONTROL *control, OUTPUT *output,
                     SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                     double *dTmp, char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
-  *dTmp = fdEqRotRate(body, iBody, body[iOrbiter].dMeanMotion,
-                      body[iOrbiter].dEccSq, control->Evolve.iEqtideModel,
-                      control->Evolve.bDiscreteRot);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
+  *dTmp        = fdEqRotRate(body, iBody, body[iOrbiter].dMeanMotion,
+                             body[iOrbiter].dEccSq, control->Evolve.iEqtideModel,
+                             control->Evolve.bDiscreteRot);
 
   if (output->bDoNeg[iBody]) {
     *dTmp *= output->dNeg;
@@ -2747,7 +2747,7 @@ void WriteEqRotRateCont(BODY *body, CONTROL *control, OUTPUT *output,
                         SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                         double *dTmp, char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
 
   // To CPL, or to CTL? That is the question
   if (control->Evolve.iEqtideModel == CPL) {
@@ -2770,7 +2770,7 @@ void WriteEqRotRateCont(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotRateDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
                             SYSTEM *system, UNITS *units, UPDATE *update,
                             int iBody, double *dTmp, char cUnit[]) {
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
 
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdCPLEqRotRateDiscrete(body[iOrbiter].dMeanMotion,
@@ -2792,7 +2792,7 @@ void WriteEqTidePower(BODY *body, CONTROL *control, OUTPUT *output,
                       SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                       double *dTmp, char cUnit[]) {
 
-  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  int iOrbiter = fiAssignTidalOrbiter(body, iBody);
   // Why is this Eq??????? XXX
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdCPLTidePowerEq(body[iBody].dTidalZ[iOrbiter], body[iBody].dEccSq,
@@ -2932,13 +2932,13 @@ void WritePowerEqtide(BODY *body, CONTROL *control, OUTPUT *output,
   //*dTmp = body[iBody].dTidalPowMan;
 
 
-  int jBody=0;
-  if (iBody==0)
+  int jBody = 0;
+  if (iBody == 0)
     jBody++;
-    /*
-  fprintf(stderr,"\niBody: %d\n",iBody);
-  fprintf(stderr,"\njBody: %d\n",jBody);
-  fprintf(stderr,"TidalZ: %lf\n",body[iBody].dTidalZ[jBody]);
+  /*
+fprintf(stderr,"\niBody: %d\n",iBody);
+fprintf(stderr,"\njBody: %d\n",jBody);
+fprintf(stderr,"TidalZ: %lf\n",body[iBody].dTidalZ[jBody]);
 */
 
   if (output->bDoNeg[iBody]) {
@@ -3433,7 +3433,7 @@ void LogBodyEqtide(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
   fprintf(fp, "----- EQTIDE PARAMETERS (%s)------\n", body[iBody].cName);
   for (iOut = iStart; iOut < OUTENDEQTIDE; iOut++) {
     if (output[iOut].iNum > 0) {
-      //fprintf(stderr,"iOut = %d.\n",iOut);
+      // fprintf(stderr,"iOut = %d.\n",iOut);
       WriteLogEntry(body, control, &output[iOut], system, update, fnWrite[iOut],
                     fp, iBody);
     }
@@ -3906,11 +3906,11 @@ void ForceBehaviorEqtide(BODY *body, MODULE *module, EVOLVE *evolve, IO *io,
 
 /**
  * Identify and return the index for the perturbing body.
- * 
- * 
-*/
+ *
+ *
+ */
 
-int fiAssignTidalPerturber(BODY *body,int iBody) {
+int fiAssignTidalPerturber(BODY *body, int iBody) {
   int iPerturber;
 
   if (!bPrimary(body, iBody)) {
@@ -3925,7 +3925,7 @@ int fiAssignTidalPerturber(BODY *body,int iBody) {
 
 /* Non-central body contains the orbital parameters */
 
-int fiAssignTidalOrbiter(BODY *body,int iBody) {
+int fiAssignTidalOrbiter(BODY *body, int iBody) {
   int iOrbiter;
 
   if (!bPrimary(body, iBody)) {
@@ -3963,8 +3963,8 @@ double fdCPLTidePower(BODY *body, int iBody) {
 
     // Does this work with DF's changes to da/dt with the synchronous case?
     // See Fleming et al., 2018
-    //fprintf(stderr,"\niBody: %d\n",iBody);
-    //fprintf(stderr,"TidalZ[%d]: %lf\n",iIndex,body[iBody].dTidalZ[iIndex]);
+    // fprintf(stderr,"\niBody: %d\n",iBody);
+    // fprintf(stderr,"TidalZ[%d]: %lf\n",iIndex,body[iBody].dTidalZ[iIndex]);
 
     dOrbPow += -body[iBody].dTidalZ[iIndex] / 8 *
                (4 * body[iBody].iTidalEpsilon[iIndex][0] +

diff --git a/src/evolve.c b/src/evolve.c
@@ -423,7 +423,8 @@ void RungeKutta4Step(BODY *body, CONTROL *control, SYSTEM *system,
       iNumEqns   = update[iBody].iNumEqns[iVar];
       for (iEqn = 0; iEqn < iNumEqns; iEqn++) {
         daDerivVar += iDir * evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
-        evolve->daDerivProc[0][iBody][iVar][iEqn] = evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
+        evolve->daDerivProc[0][iBody][iVar][iEqn] =
+              evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
       }
       evolve->daDeriv[0][iBody][iVar] = daDerivVar;
     }
@@ -465,7 +466,8 @@ void RungeKutta4Step(BODY *body, CONTROL *control, SYSTEM *system,
       iNumEqns   = update[iBody].iNumEqns[iVar];
       for (iEqn = 0; iEqn < iNumEqns; iEqn++) {
         daDerivVar += iDir * evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
-        evolve->daDerivProc[1][iBody][iVar][iEqn] = evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
+        evolve->daDerivProc[1][iBody][iVar][iEqn] =
+              evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
       }
       evolve->daDeriv[1][iBody][iVar] = daDerivVar;
     }
@@ -508,7 +510,8 @@ void RungeKutta4Step(BODY *body, CONTROL *control, SYSTEM *system,
       iNumEqns   = update[iBody].iNumEqns[iVar];
       for (iEqn = 0; iEqn < iNumEqns; iEqn++) {
         daDerivVar += iDir * evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
-        evolve->daDerivProc[2][iBody][iVar][iEqn] = evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
+        evolve->daDerivProc[2][iBody][iVar][iEqn] =
+              evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
       }
       evolve->daDeriv[2][iBody][iVar] = daDerivVar;
     }
@@ -557,7 +560,8 @@ void RungeKutta4Step(BODY *body, CONTROL *control, SYSTEM *system,
         iNumEqns                        = update[iBody].iNumEqns[iVar];
         for (iEqn = 0; iEqn < iNumEqns; iEqn++) {
           daDerivVar += iDir * evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
-          evolve->daDerivProc[3][iBody][iVar][iEqn] = evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
+          evolve->daDerivProc[3][iBody][iVar][iEqn] =
+                evolve->tmpUpdate[iBody].daDerivProc[iVar][iEqn];
         }
         evolve->daDeriv[3][iBody][iVar] = daDerivVar;
       }
@@ -574,11 +578,12 @@ void RungeKutta4Step(BODY *body, CONTROL *control, SYSTEM *system,
                                      2 * evolve->daDeriv[2][iBody][iVar] +
                                      evolve->daDeriv[3][iBody][iVar]);
       for (iEqn = 0; iEqn < update[iBody].iNumEqns[iVar]; iEqn++) {
-        update[iBody].daDerivProc[iVar][iEqn] = 1. / 6 *
-                                    (evolve->daDerivProc[0][iBody][iVar][iEqn] +
-                                     2 * evolve->daDerivProc[1][iBody][iVar][iEqn] +
-                                     2 * evolve->daDerivProc[2][iBody][iVar][iEqn] +
-                                     evolve->daDerivProc[3][iBody][iVar][iEqn]);
+        update[iBody].daDerivProc[iVar][iEqn] =
+              1. / 6 *
+              (evolve->daDerivProc[0][iBody][iVar][iEqn] +
+               2 * evolve->daDerivProc[1][iBody][iVar][iEqn] +
+               2 * evolve->daDerivProc[2][iBody][iVar][iEqn] +
+               evolve->daDerivProc[3][iBody][iVar][iEqn]);
       }
 
       if (update[iBody].iaType[iVar][0] == 0 ||
@@ -603,11 +608,9 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
             fnUpdateVariable ***fnUpdate, fnWriteOutput *fnWrite,
             fnIntegrate fnOneStep) {
   /* Master evolution routine that controls the simulation integration. */
-  int iDir, iBody, iModule, nSteps; // Dummy counting variables
-  double dDt, dFoo;                 // Next timestep, dummy variable
-  double dEqSpinRate;               // Store the equilibrium spin rate
-
-  nSteps = 0;
+  int iDir, iBody, iModule; // Dummy counting variables
+  double dDt, dFoo;         // Next timestep, dummy variable
+  double dEqSpinRate;       // Store the equilibrium spin rate
 
   if (control->Evolve.bDoForward) {
     iDir = 1;
@@ -645,6 +648,8 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
    *
    */
 
+  control->Evolve.iStepsSinceLastOutput = 0;
+  control->Evolve.iTotalSteps           = 0;
   while (control->Evolve.dTime < control->Evolve.dStopTime) {
     /* Take one step */
     fnOneStep(body, control, system, update, fnUpdate, &dDt, iDir);
@@ -678,17 +683,19 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module,
     }
 
     control->Evolve.dTime += dDt;
-    nSteps++;
+    control->Evolve.iStepsSinceLastOutput++;
 
     /* Time for Output? */
     if (control->Evolve.dTime >= control->Io.dNextOutput) {
-      control->Evolve.nSteps += nSteps;
+      control->Evolve.iTotalSteps += control->Evolve.iStepsSinceLastOutput;
       WriteOutput(body, control, files, output, system, update, fnWrite);
       // Timesteps are synchronized with the output time, so this statement is
       // sufficient
       control->Io.dNextOutput += control->Io.dOutputTime;
-      //printf("%d\n",nSteps);
-      nSteps = 0;
+      if (control->Evolve.dTime < control->Evolve.dStopTime) {
+        // Reset counter if we are not at the end of the simulation
+        control->Evolve.iStepsSinceLastOutput = 0;
+      }
     }
 
     /* Get auxiliary properties for next step -- first call

diff --git a/src/magmoc.c b/src/magmoc.c
@@ -1432,7 +1432,6 @@ void PropsAuxMagmOc(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update,
                     int iBody) {
   double dCurrentTime     = evolve->dTime;
   double dCurrentTimeStep = evolve->dTimeStep;
-  double dCurrentStepNum  = evolve->nSteps;
   double dAveMolarMassAtm;
 
   // body[iBody].dSurfTemp = body[iBody].dPotTemp;

diff --git a/src/output.c b/src/output.c
@@ -92,7 +92,7 @@ void WriteDeltaTime(BODY *body, CONTROL *control, OUTPUT *output,
 
   if (control->Evolve.bVarDt) {
     if (control->Evolve.dTime > 0) {
-      *dTmp = control->Io.dOutputTime / control->Evolve.nSteps;
+      *dTmp = control->Io.dOutputTime / control->Evolve.iStepsSinceLastOutput;
     } else {
       if (control->Io.iVerbose >= VERBINPUT && !control->Io.bDeltaTimeMessage) {
         fprintf(stderr, "INFO: DeltaTime output for first step is defined to "
@@ -2199,16 +2199,6 @@ void WriteLog(BODY *body, CONTROL *control, FILES *files, MODULE *module,
 
   /* Bodies' Properties */
   LogBody(body, control, files, module, output, system, fnWrite, fp, update);
-
-  /* Deprecated
-  if (iEnd) {
-    dTotTime = difftime(time(NULL),dStartTime);
-    fprintf(fp,"\nRuntime = %d s\n", (int)dTotTime);
-    fprintf(fp,"Total Number of Steps = %d\n",control->Evolve.nSteps);
-    if (control->Io.iVerbose >= VERBPROG)
-      printf("Runtime = %d s\n", (int)dTotTime);
-  }
-  */
   fclose(fp);
 }
 

diff --git a/src/verify.c b/src/verify.c
@@ -1109,7 +1109,6 @@ void VerifyOptions(BODY *body, CONTROL *control, FILES *files, MODULE *module,
   int iBody, iModule;
 
   control->Evolve.dTime  = 0;
-  control->Evolve.nSteps = 0;
 
   VerifyAge(body, control, options);
   VerifyNames(body, control, options);

diff --git a/src/vplanet.h b/src/vplanet.h
@@ -1734,7 +1734,8 @@ struct EVOLVE {
   double dStopTime;  /**< Integration Stop Time */
   double dTimeStep;  /**< Integration Time step */
   int bVarDt;        /**< Use Variable Timestep? */
-  int nSteps;        /**< Total Number of Steps */
+  int iTotalSteps;   /**< Total Number of Steps */
+  int iStepsSinceLastOutput;
   double dMinValue;  /**< Minimum Value for Eccentricity and Obliquity to be
                         Integrated */
   int bFirstStep;    /**< Has the First Dtep Been Taken? */