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 "NonlinearThread.h"
11 : #include "NonlinearSystem.h"
12 : #include "Problem.h"
13 : #include "FEProblem.h"
14 : #include "KernelBase.h"
15 : #include "IntegratedBCBase.h"
16 : #include "DGKernelBase.h"
17 : #include "InterfaceKernelBase.h"
18 : #include "Material.h"
19 : #include "TimeKernel.h"
20 : #include "SwapBackSentinel.h"
21 : #include "FVTimeKernel.h"
22 : #include "ComputeJacobianThread.h"
23 :
24 : #include "libmesh/threads.h"
25 :
26 3528072 : NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
27 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
28 7056144 : _nl(fe_problem.currentNonlinearSystem()),
29 3528072 : _num_cached(0),
30 3528072 : _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31 3528072 : _dg_kernels(_nl.getDGKernelWarehouse()),
32 3528072 : _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33 3528072 : _kernels(_nl.getKernelWarehouse()),
34 3528072 : _hdg_kernels(_nl.getHDGKernelWarehouse()),
35 6861205 : _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
36 7068522 : _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
37 207317 : _fe_problem.haveFV()),
38 7056144 : _should_execute_dg(false)
39 : {
40 3528072 : }
41 :
42 : // Splitting Constructor
43 1232 : NonlinearThread::NonlinearThread(NonlinearThread & x, Threads::split split)
44 : : ThreadedElementLoop<ConstElemRange>(x, split),
45 1232 : _nl(x._nl),
46 1232 : _num_cached(x._num_cached),
47 1232 : _integrated_bcs(x._integrated_bcs),
48 1232 : _dg_kernels(x._dg_kernels),
49 1232 : _interface_kernels(x._interface_kernels),
50 1232 : _kernels(x._kernels),
51 1232 : _tag_kernels(x._tag_kernels),
52 1232 : _hdg_kernels(x._hdg_kernels),
53 1232 : _has_active_objects(x._has_active_objects),
54 1232 : _should_execute_dg(x._should_execute_dg)
55 : {
56 1232 : }
57 :
58 3529276 : NonlinearThread::~NonlinearThread() {}
59 :
60 : void
61 3528809 : NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
62 : {
63 3528809 : if (_has_active_objects)
64 3404759 : ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
65 3528787 : }
66 :
67 : void
68 4208096 : NonlinearThread::subdomainChanged()
69 : {
70 : // This should come first to setup the residual objects before we do dependency determination of
71 : // material properties and variables
72 4208096 : determineObjectWarehouses();
73 :
74 4208096 : _fe_problem.subdomainSetup(_subdomain, _tid);
75 :
76 : // Update variable Dependencies
77 4208096 : std::set<MooseVariableFEBase *> needed_moose_vars;
78 4208096 : _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
79 4208096 : _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
80 4208096 : _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
81 4208096 : _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
82 :
83 : // Update FE variable coupleable vector tags
84 4208096 : std::set<TagID> needed_fe_var_vector_tags;
85 4208096 : _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
86 4208096 : _subdomain, needed_fe_var_vector_tags, _tid);
87 4208096 : _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
88 4208096 : _subdomain, needed_fe_var_vector_tags, _tid);
89 4208096 : _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
90 4208096 : _subdomain, needed_fe_var_vector_tags, _tid);
91 :
92 : // Update material dependencies
93 4208096 : std::unordered_set<unsigned int> needed_mat_props;
94 4208096 : _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
95 4208096 : _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
96 4208096 : _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
97 4208096 : _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
98 :
99 4208096 : if (_fe_problem.haveFV())
100 : {
101 : // Re-query the finite volume elemental kernels
102 141893 : _fv_kernels.clear();
103 141893 : _fe_problem.theWarehouse()
104 141893 : .query()
105 283786 : .template condition<AttribSysNum>(_nl.number())
106 141893 : .template condition<AttribSystem>("FVElementalKernel")
107 141893 : .template condition<AttribSubdomains>(_subdomain)
108 141893 : .template condition<AttribThread>(_tid)
109 141893 : .queryInto(_fv_kernels);
110 248031 : for (const auto fv_kernel : _fv_kernels)
111 : {
112 106138 : const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
113 106138 : needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
114 106138 : const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
115 106138 : needed_mat_props.insert(fv_mp_deps.begin(), fv_mp_deps.end());
116 : }
117 : }
118 :
119 : // Cache these to avoid computing them on every side
120 4208096 : _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
121 4208096 : _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
122 :
123 4208096 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
124 4208096 : _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
125 4208096 : _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
126 4208096 : }
127 :
128 : void
129 334818187 : NonlinearThread::onElement(const Elem * const elem)
130 : {
131 : // Set up Sentinel class so that, even if reinitMaterials() throws in prepareElement, we
132 : // still remember to swap back during stack unwinding.
133 334818187 : SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
134 :
135 334818187 : prepareElement(elem);
136 :
137 334818121 : if (dynamic_cast<ComputeJacobianThread *>(this))
138 46636777 : if (_nl.getScalarVariables(_tid).size() > 0)
139 158530 : _fe_problem.reinitOffDiagScalars(_tid);
140 :
141 334818121 : computeOnElement();
142 334818167 : }
143 :
144 : void
145 320614727 : NonlinearThread::computeOnElement()
146 : {
147 320614727 : if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
148 : {
149 312554351 : const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
150 1035214025 : for (const auto & kernel : kernels)
151 722659739 : compute(*kernel);
152 : }
153 :
154 320614662 : if (_fe_problem.haveFV())
155 13940796 : for (auto kernel : _fv_kernels)
156 6085869 : compute(*kernel);
157 320614662 : }
158 :
159 : void
160 104896837 : NonlinearThread::onBoundary(const Elem * const elem,
161 : const unsigned int side,
162 : const BoundaryID bnd_id,
163 : const Elem * const lower_d_elem /*=nullptr*/)
164 : {
165 104896837 : if (_ibc_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
166 : {
167 : // Set up Sentinel class so that, after we swap in reinitMaterialsFace in prepareFace, even if
168 : // one of our callees throws we remember to swap back during stack unwinding. We put our
169 : // sentinel here as opposed to in prepareFace because we certainly don't want our materials
170 : // swapped back before we proceed to residual/Jacobian computation
171 2806068 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
172 :
173 2806068 : prepareFace(elem, side, bnd_id, lower_d_elem);
174 2806068 : computeOnBoundary(bnd_id, lower_d_elem);
175 :
176 2806066 : if (lower_d_elem)
177 22928 : accumulateLower();
178 2806066 : }
179 104896835 : }
180 :
181 : void
182 2617068 : NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
183 : {
184 2617068 : const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
185 5474707 : for (const auto & bc : bcs)
186 2857641 : if (bc->shouldApply())
187 2857464 : compute(*bc);
188 2617066 : }
189 :
190 : void
191 422604 : NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
192 : {
193 422604 : if (_ik_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
194 : {
195 :
196 : // Pointer to the neighbor we are currently working on.
197 43931 : const Elem * neighbor = elem->neighbor_ptr(side);
198 :
199 43931 : if (neighbor->active())
200 : {
201 43931 : _fe_problem.reinitNeighbor(elem, side, _tid);
202 :
203 : // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
204 : // still remember to swap back during stack unwinding. Note that face, boundary, and interface
205 : // all operate with the same MaterialData object
206 43931 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
207 43931 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
208 43931 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
209 :
210 43931 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
211 43931 : _fe_problem.reinitMaterialsNeighborOnBoundary(bnd_id, neighbor->subdomain_id(), _tid);
212 :
213 : // Has to happen after face and neighbor properties have been computed. Note that we don't use
214 : // a sentinel here because FEProblem::swapBackMaterialsFace is going to handle face materials,
215 : // boundary materials, and interface materials (e.g. it queries the boundary material data
216 : // with the current element and side
217 43931 : _fe_problem.reinitMaterialsInterface(bnd_id, _tid);
218 :
219 43931 : computeOnInterface(bnd_id);
220 :
221 43931 : accumulateNeighbor();
222 43931 : }
223 : }
224 422604 : }
225 :
226 : void
227 37317 : NonlinearThread::computeOnInterface(BoundaryID bnd_id)
228 : {
229 37317 : const auto & int_ks = _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
230 99880 : for (const auto & interface_kernel : int_ks)
231 62563 : compute(*interface_kernel);
232 37317 : }
233 :
234 : void
235 3326524 : NonlinearThread::onInternalSide(const Elem * elem, unsigned int side)
236 : {
237 3326524 : if (_should_execute_dg)
238 : {
239 : // Pointer to the neighbor we are currently working on.
240 2058260 : const Elem * neighbor = elem->neighbor_ptr(side);
241 :
242 2058260 : _fe_problem.reinitElemNeighborAndLowerD(elem, side, _tid);
243 :
244 : // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
245 : // still remember to swap back during stack unwinding.
246 2058260 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
247 2058260 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
248 :
249 2058260 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
250 2058260 : _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
251 :
252 2058260 : computeOnInternalFace(neighbor);
253 :
254 2058260 : accumulateNeighborLower();
255 2058260 : }
256 3326524 : if (_subdomain_has_hdg)
257 : {
258 : // Set up Sentinel class so that, after we swap in reinitMaterialsFace in prepareFace, even if
259 : // one of our callees throws we remember to swap back during stack unwinding. We put our
260 : // sentinel here as opposed to in prepareFace because we certainly don't want our materials
261 : // swapped back before we proceed to residual/Jacobian computation
262 1268264 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
263 :
264 1268264 : prepareFace(elem, side, Moose::INVALID_BOUNDARY_ID, nullptr);
265 1268264 : computeOnInternalFace();
266 1268264 : }
267 3326524 : }
268 :
269 : void
270 106273144 : NonlinearThread::onExternalSide(const Elem * elem, unsigned int side)
271 : {
272 : // Check that we don't have any interface kernels defined on this boundary
273 106273144 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
274 :
275 106273144 : bool has_interface_kernels = false;
276 210747375 : for (const auto bid : boundary_ids)
277 104474231 : if (_ik_warehouse && _ik_warehouse->hasActiveBoundaryObjects(bid, _tid))
278 0 : has_interface_kernels = true;
279 :
280 106273144 : if (has_interface_kernels)
281 0 : mooseError("Element ",
282 0 : elem->id(),
283 : " on side ",
284 : side,
285 : " is missing a neighbor (hence identified as an external side) but "
286 : "has interface kernel(s) defined on the boundary.");
287 106273144 : }
288 :
289 : void
290 2013540 : NonlinearThread::computeOnInternalFace(const Elem * neighbor)
291 : {
292 2013540 : const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
293 4697729 : for (const auto & dg_kernel : dgks)
294 2684189 : if (dg_kernel->hasBlocks(neighbor->subdomain_id()))
295 2681675 : compute(*dg_kernel, neighbor);
296 2013540 : }
297 :
298 : void
299 646484980 : NonlinearThread::compute(KernelBase & kernel)
300 : {
301 646484980 : compute(static_cast<ResidualObject &>(kernel));
302 646484936 : }
303 :
304 : void
305 6085669 : NonlinearThread::compute(FVElementalKernel & kernel)
306 : {
307 6085669 : compute(static_cast<ResidualObject &>(kernel));
308 6085669 : }
309 :
310 : void
311 2738490 : NonlinearThread::compute(IntegratedBCBase & bc)
312 : {
313 2738490 : compute(static_cast<ResidualObject &>(bc));
314 2738490 : }
315 :
316 : void
317 2597067 : NonlinearThread::compute(DGKernelBase & dg, const Elem * /*neighbor*/)
318 : {
319 2597067 : compute(static_cast<ResidualObject &>(dg));
320 2597067 : }
321 :
322 : void
323 60287 : NonlinearThread::compute(InterfaceKernelBase & ik)
324 : {
325 60287 : compute(static_cast<ResidualObject &>(ik));
326 60287 : }
327 :
328 : void
329 288181300 : NonlinearThread::postElement(const Elem * /*elem*/)
330 : {
331 288181300 : accumulate();
332 288181300 : }
333 :
334 : void
335 3405121 : NonlinearThread::post()
336 : {
337 3405121 : clearVarsAndMaterials();
338 3405121 : }
339 :
340 : void
341 3404759 : NonlinearThread::printGeneralExecutionInformation() const
342 : {
343 3404759 : if (_fe_problem.shouldPrintExecution(_tid))
344 : {
345 10947 : const auto & console = _fe_problem.console();
346 10947 : const auto execute_on = _fe_problem.getCurrentExecuteOnFlag();
347 21894 : console << "[DBG] Beginning elemental loop to compute " + objectType() + " on " << execute_on
348 10947 : << std::endl;
349 10947 : mooseDoOnce(
350 : console << "[DBG] Execution order on each element:" << std::endl;
351 : console << "[DBG] - kernels on element quadrature points" << std::endl;
352 : console << "[DBG] - finite volume elemental kernels on element" << std::endl;
353 : console << "[DBG] - integrated boundary conditions on element side quadrature points"
354 : << std::endl;
355 : console << "[DBG] - DG kernels on element side quadrature points" << std::endl;
356 : console << "[DBG] - interface kernels on element side quadrature points" << std::endl;);
357 10947 : }
358 3404759 : }
359 :
360 : void
361 4207532 : NonlinearThread::printBlockExecutionInformation() const
362 : {
363 : // Number of objects executing is approximated by size of warehouses
364 4207532 : const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
365 4207532 : _dg_kernels.size() + _interface_kernels.size();
366 4207532 : const auto & console = _fe_problem.console();
367 4207532 : const auto block_name = _mesh.getSubdomainName(_subdomain);
368 :
369 4207532 : if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
370 : {
371 68681 : if (_blocks_exec_printed.count(_subdomain))
372 52114 : return;
373 33134 : console << "[DBG] Ordering of " + objectType() + " Objects on block " << block_name << " ("
374 16567 : << _subdomain << ")" << std::endl;
375 16567 : if (_kernels.hasActiveBlockObjects(_subdomain, _tid))
376 : {
377 16503 : console << "[DBG] Ordering of kernels:" << std::endl;
378 49509 : console << _kernels.activeObjectsToFormattedString() << std::endl;
379 : }
380 16567 : if (_fv_kernels.size())
381 : {
382 64 : console << "[DBG] Ordering of FV elemental kernels:" << std::endl;
383 : std::string fvkernels =
384 64 : std::accumulate(_fv_kernels.begin() + 1,
385 : _fv_kernels.end(),
386 64 : _fv_kernels[0]->name(),
387 48 : [](const std::string & str_out, FVElementalKernel * kernel)
388 176 : { return str_out + " " + kernel->name(); });
389 64 : console << ConsoleUtils::formatString(fvkernels, "[DBG]") << std::endl;
390 64 : }
391 16567 : if (_dg_kernels.hasActiveBlockObjects(_subdomain, _tid))
392 : {
393 5902 : console << "[DBG] Ordering of DG kernels:" << std::endl;
394 17706 : console << _dg_kernels.activeObjectsToFormattedString() << std::endl;
395 : }
396 : }
397 4138851 : else if (_fe_problem.shouldPrintExecution(_tid) && num_objects == 0 &&
398 0 : !_blocks_exec_printed.count(_subdomain))
399 0 : console << "[DBG] No Active " + objectType() + " Objects on block " << block_name << " ("
400 0 : << _subdomain << ")" << std::endl;
401 :
402 4155418 : _blocks_exec_printed.insert(_subdomain);
403 4207532 : }
404 :
405 : void
406 104896837 : NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
407 : {
408 105000785 : if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
409 103948 : (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
410 82051 : !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
411 104872485 : return;
412 :
413 24352 : const auto & console = _fe_problem.console();
414 24352 : const auto b_name = _mesh.getBoundaryName(bid);
415 48704 : console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
416 24352 : << ")" << std::endl;
417 :
418 24352 : if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
419 : {
420 21897 : console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
421 65691 : console << _integrated_bcs.activeObjectsToFormattedString() << std::endl;
422 : }
423 :
424 : // We have not checked if we have a neighbor. This could be premature for saying we are executing
425 : // interface kernels. However, we should assume the execution will happen on another side of the
426 : // same boundary
427 24352 : if (_interface_kernels.hasActiveBoundaryObjects(bid, _tid))
428 : {
429 2455 : console << "[DBG] Ordering of interface kernels:" << std::endl;
430 7365 : console << _interface_kernels.activeObjectsToFormattedString() << std::endl;
431 : }
432 :
433 24352 : _boundaries_exec_printed.insert(bid);
434 24352 : }
435 :
436 : void
437 4074332 : NonlinearThread::prepareFace(const Elem * const elem,
438 : const unsigned int side,
439 : const BoundaryID bnd_id,
440 : const Elem * const lower_d_elem)
441 : {
442 4074332 : _fe_problem.reinitElemFace(elem, side, _tid);
443 :
444 : // Needed to use lower-dimensional variables on Materials
445 4074332 : if (lower_d_elem)
446 22928 : _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
447 :
448 4074332 : if (bnd_id != Moose::INVALID_BOUNDARY_ID)
449 : {
450 2806068 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
451 2806068 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
452 : }
453 : // Currently only used by HDG
454 : else
455 1268264 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
456 4074332 : }
457 :
458 : bool
459 1272136875 : NonlinearThread::shouldComputeInternalSide(const Elem & elem, const Elem & neighbor) const
460 : {
461 : // ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide gets expensive on high
462 : // h-refinement cases so we avoid it if possible
463 1272136875 : _should_execute_dg = false;
464 1272136875 : if (_subdomain_has_dg)
465 4114145 : _should_execute_dg =
466 4114145 : ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
467 1272136875 : return _subdomain_has_hdg || _should_execute_dg;
468 : }
|
