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       12847 : AdvectionIPHDGAssemblyHelper::validParams()
      19             : {
      20       12847 :   auto params = IPHDGAssemblyHelper::validParams();
      21       51388 :   params.addRequiredParam<MaterialPropertyName>("velocity", "Velocity vector");
      22       38541 :   params.addRequiredParam<Real>("coeff",
      23             :                                 "A constant coefficient. This could be something like a density");
      24       12847 :   return params;
      25           0 : }
      26             : 
      27         303 : AdvectionIPHDGAssemblyHelper::AdvectionIPHDGAssemblyHelper(
      28             :     const MooseObject * const moose_obj,
      29             :     MooseVariableDependencyInterface * const mvdi,
      30             :     const TransientInterface * const ti,
      31             :     SystemBase & sys,
      32             :     const Assembly & assembly,
      33             :     const THREAD_ID tid,
      34             :     const std::set<SubdomainID> & block_ids,
      35         303 :     const std::set<BoundaryID> & boundary_ids)
      36             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      37         303 :     _velocity(getADMaterialProperty<RealVectorValue>("velocity")),
      38         606 :     _face_velocity(getFaceADMaterialProperty<RealVectorValue>("velocity")),
      39         909 :     _coeff(moose_obj->getParam<Real>("coeff"))
      40             : {
      41         303 : }
      42             : 
      43             : void
      44       56629 : AdvectionIPHDGAssemblyHelper::scalarVolume()
      45             : {
      46      282413 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      47             :   {
      48      225784 :     const auto adv_quant = _coeff * _u_sol[qp];
      49      225784 :     const auto qp_term = _ip_JxW[qp] * _velocity[qp] * adv_quant;
      50     1000768 :     for (const auto i : index_range(_scalar_re))
      51      774984 :       _scalar_re(i) -= _grad_scalar_phi[i][qp] * qp_term;
      52      225784 :   }
      53       56629 : }
      54             : 
      55             : ADReal
      56      802680 : AdvectionIPHDGAssemblyHelper::computeFlux(const unsigned int qp, const ADReal & face_value)
      57             : {
      58      802680 :   const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
      59      802680 :   const auto face_phi = _coeff * face_value;
      60      802680 :   const auto internal_phi = _coeff * _u_sol[qp];
      61             :   // Short form for writing upwinding. If the velocity is in the direction of the surface normal,
      62             :   // the the internal value is used, else the face value is used
      63     1605360 :   return 0.5 * vdotn * (internal_phi + face_phi) + 0.5 * abs(vdotn) * (internal_phi - face_phi);
      64      802680 : }
      65             : 
      66             : void
      67      195906 : AdvectionIPHDGAssemblyHelper::scalarFace()
      68             : {
      69      587122 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      70             :   {
      71      391216 :     const auto qp_term = _ip_JxW_face[qp] * computeFlux(qp, _lm_u_sol[qp]);
      72     1751616 :     for (const auto i : index_range(_scalar_re))
      73     1360400 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
      74      391216 :   }
      75      195906 : }
      76             : 
      77             : void
      78      195826 : AdvectionIPHDGAssemblyHelper::lmFace()
      79             : {
      80      586934 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      81             :   {
      82      391108 :     const auto qp_term = _ip_JxW_face[qp] * computeFlux(qp, _lm_u_sol[qp]);
      83     1861752 :     for (const auto i : index_range(_lm_re))
      84     1470644 :       _lm_re(i) -= _lm_phi_face[i][qp] * qp_term;
      85      391108 :   }
      86      195826 : }
      87             : 
      88             : void
      89       10204 : AdvectionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_functor)
      90             : {
      91       30560 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      92             :   {
      93       61068 :     const auto dirichlet_value = dirichlet_functor(
      94       20356 :         Moose::ElemSideQpArg{
      95       20356 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
      96       20356 :         _ti.determineState());
      97       20356 :     const auto qp_term = _ip_JxW_face[qp] * computeFlux(qp, dirichlet_value);
      98       90416 :     for (const auto i : index_range(_scalar_re))
      99       70060 :       _scalar_re(i) += _scalar_phi_face[i][qp] * qp_term;
     100       20356 :   }
     101       10204 : }
     102             : 
     103             : void
     104          52 : AdvectionIPHDGAssemblyHelper::lmOutflow()
     105             : {
     106         104 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
     107             :   {
     108             : #ifndef NDEBUG
     109             :     const auto vdotn = _face_velocity[qp] * _ip_normals[qp];
     110             :     mooseAssert(MetaPhysicL::raw_value(vdotn) >= 0, "The velocity must create outflow conditions");
     111             : #endif
     112          52 :     const auto qp_term = _ip_JxW_face[qp] * _coeff * (_lm_u_sol[qp] - _u_sol[qp]);
     113         156 :     for (const auto i : index_range(_lm_re))
     114             :       // Force the LM solution to be equivalent to the internal solution
     115         104 :       _lm_re(i) += _lm_phi_face[i][qp] * qp_term;
     116          52 :   }
     117          52 : }