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 "InterfaceKernel.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "MooseVariableFE.h"
15 : #include "SystemBase.h"
16 : #include "AuxiliarySystem.h"
17 : #include "FEProblemBase.h"
18 :
19 : #include "libmesh/quadrature.h"
20 :
21 : template <typename T>
22 : InputParameters
23 22856 : InterfaceKernelTempl<T>::validParams()
24 : {
25 22856 : InputParameters params = InterfaceKernelBase::validParams();
26 : if (std::is_same<T, Real>::value)
27 39506 : params.registerBase("InterfaceKernel");
28 : else if (std::is_same<T, RealVectorValue>::value)
29 6206 : params.registerBase("VectorInterfaceKernel");
30 : else
31 : ::mooseError("unsupported InterfaceKernelTempl specialization");
32 22856 : return params;
33 0 : }
34 :
35 : template <typename T>
36 752 : InterfaceKernelTempl<T>::InterfaceKernelTempl(const InputParameters & parameters)
37 : : InterfaceKernelBase(parameters),
38 : NeighborMooseVariableInterface<T>(this,
39 : false,
40 : Moose::VarKindType::VAR_SOLVER,
41 : std::is_same<T, Real>::value
42 : ? Moose::VarFieldType::VAR_FIELD_STANDARD
43 : : Moose::VarFieldType::VAR_FIELD_VECTOR),
44 1504 : _var(*this->mooseVariable()),
45 752 : _normals(_assembly.normals()),
46 752 : _u(_is_implicit ? _var.sln() : _var.slnOld()),
47 752 : _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
48 752 : _phi(_assembly.phiFace(_var)),
49 752 : _grad_phi(_assembly.gradPhiFace(_var)),
50 752 : _test(_var.phiFace()),
51 752 : _grad_test(_var.gradPhiFace()),
52 1504 : _neighbor_var(*getVarHelper<MooseVariableFE<T>>("neighbor_var", 0)),
53 752 : _neighbor_value(_is_implicit ? _neighbor_var.slnNeighbor() : _neighbor_var.slnOldNeighbor()),
54 752 : _grad_neighbor_value(_neighbor_var.gradSlnNeighbor()),
55 752 : _phi_neighbor(_assembly.phiFaceNeighbor(_neighbor_var)),
56 752 : _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_neighbor_var)),
57 752 : _test_neighbor(_neighbor_var.phiFaceNeighbor()),
58 752 : _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor()),
59 1504 : _same_system(_var.sys().number() == _neighbor_var.sys().number())
60 : {
61 : // Neighbor variable dependency is added by
62 : // NeighborCoupleableMooseVariableDependencyIntermediateInterface
63 752 : addMooseVariableDependency(this->mooseVariable());
64 :
65 1504 : if (!parameters.isParamValid("boundary"))
66 0 : mooseError(
67 : "In order to use an interface kernel, you must specify a boundary where it will live.");
68 :
69 1504 : if (parameters.isParamSetByUser("save_in"))
70 : {
71 26 : if (_save_in_strings.size() != _save_in_var_side.size())
72 0 : mooseError("save_in and save_in_var_side must be the same length");
73 : else
74 : {
75 78 : for (unsigned i = 0; i < _save_in_strings.size(); ++i)
76 : {
77 52 : MooseVariable * var = &_subproblem.getStandardVariable(_tid, _save_in_strings[i]);
78 :
79 52 : if (_sys.hasVariable(_save_in_strings[i]))
80 0 : mooseError("Trying to use solution variable " + _save_in_strings[i] +
81 0 : " as a save_in variable in " + name());
82 :
83 52 : if (_save_in_var_side[i] == "m")
84 : {
85 26 : if (var->feType() != _var.feType())
86 0 : mooseError(
87 0 : "Error in " + name() +
88 : ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
89 : "and the fe_type of the the primary side nonlinear "
90 : "variable this interface kernel object is acting on.");
91 26 : _primary_save_in_residual_variables.push_back(var);
92 : }
93 : else
94 : {
95 26 : if (var->feType() != _neighbor_var.feType())
96 0 : mooseError(
97 0 : "Error in " + name() +
98 : ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
99 : "and the fe_type of the the secondary side nonlinear "
100 : "variable this interface kernel object is acting on.");
101 26 : _secondary_save_in_residual_variables.push_back(var);
102 : }
103 :
104 52 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
105 52 : addMooseVariableDependency(var);
106 : }
107 : }
108 : }
109 :
110 752 : _has_primary_residuals_saved_in = _primary_save_in_residual_variables.size() > 0;
111 752 : _has_secondary_residuals_saved_in = _secondary_save_in_residual_variables.size() > 0;
112 :
113 1504 : if (parameters.isParamSetByUser("diag_save_in"))
114 : {
115 26 : if (_diag_save_in_strings.size() != _diag_save_in_var_side.size())
116 0 : mooseError("diag_save_in and diag_save_in_var_side must be the same length");
117 : else
118 : {
119 78 : for (unsigned i = 0; i < _diag_save_in_strings.size(); ++i)
120 : {
121 52 : MooseVariable * var = &_subproblem.getStandardVariable(_tid, _diag_save_in_strings[i]);
122 :
123 52 : if (_sys.hasVariable(_diag_save_in_strings[i]))
124 0 : mooseError("Trying to use solution variable " + _diag_save_in_strings[i] +
125 0 : " as a save_in variable in " + name());
126 :
127 52 : if (_diag_save_in_var_side[i] == "m")
128 : {
129 26 : if (var->feType() != _var.feType())
130 0 : mooseError(
131 0 : "Error in " + name() +
132 : ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
133 : "and the fe_type of the the primary side nonlinear "
134 : "variable this interface kernel object is acting on.");
135 26 : _primary_save_in_jacobian_variables.push_back(var);
136 : }
137 : else
138 : {
139 26 : if (var->feType() != _neighbor_var.feType())
140 0 : mooseError(
141 0 : "Error in " + name() +
142 : ". There is a mismatch between the fe_type of the save-in Auxiliary variable "
143 : "and the fe_type of the the secondary side nonlinear "
144 : "variable this interface kernel object is acting on.");
145 26 : _secondary_save_in_jacobian_variables.push_back(var);
146 : }
147 :
148 52 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
149 52 : addMooseVariableDependency(var);
150 : }
151 : }
152 : }
153 :
154 752 : _has_primary_jacobians_saved_in = _primary_save_in_jacobian_variables.size() > 0;
155 752 : _has_secondary_jacobians_saved_in = _secondary_save_in_jacobian_variables.size() > 0;
156 752 : }
157 :
158 : template <typename T>
159 : void
160 117886 : InterfaceKernelTempl<T>::computeElemNeighResidual(Moose::DGResidualType type)
161 : {
162 : bool is_elem;
163 117886 : if (type == Moose::Element)
164 58999 : is_elem = true;
165 : else
166 58887 : is_elem = false;
167 :
168 117886 : const TemplateVariableTestValue & test_space = is_elem ? _test : _test_neighbor;
169 :
170 117886 : if (is_elem)
171 58999 : prepareVectorTag(_assembly, _var.number());
172 : else
173 58887 : prepareVectorTagNeighbor(_assembly, _neighbor_var.number());
174 :
175 538884 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
176 : {
177 420998 : initQpResidual(type);
178 5292754 : for (_i = 0; _i < test_space.size(); _i++)
179 4871756 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual(type);
180 : }
181 :
182 117886 : accumulateTaggedLocalResidual();
183 :
184 : // To save the diagonal of the Jacobian
185 117886 : if (_has_primary_residuals_saved_in && is_elem)
186 64 : for (const auto & var : _primary_save_in_residual_variables)
187 64 : var->sys().solution().add_vector(_local_re, var->dofIndices());
188 117854 : else if (_has_secondary_residuals_saved_in && !is_elem)
189 64 : for (const auto & var : _secondary_save_in_residual_variables)
190 32 : var->sys().solution().add_vector(_local_re, var->dofIndicesNeighbor());
191 117886 : }
192 :
193 : template <typename T>
194 : void
195 59123 : InterfaceKernelTempl<T>::computeResidual()
196 : {
197 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
198 : // associated only with the lower-dimensional geometric entity that is the boundary between two
199 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
200 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
201 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
202 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
203 : // to execute in the former case. In the future we should remove this and add some kind of "block"
204 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
205 : // this return
206 118198 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
207 59075 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
208 48 : return;
209 :
210 59075 : precalculateResidual();
211 :
212 : // Compute the residual for this element
213 59075 : computeElemNeighResidual(Moose::Element);
214 :
215 : // Compute the residual for the neighbor
216 : // This also prevents computing a residual if the neighbor variable is auxiliary
217 59075 : if (_same_system)
218 58963 : computeElemNeighResidual(Moose::Neighbor);
219 : }
220 :
221 : template <typename T>
222 : void
223 31092 : InterfaceKernelTempl<T>::computeElemNeighJacobian(Moose::DGJacobianType type)
224 : {
225 31092 : const TemplateVariableTestValue & test_space =
226 31092 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
227 31092 : const TemplateVariableTestValue & loc_phi =
228 31092 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
229 :
230 : unsigned int ivar, jvar;
231 :
232 31092 : switch (type)
233 : {
234 8933 : case Moose::ElementElement:
235 8933 : ivar = jvar = _var.number();
236 8933 : break;
237 6633 : case Moose::ElementNeighbor:
238 6633 : ivar = _var.number(), jvar = _neighbor_var.number();
239 6633 : break;
240 6633 : case Moose::NeighborElement:
241 6633 : ivar = _neighbor_var.number(), jvar = _var.number();
242 6633 : break;
243 8893 : case Moose::NeighborNeighbor:
244 8893 : ivar = _neighbor_var.number(), jvar = _neighbor_var.number();
245 8893 : break;
246 0 : default:
247 0 : mooseError("Unknown DGJacobianType ", type);
248 : }
249 :
250 31092 : if (type == Moose::ElementElement)
251 8933 : prepareMatrixTag(_assembly, ivar, jvar);
252 : else
253 22159 : prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
254 :
255 101700 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
256 : {
257 70608 : initQpJacobian(type);
258 510320 : for (_i = 0; _i < test_space.size(); _i++)
259 4376832 : for (_j = 0; _j < loc_phi.size(); _j++)
260 3937120 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian(type);
261 : }
262 :
263 31092 : accumulateTaggedLocalMatrix();
264 :
265 : // To save the diagonal of the Jacobian
266 31092 : if (_has_primary_jacobians_saved_in && type == Moose::ElementElement)
267 : {
268 16 : auto rows = _local_ke.m();
269 16 : DenseVector<Number> diag(rows);
270 48 : for (decltype(rows) i = 0; i < rows; i++)
271 32 : diag(i) = _local_ke(i, i);
272 :
273 32 : for (const auto & var : _primary_save_in_jacobian_variables)
274 16 : var->sys().solution().add_vector(diag, var->dofIndices());
275 16 : }
276 31076 : else if (_has_secondary_jacobians_saved_in && type == Moose::NeighborNeighbor)
277 : {
278 16 : auto rows = _local_ke.m();
279 16 : DenseVector<Number> diag(rows);
280 48 : for (decltype(rows) i = 0; i < rows; i++)
281 32 : diag(i) = _local_ke(i, i);
282 :
283 32 : for (const auto & var : _secondary_save_in_jacobian_variables)
284 16 : var->sys().solution().add_vector(diag, var->dofIndicesNeighbor());
285 16 : }
286 31092 : }
287 :
288 : template <typename T>
289 : void
290 2276 : InterfaceKernelTempl<T>::computeJacobian()
291 : {
292 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
293 : // associated only with the lower-dimensional geometric entity that is the boundary between two
294 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
295 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
296 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
297 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
298 : // to execute in the former case. In the future we should remove this and add some kind of "block"
299 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
300 : // this return
301 4552 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
302 2276 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
303 0 : return;
304 :
305 2276 : precalculateJacobian();
306 :
307 2276 : computeElemNeighJacobian(Moose::ElementElement);
308 2276 : if (_same_system)
309 2260 : computeElemNeighJacobian(Moose::NeighborNeighbor);
310 : }
311 :
312 : template <typename T>
313 : void
314 0 : InterfaceKernelTempl<T>::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
315 : unsigned int jvar)
316 : {
317 0 : const TemplateVariableTestValue & test_space =
318 0 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
319 0 : const TemplateVariableTestValue & loc_phi =
320 0 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
321 :
322 : unsigned int ivar;
323 :
324 0 : if (type == Moose::ElementElement || type == Moose::ElementNeighbor)
325 0 : ivar = _var.number();
326 : else
327 0 : ivar = _neighbor_var.number();
328 :
329 0 : if (type == Moose::ElementElement)
330 0 : prepareMatrixTag(_assembly, ivar, jvar);
331 : else
332 0 : prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
333 :
334 : // Prevent calling of Jacobian computation if jvar doesn't lie in the current block
335 0 : if ((_local_ke.m() == test_space.size()) && (_local_ke.n() == loc_phi.size()))
336 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
337 : {
338 0 : initQpOffDiagJacobian(type, jvar);
339 0 : for (_i = 0; _i < test_space.size(); _i++)
340 0 : for (_j = 0; _j < loc_phi.size(); _j++)
341 0 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
342 : }
343 :
344 0 : accumulateTaggedLocalMatrix();
345 0 : }
346 :
347 : template <typename T>
348 : void
349 13439 : InterfaceKernelTempl<T>::computeElementOffDiagJacobian(unsigned int jvar)
350 : {
351 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
352 : // associated only with the lower-dimensional geometric entity that is the boundary between two
353 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
354 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
355 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
356 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
357 : // to execute in the former case. In the future we should remove this and add some kind of "block"
358 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
359 : // this return
360 26814 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
361 13375 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
362 64 : return;
363 :
364 13375 : bool is_jvar_not_interface_var = true;
365 13375 : if (jvar == _var.number())
366 : {
367 6701 : precalculateJacobian();
368 6701 : computeElemNeighJacobian(Moose::ElementElement);
369 6701 : is_jvar_not_interface_var = false;
370 : }
371 13375 : if (jvar == _neighbor_var.number() && _same_system)
372 : {
373 6677 : precalculateJacobian();
374 6677 : computeElemNeighJacobian(Moose::ElementNeighbor);
375 6677 : is_jvar_not_interface_var = false;
376 : }
377 :
378 13375 : if (is_jvar_not_interface_var)
379 : {
380 0 : precalculateOffDiagJacobian(jvar);
381 0 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
382 0 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
383 : }
384 : }
385 :
386 : template <typename T>
387 : void
388 13415 : InterfaceKernelTempl<T>::computeNeighborOffDiagJacobian(unsigned int jvar)
389 : {
390 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
391 : // associated only with the lower-dimensional geometric entity that is the boundary between two
392 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
393 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
394 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
395 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
396 : // to execute in the former case. In the future we should remove this and add some kind of "block"
397 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
398 : // this return
399 26766 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
400 13351 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
401 64 : return;
402 :
403 : // We don't care about any contribution to the neighbor Jacobian rows if it's not in the system
404 : // we are currently working with (the variable's system)
405 13351 : if (!_same_system)
406 0 : return;
407 :
408 13351 : bool is_jvar_not_interface_var = true;
409 13351 : if (jvar == _var.number())
410 : {
411 6677 : precalculateJacobian();
412 6677 : computeElemNeighJacobian(Moose::NeighborElement);
413 6677 : is_jvar_not_interface_var = false;
414 : }
415 13351 : if (jvar == _neighbor_var.number())
416 : {
417 6677 : precalculateJacobian();
418 6677 : computeElemNeighJacobian(Moose::NeighborNeighbor);
419 6677 : is_jvar_not_interface_var = false;
420 : }
421 :
422 13351 : if (is_jvar_not_interface_var)
423 : {
424 0 : precalculateOffDiagJacobian(jvar);
425 0 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
426 0 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
427 : }
428 : }
429 :
430 : template <typename T>
431 : void
432 62 : InterfaceKernelTempl<T>::computeResidualAndJacobian()
433 : {
434 62 : computeResidual();
435 :
436 62 : if (!isImplicit())
437 0 : return;
438 :
439 238 : for (const auto & [ivariable, jvariable] : _fe_problem.couplingEntries(_tid, _sys.number()))
440 : {
441 176 : if (ivariable->isFV())
442 0 : continue;
443 :
444 176 : const unsigned int ivar = ivariable->number();
445 176 : const unsigned int jvar = jvariable->number();
446 :
447 176 : prepareShapes(jvar);
448 176 : prepareNeighborShapes(jvar);
449 :
450 176 : if (_var.number() == ivar)
451 100 : computeElementOffDiagJacobian(jvar);
452 :
453 176 : if (_same_system)
454 152 : if (_neighbor_var.number() == ivar)
455 76 : computeNeighborOffDiagJacobian(jvar);
456 : }
457 : }
458 :
459 : // Explicitly instantiates the two versions of the InterfaceKernelTempl class
460 : template class InterfaceKernelTempl<Real>;
461 : template class InterfaceKernelTempl<RealVectorValue>;
|
