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       14355 : ArrayDGDiffusion::validParams()
      21             : {
      22       14355 :   InputParameters params = ArrayDGKernel::validParams();
      23       14355 :   params.addRequiredParam<MaterialPropertyName>(
      24             :       "diff", "The diffusion (or thermal conductivity or viscosity) coefficient.");
      25       14355 :   params.addParam<Real>("sigma", 4, "sigma");
      26       14355 :   params.addParam<Real>("epsilon", 1, "epsilon");
      27       14355 :   params.addClassDescription("Implements interior penalty method for array diffusion equations.");
      28       14355 :   return params;
      29           0 : }
      30             : 
      31          47 : ArrayDGDiffusion::ArrayDGDiffusion(const InputParameters & parameters)
      32             :   : ArrayDGKernel(parameters),
      33          47 :     _epsilon(getParam<Real>("epsilon")),
      34          47 :     _sigma(getParam<Real>("sigma")),
      35          47 :     _diff(getMaterialProperty<RealEigenVector>("diff")),
      36          47 :     _diff_neighbor(getNeighborMaterialProperty<RealEigenVector>("diff")),
      37          47 :     _res1(_count),
      38          94 :     _res2(_count)
      39             : {
      40          47 : }
      41             : 
      42             : void
      43       15648 : 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       15648 :   const int elem_b_order = std::max(libMesh::Order(1), _var.order());
      50             :   const Real h_elem =
      51       15648 :       _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       15648 :   _res1.noalias() = _diff[_qp].asDiagonal() * _grad_u[_qp] * _array_normals[_qp];
      55       15648 :   _res1.noalias() += _diff_neighbor[_qp].asDiagonal() * _grad_u_neighbor[_qp] * _array_normals[_qp];
      56       15648 :   _res1 *= 0.5;
      57       15648 :   _res1 -= (_u[_qp] - _u_neighbor[_qp]) * _sigma / h_elem;
      58             : 
      59       15648 :   switch (type)
      60             :   {
      61        7824 :     case Moose::Element:
      62        7824 :       _res2.noalias() = _diff[_qp].asDiagonal() * (_u[_qp] - _u_neighbor[_qp]) * _epsilon * 0.5;
      63        7824 :       break;
      64             : 
      65        7824 :     case Moose::Neighbor:
      66        7824 :       _res2.noalias() =
      67       15648 :           _diff_neighbor[_qp].asDiagonal() * (_u[_qp] - _u_neighbor[_qp]) * _epsilon * 0.5;
      68        7824 :       break;
      69             :   }
      70       15648 : }
      71             : 
      72             : void
      73       62592 : ArrayDGDiffusion::computeQpResidual(Moose::DGResidualType type, RealEigenVector & residual)
      74             : {
      75       62592 :   switch (type)
      76             :   {
      77       31296 :     case Moose::Element:
      78       31296 :       residual = -_res1 * _test[_i][_qp] + _res2 * (_grad_test[_i][_qp] * _normals[_qp]);
      79       31296 :       break;
      80             : 
      81       31296 :     case Moose::Neighbor:
      82             :       residual =
      83       31296 :           _res1 * _test_neighbor[_i][_qp] + _res2 * (_grad_test_neighbor[_i][_qp] * _normals[_qp]);
      84       31296 :       break;
      85             :   }
      86       62592 : }
      87             : 
      88             : RealEigenVector
      89      190464 : ArrayDGDiffusion::computeQpJacobian(Moose::DGJacobianType type)
      90             : {
      91      190464 :   RealEigenVector r = RealEigenVector::Zero(_count);
      92             : 
      93      190464 :   const int elem_b_order = std::max(libMesh::Order(1), _var.order());
      94             :   const Real h_elem =
      95      190464 :       _current_elem_volume / _current_side_volume * 1.0 / Utility::pow<2>(elem_b_order);
      96             : 
      97      190464 :   switch (type)
      98             :   {
      99       47616 :     case Moose::ElementElement:
     100       47616 :       r -= _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp] * 0.5 * _diff[_qp];
     101       47616 :       r += _grad_test[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi[_j][_qp] * _diff[_qp];
     102       47616 :       r += RealEigenVector::Constant(_count, _sigma / h_elem * _phi[_j][_qp] * _test[_i][_qp]);
     103       47616 :       break;
     104             : 
     105       47616 :     case Moose::ElementNeighbor:
     106       47616 :       r -= _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test[_i][_qp] * 0.5 * _diff_neighbor[_qp];
     107       47616 :       r -= _grad_test[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi_neighbor[_j][_qp] *
     108       95232 :            _diff[_qp];
     109           0 :       r -= RealEigenVector::Constant(_count,
     110       47616 :                                      _sigma / h_elem * _phi_neighbor[_j][_qp] * _test[_i][_qp]);
     111       47616 :       break;
     112             : 
     113       47616 :     case Moose::NeighborElement:
     114       47616 :       r += _grad_phi[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp] * 0.5 * _diff[_qp];
     115       47616 :       r += _grad_test_neighbor[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi[_j][_qp] *
     116       95232 :            _diff_neighbor[_qp];
     117           0 :       r -= RealEigenVector::Constant(_count,
     118       47616 :                                      _sigma / h_elem * _phi[_j][_qp] * _test_neighbor[_i][_qp]);
     119       47616 :       break;
     120             : 
     121       47616 :     case Moose::NeighborNeighbor:
     122       47616 :       r += _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp] * 0.5 *
     123       95232 :            _diff_neighbor[_qp];
     124       47616 :       r -= _grad_test_neighbor[_i][_qp] * _normals[_qp] * _epsilon * 0.5 * _phi_neighbor[_j][_qp] *
     125       95232 :            _diff_neighbor[_qp];
     126           0 :       r += RealEigenVector::Constant(
     127       47616 :           _count, _sigma / h_elem * _phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp]);
     128       47616 :       break;
     129             :   }
     130             : 
     131      380928 :   return r;
     132           0 : }