API

The API reference is organized by responsibility:

• Interfaces and numerical backends on this page.

• Physics models and solvers: transport, drying, reactive transport, Biot coupling, poroplasticity and homogenization.

• Constitutive laws: retention, permeability, tortuosity, elasticity and effective stress.

• Plasticity and material models: BBM, Drucker-Prager and coupling to a Biot medium.

PoroMechanics Module

PoroMechanics.jl

Reactive transport and poromechanics of porous media, on finite volume and finite element backends: unsaturated flow, non-isothermal drying, Biot poroelasticity and reactive transport in cementitious materials.

Architecture

Two layers, so that a physics model stays a self-contained description of its own equations and knows nothing about time stepping or assembly:

Layer 1 — Core abstractions (this module)

• AbstractPoroModel — supertype of every physics model. Material parameters live in the concrete struct; multiple dispatch selects the constitutive behavior.

• AbstractPoroSolver — supertype for time-stepping strategies (Monolithic, SNIA…).

Layer 2 — Physics models

Each model lives in its own file under src/Models/ and defines:

• A concrete <Name>Model <: AbstractPoroModel struct holding material parameters.

• The required interface methods (storage!, flux!, bcondition! for FVM models; assemble_element! for FEM models).

Jacobians are never written by hand: VoronoiFVM.jl differentiates the FVM callbacks with ForwardDiff.jl.

Backend convention

• Transport / diffusion / flow → VoronoiFVM.jl

• Coupled mechanics (Biot, BBM…) → Ferrite.jl

License

MIT — see LICENSE.

source

Model and solver types

PoroMechanics.AbstractPoroModel Type

AbstractPoroModel

Supertype for all physics models in PoroMechanics.jl.

Concrete subtypes must implement the methods appropriate for their backend:

FVM backend (VoronoiFVM.jl) — transport / diffusion

• storage!(f, u, node, model, data) — stored amount M (the solver forms ∂M/∂t)

• flux!(f, u, edge, model, data) — inter-node fluxes (Darcy, Fick, Fourier…)

• bcondition!(f, u, node, model, data) — boundary conditions

FEM backend (Ferrite.jl) — coupled mechanics

• assemble_element!(Ke, re, el, u_el, model, cv, Δt) — element stiffness & residual

A model must also expose:

• nspecies(model)::Int — number of primary unknowns

• species_names(model)::Vector{Symbol} — human-readable names for unknowns

source

PoroMechanics.AbstractPoroSolver Type

AbstractPoroSolver

Supertype for time-stepping strategies.

source

Finite volume interface

Callbacks consumed by VoronoiFVM.jl. A model implements the ones its physics needs; the Jacobian is obtained from them by automatic differentiation, never written by hand.

PoroMechanics.storage! Function

storage!(f, u, node, model::AbstractPoroModel, data)

Fill f with the stored amounts for model at node — not their time derivative: the solver forms the accumulation term from them. Must be implemented by each concrete FVM model.

source

Storage term .

source

Storage term .

source

Storage term .

source

Stored water mass, dry-air mass and entropy, per unit volume of medium.

source

storage!(f, u, node, m::NernstPlanck, data)

φ cᵢ for each ion, and zero for the potential, which makes its row algebraic: the constraint in reaction! holds at every instant rather than accumulating.

source

PoroMechanics.flux! Function

flux!(f, u, edge, model::AbstractPoroModel, data)

Fill f with the inter-node flux terms for model at edge. Must be implemented by each concrete FVM model.

source

Two-point flux  .

VoronoiFVM divides f by the edge length, so the flux is written as a difference of node values rather than as a gradient.

source

Two-point Darcy flux  .

source

Two-point flux, with the conductivity evaluated at the mean capillary pressure of the edge.

VoronoiFVM divides f by the edge length, which is why the gravity term carries an explicit edge projection while the pressure term does not. That projection is taken along gravity_axis, so the same model drives a horizontal column and a vertical one.

Note edge.coord is the coordinate matrix of the whole grid, not of the edge's two nodes: the endpoints are edge.coord[:, edge.node[1]] and edge.coord[:, edge.node[2]]. Writing edge.coord[a, 2] - edge.coord[a, 1] instead reads global nodes 1 and 2 whatever edge is being assembled. On a uniform 1D grid those happen to be adjacent and equally spaced, so the mistake returns the right number and hides; in two dimensions it silently removes gravity and the column stops draining.

source

flux!(f, u, edge, m::NernstPlanck, data)

Scharfetter-Gummel for each ion. The potential carries no flux of its own.

VoronoiFVM.fbernoulli_pm(x) returns (B(x), B(-x)) with B(x) = x/(eˣ−1), and the fitted flux of −D(∇c + z c ∇Ψ) is D (B(−δ) c₁ − B(δ) c₂) with δ = z(Ψ₁ − Ψ₂). At δ = 0 both Bernoulli factors are one and this is the plain difference.

source

PoroMechanics.bcondition! Function

bcondition!(f, u, node, model::AbstractPoroModel, data)

Apply boundary conditions for model at boundary node. Default: no-flux (Neumann zero). Override for Dirichlet or non-zero Neumann.

source

Imposed concentrations, from the model's dirichlet field.

source

Imposed pressures, from the model's dirichlet field, resolved at the current time.

source

Imposed pressures, from the model's dirichlet field.

source

bcondition!(f, u, bnode, m::NernstPlanck, data)

One boundary tuple per unknown. The potential needs at least one Dirichlet value: only its gradient enters the fluxes, so without a reference the system is singular by a constant.

source

PoroMechanics.reaction! Function

reaction!(f, u, node, model::AbstractPoroModel, data)

Fill f with volumetric reaction (source/sink) terms and algebraic constraints for model at node.

Use cases:

• Chemical reactions between species (e.g. precipitation, equilibrium)

• Algebraic constraints (e.g. electroneutrality: f[iψ] = Σ zᵢ·cᵢ)

Default: no reaction. Override for models with coupled chemistry or constraints.

source

reaction!(f, u, node, m::NernstPlanck, data)

The electroneutrality constraint Σᵢ zᵢ cᵢ = q_background, which is what determines Ψ. Zero for the ions: they have no volumetric source here.

source

Finite element interface

Callbacks consumed by Ferrite.jl.

PoroMechanics.assemble_element! Function

assemble_element!(Ke, re, el, u_el, model::AbstractPoroModel, cv, Δt)

Assemble the element stiffness matrix Ke and residual vector re for FEM models. Must be implemented by each concrete Ferrite-based model.

source

PoroMechanics.element_matrices! Function

element_matrices!(ke1, ke2, cell, model::AbstractPoroModel, cv_u, cv_p)

Fill the stationary element matrix ke1 and the storage element matrix ke2 for linear coupled FEM models (Biot poroelasticity and similar).

The global system is solved at each time step as: A · xⁿ⁺¹ = f_ext + (1/Δt) · K2 · xⁿ, with A = K1 + K2/Δt

where K1 = Σ ke1 (elastic rigidity + Darcy conductivity) and K2 = Σ ke2 (Biot coupling + compressibility storage).

Must be implemented by each concrete Ferrite-based model that uses this split.

source

PoroMechanics.facet_load! Function

facet_load!(fe, facet, model::AbstractPoroModel, fv)

Accumulate the surface load vector fe for boundary facet facet. Default: no traction (zero Neumann). Override for pressure or traction loads.

source

Model introspection

PoroMechanics.nspecies Function

nspecies(model::AbstractPoroModel) -> Int

Return the number of primary unknowns (species) solved by model.

source

PoroMechanics.species_names Function

species_names(model::AbstractPoroModel) -> Vector{Symbol}

Return the symbolic names of the primary unknowns, e.g. [:p_l, :T].

source

Backends

Glue to the solver packages. The physics lives in the constitutive layer; these only wire it up.

Finite volumes

PoroMechanics.fvm_system Function

fvm_system(model, grid; species = 1:nspecies(model), reaction = false, kwargs...)

Build a VoronoiFVM.System for model on grid, wiring storage!, flux! and bcondition! to the model by dispatch.

reaction = true also wires reaction! — for models carrying volumetric source terms or algebraic constraints such as electroneutrality. It is off by default so that models without reactions do not pay for differentiating a no-op.

Any further keyword argument is forwarded untouched to VoronoiFVM.System.

sys = fvm_system(model, grid)
tsol = solve(sys; inival, times, control)

source

Finite elements

PoroMechanics.biot_element_matrices! Function

biot_element_matrices!(ke1, ke2, m::BiotPoroelastic, cv_u, cv_p)

Element contributions for plane-strain Biot poroelasticity, with cv_u a vector-valued CellValues for the displacement and cv_p a scalar one for the pressure.

source

PoroMechanics.radial_element_matrices! Function

radial_element_matrices!(ke1, ke2, m::BiotPoroelastic, cv_u, cv_p, coords; nhoop)

Element contributions for a radially symmetric Biot problem on a 1D mesh in r, where the displacement is the scalar .

nhoop is the number of hoop directions, and it is the only thing that separates the two geometries:

nhoop	geometry	strains	volume weight
1	long cylinder, plane strain	,  ,
2	sphere	,

The discrete strain operator carries nhoop copies of alongside , so contracting with the isotropic stiffness gives

The hoop terms are what no Cartesian element produces. The factors are why the quadrature points must stay strictly inside the elements — they do — and the weight suppresses what is left near the axis.

This is the kinematics an axisymmetric Barcelona Basic Model will need, written by hand.

source

PoroMechanics.node_dof_maps Function

node_dof_maps(dh, grid, fields...) -> NamedTuple

For each requested field, a vector giving the global dof of that field at each node.

Scalar fields get one dof per node. For a vector field, pass (:u, component) to select a component: Ferrite interleaves the components, so component c of local node loc sits at dof_range(dh, :u)[dim*(loc-1) + c].

source

PoroMechanics.combine! Function

combine!(A, K1, K2, inv_dt)

Fill A with K1 + K2/Δt, in place and without touching the sparsity pattern.

Writing K1 + inv_dt .* K2 instead would be wrong in the presence of affine constraints: sparse addition prunes entries that are numerically zero, and the master–slave couplings are exactly that until apply! condenses them — so the pattern would lose the very slots the condensation needs. All three matrices must come from the same allocate_matrix, and the arithmetic then goes straight through nzval.

source

Axisymmetric elastoplasticity

Ferrite v1 has no axisymmetric element, so the kinematics are written out: the hoop strain   makes the strain a genuine 3D tensor even though the mesh is 2D, and a constitutive model that reads gets the wrong mean stress without it.

PoroMechanics.axisymmetric_strain Function

axisymmetric_strain(cv, q, ue, coords) -> SymmetricTensor{2,3}

Strain at a quadrature point of an axisymmetric element, as a genuine 3D tensor.

The hoop strain   is not a bookkeeping detail: it is a real component, and a plasticity model that reads the mean stress as gets the wrong answer if it is missing. Working in 3D from the start is what stops that class of error.

Index order is , with and zero by symmetry.

source

PoroMechanics.axisymmetric_shape_strain Function

axisymmetric_shape_strain(cv, q, i, coords) -> SymmetricTensor{2,3}

The virtual strain of shape function i, in the same 3D form.

source

PoroMechanics.assemble_axisymmetric! Function

assemble_axisymmetric!(K, f, dh, cv, mat, states, states_old, u, Δt)

Assemble the tangent K and the internal force f of an axisymmetric mechanical problem, asking mat for its response at every quadrature point.

states_old holds the converged state of the previous step and is read; states is written with the trial state of the current iterate. Keeping them apart is what makes a Newton iteration repeatable: an iteration that overwrote the history could not be taken twice from the same starting point, and a line search or a rejected step would corrupt the material.

Returns the assembled pair; the volume element carries the axisymmetric weight .

source

PoroMechanics.newton_solve! Function

newton_solve!(u, K, f, dh, cv, mat, states, states_old, ch, Δt; tol, maxiter, maxhalve, linsolve, fext)

Newton–Raphson on the equilibrium residual, with backtracking, returning the residual norm of the initial state and every accepted correction. Returns only after convergence; throws if the iteration budget or the backtracking search is exhausted. On failure, u and states describe the last accepted iterate, not a rejected trial.

Why the backtracking is not optional. An exact tangent buys a quadratic rate near the solution; it says nothing about getting there. The step on which a material first yields is where undamped Newton fails: the increment is finite, the tangent switches from elastic to elastoplastic between one iterate and the next, and the full step overshoots. Measured on the Barcelona Basic Model under isotropic compression, the first plastic step diverges outright — the residual climbs from 4.8·10⁴ and wanders for twenty-five iterations without ever descending, while the well-conditioned tangent ( ) and the purely elastic steps that precede it converge in a single iteration. The failure is globalization, not linearization.

Halving the step until the residual decreases fixes it, and costs nothing where it is not needed: away from the transition the full step is accepted at once, so the quadratic rate survives intact. maxhalve bounds the number of trial steps; exhausting it raises an error without committing a trial. maxiter bounds accepted Newton corrections, and convergence is checked after the last permitted correction as well as at the initial state.

source