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 3522101 : NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
27 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
28 7044202 : _nl(fe_problem.currentNonlinearSystem()),
29 3522101 : _num_cached(0),
30 3522101 : _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31 3522101 : _dg_kernels(_nl.getDGKernelWarehouse()),
32 3522101 : _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33 3522101 : _kernels(_nl.getKernelWarehouse()),
34 3522101 : _hdg_kernels(_nl.getHDGKernelWarehouse()),
35 6850119 : _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
36 7057395 : _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
37 207276 : _fe_problem.haveFV()),
38 7044202 : _should_execute_dg(false)
39 : {
40 3522101 : }
41 :
42 : // Splitting Constructor
43 350357 : NonlinearThread::NonlinearThread(NonlinearThread & x, Threads::split split)
44 : : ThreadedElementLoop<ConstElemRange>(x, split),
45 350357 : _nl(x._nl),
46 350357 : _num_cached(x._num_cached),
47 350357 : _integrated_bcs(x._integrated_bcs),
48 350357 : _dg_kernels(x._dg_kernels),
49 350357 : _interface_kernels(x._interface_kernels),
50 350357 : _kernels(x._kernels),
51 350357 : _tag_kernels(x._tag_kernels),
52 350357 : _hdg_kernels(x._hdg_kernels),
53 350357 : _has_active_objects(x._has_active_objects),
54 350357 : _should_execute_dg(x._should_execute_dg)
55 : {
56 350357 : }
57 :
58 3872427 : NonlinearThread::~NonlinearThread() {}
59 :
60 : void
61 3871914 : NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
62 : {
63 3871914 : if (_has_active_objects)
64 3737674 : ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
65 3871890 : }
66 :
67 : void
68 4533501 : 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 4533501 : determineObjectWarehouses();
73 :
74 4533501 : _fe_problem.subdomainSetup(_subdomain, _tid);
75 :
76 : // Update variable Dependencies
77 4533501 : std::set<MooseVariableFEBase *> needed_moose_vars;
78 4533501 : _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
79 4533501 : _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
80 4533501 : _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
81 4533501 : _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
82 :
83 : // Update FE variable coupleable vector tags
84 4533501 : std::set<TagID> needed_fe_var_vector_tags;
85 4533501 : _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
86 4533501 : _subdomain, needed_fe_var_vector_tags, _tid);
87 4533501 : _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
88 4533501 : _subdomain, needed_fe_var_vector_tags, _tid);
89 4533501 : _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
90 4533501 : _subdomain, needed_fe_var_vector_tags, _tid);
91 :
92 : // Update material dependencies
93 4533501 : std::unordered_set<unsigned int> needed_mat_props;
94 4533501 : _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
95 4533501 : _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
96 4533501 : _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
97 4533501 : _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
98 :
99 4533501 : if (_fe_problem.haveFV())
100 : {
101 : // Re-query the finite volume elemental kernels
102 148429 : _fv_kernels.clear();
103 148429 : _fe_problem.theWarehouse()
104 148429 : .query()
105 296858 : .template condition<AttribSysNum>(_nl.number())
106 148429 : .template condition<AttribSystem>("FVElementalKernel")
107 148429 : .template condition<AttribSubdomains>(_subdomain)
108 148429 : .template condition<AttribThread>(_tid)
109 148429 : .queryInto(_fv_kernels);
110 260545 : for (const auto fv_kernel : _fv_kernels)
111 : {
112 112116 : const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
113 112116 : needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
114 112116 : const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
115 112116 : 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 4533501 : _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
121 4533501 : _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
122 :
123 4533501 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
124 4533501 : _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
125 4533501 : _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
126 4533501 : }
127 :
128 : void
129 334535188 : 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 334535188 : SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
134 :
135 334535188 : prepareElement(elem);
136 :
137 334535122 : if (dynamic_cast<ComputeJacobianThread *>(this))
138 46599764 : if (_nl.getScalarVariables(_tid).size() > 0)
139 158306 : _fe_problem.reinitOffDiagScalars(_tid);
140 :
141 334535122 : computeOnElement();
142 334535166 : }
143 :
144 : void
145 320334472 : NonlinearThread::computeOnElement()
146 : {
147 320334472 : if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
148 : {
149 312274076 : const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
150 1034357582 : for (const auto & kernel : kernels)
151 722083577 : compute(*kernel);
152 : }
153 :
154 320334401 : if (_fe_problem.haveFV())
155 13940796 : for (auto kernel : _fv_kernels)
156 6085869 : compute(*kernel);
157 320334401 : }
158 :
159 : void
160 104752174 : 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 104752174 : 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 2802414 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
172 :
173 2802414 : prepareFace(elem, side, bnd_id, lower_d_elem);
174 2802414 : computeOnBoundary(bnd_id, lower_d_elem);
175 :
176 2802412 : if (lower_d_elem)
177 22928 : accumulateLower();
178 2802412 : }
179 104752172 : }
180 :
181 : void
182 2613526 : NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
183 : {
184 2613526 : const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
185 5467303 : for (const auto & bc : bcs)
186 2853779 : if (bc->shouldApply())
187 2853602 : compute(*bc);
188 2613524 : }
189 :
190 : void
191 422609 : NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
192 : {
193 422609 : 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 422609 : }
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 106128528 : NonlinearThread::onExternalSide(const Elem * elem, unsigned int side)
271 : {
272 : // Check that we don't have any interface kernels defined on this boundary
273 106128528 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
274 :
275 106128528 : bool has_interface_kernels = false;
276 210458091 : for (const auto bid : boundary_ids)
277 104329563 : if (_ik_warehouse && _ik_warehouse->hasActiveBoundaryObjects(bid, _tid))
278 0 : has_interface_kernels = true;
279 :
280 106128528 : 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 106128528 : }
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 645978357 : NonlinearThread::compute(KernelBase & kernel)
300 : {
301 645978357 : compute(static_cast<ResidualObject &>(kernel));
302 645978308 : }
303 :
304 : void
305 6085669 : NonlinearThread::compute(FVElementalKernel & kernel)
306 : {
307 6085669 : compute(static_cast<ResidualObject &>(kernel));
308 6085669 : }
309 :
310 : void
311 2734756 : NonlinearThread::compute(IntegratedBCBase & bc)
312 : {
313 2734756 : compute(static_cast<ResidualObject &>(bc));
314 2734756 : }
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 287935309 : NonlinearThread::postElement(const Elem * /*elem*/)
330 : {
331 287935309 : accumulate();
332 287935309 : }
333 :
334 : void
335 3738079 : NonlinearThread::post()
336 : {
337 3738079 : clearVarsAndMaterials();
338 3738079 : }
339 :
340 : void
341 3737674 : NonlinearThread::printGeneralExecutionInformation() const
342 : {
343 3737674 : 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 3737674 : }
359 :
360 : void
361 4532890 : NonlinearThread::printBlockExecutionInformation() const
362 : {
363 : // Number of objects executing is approximated by size of warehouses
364 4532890 : const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
365 4532890 : _dg_kernels.size() + _interface_kernels.size();
366 4532890 : const auto & console = _fe_problem.console();
367 4532890 : const auto block_name = _mesh.getSubdomainName(_subdomain);
368 :
369 4532890 : if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
370 : {
371 62777 : if (_blocks_exec_printed.count(_subdomain))
372 46210 : 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 4470113 : 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 4486680 : _blocks_exec_printed.insert(_subdomain);
403 4532890 : }
404 :
405 : void
406 104752174 : NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
407 : {
408 104849495 : if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
409 97321 : (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
410 76470 : !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
411 104728868 : return;
412 :
413 23306 : const auto & console = _fe_problem.console();
414 23306 : const auto b_name = _mesh.getBoundaryName(bid);
415 46612 : console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
416 23306 : << ")" << std::endl;
417 :
418 23306 : if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
419 : {
420 20851 : console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
421 62553 : 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 23306 : 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 23306 : _boundaries_exec_printed.insert(bid);
434 23306 : }
435 :
436 : void
437 4070678 : 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 4070678 : _fe_problem.reinitElemFace(elem, side, _tid);
443 :
444 : // Needed to use lower-dimensional variables on Materials
445 4070678 : if (lower_d_elem)
446 22928 : _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
447 :
448 4070678 : if (bnd_id != Moose::INVALID_BOUNDARY_ID)
449 : {
450 2802414 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
451 2802414 : _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 4070678 : }
457 :
458 : bool
459 1271159839 : 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 1271159839 : _should_execute_dg = false;
464 1271159839 : if (_subdomain_has_dg)
465 4114145 : _should_execute_dg =
466 4114145 : ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
467 1271159839 : return _subdomain_has_hdg || _should_execute_dg;
468 : }
|
