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       44387 : DiffusionIPHDGAssemblyHelper::validParams()
      23             : {
      24       44387 :   auto params = IPHDGAssemblyHelper::validParams();
      25       44387 :   params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
      26      133161 :   params.addParam<Real>("alpha",
      27       88774 :                         1,
      28             :                         "The stabilization coefficient required for discontinuous Galerkin "
      29             :                         "schemes.");
      30       44387 :   return params;
      31           0 : }
      32             : 
      33         802 : 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         802 :     const std::set<BoundaryID> & boundary_ids)
      42             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      43         802 :     _diff(this->getADMaterialProperty<Real>("diffusivity")),
      44         802 :     _face_diff(this->getFaceADMaterialProperty<Real>("diffusivity")),
      45        1604 :     _alpha(moose_obj->getParam<Real>("alpha"))
      46             : {
      47         802 : }
      48             : 
      49             : void
      50      566240 : DiffusionIPHDGAssemblyHelper::scalarVolume()
      51             : {
      52     2862700 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      53    10103120 :     for (const auto i : index_range(_scalar_re))
      54     7806660 :       _scalar_re(i) += _ip_JxW[qp] * (_grad_scalar_phi[i][qp] * _diff[qp] * _grad_u_sol[qp]);
      55      566240 : }
      56             : 
      57             : void
      58     2035564 : DiffusionIPHDGAssemblyHelper::scalarFace()
      59             : {
      60     2035564 :   const auto h_elem = _elem_volume / _side_area;
      61             : 
      62     9038368 :   for (const auto i : index_range(_scalar_re))
      63    21080232 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      64             :     {
      65    28154856 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
      66    42232284 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
      67    42232284 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      68    56309712 :                        (_lm_u_sol[qp] - _u_sol[qp]) * _scalar_phi_face[i][qp];
      69    42232284 :       _scalar_re(i) += _ip_JxW_face[qp] * (_lm_u_sol[qp] - _u_sol[qp]) * _face_diff[qp] *
      70    56309712 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
      71             :     }
      72     2035564 : }
      73             : 
      74             : void
      75     2035564 : DiffusionIPHDGAssemblyHelper::lmFace()
      76             : {
      77     2035564 :   const auto h_elem = _elem_volume / _side_area;
      78             : 
      79    17320140 :   for (const auto i : index_range(_lm_re))
      80    46045248 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      81             :     {
      82    61521344 :       _lm_re(i) += _ip_JxW_face[qp] * _face_diff[qp] * _grad_u_sol[qp] * _ip_normals[qp] *
      83    61521344 :                    _lm_phi_face[i][qp];
      84    92282016 :       _lm_re(i) += _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      85   123042688 :                    (_lm_u_sol[qp] - _u_sol[qp]) * _lm_phi_face[i][qp];
      86             :     }
      87     2035564 : }
      88             : 
      89             : void
      90       57532 : DiffusionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_value)
      91             : {
      92       57532 :   const auto h_elem = _elem_volume / _side_area;
      93             : 
      94      173856 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      95             :   {
      96      348972 :     const auto scalar_value = dirichlet_value(
      97      116324 :         Moose::ElemSideQpArg{
      98      116324 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
      99      116324 :         _ti.determineState());
     100             : 
     101      507632 :     for (const auto i : index_range(_u_dof_indices))
     102             :     {
     103      782616 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
     104     1173924 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
     105     1173924 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
     106     1565232 :                        (scalar_value - _u_sol[qp]) * _scalar_phi_face[i][qp];
     107     1173924 :       _scalar_re(i) += _ip_JxW_face[qp] * (scalar_value - _u_sol[qp]) * _face_diff[qp] *
     108     1565232 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
     109             :     }
     110             :   }
     111       57532 : }