KernelScalarBase

The KernelScalarBase is part of the scalar augmentation system to complement the ScalarKernel class. Its principal purpose is to add standard quadrature loops and assembly routines to handle the contributions from a single added scalar variable to a Kernel class, including the entire row of the Jacobian. These routines will assist with representing weak form terms in a partial differential equation involving scalar variables integrated over the interior of a domain. As usual, this piece of physics is referred to as the "residual" and is evaluated at integration quadrature points within that domain. To implement your own physics in MOOSE, you create your own kernel by subclassing the MOOSE KernelScalarBase class.

The Kernel scalar augmentation system supports the use of AD for residual calculations, as such there are two options for creating field-scalar coupling objects: KernelScalarBase and ADKernelScalarBase. To further understand automatic differentiation, please refer to the Automatic Differentiation page for more information.

Developers should read all sections; users can find Parameters described at the bottom.

Creation of Kernel Scalar Coupling Classes

Each Kernel object has a focus field variable or spatial variable; its job is to contribute to the residual as well as the row of the Jacobian matrix. Herein, as in the source code of Kernels, this spatial variable will be called _var. In a coupled (multi-physics) weak form, all domain integral terms containing the test function of _var are potential candidates for Kernel contributions.

The philosophy of the scalar augmentation class KernelScalarBase is to add a focus scalar variable referred to as _kappa to the Kernel object so that all terms in the coupled weak form that involve _var, _kappa, and/or their test functions can be assembled in one or multiple class instances. This philosophy is similar to how the lower dimensional variable _lambda is added to the element faces of DGKernel and IntegratedBC objects associated with the hybrid finite element method (HFEM). Documentation for that approach can be found HFEMDiffusion and HFEMDirichletBC along with the base classes DGLowerDKernel and LowerDIntegratedBC.

In a KernelScalarBase subclass, a naming scheme is established for the quadrature point methods of the two variable types: methods contributing to the test function of _kappa have "Scalar" near the front and methods contributing to the trial function of scalar variables in the Jacobian have "Scalar" at the end. The computeScalarQpResidual() function should be overridden (see Parameters for cases where the scalar should be suppressed). The computeQpResidual() function must be overridden as usual for Kernels, although it may return zero.

For non-AD objects, several contributions to the Jacobian matrix can be optionally overridden for use in Newton-type nonlinear solvers. As mentioned later, the developer should choose and document which terms (rows) of the residual and terms (rows and columns) of the Jacobian will be attributed to an instance of the developed class. These choices can be motivated by whether some terms in the weak form can be or have already been implemented within other MOOSE classes.

computeQpJacobian(): Jacobian component d-_var-residual / d-_var
computeQpOffDiagJacobian(jvar_num): off-diagonal Jacobian component d-_var-residual / d-jvar
computeQpOffDiagJacobianScalar(svar_num): off-diagonal Jacobian component d-_var-residual / d-svar
computeScalarQpJacobian(): Jacobian component d-_kappa-residual / d-_kappa
computeScalarQpOffDiagJacobian(jvar_num): off-diagonal Jacobian component d-_kappa-residual / d-jvar
computeScalarQpOffDiagJacobianScalar(svar_num): off-diagonal Jacobian component d-_kappa-residual / d-svar

Examples of some of these methods are shown below in Examples from Source Code. Loops over the coupled variables wrap around these quadrature loops. The integer for the spatial variable is jvar_num and the integer for the scalar variable is svar_num.

Also, there are some pre-calculation routines that are called within the quadrature loop once before the loop over spatial variable test and shape functions as well as before the loop over scalar components. These methods are useful for material or stabilization calculations.

initScalarQpResidual(): evaluations depending on qp but independent of test functions
initScalarQpJacobian(jvar_num): evaluations depending on qp but independent of test and shape functions
initScalarQpOffDiagJacobian(jsvar): evaluations depending on qp but independent of test and shape functions

In addition to those mentioned in the Kernel documentation, you have access to several member variables inside your KernelScalarBase class for computing the residual and Jacobian values in the above mentioned functions:

_h, _l: indices for the current test and trial scalar component respectively.
_kappa: value of the scalar variable this Kernel operates on; indexed by _h (i.e. _kappa[_h]).
_kappa_var: ID number of this scalar variable; useful to differentiate from others.
_k_order: order (number of components) of the scalar variable.

Since the test and trial "shape" functions of a scalar are "1", variables are not required for that value. Examples of the source codes below demonstrate this fact.

warning:AD global indexing required

ADKernelScalarBase only works with MOOSE configured with global AD indexing (the default).

note:Parallelization of scalar contributions

While these quadrature loops are convenient for implementation in a single object, the speed of parallel execution may be slower due to the sequential assembly needed from each element assembling to the same scalar variable _kappa. For greater speed, the developer may instead implement the terms for computeScalarQpResidual() and computeScalarQpJacobian() through a derived class of ElementIntegralUserObject as discussed at Coupling with Spatial Variables.

Examples from Source Code

As mentioned, the computeScalarQpResidual method should be overridden for both flavors of kernels, non-AD and AD. As an example, consider the scalar residual weak form term of the ScalarLMKernel class:

$(1)var element = document.getElementById("moose-equation-405d71ca-7817-412e-be54-40e3485c0da8");katex.render(" F^{(\\lambda)} \\equiv \\int_{\\Omega} \\phi^h \\;\\text{d}\\Omega - V_0 = 0 ", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

The ScalarLMKernel class is implemented using the GenericKernelScalar template class to contain both the AD and non-AD versions within the same source files. The computeScalarQpResidual method for this class is provided in Listing 1, where _value/_pp_value is equal to $var element = document.getElementById("moose-equation-30fd985e-bc2b-439e-95c8-3b5cdb3019f1");katex.render("V_0", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

Listing 1: The C++ weak-form residual statement of Eq. (1).

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

Meanwhile, the contribution to the spatial variable residual of this object is associated with Eq. (2) and implemented in Listing 2 (note that the scalar variable _kappa is termed as $var element = document.getElementById("moose-equation-d38083cd-5910-42c8-aaa5-d21f863f5169");katex.render("\\lambda^h", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ in this weak form).

$(2)var element = document.getElementById("moose-equation-3019eae5-f9fe-40fa-af86-e57379590af9");katex.render(" F^{(\\phi)}_i \\equiv \\lambda^h \\int_{\\Omega} \\varphi_i \\;\\text{d}\\Omega ", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

Listing 2: The C++ weak-form residual statement of Eq. (2).

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

This object also overrides the computeScalarQpOffDiagJacobian method to define the Jacobian term related to Eq. (1) as shown in Listing 3.

Listing 3: The C++ weak-form Jacobian for d-_kappa-residual / d-jvar.

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

Notice that there is a conditional to confirm that the coupled jvar is the focus variable _var, otherwise it returns zero. Also, this method only returns a "Real" value since this method is only called by the non-AD version of the class during Jacobian computation; an assert is used to verify this intention.

Similarly, it also overrides the computeQpOffDiagJacobianScalar method to define the Jacobian term related to Eq. (2) as shown in Listing 4.

Listing 4: The C++ weak-form Jacobian for d-_var-residual / d-svar.

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

Also notice the conditional that confirms the coupled svar is the focus scalar _kappa, otherwise it returns zero.

Depending upon the weak form and its coupling terms between spatial and scalar variables, not all of the methods listed in Creation of Kernel Scalar Coupling Classes need to be overridden.

The AD version of this object, ADScalarLMKernel, only requires the residual implementation. A solely AD source file would only need to override computeScalarQpResidual and computeQpResidual and leave all the Jacobian methods as base definitions, which return zero. See MortarScalarBase for examples of AD-only and non-AD separate classes.

As a more complicated example of the scalar augmentation system for kernels, the SolidMechanics test app contains headers, source, and test files for an alternative implementation of the "HomogenizedTotalLagrangianStressDivergence" system from the SolidMechanics module. This Kernel is designated with the suffix "S" to distinguish from the existing objects in the module. Also, there are other intermediate classes such as "TotalLagrangianStressDivergence" that are also designated with an "S" suffix. These other classes are needed since the lower class needs to also derive from KernelScalarBase. Meanwhile, they do not need the scalar_variable parameter and function identically to their original module source object; see the Parameters section for a comment about leaving this parameter blank.

The scalar augmentation system is designed such that multiple scalar variables can be coupled to an instance of the Kernel class, each focusing on one scalar from the list. This approach is similar to how SolidMechanics module classes operator on one component variable of the displacement vector field and are coupled to the other components. The developer can decide how to organize the coupling and off-diagonal Jacobian terms in a logical way and document this for the user.

Parameters

There is one required parameters the user must supply for a kernel derived from KernelScalarBase:

scalar_variable: the focus scalar variable of the kernel, for which assembly of the residual and Jacobian contributions will occur. It must be a MooseScalarVariable. This parameter may be renamed in a derived class to be more physically meaningful.

If the scalar_variable parameter is not specified, then the derived class will behave identically to a regular Kernel, namely without any scalar functionality. This feature is useful if the scalar augmentation in inserted into a class structure with several levels and not all derived members use scalar variables.

As an example, the parameter listing is shown below for the ScalarLMKernel object with the scalar_variable parameter renamed to kappa:

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [diff]
    type = Diffusion<<<{"description": "The Laplacian operator ($-\\nabla \\cdot \\nabla u$), with the weak form of $(\\nabla \\phi_i, \\nabla u_h)$.", "href": "Diffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
  []

  [ffnk]
    type = BodyForce<<<{"description": "Demonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form $(\\psi_i, -f)$.", "href": "BodyForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    function<<<{"description": "A function that describes the body force"}>>> = ffn
  []

  [sk_lm]
    type = ScalarLMKernel<<<{"description": "This class is used to enforce integral of phi = V_0 with a Lagrange multiplier approach.", "href": "ScalarLMKernel.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    kappa<<<{"description": "Primary coupled scalar variable"}>>> = lambda
    pp_name<<<{"description": "Name of the Postprocessor containing the volume of the domain."}>>> = pp
    value<<<{"description": "Given (constant) which we want the integral of the solution variable to match."}>>> = 2.666666666666666
  []
[]

Note: to avoid an error message "Variable 'kappa' does not exist in this system", the following block should be added to the input file:

[Problem<<<{"href": "../../syntax/Problem/index.html"}>>>]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

There are also some optional parameters that can be supplied to KernelScalarBase classes. They are:

compute_scalar_residuals: Whether to compute scalar residuals. This will automatically be set to false if a scalar_variable parameter is not supplied. Other cases where the user may want to set this to false is during testing when the scalar variable is an AuxVariable and not a solution variable in the system.
compute_field_residuals: Whether to compute residuals for the primal field variable. If several KernelScalarBase objects are used in the input file to compute different rows (i.e. different variables) of the global residual, then some objects can be targeted to field variable rows and others to scalar variable rows.

(framework/src/kernels/ScalarLMKernel.C)

// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(framework/src/kernels/ScalarLMKernel.C)

// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(framework/src/kernels/ScalarLMKernel.C)

// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(framework/src/kernels/ScalarLMKernel.C)

// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#include "ScalarLMKernel.h"

registerMooseObject("MooseApp", ScalarLMKernel);
registerMooseObject("MooseApp", ADScalarLMKernel);

template <bool is_ad>
InputParameters
ScalarLMKernelTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernelScalar<is_ad>::validParams();
  params.addClassDescription("This class is used to enforce integral of phi = V_0 with a "
                             "Lagrange multiplier approach.");
  params.renameCoupledVar("scalar_variable", "kappa", "Primary coupled scalar variable");
  params.addRequiredParam<PostprocessorName>(
      "pp_name", "Name of the Postprocessor containing the volume of the domain.");
  params.addRequiredParam<Real>(
      "value", "Given (constant) which we want the integral of the solution variable to match.");

  return params;
}

template <bool is_ad>
ScalarLMKernelTempl<is_ad>::ScalarLMKernelTempl(const InputParameters & parameters)
  : GenericKernelScalar<is_ad>(parameters),
    _value(this->template getParam<Real>("value")),
    _pp_value(this->getPostprocessorValue("pp_name"))
{
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeQpResidual()
{
  return _kappa[0] * _test[_i][_qp];
}

template <bool is_ad>
GenericReal<is_ad>
ScalarLMKernelTempl<is_ad>::computeScalarQpResidual()
{
  return _u[_qp] - _value / _pp_value;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpJacobian()
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(
      !is_ad,
      "In ADScalarLMKernel, computeScalarQpJacobian should not be called. Check computeJacobian "
      "implementation.");
  return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeQpOffDiagJacobianScalar(unsigned int svar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeQpOffDiagJacobianScalar should not be called. Check "
              "computeOffDiagJacobianScalar "
              "implementation.");
  if (svar == _kappa_var)
    return _test[_i][_qp];
  else
    return 0.;
}

template <bool is_ad>
Real
ScalarLMKernelTempl<is_ad>::computeScalarQpOffDiagJacobian(unsigned int jvar)
{
  // This function will never be called for the AD version. But because C++ does
  // not support an optional function declaration based on a template parameter,
  // we must keep this template for all cases.
  mooseAssert(!is_ad,
              "In ADScalarLMKernel, computeScalarQpOffDiagJacobian should not be called. Check "
              "computeOffDiagJacobian "
              "implementation.");
  if (jvar == _var.number())
    return _phi[_j][_qp];
  else
    return 0.;
}

template class ScalarLMKernelTempl<false>;
template class ScalarLMKernelTempl<true>;

(test/tests/kernels/scalar_kernel_constraint/scalar_constraint_kernel.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  nx = 2
  ny = 2
  elem_type = QUAD9
[]

[Functions]
  [exact_fn]
    type = ParsedFunction
    value = 'x*x+y*y'
  []

  [ffn]
    type = ParsedFunction
    value = -4
  []

  [bottom_bc_fn]
    type = ParsedFunction
    value = -2*y
  []

  [right_bc_fn]
    type = ParsedFunction
    value =  2*x
  []

  [top_bc_fn]
    type = ParsedFunction
    value =  2*y
  []

  [left_bc_fn]
    type = ParsedFunction
    value = -2*x
  []
[]

[Variables]
  [u]
    family = LAGRANGE
    order = SECOND
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
[]

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

  [ffnk]
    type = BodyForce
    variable = u
    function = ffn
  []

  [sk_lm]
    type = ScalarLMKernel
    variable = u
    kappa = lambda
    pp_name = pp
    value = 2.666666666666666
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [bottom]
    type = FunctionNeumannBC
    variable = u
    boundary = 'bottom'
    function = bottom_bc_fn
  []
  [right]
    type = FunctionNeumannBC
    variable = u
    boundary = 'right'
    function = right_bc_fn
  []
  [top]
    type = FunctionNeumannBC
    variable = u
    boundary = 'top'
    function = top_bc_fn
  []
  [left]
    type = FunctionNeumannBC
    variable = u
    boundary = 'left'
    function = left_bc_fn
  []
[]

[Postprocessors]
  # integrate the volume of domain since original objects set
  # int(phi)=V0, rather than int(phi-V0)=0
  [pp]
    type = FunctionElementIntegral
    function = 1
    execute_on = initial
  []
  [l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
    execute_on = 'initial timestep_end'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
    solve_type = 'NEWTON'
  []
[]

[Executioner]
  type = Steady
  nl_rel_tol = 1e-9
  l_tol = 1.e-10
  nl_max_its = 10
  # This example builds an indefinite matrix, so "-pc_type hypre -pc_hypre_type boomeramg" cannot
  # be used reliably on this problem
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  # This is a linear problem, so we don't need to recompute the
  # Jacobian. This isn't a big deal for a Steady problems, however, as
  # there is only one solve.
  solve_type = 'LINEAR'
[]

[Outputs]
  exodus = true
  hide = lambda
[]

(test/tests/kernels/scalar_kernel_constraint/scalar_constraint_kernel.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  nx = 2
  ny = 2
  elem_type = QUAD9
[]

[Functions]
  [exact_fn]
    type = ParsedFunction
    value = 'x*x+y*y'
  []

  [ffn]
    type = ParsedFunction
    value = -4
  []

  [bottom_bc_fn]
    type = ParsedFunction
    value = -2*y
  []

  [right_bc_fn]
    type = ParsedFunction
    value =  2*x
  []

  [top_bc_fn]
    type = ParsedFunction
    value =  2*y
  []

  [left_bc_fn]
    type = ParsedFunction
    value = -2*x
  []
[]

[Variables]
  [u]
    family = LAGRANGE
    order = SECOND
  []
  [lambda]
    family = SCALAR
    order = FIRST
  []
[]

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

  [ffnk]
    type = BodyForce
    variable = u
    function = ffn
  []

  [sk_lm]
    type = ScalarLMKernel
    variable = u
    kappa = lambda
    pp_name = pp
    value = 2.666666666666666
  []
[]

[Problem]
  kernel_coverage_check = false
  error_on_jacobian_nonzero_reallocation = true
[]

[BCs]
  [bottom]
    type = FunctionNeumannBC
    variable = u
    boundary = 'bottom'
    function = bottom_bc_fn
  []
  [right]
    type = FunctionNeumannBC
    variable = u
    boundary = 'right'
    function = right_bc_fn
  []
  [top]
    type = FunctionNeumannBC
    variable = u
    boundary = 'top'
    function = top_bc_fn
  []
  [left]
    type = FunctionNeumannBC
    variable = u
    boundary = 'left'
    function = left_bc_fn
  []
[]

[Postprocessors]
  # integrate the volume of domain since original objects set
  # int(phi)=V0, rather than int(phi-V0)=0
  [pp]
    type = FunctionElementIntegral
    function = 1
    execute_on = initial
  []
  [l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
    execute_on = 'initial timestep_end'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
    solve_type = 'NEWTON'
  []
[]

[Executioner]
  type = Steady
  nl_rel_tol = 1e-9
  l_tol = 1.e-10
  nl_max_its = 10
  # This example builds an indefinite matrix, so "-pc_type hypre -pc_hypre_type boomeramg" cannot
  # be used reliably on this problem
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  # This is a linear problem, so we don't need to recompute the
  # Jacobian. This isn't a big deal for a Steady problems, however, as
  # there is only one solve.
  solve_type = 'LINEAR'
[]

[Outputs]
  exodus = true
  hide = lambda
[]

Creation of Kernel Scalar Coupling Classes
Examples from Source Code
Parameters