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 "DiffusionIPHDGAssemblyHelper.h"
      11             : #include "MooseTypes.h"
      12             : #include "MooseVariableDependencyInterface.h"
      13             : #include "MooseVariableScalar.h"
      14             : #include "SystemBase.h"
      15             : #include "MooseMesh.h"
      16             : #include "MooseObject.h"
      17             : #include "MaterialPropertyInterface.h"
      18             : 
      19             : using namespace libMesh;
      20             : 
      21             : InputParameters
      22       44539 : DiffusionIPHDGAssemblyHelper::validParams()
      23             : {
      24       44539 :   auto params = IPHDGAssemblyHelper::validParams();
      25      178156 :   params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
      26       89078 :   params.addParam<Real>("alpha",
      27       89078 :                         1,
      28             :                         "The stabilization coefficient required for discontinuous Galerkin "
      29             :                         "schemes.");
      30       44539 :   return params;
      31           0 : }
      32             : 
      33         878 : DiffusionIPHDGAssemblyHelper::DiffusionIPHDGAssemblyHelper(
      34             :     const MooseObject * const moose_obj,
      35             :     MooseVariableDependencyInterface * const mvdi,
      36             :     const TransientInterface * const ti,
      37             :     SystemBase & sys,
      38             :     const Assembly & assembly,
      39             :     const THREAD_ID tid,
      40             :     const std::set<SubdomainID> & block_ids,
      41         878 :     const std::set<BoundaryID> & boundary_ids)
      42             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      43         878 :     _diff(this->getADMaterialProperty<Real>("diffusivity")),
      44        1756 :     _face_diff(this->getFaceADMaterialProperty<Real>("diffusivity")),
      45        2634 :     _alpha(moose_obj->getParam<Real>("alpha"))
      46             : {
      47         878 : }
      48             : 
      49             : void
      50      574783 : DiffusionIPHDGAssemblyHelper::scalarVolume()
      51             : {
      52     2905415 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      53    10239808 :     for (const auto i : index_range(_scalar_re))
      54     7909176 :       _scalar_re(i) += _ip_JxW[qp] * (_grad_scalar_phi[i][qp] * _diff[qp] * _grad_u_sol[qp]);
      55      574783 : }
      56             : 
      57             : void
      58     2065775 : DiffusionIPHDGAssemblyHelper::scalarFace()
      59             : {
      60     2065775 :   const auto h_elem = _elem_volume / _side_area;
      61             : 
      62     9159212 :   for (const auto i : index_range(_scalar_re))
      63    21352131 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      64             :     {
      65    28517388 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
      66    42776082 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
      67    42776082 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      68    57034776 :                        (_lm_u_sol[qp] - _u_sol[qp]) * _scalar_phi_face[i][qp];
      69    42776082 :       _scalar_re(i) += _ip_JxW_face[qp] * (_lm_u_sol[qp] - _u_sol[qp]) * _face_diff[qp] *
      70    57034776 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
      71             :     }
      72     2065775 : }
      73             : 
      74             : void
      75     2065775 : DiffusionIPHDGAssemblyHelper::lmFace()
      76             : {
      77     2065775 :   const auto h_elem = _elem_volume / _side_area;
      78             : 
      79    17578809 :   for (const auto i : index_range(_lm_re))
      80    46730622 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      81             :     {
      82    62435176 :       _lm_re(i) += _ip_JxW_face[qp] * _face_diff[qp] * _grad_u_sol[qp] * _ip_normals[qp] *
      83    62435176 :                    _lm_phi_face[i][qp];
      84    93652764 :       _lm_re(i) += _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      85   124870352 :                    (_lm_u_sol[qp] - _u_sol[qp]) * _lm_phi_face[i][qp];
      86             :     }
      87     2065775 : }
      88             : 
      89             : void
      90       59162 : DiffusionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_value)
      91             : {
      92       59162 :   const auto h_elem = _elem_volume / _side_area;
      93             : 
      94      178746 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      95             :   {
      96      358752 :     const auto scalar_value = dirichlet_value(
      97      119584 :         Moose::ElemSideQpArg{
      98      119584 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
      99      119584 :         _ti.determineState());
     100             : 
     101      520672 :     for (const auto i : index_range(_u_dof_indices))
     102             :     {
     103      802176 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
     104     1203264 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
     105     1203264 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
     106     1604352 :                        (scalar_value - _u_sol[qp]) * _scalar_phi_face[i][qp];
     107     1203264 :       _scalar_re(i) += _ip_JxW_face[qp] * (scalar_value - _u_sol[qp]) * _face_diff[qp] *
     108     1604352 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
     109             :     }
     110             :   }
     111       59162 : }