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 58221 : AdvectionIPHDGAssemblyHelper::validParams() 19 : { 20 58221 : auto params = IPHDGAssemblyHelper::validParams(); 21 232884 : params.addRequiredParam<MaterialPropertyName>("velocity", "Velocity vector"); 22 174663 : params.addParam<Real>( 23 116442 : "coeff", 1, "A constant coefficient. This could be something like a density"); 24 116442 : params.addParam<bool>("self_advection", 25 116442 : 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 58221 : 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 : }