Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ADInterfaceKernel.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "MooseVariableFE.h"
      15             : #include "SystemBase.h"
      16             : #include "ADUtils.h"
      17             : 
      18             : // libmesh includes
      19             : #include "libmesh/quadrature.h"
      20             : 
      21             : template <typename T>
      22             : InputParameters
      23      100085 : ADInterfaceKernelTempl<T>::validParams()
      24             : {
      25      100085 :   InputParameters params = InterfaceKernelBase::validParams();
      26             :   if (std::is_same<T, Real>::value)
      27       71511 :     params.registerBase("InterfaceKernel");
      28             :   else if (std::is_same<T, RealVectorValue>::value)
      29       28574 :     params.registerBase("VectorInterfaceKernel");
      30             :   else
      31             :     ::mooseError("unsupported ADInterfaceKernelTempl specialization");
      32      100085 :   return params;
      33           0 : }
      34             : 
      35             : template <typename T>
      36         124 : ADInterfaceKernelTempl<T>::ADInterfaceKernelTempl(const InputParameters & parameters)
      37             :   : InterfaceKernelBase(parameters),
      38             :     NeighborMooseVariableInterface<T>(this,
      39             :                                       false,
      40             :                                       Moose::VarKindType::VAR_SOLVER,
      41             :                                       std::is_same<T, Real>::value
      42             :                                           ? Moose::VarFieldType::VAR_FIELD_STANDARD
      43             :                                           : Moose::VarFieldType::VAR_FIELD_VECTOR),
      44         248 :     _var(*this->mooseVariable()),
      45         124 :     _normals(_assembly.normals()),
      46         124 :     _u(_var.adSln()),
      47         124 :     _grad_u(_var.adGradSln()),
      48         124 :     _ad_JxW(_assembly.adJxWFace()),
      49         124 :     _ad_coord(_assembly.adCoordTransformation()),
      50         124 :     _ad_q_point(_assembly.adQPoints()),
      51         124 :     _phi(_assembly.phiFace(_var)),
      52         124 :     _test(_var.phiFace()),
      53         124 :     _grad_test(_var.adGradPhiFace()),
      54         124 :     _neighbor_var(*getVarHelper<MooseVariableFE<T>>("neighbor_var", 0)),
      55         124 :     _neighbor_value(_neighbor_var.adSlnNeighbor()),
      56         124 :     _grad_neighbor_value(_neighbor_var.adGradSlnNeighbor()),
      57         124 :     _phi_neighbor(_assembly.phiFaceNeighbor(_neighbor_var)),
      58         124 :     _test_neighbor(_neighbor_var.phiFaceNeighbor()),
      59         248 :     _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor())
      60             : 
      61             : {
      62         124 :   _subproblem.haveADObjects(true);
      63             : 
      64         124 :   addMooseVariableDependency(this->mooseVariable());
      65             : 
      66         124 :   if (!parameters.isParamValid("boundary"))
      67           0 :     mooseError(
      68             :         "In order to use an interface kernel, you must specify a boundary where it will live.");
      69         124 : }
      70             : 
      71             : template <typename T>
      72             : void
      73        2486 : ADInterfaceKernelTempl<T>::computeElemNeighResidual(Moose::DGResidualType type)
      74             : {
      75             :   bool is_elem;
      76        2486 :   if (type == Moose::Element)
      77        1243 :     is_elem = true;
      78             :   else
      79        1243 :     is_elem = false;
      80             : 
      81        2486 :   const ADTemplateVariableTestValue<T> & test_space = is_elem ? _test : _test_neighbor;
      82             : 
      83        2486 :   if (is_elem)
      84        1243 :     prepareVectorTag(_assembly, _var.number());
      85             :   else
      86        1243 :     prepareVectorTagNeighbor(_assembly, _neighbor_var.number());
      87             : 
      88        4972 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      89             :   {
      90        2486 :     initQpResidual(type);
      91        2486 :     const auto jxw_p = _JxW[_qp] * _coord[_qp];
      92        7458 :     for (_i = 0; _i < test_space.size(); _i++)
      93        4972 :       _local_re(_i) += jxw_p * raw_value(computeQpResidual(type));
      94             :   }
      95             : 
      96        2486 :   accumulateTaggedLocalResidual();
      97        2486 : }
      98             : 
      99             : template <typename T>
     100             : void
     101        1243 : ADInterfaceKernelTempl<T>::computeResidual()
     102             : {
     103             :   // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
     104             :   // associated only with the lower-dimensional geometric entity that is the boundary between two
     105             :   // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
     106             :   // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
     107             :   // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
     108             :   // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
     109             :   // to execute in the former case. In the future we should remove this and add some kind of "block"
     110             :   // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
     111             :   // this return
     112        2486 :   if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
     113        1243 :       !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
     114           0 :     return;
     115             : 
     116        1243 :   precalculateResidual();
     117             : 
     118             :   // Compute the residual for this element
     119        1243 :   computeElemNeighResidual(Moose::Element);
     120             : 
     121             :   // Compute the residual for the neighbor
     122        1243 :   computeElemNeighResidual(Moose::Neighbor);
     123             : }
     124             : 
     125             : template <typename T>
     126             : void
     127           0 : ADInterfaceKernelTempl<T>::computeElemNeighJacobian(Moose::DGJacobianType type)
     128             : {
     129             :   mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
     130             :               "With AD you should need one call per side");
     131             : 
     132           0 :   const ADTemplateVariableTestValue<T> & test_space =
     133           0 :       (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
     134             : 
     135           0 :   std::vector<ADReal> residuals(test_space.size(), 0);
     136             : 
     137           0 :   switch (type)
     138             :   {
     139           0 :     case Moose::ElementElement:
     140           0 :       resid_type = Moose::Element;
     141           0 :       break;
     142           0 :     case Moose::ElementNeighbor:
     143           0 :       resid_type = Moose::Element;
     144           0 :       break;
     145           0 :     case Moose::NeighborElement:
     146           0 :       resid_type = Moose::Neighbor;
     147           0 :       break;
     148           0 :     case Moose::NeighborNeighbor:
     149           0 :       resid_type = Moose::Neighbor;
     150           0 :       break;
     151           0 :     default:
     152           0 :       mooseError("Unknown DGJacobianType ", type);
     153             :   }
     154             : 
     155           0 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     156             :   {
     157           0 :     initQpResidual(resid_type);
     158           0 :     const auto jxw_c = _ad_JxW[_qp] * _ad_coord[_qp];
     159           0 :     for (_i = 0; _i < test_space.size(); _i++)
     160           0 :       residuals[_i] += jxw_c * computeQpResidual(resid_type);
     161             :   }
     162             : 
     163           0 :   const bool element_var_is_var = (type == Moose::ElementElement || type == Moose::ElementNeighbor);
     164           0 :   addJacobian(_assembly,
     165             :               residuals,
     166           0 :               element_var_is_var ? _var.dofIndices() : _neighbor_var.dofIndicesNeighbor(),
     167           0 :               element_var_is_var ? _var.scalingFactor() : _neighbor_var.scalingFactor());
     168           0 : }
     169             : 
     170             : template <typename T>
     171             : void
     172           0 : ADInterfaceKernelTempl<T>::computeJacobian()
     173             : {
     174             :   // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
     175             :   // associated only with the lower-dimensional geometric entity that is the boundary between two
     176             :   // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
     177             :   // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
     178             :   // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
     179             :   // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
     180             :   // to execute in the former case. In the future we should remove this and add some kind of "block"
     181             :   // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
     182             :   // this return
     183           0 :   if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
     184           0 :       !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
     185           0 :     return;
     186             : 
     187           0 :   precalculateJacobian();
     188             : 
     189           0 :   computeElemNeighJacobian(Moose::ElementElement);
     190           0 :   computeElemNeighJacobian(Moose::NeighborNeighbor);
     191             : }
     192             : 
     193             : template <typename T>
     194             : void
     195         298 : ADInterfaceKernelTempl<T>::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int)
     196             : {
     197             :   mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
     198             :               "With AD you should need one call per side");
     199             : 
     200         298 :   const ADTemplateVariableTestValue<T> & test_space =
     201             :       type == Moose::ElementElement ? _test : _test_neighbor;
     202             : 
     203         298 :   if (type == Moose::ElementElement)
     204         149 :     resid_type = Moose::Element;
     205             :   else
     206         149 :     resid_type = Moose::Neighbor;
     207             : 
     208         298 :   std::vector<ADReal> residuals(test_space.size(), 0);
     209             : 
     210         596 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     211             :   {
     212         298 :     initQpResidual(resid_type);
     213         298 :     const auto jxw_c = _ad_JxW[_qp] * _ad_coord[_qp];
     214         894 :     for (_i = 0; _i < test_space.size(); _i++)
     215         596 :       residuals[_i] += jxw_c * computeQpResidual(resid_type);
     216             :   }
     217             : 
     218             :   // We assert earlier that the type cannot be Moose::ElementNeighbor (nor Moose::NeighborElement)
     219         894 :   addJacobian(_assembly,
     220             :               residuals,
     221         149 :               type == Moose::ElementElement ? _var.dofIndices()
     222         149 :                                             : _neighbor_var.dofIndicesNeighbor(),
     223         298 :               type == Moose::ElementElement ? _var.scalingFactor() : _neighbor_var.scalingFactor());
     224         298 : }
     225             : 
     226             : template <typename T>
     227             : void
     228         343 : ADInterfaceKernelTempl<T>::computeElementOffDiagJacobian(unsigned int jvar)
     229             : {
     230             :   // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
     231             :   // associated only with the lower-dimensional geometric entity that is the boundary between two
     232             :   // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
     233             :   // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
     234             :   // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
     235             :   // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
     236             :   // to execute in the former case. In the future we should remove this and add some kind of "block"
     237             :   // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
     238             :   // this return
     239         686 :   if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
     240         343 :       !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
     241           0 :     return;
     242             : 
     243         343 :   if (jvar != _var.number())
     244             :     // We only need to do these computations a single time because AD computes all the derivatives
     245             :     // at once
     246         194 :     return;
     247             : 
     248         149 :   precalculateOffDiagJacobian(jvar);
     249             : 
     250             :   // Again AD does Jacobians all at once so we only need to call with ElementElement
     251         149 :   computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
     252             : }
     253             : 
     254             : template <typename T>
     255             : void
     256         343 : ADInterfaceKernelTempl<T>::computeNeighborOffDiagJacobian(unsigned int jvar)
     257             : {
     258             :   // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
     259             :   // associated only with the lower-dimensional geometric entity that is the boundary between two
     260             :   // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
     261             :   // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
     262             :   // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
     263             :   // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
     264             :   // to execute in the former case. In the future we should remove this and add some kind of "block"
     265             :   // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
     266             :   // this return
     267         686 :   if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
     268         343 :       !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
     269           0 :     return;
     270             : 
     271         343 :   if (jvar != _neighbor_var.number())
     272             :     // We only need to do these computations a single time because AD computes all the derivatives
     273             :     // at once
     274         194 :     return;
     275             : 
     276         149 :   precalculateOffDiagJacobian(jvar);
     277             : 
     278             :   // Again AD does Jacobians all at once so we only need to call with NeighborNeighbor
     279         149 :   computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
     280             : }
     281             : 
     282             : // Explicitly instantiates the two versions of the ADInterfaceKernelTempl class
     283             : template class ADInterfaceKernelTempl<Real>;
     284             : template class ADInterfaceKernelTempl<RealVectorValue>;