Skip to content

Commit

Permalink
Added new test for hydrogen loss.
Browse files Browse the repository at this point in the history 
Fixed several floating point errors.
  • Loading branch information
Rory Barnes committed Apr 16, 2024
1 parent b6ef1d0 commit a5a17ff
Show file tree
Hide file tree
Showing 9 changed files with 216 additions and 47 deletions.
51 changes: 30 additions & 21 deletions src/atmesc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3356,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 @@ -4303,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
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
24 changes: 24 additions & 0 deletions tests/Atmesc/HydELimNoXUVLopez12/planet.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
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
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
119 changes: 119 additions & 0 deletions tests/Atmesc/HydELimNoXUVLopez12/test_HydELimNoXUVLopez12.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
from benchmark import Benchmark, benchmark
import astropy.units as u

@benchmark(
{
"log.initial.system.Age": {"value": 3.155760e+13, "unit": u.sec},
"log.initial.system.Time": {"value": 0.000000, "unit": u.sec},
"log.initial.system.TotAngMom": {"value": 6.108249e+36, "unit": (u.kg * u.m ** 2) / u.sec},
"log.initial.system.TotEnergy": {"value": 1.948502e+32, "unit": u.Joule},
"log.initial.system.PotEnergy": {"value": -2.725202e+31, "unit": u.Joule},
"log.initial.system.KinEnergy": {"value": 2.221022e+32, "unit": u.Joule},
"log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec},
"log.initial.planet.Mass": {"value": 2.000000, "unit": u.Mearth},
"log.initial.planet.Radius": {"value": 2.096446e+08, "unit": u.m},
"log.initial.planet.RadGyra": {"value": 0.400000},
"log.initial.planet.BodyType": {"value": 0.000000},
"log.initial.planet.Density": {"value": 0.309474, "unit": u.kg / u.m ** 3},
"log.initial.planet.HZLimitDryRunaway": {"value": -1.000000, "unit": u.m},
"log.initial.planet.HZLimRecVenus": {"value": -1.000000},
"log.initial.planet.HZLimRunaway": {"value": -1.000000},
"log.initial.planet.HZLimMoistGreenhouse": {"value": -1.000000},
"log.initial.planet.HZLimMaxGreenhouse": {"value": -1.000000},
"log.initial.planet.HZLimEarlyMars": {"value": -1.000000},
"log.initial.planet.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3},
"log.initial.planet.MeanMotion": {"value": -1.000000, "unit": 1 / u.sec},
"log.initial.planet.OrbPeriod": {"value": -1.000000, "unit": u.sec},
"log.initial.planet.SemiMajorAxis": {"value": -1.000000, "unit": u.m},
"log.initial.planet.LXUVTot": {"value": -1.000000, "unit": u.kg / u.sec ** 3},
"log.initial.planet.SurfWaterMass": {"value": 0.000000, "unit": u.kg},
"log.initial.planet.EnvelopeMass": {"value": 1.000000, "unit": u.Mearth},
"log.initial.planet.OxygenMass": {"value": 0.000000, "unit": u.kg},
"log.initial.planet.RGLimit": {"value": 0.000000, "unit": u.m},
"log.initial.planet.XO": {"value": 0.000000},
"log.initial.planet.EtaO": {"value": 0.000000},
"log.initial.planet.PlanetRadius": {"value": 32.869442, "unit": u.Rearth},
"log.initial.planet.OxygenMantleMass": {"value": 0.000000, "unit": u.kg},
"log.initial.planet.RadXUV": {"value": -1.000000, "unit": u.m},
"log.initial.planet.RadSolid": {"value": -1.000000, "unit": u.m},
"log.initial.planet.PresXUV": {"value": 5.000000},
"log.initial.planet.ScaleHeight": {"value": -1.000000, "unit": u.m},
"log.initial.planet.ThermTemp": {"value": 400.000000, "unit": u.K},
"log.initial.planet.AtmGasConst": {"value": 4124.000000},
"log.initial.planet.PresSurf": {"value": -1.000000, "unit": u.Pa},
"log.initial.planet.DEnvMassDt": {"value": 0.000000, "unit": u.kg / u.sec},
"log.initial.planet.FXUV": {"value": 0.000000, "unit": u.W / u.m ** 2},
"log.initial.planet.AtmXAbsEffH2O": {"value": 0.300000},
"log.initial.planet.RocheRadius": {"value": 1626.273816, "unit": u.Rearth},
"log.initial.planet.BondiRadius": {"value": 21.602817, "unit": u.Rearth},
"log.initial.planet.HEscapeRegime": {"value": 3.000000},
"log.initial.planet.RRCriticalFlux": {"value": 0.012799, "unit": u.W / u.m ** 2},
"log.initial.planet.CrossoverMass": {"value": 0.000000, "unit": u.kg},
"log.initial.planet.WaterEscapeRegime": {"value": 8.000000},
"log.initial.planet.FXUVCRITDRAG": {"value": 3.000221e-05, "unit": u.W / u.m ** 2},
"log.initial.planet.HREFFLUX": {"value": 0.000000, "unit": 1 / u.m ** 2 / u.sec},
"log.initial.planet.XO2": {"value": 0.000000},
"log.initial.planet.XH2O": {"value": 0.000000},
"log.initial.planet.HDiffFlux": {"value": 3.512237e+14, "unit": 1 / u.m ** 2 / u.sec},
"log.initial.planet.HRefODragMod": {"value": 1.000000},
"log.initial.planet.KTide": {"value": 0.969687},
"log.initial.planet.RGDuration": {"value": 1.00000e+06, "unit": u.yr},
"log.final.system.Age": {"value": 6.311520e+13, "unit": u.sec},
"log.final.system.Time": {"value": 3.155760e+13, "unit": u.sec},
"log.final.system.TotAngMom": {"value": 6.108249e+36, "unit": (u.kg * u.m ** 2) / u.sec},
"log.final.system.TotEnergy": {"value": 1.948502e+32, "unit": u.Joule},
"log.final.system.PotEnergy": {"value": -2.725202e+31, "unit": u.Joule},
"log.final.system.KinEnergy": {"value": 2.221022e+32, "unit": u.Joule},
"log.final.system.DeltaTime": {"value": 3.155760e+13, "unit": u.sec},
"log.final.planet.Mass": {"value": 2.000000, "unit": u.Mearth},
"log.final.planet.Radius": {"value": 2.096446e+08, "unit": u.m},
"log.final.planet.RadGyra": {"value": 0.400000},
"log.final.planet.BodyType": {"value": 0.000000},
"log.final.planet.Density": {"value": 0.309474, "unit": u.kg / u.m ** 3},
"log.final.planet.HZLimitDryRunaway": {"value": -1.000000, "unit": u.m},
"log.final.planet.HZLimRecVenus": {"value": -1.000000},
"log.final.planet.HZLimRunaway": {"value": -1.000000},
"log.final.planet.HZLimMoistGreenhouse": {"value": -1.000000},
"log.final.planet.HZLimMaxGreenhouse": {"value": -1.000000},
"log.final.planet.HZLimEarlyMars": {"value": -1.000000},
"log.final.planet.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3},
"log.final.planet.MeanMotion": {"value": -1.000000, "unit": 1 / u.sec},
"log.final.planet.OrbPeriod": {"value": -1.000000, "unit": u.sec},
"log.final.planet.SemiMajorAxis": {"value": -1.000000, "unit": u.m},
"log.final.planet.LXUVTot": {"value": -1.000000, "unit": u.kg / u.sec ** 3},
"log.final.planet.SurfWaterMass": {"value": 0.000000, "unit": u.kg},
"log.final.planet.EnvelopeMass": {"value": 1.000000, "unit": u.Mearth},
"log.final.planet.OxygenMass": {"value": 0.000000, "unit": u.kg},
"log.final.planet.RGLimit": {"value": 0.000000, "unit": u.m},
"log.final.planet.XO": {"value": 0.000000},
"log.final.planet.EtaO": {"value": 0.000000},
"log.final.planet.PlanetRadius": {"value": 32.869442, "unit": u.Rearth},
"log.final.planet.OxygenMantleMass": {"value": 0.000000, "unit": u.kg},
"log.final.planet.RadXUV": {"value": -1.000000, "unit": u.m},
"log.final.planet.RadSolid": {"value": -1.000000, "unit": u.m},
"log.final.planet.PresXUV": {"value": 5.000000},
"log.final.planet.ScaleHeight": {"value": -1.000000, "unit": u.m},
"log.final.planet.ThermTemp": {"value": 400.000000, "unit": u.K},
"log.final.planet.AtmGasConst": {"value": 4124.000000},
"log.final.planet.PresSurf": {"value": -1.000000, "unit": u.Pa},
"log.final.planet.DEnvMassDt": {"value": 0.000000, "unit": u.kg / u.sec},
"log.final.planet.FXUV": {"value": 0.000000, "unit": u.W / u.m ** 2},
"log.final.planet.AtmXAbsEffH2O": {"value": 0.300000},
"log.final.planet.RocheRadius": {"value": 1626.273816, "unit": u.Rearth},
"log.final.planet.BondiRadius": {"value": 21.602817, "unit": u.Rearth},
"log.final.planet.HEscapeRegime": {"value": 3.000000},
"log.final.planet.RRCriticalFlux": {"value": 0.012799, "unit": u.W / u.m ** 2},
"log.final.planet.CrossoverMass": {"value": 0.000000, "unit": u.kg},
"log.final.planet.WaterEscapeRegime": {"value": 8.000000},
"log.final.planet.FXUVCRITDRAG": {"value": 3.000221e-05, "unit": u.W / u.m ** 2},
"log.final.planet.HREFFLUX": {"value": 0.000000, "unit": 1 / u.m ** 2 / u.sec},
"log.final.planet.XO2": {"value": 0.000000},
"log.final.planet.XH2O": {"value": 0.000000},
"log.final.planet.HDiffFlux": {"value": 3.512237e+14, "unit": 1 / u.m ** 2 / u.sec},
"log.final.planet.HRefODragMod": {"value": 1.000000},
"log.final.planet.KTide": {"value": 0.969687},
"log.final.planet.RGDuration": {"value": 1.00000e+06, "unit": u.yr},
}
)
class Test_HydELimConstXUVLopez12(Benchmark):
pass
16 changes: 16 additions & 0 deletions tests/Atmesc/HydELimNoXUVLopez12/vpl.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
sSystemName system
iVerbose 5
bOverwrite 1
saBodyFiles planet.in
sUnitMass solar
sUnitLength AU
sUnitTime YEARS
sUnitAngle d
bDoLog 1
iDigits 6
dMinValue 1e-10
bDoForward 1
bVarDt 1
dEta 0.01
dStopTime 1e6
dOutputTime 1e6
2 changes: 1 addition & 1 deletion tests/AtmescSpinbody/NBodyAtmEsc/vpl.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,6 @@ iVerbose 5 # Verbosity level max=5
bDoForward 1 # Perform a forward evolution?
bVarDt 0 # Use variable timestepping?
#dEta .1
dTimeStep 0.0005 # Timestep in years
dTimeStep 0.0001 # Timestep in years
dStopTime 100 # Stop time for evolution
dOutputTime 100 # Output timesteps (assuming in body files)
5 changes: 0 additions & 5 deletions tests/floatingpoint.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,6 @@ def Main():
if last_line != "Simulation completed.\n":
tot_fail += 1
print("Fail", flush=True)
print("\n------------\n")
with open(outfile, "r") as f:
print(f.read())
print("\n------------\n")

else:
print("Pass", flush=True)
os.chdir("../../")
Expand Down

0 comments on commit a5a17ff

Please sign in to comment.