ArrayMooseVariable

Overview

An array variable is define as a set of standard field variables with the same finite element family and order. Each standard variable of an array variable is referred to as a component of the array variable. An array kernel is a MOOSE kernel operating on an array variable and assembles the residuals and Jacobians for all the components of the array variable. A sample array kernel can be found at ArrayDiffusion. The purpose of having array kernels is to reduce the number of kernels when the number of components of an array variable is large (potentially hundreds or thousands) and to avoid duplicated operations otherwise with lots of standard kernels. Array kernel can be useful for radiation transport where an extra independent direction variable can result into large number of standard variables. Similarly as array kernel, we have array initial conditions, array boundary conditions, array DG kernels, array interface kernels, array constraints, etc.

The design:

• to use variable groups in libMesh to group components in an array variable together.

• to use template to avoid code duplication with standard and vector variables.

• to use Eigen::Matrix to hold local dofs, solutions on quadrature points for an array variable to ease the local operations.

• to use dense matrices or vectors as standard or vector variables with proper sizes for holding the local Jacobians and residuals that are to be assembled into a global Jacobian matrix and a global residual vector.

The following map is useful for understanding the template:

The three rows correspond to standard, vector and array variables. OutputType is the data type used for templating. RealEigenVector is a typedef in (moose/framework/include/utils/MooseTypes.h) as Eigen::Matrix<Real, Eigen::Dynamic, 1>. OutputShape is for the type of shape functions and OutputData is the type of basis function expansion coefficients that are stored in the moose array variable grabbed from the solution vector.

Kernels for Array Variables

Since array variables have a different OutputData type, standard Kernels cannot be used on array variables. Kernels for array variables, or ArrayKernels, must derive from ArrayKernel which have different virtual functions for computeQpResidual, computeQpJacobian, and computeQpOffDiagJacobian. The declaration of these functions are below:

virtual void computeQpResidual(RealEigenVector & residual) = 0;
virtual RealEigenVector computeQpJacobian();
virtual RealEigenMatrix computeQpOffDiagJacobian(const MooseVariableFEBase & jvar);

When defining a computeQpResidual in a derived class, this function must define the residual in the input argument (residual). This input is already properly sized when called in ArrayKernel.C. computeQpJacobian must return a vector defining the on-diagonal terms of the Jacobian. computeQpOffDiagJacobian must return a matrix with number of rows equal to the number of components and number of columns being the number of components in jvar.

Using ArrayDiffusion as an example. The computeQpResidual function has

void
ArrayDiffusion::computeQpResidual(RealEigenVector & residual)
{
  // WARNING: the noalias() syntax is an Eigen optimization tactic, it avoids creating
  // a temporary object for the matrix multiplication on the right-hand-side. However,
  // it should be used with caution because it could cause unintended results,
  // developers should NOT use it if the vector on the left-hand-side appears on the
  // right-hand-side, for instance:
  //   vector = matrix * vector;
  // See http://eigen.tuxfamily.org/dox/group__TopicAliasing.html for more details.
  if (_d)
    residual.noalias() = (*_d)[_qp] * _grad_u[_qp] * _array_grad_test[_i][_qp];
  else if (_d_array)
    residual.noalias() = (*_d_array)[_qp].asDiagonal() * _grad_u[_qp] * _array_grad_test[_i][_qp];
  else
    residual.noalias() = (*_d_2d_array)[_qp] * _grad_u[_qp] * _array_grad_test[_i][_qp];
}

where _grad_u[_qp] is an Eigen::Matrix with number of rows being equal to the number of components of the array variable and number of columns being LIBMESH_DIM. _array_grad_test[_i][_qp] is an Eigen::Map of classic _grad_test[_i][_qp] which is in type of Gradient. Thanks to the Eigen matrix arithmetic operators, we can have a simple multiplication expression here. _d is a pointer of a material property of Real type for scalar diffusion coefficient. Here we assume the diffusion coefficient is the same for all components. _d_array is a pointer to a RealEigenVector material property, here we assume there is a diffusion coefficient for each component with no coupling. _d_2d_array is a pointer to a RealEigenMatrix material property, where the diffusion coefficient is represented as dense matrix. See ArrayDiffusion for more details.

Correspondingly the computeQpJacobian has

RealEigenVector
ArrayDiffusion::computeQpJacobian()
{
  if (_d)
    return RealEigenVector::Constant(_var.count(),
                                     _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d)[_qp]);
  else if (_d_array)
    return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_array)[_qp];
  else
    return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_2d_array)[_qp].diagonal();
}

It is noted that only the diagonal entries of the diffusion coefficients are used in the fully-coupled case because computeQpJacobian is supposed to only assemble the block-diagonal part of the Jacobian. The full local Jacobian is assembled in function computeQpOffDiagJacobian, where when the off-diagonal variable is the array variable, we have

RealEigenMatrix
ArrayDiffusion::computeQpOffDiagJacobian(const MooseVariableFEBase & jvar)
{
  if (jvar.number() == _var.number() && _d_2d_array)
    return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_2d_array)[_qp];
  else
    return ArrayKernel::computeQpOffDiagJacobian(jvar);
}

The retuned value is in type of an Eigen matrix with number of rows and columns equal to the number of components.

ArrayKernel also has virtual functions including:

virtual void initQpResidual();
virtual void initQpJacobian();
virtual void initQpOffDiagJacobian(const MooseVariableFEBase & jvar);

Which are functions that are called inside the quadrature point loop, but outside the test/shape function loop. This is useful to perform operations that depend on position (quadrature point) but do not depend on the test/shape function. For instance, ArrayDiffusion uses initQpResidual to check if the size of the vector or matrix diffusion coefficient matches the number of components in the variable:

void
ArrayDiffusion::initQpResidual()
{
  if (_d_array)
  {
    mooseAssert((*_d_array)[_qp].size() == _var.count(),
                "diffusion_coefficient size is inconsistent with the number of components of array "
                "variable");
  }
  else if (_d_2d_array)
  {
    mooseAssert((*_d_2d_array)[_qp].cols() == _var.count(),
                "diffusion_coefficient size is inconsistent with the number of components of array "
                "variable");
    mooseAssert((*_d_2d_array)[_qp].rows() == _var.count(),
                "diffusion_coefficient size is inconsistent with the number of components of array "
                "variable");
  }
}

Future work:

• To change the current dof ordering for elemental variables so that we can avoid bunch of if statements with isNodal() in MooseVariableFE.C, (refer to libMesh Issue 2114).

• To use Eigen::Map for faster solution vector access.

• To implement ArrayInterfaceKernel and ArrayConstraints.

Useful Eigen API

Linear algebra is very easy using Eigen, it has a lot of API for matrix arithmetic and manipulation that simplifies code and helps developers avoid writing their own loops. For a full list of functions, visit the Eigen matrix doxygen, which is relevant to the RealEigenVector and RealEigenMatrix types in MOOSE. Below is a list of commonly used functions for residual and jacobian evaluations. For exposition here is a definition of a vector and matrix:

• Set all elements to same value:

• Change matrix representation:

• Element-wise operations:

Eigen Tips and Tricks (Advanced)

Eigen has some unique features that, when used properly, can significantly impact performance. Here are some recommendations that can improve code performance.

• Aliasing is a technique in Eigen that constructs temporary objects when performing matrix multiplications, this is to avoid overriding data that needs to be used later in the computation. For instance vec = mat * vec will create a temporary vector for mat * vec then assign it to vec at the end. However, vec2 = mat * vec1 does not need this temporary object and assign the result to vec2 directly, this aliasing can be avoided by doing vec2.noalias(). The noalias() function should be used with extreme caution since it can cause erroneous results.

• Eigen uses what's known as expression templates, enabling operations to be known at compile time. This allows multiple operations to occur in a single element loop, providing more compiler optimization and improved cache efficiency. With this in mind, it is often better to write multiple Eigen operations in a single line or assignment. For instance, with the following syntax:

  a = 3*b + 4*c + 5*d
  

  Eigen will compile to a single for loop:

  for (unsigned int i = 0; i < b.size(); ++i)
    a[i] = 3*b[i] + 4*c[i] + 5*d[i];
  

• Eigen also has a interface for accessing raw buffers using its Map class.

For a better understanding of Eigen and using it to its full potential, it is highly recommended to go through Eigen's tutorials.

Input Parameters