www.mooseframework.org
InterfaceKernel.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "InterfaceKernel.h"
11 
12 // MOOSE includes
13 #include "Assembly.h"
14 #include "MooseVariableFE.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/quadrature.h"
18 
19 template <typename T>
22 {
24  if (std::is_same<T, Real>::value)
25  params.registerBase("InterfaceKernel");
26  else if (std::is_same<T, RealVectorValue>::value)
27  params.registerBase("VectorInterfaceKernel");
28  else
29  ::mooseError("unsupported InterfaceKernelTempl specialization");
30  return params;
31 }
32 
33 template <typename T>
35  : InterfaceKernelBase(parameters),
37  false,
39  std::is_same<T, Real>::value
42  _var(*this->mooseVariable()),
43  _normals(_assembly.normals()),
44  _u(_is_implicit ? _var.sln() : _var.slnOld()),
45  _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
46  _phi(_assembly.phiFace(_var)),
47  _grad_phi(_assembly.gradPhiFace(_var)),
48  _test(_var.phiFace()),
49  _grad_test(_var.gradPhiFace()),
50  _neighbor_var(*getVarHelper<MooseVariableFE<T>>("neighbor_var", 0)),
51  _neighbor_value(_is_implicit ? _neighbor_var.slnNeighbor() : _neighbor_var.slnOldNeighbor()),
52  _grad_neighbor_value(_neighbor_var.gradSlnNeighbor()),
53  _phi_neighbor(_assembly.phiFaceNeighbor(_neighbor_var)),
54  _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_neighbor_var)),
55  _test_neighbor(_neighbor_var.phiFaceNeighbor()),
56  _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor())
57 
58 {
60 
61  if (!parameters.isParamValid("boundary"))
62  mooseError(
63  "In order to use an interface kernel, you must specify a boundary where it will live.");
64 
65  if (parameters.isParamSetByUser("save_in"))
66  {
67  if (_save_in_strings.size() != _save_in_var_side.size())
68  mooseError("save_in and save_in_var_side must be the same length");
69  else
70  {
71  for (unsigned i = 0; i < _save_in_strings.size(); ++i)
72  {
74 
76  mooseError("Trying to use solution variable " + _save_in_strings[i] +
77  " as a save_in variable in " + name());
78 
79  if (_save_in_var_side[i] == "m")
80  {
81  if (var->feType() != _var.feType())
82  mooseError(
83  "Error in " + name() +
84  ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
85  "and the fe_type of the the primary side nonlinear "
86  "variable this interface kernel object is acting on.");
88  }
89  else
90  {
91  if (var->feType() != _neighbor_var.feType())
92  mooseError(
93  "Error in " + name() +
94  ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
95  "and the fe_type of the the secondary side nonlinear "
96  "variable this interface kernel object is acting on.");
98  }
99 
102  }
103  }
104  }
105 
108 
109  if (parameters.isParamSetByUser("diag_save_in"))
110  {
112  mooseError("diag_save_in and diag_save_in_var_side must be the same length");
113  else
114  {
115  for (unsigned i = 0; i < _diag_save_in_strings.size(); ++i)
116  {
118 
120  mooseError("Trying to use solution variable " + _diag_save_in_strings[i] +
121  " as a save_in variable in " + name());
122 
123  if (_diag_save_in_var_side[i] == "m")
124  {
125  if (var->feType() != _var.feType())
126  mooseError(
127  "Error in " + name() +
128  ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
129  "and the fe_type of the the primary side nonlinear "
130  "variable this interface kernel object is acting on.");
132  }
133  else
134  {
135  if (var->feType() != _neighbor_var.feType())
136  mooseError(
137  "Error in " + name() +
138  ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
139  "and the fe_type of the the secondary side nonlinear "
140  "variable this interface kernel object is acting on.");
142  }
143 
146  }
147  }
148  }
149 
152 }
153 
154 template <typename T>
155 void
157 {
158  bool is_elem;
159  if (type == Moose::Element)
160  is_elem = true;
161  else
162  is_elem = false;
163 
164  const TemplateVariableTestValue & test_space = is_elem ? _test : _test_neighbor;
165 
166  if (is_elem)
167  prepareVectorTag(_assembly, _var.number());
168  else
169  prepareVectorTagNeighbor(_assembly, _neighbor_var.number());
170 
171  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
172  {
173  initQpResidual(type);
174  for (_i = 0; _i < test_space.size(); _i++)
175  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual(type);
176  }
177 
178  accumulateTaggedLocalResidual();
179 
180  if (_has_primary_residuals_saved_in && is_elem)
181  {
182  Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
183  for (const auto & var : _primary_save_in_residual_variables)
184  {
185  var->sys().solution().add_vector(_local_re, var->dofIndices());
186  }
187  }
188  else if (_has_secondary_residuals_saved_in && !is_elem)
189  {
190  Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
191  for (const auto & var : _secondary_save_in_residual_variables)
192  var->sys().solution().add_vector(_local_re, var->dofIndicesNeighbor());
193  }
194 }
195 
196 template <typename T>
197 void
199 {
200  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
201  // associated only with the lower-dimensional geometric entity that is the boundary between two
202  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
203  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
204  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
205  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
206  // to execute in the former case. In the future we should remove this and add some kind of "block"
207  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
208  // this return
209  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
210  !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
211  return;
212 
213  precalculateResidual();
214 
215  // Compute the residual for this element
216  computeElemNeighResidual(Moose::Element);
217 
218  // Compute the residual for the neighbor
219  computeElemNeighResidual(Moose::Neighbor);
220 }
221 
222 template <typename T>
223 void
225 {
226  const TemplateVariableTestValue & test_space =
227  (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
228  const TemplateVariableTestValue & loc_phi =
229  (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
230 
231  unsigned int ivar, jvar;
232 
233  switch (type)
234  {
236  ivar = jvar = _var.number();
237  break;
239  ivar = _var.number(), jvar = _neighbor_var.number();
240  break;
242  ivar = _neighbor_var.number(), jvar = _var.number();
243  break;
245  ivar = _neighbor_var.number(), jvar = _neighbor_var.number();
246  break;
247  default:
248  mooseError("Unknown DGJacobianType ", type);
249  }
250 
251  if (type == Moose::ElementElement)
252  prepareMatrixTag(_assembly, ivar, jvar);
253  else
254  prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
255 
256  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
257  {
258  initQpJacobian(type);
259  for (_i = 0; _i < test_space.size(); _i++)
260  for (_j = 0; _j < loc_phi.size(); _j++)
261  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian(type);
262  }
263 
264  accumulateTaggedLocalMatrix();
265 
266  if (_has_primary_jacobians_saved_in && type == Moose::ElementElement)
267  {
268  auto rows = _local_ke.m();
269  DenseVector<Number> diag(rows);
270  for (decltype(rows) i = 0; i < rows; i++)
271  diag(i) = _local_ke(i, i);
272 
273  Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
274  for (const auto & var : _primary_save_in_jacobian_variables)
275  var->sys().solution().add_vector(diag, var->dofIndices());
276  }
277  else if (_has_secondary_jacobians_saved_in && type == Moose::NeighborNeighbor)
278  {
279  auto rows = _local_ke.m();
280  DenseVector<Number> diag(rows);
281  for (decltype(rows) i = 0; i < rows; i++)
282  diag(i) = _local_ke(i, i);
283 
284  Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
285  for (const auto & var : _secondary_save_in_jacobian_variables)
286  var->sys().solution().add_vector(diag, var->dofIndicesNeighbor());
287  }
288 }
289 
290 template <typename T>
291 void
293 {
294  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
295  // associated only with the lower-dimensional geometric entity that is the boundary between two
296  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
297  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
298  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
299  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
300  // to execute in the former case. In the future we should remove this and add some kind of "block"
301  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
302  // this return
303  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
304  !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
305  return;
306 
307  precalculateJacobian();
308 
309  computeElemNeighJacobian(Moose::ElementElement);
310  computeElemNeighJacobian(Moose::NeighborNeighbor);
311 }
312 
313 template <typename T>
314 void
316  unsigned int jvar)
317 {
318  const TemplateVariableTestValue & test_space =
319  (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
320  const TemplateVariableTestValue & loc_phi =
321  (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
322 
323  unsigned int ivar;
324 
325  if (type == Moose::ElementElement || type == Moose::ElementNeighbor)
326  ivar = _var.number();
327  else
328  ivar = _neighbor_var.number();
329 
330  if (type == Moose::ElementElement)
331  prepareMatrixTag(_assembly, ivar, jvar);
332  else
333  prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
334 
335  // Prevent calling of Jacobian computation if jvar doesn't lie in the current block
336  if ((_local_ke.m() == test_space.size()) && (_local_ke.n() == loc_phi.size()))
337  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
338  {
339  initQpOffDiagJacobian(type, jvar);
340  for (_i = 0; _i < test_space.size(); _i++)
341  for (_j = 0; _j < loc_phi.size(); _j++)
342  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
343  }
344 
345  accumulateTaggedLocalMatrix();
346 }
347 
348 template <typename T>
349 void
351 {
352  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
353  // associated only with the lower-dimensional geometric entity that is the boundary between two
354  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
355  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
356  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
357  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
358  // to execute in the former case. In the future we should remove this and add some kind of "block"
359  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
360  // this return
361  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
362  !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
363  return;
364 
365  bool is_jvar_not_interface_var = true;
366  if (jvar == _var.number())
367  {
368  precalculateJacobian();
369  computeElemNeighJacobian(Moose::ElementElement);
370  is_jvar_not_interface_var = false;
371  }
372  if (jvar == _neighbor_var.number())
373  {
374  precalculateJacobian();
375  computeElemNeighJacobian(Moose::ElementNeighbor);
376  is_jvar_not_interface_var = false;
377  }
378 
379  if (is_jvar_not_interface_var)
380  {
381  precalculateOffDiagJacobian(jvar);
382  computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
383  computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
384  }
385 }
386 
387 template <typename T>
388 void
390 {
391  // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
392  // associated only with the lower-dimensional geometric entity that is the boundary between two
393  // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
394  // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
395  // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
396  // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
397  // to execute in the former case. In the future we should remove this and add some kind of "block"
398  // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
399  // this return
400  if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
401  !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
402  return;
403 
404  bool is_jvar_not_interface_var = true;
405  if (jvar == _var.number())
406  {
407  precalculateJacobian();
408  computeElemNeighJacobian(Moose::NeighborElement);
409  is_jvar_not_interface_var = false;
410  }
411  if (jvar == _neighbor_var.number())
412  {
413  precalculateJacobian();
414  computeElemNeighJacobian(Moose::NeighborNeighbor);
415  is_jvar_not_interface_var = false;
416  }
417 
418  if (is_jvar_not_interface_var)
419  {
420  precalculateOffDiagJacobian(jvar);
421  computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
422  computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
423  }
424 }
425 
426 template <typename T>
427 void
429 {
430  computeResidual();
431 
432  if (!isImplicit())
433  return;
434 
435  for (const auto & [ivariable, jvariable] : _fe_problem.couplingEntries(_tid, _sys.number()))
436  {
437  if (ivariable->isFV())
438  continue;
439 
440  const unsigned int ivar = ivariable->number();
441  const unsigned int jvar = jvariable->number();
442 
443  prepareShapes(jvar);
444  prepareNeighborShapes(jvar);
445 
446  if (_var.number() == ivar)
447  computeElementOffDiagJacobian(jvar);
448 
449  if (_neighbor_var.number() == ivar)
450  computeNeighborOffDiagJacobian(jvar);
451  }
452 }
453 
454 // Explicitly instantiates the two versions of the InterfaceKernelTempl class
455 template class InterfaceKernelTempl<Real>;
VarFieldType
Definition: MooseTypes.h:634
std::vector< MooseVariableFEBase * > _secondary_save_in_residual_variables
The aux variables to save the secondary contributions to.
static InputParameters validParams()
MultiMooseEnum _diag_save_in_var_side
MultiMooseEnum specifying whether jacobian save-in aux variables correspond to primary or secondary s...
Class for stuff related to variables.
Definition: Adaptivity.h:31
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
InterfaceKernel and VectorInterfaceKernel is responsible for interfacing physics across subdomains...
InterfaceKernelTempl(const InputParameters &parameters)
unsigned int size() const
Return the number of active items in the MultiMooseEnum.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
std::vector< AuxVariableName > _diag_save_in_strings
The names of the aux variables that will be used to save-in jacobians (includes both primary and seco...
DGResidualType
Definition: MooseTypes.h:656
THREAD_ID _tid
The thread ID for this kernel.
MultiMooseEnum _save_in_var_side
MultiMooseEnum specifying whether residual save-in aux variables correspond to primary or secondary s...
virtual const std::string & name() const
Get the name of the class.
Definition: MooseBase.h:56
const FEType & feType() const
Get the type of finite element object.
bool _has_primary_jacobians_saved_in
Whether there are primary jacobian aux variables.
SystemBase & _sys
Reference to the EquationSystem object.
void registerBase(const std::string &value)
This method must be called from every base "Moose System" to create linkage with the Action System...
MooseVariableFE< T > * mooseVariable() const
Enhances MooseVariableInterface interface provide values from neighbor elements.
const MooseVariableFE< T > & _neighbor_var
Coupled neighbor variable.
virtual void computeElementOffDiagJacobian(unsigned int jvar) override
Selects the correct Jacobian type and routine to call for the primary variable jacobian.
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
virtual void computeElemNeighResidual(Moose::DGResidualType type)
Using the passed DGResidual type, selects the correct test function space and residual block...
virtual void computeResidualAndJacobian() override
Computes the residual and Jacobian for the current side.
std::vector< MooseVariableFEBase * > _primary_save_in_jacobian_variables
The aux variables to save the diagonal Jacobian contributions of the primary variables to...
SubProblem & _subproblem
Reference to this kernel&#39;s SubProblem.
VarKindType
Framework-wide stuff.
Definition: MooseTypes.h:627
virtual void computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, unsigned int jvar)
Using the passed DGJacobian type, selects the correct test function and trial function spaces and jac...
virtual void computeJacobian() override
Computes the jacobian for the current side.
virtual MooseVariable & getStandardVariable(const THREAD_ID tid, const std::string &var_name)=0
Returns the variable reference for requested MooseVariable which may be in any system.
std::vector< MooseVariableFEBase * > _primary_save_in_residual_variables
The aux variables to save the primary residual contributions to.
virtual bool hasVariable(const std::string &var_name) const
Query a system for a variable.
Definition: SystemBase.C:800
DGJacobianType
Definition: MooseTypes.h:662
bool isParamSetByUser(const std::string &name) const
Method returns true if the parameter was by the user.
virtual void computeElemNeighJacobian(Moose::DGJacobianType type)
Using the passed DGJacobian type, selects the correct test function and trial function spaces and jac...
void addMooseVariableDependency(MooseVariableFieldBase *var)
Call this function to add the passed in MooseVariableFieldBase as a variable that this object depends...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void addVariableToZeroOnResidual(std::string var_name)
Adds this variable to the list of variables to be zeroed during each residual evaluation.
Definition: SystemBase.C:163
std::vector< MooseVariableFEBase * > _secondary_save_in_jacobian_variables
The aux variables to save the diagonal Jacobian contributions of the secondary variables to...
bool _has_secondary_residuals_saved_in
Whether there are secondary residual aux variables.
std::vector< AuxVariableName > _save_in_strings
The names of the aux variables that will be used to save-in residuals (includes both primary and seco...
bool _has_primary_residuals_saved_in
Whether there are primary residual aux variables.
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
virtual void addVariableToZeroOnJacobian(std::string var_name)
Adds this variable to the list of variables to be zeroed during each Jacobian evaluation.
Definition: SystemBase.C:169
const InputParameters & parameters() const
Get the parameters of the object.
InterfaceKernelBase is the base class for all InterfaceKernel type classes.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
bool _has_secondary_jacobians_saved_in
Whether there are secondary jacobian aux variables.
MooseVariableFE< T > & _var
The primary side MooseVariable.
SystemBase & sys()
Get the system this variable is part of.
virtual void computeResidual() override
Computes the residual for the current side.
static InputParameters validParams()
virtual void computeNeighborOffDiagJacobian(unsigned int jvar) override
Selects the correct Jacobian type and routine to call for the secondary variable jacobian.
bool isParamValid(const std::string &name) const
This method returns parameters that have been initialized in one fashion or another, i.e.