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 "ArrayDGKernel.h"
      11             : #include "Assembly.h"
      12             : #include "MooseVariable.h"
      13             : #include "Problem.h"
      14             : #include "SubProblem.h"
      15             : #include "SystemBase.h"
      16             : #include "MaterialData.h"
      17             : #include "ParallelUniqueId.h"
      18             : 
      19             : #include "libmesh/dof_map.h"
      20             : #include "libmesh/dense_vector.h"
      21             : #include "libmesh/numeric_vector.h"
      22             : #include "libmesh/dense_subvector.h"
      23             : #include "libmesh/libmesh_common.h"
      24             : #include "libmesh/quadrature.h"
      25             : 
      26             : InputParameters
      27       43142 : ArrayDGKernel::validParams()
      28             : {
      29       43142 :   InputParameters params = DGKernelBase::validParams();
      30       43142 :   return params;
      31             : }
      32             : 
      33         182 : ArrayDGKernel::ArrayDGKernel(const InputParameters & parameters)
      34             :   : DGKernelBase(parameters),
      35             :     NeighborMooseVariableInterface(
      36             :         this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_ARRAY),
      37         364 :     _var(*mooseVariable()),
      38         182 :     _u(_is_implicit ? _var.sln() : _var.slnOld()),
      39         182 :     _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
      40             : 
      41         182 :     _phi(_var.phiFace()),
      42         182 :     _grad_phi(_var.gradPhiFace()),
      43         182 :     _array_grad_phi(_var.arrayGradPhiFace()),
      44             : 
      45         182 :     _test(_var.phiFace()),
      46         182 :     _grad_test(_var.gradPhiFace()),
      47         182 :     _array_grad_test(_var.arrayGradPhiFace()),
      48             : 
      49         182 :     _phi_neighbor(_var.phiFaceNeighbor()),
      50         182 :     _grad_phi_neighbor(_var.gradPhiFaceNeighbor()),
      51         182 :     _array_grad_phi_neighbor(_var.arrayGradPhiFaceNeighbor()),
      52             : 
      53         182 :     _test_neighbor(_var.phiFaceNeighbor()),
      54         182 :     _grad_test_neighbor(_var.gradPhiFaceNeighbor()),
      55         182 :     _array_grad_test_neighbor(_var.arrayGradPhiFaceNeighbor()),
      56             : 
      57         182 :     _u_neighbor(_is_implicit ? _var.slnNeighbor() : _var.slnOldNeighbor()),
      58         182 :     _grad_u_neighbor(_is_implicit ? _var.gradSlnNeighbor() : _var.gradSlnOldNeighbor()),
      59             : 
      60         182 :     _array_normals(_assembly.mappedNormals()),
      61         182 :     _count(_var.count()),
      62             : 
      63         364 :     _work_vector(_count)
      64             : {
      65         182 :   addMooseVariableDependency(mooseVariable());
      66             : 
      67         182 :   _save_in.resize(_save_in_strings.size());
      68         182 :   _diag_save_in.resize(_diag_save_in_strings.size());
      69             : 
      70         195 :   for (unsigned int i = 0; i < _save_in_strings.size(); i++)
      71             :   {
      72          13 :     MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
      73          13 :                                                          _save_in_strings[i],
      74             :                                                          Moose::VarKindType::VAR_AUXILIARY,
      75             :                                                          Moose::VarFieldType::VAR_FIELD_ARRAY);
      76             : 
      77          13 :     if (_sys.hasVariable(_save_in_strings[i]))
      78           0 :       mooseError("Trying to use solution variable " + _save_in_strings[i] +
      79           0 :                  " as a save_in variable in " + name());
      80             : 
      81          13 :     if (var->feType() != _var.feType())
      82           0 :       paramError(
      83             :           "save_in",
      84             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
      85           0 :           moose::internal::incompatVarMsg(*var, _var));
      86             : 
      87          13 :     _save_in[i] = var;
      88          13 :     var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
      89          13 :     addMooseVariableDependency(var);
      90             :   }
      91             : 
      92         182 :   _has_save_in = _save_in.size() > 0;
      93             : 
      94         195 :   for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
      95             :   {
      96          13 :     MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
      97          13 :                                                          _diag_save_in_strings[i],
      98             :                                                          Moose::VarKindType::VAR_AUXILIARY,
      99             :                                                          Moose::VarFieldType::VAR_FIELD_ARRAY);
     100             : 
     101          13 :     if (_sys.hasVariable(_diag_save_in_strings[i]))
     102           0 :       mooseError("Trying to use solution variable " + _diag_save_in_strings[i] +
     103           0 :                  " as a diag_save_in variable in " + name());
     104             : 
     105          13 :     if (var->feType() != _var.feType())
     106           0 :       paramError(
     107             :           "diag_save_in",
     108             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
     109           0 :           moose::internal::incompatVarMsg(*var, _var));
     110             : 
     111          13 :     _diag_save_in[i] = var;
     112          13 :     var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
     113          13 :     addMooseVariableDependency(var);
     114             :   }
     115             : 
     116         182 :   _has_diag_save_in = _diag_save_in.size() > 0;
     117         182 : }
     118             : 
     119             : void
     120       18552 : ArrayDGKernel::computeElemNeighResidual(Moose::DGResidualType type)
     121             : {
     122             :   bool is_elem;
     123       18552 :   if (type == Moose::Element)
     124        9276 :     is_elem = true;
     125             :   else
     126        9276 :     is_elem = false;
     127             : 
     128       18552 :   const ArrayVariableTestValue & test_space = is_elem ? _test : _test_neighbor;
     129             : 
     130       18552 :   if (is_elem)
     131        9276 :     prepareVectorTag(_assembly, _var.number());
     132             :   else
     133        9276 :     prepareVectorTagNeighbor(_assembly, _var.number());
     134             : 
     135       77496 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     136             :   {
     137       58944 :     initQpResidual(type);
     138      564096 :     for (_i = 0; _i < test_space.size(); _i++)
     139             :     {
     140      505152 :       _work_vector.setZero();
     141      505152 :       computeQpResidual(type, _work_vector);
     142             :       mooseAssert(_work_vector.size() == _count,
     143             :                   "Size of local residual is not equal to the number of array variable compoments");
     144      505152 :       _work_vector *= _JxW[_qp] * _coord[_qp];
     145      505152 :       _assembly.saveLocalArrayResidual(_local_re, _i, test_space.size(), _work_vector);
     146             :     }
     147             :   }
     148             : 
     149       18552 :   accumulateTaggedLocalResidual();
     150             : 
     151       18552 :   if (_has_save_in)
     152             :   {
     153        1152 :     Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
     154        2304 :     for (const auto & var : _save_in)
     155             :     {
     156        1152 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     157        1152 :       if (!avar)
     158           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     159             : 
     160        1152 :       if (is_elem)
     161         576 :         avar->addSolution(_local_re);
     162             :       else
     163         576 :         avar->addSolutionNeighbor(_local_re);
     164             :     }
     165        1152 :   }
     166       18552 : }
     167             : 
     168             : void
     169         768 : ArrayDGKernel::computeElemNeighJacobian(Moose::DGJacobianType type)
     170             : {
     171         768 :   const ArrayVariableTestValue & test_space =
     172         768 :       (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
     173        1344 :   const ArrayVariableTestValue & loc_phi =
     174         576 :       (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
     175             : 
     176         768 :   if (type == Moose::ElementElement)
     177         192 :     prepareMatrixTag(_assembly, _var.number(), _var.number());
     178             :   else
     179         576 :     prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), type);
     180             : 
     181        3840 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     182             :   {
     183        3072 :     initQpJacobian(type);
     184       33792 :     for (_i = 0; _i < test_space.size(); _i++)
     185      337920 :       for (_j = 0; _j < loc_phi.size(); _j++)
     186             :       {
     187      307200 :         RealEigenVector v = _JxW[_qp] * _coord[_qp] * computeQpJacobian(type);
     188      307200 :         _assembly.saveDiagLocalArrayJacobian(
     189      614400 :             _local_ke, _i, test_space.size(), _j, loc_phi.size(), _var.number(), v);
     190      307200 :       }
     191             :   }
     192             : 
     193         768 :   accumulateTaggedLocalMatrix();
     194             : 
     195         768 :   if (_has_diag_save_in && (type == Moose::ElementElement || type == Moose::NeighborNeighbor))
     196             :   {
     197           0 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     198           0 :     Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
     199           0 :     for (const auto & var : _diag_save_in)
     200             :     {
     201           0 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     202           0 :       if (!avar)
     203           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     204             : 
     205           0 :       if (type == Moose::ElementElement)
     206           0 :         avar->addSolution(diag);
     207             :       else
     208           0 :         avar->addSolutionNeighbor(diag);
     209             :     }
     210           0 :   }
     211         768 : }
     212             : 
     213     2121600 : RealEigenVector ArrayDGKernel::computeQpJacobian(Moose::DGJacobianType)
     214             : {
     215     4243200 :   return RealEigenVector::Zero(_count);
     216             : }
     217             : 
     218             : void
     219        1488 : ArrayDGKernel::computeOffDiagJacobian(const unsigned int jvar_num)
     220             : {
     221        1488 :   if (!excludeBoundary())
     222             :   {
     223        1488 :     const auto & jvar = getVariable(jvar_num);
     224             : 
     225        1488 :     precalculateOffDiagJacobian(jvar_num);
     226             : 
     227             :     // Compute element-element Jacobian
     228        1488 :     computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
     229             : 
     230             :     // Compute element-neighbor Jacobian
     231        1488 :     computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
     232             : 
     233             :     // Compute neighbor-element Jacobian
     234        1488 :     computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
     235             : 
     236             :     // Compute neighbor-neighbor Jacobian
     237        1488 :     computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
     238             :   }
     239        1488 : }
     240             : 
     241             : void
     242       13632 : ArrayDGKernel::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
     243             :                                                const MooseVariableFEBase & jvar)
     244             : {
     245       13632 :   const ArrayVariableTestValue & test_space =
     246       13632 :       (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
     247             : 
     248       13632 :   if (type == Moose::ElementElement)
     249        3408 :     prepareMatrixTag(_assembly, _var.number(), jvar.number());
     250             :   else
     251       10224 :     prepareMatrixTagNeighbor(_assembly, _var.number(), jvar.number(), type);
     252             : 
     253       13632 :   if (_local_ke.n() == 0 || _local_ke.m() == 0)
     254        3840 :     return;
     255             : 
     256        9792 :   if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
     257             :   {
     258           0 :     const auto & jv0 = static_cast<const MooseVariable &>(jvar);
     259             :     const VariableTestValue & loc_phi =
     260           0 :         (type == Moose::ElementElement || type == Moose::NeighborElement) ? jv0.phiFace()
     261           0 :                                                                           : jv0.phiFaceNeighbor();
     262             : 
     263           0 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     264             :     {
     265           0 :       initQpOffDiagJacobian(type, jvar);
     266           0 :       for (_i = 0; _i < test_space.size(); _i++)
     267           0 :         for (_j = 0; _j < loc_phi.size(); _j++)
     268             :         {
     269           0 :           RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
     270           0 :           _assembly.saveFullLocalArrayJacobian(_local_ke,
     271             :                                                _i,
     272             :                                                test_space.size(),
     273             :                                                _j,
     274             :                                                loc_phi.size(),
     275           0 :                                                _var.number(),
     276             :                                                jvar.number(),
     277             :                                                v);
     278           0 :         }
     279             :     }
     280             :   }
     281        9792 :   else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
     282             :   {
     283        9792 :     const auto & jv1 = static_cast<const ArrayMooseVariable &>(jvar);
     284             :     const ArrayVariableTestValue & loc_phi =
     285        9792 :         (type == Moose::ElementElement || type == Moose::NeighborElement) ? jv1.phiFace()
     286       14688 :                                                                           : jv1.phiFaceNeighbor();
     287             : 
     288       37440 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     289             :     {
     290       27648 :       initQpOffDiagJacobian(type, jvar);
     291      242304 :       for (_i = 0; _i < test_space.size(); _i++)
     292     2219520 :         for (_j = 0; _j < loc_phi.size(); _j++)
     293             :         {
     294     2004864 :           RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
     295     4009728 :           _assembly.saveFullLocalArrayJacobian(_local_ke,
     296             :                                                _i,
     297             :                                                test_space.size(),
     298             :                                                _j,
     299             :                                                loc_phi.size(),
     300     2004864 :                                                _var.number(),
     301             :                                                jvar.number(),
     302             :                                                v);
     303     2004864 :         }
     304             :     }
     305             :   }
     306             :   else
     307           0 :     mooseError("Vector variable cannot be coupled into array DG kernel currently");
     308             : 
     309        9792 :   accumulateTaggedLocalMatrix();
     310             : 
     311       10560 :   if (_has_diag_save_in && (type == Moose::ElementElement || type == Moose::NeighborNeighbor) &&
     312         768 :       _var.number() == jvar.number())
     313             :   {
     314         768 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     315         768 :     Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
     316        1536 :     for (const auto & var : _diag_save_in)
     317             :     {
     318         768 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     319         768 :       if (!avar)
     320           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     321             : 
     322         768 :       if (type == Moose::ElementElement)
     323         384 :         avar->addSolution(diag);
     324             :       else
     325         384 :         avar->addSolutionNeighbor(diag);
     326             :     }
     327         768 :   }
     328             : }
     329             : 
     330             : RealEigenMatrix
     331     2004864 : ArrayDGKernel::computeQpOffDiagJacobian(Moose::DGJacobianType type,
     332             :                                         const MooseVariableFEBase & jvar)
     333             : {
     334     2004864 :   if (jvar.number() == _var.number())
     335     4009728 :     return computeQpJacobian(type).asDiagonal();
     336             :   else
     337           0 :     return RealEigenMatrix::Zero(_var.count(), jvar.count());
     338             : }