Rheology and phases
Phase ratios, buoyancy, viscosity, melt fraction, dissolved volatiles, the yield function and plastic gradients, and marker-chain topography correction.
JustRelax.JustRelax2D.fn_ratio Method
fn_ratio(fn::F, rheology::NTuple{N, AbstractMaterialParamsStruct}, ratio) where {N, F}Average the function fn over the material phases in rheology using the phase ratios ratio.
JustRelax.JustRelax2D.update_phase_ratios_2D! Method
update_phase_ratios_2D!(
phase_ratios::JustPIC.PhaseRatios, phase_arrays::NTuple{N, AbstractMatrix}, xci, xvi
)JustRelax routine based on JustPIC.update_phase_ratios!. Update the center, vertex and velocity-face phase ratios in phase_ratios from the 2-D phase_arrays, given the cell-center coordinates xci and vertex coordinates xvi. The phase arrays need to be AbstractArrays and have values between 0 and 1.
#Example:
nx, ny = 100, 100
phase_1 = zeros(nx, ny)
phase_1[User_criterion .== true] .= 1.0
phase_2 = zeros(nx, ny)
phase_2[User_criterion .== false] .= 1.0
phase_arrays = (phase_1, phase_2)
# Advect both phase arrays and update phase ratios
update_phase_ratios_2D!(phase_ratios, phase_arrays, xci, xvi)JustRelax.JustRelax2D.update_phase_ratios_3D! Method
update_phase_ratios_3D!(
phase_ratios::JustPIC.PhaseRatios, phase_arrays::NTuple{N, AbstractArray}, xci, xvi
)JustRelax routine based on JustPIC.update_phase_ratios!. Update the center, vertex, velocity-face and shear-stress-midpoint phase ratios in phase_ratios from the 3-D phase_arrays, given the cell-center coordinates xci and vertex coordinates xvi. The phase arrays need to be AbstractArrays and have values between 0 and 1.
#Example:
nx, ny, nz = 100, 100, 100
phase_1 = zeros(nx, ny, nz)
phase_1[User_criterion .== true] .= 1.0
phase_2 = zeros(nx, ny, nz)
phase_2[User_criterion .== false] .= 1.0
phase_arrays = (phase_1, phase_2)
# Advect both phase arrays and update phase ratios
update_phase_ratios_3D!(phase_ratios, phase_arrays, xci, xvi)JustRelax.JustRelax2D.compute_buoyancy Method
compute_buoyancy(rheology, args, phase_ratios)Compute the buoyancy forces based on the given rheology, arguments, and phase ratios.
Arguments
rheology: The rheology used to compute the buoyancy forces.args: Additional arguments required by the rheology.phase_ratios: The ratios of the different phases.
JustRelax.JustRelax2D.compute_buoyancy Method
compute_buoyancy(rheology, args)Compute the buoyancy forces based on the given rheology and arguments.
Arguments
rheology: The rheology used to compute the buoyancy forces.args: Additional arguments required for the computation.
JustRelax.JustRelax2D.compute_buoyancy Method
compute_buoyancy(rheology::MaterialParams, args, phase_ratios)Compute the buoyancy forces for a given set of material parameters, arguments, and phase ratios.
Arguments
rheology: The material parameters.args: The arguments.phase_ratios: The phase ratios.
JustRelax.JustRelax2D.compute_buoyancy Method
compute_buoyancy(rheology::MaterialParams, args)Compute the buoyancy forces based on the given rheology parameters and arguments.
Arguments
rheology::MaterialParams: The material parameters for the rheology.args: The arguments for the computation.
JustRelax.JustRelax2D.compute_ρg! Method
compute_ρg!(ρg, rheology, args)Calculate the buoyance forces ρg for the given GeoParams.jl rheology object and correspondent arguments args.
JustRelax.JustRelax2D.compute_ρg! Method
compute_ρg!(ρg, phase_ratios, rheology, args; air_phase = 0)Calculate the buoyance forces ρg for the given GeoParams.jl rheology object and correspondent arguments args. The phase_ratios are used to compute the density of the composite rheology.
A non-zero air_phase drops that phase from the average and renormalizes over the remaining ones, so ρg is the buoyancy of the rock alone rather than of the rock-air mixture filling the cell. Solvers that weight ρg by a rock volume fraction need this: with the mixture density the rock fraction is counted twice.
JustRelax.JustRelax2D.compute_viscosity! Method
compute_viscosity!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place by evaluating rheology at the strain-rate invariant (equivalent to compute_viscosity_εII!; see there for the arguments, and compute_viscosity_τII! for the stress-invariant alternative).
JustRelax.JustRelax2D.compute_viscosity_εII! Method
compute_viscosity_εII!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place from the second invariant of the strain rate (εII); otherwise identical to compute_viscosity_τII!.
JustRelax.JustRelax2D.compute_viscosity_τII! Method
compute_viscosity_τII!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place from the second invariant of the deviatoric stress (τII), evaluating rheology (a single GeoParams.MaterialParams, or one per phase when phase_ratios is given) at each cell and relaxing towards the new value with factor relaxation (1.0 = no damping). cutoff = (ηmin, ηmax) clamps the result. air_phase (multi-phase form only) excludes that phase from the update.
See also compute_viscosity_εII! for the strain-rate-invariant convention, and compute_viscosity! for the rheology-driven default (εII).
JustRelax.JustRelax2D.viscosity_phase_ratio Method
viscosity_phase_ratio(air_phase, ratio)Phase ratio to average viscosity over, with air_phase dropped and the remaining phases renormalized. A cell holding nothing but air keeps its own ratio: averaging over no phase at all would make the harmonic mean Inf, which then spreads through ητ into neighbouring cells that do carry rock.
JustRelax.JustRelax2D.compute_melt_fraction! Method
compute_melt_fraction!(ϕ, rheology, args)
compute_melt_fraction!(ϕ, phase_ratios::JustPIC.PhaseRatios, rheology, args)
compute_melt_fraction!(ϕ, dϕdT, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill the melt-fraction array ϕ from the GeoParams melting parameterisation of rheology, with args supplying the state variables it needs (typically P and T, as scalars or index-matched arrays). Given phase_ratios, the melt fraction is averaged over the phases present in each cell.
The five-argument form also fills dϕdT, the temperature derivative of the melt fraction, in the same pass. dϕdT is what activates GeoParams.Latent_HeatCapacity: it contributes Q_L * dϕdT to Cp, and defaults to zero when absent from the args handed to the thermal kernels. ϕ is clamped to [0, 1] and dϕdT is not, so a melting law with apply_bounds = false can return a flat ϕ alongside a nonzero dϕdT.
There is no fused single-phase form: pair the three-argument method with compute_melt_fraction_derivative!, at the cost of a second pass over the grid.
JustRelax.JustRelax2D.compute_melt_fraction_derivative! Method
compute_melt_fraction_derivative!(dϕdT, rheology, args)
compute_melt_fraction_derivative!(dϕdT, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill dϕdT, the temperature derivative of the melt fraction, without touching ϕ. Use this when ϕ is advected on particles or otherwise not recomputed on the grid; when it is, the fused compute_melt_fraction!(ϕ, dϕdT, …) does both in one pass.
In-place GeoParams equivalent: compute_dϕdT!.
JustRelax.JustRelax2D.compute_dissolved_volatiles! Method
compute_dissolved_volatiles!(mH2O, mCO2, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill the dissolved H2O and CO2 mass-fraction arrays from the GeoParams solubility closures (Liu2005_Solubility, Mafic_Solubility). Mirrors compute_melt_fraction! but writes two arrays, because compute_dissolved returns the (m_h2o, m_co2) pair. args supplies P, T, and the CO2 mole fraction of the gas X_co2 (scalars or index-matched arrays).
JustRelax.JustRelax2D.update_phases_given_markerchain! Method
update_phases_given_markerchain!(phase, chain::MarkerChain, particles::Particles, origin, di, air_phase, args = ())Deactivate the particles that end up on the wrong side of the free surface tracked by chain: air-phase particles below it and rock particles above it. Their coordinates and every field in args are set to NaN and their index entry to false, so that particle injection re-seeds those cells from their neighbours.
origin and di are the origin and grid spacing of the particle grid, and air_phase the phase index standing for air.
JustRelax.JustRelax2D.compute_plastic_gradients_phase Method
compute_plastic_gradients_phase(rheology, phase, τij; P, τII, EII, …) -> (dQdτ, dQdP, dFdP)
compute_plastic_gradients_phase(rheology, ratio, τij; P, τII, EII, …) -> (dQdτ, dQdP, dFdP)Return the plastic flow gradients required by the return-mapping update:
dQdτ::NTuple— gradient of the plastic potentialQwith respect to the deviatoric stress tensorτij, in tensor convention (shear slots already halved, i.e.ε_pl_xy = λ * dQdτ[xy]directly without a factor of 1/2).dQdP— scalar gradient ofQwith respect to pressure (drives volumetric plastic strain rateε_vol_pl = -λ · dQdP).dFdP— scalar gradient of the yield functionFwith respect to pressure; combined withdQdPit enters the volume-closure termK·dt·dFdP·dQdPof the λ denominator.
τij is the deviatoric stress at which to evaluate the gradients (typically the trial stress). The slot ordering matches @stress(stokes):
2D:
(xx, yy, xy)3D:
(xx, yy, zz, yz, xz, xy)
The second positional argument selects the phase the same way as compute_yieldfunction_phase (single phase::Integer or per-phase ratio). For composite rheologies without a plastic element, all three return values are zero.
JustRelax.JustRelax2D.compute_yieldfunction_phase Method
compute_yieldfunction_phase(rheology, phase; P, τII, EII, …)
compute_yieldfunction_phase(rheology, ratio; P, τII, EII, …)Evaluate the plastic yield function F of the CompositeRheology for a single phase or for a phase-weighted mixture at a given stress state.
rheology is the tuple of MaterialParams for all phases. The second positional argument selects the phase:
phase::Integerpicks the rheology of a single phase.ratio::NTuple{N}/ratio::SVector{N}provides per-phase volume fractions; phases with zero weight are skipped and the remaining contributions are summed (linear blend ofF).
Keyword arguments are forwarded to GeoParams.compute_yieldfunction of the plastic primitive (typically P, τII, EII). For composite rheologies without a plastic element, F = τII is returned so a weighted sum is not artificially driven to zero by elastic-only phases.
JustRelax.JustRelax3D.fn_ratio Method
fn_ratio(fn::F, rheology::NTuple{N, AbstractMaterialParamsStruct}, ratio) where {N, F}Average the function fn over the material phases in rheology using the phase ratios ratio.
JustRelax.JustRelax3D.update_phase_ratios_2D! Method
update_phase_ratios_2D!(
phase_ratios::JustPIC.PhaseRatios, phase_arrays::NTuple{N, AbstractMatrix}, xci, xvi
)JustRelax routine based on JustPIC.update_phase_ratios!. Update the center, vertex and velocity-face phase ratios in phase_ratios from the 2-D phase_arrays, given the cell-center coordinates xci and vertex coordinates xvi. The phase arrays need to be AbstractArrays and have values between 0 and 1.
#Example:
nx, ny = 100, 100
phase_1 = zeros(nx, ny)
phase_1[User_criterion .== true] .= 1.0
phase_2 = zeros(nx, ny)
phase_2[User_criterion .== false] .= 1.0
phase_arrays = (phase_1, phase_2)
# Advect both phase arrays and update phase ratios
update_phase_ratios_2D!(phase_ratios, phase_arrays, xci, xvi)JustRelax.JustRelax3D.update_phase_ratios_3D! Method
update_phase_ratios_3D!(
phase_ratios::JustPIC.PhaseRatios, phase_arrays::NTuple{N, AbstractArray}, xci, xvi
)JustRelax routine based on JustPIC.update_phase_ratios!. Update the center, vertex, velocity-face and shear-stress-midpoint phase ratios in phase_ratios from the 3-D phase_arrays, given the cell-center coordinates xci and vertex coordinates xvi. The phase arrays need to be AbstractArrays and have values between 0 and 1.
#Example:
nx, ny, nz = 100, 100, 100
phase_1 = zeros(nx, ny, nz)
phase_1[User_criterion .== true] .= 1.0
phase_2 = zeros(nx, ny, nz)
phase_2[User_criterion .== false] .= 1.0
phase_arrays = (phase_1, phase_2)
# Advect both phase arrays and update phase ratios
update_phase_ratios_3D!(phase_ratios, phase_arrays, xci, xvi)JustRelax.JustRelax3D.compute_buoyancy Method
compute_buoyancy(rheology, args, phase_ratios)Compute the buoyancy forces based on the given rheology, arguments, and phase ratios.
Arguments
rheology: The rheology used to compute the buoyancy forces.args: Additional arguments required by the rheology.phase_ratios: The ratios of the different phases.
JustRelax.JustRelax3D.compute_buoyancy Method
compute_buoyancy(rheology, args)Compute the buoyancy forces based on the given rheology and arguments.
Arguments
rheology: The rheology used to compute the buoyancy forces.args: Additional arguments required for the computation.
JustRelax.JustRelax3D.compute_buoyancy Method
compute_buoyancy(rheology::MaterialParams, args, phase_ratios)Compute the buoyancy forces for a given set of material parameters, arguments, and phase ratios.
Arguments
rheology: The material parameters.args: The arguments.phase_ratios: The phase ratios.
JustRelax.JustRelax3D.compute_buoyancy Method
compute_buoyancy(rheology::MaterialParams, args)Compute the buoyancy forces based on the given rheology parameters and arguments.
Arguments
rheology::MaterialParams: The material parameters for the rheology.args: The arguments for the computation.
JustRelax.JustRelax3D.compute_ρg! Method
compute_ρg!(ρg, rheology, args)Calculate the buoyance forces ρg for the given GeoParams.jl rheology object and correspondent arguments args.
JustRelax.JustRelax3D.compute_ρg! Method
compute_ρg!(ρg, phase_ratios, rheology, args; air_phase = 0)Calculate the buoyance forces ρg for the given GeoParams.jl rheology object and correspondent arguments args. The phase_ratios are used to compute the density of the composite rheology.
A non-zero air_phase drops that phase from the average and renormalizes over the remaining ones, so ρg is the buoyancy of the rock alone rather than of the rock-air mixture filling the cell. Solvers that weight ρg by a rock volume fraction need this: with the mixture density the rock fraction is counted twice.
JustRelax.JustRelax3D.compute_viscosity! Method
compute_viscosity!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place by evaluating rheology at the strain-rate invariant (equivalent to compute_viscosity_εII!; see there for the arguments, and compute_viscosity_τII! for the stress-invariant alternative).
JustRelax.JustRelax3D.compute_viscosity_εII! Method
compute_viscosity_εII!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place from the second invariant of the strain rate (εII); otherwise identical to compute_viscosity_τII!.
JustRelax.JustRelax3D.compute_viscosity_τII! Method
compute_viscosity_τII!(stokes::StokesArrays, [phase_ratios,] args, rheology, cutoff; air_phase=0, relaxation=1.0)Update stokes.viscosity.η in place from the second invariant of the deviatoric stress (τII), evaluating rheology (a single GeoParams.MaterialParams, or one per phase when phase_ratios is given) at each cell and relaxing towards the new value with factor relaxation (1.0 = no damping). cutoff = (ηmin, ηmax) clamps the result. air_phase (multi-phase form only) excludes that phase from the update.
See also compute_viscosity_εII! for the strain-rate-invariant convention, and compute_viscosity! for the rheology-driven default (εII).
JustRelax.JustRelax3D.viscosity_phase_ratio Method
viscosity_phase_ratio(air_phase, ratio)Phase ratio to average viscosity over, with air_phase dropped and the remaining phases renormalized. A cell holding nothing but air keeps its own ratio: averaging over no phase at all would make the harmonic mean Inf, which then spreads through ητ into neighbouring cells that do carry rock.
JustRelax.JustRelax3D.compute_melt_fraction! Method
compute_melt_fraction!(ϕ, rheology, args)
compute_melt_fraction!(ϕ, phase_ratios::JustPIC.PhaseRatios, rheology, args)
compute_melt_fraction!(ϕ, dϕdT, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill the melt-fraction array ϕ from the GeoParams melting parameterisation of rheology, with args supplying the state variables it needs (typically P and T, as scalars or index-matched arrays). Given phase_ratios, the melt fraction is averaged over the phases present in each cell.
The five-argument form also fills dϕdT, the temperature derivative of the melt fraction, in the same pass. dϕdT is what activates GeoParams.Latent_HeatCapacity: it contributes Q_L * dϕdT to Cp, and defaults to zero when absent from the args handed to the thermal kernels. ϕ is clamped to [0, 1] and dϕdT is not, so a melting law with apply_bounds = false can return a flat ϕ alongside a nonzero dϕdT.
There is no fused single-phase form: pair the three-argument method with compute_melt_fraction_derivative!, at the cost of a second pass over the grid.
JustRelax.JustRelax3D.compute_melt_fraction_derivative! Method
compute_melt_fraction_derivative!(dϕdT, rheology, args)
compute_melt_fraction_derivative!(dϕdT, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill dϕdT, the temperature derivative of the melt fraction, without touching ϕ. Use this when ϕ is advected on particles or otherwise not recomputed on the grid; when it is, the fused compute_melt_fraction!(ϕ, dϕdT, …) does both in one pass.
In-place GeoParams equivalent: compute_dϕdT!.
JustRelax.JustRelax3D.compute_dissolved_volatiles! Method
compute_dissolved_volatiles!(mH2O, mCO2, phase_ratios::JustPIC.PhaseRatios, rheology, args)Fill the dissolved H2O and CO2 mass-fraction arrays from the GeoParams solubility closures (Liu2005_Solubility, Mafic_Solubility). Mirrors compute_melt_fraction! but writes two arrays, because compute_dissolved returns the (m_h2o, m_co2) pair. args supplies P, T, and the CO2 mole fraction of the gas X_co2 (scalars or index-matched arrays).
JustRelax.JustRelax3D.update_phases_given_markerchain! Method
update_phases_given_markerchain!(phase, chain::MarkerChain, particles::Particles, origin, di, air_phase, args = ())Deactivate the particles that end up on the wrong side of the free surface tracked by chain: air-phase particles below it and rock particles above it. Their coordinates and every field in args are set to NaN and their index entry to false, so that particle injection re-seeds those cells from their neighbours.
origin and di are the origin and grid spacing of the particle grid, and air_phase the phase index standing for air.
JustRelax.JustRelax3D.compute_plastic_gradients_phase Method
compute_plastic_gradients_phase(rheology, phase, τij; P, τII, EII, …) -> (dQdτ, dQdP, dFdP)
compute_plastic_gradients_phase(rheology, ratio, τij; P, τII, EII, …) -> (dQdτ, dQdP, dFdP)Return the plastic flow gradients required by the return-mapping update:
dQdτ::NTuple— gradient of the plastic potentialQwith respect to the deviatoric stress tensorτij, in tensor convention (shear slots already halved, i.e.ε_pl_xy = λ * dQdτ[xy]directly without a factor of 1/2).dQdP— scalar gradient ofQwith respect to pressure (drives volumetric plastic strain rateε_vol_pl = -λ · dQdP).dFdP— scalar gradient of the yield functionFwith respect to pressure; combined withdQdPit enters the volume-closure termK·dt·dFdP·dQdPof the λ denominator.
τij is the deviatoric stress at which to evaluate the gradients (typically the trial stress). The slot ordering matches @stress(stokes):
2D:
(xx, yy, xy)3D:
(xx, yy, zz, yz, xz, xy)
The second positional argument selects the phase the same way as compute_yieldfunction_phase (single phase::Integer or per-phase ratio). For composite rheologies without a plastic element, all three return values are zero.
JustRelax.JustRelax3D.compute_yieldfunction_phase Method
compute_yieldfunction_phase(rheology, phase; P, τII, EII, …)
compute_yieldfunction_phase(rheology, ratio; P, τII, EII, …)Evaluate the plastic yield function F of the CompositeRheology for a single phase or for a phase-weighted mixture at a given stress state.
rheology is the tuple of MaterialParams for all phases. The second positional argument selects the phase:
phase::Integerpicks the rheology of a single phase.ratio::NTuple{N}/ratio::SVector{N}provides per-phase volume fractions; phases with zero weight are skipped and the remaining contributions are summed (linear blend ofF).
Keyword arguments are forwarded to GeoParams.compute_yieldfunction of the plastic primitive (typically P, τII, EII). For composite rheologies without a plastic element, F = τII is returned so a weighted sum is not artificially driven to zero by elastic-only phases.