Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Steps #263

Merged
merged 2 commits into from
Feb 5, 2024
Merged

Steps #263

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ae3a1e1
Changed output DeltaTime to be average timestep over previous output.
Feb 5, 2024
d15e244
Updated testing for better performance when recreating all tests. Als…
RoryBarnes Feb 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
324 changes: 185 additions & 139 deletions src/atmesc.c
View file
Open in desktop

Large diffs are not rendered by default.

56 changes: 28 additions & 28 deletions src/eqtide.c
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading