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       10295 : DiffusionIPHDGAssemblyHelper::validParams()
      23             : {
      24       10295 :   auto params = IPHDGAssemblyHelper::validParams();
      25       41180 :   params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
      26       30885 :   params.addRequiredParam<Real>("alpha",
      27             :                                 "The stabilization coefficient required for discontinuous Galerkin "
      28             :                                 "schemes.");
      29       10295 :   return params;
      30           0 : }
      31             : 
      32         562 : DiffusionIPHDGAssemblyHelper::DiffusionIPHDGAssemblyHelper(
      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         562 :     const std::set<BoundaryID> & boundary_ids)
      41             :   : IPHDGAssemblyHelper(moose_obj, mvdi, ti, sys, assembly, tid, block_ids, boundary_ids),
      42         562 :     _diff(this->getADMaterialProperty<Real>("diffusivity")),
      43        1124 :     _face_diff(this->getFaceADMaterialProperty<Real>("diffusivity")),
      44        1686 :     _alpha(moose_obj->getParam<Real>("alpha"))
      45             : {
      46         562 : }
      47             : 
      48             : void
      49      100501 : DiffusionIPHDGAssemblyHelper::scalarVolume()
      50             : {
      51      525005 :   for (const auto qp : make_range(_ip_qrule->n_points()))
      52     1796296 :     for (const auto i : index_range(_scalar_re))
      53     1371792 :       _scalar_re(i) += _ip_JxW[qp] * (_grad_scalar_phi[i][qp] * _diff[qp] * _grad_u_sol[qp]);
      54      100501 : }
      55             : 
      56             : void
      57      363318 : DiffusionIPHDGAssemblyHelper::scalarFace()
      58             : {
      59      363318 :   const auto h_elem = _ip_current_elem->hmax();
      60             : 
      61     1547016 :   for (const auto i : index_range(_scalar_re))
      62     3602394 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      63             :     {
      64     4837392 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
      65     7256088 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
      66     7256088 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      67     9674784 :                        (_lm_u_sol[qp] - _u_sol[qp]) * _scalar_phi_face[i][qp];
      68     7256088 :       _scalar_re(i) += _ip_JxW_face[qp] * (_lm_u_sol[qp] - _u_sol[qp]) * _face_diff[qp] *
      69     9674784 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
      70             :     }
      71      363318 : }
      72             : 
      73             : void
      74      363318 : DiffusionIPHDGAssemblyHelper::lmFace()
      75             : {
      76      363318 :   const auto h_elem = _ip_current_elem->hmax();
      77             : 
      78     2442048 :   for (const auto i : index_range(_lm_re))
      79     6372990 :     for (const auto qp : make_range(_ip_qrule_face->n_points()))
      80             :     {
      81     8588520 :       _lm_re(i) += _ip_JxW_face[qp] * _face_diff[qp] * _grad_u_sol[qp] * _ip_normals[qp] *
      82     8588520 :                    _lm_phi_face[i][qp];
      83    12882780 :       _lm_re(i) += _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
      84    17177040 :                    (_lm_u_sol[qp] - _u_sol[qp]) * _lm_phi_face[i][qp];
      85             :     }
      86      363318 : }
      87             : 
      88             : void
      89       19236 : DiffusionIPHDGAssemblyHelper::scalarDirichlet(const Moose::Functor<Real> & dirichlet_value)
      90             : {
      91       19236 :   const auto h_elem = _ip_current_elem->hmax();
      92             : 
      93       58608 :   for (const auto qp : make_range(_ip_qrule_face->n_points()))
      94             :   {
      95      118116 :     const auto scalar_value = dirichlet_value(
      96       39372 :         Moose::ElemSideQpArg{
      97       39372 :             _ip_current_elem, _ip_current_side, qp, _ip_qrule_face, _ip_q_point_face[qp]},
      98       39372 :         _ti.determineState());
      99             : 
     100      166560 :     for (const auto i : index_range(_u_dof_indices))
     101             :     {
     102      254376 :       _scalar_re(i) -= _ip_JxW_face[qp] * _face_diff[qp] * _scalar_phi_face[i][qp] *
     103      381564 :                        (_grad_u_sol[qp] * _ip_normals[qp]);
     104      381564 :       _scalar_re(i) -= _ip_JxW_face[qp] * _alpha / h_elem * _face_diff[qp] *
     105      508752 :                        (scalar_value - _u_sol[qp]) * _scalar_phi_face[i][qp];
     106      381564 :       _scalar_re(i) += _ip_JxW_face[qp] * (scalar_value - _u_sol[qp]) * _face_diff[qp] *
     107      508752 :                        _grad_scalar_phi_face[i][qp] * _ip_normals[qp];
     108             :     }
     109             :   }
     110       19236 : }