Skip to content

Commit

Permalink
Merge pull request #289 from VirtualPlanetaryLaboratory/ConstXUV
Browse files Browse the repository at this point in the history 
Const xuv
  • Loading branch information
@RoryBarnes
RoryBarnes authored Apr 16, 2024
2 parents 95b0b1b + 192afbd commit 2fadd16
Show file tree
Hide file tree
Showing 62 changed files with 24,082 additions and 12,699 deletions.
140 changes: 75 additions & 65 deletions src/atmesc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2108,43 +2108,37 @@ void VerifyAtmEsc(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
}

if (body[iBody].iPlanetRadiusModel == ATMESC_LEHMER17) {
if (body[iBody].dEnvelopeMass >= 0.5 * body[iBody].dMass) {
fprintf(stderr,
"ERROR: %s's Envelope mass is greater than 50%% of its total "
"mass, which ",
if (body[0].bStellar) {
if (body[iBody].dEnvelopeMass > 0.5 * body[iBody].dMass) {
fprintf(stderr,
"ERROR: %s's Envelope mass is greater than 50%% of its total "
"mass, which ",
body[iBody].cName);
fprintf(
stderr,
"is not allowed for the Lehmer-Catling (2017) envelope model.\n");
DoubleLineExit(files->Infile[iBody + 1].cIn, files->Infile[iBody + 1].cIn,
options[OPT_ENVELOPEMASS].iLine[iBody + 1],
options[OPT_MASS].iLine[iBody + 1]);
}
if (body[iBody].dEnvelopeMass >= 0.1 * body[iBody].dMass) {
fprintf(
stderr,
"WARNING: Envelope masses more than 10%% of the total mass are not "
"recommended for the Lehmer-Catling (2017) envelope model. %s's "
"envelope ",
body[iBody].cName);
fprintf(
stderr,
"is not allowed for the Lehmer-Catling (2017) envelope model.\n");
DoubleLineExit(files->Infile[iBody + 1].cIn, files->Infile[iBody + 1].cIn,
options[OPT_ENVELOPEMASS].iLine[iBody + 1],
options[OPT_ENVELOPEMASS].iLine[iBody + 1]);
}
if (body[iBody].dEnvelopeMass >= 0.1 * body[iBody].dMass) {
fprintf(
stderr,
"WARNING: Envelope masses more than 10%% of the total mass are not "
"recommended for the Lehmer-Catling (2017) envelope model. %s's "
"envelope ",
body[iBody].cName);
fprintf(stderr, "mass exceeds this threshold.\n");
}
fprintf(stderr, "mass exceeds this threshold.\n");
}

// Get thermal temperature
if (body[iBody].bAutoThermTemp) {
body[iBody].dThermTemp = fdThermalTemp(body, iBody);
// Calculate auxiliary properties
body[iBody].dRadSolid = fdMassToRad_LehmerCatling17(body[iBody].dMass - body[iBody].dEnvelopeMass);
AuxPropsLehmer17(body,iBody);
} else {
fprintf(stderr,
"ERROR: The Lehmer & Catling (2017) model requires a star.\n");
exit(EXIT_INPUT);
}
body[iBody].dRadSolid = fdMassToRad_LehmerCatling17(
body[iBody].dMass - body[iBody].dEnvelopeMass);
body[iBody].dGravAccel = BIGG *
(body[iBody].dMass - body[iBody].dEnvelopeMass) /
(body[iBody].dRadSolid * body[iBody].dRadSolid);
body[iBody].dScaleHeight = body[iBody].dAtmGasConst *
body[iBody].dThermTemp / body[iBody].dGravAccel;
body[iBody].dPresSurf =
fdLehmerPres(body[iBody].dEnvelopeMass, body[iBody].dGravAccel,
body[iBody].dRadSolid);
body[iBody].dRadXUV = fdLehmerRadius(body, iBody);
} else {
int iCol, bError = 0;
for (iCol = 0; iCol < files->Outfile[iBody].iNumCols; iCol++) {
Expand Down Expand Up @@ -2830,10 +2824,14 @@ void WriteRGLimit(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
double flux = fdHZRG14(body, iBody);

// Convert to semi-major axis *at current eccentricity!*
*dTmp = pow(4 * PI * flux /
(body[0].dLuminosity *
pow((1 - body[iBody].dEcc * body[iBody].dEcc), 0.5)),
-0.5);
if (body[0].dLuminosity == 0) {
*dTmp = -1;
} else {
*dTmp = pow(4 * PI * flux /
(body[0].dLuminosity *
pow((1 - body[iBody].dEcc * body[iBody].dEcc), 0.5)),
-0.5);
}

if (output->bDoNeg[iBody]) {
*dTmp *= output->dNeg;
Expand Down Expand Up @@ -3358,10 +3356,14 @@ void WriteFXUVCRITDRAG(BODY *body, CONTROL *control, OUTPUT *output,
double GPotential =
(BIGG * body[iBody].dMass * PROTONMASS) / (body[iBody].dRadius);

*dTmp = ((4 * BDIFF * pow(GPotential, 2)) /
(body[iBody].dAtmXAbsEffH2O * KBOLTZ * body[iBody].dFlowTemp *
body[iBody].dRadius)) *
(QOH - 1) * (1 - XO);
if (body[iBody].dAtmXAbsEffH2O > 0 && body[iBody].dFlowTemp > 0 && body[iBody].dRadius > 0) {
*dTmp = ((4 * BDIFF * pow(GPotential, 2)) /
(body[iBody].dAtmXAbsEffH2O * KBOLTZ * body[iBody].dFlowTemp *
body[iBody].dRadius)) *
(QOH - 1) * (1 - XO);
} else {
*dTmp = -1;
}
if (output->bDoNeg[iBody]) {
*dTmp *= output->dNeg;
strcpy(cUnit, output->cNeg);
Expand Down Expand Up @@ -4305,24 +4307,29 @@ double fdXUVEfficiencyBolmont2016(double dFXUV) {
double c3 = -0.89880083;

// Convert to erg/cm^2/s and take the log
double x = log10(dFXUV * 1.e3);
double y;

// Piecewise polynomial fit and handle edge cases
if ((x >= -2) && (x < -1)) {
y = pow(10, a0 * x * x + a1 * x + a2);
} else if ((x >= -1) && (x < 0)) {
y = pow(10, b0 * x * x * x + b1 * x * x + b2 * x + b3);
} else if ((x >= 0) && (x <= 5)) {
y = pow(10, c0 * x * x * x + c1 * x * x + c2 * x + c3);
} else if (x < -2) { // Lower flux bound
y = 0.001;
} else if (x > 5) { // Upper flux bound
y = 0.01;
} else { // Base case that never happens but DPF likes code symmetry
y = 0.1;
double dWaterEscapeEfficiency;

if (dFXUV > 0) {
double x = log10(dFXUV * 1.e3);

// Piecewise polynomial fit and handle edge cases
if ((x >= -2) && (x < -1)) {
dWaterEscapeEfficiency = pow(10, a0 * x * x + a1 * x + a2);
} else if ((x >= -1) && (x < 0)) {
dWaterEscapeEfficiency = pow(10, b0 * x * x * x + b1 * x * x + b2 * x + b3);
} else if ((x >= 0) && (x <= 5)) {
dWaterEscapeEfficiency = pow(10, c0 * x * x * x + c1 * x * x + c2 * x + c3);
} else if (x < -2) { // Lower flux bound
dWaterEscapeEfficiency = 0.001;
} else if (x > 5) { // Upper flux bound
dWaterEscapeEfficiency = 0.01;
} else { // Base case that never happens but DPF likes code symmetry
dWaterEscapeEfficiency = 0.1;
}
} else {
dWaterEscapeEfficiency = 0; // No escape if F_XUV = 0
}
return y;
return dWaterEscapeEfficiency;
}

/**
Expand Down Expand Up @@ -4400,14 +4407,17 @@ double fdRocheRadius(BODY *body, int iBody) {
@return Body's Bondi radius
*/
double fdBondiRadius(BODY *body, int iBody) {

double dBondiRadius;
// Compute sound speed in planet's atmosphere assuming a diatomic H atmosphere
// assuming body 0 is the star as it should be when using atmesc
double cs = fdEqH2AtmosphereSoundSpeed(body[0].dTemperature, body[0].dRadius,
// assuming body 0 is the star
if (body[0].bStellar) {
double dSoundSpeed = fdEqH2AtmosphereSoundSpeed(body[0].dTemperature, body[0].dRadius,
body[iBody].dSemi);
double rb = BIGG * body[iBody].dMass / (2.0 * cs * cs);

return rb;
dBondiRadius = BIGG * body[iBody].dMass / (2.0 * dSoundSpeed * dSoundSpeed);
} else {
dBondiRadius = -1;
}
return dBondiRadius;
}

/**
Expand Down
19 changes: 12 additions & 7 deletions src/body.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,10 +61,10 @@ double fdDPerDt(double dRotRate, double dDrotrateDt) {
@param dDeriv The value of the variable's time derivative
@return The timescale of the variable's change: |x/(dx/dt)|. If the
derivative is 0, return 0.
derivative is <1e-100, return 0.
*/
double fdTimescale(double dVar, double dDeriv) {
if (dDeriv != 0) {
if (dDeriv > 1e-100) {
return fabs(dVar / dDeriv);
} else {
return 0;
Expand All @@ -84,16 +84,16 @@ controlled by multiple processes.
@return The timescale of the variable's change: |x/Sum(dx/dt)|
*/
double fdTimescaleMulti(double dVar, double *dDeriv, int iNum) {
double dTime;
double dTime,dTotalDerivative;
int iPert;

dTime = 0;
dTotalDerivative = 0;
for (iPert = 0; iPert < iNum; iPert++) {
if (dDeriv[iPert] != 0) {
dTime += dDeriv[iPert]; // Note that here dTime is actullay the rate
dTotalDerivative += dDeriv[iPert];
}
dTime = fabs(dVar / dTime);
}
dTime = fabs(dVar / dTotalDerivative);
return dTime;
}

Expand Down Expand Up @@ -580,10 +580,15 @@ double CalcDynEllipEq(BODY *body, int iBody) {
double fdLehmerRadius(BODY *body, int iBody) {
double dRadXUV, dRoche;

dRadXUV = body[iBody].dRadSolid * body[iBody].dRadSolid /
// Set floor for surface pressure to prevent overflow error
if (body[iBody].dPresSurf > 1e-100) {
dRadXUV = body[iBody].dRadSolid * body[iBody].dRadSolid /
(body[iBody].dScaleHeight *
log(body[iBody].dPresXUV / body[iBody].dPresSurf) +
body[iBody].dRadSolid);
} else {
dRadXUV = body[iBody].dRadSolid;
}
dRoche = fdRocheRadius(body, iBody);
// printf("%lf %lf %lf %lf
// %lf\n",body[iBody].dPresXUV,body[iBody].dPresSurf,body[iBody].dGravAccel,body[iBody].dEnvelopeMass,dRadXUV);
Expand Down
3 changes: 2 additions & 1 deletion src/stellar.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -463,7 +463,7 @@ void InitializeOptionsStellar(OPTIONS *options, fnReadOption fnRead[]) {
"gigayears.");

sprintf(options[OPT_STELLARMODEL].cName, "sStellarModel");
sprintf(options[OPT_STELLARMODEL].cDescr, "Luminosity evolution model");
sprintf(options[OPT_STELLARMODEL].cDescr, "Stellar evolution model");
sprintf(options[OPT_STELLARMODEL].cDefault, "BARAFFE");
sprintf(options[OPT_STELLARMODEL].cValues, "BARAFFE PROXIMA SINEWAVE NONE");
options[OPT_STELLARMODEL].iType = 3;
Expand Down Expand Up @@ -1524,6 +1524,7 @@ double fdTemperature(BODY *body, SYSTEM *system, int *iaBody) {
}
}
if (body[iaBody[0]].iStellarModel == STELLAR_MODEL_SINEWAVE) {
printf("%lf\n",body[iaBody[0]].dLuminosity);
foo = fdEffectiveTemperature(body,iaBody[0]);
return foo;
}
Expand Down
24 changes: 12 additions & 12 deletions src/vplanet.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,18 +28,18 @@ We need this wrapper so we can call `main_impl` from Python.
*/
int main_impl(int argc, char *argv[]) {
#ifdef DEBUG
#ifdef __x86_64__
// feenableexcept(FE_INVALID | FE_OVERFLOW);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO);
//_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW);
fprintf(stderr, "INFO: Floating point trapping enabled.\n");
#else
fprintf(stderr,
"WARNING: Floating point trapping only enabled for x86 "
"architectures.\n");
#endif
#ifdef __x86_64__
// feenableexcept(FE_INVALID | FE_OVERFLOW);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO);
//_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW);
fprintf(stderr, "INFO: Floating point trapping enabled.\n");
#else
fprintf(stderr,
"WARNING: Floating point trapping only enabled for x86 "
"architectures.\n");
#endif
#endif

// struct timeval start, end;
Expand Down
25 changes: 25 additions & 0 deletions tests/Atmesc/HydELimConstXUVLopez12/planet.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
sName planet # Body's name
saModules atmesc # Active modules

# Physical Parameters
dMass -2.0 # Mass, negative -> Earth masses
sPlanetRadiusModel lopez12 # Mass-radius model
dRadGyra 0.4 # Radius of gyration; ang. mom. coeff.
dAge 1.0e6 # Age [yr]

# ATMESC Parameters
dFXUV -100 # Incident XUV flux (constant)
dXFrac 1.0 # X-Ray/XUV absorption radius in planet radii
dAtmXAbsEffH 0.1 # H X-ray/XUV absorption efficiency (epsilon)
dSurfWaterMass 0.0 # Initial water mass, negative -> Earth oceans
dEnvelopeMass -1.0 # Initial H envelope mass, negative -> Earth Mass
bInstantO2Sink 0 # Is Oxygen instantly absorbed by the surface?
bHaltSurfaceDesiccated 0 # Halt when dry?
bHaltEnvelopeGone 0 # Halt when H enevlope evaporated?
dJeansTime -12.0 # Time when flow transitions to ballistic escape (Gyr)
bUseEnergyLimited 1 # Is the flow energy-limited?
bUseRRLimited 0 # Is the flow radiation/recombination-limited?
bUseBondiLimited 0 # Is the flow Bondi-limited?
bAtmEscAuto 0 # Should atmesc decide the escape regime?

saOutputOrder Time -Mass -EnvelopeMass -PlanetRadius -BondiRadius -RocheRadius DEnvMassDt
Loading

0 comments on commit 2fadd16

Please sign in to comment.