Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "MassFluxPenaltyIPHDG.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariableFE.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", MassFluxPenaltyIPHDG);
      16             : 
      17             : InputParameters
      18          80 : MassFluxPenaltyIPHDG::validParams()
      19             : {
      20          80 :   InputParameters params = HDGKernel::validParams();
      21         160 :   params.addRequiredParam<NonlinearVariableName>("u", "The x-velocity");
      22         160 :   params.addRequiredParam<NonlinearVariableName>("v", "The y-velocity");
      23         160 :   params.addRequiredParam<NonlinearVariableName>("u_face", "The x-velocity on the face");
      24         160 :   params.addRequiredParam<NonlinearVariableName>("v_face", "The y-velocity on the face");
      25         160 :   params.addRequiredRangeCheckedParam<unsigned short>(
      26             :       "component", "0<=component<=1", "The velocity component this object is being applied to");
      27         160 :   params.addParam<Real>("gamma", 1, "The penalty to multiply the jump");
      28          80 :   params.addClassDescription("introduces a jump correction on internal faces for grad-div "
      29             :                              "stabilization for discontinuous Galerkin methods. Because this is "
      30             :                              "derived from HDGKernel this class executes twice per face.");
      31          80 :   return params;
      32           0 : }
      33             : 
      34          40 : MassFluxPenaltyIPHDG::MassFluxPenaltyIPHDG(const InputParameters & parameters)
      35             :   : HDGKernel(parameters),
      36          40 :     _vel_x_var(_sys.getFieldVariable<Real>(_tid, getParam<NonlinearVariableName>("u"))),
      37          80 :     _vel_y_var(_sys.getFieldVariable<Real>(_tid, getParam<NonlinearVariableName>("v"))),
      38          80 :     _vel_x_face_var(_sys.getFieldVariable<Real>(_tid, getParam<NonlinearVariableName>("u_face"))),
      39          80 :     _vel_y_face_var(_sys.getFieldVariable<Real>(_tid, getParam<NonlinearVariableName>("v_face"))),
      40          40 :     _vel_x(_vel_x_var.adSln()),
      41          40 :     _vel_y(_vel_y_var.adSln()),
      42          40 :     _vel_x_face(_vel_x_face_var.adSln()),
      43          40 :     _vel_y_face(_vel_y_face_var.adSln()),
      44          40 :     _vel_x_phi(_vel_x_var.phiFace()),
      45          40 :     _vel_y_phi(_vel_y_var.phiFace()),
      46          40 :     _vel_x_face_phi(_vel_x_face_var.phiFace()),
      47          40 :     _vel_y_face_phi(_vel_y_face_var.phiFace()),
      48          80 :     _comp(getParam<unsigned short>("component")),
      49         120 :     _gamma(getParam<Real>("gamma"))
      50             : {
      51          40 :   if (_mesh.dimension() > 2)
      52           0 :     mooseError("Only two-dimensional velocity is currently implemented");
      53          40 : }
      54             : 
      55             : template <typename T>
      56             : void
      57     1224704 : MassFluxPenaltyIPHDG::computeOnSideHelper(std::vector<T> & residuals,
      58             :                                           const MooseArray<std::vector<Real>> & test,
      59             :                                           const Real sign)
      60             : {
      61     1224704 :   residuals.resize(test.size());
      62    10409984 :   for (auto & r : residuals)
      63     9185280 :     r = 0;
      64             : 
      65     3367936 :   auto return_residual = [this]() -> T
      66             :   {
      67             :     if constexpr (std::is_same<T, Real>::value)
      68     2525952 :       return MetaPhysicL::raw_value(computeQpResidualOnSide());
      69             :     else
      70     1148160 :       return computeQpResidualOnSide();
      71             :   };
      72             : 
      73     4898816 :   for (_qp = 0; _qp < _qrule_face->n_points(); ++_qp)
      74             :   {
      75     4822272 :     const auto qpres = _JxW_face[_qp] * _coord[_qp] * return_residual() * sign;
      76    31229952 :     for (const auto i : index_range(test))
      77    36167040 :       residuals[i] += qpres * test[i][_qp];
      78             :   }
      79     1224704 : }
      80             : 
      81             : void
      82      420992 : MassFluxPenaltyIPHDG::computeResidualOnSide()
      83             : {
      84      420992 :   const auto & var = _comp == 0 ? _vel_x_var : _vel_y_var;
      85      420992 :   const auto & face_var = _comp == 0 ? _vel_x_face_var : _vel_y_face_var;
      86      420992 :   const auto & test = _comp == 0 ? _vel_x_phi : _vel_y_phi;
      87      420992 :   const auto & face_test = _comp == 0 ? _vel_x_face_phi : _vel_y_face_phi;
      88      420992 :   computeOnSideHelper(_residuals, test, 1);
      89      420992 :   addResiduals(_assembly, _residuals, var.dofIndices(), var.scalingFactor());
      90      420992 :   computeOnSideHelper(_residuals, face_test, -1);
      91      420992 :   addResiduals(_assembly, _residuals, face_var.dofIndices(), face_var.scalingFactor());
      92      420992 : }
      93             : 
      94             : void
      95      191360 : MassFluxPenaltyIPHDG::computeJacobianOnSide()
      96             : {
      97      191360 :   const auto & var = _comp == 0 ? _vel_x_var : _vel_y_var;
      98      191360 :   const auto & face_var = _comp == 0 ? _vel_x_face_var : _vel_y_face_var;
      99      191360 :   const auto & test = _comp == 0 ? _vel_x_phi : _vel_y_phi;
     100      191360 :   const auto & face_test = _comp == 0 ? _vel_x_face_phi : _vel_y_face_phi;
     101      191360 :   computeOnSideHelper(_ad_residuals, test, 1);
     102      191360 :   addJacobian(_assembly, _ad_residuals, var.dofIndices(), var.scalingFactor());
     103      191360 :   computeOnSideHelper(_ad_residuals, face_test, -1);
     104      191360 :   addJacobian(_assembly, _ad_residuals, face_var.dofIndices(), face_var.scalingFactor());
     105      191360 : }
     106             : 
     107             : ADReal
     108     3674112 : MassFluxPenaltyIPHDG::computeQpResidualOnSide()
     109             : {
     110             :   const ADRealVectorValue soln_jump(
     111     7348224 :       _vel_x[_qp] - _vel_x_face[_qp], _vel_y[_qp] - _vel_y_face[_qp], 0);
     112             : 
     113     7348224 :   return _gamma * soln_jump * _normals[_qp] * _normals[_qp](_comp);
     114             : }