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 3535193 : NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
27 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
28 7070386 : _nl(fe_problem.currentNonlinearSystem()),
29 3535193 : _num_cached(0),
30 3535193 : _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31 3535193 : _dg_kernels(_nl.getDGKernelWarehouse()),
32 3535193 : _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33 3535193 : _kernels(_nl.getKernelWarehouse()),
34 3535193 : _hdg_kernels(_nl.getHDGKernelWarehouse()),
35 6875526 : _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
36 7082878 : _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
37 207352 : _fe_problem.haveFV()),
38 7070386 : _should_execute_dg(false)
39 : {
40 3535193 : }
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 3536397 : NonlinearThread::~NonlinearThread() {}
59 :
60 : void
61 3535930 : NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
62 : {
63 3535930 : if (_has_active_objects)
64 3411861 : ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
65 3535908 : }
66 :
67 : void
68 4215672 : 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 4215672 : determineObjectWarehouses();
73 :
74 4215672 : _fe_problem.subdomainSetup(_subdomain, _tid);
75 :
76 : // Update variable Dependencies
77 4215672 : std::set<MooseVariableFEBase *> needed_moose_vars;
78 4215672 : _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
79 4215672 : _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
80 4215672 : _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
81 4215672 : _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
82 :
83 : // Update FE variable coupleable vector tags
84 4215672 : std::set<TagID> needed_fe_var_vector_tags;
85 4215672 : _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
86 4215672 : _subdomain, needed_fe_var_vector_tags, _tid);
87 4215672 : _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
88 4215672 : _subdomain, needed_fe_var_vector_tags, _tid);
89 4215672 : _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
90 4215672 : _subdomain, needed_fe_var_vector_tags, _tid);
91 :
92 : // Update material dependencies
93 4215672 : std::unordered_set<unsigned int> needed_mat_props;
94 4215672 : _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
95 4215672 : _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
96 4215672 : _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
97 4215672 : _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
98 :
99 4215672 : if (_fe_problem.haveFV())
100 : {
101 : // Re-query the finite volume elemental kernels
102 141909 : _fv_kernels.clear();
103 141909 : _fe_problem.theWarehouse()
104 141909 : .query()
105 283818 : .template condition<AttribSysNum>(_nl.number())
106 141909 : .template condition<AttribSystem>("FVElementalKernel")
107 141909 : .template condition<AttribSubdomains>(_subdomain)
108 141909 : .template condition<AttribThread>(_tid)
109 141909 : .queryInto(_fv_kernels);
110 248063 : for (const auto fv_kernel : _fv_kernels)
111 : {
112 106154 : const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
113 106154 : needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
114 106154 : const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
115 106154 : 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 4215672 : _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
121 4215672 : _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
122 :
123 4215672 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
124 4215672 : _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
125 4215672 : _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
126 4215672 : }
127 :
128 : void
129 335565545 : 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 335565545 : SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
134 :
135 335565545 : prepareElement(elem);
136 :
137 335565479 : if (dynamic_cast<ComputeJacobianThread *>(this))
138 46743198 : if (_nl.getScalarVariables(_tid).size() > 0)
139 158554 : _fe_problem.reinitOffDiagScalars(_tid);
140 :
141 335565479 : computeOnElement();
142 335565525 : }
143 :
144 : void
145 321276612 : NonlinearThread::computeOnElement()
146 : {
147 321276612 : if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
148 : {
149 313199705 : const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
150 1037248987 : for (const auto & kernel : kernels)
151 724049347 : compute(*kernel);
152 : }
153 :
154 321276547 : if (_fe_problem.haveFV())
155 13945776 : for (auto kernel : _fv_kernels)
156 6088359 : compute(*kernel);
157 321276547 : }
158 :
159 : void
160 105114596 : 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 105114596 : 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 2868825 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
172 :
173 2868825 : prepareFace(elem, side, bnd_id, lower_d_elem);
174 2868825 : computeOnBoundary(bnd_id, lower_d_elem);
175 :
176 2868823 : if (lower_d_elem)
177 22916 : accumulateLower();
178 2868823 : }
179 105114594 : }
180 :
181 : void
182 2658220 : NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
183 : {
184 2658220 : const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
185 5556587 : for (const auto & bc : bcs)
186 2898369 : if (bc->shouldApply())
187 2898192 : compute(*bc);
188 2658218 : }
189 :
190 : void
191 422500 : NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
192 : {
193 422500 : if (_ik_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
194 : {
195 :
196 : // Pointer to the neighbor we are currently working on.
197 43901 : const Elem * neighbor = elem->neighbor_ptr(side);
198 :
199 43901 : if (neighbor->active())
200 : {
201 43901 : _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 43901 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
207 43901 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
208 43901 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
209 :
210 43901 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
211 43901 : _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 43901 : _fe_problem.reinitMaterialsInterface(bnd_id, _tid);
218 :
219 43901 : computeOnInterface(bnd_id);
220 :
221 43901 : accumulateNeighbor();
222 43901 : }
223 : }
224 422500 : }
225 :
226 : void
227 37287 : NonlinearThread::computeOnInterface(BoundaryID bnd_id)
228 : {
229 37287 : const auto & int_ks = _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
230 99790 : for (const auto & interface_kernel : int_ks)
231 62503 : compute(*interface_kernel);
232 37287 : }
233 :
234 : void
235 3322220 : NonlinearThread::onInternalSide(const Elem * elem, unsigned int side)
236 : {
237 3322220 : if (_should_execute_dg)
238 : {
239 : // Pointer to the neighbor we are currently working on.
240 2057060 : const Elem * neighbor = elem->neighbor_ptr(side);
241 :
242 2057060 : _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 2057060 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
247 2057060 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
248 :
249 2057060 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
250 2057060 : _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
251 :
252 2057060 : computeOnInternalFace(neighbor);
253 :
254 2057060 : accumulateNeighborLower();
255 2057060 : }
256 3322220 : 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 1265160 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
263 :
264 1265160 : prepareFace(elem, side, Moose::INVALID_BOUNDARY_ID, nullptr);
265 1265160 : computeOnInternalFace();
266 1265160 : }
267 3322220 : }
268 :
269 : void
270 106496894 : NonlinearThread::onExternalSide(const Elem * elem, unsigned int side)
271 : {
272 : // Check that we don't have any interface kernels defined on this boundary
273 106496894 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
274 :
275 106496894 : bool has_interface_kernels = false;
276 211188988 : for (const auto bid : boundary_ids)
277 104692094 : if (_ik_warehouse && _ik_warehouse->hasActiveBoundaryObjects(bid, _tid))
278 0 : has_interface_kernels = true;
279 :
280 106496894 : 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 106496894 : }
288 :
289 : void
290 2012556 : NonlinearThread::computeOnInternalFace(const Elem * neighbor)
291 : {
292 2012556 : const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
293 4694789 : for (const auto & dg_kernel : dgks)
294 2682233 : if (dg_kernel->hasBlocks(neighbor->subdomain_id()))
295 2679719 : compute(*dg_kernel, neighbor);
296 2012556 : }
297 :
298 : void
299 647835010 : NonlinearThread::compute(KernelBase & kernel)
300 : {
301 647835010 : compute(static_cast<ResidualObject &>(kernel));
302 647834966 : }
303 :
304 : void
305 6088159 : NonlinearThread::compute(FVElementalKernel & kernel)
306 : {
307 6088159 : compute(static_cast<ResidualObject &>(kernel));
308 6088159 : }
309 :
310 : void
311 2779262 : NonlinearThread::compute(IntegratedBCBase & bc)
312 : {
313 2779262 : compute(static_cast<ResidualObject &>(bc));
314 2779262 : }
315 :
316 : void
317 2595111 : NonlinearThread::compute(DGKernelBase & dg, const Elem * /*neighbor*/)
318 : {
319 2595111 : compute(static_cast<ResidualObject &>(dg));
320 2595111 : }
321 :
322 : void
323 60227 : NonlinearThread::compute(InterfaceKernelBase & ik)
324 : {
325 60227 : compute(static_cast<ResidualObject &>(ik));
326 60227 : }
327 :
328 : void
329 288822237 : NonlinearThread::postElement(const Elem * /*elem*/)
330 : {
331 288822237 : accumulate();
332 288822237 : }
333 :
334 : void
335 3412223 : NonlinearThread::post()
336 : {
337 3412223 : clearVarsAndMaterials();
338 3412223 : }
339 :
340 : void
341 3411861 : NonlinearThread::printGeneralExecutionInformation() const
342 : {
343 3411861 : if (_fe_problem.shouldPrintExecution(_tid))
344 : {
345 10935 : const auto & console = _fe_problem.console();
346 10935 : const auto execute_on = _fe_problem.getCurrentExecuteOnFlag();
347 21870 : console << "[DBG] Beginning elemental loop to compute " + objectType() + " on " << execute_on
348 10935 : << std::endl;
349 10935 : 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 10935 : }
358 3411861 : }
359 :
360 : void
361 4215108 : NonlinearThread::printBlockExecutionInformation() const
362 : {
363 : // Number of objects executing is approximated by size of warehouses
364 4215108 : const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
365 4215108 : _dg_kernels.size() + _interface_kernels.size();
366 4215108 : const auto & console = _fe_problem.console();
367 4215108 : const auto block_name = _mesh.getSubdomainName(_subdomain);
368 :
369 4215108 : if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
370 : {
371 68597 : if (_blocks_exec_printed.count(_subdomain))
372 52051 : return;
373 33092 : console << "[DBG] Ordering of " + objectType() + " Objects on block " << block_name << " ("
374 16546 : << _subdomain << ")" << std::endl;
375 16546 : if (_kernels.hasActiveBlockObjects(_subdomain, _tid))
376 : {
377 16482 : console << "[DBG] Ordering of kernels:" << std::endl;
378 49446 : console << _kernels.activeObjectsToFormattedString() << std::endl;
379 : }
380 16546 : 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 16546 : if (_dg_kernels.hasActiveBlockObjects(_subdomain, _tid))
392 : {
393 5893 : console << "[DBG] Ordering of DG kernels:" << std::endl;
394 17679 : console << _dg_kernels.activeObjectsToFormattedString() << std::endl;
395 : }
396 : }
397 4146511 : 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 4163057 : _blocks_exec_printed.insert(_subdomain);
403 4215108 : }
404 :
405 : void
406 105114596 : NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
407 : {
408 105218475 : if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
409 103879 : (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
410 82018 : !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
411 105090283 : return;
412 :
413 24313 : const auto & console = _fe_problem.console();
414 24313 : const auto b_name = _mesh.getBoundaryName(bid);
415 48626 : console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
416 24313 : << ")" << std::endl;
417 :
418 24313 : if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
419 : {
420 21861 : console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
421 65583 : 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 24313 : if (_interface_kernels.hasActiveBoundaryObjects(bid, _tid))
428 : {
429 2452 : console << "[DBG] Ordering of interface kernels:" << std::endl;
430 7356 : console << _interface_kernels.activeObjectsToFormattedString() << std::endl;
431 : }
432 :
433 24313 : _boundaries_exec_printed.insert(bid);
434 24313 : }
435 :
436 : void
437 4133985 : 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 4133985 : _fe_problem.reinitElemFace(elem, side, _tid);
443 :
444 : // Needed to use lower-dimensional variables on Materials
445 4133985 : if (lower_d_elem)
446 22916 : _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
447 :
448 4133985 : if (bnd_id != Moose::INVALID_BOUNDARY_ID)
449 : {
450 2868825 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
451 2868825 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
452 : }
453 : // Currently only used by HDG
454 : else
455 1265160 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
456 4133985 : }
457 :
458 : bool
459 1274512071 : 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 1274512071 : _should_execute_dg = false;
464 1274512071 : if (_subdomain_has_dg)
465 4111745 : _should_execute_dg =
466 4111745 : ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
467 1274512071 : return _subdomain_has_hdg || _should_execute_dg;
468 : }
|
