Idaho National Laboratory

www.inl.gov

MOOSE Introduction

Multi-physics Object Oriented Simulation Environment

Governing Equations

Conservation of Mass:

\nabla \cdot \vec{u} = 0

Conservation of Energy:

C\left( \frac{\partial T}{\partial t} + \epsilon \vec{u}\cdot\nabla T \right) - \nabla\cdot k \nabla T = 0

Darcy's Law:

\vec{u} = -\frac{\mathbf{K}}{\mu} (\nabla p - \rho \vec{g})

where $\vec{u}$ is the fluid velocity, $\epsilon$ is porosity, $\mathbf{K}$ is the permeability tensor, $\mu$ is fluid viscosity, $p$ is the pressure, $\rho$ is the density, $\c_p$ is the specific heat, $\vec{g}$ is the gravity vector, and $T$ is the temperature.

Assuming that $\vec{g}=0$ and imposing the divergence-free condition of Eq. (1) to Eq. (3) leads to the following system of two equations in the unknowns $p$ and $T$ :

-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p = 0

\rho c_p \left( \frac{\partial T}{\partial t} + \epsilon \vec{u}\cdot\nabla T \right) - \nabla \cdot k \nabla T = 0

The parameters $\rho$ , $c_p$ , and $k$ are the porosity-dependent density, specific heat capacity, and thermal conductivity of the combined fluid/solid medium, defined by:

\rho \equiv \epsilon \rho_f + (1-\epsilon) \rho_s \\ \rho c_p \equiv \epsilon \rho_f {c_p}_f + (1-\epsilon) \rho_s {c_p}_s \\ k \equiv \epsilon k_f + (1-\epsilon) k_s

where $\epsilon$ is the porosity, $c_p$ is the specific heat, and the subscripts $f$ and $s$ refer to fluid and solid, respectively.

Material Properties

Property	Value	Units
Viscosity of water, $\mu_f$	$7.98\times10^{-4}$	$\textrm{P}\cdot\textrm{s}$
Density of water, $\rho_f$	995.7	$\textrm{kg}/\textrm{m}^3$
Density of steel, $\rho_s$	8000	$\textrm{kg}/\textrm{m}^3$
Thermal conductivity of water, $k_f$	0.6	$\textrm{W}/\textrm{m}\,\textrm{K}$
Thermal conductivity of steel, $k_s$	18	$\textrm{W}/\textrm{m}\,\textrm{K}$
Specific heat capacity of water, $c_p{_f}$	4181.3	$\textrm{J}/(\textrm{kg}\,\textrm{K})$
Specific heat capacity of steel, $c_p{_s}$	466	$\textrm{J}/(\textrm{kg}\,\textrm{K})$

Step 6: Equation Coupling

Solve the pressure and temperature in a coupled system of equations by adding the advection term to the heat equation.

-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p = 0 \\ \rho c_p \left( \frac{\partial T}{\partial t} + \epsilon \vec{u}\cdot\nabla T \right) - \nabla \cdot k \nabla T = 0

First, consider the steady-state diffusion equation on the domain $\Omega$ : find $u$ such that

-\nabla \cdot \nabla u = 0 \in \Omega,

where $u = 4000$ on the left, $u = 0$ on the right and with $\nabla u \cdot \hat{n} = 0$ on the remaining boundaries.

The weak form of this equation, in inner-product notation, is given by:

\left(\nabla \psi_i, \nabla u_h\right) = 0 \quad \forall \psi_i,

where $\psi_i$ are the test functions and $u_h$ is the finite element solution.

Polynomial Fitting

To introduce the concept of FEM, consider a polynomial fitting exercise. When fitting a polynomial there is a known set of points as well as a set of coefficients that are unkown for a function that has the form:

f(x) = a + bx + cx^2 + \dots,

where $a$ , $b$ and $c$ are scalar coefficients and $1$ , $x$ , $x^2$ are "basis functions". Thus, the problem is to find $a$ , $b$ , $c$ , etc. such that $f(x)$ passes through the points given.

More generally,

f(x) = \sum_{i=0}^d c_i x^i,

where the $c_i$ are coefficients to be determined. $f(x)$ is unique and interpolary if $d+1$ is the same as the number of points needed to fit. This defines a linear system that must be solved to find the coefficients.

Polynomial Example

Define a set of points:

(x_1, y_1) = (1,4) \\ (x_2, y_2) = (3,1) \\ (x_3, y_3) = (4,2)

Substitute $(x_i, y_i)$ data into the model:

y_i = a + bx_i + cx_i^2, i=1,2,3.

This leads to the following linear system for $a$ , $b$ , and $c$ :

\begin{bmatrix} 1 & 1 & 1 \\ 1 & 3 & 9 \\ 1 & 4 & 16 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} 4 \\ 1 \\ 2 \end{bmatrix}

Solving for the coefficients results in:

\begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} 8 \\ \frac{29}{6} \\ \frac{5}{6} \end{bmatrix}

These coefficients define the solution function:

f(x) = 8 - \frac{29}{6} x + \frac{5}{6} x^2

The solution is the function, not the coefficients.

The coefficients are meaningless, they are just numbers used to define a function.

The solution is not the coefficients, but rather the function created when they are multiplied by their respective basis functions and summed.

The function $f(x)$ does go through the points given, but it is also defined everywhere in between.

$f(x)$ can be evaluated at the point $x=2$ , for example, by computing:

f(2) = \sum_{i=0}^2 c_i 2^i = \sum_{i=0}^2 c_i g_i(2),

where the $c_i$ correspond to the coefficients in the solution vector, and the $g_i$ are the respective functions.

FEM can be used to solve both linear and nonlinear PDEs

FEM is a method for numerically approximating the solution to partial differential equations (PDEs). FEM is widely applicable for a large range of PDEs and domains.

Example PDEs: Have you seen them before? Are they linear/nonlinear? Coupled?

-\nabla \cdot \nabla u = q

\dfrac{\partial u_x}{\partial t} -\nabla\cdot \mu \nabla u_x + \vec{u} \cdot \nabla u_x = 0

C \dfrac{\partial T}{\partial t} -\nabla \cdot k(T) \nabla T = q(\psi) \\ \dfrac{1}{v} \dfrac{\partial \psi}{\partial t} + \Omega \cdot \nabla \psi + \Sigma(T) \psi = Q(\psi, T)

Weak Form

Using FEM to find the solution to a PDE starts with forming a "weighted residual" or "variational statement" or "weak form", this processes if referred to here as generating a weak form.

The weak form provides flexibility, both mathematically and numerically and it is needed by MOOSE to solve a problem.

Generating a weak form generally involves these steps:

Write down strong form of PDE.
Rearrange terms so that zero is on the right of the equals sign.
Multiply the whole equation by a "test" function $\psi$ .
Integrate the whole equation over the domain $\Omega$ .
Integrate by parts and use the divergence theorem to get the desired derivative order on your functions and simultaneously generate boundary integrals.

note: Exercise

Obtain the weak form for the equations listed on the previous slide and the shape functions.

Integration by Parts and Divergence Theorem

Suppose $\varphi$ is a scalar function, $\vec{v}$ is a vector function, and both are continuously differentiable functions, then the product rule states:

\nabla\cdot(\varphi\vec{v}) = \varphi(\nabla\cdot\vec{v}) + \vec{v}\cdot(\nabla \varphi)

The function can be integrated over the volume $V$ and rearranged as:

\int \varphi(\nabla\cdot\vec{v})\,dV = \int \nabla\cdot(\varphi\vec{v})\,dV - \int \vec{v}\cdot(\nabla\varphi)\,dV

The divergence theorem transforms a volume integral into a surface integral on surface $s$ :

\int \nabla \cdot (\varphi\vec{v})\,\text{d}V = \int \varphi\vec{v}\cdot\hat{n}\,\text{d}s,

where $\hat{n}$ is the outward normal vector for surface $s$ . Combining Eq. (7) and Eq. (8) yield:

\boxed{\int \varphi (\nabla\cdot\vec{v})\,dV = \int \varphi\vec{v}\cdot\hat{n}\,ds - \int\vec{v}\cdot(\nabla\varphi)\,dV}

Example: Advection-Diffusion

(1) Write the strong form of the equation:

-\nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u = f

(2) Rearrange to get zero on the right-hand side:

-\nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u - f = 0

(3) Multiply by the test function $\psi$ :

-\psi \left(\nabla\cdot k\nabla u\right) + \psi\left(\vec{\beta} \cdot \nabla u\right) - \psi f = 0

(4) Integrate over the domain $\Omega$ :

{-\int_\Omega\psi \left(\nabla\cdot k\nabla u\right)} + \int_\Omega\psi\left(\vec{\beta} \cdot \nabla u\right) - \displaystyle\int_\Omega\psi f = 0

(5) Integrate by parts and apply the divergence theorem, by using Eq. (9) on the left-most term of the PDE:

\int_\Omega\nabla\psi\cdot k\nabla u - \int_{\partial\Omega} \psi \left(k\nabla u \cdot \hat{n}\right) + \int_\Omega\psi\left(\vec{\beta} \cdot \nabla u\right) - \int_\Omega\psi f = 0

Write in inner product notation. Each term of the equation will inherit from an existing MOOSE type as shown below.

\underbrace{\left(\nabla\psi, k\nabla u \right)}_{Kernel} - \underbrace{\langle\psi, k\nabla u\cdot \hat{n} \rangle}_{BoundaryCondition} + \underbrace{\left(\psi, \vec{\beta} \cdot \nabla u\right)}_{Kernel} - \underbrace{\left(\psi, f\right)}_{Kernel} = 0

Corresponding MOOSE input file blocks

\underbrace{\left(\nabla\psi, k\nabla u \right)}_{Kernel} - \underbrace{\langle\psi, k\nabla u\cdot \hat{n} \rangle}_{BoundaryCondition} + \underbrace{\left(\psi, \vec{\beta} \cdot \nabla u\right)}_{Kernel} - \underbrace{\left(\psi, f\right)}_{Kernel} = 0

[Kernels]
  [diff]
    type = Diffusion
    variable = u
  []
[]

(test/tests/kernels/2d_diffusion/2d_diffusion_neumannbc_test.i)

$\quad$

[BCs]
  [right]
    type = NeumannBC
    variable = u
    boundary = 1
    value = 1
  []
[]

(test/tests/kernels/2d_diffusion/2d_diffusion_neumannbc_test.i)

$\quad$

[Kernels]
  [adv_u]
    implicit = false
    type = ConservativeAdvection
    variable = u
    velocity = '1 0 0'
  []
[]

(test/tests/dgkernels/1d_advection_dg/1d_advection_dg.i)

$\quad$

[Kernels]
  [ffn]
    type = BodyForce
    variable = u
    function = f_fn
  []
[]

(test/tests/bcs/nodal_normals/circle_tris.i)

Shape Functions

The discretized expansion of $u$ takes on the following form:

u \approx u_h = \sum_{j=1}^N u_j \phi_j,

where $\phi_j$ are the "basis functions", which form the basis for the "trial function", $u_h$ . $N$ is the total number of functions for the discretized domain.

The gradient of $u$ can be expanded similarly:

\nabla u \approx \nabla u_h = \sum_{j=1}^N u_j \nabla \phi_j

In the Galerkin finite element method, the same basis functions are used for both the trial and test functions:

\psi = \{\phi_i\}_{i=1}^N

Substituting these expansions back into the example weak form (Eq. (10)) yields:

\left(\nabla\psi_i, k\nabla u_h \right) - \langle\psi_i, k\nabla u_h\cdot \hat{n} \rangle + \left(\psi_i, \vec{\beta} \cdot \nabla u_h\right) - \left(\psi_i, f\right) = 0, \quad i=1,\ldots,N

The left-hand side of the equation above is referred to as the $i^{th}$ component of the "residual vector," $\vec{R}_i(u_h)$ .

Shape Functions are the functions that get multiplied by coefficients and summed to form the solution.

Individual shape functions are restrictions of the global basis functions to individual elements.

They are analogous to the $x^n$ functions from polynomial fitting (in fact, you can use those as shape functions).

Typical shape function families: Lagrange, Hermite, Hierarchic, Monomial, Clough-Toucher

Lagrange shape functions are the most common, which are interpolatory at the nodes, i.e., the coefficients correspond to the values of the functions at the nodes.

2D Lagrange Shape Functions

Example bi-quadratic basis functions defined on the Quad9 element:

$\psi_0$ is associated to a "corner" node, it is zero on the opposite edges.
$\psi_4$ is associated to a "mid-edge" node, it is zero on all other edges.
$\psi_8$ is associated to the "center" node, it is symmetric and $\geq 0$ on the element.

$\psi_0$

$\psi_4$

$\psi_8$

Numerical Integration

The only remaining non-discretized parts of the weak form are the integrals. First, split the domain integral into a sum of integrals over elements:

\int_{\Omega} f(\vec{x}) \;\text{d}\vec{x} = \sum_e \int_{\Omega_e} f(\vec{x}) \;\text{d}\vec{x}

Through a change of variables, the element integrals are mapped to integrals over the "reference" elements $\hat{\Omega}_e$ .

\sum_e \int_{\Omega_e} f(\vec{x}) \;\text{d}\vec{x} = \sum_e \int_{\hat{\Omega}_e} f(\vec{\xi}) \left|\mathcal{J}_e\right| \;\text{d}\vec{\xi},

where $\mathcal{J}_e$ is the Jacobian of the map from the physical element to the reference element.

Quadrature

Quadrature, typically "Gaussian quadrature", is used to approximate the reference element integrals numerically.

\int_{\Omega} f(\vec{x})\;d\vec{x} \approx \sum_q w_q f(\vec{x}_q),

where $w_q$ is the weight function at quadrature point $q$ .

Under certain common situations, the quadrature approximation is exact. For example, in 1 dimension, Gaussian Quadrature can exactly integrate polynomials of order $2n-1$ with $n$ quadrature points.

Applying the quadrature to Eq. (12) we can simply compute:

\sum_e \int_{\hat{\Omega}_e} f(\vec{\xi}) \left|\mathcal{J}_e\right| \;\text{d}\vec{\xi} \approx \sum_e \sum_{q} w_{q} f( \vec{x}_{q}) \left|\mathcal{J}_e(\vec{x}_{q})\right|,

where $\vec{x}_{q}$ is the spatial location of the $q^\textrm{th}$ quadrature point and $w_{q}$ is its associated weight.

MOOSE handles multiplication by the Jacobian ( $\mathcal{J}_e$ ) and the weight ( $w_{q}$ ) automatically, thus your Kernel object is only responsible for computing the $f(\vec{x}_{q})$ part of the integrand.

Sampling $u_h$ at the quadrature points yields:

\begin{aligned} u(\vec{x}_{q}) &\approx u_h(\vec{x}_{q}) = \sum u_j \phi_j(\vec{x}_{q}) \\ \nabla u (\vec{x}_{q}) &\approx \nabla u_h(\vec{x}_{q}) = \sum u_j \nabla \phi_j(\vec{x}_{q}) \end{aligned}

Thus, the weak form of Eq. (11) becomes:

\begin{aligned} \vec{R}_i(u_h) &= \sum_e \sum_{q} w_{q} \left|\mathcal{J}_e\right|\underbrace{\left[ \nabla\psi_i\cdot k \nabla u_h + \psi_i \left(\vec\beta\cdot \nabla u_h \right) - \psi_i f \right](\vec{x}_{q})}_{\textrm{Kernel Object(s)}} \\ &- \sum_f \sum_{q_{face}} w_{q_{face}} \left|\mathcal{J}_f\right|\underbrace{\left[\psi_i k \nabla u_h \cdot \vec{n} \right](\vec x_{q_{face}})}_{\textrm{BoundaryCondition Object(s)}} \end{aligned}

The second sum is over boundary faces, $f$ . MOOSE Kernel or BoundaryCondition objects provide each of the terms in square brackets (evaluated at $\vec{x}_{q}$ or $\vec x_{q_{face}}$ as necessary), respectively.

Newton's Method

Newton's method is a "root finding" method with good convergence properties, in "update form", for finding roots of a scalar equation it is defined as: $f(x)=0$ , $f(x): \mathbb{R} \rightarrow \mathbb{R}$ is given by

\begin{aligned} f'(x_n) \delta x_{n+1} &= -f(x_n) \\ x_{n+1} &= x_n + \delta x_{n+1} \end{aligned}

Newton's Method in MOOSE

The residual, $\vec{R}_i(u_h)$ , as defined by Eq. (13) is a nonlinear system of equations,

\vec{R}_i(u_h)=0, \qquad i=1,\ldots, N,

that is used to solve for the coefficients $u_j, j=1,\ldots,N$ .

For this system of nonlinear equations Newton's method is defined as:

\begin{aligned} \mathbf{J}(\vec{u}_n) \delta\vec{u}_{n+1} &= -\vec{R}(\vec{u}_n) \\ \vec{u}_{n+1} &= \vec{u}_n + \delta\vec{u}_{n+1} \end{aligned}

where $\mathbf{J}(\vec{u}_n)$ is the Jacobian matrix evaluated at the current iterate:

J_{ij}(\vec{u}_n) = \frac{\partial \vec{R}_i(\vec{u}_n)}{\partial u_j}

JFNK

Uses a Krylov subspace-based linear solver

\begin{aligned} K_j = \{\mathbf{r}, \mathbf{J}\mathbf{r},\mathbf{J}^2\mathbf{r},...,\mathbf{J}^{j-1}\mathbf{r}\} \end{aligned}

The action of the Jacobian is approximated by:

\begin{aligned} \mathbf{J}(\mathbf{u})\mathbf{v} \approx \frac{\mathbf{R}(\mathbf{u}+\epsilon \mathbf{v})-\mathbf{R}(\mathbf{u})}{\epsilon} \end{aligned}

The Kernel method computeQpResidual is called to compute $\vec{R}(\vec{u}_n)$ during the nonlinear step (Eq. (14)).

During each linear step of JFNK, the computeQpResidual method is called to approximate the action of the Jacobian on the Krylov vector.

PJFNK

The action of the preconditioned Jacobian is approximated by:

\begin{aligned} \mathbf{JP}^{-1}\mathbf{v} \approx \frac{\mathbf{R}(\mathbf{u}+\epsilon \mathbf{P}^{-1}\mathbf{v})-\mathbf{R}(\mathbf{u})}{\epsilon} \end{aligned}

The Kernel method computeQpResidual is called to compute $\vec{R}(\vec{u}_n)$ during the nonlinear step (Eq. (14)).

During each linear step of PJFNK, the computeQpResidual method is called to approximate the action of the Jacobian on the Krylov vector. The computeQpJacobian and computeQpOffDiagJacobian methods are used to compute values for the preconditioning matrix.

NEWTON

The Kernel method computeQpResidual is called to compute $\vec{R}(\vec{u}_n)$ during the nonlinear step (Eq. (14)).

The computeQpJacobian and computeQpOffDiagJacobian methods are used to compute the preconditioning matrix. It is assumed that the preconditioning matrix is the Jacobian matrix, thus the residual and Jacobian calculations are able to remain constant during linear iterations.

Relies on two techniques:

chain rule

\dfrac{\partial f(g(x))}{\partial x} = \dfrac{\partial f(g(x))}{\partial g(x)} \dfrac{\partial g(x)}{\partial x}

operator overloading

\bold{x} = (x, \dfrac{\partial x}{\partial x}, \dfrac{\partial x}{\partial y}) = (x, 1, 0) \\ \bold{y} = (y, \dfrac{\partial y}{\partial x}, \dfrac{\partial y}{\partial y}) = (y, 0, 1) \\ x + y = (x + y, 1, 1) \\ x * y = (x * y, y, x)

One thing to note is that the derivatives are with regards to the degrees of freedom, but the residual is computed at quadrature points! There are therefore often several non-zero coefficients even for simply $\_u[\_qp]$ (value of variable u at a quadrature point).

FEM Derivative Identities

The following relationships are useful when computing Jacobian terms.

\frac{\partial u_h}{\partial u_j} = \sum_k\frac{\partial }{\partial u_j}\left(u_k \phi_k\right) = \phi_j

\frac{\partial \left(\nabla u_h\right)}{\partial u_j} = \sum_k \frac{\partial }{\partial u_j}\left(u_k \nabla \phi_k\right) = \nabla \phi_j

Newton for a Simple Equation

Again, consider the advection-diffusion equation with nonlinear $k$ , $\vec{\beta}$ , and $f$ :

\begin{aligned} - \nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u = f \end{aligned}

Thus, the $i^{th}$ component of the residual vector is:

\begin{aligned} \vec{R}_i(u_h) = \left(\nabla\psi_i, k\nabla u_h \right) - \langle\psi_i, k\nabla u_h\cdot \hat{n} \rangle + \left(\psi_i, \vec{\beta} \cdot \nabla u_h\right) - \left(\psi_i, f\right) \end{aligned}

\begin{aligned} \vec{R}_i(u_h) = \left(\nabla\psi_i, k\nabla u_h \right) - \langle\psi_i, k\nabla u_h\cdot \hat{n} \rangle + \left(\psi_i, \vec{\beta} \cdot \nabla u_h\right) - \left(\psi_i, f\right) \end{aligned}

Using the previously-defined rules in Eq. (18) and Eq. (19) for $\frac{\partial u_h}{\partial u_j}$ and $\frac{\partial \left(\nabla u_h\right)}{\partial u_j}$ , the $(i,j)$ entry of the Jacobian is then:

\begin{aligned} J_{ij}(u_h) &= \left(\nabla\psi_i, \frac{\partial k}{\partial u_j}\nabla u_h \right) + \left(\nabla\psi_i, k \nabla \phi_j \right) - \left \langle\psi_i, \frac{\partial k}{\partial u_j}\nabla u_h\cdot \hat{n} \right\rangle \\&- \left \langle\psi_i, k\nabla \phi_j\cdot \hat{n} \right\rangle + \left(\psi_i, \frac{\partial \vec{\beta}}{\partial u_j} \cdot\nabla u_h\right) + \left(\psi_i, \vec{\beta} \cdot \nabla \phi_j\right) - \left(\psi_i, \frac{\partial f}{\partial u_j}\right) \end{aligned}

That even for this "simple" equation, the Jacobian entries are nontrivial: they depend on the partial derivatives of $k$ , $\vec{\beta}$ , and $f$ , which may be difficult or time-consuming to compute analytically.

In a multiphysics setting with many coupled equations and complicated material properties, the Jacobian might be extremely difficult to determine.

Intrinsic Type	Variant(s)
bool	\|
char	unsigned
int	unsigned, long, short
float	\|
double	long
void	\|

Purpose	Symbols
Math	`+ - * / % += -= /= %= ++ --`
Comparison	`< > <= >= != ==`
Logical Comparison	`&& \|\| !`
Memory	`* & new delete sizeof`
Assignment	`=`
Member Access	`-> .`
Name Resolution	`::`

Kernel Object

To implement the Darcy pressure equation, a Kernel object is needed to add the coefficient to diffusion equation.

-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p = 0,

where $\textbf{K}$ is the permeability tensor and $\mu$ is the fluid viscosity.

A Kernel is C++ class, which inherits from MooseObject that is used by MOOSE for coding volume integrals of a partial differential equation (PDE).

Kernel Object Members

_u, _grad_u
Value and gradient of the variable this Kernel is operating on

_test, _grad_test
Value ( $\psi$ ) and gradient ( $\nabla \psi$ ) of the test functions at the quadrature points

_phi, _grad_phi
Value ( $\phi$ ) and gradient ( $\nabla \phi$ ) of the trial functions at the quadrature points

_q_point
Coordinates of the current quadrature point

_i, _j
Current index for test and trial functions, respectively

_qp
Current quadrature point index

Base	Override	Use
Kernel ADKernel	computeQpResidual	Use when the term in the PDE is multiplied by both the test function and the gradient of the test function (`_test` and `_grad_test` must be applied)
KernelValue ADKernelValue	precomputeQpResidual	Use when the term computed in the PDE is only multiplied by the test function (do not use `_test` in the override, it is applied automatically)
KernelGrad ADKernelGrad	precomputeQpResidual	Use when the term computed in the PDE is only multiplied by the gradient of the test function (do not use `_grad_test` in the override, it is applied automatically)

Diffusion

Recall the steady-state diffusion equation on the 3D domain $\Omega$ :

-\nabla \cdot \nabla u = 0 \in \Omega

The weak form of this equation includes a volume integral, which in inner-product notation, is given by:

\left(\nabla \psi_i, \nabla u_h\right) = 0 \quad\forall \psi_i,

where $\psi_i$ are the test functions and $u_h$ is the finite element solution.

Problem Statement

Given a domain $\Omega$ , find $u$ such that:

-\nabla \cdot \left( k(u) \nabla u \right) + \kappa u = 0 \in \Omega,

and

\quad k(u) \nabla u \cdot \hat{n} = \sigma \in \partial \Omega,

where

k(u) \equiv \frac{1}{\sqrt{1 + |\nabla u|^2}} \\ \kappa=1 \\ \sigma=0.2

Extension	Description
.dat	Tecplot ASCII file
.e, .exd	Sandia's ExodusII format
.fro	ACDL's surface triangulation file
.gmv	LANL's GMV (General Mesh Viewer) format
.mat	Matlab triangular ASCII file (read only)
.msh	GMSH ASCII file
.n, .nem	Sandia's Nemesis format
.plt	Tecplot binary file (write only)
.node, .ele; .poly	TetGen ASCII file (read; write)
.inp	Abaqus .inp format (read only)
.ucd	AVS's ASCII UCD format
.unv	I-deas Universal format
.xda, .xdr	libMesh formats
.vtk, .pvtu	Visualization Toolkit

Distributed Mesh

Changing the type to distributed when running in parallel operates such that only the portion of the mesh owned by a processor is stored on that processor.

parallel_type = distributed

If the mesh is too large to read in on a single processor, it can be split prior to the simulation.

Copy the mesh to a large memory machine
Use the --split-mesh option to split the mesh into $n$ pieces
Run the executable with --use-split

Short-cut	Sub-block ("type=")	Description
console	Console	Writes to the screen and optionally a file
exodus	Exodus	The most common,well supported, and controllable output type
vtk	VTK	Visualization Toolkit format, requires `--enable-vtk` when building libMesh
gmv	GMV	General Mesh Viewer format
nemesis	Nemesis	Parallel ExodusII format
tecplot	Tecplot	Requires `--enable-tecplot` when building libMesh
xda	XDA	libMesh internal format (ascii)
xdr	XDR	libMesh internal format (binary)
csv	CSV	Comma separated scalar values
gnuplot	GNUPlot	Only support scalar outputs
checkpoint	Checkpoint	MOOSE internal format used for restart and recovery

Type	AuxKernel	Variable Name(s)
Real	`MaterialRealAux`	prop
RealVectorValue	`MaterialRealVectorValueAux`	prop_1, prop_2, and prop_3
RealTensorValue	`MaterialRealTensorValueAux`	prop_11, prop_12, prop_13, prop_21, etc.

Function System

A system for defining analytic expressions based on the spatial location ( $x$ , $y$ , $z$ ) and time, $t$ .

DarcyVelocity AuxKernel

The primary unknown ("nonlinear variable") is the pressure

Once the pressure is computed, the AuxiliarySystem can compute and output the velocity field using the coupled pressure variable and the permeability and viscosity properties.

Auxiliary variables come in two flavors: Nodal and Elemental.

Nodal auxiliary variables cannot couple to gradients of nonlinear variables since gradients of $C_0$ continuous variables are not well-defined at the nodes.

Elemental auxiliary variables can couple to gradients of nonlinear variables since evaluation occurs in the element interiors.

Option	Definition
`dt`	Starting time step size
`num_steps`	Number of time steps
`start_time`	The start time of the simulation
`end_time`	The end time of the simulation
`scheme`	Time integration scheme (discussed next)

Option	Definition
`steady_state_detection`	Whether to try and detect achievement of steady-state (Default = `false`)
`steady_state_tolerance`	Used for determining a steady-state; Compared against the difference in solution vectors between current and old time steps (Default = `1e-8`)

Option	Definition
`l_tol`	Linear Tolerance (default: 1e-5)
`l_max_its`	Max Linear Iterations (default: 10000)
`nl_rel_tol`	Nonlinear Relative Tolerance (default: 1e-8)
`nl_abs_tol`	Nonlinear Absolute Tolerance (default: 1e-50)
`nl_max_its`	Max Nonlinear Iterations (default: 50)

Base	Override	Use
TimeKernel ADTimeKernel	computeQpResidual	Use when the time term in the PDE is multiplied by both the test function and the gradient of the test function (`_test` and `_grad_test` must be applied)
TimeKernelValue ADTimeKernelValue	precomputeQpResidual	Use when the time term computed in the PDE is only multiplied by the test function (do not use `_test` in the override, it is applied automatically)
TimeKernelGrad ADTimeKernelGrad	precomputeQpResidual	Use when the time term computed in the PDE is only multiplied by the gradient of the test function (do not use `_grad_test` in the override, it is applied automatically)

Convergence Rates

Consider the test problem:

\begin{array}{rl} \frac{\partial u}{\partial t} - \nabla^2 u &= f \\ u(t=0)&= u_0 \\ \left. u \right|_{\partial \Omega} &= u_D \end{array}

for $t=(0,T]$ , and $\Omega=(-1,1)^2$ , $f$ is chosen so the exact solution is given by $u = t^3 (x^2 + y^2)$ and $u_0$ and $u_D$ are the initial and Dirichlet boundary conditions corresponding to this exact solution.

TimeSequenceStepper

Provide a vector of time points using parameter time_sequence, the object simply moves through these time points.

The $t_{start}$ and $t_{end}$ parameters are automatically added to the sequence.

Only time points satisfying $t_{start} < t <t_{end}$ are considered.

If a solve fails at step $n$ an additional time point $t_{new} = \frac{1}{2}(t_{n+1}+t_n)$ is inserted and the step is resolved.

ADIntegratedBC Object Members

_u, _grad_u
Value and gradient of the variable this Kernel is operating on

_test, _grad_test
Value ( $\psi$ ) and gradient ( $\nabla \psi$ ) of the test functions at the quadrature points

_phi, _grad_phi
Value ( $\phi$ ) and gradient ( $\nabla \phi$ ) of the trial functions at the quadrature points

_i, _j, _qp
Current index for test function, trial functions, and quadrature point, respectively

_normals:
Outward normal vector for boundary element

_boundary_id
The boundary ID

_current_elem, _current_side
A pointer to the element and index to the current boundary side

Integrated BCs

Integrated BCs (including Neumann BCs) are actually integrated over the external face of an element.

\left\{ \begin{array}{rl} (\nabla u, \nabla \psi_i) - (f, \psi_i) - \langle \nabla u\cdot \hat{\boldsymbol n}, \psi_i\rangle &= 0 \quad \forall i \\ \nabla u \cdot \hat{\boldsymbol n} &= g_1\quad \text{on} \quad\partial\Omega \end{array} \right.

becomes:

(\nabla u, \nabla \psi_i) - (f, \psi_i) - \langle g_1, \psi_i\rangle = 0 \quad \forall i

If $\nabla u \cdot \hat{\boldsymbol n} = 0$ , then the boundary integral is zero ("natural boundary condition").

The pressure and heat equations have been developed independently up to this point, in this step, they are coupled together.

-\nabla \cdot \frac{\mathbf{K}}{\mu} \nabla p = 0 \\ \rho c_p \frac{\partial T}{\partial t} + \underbrace{\rho c_p \epsilon \vec{u}\cdot\nabla T}_{\textrm{DarcyAdvection}} - \nabla \cdot k \nabla T = 0

Objects have been created for everything except the $\vec{u}\cdot\nabla T$ term; a Kernel, DarcyAdvection, will be developed for this term.
A more sophisticated Material object will be created that includes temperature dependence.

h-Adaptivity

$h$ -adaptivity is a method of automatically refining/coarsening the mesh in regions of high/low estimated solution error.

Concentrate degrees of freedom (DOFs) where the error is highest, while reducing DOFs where the solution is already well-captured.

Compute a measure of error using an Indicator object
Mark an element for refinement or coarsening based on the error using a Marker object

Mesh adaptivity can be employed in both Steady and Transient executioners.

Indicator Objects

Indicators report an amount of "error" for each element, built-in Indicators include:

GradientJumpIndicator
Jump in the gradient of a variable across element edges. A good "curvature" indicator that works well over a wide range of problems.

FluxJumpIndicator
Similar to GradientJump, except that a scalar coefficient (e.g. thermal conductivity) can be provided to produce a physical "flux" quantity.

LaplacianJumpIndicator
Jump in the second derivative of a variable. Only useful for $C^1$ shape functions.

AnalyticIndicator
Computes the difference between the finite element solution and a user-supplied Function representing the analytic solution to the problem.

Step 8: Input File

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 30
    ny = 3
    xmax = 0.304 # Length of test chamber
    ymax = 0.0257 # Test chamber radius
  []
  coord_type = RZ
  rz_coord_axis = X
  uniform_refine = 2
[]

[Variables]
  [pressure]
  []
  [temperature]
    initial_condition = 300 # Start at room temperature
  []
[]

[AuxVariables]
  [velocity]
    order = CONSTANT
    family = MONOMIAL_VEC
  []
[]

[Kernels]
  [darcy_pressure]
    type = DarcyPressure
    variable = pressure
  []
  [heat_conduction]
    type = ADHeatConduction
    variable = temperature
  []
  [heat_conduction_time_derivative]
    type = ADHeatConductionTimeDerivative
    variable = temperature
  []
  [heat_convection]
    type = DarcyAdvection
    variable = temperature
    pressure = pressure
  []
[]

[AuxKernels]
  [velocity]
    type = DarcyVelocity
    variable = velocity
    execute_on = timestep_end
    pressure = pressure
  []
[]

[BCs]
  [inlet]
    type = DirichletBC
    variable = pressure
    boundary = left
    value = 4000 # (Pa) From Figure 2 from paper.  First data point for 1mm spheres.
  []
  [outlet]
    type = DirichletBC
    variable = pressure
    boundary = right
    value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
  []
  [inlet_temperature]
    type = FunctionDirichletBC
    variable = temperature
    boundary = left
    function = 'if(t<0,350+50*t,350)'
  []
  [outlet_temperature]
    type = HeatConductionOutflow
    variable = temperature
    boundary = right
  []
[]

[Materials]
  [column]
    type = PackedColumn
    radius = 1
    temperature = temperature
    porosity = '0.25952 + 0.7*y/0.0257'
  []
[]

[Postprocessors]
  [average_temperature]
    type = ElementAverageValue
    variable = temperature
  []
  [outlet_heat_flux]
    type = ADSideDiffusiveFluxIntegral
    variable = temperature
    boundary = right
    diffusivity = thermal_conductivity
  []
[]

[VectorPostprocessors]
  [temperature_sample]
    type = LineValueSampler
    num_points = 500
    start_point = '0.1 0      0'
    end_point =   '0.1 0.0257 0'
    variable = temperature
    sort_by = y
  []
[]

[Problem]
  type = FEProblem
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  automatic_scaling = true

  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'

  end_time = 100
  dt = 0.25
  start_time = -1

  steady_state_tolerance = 1e-5
  steady_state_detection = true

  [TimeStepper]
    type = FunctionDT
    function = 'if(t<0,0.1,0.25)'
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

(tutorials/darcy_thermo_mech/step08_postprocessors/problems/step8.i)

Mechanics

Compute the elastic and thermal strain if the tube is only allowed to expand along the axial (y) direction.

\begin{aligned} \nabla \cdot (\boldsymbol{\sigma} + \boldsymbol{\sigma}_0) + \boldsymbol{b} =& \boldsymbol{0} \;\mathrm{in}\;\Omega \\ \boldsymbol{u} =& \boldsymbol{g}\;\mathrm{in}\;\Gamma_{ \boldsymbol{g}} \\ \boldsymbol{\sigma} \cdot \boldsymbol{n}=&\boldsymbol{t}\;\mathrm{in}\;\Gamma_{ \boldsymbol{t}} \end{aligned}

where $\boldsymbol{\sigma}$ is the Cauchy stress tensor, $\boldsymbol{\sigma}_0$ is an additional source of stress (such as pore pressure), $\boldsymbol{u}$ is the displacement vector, $\boldsymbol{b}$ is the body force, $\boldsymbol{n}$ is the unit normal to the boundary, $\boldsymbol{g}$ is the prescribed displacement on the boundary and $\boldsymbol{t}$ is the prescribed traction on the boundary.

Optimization

The MOOSE optimization module provides functionality for solving inverse optimization problems in MOOSE. It is based on PDE constrained optimization using the PETSc TAO optimization solver.

Figure 1: Optimization cycle example for parameterizing an internal heat source distribution $q_v$ to match the simulated and experimental temperature field, $T$ and $\widetilde{T}$ , respectively.

PackedColumn.h

#pragma once

#include "Material.h"

// A helper class from MOOSE that linear interpolates x,y data
#include "LinearInterpolation.h"

/**
 * Material-derived objects override the computeQpProperties()
 * function.  They must declare and compute material properties for
 * use by other objects in the calculation such as Kernels and
 * BoundaryConditions.
 */
class PackedColumn : public Material
{
public:
  static InputParameters validParams();

  PackedColumn(const InputParameters & parameters);

protected:
  /**
   * Necessary override.  This is where the values of the properties
   * are computed.
   */
  virtual void computeQpProperties() override;

  /**
   * Helper function for reading CSV data for use in an interpolator object.
   */
  bool initInputData(const std::string & param_name, ADLinearInterpolation & interp);

  /// The radius of the spheres in the column
  const Function & _input_radius;

  /// The input porosity
  const Function & _input_porosity;

  /// Temperature
  const ADVariableValue & _temperature;

  /// Compute permeability based on the radius (mm)
  LinearInterpolation _permeability_interpolation;

  /// Fluid viscosity
  bool _use_fluid_mu_interp;
  const Real & _fluid_mu;
  ADLinearInterpolation _fluid_mu_interpolation;

  /// Fluid thermal conductivity
  bool _use_fluid_k_interp = false;
  const Real & _fluid_k;
  ADLinearInterpolation _fluid_k_interpolation;

  /// Fluid density
  bool _use_fluid_rho_interp = false;
  const Real & _fluid_rho;
  ADLinearInterpolation _fluid_rho_interpolation;

  /// Fluid specific heat
  bool _use_fluid_cp_interp;
  const Real & _fluid_cp;
  ADLinearInterpolation _fluid_cp_interpolation;

  /// Fluid thermal expansion coefficient
  bool _use_fluid_cte_interp;
  const Real & _fluid_cte;
  ADLinearInterpolation _fluid_cte_interpolation;

  /// Solid thermal conductivity
  bool _use_solid_k_interp = false;
  const Real & _solid_k;
  ADLinearInterpolation _solid_k_interpolation;

  /// Solid density
  bool _use_solid_rho_interp = false;
  const Real & _solid_rho;
  ADLinearInterpolation _solid_rho_interpolation;

  /// Solid specific heat
  bool _use_solid_cp_interp;
  const Real & _solid_cp;
  ADLinearInterpolation _solid_cp_interpolation;

  /// Solid thermal expansion coefficient
  bool _use_solid_cte_interp;
  const Real & _solid_cte;
  ADLinearInterpolation _solid_cte_interpolation;

  /// The permeability (K)
  ADMaterialProperty<Real> & _permeability;

  /// The porosity (eps)
  ADMaterialProperty<Real> & _porosity;

  /// The viscosity of the fluid (mu)
  ADMaterialProperty<Real> & _viscosity;

  /// The bulk thermal conductivity
  ADMaterialProperty<Real> & _thermal_conductivity;

  /// The bulk heat capacity
  ADMaterialProperty<Real> & _specific_heat;

  /// The bulk density
  ADMaterialProperty<Real> & _density;

  /// The bulk thermal expansion coefficient
  ADMaterialProperty<Real> & _thermal_expansion;
};

(tutorials/darcy_thermo_mech/step09_mechanics/include/materials/PackedColumn.h)

PackedColumn.C

#include "PackedColumn.h"
#include "Function.h"
#include "DelimitedFileReader.h"

registerMooseObject("DarcyThermoMechApp", PackedColumn);

InputParameters
PackedColumn::validParams()
{
  InputParameters params = Material::validParams();
  params.addRequiredCoupledVar("temperature", "The temperature (C) of the fluid.");

  // Add a parameter to get the radius of the spheres in the column
  // (used later to interpolate permeability).
  params.addParam<FunctionName>("radius",
                                "1.0",
                                "The radius of the steel spheres (mm) that are packed in the "
                                "column for computing permeability.");

  // http://en.wikipedia.org/wiki/Close-packing_of_equal_spheres
  params.addParam<FunctionName>(
      "porosity", 0.25952, "Porosity of porous media, default is for closed packed spheres.");

  // Fluid properties
  params.addParam<Real>(
      "fluid_viscosity", 1.002e-3, "Fluid viscosity (Pa s); default is for water at 20C).");
  params.addParam<FileName>(
      "fluid_viscosity_file",
      "The name of a file containing the fluid viscosity (Pa-s) as a function of temperature "
      "(C); if provided the constant value is ignored.");

  params.addParam<Real>("fluid_thermal_conductivity",
                        0.59803,
                        "Fluid thermal conductivity (W/(mK); default is for water at 20C).");
  params.addParam<FileName>(
      "fluid_thermal_conductivity_file",
      "The name of a file containing fluid thermal conductivity (W/(mK)) as a function of "
      "temperature (C); if provided the constant value is ignored.");

  params.addParam<Real>(
      "fluid_density", 998.21, "Fluid density (kg/m^3); default is for water at 20C).");
  params.addParam<FileName>("fluid_density_file",
                            "The name of a file containing fluid density (kg/m^3) as a function "
                            "of temperature (C); if provided the constant value is ignored.");

  params.addParam<Real>(
      "fluid_specific_heat", 4157.0, "Fluid specific heat (J/(kgK); default is for water at 20C).");
  params.addParam<FileName>(
      "fluid_specific_heat_file",
      "The name of a file containing fluid specific heat (J/(kgK) as a function of temperature "
      "(C); if provided the constant value is ignored.");

  params.addParam<Real>("fluid_thermal_expansion",
                        2.07e-4,
                        "Fluid thermal expansion coefficient (1/K); default is for water at 20C).");
  params.addParam<FileName>("fluid_thermal_expansion_file",
                            "The name of a file containing fluid thermal expansion coefficient "
                            "(1/K) as a function of temperature "
                            "(C); if provided the constant value is ignored.");

  // Solid properties
  // https://en.wikipedia.org/wiki/Stainless_steel#Properties
  params.addParam<Real>("solid_thermal_conductivity",
                        15.0,
                        "Solid thermal conductivity (W/(mK); default is for AISI/ASTIM 304 "
                        "stainless steel at 20C).");
  params.addParam<FileName>(
      "solid_thermal_conductivity_file",
      "The name of a file containing solid thermal conductivity (W/(mK)) as a function of "
      "temperature (C); if provided the constant value is ignored.");

  params.addParam<Real>(
      "solid_density",
      7900,
      "Solid density (kg/m^3); default is for AISI/ASTIM 304 stainless steel at 20C).");
  params.addParam<FileName>("solid_density_file",
                            "The name of a file containing solid density (kg/m^3) as a function "
                            "of temperature (C); if provided the constant value is ignored.");

  params.addParam<Real>(
      "solid_specific_heat",
      500,
      "Solid specific heat (J/(kgK); default is for AISI/ASTIM 304 stainless steel at 20C).");
  params.addParam<FileName>(
      "solid_specific_heat_file",
      "The name of a file containing solid specific heat (J/(kgK) as a function of temperature "
      "(C); if provided the constant value is ignored.");

  params.addParam<Real>("solid_thermal_expansion",
                        17.3e-6,
                        "Solid thermal expansion coefficient (1/K); default is for water at 20C).");
  params.addParam<FileName>("solid_thermal_expansion_file",
                            "The name of a file containing solid thermal expansion coefficient "
                            "(1/K) as a function of temperature "
                            "(C); if provided the constant value is ignored.");
  return params;
}

PackedColumn::PackedColumn(const InputParameters & parameters)
  : Material(parameters),
    // Get the one parameter from the input file
    _input_radius(getFunction("radius")),
    _input_porosity(getFunction("porosity")),
    _temperature(adCoupledValue("temperature")),

    // Fluid
    _fluid_mu(getParam<Real>("fluid_viscosity")),
    _fluid_k(getParam<Real>("fluid_thermal_conductivity")),
    _fluid_rho(getParam<Real>("fluid_density")),
    _fluid_cp(getParam<Real>("fluid_specific_heat")),
    _fluid_cte(getParam<Real>("fluid_thermal_expansion")),

    // Solid
    _solid_k(getParam<Real>("solid_thermal_conductivity")),
    _solid_rho(getParam<Real>("solid_density")),
    _solid_cp(getParam<Real>("solid_specific_heat")),
    _solid_cte(getParam<Real>("solid_thermal_expansion")),

    // Material Properties being produced by this object
    _permeability(declareADProperty<Real>("permeability")),
    _porosity(declareADProperty<Real>("porosity")),
    _viscosity(declareADProperty<Real>("viscosity")),
    _thermal_conductivity(declareADProperty<Real>("thermal_conductivity")),
    _specific_heat(declareADProperty<Real>("specific_heat")),
    _density(declareADProperty<Real>("density")),
    _thermal_expansion(declareADProperty<Real>("thermal_expansion"))
{
  // Set data for permeability interpolation
  std::vector<Real> sphere_sizes = {1, 3};
  std::vector<Real> permeability = {0.8451e-9, 8.968e-9};
  _permeability_interpolation.setData(sphere_sizes, permeability);

  // Fluid viscosity, thermal conductivity, density, and specific heat
  _use_fluid_mu_interp = initInputData("fluid_viscosity_file", _fluid_mu_interpolation);
  _use_fluid_k_interp = initInputData("fluid_thermal_conductivity_file", _fluid_k_interpolation);
  _use_fluid_rho_interp = initInputData("fluid_density_file", _fluid_rho_interpolation);
  _use_fluid_cp_interp = initInputData("fluid_specific_heat_file", _fluid_cp_interpolation);
  _use_fluid_cte_interp = initInputData("fluid_thermal_expansion_file", _fluid_cte_interpolation);

  // Solid thermal conductivity, density, and specific heat
  _use_solid_k_interp = initInputData("solid_thermal_conductivity_file", _solid_k_interpolation);
  _use_solid_rho_interp = initInputData("solid_density_file", _solid_rho_interpolation);
  _use_solid_cp_interp = initInputData("solid_specific_heat_file", _solid_cp_interpolation);
  _use_solid_cte_interp = initInputData("solid_thermal_expansion_file", _solid_cte_interpolation);
}

void
PackedColumn::computeQpProperties()
{
  // Current temperature
  ADReal temp = _temperature[_qp] - 273.15;

  // Permeability
  Real radius_value = _input_radius.value(_t, _q_point[_qp]);
  mooseAssert(radius_value >= 1 && radius_value <= 3,
              "The radius range must be in the range [1, 3], but " << radius_value << " provided.");
  _permeability[_qp] = _permeability_interpolation.sample(radius_value);

  // Porosity
  Real porosity_value = _input_porosity.value(_t, _q_point[_qp]);
  mooseAssert(porosity_value > 0 && porosity_value <= 1,
              "The porosity range must be in the range (0, 1], but " << porosity_value
                                                                     << " provided.");
  _porosity[_qp] = porosity_value;

  // Fluid properties
  _viscosity[_qp] =
      _use_fluid_mu_interp ? _fluid_mu_interpolation.sample(raw_value(temp)) : _fluid_mu;
  ADReal fluid_k = _use_fluid_k_interp ? _fluid_k_interpolation.sample(raw_value(temp)) : _fluid_k;
  ADReal fluid_rho =
      _use_fluid_rho_interp ? _fluid_rho_interpolation.sample(raw_value(temp)) : _fluid_rho;
  ADReal fluid_cp =
      _use_fluid_cp_interp ? _fluid_cp_interpolation.sample(raw_value(temp)) : _fluid_cp;
  ADReal fluid_cte =
      _use_fluid_cte_interp ? _fluid_cte_interpolation.sample(raw_value(temp)) : _fluid_cte;

  // Solid properties
  ADReal solid_k = _use_solid_k_interp ? _solid_k_interpolation.sample(raw_value(temp)) : _solid_k;
  ADReal solid_rho =
      _use_solid_rho_interp ? _solid_rho_interpolation.sample(raw_value(temp)) : _solid_rho;
  ADReal solid_cp =
      _use_solid_cp_interp ? _solid_cp_interpolation.sample(raw_value(temp)) : _solid_cp;
  ADReal solid_cte =
      _use_solid_cte_interp ? _solid_cte_interpolation.sample(raw_value(temp)) : _solid_cte;

  // Compute the heat conduction material properties as a linear combination of
  // the material properties for fluid and steel.
  _thermal_conductivity[_qp] = _porosity[_qp] * fluid_k + (1.0 - _porosity[_qp]) * solid_k;
  _density[_qp] = _porosity[_qp] * fluid_rho + (1.0 - _porosity[_qp]) * solid_rho;
  _specific_heat[_qp] = _porosity[_qp] * fluid_cp + (1.0 - _porosity[_qp]) * solid_cp;
  _thermal_expansion[_qp] = _porosity[_qp] * fluid_cte + (1.0 - _porosity[_qp]) * solid_cte;
}

bool
PackedColumn::initInputData(const std::string & param_name, ADLinearInterpolation & interp)
{
  if (isParamValid(param_name))
  {
    const std::string & filename = getParam<FileName>(param_name);
    MooseUtils::DelimitedFileReader reader(filename, &_communicator);
    reader.setComment("#");
    reader.read();
    interp.setData(reader.getData(0), reader.getData(1));
    return true;
  }
  return false;
}

(tutorials/darcy_thermo_mech/step09_mechanics/src/materials/PackedColumn.C)

Step 9: Input File

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Mesh]
  [generate]
    type = GeneratedMeshGenerator
    dim = 2
    ny = 200
    nx = 10
    ymax = 0.304 # Length of test chamber
    xmax = 0.0257 # Test chamber radius
  []
  [bottom]
    type = SubdomainBoundingBoxGenerator
    input = generate
    location = inside
    bottom_left = '0 0 0'
    top_right = '0.01285 0.304 0'
    block_id = 1
  []
  coord_type = RZ
[]

[Variables]
  [pressure]
  []
  [temperature]
    initial_condition = 300 # Start at room temperature
  []
[]

[AuxVariables]
  [velocity]
    order = CONSTANT
    family = MONOMIAL_VEC
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    # This block adds all of the proper Kernels, strain calculators, and Variables
    # for SolidMechanics in the correct coordinate system (autodetected)
    add_variables = true
    strain = FINITE
    eigenstrain_names = eigenstrain
    use_automatic_differentiation = true
    generate_output = 'vonmises_stress elastic_strain_xx elastic_strain_yy strain_xx strain_yy'
  []
[]

[Kernels]
  [darcy_pressure]
    type = DarcyPressure
    variable = pressure
  []
  [heat_conduction]
    type = ADHeatConduction
    variable = temperature
  []
  [heat_conduction_time_derivative]
    type = ADHeatConductionTimeDerivative
    variable = temperature
  []
  [heat_convection]
    type = DarcyAdvection
    variable = temperature
    pressure = pressure
  []
[]

[AuxKernels]
  [velocity]
    type = DarcyVelocity
    variable = velocity
    execute_on = timestep_end
    pressure = pressure
  []
[]

[BCs]
  [inlet]
    type = DirichletBC
    variable = pressure
    boundary = bottom
    value = 4000 # (Pa) From Figure 2 from paper.  First data point for 1mm spheres.
  []
  [outlet]
    type = DirichletBC
    variable = pressure
    boundary = top
    value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
  []
  [inlet_temperature]
    type = FunctionDirichletBC
    variable = temperature
    boundary = bottom
    function = 'if(t<0,350+50*t,350)'
  []
  [outlet_temperature]
    type = HeatConductionOutflow
    variable = temperature
    boundary = top
  []
  [hold_inlet]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  []
  [hold_center]
    type = DirichletBC
    variable = disp_r
    boundary = left
    value = 0
  []
  [hold_outside]
    type = DirichletBC
    variable = disp_r
    boundary = right
    value = 0
  []
[]

[Materials]
  viscosity_file = data/water_viscosity.csv
  density_file = data/water_density.csv
  thermal_conductivity_file = data/water_thermal_conductivity.csv
  specific_heat_file = data/water_specific_heat.csv
  thermal_expansion_file = data/water_thermal_expansion.csv

  [column_top]
    type = PackedColumn
    block = 0
    temperature = temperature
    radius = 1.15
    fluid_viscosity_file = ${viscosity_file}
    fluid_density_file = ${density_file}
    fluid_thermal_conductivity_file = ${thermal_conductivity_file}
    fluid_specific_heat_file = ${specific_heat_file}
    fluid_thermal_expansion_file = ${thermal_expansion_file}
  []
  [column_bottom]
    type = PackedColumn
    block = 1
    temperature = temperature
    radius = 1
    fluid_viscosity_file = ${viscosity_file}
    fluid_density_file = ${density_file}
    fluid_thermal_conductivity_file = ${thermal_conductivity_file}
    fluid_specific_heat_file = ${specific_heat_file}
    fluid_thermal_expansion_file = ${thermal_expansion_file}
  []

  [elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 200e9 # (Pa) from wikipedia
    poissons_ratio = .3 # from wikipedia

  []
  [elastic_stress]
    type = ADComputeFiniteStrainElasticStress
  []
  [thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    stress_free_temperature = 300
    eigenstrain_name = eigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1e-5 # TM modules doesn't support material property, but it will
  []
[]

[Postprocessors]
  [average_temperature]
    type = ElementAverageValue
    variable = temperature
  []
[]

[Problem]
  type = FEProblem
[]

[Executioner]
  type = Transient
  start_time = -1
  end_time = 200
  steady_state_tolerance = 1e-7
  steady_state_detection = true
  dt = 0.25
  solve_type = PJFNK
  automatic_scaling = true
  compute_scaling_once = false
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  #petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  #petsc_options_value = 'hypre boomeramg 500'
  line_search = none
  [TimeStepper]
    type = FunctionDT
    function = 'if(t<0,0.1,0.25)'
  []
[]

[Outputs]
  [out]
    type = Exodus
    elemental_as_nodal = true
  []
[]

(tutorials/darcy_thermo_mech/step09_mechanics/problems/step9.i)

Coupling terminology

Loosely-Coupled
- Each physics solved with a separate linear/nonlinear solve.
- Data exchange once per timestep (typically)
Tightly-Coupled / Picard
- Each physics solved with a separate linear/nonlinear solve.
- Data is exchanged and physics re-solved until “convergence”
Fully-Coupled
- All physics solved for in one linear/nonlinear solve

Example scheme (implicit-explicit)

\text{solve }M(u_n, v_n) u_{n+1/2} = 0\\ \text{then }N(u_{n+1/2}, v_n) v_{n+1} = 0\\

\text{solve }M(u_{n,i}, v_{n,i}) u_{n,i+1} = 0\\ \text{then }N(u_{n,i+1}, v_{n,i}) v_{n,i+1} = 0\\ \text{then }M(u_{n,i+1}, v_{n,i+1}) u_{n,i+2} = 0\\ \text{etc }

\text{solve }\begin{bmatrix}M(u_n, v_n) \\ N(u_n, v_n)\end{bmatrix} \begin{bmatrix}u_n v_n\end{bmatrix} = \begin{bmatrix}0 0\end{bmatrix} \\

SetupDarcySimulation.C

#include "SetupDarcySimulation.h"

// MOOSE includes
#include "FEProblem.h"
#include "AuxiliarySystem.h"

// libMesh includes
#include "libmesh/fe.h"

registerMooseAction("DarcyThermoMechApp", SetupDarcySimulation, "add_aux_variable");
registerMooseAction("DarcyThermoMechApp", SetupDarcySimulation, "add_aux_kernel");
registerMooseAction("DarcyThermoMechApp", SetupDarcySimulation, "add_kernel");

InputParameters
SetupDarcySimulation::validParams()
{
  InputParameters params = Action::validParams();
  params.addParam<VariableName>("pressure", "pressure", "The pressure variable.");
  params.addParam<VariableName>("temperature", "temperature", "The temperature variable.");
  params.addParam<bool>(
      "compute_velocity", true, "Enable the auxiliary calculation of velocity from pressure.");
  params.addParam<bool>("compute_pressure", true, "Enable the computation of pressure.");
  params.addParam<bool>("compute_temperature", true, "Enable the computation of temperature.");
  return params;
}

SetupDarcySimulation::SetupDarcySimulation(const InputParameters & parameters)
  : Action(parameters),
    _compute_velocity(getParam<bool>("compute_velocity")),
    _compute_pressure(getParam<bool>("compute_pressure")),
    _compute_temperature(getParam<bool>("compute_temperature"))
{
}

void
SetupDarcySimulation::act()
{
  // Actual names of input variables
  const std::string pressure = getParam<VariableName>("pressure");
  const std::string temperature = getParam<VariableName>("temperature");

  // Velocity AuxVariables
  if (_compute_velocity && _current_task == "add_aux_variable")
  {
    InputParameters var_params = _factory.getValidParams("VectorMooseVariable");
    var_params.set<MooseEnum>("family") = "MONOMIAL_VEC";
    var_params.set<MooseEnum>("order") = "CONSTANT";

    _problem->addAuxVariable("VectorMooseVariable", "velocity", var_params);
  }

  // Velocity AuxKernels
  else if (_compute_velocity && _current_task == "add_aux_kernel")
  {
    InputParameters params = _factory.getValidParams("DarcyVelocity");
    params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_END;
    params.set<std::vector<VariableName>>("pressure") = {pressure};

    params.set<AuxVariableName>("variable") = "velocity";
    _problem->addAuxKernel("DarcyVelocity", "velocity", params);
  }

  // Kernels
  else if (_current_task == "add_kernel")
  {
    // Flags for aux variables
    const bool is_pressure_aux = _problem->getAuxiliarySystem().hasVariable(pressure);
    const bool is_temperature_aux = _problem->getAuxiliarySystem().hasVariable(temperature);

    // Pressure
    if (_compute_pressure && !is_pressure_aux)
    {
      InputParameters params = _factory.getValidParams("DarcyPressure");
      params.set<NonlinearVariableName>("variable") = pressure;
      _problem->addKernel("DarcyPressure", "darcy_pressure", params);
    }

    // Temperature
    if (_compute_temperature && !is_temperature_aux)
    {
      {
        InputParameters params = _factory.getValidParams("ADHeatConduction");
        params.set<NonlinearVariableName>("variable") = temperature;
        _problem->addKernel("ADHeatConduction", "heat_conduction", params);
      }

      {
        InputParameters params = _factory.getValidParams("ADHeatConductionTimeDerivative");
        params.set<NonlinearVariableName>("variable") = temperature;
        _problem->addKernel("ADHeatConductionTimeDerivative", "heat_conduction_time", params);
      }

      {
        InputParameters params = _factory.getValidParams("DarcyAdvection");
        params.set<NonlinearVariableName>("variable") = temperature;
        params.set<std::vector<VariableName>>("pressure") = {pressure};
        _problem->addKernel("DarcyAdvection", "heat_advection", params);
      }
    }
  }
}

(tutorials/darcy_thermo_mech/step11_action/src/actions/SetupDarcySimulation.C)

Step 11: Input File

[GlobalParams]
  displacements = 'disp_r disp_z'
[]

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    ny = 200
    nx = 10
    ymax = 0.304 # Length of test chamber
    xmax = 0.0257 # Test chamber radius
  []
  coord_type = RZ
[]

[Variables]
  [pressure]
  []
  [temperature]
    initial_condition = 300 # Start at room temperature
  []
[]

[DarcyThermoMech]
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    # This block adds all of the proper Kernels, strain calculators, and Variables
    # for SolidMechanics in the correct coordinate system (autodetected)
    add_variables = true
    strain = FINITE
    eigenstrain_names = eigenstrain
    use_automatic_differentiation = true
    generate_output = 'vonmises_stress elastic_strain_xx elastic_strain_yy strain_xx strain_yy'
  []
[]

[BCs]
  [inlet]
    type = DirichletBC
    variable = pressure
    boundary = bottom
    value = 4000 # (Pa) From Figure 2 from paper.  First data point for 1mm spheres.
  []
  [outlet]
    type = DirichletBC
    variable = pressure
    boundary = top
    value = 0 # (Pa) Gives the correct pressure drop from Figure 2 for 1mm spheres
  []
  [inlet_temperature]
    type = FunctionDirichletBC
    variable = temperature
    boundary = bottom
    function = 'if(t<0,350+50*t,350)'
  []
  [outlet_temperature]
    type = HeatConductionOutflow
    variable = temperature
    boundary = top
  []
  [hold_inlet]
    type = DirichletBC
    variable = disp_z
    boundary = bottom
    value = 0
  []
  [hold_center]
    type = DirichletBC
    variable = disp_r
    boundary = left
    value = 0
  []
  [hold_outside]
    type = DirichletBC
    variable = disp_r
    boundary = right
    value = 0
  []
[]

[Materials]
  viscosity_file = data/water_viscosity.csv
  density_file = data/water_density.csv
  thermal_conductivity_file = data/water_thermal_conductivity.csv
  specific_heat_file = data/water_specific_heat.csv
  thermal_expansion_file = data/water_thermal_expansion.csv

  [column]
    type = PackedColumn
    block = 0
    temperature = temperature
    radius = 1.15
    fluid_viscosity_file = ${viscosity_file}
    fluid_density_file = ${density_file}
    fluid_thermal_conductivity_file = ${thermal_conductivity_file}
    fluid_specific_heat_file = ${specific_heat_file}
    fluid_thermal_expansion_file = ${thermal_expansion_file}
  []

  [elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 200e9 # (Pa) from wikipedia
    poissons_ratio = .3 # from wikipedia

  []
  [elastic_stress]
    type = ADComputeFiniteStrainElasticStress
  []
  [thermal_strain]
    type = ADComputeThermalExpansionEigenstrain
    stress_free_temperature = 300
    eigenstrain_name = eigenstrain
    temperature = temperature
    thermal_expansion_coeff = 1e-5
  []
[]

[Postprocessors]
  [average_temperature]
    type = ElementAverageValue
    variable = temperature
  []
[]

[Problem]
  type = FEProblem
[]

[Executioner]
  type = Transient
  start_time = -1
  end_time = 200
  steady_state_tolerance = 1e-7
  steady_state_detection = true
  dt = 0.25
  solve_type = PJFNK
  automatic_scaling = true
  compute_scaling_once = false
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 500'
  line_search = none
  [TimeStepper]
    type = FunctionDT
    function = 'if(t<0,0.1,0.25)'
  []
[]

[Outputs]
  [out]
    type = Exodus
    elemental_as_nodal = true
  []
[]

(tutorials/darcy_thermo_mech/step11_action/problems/step11.i)

Error Analysis

To compare two solutions (or a solution and an analytical solution) $f_1$ and $f_2$ , the following expressions are frequently used:

||f_1-f_2||^2_{L_2(\Omega)} = \int_{\Omega} (f_1-f_2)^2 \;\text{d}\Omega \\ ||f_1-f_2||^2_{H_{1,\text{semi}}(\Omega)} = \int_{\Omega} \left|\nabla \left(f_1-f_2\right)\right|^2 \;\text{d}\Omega

From finite element theory, the convergence rates are known for these quantities on successively refined grids. These two calculations are computed in MOOSE by utilizing the ElementL2Error or ElementH1SemiError postprocessor objects, respectively.

Approximating the Solution

First, domain $\Omega$ is split into $N_c$ cells ( $\Omega_{i},~i=1,...,N_c$ ).

In MOOSE, the cell-centered finite volume approximation (Moukalled et al. (2016) and Jasak (1996)) is used due to the fact that it supports both unstructured and structured meshes:

Constant over the cells (value at cell centroid): $u \approx u_C$ if $r \in \Omega_C$
Constant over the faces (value at face centroid): $u \approx u_f$ if $r \in \partial\Omega_f$

Error analysis is done using Taylor expansions, we prefer second-order discretization/interpolation schemes.

FVM using an Advection-Diffusion-Reaction Problem

(1) Write the strong form of the equation:

-\nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u + cu - f = 0

(2) Integrate the equation over the domain $\Omega$ :

\int_\Omega \left(-\nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u +cu -f \right) dV = 0

(3) Split the integral into different terms:

-\int_\Omega \nabla\cdot k\nabla u dV + \int_\Omega \vec{\beta} \cdot \nabla u dV + \int_\Omega cu dV -\int_\Omega f dV = 0

Approximating the Reaction and Source Terms

Split into cell-wise integrals:

\int_\Omega cu ~dV = \sum_i^{N_c}\int_{\Omega_i} c u ~dV, \quad \text{and} \quad \int_\Omega f ~dV = \sum_i^{N_c}\int_{\Omega_i} f ~dV

Use the finite volume approximation (on cell C):

\int_{\Omega_C} c u ~dV\approx c_Cu_CV_C, \quad \text{and} \quad \int_{\Omega_C} f ~dV \approx f_C V_C

Approximating the Diffusion Term

Othogonal scenario: ( $S_f$ and $\delta_{NP}$ are parallel)

Non-othogonal scenario:

Split into cell-wise integrals and use the divergence theorem:

\int_\Omega\nabla\cdot k\nabla u dV = \sum_i^{N_c}\int_{\Omega_i}\nabla\cdot k\nabla u dV = \sum_i^{N_c}\int_{\partial\Omega_i} k\nabla u \cdot \hat{n} dS

On an internal cell (let's say on cell C):

\int_{\partial\Omega_C} k\nabla u \cdot \hat{n} dS = \sum_f^{N_{f,C}}\int_{\partial\Omega_{C,f}}k\nabla u \cdot \hat{n} dS \approx \sum_f^{N_{f,C}}k_f\left(\nabla u\right)_f \cdot \hat{n}_f |S_f|

$\left(\nabla u\right)_f \cdot \hat{n}$ for orthogonal scenario: $\frac{u_N-u_C}{\delta_{CN}}$ (cheap and accurate)
$\left(\nabla u\right)_f \cdot \hat{n}$ for nonorthogonal scenario: $\frac{u_N-u_C}{\delta_{CN}}|\Delta_f|+\hat{k}_f\cdot \left(\nabla u\right)_f$
- The surface normal ( $\hat{n}$ ) is expanded into a $\delta_{CN}$ -parallel ( $\hat{\Delta}_f$ ) and a remaining vector ( $\hat{k}_f$ ): $\hat{n}_f=\hat{\Delta}_f + \hat{k}_f$
- $\left(\nabla u\right)_f$ on the right hand side is computed using the interpolation of cell gradients (next slide)

Interpolation and Gradient Estimation

Basic interpolation scheme:

u_f = \left(\frac{\delta_{fN}}{\delta_{CN}}\right) u_C + \left(1- \frac{\delta_{fN}}{\delta_{CN}}\right) u_N

Cell gradients can be estimated using the Green-Gauss theorem:

\int_{\Omega_i}\nabla u ~dV \approx (\nabla u)_i V_i

At the same time:

\int_{\Omega_i}\nabla u ~dV = \int_{\partial\Omega_i}u \hat{n} dS \approx \sum_f u_f \hat{n}_f |S_f|

From the two equations:

(\nabla u)_i = \frac{1}{V_i} \sum_f u_f \hat{n}_f |S_f|

(Alternative approaches: least-squares, node-based interpolation, skewness-corrected interpolation)

Approximating the Advection Term

Split to integral to cell-wise integrals and use the divergence theorem:

\int_\Omega \vec{\beta} \cdot \nabla u~ dV = \sum_i^N\int_{\Omega_i}\vec{\beta} \cdot \nabla u~ dV = \sum_i^N\int_{\partial\Omega_i} (\vec{\beta}~u) \cdot \hat{n}~ dS

On an internal cell (let's say on cell C), assuming that $\beta$ is constant:

\int_{\partial\Omega_i} (\vec{\beta}~u) \cdot \hat{n}~ dS = \sum_f^{N_{f,C}} \int_{\partial\Omega_{C,f}} (\vec{\beta}~u) \cdot \hat{n}~ dS \approx \sum_f^{N_{f,C}} (\vec{\beta}_f~ u_f) \cdot \hat{n}_f |S_f|

Common interpolation method for $u_f$ :

Upwind
(Weighted) Average
TVD limited schemes: van Leer, Minmod

MOOSE Workshop