From 551334f328ae4c293bc7279113b8c0f05e598a36 Mon Sep 17 00:00:00 2001 From: Andrew E Slaughter Date: Wed, 13 Mar 2019 11:15:38 -0600 Subject: [PATCH] Create a generalized stabilization ADKernel (closes #13044) --- framework/include/kernels/ADKernelSUPG.h | 20 +- .../include/kernels/ADKernelStabilized.h | 54 ++++++ framework/src/kernels/ADKernelSUPG.C | 151 +-------------- framework/src/kernels/ADKernelStabilized.C | 176 ++++++++++++++++++ 4 files changed, 241 insertions(+), 160 deletions(-) create mode 100644 framework/include/kernels/ADKernelStabilized.h create mode 100644 framework/src/kernels/ADKernelStabilized.C diff --git a/framework/include/kernels/ADKernelSUPG.h b/framework/include/kernels/ADKernelSUPG.h index 1be204569d6b..f4c1ca226dbf 100644 --- a/framework/include/kernels/ADKernelSUPG.h +++ b/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::_velocity; \ using ADKernelSUPGTempl::_tau #define usingKernelSUPGMembers usingTemplKernelSUPGMembers(Real) @@ -31,28 +31,18 @@ declareADValidParams(ADKernelSUPG); declareADValidParams(ADVectorKernelSUPG); template -class ADKernelSUPGTempl : public ADKernelTempl +class ADKernelSUPGTempl : public ADKernelStabilizedTempl { 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::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 */ diff --git a/framework/include/kernels/ADKernelStabilized.h b/framework/include/kernels/ADKernelStabilized.h new file mode 100644 index 000000000000..2090d0c8f854 --- /dev/null +++ b/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 +class ADKernelStabilizedTempl; + +template +using ADKernelStabilized = ADKernelStabilizedTempl; +template +using ADVectorKernelStabilized = ADKernelStabilizedTempl; + +declareADValidParams(ADKernelStabilized); +declareADValidParams(ADVectorKernelStabilized); + +template +class ADKernelStabilizedTempl : public ADKernelTempl +{ +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::type>::OutputValue + precomputeQpStrongResidual() = 0; + + virtual ADRealVectorValue computeQpStabilization() = 0; + + virtual ADResidual computeQpResidual() override final; + + usingTemplKernelMembers(T); +}; + +#endif /* ADKERNELSTABILIZED_H */ diff --git a/framework/src/kernels/ADKernelSUPG.C b/framework/src/kernels/ADKernelSUPG.C index 5766de7a4946..fd642566fb56 100644 --- a/framework/src/kernels/ADKernelSUPG.C +++ b/framework/src/kernels/ADKernelSUPG.C @@ -16,164 +16,25 @@ #define PGParams \ params.addParam( \ - "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 ADKernelSUPGTempl::ADKernelSUPGTempl(const InputParameters & parameters) - : ADKernelTempl(parameters), + : ADKernelStabilizedTempl(parameters), _tau(adGetADMaterialProperty("tau_name")), _velocity(adCoupledVectorValue("velocity")) { } template -void -ADKernelSUPGTempl::computeResidual() +ADRealVectorValue +ADKernelSUPGTempl::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::computeResidual() -{ -} - -template <> -void -ADKernelSUPGTempl::computeResidual() -{ -} - -template -void -ADKernelSUPGTempl::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 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::computeJacobian() -{ -} - -template <> -void -ADKernelSUPGTempl::computeJacobian() -{ -} - -template -void -ADKernelSUPGTempl::computeADOffDiagJacobian() -{ - std::vector 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> & 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::computeADOffDiagJacobian() -{ -} - -template <> -void -ADKernelSUPGTempl::computeADOffDiagJacobian() -{ -} - -template -ADResidual -ADKernelSUPGTempl::computeQpResidual() -{ - mooseError("Override precomputeQpStrongResidual() in your ADKernelSUPG derived class!"); + return _velocity[_qp] * _tau[_qp]; } template class ADKernelSUPGTempl; diff --git a/framework/src/kernels/ADKernelStabilized.C b/framework/src/kernels/ADKernelStabilized.C new file mode 100644 index 000000000000..293f059c9188 --- /dev/null +++ b/framework/src/kernels/ADKernelStabilized.C @@ -0,0 +1,176 @@ +//* 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 "ADKernelStabilized.h" +#include "MathUtils.h" +#include "Assembly.h" + +// libmesh includes +#include "libmesh/threads.h" + +defineADValidParams(ADKernelStabilized, ADKernel, ); +defineADValidParams(ADVectorKernelStabilized, ADVectorKernel, ); + +template +ADKernelStabilizedTempl::ADKernelStabilizedTempl( + const InputParameters & parameters) + : ADKernelTempl(parameters) +{ +} + +template +void +ADKernelStabilizedTempl::computeResidual() +{ + 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]; + for (_i = 0; _i < n_test; _i++) // target for auto vectorization + _local_re(_i) += _grad_test[_i][_qp] * computeQpStabilization() * 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 +ADKernelStabilizedTempl::computeResidual() +{ +} + +template <> +void +ADKernelStabilizedTempl::computeResidual() +{ +} + +template +void +ADKernelStabilizedTempl::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]; + for (_i = 0; _i < _grad_test.size(); _i++) + { + const auto residual = _grad_test[_i][_qp] * computeQpStabilization() * 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 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 +ADKernelStabilizedTempl::computeJacobian() +{ +} + +template <> +void +ADKernelStabilizedTempl::computeJacobian() +{ +} + +template +void +ADKernelStabilizedTempl::computeADOffDiagJacobian() +{ + std::vector residuals(_grad_test.size(), 0); + + precalculateResidual(); + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + const auto value = precomputeQpStrongResidual() * _ad_JxW[_qp] * _ad_coord[_qp]; + for (_i = 0; _i < _grad_test.size(); _i++) + residuals[_i] += _grad_test[_i][_qp] * computeQpStabilization() * value; + } + + std::vector> & 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 +ADKernelStabilizedTempl::computeADOffDiagJacobian() +{ +} + +template <> +void +ADKernelStabilizedTempl::computeADOffDiagJacobian() +{ +} + +template +ADResidual +ADKernelStabilizedTempl::computeQpResidual() +{ + mooseError("Override precomputeQpStrongResidual() in your ADKernelStabilized derived class!"); +} + +template class ADKernelStabilizedTempl; +template class ADKernelStabilizedTempl; +template class ADKernelStabilizedTempl; +template class ADKernelStabilizedTempl;