Intra-variable Jacobian (#32212)

Author: dschwen

framework/include/nodalkernels/ADArrayNodalKernel.h

@@ -66,11 +66,16 @@ class ADArrayNodalKernel : public NodalKernelBase, public MooseVariableInterface
 
   /**
    * Dummy method so we can make derived generic classes that template on <bool is_ad> */
-  virtual RealEigenVector computeQpJacobian()
+  virtual void computeQpJacobian()
   {
     mooseError("I'm an AD object, so computeQpJacobian should never be called");
   }
 
+  void setJacobian(unsigned int, unsigned int, Real)
+  {
+    mooseError("I'm an AD object, so setJacobian should never be called");
+  }
+
   /// variable this works on
   MooseVariableFE<RealEigenVector> & _var;
 

framework/include/nodalkernels/ArrayNodalKernel.h

@@ -65,20 +65,21 @@ class ArrayNodalKernel : public NodalKernelBase, public MooseVariableInterface<R
   virtual void computeQpResidual(RealEigenVector & residual) = 0;
 
   /**
-   * The user can override this function to compute the "on-diagonal"
-   * Jacobian contribution.  If not overriden,
-   * returns an array of 1s.
-   * @return On-diagonal Jacobian contribution
+   * The user can override this function to compute the intra-variable
+   * off-diagonal Jacobian contribution (the coupling between array
+   * components of the variable itself). Use setJacobian in this function.
    */
-  virtual RealEigenVector computeQpJacobian();
+  virtual void computeQpJacobian() {}
 
   /**
-   * This is the virtual that derived classes should override for
-   * computing an off-diagonal jacobian component.
-   * @param Coupled variable for the off-diagonal Jacobian contribution
-   * @return Off-diagonal Jacobian contribution
+   * The user can override this function to compute the inter-variable
+   * off-diagonal Jacobian contribution (the coupling between array
+   * components of of the kernel variable and the coupled variable jvar).
+   * Use setJacobian in this function.
    */
-  virtual RealEigenMatrix computeQpOffDiagJacobian(const MooseVariableFieldBase & jvar);
+  virtual void computeQpOffDiagJacobian(unsigned int /*jvar*/) {}
+
+  void setJacobian(unsigned int i, unsigned int j, Real value);
 
   /// variable this works on
   MooseVariableFE<RealEigenVector> & _var;
@@ -92,4 +93,11 @@ class ArrayNodalKernel : public NodalKernelBase, public MooseVariableInterface<R
 private:
   /// Work vector for residual
   RealEigenVector _work_vector;
+
+  /// DOF indices
+  const std::vector<dof_id_type> * _ivar_indices;
+  const std::vector<dof_id_type> * _jvar_indices;
+
+  /// scaling factors
+  const std::vector<Real> & _scaling;
 };

framework/include/nodalkernels/ArrayReactionNodalKernel.h

@@ -24,7 +24,7 @@ class ArrayReactionNodalKernelTempl : public GenericArrayNodalKernel<is_ad>
 
 protected:
   virtual void computeQpResidual(GenericRealEigenVector<is_ad> & residual) override;
-  virtual RealEigenVector computeQpJacobian() override;
+  virtual void computeQpJacobian() override;
 
   /// rate coefficient
   const RealEigenVector & _coeff;

framework/include/nodalkernels/ArrayTimeDerivativeNodalKernel.h

@@ -23,7 +23,7 @@ class ArrayTimeDerivativeNodalKernel : public ArrayNodalKernel
 
 protected:
   virtual void computeQpResidual(RealEigenVector & residual) override;
-  virtual RealEigenVector computeQpJacobian() override;
+  virtual void computeQpJacobian() override;
 
   /// Time derivative of u
   const ArrayVariableValue & _u_dot;

framework/include/nodalkernels/GenericArrayNodalKernel.h

@@ -35,4 +35,5 @@ class GenericArrayNodalKernel<true> : public ADArrayNodalKernel
   using GenericArrayNodalKernel<is_ad>::_count;                                                    \
   using GenericArrayNodalKernel<is_ad>::_u;                                                        \
   using GenericArrayNodalKernel<is_ad>::_qp;                                                       \
+  using GenericArrayNodalKernel<is_ad>::setJacobian;                                               \
   using GenericArrayNodalKernel<is_ad>::paramError

framework/src/nodalkernels/ArrayNodalKernel.C

@@ -29,7 +29,8 @@ ArrayNodalKernel::ArrayNodalKernel(const InputParameters & parameters)
     _var(*mooseVariable()),
     _u(_var.dofValues()),
     _count(_var.count()),
-    _work_vector(_count)
+    _work_vector(_count),
+    _scaling(_var.arrayScalingFactor())
 {
   addMooseVariableDependency(mooseVariable());
 }
@@ -41,7 +42,7 @@ ArrayNodalKernel::computeResidual()
     return;
   _qp = 0;
   computeQpResidual(_work_vector);
-  addResiduals(_assembly, _work_vector, _var.dofIndices(), _var.arrayScalingFactor());
+  addResiduals(_assembly, _work_vector, _var.dofIndices(), _scaling);
 }
 
 void
@@ -50,12 +51,12 @@ ArrayNodalKernel::computeJacobian()
   if (!_var.isNodalDefined())
     return;
   _qp = 0;
-  const auto jacobian = computeQpJacobian();
-  const auto & dof_indices = _var.dofIndices();
-  mooseAssert(dof_indices.size() == _count, "The number of dofs should be equal to count");
-  for (const auto i : make_range(_count))
-    addJacobianElement(
-        _assembly, jacobian(i), dof_indices[i], dof_indices[i], _var.arrayScalingFactor()[i]);
+
+  _ivar_indices = &_var.dofIndices();
+  _jvar_indices = &_var.dofIndices();
+  mooseAssert(_ivar_indices->size() == _count, "The number of dofs should be equal to count");
+
+  computeQpJacobian();
 }
 
 void
@@ -69,31 +70,18 @@ ArrayNodalKernel::computeOffDiagJacobian(const unsigned int jvar_num)
   else
   {
     const auto & jvar = getVariable(jvar_num);
-    const auto jacobian = computeQpOffDiagJacobian(jvar);
-    const auto & ivar_indices = _var.dofIndices();
-    const auto & jvar_indices = jvar.dofIndices();
+    _ivar_indices = &_var.dofIndices();
+    _jvar_indices = &jvar.dofIndices();
+    mooseAssert(_ivar_indices->size() == _count, "The number of dofs should be equal to count");
+    mooseAssert(_jvar_indices->size() == jvar.count(),
+                "The number of dofs should be equal to count");
 
-    mooseAssert(ivar_indices.size() == _count, "The number of dofs should be equal to count");
-    mooseAssert(jvar_indices.size() == jvar.count(), "The number of dofs should be equal to count");
-
-    for (const auto i : make_range(_var.count()))
-      for (const auto j : make_range(jvar.count()))
-        addJacobianElement(_assembly,
-                           jacobian(i, j),
-                           ivar_indices[i],
-                           jvar_indices[j],
-                           _var.arrayScalingFactor()[i]);
+    computeQpOffDiagJacobian(jvar_num);
   }
 }
 
-RealEigenVector
-ArrayNodalKernel::computeQpJacobian()
-{
-  return RealEigenVector::Zero(_count);
-}
-
-RealEigenMatrix
-ArrayNodalKernel::computeQpOffDiagJacobian(const MooseVariableFieldBase & jvar)
+void
+ArrayNodalKernel::setJacobian(unsigned int i, unsigned int j, Real value)
 {
-  return RealEigenMatrix::Zero(_count, jvar.count());
+  addJacobianElement(_assembly, value, (*_ivar_indices)[i], (*_jvar_indices)[j], _scaling[i]);
 }

framework/src/nodalkernels/ArrayReactionNodalKernel.C

@@ -43,10 +43,11 @@ ArrayReactionNodalKernelTempl<is_ad>::computeQpResidual(GenericRealEigenVector<i
 }
 
 template <bool is_ad>
-RealEigenVector
+void
 ArrayReactionNodalKernelTempl<is_ad>::computeQpJacobian()
 {
-  return _coeff;
+  for (const auto i : index_range(_coeff))
+    setJacobian(i, i, _coeff(i));
 }
 
 template class ArrayReactionNodalKernelTempl<false>;

framework/src/nodalkernels/ArrayTimeDerivativeNodalKernel.C

@@ -34,8 +34,9 @@ ArrayTimeDerivativeNodalKernel::computeQpResidual(RealEigenVector & residual)
   residual = _u_dot[_qp];
 }
 
-RealEigenVector
+void
 ArrayTimeDerivativeNodalKernel::computeQpJacobian()
 {
-  return RealEigenVector::Constant(_count, _du_dot_du[_qp]);
+  for (const auto i : make_range(_count))
+    setJacobian(i, i, _du_dot_du[_qp]);
 }