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