Skip to content

Commit

Permalink
Merge pull request #263 from VirtualPlanetaryLaboratory/Steps
Browse files Browse the repository at this point in the history 
Steps
  • Loading branch information
@RoryBarnes
RoryBarnes authored Feb 5, 2024
2 parents a0b4594 + d15e244 commit aefaed8
Show file tree
Hide file tree
Showing 56 changed files with 3,924 additions and 1,974 deletions.
324 changes: 185 additions & 139 deletions src/atmesc.c
View file

Large diffs are not rendered by default.

56 changes: 28 additions & 28 deletions src/eqtide.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand Down Expand Up @@ -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++) {
Expand Down Expand Up @@ -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));

Expand All @@ -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) {
Expand All @@ -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));
Expand All @@ -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;
Expand All @@ -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) {
Expand All @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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]) {
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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)) {
Expand All @@ -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)) {
Expand Down Expand Up @@ -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] +
Expand Down
43 changes: 25 additions & 18 deletions src/evolve.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand All @@ -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 ||
Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion src/magmoc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
12 changes: 1 addition & 11 deletions src/output.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 "
Expand Down Expand Up @@ -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);
}

Expand Down
1 change: 0 additions & 1 deletion src/verify.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
3 changes: 2 additions & 1 deletion src/vplanet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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? */
Expand Down
Loading

0 comments on commit aefaed8

Please sign in to comment.