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.

source

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)

source

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)

source

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.

source

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.

source

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.

source

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.

source

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.

source

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.

source

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

source

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!.

source

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

source

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.

source

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.

source

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!.

source

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

source

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. The topography is linearly interpolated between chain.cell_vertices and chain.h_vertices. 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.

source

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 potential Q with 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 of Q with respect to pressure (drives volumetric plastic strain rate ε_vol_pl = -λ · dQdP).

• dFdP — scalar gradient of the yield function F with respect to pressure; combined with dQdP it enters the volume-closure term K·dt·dFdP·dQdP of 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.

source

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::Integer picks 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 of F).

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.

source

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.

source

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)

source

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)

source

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.

source

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.

source

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.

source

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.

source

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.

source

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.

source

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

source

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!.

source

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

source

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.

source

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.

source

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!.

source

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

source

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. The topography is linearly interpolated between chain.cell_vertices and chain.h_vertices. 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.

source

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 potential Q with 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 of Q with respect to pressure (drives volumetric plastic strain rate ε_vol_pl = -λ · dQdP).

• dFdP — scalar gradient of the yield function F with respect to pressure; combined with dQdP it enters the volume-closure term K·dt·dFdP·dQdP of 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.

source

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::Integer picks 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 of F).

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.

source