Kernels System

A "Kernel" is a piece of physics. It can represent one or more operators or terms in the weak form of a partial differential equation. With all terms on the left-hand-side, their sum is referred to as the "residual". The residual is evaluated at several integration quadrature points over the problem domain. To implement your own physics in MOOSE, you create your own kernel by subclassing the MOOSE Kernel class.

The Kernel system supports the use of automatic differentiation (AD) for residual calculations, as such there are two options for creating Kernel objects: Kernel and ADKernel. To further understand automatic differentiation, please refer to the Automatic Differentiation page for more information.

In a Kernel subclass the computeQpResidual() function must be overridden. This is where you implement your PDE weak form terms. For non-AD objects the following member functions can optionally be overridden:

  • computeQpJacobian()

  • computeQpOffDiagJacobian()

These two functions provide extra information that can help the numerical solver(s) converge faster and better.

Inside your Kernel class, you have access to several member variables for computing the residual and Jacobian values in the above mentioned functions:

  • _i, _j: indices for the current test and trial shape functions respectively.

  • _qp: current quadrature point index.

  • _u, _grad_u: value and gradient of the variable this Kernel operates on; indexed by _qp (i.e. _u[_qp]).

  • _test, _grad_test: value () and gradient () of the test functions at the q-points; indexed by _i and then _qp (i.e., _test[_i][_qp]).

  • _phi, _grad_phi: value () and gradient () of the trial functions at the q-points; indexed by _j and then _qp (i.e., _phi[_j][_qp]).

  • _q_point: XYZ coordinates of the current quadrature point.

  • _current_elem: pointer to the current element being operated on.

Optimized Kernel Objects

Depending on the residual calculation being performed it is sometimes possible to optimize the calculation of the residual by precomputing values during the finite element assembly of the residual vector. The following table details the various Kernel base classes that can be used for as base classes to improve performance.

BaseOverrideUse
Kernel
ADKernel
computeQpResidualUse when the term in the partial differential equation (PDE) is multiplied by both the test function and the gradient of the test function (_test and _grad_test must be applied)
KernelValue
ADKernelValue
precomputeQpResidualUse 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
precomputeQpResidualUse 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)

Custom Kernel Creation

To create a custom kernel, you can follow the pattern of the Diffusion or ADDiffusion objects implemented and included in the MOOSE framework. Additionally, Example 2 in MOOSE provides a step-by-step overview of creating your own custom kernel. The following describes that calculation of the diffusion term of a PDE.

The strong-form of the diffusion equation is defined on a 3-D domain as: find such that

(1)

where is defined as the boundary on which the value of is fixed to a known constant , is defined as the boundary on which the flux across the boundary is fixed to a known constant , and is the boundary outward normal.

The weak form is generated by multiplying by a test function () and integrating over the domain (using inner-product notation):

and then integrating by parts which gives the weak form:

(2)

where is known as the trial function that defines the finite element discretization, , with being the basis functions.

The Jacobian, which is the derivative of Eq. (2) with respect to , is defined as:

(3)

As mentioned, the computeQpResidual method must be overridden for both flavors of kernels non-AD and AD. The computeQpResidual method for the non-AD version, Diffusion, is provided in Listing 1.

Listing 1: The C++ weak-form residual statement of Eq. (2) as implemented in the Diffusion kernel.

Real
Diffusion::computeQpResidual()
{
  return _grad_u[_qp] * _grad_test[_i][_qp];
}
(moose/framework/src/kernels/Diffusion.C)

This object also overrides the computeQpJacobian method to define Jacobian term of (moose/test/tests/test_harness/duplicate_outputs_analyzejacobian) as shown in Listing 2.

Listing 2: The C++ weak-form Jacobian statement of (moose/test/tests/test_harness/duplicate_outputs_analyzejacobian) as implemented in the Diffusion kernel.

Real
Diffusion::computeQpJacobian()
{
  return _grad_phi[_j][_qp] * _grad_test[_i][_qp];
}
(moose/framework/src/kernels/Diffusion.C)

The AD version of this object, ADDiffusion, relies on an optimized kernel object (see Optimized Kernel Objects), as such it overrides precomputeQpResidual as follows.

Listing 3: The C++ pre-computed portions of the weak-form residual statement of Eq. (2) as implemented in the ADDiffusion kernel.

ADDiffusion::precomputeQpResidual()
{
  return _grad_u[_qp];
}
(moose/framework/src/kernels/ADDiffusion.C)

Time Derivative Kernels

You can create a time-derivative term/kernel by subclassing TimeKernel instead of Kernel. For example, the residual contribution for a time derivative term is:

where is the finite element solution, and

(4)

because you can interchange the order of differentiation and summation.

In the equation above, is the time derivative of the th finite element coefficient of . While the exact form of this derivative depends on the time stepping scheme, without much loss of generality, we can assume the following form for the time derivative:

for some constants , which depend on and the timestepping method.

The derivative of equation Eq. (4) with respect to is then:

So that the Jacobian term for equation Eq. (4) is

where is what we call du_dot_du in MOOSE.

Therefore the computeQpResidual() function for our time-derivative term kernel looks like:


return _test[_i][_qp] * _u_dot[_qp];

And the corresponding computeQpJacobian() is:


return _test[_i][_qp] * _phi[_j][_qp] * _du_dot_du[_qp];

Coupling with Scalar Variables

If the weak form has contributions from scalar variables, then this contribution can be treated similarly as coupling from other spatial variables. See the Coupleable interface for how to obtain the variable values. Residual contributions are simply added to the computeQpResidual() function. Jacobian terms from the test spatial variable and incremental scalar variable are added by overriding the function computeQpOffDiagJacobianScalar().

Contributions to the scalar variable weak equation (test scalar variable terms) are not natively treated by the Kernel class. Inclusion of these residual and Jacobian contributions are discussed within ScalarKernels and specifically KernelScalarBase.

Further Kernel Documentation

Several specialized kernel types exist in MOOSE each with useful functionality. Details for each are in the sections below.

Available Objects

  • Moose App
  • ADBodyForceDemonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form .
  • ADCoefReactionImplements the residual term (p*u, test)
  • ADConservativeAdvectionConservative form of which in its weak form is given by: .
  • ADCoupledForceImplements a source term proportional to the value of a coupled variable. Weak form: .
  • ADCoupledTimeDerivativeTime derivative Kernel that acts on a coupled variable. Weak form: .
  • ADDiffusionSame as Diffusion in terms of physics/residual, but the Jacobian is computed using forward automatic differentiation
  • ADMatBodyForceKernel that defines a body force modified by a material property
  • ADMatCoupledForceKernel representing the contribution of the PDE term , where is a material property coefficient, is a coupled scalar field variable, and Jacobian derivatives are calculated using automatic differentiation.
  • ADMatDiffusionDiffusion equation kernel that takes an isotropic diffusivity from a material property
  • ADMatReactionKernel representing the contribution of the PDE term , where is a reaction rate material property, is a scalar variable (nonlinear or coupled), and whose Jacobian contribution is calculated using automatic differentiation.
  • ADMaterialPropertyValueResidual term (u - prop) to set variable u equal to a given material property prop
  • ADReactionImplements a simple consuming reaction term with weak form .
  • ADScalarLMKernelThis class is used to enforce integral of phi = V_0 with a Lagrange multiplier approach.
  • ADTimeDerivativeThe time derivative operator with the weak form of .
  • ADVectorDiffusionThe Laplacian operator (), with the weak form of . The Jacobian is computed using automatic differentiation
  • ADVectorTimeDerivativeThe time derivative operator with the weak form of .
  • AnisotropicDiffusionAnisotropic diffusion kernel with weak form given by .
  • ArrayBodyForceApplies body forces specified with functions to an array variable.
  • ArrayCoupledTimeDerivativeTime derivative Array Kernel that acts on a coupled variable. Weak form: . The coupled variable and the variable must have the same dimensionality
  • ArrayDiffusionThe array Laplacian operator (), with the weak form of .
  • ArrayReactionThe array reaction operator with the weak form of .
  • ArrayTimeDerivativeArray time derivative operator with the weak form of .
  • BodyForceDemonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form .
  • CoefReactionImplements the residual term (p*u, test)
  • CoefTimeDerivativeThe time derivative operator with the weak form of .
  • ConservativeAdvectionConservative form of which in its weak form is given by: .
  • CoupledForceImplements a source term proportional to the value of a coupled variable. Weak form: .
  • CoupledTimeDerivativeTime derivative Kernel that acts on a coupled variable. Weak form: .
  • DiffusionThe Laplacian operator (), with the weak form of .
  • DivFieldThe divergence operator optionally scaled by a constant scalar coefficient. Weak form: .
  • FunctionDiffusionDiffusion with a function coefficient.
  • GradFieldThe gradient operator optionally scaled by a constant scalar coefficient. Weak form: .
  • MassEigenKernelAn eigenkernel with weak form where is the eigenvalue.
  • MassLumpedTimeDerivativeLumped formulation of the time derivative . Its corresponding weak form is where denotes the time derivative of the solution coefficient associated with node .
  • MassMatrixComputes a finite element mass matrix
  • MatBodyForceKernel that defines a body force modified by a material property
  • MatCoupledForceImplements a forcing term RHS of the form PDE = RHS, where RHS = Sum_j c_j * m_j * v_j. c_j, m_j, and v_j are provided as real coefficients, material properties, and coupled variables, respectively.
  • MatDiffusionDiffusion equation Kernel that takes an isotropic Diffusivity from a material property
  • MatReactionKernel to add -L*v, where L=reaction rate, v=variable
  • MaterialDerivativeRankFourTestKernelClass used for testing derivatives of a rank four tensor material property.
  • MaterialDerivativeRankTwoTestKernelClass used for testing derivatives of a rank two tensor material property.
  • MaterialDerivativeTestKernelClass used for testing derivatives of a scalar material property.
  • MaterialPropertyValueResidual term (u - prop) to set variable u equal to a given material property prop
  • NullKernelKernel that sets a zero residual.
  • ReactionImplements a simple consuming reaction term with weak form .
  • ScalarLMKernelThis class is used to enforce integral of phi = V_0 with a Lagrange multiplier approach.
  • ScalarLagrangeMultiplierThis class is used to enforce integral of phi = V_0 with a Lagrange multiplier approach.
  • TimeDerivativeThe time derivative operator with the weak form of .
  • UserForcingFunctionDemonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form .
  • VectorBodyForceDemonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form .
  • VectorCoupledTimeDerivativeTime derivative Kernel that acts on a coupled vector variable. Weak form: .
  • VectorDiffusionThe Laplacian operator (), with the weak form of .
  • VectorFunctionReactionKernel representing the contribution of the PDE term , where is a function coefficient and is a vector variable.
  • VectorTimeDerivativeThe time derivative operator with the weak form of .
  • Zapdos App
  • AccelerationByAveragingAn acceleration scheme based on averaging a density over a periodic cycle
  • AxisymmetricCurlZThe Z-component of an axisymmetric curl.
  • ChargeSourceMoles_KVUsed for adding charged sources to Poisson’s equation. This kernel assumes that densities are measured in units of mol/m as opposed to #/m
  • CoeffDiffusionGeneric diffusion term (densities must be in logarithmic form), where the Jacobian is computed using forward automatic differentiation.
  • CoeffDiffusionForShootMethodThe derivative of the generic diffusion term used to calculate the sensitivity value for the shoothing method.(Densities must be in logarithmic form)
  • CoeffDiffusionLinGeneric linear diffusion term (Values are NOT in logarithmic form), where the Jacobian is computed using forward automatic differentiation.
  • DriftDiffusionGeneric drift-diffusion equation that contains both an electric field driven advection term and a diffusion term (Densities must be in logarithmic form)
  • EEDFReactionLogForShootMethodThe derivative of an EEDF reaction term used to calculate the sensitivity variable for the shoothing method.(Densities must be in logarithmic form)
  • EFieldAdvectionGeneric electric field driven advection term. (Densities must be in logarithmic form.)
  • EFieldArtDiffGeneric artificial electric field driven advection term (Densities must be in logarithmic form)
  • EFieldMagnitudeSourceElectric field magnitude term based on the electrostatic approximation
  • ElectronEnergyLossFromElasticElectron energy loss term for elastic collisions using Townsend coefficient (Densities must be in logarithmic form)
  • ElectronEnergyLossFromExcitationElectron energy loss term for inelastic excitation collisions using Townsend coefficient, the energy lost in Volts in a single excitation collision (Densities must be in logarithmic form)
  • ElectronEnergyLossFromIonizationElectron energy loss term for inelastic ionization collisions using Townsend coefficients, the energy lost in Volts in a single ionization collision (Densities must be in logarithmic form)
  • ElectronEnergyTermElasticRateElectron energy loss term for elastic collisions using reaction rate coefficients (Densities must be in logarithmic form)
  • ElectronEnergyTermRateElectron energy loss term for inelastic collisions using reaction rate coefficients. Threshold energy is the energy lost in Volts in a single collision (Densities must be in logarithmic form)
  • ElectronTimeDerivativeGeneric accumulation term for variables in logarithmic form.
  • ElectronsFromIonizationRate of production of electrons from ionization using Townsend coefficients (Electron density must be in logarithmic form)
  • ExcitationReactionRate of production of metastables from excitation using Townsend coefficients (Densities must be in logarithmic form)
  • IonsFromIonizationRate of production of ions from ionization using Townsend coefficients (Ion density must be in logarithmic form)
  • JouleHeatingJoule heating term for electrons (densities must be in logarithmic form), where the Jacobian is computed using forward automatic differentiation.
  • LogStabilizationMolesKernel stabilizes solution variable u in places where u → 0; b is the offset valuespecified by the user. A typical value for b is 20.
  • ProductAABBRxnGeneric second order reaction source term in which two molecules of v are produced from two molecules of u (Densities must be in logarithmic form)
  • ProductFirstOrderRxnGeneric first order reaction source term for u (v is the reactant and densities must be in logarithmic form)
  • ReactantAARxnGeneric second order reaction sink term for u in which twomolecules of u are consumed(Densities must be in logarithmic form)
  • ReactantFirstOrderRxnGeneric first order reaction sink term for u (u is the reactant)(Densities must be in logarithmic form)
  • ReactionSecondOrderLogForShootMethodThe derivative of a second order reaction term used to calculate the sensitivity variable for the shoothing method. (Densities must be in logarithmic form)
  • ReactionThirdOrderLogForShootMethodThe derivative of a third order reaction term used to calculate the sensitivity variable for the shoothing method. (Densities must be in logarithmic form)
  • ScaledReactionThe multiple of a given variable (Used for calculating the effective ion potential for a given collision frequency)
  • ShootMethodLogAn acceleration scheme based on the shooting method
  • TM0CylindricalThe axisymmetric wave equation for the azimuthal component of the magnetizing field.
  • TM0CylindricalErThe axisymmetric wave equation for the radial component of the electric field.
  • TM0CylindricalEzThe axisymmetric wave equation for the axial component of the electric field.
  • UserSourceUser defined source term
  • Crane App
  • ADEEDFElasticLog
  • ADEEDFElasticTownsendLog
  • ADEEDFEnergyLog
  • ADEEDFEnergyTownsendLogAdds the change in enthalpy from a chemical reaction to the electron and/or gas temperature equations.
  • ADEEDFReactionLog
  • ADEEDFReactionTownsendLog
  • ADReactionSecondOrderLog
  • ADReactionThirdOrderLog
  • EEDFElasticLog
  • EEDFElasticTownsendLog
  • EEDFEnergyLog
  • EEDFEnergyTownsendLog
  • EEDFReactionLog
  • EEDFReactionTownsendLog
  • LogStabilization
  • ReactionFirstOrder
  • ReactionFirstOrderEnergy
  • ReactionFirstOrderEnergyLog
  • ReactionFirstOrderLog
  • ReactionSecondOrder
  • ReactionSecondOrderEnergy
  • ReactionSecondOrderEnergyLog
  • ReactionSecondOrderLog
  • ReactionThirdOrder
  • ReactionThirdOrderEnergy
  • ReactionThirdOrderEnergyLog
  • ReactionThirdOrderLog
  • TimeDerivativeLogThe time derivative operator with the weak form of .
  • Squirrel App
  • ConservativeTemperatureAdvectionConservative form of which in its weak form is given by: .
  • CtrlConservativeAdvection
  • CtrlConservativeTemperatureAdvection
  • MatINSTemperatureTimeDerivativeThe time derivative operator with the weak form of .
  • NonConservativeAdvection
  • PotentialAdvection
  • VelocityFunctionConservativeAdvection
  • VelocityFunctionTemperatureAdvection