Reference

HydroTrixi

HydroTrixi.jl is an adaptive discontinuous spectral-element framework for hydrologic problems based on the Trixi.jl and SciML ecosystems.

examples_dir()

Return the absolute path to the package examples/ directory.

default_algorithm(semi)

Return a recommended OrdinaryDiffEq time-integration algorithm for the semidiscretization semi, suitable for passing to SciMLBase.solve.

The chosen algorithm depends on the type of semi:

• For a generic Trixi.AbstractSemidiscretization, returns Tsit5(), a 5(4) adaptive explicit Runge-Kutta method that is a sensible default for non-stiff problems.
• For a SemidiscretizationImplicit, returns Rodas5P(autodiff = AutoFiniteDiff()), a stiffly accurate linearly implicit Rosenbrock method.

compute_eoc(errors; refinement_factor=2.0)

Compute the estimated orders of convergence between successive entries of errors. The first entry of the returned vector is NaN because no coarser reference value is available.

HydrologicProblem(; equations, initial_condition, boundary_conditions, domain, tspan,
                    source_terms = nothing, state_to_evolved = nothing,
                    evolved_to_state = nothing)

A container for reusable hydrologic PDE problem data. It stores the governing equations, an initial_condition, boundary_conditions, the spatial domain, the time interval tspan, optional source_terms, and optional state-to-evolved and evolved-to-state maps for mixed-form problems.

The domain must be a pair (x_min, x_max) of coordinate tuples with matching dimension. The type parameter NDIMS records that dimension for dispatch and introspection.

RichardsEquation1D(; soil_model)

A one-dimensional Richards equation model written in terms of the pressure head $\psi$ for use with Trixi.jl's parabolic spatial discretization infrastructure. It represents the spatial operator in either the pressure-head form

\[\frac{\partial \theta}{\partial \psi} \frac{\partial \psi}{\partial t} = \frac{\partial}{\partial z} \left( K(\psi) \left( \frac{\partial \psi}{\partial z} - 1 \right) \right),\]

or the mixed form

\[\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z} \left( K(\psi) \left( \frac{\partial \psi}{\partial z} - 1 \right) \right),\]

where the water retention map between $\theta$ and $\psi$ is supplied by a separate semidiscretization wrapper. The hydraulic conductivity is supplied through soil_model, with hydraulic_conductivity(psi, equations) dispatching on the model type parameter SoilModel. If soil_model is omitted, it defaults to a Haverkamp model parameterized with the Celia (1990) reference values reported in Ireson et al. (2023), Eq. (25), in SI units (lengths in m, time in seconds).

water_content(u, equations::RichardsEquation1D)

Return the volumetric water content associated with the pressure head state u under the Richards equation model equations.

water_content_timederivative

Analysis integral for the time derivative of the total water content, $\mathrm{d} \int \theta \,\mathrm{d}x / \mathrm{d}t$.

For the constitutive mixed form, this integrates the evolved water-content derivative. For the pressure-head form, this integrates $C(\psi) \partial \psi / \partial t$ using water_capacity.

mass_bias
mass_bias(u_ode, semi::SemidiscretizationImplicit)
mass_bias(u_ode, semi::SemidiscretizationImplicit, initial_water_content)

Return the mass balance bias

\[\epsilon_B = \int_{t_0}^{t_M} \left(q(t, 0) - q(t, z_N)\right)\,\mathrm{d}t - \int_0^L \left(\theta(t_M, z) - \theta(t_0, z)\right)\,\mathrm{d}z.\]

This diagnostic requires PassiveVariablesBoundaryFlux1D, which stores the time-integrated boundary fluxes in the DG surface orientation. The initial total water content is recorded when the implicit coefficients are initialized and is used by the analysis callback.

mass_bias_history(sol; initial_water_content = nothing)
mass_bias_history(analysis_path::AbstractString;
                  time_column = "time", mass_bias_column = "mass_bias")

Return the saved times and corresponding mass_bias values for sol.

The solution must use a SemidiscretizationImplicit with PassiveVariablesBoundaryFlux1D. If no initial water content has been recorded for the semidiscretization and the first saved state is not at the initial time, pass the initial total water content explicitly with initial_water_content.

When analysis_path is provided, read the time and mass-bias columns from a Trixi analysis file written by AnalysisCallback(save_analysis = true).

water_capacity(u, equations::RichardsEquation1D)

Return the water capacity $C(\psi) = \mathrm{d}\theta / \mathrm{d}\psi$ associated with the pressure head state u under the Richards equation model equations.

pressure_head(u)
pressure_head(u, equations::RichardsEquation1D)

Return the pressure head $\psi$ stored in state u.

For scalar states, u is returned directly. For vector-like states, the first component is interpreted as pressure head.

pressure_head_from_water_content(theta, equations::RichardsEquation1D)

Return the pressure head associated with water content theta for the retention curve stored in equations.

Haverkamp(; saturated_hydraulic_conductivity, alpha, beta, A, gamma, theta_s, theta_r)

A Haverkamp constitutive model for Richards' equation, written in the form used by Celia et al. (1990) and Ireson et al. (2023). For $\psi < 0$, the effective saturation and hydraulic conductivity are

\[S_e(\psi) = \frac{\alpha}{\alpha + |\psi|^\beta}, \qquad K(\psi) = K_s\,\frac{A}{A + |\psi|^\gamma},\]

with saturated values recovered for $\psi \ge 0$. The parameters $\alpha$ and $A$ carry units of $[L]^\beta$ and $[L]^\gamma$ respectively, where $[L]$ is the length unit chosen by the caller for $\psi$. The associated water content is

\[\theta(\psi) = \theta_r + (\theta_s - \theta_r) S_e(\psi).\]

References

• Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P. J., Vachaud, G. (1977). A comparison of numerical simulation models for one-dimensional infiltration. Soil Science Society of America Journal, 41(2), 285-294. DOI: 10.2136/sssaj1977.03615995004100020024x
• Celia, M. A., Bouloutas, E. T., Zarba, R. L. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 26(7), 1483-1496. DOI: 10.1029/WR026i007p01483
• Ireson, A. M., Spiteri, R. J., Clark, M. P., Mathias, S. A. (2023). A simple, efficient, mass-conservative approach to solving Richards' equation (openRE, v1.0). Geoscientific Model Development, 16, 659-677. DOI: 10.5194/gmd-16-659-2023

VanGenuchten(; saturated_hydraulic_conductivity, alpha, n, theta_s, theta_r,
               m = 1 - 1 / n, pore_connectivity = 1 / 2)

A van Genuchten-Mualem soil-hydraulic model, written in terms of the effective saturation following Ireson et al. (2023). For $\psi < 0$ the effective saturation is

\[S_e(\psi) = \left(1 + (\alpha |\psi|)^n\right)^{-m},\]

with $S_e(\psi) = 1$ for $\psi \ge 0$. The hydraulic conductivity is

\[K(\psi) = K_s S_e(\psi)^l \left(1 - \left(1 - S_e(\psi)^{1 / m}\right)^m\right)^2,\]

where $K_s$ is the saturated_hydraulic_conductivity parameter and $l$ is the pore_connectivity parameter. The volumetric water content $\theta(\psi) = \theta_r + (\theta_s - \theta_r) S_e(\psi)$ follows from the generic water_content relation using the residual and saturated water contents theta_r and theta_s.

References

• van Genuchten, M. Th. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5), 892-898. DOI: 10.2136/sssaj1980.03615995004400050002x
• Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3), 513-522. DOI: 10.1029/WR012i003p00513

HydrologicProblemCelia1990(; tspan = (0.0, 360.0))

Return the one-dimensional Richards-equation infiltration problem introduced by Celia, Bouloutas, and Zarba (1990) as a HydrologicProblem for HydroTrixi semidiscretizations.

The problem uses the Haverkamp constitutive laws in the form reported by Ireson et al. (2023), Eq. (25):

\[\theta(\psi) = \theta_r + (\theta_s - \theta_r) \frac{\alpha}{\alpha + |\psi|^\beta}, \qquad K(\psi) = K_s\,\frac{A}{A + |\psi|^\gamma},\]

for $\psi < 0$, with saturated values $\theta = \theta_s$ and $K = K_s$ for $\psi \ge 0$. The setup uses depth $z$ in metres, positive downward on $z \in [0, 0.4]$, and time in seconds on $t \in [0, 360]$. The pressure head is initialized at $-0.615$ m, with a Dirichlet boundary value of $-0.207$ m at the top boundary (x_neg) and $-0.615$ m at the bottom boundary (x_pos).

The returned problem setup contains the fields equations, state_to_evolved, evolved_to_state, initial_condition, boundary_conditions, domain, and tspan.

References

• Celia, M. A., Bouloutas, E. T., Zarba, R. L. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 26(7), 1483-1496. DOI: 10.1029/WR026i007p01483
• Ireson, A. M., Spiteri, R. J., Clark, M. P., Mathias, S. A. (2023). A simple, efficient, mass-conservative approach to solving Richards' equation (openRE, v1.0). Geoscientific Model Development, 16, 659-677. DOI: 10.5194/gmd-16-659-2023

HydrologicProblemRichardsManufacturedSolution(; tspan = (0.0, 120.0),
                                                soil_model = default_soil_model())

Return a one-dimensional Richards equation manufactured-solution problem. The pressure head is given by richards_manufactured_solution, the boundary conditions are Dirichlet values from that profile, and the source term is source_terms_richards_manufactured_solution. The default setup uses the same Haverkamp parameters as HydrologicProblemCelia1990.

The problem uses depth $z$ in metres on $z \in [0, 0.2]$ and time in seconds on $t \in [0, 120]$ by default. It is intended for regression and convergence checks of mixed and pressure-head forms of the Richards equation.

References

• List, F., Radu, F. A. (2016). A study on iterative methods for solving the Richards equation. Computational Geosciences, 20, 341-353. DOI: 10.1007/s10596-016-9566-3
• Keita, S., Beljadid, A., Bourgault, Y. (2021). Implicit and semi-implicit second-order time stepping methods for the Richards equation. arXiv:2105.05224

richards_manufactured_solution(x, t, equations)

Return the manufactured pressure head profile used by HydrologicProblemRichardsManufacturedSolution. In SI units, the pressure head is

\[\psi(z, t) = 0.204 \tanh\left(\frac{1}{2}\left(100z + \frac{t}{12} - 15\right)\right) - 0.411.\]

source_terms_richards_manufactured_solution(u, gradients, x, t, equations)

Return the storage source term that makes richards_manufactured_solution solve the one-dimensional Richards equation with Haverkamp constitutive laws,

\[s = C(\psi)\frac{\partial \psi}{\partial t} - \frac{\partial}{\partial z} \left(K(\psi)\left(\frac{\partial \psi}{\partial z} - 1\right)\right).\]

HydrologicProblemRichardsClosedColumn(; soil_model, domain, tspan, base_head, amplitude)

Return a one-dimensional Richards-equation redistribution problem on a closed column with homogeneous Neumann boundary conditions at both ends. The setup is intended for mass-conservation studies because the net boundary flux vanishes, so the total water content should remain constant up to time-integration and nonlinear-solve tolerances.

The problem uses the mixed Richards formulation with state-to-evolved map $\theta = \Theta(\psi)$. The pressure head is initialized as

\[\psi(z, 0) = \psi_{\text{base}} + z + A \left( 1 - \cos\left(\frac{2\pi (z - z_{\min})}{L}\right) \right),\]

where $\psi_{\text{base}}$ is base_head, $A$ is amplitude, and $L = z_{\max} - z_{\min}$. This profile satisfies $\partial \psi / \partial z = 1$ at both boundaries, and therefore the flux

\[q = K(\psi) \left( \frac{\partial \psi}{\partial z} - 1 \right)\]

vanishes at $z = z_{\min}$ and $z = z_{\max}$ at the initial time, and should remain zero for all time due to the homogeneous Neumann boundary conditions.

The returned problem setup contains the fields equations, state_to_evolved, evolved_to_state, initial_condition, boundary_conditions, domain, and tspan.

BoundaryConditionDirichletPenalty(boundary_value_function; penalty=1.0)

Dirichlet boundary condition with a configurable boundary penalty of the form $\tau (u_{\text{inner}} - u_b)$ added to the boundary flux on the divergence pass for parabolic diffusion operators. The boundary_value_function is called as boundary_value_function(x, t, equations), and penalty may be either a scalar or a function called as penalty(x, t, equations).

References

• Arnold, D. N. (1982). An interior penalty finite element method with discontinuous elements. SIAM Journal on Numerical Analysis, 19(4), 742-760. DOI: 10.1137/0719052
• Arnold, D. N., Brezzi, F., Cockburn, B., Marini, L. D. (2002). Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5), 1749-1779. DOI: 10.1137/S0036142901384162

MixedForm()

Select the mixed form of the Richards equation in SemidiscretizationImplicit. The evolved variable is water content, and pressure head is recovered through the water retention relation.

PressureHeadForm()
PressureHeadForm(; transfer_variables = pressure_head)

Select the pressure-head form of the Richards equation in SemidiscretizationImplicit. The evolved state is pressure head, and the temporal operator divides the spatial residual by the nodal water capacity $C(\psi)$. The optional transfer_variables map controls which variable is transferred during adaptive mesh refinement. The default transfers pressure head directly. Use transfer_variables = water_content to transfer water content and reconstruct pressure head after each mesh update.

SemidiscretizationImplicit{Semidiscretization, TemporalOperator, PassiveVariables}

A semidiscretization wrapper that augments a spatial semidiscretization with a temporal operator, so that the semi-discrete problem is not restricted to the explicit form $\mathrm{d}\mathbf{u}/\mathrm{d}t = R(\mathbf{u}, t)$. For example, TemporalOperatorConstitutive enforces a map between state and evolved variables, while TemporalOperatorCapacity applies a nodal capacity function to the spatial residual. Passive scalar variables can be appended to the ODE state for diagnostics that are integrated by the time integrator but do not affect the physical residual.

TemporalOperatorConstitutive(state_to_evolved; evolved_to_state = nothing)

Temporal operator for a SemidiscretizationImplicit that takes the form

\[\frac{\mathrm{d}\mathbf{u}_\mathrm{evolved}}{\mathrm{d} t} = R(\mathbf{u}_\mathrm{state}, t), \qquad 0 = \mathbf{u}_\mathrm{evolved} - S(\mathbf{u}_\mathrm{state}),\]

where $R$ is the operator associated with the spatial discretization semi_base, and $S$ is state_to_evolved, which maps the state variable(s) to the evolved variable(s). The optional evolved_to_state inverse is required by AMR to reconstruct the algebraic state after transferring the evolved block. In the case of the Richards equation, for example, the state variable is the pressure head, and the evolved variable is the water content, which are related by the water retention curve (i.e., water_content). HydroTrixi.jl then uses SciML's DAE solvers to integrate the resulting system of equations based on the following mass-matrix form: ```math \begin{bmatrix} I & 0 \ 0 & 0 \end{bmatrix} \frac{\mathrm{d}}{\mathrm{d} t} \begin{bmatrix} \mathbf{u}\mathrm{evolved} \ \mathbf{u}\mathrm{state} \end{bmatrix} = \begin{bmatrix} R(\mathbf{u}\mathrm{state}, t) \ \mathbf{u}\mathrm{evolved} - S(\mathbf{u}_\mathrm{state}) \end{bmatrix}.

TemporalOperatorCapacity(capacity_function; transfer_variables, transfer_to_state)

Temporal operator for a SemidiscretizationImplicit that stores the state variable directly and applies a nodal capacity function to the spatial residual,

\[\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = C(\mathbf{u})^{-1} R(\mathbf{u}, t).\]

For the pressure-head form of the Richards equation, the state variable is $\psi$ and $C(\psi) = \mathrm{d}\theta / \mathrm{d}\psi$. The capacity must be strictly positive at all nodal states. The optional AMR transfer maps convert the state to the transferred variable before mesh adaptation and reconstruct the state afterwards.

NoPassiveVariables()

Passive-variable configuration for implicit semidiscretizations without appended diagnostic variables.

PassiveVariables(nvariables, initial_condition, rhs!)

Append nvariables passive scalar variables to a SemidiscretizationImplicit. The passive variables are integrated by the time integrator but do not affect the physical residual. The function initial_condition(semi_base, t) returns the initial passive values, and rhs!(dq, u_physical, du_physical, semi, t) fills their time derivatives.

PassiveVariablesBoundaryFlux1D()

Append two passive scalar variables that store the time-integrated numerical boundary fluxes at the negative and positive boundaries of a one-dimensional scalar problem. The stored fluxes use the same sign convention as the parabolic flux in the semidiscrete operator.

passive_variable_view(u_ode, semi::SemidiscretizationImplicit)

Return a view of the passive scalar variables appended to the ODE state of semi.

passive_variables(u_ode, semi::SemidiscretizationImplicit)

Return a copy of the passive scalar variables appended to the ODE state of semi.

boundary_flux_integrals(u_ode, semi::SemidiscretizationImplicit)

Return the integrated negative- and positive-boundary fluxes stored by PassiveVariablesBoundaryFlux1D as a named tuple with fields x_neg and x_pos.

plot_solution_1d(sol; kwargs...)

Plot the final one-dimensional solution profile stored in sol, save it to output_path, and return the CairoMakie.Figure.

This method is provided by HydroTrixiVisualizationExt and becomes available when CairoMakie and LaTeXStrings are loaded.

animate_solution_1d(sol; kwargs...)
animate_solution_1d(ode::SciMLBase.ODEProblem; callback, kwargs...)

Animate a one-dimensional solution and save it to output_path.

The first method renders frames from saved states in sol. The second method advances an ODEProblem to the requested frame times while applying callback, which is useful for animations with adaptive meshes.

These methods are provided by HydroTrixiVisualizationExt and become available when CairoMakie and LaTeXStrings are loaded.

plot_convergence_1d(ndofs, l2_errors, linf_errors; kwargs...)

Plot one-dimensional convergence data against the number of degrees of freedom, save the figure to output_path, and return the CairoMakie.Figure.

This method is provided by HydroTrixiVisualizationExt and becomes available when CairoMakie and LaTeXStrings are loaded.

plot_mass_bias(sol; kwargs...)
plot_mass_bias(analysis_path::AbstractString;
               time_column = "time", mass_bias_column = "mass_bias", kwargs...)
plot_mass_bias(sources; labels, kwargs...)

Plot the saved mass bias time history stored in sol, in a Trixi analysis file, or in multiple sources, save it to output_path, and return the CairoMakie.Figure. By default, the mass bias is scaled by an automatically selected power of ten and the corresponding exponent is shown above the top-left corner of the axis. Pass absolute = true, yscale = log10 to plot $|\epsilon_b|$ on a logarithmic axis.

This method is provided by HydroTrixiVisualizationExt and becomes available when CairoMakie and LaTeXStrings are loaded.