Create models of rotating and non-rotating planets by solving the differential
equations for hydrostatic equilibrium, and create initial conditions for e.g.
smoothed particle hydrodynamics (SPH) simulations by placing particles that
precisely match the planet's profiles.
See the tutorial.ipynb
notebook for an introductory tutorial and examples
(https://github.com/srbonilla/WoMa), with additional documentation in this
file, below. Please see the class and function docstrings for full details.
Presented in Ruiz-Bonilla et al. (2020, MNRAS 500:3).
Includes SEAGen (https://github.com/jkeger/seagen;
Kegerreis et al. 2019, MNRAS 487:4)
with modifications for spinning planets.
Sergio Ruiz-Bonilla: sergio.ruiz-bonilla@durham.ac.uk
Jacob Kegerreis: jacob.kegerreis@durham.ac.uk
Tom Sandnes: thomas.d.sandnes@durham.ac.uk
If you find any bugs, potential improvements, or features worth adding, then
please let us know!
Files
tutorial.ipynb
Jupyter notebook tutorial and examples.woma/
Code directory.
main.py
The main program classes and functions.data/
Data for equation of state (EoS) tables.eos/
EoS and temperature-density relations.spherical_funcs/
Functions for spherical profiles.spin_funcs/
Functions for spinning profiles.misc/
Miscellaneous functions and constants.
README.md
This file. General info and documentation.LICENSE.txt
GNU general public license v3+.setup.py
, setup.cfg
, MANIFEST.in
PyPI package files.
Installation and requirements
- Install the package with
pip install woma
, or downloaded directly from https://pypi.org/project/woma/ - Python 3 (tested with 3.6.9)
+ seagen>=1.5, numpy, numba>=0.50.1, h5py
Documentation
See the tutorial.ipynb
notebook for a tutorial and examples.
This documentation summarises the different primary options available.
Full documentation is provided in the class and function docstrings,
particularly for the various utilities (e.g. equations of state functions)
that are not the focus of the primary initial-conditions topics
in the tutorial and here.
Most functions take an optional verbosity
argument
that controls the amount of printed information.
Set 0
for no printing, 1
for standard output, or 2
for extra details.
Contents
- Spherical profiles
- Spinning profiles
- Particle placement
- (Bonus) Impact initial conditions
Notation etc.
- Formatted with black.
- Arrays are explicitly labelled with a prefix
A1_
, or An_
for an
n
-dimensional array. numba
provides great speed increases, but places some odd requirements and
limitations that lead us to some awkward or ugly but functional approaches...
1. Spherical profiles
All profiles require the temperature and either the pressure or density at the
surface: T_s
and P_s
or rho_s
,
plus the material and temperature-density relation for each layer (see below):
A1_mat_layer
and A1_T_rho_type
.
For a fixed-entropy relation, only the surface density is needed.
The optional parameter num_prof
sets the number of profile integration steps.
Default 1000.
Equations of state (EoS)
It is important to check that the EoS you use are appropriate
for your application.
So far, we have implemented several Tillotson, ANEOS, SESAME,
and Hubbard & MacFarlane (1980) materials, with more on the way.
Custom materials in SESAME-style tables can also be provided.
Materials are set for each layer with A1_mat_layer
using the following material names,
which are converted internally into material IDs,
set by a base type ID (multiplied by 100) plus a minor type:
- Tillotson (Melosh 1989; Benz & Asphaug 1999):
1
- Iron:
Til_iron
: 100
- Granite:
Til_granite
: 101
- Water:
Til_water
: 102
- Basalt:
Til_basalt
: 103
- Hubbard & MacFarlane (1980):
2
- Hydrogen-helium atmosphere:
HM80_HHe
: 200
- Ice H20-CH4-NH3 mix:
HM80_ice
: 201
- Rock SiO2-MgO-FeS-FeO mix:
HM80_rock
: 202
- SESAME (Lyon & Johnson 1992) and others in similar formats:
3
- Iron (2140):
SESAME_iron
: 300
- Basalt (7530):
SESAME_basalt
: 301
- Water (7154):
SESAME_water
: 302
- Senft & Stewart (2008) water:
SS08_water
: 303
- Haldemann, J. et al. (2020) AQUA:
AQUA
: 304
- Chabrier, G. et al. (2019) Hydrogen:
CMS19_H
: 305
- Chabrier, G. et al. (2019) Helium:
CMS19_He
: 306
- Chabrier & Debras (2021) H/He mixture Y=0.245 (Jupiter):
CD21_HHe
: 307
- ANEOS (in SESAME-style tables):
4
- Forsterite (Stewart et al. 2019):
ANEOS_forsterite
: 400
- Iron (Stewart, zenodo.org/record/3866507):
ANEOS_iron
: 401
- Fe85Si15 (Stewart, zenodo.org/record/3866550):
ANEOS_Fe85Si15
: 402
- Custom (in SESAME-style tables):
9
- User-provided custom material(s):
900
, 901
, ..., 904
These are defined in woma/misc/glob_vars.py
,
including the file paths for custom tables.
Temperature-density relations
These relations are set for each layer with A1_T_rho_type
:
"adiabatic"
: Adiabatic, only available for some EoS."power=a"
where a
is a float: A power law T ~ rho^a
.
So e.g. "power=0"
for isothermal."entropy=s"
where s
is a float: A fixed specific entropy (J K^-1 kg^-1).
Similar to adiabatic, but uses this entropy directly instead of deriving it
from the temperature and density or pressure.
Profile generation
There are several options for which additional parameters are set and which
unknowns are found, depending on the number of layers in the planet.
Most of these functions are simple iterative bisection searches over the unknown
parameter(s) to find a valid planet in hydrostatic equilibrium that satisfies
the set attribute values.
The additional function arguments like R_max
set things like the upper bound
for an iteration so usually do not need to be precise.
Optional arguments for these functions (in addition to the verbosity
) set:
tol
: The tolerance for finding unknown parameters as a fractional
difference between two consecutive iterations. Default usually 0.001,
depending on the method.num_attempt
: The maximum number of iteration attempts if the tolerance has
still not been reached. Default usually 40, depending on the method.
If the outer radii or masses of some but not all layers are required as inputs,
then the unknown elements in the input arrays can be left as None
, e.g.:
A1_R_layer = [3.1415, None]
or A1_M_layer = [None, 1.23, 4.56]
.
1 layer
gen_prof_L1_find_R_given_M()
, requires:
- Total mass:
self.M
- Maximum radius:
R_max
gen_prof_L1_find_M_given_R()
, requires:
- Total radius:
self.R
- Maximum mass:
M_max
2 layers
gen_prof_L2_find_R1_given_M_R()
, requires:
- Total radius:
self.R
- Total mass:
self.M
gen_prof_L2_find_M_given_R_R1()
, requires:
- Total radius:
self.R
- Layer 1 outer radius:
self.A1_R_layer[0]
- Maximum mass:
M_max
gen_prof_L2_find_R_given_M_R1()
, requires:
- Total mass:
self.M
- Layer 1 outer radius:
self.A1_R_layer[0]
- Maximum radius:
R_max
gen_prof_L2_find_R_R1_given_M1_M2()
, requires:
- Layer 1 and 2 masses:
self.A1_M_layer
- Minimum and maximum radii:
R_min
, R_max
3 layers
gen_prof_L3_find_M_given_R_R1_R2()
, requires:
- Layer 1, 2, and 3 outer radii:
self.A1_R_layer
gen_prof_L3_find_R1_given_M_R_R2()
, requires:
- Total mass:
self.M
- Total radius:
self.R
- Layer 2 outer radius:
self.A1_R_layer[1]
gen_prof_L3_find_R2_given_M_R_R1()
, requires:
- Total mass:
self.M
- Total radius:
self.R
- Layer 1 outer radius:
self.A1_R_layer[0]
gen_prof_L3_find_R_given_M_R1_R2()
, requires:
- Total mass:
self.M
- Layer 1 and 2 outer radii:
self.A1_R_layer[0]
, [1]
gen_prof_L3_find_R1_R2_given_M_R_I()
, requires:
Adding layers
gen_prof_given_inner_prof()
:
After generating an initial planet, a new layer can be added on top by
integrating outwards. Requires:
- Name of the material in the new layer:
mat
- Temperature-density relation in the new layer:
T_rho_type
- Minimum density at which the new layer will stop:
rho_min
- Minimum pressure at which the new layer will stop:
P_min
Additional parameters
See the main.py
and other docstrings for full details.
The class's num_prof
parameter sets the number of radial profile steps,
while the profile generating functions take arguments like tol
and/or num_attempt
that control the convergence criterion and
maximum number of iterations to do to find the unknown parameters.
2. Spinning profiles
See tutorial.ipynb
for the main usage:
spherical_planet = woma.Planet( . . . )
spin_planet = woma.SpinPlanet(
planet = spherical_planet,
period = 24,
)
The output attributes available from the spin_planet
object
are documented in the SpinPlanet
class docstring in woma/main.py
.
The primary outputs are the arrays of properties of the nested spheroids,
including their equatorial and polar radii (semi-major and semi-minor axes),
A1_R
and A1_Z
,
and for example their masses, densities, pressures, and temperatures,
A1_m
, A1_rho
, A1_P
, and A1_T
.
Additional parameters are similar to those in the spherical case mentioned
above. See the docstrings for full details.
3. Particle placement
See tutorial.ipynb
for the main usage:
planet = woma.Planet( . . . )
planet = woma.SpinPlanet( . . . )
N_particles = 1e7
N_ngb = 48
particles = woma.ParticleSet(planet, N_particles, N_ngb=N_ngb)
The output attributes available from the particles
object
are documented in the ParticlePlanet
class docstring in woma/main.py
.
Particles can be saved to a SWIFT-style HDF5 file with the save()
method.
See its docstring and save_particle_data()
in woma/misc/io.py
for details.
The specific entropies of the particles can be added to the object using the
set_material_entropies()
or calculate_entropies()
methods.
4. Impact initial conditions
One of the motivations for WoMa's development was to create initial conditions
for the modelling of planetary giant impacts. (Check out the open-source SWIFT
code at www.swiftsim.com and see e.g. http://icc.dur.ac.uk/giant_impacts/ for
more info about these applications.)
Therefore, we include here some simple utilities for setting up the initial
conditions for an impact scenario between two planets, as described in Appendix
A of Kegerreis et al. (2020), ApJ 897:161.
See the tutorial.ipynb
notebook for an example and the docstrings of
impact_pos_vel_b_v_c_r()
and impact_pos_vel_b_v_c_t()
in
woma/misc/utils.py
for the full documentation.