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 "DGLowerDKernel.h"
      11             : 
      12             : #include "Assembly.h"
      13             : #include "MooseVariable.h"
      14             : #include "Problem.h"
      15             : #include "SubProblem.h"
      16             : #include "SystemBase.h"
      17             : #include "MaterialData.h"
      18             : #include "ParallelUniqueId.h"
      19             : 
      20             : #include "libmesh/dof_map.h"
      21             : #include "libmesh/dense_vector.h"
      22             : #include "libmesh/numeric_vector.h"
      23             : #include "libmesh/dense_subvector.h"
      24             : #include "libmesh/libmesh_common.h"
      25             : #include "libmesh/quadrature.h"
      26             : 
      27             : InputParameters
      28       28734 : DGLowerDKernel::validParams()
      29             : {
      30       28734 :   InputParameters params = DGKernel::validParams();
      31       28734 :   params.addRequiredCoupledVar("lowerd_variable", "Lagrange multiplier");
      32       28734 :   return params;
      33           0 : }
      34             : 
      35         105 : DGLowerDKernel::DGLowerDKernel(const InputParameters & parameters)
      36             :   : DGKernel(parameters),
      37         105 :     _lowerd_var(*getVar("lowerd_variable", 0)),
      38         105 :     _lambda(_is_implicit ? _lowerd_var.slnLower() : _lowerd_var.slnLowerOld()),
      39             : 
      40         105 :     _phi_lambda(_lowerd_var.phiLower()),
      41         210 :     _test_lambda(_lowerd_var.phiLower())
      42             : {
      43         105 :   const auto & lower_domains = _lowerd_var.activeSubdomains();
      44         223 :   for (const auto & id : _mesh.interiorLowerDBlocks())
      45         118 :     if (lower_domains.count(id) == 0)
      46           0 :       mooseDocumentedError(
      47             :           "moose",
      48             :           29151,
      49           0 :           "'lowerd_variable' must be defined on the interior lower-dimensional subdomain '" +
      50           0 :               _mesh.getSubdomainName(id) +
      51             :               "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
      52             :               "restrictive.");
      53             : 
      54         249 :   for (const auto & id : _var.activeSubdomains())
      55         144 :     if (_mesh.interiorLowerDBlocks().count(id) > 0)
      56           0 :       paramError("variable",
      57           0 :                  "Must not be defined on the interior lower-dimensional subdomain'" +
      58           0 :                      _mesh.getSubdomainName(id) + "'");
      59             : 
      60             :   // Note: the above two conditions also ensure that the variable and lower-d variable are
      61             :   // different.
      62         105 : }
      63             : 
      64             : void
      65        2528 : DGLowerDKernel::computeResidual()
      66             : {
      67        2528 :   if (!excludeBoundary())
      68             :   {
      69        2528 :     precalculateResidual();
      70             : 
      71             :     // Compute the residual for this element
      72        2528 :     computeElemNeighResidual(Moose::Element);
      73             : 
      74             :     // Compute the residual for the neighbor
      75        2528 :     computeElemNeighResidual(Moose::Neighbor);
      76             : 
      77        2528 :     computeLowerDResidual();
      78             :   }
      79        2528 : }
      80             : 
      81             : void
      82        2528 : DGLowerDKernel::computeLowerDResidual()
      83             : {
      84        2528 :   prepareVectorTagLower(_assembly, _lowerd_var.number());
      85             : 
      86       27104 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      87             :   {
      88       24576 :     initLowerDQpResidual();
      89       49152 :     for (_i = 0; _i < _test_lambda.size(); _i++)
      90       24576 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeLowerDQpResidual();
      91             :   }
      92             : 
      93        2528 :   accumulateTaggedLocalResidual();
      94        2528 : }
      95             : 
      96             : void
      97           0 : DGLowerDKernel::computeJacobian()
      98             : {
      99           0 :   if (!excludeBoundary())
     100             :   {
     101           0 :     precalculateJacobian();
     102             : 
     103             :     // Compute element-element Jacobian
     104           0 :     computeElemNeighJacobian(Moose::ElementElement);
     105             : 
     106             :     // Compute element-neighbor Jacobian
     107           0 :     computeElemNeighJacobian(Moose::ElementNeighbor);
     108             : 
     109             :     // Compute neighbor-element Jacobian
     110           0 :     computeElemNeighJacobian(Moose::NeighborElement);
     111             : 
     112             :     // Compute neighbor-neighbor Jacobian
     113           0 :     computeElemNeighJacobian(Moose::NeighborNeighbor);
     114             : 
     115             :     // Compute the other five pieces of Jacobian related with lower-d variable
     116           0 :     computeLowerDJacobian(Moose::LowerLower);
     117             :   }
     118           0 : }
     119             : 
     120             : void
     121        6320 : DGLowerDKernel::computeLowerDJacobian(Moose::ConstraintJacobianType type)
     122             : {
     123             :   mooseAssert(type != Moose::PrimaryPrimary && type != Moose::PrimarySecondary &&
     124             :                   type != Moose::SecondaryPrimary && type != Moose::SecondarySecondary,
     125             :               "Jacobian types without lower should be handled in computeElemNeighJacobian");
     126             : 
     127       11376 :   const auto & test_space =
     128        5056 :       (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
     129             :           ? _test_lambda
     130             :           : (type == Moose::PrimaryLower ? _test : _test_neighbor);
     131             :   const auto ivar =
     132        5056 :       (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
     133        8848 :           ? _lowerd_var.number()
     134        2528 :           : _var.number();
     135             : 
     136       11376 :   const auto & loc_phi =
     137        5056 :       (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
     138             :           ? _phi_lambda
     139             :           : (type == Moose::LowerPrimary ? _phi : _phi_neighbor);
     140             :   const auto jvar =
     141        5056 :       (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
     142        8848 :           ? _lowerd_var.number()
     143        2528 :           : _var.number();
     144             : 
     145        6320 :   prepareMatrixTagLower(_assembly, ivar, jvar, type);
     146             : 
     147        6320 :   if (_local_ke.n() == 0 || _local_ke.m() == 0)
     148           0 :     return;
     149             : 
     150       67760 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     151             :   {
     152       61440 :     initLowerDQpJacobian(type);
     153      538624 :     for (_i = 0; _i < test_space.size(); _i++)
     154     1370112 :       for (_j = 0; _j < loc_phi.size(); _j++)
     155      892928 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeLowerDQpJacobian(type);
     156             :   }
     157             : 
     158        6320 :   accumulateTaggedLocalMatrix();
     159             : }
     160             : 
     161             : void
     162        2528 : DGLowerDKernel::computeOffDiagJacobian(const unsigned int jvar_num)
     163             : {
     164        2528 :   if (!excludeBoundary())
     165             :   {
     166        2528 :     precalculateOffDiagJacobian(jvar_num);
     167             : 
     168        2528 :     if (jvar_num == variable().number())
     169             :     {
     170        1264 :       computeElemNeighJacobian(Moose::ElementElement);
     171        1264 :       computeElemNeighJacobian(Moose::ElementNeighbor);
     172        1264 :       computeElemNeighJacobian(Moose::NeighborElement);
     173        1264 :       computeElemNeighJacobian(Moose::NeighborNeighbor);
     174        1264 :       computeLowerDJacobian(Moose::LowerSecondary);
     175        1264 :       computeLowerDJacobian(Moose::LowerPrimary);
     176             :     }
     177        1264 :     else if (jvar_num == _lowerd_var.number())
     178             :     {
     179        1264 :       computeLowerDJacobian(Moose::LowerLower);
     180        1264 :       computeLowerDJacobian(Moose::SecondaryLower);
     181        1264 :       computeLowerDJacobian(Moose::PrimaryLower);
     182             :     }
     183             :     else
     184             :     {
     185           0 :       const auto & jvar = getVariable(jvar_num);
     186             : 
     187             :       // Compute element-element Jacobian
     188           0 :       computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
     189             : 
     190             :       // Compute element-neighbor Jacobian
     191           0 :       computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
     192             : 
     193             :       // Compute neighbor-element Jacobian
     194           0 :       computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
     195             : 
     196             :       // Compute neighbor-neighbor Jacobian
     197           0 :       computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
     198             : 
     199             :       // Compute the other five pieces of Jacobian related with lower-d variable
     200           0 :       computeOffDiagLowerDJacobian(Moose::LowerLower, jvar);
     201           0 :       computeOffDiagLowerDJacobian(Moose::LowerSecondary, jvar);
     202           0 :       computeOffDiagLowerDJacobian(Moose::LowerPrimary, jvar);
     203           0 :       computeOffDiagLowerDJacobian(Moose::SecondaryLower, jvar);
     204           0 :       computeOffDiagLowerDJacobian(Moose::PrimaryLower, jvar);
     205             :     }
     206             :   }
     207        2528 : }
     208             : 
     209             : void
     210           0 : DGLowerDKernel::computeOffDiagLowerDJacobian(Moose::ConstraintJacobianType type,
     211             :                                              const MooseVariableFEBase & jvar)
     212             : {
     213             :   mooseAssert(type != Moose::PrimaryPrimary && type != Moose::PrimarySecondary &&
     214             :                   type != Moose::SecondaryPrimary && type != Moose::SecondarySecondary,
     215             :               "Jacobian types without lower should be handled in computeOffDiagElemNeighJacobian");
     216             : 
     217           0 :   const auto & test_space =
     218           0 :       (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
     219             :           ? _test_lambda
     220             :           : (type == Moose::PrimaryLower ? _test : _test_neighbor);
     221             :   const auto ivar =
     222           0 :       (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
     223           0 :           ? _lowerd_var.number()
     224           0 :           : _var.number();
     225             : 
     226           0 :   prepareMatrixTagLower(_assembly, ivar, jvar.number(), type);
     227           0 :   if (_local_ke.n() == 0 || _local_ke.m() == 0)
     228           0 :     return;
     229             : 
     230           0 :   if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
     231             :   {
     232           0 :     const auto & jv0 = static_cast<const MooseVariable &>(jvar);
     233             :     const VariableTestValue & loc_phi =
     234           0 :         (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
     235           0 :             ? jv0.phiLower()
     236           0 :             : (type == Moose::LowerPrimary ? jv0.phiFace() : jv0.phiFaceNeighbor());
     237             : 
     238           0 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     239             :     {
     240           0 :       initLowerDQpOffDiagJacobian(type, jvar);
     241           0 :       for (_i = 0; _i < test_space.size(); _i++)
     242           0 :         for (_j = 0; _j < loc_phi.size(); _j++)
     243           0 :           _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeLowerDQpOffDiagJacobian(type, jvar);
     244             :     }
     245             :   }
     246           0 :   else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
     247           0 :     mooseError("Array variable cannot be coupled into DG kernel currently");
     248             :   else
     249           0 :     mooseError("Vector variable cannot be coupled into DG kernel currently");
     250             : 
     251           0 :   accumulateTaggedLocalMatrix();
     252             : }