diff --git a/Makefile b/Makefile index 8e36d81f8..a0fea5c7c 100644 --- a/Makefile +++ b/Makefile @@ -30,7 +30,7 @@ legacy: @echo "==========================================================================================================" debug: - -gcc -g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -gcc -g -D DEBUG -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" debug_no_AE: -gcc -g -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" @@ -52,7 +52,7 @@ cpp: g++ -o bin/vplanet src/*.c -lm -O3 -fopenmp -fpermissive -w -DGITVERSION=\"$(GITVERSION)\" warnings: - clang -Weverything src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" + -gcc -g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" parallel: gcc -o bin/vplanet src/*.c -lm -O3 -fopenmp -DGITVERSION=\"$(GITVERSION)\" diff --git a/examples/CosmicShoreline/CosmicShoreline.png b/examples/CosmicShoreline/CosmicShoreline.png new file mode 100644 index 000000000..92370e185 Binary files /dev/null and b/examples/CosmicShoreline/CosmicShoreline.png differ diff --git a/examples/CosmicShoreline/README.rst b/examples/CosmicShoreline/README.rst new file mode 100644 index 000000000..ff41549f3 --- /dev/null +++ b/examples/CosmicShoreline/README.rst @@ -0,0 +1,37 @@ +The "Cosmic Shoreline" +========== + +Overview +-------- + +A reproduction of Zahnle & Catling (2018)'s cosmic shoreline, +a proposed boundary between planet with and without atmosphere. + +=================== ============ +**Date** 09/12/2024 +**Author** Rory Barnes +**Modules** Atmesc, STELLAR +**Approx. runtime** 10 seonds +=================== ============ + + + + +To run this example +------------------- + +.. code-block:: bash + + python makeplot.py + + +Expected output +--------------- + +.. figure:: CosmicShoreline.png + :width: 100% + :align: center + +The blue line is the cosmic shoreline and the blue points +represent the Solar System planets. The Sun's XUV luminosity +is assumed to follow the Ribas et al. (2005) model. diff --git a/examples/CosmicShoreline/earth.in b/examples/CosmicShoreline/earth.in new file mode 100644 index 000000000..d0374ceb5 --- /dev/null +++ b/examples/CosmicShoreline/earth.in @@ -0,0 +1,6 @@ +sName Earth +saModules atmesc +dMass 3.0018090452e-06 +dRadius -1. +dSemi 1.00000321 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/george.in b/examples/CosmicShoreline/george.in new file mode 100644 index 000000000..618508379 --- /dev/null +++ b/examples/CosmicShoreline/george.in @@ -0,0 +1,6 @@ +sName George +saModules atmesc +dMass -14.54 +dRadius -4.007 +dSemi 19.19203990 +saOutputOrder Time -fXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/jupiter.in b/examples/CosmicShoreline/jupiter.in new file mode 100644 index 000000000..d3c34bc2c --- /dev/null +++ b/examples/CosmicShoreline/jupiter.in @@ -0,0 +1,6 @@ +sName Jupiter +saModules atmesc +dMass -317.83 +dRadius -11.2 +dSemi 5.20880408 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/makeplot.py b/examples/CosmicShoreline/makeplot.py new file mode 100644 index 000000000..e74a34504 --- /dev/null +++ b/examples/CosmicShoreline/makeplot.py @@ -0,0 +1,73 @@ +import os +import pathlib +import subprocess +import sys + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import vplot as vpl +from matplotlib.ticker import MaxNLocator +from mpl_toolkits.axes_grid1 import make_axes_locatable + +import vplanet + +path = pathlib.Path(__file__).parents[0].absolute() +sys.path.insert(1, str(path.parents[0])) + +output = vplanet.run(units = False) + +# Plot! +fig = plt.figure(figsize=(8.5, 6)) + +fxuv_earth = output.log.final.Earth.CumulativeFXUV + +fxuv = [] +fxuv.append(output.log.final.Mercury.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Venus.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Earth.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Mars.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Jupiter.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Saturn.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.George.CumulativeFXUV/fxuv_earth) +fxuv.append(output.log.final.Neptune.CumulativeFXUV/fxuv_earth) + +escvel = [] +escvel.append(output.log.final.Mercury.EscapeVelocity/1e3) +escvel.append(output.log.final.Venus.EscapeVelocity/1e3) +escvel.append(output.log.final.Earth.EscapeVelocity/1e3) +escvel.append(output.log.final.Mars.EscapeVelocity/1e3) +escvel.append(output.log.final.Jupiter.EscapeVelocity/1e3) +escvel.append(output.log.final.Saturn.EscapeVelocity/1e3) +escvel.append(output.log.final.George.EscapeVelocity/1e3) +escvel.append(output.log.final.Neptune.EscapeVelocity/1e3) + +shorelinex = [] +shorelinex.append(0.2) +shorelinex.append(60) + +shoreliney = [] +shoreliney.append(1e-6) +shoreliney.append(1e4) + +plt.xlabel('Escape Velocity [km/s]') +plt.ylabel('Normalized Cumulative XUV Flux') +plt.plot(shorelinex,shoreliney,color=vpl.colors.pale_blue) +plt.plot(escvel,fxuv,'o',color='k') +plt.xscale('log') +plt.yscale('log') +plt.ylim(1e-6,1e4) +plt.xlim(0.1,200) + +plt.annotate('Mercury',(2.7,9)) +plt.annotate('Venus',(10.2,2.6)) +plt.annotate('Earth',(11.5,0.5)) +plt.annotate('Mars',(5,0.2)) +plt.annotate('Jupiter',(60,0.05)) +plt.annotate('Saturn',(37,0.014)) +plt.annotate('Uranus',(22,0.0034)) +plt.annotate('Neptune',(24,0.0006)) + + +# Save figure +fig.savefig(path / f"CosmicShoreline.png", bbox_inches="tight", dpi=200) diff --git a/examples/CosmicShoreline/mars.in b/examples/CosmicShoreline/mars.in new file mode 100644 index 000000000..67b9ead8e --- /dev/null +++ b/examples/CosmicShoreline/mars.in @@ -0,0 +1,6 @@ +sName Mars +saModules atmesc +dMass -0.107 +dRadius -0.53202 +dSemi 1.52366290 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/mercury.in b/examples/CosmicShoreline/mercury.in new file mode 100644 index 000000000..d7a13d954 --- /dev/null +++ b/examples/CosmicShoreline/mercury.in @@ -0,0 +1,6 @@ +sName Mercury +saModules atmesc +dMass -0.0553 +dRadius -0.383 +dSemi 0.38709897 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/neptune.in b/examples/CosmicShoreline/neptune.in new file mode 100644 index 000000000..e7a477588 --- /dev/null +++ b/examples/CosmicShoreline/neptune.in @@ -0,0 +1,6 @@ +sName Neptune +saModules atmesc +dMass -17.15 +dRadius -3.883 +dSemi 30.07050641 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/saturn.in b/examples/CosmicShoreline/saturn.in new file mode 100644 index 000000000..979218df2 --- /dev/null +++ b/examples/CosmicShoreline/saturn.in @@ -0,0 +1,6 @@ +sName Saturn +saModules atmesc +dMass -95.16 +dRadius -9.449 +dSemi 9.53999265 +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/sun.in b/examples/CosmicShoreline/sun.in new file mode 100644 index 000000000..80e98f955 --- /dev/null +++ b/examples/CosmicShoreline/sun.in @@ -0,0 +1,18 @@ +# The host star +sName Sun # Body's name +saModules stellar # Modules to apply, exact spelling required + +# Output +saOutputOrder Age -Luminosity -LXUVStellar -Radius Temperature + +# Physical Parameters +dMass 1.00 # Mass, solar masses +dAge 5e7 +dObliquity 0 +dRotPeriod -30 # Rotation period, negative -> days +dRadGyra 0.5 # Radius of gyration (moment of inertia constant) + +# STELLAR Parameters +sStellarModel baraffe # Stellar evolution model: `baraffe` or `none` +dSatXUVFrac 1.e-3 # Saturation level of the XUV luminosity +dSatXUVTime -0.1 diff --git a/examples/CosmicShoreline/venus.in b/examples/CosmicShoreline/venus.in new file mode 100644 index 000000000..7fef281d6 --- /dev/null +++ b/examples/CosmicShoreline/venus.in @@ -0,0 +1,29 @@ +# Planet a parameters +sName Venus # Body's name +saModules atmesc # Modules to apply, exact spelling required + +# Physical Properties +dMass -0.815 # Mass, negative -> Earth masses +dRadius -0.9499 # Radius, negative -> Earth radii +dRotPeriod -243. # Rotation period, negative -> days +dObliquity 180. # Retrograde rotation +dRadGyra 0.5 # Radius of gyration (moment of inertia constant) + +# ATMESC Properties +dXFrac 1.0 # X-Ray/XUV absorption radius (fraction of planet radius) +dSurfWaterMass -3.0 # Initial surface water (Earth oceans) +dEnvelopeMass 0 # Initial envelope mass (Earth masses) +bHaltSurfaceDesiccated 0 # Halt when dry? +bHaltEnvelopeGone 0 # Halt when evaporated? +dMinSurfWaterMass -1.e-5 # Planet is desiccated when water content drops below this (Earth oceans) +sWaterLossModel lbexact +sPlanetRadiusModel none +bInstantO2Sink 1 +sAtmXAbsEffH2OModel bolmont16 + +# Orbital Properties +dSemi -0.723 # Semi-major axis, negative -> AU +dEcc 0.006772 # Eccentricity + +# Output +saOutputOrder Time -FXUV -CumulativeFXUV diff --git a/examples/CosmicShoreline/vpl.in b/examples/CosmicShoreline/vpl.in new file mode 100644 index 000000000..c53e2fbce --- /dev/null +++ b/examples/CosmicShoreline/vpl.in @@ -0,0 +1,25 @@ +sSystemName solarsystem # System Name +iVerbose 5 # Verbosity level +bOverwrite 1 # Allow file overwrites? + +# List of "body files" that contain body-specific parameters +saBodyFiles sun.in mercury.in venus.in earth.in mars.in $ + jupiter.in saturn.in george.in neptune.in + +# Input/Output Units +sUnitMass solar # Options: gram, kg, Earth, Neptune, Jupiter, solar +sUnitLength aU # Options: cm, m, km, Earth, Jupiter, solar, AU +sUnitTime YEARS # Options: sec, day, year, Myr, Gyr +sUnitAngle d # Options: deg, rad + +# Input/Output +bDoLog 1 # Write a log file? +iDigits 6 # Maximum number of digits to right of decimal +dMinValue 1e-10 # Minimum value of eccentricity/obliquity + +# Evolution Parameters +bDoForward 1 # Perform a forward evolution? +bVarDt 1 # Use variable timestepping? +dEta 0.01 # Coefficient for variable timestepping +dStopTime 4.5e9 # Stop time for evolution +dOutputTime 1e8 # Output timesteps (assuming in body files) diff --git a/src/atmesc.c b/src/atmesc.c index 979afe37e..69dfa656d 100644 --- a/src/atmesc.c +++ b/src/atmesc.c @@ -2002,6 +2002,17 @@ void fnPropsAuxAtmEsc(BODY *body, EVOLVE *evolve, IO *io, UPDATE *update, } +void fvCumulativeAtmEsc(BODY *body,EVOLVE *evolve,SYSTEM *system,double dDt,int iBody) { + if (evolve->bFirstStep) { + body[iBody].dFXUVCumulative = 0; + body[iBody].dFXUVLast = fdXUVFlux(body,iBody); + } else { + body[iBody].dFXUV = fdXUVFlux(body,iBody); + body[iBody].dFXUVCumulative += fdTrapezoidalArea(body[iBody].dFXUV,body[iBody].dFXUVLast,dDt); + body[iBody].dFXUVLast = body[iBody].dFXUV; + } +} + /** Assigns functions returning the time-derivatives of each variable to the magical matrix of function pointers. @@ -3454,6 +3465,20 @@ void WriteHRefODragMod(BODY *body, CONTROL *control, OUTPUT *output, fvFormattedString(cUnit, ""); } +void WriteCumulativeFXUV(BODY *body, CONTROL *control, OUTPUT *output, + SYSTEM *system, UNITS *units, UPDATE *update, int iBody, + double *dTmp, char **cUnit) { + *dTmp = body[iBody].dFXUVCumulative; + + if (output->bDoNeg[iBody]) { + *dTmp *= output->dNeg; + fvFormattedString(cUnit, output->cNeg); + } else { + *dTmp *= fdUnitsEnergyFlux(units->iTime,units->iMass,units->iLength); + fsUnitsEnergyFlux(units,cUnit); + } +} + /** Logs the molecular oxygen mixing ratio. @@ -3693,6 +3718,16 @@ void InitializeOutputAtmEsc(OUTPUT *output, fnWriteOutput fnWrite[]) { output[OUT_FXUV].iModuleBit = ATMESC; fnWrite[OUT_FXUV] = &WriteFXUV; + fvFormattedString(&output[OUT_CUMULATIVEFXUV].cName, "CumulativeFXUV"); + fvFormattedString(&output[OUT_CUMULATIVEFXUV].cDescr, "Cumulative XUV flux"); + fvFormattedString(&output[OUT_CUMULATIVEFXUV].cNeg, "W/m^2"); + output[OUT_CUMULATIVEFXUV].bNeg = 1; + output[OUT_CUMULATIVEFXUV].dNeg = 1; + output[OUT_CUMULATIVEFXUV].iNum = 1; + output[OUT_CUMULATIVEFXUV].iModuleBit = ATMESC; + fnWrite[OUT_CUMULATIVEFXUV] = &WriteCumulativeFXUV; + + fvFormattedString(&output[OUT_HESCAPEREGIME].cName, "HEscapeRegime"); fvFormattedString(&output[OUT_HESCAPEREGIME].cDescr, "Integer flag for H envelope escape regime"); diff --git a/src/atmesc.h b/src/atmesc.h index bb4392110..89a0e740a 100644 --- a/src/atmesc.h +++ b/src/atmesc.h @@ -171,6 +171,7 @@ void FinalizeUpdateMassAtmEsc(BODY *, UPDATE *, int *, int, int, int); #define OUT_HDIFFFLUX 1238 /**< The diffusion limited flux (the true flux in the diffusion regime) */ #define OUT_HREFODRAGMOD 1239 /**< Multiply by H ref flux to get H flux with drag of oxgyen */ #define OUT_KTIDE 1240 /**< Gravitational enhancement of mass loss */ +#define OUT_CUMULATIVEFXUV 1250 void InitializeOutputAtmEsc(OUTPUT *, fnWriteOutput[]); void InitializeOutputFunctionAtmEsc(OUTPUT *, int, int); @@ -212,6 +213,8 @@ void LogBodyAtmEsc(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UPDATE *, void fnForceBehaviorAtmEsc(BODY *, MODULE *, EVOLVE *, IO *, SYSTEM *, UPDATE *, fnUpdateVariable ***, int, int); void fnPropsAuxAtmEsc(BODY *, EVOLVE *, IO *, UPDATE *, int); +void fvCumulativeAtmEsc(BODY *,EVOLVE *,SYSTEM *,double,int); + double fdDSurfaceWaterMassDt(BODY *, SYSTEM *, int *); double fdDEnvelopeMassDt(BODY *, SYSTEM *, int *); double fdDEnvelopeMassDtBondiLimited(BODY *, SYSTEM *, int *); diff --git a/src/body.c b/src/body.c index 3c31b2757..59e9b93f3 100644 --- a/src/body.c +++ b/src/body.c @@ -1528,4 +1528,11 @@ void fdHabitableZoneKopparapu2013(BODY *body, int iNumBodies, double fdEffectiveTemperature(BODY *body,int iBody) { double dTeff = pow((body[iBody].dLuminosity/(4*PI*SIGMA*body[iBody].dRadius*body[iBody].dRadius)),0.25); return dTeff; +} + +double fdEscapeVelocity(BODY *body,int iBody) { + double dEscVel; + + dEscVel = sqrt(2*BIGG*body[iBody].dMass/body[iBody].dRadius); + return dEscVel; } \ No newline at end of file diff --git a/src/body.h b/src/body.h index f50c57b24..fa3ecf016 100644 --- a/src/body.h +++ b/src/body.h @@ -110,6 +110,7 @@ double CalcDynEllipEq(BODY *, int); void fdHabitableZoneKopparapu2013(BODY *, int, double *); double fdEffectiveTemperature(BODY*,int); +double fdEscapeVelocity(BODY *,int); // RB: Move diff --git a/src/evolve.c b/src/evolve.c index 832eb495c..e4cf87efa 100644 --- a/src/evolve.c +++ b/src/evolve.c @@ -13,6 +13,11 @@ #include "vplanet.h" +double fdTrapezoidalArea(double dY1,double dY2,double dDeltaX) { + double dArea = 0.5*dDeltaX*(dY1 + dY2); + return dArea; +} + void PropsAuxGeneral(BODY *body, CONTROL *control) { /* Recompute the mean motion, necessary for most modules */ int iBody; // Dummy counting variable @@ -50,6 +55,17 @@ void PropertiesAuxiliary(BODY *body, CONTROL *control, SYSTEM *system, } } +void fvCumulative(BODY *body,CONTROL *control,SYSTEM *system,double dDt) { + int iBody; + + for (iBody=0;iBodyEvolve.iNumBodies;iBody++) { + if (body[iBody].bAtmEsc) { + fvCumulativeAtmEsc(body,&control->Evolve,system,dDt,iBody); + } + } + +} + void CalculateDerivatives(BODY *body, SYSTEM *system, UPDATE *update, fnUpdateVariable ***fnUpdate, int iNumBodies) { int iBody, iVar, iEqn; @@ -669,6 +685,8 @@ void Evolve(BODY *body, CONTROL *control, FILES *files, MODULE *module, } } + fvCumulative(body,control,system,dDt); + fdGetUpdateInfo(body, control, system, update, fnUpdate); /* Halt? */ diff --git a/src/evolve.h b/src/evolve.h index 641e4dedc..6eec7a40b 100644 --- a/src/evolve.h +++ b/src/evolve.h @@ -16,6 +16,7 @@ /* @cond DOXYGEN_OVERRIDE */ +double fdTrapezoidalArea(double,double,double); void PropertiesAuxiliary(BODY *, CONTROL *, SYSTEM *, UPDATE *); void fdGetUpdateInfo(BODY *, CONTROL *, SYSTEM *, UPDATE *, fnUpdateVariable ***); diff --git a/src/output.c b/src/output.c index 846d5a821..354386f2c 100644 --- a/src/output.c +++ b/src/output.c @@ -128,6 +128,23 @@ void WriteDensity(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, } } +/* + E +*/ + +void WriteEscapeVelocity(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, + UNITS *units, UPDATE *update, int iBody, double *dTmp, + char **cUnit) { + + *dTmp = fdEscapeVelocity(body,iBody); + if (output->bDoNeg[iBody]) { + *dTmp *= output->dNeg; + fvFormattedString(cUnit, output->cNeg); + } else { + *dTmp *= fdUnitsTime(units->iTime) / fdUnitsLength(units->iLength); + fsUnitsVel(units, cUnit); + } +} /* * H @@ -1305,6 +1322,19 @@ void InitializeOutputGeneral(OUTPUT *output, fnWriteOutput fnWrite[]) { output[OUT_DENSITY].iModuleBit = 1; fnWrite[OUT_DENSITY] = &WriteDensity; +/* + E +*/ + + fvFormattedString(&output[OUT_ESCAPEVELOCITY].cName, "EscapeVelocity"); + fvFormattedString(&output[OUT_ESCAPEVELOCITY].cDescr, "Escape Velocity"); + fvFormattedString(&output[OUT_ESCAPEVELOCITY].cNeg, "m/s"); + output[OUT_ESCAPEVELOCITY].bNeg = 1; + output[OUT_ESCAPEVELOCITY].dNeg = 1; + output[OUT_ESCAPEVELOCITY].iNum = 1; + output[OUT_ESCAPEVELOCITY].iModuleBit = 1; + fnWrite[OUT_ESCAPEVELOCITY] = &WriteEscapeVelocity; + /* * H */ diff --git a/src/output.h b/src/output.h index 369fea3a5..f5c5d3a97 100644 --- a/src/output.h +++ b/src/output.h @@ -108,6 +108,7 @@ #define OUT_LOSTENG 691 #define OUT_LOSTANGMOM 692 +#define OUT_ESCAPEVELOCITY 693 /* @cond DOXYGEN_OVERRIDE */ diff --git a/src/vplanet.h b/src/vplanet.h index de3ed3395..ff52fb16a 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -262,6 +262,8 @@ struct BODY { double dThermTemp; /**< Thermosphere's temperature in Lehmer-Catling model */ double dAtmGasConst; /**< Atmosphere's gas constant in Lehmer-Catling model */ double dFXUV; /**< XUV Flux at planet's atmosphere */ + double dFXUVLast; /**< XUV at planet during last time step */ + double dFXUVCumulative; double dJeansTime; /**< Jeans timescale for atmospheric escape */ double dFlowTemp; /**< Temperature of the hydrodynamic flow */ double dRocheRadius; /**< Radius of the Roche lobe */