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 3519717 : NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
27 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
28 7039434 : _nl(fe_problem.currentNonlinearSystem()),
29 3519717 : _num_cached(0),
30 3519717 : _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31 3519717 : _dg_kernels(_nl.getDGKernelWarehouse()),
32 3519717 : _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33 3519717 : _kernels(_nl.getKernelWarehouse()),
34 3519717 : _hdg_kernels(_nl.getHDGKernelWarehouse()),
35 6845380 : _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
36 7052713 : _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
37 207333 : _fe_problem.haveFV()),
38 7039434 : _should_execute_dg(false)
39 : {
40 3519717 : }
41 :
42 : // Splitting Constructor
43 350072 : NonlinearThread::NonlinearThread(NonlinearThread & x, Threads::split split)
44 : : ThreadedElementLoop<ConstElemRange>(x, split),
45 350072 : _nl(x._nl),
46 350072 : _num_cached(x._num_cached),
47 350072 : _integrated_bcs(x._integrated_bcs),
48 350072 : _dg_kernels(x._dg_kernels),
49 350072 : _interface_kernels(x._interface_kernels),
50 350072 : _kernels(x._kernels),
51 350072 : _tag_kernels(x._tag_kernels),
52 350072 : _hdg_kernels(x._hdg_kernels),
53 350072 : _has_active_objects(x._has_active_objects),
54 350072 : _should_execute_dg(x._should_execute_dg)
55 : {
56 350072 : }
57 :
58 3869758 : NonlinearThread::~NonlinearThread() {}
59 :
60 : void
61 3869245 : NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
62 : {
63 3869245 : if (_has_active_objects)
64 3734966 : ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
65 3869221 : }
66 :
67 : void
68 4530726 : 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 4530726 : determineObjectWarehouses();
73 :
74 4530726 : _fe_problem.subdomainSetup(_subdomain, _tid);
75 :
76 : // Update variable Dependencies
77 4530726 : std::set<MooseVariableFEBase *> needed_moose_vars;
78 4530726 : _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
79 4530726 : _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
80 4530726 : _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
81 4530726 : _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
82 :
83 : // Update FE variable coupleable vector tags
84 4530726 : std::set<TagID> needed_fe_var_vector_tags;
85 4530726 : _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
86 4530726 : _subdomain, needed_fe_var_vector_tags, _tid);
87 4530726 : _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
88 4530726 : _subdomain, needed_fe_var_vector_tags, _tid);
89 4530726 : _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
90 4530726 : _subdomain, needed_fe_var_vector_tags, _tid);
91 :
92 : // Update material dependencies
93 4530726 : std::unordered_set<unsigned int> needed_mat_props;
94 4530726 : _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
95 4530726 : _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
96 4530726 : _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
97 4530726 : _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
98 :
99 4530726 : if (_fe_problem.haveFV())
100 : {
101 : // Re-query the finite volume elemental kernels
102 148457 : _fv_kernels.clear();
103 148457 : _fe_problem.theWarehouse()
104 148457 : .query()
105 296914 : .template condition<AttribSysNum>(_nl.number())
106 148457 : .template condition<AttribSystem>("FVElementalKernel")
107 148457 : .template condition<AttribSubdomains>(_subdomain)
108 148457 : .template condition<AttribThread>(_tid)
109 148457 : .queryInto(_fv_kernels);
110 260601 : for (const auto fv_kernel : _fv_kernels)
111 : {
112 112144 : const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
113 112144 : needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
114 112144 : const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
115 112144 : 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 4530726 : _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
121 4530726 : _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
122 :
123 4530726 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
124 4530726 : _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
125 4530726 : _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
126 4530726 : }
127 :
128 : void
129 334353147 : 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 334353147 : SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
134 :
135 334353147 : prepareElement(elem);
136 :
137 334353081 : if (dynamic_cast<ComputeJacobianThread *>(this))
138 46593365 : if (_nl.getScalarVariables(_tid).size() > 0)
139 158338 : _fe_problem.reinitOffDiagScalars(_tid);
140 :
141 334353081 : computeOnElement();
142 334353125 : }
143 :
144 : void
145 320151743 : NonlinearThread::computeOnElement()
146 : {
147 320151743 : if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
148 : {
149 312088647 : const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
150 1033774934 : for (const auto & kernel : kernels)
151 721686355 : compute(*kernel);
152 : }
153 :
154 320151675 : if (_fe_problem.haveFV())
155 13946270 : for (auto kernel : _fv_kernels)
156 6088606 : compute(*kernel);
157 320151675 : }
158 :
159 : void
160 104681061 : 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 104681061 : 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 2802580 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
172 :
173 2802580 : prepareFace(elem, side, bnd_id, lower_d_elem);
174 2802580 : computeOnBoundary(bnd_id, lower_d_elem);
175 :
176 2802578 : if (lower_d_elem)
177 22916 : accumulateLower();
178 2802578 : }
179 104681059 : }
180 :
181 : void
182 2613612 : NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
183 : {
184 2613612 : const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
185 5467255 : for (const auto & bc : bcs)
186 2853645 : if (bc->shouldApply())
187 2853468 : compute(*bc);
188 2613610 : }
189 :
190 : void
191 422505 : NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
192 : {
193 422505 : 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 422505 : }
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 3319760 : NonlinearThread::onInternalSide(const Elem * elem, unsigned int side)
236 : {
237 3319760 : if (_should_execute_dg)
238 : {
239 : // Pointer to the neighbor we are currently working on.
240 2057704 : const Elem * neighbor = elem->neighbor_ptr(side);
241 :
242 2057704 : _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 2057704 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
247 2057704 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
248 :
249 2057704 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
250 2057704 : _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
251 :
252 2057704 : computeOnInternalFace(neighbor);
253 :
254 2057704 : accumulateNeighborLower();
255 2057704 : }
256 3319760 : 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 1262056 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
263 :
264 1262056 : prepareFace(elem, side, Moose::INVALID_BOUNDARY_ID, nullptr);
265 1262056 : computeOnInternalFace();
266 1262056 : }
267 3319760 : }
268 :
269 : void
270 106057440 : NonlinearThread::onExternalSide(const Elem * elem, unsigned int side)
271 : {
272 : // Check that we don't have any interface kernels defined on this boundary
273 106057440 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
274 :
275 106057440 : bool has_interface_kernels = false;
276 210315994 : for (const auto bid : boundary_ids)
277 104258554 : if (_ik_warehouse && _ik_warehouse->hasActiveBoundaryObjects(bid, _tid))
278 0 : has_interface_kernels = true;
279 :
280 106057440 : 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 106057440 : }
288 :
289 : void
290 2012984 : NonlinearThread::computeOnInternalFace(const Elem * neighbor)
291 : {
292 2012984 : const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
293 4696077 : for (const auto & dg_kernel : dgks)
294 2683093 : if (dg_kernel->hasBlocks(neighbor->subdomain_id()))
295 2680579 : compute(*dg_kernel, neighbor);
296 2012984 : }
297 :
298 : void
299 645601546 : NonlinearThread::compute(KernelBase & kernel)
300 : {
301 645601546 : compute(static_cast<ResidualObject &>(kernel));
302 645601499 : }
303 :
304 : void
305 6088406 : NonlinearThread::compute(FVElementalKernel & kernel)
306 : {
307 6088406 : compute(static_cast<ResidualObject &>(kernel));
308 6088406 : }
309 :
310 : void
311 2734658 : NonlinearThread::compute(IntegratedBCBase & bc)
312 : {
313 2734658 : compute(static_cast<ResidualObject &>(bc));
314 2734658 : }
315 :
316 : void
317 2595971 : NonlinearThread::compute(DGKernelBase & dg, const Elem * /*neighbor*/)
318 : {
319 2595971 : compute(static_cast<ResidualObject &>(dg));
320 2595971 : }
321 :
322 : void
323 60227 : NonlinearThread::compute(InterfaceKernelBase & ik)
324 : {
325 60227 : compute(static_cast<ResidualObject &>(ik));
326 60227 : }
327 :
328 : void
329 287759669 : NonlinearThread::postElement(const Elem * /*elem*/)
330 : {
331 287759669 : accumulate();
332 287759669 : }
333 :
334 : void
335 3735374 : NonlinearThread::post()
336 : {
337 3735374 : clearVarsAndMaterials();
338 3735374 : }
339 :
340 : void
341 3734966 : NonlinearThread::printGeneralExecutionInformation() const
342 : {
343 3734966 : 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 3734966 : }
359 :
360 : void
361 4530115 : NonlinearThread::printBlockExecutionInformation() const
362 : {
363 : // Number of objects executing is approximated by size of warehouses
364 4530115 : const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
365 4530115 : _dg_kernels.size() + _interface_kernels.size();
366 4530115 : const auto & console = _fe_problem.console();
367 4530115 : const auto block_name = _mesh.getSubdomainName(_subdomain);
368 :
369 4530115 : if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
370 : {
371 62693 : if (_blocks_exec_printed.count(_subdomain))
372 46147 : 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 4467422 : 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 4483968 : _blocks_exec_printed.insert(_subdomain);
403 4530115 : }
404 :
405 : void
406 104681061 : NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
407 : {
408 104778313 : if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
409 97252 : (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
410 76437 : !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
411 104657794 : return;
412 :
413 23267 : const auto & console = _fe_problem.console();
414 23267 : const auto b_name = _mesh.getBoundaryName(bid);
415 46534 : console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
416 23267 : << ")" << std::endl;
417 :
418 23267 : if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
419 : {
420 20815 : console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
421 62445 : 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 23267 : 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 23267 : _boundaries_exec_printed.insert(bid);
434 23267 : }
435 :
436 : void
437 4064636 : 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 4064636 : _fe_problem.reinitElemFace(elem, side, _tid);
443 :
444 : // Needed to use lower-dimensional variables on Materials
445 4064636 : if (lower_d_elem)
446 22916 : _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
447 :
448 4064636 : if (bnd_id != Moose::INVALID_BOUNDARY_ID)
449 : {
450 2802580 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
451 2802580 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
452 : }
453 : // Currently only used by HDG
454 : else
455 1262056 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
456 4064636 : }
457 :
458 : bool
459 1270513601 : 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 1270513601 : _should_execute_dg = false;
464 1270513601 : if (_subdomain_has_dg)
465 4113033 : _should_execute_dg =
466 4113033 : ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
467 1270513601 : return _subdomain_has_hdg || _should_execute_dg;
468 : }
|
