Erebus.jl

Compute 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]

Add 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

Apply 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

Apply 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

Apply 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

Assemble 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

Assemble 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

Assemble 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

Backtrack 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

Compute 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]

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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).

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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)

Compute 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)

Compute 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

Compute 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

Compute 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

Compute 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

Assess 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

Compute 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

Compute 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

Compute 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

Compute 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

Compute 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

Define 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

Calculate 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]

Compute 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

Compute 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

Decide 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

Decide 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

Compute 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

Compute 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

Compute 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]

Compute 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

Get 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

Get 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]]

Initialize 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

Interpolate 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

Interpolate 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

Interpolate 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

Compute 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]

Compute 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]

Interpolate 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

Interpolate 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

Interpolate 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

Interpolate 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

Interpolate 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

Move 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

Parse command line arguments and feed them to the main function.

parse_commandline()

Details:

- nothing

Returns

- parsed_args: parsed command line arguments

Perform 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

Perform 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

Compare 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

Process 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

Process 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

Recompute 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

Reduce 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

Add 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

Reset 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

Reset 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

Runs the simulation with the given parameters.

run_simulation()

Details

- nothing

Returns

- nothing

Convert seconds to Ma (millions of years).

s_to_Ma(s)

Details

- s: period in seconds

Returns

- Ma: period in millions of years

Save 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

Set 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

Set 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

Set 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

Set 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

Set 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

Set 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

Set 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

Set 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

Set 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]

Set 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

Main simulation loop: run calculations with timestepping.

simulation_loop(output_path)

Details

- output_path: Absolute path where to save simulation output files

Returns

- nothing

Apply 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

Compute 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

Update 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

Update 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

Update 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

Update 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

Update 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

Compute 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

Compute 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