Merge pull request #13057 from aeslaughter/ad-kernel-stabilized-13044

Create a generalized stabilization ADKernel

framework/include/kernels/ADKernelSUPG.h

@@ -10,10 +10,10 @@
 #ifndef ADKERNELSUPG_H
 #define ADKERNELSUPG_H
 
-#include "ADKernel.h"
+#include "ADKernelStabilized.h"
 
 #define usingTemplKernelSUPGMembers(type)                                                          \
-  usingTemplKernelMembers(type);                                                                   \
+  usingTemplKernelStabilizedMembers(type);                                                         \
   using ADKernelSUPGTempl<type, compute_stage>::_velocity;                                         \
   using ADKernelSUPGTempl<type, compute_stage>::_tau
 #define usingKernelSUPGMembers usingTemplKernelSUPGMembers(Real)
@@ -31,28 +31,18 @@ declareADValidParams(ADKernelSUPG);
 declareADValidParams(ADVectorKernelSUPG);
 
 template <typename T, ComputeStage compute_stage>
-class ADKernelSUPGTempl : public ADKernelTempl<T, compute_stage>
+class ADKernelSUPGTempl : public ADKernelStabilizedTempl<T, compute_stage>
 {
 public:
   ADKernelSUPGTempl(const InputParameters & parameters);
 
-  virtual void computeResidual() override;
-  virtual void computeJacobian() override;
-  virtual void computeADOffDiagJacobian() override;
-
 protected:
-  /**
-   * Called before forming the residual for an element
-   */
-  virtual typename OutputTools<typename Moose::ValueType<T, compute_stage>::type>::OutputValue
-  precomputeQpStrongResidual() = 0;
-
-  virtual ADResidual computeQpResidual() override final;
+  ADRealVectorValue virtual computeQpStabilization() override;
 
   const ADMaterialProperty(Real) & _tau;
   const ADVectorVariableValue & _velocity;
 
-  usingTemplKernelMembers(T);
+  usingTemplKernelStabilizedMembers(T);
 };
 
 #endif /* ADKERNELSUPG_H */

framework/include/kernels/ADKernelStabilized.h

@@ -0,0 +1,54 @@
+//* 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
+
+#ifndef ADKERNELSTABILIZED_H
+#define ADKERNELSTABILIZED_H
+
+#include "ADKernel.h"
+
+#define usingTemplKernelStabilizedMembers(type) usingTemplKernelMembers(type)
+#define usingKernelStabilizedMembers usingTemplKernelStabilizedMembers(Real)
+#define usingVectorKernelStabilizedMembers usingTemplKernelStabilizedMembers(RealVectorValue)
+
+template <typename, ComputeStage>
+class ADKernelStabilizedTempl;
+
+template <ComputeStage compute_stage>
+using ADKernelStabilized = ADKernelStabilizedTempl<Real, compute_stage>;
+template <ComputeStage compute_stage>
+using ADVectorKernelStabilized = ADKernelStabilizedTempl<RealVectorValue, compute_stage>;
+
+declareADValidParams(ADKernelStabilized);
+declareADValidParams(ADVectorKernelStabilized);
+
+template <typename T, ComputeStage compute_stage>
+class ADKernelStabilizedTempl : public ADKernelTempl<T, compute_stage>
+{
+public:
+  ADKernelStabilizedTempl(const InputParameters & parameters);
+
+  virtual void computeResidual() override;
+  virtual void computeJacobian() override;
+  virtual void computeADOffDiagJacobian() override;
+
+protected:
+  /**
+   * Called before forming the residual for an element
+   */
+  virtual typename OutputTools<typename Moose::ValueType<T, compute_stage>::type>::OutputValue
+  precomputeQpStrongResidual() = 0;
+
+  virtual ADRealVectorValue computeQpStabilization() = 0;
+
+  virtual ADResidual computeQpResidual() override final;
+
+  usingTemplKernelMembers(T);
+};
+
+#endif /* ADKERNELSTABILIZED_H */

framework/src/kernels/ADKernelSUPG.C

@@ -16,164 +16,25 @@
 
 #define PGParams                                                                                   \
   params.addParam<MaterialPropertyName>(                                                           \
-      "tau_name", "tau", "The name of the stabilization parameter tau");                           \
-  params.addRequiredCoupledVar("velocity", "The velocity variable")
+      "tau_name", "tau", "The name of the stabilization parameter tau.");                          \
+  params.addRequiredCoupledVar("velocity", "The velocity variable.")
 
 defineADValidParams(ADKernelSUPG, ADKernel, PGParams;);
 defineADValidParams(ADVectorKernelSUPG, ADVectorKernel, PGParams;);
 
 template <typename T, ComputeStage compute_stage>
 ADKernelSUPGTempl<T, compute_stage>::ADKernelSUPGTempl(const InputParameters & parameters)
-  : ADKernelTempl<T, compute_stage>(parameters),
+  : ADKernelStabilizedTempl<T, compute_stage>(parameters),
     _tau(adGetADMaterialProperty<Real>("tau_name")),
     _velocity(adCoupledVectorValue("velocity"))
 {
 }
 
 template <typename T, ComputeStage compute_stage>
-void
-ADKernelSUPGTempl<T, compute_stage>::computeResidual()
+ADRealVectorValue
+ADKernelSUPGTempl<T, compute_stage>::computeQpStabilization()
 {
-  prepareVectorTag(_assembly, _var.number());
-
-  precalculateResidual();
-  const unsigned int n_test = _grad_test.size();
-  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
-  {
-    const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp] * _tau[_qp];
-    for (_i = 0; _i < n_test; _i++) // target for auto vectorization
-      _local_re(_i) += _grad_test[_i][_qp] * _velocity[_qp] * value;
-  }
-
-  accumulateTaggedLocalResidual();
-
-  if (_has_save_in)
-  {
-    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
-    for (unsigned int i = 0; i < _save_in.size(); i++)
-      _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
-  }
-}
-
-template <>
-void
-ADKernelSUPGTempl<Real, JACOBIAN>::computeResidual()
-{
-}
-
-template <>
-void
-ADKernelSUPGTempl<RealVectorValue, JACOBIAN>::computeResidual()
-{
-}
-
-template <typename T, ComputeStage compute_stage>
-void
-ADKernelSUPGTempl<T, compute_stage>::computeJacobian()
-{
-  prepareMatrixTag(_assembly, _var.number(), _var.number());
-
-  size_t ad_offset = _var.number() * _sys.getMaxVarNDofsPerElem();
-
-  precalculateResidual();
-  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
-  {
-    // This will also compute the derivative with respect to all dofs
-    const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp] * _tau[_qp];
-    for (_i = 0; _i < _grad_test.size(); _i++)
-    {
-      const auto residual = _grad_test[_i][_qp] * _velocity[_qp] * value;
-      for (_j = 0; _j < _var.phiSize(); _j++)
-        _local_ke(_i, _j) += residual.derivatives()[ad_offset + _j];
-    }
-  }
-  accumulateTaggedLocalMatrix();
-
-  if (_has_diag_save_in)
-  {
-    unsigned int rows = _local_ke.m();
-    DenseVector<Number> diag(rows);
-    for (unsigned int i = 0; i < rows; i++)
-      diag(i) = _local_ke(i, i);
-
-    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
-    for (unsigned int i = 0; i < _diag_save_in.size(); i++)
-      _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
-  }
-}
-
-template <>
-void
-ADKernelSUPGTempl<Real, RESIDUAL>::computeJacobian()
-{
-}
-
-template <>
-void
-ADKernelSUPGTempl<RealVectorValue, RESIDUAL>::computeJacobian()
-{
-}
-
-template <typename T, ComputeStage compute_stage>
-void
-ADKernelSUPGTempl<T, compute_stage>::computeADOffDiagJacobian()
-{
-  std::vector<DualReal> residuals(_grad_test.size(), 0);
-
-  precalculateResidual();
-  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
-  {
-    const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp] * _tau[_qp];
-    for (_i = 0; _i < _grad_test.size(); _i++)
-      residuals[_i] += _grad_test[_i][_qp] * _velocity[_qp] * value;
-  }
-
-  std::vector<std::pair<MooseVariableFEBase *, MooseVariableFEBase *>> & ce =
-      _assembly.couplingEntries();
-  for (const auto & it : ce)
-  {
-    MooseVariableFEBase & ivariable = *(it.first);
-    MooseVariableFEBase & jvariable = *(it.second);
-
-    unsigned int ivar = ivariable.number();
-    unsigned int jvar = jvariable.number();
-
-    if (ivar != _var.number())
-      continue;
-
-    size_t ad_offset = jvar * _sys.getMaxVarNDofsPerElem();
-
-    prepareMatrixTag(_assembly, ivar, jvar);
-
-    if (_local_ke.m() != _grad_test.size() || _local_ke.n() != jvariable.phiSize())
-      continue;
-
-    precalculateResidual();
-    for (_i = 0; _i < _grad_test.size(); _i++)
-      for (_j = 0; _j < jvariable.phiSize(); _j++)
-        _local_ke(_i, _j) += residuals[_i].derivatives()[ad_offset + _j];
-
-    accumulateTaggedLocalMatrix();
-  }
-}
-
-template <>
-void
-ADKernelSUPGTempl<Real, RESIDUAL>::computeADOffDiagJacobian()
-{
-}
-
-template <>
-void
-ADKernelSUPGTempl<RealVectorValue, RESIDUAL>::computeADOffDiagJacobian()
-{
-}
-
-template <typename T, ComputeStage compute_stage>
-ADResidual
-ADKernelSUPGTempl<T, compute_stage>::computeQpResidual()
-{
-  mooseError("Override precomputeQpStrongResidual() in your ADKernelSUPG derived class!");
+  return _velocity[_qp] * _tau[_qp];
 }
 
 template class ADKernelSUPGTempl<Real, RESIDUAL>;