virtualplanetarylaboratory / vplanet Goto Github PK
View Code? Open in Web Editor NEWThe Virtual Planet Simulator
License: MIT License
The Virtual Planet Simulator
License: MIT License
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.
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...
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.
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:
0.0 0.3 1.0
and so on.
I found two issues with the density output using the STELLAR module:
infiles and log:
star.txt
vpl.txt
Trappist.log
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.
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.
A halt needs to be added to ensure POISE is not used at high eccentricity.
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.
The AtmEsc switch bStopWaterLossInHZ is currently broken. Whether set or not, water will be photolyzed and hydrogen will escape in the HZ.
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?
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.
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.
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.
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.
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.
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!
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.
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.
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.
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.
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.
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.
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
.
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 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.
v1.0.2 is currently not working with shared axis labels. See VirtualPlanetaryLaboratory/vplot#8.
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.
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.
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.
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
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.
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
Add a file explaining how users can submit PRs with proposed additions to VPLANET.
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.
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.
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!
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?
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.
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)\"
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).
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.