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 3867237 : NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
27 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
28 7734474 : _nl(fe_problem.currentNonlinearSystem()),
29 3867237 : _num_cached(0),
30 3867237 : _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31 3867237 : _dg_kernels(_nl.getDGKernelWarehouse()),
32 3867237 : _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33 3867237 : _kernels(_nl.getKernelWarehouse()),
34 3867237 : _hdg_kernels(_nl.getHDGKernelWarehouse()),
35 7523475 : _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
36 7796867 : _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
37 273392 : _fe_problem.haveFV()),
38 7734474 : _should_execute_dg(false)
39 : {
40 3867237 : }
41 :
42 : // Splitting Constructor
43 339925 : NonlinearThread::NonlinearThread(NonlinearThread & x, Threads::split split)
44 : : ThreadedElementLoop<ConstElemRange>(x, split),
45 339925 : _nl(x._nl),
46 339925 : _num_cached(x._num_cached),
47 339925 : _integrated_bcs(x._integrated_bcs),
48 339925 : _dg_kernels(x._dg_kernels),
49 339925 : _interface_kernels(x._interface_kernels),
50 339925 : _kernels(x._kernels),
51 339925 : _tag_kernels(x._tag_kernels),
52 339925 : _hdg_kernels(x._hdg_kernels),
53 339925 : _has_active_objects(x._has_active_objects),
54 339925 : _should_execute_dg(x._should_execute_dg)
55 : {
56 339925 : }
57 :
58 4207119 : NonlinearThread::~NonlinearThread() {}
59 :
60 : void
61 4206501 : NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
62 : {
63 4206501 : if (_has_active_objects)
64 4106276 : ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
65 4206467 : }
66 :
67 : void
68 5012461 : 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 5012461 : determineObjectWarehouses();
73 :
74 5012461 : _fe_problem.subdomainSetup(_subdomain, _tid);
75 :
76 : // Update variable Dependencies
77 5012461 : std::set<MooseVariableFEBase *> needed_moose_vars;
78 5012461 : _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
79 5012461 : _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
80 5012461 : _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
81 5012461 : _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
82 :
83 : // Update FE variable coupleable vector tags
84 5012461 : std::set<TagID> needed_fe_var_vector_tags;
85 5012461 : _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
86 5012461 : _subdomain, needed_fe_var_vector_tags, _tid);
87 5012461 : _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
88 5012461 : _subdomain, needed_fe_var_vector_tags, _tid);
89 5012461 : _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
90 5012461 : _subdomain, needed_fe_var_vector_tags, _tid);
91 :
92 : // Update material dependencies
93 5012461 : std::unordered_set<unsigned int> needed_mat_props;
94 5012461 : _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
95 5012461 : _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
96 5012461 : _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
97 5012461 : _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
98 :
99 5012461 : if (_fe_problem.haveFV())
100 : {
101 : // Re-query the finite volume elemental kernels
102 250001 : _fv_kernels.clear();
103 250001 : _fe_problem.theWarehouse()
104 250001 : .query()
105 500002 : .template condition<AttribSysNum>(_nl.number())
106 250001 : .template condition<AttribSystem>("FVElementalKernel")
107 250001 : .template condition<AttribSubdomains>(_subdomain)
108 250001 : .template condition<AttribThread>(_tid)
109 250001 : .queryInto(_fv_kernels);
110 542004 : for (const auto fv_kernel : _fv_kernels)
111 : {
112 292003 : const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
113 292003 : needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
114 292003 : const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
115 292003 : 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 5012461 : _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
121 5012461 : _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
122 :
123 5012461 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
124 5012461 : _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
125 5012461 : _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
126 5012461 : }
127 :
128 : void
129 384803259 : 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 384803259 : SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
134 :
135 384803259 : prepareElement(elem);
136 :
137 384803170 : if (dynamic_cast<ComputeJacobianThread *>(this))
138 55438562 : if (_nl.getScalarVariables(_tid).size() > 0)
139 177778 : _fe_problem.reinitOffDiagScalars(_tid);
140 :
141 384803170 : computeOnElement();
142 384803229 : }
143 :
144 : void
145 366795491 : NonlinearThread::computeOnElement()
146 : {
147 366795491 : if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
148 : {
149 353153441 : const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
150 1170294630 : for (const auto & kernel : kernels)
151 817141268 : compute(*kernel);
152 : }
153 :
154 366795412 : if (_fe_problem.haveFV())
155 28601554 : for (auto kernel : _fv_kernels)
156 15152371 : compute(*kernel);
157 366795412 : }
158 :
159 : void
160 117200596 : 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 117200596 : 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 3239737 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
172 :
173 3239737 : prepareFace(elem, side, bnd_id, lower_d_elem);
174 3239737 : computeOnBoundary(bnd_id, lower_d_elem);
175 :
176 3239733 : if (lower_d_elem)
177 : {
178 25484 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179 25484 : accumulateLower();
180 25484 : }
181 3239733 : }
182 117200592 : }
183 :
184 : void
185 3019352 : NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
186 : {
187 3019352 : const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
188 6349263 : for (const auto & bc : bcs)
189 3329915 : if (bc->shouldApply())
190 3329717 : compute(*bc);
191 3019348 : }
192 :
193 : void
194 478822 : NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
195 : {
196 478822 : if (_ik_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
197 : {
198 :
199 : // Pointer to the neighbor we are currently working on.
200 54438 : const Elem * neighbor = elem->neighbor_ptr(side);
201 :
202 54438 : if (neighbor->active())
203 : {
204 54438 : _fe_problem.reinitNeighbor(elem, side, _tid);
205 :
206 : // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
207 : // still remember to swap back during stack unwinding. Note that face, boundary, and interface
208 : // all operate with the same MaterialData object
209 54438 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
210 54438 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
211 54438 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
212 :
213 54438 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
214 54438 : _fe_problem.reinitMaterialsNeighborOnBoundary(bnd_id, neighbor->subdomain_id(), _tid);
215 :
216 : // Has to happen after face and neighbor properties have been computed. Note that we don't use
217 : // a sentinel here because FEProblem::swapBackMaterialsFace is going to handle face materials,
218 : // boundary materials, and interface materials (e.g. it queries the boundary material data
219 : // with the current element and side
220 54438 : _fe_problem.reinitMaterialsInterface(bnd_id, _tid);
221 :
222 54438 : computeOnInterface(bnd_id);
223 :
224 : {
225 54438 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
226 54438 : accumulateNeighbor();
227 54438 : }
228 54438 : }
229 : }
230 478822 : }
231 :
232 : void
233 46962 : NonlinearThread::computeOnInterface(BoundaryID bnd_id)
234 : {
235 46962 : const auto & int_ks = _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
236 126877 : for (const auto & interface_kernel : int_ks)
237 79915 : compute(*interface_kernel);
238 46962 : }
239 :
240 : void
241 6575889 : NonlinearThread::onInternalSide(const Elem * elem, unsigned int side)
242 : {
243 6575889 : if (_should_execute_dg)
244 : {
245 : // Pointer to the neighbor we are currently working on.
246 2419643 : const Elem * neighbor = elem->neighbor_ptr(side);
247 :
248 2419643 : _fe_problem.reinitElemNeighborAndLowerD(elem, side, _tid);
249 :
250 : // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
251 : // still remember to swap back during stack unwinding.
252 2419643 : SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
253 2419643 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
254 :
255 2419643 : SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
256 2419643 : _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
257 :
258 2419643 : computeOnInternalFace(neighbor);
259 :
260 : {
261 2419643 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
262 2419643 : accumulateNeighborLower();
263 2419643 : }
264 2419643 : }
265 6575889 : if (_subdomain_has_hdg)
266 : {
267 : // Set up Sentinel class so that, after we swap in reinitMaterialsFace in prepareFace, even if
268 : // one of our callees throws we remember to swap back during stack unwinding. We put our
269 : // sentinel here as opposed to in prepareFace because we certainly don't want our materials
270 : // swapped back before we proceed to residual/Jacobian computation
271 4156246 : SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
272 :
273 4156246 : prepareFace(elem, side, Moose::INVALID_BOUNDARY_ID, nullptr);
274 4156246 : computeOnInternalFace();
275 4156246 : }
276 6575889 : }
277 :
278 : void
279 2369063 : NonlinearThread::computeOnInternalFace(const Elem * neighbor)
280 : {
281 2369063 : const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
282 5583260 : for (const auto & dg_kernel : dgks)
283 3214197 : if (dg_kernel->hasBlocks(neighbor->subdomain_id()))
284 3211369 : compute(*dg_kernel, neighbor);
285 2369063 : }
286 :
287 : void
288 729565331 : NonlinearThread::compute(KernelBase & kernel)
289 : {
290 729565331 : compute(static_cast<ResidualObject &>(kernel));
291 729565277 : }
292 :
293 : void
294 15152371 : NonlinearThread::compute(FVElementalKernel & kernel)
295 : {
296 15152371 : compute(static_cast<ResidualObject &>(kernel));
297 15152371 : }
298 :
299 : void
300 3192621 : NonlinearThread::compute(IntegratedBCBase & bc)
301 : {
302 3192621 : compute(static_cast<ResidualObject &>(bc));
303 3192621 : }
304 :
305 : void
306 3107731 : NonlinearThread::compute(DGKernelBase & dg, const Elem * /*neighbor*/)
307 : {
308 3107731 : compute(static_cast<ResidualObject &>(dg));
309 3107731 : }
310 :
311 : void
312 76885 : NonlinearThread::compute(InterfaceKernelBase & ik)
313 : {
314 76885 : compute(static_cast<ResidualObject &>(ik));
315 76885 : }
316 :
317 : void
318 329364554 : NonlinearThread::postElement(const Elem * /*elem*/)
319 : {
320 329364554 : accumulate();
321 329364554 : }
322 :
323 : void
324 4106765 : NonlinearThread::post()
325 : {
326 4106765 : clearVarsAndMaterials();
327 4106765 : }
328 :
329 : void
330 4106276 : NonlinearThread::printGeneralExecutionInformation() const
331 : {
332 4106276 : if (_fe_problem.shouldPrintExecution(_tid))
333 : {
334 9373 : const auto & console = _fe_problem.console();
335 9373 : const auto execute_on = _fe_problem.getCurrentExecuteOnFlag();
336 18746 : console << "[DBG] Beginning elemental loop to compute " + objectType() + " on " << execute_on
337 9373 : << std::endl;
338 9373 : mooseDoOnce(
339 : console << "[DBG] Execution order on each element:" << std::endl;
340 : console << "[DBG] - kernels on element quadrature points" << std::endl;
341 : console << "[DBG] - finite volume elemental kernels on element" << std::endl;
342 : console << "[DBG] - integrated boundary conditions on element side quadrature points"
343 : << std::endl;
344 : console << "[DBG] - DG kernels on element side quadrature points" << std::endl;
345 : console << "[DBG] - interface kernels on element side quadrature points" << std::endl;);
346 9373 : }
347 4106276 : }
348 :
349 : void
350 5011788 : NonlinearThread::printBlockExecutionInformation() const
351 : {
352 : // Number of objects executing is approximated by size of warehouses
353 5011788 : const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
354 5011788 : _dg_kernels.size() + _interface_kernels.size();
355 5011788 : const auto & console = _fe_problem.console();
356 5011788 : const auto block_name = _mesh.getSubdomainName(_subdomain);
357 :
358 5011788 : if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
359 : {
360 78969 : if (_blocks_exec_printed.count(_subdomain))
361 62472 : return;
362 32994 : console << "[DBG] Ordering of " + objectType() + " Objects on block " << block_name << " ("
363 16497 : << _subdomain << ")" << std::endl;
364 16497 : if (_kernels.hasActiveBlockObjects(_subdomain, _tid))
365 : {
366 16385 : console << "[DBG] Ordering of kernels:" << std::endl;
367 49155 : console << _kernels.activeObjectsToFormattedString() << std::endl;
368 : }
369 16497 : if (_fv_kernels.size())
370 : {
371 112 : console << "[DBG] Ordering of FV elemental kernels:" << std::endl;
372 : std::string fvkernels =
373 112 : std::accumulate(_fv_kernels.begin() + 1,
374 : _fv_kernels.end(),
375 112 : _fv_kernels[0]->name(),
376 84 : [](const std::string & str_out, FVElementalKernel * kernel)
377 308 : { return str_out + " " + kernel->name(); });
378 112 : console << ConsoleUtils::formatString(fvkernels, "[DBG]") << std::endl;
379 112 : }
380 16497 : if (_dg_kernels.hasActiveBlockObjects(_subdomain, _tid))
381 : {
382 7424 : console << "[DBG] Ordering of DG kernels:" << std::endl;
383 22272 : console << _dg_kernels.activeObjectsToFormattedString() << std::endl;
384 : }
385 : }
386 4932819 : else if (_fe_problem.shouldPrintExecution(_tid) && num_objects == 0 &&
387 0 : !_blocks_exec_printed.count(_subdomain))
388 0 : console << "[DBG] No Active " + objectType() + " Objects on block " << block_name << " ("
389 0 : << _subdomain << ")" << std::endl;
390 :
391 4949316 : _blocks_exec_printed.insert(_subdomain);
392 5011788 : }
393 :
394 : void
395 117200596 : NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
396 : {
397 117303497 : if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
398 102901 : (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
399 80819 : !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
400 117175298 : return;
401 :
402 25298 : const auto & console = _fe_problem.console();
403 25298 : const auto b_name = _mesh.getBoundaryName(bid);
404 50596 : console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
405 25298 : << ")" << std::endl;
406 :
407 25298 : if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
408 : {
409 22082 : console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
410 66246 : console << _integrated_bcs.activeObjectsToFormattedString() << std::endl;
411 : }
412 :
413 : // We have not checked if we have a neighbor. This could be premature for saying we are executing
414 : // interface kernels. However, we should assume the execution will happen on another side of the
415 : // same boundary
416 25298 : if (_interface_kernels.hasActiveBoundaryObjects(bid, _tid))
417 : {
418 3216 : console << "[DBG] Ordering of interface kernels:" << std::endl;
419 9648 : console << _interface_kernels.activeObjectsToFormattedString() << std::endl;
420 : }
421 :
422 25298 : _boundaries_exec_printed.insert(bid);
423 25298 : }
424 :
425 : void
426 7395983 : NonlinearThread::prepareFace(const Elem * const elem,
427 : const unsigned int side,
428 : const BoundaryID bnd_id,
429 : const Elem * const lower_d_elem)
430 : {
431 7395983 : _fe_problem.reinitElemFace(elem, side, _tid);
432 :
433 : // Needed to use lower-dimensional variables on Materials
434 7395983 : if (lower_d_elem)
435 25484 : _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
436 :
437 7395983 : if (bnd_id != Moose::INVALID_BOUNDARY_ID)
438 : {
439 3239737 : _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
440 3239737 : _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
441 : }
442 : // Currently only used by HDG
443 : else
444 4156246 : _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
445 7395983 : }
446 :
447 : bool
448 1454219041 : NonlinearThread::shouldComputeInternalSide(const Elem & elem, const Elem & neighbor) const
449 : {
450 : // ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide gets expensive on high
451 : // h-refinement cases so we avoid it if possible
452 1454219041 : _should_execute_dg = false;
453 1454219041 : if (_subdomain_has_dg)
454 4836691 : _should_execute_dg =
455 4836691 : ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
456 1454219041 : return _subdomain_has_hdg || _should_execute_dg;
457 : }
|
