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       60141 : AdvectionIPHDGAssemblyHelper::validParams()
      19             : {
      20       60141 :   auto params = IPHDGAssemblyHelper::validParams();
      21      240564 :   params.addRequiredParam<MaterialPropertyName>("velocity", "Velocity vector");
      22      180423 :   params.addParam<Real>(
      23      120282 :       "coeff", 1, "A constant coefficient. This could be something like a density");
      24      120282 :   params.addParam<bool>("self_advection",
      25      120282 :                         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       60141 :   return params;
      30           0 : }
      31             : 
      32         582 : 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         582 :     const std::set<BoundaryID> & boundary_ids)
      41             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      42         582 :     _velocity(getADMaterialProperty<RealVectorValue>("velocity")),
      43        1164 :     _face_velocity(getFaceADMaterialProperty<RealVectorValue>("velocity")),
      44        1164 :     _coeff(moose_obj->getParam<Real>("coeff")),
      45        1746 :     _self_advection(moose_obj->getParam<bool>("self_advection"))
      46             : {
      47         582 : }
      48             : 
      49             : void
      50      518553 : AdvectionIPHDGAssemblyHelper::scalarVolume()
      51             : {
      52     2591901 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      53             :   {
      54     2073348 :     ADReal adv_quant = _coeff;
      55     2073348 :     if (_self_advection)
      56     2073348 :       adv_quant *= _u_sol[qp];
      57     2073348 :     const auto qp_term = _ip_JxW[qp] * _velocity[qp] * adv_quant;
      58     9209916 :     for (const auto i : index_range(_scalar_re))
      59     7136568 :       _scalar_re(i) -= _grad_scalar_phi[i][qp] * qp_term;
      60     2073348 :   }
      61      518553 : }
      62             : 
      63             : void
      64     1851129 : AdvectionIPHDGAssemblyHelper::scalarFace()
      65             : {
      66     5552694 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      67             :   {
      68     3701565 :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      69     3701565 :     ADReal adv_quant = _coeff;
      70     3701565 :     if (_self_advection)
      71     3701565 :       adv_quant *= (MetaPhysicL::raw_value(vdotn) >= 0 ? _u_sol[qp] : _lm_u_sol[qp]);
      72     3701565 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
      73    16597611 :     for (const auto i : index_range(_scalar_re))
      74    12896046 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
      75     3701565 :   }
      76     1851129 : }
      77             : 
      78             : void
      79     1851003 : AdvectionIPHDGAssemblyHelper::lmFace()
      80             : {
      81     5552379 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      82             :   {
      83     3701376 :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      84     3701376 :     ADReal adv_quant = _coeff;
      85     3701376 :     if (_self_advection)
      86     3701376 :       adv_quant *= (MetaPhysicL::raw_value(vdotn) >= 0 ? _u_sol[qp] : _lm_u_sol[qp]);
      87     3701376 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
      88    31283784 :     for (const auto i : index_range(_lm_re))
      89    27582408 :       _lm_re(i) -= _lm_phi_face[i][qp] * qp_term;
      90     3701376 :   }
      91     1851003 : }
      92             : 
      93             : void
      94       47439 : AdvectionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_functor)
      95             : {
      96      142254 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      97             :   {
      98       94815 :     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      284445 :     const auto dirichlet_value = dirichlet_functor(
     101       94815 :         Moose::ElemSideQpArg{
     102       94815 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
     103       94815 :         _ti.determineState());
     104       94815 :     const auto adv_quant = dirichlet_value * _coeff;
     105       94815 :     const auto qp_term = _ip_JxW_face[qp] * vdotn * adv_quant;
     106      421497 :     for (const auto i : index_range(_scalar_re))
     107      326682 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
     108       94815 :   }
     109       47439 : }
     110             : 
     111             : void
     112         126 : AdvectionIPHDGAssemblyHelper::lmOutflow()
     113             : {
     114         315 :   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         189 :     const auto qp_term = _ip_JxW_face[qp] * _coeff * (_lm_u_sol[qp] - _u_sol[qp]);
     122        1071 :     for (const auto i : index_range(_lm_re))
     123             :       // Force the LM solution to be equivalent to the internal solution
     124         882 :       _lm_re(i) += _lm_phi_face[i][qp] * qp_term;
     125         189 :   }
     126         126 : }