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 "AdvectionIPHDGAssemblyHelper.h"
      11             : #include "MooseTypes.h"
      12             : #include "MooseObject.h"
      13             : #include "MaterialPropertyInterface.h"
      14             : 
      15             : using namespace libMesh;
      16             : 
      17             : InputParameters
      18       43800 : AdvectionIPHDGAssemblyHelper::validParams()
      19             : {
      20       43800 :   auto params = IPHDGAssemblyHelper::validParams();
      21       43800 :   params.addRequiredParam<MaterialPropertyName>("velocity", "Velocity vector");
      22      131400 :   params.addParam<Real>(
      23       87600 :       "coeff", 1, "A constant coefficient. This could be something like a density");
      24      131400 :   params.addParam<bool>("self_advection",
      25       87600 :                         true,
      26             :                         "Whether this kernel should advect its variables, e.g. its "
      27             :                         "variable/side_variable pair. If false, we will advect "
      28             :                         "unity (possibly multiplied by the 'coeff' parameter");
      29       43800 :   return params;
      30           0 : }
      31             : 
      32         504 : AdvectionIPHDGAssemblyHelper::AdvectionIPHDGAssemblyHelper(
      33             :     const MooseObject * const moose_obj,
      34             :     MooseVariableDependencyInterface * const mvdi,
      35             :     const TransientInterface * const ti,
      36             :     SystemBase & sys,
      37             :     const Assembly & assembly,
      38             :     const THREAD_ID tid,
      39             :     const std::set<SubdomainID> & block_ids,
      40         504 :     const std::set<BoundaryID> & boundary_ids)
      41             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      42         504 :     _velocity(getADMaterialProperty<RealVectorValue>("velocity")),
      43         504 :     _face_velocity(getFaceADMaterialProperty<RealVectorValue>("velocity")),
      44         504 :     _coeff(moose_obj->getParam<Real>("coeff")),
      45        1008 :     _self_advection(moose_obj->getParam<bool>("self_advection"))
      46             : {
      47         504 : }
      48             : 
      49             : void
      50      516114 : AdvectionIPHDGAssemblyHelper::scalarVolume()
      51             : {
      52     2580030 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      53             :   {
      54     2063916 :     ADReal adv_quant = _coeff;
      55     2063916 :     if (_self_advection)
      56     2063916 :       adv_quant *= _u_sol[qp];
      57     2063916 :     const auto qp_term = _ip_JxW[qp] * _velocity[qp] * adv_quant;
      58     9172404 :     for (const auto i : index_range(_scalar_re))
      59     7108488 :       _scalar_re(i) -= _grad_scalar_phi[i][qp] * qp_term;
      60     2063916 :   }
      61      516114 : }
      62             : 
      63             : void
      64     1844649 : AdvectionIPHDGAssemblyHelper::scalarFace()
      65             : {
      66     5533434 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      67             :   {
      68     3688785 :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      69     3688785 :     ADReal adv_quant = _coeff;
      70     3688785 :     if (_self_advection)
      71     3688785 :       adv_quant *= (MetaPhysicL::raw_value(vdotn) >= 0 ? _u_sol[qp] : _lm_u_sol[qp]);
      72     3688785 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
      73    16546851 :     for (const auto i : index_range(_scalar_re))
      74    12858066 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
      75     3688785 :   }
      76     1844649 : }
      77             : 
      78             : void
      79     1844622 : AdvectionIPHDGAssemblyHelper::lmFace()
      80             : {
      81     5533380 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      82             :   {
      83     3688758 :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      84     3688758 :     ADReal adv_quant = _coeff;
      85     3688758 :     if (_self_advection)
      86     3688758 :       adv_quant *= (MetaPhysicL::raw_value(vdotn) >= 0 ? _u_sol[qp] : _lm_u_sol[qp]);
      87     3688758 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
      88    31196034 :     for (const auto i : index_range(_lm_re))
      89    27507276 :       _lm_re(i) -= _lm_phi_face[i][qp] * qp_term;
      90     3688758 :   }
      91     1844622 : }
      92             : 
      93             : void
      94       47403 : AdvectionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_functor)
      95             : {
      96      142182 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      97             :   {
      98       94779 :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      99             :     mooseAssert(_self_advection, "This shouldn't be called if we are not self-advecting");
     100      284337 :     const auto dirichlet_value = dirichlet_functor(
     101       94779 :         Moose::ElemSideQpArg{
     102       94779 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
     103       94779 :         _ti.determineState());
     104       94779 :     const auto adv_quant = dirichlet_value * _coeff;
     105       94779 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
     106      421425 :     for (const auto i : index_range(_scalar_re))
     107      326646 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
     108       94779 :   }
     109       47403 : }
     110             : 
     111             : void
     112          27 : AdvectionIPHDGAssemblyHelper::lmOutflow()
     113             : {
     114          54 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
     115             :   {
     116             : #ifndef NDEBUG
     117             :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
     118             :     mooseAssert(MetaPhysicL::raw_value(vdotn) >= 0, "The velocity must create outflow conditions");
     119             :     mooseAssert(_self_advection, "This shouldn't be called if we are not self-advecting");
     120             : #endif
     121          27 :     const auto qp_term = _ip_JxW_face[qp] * _coeff * (_lm_u_sol[qp] - _u_sol[qp]);
     122          81 :     for (const auto i : index_range(_lm_re))
     123             :       // Force the LM solution to be equivalent to the internal solution
     124          54 :       _lm_re(i) += _lm_phi_face[i][qp] * qp_term;
     125          27 :   }
     126          27 : }