Previousarrow_backNextarrow_forward

Step 10: Develop an AuxKernel Object

This step introduces the auxiliary system in MOOSE, which is based on the AuxKernel class and operates on AuxVariables. These two components are analogous to Kernel objects operating on Variables, e.g., a DarcyPressure object operating on a pressure variable. However, AuxVariables do not directly affect the NonlinearSystem—the various methods that contribute to solving PDEs by the Galerkin FEM.

To demonstrate the basic use of this system, a new VectorAuxKernel object that computes the volumetric flux (velocity) associated with a given pressure gradient will be developed in accordance with Darcy's Law. A simple regression test for this new class will be created. It will also be implemented in the pressure vessel model of pressure_diffusion.i.

Auxiliary Variables

In Step 5, it was mentioned that kernels compute residual contributions from the volumetric terms in a PDE. The weak form of these terms exposes the primary dependent variable of a given problem, i.e., that which is multiplied by the test function gradient , and which can only be determined by enforcing Dirichlet boundary conditions and solving the resulting system of equations. For example, the primary variable of Darcy's equation is pressure. In the case of solid mechanics (referred to as Solid Mechanics in MOOSE), it is displacement. In MOOSE lingo, these types of variables are referred to as nonlinear variables.

In contrast to nonlinear variables, auxiliary variables are never the primary dependent. These variables can always be computed explicitly based on given information during any part of a solve procedure. But, similar to nonlinear variables, they are continuous and possibly differentiable over individual finite elements and vary based on an assumed shape function. Often, they depend on nonlinear variables and/or vice versa. For example, one auxiliary variable in Darcy's equation is the volumetric flux, which depends on the gradient of a scalar pressure field. In the case of solid mechanics, stress and strain, which both depend on the gradient of a displacement field, are considered as auxiliary variables. In structural dynamics, acceleration and velocity are considered to be auxiliary, but the displacement at any given moment depends on their values. In other cases still, auxiliary variables may have no relation to the nonlinear ones.

Auxiliary variables are declared in the [AuxVariables] block and have access to the same parameters that their nonlinear counterparts do, e.g., "order" and "family". They are designated as one of the following two types depending on their input for the "family" parameter:

  • Elemental (family = MONOMIAL or family = MONOMIAL_VEC)

  • Nodal (family = LAGRANGE or family = LAGRANGE_VEC)

It is necessary to distinguish AuxVariables based on FE shape function family because their selection dictates how they are treated by AuxKernel objects, which are discussed in the next section.

commentnote:Variable coupling rules

Elemental variables can couple to any other type of variable, including nonlinear variables and other auxiliary ones. Nodal variables, however, can only couple to nonlinear variables and other nodal auxiliary ones, but not elemental ones.

Auxiliary Kernel Objects

The AuxKernels System in MOOSE is one for computing and setting the explicitly known values of auxiliary variables. The main advantage of AuxKernel (or VectorAuxKernel) objects is that the variables they operate on are treated like any other MOOSE variable and can be used to decouple systems of equations or for various postprocessing tasks. These objects have access to the following variable members:

  • _u, _grad_u
    Value and gradient of the variable being operated on.

  • _q_point
    Coordinates of the current quadrature point. This member is only valid when operating elemental variables.

  • _qp
    Current quadrature point or node index.

  • _current_elem
    Pointer to the current element being operated on (elemental only).

  • _current_node
    Pointer to the current node being operate on (nodal only).

One major difference between AuxKernel objects and Kernel ones is that the former do not compute residuals and, therefore, do not involve test functions . Also, all AuxKernel objects must override the computeValue() method as opposed to computeQpResidual() (or any of the other virtual function members discussed in the Kernel Objects section).

Every AuxKernel object will add a parameter "variable" taken from the validParams() method of the base class. This parameter must provide a string identifying the AuxVariable object to write the results to. When the variable type is elemental, the computeValue() method averages the values over the element quadrature points (QPs) weighted by their Gaussian quadratures in proportion to the total volume, area, or length of the element. When the variable type is nodal, the values are computed directly at nodes and the _qp member is actually a node index.

commentnote:Only elemental variables can couple to material properties.

AuxKernel objects operating on elemental variables can easily couple to material properties because both are defined at the Gauss QPs, which will be accessed by the computeValue() method accordingly. If the variable is nodal, however, the QPs cannot be accessed and so neither can the material properties.

Demonstration

The volumetric flux of a fluid through a porous medium is governed by Darcy's Law as it was given in Problem Statement, Eq. (3). Assuming a negligible gravitational field and a homogeneous medium, the law can be expressed as

(1)

Assuming every term in Eq. (1) is known, including the pressure variable , it can be developed as a VectorAuxKernel object. This new MooseObject shall accept an input for the name of the nonlinear variable holding the values for and retrieve the and terms from the available material properties.

In a space, such as the axisymmetric pipe model developed in pressure_diffusion.i, the pressure gradient has two components, e.g., . According to Eq. (1), the -component of the volumetric flux through a distance must be

(2)

If is constant along the -direction (i.e., if ), then the pressure throughout the plane can be described by a univariate function of ; independent of . Assuming that takes the form of a linear Lagrange polynomial, this function is given by

(3)

where and are the pressures at the domain boundaries: and , respectively. Differentiating Eq. (3) and substituting it in Eq. (2) leads to the following:

(4)

Thus, for any model that satisfies the foregoing assumptions, especially so that , is a constant monomial that depends only on , , , , and .

In the pipe model of pressure_diffusion.i, and are the only boundary conditions (BCs) applied. And since there are no pressure fluxes through any of the boundaries, the model indeed satisfies . Therefore, Eq. (4) can be expected to hold for this problem because the pressure variable uses the default parameters: order = FIRST and family = LAGRANGE. To verify that a new class designed to evaluate Eq. (1) is working properly, it shall be confirmed that the L2-norm of the result in all elements is the value obtained by substituting , , , , and . In addition, a simpler model that is also expected to satisfy Eq. (4) shall be developed as a regression test.

Source Code

To evaluate Eq. (1), a new VectorAuxKernel object can be created and it shall be called DarcyVelocity. (The word "velocity" is used since has units of distance-per-time, although it does not actually represent the flow velocity of a fluid; rather, it represents discharge-per-area). Start by making the directories to store files for objects that are part of the AuxKernels System:

cd ~/projects/babbler
mkdir include/auxkernels src/auxkernels

In include/auxkernels, create a file name DarcyVelocity.h and add the code given in Listing 1. Here, the AuxKernel.h header file was included because it provides the VectorAuxKernel base class that needs to be inherited. The validParams() and constructor methods were declared as public members. In the protected member declarations, the computeValue() method from the base class was overridden. This method returns a RealVectorValue data type—a tuple of values corresponding to the components of a vector in . A reference to a variable _grad_p was declared and will be used to couple to the pressure (set by a DarcyPressure object) and compute . The VariableGradient type will numerically differentiate the input over the field before writing its point values to the _qp indices of _grad_p. Finally, two ADMaterialProperty<Real> variables, _permeability and _viscosity, were declared and will be used to consume the values for and (set by a PackedColumn object), which is exactly what the DarcyPressure class was made to do in the previous step.

Listing 1: Header file for the DarcyVelocity class.


#pragma once

#include "AuxKernel.h"

/**
 * Auxiliary kernel responsible for computing the Darcy velocity (discharge per unit area) vector
 * given permeability, viscosity, and the pressure gradient of the medium.
 */
class DarcyVelocity : public VectorAuxKernel
{
public:
  static InputParameters validParams();

  DarcyVelocity(const InputParameters & parameters);

protected:
  /// AuxKernels MUST override computeValue(), which is called on every Gauss QP for elemental
  /// AuxVariables. For nodal AuxVariables, it is called on every node instead and the _qp index
  /// automatically refers to those nodes.
  virtual RealVectorValue computeValue() override;

  /// The gradient of a coupled variable, i.e., pressure
  const VariableGradient & _grad_p;

  /// The material properties which hold the values for K and mu
  const ADMaterialProperty<Real> & _permeability;
  const ADMaterialProperty<Real> & _viscosity;
};

In src/auxkernels, create a file name DarcyVelocity.C and add the code given in Listing 2. For the dependencies, the header file for the MetaPhysicL namespace was included in addition to DarcyVelocity.h, because it provides a simple method for accessing the raw values of AD types, which are usually a complex set of containers storing a variable and all of its derivatives. The method will be used to handle the ADMaterialProperty data.

The only input that needs to be appended to those defined by VectorAuxKernel::validParams() is the name of the nonlinear variable containing the values for . Thus, the addRequiredCoupledVar() method was used to define a parameter "pressure". This InputParameters member is like addRequiredParam(), but the only type of data it accepts is a string identifying the variable object, for which, if not found, will cause an error.

In the constructor definition, all components of the gradient of the variable provided to the "pressure" parameter were retrieved by the coupledGradient() method from the Coupleable class—a member of the Interfaces System in MOOSE. The data returned is then used to set _grad_p. Next, two material properties by the names "permeability" and "viscosity" are retrieved from the available material properties and used to set _permeability and _viscosity, respectively. Finally, the required computeValue() method is defined with the right-hand side of Eq. (1) discretized over all QPs. Here, the derivatives of the _permeability and _viscosity variables are not of concern (and are zero at all QPs anyways), but the MetaPhysicL::raw_value() function will ensure that the value is returned when the AD data is passed to it.

Listing 2: Source file for the DarcyVelocity class.


#include "DarcyVelocity.h"

#include "metaphysicl/raw_type.h"

registerMooseObject("BabblerApp", DarcyVelocity);

InputParameters
DarcyVelocity::validParams()
{
  InputParameters params = VectorAuxKernel::validParams();
  params.addClassDescription("Compute volumetric flux (which has units of velocity) for a given "
                             "pressure gradient, permeability, and viscosity using Darcy's Law: "
                             "$\\vec{u} = -\\frac{\\mathbf{K}}{\\mu} \\nabla p$");

  // Required parameter for specifying the MooseVariable to couple to, i.e., pressure
  params.addRequiredCoupledVar("pressure", "The pressure field.");

  return params;
}

DarcyVelocity::DarcyVelocity(const InputParameters & parameters)
  : VectorAuxKernel(parameters),
    _grad_p(coupledGradient("pressure")),

    // Note that only AuxKernels operating on elemental AuxVariables can consume a MaterialProperty
    // reference, since they are defined at the Gauss QPs within the element
    _permeability(getADMaterialProperty<Real>("permeability")),
    _viscosity(getADMaterialProperty<Real>("viscosity"))
{
}

RealVectorValue
DarcyVelocity::computeValue()
{
  // Access the gradient of the pressure at the QP. The MetaPhysicL::raw_value() method will return
  // the computed value from an automatically differentiable type, like ADMaterialProperty, as
  // opposed to any of its gradients WRT the spatial domain, which is what we want in this case.
  return -MetaPhysicL::raw_value(_permeability[_qp] / _viscosity[_qp]) * _grad_p[_qp];
}

Be sure to recompile the application before proceeding:

cd ~/projects/babbler
make -j4

Input File

With the new DarcyVelocity class available, it is now possible to compute the volumetric flux of the fluid in the pressure vessel model. First, recall that earlier in the Demonstration section it was noted that is given by a constant monomial expression, i.e., Eq. (4). It could also be shown that is a vector of constant monomial expressions when is approximated by a linear Lagrange polynomial, even if . It therefore makes sense to add the [AuxVariables] block to pressure_diffusion.i using the following syntax:

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [velocity]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT # Since "pressure" is approximated linearly, its gradient must be constant
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL_VEC # A monomial interpolation means this is an elemental AuxVariable
  []
[]

The [velocity] block creates an elemental variable with order = CONSTANT and family = MONOMIAL_VEC, where MONOMIAL_VEC types are just a vector of MONOMIAL ones. This is good, because that's exactly what should be and it is a requirement that variables associated with DarcyVelocity objects be RealVectorValue objects (or something similar).

Proceed with adding the following [AuxKernels] block to pressure_diffusion.i:

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [velocity]
    type = DarcyVelocity
    variable = velocity # Store volumetric flux vector in "velocity" variable from above
    pressure = pressure # Couple to the "pressure" variable from above
    execute_on = TIMESTEP_END # Perform calculation at the end of the solve step - after Kernels run
  []
[]

Here, a DarcyVelocity object was created. The associated variable that stores the results was identified along with the nonlinear variable associated with the DarcyPressure object. The "execute_on" parameter is made available through the base class and is used to control when the object runs via the ExecFlagEnum methods. For this demonstration, the velocity is calculated based on the current steady-state solution. Setting execute_on = TIMESTEP_END ensures that the DarcyVelocity object is constructed only after a DarcyPressure object has computed the current solution .

Now, execute the application:

cd ~/projects/babbler/problems
../babbler-opt -i pressure_diffusion.i

If the program ran successfully, a formal test of the DarcyVelocity class is in order. Start by creating a directory to store the test files:

cd ~/projects/babbler
mkdir test/tests/auxkernels/darcy_velocity

In this folder, create a file named darcy_velocity_test.i and add the inputs given in Listing 3. The model created in this file is the same as in darcy_pressure_test.i, except with [AuxVariables] and [AuxKernels] blocks added using the same syntax just demonstrated for pressure_diffusion.i.

Notice in Listing 3 that the nonlinear variable pressure is created with order = FIRST and family = LAGRANGE, even though they are the defaults. This is to emphasize a major assumption behind the test, as it is a requirement if the resulting velocity vector is to be composed of constant monomials. Also notice that l_tol = 1e-07 is set in the [Executioner] block. This forces the algebraic solvers to converge on an absolute residual that is at least less than the initial residual before the PJFNK type solver attempts to compute a new Jacobian or, rather, its action on the solution vector (see Jacobian Definition). Using a tighter relative linear tolerance improves the overall accuracy of the solution. For testing purposes, it needs to be verified that the new DarcyVelocity class produces the expected result . However, it was found that the default linear tolerance of yields very slight deviations in between elements.

Listing 3: Input file to test the DarcyVelocity class with an Exodiff object.

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 10
[]

[Variables<<<{"href": "../../../syntax/Variables/index.html"}>>>]
  [pressure]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = LAGRANGE
  []
[]

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [velocity]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL_VEC
  []
[]

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [diffusion]
    type = DarcyPressure
    variable = pressure
  []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [velocity]
    type = DarcyVelocity
    variable = velocity
    pressure = pressure
    execute_on = TIMESTEP_END
  []
[]

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [column]
    type = PackedColumn
  []
[]

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [left]
    type = ADDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/ADDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 0
  []
  [right]
    type = ADDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/ADDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pressure
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 1
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Steady
  solve_type = PJFNK
  l_tol = 1e-07 # tighter tolerance to acheive a numerically constant velocity field on 64-bit procs
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]

Run the input file and confirm that the linear solver always iterates until it computes an absolute residual that is at least less than the one computed on the zero-th iteration:

cd ~/projects/babbler/test/tests/auxkernels/darcy_velocity
../../../../babbler-opt -i darcy_velocity_test.i

Results

Use PEACOCK to visualize the solution of the pressure_diffusion.i model:

cd ~/projects/babbler/problems
peacock -r pressure_diffusion_out.e

Use the dropdown menu to select the "velocity_" variable and render its contours. Also check that the selected component in the adjacent dropdown is "Magnitude" and confirm that the result resembles the one shown in Figure 1.

Figure 1: Pressure vessel model results for indicating differences smaller than between any two elements.

From Figure 1, it can be concluded that the DarcyVelocity object produced the expected results of at all points in the mesh. Still, the renderer seems to be registering slightly different values between elements as indicated by the many different colors displayed. For practical purposes, these differences are negligible, but they can be neutralized by using tighter solver tolerances as discussed in the Input File section. To test this hypothesis, there's no need to modify the input file as any parameter can be overridden from the terminal, e.g., the following will rerun it with a relative linear tolerance of :

cd ~/projects/babbler/problems
../babbler-opt -i pressure_diffusion.i Executioner/l_tol=1e-16

Next, rerun peacock -r on the new results and observe the uniformity in the contours for , as illustrated in Figure 2. Note that the contour mappings were rounded off to 15 decimal digits. The defaults for the "Min" and "Max" fields in the ExodusViewer tab of PEACOCK might have to be trimmed to this many significant digits to achieve a solid contour like the one shown.

Figure 2: Pressure vessel model results for using a relative linear tolerance of . Every element shows exactly the predicted velocity magnitude of .

schooltip:Be aware of the many different types of controllable tolerances.

There is usually a variety of parameters available to Executioner objects that are used to set solver tolerances. Sometimes, even slight changes to them can have profound numerical effects, and certain types of tolerances supersede others. For example, the same desired effect of constantness (up to 16 decimal digits—the available precision for IEEE 754 decimal64 formatted floats (referred to as double precision in C++)) can be achieved by setting nl_rel_tol = 1e-12, but more nonlinear iterations are required in this case.

For more information about common solver parameters, please see the Steady or Transient page.

Now, use PEACOCK to inspect the results of the DarcyVelocity test:

cd ~/projects/babbler/test/tests/auxkernels/darcy_velocity
peacock -r darcy_velocity_test_out.e

Again, select the contours for the velocity magnitude and confirm that a solid contour is rendered up to 15 significant digits, as illustrated in Figure 3. This model was created by Listing 3 and uses , , and . And since the PackedColumn object used the defaults for the "diameter" and "viscosity" parameters, and . Substituting these values in Eq. (4) yields

Thus, the figure indicates that the darcy_velocity_test.i model must be correct since .

Figure 3: Rendering of the results for produced by the DarcyVelocity test. Every element shows exactly the predicted velocity magnitude of .

Test

Since the results of the darcy_velocity_test.i input file have been deemed good, the ExodusII output can now become a certified gold file:

cd ~/projects/babbler/test/tests/auxkernels/darcy_velocity
mkdir gold
mv darcy_velocity_test_out.e gold
createTASK:Create the specification for the DarcyVelocity test.

With the darcy_velocity_test.i file written and the darcy_velocity_test_out.e file moved to gold, a test specification file named tests can be generated in test/tests/auxkernels/darcy_velocity and is left up to the reader to determine the syntax. If necessary, please review Step 8.

Be sure to write the tests file as instructed above. Finally, run TestHarness:

cd ~/projects/babbler
./run_tests -j4

If the tests passed, the terminal output should look something like that shown below.


test:kernels/darcy_pressure.test .......................................................................... OK
test:materials/packed_column.test ......................................................................... OK
test:kernels/simple_diffusion.test ........................................................................ OK
test:auxkernels/darcy_velocity.test ....................................................................... OK
--------------------------------------------------------------------------------------------------------------
Ran 4 tests in 0.4 seconds. Average test time 0.1 seconds, maximum test time 0.1 seconds.
4 passed, 0 skipped, 0 pending, 0 failed

Commit

Add all of the changes made in this step to the git tracker:

cd ~/projects/babbler
git add *
schooltip:Keep tabs on local changes.

As changes to the repository are made, it can be helpful to frequently check the outputs from git status and/or git diff. Once a file is in good shape, it can be staged using a git add and further changes can be made before committing, but be careful when using wildcards!

Now, commit and push the changes to the remote repository:

git commit -m "developed auxkernel for computing the velocity associated with a pressure gradient in accordance with Darcy's law"
git push
Previousarrow_backNextarrow_forward