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