Erebus.jl
Documentation for Erebus.jl
Erebus.Q_radiogenic
Erebus.add_vrk4
Erebus.apply_insulating_boundary_conditions!
Erebus.apply_subgrid_stress_diffusion!
Erebus.apply_subgrid_temperature_diffusion!
Erebus.assemble_gravitational_lse!
Erebus.assemble_hydromechanical_lse!
Erebus.assemble_thermal_lse!
Erebus.backtrace_pressures_rk4!
Erebus.calculate_radioactive_heating
Erebus.compute_Aϕ!
Erebus.compute_adiabatic_heating!
Erebus.compute_basic_node_properties!
Erebus.compute_displacement_timestep
Erebus.compute_fluid_velocities!
Erebus.compute_gibbs_free_energy
Erebus.compute_gravity_solution!
Erebus.compute_kfluidm
Erebus.compute_ksolidm
Erebus.compute_marker_properties!
Erebus.compute_molarfraction!
Erebus.compute_nodal_adjustment!
Erebus.compute_p_node_properties!
Erebus.compute_reaction_constant
Erebus.compute_relative_enthalpy
Erebus.compute_rhocpfluidm
Erebus.compute_rotation_rate!
Erebus.compute_shear_heating!
Erebus.compute_stress_strainrate!
Erebus.compute_thermochemical_iteration_outcome
Erebus.compute_thermodynamic_xfer!
Erebus.compute_velocities!
Erebus.compute_vx_node_properties!
Erebus.compute_vy_node_properties!
Erebus.compute_Δtreaction
Erebus.define_markers!
Erebus.distance
Erebus.dot4
Erebus.etatotal_rocks
Erebus.finalize_plastic_iteration_pass!
Erebus.finalize_thermochemical_iteration_pass
Erebus.fix
Erebus.fix_distances
Erebus.fix_weights
Erebus.get_viscosities_stresses_density_gradients!
Erebus.grid_average
Erebus.grid_vector
Erebus.initialize_pardiso!
Erebus.interpolate_add_to_grid!
Erebus.interpolate_add_to_marker!
Erebus.interpolate_to_marker!
Erebus.kphi
Erebus.ktotal
Erebus.marker_to_basic_nodes!
Erebus.marker_to_p_nodes!
Erebus.marker_to_vx_nodes!
Erebus.marker_to_vy_nodes!
Erebus.molarfraction_marker_to_p_nodes!
Erebus.move_markers_rk4!
Erebus.parse_commandline
Erebus.perform_thermal_iterations!
Erebus.perform_thermochemical_reaction!
Erebus.positive_max!
Erebus.process_gravitational_solution!
Erebus.process_hydromechanical_solution!
Erebus.recompute_bulk_viscosity!
Erebus.reduce_add_3darray!
Erebus.replenish_markers!
Erebus.reset_interpolated_properties!
Erebus.reset_thermochemical_properties!
Erebus.run_simulation
Erebus.s_to_Ma
Erebus.save_state
Erebus.setup_dynamic_simulation_parameters
Erebus.setup_gravitational_lse
Erebus.setup_hydromechanical_lse
Erebus.setup_interpolated_properties
Erebus.setup_marker_geometry_helpers
Erebus.setup_marker_properties
Erebus.setup_marker_properties_helpers
Erebus.setup_staggered_grid_properties
Erebus.setup_staggered_grid_properties_helpers
Erebus.setup_thermal_lse
Erebus.simulation_loop
Erebus.symmetrize_p_node_observables!
Erebus.total
Erebus.update_marker_population_geometry!
Erebus.update_marker_porosity!
Erebus.update_marker_stress!
Erebus.update_marker_temperature!
Erebus.update_marker_viscosity!
Erebus.update_p_nodes_melt_composition!
Erebus.ηᶠcur_inv_kᵠ
Module Functions
Erebus.Q_radiogenic
— MethodCompute radiogenic heat production of isotope mixture.
Q_radiogenic(f, ratio, E, tau, time)
Details
- f: fraction of radioactive matter [atoms/kg]
- ratio: initial ratio of radioactive to non-radioactive isotopes
- E: heat energy [J]
- tau: exp decay mean lifetime ``\tau=\frac{t_{1/2}}{\log{2}}`` [s]
- time: time elapsed since start of radioactive decay [s]
Returns
- Q: radiogenic heat production [W/kg]
Erebus.add_vrk4
— MethodAdd a RK4 stage velocity to the RK4 velocity vector.
add_vrk4(vrk4, v, rk)
Details:
- vrk4: current RK4 velocity vector
- v: RK4 velocity of stage `rk` to be added to velocity vector
- rk: RK4 stage number
Returns
- vrk4: updated RK4 velocity vector
Erebus.apply_insulating_boundary_conditions!
— MethodApply insulating boundary conditions to given array:
[x x x x x x; x a b c d x; x e f g h x; x x x x x x]
becomes
[a a b c d d; a a b c d d; e e f g h h; e e f g h h]
Details
- t: array to apply insulating boundary conditions to
Returns
- nothing
Erebus.apply_subgrid_stress_diffusion!
— MethodApply xy (basic grid) and xx (P grid) subgrid stress diffusion to markers.
apply_subgrid_stress_diffusion!(
xm,
ym,
tm,
inv_gggtotalm,
sxxm,
sxym,
SXX0,
SXY0,
DSXX,
DSXY,
SXXSUM,
SXYSUM,
WTPSUM,
WTSUM,
dt,
marknum
)
Details
- xm: marker x-coordinates
- ym: marker y-coordinates
- tm: marker type
- inv_gggtotalm: inverse of total shear modulus of markers
- sxxm: marker σ′xx [Pa]
- sxym: marker σxy [Pa]
- SXX0: σ₀′xx at P nodes [1/s]
- SXY0: σ₀xy at basic nodes [1/s]
- DSXX: stress change Δσ′xx at P nodes
- DSXY: stress change Δσxy at basic nodes
- SXXSUM: interpolation of SXX at P nodes
- SXYSUM: interpolation of SXY at basic nodes
- WTPSUM: interpolation weights at P nodes
- WTSUM: interpolation weights at basic nodes
- dt: computational time step
- marknum: total number of markers in use
Returns
- nothing
Erebus.apply_subgrid_temperature_diffusion!
— MethodApply subgrid temperature diffusion to markers.
apply_subgrid_temperature_diffusion!(
xm,
ym,
tm,
tkm,
phim,
tk1,
DT,
TKSUM,
RHOCPSUM,
dt,
marknum,
mode
)
Details
- xm: marker x-coordinates
- ym: marker y-coordinates
- tm: marker type
- tkm: marker temperature
- phim: marker porosity
- tk1: current temperature at P nodes
- DT: temperature change at P nodes
- TKSUM: interpolation of TK at P nodes
- RHOCPSUM: interpolation of RHOCP at P nodes
- dt: computational time step
- marknum: total number of markers in use
- mode: marker property computation mode
- 1: dynamic, based on (Touloukian, 1970; Hobbs, 1974;
Travis and Schubert, 2005)
- 2: constant parameter rhocpfluidm
Returns
- nothing
Erebus.assemble_gravitational_lse!
— MethodAssemble the LHS sparse coefficient matrix and fill RHS coefficient vector of the Poisson equation to be solved for the gravitational potential Φ.
assemble_gravitational_lse!(RHO, RP)
Details
- RHO: density at P nodes
- RP: right hand side coefficient vector
Returns
- LP: LHS sparse coefficient matrix
Erebus.assemble_hydromechanical_lse!
— MethodAssemble hydromechanical system of equations.
assemble_hydromechanical_lse!(
ETA,
ETAP,
GGG,
GGGP,
SXY0,
SXX0,
RHOX,
RHOY,
RHOFX,
RHOFY,
RX,
RY,
ETAPHI,
BETAPHI,
PHI,
gx,
gy,
pr0,
pf0,
DMP,
dt,
R
)
Details
- ETA: viscosity at basic nodes
- ETAP: viscosity at P nodes
- GGG: shear modulus at basic nodes
- GGGP: shear modulus at P nodes
- SXY0: previous XY stress at basic nodes
- SXX0: previous XX stress at P nodes
- RHOX: total density at Vx nodes
- RHOY: total density at Vy nodes
- RHOFX: fluid density at Vx nodes
- RHOFY: fluid density at Vy nodes
- RX: ηfluid/Kϕ at Vx nodes
- RY: ηfluid/Kϕ at Vy nodes
- ETAPHI: bulk viscosity at P nodes
- BETAPHI: bulk compressibility at P nodes
- PHI: porosity at P nodes
- gx: x gravitational acceleration at Vx nodes
- gy: y gravitational acceleration at Vy nodes
- pr0: previous total pressure at P nodes
- pf0: previous fluid pressure at P nodes
- DMP: mass transfer term at P nodes
- dt: time step
- R: vector to store RHS coefficients
Returns
- L: LHS coefficient matrix
Erebus.assemble_thermal_lse!
— MethodAssemble the LHS sparse coefficient matrix and fill RHS coefficient vector of the energy conservation (heat) equation.
assemble_thermal_lse!(
tk1,
RHOCP,
KX,
KY,
HR,
HA,
HS,
DHP,
RT,
dt
)
Details
- tk1: current temperature at P nodes
- RHOCP: volumetric heat capacity at P nodes
- KX: thermal conductivity at Vx nodes
- KY: thermal conductivity at Vy nodes
- HR: radioactive heating at P nodes
- HA: adiabatic heating at P nodes
- HS: shear heating at P nodes
- DHP: latent heating (HL) at P nodes
- RT: thermal RHS coefficient vector
- dt: current time step length
Returns
- LT: LHS sparse coefficient matrix
Erebus.backtrace_pressures_rk4!
— MethodBacktrack pressure nodes using classic Runge-Kutta integration (RK4) to update total, solid, and fluid pressure under consideration of two-phase flow velocities.
backtrace_pressures_rk4!(
pr,
pr0,
ps,
ps0,
pf,
pf0,
vx,
vy,
vxf,
vyf,
dt
)
Details
- pr: total pressure at P nodes
- pr0: previous time step total pressure at P nodes
- ps: solid pressure at P nodes
- ps0: previous time step solid pressure at P nodes
- pf: fluid pressure at P nodes
- pf0: previous time step fluid pressure at P nodes
- vx: solid vx-velocity at Vx nodes
- vy: solid vy-velocity at Vy nodes
- vxf: fluid vx-velocity at Vx nodes
- vyf: fluid vy-velocity at Vy nodes
- dt: computational time step
Returns
- nothing
Erebus.calculate_radioactive_heating
— MethodCompute radiogenic heat production of 26Al and 60Fe isotopes.
calculate_radioactive_heating(al, fe, timesum)
Details
- al: true if radioactive isotope 26Al is present
- fe: true if radioactive isotope 60Fe is present
- timesum: time elapsed since initial conditions at start of simulation
Returns
- hrsolidm: radiogenic heat production of 26Al [W/m^3]
- hrfluidm: radiogenic heat production of 60Fe [W/m^3]
Erebus.compute_Aϕ!
— MethodCompute porosity coefficient Aϕ = Dln[(1-ϕ)/ϕ]/Dt
compute_Aϕ!(
APHI,
ETAPHI,
BETAPHI,
PHI,
pr,
pf,
pr0,
pf0,
dt
)
Details
In
- ETAPHI: bulk viscosity at P Nodes
- BETAPHI: bulk compressibility at P nodes
- PHI: porosity at P Nodes
- pr: total pressure at P nodes
- pf: fluid pressure at P nodes
- pr0: previous step total pressure at P nodes
- pf0: previous step fluid pressure at P nodes
- dt: time step
Out
- APHI: porosity coefficient at P nodes
Returns
- aphimax: maximum absolute porosity coefficient
Erebus.compute_adiabatic_heating!
— MethodCompute adiabatic heating based on basic (temperature) and P grids.
compute_adiabatic_heating!(
HA,
tk1,
ALPHA,
ALPHAF,
PHI,
vx,
vy,
vxf,
vyf,
ps,
pf
)
Details
- HA: adiabatic heating at P nodes
- tk1: previous temperature at P nodes
- ALPHA: thermal expansion coefficient at P nodes
- ALPHAF: fluid thermal expansion coefficient at P nodes
- PHI: porosity at P nodes
- vx: solid vx-velocity at Vx nodes
- vy: solid vy-velocity at Vy nodes
- vxf: fluid vx-velocity at Vx nodes
- vyf: fluid vy-velocity at Vy nodes
- ps: solid pressure at P nodes
- pf: fluid pressure at P nodes
Returns
- nothing
Erebus.compute_basic_node_properties!
— MethodCompute properties of basic nodes based on interpolation arrays.
compute_basic_node_properties!(
ETA0SUM,
ETASUM,
GGGSUM,
SXYSUM,
COHSUM,
TENSUM,
FRISUM,
WTSUM,
ETA0,
ETA,
GGG,
SXY0,
COH,
TEN,
FRI,
YNY
)
Details
- ETA0SUM: ETA0 interpolation array
- ETASUM: ETA interpolation array
- GGGSUM: GGG interpolation array
- SXYSUM: SXY interpolation array
- COHSUM: COH interpolation array
- TENSUM: TEN interpolation array
- FRISUM: FRI interpolation array
- WTSUM: WT interpolation array
- ETA0: ETA0 basic node array
- ETA: ETA basic node array
- GGG: GGG basic node array
- SXY0: SXY basic node array
- COH: COH basic node array
- TEN: TEN basic node array
- FRI: FRI basic node array
- YNY: YNY basic node array
Returns
- nothing
Erebus.compute_displacement_timestep
— MethodCompute velocity-/displacement-limited time step.
compute_displacement_timestep(vx, vy, vxf, vyf, dt, aphimax)
Details
- vx: solid vx velocity at Vx nodes
- vy: solid vy velocity at Vy nodes
- vxf: fluid vx velocity at Vx nodes
- vyf: fluid vy velocity at Vy nodes
- dt: current time step
- aphimax: maximum observed porosity coefficient
Returns
- dt: displacement time step
Erebus.compute_fluid_velocities!
— MethodCompute current fluid velocities.
compute_fluid_velocities!(
PHIX,
PHIY,
qxD,
qyD,
vx,
vy,
vxf,
vyf
)
Details
In
- PHIX: porosity at Vx nodes
- PHIY: porosity at Vy nodes
- qxD: qx-Darcy flux at Vx nodes
- qyD: qy-Darcy flux at Vy nodes
- vx: solid velocity at Vx nodes
- vy: solid velocity at Vy nodes
Out
- vxf: fluid vx velocity at Vx nodes
- vyf: fluid vy velocity at Vy nodes
Returns
- nothing
Erebus.compute_gibbs_free_energy
— MethodCompute molar Gibbs free energy for single dehydration reaction Wsilicate = Dsilicate + H₂O (16.144)
compute_gibbs_free_energy(T, pf, XDˢ, XWˢ, Δt, Δtr)
Details
- T: temperature
- pf: fluid pressure
- XDˢ: molar fraction of dry solid
- XWˢ: molar fraction of wet solid
- Δt: timestep size
- Δtr: total reaction time Δtreaction
Returns
- ΔGWD: molar Gibbs free energy for single dehydration reaction (16.165a/b).
Erebus.compute_gravity_solution!
— MethodCompute gravity solution in P nodes to obtain gravitational accelerations gx for Vx nodes, gy for Vy nodes.
compute_gravity_solution!(SP, RP, RHO, FI, gx, gy)
Details
- SP: solution vector
- RP: right hand side vector
- RHO: density at P nodes
- FI: gravity potential at P nodes
- gx: x gravitational acceleration at Vx nodes
- gy: y gravitational acceleration at Vy nodes
Returns
- nothing
Erebus.compute_kfluidm
— MethodCompute thermal conductivity of H₂O (fluid phase) based on temperature.
compute_kfluidm(T, mode)
Details
- T: temperature [K]
- mode:
- 1: dynamic, based on (Touloukian, 1970; Hobbs, 1974;
Grimm & Mcsween, 1989; Bland & Travis, 2017)
- 9: constant parameter ksolidm
Returns
- kᶠ: thermal conductivity of fluid
Erebus.compute_ksolidm
— MethodCompute thermal conductivity of silicate (solid phase) based on temperature.
compute_ksolidm(T, mode)
Details
- T: temperature [K]
- mode:
- 1: dynamic, based on (Gerya, 2019)
- 9: constant parameter ksolidm
Returns
- kᶠ: thermal conductivity of solid
Erebus.compute_marker_properties!
— MethodCompute properties of given marker and save them to corresponding arrays.
compute_marker_properties!(
m,
tm,
tkm,
rhototalm,
rhocptotalm,
etatotalm,
hrtotalm,
ktotalm,
tkm_rhocptotalm,
etafluidcur_inv_kphim,
hrsolidm,
hrfluidm,
phim,
XWˢm₀,
mode
)
Details
- m: marker number
- tm: type of markers
- tkm: temperature of markers
- rhototalm: total density of markers
- rhocptotalm: total volumetric heat capacity of markers
- etatotalm: total viscosity of markers
- hrtotalm: total radiogenic heat production of markers
- ktotalm: total thermal conductivity of markers
- tkm_rhocptotalm: total thermal energy of markers
- etafluidcur_inv_kphim: (fluid viscosity)/permeability of markers
- hrsolidm: vector of radiogenic heat production of solid materials
- hrfluidm: vector of radiogenic heat production of fluid materials
- phim: porosity of markers
- XWˢm₀: previous wet solid molar fraction of markers
- mode: marker property computation mode
- 1: dynamic, based on (Touloukian, 1970; Hobbs, 1974;
Travis and Schubert, 2005)
- 2: constant parameter rhocpfluidm
Returns
- nothing
Erebus.compute_molarfraction!
— MethodCompute wet solid molar fraction at P nodes based on interpolation arrays.
compute_molarfraction!(XWSSUM, WTPSUM, XWS)
Details
- XWSSUM: XWX interpolation array
- WTPSUM: WTP interpolation array
- XWS: wet solid molar fraction at P nodes
Returns
- nothing
Erebus.compute_nodal_adjustment!
— MethodCompute nodal adjustment and return plastic iterations completeness status.
compute_nodal_adjustment!(
ETA,
ETA0,
ETA5,
GGG,
SXX,
SXY,
pr,
pf,
COH,
TEN,
FRI,
YNY,
YNY5,
YERRNOD,
DSY,
dt,
iplast
)
Details
- ETA: viscoplastic viscosity at basic nodes
- ETA0: previous time step viscoplastic viscosity at basic nodes
- ETA5: plastic iterations viscoplastic viscosity at basic nodes
- GGG: shear modulus at basic nodes
- SXX: σ′xx at P nodes
- SXY: σxy at basic nodes
- pr: total pressure at P nodes
- pf: fluid pressure at P nodes
- COH: compressive strength at basic nodes
- TEN: tensile strength at basic nodes
- FRI: friction at basic nodes
- YNY: plastic yielding status at basic nodes
- YNY5: plastic iterations plastic yielding status at basic nodes
- YERRNOD: vector of summed yielding errors of nodes over plastic iterations
- DSY: (SIIB-syield) at basic nodes
- dt: time set
- iplast: plastic iteration step
Returns
- plastic_iterations_complete: true if plastic iterations complete
Erebus.compute_p_node_properties!
— MethodCompute properties of P nodes based on interpolation arrays.
compute_p_node_properties!(
RHOSUM,
RHOCPSUM,
ALPHASUM,
ALPHAFSUM,
HRSUM,
GGGPSUM,
SXXSUM,
TKSUM,
PHISUM,
WTPSUM,
RHO,
RHOCP,
ALPHA,
ALPHAF,
HR,
GGGP,
SXX0,
tk1,
PHI,
BETAPHI
)
Details
- GGGPSUM: GGGP interpolation array
- SXX0SUM: SXX0 interpolation array
- RHOSUM: RHO interpolation array
- RHOCPSUM: RHOCP interpolation array
- ALPHASUM: ALPHA interpolation array
- ALPHAFSUM: ALPHAF interpolation array
- HRSUM: HR interpolation array
- PHISUM: PHI interpolation array
- TKSUM: TK interpolation array
- WTPSUM: WTP interpolation array
- GGGP: GGGP P node array
- SXX0: SXX0 P node array
- RHO: RHO P node array
- RHOCP: RHOCP P node array
- ALPHA: ALPHA P node array
- ALPHAF: ALPHAF P node array
- HR: HR P node array
- PHI: PHI P node array
- BETAPHI: BETAPHI P node array
- tk1: tk1 P node array
Returns
- nothing
Erebus.compute_reaction_constant
— MethodCompute dehydration reaction constant (16.151).
compute_reaction_constant(T, pf, ΔGWD)
Details
- T: temperature
- pf: fluid pressure
- ΔGWD: Gibbs free energy for dehydration reaction
Returns
- KWD: dehydration reaction constant (16.151)
Erebus.compute_relative_enthalpy
— MethodCompute relative enthalpy of system for single dehydration reaction Wsilicate = Dsilicate + H₂O (16.144).
compute_relative_enthalpy(Xsolid, XWsolid)
Details
- T: temperature
- pf: fluid pressure
- XDsolid: molar fraction of dry solid
- XWsolid: molar fraction of wet solid
- Δt: timestep size
Returns
- Hᵗ: relative enthalpy of system for single dehydration reaction (16.163)
Erebus.compute_rhocpfluidm
— MethodCompute volumetric isobaric heat capacity of H₂O (fluid phase) based on temperature.
compute_rhocpfluidm(T, mode)
Details
- T: temperature [K]
- mode: marker property computation mode
- 1: dynamic, based on (Touloukian, 1970; Hobbs, 1974;
Travis and Schubert, 2005)
- 2: constant parameter rhocpfluidm
Returns
- ρᶠCₚᶠ: volumetric isobaric heat capacity of fluid
Erebus.compute_rotation_rate!
— MethodCompute rotatation rate in basic nodes based on velocity derivatives at Vx and Vy nodes.
compute_rotation_rate!(vx, vy, wyx)
Details
- vx: solid vx-velocity at Vx nodes
- vy: solid vy-velocity at Vy nodes
- wyx: rotation rate at basic nodes
Returns
- nothing
Erebus.compute_shear_heating!
— MethodCompute shear heating based on basic (temperature) and P grids.
compute_shear_heating!(
HS,
ETA,
SXY,
ETAP,
SXX,
RX,
RY,
qxD,
qyD,
PHI,
ETAPHI,
pr,
pf
)
Details
- HS: shear heating
- ETA: viscoplastic viscosity at basic nodes
- SXY: σ₀xy XY stress at basic nodes
- ETAP: viscosity at P nodes
- SXX: σ₀xy XY stress at basic nodes
- RX: ηfluid/Kϕ at Vx nodes
- RY: ηfluid/Kϕ at Vy nodes
- qxD: qx-Darcy flux at Vx nodes
- qyD: qy-Darcy flux at Vy nodes
- PHI: porosity at P nodes
- ETAPHI: bulk viscosity at P nodes
- pr: total pressure at P nodes
- pf: fluid pressure at P nodes
Returns
- nothing
Erebus.compute_stress_strainrate!
— MethodCompute stress, stress change, and strain rate components.
compute_stress_strainrate!(
vx,
vy,
ETA,
GGG,
ETAP,
GGGP,
SXX0,
SXY0,
EXX,
EXY,
SXX,
SXY,
DSXX,
DSXY,
EII,
SII,
dt
)
Details
In
- vx: solid vx velocity at Vx nodes
- vy: solid vy velocity at Vy nodes
- ETA: viscosity at basic nodes
- GGG: shear modulus at basic nodes
- ETAP: viscosity at P nodes
- GGGP: shear modulus at P nodes
- SXX0: previous time step σ₀′xx at P nodes
- SXY0: previous time step σ₀xy at basic nodes
- dt: computational time step
Out
- EXX: ϵxx at P nodes
- EXY: ϵxy at basic nodes
- SXX: σ′xx P nodes
- SXY: σxy at basic nodes
- DSXX: stress change Δσ′xx at P nodes
- DSXY: stress change Δσxy at basic nodes
- EII: second strain rate invariant ϵᴵᴵ at P nodes
- SII: second stress invariant σᴵᴵ at P nodes
Returns
- nothing
Erebus.compute_thermochemical_iteration_outcome
— MethodAssess outcome of thermochemical iteration and return thermochemical iterations completeness status.
compute_thermochemical_iteration_outcome(
DMP,
pf,
pf0,
titer
)
Details:
- DMP: mass transfer term at P nodes
- pf: fluid pressure at P nodes
- pf0: previous time step fluid pressure at P nodes
- titer: current thermochemical iteration counter
Returns
- dt: adjusted next time step
Erebus.compute_thermodynamic_xfer!
— MethodCompute thermodynamic properties at P nodes based on interpolation arrays.
compute_thermodynamic_xfer!(
DMPSUM,
DHPSUM,
WTPSUM,
DMP,
DHP
)
Details
- DMPSUM: DMP interpolation array
- DHPSUM: DHP interpolation array
- WTPSUM: WTP interpolation array
- DMP: mass transfer term at P nodes
- DHP: enthalpy transfer/latent heating term at P nodes
Returns
- nothing
Erebus.compute_velocities!
— MethodCompute solid velocities, fluid velocities at P nodes.
compute_velocities!(vx, vy, vxf, vyf, vxp, vyp, vxpf, vypf)
Details
- vx: solid vx-velocity at Vx nodes
- vy: solid vy-velocity at Vy nodes
- vxf: fluid vx-velocity at Vx nodes
- vyf: fluid vy-velocity at Vy nodes
- vxp: solid vx-velocity at P nodes
- vyp: solid vy-velocity at P nodes
- vxpf: fluid vx-velocity at P nodes
- vypf: fluid vy-velocity at P nodes
Returns
- nothing
Erebus.compute_vx_node_properties!
— MethodCompute properties of Vx nodes based on interpolation arrays.
compute_vx_node_properties!(
RHOXSUM,
RHOFXSUM,
KXSUM,
PHIXSUM,
RXSUM,
WTXSUM,
RHOX,
RHOFX,
KX,
PHIX,
RX
)
Details
- RHOXSUM: RHOX interpolation array
- RHOFXSUM: RHOFX interpolation array
- KXSUM: KX interpolation array
- PHIXSUM: PHIX interpolation array
- RXSUM: RX interpolation array
- WTXSUM: WTX interpolation array
- RHOX: RHOX Vx node array
- RHOFX: RHOFX Vx node array
- KX: KX Vx node array
- PHIX: PHIX Vx node array
- RX: RX Vx node array
Returns
- nothing
Erebus.compute_vy_node_properties!
— MethodCompute properties of Vy nodes based on interpolation arrays.
compute_vy_node_properties!(
RHOYSUM,
RHOFYSUM,
KYSUM,
PHIYSUM,
RYSUM,
WTYSUM,
RHOY,
RHOFY,
KY,
PHIY,
RY
)
Details
- RHOYSUM: RHOY interpolation array
- RHOFYSUM: RHOFY interpolation array
- KYSUM: KY interpolation array
- PHIYSUM: PHIY interpolation array
- RYSUM: RY interpolation array
- WTYSUM: WTY interpolation array
- RHOY: RHOY Vy node array
- RHOFY: RHOFY Vy node array
- KY: KY Vy node array
- PHIY: PHIY Vy node array
- RY: RY Vy node array
Returns
- nothing
Erebus.compute_Δtreaction
— MethodCompute dehydration reaction time Δtreaction based on temperature and porosity according to selected method:
compute_Δtreaction(T, ϕ, mode)
Details
- T: temperature
- ϕ: porosity
- mode:
- 1: Gaussian form reaction rate coefficient, based on (Martin & Fyfe,
1970; Emmanuel & Berkowitz, 2006; Iyer et al., 2012).
- 2: pseudo-Arrhenius form reaction rate coefficient, based on
(Bland & Travis, 2017).
- 3: Arrhenius form reaction coefficent, based on (Travis et al., 2018).
- 9: constant parameter Δtreaction
Returns
- Δtreaction: dehydration reaction time
Erebus.define_markers!
— MethodDefine initial set of markers according to model parameters
define_markers!(
xm,
ym,
tm,
phim,
etavpm,
rhototalm,
rhocptotalm,
etatotalm,
hrtotalm,
ktotalm,
tkm,
inv_gggtotalm,
fricttotalm,
cohestotalm,
tenstotalm,
rhofluidcur,
alphasolidcur,
alphafluidcur,
XWsolidm0;
randomized
)
Details
- xm: x coordinates of markers
- ym: y coordinates of markers
- tm: material type of markers
- phim: porosity of markers
- etavpm: matrix viscosity of markers
- rhototalm: total density of markers
- rhocptotalm: total volumetric heat capacity of markers
- etatotalm: total viscosity of markers
- hrtotalm: total radiogenic heat production of markers
- ktotalm: total thermal conductivity of markers
- tkm: temperature of markers
- inv_gggtotalm: inverse of total shear modulus of markers
- fricttotalm: total friction coefficient of markers
- cohestotalm: total compressive strength of markers
- tenstotalm: total tensile strength of markers
- rhofluidcur: fluid density of markers
- alphasolidcur: solid thermal expansion coefficient of markers
- alphafluidcur: fluid thermal expansion coefficient of markers
- XWsolidm0: previous wet solid molar fraction of markers
- randomized: uniformly random-distribute marker x/y positions within cells
and randomly set initial marker porosity
Returns
- nothing
Erebus.distance
— MethodCalculate Euclidean distance between two point coordinates.
distance(x1, y1, x2, y2)
Details
- x1: x-coordinate of point 1 [m]
- y1: y-coordinate of point 1 [m]
- x2: x-coordinate of point 2 [m]
- y2: y-coordinate of point 2 [m]
Returns
- Euclidean distance between point 1 and point 2 [m]
Erebus.dot4
— MethodCompute inner product of two 4-vectors of reals.
dot4(v1, v2)
Details
- v1: first 4-vector
- v2: second 4-vector
Returns
- inner product of v1 and v2
Erebus.etatotal_rocks
— MethodCompute total rocky marker viscosity based on temperature and material type.
etatotal_rocks(tkmm, tmm)
Details
- tkmm: marker temperature [K]
- tmm: marker type [1, 2]
Returns
- etatotal: rocky marker temperature-dependent total viscosity
Erebus.finalize_plastic_iteration_pass!
— MethodDecide next pass plastic iteration time step, viscoplastic viscosity, and basic node yielding status.
finalize_plastic_iteration_pass!(
ETA,
ETA5,
ETA00,
YNY,
YNY5,
YNY00,
YNY_inv_ETA,
dt,
iplast
)
Details:
- ETA: viscoplastic viscosity at basic nodes
- ETA5: plastic iterations viscoplastic viscosity at basic nodes
- ETA00: previous time step viscoplastic viscosity at basic nodes
- YNY: plastic yielding status at basic nodes
- YNY5: plastic iterations plastic yielding status at basic nodes
- YNY00: previous time step plastic yielding status at basic nodes
- YNY_inv_ETA: inverse of plastic viscosity at yielding basic nodes
- dt: current time step
- iplast: current plastic iteration counter
Returns
- dt: adjusted next time step
Erebus.finalize_thermochemical_iteration_pass
— MethodDecide next pass thermochemical iteration time step.
finalize_thermochemical_iteration_pass(
maxDTcurrent,
dt,
titer
)
Details:
- maxDTcurrent: maximum temperature difference between current and
previous time step
- dt: current time step duration
- titer: current thermochemical iteration counter
Returns
- dt: adjusted next time step
Erebus.fix
— MethodCompute top and left grid nodes indices (i, j) and x- and y-distances to that grid node (i, j) for given (x, y) position and grid axes.
fix(x, y, x_axis, y_axis, dx, dy, jmin, jmax, imin, imax)
Details
- x: x-position [m]
- y: y-position [m]
- x_axis: x-grid reference axis array [m]
- y_axis: y-grid reference axis array [m]
- dx: x-grid axis mesh width [m]
- dy: y-grid axis mesh width [m]
- jmin: minimum assignable index on x-grid axis (basic/Vx/Vy/P)
- jmax: maximum assignable index on x-grid axis (basic/Vx/Vy/P)
- imin: minimum assignable index on y-grid axis (basic/Vx/Vy/P)
- imax: maximum assignable index on y-grid axis (basic/Vx/Vy/P)
Returns
- i: top (with reference to y) node index on y-grid axis
- j: left (with reference to x) node index on x-grid axis
Erebus.fix_distances
— MethodCompute top and left grid nodes indices (i, j) and x- and y-distances to that grid node (i, j) for given (x, y) position and grid axes.
fix_distances(
x,
y,
x_axis,
y_axis,
dx,
dy,
jmin,
jmax,
imin,
imax
)
Details
- x: x-position [m]
- y: y-position [m]
- x_axis: x-grid reference axis array [m]
- y_axis: y-grid reference axis array [m]
- dx: x-grid axis mesh width [m]
- dy: y-grid axis mesh width [m]
- jmin: minimum assignable index on x-grid axis (basic/Vx/Vy/P)
- jmax: maximum assignable index on x-grid axis (basic/Vx/Vy/P)
- imin: minimum assignable index on y-grid axis (basic/Vx/Vy/P)
- imax: maximum assignable index on y-grid axis (basic/Vx/Vy/P)
Returns
- i: top (with reference to y) node index on y-grid axis
- j: left (with reference to x) node index on x-grid axis
- dxmj: x-distance from (x, y) point to (i, j) node
- dymi: y-distance from (x, y) point to (i, j) node
Erebus.fix_weights
— MethodCompute top and left grid nodes indices and bilinear interpolation weigths to nearest four grid nodes for given (x, y) position and grid axes.
fix_weights(
x,
y,
x_axis,
y_axis,
dx,
dy,
jmin,
jmax,
imin,
imax
)
Details
- x: x-position [m]
- y: y-position [m]
- x_axis: x-grid reference axis array [m]
- y_axis: y-grid reference axis array [m]
- dx: x-grid axis mesh width [m]
- dy: y-grid axis mesh width [m]
- jmin: minimum assignable index on x-grid axis (basic/Vx/Vy/P)
- jmax: maximum assignable index on x-grid axis (basic/Vx/Vy/P)
- imin: minimum assignable index on y-grid axis (basic/Vx/Vy/P)
- imax: maximum assignable index on y-grid axis (basic/Vx/Vy/P)
Returns
- i: top (with reference to y) node index on y-grid axis
- j: left (with reference to x) node index on x-grid axis
- bilinear_weights: vector of 4 bilinear interpolation weights to
nearest four grid nodes:
[wtmij : i , j node,
wtmi1j : i+1, j node,
wtmij1 : i , j+1 node,
wtmi1j1: i+1, j+1 node]
Erebus.get_viscosities_stresses_density_gradients!
— MethodCompute viscosities, stresses, and density gradients for hydromechanical solver.
get_viscosities_stresses_density_gradients!(
ETA,
ETAP,
GGG,
GGGP,
SXY0,
SXX0,
RHOX,
RHOY,
dt,
ETAcomp,
ETAPcomp,
SXYcomp,
SXXcomp,
SYYcomp,
dRHOXdx,
dRHOXdy,
dRHOYdx,
dRHOYdy
)
Details
In
- ETA: viscoplastic viscosity at basic nodes
- ETAP: viscosity at P nodes
- GGG: shear modulus at basic nodes
- GGGP: shear modulus at P nodes
- SXY0: σ₀xy XY stress at basic nodes
- SXX0:σ₀xy XY stress at basic nodes
- RHOX: density at Vx nodes
- RHOY: density at Vy nodes
- dt: time step
Out
- ETAcomp: computational viscosity at basic nodes
- ETAPcomp: computational viscosity at P nodes
- SXYcomp: previous XY stresses at basic nodes
- SXXcomp: previous XX stresses at P nodes
- SYYcomp: previous YY stresses at P nodes
- dRHOXdx: density gradient at Vx nodes in x direction
- dRHOXdy: density gradient at Vx nodes in y direction
- dRHOYdx: density gradient at Vy nodes in x direction
- dRHOYdy: density gradient at Vy nodes in y direction
Returns
- nothing
Erebus.grid_average
— MethodGet the arithmetic average of values from a grid 4-stencil anchored at (i, j).
grid_average(i, j, grid)
Details
- i: top left grid node column index
- j: top left grid node row index
- grid: data from which to get the average value
Returns
- grid_average: (grid[i, j]+grid[i+1, j]+grid[i, j+1]+grid[i+1, j+1]) / 4
Erebus.grid_vector
— MethodGet a 4-vector of values from a grid 4-stencil anchored at (i, j) in column-major order.
grid_vector(i, j, grid)
Details
- i: top left grid node column index
- j: top left grid node row index
- grid: data from which to build Vector
Returns
- grid_vector: 4-vector of values
[grid[i, j], grid[i+1, j], grid[i, j+1], grid[i+1, j+1]]
Erebus.initialize_pardiso!
— MethodInitialize iparm
parameters of Pardiso MKL solver.
initialize_pardiso!(pardiso_solver, iparms_dict)
Details
- ps: Instance of pardiso solver
- iparms_dict: dictionary of iparm parameters
Returns
- nothing
Erebus.interpolate_add_to_grid!
— MethodInterpolate a property to neareast four nodes on a given grid location using given bilinear interpolation weights.
Details
- i: top (with reference to y) node index on vertical y-grid axis
- j: left (with reference to x) node index on horizontal x-grid axis
- weights: vector of 4 bilinear interpolation weights to
nearest four grid nodes:
[wtmij : i , j node,
wtmi1j : i+1, j node,
wtmij1 : i , j+1 node,
wtmi1j1: i+1, j+1 node]
- property: property to be interpolated to grid using weights
- grid: threaded grid array on which to interpolate property
Returns
- nothing
Erebus.interpolate_add_to_marker!
— MethodInterpolate a property from nearest the four nodes on a given grid to a marker and add it to the markers property.
Details
- m: number of marker to interpolate to
- i: top (with reference to y) node index on vertical y-grid axis
- j: left (with reference to x) node index on horizontal x-grid axis
- weights: vector of 4 bilinear interpolation weights to
nearest four grid nodes:
[wtmij : i , j node,
wtmi1j : i+1, j node,
wtmij1 : i , j+1 node,
wtmi1j1: i+1, j+1 node]
- marker_property: marker property array into which to interpolate and add
- property_grid: grid whose property is to be interpolated to marker
- m: marker
Erebus.interpolate_to_marker!
— MethodInterpolate a property from nearest the four nodes on a given grid to a marker.
Details
- m: number of marker to interpolate to
- i: top (with reference to y) node index on vertical y-grid axis
- j: left (with reference to x) node index on horizontal x-grid axis
- weights: vector of 4 bilinear interpolation weights to
nearest four grid nodes:
[wtmij : i , j node,
wtmi1j : i+1, j node,
wtmij1 : i , j+1 node,
wtmi1j1: i+1, j+1 node]
- marker_property: marker property array into which to interpolate
- property_grid: grid whose property is to be interpolated to marker
- m: marker
Erebus.kphi
— MethodCompute porosity-dependent permeability (eqn 16.64 in Gerya (2019)).
kphi(kphim0m, phimm)
Details
- kphim0m: standard (reference) permeability (of marker type) [m^2]
- phimm: actual (marker) porosity
Returns
- kphim: empirical porosity-dependent permeability [m^2]
Erebus.ktotal
— MethodCompute total thermal conductivity of two-phase material.
ktotal(ksolid, kfluid, phi)
Details
- ksolid: solid thermal conductivity [W/m/K]
- kfluid: fluid thermal conductivity [W/m/K]
- phi: fraction of solid
Returns
- ktotal: total thermal conductivity of mixed phase [W/m/K]
Erebus.marker_to_basic_nodes!
— MethodInterpolate selected marker properties to basic nodes.
marker_to_basic_nodes!(
m,
xmm,
ymm,
etatotalm,
etavpm,
inv_gggtotalm,
sxym,
cohestotalm,
tenstotalm,
fricttotalm,
ETA0SUM,
ETASUM,
GGGSUM,
SXYSUM,
COHSUM,
TENSUM,
FRISUM,
WTSUM
)
Details
- m: marker number
- xmm: marker's x-position [m]
- ymm: marker's y-position [m]
- etatotalm: total viscosity of markers
- etavpm: matrix viscosity of markers
- inv_gggtotalm: inverse of total shear modulus of markers
- sxym: marker σxy [Pa]
- cohestotalm: total compressive strength of markers
- tenstotalm: total tensile strength of markers
- fricttotalm: total friction coefficient of markers
- ETA0SUM: viscous viscosity interpolated to basic nodes
- ETASUM: viscoplastic viscosity interpolated to basic nodes
- GGGSUM: shear modulus interpolated to basic nodes
- SXYSUM: σxy shear stress interpolated to basic nodes
- COHSUM: compressive strength interpolated to basic nodes
- TENSUM: tensile strength interpolated to basic nodes
- FRISUM: friction interpolated to basic nodes
- WTSUM: weight array for bilinear interpolation to basic nodes
Returns
- nothing
Erebus.marker_to_p_nodes!
— MethodInterpolate selected marker properties to P nodes.
marker_to_p_nodes!(
m,
xmm,
ymm,
inv_gggtotalm,
sxxm,
rhototalm,
rhocptotalm,
alphasolidcur,
alphafluidcur,
hrtotalm,
phim,
tkm_rhocptotalm,
GGGPSUM,
SXXSUM,
RHOSUM,
RHOCPSUM,
ALPHASUM,
ALPHAFSUM,
HRSUM,
PHISUM,
TKSUM,
WTPSUM
)
Details
- m: marker number
- xmm: marker x-position [m]
- ymm: marker y-position [m]
- inv_gggtotalm: inverse of total shear modulus of markers
- sxxm: marker σ′xx [Pa]
- rhototalm: total density of markers
- rhocptotalm: total volumetric heat capacity of markers
- alphasolidcur: solid thermal expansion coefficient of markers
- alphafluidcur: fluid thermal expansion coefficient of markers
- hrtotalm: total radiogenic heat production of markers
- phim: marker porosity
- tkm_rhocptotalm: total thermal energy of markers
- GGGPSUM: shear modulus interpolated to P nodes
- SXXSUM: σ'xx interpolated to P nodes
- RHOSUM: density interpolated to P nodes
- RHOCPSUM: volumetric heat capacity interpolated to P nodes
- ALPHASUM: thermal expansion interpolated to P nodes
- ALPHAFSUM: fluid thermal expansion interpolated to P nodes
- HRSUM: radioactive heating interpolated to P nodes
- PHISUM: porosity interpolated to P nodes
- TKSUM: heat capacity interpolated to P nodes
- WTPSUM: weight for bilinear interpolation to P nodes
Returns
- nothing
Erebus.marker_to_vx_nodes!
— MethodInterpolate selected marker properties to Vx nodes.
marker_to_vx_nodes!(
m,
xmm,
ymm,
rhototalm,
rhofluidcur,
ktotalm,
phim,
etafluidcur_inv_kphim,
RHOXSUM,
RHOFXSUM,
KXSUM,
PHIXSUM,
RXSUM,
WTXSUM
)
Details
- m: marker number
- xmm: marker's x-position [m]
- ymm: marker's y-position [m]
- rhototalm: total density of markers
- rhofluidcur: fluid density of markers
- ktotalm: total thermal conductivity of markers
- phim: porosity of markers
- etafluidcur_inv_kphim: fluid viscosity over permeability of markers
- RHOXSUM: density interpolated to Vx nodes
- RHOFXSUM: fluid density interpolated to Vx nodes
- KXSUM: thermal conductivity interpolated to Vx nodes
- PHIXSUM: porosity interpolated to Vx nodes
- RXSUM: ηfluid/kϕ interpolated to Vx nodes
- WTXSUM: weight for bilinear interpolation to Vx nodes
Returns
- nothing
Erebus.marker_to_vy_nodes!
— MethodInterpolate selected marker properties to Vy nodes.
marker_to_vy_nodes!(
m,
xmm,
ymm,
rhototalm,
rhofluidcur,
ktotalm,
phim,
etafluidcur_inv_kphim,
RHOYSUM,
RHOFYSUM,
KYSUM,
PHIYSUM,
RYSUM,
WTYSUM
)
Details
- m: marker number
- xmm: marker's x-position [m]
- ymm: marker's y-position [m]
- rhototalm: total density of markers
- rhofluidcur: fluid density of markers
- ktotalm: total thermal conductivity of markers
- phim: porosity of markers
- etafluidcur_inv_kphim: fluid viscosity over permeability of markers
- RHOYSUM: density interpolated to Vy nodes
- RHOFYSUM: fluid density interpolated to Vy nodes
- KYSUM: thermal conductivity interpolated to Vy nodes
- PHIYSUM: porosity interpolated to Vy nodes
- RYSUM: ηfluid/kϕ interpolated to Vy nodes
- WTYSUM: weight for bilinear interpolation to Vy nodes
Returns
- nothing
Erebus.molarfraction_marker_to_p_nodes!
— MethodInterpolate marker wet silicate (solid) fraction to P nodes.
molarfraction_marker_to_p_nodes!(
m,
xmm,
ymm,
XWsolidm0,
XWSSUM,
WTPSUM
)
Details
- m: marker number
- xmm: marker x-position [m]
- ymm: marker y-position [m]
- XWsolidm0: previous marker wet silicate (solid) molar fraction
- XWSSUM: wet silicate (solid) molar fraction interpolated to P nodes
- WTPSUM: weight for bilinear interpolation to P nodes
Returns
- nothing
Erebus.move_markers_rk4!
— MethodMove markers using classic Runge-Kutta integration (RK4) taking into account two-phase flow and solid-fluid temperatures.
move_markers_rk4!(
xm,
ym,
tm,
tkm,
phim,
sxym,
sxxm,
vx,
vy,
vxf,
vyf,
wyx,
tk2,
marknum,
dt,
mode
)
Details
- xm: x-coordinate of markers
- ym: y-coordinate of markers
- tm: type of markers
- tkm: temperature of markers
- phim: porosity of markers
- sxxm: marker σ′xx
- sxym: marker σxy
- vx: solid vx-velocity at Vx nodes
- vy: solid vy-velocity at Vy nodes
- vxf: fluid vx-velocity at Vx nodes
- vyf: fluid vy-velocity at Vy nodes
- wyx: rotation rate at basic nodes
- tk2: next temperature at P nodes
- marknum: number of markers in use
- dt: computational time step
- mode: marker property computation mode
- 1: dynamic, based on (Touloukian, 1970; Hobbs, 1974;
Travis and Schubert, 2005)
- 2: constant parameter rhocpfluidm
Returns
- nothing
Erebus.parse_commandline
— MethodParse command line arguments and feed them to the main function.
parse_commandline()
Details:
- nothing
Returns
- parsed_args: parsed command line arguments
Erebus.perform_thermal_iterations!
— MethodPerform thermal iterations to time step thermal field at P nodes.
perform_thermal_iterations!(
tk0,
tk1,
tk2,
DT,
DT0,
RHOCP,
KX,
KY,
HR,
HA,
HS,
DHP,
RT,
ST,
dt
)
Details
- tk0: previous temperature at P nodes
- tk1: current temperature at P nodes
- tk2: next temperature at P nodes
- DT: calculated temperature difference at P nodes
- DT0: previous calculated temperature difference at P nodes
- RHOCP: volumetric heat capacity at P nodes
- KX: thermal conductivity at Vx nodes
- KY: thermal conductivity at Vy nodes
- HR: radioactive heating at P nodes
- HA: adiabatic heating at P nodes
- HS: shear heating at P nodes
- DHP: latent heating (HL) at P nodes
- RT: thermal RHS coefficient vector
- ST: thermal solution vector
- dt: computational time step
Returns
- nothing
Erebus.perform_thermochemical_reaction!
— MethodPerform hydrothermomechanical iterations to time step thermal field at P nodes.
perform_thermochemical_reaction!(
DMP,
DHP,
DMPSUM,
DHPSUM,
WTPSUM,
pf,
tk2,
tm,
xm,
ym,
XWˢm₀,
XWˢm,
phim,
phinewm,
pfm₀,
marknum,
Δt,
timestep,
titer
)
Details
- DMP: mass transfer term at P nodes
- DHP: enthalpy transfer/latent heating term at P nodes
- DMPSUM: interpolation of DMP (mass transfer term) at P nodes
- DHPSUM: interpolation of DHP (enthalpy transfer term) at P nodes
- WTPSUM: interpolation weights at P nodes
- pf: fluid pressure at P nodes
- tk2: next temperature at P nodes
- tm: type of markers
- xm: x-coordinate of markers
- ym: y-coordinate of markers
- XWˢm₀: previous marker wet silicate (solid) fraction
- XWˢm: current marker wet silicate (solid) fraction
- phim: current marker porosity
- phinewm: next generation marker porosity
- pfm₀: previous marker fluid pressure
- marknum: current total number of markers
- Δt: current time step length
- timestep: current time step
- titer: current thermochemical iteration number
Returns
- nothing
Erebus.positive_max!
— MethodCompare two arrays of identical sizes element-wise and fill a third array with the larger value if positive and zero otherwise.
positive_max!(A, B, C)
Details
- A: first array
- B: second array
- C: result array
Returns
- nothing
Erebus.process_gravitational_solution!
— MethodProcess gravitational potential solution vector to output physical observables.
process_gravitational_solution!(SP, FI, gx, gy)
Details
- FI: gravitational potential
- gx: x-component of gravitational acceleration
- gy: y-component of gravitational acceleration
Returns
- nothing
Erebus.process_hydromechanical_solution!
— MethodProcess hydromechanical solution vector to output physical observables.
process_hydromechanical_solution!(
S,
vx,
vy,
pr,
qxD,
qyD,
pf
)
Details
- S: hydromechanical solution vector
- vx: solid velocity at Vx nodes
- vy: solid velocity at Vy nodes
- pr: total pressure at P nodes
- qxD: qx-Darcy flux at Vx nodes
- qyD: qy-Darcy flux at Vy nodes
- pf: fluid pressure at P nodes
Returns
- nothing
Erebus.recompute_bulk_viscosity!
— MethodRecompute bulk viscosity at P nodes.
Details
- ETA: viscoplastic viscosity at basic nodes
- ETAP: viscosity at P nodes
- ETAPHI: bulk viscosity at P nodes
- PHI: porosity at P nodes
- etaphikoef: coefficient: shear viscosity -> compaction viscosity
Returns
- nothing
Erebus.reduce_add_3darray!
— MethodReduce a 3D (i, j, k) array along its third (k) axis by addition and write the result into (i, j, 1) without reallocating the array's memory.
reduce_add_3darray!(A)
Details
- A: 3D array [i, j, k]
Returns
- nothing
Erebus.replenish_markers!
— MethodAdd markers to populate currently sparsely filled grid areas.
replenish_markers!(
xm,
ym,
tm,
tkm,
phim,
sxxm,
sxym,
etavpm,
phinewm,
pfm0,
XWsolidm,
XWsolidm0,
rhototalm,
rhocptotalm,
etatotalm,
hrtotalm,
ktotalm,
inv_gggtotalm,
fricttotalm,
cohestotalm,
tenstotalm,
rhofluidcur,
alphasolidcur,
alphafluidcur,
tkm_rhocptotalm,
etafluidcur_inv_kphim,
mdis,
mnum;
randomized
)
Details
- xm: horizontal x-position of markers
- ym: vertical y-position of markers
- tm: type of markers
- tkm: temperature of markers
- phim: porosity of markers
- sxxm: marker σ′xx of markers
- sxym: σxy of markers
- etavpm: viscoplastic viscosity of markers
- phinewm: reacted porosity of markers
- pfm0: previous fluid pressure of markers
- XWsolidm: melt molar fraction of markers
- XWsolidm0: previous melt molar fraction of markers
- rhototalm: total density of markers
- rhocptotalm : total volumetric heat capacity of markers
- etatotalm: total viscosity of markers
- hrtotalm: total radiogenic heat production of markers
- ktotalm: total thermal conductivity of markers
- inv_gggtotalm: inverse of total shear modulus of markers
- fricttotalm: total friction coefficient of markers
- cohestotalm: total compressive strength of markers
- tenstotalm: total tensile strength of markers
- rhofluidcur: fluid density of markers
- alphasolidcur: solid thermal expansion coefficient of markers
- alphafluidcur: fluid thermal expansion coefficient of markers
- tkm_rhocptotalm: total thermal energy of markers
- etafluidcur_inv_kphim: fluid viscosity over permeability of markers
- mdis: minimum distance of marker launch anchor points to nearest marker
- mnum: number of marker nearest to marker launch anchor positions
Returns
- marknum: updated number of markers in use
Erebus.reset_interpolated_properties!
— MethodReset properties to be interpolated from markers to staggered grid.
reset_interpolated_properties!(
ETA0SUM,
ETASUM,
GGGSUM,
SXYSUM,
COHSUM,
TENSUM,
FRISUM,
WTSUM,
RHOXSUM,
RHOFXSUM,
KXSUM,
PHIXSUM,
RXSUM,
WTXSUM,
RHOYSUM,
RHOFYSUM,
KYSUM,
PHIYSUM,
RYSUM,
WTYSUM,
RHOSUM,
RHOCPSUM,
ALPHASUM,
ALPHAFSUM,
HRSUM,
GGGPSUM,
SXXSUM,
TKSUM,
PHISUM,
WTPSUM
)
Details
- ETA0SUM: interpolation of ETA0 at basic nodes
- ETASUM: interpolation of ETA at basic nodes
- GGGSUM: interpolation of GGG at basic nodes
- SXYSUM: interpolation of SXY at basic nodes
- COHSUM: interpolation of COH at basic nodes
- TENSUM: interpolation of TEN at basic nodes
- FRISUM: interpolation of FRI at basic nodes
- WTSUM: interpolation weights at basic nodes
- RHOXSUM: interpolation of RHOX at Vx nodes
- RHOFXSUM: interpolation of RHOFX at Vx nodes
- KXSUM: interpolation of KX at Vx nodes
- PHIXSUM: interpolation of PHIX at Vx nodes
- RXSUM: interpolation of RX at Vx nodes
- WTXSUM: interpolation weights at Vx nodes
- RHOYSUM: interpolation of RHOY at Vy nodes
- RHOFYSUM: interpolation of RHOFY at Vy nodes
- KYSUM: interpolation of KY at Vy nodes
- PHIYSUM: interpolation of PHIX at Vy nodes
- RYSUM: interpolation of RY at Vy nodes
- WTYSUM: interpolation weights at Vy nodes
- RHOSUM: interpolation of RHO at P nodes
- RHOCPSUM: interpolation of RHOCP at P nodes
- ALPHASUM: interpolation of ALPHA at P nodes
- ALPHAFSUM: interpolation of ALPHAF at P nodes
- HRSUM: interpolation of HR at P nodes
- GGGPSUM: interpolation of GGGP at P nodes
- SXXSUM: interpolation of SXX at P nodes
- TKSUM: interpolation of TK at P nodes
- PHISUM: interpolation of PHI at P nodes
- WTPSUM: interpolation weights at P nodes
Returns
- nothing
Erebus.reset_thermochemical_properties!
— MethodReset thermochemical properties to be interpolated to staggered grid.
reset_thermochemical_properties!(DMPSUM, DHPSUM, WTPSUM)
Details
- DMPSUM: interpolation of DMP at P nodes
- DHPSUM: interpolation of DHP at P nodes
- WTPSUM: interpolation weights at P nodes
Returns
- nothing
Erebus.run_simulation
— MethodRuns the simulation with the given parameters.
run_simulation()
Details
- nothing
Returns
- nothing
Erebus.s_to_Ma
— MethodConvert seconds to Ma (millions of years).
s_to_Ma(s)
Details
- s: period in seconds
Returns
- Ma: period in millions of years
Erebus.save_state
— MethodSave simulation state to JLD2 output file named after current timestep.
save_state(
output_path,
timestep,
dt,
timesum,
marknum,
ETA,
ETA0,
GGG,
EXY,
SXY,
SXY0,
wyx,
COH,
TEN,
FRI,
YNY,
RHOX,
RHOFX,
KX,
PHIX,
vx,
vxf,
RX,
qxD,
gx,
RHOY,
RHOFY,
KY,
PHIY,
vy,
vyf,
RY,
qyD,
gy,
RHO,
RHOCP,
ALPHA,
ALPHAF,
HR,
HA,
HS,
ETAP,
GGGP,
EXX,
SXX,
SXX0,
tk1,
tk2,
vxp,
vyp,
vxpf,
vypf,
pr,
pf,
ps,
pr0,
pf0,
ps0,
ETAPHI,
BETAPHI,
PHI,
APHI,
FI,
ETA5,
ETA00,
YNY5,
YNY00,
YNY_inv_ETA,
DSXY,
EII,
SII,
DSXX,
DMP,
DHP,
XWS,
XWsolidm0,
xm,
ym,
tm,
tkm,
sxxm,
sxym,
etavpm,
phim,
rhototalm,
rhocptotalm,
etatotalm,
hrtotalm,
ktotalm,
tkm_rhocptotalm,
etafluidcur_inv_kphim,
inv_gggtotalm,
fricttotalm,
cohestotalm,
tenstotalm,
rhofluidcur,
alphasolidcur,
alphafluidcur
)
Details
- output_path: absolute path to output directory
- timestep: current time step number
- dt: time step
- timesum: total simulation time
- marknum: number of markers
- ETA... : simulation state variables
Returns
- nothing
Erebus.setup_dynamic_simulation_parameters
— MethodSet up and initialize dynamic simulation parameters.
setup_dynamic_simulation_parameters()
Details
- nothing
Returns
- timestep: simulation starting time step count
- dt: simulation initial computational time step [s]
- timesum: simulation starting time [s]
- marknum: initial number of markers
- hrsolidm: initial radiogenic heat production solid phase
- hrfluidm: initial radiogenic heat production fluid phase
- YERRNOD: vector of summed yielding errors of nodes over plastic iterations
Erebus.setup_gravitational_lse
— MethodSet up gravitational linear system of equations structures.
setup_gravitational_lse()
Details
- nothing
Returns
- RP: gravitational linear system of equations: RHS vector
- SP: gravitational linear system of equations: solution vector
Erebus.setup_hydromechanical_lse
— MethodSet up hydromechanical linear system of equations structures.
setup_hydromechanical_lse()
Details
- nothing
Returns
- R: hydromechanical linear system of equations: RHS vector
- S: hydromechanical linear system of equations: solution vector
Erebus.setup_interpolated_properties
— MethodSet up properties to be interpolated from markers to staggered grid.
setup_interpolated_properties()
Details
- nothing
Returns
- ETA0SUM: interpolation of ETA0 at basic nodes
- ETASUM: interpolation of ETA at basic nodes
- GGGSUM: interpolation of GGG at basic nodes
- SXYSUM: interpolation of SXY at basic nodes
- COHSUM: interpolation of COH at basic nodes
- TENSUM: interpolation of TEN at basic nodes
- FRISUM: interpolation of FRI at basic nodes
- WTSUM: interpolation weights at basic nodes
- RHOXSUM: interpolation of RHOX at Vx nodes
- RHOFXSUM: interpolation of RHOFX at Vx nodes
- KXSUM: interpolation of KX at Vx nodes
- PHIXSUM: interpolation of PHIX at Vx nodes
- RXSUM: interpolation of RX at Vx nodes
- WTXSUM: interpolation weights at Vx nodes
- RHOYSUM: interpolation of RHOY at Vy nodes
- RHOFYSUM: interpolation of RHOFY at Vy nodes
- KYSUM: interpolation of KY at Vy nodes
- PHIYSUM: interpolation of PHIX at Vy nodes
- RYSUM: interpolation of RY at Vy nodes
- WTYSUM: interpolation weights at Vy nodes
- RHOSUM: interpolation of RHO at P nodes
- RHOCPSUM: interpolation of RHOCP at P nodes
- ALPHASUM: interpolation of ALPHA at P nodes
- ALPHAFSUM: interpolation of ALPHAF at P nodes
- HRSUM: interpolation of HR at P nodes
- GGGPSUM: interpolation of GGGP at P nodes
- SXXSUM: interpolation of SXX at P nodes
- TKSUM: interpolation of TK at P nodes
- PHISUM: interpolation of PHI at P nodes
- DMPSUM: interpolation of DMP at P nodes
- DHPSUM: interpolation of DHP at P nodes
- XWSSUM: interpolation of XWS at P nodes
- WTPSUM: interpolation weights at P nodes
Erebus.setup_marker_geometry_helpers
— MethodSet up additional marker geometry helpers to facilitate marker handling.
setup_marker_geometry_helpers()
Details
- nothing
Returns
- mdis: minimum distance of marker launch anchor points to nearest marker
- mnum: number of marker nearest to marker launch anchor positions
Erebus.setup_marker_properties
— MethodSet up geodesic and physical properties of the set of markers.
setup_marker_properties(marknum; randomized)
Details
- randomized: fill in random values for grid properties instead of zeros
Returns
- xm : horizontal marker coordinate [m]
- ym : vertical marker coordinate [m]
- tm : marker material type
- tkm : marker temperature [K]
- sxxm : marker σ′xx [Pa]
- sxym : marker σxy [Pa]
- etavpm : marker viscoplastic viscosity [Pa]
- phim : marker porosity
Erebus.setup_marker_properties_helpers
— MethodSet up additional helper marker properties to facility comptuations.
setup_marker_properties_helpers(marknum; randomized)
Details
- randomized: fill in random values for grid properties instead of zeros
Returns
- rhototalm: total density of markers
- rhocptotalm : total volumetric heat capacity of markers
- etatotalm: total viscosity of markers
- hrtotalm: total radiogenic heat production of markers
- ktotalm: total thermal conductivity of markers
- tkm_rhocptotalm: total thermal energy of markers
- etafluidcur_inv_kphim: fluid viscosity over permeability of markers
- inv_gggtotalm: inverse of total shear modulus of markers
- fricttotalm: total friction coefficient of markers
- cohestotalm: total compressive strength of markers
- tenstotalm: total tensile strength of markers
- rhofluidcur: fluid density of markers
- alphasolidcur: solid thermal expansion coefficient of markers
- alphafluidcur: fluid thermal expansion coefficient of markers
Erebus.setup_staggered_grid_properties
— MethodSet up staggered grid properties for basic, Vx, Vy, and P nodes.
setup_staggered_grid_properties(; randomized)
Details
- randomized: fill in random values for grid properties instead of zeros
Returns
- ETA : viscoplastic viscosity at basic nodes [Pa*s]
- ETA0 : viscous viscosity at basic nodes [Pa*s]
- GGG : shear modulus at basic nodes [Pa]
- EXY : ϵxy at basic nodes [1/s]
- SXY : σxy at basic nodes [1/s]
- SXY0 : σ₀xy at basic nodes [1/s]
- wyx : rotation rate at basic nodes [1/s]
- COH : compressive strength at basic nodes [Pa]
- TEN : tensile strength at basic nodes [Pa]
- FRI : friction at basic nodes
- YNY : plastic yielding node property at basic nodes
- RHOX : density at Vx nodes [kg/m^3]
- RHOFX : fluid density at Vx nodes [kg/m^3]
- KX : thermal conductivity at Vx nodes [W/m/K]
- PHIX : porosity at Vx nodes
- vx : solid vx-velocity at Vx nodes [m/s]
- vxf : fluid vx-velocity at Vx nodes [m/s]
- RX : etafluid/kphi ratio at Vx nodes [kg m⁻³s⁻¹]
- qxD : qx-darcy flux at Vx nodes [m/s]
- gx : gx-gravity at Vx nodes [m/s^2]
- RHOY : density at Vy nodes [kg/m^3]
- RHOFY : fluid density at Vy nodes [kg/m^3]
- KY : thermal conductivity at Vy nodes [W/m/K]
- PHIY : porosity at Vy nodes
- vy : solid vy-velocity at Vy nodes [m/s]
- vyf : fluid vy-velocity at Vy nodes [m/s]
- RY : etafluid/kphi ratio at Vy nodes [kg m⁻³s⁻¹]
- qyD : qy-Darcy flux at Vy nodes [m/s]
- gy : gy-gravity at Vy nodes [m/s^2]
- RHO : density at P nodes [kg/m^3]
- RHOCP : volumetric heat capacity at P nodes [J/m^3/K]
- ALPHA : thermal expansion at P nodes [K⁻¹]
- ALPHAF : fluid thermal expansion at P nodes [K⁻¹]
- HR : radioactive heating at P nodes [W/m^3]
- HA : adiabatic heating at P nodes [W/m^3]
- HS : shear heating at P nodes [W/m^3]
- ETAP : viscosity at P nodes [Pa*s]
- GGGP : shear modulus at P nodes [Pa]
- EXX : ϵxx at P nodes [1/s]
- SXX : σ′xx at P nodes [1/s]
- SXX0 : σ₀′xx at P nodes [1/s]
- tk1 : current temperature at P nodes [K]
- tk2 : next temperature at P nodes [K]
- DT: temperature difference at P nodes [K]
- DT0: previous temperature difference at P nodes [K]
- vxp : solid vx in pressure nodes at P nodes [m/s]
- vyp : solid vy in pressure nodes at P nodes [m/s]
- vxpf : fluid vx in pressure nodes at P nodes [m/s]
- vypf : fluid vy in pressure nodes at P nodes [m/s]
- pr : total pressure at P nodes [Pa]
- pf : fluid pressure at P nodes [Pa]
- ps : solid pressure at P nodes [Pa]
- pr0 : previous total pressure at P nodes [Pa]
- pf0 : previous fluid pressure at P nodes [Pa]
- ps0 : previous solid pressure at P nodes [Pa]
- ETAPHI : bulk viscosity at P nodes [Pa*s]
- BETAPHI : bulk compresibility at P nodes [Pa*s]
- PHI : porosity at P nodes
- APHI : Dln at P nodes [(1-ϕ)/ϕ]/Dt
- FI : gravity potential at P nodes [J/kg]
- DMP: mass transfer term at P nodes
- DHP: enthalpy transfer/latent heating term at P nodes
- XWS: wet solid fraction at P nodes
Erebus.setup_staggered_grid_properties_helpers
— MethodSet up additional helper staggered grid properties to facilitate computations.
setup_staggered_grid_properties_helpers(; randomized)
Details
- randomized: fill in random values for grid properties instead of zeros
Returns
- ETA5: plastic iterations viscoplastic viscosity at basic nodes [Pa⋅s]
- ETA00: previous viscous viscosity at basic nodes [Pa⋅s]
- YNY5: plastic iterations plastic yielding node property at basic nodes
- YNY00: previous plastic yielding node property at basic nodes
- DSXY: stress change Δσxy at basic nodes [Pa]
- DSY: (SIIB-syield) at basic nodes
- EII :second strain rate invariant at P nodes [1/s]
- SII :second stress invariant at P nodes [Pa]
- DSXX :stress change Δσ′xx at P nodes [Pa]
- tk0: previous temperature at P nodes [K]
Erebus.setup_thermal_lse
— MethodSet up thermal linear system of equations structures.
setup_thermal_lse()
Details
- nothing
Returns
- RT: thermal linear system of equations: RHS vector
- ST: thermal linear system of equations: solution vector
Erebus.simulation_loop
— MethodMain simulation loop: run calculations with timestepping.
simulation_loop(output_path)
Details
- output_path: Absolute path where to save simulation output files
Returns
- nothing
Erebus.symmetrize_p_node_observables!
— MethodApply symmetry to P node observables.
symmetrize_p_node_observables!(SXX, APHI, PHI, pr, pf, ps)
Details
- SXX: σ′xx at P nodes
- APHI: Aϕ = Dln[(1-ϕ)/ϕ]/Dt at P nodes
- PHI: porosity at P nodes
- pr: total pressure at P nodes
- pf: fluid pressure at P nodes
- ps: solid pressure at P nodes
Returns
- nothing
Erebus.total
— MethodCompute convex combination of fluid and solid properties to get total property.
total(solid, fluid, ϕ)
Details
- fluid: fluid properties
- solid: solid properties
- ϕ: porosity (fraction of fluid)
Returns
- total: computed total property
Erebus.update_marker_population_geometry!
— MethodUpdate marker population geometry status given a marker number and nearest top/left marker grid point.
update_marker_population_geometry!(
m,
i,
j,
xm,
ym,
mdis,
mnum
)
Details
- m: marker number
- i: vertical index of top/left marker grid point
- j: horizontal index of top/left marker grid point
- xm: horizontal x-position of markers
- ym: vertical y-position of markers
- mdis: minimum distance of marker launch anchor points to nearest marker
- mnum: number of marker nearest to marker launch anchor positions
Returns
- nothing
Erebus.update_marker_porosity!
— MethodUpdate marker porosity for compaction based on Dln[(1-ϕ)/ϕ]/Dt at P grid.
update_marker_porosity!(xm, ym, tm, phim, APHI, dt, marknum)
Details
- xm: x-coordinates of markers
- ym: y-coordinates of markers
- tm: marker type
- phim: marker porosity
- APHI: Dln[(1-ϕ)/ϕ]/Dt at P nodes
- dt: computational time step
- marknum: total number of markers in use
Returns
- nothing
Erebus.update_marker_stress!
— MethodUpdate marker stress based on xy (basic grid) and xx (P grid) stress changes.
update_marker_stress!(
xm,
ym,
sxxm,
sxym,
DSXX,
DSXY,
marknum
)
Details
- xm: x-coordinates of markers
- ym: y-coordinates of markers
- sxxm: marker σ′xx [Pa]
- sxym: marker σxy [Pa]
- DSXX: stress change Δσ′xx at P nodes
- DSXY: stress change Δσxy at basic nodes
- marknum: total number of markers in use
Returns
- nothing
Erebus.update_marker_temperature!
— MethodUpdate marker temperature based on P grid temperature changes.
update_marker_temperature!(
xm,
ym,
tkm,
DT,
tk2,
timestep,
marknum
)
Details
- xm: x-coordinates of markers
- ym: y-coordinates of markers
- tkm: marker temperature
- DT: temperature change at P nodes
- tk2: next temperature at P nodes
- timestep: current time step
- marknum: total number of markers in use
Returns
- nothing
Erebus.update_marker_viscosity!
— MethodUpdate marker viscoplastic viscosity based on basic node yielding status and marker temperature- and material-based viscosity.
update_marker_viscosity!(
m,
xm,
ym,
tm,
tkm,
etatotalm,
etavpm,
YNY,
YNY_inv_ETA
)
Details
- m: marker number
- xm: x-coordinates of markers
- ym: y-coordinates of markers
- tm: type of markers
- tkm: temperature of markers
- etatotalm: total viscosity of markers
- etavpm: matrix viscosity of markers
- YNY: plastic yielding node property at basic nodes
- YNY_inv_ETA: inverse viscoplastic viscosity at yielding basic nodes
Returns
-nothing
Erebus.update_p_nodes_melt_composition!
— MethodCompute P nodes wet silicate (solid) fraction from markers.
update_p_nodes_melt_composition!(
xm,
ym,
XWsolidm0,
XWS,
XWSSUM,
WTPSUM,
marknum
)
Details
- xm: x-coordinate of markers
- ym: y-coordinate of markers
- XWsolidm0: previous marker wet silicate (solid) molar fraction
- XWS: wet silicate (solid) molar fraction at P nodes
- XWSSUM: wet silicate (solid) molar fraction interpolated to P nodes
- WTPSUM: weight for bilinear interpolation to P nodes
- marknum: current total number of markers in use
Returns
- nothing
Erebus.ηᶠcur_inv_kᵠ
— MethodCompute inverse of porosity-dependent permeability (eqn 16.64 in Gerya (2019)) times current fluid viscosity.
ηᶠcur_inv_kᵠ(kϕᵣ, ϕ, ηᶠcur)
Details
- kϕᵣ: reference permeability [m^2]
- ϕ: current porosity
- ηᶠcur: current fluid viscosity [Pa s]
Returns
- etafluidcur_inv_kphi: inverse empirical porosity-dependent permeability
times current fluid viscosity