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 "ArrayDGDiffusion.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariableFE.h"
      14             : 
      15             : #include "libmesh/utility.h"
      16             : 
      17             : registerMooseObject("MooseApp", ArrayDGDiffusion);
      18             : 
      19             : InputParameters
      20        3145 : ArrayDGDiffusion::validParams()
      21             : {
      22        3145 :   InputParameters params = ArrayDGKernel::validParams();
      23       12580 :   params.addRequiredParam<MaterialPropertyName>(
      24             :       "diff", "The diffusion (or thermal conductivity or viscosity) coefficient.");
      25       12580 :   params.addParam<Real>("sigma", 4, "sigma");
      26       12580 :   params.addParam<Real>("epsilon", 1, "epsilon");
      27        3145 :   params.addClassDescription("Implements interior penalty method for array diffusion equations.");
      28        3145 :   return params;
      29           0 : }
      30             : 
      31          44 : ArrayDGDiffusion::ArrayDGDiffusion(const InputParameters & parameters)
      32             :   : ArrayDGKernel(parameters),
      33          44 :     _epsilon(getParam<Real>("epsilon")),
      34          88 :     _sigma(getParam<Real>("sigma")),
      35          88 :     _diff(getMaterialProperty<RealEigenVector>("diff")),
      36          88 :     _diff_neighbor(getNeighborMaterialProperty<RealEigenVector>("diff")),
      37          44 :     _res1(_count),
      38          88 :     _res2(_count)
      39             : {
      40          44 : }
      41             : 
      42             : void
      43       11904 : ArrayDGDiffusion::initQpResidual(Moose::DGResidualType type)
      44             : {
      45             :   mooseAssert(_diff[_qp].size() == _count && _diff_neighbor[_qp].size() == _count,
      46             :               "'diff' size is inconsistent with the number of components of array "
      47             :               "variable");
      48             : 
      49       11904 :   const int elem_b_order = std::max(libMesh::Order(1), _var.order());
      50             :   const Real h_elem =
      51       11904 :       _current_elem_volume / _current_side_volume * 1.0 / Utility::pow<2>(elem_b_order);
      52             : 
      53             :   // WARNING: use noalias() syntax with caution. See ArrayDiffusion.C for more details.
      54       11904 :   _res1.noalias() = _diff[_qp].asDiagonal() * _grad_u[_qp] * _array_normals[_qp];
      55       11904 :   _res1.noalias() += _diff_neighbor[_qp].asDiagonal() * _grad_u_neighbor[_qp] * _array_normals[_qp];
      56       11904 :   _res1 *= 0.5;
      57       11904 :   _res1 -= (_u[_qp] - _u_neighbor[_qp]) * _sigma / h_elem;
      58             : 
      59       11904 :   switch (type)
      60             :   {
      61        5952 :     case Moose::Element:
      62        5952 :       _res2.noalias() = _diff[_qp].asDiagonal() * (_u[_qp] - _u_neighbor[_qp]) * _epsilon * 0.5;
      63        5952 :       break;
      64             : 
      65        5952 :     case Moose::Neighbor:
      66        5952 :       _res2.noalias() =
      67       11904 :           _diff_neighbor[_qp].asDiagonal() * (_u[_qp] - _u_neighbor[_qp]) * _epsilon * 0.5;
      68        5952 :       break;
      69             :   }
      70       11904 : }
      71             : 
      72             : void
      73       47616 : ArrayDGDiffusion::computeQpResidual(Moose::DGResidualType type, RealEigenVector & residual)
      74             : {
      75       47616 :   switch (type)
      76             :   {
      77       23808 :     case Moose::Element:
      78       23808 :       residual = -_res1 * _test[_i][_qp] + _res2 * (_grad_test[_i][_qp] * _normals[_qp]);
      79       23808 :       break;
      80             : 
      81       23808 :     case Moose::Neighbor:
      82             :       residual =
      83       23808 :           _res1 * _test_neighbor[_i][_qp] + _res2 * (_grad_test_neighbor[_i][_qp] * _normals[_qp]);
      84       23808 :       break;
      85             :   }
      86       47616 : }
      87             : 
      88             : RealEigenVector
      89      172032 : ArrayDGDiffusion::computeQpJacobian(Moose::DGJacobianType type)
      90             : {
      91      172032 :   RealEigenVector r = RealEigenVector::Zero(_count);
      92             : 
      93      172032 :   const int elem_b_order = std::max(libMesh::Order(1), _var.order());
      94             :   const Real h_elem =
      95      172032 :       _current_elem_volume / _current_side_volume * 1.0 / Utility::pow<2>(elem_b_order);
      96             : 
      97      172032 :   switch (type)
      98             :   {
      99       43008 :     case Moose::ElementElement:
     100       43008 :       r -= _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp] * 0.5 * _diff[_qp];
     101       43008 :       r += _grad_test[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi[_j][_qp] * _diff[_qp];
     102       43008 :       r += RealEigenVector::Constant(_count, _sigma / h_elem * _phi[_j][_qp] * _test[_i][_qp]);
     103       43008 :       break;
     104             : 
     105       43008 :     case Moose::ElementNeighbor:
     106       43008 :       r -= _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test[_i][_qp] * 0.5 * _diff_neighbor[_qp];
     107       43008 :       r -= _grad_test[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi_neighbor[_j][_qp] *
     108       86016 :            _diff[_qp];
     109           0 :       r -= RealEigenVector::Constant(_count,
     110       43008 :                                      _sigma / h_elem * _phi_neighbor[_j][_qp] * _test[_i][_qp]);
     111       43008 :       break;
     112             : 
     113       43008 :     case Moose::NeighborElement:
     114       43008 :       r += _grad_phi[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp] * 0.5 * _diff[_qp];
     115       43008 :       r += _grad_test_neighbor[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi[_j][_qp] *
     116       86016 :            _diff_neighbor[_qp];
     117           0 :       r -= RealEigenVector::Constant(_count,
     118       43008 :                                      _sigma / h_elem * _phi[_j][_qp] * _test_neighbor[_i][_qp]);
     119       43008 :       break;
     120             : 
     121       43008 :     case Moose::NeighborNeighbor:
     122       43008 :       r += _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp] * 0.5 *
     123       86016 :            _diff_neighbor[_qp];
     124       43008 :       r -= _grad_test_neighbor[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi_neighbor[_j][_qp] *
     125       86016 :            _diff_neighbor[_qp];
     126           0 :       r += RealEigenVector::Constant(
     127       43008 :           _count, _sigma / h_elem * _phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp]);
     128       43008 :       break;
     129             :   }
     130             : 
     131      344064 :   return r;
     132           0 : }