Index

Every exported function and type in one flat list. For the same functions grouped by topic, see the other pages under API reference in the sidebar.

JustRelax.JustRelax Module

julia

JustRelax

Matrix-free, accelerated pseudo-transient solvers for geodynamic multi-physics problems.

The root module holds what is shared across dimensions and backends: the backend types (CPUBackend, CUDABackend, AMDGPUBackend) and their traits, the staggered-grid Geometry, the MPI topology handle IGG, the field containers, and the boundary-condition types. Solvers and kernels live in the dimension-specific submodules JustRelax.JustRelax2D and JustRelax.JustRelax3D, one of which a model loads alongside this module; I/O lives in JustRelax.DataIO.

Material properties come from GeoParams.jl, particle advection from JustPIC.jl, kernel parallelism from ParallelStencil.jl, and distributed-memory decomposition from ImplicitGlobalGrid.jl.

JustRelax.AMDGPUBackend Type

julia

AMDGPUBackend

Backend tag selecting ROCArrays, running on an AMD GPU. Requires AMDGPU.jl to be loaded (see Selecting the backend).

JustRelax.AbstractBackend Type

julia

AbstractBackend

Supertype for backend tags (CPUBackend, CUDABackend, AMDGPUBackend) selecting which array type and device a model runs on.

JustRelax.AbstractBoundaryConditions Type

julia

AbstractBoundaryConditions

Supertype for all boundary condition types, e.g. TemperatureBoundaryConditions.

JustRelax.AbstractFlowBoundaryConditions Type

julia

AbstractFlowBoundaryConditions

Supertype for velocity/displacement boundary condition types (VelocityBoundaryConditions, DisplacementBoundaryConditions).

JustRelax.BackendTrait Type

julia

BackendTrait

Supertype for the traits backend(x) returns to dispatch solver code on the array type backing x, independent of x's own type hierarchy.

JustRelax.CPUBackend Type

julia

CPUBackend

Backend tag selecting plain Arrays, running on the CPU via ParallelStencil's Threads target. The default backend.

JustRelax.CPUBackendTrait Type

julia

CPUBackendTrait

Trait returned by backend(x) when x is backed by a plain Array.

JustRelax.CUDABackend Type

julia

CUDABackend

Backend tag selecting CuArrays, running on an Nvidia GPU. Only defined once CUDA.jl is loaded (see Selecting the backend).

JustRelax.DYREL Type

julia

struct DYREL{T, F}

Structure containing parameters and arrays for the DYREL (Dynamic Relaxation) solver.

Fields

γ_eff: Effective penalty parameter.
Dx, Dy, Dz: Diagonal preconditioners for velocity updates in x, y, (and z) directions.
λmaxVx, λmaxVy, λmaxVz: Maximum eigenvalues for stability calculation.
dVxdτ, dVydτ, dVzdτ: Pseudo-time step related damping terms.
dτVx, dτVy, dτVz: Pseudo-time steps for velocity fields.
dVx, dVy, dVz: Velocity increments for the current iteration.
βVx, βVy, βVz: Damping coefficients for momentum equation.
cVx, cVy, cVz: Damping coefficients related to dynamic relaxation.
αVx, αVy, αVz: Scaling factors for damping.
ηb: Bulk viscosity field.
P_num: Numerical pressure scratch field.
Rx0, Ry0, Rz0: Velocity residual history scratch fields.
CFL: Courant-Friedrichs-Lewy number.
γfact: Penalty scaling factor.
ϵ: General convergence tolerance.
ϵ_vel: Velocity convergence tolerance.
c_fact: Damping scaling factor.

JustRelax.DisplacementBoundaryConditions Type

julia

DisplacementBoundaryConditions(; no_slip, free_slip, periodic, free_surface=false)

Define 2D or 3D boundary conditions for the displacement field. Each face is controlled independently with a named tuple. Use exactly four faces in 2D (left, right, top, bot) or six in 3D (also front, back).

A face can carry at most one of no_slip, free_slip, and periodic; a face where all three are false is left for the caller to prescribe explicitly. Periodic faces must be enabled in pairs (left/right, front/back, or bot/top), and a periodic top face is incompatible with free_surface=true.

JustRelax.Geometry Type

julia

struct Geometry{nDim,V,D,T}

A staggered Cartesian grid in nDim dimensions.

Geometry stores the domain size, origin, cell spacing, cell-centered coordinates, vertex coordinates, and the staggered velocity grids used throughout JustRelax.

JustRelax.Geometry Method

julia

Geometry(TA, xvi::Vararg{<:AbstractVector, nDim})
Geometry(xvi::NTuple{nDim, <:AbstractVector})

Build a staggered grid from explicit vertex coordinates along each dimension.

This constructor is useful for refined or otherwise nonuniform meshes. Cell-centered coordinates, local spacings, and staggered velocity grids are derived from the supplied vertex coordinates. TA can be used to move the generated arrays to a target array type.

Arguments

TA: Array constructor used to materialize the coordinate arrays, for example Array or a backend-specific array type.
xvi: One vertex-coordinate vector per dimension.

JustRelax.Geometry Method

julia

Geometry(ni, li; origin = ntuple(_ -> 0.0, Val(nDim)))

Build a uniform staggered grid with ni cells and physical domain lengths li.

When ImplicitGlobalGrid has been initialized, the grid spacing is computed from the global grid dimensions and the returned coordinates correspond to the local MPI subdomain. Otherwise a serial grid covering the full domain is created.

Arguments

ni: Number of local grid cells in each direction.
li: Physical domain length in each direction.

Keywords

origin: Lower-left or lower-front corner of the domain.

Returns

A Geometry with cell-centered coordinates xci, vertex coordinates xvi, and staggered velocity coordinates xi_vel.

JustRelax.GeometryAnnulus Type

julia

GeometryAnnulus{nDim,V,D,T} <: AbstractGrid{nDim,V,D,T}

A staggered annular grid whose first two coordinate directions are (θ, r): angular position followed by radius. Accordingly, xci[1] and xvi[1] contain angular coordinates in radians, while xci[2] and xvi[2] contain radial coordinates.

The grid stores its cell counts, coordinate extents, origin, spacings, cell-centered coordinates, vertex coordinates, and staggered velocity grids.

JustRelax.GeometryAnnulus Method

julia

GeometryAnnulus(TA, θv, rv)
GeometryAnnulus((θv, rv))

Build a staggered annular grid from explicit vertex-coordinate vectors in (θ, r) order. θv contains angular vertices and rv contains radial vertices. Cell-centered coordinates and local spacings are derived from these vectors. All values in θv must be given in radians.

TA is the array constructor used to materialize the generated coordinate arrays, such as Array or a backend-specific array type. The tuple form uses Array.

JustRelax.GeometryAnnulus Method

julia

GeometryAnnulus(ni, li; origin = ntuple(_ -> 0.0, Val(nDim)))

Build a uniform staggered annular grid in (θ, r) coordinate order.

For a two-dimensional annulus, ni = (nθ, nr), li = (lθ, lr), and origin = (θ₀, r₀). Thus, the first entries describe the angular direction and the second entries describe the radial direction. Both lθ and θ₀ must be given in radians. When ImplicitGlobalGrid is initialized, the coordinates cover the local MPI subdomain; otherwise they cover the full domain.

Arguments

ni: Number of cells in each coordinate direction, ordered as (nθ, nr).
li: Coordinate extents, ordered as angular extent then radial extent.

Keywords

origin: Lower coordinate bounds, ordered as (θ₀, r₀).

Returns

A GeometryAnnulus with cell-centered, vertex, and staggered velocity coordinates in (θ, r) order.

JustRelax.IGG Type

julia

IGG(me, dims, nprocs, coords, comm_cart)

Container for the Cartesian MPI topology returned by ImplicitGlobalGrid.init_global_grid.

This is typically created as:

julia

igg = IGG(init_global_grid(nx, ny, nz; init_MPI = true)...)

and then passed around so code can access the current rank, Cartesian coordinates, and communicator associated with the distributed grid decomposition.

JustRelax.NonCPUBackendTrait Type

julia

NonCPUBackendTrait

Trait returned by backend(x) for any AbstractArray other than Array that isn't a CuArray/ROCArray — the generic fallback for GPU array types without a dedicated trait.

JustRelax.PTStokesCoeffs Type

julia

PTStokesCoeffs(li, di; ϵ_rel=1e-6, ϵ_abs=1e-12, Re=3π, CFL=0.9/√2.1, r=0.7)

Pseudo-transient damping coefficients for the Stokes solver, derived from the domain size li, grid spacing di, Reynolds number Re and bulk-to-shear damping ratio r following Räss et al. (2022). Passed as pt_stokes to solve!.

ηdτ / ητ is the local velocity pseudo-time step, while θ_dτ controls the stress and pressure updates. In the 2D free-surface-stabilized velocity kernel, the vertical pseudo-time step also includes the local diagonal -dt * ∂y(ρg) introduced by stabilization.

Keyword arguments

ϵ_rel, ϵ_abs: relative/absolute convergence tolerances.
Re: Reynolds number.
CFL: Courant-Friedrichs-Lewy number bounding the pseudo-time step.
r: ratio of the damping coefficients for the bulk and shear rheology.

JustRelax.RockRatio Type

julia

RockRatio{T, N} <: AbstractMask

A struct representing the rock ratio mask used in variational Stokes solvers. It contains arrays for center and vertex values, as well as velocity components and shear components.

JustRelax.StressParticles Type

julia

StressParticles{backend, nNormal, nShear, T}

Particle-borne deviatoric stress and vorticity: the normal components τ_normal, the shear components τ_shear, and the vorticity components ω, each a tuple of particle cell arrays. Carrying the old stress on the particles instead of on the grid keeps it attached to the material as it advects and rotates.

Build one from the particles it follows with StressParticles(particles), advance it with rotate_stress!, and write it back onto stokes.τ_o with stress2grid!.

JustRelax.TemperatureBoundaryConditions Type

julia

TemperatureBoundaryConditions(; no_flux, constant_flux, constant_value, periodic, dirichlet)

Create thermal boundary conditions for 2D or 3D temperature fields.

Boundary tuples use left, right, top, and bot in 2D. In 3D they also use front and back. Omitted faces are filled with false, and the dimensionality is inferred from the longest boundary tuple that is passed: pass a complete six-face tuple for a 3D set, since the defaults are four-face 2D tuples. Tuples with any other number of faces are rejected.

The face values have the following meaning:

no_flux: true copies the adjacent interior temperature into the ghost layer.
constant_value: numeric values prescribe the boundary temperature through the ghost value Tghost = 2 * value - Tinterior.
constant_flux: numeric values prescribe heat fluxes in the pseudo-transient diffusion flux kernels.
periodic: true copies the opposite interior temperature into the ghost layer.
false: leaves that boundary inactive for the corresponding condition.

Periodic faces must be enabled in pairs (left/right, front/back, or bot/top) and cannot also use no_flux, constant_flux, or constant_value.

dirichlet accepts the mask-based Dirichlet forms supported by Dirichlet, for example (; constant = value, mask = mask).

Examples

julia

TemperatureBoundaryConditions(;
    no_flux = (left = true, right = true, top = false, bot = false),
    constant_value = (top = 273.0, bot = 1573.0),
)

TemperatureBoundaryConditions(;
    no_flux = (left = false, right = false, top = false, bot = false),
    constant_value = (top = 273.0, bot = 1573.0),
    periodic = (left = true, right = true, top = false, bot = false),
)

JustRelax.VelocityBoundaryConditions Type

julia

VelocityBoundaryConditions(; no_slip, free_slip, periodic, free_surface=false)

Define 2D or 3D boundary conditions for the velocity field. Face names are left, right, top, and bot in 2D, with front and back added in 3D.

no_slip, free_slip, and periodic are mutually exclusive on each face. Faces where all three are false are not modified by flow_bcs!, which is how a prescribed velocity field is imposed. Periodic faces must be enabled in pairs (left/right, front/back, or bot/top), and a periodic top face is incompatible with free_surface=true.

JustRelax.PTArray Method

julia

PTArray()
PTArray(::Type{<:AbstractBackend})

The array type associated with a backend tag: Array for CPUBackend, CuArray for CUDABackend, ROCArray for AMDGPUBackend. PTArray() (no argument) is the CPU default.

JustRelax.check_flow_bcs Method

julia

check_flow_bcs(no_slip, free_slip, periodic, free_surface)

Throw if flow boundary conditions conflict or if a periodic direction is not paired. A boundary flagged as neither no_slip, free_slip, nor periodic is left untouched by flow_bcs!, which is how a prescribed velocity field is imposed: the caller writes the boundary and ghost values itself.

JustRelax.check_periodic_bcs Method

julia

check_periodic_bcs(stokes, bcs, igg, di_center)

Throw unless the periodic directions of bcs are ones the solver can actually solve.

The momentum row of a periodic seam only exists if stokes was allocated for it, which is what StokesArrays(backend, ni, bcs) does; a stokes built without the boundary conditions leaves that row out and the seam velocity is then frozen at its initial value. The other two conditions mark combinations that are not implemented rather than ones that are wrong in principle.

JustRelax.lazy_grid Method

julia

lazy_grid(di, ni, Li; origin = ntuple(_ -> zero(T1), Val(N)))

Create cell-centered and vertex coordinates for a serial uniform grid.

di gives the spacing in each direction, ni the number of cells, and Li the physical lengths of the domain.

JustRelax.lazy_grid_MPI Method

julia

lazy_grid_MPI(di, ni; origin = ntuple(_ -> zero(T1), Val(N)))

Create local cell-centered and vertex coordinates for a uniform grid distributed with ImplicitGlobalGrid.

The returned coordinates are shifted by origin and correspond to the local MPI rank.

JustRelax.legacy_uniform_grid Method

julia

legacy_uniform_grid(ni, di)

Construct a uniform Geometry from grid sizes ni and cell spacings di.

This helper preserves the older API used by some solver code. In MPI mode the physical domain lengths are reconstructed from the global grid dimensions, so the resulting geometry matches the full distributed domain rather than only the local chunk.

JustRelax.periodic_dims Function

julia

periodic_dims(bcs::AbstractFlowBoundaryConditions)
periodic_dims(stokes::StokesArrays)

Periodicity of each spatial direction, as an N-tuple of Bools ordered (x, y[, z]).

A boundary-condition object reports what the caller asked for. A StokesArrays reports what its momentum residuals were allocated for: in a periodic direction the two coincident boundary faces are a single unknown, so that direction carries one residual row per cell rather than one per interior face. check_periodic_bcs requires the two to agree.

JustRelax.plot_field Function

julia

plot_field(data, index, grid; kwargs...)

Heatmap of one field of a 2D cell array, e.g. the volume fraction of a single phase, and return the Makie.Figure. Requires a Makie backend (e.g. using GLMakie or using CairoMakie) to be loaded.

Arguments

data: cell array holding the field, e.g. phase_ratios.center.
index: which entry of each cell to plot, e.g. the phase number.
grid: tuple of coordinates matching data, as LinRanges. The xci and xvi of a uniform Geometry qualify; a nonuniform grid stores its coordinates as arrays and is not accepted.

Keyword arguments

colormap: colormap of the heatmap. Default: :roma.
title: plot title. Default: "Field Plot".
filename: if given, the figure is also written to this path. Default: nothing.
resolution: figure size in pixels. Default: (1200, 1000).
labelsize, titlesize: font sizes of the axis labels and title. Default: 35, 50.
conversion: divisor applied to the coordinates. Default: 1.0e3, i.e. m to km.
units: unit shown in the axis labels. Default: :km.

Example

julia

f = plot_field(phase_ratios.center, 2, xci; title = "Phase 2 fraction")

JustRelax.plot_particles Function

julia

plot_particles(particles, pPhases; kwargs...)

Scatter the particle positions of a 2D model, colored by phase, and return the Makie.Figure. Requires a Makie backend (e.g. using GLMakie or using CairoMakie) to be loaded.

Arguments

particles: JustPIC.Particles whose coordinates are plotted.
pPhases: per-particle phase, as a cell array.

Keyword arguments

chain: JustPIC.MarkerChain drawn on top of the particles, e.g. a free surface. Default: nothing.
clrmap: colormap for the phases. Default: :roma.
title: plot title. Default: "Particle Position".
filename: if given, the figure is also written to this path. Default: nothing.
resolution: figure size in pixels. Default: (1200, 1200).
labelsize, titlesize: font sizes of the axis labels and title. Default: 35, 50.
linecolor, markersize: color and marker size of the marker chain. Default: :black, 1.0.
conversion: divisor applied to the coordinates. Default: 1.0e3, i.e. m to km.
units: unit shown in the axis labels. Default: :km.

JustRelax.reject_periodic_bcs Method

julia

reject_periodic_bcs(bcs, solver)

Throw if bcs asks for a periodic direction, naming the solver that does not implement it.

Only the DYREL solvers write the momentum row for a periodic seam: solve_DYREL! in 2D and 3D, and solve_VariationalDYREL! in 2D. Every other path would leave that row unsolved and silently return a field pinned to its initial guess at the seam.

JustRelax.unwrap Method

julia

unwrap(x::StressParticles)

Flatten x into a single tuple (τ_normal..., τ_shear..., ω...) of its underlying particle cell arrays.

JustRelax.velocity_grids Method

julia

velocity_grids(xci, xvi, di)

Build staggered velocity coordinates from cell-centered and vertex grids.

For each velocity component, the coordinate along that component lives on vertices, while the transverse directions are extended with one ghost point on either side. Both uniform spacings and nonuniform spacing vectors are supported in 2D and 3D.

Arguments

xci: Cell-centered coordinates in each direction.
xvi: Vertex coordinates in each direction.
di: Cell spacing as either scalars for a uniform grid or vectors for a nonuniform grid.

JustRelax.versioninfo Function

julia

versioninfo(io::IO=stdout; verbose::Bool=false)

Print information about the version of JustRelax in use. The output includes:

JustRelax version and installation method
Git commit information (if available)
Platform information
Julia version
Key dependencies (verbose mode)
Environment variables (verbose mode)

The output is controlled with boolean keyword arguments:

verbose: print all additional information including dependencies and environment

See also: Base.versioninfo().

JustRelax.x_g Method

julia

x_g(idx, dxi, nxi_or_A)

Global x coordinate of local index idx on a grid with spacing dxi, given either the local axis length nxi or an array A sized along that axis. Accounts for the MPI domain decomposition and periodicity set up by init_global_grid (see IGG).

JustRelax.y_g Method

julia

y_g(idx, dxi, nxi_or_A)

Global y coordinate of local index idx, analogous to x_g.

JustRelax.z_g Method

julia

z_g(idx, dxi, nxi_or_A)

Global z coordinate of local index idx, analogous to x_g.

JustRelax.@dx Macro

julia

@dx(dx, I)

x spacing at grid index I, reading from dx — see @dxi.

JustRelax.@dxi Macro

julia

@dxi(dxi, I, J[, K])

Cell spacing (dx, dy[, dz]) at grid index (I, J[, K]), reading from dxi — a tuple of per-axis spacing that may be a constant Number (uniform grid) or an AbstractVector (non-uniform grid, indexed at I/J/K).

JustRelax.@dy Macro

julia

@dy(dy, J)

y spacing at grid index J, reading from dy — see @dxi.

JustRelax.@dz Macro

julia

@dz(dz, K)

z spacing at grid index K, reading from dz — see @dxi.

JustRelax.JustRelax2D Module

julia

JustRelax.JustRelax2D

Two-dimensional solvers, kernels, and constructors, on the CPU backend.

The submodule is loaded with using JustRelax.JustRelax2D, and its ParallelStencil environment is initialized for two dimensions when it loads. Loading CUDA.jl or AMDGPU.jl before it adds the matching GPU methods through a package extension; the entry points and their signatures stay the same, and the backend is chosen by the array type of the containers passed in. See JustRelax.JustRelax3D for the three-dimensional counterpart.

JustRelax.JustRelax2D.DYREL Function

julia

DYREL(ni::NTuple{N, Integer}[, periodic]; ϵ=1e-6, ϵ_vel=1e-6, CFL=0.99, c_fact=0.5, γfact=20.0) where N

Creates a new DYREL struct with fields initialized to zero.

Arguments

ni: Tuple containing the grid dimensions (nx, ny) for 2D or (nx, ny, nz) for 3D.
periodic: N-tuple marking the periodic directions, which each carry one extra momentum row (see momentum_rows). Defaults to all-false. The StokesArrays method below reads it off the containers instead, so the two cannot disagree.

Keyword arguments

ϵ: General convergence tolerance. Default: 1.0e-6.
ϵ_vel: Velocity convergence tolerance. Default: 1.0e-6.
CFL: Courant-Friedrichs-Lewy number. Default: 0.99.
c_fact: Damping scaling factor. Default: 0.5.
γfact: Penalty scaling factor. Default: 20.0.

JustRelax.JustRelax2D.DYREL! Method

julia

DYREL!(dyrel::JustRelax.DYREL, stokes::JustRelax.StokesArrays, rheology, phase_ratios, di, dt; CFL=dyrel.CFL, γfact=dyrel.γfact)

Updates the fields of the DYREL struct in-place for the current time step.

This function recomputes:

Bulk viscosity and penalty parameter γ_eff.
Gershgorin estimates for eigenvalues and preconditioners.
Damping coefficients.

Arguments

dyrel: JustRelax.DYREL struct to modify.
stokes: JustRelax.StokesArrays containing current simulation state.
rheology, phase_ratios: Material properties.
di: Grid spacing.
dt: Current time step.
CFL: Courant number (default: the value stored in dyrel).
γfact: Penalty factor (default: the value stored in dyrel).

Returns nothing.

JustRelax.JustRelax2D.DYREL Method

julia

DYREL(stokes, rheology, phase_ratios, di, dt; ϵ=1e-6, ϵ_vel=1e-6, CFL=0.99, c_fact=0.5, γfact=20.0)

Constructs and initializes a DYREL object based on existing Stokes fields.

This function:

Allocates zero-initialized arrays using grid dimensions from stokes.
Computes initial bulk viscosity and penalty parameters.
Computes Gershgorin estimates for eigenvalues and preconditioners.
Updates damping coefficients.

Arguments

stokes: JustRelax.StokesArrays struct.
rheology: Material properties.
phase_ratios: Phase fraction information.
di: Grid spacing tuple.
dt: Time step.

Keyword arguments

ϵ: General convergence tolerance. Default: 1.0e-6.
ϵ_vel: Velocity convergence tolerance. Default: 1.0e-6.
CFL: Courant-Friedrichs-Lewy number. Default: 0.99.
c_fact: Damping scaling factor. Default: 0.5.
γfact: Factor for the penalty parameter calculation. Default: 20.0.

JustRelax.JustRelax2D.Displacement Method

julia

Displacement(nx::Integer, ny::Integer, nz::Integer)

Create the displacement arrays for the Stokes solver in 3D.

Fields

Ux: Displacement in x direction at their staggered location
Uy: Displacement in y direction at their staggered location
Uz: Displacement in z direction at their staggered location

JustRelax.JustRelax2D.Displacement Method

julia

Displacement(nx::Integer, ny::Integer)

Create the displacement arrays for the Stokes solver in 2D.

Fields

Ux: Displacement in x direction at their staggered location
Uy: Displacement in y direction at their staggered location

JustRelax.JustRelax2D.PTThermalCoeffs Method

julia

PTThermalCoeffs(rheology, phase_ratios, args, dt, ni, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from a multi-phase rheology.

The coefficient arrays are evaluated on the thermal grid of size ni using the phase ratios in phase_ratios.center together with the thermodynamic state stored in args.

JustRelax.JustRelax2D.PTThermalCoeffs Method

julia

PTThermalCoeffs(K, ρCp, dt, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from conductivity and volumetric heat-capacity arrays.

di and li are the grid spacing and domain-length tuples used to estimate the local pseudo-time step and relaxation factor stored in the returned JustRelax.PTThermalCoeffs.

JustRelax.JustRelax2D.PTThermalCoeffs Method

julia

PTThermalCoeffs(rheology, args, dt, ni, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from a single rheology state.

This method evaluates conductivity and ρCp from rheology and the cellwise fields in args, then stores the resulting relaxation arrays in a JustRelax.PTThermalCoeffs object.

JustRelax.JustRelax2D.PrincipalStress Method

julia

PrincipalStress(ni::NTuple{N, Integer}) where {N}

Create the principal stress arrays for the Stokes solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

Fields

σ1: First principal stress
σ2: Second principal stress
σ3: Third principal stress (only in 3D). In 2D it is a placeholder array of size (2, 1, 1).

JustRelax.JustRelax2D.Residual Method

julia

Residual(nx::Integer, ny::Integer[, nz::Integer])
Residual(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool})

Create the residual arrays for the Stokes solver.

Fields

Rx, Ry[, Rz]: Residuals for the momentum equations
RP: Residual for the continuity equation

periodic marks the directions whose two boundary faces are the same plane; each of those gains one momentum row (see momentum_rows). It defaults to all-false.

JustRelax.JustRelax2D.RockRatio Method

julia

RockRatio(nx, ny, nz)

Create a RockRatio object for a 3D grid with dimensions nx x ny x nz on a staggered grid.

JustRelax.JustRelax2D.RockRatio Method

julia

RockRatio(nx, ny)

Create a RockRatio object for a 2D grid with dimensions nx x ny on a staggered grid.

JustRelax.JustRelax2D.RockRatio Method

julia

RockRatio(backend, ni)
RockRatio(backend, ni...)

Allocate a RockRatio on backend for a staggered grid of ni cells, with every volume fraction initialized to zero.

backend is the same backend type passed to the other allocators (CPUBackend, CUDABackend, AMDGPUBackend), and ni is (nx, ny) in 2D or (nx, ny, nz) in 3D. Fill the fractions with update_rock_ratio!, or with JustPIC.compute_rock_fraction! when the liquid domain is bounded by a marker chain.

JustRelax.JustRelax2D.StokesArrays Method

julia

StokesArrays(ni::NTuple{N,Integer}) where {N}
StokesArrays(ni::NTuple{N,Integer}, bcs::AbstractFlowBoundaryConditions)
StokesArrays(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool})

Create the Stokes arrays object in 2D or 3D.

Passing the flow boundary conditions sizes the momentum residuals for the periodic directions they declare, which is what gives a periodic seam a momentum row; without them every direction is taken to be non-periodic. The solvers check the two against each other, so the boundary conditions have to be built first.

Fields

P: Pressure field
P0: Previous pressure field
∇V: Velocity gradient
V: Velocity fields
Q: Volumetric source/sink term e.g. ΔV/V_tot [m³/m³]
U: Displacement fields
ω: Vorticity field
τ: Stress tensors
τ_o: Old stress tensors
ε: Strain rate tensors
ε_pl: Plastic strain rate tensors
EII_pl: Second invariant of the accumulated plastic strain
viscosity: Viscosity fields
R: Residual fields
Δε: Strain increment tensor
∇U: Displacement gradient
λ : plastic multiplier @ centers
λv : plastic multiplier @ vertices
λv_yz, λv_xz, λv_xy : 3D plastic multiplier on shear staggered grids
ΔPψ : pressure correction in dilatant case

JustRelax.JustRelax2D.StressParticles Method

julia

StressParticles(particles::Particles)

Allocate the stress and vorticity cell arrays that follow particles, on the same backend and with the same per-cell capacity. Two normal and one shear component in 2-D, three of each in 3-D.

JustRelax.JustRelax2D.SymmetricTensor Method

julia

SymmetricTensor(nx::Integer, ny::Integer, nz::Integer)

Create the symmetric tensor arrays for the Stokes solver in 3D.

Fields

xx: xx component of the tensor at cell centers
yy: yy component of the tensor at cell centers
zz: zz component of the tensor at cell centers
xx_v: xx component of the tensor at vertices
yy_v: yy component of the tensor at vertices
zz_v: zz component of the tensor at vertices
xy: xy component of the tensor at vertices
yz: yz component of the tensor at vertices
xz: xz component of the tensor at vertices
yz_c: yz component of the tensor at cell centers
xz_c: xz component of the tensor at cell centers
xy_c: xy component of the tensor at cell centers
II: second invariant of the tensor at cell centers

JustRelax.JustRelax2D.SymmetricTensor Method

julia

SymmetricTensor(nx::Integer, ny::Integer)

Create the symmetric tensor arrays for the Stokes solver in 2D.

Fields

xx: xx component of the tensor at cell centers
yy: yy component of the tensor at cell centers
xx_v: xx component of the tensor at vertices
yy_v: yy component of the tensor at vertices
xy: xy component of the tensor at vertices
xy_c: xy component of the tensor at cell centers
II: second invariant of the tensor at cell centers

JustRelax.JustRelax2D.ThermalArrays Method

julia

ThermalArrays(nx::Integer, ny::Integer, nz::Integer)

3D constructor for the thermal arrays for the heat diffusion solver with the extents given by nx, ny and nz.

Fields

T: Temperature at cell centers with one ghost node on every boundary (nx + 2, ny + 2, nz + 2)
Told: Temperature at previous time step at cell centers with ghost nodes (nx + 2, ny + 2, nz + 2)
ΔT: Temperature change at cell centers with ghost nodes (nx + 2, ny + 2, nz + 2)
adiabatic: Adiabatic term α (u ⋅ ∇P) at cell centers (nx, ny, nz)
dT_dt: Time derivative of temperature at cell centers (nx, ny, nz)
qTx: Conductive heat flux in x direction on cell faces (nx + 1, ny, nz)
qTy: Conductive heat flux in y direction on cell faces (nx, ny + 1, nz)
qTz: Conductive heat flux in z direction on cell faces (nx, ny, nz + 1)
qTx2: Conductive heat flux in x direction on cell faces for second order scheme (nx + 1, ny, nz)
qTy2: Conductive heat flux in y direction on cell faces for second order scheme (nx, ny + 1, nz)
qTz2: Conductive heat flux in z direction on cell faces for second order scheme (nx, ny, nz + 1)
H: Source terms at cell centers (nx, ny, nz)
shear_heating: Shear heating terms at cell centers (nx, ny, nz)
ResT: Residual of the temperature equation at cell centers (nx, ny, nz)

JustRelax.JustRelax2D.ThermalArrays Method

julia

ThermalArrays(nx::Integer, ny::Integer)

2D constructor for the thermal arrays for the heat diffusion solver with the extents given by nx and ny.

Fields

T: Temperature at cell centers with one ghost node on every boundary (nx + 2, ny + 2)
Told: Temperature at previous time step at cell centers with ghost nodes (nx + 2, ny + 2)
ΔT: Temperature change at cell centers with ghost nodes (nx + 2, ny + 2)
adiabatic: Adiabatic term α (u ⋅ ∇P) at cell centers (nx, ny)
dT_dt: Time derivative of temperature at cell centers (nx, ny)
qTx: Conductive heat flux in x direction on cell faces (nx + 1, ny)
qTy: Conductive heat flux in y direction on cell faces (nx, ny + 1)
qTx2: Conductive heat flux in x direction on cell faces for second order scheme (nx + 1, ny)
qTy2: Conductive heat flux in y direction on cell faces for second order scheme (nx, ny + 1)
H: Source terms at cell centers (nx, ny)
shear_heating: Shear heating terms at cell centers (nx, ny)
ResT: Residual of the temperature equation at cell centers (nx, ny)

JustRelax.JustRelax2D.ThermalArrays Method

julia

ThermalArrays(ni::NTuple{N, Integer}) where {N}

Create the thermal arrays for the heat diffusion solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

JustRelax.JustRelax2D.ThermalArrays Method

julia

ThermalArrays(::Backend, ni::NTuple{N, Integer}) where {N}

Internal entry point function for the ThermalArrays constructor. This allows for dispatching on the backend type and then calling the main constructor with the dimensions.

JustRelax.JustRelax2D.Velocity Method

julia

Velocity(nx::Integer, ny::Integer, nz::Integer)

Create the velocity arrays for the Stokes solver in 3D.

Fields

Vx: Velocity in x direction (nx + 1, ny + 2, nz + 2)
Vy: Velocity in y direction (nx + 2, ny + 1, nz + 2)
Vz: Velocity in z direction (nx + 2, ny + 2, nz + 1)

JustRelax.JustRelax2D.Velocity Method

julia

Velocity(nx::Integer, ny::Integer)

Create the velocity arrays for the Stokes solver in 2D.

Fields

Vx: Velocity in x direction (nx + 1, ny + 2)
Vy: Velocity in y direction (nx + 2, ny + 1)

JustRelax.JustRelax2D.Viscosity Method

julia

Viscosity(ni::NTuple{N, Integer}) where {N}

Create the viscosity arrays for the Stokes solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

Fields

η: Viscosity at cell centers
ηv: Viscosity at vertices
η_vep: Viscosity for visco-elastic-plastic rheology
ητ: Pseudo-transient viscosity for stress update

JustRelax.JustRelax2D.Vorticity Method

julia

Vorticity(nx::Integer, ny::Integer, nz::Integer)

Create the vorticity arrays for the Stokes solver in 3D.

Fields

yz: Vorticity component yz at their staggered location
xz: Vorticity component xz at their staggered location
xy: Vorticity component xy at their staggered location
yz_c: Vorticity component yz at the cell centers
xz_c: Vorticity component xz at the cell centers
xy_c: Vorticity component xy at the cell centers

JustRelax.JustRelax2D.Vorticity Method

julia

Vorticity(nx::Integer, ny::Integer)

Create the vorticity arrays for the Stokes solver in 2D.

Fields

xy: Vorticity component xy at vertices

JustRelax.JustRelax2D.WENO5 Method

julia

WENO5(backend, method::Val, ni::NTuple)
WENO5(method::Val, ni::NTuple)

Allocate a WENO5 advection scheme state for a grid of size ni. method selects the smoothness-indicator variant: Val(1) for Jiang-Shu, Val(2) for the Z-variant.

JustRelax.JustRelax2D.WENO_advection! Method

julia

WENO_advection!(u, Vxi, weno, di, ni, dt)

Perform the advection step of the Weighted Essentially Non-Oscillatory (WENO) scheme for the solution of hyperbolic partial differential equations.

Arguments

u: field to be advected.
Vxi: velocity field.
weno: structure containing the WENO scheme parameters and temporary variables.
di: grid spacing.
ni: number of grid points.
dt: time step.

Description

The function approximates the advected fluxes using the WENO scheme and use a strong-stability preserving (SSP) Runge-Kutta method of order 3 for the time integration.

JustRelax.JustRelax2D._bc_coordinate Method

julia

_bc_coordinate(A, x)

Materialize coordinate vector x on the same device as array A, so that background-field kernels can index it. collect is applied first because a lazy range cannot be transferred to a device array directly.

JustRelax.JustRelax2D._clamped_index Method

julia

clamped_indices(ni, [periodic,] I...)

Cell indices of the stencil around vertex I on a grid of ni cells, kept in range.

Outside a periodic direction an index that falls off the grid is clamped onto its in-range neighbour, so a vertex-centred average degenerates to the one-sided average of the cells that do exist. A direction listed in periodic has no edge: the cell on the far side of the seam is a real neighbour, so the index wraps onto it and both copies of the seam plane see the same stencil. Omitting periodic clamps every direction.

JustRelax.JustRelax2D._heatdiffusion_PT! Method

julia

_heatdiffusion_PT!(thermal, pt_thermal, thermal_bc, K, ρCp, dt, grid;
    igg, b_width, iterMax, nout, verbose)

Solve the heat equation with pseudo-transient iterations using precomputed, cell-centered material properties.

K is the thermal conductivity field and ρCp is the volumetric heat-capacity field on the thermal grid. pt_thermal supplies the pseudo-transient coefficients, thermal_bc applies the temperature boundary conditions after each update, and grid provides metric terms and halo layout.

Returns a named tuple containing the iteration numbers and residual norms sampled every nout iterations.

JustRelax.JustRelax2D._heatdiffusion_PT! Method

julia

_heatdiffusion_PT!(thermal, pt_thermal, thermal_bc, rheology, args, dt, grid;
    igg, phase, stokes, b_width, iterMax, nout, verbose)

Solve the heat equation with pseudo-transient iterations using thermal properties derived from rheology.

args is a named tuple of thermodynamic fields, typically including T and P. Entries sized like the thermal cell centers are read as-is; larger entries are offset by one to skip their ghost nodes. When phase is provided, pseudo-transient coefficients are recomputed from the local phase ratios each iteration. When stokes is provided, thermal.adiabatic is refreshed before the iteration loop to include the adiabatic heating contribution.

Returns a named tuple containing the sampled iteration counts and residual history.

JustRelax.JustRelax2D._update_rock_ratio! Method

julia

_update_rock_ratio!(ϕ, ratio, air_phase)

Inner kernel of update_rock_ratio that clamps the computed rock ratio to the range [0, 1] for the given ratio and air_phase.

JustRelax.JustRelax2D.accumulate_tensor! Method

julia

accumulate_tensor!(II, A::JustRelax.SymmetricTensor, dt)

Accumulate the second invariant of the symmetric tensor A over a time step: II[I] += dt * A_II[I].

The tensor components are gathered onto the cell centers where II lives, so A may hold its shear components on the vertices. Used to integrate the deviatoric plastic strain rate ε_pl into EII_pl; the volumetric counterpart is accumulate_vol!.

JustRelax.JustRelax2D.accumulate_vol! Method

julia

accumulate_vol!(EVol_pl, ε_vol_pl, dt)

Accumulate the volumetric plastic strain over a time step: EVol_pl[I] += dt * ε_vol_pl[I].

ε_vol_pl is the volumetric plastic strain rate at cell centers (= λ · (-dQ/dP), set inside the stress kernel), and is the volumetric counterpart of the deviatoric ε_pl. EVol_pl is the running invariant accumulated through time and is distinct from EII_pl (which integrates the second invariant of the deviatoric plastic strain rate via accumulate_tensor!).

JustRelax.JustRelax2D.adiabatic_heating! Method

julia

adiabatic_heating!(thermal, stokes, rheology, phases, _dt, grid)

Fill thermal.adiabatic with the adiabatic heating term inferred from the pressure change between stokes.P0 and stokes.P.

The kernels average the local thermal expansivity over the temperature nodes and scale the pressure increment by inv(dt), passed here as _dt. When stokes is nothing, the no-op overloads leave the field unchanged.

JustRelax.JustRelax2D.allzero Method

julia

allzero(x::Vararg{T,N}) where {T,N}

Check if all elements in x are zero.

Arguments

x::Vararg{T,N}: The input array.

Returns

Bool: true if all elements in x are zero, false otherwise.

JustRelax.JustRelax2D.apply_free_surface_diagonal! Method

julia

apply_free_surface_diagonal!(Dn, λmaxVn, ρgn, di_center, dt)

Add the free-surface diagonal -dt * ∂n(ρg) to the normal DYREL preconditioner and its Gershgorin row bound. The normal direction is y in 2D and z in 3D. Passing dt = 0 is a no-op.

JustRelax.JustRelax2D.assign! Method

julia

assign!(B::AbstractArray{T,N}, A::AbstractArray{T,N}) where {T,N}

Assigns the values of array A to array B in parallel.

Arguments

B::AbstractArray{T,N}: The destination array.
A::AbstractArray{T,N}: The source array.

JustRelax.JustRelax2D.center2vertex! Method

julia

center2vertex!(vertex, center)
center2vertex!(vertex_yz, vertex_xz, vertex_xy, center_yz, center_xz, center_xy)

Interpolates the values at the cell center(s) onto vertex points. The 6-argument method interpolates the three shear-stress/strain-rate components of a 3D SymmetricTensor onto their respective face vertices.

JustRelax.JustRelax2D.compute_P! Method

compute_P!(P, P0, RP, ∇V, Q, ΔT, η, rheology::NTuple{N,MaterialParams}, phase_ratio::C, dt, r, θ_dτ)

Compute the pressure field P and the residual RP for the compressible case. This function introduces thermal stresses after the implementation of Kiss et al. (2023).

Arguments

P: pressure field
RP: residual field
∇V: divergence of the velocity field
Q: volumetric source/sink term which should have the properties of dV/V_tot [m³/m³] normalized per cell, default is zero.
ΔT: temperature difference on the cell center, to account for thermal stresses. The thermal expansivity α is computed from the material parameters.
η: viscosity field
rheology: material parameters
phase_ratio: phase field
melt_fraction: melt fraction field, used for the thermal expansion coefficient if provided
dt: time step
r: relaxation parameter for the pressure update
θ_dτ: numerical parameter for the pressure update

JustRelax.JustRelax2D.compute_V! Method

julia

compute_V!(Vx, Vy, Vz, Rx, Ry, Rz, P, fx, fy, fz, τxx, τyy, τzz, τyz, τxz, τxy, ητ, ηdτ, ϕ, _dx, _dy, _dz)

Compute the 3D velocity field V from the pressure P, stress components τ, body forces f, and other parameters, with the rock ratio ϕ and grid spacing _dx, _dy, _dz.

JustRelax.JustRelax2D.compute_V! Method

julia

compute_V!(Vx, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ, _dx, _dy, dt)

Compute the velocity field V with the timestep dt from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ, grid spacing _dx, _dy, and time step dt.

JustRelax.JustRelax2D.compute_V! Method

julia

compute_V!(Vx, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ, _dx, _dy)

Compute the velocity field V from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

Each active face update uses its own liquid face mass. When the density-gradient free-surface correction is enabled, it is included implicitly in that local diagonal; inactive faces receive zero velocity and zero residual.

JustRelax.JustRelax2D.compute_Vx! Method

julia

compute_Vx!(Vx, Rx, P, τxx, τxy, ηdτ, ρgx, ητ, ϕ, _dx, _dy)

Compute the x-component of the velocity field Vx from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

JustRelax.JustRelax2D.compute_Vy! Method

julia

compute_Vy!(Vy, Vx_on_Vy, Ry, P, τyy, τxy, ηdτ, ρgy, ητ, ϕ, _dx, _dy, dt)

Compute the y-component of the velocity field Vy from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ, grid spacing _dx, _dy, and time step dt.

JustRelax.JustRelax2D.compute_air_ratio Method

julia

compute_air_ratio(phase_ratio, air_phase, inds...)

Compute the air ratio at the given indices based on the phase_ratio and air_phase.

JustRelax.JustRelax2D.compute_bulk_viscosity_and_penalty! Method

julia

compute_bulk_viscosity_and_penalty!(dyrel, stokes, rheology, phase_ratios, γfact, dt)

Computes the bulk viscosity ηb and the effective penalty parameter γ_eff.

Bulk Viscosity (ηb): Computed based on the bulk modulus of the material phases.

If Kb is infinite (incompressible), ηb defaults to γfact * η_mean.
Otherwise ηb = Kb * dt.

Penalty Parameter (γ_eff): A combination of numerical (γ_num) and physical (γ_phy) penalty terms.

γ_num = γfact * η_mean, with η_mean the mean of the finite viscosities
γ_phy = Kb * dt (or γ_num where Kb is infinite)
γ_eff = (γ_phy * γ_num) / (γ_phy + γ_num)

Arguments

dyrel: JustRelax.DYREL struct to update.
stokes: JustRelax.StokesArrays.
rheology: Material properties.
phase_ratios: Phase fraction information.
γfact: Numerical factor for penalty parameter (default: 20.0).
dt: Time step.

This function parallelizes the computation across grid cells.

JustRelax.JustRelax2D.compute_buoyancy Method

julia

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

julia

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

julia

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

julia

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_dissolved_volatiles! Method

julia

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.compute_dt Method

julia

compute_dt(S::JustRelax.StokesArrays, args...)

Compute the time step dt for the simulation.

JustRelax.JustRelax2D.compute_lithostatic_pressure! Method

julia

compute_lithostatic_pressure!(P, ρg, dz)
compute_lithostatic_pressure!(P, ρg, dz, igg::IGG)

Integrate the vertical component of the buoyancy force ρg down the columns of the cell-centered pressure P. The vertical direction is the last dimension of P and points upwards, so the last entry of a column is the shallowest cell.

dz is either a constant cell height or a vector holding the height of every cell along the vertical direction. Since P is cell-centered, cell j carries the weight of all the cells above it plus half of its own,

julia

P[j] = Σ_{k>j} ρg[k] * dz[k] + ρg[j] * dz[j] / 2

The three-argument method integrates each column within the local subdomain, and throws if the vertical direction is split across MPI ranks. Pass the MPI topology igg to also collect the weight of the cells held by the ranks stacked above the local subdomain; ρg must then be up to date on the halo cells, and the resulting P agrees on the cells that neighboring ranks share.

JustRelax.JustRelax2D.compute_maxloc! Method

julia

compute_maxloc!(B, A; window)

Compute the maximum value of A in the window = (width_x, width_y, width_z) and store the result in B.

JustRelax.JustRelax2D.compute_melt_fraction! Method

julia

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

julia

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_plastic_gradients_phase Method

julia

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.

JustRelax.JustRelax2D.compute_principal_stresses! Method

julia

compute_principal_stresses!(stokes, σ::PrincipalStress)

In-place version of compute_principal_stresses, writing into a pre-allocated σ.

JustRelax.JustRelax2D.compute_principal_stresses Method

julia

compute_principal_stresses(backend, stokes::StokesArrays)

Compute the principal deviatoric stresses (eigenvalues and eigenvectors of the stress tensor) at cell centers from stokes, returning a new PrincipalStress.

JustRelax.JustRelax2D.compute_rock_ratio Method

julia

compute_rock_ratio(phase_ratio, air_phase, inds...)

Compute the rock ratio at the given indices based on the phase_ratio and air_phase.

JustRelax.JustRelax2D.compute_shear_heating! Method

julia

compute_shear_heating!(thermal, stokes, rheology, dt)

Populate thermal.shear_heating from the current deviatoric stress and strain-rate fields stored in stokes.

The kernel evaluates the elastic strain-rate contribution from stokes.τ, stokes.τ_o, the shear modulus in rheology, and the pseudo-time step dt. The resulting volumetric heating term is written in place on thermal cell centers.

When phase ratios are passed as an extra positional argument, phase-weighted material properties are used instead of a single rheology state.

JustRelax.JustRelax2D.compute_strain_rate! Method

julia

compute_strain_rate!(εxx, εyy, εzz, εyz, εxz, εxy, ∇V, Vx, Vy, Vz, ϕ, _dx, _dy, _dz)

Compute the 3D components of the strain rate tensor ε from the velocity field V and its divergence ∇V, taking into account the rock ratio ϕ and grid spacing _dx, _dy, _dz.

JustRelax.JustRelax2D.compute_strain_rate! Method

julia

compute_strain_rate!(εxx, εyy, εxy, ∇V, Vx, Vy, ϕ, _dx, _dy)

Compute the components of the strain rate tensor ε from the velocity field V and its divergence ∇V, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

ϕ selects which entries are computed at all; it does not scale them. The rock fraction enters the momentum balance once, where the stress divergence is taken (d_xa(τxx, ϕ.center), d_yi(τxy, ϕ.vertex)). Scaling ε by it as well would weight the deviatoric term by ϕ² while pressure and buoyancy carry ϕ, which softens partially filled cells relative to the rest of the momentum equation.

JustRelax.JustRelax2D.compute_strain_rate_from_increment! Method

julia

compute_strain_rate_from_increment!(εxx, εyy, εxy, Δεxx, Δεyy, Δεxy, ϕ, _dt)

Compute the components of the strain rate tensor ε from the strain increments Δε, taking into account the rock ratio ϕ and time step _dt.

As in compute_strain_rate!, ϕ only selects which entries are computed; the rock fraction is applied once, in the stress divergence.

JustRelax.JustRelax2D.compute_viscosity! Method

julia

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

julia

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

julia

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.compute_yieldfunction_phase Method

julia

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.

JustRelax.JustRelax2D.compute_α Method

julia

compute_α(rheology, phase)

Return the thermal expansivity α used by the adiabatic heating kernels.

phase can be a single phase index, nothing, or a phase-ratio vector. In the latter case the result is phase-weighted.

JustRelax.JustRelax2D.compute_ρg! Method

julia

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

julia

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_∇V! Method

julia

compute_∇V!(∇V, V, ϕ, _di)

Compute the divergence of the velocity field V and store it in ∇V, taking into account the rock ratio ϕ and grid spacing _di.

The divergence is evaluated only for retained pressure cells. The local continuity degree of freedom is eliminated when one of its surrounding staggered velocity faces is disconnected:

julia

             Vy[i, j+1]
                   o
                   |
    Vx[i, j]  o--- p[i,j] ---o  Vx[i+1, j]
                   |
                   o
             Vy[i, j]

This is null-space elimination, not an air velocity equation.

JustRelax.JustRelax2D.continuation_linear Method

julia

continuation_linear(x_new, x_old, ν)

Do a continuation step (1-ν)*x_old + ν*x_new with damping parameter ν

JustRelax.JustRelax2D.continuation_log Method

julia

continuation_log(x_new, x_old, ν)

Do a continuation step exp((1-ν)*log(x_old) + ν*log(x_new)) with damping parameter ν

JustRelax.JustRelax2D.displacement2velocity! Method

julia

displacement2velocity!(stokes::StokesArrays, dt)
displacement2velocity!(stokes::StokesArrays, dt, flow_bcs::AbstractFlowBoundaryConditions)

Set stokes.V (velocity) to stokes.U / dt (displacement divided by the time step), in place — the inverse of velocity2displacement!. The 3-argument method is a no-op when flow_bcs isa VelocityBoundaryConditions (velocity is already the primary variable).

JustRelax.JustRelax2D.flow_bcs! Method

julia

flow_bcs!(stokes, bcs::VelocityBoundaryConditions)
flow_bcs!(stokes, bcs::DisplacementBoundaryConditions)
flow_bcs!(bcs, Vx, Vy[, Vz])

Apply no-slip, free-slip, and periodic flow boundary conditions to staggered velocity or displacement arrays. The array form accepts the boundary-condition object first; the stokes form accepts it second. Boundary updates are executed through ParallelStencil kernels on the selected backend.

The three conditions are applied in the order no-slip, free-slip, periodic, so a face carrying more than one of them would keep only the last. The constructors reject such combinations. Faces where all three are false are left untouched. Periodic conditions match normal components at paired boundary planes and copy opposite interior values into tangential ghost planes.

JustRelax.JustRelax2D.fn_ratio Method

julia

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.getindex_NamedTuple Method

julia

getindex_NamedTuple(args::NamedTuple, [sz_min::NTuple], I...)

Sample every array in args at I, offsetting entries larger than sz_min by one to skip their ghost nodes. sz_min is the size of the grid I indexes over; pass it whenever it is known, as inferring it from args reads the ghost node as a cell center when args carries no cell-sized entry.

JustRelax.JustRelax2D.heatdiffusion_PT! Method

julia

heatdiffusion_PT!(thermal, args...; kwargs...)

Dispatch pseudo-transient thermal diffusion to the backend associated with thermal.

See the _heatdiffusion_PT! methods below for the supported argument groups: constant K and ρCp fields, or rheology-driven properties with optional phase ratios and Stokes fields for adiabatic heating.

JustRelax.JustRelax2D.interp_Vx_on_Vy! Method

julia

interp_Vx_on_Vy!(Vx_on_Vy, Vx)

Interpolates the values of Vx onto the grid points of Vy.

Arguments

Vx_on_Vy::AbstractArray: Vx at Vy grid points.
Vx::AbstractArray: Vx at its staggered grid points.

JustRelax.JustRelax2D.isvalid_c Method

julia

isvalid_c(ϕ::JustRelax.RockRatio, inds...)

Check whether the 3D pressure degree of freedom is connected to liquid. As in 2D, the row is retained when the centre carries liquid and all six adjacent velocity faces are active.

JustRelax.JustRelax2D.isvalid_c Method

julia

isvalid_c(ϕ::JustRelax.RockRatio, inds...)

Check whether the 2D pressure degree of freedom is connected to liquid.

The cell-centred pressure row is retained when ϕ.center[i,j] carries liquid and all four adjacent velocity faces are active:

julia

             Vy[i, j+1]
                   o
                   |
    Vx[i, j]  o--- p[i,j] ---o  Vx[i+1, j]
                   |
                   o
             Vy[i, j]

This is the local null-space elimination of the matrix-free reduced system. ∇V is built from these four faces, and the pressure of the cell acts back on them through the weighted gradient, so a row that keeps an inactive face is a row whose divergence no free velocity can relieve: the Powell-Hestenes penalty then drives its pressure without bound.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax2D.isvalid_v Method

julia

isvalid_v(ϕ::JustRelax.RockRatio, inds...)

Check if ϕ.vertex[inds...] is a not a nullspace in 3D.

JustRelax.JustRelax2D.isvalid_v Method

julia

isvalid_v(ϕ::JustRelax.RockRatio, inds...)

Check if ϕ.vertex[inds...] is a not a nullspace 2D.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax2D.isvalid_velocity Method

julia

isvalid_velocity(ϕ::JustRelax.RockRatio, inds...)

Check if the velocity components at the given indices are not nullspaces in 3D.

JustRelax.JustRelax2D.isvalid_velocity Method

julia

isvalid_velocity(ϕ::JustRelax.RockRatio, inds...)

Check if the velocity components at the given indices are not nullspaces in 2D.

JustRelax.JustRelax2D.isvalid_vx Method

julia

isvalid_vx(ϕ::JustRelax.RockRatio, inds...)

Check whether the x-velocity degree of freedom at inds is not a nullspace.

The face is retained when its own control volume carries liquid, ϕ.Vx > 0. A Vx/Vy control volume straddles two cell halves, so it runs dry before the cells it separates do; the face is then a rigid lid on the cut cell behind it. That is what pairs with isvalid_c: a face with no liquid contributes to the divergence of both cells it separates without being able to relieve it, so the two rules have to draw the boundary in the same place.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax2D.isvalid_vy Method

julia

isvalid_vy(ϕ::JustRelax.RockRatio, inds...)

Check whether the y-velocity degree of freedom at inds is not a nullspace. Retained when ϕ.Vy > 0; see isvalid_vx.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax2D.isvalid_vz Method

julia

isvalid_vz(ϕ::JustRelax.RockRatio, inds...)

Check whether the z-velocity degree of freedom at inds is not a nullspace. Retained when ϕ.Vz > 0; see isvalid_vx.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax2D.isvalid_xy Method

julia

isvalid_xy(ϕ, inds...)

Check if the xy shear component at the given indices is not a nullspace.

JustRelax.JustRelax2D.isvalid_xz Method

julia

isvalid_xz(ϕ, inds...)

Check if the xz shear component at the given indices is not a nullspace.

JustRelax.JustRelax2D.isvalid_yz Method

julia

isvalid_yz(ϕ, inds...)

Check if the yz shear component at the given indices is not a nullspace.

JustRelax.JustRelax2D.maximum_mpi Method

julia

maximum_mpi(A)

Compute the maximum value of array A across all MPI processes.

JustRelax.JustRelax2D.mean_mpi Method

julia

mean_mpi(A)

Compute the mean of array A across all MPI processes.

JustRelax.JustRelax2D.minimum_mpi Method

julia

minimum_mpi(A)

Compute the minimum value of array A across all MPI processes.

JustRelax.JustRelax2D.momentum_rows Method

julia

momentum_rows(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool}, d)

Shape of the momentum residual of direction d on a grid of ni cells.

A non-periodic direction has ni[d] - 1 interior faces to solve, the two boundary faces being prescribed. A periodic direction has ni[d]: its two boundary faces are the same plane, so they form one extra unknown, stored as the last row of the residual and mapped to the upper face.

JustRelax.JustRelax2D.multi_copy! Method

julia

multi_copy!(dst::NTuple{N, T}, src::NTuple{N, T}) where {N, T}

Copy data from the tuple of arrays src into the tuple of arrays dst in parallel.

JustRelax.JustRelax2D.norm_mpi Method

julia

norm_mpi(A)

Compute the L2 norm of array A across all MPI processes.

JustRelax.JustRelax2D.pureshear_bc! Method

julia

pureshear_bc!(stokes, xci, xvi, εbg)
pureshear_bc!(stokes, xci, xvi, εbg, backend)

Initialize a pure-shear background velocity field on the staggered grids. xci contains cell-center coordinates and xvi contains velocity-grid coordinates; each component is built from the vertex coordinates of its own direction. In 2D the kernels set Vx = εbg*x and Vy = -εbg*y. In 3D they set Vx = εbg*x, Vy = εbg*y, and Vz = -εbg*z. Ghost layers are left untouched so that subsequent flow boundary-condition and halo updates can set them consistently.

All field updates are performed by ParallelStencil kernels on the backend of stokes. The five-argument form remains available for compatibility; its backend argument is redundant, as the backend is inferred from stokes.

JustRelax.JustRelax2D.rotate_stress! Method

julia

rotate_stress!(τ_particles::StressParticles, stokes, particles, dt)

Interpolate the current deviatoric stress stokes.τ and vorticity stokes.ω onto the particles and rotate the particle stress over dt. stokes.ω must hold the vorticity of the current velocity field. Use stress2grid! afterwards to map the rotated stress back onto stokes.τ_o.

JustRelax.JustRelax2D.rotate_stress_particles! Method

julia

rotate_stress_particles!(τ::NTuple, ω::NTuple, particles::Particles, dt; method = :matrix)

Rotate the deviatoric stress carried by each active particle over dt with the local vorticity, using GeoParams' elastic stress rotation. τ holds the stress components and ω the vorticity components, as particle cell arrays. method is accepted for call-site compatibility and does not select an algorithm.

JustRelax.JustRelax2D.shear2center! Method

julia

shear2center!(A::SymmetricTensor)

Interpolate the shear components of A onto the cell centers, in place.

JustRelax.JustRelax2D.simpleshear_bc! Method

julia

simpleshear_bc!(stokes, xci, xvi, γbg)
simpleshear_bc!(stokes, xci, xvi, γbg, backend)

Initialize an xy simple-shear background velocity field on the staggered grids. The imposed field is Vx = γbg * y; the other velocity components are set to zero. xci contains cell-center coordinates and xvi contains velocity-grid coordinates. Ghost layers are left untouched.

All field updates are performed by ParallelStencil kernels on the backend of stokes. The five-argument form remains available for compatibility; its backend argument is redundant, as the backend is inferred from stokes.

JustRelax.JustRelax2D.solve! Method

julia

solve!(stokes::StokesArrays, pt_stokes, grid, flow_bcs, ρg, phase_ratios, rheology, args, dt, igg; kwargs...)

Solve the 2D viscoelastoplastic Stokes equations to pseudo-transient convergence, updating stokes (velocity, pressure, stress, viscosity) in place for one physical time step dt. This is the general, multi-phase call form used by most models; dispatch on the type of the sixth argument also accepts, for simpler/benchmark setups:

a single-phase rheology::GeoParams.MaterialParams in place of phase_ratios/rheology (drop the phase_ratios argument), or
constant K (bulk modulus) or K, G (bulk and shear modulus) fields in place of phase_ratios/rheology/args, for linear (visco)elastic problems with no material rheology.

grid may also be replaced by the grid spacing di alone (a NTuple/NamedTuple).

Arguments

stokes: solver state allocated with StokesArrays.
pt_stokes: pseudo-transient coefficients, e.g. from PTStokesCoeffs.
grid: the model Geometry.
flow_bcs: velocity/displacement boundary conditions.
ρg: buoyancy forcing (ρgx, ρgy).
phase_ratios: per-cell/per-node phase fractions (a JustPIC.PhaseRatios).
rheology: one GeoParams.MaterialParams per phase.
args: auxiliary fields (e.g. temperature T, pressure P) used by the constitutive updates.
dt: physical time step.
igg: the distributed-grid context (IGG).

Keyword arguments

iterMax = 50e3, iterMin = 100: min/max pseudo-transient iterations.
nout = 500: check convergence every nout iterations.
viscosity_cutoff = (-Inf, Inf): clamp bounds for the effective viscosity.
viscosity_relaxation = 1e-2: relaxation factor for nonlinear viscosity updates.
λ_relaxation = 0.2: relaxation factor for the plastic multiplier.
strain_increment = false: accumulate strain increments instead of overwriting.
free_surface = false: enable free-surface stabilization.
b_width = (4, 4, 0): halo width used when overlapping communication and computation.
verbose = true: print convergence progress.

Dispatches on the CPU/CUDA/AMDGPU backend selected by stokes.

JustRelax.JustRelax2D.solve_DYREL! Method

julia

solve_DYREL!(
    stokes, ρg, dyrel, flow_bcs, phase_ratios, rheology, args, grid, dt, igg;
    kwargs...,
)

Solve the Stokes system with the self-tuned dynamic relaxation (DYREL) method.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
ρg: buoyancy forces arrays.
dyrel: DYREL-specific parameters and fields.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
grid: Geometry object carrying grid spacing and staggered-grid coordinates. A legacy 2D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
λ_relaxation_DR: Relaxation factor for dynamic relaxation. Default: 1.
λ_relaxation_PH: Relaxation factor for Powell-Hestenes iterations. Default: 1.
pressure_relaxation: Relaxation factor for the Powell-Hestenes pressure update. Default: 1.
iterMax_PH: Maximum number of Powell-Hestenes passes. Default: 1.0e3.
iterMax_DR: Maximum number of iterations for each dynamic-relaxation solve. Default: 50.0e3.
iterMax: Compatibility alias for iterMax_DR; used when iterMax_DR is not given.
total_iterMax: Maximum number of total dynamic-relaxation iterations. Default: 50.0e3.
nout: Output frequency for residuals. Default: 100.
rel_drop: Relative residual drop tolerance. Default: 1.0e-2.
verbose_PH: Print Powell-Hestenes iteration info. Default: true.
verbose_DR: Print Dynamic Relaxation iteration info. Default: true.
linear_viscosity: Whether to use linear viscosity. Default: false.
free_surface: Include the density-gradient free-surface stabilization term. Default: false.
update_material: Recompute viscosity and buoyancy from rheology. Set to false when those fields are prescribed by the caller. Default: true.

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax2D.solve_VariationalDYREL! Method

julia

solve_VariationalDYREL!(stokes, ρg, dyrel, flow_bcs, phase_ratios, ϕ,
    rheology, args, grid, dt, igg; kwargs...)

Solve the 2D variational Stokes problem with DYREL relaxation and the RockRatio volume weights. This is a separate entry point from solve_DYREL!; the latter remains the standard, unweighted DYREL solver.

Center fractions weight pressure and normal stress, vertex fractions weight shear stress, and face fractions weight momentum rows. Rows whose volume fraction vanishes are eliminated rather than solved with air properties.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
ρg: buoyancy forces arrays.
dyrel: DYREL-specific parameters and fields, built with the same ϕ.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
ϕ: JustRelax.RockRatio carrying the cell, vertex and face volume fractions.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
grid: Geometry{2} object carrying grid spacing and staggered-grid coordinates. A legacy 2D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

air_phase: Phase index excluded from material averages; 0 disables the correction. Default: 0.
viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
λ_relaxation_DR: Relaxation factor for dynamic relaxation. Default: 1.
λ_relaxation_PH: Relaxation factor for Powell-Hestenes iterations. Default: 1.
pressure_relaxation: Relaxation factor for the Powell-Hestenes pressure update. Default: 1.
iterMax_PH: Maximum number of Powell-Hestenes passes. Default: 1.0e3.
iterMax_DR: Maximum number of iterations for each dynamic-relaxation solve. Default: 50.0e3.
iterMax: Alias for iterMax_DR; used when iterMax_DR is not given.
total_iterMax: Maximum number of total dynamic-relaxation iterations. Default: 50.0e3.
nout: Output frequency for residuals. Default: 100.
rel_drop: Relative residual drop tolerance. Default: 1.0e-2.
verbose_PH: Print Powell-Hestenes iteration info. Default: true.
verbose_DR: Print Dynamic Relaxation iteration info. Default: true.
linear_viscosity: Whether to use linear viscosity. Default: false.
free_surface: Include the density-gradient free-surface stabilization term. Default: false.

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax2D.solve_VariationalStokes! Method

julia

solve_VariationalStokes!(backend::BackendTrait, stokes::JustRelax.StokesArrays, args...; kwargs)

Stokes solver entry point for variational Stokes solvers. This function dispatches to the appropriate implementation based on the backend provided in the function call.

JustRelax.JustRelax2D.solve_VariationalStokes! Method

julia

solve_VariationalStokes!(stokes::JustRelax.StokesArrays, args...; kwargs...)

Solve the 2D volume-fraction variational Stokes problem with matrix-free pseudo-transient iterations.

ϕ carries liquid weights at pressure cells, stress vertices, and staggered velocity faces. A pressure degree of freedom is retained only when its cell and all four surrounding velocity faces are connected to liquid:

julia

             Vy[i, j+1]
                   o
                   |
    Vx[i, j]  o--- p[i,j] ---o  Vx[i+1, j]
                   |
                   o
             Vy[i, j]

   inactive face => pressure row and disconnected velocity row eliminated

Zero-weight rows are written as zero instead of being solved with air material properties. Positive sliver fractions remain active; their velocity diagonal uses the bounded face mass max(ϕ_face, 0.1).

The free_surface keyword enables the density-gradient correction in the vertical momentum row. In this solver it is included implicitly in the local face diagonal, so the physical timestep does not create an explicit feedback instability.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
pt_stokes: Pseudo-transient coefficients, from PTStokesCoeffs.
grid: Geometry{2} object carrying grid spacing and staggered-grid coordinates. A legacy 2D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
ρg: buoyancy forces arrays.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
ϕ: JustRelax.RockRatio carrying the cell, vertex and face volume fractions.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

air_phase: Phase index excluded from material averages; 0 disables the correction. Default: 0.
viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
λ_relaxation: Relaxation factor for the plastic multiplier. Default: 0.2.
strain_increment: Solve for displacement increments alongside velocity. Default: false.
iterMax: Maximum number of pseudo-transient iterations. Default: 50.0e3.
iterMin: Minimum number of pseudo-transient iterations. Default: 1.0e2.
nout: Output frequency for residuals. Default: 500.
verbose: Print iteration info. Default: true.
free_surface: Include the density-gradient free-surface stabilization term. Default: false.
b_width: Halo width used to overlap communication with computation. Default: (4, 4, 0).

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax2D.stress2grid! Method

julia

stress2grid!(stokes, τ_particles::StressParticles, particles)

Interpolate the particle stress in τ_particles back onto the old-stress fields stokes.τ_o: normal components onto the cell centers, and shear components onto the vertices in 2D or onto the cell centers and edges in 3D, matching where the stress kernels read them from. Counterpart of rotate_stress!, and the step that hands the rotated stress to the next Stokes solve.

JustRelax.JustRelax2D.subgrid_characteristic_time! Method

julia

subgrid_characteristic_time!(subgrid_arrays, particles, dt₀, phases, rheology, thermal::ThermalArrays, stokes::StokesArrays[, di])

Compute, per cell, the characteristic thermal diffusion timescale dt₀ = ρCp / (2 K Σ dxi⁻²) used for JustPIC's subgrid-diffusion correction of particle temperature, evaluating rheology's density/heat-capacity/conductivity at the local phase (from phases, either a JustPIC.PhaseRatios or an integer phase-id array) and temperature/pressure. dt₀ must be ghosted, with size size(stokes.P) .+ 2; values are written to its interior region.

JustRelax.JustRelax2D.sum_mpi Method

julia

sum_mpi(A)

Compute the sum of array A across all MPI processes.

JustRelax.JustRelax2D.take Method

julia

take(fldr::String)

Create folder fldr if it does not exist.

JustRelax.JustRelax2D.tensor_invariant! Method

julia

tensor_invariant!(A::JustRelax.SymmetricTensor)

Compute the tensor invariant of the given symmetric tensor A.

Arguments

A::JustRelax.SymmetricTensor: The input symmetric tensor.

JustRelax.JustRelax2D.thermal_bcs! Method

julia

thermal_bcs!(thermal, bcs::TemperatureBoundaryConditions)
thermal_bcs!(T::AbstractArray, bcs::TemperatureBoundaryConditions)

Apply thermal ghost-cell boundary conditions to a temperature field.

thermal_bcs! applies the scalar temperature conditions stored in bcs:

constant_value faces are applied first using Tghost = 2 * value - Tinterior.
no_flux faces are applied next by copying the adjacent interior temperature.
periodic faces are applied last by copying the opposite interior temperature into the ghost layer.

Faces set to false are ignored. Periodic faces must be paired by direction and cannot also carry another thermal condition. Prescribed constant_flux values are not applied here; they are consumed by the pseudo-transient heat-diffusion compute_flux! kernels.

JustRelax.JustRelax2D.unpack Method

julia

unpack(x::T)

Generated function to unpack the fields of the struct x into a tuple.

JustRelax.JustRelax2D.update_T Method

julia

update_T(::Nothing, b_width, thermal, ρCp, pt_thermal, dirichlet, _dt, _di, ni)
update_T(::Nothing, b_width, thermal, rheology, phase, pt_thermal, dirichlet, _dt, _di, ni, args)

Launch the pseudo-transient temperature update kernel over the active thermal domain.

These wrappers select the appropriate kernel overload depending on whether the solver works with precomputed ρCp fields or rheology-derived properties.

JustRelax.JustRelax2D.update_dτV_α_β! Method

julia

update_dτV_α_β!(dτV, βV, αV, cV, λmaxV, CFL_v)

Generic N-dimensional version (2D and 3D) for updating pseudo-time step and acceleration parameters.

Computes the pseudo-time step dτV from the maximum eigenvalue estimate λmaxV and CFL number, then updates the damping parameters βV and acceleration parameters αV accordingly.

Arguments

dτV::NTuple{N, AbstractArray{T, N}}: Tuple of pseudo-time steps for each velocity component
βV::NTuple{N, AbstractArray{T, N}}: Tuple of damping parameters β for each velocity component
αV::NTuple{N, AbstractArray{T, N}}: Tuple of acceleration parameters α for each velocity component
cV::NTuple{N, AbstractArray{T, N}}: Tuple of preconditioner diagonal entries for each velocity component
λmaxV::NTuple{N, AbstractArray{T, N}}: Tuple of maximum eigenvalue estimates for each velocity component
CFL_v::Real: CFL number for velocity

JustRelax.JustRelax2D.update_phase_ratios_2D! Method

julia

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:

julia

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

julia

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:

julia

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

julia

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.

JustRelax.JustRelax2D.update_pt_thermal_arrays! Method

julia

update_pt_thermal_arrays!(pt_thermal, phase_ratios, rheology, args, _dt)

Recompute the pseudo-transient thermal coefficient arrays stored in pt_thermal from phase-weighted material properties.

This helper is used by the pseudo-transient thermal solver when the local phase mixture changes over time.

JustRelax.JustRelax2D.update_rock_ratio! Method

julia

update_rock_ratio!(ϕ::JustRelax.RockRatio, phase_ratios, air_phase)

Update the rock ratio ϕ based on the provided phase_ratios and air_phase.

Arguments

ϕ::JustRelax.RockRatio: The rock ratio object to be updated.
phase_ratios: The ratios of different phases present.
air_phase: The phase representing air.

JustRelax.JustRelax2D.update_rock_ratio! Method

julia

update_rock_ratio!(ϕ::JustRelax.RockRatio, phase_ratios, air_phase)

Update the rock ratio ϕ for a 3D grid based on the provided phase_ratios and air_phase.

Arguments

ϕ::JustRelax.RockRatio: The rock ratio object to be updated.
phase_ratios: The ratios of different phases present.
air_phase: The phase representing air.

JustRelax.JustRelax2D.update_rock_ratio_cv! Method

julia

update_rock_ratio_cv!(ϕ, ratio_center, ratio_vertex, air_phase)

Update the rock ratio for both center and vertex values based on the provided ratio_center, ratio_vertex, and air_phase.

JustRelax.JustRelax2D.update_thermal_coeffs! Method

julia

update_thermal_coeffs!(pt_thermal, rheology, phase_ratios, args, dt)
update_thermal_coeffs!(pt_thermal, rheology, args, dt)

Refresh the pseudo-transient coefficient arrays stored in pt_thermal.

Use this when conductivity or volumetric heat capacity changes during the solve, for example because temperature, pressure, or phase ratios evolve. The update is performed in place on pt_thermal.θr_dτ and pt_thermal.dτ_ρ.

JustRelax.JustRelax2D.update_α_β! Method

julia

update_α_β!(βV, αV, dτV, cV)

Generic N-dimensional version (2D and 3D) of the acceleration parameters update.

Updates the damping parameters βV and αV for each velocity component based on the pseudo-time step dτV and the preconditioner diagonal cV.

Arguments

βV::NTuple{N, AbstractArray{T, N}}: Tuple of damping parameters β for each velocity component
αV::NTuple{N, AbstractArray{T, N}}: Tuple of acceleration parameters α for each velocity component
dτV::NTuple{N, AbstractArray{T, N}}: Tuple of pseudo-time steps for each velocity component
cV::NTuple{N, AbstractArray{T, N}}: Tuple of preconditioner diagonal entries for each velocity component

JustRelax.JustRelax2D.velocity2center! Method

julia

velocity2center!(Vx_c, Vy_c, Vx, Vy)

In-place interpolation of the velocity field Vx, Vy from a staggered grid with ghost nodes onto the pre-allocated Vx_c, Vy_c 2D arrays located at the cell centers.

JustRelax.JustRelax2D.velocity2center! Method

julia

velocity2center!(Vx_c, Vy_c, Vz_c, Vx, Vy, Vz)

In-place interpolation of the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the pre-allocated Vx_c, Vy_c, Vz_c 3D arrays located at the cell centers.

JustRelax.JustRelax2D.velocity2displacement! Method

julia

velocity2displacement!(stokes::StokesArrays, dt)

Set stokes.U (displacement) to stokes.V * dt (velocity times the time step), in place.

JustRelax.JustRelax2D.velocity2vertex! Method

julia

velocity2vertex!(Vx_v, Vy_v, Vx, Vy)

In-place interpolation of the velocity field Vx, Vy from a staggered grid with ghost nodes onto the pre-allocated Vx_v, Vy_v 2D arrays located at the grid vertices.

JustRelax.JustRelax2D.velocity2vertex! Method

julia

velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz)

In-place interpolation of the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the pre-allocated Vx_d, Vy_d, Vz_d 3D arrays located at the grid vertices.

JustRelax.JustRelax2D.velocity2vertex Method

julia

velocity2vertex(Vx, Vy, Vz)

Interpolate the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the grid vertices.

JustRelax.JustRelax2D.vertex2center! Method

julia

vertex2center!(center, vertex; ghost_x = false, ghost_y = false, ghost_z = false)

Interpolates the values at the vertex onto center points.

center may carry a ring of ghost nodes that vertex does not; each ghost_* keyword shifts the write index by one along that dimension, so the interpolated block lands in the interior of center instead of overwriting its ghost nodes. ghost_z is ignored in 2D.

JustRelax.JustRelax2D.viscosity_phase_ratio Method

julia

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.@add Macro

julia

@add(I, args...)

Add I to the scalars in args

JustRelax.JustRelax2D.@allocate Macro

julia

@allocate(ni...)

Convenience macro to allocate a PTArray of size ni... with undef values.

JustRelax.JustRelax2D.@copy Macro

julia

copy(B, A)

convenience macro to copy data from the array A into array B

JustRelax.JustRelax2D.@displacement Macro

julia

@displacement(U)

Unpacks the displacement arrays U from the StokesArrays A.

JustRelax.JustRelax2D.@idx Macro

julia

@idx(args...)

Make a linear range from 1 to args[i], with i ∈ [1, ..., n]

JustRelax.JustRelax2D.@normal Macro

julia

@normal(A)

Unpacks the normal components of the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@plastic_strain Macro

julia

@plastic_strain(A)

Unpacks the plastic strain rate tensor ε_pl from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@qT Macro

julia

@qT(V)

Unpacks the flux arrays qT_i from the ThermalArrays A.

JustRelax.JustRelax2D.@qT2 Macro

julia

@qT2(V)

Unpacks the flux arrays qT2_i from the ThermalArrays A.

JustRelax.JustRelax2D.@residuals Macro

julia

@residuals(A)

Unpacks the momentum residuals from A.

JustRelax.JustRelax2D.@shear Macro

julia

@shear(A)

Unpacks the shear components of the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@shear_center Macro

julia

@shear_center(A)

Unpacks the shear components of the symmetric tensor A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@strain Macro

julia

@strain(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@strain_center Macro

julia

@strain_center(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@strain_increment Macro

julia

@strain_increment(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@stress Macro

julia

@stress(A)

Unpacks the deviatoric stress tensor τ from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@stress_center Macro

julia

@stress_center(A)

Unpacks the deviatoric stress tensor τ from the StokesArrays A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@tensor Macro

julia

@tensor(A)

Unpacks the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@tensor_center Macro

julia

@tensor_center(A)

Unpacks the symmetric tensor A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@tensor_vertex Macro

julia

tensor_vertex(A)

Unpacks the symmetric tensor A, where its components are defined in the vertices of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax2D.@tuple Macro

julia

@tuple(A)

Convenience maktro to unpack the fields of the struct A into a tuple. Works with Velocity and SymmetricTensor structs.

JustRelax.JustRelax2D.@unpack Macro

julia

@unpack(x)

Convenience macro to unpack the fields of the struct x into a tuple.

JustRelax.JustRelax2D.@velocity Macro

julia

@velocity(V)

Unpacks the velocity arrays V from the StokesArrays A.

JustRelax.JustRelax3D Module

julia

JustRelax.JustRelax3D

Three-dimensional solvers, kernels, and constructors, on the CPU backend.

The submodule is loaded with using JustRelax.JustRelax3D, and its ParallelStencil environment is initialized for three dimensions when it loads. Loading CUDA.jl or AMDGPU.jl before it adds the matching GPU methods through a package extension; the entry points and their signatures stay the same, and the backend is chosen by the array type of the containers passed in. See JustRelax.JustRelax2D for the two-dimensional counterpart.

JustRelax.JustRelax3D.DYREL Function

julia

DYREL(ni::NTuple{N, Integer}[, periodic]; ϵ=1e-6, ϵ_vel=1e-6, CFL=0.99, c_fact=0.5, γfact=20.0) where N

Creates a new DYREL struct with fields initialized to zero.

Arguments

ni: Tuple containing the grid dimensions (nx, ny) for 2D or (nx, ny, nz) for 3D.
periodic: N-tuple marking the periodic directions, which each carry one extra momentum row (see momentum_rows). Defaults to all-false. The StokesArrays method below reads it off the containers instead, so the two cannot disagree.

Keyword arguments

ϵ: General convergence tolerance. Default: 1.0e-6.
ϵ_vel: Velocity convergence tolerance. Default: 1.0e-6.
CFL: Courant-Friedrichs-Lewy number. Default: 0.99.
c_fact: Damping scaling factor. Default: 0.5.
γfact: Penalty scaling factor. Default: 20.0.

JustRelax.JustRelax3D.DYREL! Method

julia

DYREL!(dyrel::JustRelax.DYREL, stokes::JustRelax.StokesArrays, rheology, phase_ratios, di, dt; CFL=dyrel.CFL, γfact=dyrel.γfact)

Updates the fields of the DYREL struct in-place for the current time step.

This function recomputes:

Bulk viscosity and penalty parameter γ_eff.
Gershgorin estimates for eigenvalues and preconditioners.
Damping coefficients.

Arguments

dyrel: JustRelax.DYREL struct to modify.
stokes: JustRelax.StokesArrays containing current simulation state.
rheology, phase_ratios: Material properties.
di: Grid spacing.
dt: Current time step.
CFL: Courant number (default: the value stored in dyrel).
γfact: Penalty factor (default: the value stored in dyrel).

Returns nothing.

JustRelax.JustRelax3D.DYREL Method

julia

DYREL(stokes, rheology, phase_ratios, di, dt; ϵ=1e-6, ϵ_vel=1e-6, CFL=0.99, c_fact=0.5, γfact=20.0)

Constructs and initializes a DYREL object based on existing Stokes fields.

This function:

Allocates zero-initialized arrays using grid dimensions from stokes.
Computes initial bulk viscosity and penalty parameters.
Computes Gershgorin estimates for eigenvalues and preconditioners.
Updates damping coefficients.

Arguments

stokes: JustRelax.StokesArrays struct.
rheology: Material properties.
phase_ratios: Phase fraction information.
di: Grid spacing tuple.
dt: Time step.

Keyword arguments

ϵ: General convergence tolerance. Default: 1.0e-6.
ϵ_vel: Velocity convergence tolerance. Default: 1.0e-6.
CFL: Courant-Friedrichs-Lewy number. Default: 0.99.
c_fact: Damping scaling factor. Default: 0.5.
γfact: Factor for the penalty parameter calculation. Default: 20.0.

JustRelax.JustRelax3D.Displacement Method

julia

Displacement(nx::Integer, ny::Integer, nz::Integer)

Create the displacement arrays for the Stokes solver in 3D.

Fields

Ux: Displacement in x direction at their staggered location
Uy: Displacement in y direction at their staggered location
Uz: Displacement in z direction at their staggered location

JustRelax.JustRelax3D.Displacement Method

julia

Displacement(nx::Integer, ny::Integer)

Create the displacement arrays for the Stokes solver in 2D.

Fields

Ux: Displacement in x direction at their staggered location
Uy: Displacement in y direction at their staggered location

JustRelax.JustRelax3D.PTThermalCoeffs Method

julia

PTThermalCoeffs(rheology, phase_ratios, args, dt, ni, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from a multi-phase rheology.

The coefficient arrays are evaluated on the thermal grid of size ni using the phase ratios in phase_ratios.center together with the thermodynamic state stored in args.

JustRelax.JustRelax3D.PTThermalCoeffs Method

julia

PTThermalCoeffs(K, ρCp, dt, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from conductivity and volumetric heat-capacity arrays.

di and li are the grid spacing and domain-length tuples used to estimate the local pseudo-time step and relaxation factor stored in the returned JustRelax.PTThermalCoeffs.

JustRelax.JustRelax3D.PTThermalCoeffs Method

julia

PTThermalCoeffs(rheology, args, dt, ni, di, li; ϵ = 1.0e-8, CFL = 0.9 / √3)

Construct pseudo-transient thermal coefficients from a single rheology state.

This method evaluates conductivity and ρCp from rheology and the cellwise fields in args, then stores the resulting relaxation arrays in a JustRelax.PTThermalCoeffs object.

JustRelax.JustRelax3D.PrincipalStress Method

julia

PrincipalStress(ni::NTuple{N, Integer}) where {N}

Create the principal stress arrays for the Stokes solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

Fields

σ1: First principal stress
σ2: Second principal stress
σ3: Third principal stress (only in 3D). In 2D it is a placeholder array of size (2, 1, 1).

JustRelax.JustRelax3D.Residual Method

julia

Residual(nx::Integer, ny::Integer[, nz::Integer])
Residual(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool})

Create the residual arrays for the Stokes solver.

Fields

Rx, Ry[, Rz]: Residuals for the momentum equations
RP: Residual for the continuity equation

periodic marks the directions whose two boundary faces are the same plane; each of those gains one momentum row (see momentum_rows). It defaults to all-false.

JustRelax.JustRelax3D.RockRatio Method

julia

RockRatio(nx, ny, nz)

Create a RockRatio object for a 3D grid with dimensions nx x ny x nz on a staggered grid.

JustRelax.JustRelax3D.RockRatio Method

julia

RockRatio(nx, ny)

Create a RockRatio object for a 2D grid with dimensions nx x ny on a staggered grid.

JustRelax.JustRelax3D.RockRatio Method

julia

RockRatio(backend, ni)
RockRatio(backend, ni...)

Allocate a RockRatio on backend for a staggered grid of ni cells, with every volume fraction initialized to zero.

backend is the same backend type passed to the other allocators (CPUBackend, CUDABackend, AMDGPUBackend), and ni is (nx, ny) in 2D or (nx, ny, nz) in 3D. Fill the fractions with update_rock_ratio!, or with JustPIC.compute_rock_fraction! when the liquid domain is bounded by a marker chain.

JustRelax.JustRelax3D.StokesArrays Method

julia

StokesArrays(ni::NTuple{N,Integer}) where {N}
StokesArrays(ni::NTuple{N,Integer}, bcs::AbstractFlowBoundaryConditions)
StokesArrays(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool})

Create the Stokes arrays object in 2D or 3D.

Passing the flow boundary conditions sizes the momentum residuals for the periodic directions they declare, which is what gives a periodic seam a momentum row; without them every direction is taken to be non-periodic. The solvers check the two against each other, so the boundary conditions have to be built first.

Fields

P: Pressure field
P0: Previous pressure field
∇V: Velocity gradient
V: Velocity fields
Q: Volumetric source/sink term e.g. ΔV/V_tot [m³/m³]
U: Displacement fields
ω: Vorticity field
τ: Stress tensors
τ_o: Old stress tensors
ε: Strain rate tensors
ε_pl: Plastic strain rate tensors
EII_pl: Second invariant of the accumulated plastic strain
viscosity: Viscosity fields
R: Residual fields
Δε: Strain increment tensor
∇U: Displacement gradient
λ : plastic multiplier @ centers
λv : plastic multiplier @ vertices
λv_yz, λv_xz, λv_xy : 3D plastic multiplier on shear staggered grids
ΔPψ : pressure correction in dilatant case

JustRelax.JustRelax3D.StressParticles Method

julia

StressParticles(particles::Particles)

Allocate the stress and vorticity cell arrays that follow particles, on the same backend and with the same per-cell capacity. Two normal and one shear component in 2-D, three of each in 3-D.

JustRelax.JustRelax3D.SymmetricTensor Method

julia

SymmetricTensor(nx::Integer, ny::Integer, nz::Integer)

Create the symmetric tensor arrays for the Stokes solver in 3D.

Fields

xx: xx component of the tensor at cell centers
yy: yy component of the tensor at cell centers
zz: zz component of the tensor at cell centers
xx_v: xx component of the tensor at vertices
yy_v: yy component of the tensor at vertices
zz_v: zz component of the tensor at vertices
xy: xy component of the tensor at vertices
yz: yz component of the tensor at vertices
xz: xz component of the tensor at vertices
yz_c: yz component of the tensor at cell centers
xz_c: xz component of the tensor at cell centers
xy_c: xy component of the tensor at cell centers
II: second invariant of the tensor at cell centers

JustRelax.JustRelax3D.SymmetricTensor Method

julia

SymmetricTensor(nx::Integer, ny::Integer)

Create the symmetric tensor arrays for the Stokes solver in 2D.

Fields

xx: xx component of the tensor at cell centers
yy: yy component of the tensor at cell centers
xx_v: xx component of the tensor at vertices
yy_v: yy component of the tensor at vertices
xy: xy component of the tensor at vertices
xy_c: xy component of the tensor at cell centers
II: second invariant of the tensor at cell centers

JustRelax.JustRelax3D.ThermalArrays Method

julia

ThermalArrays(nx::Integer, ny::Integer, nz::Integer)

3D constructor for the thermal arrays for the heat diffusion solver with the extents given by nx, ny and nz.

Fields

T: Temperature at cell centers with one ghost node on every boundary (nx + 2, ny + 2, nz + 2)
Told: Temperature at previous time step at cell centers with ghost nodes (nx + 2, ny + 2, nz + 2)
ΔT: Temperature change at cell centers with ghost nodes (nx + 2, ny + 2, nz + 2)
adiabatic: Adiabatic term α (u ⋅ ∇P) at cell centers (nx, ny, nz)
dT_dt: Time derivative of temperature at cell centers (nx, ny, nz)
qTx: Conductive heat flux in x direction on cell faces (nx + 1, ny, nz)
qTy: Conductive heat flux in y direction on cell faces (nx, ny + 1, nz)
qTz: Conductive heat flux in z direction on cell faces (nx, ny, nz + 1)
qTx2: Conductive heat flux in x direction on cell faces for second order scheme (nx + 1, ny, nz)
qTy2: Conductive heat flux in y direction on cell faces for second order scheme (nx, ny + 1, nz)
qTz2: Conductive heat flux in z direction on cell faces for second order scheme (nx, ny, nz + 1)
H: Source terms at cell centers (nx, ny, nz)
shear_heating: Shear heating terms at cell centers (nx, ny, nz)
ResT: Residual of the temperature equation at cell centers (nx, ny, nz)

JustRelax.JustRelax3D.ThermalArrays Method

julia

ThermalArrays(nx::Integer, ny::Integer)

2D constructor for the thermal arrays for the heat diffusion solver with the extents given by nx and ny.

Fields

T: Temperature at cell centers with one ghost node on every boundary (nx + 2, ny + 2)
Told: Temperature at previous time step at cell centers with ghost nodes (nx + 2, ny + 2)
ΔT: Temperature change at cell centers with ghost nodes (nx + 2, ny + 2)
adiabatic: Adiabatic term α (u ⋅ ∇P) at cell centers (nx, ny)
dT_dt: Time derivative of temperature at cell centers (nx, ny)
qTx: Conductive heat flux in x direction on cell faces (nx + 1, ny)
qTy: Conductive heat flux in y direction on cell faces (nx, ny + 1)
qTx2: Conductive heat flux in x direction on cell faces for second order scheme (nx + 1, ny)
qTy2: Conductive heat flux in y direction on cell faces for second order scheme (nx, ny + 1)
H: Source terms at cell centers (nx, ny)
shear_heating: Shear heating terms at cell centers (nx, ny)
ResT: Residual of the temperature equation at cell centers (nx, ny)

JustRelax.JustRelax3D.ThermalArrays Method

julia

ThermalArrays(ni::NTuple{N, Integer}) where {N}

Create the thermal arrays for the heat diffusion solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

JustRelax.JustRelax3D.ThermalArrays Method

julia

ThermalArrays(::Backend, ni::NTuple{N, Integer}) where {N}

Internal entry point function for the ThermalArrays constructor. This allows for dispatching on the backend type and then calling the main constructor with the dimensions.

JustRelax.JustRelax3D.Velocity Method

julia

Velocity(nx::Integer, ny::Integer, nz::Integer)

Create the velocity arrays for the Stokes solver in 3D.

Fields

Vx: Velocity in x direction (nx + 1, ny + 2, nz + 2)
Vy: Velocity in y direction (nx + 2, ny + 1, nz + 2)
Vz: Velocity in z direction (nx + 2, ny + 2, nz + 1)

JustRelax.JustRelax3D.Velocity Method

julia

Velocity(nx::Integer, ny::Integer)

Create the velocity arrays for the Stokes solver in 2D.

Fields

Vx: Velocity in x direction (nx + 1, ny + 2)
Vy: Velocity in y direction (nx + 2, ny + 1)

JustRelax.JustRelax3D.Viscosity Method

julia

Viscosity(ni::NTuple{N, Integer}) where {N}

Create the viscosity arrays for the Stokes solver in 2D or 3D with the extents given by ni (nx x ny or `nx x ny x nz``).

Fields

η: Viscosity at cell centers
ηv: Viscosity at vertices
η_vep: Viscosity for visco-elastic-plastic rheology
ητ: Pseudo-transient viscosity for stress update

JustRelax.JustRelax3D.Vorticity Method

julia

Vorticity(nx::Integer, ny::Integer, nz::Integer)

Create the vorticity arrays for the Stokes solver in 3D.

Fields

yz: Vorticity component yz at their staggered location
xz: Vorticity component xz at their staggered location
xy: Vorticity component xy at their staggered location
yz_c: Vorticity component yz at the cell centers
xz_c: Vorticity component xz at the cell centers
xy_c: Vorticity component xy at the cell centers

JustRelax.JustRelax3D.Vorticity Method

julia

Vorticity(nx::Integer, ny::Integer)

Create the vorticity arrays for the Stokes solver in 2D.

Fields

xy: Vorticity component xy at vertices

JustRelax.JustRelax3D.WENO5 Method

julia

WENO5(backend, method::Val, ni::NTuple)
WENO5(method::Val, ni::NTuple)

Allocate a WENO5 advection scheme state for a grid of size ni. method selects the smoothness-indicator variant: Val(1) for Jiang-Shu, Val(2) for the Z-variant.

JustRelax.JustRelax3D.WENO_advection! Method

julia

WENO_advection!(u, Vxi, weno, di, ni, dt)

Perform the advection step of the Weighted Essentially Non-Oscillatory (WENO) scheme for the solution of hyperbolic partial differential equations.

Arguments

u: field to be advected.
Vxi: velocity field.
weno: structure containing the WENO scheme parameters and temporary variables.
di: grid spacing.
ni: number of grid points.
dt: time step.

Description

The function approximates the advected fluxes using the WENO scheme and use a strong-stability preserving (SSP) Runge-Kutta method of order 3 for the time integration.

JustRelax.JustRelax3D._bc_coordinate Method

julia

_bc_coordinate(A, x)

Materialize coordinate vector x on the same device as array A, so that background-field kernels can index it. collect is applied first because a lazy range cannot be transferred to a device array directly.

JustRelax.JustRelax3D._clamped_index Method

julia

clamped_indices(ni, [periodic,] I...)

Cell indices of the stencil around vertex I on a grid of ni cells, kept in range.

Outside a periodic direction an index that falls off the grid is clamped onto its in-range neighbour, so a vertex-centred average degenerates to the one-sided average of the cells that do exist. A direction listed in periodic has no edge: the cell on the far side of the seam is a real neighbour, so the index wraps onto it and both copies of the seam plane see the same stencil. Omitting periodic clamps every direction.

JustRelax.JustRelax3D._heatdiffusion_PT! Method

julia

_heatdiffusion_PT!(thermal, pt_thermal, thermal_bc, K, ρCp, dt, grid;
    igg, b_width, iterMax, nout, verbose)

Solve the heat equation with pseudo-transient iterations using precomputed, cell-centered material properties.

K is the thermal conductivity field and ρCp is the volumetric heat-capacity field on the thermal grid. pt_thermal supplies the pseudo-transient coefficients, thermal_bc applies the temperature boundary conditions after each update, and grid provides metric terms and halo layout.

Returns a named tuple containing the iteration numbers and residual norms sampled every nout iterations.

JustRelax.JustRelax3D._heatdiffusion_PT! Method

julia

_heatdiffusion_PT!(thermal, pt_thermal, thermal_bc, rheology, args, dt, grid;
    igg, phase, stokes, b_width, iterMax, nout, verbose)

Solve the heat equation with pseudo-transient iterations using thermal properties derived from rheology.

args is a named tuple of thermodynamic fields, typically including T and P. Entries sized like the thermal cell centers are read as-is; larger entries are offset by one to skip their ghost nodes. When phase is provided, pseudo-transient coefficients are recomputed from the local phase ratios each iteration. When stokes is provided, thermal.adiabatic is refreshed before the iteration loop to include the adiabatic heating contribution.

Returns a named tuple containing the sampled iteration counts and residual history.

JustRelax.JustRelax3D._update_rock_ratio! Method

julia

_update_rock_ratio!(ϕ, ratio, air_phase)

Inner kernel of update_rock_ratio that clamps the computed rock ratio to the range [0, 1] for the given ratio and air_phase.

JustRelax.JustRelax3D.accumulate_tensor! Method

julia

accumulate_tensor!(II, A::JustRelax.SymmetricTensor, dt)

Accumulate the second invariant of the symmetric tensor A over a time step: II[I] += dt * A_II[I].

The tensor components are gathered onto the cell centers where II lives, so A may hold its shear components on the vertices. Used to integrate the deviatoric plastic strain rate ε_pl into EII_pl; the volumetric counterpart is accumulate_vol!.

JustRelax.JustRelax3D.accumulate_vol! Method

julia

accumulate_vol!(EVol_pl, ε_vol_pl, dt)

Accumulate the volumetric plastic strain over a time step: EVol_pl[I] += dt * ε_vol_pl[I].

ε_vol_pl is the volumetric plastic strain rate at cell centers (= λ · (-dQ/dP), set inside the stress kernel), and is the volumetric counterpart of the deviatoric ε_pl. EVol_pl is the running invariant accumulated through time and is distinct from EII_pl (which integrates the second invariant of the deviatoric plastic strain rate via accumulate_tensor!).

JustRelax.JustRelax3D.adiabatic_heating! Method

julia

adiabatic_heating!(thermal, stokes, rheology, phases, _dt, grid)

Fill thermal.adiabatic with the adiabatic heating term inferred from the pressure change between stokes.P0 and stokes.P.

The kernels average the local thermal expansivity over the temperature nodes and scale the pressure increment by inv(dt), passed here as _dt. When stokes is nothing, the no-op overloads leave the field unchanged.

JustRelax.JustRelax3D.allzero Method

julia

allzero(x::Vararg{T,N}) where {T,N}

Check if all elements in x are zero.

Arguments

x::Vararg{T,N}: The input array.

Returns

Bool: true if all elements in x are zero, false otherwise.

JustRelax.JustRelax3D.apply_free_surface_diagonal! Method

julia

apply_free_surface_diagonal!(Dn, λmaxVn, ρgn, di_center, dt)

Add the free-surface diagonal -dt * ∂n(ρg) to the normal DYREL preconditioner and its Gershgorin row bound. The normal direction is y in 2D and z in 3D. Passing dt = 0 is a no-op.

JustRelax.JustRelax3D.assign! Method

julia

assign!(B::AbstractArray{T,N}, A::AbstractArray{T,N}) where {T,N}

Assigns the values of array A to array B in parallel.

Arguments

B::AbstractArray{T,N}: The destination array.
A::AbstractArray{T,N}: The source array.

JustRelax.JustRelax3D.center2vertex! Method

julia

center2vertex!(vertex, center)
center2vertex!(vertex_yz, vertex_xz, vertex_xy, center_yz, center_xz, center_xy)

Interpolates the values at the cell center(s) onto vertex points. The 6-argument method interpolates the three shear-stress/strain-rate components of a 3D SymmetricTensor onto their respective face vertices.

JustRelax.JustRelax3D.compute_P! Method

compute_P!(P, P0, RP, ∇V, Q, ΔT, η, rheology::NTuple{N,MaterialParams}, phase_ratio::C, dt, r, θ_dτ)

Compute the pressure field P and the residual RP for the compressible case. This function introduces thermal stresses after the implementation of Kiss et al. (2023).

Arguments

P: pressure field
RP: residual field
∇V: divergence of the velocity field
Q: volumetric source/sink term which should have the properties of dV/V_tot [m³/m³] normalized per cell, default is zero.
ΔT: temperature difference on the cell center, to account for thermal stresses. The thermal expansivity α is computed from the material parameters.
η: viscosity field
rheology: material parameters
phase_ratio: phase field
melt_fraction: melt fraction field, used for the thermal expansion coefficient if provided
dt: time step
r: relaxation parameter for the pressure update
θ_dτ: numerical parameter for the pressure update

JustRelax.JustRelax3D.compute_V! Method

julia

compute_V!(Vx, Vy, Vz, Rx, Ry, Rz, P, fx, fy, fz, τxx, τyy, τzz, τyz, τxz, τxy, ητ, ηdτ, ϕ, _dx, _dy, _dz)

Compute the 3D velocity field V from the pressure P, stress components τ, body forces f, and other parameters, with the rock ratio ϕ and grid spacing _dx, _dy, _dz.

JustRelax.JustRelax3D.compute_V! Method

julia

compute_V!(Vx, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ, _dx, _dy, dt)

Compute the velocity field V with the timestep dt from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ, grid spacing _dx, _dy, and time step dt.

JustRelax.JustRelax3D.compute_V! Method

julia

compute_V!(Vx, Vy, Rx, Ry, P, τxx, τyy, τxy, ηdτ, ρgx, ρgy, ητ, ϕ, _dx, _dy)

Compute the velocity field V from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

Each active face update uses its own liquid face mass. When the density-gradient free-surface correction is enabled, it is included implicitly in that local diagonal; inactive faces receive zero velocity and zero residual.

JustRelax.JustRelax3D.compute_Vx! Method

julia

compute_Vx!(Vx, Rx, P, τxx, τxy, ηdτ, ρgx, ητ, ϕ, _dx, _dy)

Compute the x-component of the velocity field Vx from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

JustRelax.JustRelax3D.compute_Vy! Method

julia

compute_Vy!(Vy, Vx_on_Vy, Ry, P, τyy, τxy, ηdτ, ρgy, ητ, ϕ, _dx, _dy, dt)

Compute the y-component of the velocity field Vy from the pressure P, stress components τ, and other parameters, taking into account the rock ratio ϕ, grid spacing _dx, _dy, and time step dt.

JustRelax.JustRelax3D.compute_air_ratio Method

julia

compute_air_ratio(phase_ratio, air_phase, inds...)

Compute the air ratio at the given indices based on the phase_ratio and air_phase.

JustRelax.JustRelax3D.compute_bulk_viscosity_and_penalty! Method

julia

compute_bulk_viscosity_and_penalty!(dyrel, stokes, rheology, phase_ratios, γfact, dt)

Computes the bulk viscosity ηb and the effective penalty parameter γ_eff.

Bulk Viscosity (ηb): Computed based on the bulk modulus of the material phases.

If Kb is infinite (incompressible), ηb defaults to γfact * η_mean.
Otherwise ηb = Kb * dt.

Penalty Parameter (γ_eff): A combination of numerical (γ_num) and physical (γ_phy) penalty terms.

γ_num = γfact * η_mean, with η_mean the mean of the finite viscosities
γ_phy = Kb * dt (or γ_num where Kb is infinite)
γ_eff = (γ_phy * γ_num) / (γ_phy + γ_num)

Arguments

dyrel: JustRelax.DYREL struct to update.
stokes: JustRelax.StokesArrays.
rheology: Material properties.
phase_ratios: Phase fraction information.
γfact: Numerical factor for penalty parameter (default: 20.0).
dt: Time step.

This function parallelizes the computation across grid cells.

JustRelax.JustRelax3D.compute_buoyancy Method

julia

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

julia

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

julia

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

julia

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_dissolved_volatiles! Method

julia

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.compute_dt Method

julia

compute_dt(S::JustRelax.StokesArrays, args...)

Compute the time step dt for the simulation.

JustRelax.JustRelax3D.compute_lithostatic_pressure! Method

julia

compute_lithostatic_pressure!(P, ρg, dz)
compute_lithostatic_pressure!(P, ρg, dz, igg::IGG)

Integrate the vertical component of the buoyancy force ρg down the columns of the cell-centered pressure P. The vertical direction is the last dimension of P and points upwards, so the last entry of a column is the shallowest cell.

dz is either a constant cell height or a vector holding the height of every cell along the vertical direction. Since P is cell-centered, cell j carries the weight of all the cells above it plus half of its own,

julia

P[j] = Σ_{k>j} ρg[k] * dz[k] + ρg[j] * dz[j] / 2

The three-argument method integrates each column within the local subdomain, and throws if the vertical direction is split across MPI ranks. Pass the MPI topology igg to also collect the weight of the cells held by the ranks stacked above the local subdomain; ρg must then be up to date on the halo cells, and the resulting P agrees on the cells that neighboring ranks share.

JustRelax.JustRelax3D.compute_maxloc! Method

julia

compute_maxloc!(B, A; window)

Compute the maximum value of A in the window = (width_x, width_y, width_z) and store the result in B.

JustRelax.JustRelax3D.compute_melt_fraction! Method

julia

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

julia

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_plastic_gradients_phase Method

julia

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.

JustRelax.JustRelax3D.compute_principal_stresses! Method

julia

compute_principal_stresses!(stokes, σ::PrincipalStress)

In-place version of compute_principal_stresses, writing into a pre-allocated σ.

JustRelax.JustRelax3D.compute_principal_stresses Method

julia

compute_principal_stresses(backend, stokes::StokesArrays)

Compute the principal deviatoric stresses (eigenvalues and eigenvectors of the stress tensor) at cell centers from stokes, returning a new PrincipalStress.

JustRelax.JustRelax3D.compute_rock_ratio Method

julia

compute_rock_ratio(phase_ratio, air_phase, inds...)

Compute the rock ratio at the given indices based on the phase_ratio and air_phase.

JustRelax.JustRelax3D.compute_shear_heating! Method

julia

compute_shear_heating!(thermal, stokes, rheology, dt)

Populate thermal.shear_heating from the current deviatoric stress and strain-rate fields stored in stokes.

The kernel evaluates the elastic strain-rate contribution from stokes.τ, stokes.τ_o, the shear modulus in rheology, and the pseudo-time step dt. The resulting volumetric heating term is written in place on thermal cell centers.

When phase ratios are passed as an extra positional argument, phase-weighted material properties are used instead of a single rheology state.

JustRelax.JustRelax3D.compute_strain_rate! Method

julia

compute_strain_rate!(εxx, εyy, εzz, εyz, εxz, εxy, ∇V, Vx, Vy, Vz, ϕ, _dx, _dy, _dz)

Compute the 3D components of the strain rate tensor ε from the velocity field V and its divergence ∇V, taking into account the rock ratio ϕ and grid spacing _dx, _dy, _dz.

JustRelax.JustRelax3D.compute_strain_rate! Method

julia

compute_strain_rate!(εxx, εyy, εxy, ∇V, Vx, Vy, ϕ, _dx, _dy)

Compute the components of the strain rate tensor ε from the velocity field V and its divergence ∇V, taking into account the rock ratio ϕ and grid spacing _dx, _dy.

ϕ selects which entries are computed at all; it does not scale them. The rock fraction enters the momentum balance once, where the stress divergence is taken (d_xa(τxx, ϕ.center), d_yi(τxy, ϕ.vertex)). Scaling ε by it as well would weight the deviatoric term by ϕ² while pressure and buoyancy carry ϕ, which softens partially filled cells relative to the rest of the momentum equation.

JustRelax.JustRelax3D.compute_strain_rate_from_increment! Method

julia

compute_strain_rate_from_increment!(εxx, εyy, εxy, Δεxx, Δεyy, Δεxy, ϕ, _dt)

Compute the components of the strain rate tensor ε from the strain increments Δε, taking into account the rock ratio ϕ and time step _dt.

As in compute_strain_rate!, ϕ only selects which entries are computed; the rock fraction is applied once, in the stress divergence.

JustRelax.JustRelax3D.compute_viscosity! Method

julia

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

julia

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

julia

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.compute_yieldfunction_phase Method

julia

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.

JustRelax.JustRelax3D.compute_α Method

julia

compute_α(rheology, phase)

Return the thermal expansivity α used by the adiabatic heating kernels.

phase can be a single phase index, nothing, or a phase-ratio vector. In the latter case the result is phase-weighted.

JustRelax.JustRelax3D.compute_ρg! Method

julia

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

julia

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_∇V! Method

julia

compute_∇V!(∇V, V, ϕ, _di)

Compute the divergence of the velocity field V and store it in ∇V, taking into account the rock ratio ϕ and grid spacing _di.

The divergence is evaluated only for retained pressure cells. The local continuity degree of freedom is eliminated when one of its surrounding staggered velocity faces is disconnected:

julia

             Vy[i, j+1]
                   o
                   |
    Vx[i, j]  o--- p[i,j] ---o  Vx[i+1, j]
                   |
                   o
             Vy[i, j]

This is null-space elimination, not an air velocity equation.

JustRelax.JustRelax3D.continuation_linear Method

julia

continuation_linear(x_new, x_old, ν)

Do a continuation step (1-ν)*x_old + ν*x_new with damping parameter ν

JustRelax.JustRelax3D.continuation_log Method

julia

continuation_log(x_new, x_old, ν)

Do a continuation step exp((1-ν)*log(x_old) + ν*log(x_new)) with damping parameter ν

JustRelax.JustRelax3D.displacement2velocity! Method

julia

displacement2velocity!(stokes::StokesArrays, dt)
displacement2velocity!(stokes::StokesArrays, dt, flow_bcs::AbstractFlowBoundaryConditions)

Set stokes.V (velocity) to stokes.U / dt (displacement divided by the time step), in place — the inverse of velocity2displacement!. The 3-argument method is a no-op when flow_bcs isa VelocityBoundaryConditions (velocity is already the primary variable).

JustRelax.JustRelax3D.flow_bcs! Method

julia

flow_bcs!(stokes, bcs::VelocityBoundaryConditions)
flow_bcs!(stokes, bcs::DisplacementBoundaryConditions)
flow_bcs!(bcs, Vx, Vy[, Vz])

Apply no-slip, free-slip, and periodic flow boundary conditions to staggered velocity or displacement arrays. The array form accepts the boundary-condition object first; the stokes form accepts it second. Boundary updates are executed through ParallelStencil kernels on the selected backend.

The three conditions are applied in the order no-slip, free-slip, periodic, so a face carrying more than one of them would keep only the last. The constructors reject such combinations. Faces where all three are false are left untouched. Periodic conditions match normal components at paired boundary planes and copy opposite interior values into tangential ghost planes.

JustRelax.JustRelax3D.fn_ratio Method

julia

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.getindex_NamedTuple Method

julia

getindex_NamedTuple(args::NamedTuple, [sz_min::NTuple], I...)

Sample every array in args at I, offsetting entries larger than sz_min by one to skip their ghost nodes. sz_min is the size of the grid I indexes over; pass it whenever it is known, as inferring it from args reads the ghost node as a cell center when args carries no cell-sized entry.

JustRelax.JustRelax3D.heatdiffusion_PT! Method

julia

heatdiffusion_PT!(thermal, args...; kwargs...)

Dispatch pseudo-transient thermal diffusion to the backend associated with thermal.

See the _heatdiffusion_PT! methods below for the supported argument groups: constant K and ρCp fields, or rheology-driven properties with optional phase ratios and Stokes fields for adiabatic heating.

JustRelax.JustRelax3D.interp_Vx_on_Vy! Method

julia

interp_Vx_on_Vy!(Vx_on_Vy, Vx)

Interpolates the values of Vx onto the grid points of Vy.

Arguments

Vx_on_Vy::AbstractArray: Vx at Vy grid points.
Vx::AbstractArray: Vx at its staggered grid points.

JustRelax.JustRelax3D.isvalid_c Method

julia

isvalid_c(ϕ::JustRelax.RockRatio, inds...)

Check whether the 3D pressure degree of freedom is connected to liquid. As in 2D, the row is retained when the centre carries liquid and all six adjacent velocity faces are active.

JustRelax.JustRelax3D.isvalid_c Method

julia

isvalid_c(ϕ::JustRelax.RockRatio, inds...)

Check whether the 2D pressure degree of freedom is connected to liquid.

The cell-centred pressure row is retained when ϕ.center[i,j] carries liquid and all four adjacent velocity faces are active:

julia

             Vy[i, j+1]
                   o
                   |
    Vx[i, j]  o--- p[i,j] ---o  Vx[i+1, j]
                   |
                   o
             Vy[i, j]

This is the local null-space elimination of the matrix-free reduced system. ∇V is built from these four faces, and the pressure of the cell acts back on them through the weighted gradient, so a row that keeps an inactive face is a row whose divergence no free velocity can relieve: the Powell-Hestenes penalty then drives its pressure without bound.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax3D.isvalid_v Method

julia

isvalid_v(ϕ::JustRelax.RockRatio, inds...)

Check if ϕ.vertex[inds...] is a not a nullspace in 3D.

JustRelax.JustRelax3D.isvalid_v Method

julia

isvalid_v(ϕ::JustRelax.RockRatio, inds...)

Check if ϕ.vertex[inds...] is a not a nullspace 2D.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax3D.isvalid_velocity Method

julia

isvalid_velocity(ϕ::JustRelax.RockRatio, inds...)

Check if the velocity components at the given indices are not nullspaces in 3D.

JustRelax.JustRelax3D.isvalid_velocity Method

julia

isvalid_velocity(ϕ::JustRelax.RockRatio, inds...)

Check if the velocity components at the given indices are not nullspaces in 2D.

JustRelax.JustRelax3D.isvalid_vx Method

julia

isvalid_vx(ϕ::JustRelax.RockRatio, inds...)

Check whether the x-velocity degree of freedom at inds is not a nullspace.

The face is retained when its own control volume carries liquid, ϕ.Vx > 0. A Vx/Vy control volume straddles two cell halves, so it runs dry before the cells it separates do; the face is then a rigid lid on the cut cell behind it. That is what pairs with isvalid_c: a face with no liquid contributes to the divergence of both cells it separates without being able to relieve it, so the two rules have to draw the boundary in the same place.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax3D.isvalid_vy Method

julia

isvalid_vy(ϕ::JustRelax.RockRatio, inds...)

Check whether the y-velocity degree of freedom at inds is not a nullspace. Retained when ϕ.Vy > 0; see isvalid_vx.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax3D.isvalid_vz Method

julia

isvalid_vz(ϕ::JustRelax.RockRatio, inds...)

Check whether the z-velocity degree of freedom at inds is not a nullspace. Retained when ϕ.Vz > 0; see isvalid_vx.

Arguments

ϕ::JustRelax.RockRatio: The RockRatio object to check against.
inds: Cartesian indices to check.

JustRelax.JustRelax3D.isvalid_xy Method

julia

isvalid_xy(ϕ, inds...)

Check if the xy shear component at the given indices is not a nullspace.

JustRelax.JustRelax3D.isvalid_xz Method

julia

isvalid_xz(ϕ, inds...)

Check if the xz shear component at the given indices is not a nullspace.

JustRelax.JustRelax3D.isvalid_yz Method

julia

isvalid_yz(ϕ, inds...)

Check if the yz shear component at the given indices is not a nullspace.

JustRelax.JustRelax3D.maximum_mpi Method

julia

maximum_mpi(A)

Compute the maximum value of array A across all MPI processes.

JustRelax.JustRelax3D.mean_mpi Method

julia

mean_mpi(A)

Compute the mean of array A across all MPI processes.

JustRelax.JustRelax3D.minimum_mpi Method

julia

minimum_mpi(A)

Compute the minimum value of array A across all MPI processes.

JustRelax.JustRelax3D.momentum_rows Method

julia

momentum_rows(ni::NTuple{N,Integer}, periodic::NTuple{N,Bool}, d)

Shape of the momentum residual of direction d on a grid of ni cells.

A non-periodic direction has ni[d] - 1 interior faces to solve, the two boundary faces being prescribed. A periodic direction has ni[d]: its two boundary faces are the same plane, so they form one extra unknown, stored as the last row of the residual and mapped to the upper face.

JustRelax.JustRelax3D.multi_copy! Method

julia

multi_copy!(dst::NTuple{N, T}, src::NTuple{N, T}) where {N, T}

Copy data from the tuple of arrays src into the tuple of arrays dst in parallel.

JustRelax.JustRelax3D.norm_mpi Method

julia

norm_mpi(A)

Compute the L2 norm of array A across all MPI processes.

JustRelax.JustRelax3D.pureshear_bc! Method

julia

pureshear_bc!(stokes, xci, xvi, εbg)
pureshear_bc!(stokes, xci, xvi, εbg, backend)

Initialize a pure-shear background velocity field on the staggered grids. xci contains cell-center coordinates and xvi contains velocity-grid coordinates; each component is built from the vertex coordinates of its own direction. In 2D the kernels set Vx = εbg*x and Vy = -εbg*y. In 3D they set Vx = εbg*x, Vy = εbg*y, and Vz = -εbg*z. Ghost layers are left untouched so that subsequent flow boundary-condition and halo updates can set them consistently.

All field updates are performed by ParallelStencil kernels on the backend of stokes. The five-argument form remains available for compatibility; its backend argument is redundant, as the backend is inferred from stokes.

JustRelax.JustRelax3D.rotate_stress! Method

julia

rotate_stress!(τ_particles::StressParticles, stokes, particles, dt)

Interpolate the current deviatoric stress stokes.τ and vorticity stokes.ω onto the particles and rotate the particle stress over dt. stokes.ω must hold the vorticity of the current velocity field. Use stress2grid! afterwards to map the rotated stress back onto stokes.τ_o.

JustRelax.JustRelax3D.rotate_stress_particles! Method

julia

rotate_stress_particles!(τ::NTuple, ω::NTuple, particles::Particles, dt; method = :matrix)

Rotate the deviatoric stress carried by each active particle over dt with the local vorticity, using GeoParams' elastic stress rotation. τ holds the stress components and ω the vorticity components, as particle cell arrays. method is accepted for call-site compatibility and does not select an algorithm.

JustRelax.JustRelax3D.shear2center! Method

julia

shear2center!(A::SymmetricTensor)

Interpolate the shear components of A onto the cell centers, in place.

JustRelax.JustRelax3D.simpleshear_bc! Method

julia

simpleshear_bc!(stokes, xci, xvi, γbg)
simpleshear_bc!(stokes, xci, xvi, γbg, backend)

Initialize an xy simple-shear background velocity field on the staggered grids. The imposed field is Vx = γbg * y; the other velocity components are set to zero. xci contains cell-center coordinates and xvi contains velocity-grid coordinates. Ghost layers are left untouched.

All field updates are performed by ParallelStencil kernels on the backend of stokes. The five-argument form remains available for compatibility; its backend argument is redundant, as the backend is inferred from stokes.

JustRelax.JustRelax3D.solve! Method

julia

solve!(stokes::StokesArrays, pt_stokes, grid, flow_bcs, ρg, phase_ratios, rheology, args, dt, igg; kwargs...)

Solve the 3D viscoelastoplastic Stokes equations to pseudo-transient convergence, updating stokes (velocity, pressure, stress, viscosity) in place for one physical time step dt. This is the general, multi-phase call form used by most models; dispatch on the type of the sixth argument also accepts, for simpler/benchmark setups:

a single-phase rheology::GeoParams.MaterialParams in place of phase_ratios/rheology (drop the phase_ratios argument), or
constant K, G (bulk and shear modulus) fields in place of phase_ratios/rheology/args, for linear viscoelastic problems with no material rheology.

grid may also be replaced by the grid spacing di alone (a NTuple/NamedTuple).

Arguments

stokes: solver state allocated with StokesArrays.
pt_stokes: pseudo-transient coefficients, e.g. from PTStokesCoeffs.
grid: the model Geometry.
flow_bcs: velocity/displacement boundary conditions.
ρg: buoyancy forcing (ρgx, ρgy, ρgz).
phase_ratios: per-cell/per-node phase fractions (a JustPIC.PhaseRatios).
rheology: one GeoParams.MaterialParams per phase.
args: auxiliary fields (e.g. temperature T, pressure P) used by the constitutive updates.
dt: physical time step.
igg: the distributed-grid context (IGG).

Keyword arguments

iterMax = 10e3: maximum pseudo-transient iterations.
nout = 500: check convergence every nout iterations.
viscosity_relaxation = 1e-2: relaxation factor for nonlinear viscosity updates.
viscosity_cutoff = (-Inf, Inf): clamp bounds for the effective viscosity.
λ_relaxation = 0.2: relaxation factor for the plastic multiplier.
b_width = (4, 4, 4): halo width used when overlapping communication and computation.
verbose = true: print convergence progress.

Dispatches on the CPU/CUDA/AMDGPU backend selected by stokes.

JustRelax.JustRelax3D.solve_DYREL! Method

julia

solve_DYREL!(
    stokes, ρg, dyrel, flow_bcs, phase_ratios, rheology, args, grid, dt, igg;
    kwargs...,
)

Solve the Stokes system with the self-tuned dynamic relaxation (DYREL) method.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
ρg: buoyancy forces arrays.
dyrel: DYREL-specific parameters and fields.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
grid: Geometry object carrying grid spacing and staggered-grid coordinates. A legacy 2D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
λ_relaxation_DR: Relaxation factor for dynamic relaxation. Default: 1.
λ_relaxation_PH: Relaxation factor for Powell-Hestenes iterations. Default: 1.
pressure_relaxation: Relaxation factor for the Powell-Hestenes pressure update. Default: 1.
iterMax_PH: Maximum number of Powell-Hestenes passes. Default: 1.0e3.
iterMax_DR: Maximum number of iterations for each dynamic-relaxation solve. Default: 50.0e3.
iterMax: Compatibility alias for iterMax_DR; used when iterMax_DR is not given.
total_iterMax: Maximum number of total dynamic-relaxation iterations. Default: 50.0e3.
nout: Output frequency for residuals. Default: 100.
rel_drop: Relative residual drop tolerance. Default: 1.0e-2.
verbose_PH: Print Powell-Hestenes iteration info. Default: true.
verbose_DR: Print Dynamic Relaxation iteration info. Default: true.
linear_viscosity: Whether to use linear viscosity. Default: false.
free_surface: Include the density-gradient free-surface stabilization term. Default: false.
update_material: Recompute viscosity and buoyancy from rheology. Set to false when those fields are prescribed by the caller. Default: true.

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax3D.solve_VariationalDYREL! Method

julia

solve_VariationalDYREL!(stokes, ρg, dyrel, flow_bcs, phase_ratios, ϕ,
    rheology, args, grid, dt, igg; kwargs...)

Solve the 2D variational Stokes problem with DYREL relaxation and the RockRatio volume weights. This is a separate entry point from solve_DYREL!; the latter remains the standard, unweighted DYREL solver.

Center fractions weight pressure and normal stress, vertex fractions weight shear stress, and face fractions weight momentum rows. Rows whose volume fraction vanishes are eliminated rather than solved with air properties.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
ρg: buoyancy forces arrays.
dyrel: DYREL-specific parameters and fields, built with the same ϕ.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
ϕ: JustRelax.RockRatio carrying the cell, vertex and face volume fractions.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
grid: Geometry{2} object carrying grid spacing and staggered-grid coordinates. A legacy 2D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

air_phase: Phase index excluded from material averages; 0 disables the correction. Default: 0.
viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
λ_relaxation_DR: Relaxation factor for dynamic relaxation. Default: 1.
λ_relaxation_PH: Relaxation factor for Powell-Hestenes iterations. Default: 1.
pressure_relaxation: Relaxation factor for the Powell-Hestenes pressure update. Default: 1.
iterMax_PH: Maximum number of Powell-Hestenes passes. Default: 1.0e3.
iterMax_DR: Maximum number of iterations for each dynamic-relaxation solve. Default: 50.0e3.
iterMax: Alias for iterMax_DR; used when iterMax_DR is not given.
total_iterMax: Maximum number of total dynamic-relaxation iterations. Default: 50.0e3.
nout: Output frequency for residuals. Default: 100.
rel_drop: Relative residual drop tolerance. Default: 1.0e-2.
verbose_PH: Print Powell-Hestenes iteration info. Default: true.
verbose_DR: Print Dynamic Relaxation iteration info. Default: true.
linear_viscosity: Whether to use linear viscosity. Default: false.
free_surface: Include the density-gradient free-surface stabilization term. Default: false.

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax3D.solve_VariationalStokes! Method

julia

solve_VariationalStokes!(stokes::JustRelax.StokesArrays, args...; kwargs...)

Solve the 3D volume-fraction variational Stokes problem with matrix-free pseudo-transient iterations.

ϕ carries liquid weights at pressure cells, stress vertices and staggered velocity faces; zero-weight rows are written as zero instead of being solved with air material properties. See the 2D method for the pressure/velocity connectivity rule the weights encode.

Arguments (in the following order)

stokes: JustRelax.StokesArrays containing the simulation fields.
pt_stokes: Pseudo-transient coefficients, from PTStokesCoeffs.
grid: Geometry{3} object carrying grid spacing and staggered-grid coordinates. A legacy 3D spacing tuple or named tuple is also accepted and converted to a uniform Geometry.
flow_bcs: AbstractFlowBoundaryConditions defining velocity boundary conditions.
ρg: buoyancy forces arrays.
phase_ratios: JustPIC.PhaseRatios for material phase tracking.
ϕ: JustRelax.RockRatio carrying the cell, vertex and face volume fractions.
rheology: Material properties and rheological laws.
args: Tuple of additional arguments needed to update viscosity, stress, and buoyancy forces.
dt: Time step.
igg: IGG object for global grid information (MPI).

Keyword Arguments

air_phase: Phase index excluded from material averages; 0 disables the correction. Default: 0.
viscosity_cutoff: Limits for viscosity (min, max). Default: (-Inf, Inf).
viscosity_relaxation: Relaxation factor for viscosity updates. Default: 1.0e-2.
iterMax: Maximum number of pseudo-transient iterations. Default: 10.0e3.
nout: Output frequency for residuals. Default: 500.
verbose: Print iteration info. Default: true.
b_width: Halo width used to overlap communication with computation. Default: (4, 4, 4).

Options may be passed either as plain keywords or bundled as a single kwargs = (; ...) NamedTuple.

JustRelax.JustRelax3D.stress2grid! Method

julia

stress2grid!(stokes, τ_particles::StressParticles, particles)

Interpolate the particle stress in τ_particles back onto the old-stress fields stokes.τ_o: normal components onto the cell centers, and shear components onto the vertices in 2D or onto the cell centers and edges in 3D, matching where the stress kernels read them from. Counterpart of rotate_stress!, and the step that hands the rotated stress to the next Stokes solve.

JustRelax.JustRelax3D.subgrid_characteristic_time! Method

julia

subgrid_characteristic_time!(subgrid_arrays, particles, dt₀, phases, rheology, thermal::ThermalArrays, stokes::StokesArrays[, di])

Compute, per cell, the characteristic thermal diffusion timescale dt₀ = ρCp / (2 K Σ dxi⁻²) used for JustPIC's subgrid-diffusion correction of particle temperature, evaluating rheology's density/heat-capacity/conductivity at the local phase (from phases, either a JustPIC.PhaseRatios or an integer phase-id array) and temperature/pressure. dt₀ must be ghosted, with size size(stokes.P) .+ 2; values are written to its interior region.

JustRelax.JustRelax3D.sum_mpi Method

julia

sum_mpi(A)

Compute the sum of array A across all MPI processes.

JustRelax.JustRelax3D.take Method

julia

take(fldr::String)

Create folder fldr if it does not exist.

JustRelax.JustRelax3D.tensor_invariant! Method

julia

tensor_invariant!(A::JustRelax.SymmetricTensor)

Compute the tensor invariant of the given symmetric tensor A.

Arguments

A::JustRelax.SymmetricTensor: The input symmetric tensor.

JustRelax.JustRelax3D.thermal_bcs! Method

julia

thermal_bcs!(thermal, bcs::TemperatureBoundaryConditions)
thermal_bcs!(T::AbstractArray, bcs::TemperatureBoundaryConditions)

Apply thermal ghost-cell boundary conditions to a temperature field.

thermal_bcs! applies the scalar temperature conditions stored in bcs:

constant_value faces are applied first using Tghost = 2 * value - Tinterior.
no_flux faces are applied next by copying the adjacent interior temperature.
periodic faces are applied last by copying the opposite interior temperature into the ghost layer.

Faces set to false are ignored. Periodic faces must be paired by direction and cannot also carry another thermal condition. Prescribed constant_flux values are not applied here; they are consumed by the pseudo-transient heat-diffusion compute_flux! kernels.

JustRelax.JustRelax3D.unpack Method

julia

unpack(x::T)

Generated function to unpack the fields of the struct x into a tuple.

JustRelax.JustRelax3D.update_T Method

julia

update_T(::Nothing, b_width, thermal, ρCp, pt_thermal, dirichlet, _dt, _di, ni)
update_T(::Nothing, b_width, thermal, rheology, phase, pt_thermal, dirichlet, _dt, _di, ni, args)

Launch the pseudo-transient temperature update kernel over the active thermal domain.

These wrappers select the appropriate kernel overload depending on whether the solver works with precomputed ρCp fields or rheology-derived properties.

JustRelax.JustRelax3D.update_dτV_α_β! Method

julia

update_dτV_α_β!(dτV, βV, αV, cV, λmaxV, CFL_v)

Generic N-dimensional version (2D and 3D) for updating pseudo-time step and acceleration parameters.

Computes the pseudo-time step dτV from the maximum eigenvalue estimate λmaxV and CFL number, then updates the damping parameters βV and acceleration parameters αV accordingly.

Arguments

dτV::NTuple{N, AbstractArray{T, N}}: Tuple of pseudo-time steps for each velocity component
βV::NTuple{N, AbstractArray{T, N}}: Tuple of damping parameters β for each velocity component
αV::NTuple{N, AbstractArray{T, N}}: Tuple of acceleration parameters α for each velocity component
cV::NTuple{N, AbstractArray{T, N}}: Tuple of preconditioner diagonal entries for each velocity component
λmaxV::NTuple{N, AbstractArray{T, N}}: Tuple of maximum eigenvalue estimates for each velocity component
CFL_v::Real: CFL number for velocity

JustRelax.JustRelax3D.update_phase_ratios_2D! Method

julia

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:

julia

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

julia

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:

julia

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

julia

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.

JustRelax.JustRelax3D.update_pt_thermal_arrays! Method

julia

update_pt_thermal_arrays!(pt_thermal, phase_ratios, rheology, args, _dt)

Recompute the pseudo-transient thermal coefficient arrays stored in pt_thermal from phase-weighted material properties.

This helper is used by the pseudo-transient thermal solver when the local phase mixture changes over time.

JustRelax.JustRelax3D.update_rock_ratio! Method

julia

update_rock_ratio!(ϕ::JustRelax.RockRatio, phase_ratios, air_phase)

Update the rock ratio ϕ based on the provided phase_ratios and air_phase.

Arguments

ϕ::JustRelax.RockRatio: The rock ratio object to be updated.
phase_ratios: The ratios of different phases present.
air_phase: The phase representing air.

JustRelax.JustRelax3D.update_rock_ratio! Method

julia

update_rock_ratio!(ϕ::JustRelax.RockRatio, phase_ratios, air_phase)

Update the rock ratio ϕ for a 3D grid based on the provided phase_ratios and air_phase.

Arguments

ϕ::JustRelax.RockRatio: The rock ratio object to be updated.
phase_ratios: The ratios of different phases present.
air_phase: The phase representing air.

JustRelax.JustRelax3D.update_rock_ratio_cv! Method

julia

update_rock_ratio_cv!(ϕ, ratio_center, ratio_vertex, air_phase)

Update the rock ratio for both center and vertex values based on the provided ratio_center, ratio_vertex, and air_phase.

JustRelax.JustRelax3D.update_thermal_coeffs! Method

julia

update_thermal_coeffs!(pt_thermal, rheology, phase_ratios, args, dt)
update_thermal_coeffs!(pt_thermal, rheology, args, dt)

Refresh the pseudo-transient coefficient arrays stored in pt_thermal.

Use this when conductivity or volumetric heat capacity changes during the solve, for example because temperature, pressure, or phase ratios evolve. The update is performed in place on pt_thermal.θr_dτ and pt_thermal.dτ_ρ.

JustRelax.JustRelax3D.update_α_β! Method

julia

update_α_β!(βV, αV, dτV, cV)

Generic N-dimensional version (2D and 3D) of the acceleration parameters update.

Updates the damping parameters βV and αV for each velocity component based on the pseudo-time step dτV and the preconditioner diagonal cV.

Arguments

βV::NTuple{N, AbstractArray{T, N}}: Tuple of damping parameters β for each velocity component
αV::NTuple{N, AbstractArray{T, N}}: Tuple of acceleration parameters α for each velocity component
dτV::NTuple{N, AbstractArray{T, N}}: Tuple of pseudo-time steps for each velocity component
cV::NTuple{N, AbstractArray{T, N}}: Tuple of preconditioner diagonal entries for each velocity component

JustRelax.JustRelax3D.velocity2center! Method

julia

velocity2center!(Vx_c, Vy_c, Vx, Vy)

In-place interpolation of the velocity field Vx, Vy from a staggered grid with ghost nodes onto the pre-allocated Vx_c, Vy_c 2D arrays located at the cell centers.

JustRelax.JustRelax3D.velocity2center! Method

julia

velocity2center!(Vx_c, Vy_c, Vz_c, Vx, Vy, Vz)

In-place interpolation of the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the pre-allocated Vx_c, Vy_c, Vz_c 3D arrays located at the cell centers.

JustRelax.JustRelax3D.velocity2displacement! Method

julia

velocity2displacement!(stokes::StokesArrays, dt)

Set stokes.U (displacement) to stokes.V * dt (velocity times the time step), in place.

JustRelax.JustRelax3D.velocity2vertex! Method

julia

velocity2vertex!(Vx_v, Vy_v, Vx, Vy)

In-place interpolation of the velocity field Vx, Vy from a staggered grid with ghost nodes onto the pre-allocated Vx_v, Vy_v 2D arrays located at the grid vertices.

JustRelax.JustRelax3D.velocity2vertex! Method

julia

velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz)

In-place interpolation of the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the pre-allocated Vx_d, Vy_d, Vz_d 3D arrays located at the grid vertices.

JustRelax.JustRelax3D.velocity2vertex Method

julia

velocity2vertex(Vx, Vy, Vz)

Interpolate the velocity field Vx, Vy, Vz from a staggered grid with ghost nodes onto the grid vertices.

JustRelax.JustRelax3D.vertex2center! Method

julia

vertex2center!(center, vertex; ghost_x = false, ghost_y = false, ghost_z = false)

Interpolates the values at the vertex onto center points.

center may carry a ring of ghost nodes that vertex does not; each ghost_* keyword shifts the write index by one along that dimension, so the interpolated block lands in the interior of center instead of overwriting its ghost nodes. ghost_z is ignored in 2D.

JustRelax.JustRelax3D.viscosity_phase_ratio Method

julia

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.@add Macro

julia

@add(I, args...)

Add I to the scalars in args

JustRelax.JustRelax3D.@allocate Macro

julia

@allocate(ni...)

Convenience macro to allocate a PTArray of size ni... with undef values.

JustRelax.JustRelax3D.@copy Macro

julia

copy(B, A)

convenience macro to copy data from the array A into array B

JustRelax.JustRelax3D.@displacement Macro

julia

@displacement(U)

Unpacks the displacement arrays U from the StokesArrays A.

JustRelax.JustRelax3D.@idx Macro

julia

@idx(args...)

Make a linear range from 1 to args[i], with i ∈ [1, ..., n]

JustRelax.JustRelax3D.@normal Macro

julia

@normal(A)

Unpacks the normal components of the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@plastic_strain Macro

julia

@plastic_strain(A)

Unpacks the plastic strain rate tensor ε_pl from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@qT Macro

julia

@qT(V)

Unpacks the flux arrays qT_i from the ThermalArrays A.

JustRelax.JustRelax3D.@qT2 Macro

julia

@qT2(V)

Unpacks the flux arrays qT2_i from the ThermalArrays A.

JustRelax.JustRelax3D.@residuals Macro

julia

@residuals(A)

Unpacks the momentum residuals from A.

JustRelax.JustRelax3D.@shear Macro

julia

@shear(A)

Unpacks the shear components of the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@shear_center Macro

julia

@shear_center(A)

Unpacks the shear components of the symmetric tensor A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@strain Macro

julia

@strain(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@strain_center Macro

julia

@strain_center(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@strain_increment Macro

julia

@strain_increment(A)

Unpacks the strain rate tensor ε from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@stress Macro

julia

@stress(A)

Unpacks the deviatoric stress tensor τ from the StokesArrays A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@stress_center Macro

julia

@stress_center(A)

Unpacks the deviatoric stress tensor τ from the StokesArrays A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@tensor Macro

julia

@tensor(A)

Unpacks the symmetric tensor A, where its components are defined in the staggered grid. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@tensor_center Macro

julia

@tensor_center(A)

Unpacks the symmetric tensor A, where its components are defined in the center of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@tensor_vertex Macro

julia

tensor_vertex(A)

Unpacks the symmetric tensor A, where its components are defined in the vertices of the grid cells. Shear components are unpack following Voigt's notation.

JustRelax.JustRelax3D.@tuple Macro

julia

@tuple(A)

Convenience maktro to unpack the fields of the struct A into a tuple. Works with Velocity and SymmetricTensor structs.

JustRelax.JustRelax3D.@unpack Macro

julia

@unpack(x)

Convenience macro to unpack the fields of the struct x into a tuple.

JustRelax.JustRelax3D.@velocity Macro

julia

@velocity(V)

Unpacks the velocity arrays V from the StokesArrays A.

JustRelax.DataIO Module

julia

JustRelax.DataIO

Checkpointing and output writing.

checkpointing_hdf5/load_checkpoint_hdf5 and checkpointing_jld2/load_checkpoint_jld2 save and restore the model state, while save_vtk, save_particles, and save_marker_chain write fields, particles, and marker chains for ParaView. Every routine moves the data to the CPU before writing, whatever backend the fields live on. The submodule is loaded together with JustRelax; its names are reached as JustRelax.DataIO.f or by using JustRelax.DataIO.

JustRelax.DataIO.VTKDataSeries Type

julia

VTKDataSeries(full_name::String, xi)

A ParaView collection (.pvd) accumulating a time series of VTK snapshots on the grid xi, written under full_name. Add each time step with Base.append!.

JustRelax.DataIO.add_field! Method

julia

add_field!(vtk, name, array, npoints, ncells, precision)

Write array to vtk as point data or cell data, whichever its size matches. A vertex grid holds one cell fewer than it has nodes per dimension, so the two sizes are always distinguishable.

JustRelax.DataIO.append! Method

julia

append!(data_series::VTKDataSeries, data::NamedTuple, time_step, seconds)

Write data as a new VTK snapshot on data_series's grid and add it to the time series at time seconds, under the name time_step.

JustRelax.DataIO.center_coordinates Method

julia

center_coordinates(grid::Geometry)

The cell-center coordinate vectors of grid (grid.xci), collected into plain Vectors for serialization.

JustRelax.DataIO.checkpointing_hdf5 Method

julia

checkpointing_hdf5(dst, stokes, T, time, timestep; precision = Float32)

Save the state of the model at time in dst as the HDF5 file checkpoint.h5, so that the run can be restarted from it with load_checkpoint_hdf5.

The file holds the velocity components of stokes (Vx, Vy, and Vz in 3D), its pressure and viscosity, the temperature T, and time/timestep. Fields are transferred to the CPU before writing, whatever the backend. The file is written to a temporary directory first and moved into place, so an interrupted call leaves any previous checkpoint intact.

Arguments

dst: Directory the checkpoint is written to; created if it does not exist.
stokes: JustRelax.StokesArrays holding the velocity, pressure, and viscosity fields.
T: Temperature field.
time: Simulation time.
timestep: Time step.

Keyword arguments

precision: element type the arrays are converted to before writing.

JustRelax.DataIO.checkpointing_jld2 Method

julia

checkpointing_jld2(dst, stokes, [thermal,] time, timestep[, igg]; kwargs...)

Save necessary data in dst as a jld2 file to restart the model from the state at time. If run in parallel, the file will be named after the corresponidng rank e.g. checkpoint0000.jld2 and thus can be loaded by the processor while restarting the simulation. If you want to restart your simulation from the checkpoint you can use load() and specify the MPI rank by providing a dollar sign and the rank number.

Arguments

dst: The destination directory where the checkpoint file will be saved.
stokes: The stokes flow variables to be saved.
thermal: (Optional) The thermal variables to be saved.
time: The current simulation time.
timestep: The current timestep.
igg: (Optional) The IGG struct for parallel runs.

Keyword Arguments

kwargs...: Additional variables to be saved in the checkpoint file. These will be added to the base checkpoint data. Example
julia

JustRelax.DataIO.load_checkpoint_hdf5 Method

julia

load_checkpoint_hdf5(file_path)

Load the state of the simulation from an .h5 file.

Arguments

file_path: The path to the .h5 file.

Returns

P: The loaded state of the pressure variable.
T: The loaded state of the temperature variable.
Vx: The loaded state of the x-component of the velocity variable.
Vy: The loaded state of the y-component of the velocity variable.
Vz: The loaded state of the z-component of the velocity variable, or nothing for a 2D checkpoint.
η: The loaded state of the viscosity variable.
t: The loaded simulation time.
dt: The loaded time step.

All arrays are returned on the CPU, at the precision they were written with; move them to the device with PTArray(backend)(A) before copying them back into a StokesArrays. See checkpointing_hdf5 for the writing side.

Example

julia

file_path = joinpath("path/to/your/output", "checkpoint.h5")
P, T, Vx, Vy, Vz, η, t, dt = load_checkpoint_hdf5(file_path)

JustRelax.DataIO.load_checkpoint_jld2 Method

julia

load_checkpoint_jld2(file_path)

Load the state of the simulation from a .jld2 file.

Arguments

file_path: The path to the .jld2 file.

Returns

stokes: The loaded state of the stokes variable.
thermal: The loaded state of the thermal variable. Can be nothing if not present in the file.
time: The loaded simulation time.
timestep: The loaded time step.
igg: The IGG struct if needed for parallel runs.

Example

julia

stokes, thermal, time, timestep = load_checkpoint_jld2("path/to/checkpoint.jld2", igg)

or without thermal

julia

stokes, _, time, timestep = load_checkpoint_jld2("path/to/checkpoint.jld2")

JustRelax.DataIO.metadata Method

julia

metadata(src, dst, files...)

Copy files..., Manifest.toml, and Project.toml from src to dst

JustRelax.DataIO.pack_velocity Function

julia

pack_velocity(velocity::Tuple, precision, slices = nothing)

Pack the velocity components into the (3, size...) array expected for a VTK vector attribute. Readers copy three values per tuple regardless of the declared NumberOfComponents, so a two-component array leaves the third component undefined and ParaView orients its glyphs from uninitialized memory. Components beyond N are written as zeros.

slices optionally restricts each component to a sub-range (used to trim the ghost layer of MPI-distributed arrays).

JustRelax.DataIO.save_data Method

julia

save_data(file, data, precision)
save_data(file, grid::Geometry)

Write data (converted to precision) into the open HDF5 file under its own variable name. The Geometry method instead writes the cell-center/vertex coordinate vectors (Xc/Yc[/Zc], Xv/Yv[/Zv]). Used internally by save_hdf5.

JustRelax.DataIO.save_hdf5 Method

julia

save_hdf5(dst, fname, data...)

Save each entry of data as the fname.h5 HDF5 file in the folder dst, creating the folder if it does not exist.

JustRelax.DataIO.save_hdf5 Method

julia

save_hdf5(fname, data...; precision = Float32)

Save each entry of data as the fname.h5 HDF5 file, one variable per entry, named after the variable passed at the call site.

Keyword arguments

precision: element type the arrays are converted to before writing.

JustRelax.DataIO.save_marker_chain Method

julia

save_marker_chain(fname::String, chain::MarkerChain; conversion=1.0e3, pvd=nothing, t=0.0)

Save a vector of points as a line in a VTK file.

Arguments

fname::String: The name of the VTK file to save. The extension .vtk will be appended to the name.
chain::MarkerChain: Marker chain object from JustPIC.jl.
conversion: Conversion factor for coordinates (default: 1.0e3)
pvd::Union{Nothing, String}: Optional ParaView collection filename for time series
t::Number: Time value (default: 0.0)

JustRelax.DataIO.save_particles Method

julia

save_particles(particles::Particles{B, 2}, pPhases; conversion = 1e3, fname::String = "./particles", pvd=nothing, t=0.0) where B

Save particle data and their material phase to a VTK file.

Arguments

particles::Particles{B, 2}: The particle data, where B is the type of the particle coordinates.
pPhases: The phases of the particles.
conversion: A conversion factor for the particle coordinates (default is 1e3).
fname::String: The name of the VTK file to save (default is "./particles").
pvd::Union{Nothing, String}: Optional ParaView collection filename for time series
t::Number: Time value (default: 0.0)

JustRelax.DataIO.save_particles Method

julia

save_particles(particles, igg::IGG; pPhases=nothing, conversion=1e3, fname="./particles", pvd=nothing, t=0.0, precision=Float32)

Parallel (MPI) counterpart of save_particles: each rank writes its own active particles as an unstructured .vtu piece, tied together by <fname>.pvtu. Works for 2D and 3D particles. pPhases (a CellArray of phase ids) is written as the phase point field when given, otherwise a constant is used. If pvd is given, the .pvtu datasets are appended to <pvd>.pvd at time t to build a time series (only rank 0 touches the collection).

JustRelax.DataIO.save_particles Method

julia

save_particles(particles::Particles{B, 2}; conversion = 1e3, fname::String = "./particles", pvd=nothing, t=0.0) where B

Save particle data to a VTK file.

Arguments

particles::Particles{B, 2}: The particle data, where B is the type of the particle coordinates.
conversion: A conversion factor for the particle coordinates (default is 1e3).
fname::String: The name of the VTK file to save (default is "./particles").
pvd::Union{Nothing, String}: Optional ParaView collection filename for time series
t::Number: Time value (default: 0.0)

JustRelax.DataIO.save_pvtk Method

julia

save_pvtk(fname, di::NTuple{N}, data_v::NamedTuple, data_c::NamedTuple, velocity::Tuple, igg::IGG; t=nothing, precision=Float32, pvd=nothing)

Parallel (MPI) counterpart of the serial save_vtk for an ImplicitGlobalGrid-distributed grid (requires ImplicitGlobalGrid ≥ 0.17). Writes vertex fields data_v + velocity to <fname>_vertex.pvti and cell fields data_c to <fname>_center.pvti, one .vti piece per rank. di is the global grid spacing (e.g. grid.di.center). Ranks overlap by one ghost layer, so update_halo! before writing. If pvd is given, the datasets are appended to <pvd>_vertex.pvd / <pvd>_center.pvd at time t.

JustRelax.DataIO.save_vtk Method

julia

save_vtk(fname::String, xvi, xci, data_v::NamedTuple, data_c::NamedTuple, velocity; t=0, pvd=nothing)

Save vertex and cell data to a single VTK file. The file holds the grid spanned by the vertices xvi; data_v and velocity are written as point data and data_c as cell data of that same grid.

Arguments

fname::String: The filename for the VTK file (without extension)
xvi: Vertex coordinates (tuple of coordinate arrays)
xci: Cell center coordinates (tuple of coordinate arrays); must have one entry fewer per dimension than xvi
data_v::NamedTuple: Data defined at vertices
data_c::NamedTuple: Data defined at cell centers. Fields of data_v and data_c are written as point or cell data according to their size, so a cell-centered field passed in data_v still lands on the cells
velocity::Tuple: Velocity components, each an array defined at the vertices
t::Number: Time value (default: 0)
pvd::Union{Nothing, String}: Optional ParaView collection filename. If provided, the VTK file will be added to a time series collection. WriteVTK.jl automatically handles creating new collections or appending to existing ones.

Examples

julia

# Basic usage
save_vtk("output", xvi, xci, data_v, data_c, velocity; t=1.0)

# With ParaView collection for time series
save_vtk("timestep_001", xvi, xci, data_v, data_c, velocity; t=1.0, pvd="simulation")
save_vtk("timestep_002", xvi, xci, data_v, data_c, velocity; t=2.0, pvd="simulation")
# This creates simulation.pvd containing the time series

# Time series example
times = 0:0.1:10
for (i, t) in enumerate(times)
    fname = "timestep_$(lpad(i, 3, '0'))"
    save_vtk(fname, xvi, xci, data_v, data_c, velocity; t=t, pvd="full_simulation")
end

JustRelax.DataIO.save_vtk Method

julia

save_vtk(fname::String, xci, data_c::NamedTuple, velocity; t=nothing, pvd=nothing)

Save VTK data with cell-centered data and velocity field.

Arguments

fname::String: The filename for the VTK file (without extension)
xci: Cell center coordinates (tuple of coordinate arrays)
data_c::NamedTuple: Data defined at cell centers
velocity::Tuple: Velocity components, each an array defined on the grid nodes
t::Number: Time value (default: nothing)
pvd::Union{Nothing, String}: Optional ParaView collection filename. If provided, the VTK file will be added to a time series collection. WriteVTK.jl automatically handles creating new collections or appending to existing ones.

Examples

julia

# Basic usage
save_vtk("output", xci, data_c, velocity; t=1.0)

# With ParaView collection
save_vtk("timestep_001", xci, data_c, velocity; t=1.0, pvd="simulation")

JustRelax.DataIO.vertex_coordinates Method

julia

vertex_coordinates(grid::Geometry)

The cell-vertex coordinate vectors of grid (grid.xvi), collected into plain Vectors for serialization.

The 2D variational-Stokes entry point is JustRelax.JustRelax2D.solve_VariationalStokes!. Its volume-fraction container and update routines are documented on the 2D variational Stokes page.

The dynamic-relaxation counterpart is JustRelax.JustRelax2D.solve_VariationalDYREL!. It requires a 2D RockRatio in both the DYREL constructor and solver call; see DYREL.