virtualplanetarylaboratory/vplanet

The Virtual Planet Simulator

VPLanet: The Virtual Planet Simulator

Overview

VPLanet is software to simulate planetary system evolution, with a focus on habitability. Physical models, typically consisting of ordinary differential equations, are coupled together to simulate evolution, from planetary cores to passing stars, for the age of a system. We strive for full transparency and reproducibility in our software, and this repository contains 1) the source code, 2) extensive documentation, 3) scripts and files to generate published figures and perform parameter sweeps, and 4) scripts to validate the current release. We can't claim we found life beyond the Earth with closed-source or unreliable software!

To get started, ensure you have clang/gcc installed and follow the Installation Guide. You can also watch videos on our YouTube channel on how to install and run VPLanet, as well as updates on recent results .

Modules

VPLanet currently consists of 13 functioning "modules," each containing a set of equations that simulates a specific physical process:

AtmEsc: Roche lobe overflow and thermal escape (energy-limited and radiation-recombination-limited) of an atmosphere, including water photolyzation, hydrogen escape, oxygen escape, and oxygen build-up.

Binary: Orbital evolution of a single circumbinary planet.

DistOrb: 2nd and 4th order semi-analytic models of orbital evolution outside of resonance.

DistRot: Evolution of a world's rotational axis due to orbital evolution and the stellar torque.

EqTide: Tidal evolution in the equilibrium tide framework.

Flare: Flare frequency distribution and flare XUV luminosity evolution in low-mass stars.

GalHabit: Evolution of a wide orbit due to the galactic tide and impulses from passing stars (including radial migration).

MagmOc: Thermal and geochemical evolution of a magma ocean.

POISE: Energy balance climate model including dynamic ice sheets and lithospheric compression/rebound.

RadHeat: Radiogenic heating in a world's core, mantle, and crust.

SpiNBody: N-body integrator for the evolution of a system of massive particles.

Stellar: Evolution of a star's bolometric and XUV luminosity, temperature, radius, and mass concentration. Also includes magnetic braking and stellar wind spin-down.

ThermInt: Thermal interior evolution, including magnetic fields, for planets undergoing plate tectonics or stagnant lid evolution.

Many of these modules can be combined together to simulate numerous phenomena and feedback loops in planetary systems.

Resources

The examples/ directory contains input files and scripts for generating the figures in Barnes et al. (2020) and subsequent publications. The Manual/ directory contains the pdf of Barnes et al. (2020) that describes the physics of the first 11 modules, validates the software against observations and/or past results, and includes figures from the examples/ directory.

An ecosystem of support software is also publicly available in other repositories of the Virtual Planetary Laboratory. The vplot package is both a command line tool to quickly plot the evolution of a single simulation and a Pythom module for generating publication-worthy figures. The VSPACE script generates input files for a parameter space sweep, which can then be performed on an arbitrary number of cores with MultiPlanet. For large parameter sweeps, an enormous amount of data can be generated, which can slow analyses. To overcome this barrier, the BigPlanet code can both compress datasets into HDF5 format, including statistics of an integration, and tools to facilitate plotting. These three scripts can be executed from the command line to seamlessly perform parameter sweeps. These Python scripts are optimized for anaconda distributions.

The "pip-install" badge indicates if the latest executables are available for installation via pip for the Python distributions and operating systems for which we perform unit tests (see below). (Note: This badge is currently broken due to an unknown GitHub issue, so if you pip install, you may want to check your VPLanet version, which is printed near the top of the log file that is generated when you perform a simulation.)

Code Integrity

We are committed to maintaining a stable tool for scientists to analyze planetary system. Behind the scenes, the VPLanet team maintains code integrity through automatic checks at every merge into the main branch. You can see the status of these checks via the badges above. Currently we perform "Unit Tests" for the initial and final conditions across an orthogonal set of "Test Sims" (simulations), with the numbers of tests for each shown via badges. We perform the tests across all permutations of operating systems and Python version shown by the badges.

We also check for memory violations via valgrind's memcheck tool ("memcheck") and address sanitizer ("sanitizer"), check for overflow, invalid operation, and divide-by-zero floating point exceptions ("floating-point"), and test if all the examples work across the operating system and Python versions listed after the "examples" badge.

The percentage of the lines of code that are executed by the unit tests is shown with the "codecov" badge, with details available at our Codecov account.

Community

VPLanet is a community project. We're happy to take pull requests; if you want to create one, please issue it to the dev branch. The documentation includes tutorials on adding new features and modules. VPLanet is a platform for planetary science that can grow exponentially, either by adding new physics or by adding competing models for clean comparisons.

A list of additional GitHub repositories with VPLanet examples can be found here.

If you have questions or are running into issues, you can read or post to a Discussion.

If you believe you have encountered a bug, please raise an issue using the Issues tab at the top of this page.

Acknowledgments

If you use this code to generate results used in any publication or conference contribution, please cite Barnes, R. et al. (2020), PASP, 132, 24502.

VPLanet development has been supported by NASA grants NNA13AA93A, NNX15AN35G, 80NSSC17K048, 13-13NAI7_0024, 80NSSC20K0229, and 80NSSC20K0261. We also acknowledge support from the University of Washington and the Carnegie Institute for Science.

Enjoy!

© 2018-2024 The VPLanet Team.

vplanet's Issues

Derivative values in the log file and .forward files don't match

Currently the initial derivatives are not calculated until after the log file is written. The values in the first line .forward files are correct.

the line numbers in quickstart's input files interfere when new users try to copy&paste its content.

Most new VPLanet users (like myself) will probably start with the quickstart in the docs. The example input files like vpl.in there contains line numbers, which get in the way when people copy&paste the file content to try out their first simulation.

I suggest to remove the file numbers from the quickstart doc. Happy to do this fix if it finds approval!

Add a CONTRIBUTING file

Add a file explaining how users can submit PRs with proposed additions to VPLANET.

Simulation of nested orbiting objects

While adjusting a planet's inputs to take into account any orbiting bodies works in many situations, it does not take into account any tidal forces between the planet and the orbiting body or the orbiting body and the stellar body the planet itself is orbiting around. Basically, I'm looking for a way to easily create an array of orbiting bodies around any object that is already in its own orbit.

Water Loss in HZ is broken

The AtmEsc switch bStopWaterLossInHZ is currently broken. Whether set or not, water will be photolyzed and hydrogen will escape in the HZ.

Two functions for orbit conversions

Spinbody and distorb both have functions to convert between orbital elements and Cartesian coordinates. These should be combined into 1 and moved into system.c.

EarthClimate/makeplot.py is failing

This seems to be sensitive to the version of some package, but I can't figure out which package(s) are the issue.

In examples/EarthClimate, python makeplot.py fails with this error:

Traceback (most recent call last):
File "/home/rdeitrick/vplanet_test_sims/EarthClimate_copy0/makeplot.py", line 415, in
comp2huybers("Earth", show=False)
File "/home/rdeitrick/vplanet_test_sims/EarthClimate_copy0/makeplot.py", line 108, in comp2huybers
c = plt.contourf(
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/pyplot.py", line 2492, in contourf
__ret = gca().contourf(
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/init.py", line 1414, in inner
return func(ax, *map(sanitize_sequence, args), **kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/axes/_axes.py", line 6319, in contourf
contours = mcontour.QuadContourSet(self, *args, **kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/contour.py", line 812, in init
kwargs = self._process_args(*args, **kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/contour.py", line 1446, in _process_args
x, y, z = self._contour_args(args, kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/contour.py", line 1490, in _contour_args
z = ma.masked_invalid(z, copy=False)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/numpy/ma/core.py", line 2360, in masked_invalid
res = masked_where(~(np.isfinite(a)), a, copy=copy)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/numpy/ma/core.py", line 1936, in masked_where
cond = mask_or(cond, a._mask)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/numpy/ma/core.py", line 1746, in mask_or
return make_mask(m1, copy=copy, shrink=shrink, dtype=dtype)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/numpy/ma/core.py", line 1640, in make_mask
result = _shrink_mask(result)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/numpy/ma/core.py", line 1549, in _shrink_mask
if m.dtype.names is None and not m.any():
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/astropy/units/quantity.py", line 2001, in any
raise TypeError(
TypeError: cannot evaluate truth value of quantities. Evaluate array with q.value.any(...)

The weird thing is that the script still works in an old environment, but not when I create a new environment, even if I specify the same versions of matplotlib and astropy. I suspect another package is the issue, but can't figure out which one.

I've tried adding .value to unitful quantities inside the contourf call. This removes the above issue but generates a new issue with vplot:

Traceback (most recent call last):
File "/home/rdeitrick/vplanet_test_sims/EarthClimate_copy0/makeplot.py", line 415, in
comp2huybers("Earth", show=False)
File "/home/rdeitrick/vplanet_test_sims/EarthClimate_copy0/makeplot.py", line 255, in comp2huybers
plt.savefig(path / f"EarthClimateMilankovitch.{ext}")
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/pyplot.py", line 1119, in savefig
res = fig.savefig(*args, **kwargs) # type: ignore[func-returns-value]
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/figure.py", line 3390, in savefig
self.canvas.print_figure(fname, **kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/backend_bases.py", line 2187, in print_figure
result = print_method(
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/backend_bases.py", line 2043, in
print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/backends/backend_pdf.py", line 2807, in print_pdf
self.figure.draw(renderer)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/matplotlib/artist.py", line 39, in draw_wrapper
return draw(artist, renderer, *args, **kwargs)
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/vplot-1.0.5-py3.9.egg/vplot/figure.py", line 353, in draw
self._add_labels()
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/vplot-1.0.5-py3.9.egg/vplot/figure.py", line 179, in _add_labels
unit, _, label, physical_type = _get_array_info(
File "/home/rdeitrick/anaconda3/envs/test_a6/lib/python3.9/site-packages/vplot-1.0.5-py3.9.egg/vplot/figure.py", line 33, in _get_array_info
physical_type = array.unit.physical_type.title()
AttributeError: 'PhysicalType' object has no attribute 'title'

So far, the only way I've gotten around this entirely is to remove vplot completely, and add .value to all quantities getting passed to plotting functions... but there must be a way of getting this working properly with vplot, surely?

Does Virtual Planetary work for wsn

thanks for you efforts,

i have a question, Does Virtual Planetary work for Wireless sensor networks to save energy and extend lifetime of the networks?

imp module in setup.py

The imp module (imported by setup.py) is now fully removed from Python 3.12. The good news is that setup.py doesn't actually seem to be using it! That import imp line can probably be removed. But someone should double check that in case I am missing something.

Initial ice sheets heights shouldn't be negative!

The option dinitIceSheetHeight in POISE allows negative arguments, which shouldn't be allowed. It'd be better to change it so that if the user inputs a negative value, it defaults to meters since the ice sheet height is positive definite.

Problems with distrot near coplanar retrograde obliquity

I have found that beginning near 180 deg obliquity causes distrot to hit a singularity described by O. Neron de Surgy and J. Laskar (1997). Attached is an example trial near this condition. Also the orbital integration can be halted due to a large max eccentricity when this occurs. Given the integration is variable step, it seems odd that there should be such a growth in eccentricity for the coplanar retrograde case.
Inc[180].zip

Read in input files only once

Right now the code opens and read the input files for every option. It'd be better to open them once and read the entire file into an array of chars, one line per element.

STELLAR module: issues with output density

I found two issues with the density output using the STELLAR module:

1. Density is listed under 'system properties' in the output log file instead of 'stellar parameters'
2. The log file says that the density is in solar units, but it's actually a factor of ~12 off from solar (see plot)

plot:

infiles and log:
star.txt
vpl.txt
Trappist.log

Implement the LSODA ODE solver in place of RK4

VPLanet would probably be ~10x faster, and a more flexible if it used a production quality ODE solver. I think LSODA (this is link to C version) is a good option. It automatically detects stiff ODEs, and will switch methods to be as efficient as possible. This would allow VPlanet to implement modules with stiff ODEs, like chemical reactions in the atmosphere. It looks like the POISE module already introduces stiff ODEs, which RK4 has a hard time dealing with.

CTL Model can fail for 1 star, 2 planet systems with DISTORB

In 1 star, 2 planet systems where EQTIDE and DISTORB are active, and if tides are significant on the outer planet, the evolution can fail to conserve energy. The evolution works if EQTIDE is not active on the the outer perturbing body, however.

From David Fleming.

mem-check: dirty

The MagmOc module has not been swept for memory leaks yet. The tests pass and the results are consistent with previously published results, but please use cautiously.

Maximum number of planets?

Is there a limit to the number of simulated planets per system?

I'm getting strange behavior with changing units and unrecognized keywords when using many planets with atmesc, even though the input files are (nearly) identical :

WARNING: dEnvelopeMass < 0 in file p_0.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_1.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_2.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_3.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_4.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_5.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_6.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_7.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_8.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_9.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_10.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_11.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_12.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_13.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_14.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_15.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_16.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_17.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_18.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_19.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_20.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_21.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_22.in, units assumed to be Earth.
WARNING: dEnvelopeMass < 0 in file p_23.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_24.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_25.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_26.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_27.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_28.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_29.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_30.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_31.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_32.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_33.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_34.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_35.in, units assumed to be .in.
WARNING: dEnvelopeMass < 0 in file p_36.in, units assumed to be .in.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_28.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_28.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_29.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_29.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_30.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_30.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_31.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_31.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_32.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_32.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_33.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_33.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_34.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_34.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_35.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_35.in, line 13.
ERROR: Unrecognized option "bHaltEnvelopeGone" in p_36.in, line 11.
ERROR: Unrecognized option "dAtmXAbsEffH" in p_36.in, line 13.

infiles.zip

Changing Stellar Age

Hej, I'm currently doing my master's thesis on early terrestrial evolution, and i am using the vplanet code, specifically the magmoc_earth example.
I've been chancing a few variables in the code to see how that migth affect the evolution of the planet, and one of those variables is the dSemi - I'm keeping everything else constant and in mostly in the default values.
I came across an issue with my data as the Sun's age in this example is sat to 500 kyr (very young), and since the Sun was more volatile at the age, I changed the value to be more appropriate to when the formation/accretion of Earth happened (at about 60Myr).
Now I have an issue when my planet is at dSemi 1 (au), and after running the program get a message of 'core dumped' as a result.
Have I made a mistake somewhere? Or is what I'm getting now a more realistic view/approach to when the accretion of earth happened?
Also, the solidification time is lessened by about 200kyr compared to the Hamano plot used in the example

.

Spurious behavior forced rotation and circular orbits?

If a user sets bForceEqSpin to 1, dEccentricity to 0, and a non-zero value of obliquity, then the obliquity evolves secularly to 180 degrees. I don't know if this is physical or not, but looks suspicious.

Thanks to John Ahlers for pointing this out.

Orbital and\or rotational evolution with mass loss

Congratulations on the first release of this great program!

I have a question: is it possible to simulate orbital and rotational evolution taking into account mass loss\grow in stars and planets(in the isotropic approximation)?

Thank you very much for your attention!

Orbital parameters for central body are nonsense in log file

Central bodies do not have orbital parameters (they can be orbited by many bodies), so orbital parameters are reported as -1. There might be a better way to do this...

Issue with composition

in magmoc.h file line 76, the mass fraction of MgO in the mantle is defined as 0.0367, which is very low for a mantle peridotite. I think the MASSFRACMGO is meant to be 0.367 instead? Which is more consistent with papers.

I have no real experience programming with c, and cannot get the 'makefile' to compile with the new values.

extra q: is there any way I can make the core a parameter option in the body files also?
I'm doing my master's in early terrestrial mo evolution, and would like the opportunity to change the core size. It can be relative to Earth, like with planet mass and radius, if that will make it easier for me.
Thank you for your troubles.

iVerbose only works in primary file

Tested both in my own personal branch and vplanet-main, the iVerbose input option only works if placed in the primary file (vpl.in). If placed elsewhere, then an error is output stating for example

ERROR: Unrecognized option "iVerbose" in star.in

I currently don't have the time to fix it myself, which why I'm posting here as a reminder to everyone.

Memcheck-dirty

Currently the derivatives are not calculated prior to logging so their initial values may not be correct. The first line in the output (.forward) files is correct, so use those values.

Add power/surface energy flux due to secular cooling in thermint

Currently the code does not include the power or surface energy flux generated by non-radiogenic and non-tidal sources. Thus the sum of RadPowerTotal and PowerEqtide do not sum to PowerTotal, which can be confusing. New output options that report this energy source should be added. Thanks to Stephen Kane for identifying this issue.

Enumerate option/output macros

Assigning the options and output parameters identifiers (OPT_, and OUT_) via #define is poor style and hinders debugging. The should be assigned via enum. This would, however, add the complication that they could not be part of the OPTIONS or OUTPUT structs and require considerable editing of the entire code.

Min/Max Values for Gaussian Distributions with VSPACE

If you set one parameter's minimum and/or maximum values for a Gaussian sampling distribution with VSPACE, those limits will be imposed on other parameters with a Gaussian distribution, too. As a quick hack, set all Gaussian min/max values and VSPACE will work.

Problem with DeltaT in GalHabit

I have found that time step DeltaT goes down as 1/n for each new calculation in a simulation when running GalHabit module. For example, the time step went from 1000yrs in the beginning down to 0.2 at the end of a 5Gyr simulation.

Derivative Values in Log Files and .forward Files Don

make legacy fails

When running make legacy there is the following error with declared variables in POISE:

gcc -o bin/vplanet src/*.c -lm -DGITVERSION=\"v1.0.0-1643-gc2fd79de8915981786f11adaec39f9965517a557\"
src/poise.c: In function ‘fvPoiseSeasonalInitialize’:
src/poise.c:6792:3: error: ‘for’ loop initial declarations are only allowed in C99 mode
   for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
   ^
src/poise.c:6792:3: note: use option -std=c99 or -std=gnu99 to compile your code
src/poise.c: In function ‘fvFluxesByLatitude’:
src/poise.c:6886:2: error: ‘for’ loop initial declarations are only allowed in C99 mode
  for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
  ^
src/poise.c:6902:5: error: ‘for’ loop initial declarations are only allowed in C99 mode
     for (int jLat = 0; jLat < body[iBody].iNumLats; jLat++) {
     ^
src/poise.c: In function ‘fvCalculateSeaIce’:
src/poise.c:7074:3: error: ‘for’ loop initial declarations are only allowed in C99 mode
   for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
   ^
make: [legacy] Error 1 (ignored)

this apparently happens on machines where the standard is not set to C99. Either this issue should be fixed in POISE or the makefile should specify std

-gcc -o bin/vplanet src/*.c -std=c99 -lm -DGITVERSION=\"$(GITVERSION)\"

Printing of output parameters should be limited to those that are requested

When the code outputs parameters, it loops over all possible values, MODULEOPTEND, which is inefficient. Instead the code should record the number of outputs per body and just loop over them.

Issues with ice sheet discretization in POISE

The ice sheet flow equation is not properly discretized to conserved mass. Admittedly, this is a small error compared to other aspects of the EBM, but could still be improved.

• First, the array daYBoundary should probably be renamed for clarity. Something like daDeltaYBoundary. We'll also need to introduce a second array, daDeltaYCenter. daDeltaYBoundary will be Delta y, centered on the boundaries, i.e., the difference in y coordinate between cell centers. daDeltaYCenter will then be centered on the cell centers, i.e, the difference in y coordinate between cell boundaries.

• Second, the mass conservation part. Everywhere that terms involving Delta y^2 appear, these need to use a boundary value and a centered value. Terms that look like:
  (body[iBody].daIceSheetDiff[iLat] / (body[iBody].daYBoundary[iLat] * body[iBody].daYBoundary[iLat])
  should instead look like:
  (body[iBody].daIceSheetDiff[iLat] / (body[iBody].daDeltaYBoundary[iLat] * body[iBody].daDeltaYCenter[iLat])
  etc. I think it will be relatively straightforward to replace the second reference to the boundary value with the center value in each term. This goes for every value of daIceSheetMat and daIcePropsTmp.

POISE should halt if rotational frequency is close to pericenter angular frequency

A halt needs to be added to ensure POISE is not used at high eccentricity.

Incorrect shared axis labels in examples

v1.0.2 is currently not working with shared axis labels. See VirtualPlanetaryLaboratory/vplot#8.

Move get_args.py into vplot

The examples assume the file get_args.py is available. If you move an example directory to another location, the makeplot.py file will not work unless you move get_args.py with it.

The long-term solution is probably to move the get_args function into vplot.

Issue with vplot on Window's directory structure

Hi, VPlanet team! Looking forward to this week's workshop.

I was testing out my installation and found a small bug when running
vplot after a run on a Windows 10 install (conda install w/ Python 3.8).

I received this traceback (<PATH> is just my local system path to where I installed vplanet).

Traceback (most recent call last):
  File "C:\<PATH>\vpl\Scripts\vplot-script.py", line 33, in <module>
    sys.exit(load_entry_point('vplot', 'console_scripts', 'vplot')())
  File "C:\<PATH>\vpl\lib\site-packages\vplot-1.0.2-py3.8.egg\vplot\command_line.py", line 37, in _entry_point
    auto_plot(**args.__dict__)
  File "C:\<PATH>\vpl\lib\site-packages\vplot-1.0.2-py3.8.egg\vplot\auto_plot.py", line 66, in auto_plot
    output = get_output(sysname=sysname, path=path)
  File "<PATH>\vplanet\vplanet\output.py", line 300, in get_output
    log = get_log(sysname=sysname, path=path, units=units)
  File "<PATH>\vplanet\vplanet\log.py", line 308, in get_log
    raise ValueError("Error processing line {} of {}: ".format(i, lf) + str(e))
ValueError: Error processing line 3 of solarsystem.log: too many values to unpack (expected 2)

The issue is with this line in the solarsystem.log:
Executable: C:\<PATH>\vpl\Scripts\vplanet
It is getting confused with the 2nd ":". By removing that from the text file it works just fine.

Yay for cross-platform directory non-conformity!

Semi-related: I am also getting an "Unknown" for that file's version tag, not sure if that is as intended or not.

Version: Unknown

Cheers

Forced orbital/rotational evolution with POISE

POISE can include arbitrary oscillations of eccentricity and obliquity, as well as axial precession self-consistently from the star. However, these properties are calculated in ForceBehaviorPoise, preventing them from determining the timestep. They should become primary variables (exclusive to POISE to prevent overlap with HECC, KECC, XOBL, YOBL, and ZOBL) that are part of The Matrix to ensure self-consistent evolution. Right now a user could set a timestep that is too large and the resultant climate would be wrong.

Memcheck-dirty

There are cases of uninitialized memory with the habitable zone calculations, but the tests are passing. If you use HZ limits, double check that they appear correct.

Better reporting of tidal parameters

If a user selects the DB15 model and includes rotational properties, a Warning should be printed.

The help message for sTideModel does not include the DB15 option.

IceLatitudes Issue With Unit Conversion (POISE)

When calculating the latitude for the ice sheets on when there is no ice present, the code sets it to 100 deg, however it interprets it as 100rad and when converted to degrees it shows a value of 5729.5779514719952203 (equivalent to pi/180 * 100).

Tests for outputs that begin with a number

The maketest script skips any output parameter that begins with a number. This may be due to an issue in how vplot gathers outputs. Probably the thing to do is change the output names to begin with a letter.

Add output column headers to .forward files

Gijs Mulders suggested that it would be useful to include commented-out headers in output files. These will be ignored by most python I/O functions, like numpy's genfromtxt of pandas' read_tables, and can make output more human-readable.

It could look like this:

Time Eccentricity Semi-major Axis

0.0    0.3                1.0

and so on.

memcheck dirty

The upgrade to v2.0 has created a few memory leaks. Tests are passing, but if you push the code into something that is not currently tested, please use caution.