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 101402 : InterfaceKernelTempl<T>::validParams()
24 : {
25 101402 : InputParameters params = InterfaceKernelBase::validParams();
26 : if (std::is_same<T, Real>::value)
27 87095 : params.registerBase("InterfaceKernel");
28 : else if (std::is_same<T, RealVectorValue>::value)
29 14307 : params.registerBase("VectorInterfaceKernel");
30 : else
31 : ::mooseError("unsupported InterfaceKernelTempl specialization");
32 101402 : return params;
33 0 : }
34 :
35 : template <typename T>
36 807 : 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 1614 : _var(*this->mooseVariable()),
45 807 : _normals(_assembly.normals()),
46 807 : _u(_is_implicit ? _var.sln() : _var.slnOld()),
47 807 : _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
48 807 : _phi(_assembly.phiFace(_var)),
49 807 : _grad_phi(_assembly.gradPhiFace(_var)),
50 807 : _test(_var.phiFace()),
51 807 : _grad_test(_var.gradPhiFace()),
52 807 : _neighbor_var(*getVarHelper<MooseVariableFE<T>>("neighbor_var", 0)),
53 807 : _neighbor_value(_is_implicit ? _neighbor_var.slnNeighbor() : _neighbor_var.slnOldNeighbor()),
54 807 : _grad_neighbor_value(_neighbor_var.gradSlnNeighbor()),
55 807 : _phi_neighbor(_assembly.phiFaceNeighbor(_neighbor_var)),
56 807 : _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_neighbor_var)),
57 807 : _test_neighbor(_neighbor_var.phiFaceNeighbor()),
58 807 : _grad_test_neighbor(_neighbor_var.gradPhiFaceNeighbor()),
59 1614 : _same_system(_var.sys().number() == _neighbor_var.sys().number())
60 : {
61 : // Neighbor variable dependency is added by
62 : // NeighborCoupleableMooseVariableDependencyIntermediateInterface
63 807 : addMooseVariableDependency(this->mooseVariable());
64 :
65 807 : 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 807 : 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 807 : _has_primary_residuals_saved_in = _primary_save_in_residual_variables.size() > 0;
111 807 : _has_secondary_residuals_saved_in = _secondary_save_in_residual_variables.size() > 0;
112 :
113 807 : 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 807 : _has_primary_jacobians_saved_in = _primary_save_in_jacobian_variables.size() > 0;
155 807 : _has_secondary_jacobians_saved_in = _secondary_save_in_jacobian_variables.size() > 0;
156 807 : }
157 :
158 : template <typename T>
159 : void
160 146690 : InterfaceKernelTempl<T>::computeElemNeighResidual(Moose::DGResidualType type)
161 : {
162 : bool is_elem;
163 146690 : if (type == Moose::Element)
164 73401 : is_elem = true;
165 : else
166 73289 : is_elem = false;
167 :
168 146690 : const TemplateVariableTestValue & test_space = is_elem ? _test : _test_neighbor;
169 :
170 146690 : if (is_elem)
171 73401 : prepareVectorTag(_assembly, _var.number());
172 : else
173 73289 : prepareVectorTagNeighbor(_assembly, _neighbor_var.number());
174 :
175 682868 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
176 : {
177 536178 : initQpResidual(type);
178 6856198 : for (_i = 0; _i < test_space.size(); _i++)
179 6320020 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual(type);
180 : }
181 :
182 146690 : accumulateTaggedLocalResidual();
183 :
184 : // To save the diagonal of the Jacobian
185 146690 : if (_has_primary_residuals_saved_in && is_elem)
186 : {
187 32 : Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
188 64 : for (const auto & var : _primary_save_in_residual_variables)
189 : {
190 32 : var->sys().solution().add_vector(_local_re, var->dofIndices());
191 : }
192 32 : }
193 146658 : else if (_has_secondary_residuals_saved_in && !is_elem)
194 : {
195 32 : Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
196 64 : for (const auto & var : _secondary_save_in_residual_variables)
197 32 : var->sys().solution().add_vector(_local_re, var->dofIndicesNeighbor());
198 32 : }
199 146690 : }
200 :
201 : template <typename T>
202 : void
203 73525 : InterfaceKernelTempl<T>::computeResidual()
204 : {
205 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
206 : // associated only with the lower-dimensional geometric entity that is the boundary between two
207 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
208 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
209 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
210 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
211 : // to execute in the former case. In the future we should remove this and add some kind of "block"
212 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
213 : // this return
214 147002 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
215 73477 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
216 48 : return;
217 :
218 73477 : precalculateResidual();
219 :
220 : // Compute the residual for this element
221 73477 : computeElemNeighResidual(Moose::Element);
222 :
223 : // Compute the residual for the neighbor
224 : // This also prevents computing a residual if the neighbor variable is auxiliary
225 73477 : if (_same_system)
226 73365 : computeElemNeighResidual(Moose::Neighbor);
227 : }
228 :
229 : template <typename T>
230 : void
231 32576 : InterfaceKernelTempl<T>::computeElemNeighJacobian(Moose::DGJacobianType type)
232 : {
233 32576 : const TemplateVariableTestValue & test_space =
234 32576 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
235 32576 : const TemplateVariableTestValue & loc_phi =
236 32576 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
237 :
238 : unsigned int ivar, jvar;
239 :
240 32576 : switch (type)
241 : {
242 9644 : case Moose::ElementElement:
243 9644 : ivar = jvar = _var.number();
244 9644 : break;
245 6664 : case Moose::ElementNeighbor:
246 6664 : ivar = _var.number(), jvar = _neighbor_var.number();
247 6664 : break;
248 6664 : case Moose::NeighborElement:
249 6664 : ivar = _neighbor_var.number(), jvar = _var.number();
250 6664 : break;
251 9604 : case Moose::NeighborNeighbor:
252 9604 : ivar = _neighbor_var.number(), jvar = _neighbor_var.number();
253 9604 : break;
254 0 : default:
255 0 : mooseError("Unknown DGJacobianType ", type);
256 : }
257 :
258 32576 : if (type == Moose::ElementElement)
259 9644 : prepareMatrixTag(_assembly, ivar, jvar);
260 : else
261 22932 : prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
262 :
263 109052 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
264 : {
265 76476 : initQpJacobian(type);
266 589860 : for (_i = 0; _i < test_space.size(); _i++)
267 5446472 : for (_j = 0; _j < loc_phi.size(); _j++)
268 4933088 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian(type);
269 : }
270 :
271 32576 : accumulateTaggedLocalMatrix();
272 :
273 : // To save the diagonal of the Jacobian
274 32576 : if (_has_primary_jacobians_saved_in && type == Moose::ElementElement)
275 : {
276 16 : auto rows = _local_ke.m();
277 16 : DenseVector<Number> diag(rows);
278 48 : for (decltype(rows) i = 0; i < rows; i++)
279 32 : diag(i) = _local_ke(i, i);
280 :
281 16 : Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
282 32 : for (const auto & var : _primary_save_in_jacobian_variables)
283 16 : var->sys().solution().add_vector(diag, var->dofIndices());
284 16 : }
285 32560 : else if (_has_secondary_jacobians_saved_in && type == Moose::NeighborNeighbor)
286 : {
287 16 : auto rows = _local_ke.m();
288 16 : DenseVector<Number> diag(rows);
289 48 : for (decltype(rows) i = 0; i < rows; i++)
290 32 : diag(i) = _local_ke(i, i);
291 :
292 16 : Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
293 32 : for (const auto & var : _secondary_save_in_jacobian_variables)
294 16 : var->sys().solution().add_vector(diag, var->dofIndicesNeighbor());
295 16 : }
296 32576 : }
297 :
298 : template <typename T>
299 : void
300 2956 : InterfaceKernelTempl<T>::computeJacobian()
301 : {
302 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
303 : // associated only with the lower-dimensional geometric entity that is the boundary between two
304 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
305 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
306 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
307 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
308 : // to execute in the former case. In the future we should remove this and add some kind of "block"
309 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
310 : // this return
311 5912 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
312 2956 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
313 0 : return;
314 :
315 2956 : precalculateJacobian();
316 :
317 2956 : computeElemNeighJacobian(Moose::ElementElement);
318 2956 : if (_same_system)
319 2940 : computeElemNeighJacobian(Moose::NeighborNeighbor);
320 : }
321 :
322 : template <typename T>
323 : void
324 0 : InterfaceKernelTempl<T>::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
325 : unsigned int jvar)
326 : {
327 0 : const TemplateVariableTestValue & test_space =
328 0 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
329 0 : const TemplateVariableTestValue & loc_phi =
330 0 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
331 :
332 : unsigned int ivar;
333 :
334 0 : if (type == Moose::ElementElement || type == Moose::ElementNeighbor)
335 0 : ivar = _var.number();
336 : else
337 0 : ivar = _neighbor_var.number();
338 :
339 0 : if (type == Moose::ElementElement)
340 0 : prepareMatrixTag(_assembly, ivar, jvar);
341 : else
342 0 : prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
343 :
344 : // Prevent calling of Jacobian computation if jvar doesn't lie in the current block
345 0 : if ((_local_ke.m() == test_space.size()) && (_local_ke.n() == loc_phi.size()))
346 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
347 : {
348 0 : initQpOffDiagJacobian(type, jvar);
349 0 : for (_i = 0; _i < test_space.size(); _i++)
350 0 : for (_j = 0; _j < loc_phi.size(); _j++)
351 0 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
352 : }
353 :
354 0 : accumulateTaggedLocalMatrix();
355 0 : }
356 :
357 : template <typename T>
358 : void
359 13497 : InterfaceKernelTempl<T>::computeElementOffDiagJacobian(unsigned int jvar)
360 : {
361 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
362 : // associated only with the lower-dimensional geometric entity that is the boundary between two
363 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
364 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
365 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
366 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
367 : // to execute in the former case. In the future we should remove this and add some kind of "block"
368 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
369 : // this return
370 26930 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
371 13433 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
372 64 : return;
373 :
374 13433 : bool is_jvar_not_interface_var = true;
375 13433 : if (jvar == _var.number())
376 : {
377 6732 : precalculateJacobian();
378 6732 : computeElemNeighJacobian(Moose::ElementElement);
379 6732 : is_jvar_not_interface_var = false;
380 : }
381 13433 : if (jvar == _neighbor_var.number() && _same_system)
382 : {
383 6708 : precalculateJacobian();
384 6708 : computeElemNeighJacobian(Moose::ElementNeighbor);
385 6708 : is_jvar_not_interface_var = false;
386 : }
387 :
388 13433 : if (is_jvar_not_interface_var)
389 : {
390 0 : precalculateOffDiagJacobian(jvar);
391 0 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
392 0 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
393 : }
394 : }
395 :
396 : template <typename T>
397 : void
398 13473 : InterfaceKernelTempl<T>::computeNeighborOffDiagJacobian(unsigned int jvar)
399 : {
400 : // in the gmsh mesh format (at least in the version 2 format) the "sideset" physical entities are
401 : // associated only with the lower-dimensional geometric entity that is the boundary between two
402 : // higher-dimensional element faces. It does not have a sidedness to it like the exodus format
403 : // does. Consequently we may naively try to execute an interface kernel twice, one time where _var
404 : // has dofs on _current_elem *AND* _neighbor_var has dofs on _neighbor_elem, and the other time
405 : // where _var has dofs on _neighbor_elem and _neighbor_var has dofs on _current_elem. We only want
406 : // to execute in the former case. In the future we should remove this and add some kind of "block"
407 : // awareness to interface kernels to avoid all the unnecessary reinit that happens before we hit
408 : // this return
409 26882 : if (!_var.activeOnSubdomain(_current_elem->subdomain_id()) ||
410 13409 : !_neighbor_var.activeOnSubdomain(_neighbor_elem->subdomain_id()))
411 64 : return;
412 :
413 : // We don't care about any contribution to the neighbor Jacobian rows if it's not in the system
414 : // we are currently working with (the variable's system)
415 13409 : if (!_same_system)
416 0 : return;
417 :
418 13409 : bool is_jvar_not_interface_var = true;
419 13409 : if (jvar == _var.number())
420 : {
421 6708 : precalculateJacobian();
422 6708 : computeElemNeighJacobian(Moose::NeighborElement);
423 6708 : is_jvar_not_interface_var = false;
424 : }
425 13409 : if (jvar == _neighbor_var.number())
426 : {
427 6708 : precalculateJacobian();
428 6708 : computeElemNeighJacobian(Moose::NeighborNeighbor);
429 6708 : is_jvar_not_interface_var = false;
430 : }
431 :
432 13409 : if (is_jvar_not_interface_var)
433 : {
434 0 : precalculateOffDiagJacobian(jvar);
435 0 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
436 0 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
437 : }
438 : }
439 :
440 : template <typename T>
441 : void
442 62 : InterfaceKernelTempl<T>::computeResidualAndJacobian()
443 : {
444 62 : computeResidual();
445 :
446 62 : if (!isImplicit())
447 0 : return;
448 :
449 238 : for (const auto & [ivariable, jvariable] : _fe_problem.couplingEntries(_tid, _sys.number()))
450 : {
451 176 : if (ivariable->isFV())
452 0 : continue;
453 :
454 176 : const unsigned int ivar = ivariable->number();
455 176 : const unsigned int jvar = jvariable->number();
456 :
457 176 : prepareShapes(jvar);
458 176 : prepareNeighborShapes(jvar);
459 :
460 176 : if (_var.number() == ivar)
461 100 : computeElementOffDiagJacobian(jvar);
462 :
463 176 : if (_same_system)
464 152 : if (_neighbor_var.number() == ivar)
465 76 : computeNeighborOffDiagJacobian(jvar);
466 : }
467 : }
468 :
469 : // Explicitly instantiates the two versions of the InterfaceKernelTempl class
470 : template class InterfaceKernelTempl<Real>;
471 : template class InterfaceKernelTempl<RealVectorValue>;
|
