Adding array kernel for computing time derivates of array variables #…

…14057

framework/doc/content/source/kernels/ArrayTimeDerivative.md

@@ -0,0 +1,38 @@
+# ArrayTimeDerivative
+
+## Description
+
+This array kernel implements the following piece of a weak form:
+\begin{equation}
+(\vec{u}^\ast, \mathbf{T} \dot{\vec{u}}),
+\end{equation}
+where $\vec{u}^\ast$ is the test function, $\dot{\vec{u}}$ is time derivative of the array of finite element solutions ($\dot{\vec{u}} = \left[\frac{\partial u_1}{\partial t},\frac{\partial u_2}{\partial t},...\right]^T$), and $\mathbf{T}$ is a matrix of the time derivative coefficients ($(\mathbf{T})_{n,m} = T_{n,m}$).
+
+Similarly as showed in [ArrayDiffusion.md], we can rearrange it into
+\begin{equation}
+(\vec{u}^\ast, \mathbf{T} \dot{\vec{u}}) = \sum_{e} \sum_{i=1}^{N_{\text{dof}}} \sum_{\text{qp}=1}^{N_{qp}} (|J|w)_{\text{qp}} \vec{w}_p\vec{u}_i^\ast \underline{\mathbf{T}_{\text{qp}} \dot{\vec{u}}_{\text{qp}} b_{i,\text{qp}}},
+\end{equation}
+where the underlined term is the vector provided by [ArrayTimeDerivative::computeQpResidual](ArrayTimeDerivative.C).
+Detailed explanations on the notations can be found in [ArrayDiffusion.md].
+
+In general, the reaction coefficient $\mathbf{T}$ is a square matrix with the size of the number of components.
+When it is a diagonal matrix, it can be represented by a vector.
+In such a case, the components are not coupled with this array time derivative kernel.
+If all elements of the time derivative coefficient vector are the same, we can use a scalar reaction coefficient.
+Thus this kernel gives users an option to set the coefficient to a scalar, vector, or matrix material property, corresponding to scalar, diagonal matrix, and full matrix, respectively.
+
+The local Jacobian can be found in the following equation:
+\begin{equation}
+J_{n,m,i,j} = \sum_{e} \sum_{i=1}^{N_{\text{dof}}} \sum_{j=1}^{N_{\text{dof}}} \sum_{\text{qp}=1}^{N_{qp}} (|J|w)_{\text{qp}} \vec{w}_p u_{n,i}^\ast \underline{T_{n,m,\text{qp}} b_{j,\text{qp}} b_{i,\text{qp}} \frac{\partial \dot{u}_{m,j}}{\partial u_{m,j}}},
+\end{equation}
+where $n$ and $m$ are the component row and column, respectively. The underlined part is the local Jacobian evaluated by [ArrayTimeDerivative::computeQpJacobian](ArrayTimeDerivative.C) and [ArrayTimeDerivative::computeQpOffDiagJacobian](ArrayTimeDerivative.C).
+
+## Example Input Syntax
+
+!listing tests/kernels/array_kernels/array_diffusion_reaction_transient.i block=Kernels
+
+!syntax parameters /Kernels/ArrayTimeDerivative
+
+!syntax inputs /Kernels/ArrayTimeDerivative
+
+!syntax children /Kernels/ArrayTimeDerivative

framework/include/kernels/ArrayTimeDerivative.h

@@ -0,0 +1,38 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#pragma once
+
+#include "ArrayTimeKernel.h"
+
+// Forward Declaration
+class ArrayTimeDerivative;
+
+template <>
+InputParameters validParams<ArrayTimeDerivative>();
+
+class ArrayTimeDerivative : public ArrayTimeKernel
+{
+public:
+  ArrayTimeDerivative(const InputParameters & parameters);
+
+protected:
+  virtual RealEigenVector computeQpResidual() override;
+  virtual RealEigenVector computeQpJacobian() override;
+  virtual RealEigenMatrix computeQpOffDiagJacobian(MooseVariableFEBase & jvar) override;
+
+  /// time derivative coefficient type
+  unsigned int _coeff_type;
+  /// scalar time derivative coefficient
+  const MaterialProperty<Real> * _coeff;
+  /// array time derivative coefficient
+  const MaterialProperty<RealEigenVector> * _coeff_array;
+  /// matrix time derivative coefficient
+  const MaterialProperty<RealEigenMatrix> * _coeff_2d_array;
+};

framework/include/kernels/ArrayTimeKernel.h

@@ -0,0 +1,39 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#pragma once
+
+#include "ArrayKernel.h"
+
+// Forward Declaration
+class ArrayTimeKernel;
+
+template <>
+InputParameters validParams<ArrayTimeKernel>();
+
+/**
+ * All array time kernels should inherit from this class
+ *
+ */
+class ArrayTimeKernel : public ArrayKernel
+{
+public:
+  ArrayTimeKernel(const InputParameters & parameters);
+
+protected:
+  /// Time derivative of \vec{u}
+  const ArrayVariableValue & _u_dot;
+
+  /**
+   * Derivative of u_dot with respect to u. This value is only dependent on the
+   * time integration scheme, which is the same for every componenet. Therefore,
+   * _du_dot_du is VariableValue instead of an ArrayVariableValue
+   */
+  const VariableValue & _du_dot_du;
+};

framework/src/kernels/ArrayTimeDerivative.C

@@ -0,0 +1,91 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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 "ArrayTimeDerivative.h"
+
+registerMooseObject("MooseApp", ArrayTimeDerivative);
+
+template <>
+InputParameters
+validParams<ArrayTimeDerivative>()
+{
+  InputParameters params = validParams<ArrayTimeKernel>();
+  params.addClassDescription("Array time derivative operator with the weak form of $(\\psi_i, "
+                             "\\frac{\\partial u_h}{\\partial t})$.");
+  params.addParam<MaterialPropertyName>("time_derivative_coefficient",
+                                        "The name of the time derivative coefficient. "
+                                        "Can be scalar, vector, or matrix");
+  return params;
+}
+
+ArrayTimeDerivative::ArrayTimeDerivative(const InputParameters & parameters)
+  : ArrayTimeKernel(parameters)
+{
+  if (hasMaterialProperty<Real>("time_derivative_coefficient"))
+  {
+    _coeff_type = 0;
+    _coeff = &getMaterialProperty<Real>("time_derivative_coefficient");
+  }
+  else if (hasMaterialProperty<RealEigenVector>("time_derivative_coefficient"))
+  {
+    _coeff_type = 1;
+    _coeff_array = &getMaterialProperty<RealEigenVector>("time_derivative_coefficient");
+    mooseAssert((*_coeff_array)[_qp].size() == _var.count(),
+                "time_derivative_coefficient size is inconsistent with the number of components "
+                "in array variable");
+  }
+  else if (hasMaterialProperty<RealEigenMatrix>("time_derivative_coefficient"))
+  {
+    _coeff_type = 2;
+    _coeff_2d_array = &getMaterialProperty<RealEigenMatrix>("time_derivative_coefficient");
+    mooseAssert((*_coeff_2d_array)[_qp].cols() == _var.count(),
+                "time_derivative_coefficient size is inconsistent with the number of components "
+                "in array variable");
+    mooseAssert((*_coeff_2d_array)[_qp].rows() == _var.count(),
+                "time_derivative_coefficient size is inconsistent with the number of components "
+                "in array variable");
+  }
+  else
+  {
+    MaterialPropertyName mat = getParam<MaterialPropertyName>("time_derivative_coefficient");
+    mooseError("Property " + mat + " is of unsupported type for ArrayTimeDerivative");
+  }
+}
+
+RealEigenVector
+ArrayTimeDerivative::computeQpResidual()
+{
+  if (_coeff_type == 0)
+    return (*_coeff)[_qp] * _u_dot[_qp] * _test[_i][_qp];
+  else if (_coeff_type == 1)
+    return ((*_coeff_array)[_qp].array() * _u_dot[_qp].array()) * _test[_i][_qp];
+  else
+    return (*_coeff_2d_array)[_qp] * _u_dot[_qp] * _test[_i][_qp];
+}
+
+RealEigenVector
+ArrayTimeDerivative::computeQpJacobian()
+{
+  Real tmp = _test[_i][_qp] * _phi[_j][_qp] * _du_dot_du[_qp];
+  if (_coeff_type == 0)
+    return RealEigenVector::Constant(_var.count(), tmp * (*_coeff)[_qp]);
+  else if (_coeff_type == 1)
+    return tmp * (*_coeff_array)[_qp];
+  else
+    return tmp * (*_coeff_2d_array)[_qp].diagonal();
+}
+
+RealEigenMatrix
+ArrayTimeDerivative::computeQpOffDiagJacobian(MooseVariableFEBase & jvar)
+{
+  if (jvar.number() == _var.number() && _coeff_type == 2)
+    return _phi[_j][_qp] * _test[_i][_qp] * _du_dot_du[_qp] * (*_coeff_2d_array)[_qp];
+  else
+    return ArrayKernel::computeQpOffDiagJacobian(jvar);
+}

framework/src/kernels/ArrayTimeKernel.C

@@ -0,0 +1,27 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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 "ArrayTimeKernel.h"
+
+template <>
+InputParameters
+validParams<ArrayTimeKernel>()
+{
+  InputParameters params = validParams<ArrayKernel>();
+
+  params.set<MultiMooseEnum>("vector_tags") = "time";
+  params.set<MultiMooseEnum>("matrix_tags") = "system time";
+
+  return params;
+}
+
+ArrayTimeKernel::ArrayTimeKernel(const InputParameters & parameters)
+  : ArrayKernel(parameters), _u_dot(_var.uDot()), _du_dot_du(_var.duDotDu())
+{
+}

test/tests/kernels/array_kernels/array_diffusion_reaction_transient.i

@@ -0,0 +1,91 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 4
+  ny = 4
+[]
+
+[Variables]
+  [u]
+    order = FIRST
+    family = LAGRANGE
+    components = 2
+  []
+[]
+
+[Kernels]
+  [dudt]
+    type = ArrayTimeDerivative
+    variable = u
+    time_derivative_coefficient = tc
+  []
+  [diff]
+    type = ArrayDiffusion
+    variable = u
+    diffusion_coefficient = dc
+  []
+  [reaction]
+    type = ArrayReaction
+    variable = u
+    reaction_coefficient_type = full
+    reaction_coefficient = rc
+  []
+[]
+
+[BCs]
+  [left]
+    type = ArrayDirichletBC
+    variable = u
+    boundary = 1
+    values = '0 0'
+  []
+
+  [right]
+    type = ArrayDirichletBC
+    variable = u
+    boundary = 2
+    values = '1 2'
+  []
+[]
+
+[Materials]
+  [tc]
+    type = GenericConstantArray
+    prop_name = tc
+    prop_value = '1 1'
+  []
+  [dc]
+    type = GenericConstantArray
+    prop_name = dc
+    prop_value = '1 1'
+  []
+  [rc]
+    type = GenericConstant2DArray
+    prop_name = rc
+    prop_value = '1 0; -0.1 1'
+  []
+[]
+
+[Postprocessors]
+  [intu0]
+    type = ElementIntegralArrayVariablePostprocessor
+    variable = u
+    component = 0
+  []
+  [intu1]
+    type = ElementIntegralArrayVariablePostprocessor
+    variable = u
+    component = 1
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+  dt = 0.1
+  num_steps = 10
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/kernels/array_kernels/tests

@@ -1,5 +1,5 @@
 [Tests]
-  issues = '#6881'
+  issues = '#6881 #14057'
   design = '/variables/ArrayMooseVariable.md'
   [./test]
     type = 'Exodiff'
@@ -45,4 +45,10 @@
     requirement = 'MOOSE shall provide residual save-in and Jacobian diagonal save-in on an array variable.'
     prereq = test_save_in
   [../]
+  [./test_diffusion_reaction_transient]
+    type = 'Exodiff'
+    input = 'array_diffusion_reaction_transient.i'
+    exodiff = 'array_diffusion_reaction_transient_out.e'
+    requirement = 'MOOSE shall provide array time derivative kernels on an array variable.'
+  [../]
 []