www.mooseframework.org
ComputeJacobianThread.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "ComputeJacobianThread.h"
11 
12 #include "DGKernelBase.h"
13 #include "FEProblem.h"
14 #include "IntegratedBCBase.h"
15 #include "InterfaceKernel.h"
16 #include "MooseVariableFE.h"
17 #include "NonlinearSystem.h"
18 #include "NonlocalIntegratedBC.h"
19 #include "NonlocalKernel.h"
20 #include "SwapBackSentinel.h"
21 #include "TimeDerivative.h"
22 
23 #include "libmesh/threads.h"
24 
26  const std::set<TagID> & tags)
27  : ThreadedElementLoop<ConstElemRange>(fe_problem),
28  _nl(fe_problem.getNonlinearSystemBase()),
29  _num_cached(0),
30  _integrated_bcs(_nl.getIntegratedBCWarehouse()),
31  _dg_kernels(_nl.getDGKernelWarehouse()),
32  _interface_kernels(_nl.getInterfaceKernelWarehouse()),
33  _kernels(_nl.getKernelWarehouse()),
34  _ad_jacobian_kernels(_nl.getADJacobianKernelWarehouse()),
35  _tags(tags)
36 {
37 }
38 
39 // Splitting Constructor
41  : ThreadedElementLoop<ConstElemRange>(x, split),
42  _nl(x._nl),
43  _num_cached(x._num_cached),
44  _integrated_bcs(x._integrated_bcs),
45  _dg_kernels(x._dg_kernels),
46  _interface_kernels(x._interface_kernels),
47  _kernels(x._kernels),
48  _warehouse(x._warehouse),
49  _ad_jacobian_kernels(x._ad_jacobian_kernels),
50  _tags(x._tags)
51 {
52 }
53 
55 
56 void
58 {
60  {
61  auto & kernels = _warehouse->getActiveBlockObjects(_subdomain, _tid);
62  for (const auto & kernel : kernels)
63  if (kernel->isImplicit())
64  {
65  kernel->subProblem().prepareShapes(kernel->variable().number(), _tid);
66  kernel->computeJacobian();
69  {
70  std::shared_ptr<NonlocalKernel> nonlocal_kernel =
72  if (nonlocal_kernel)
73  kernel->computeNonlocalJacobian();
74  }
75  }
76  }
78  {
80  for (const auto & kernel : kernels)
81  if (kernel->isImplicit())
82  kernel->computeJacobian();
83  }
84 }
85 
86 void
88 {
89  const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
90  for (const auto & bc : bcs)
91  if (bc->shouldApply() && bc->isImplicit())
92  {
93  bc->subProblem().prepareFaceShapes(bc->variable().number(), _tid);
94  bc->computeJacobian();
97  {
98  std::shared_ptr<NonlocalIntegratedBC> nonlocal_integrated_bc =
100  if (nonlocal_integrated_bc)
102  }
103  }
104 }
105 
106 void
108 {
109  // No need to call hasActiveObjects, this is done in the calling method (see onInternalSide)
110  const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
111  for (const auto & dg : dgks)
112  if (dg->isImplicit())
113  {
114  dg->subProblem().prepareFaceShapes(dg->variable().number(), _tid);
115  dg->subProblem().prepareNeighborShapes(dg->variable().number(), _tid);
116  if (dg->hasBlocks(neighbor->subdomain_id()))
117  dg->computeJacobian();
118  }
119 }
120 
121 void
123 {
124  // No need to call hasActiveObjects, this is done in the calling method (see onInterface)
125  const auto & intks = _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
126  for (const auto & intk : intks)
127  if (intk->isImplicit())
128  {
129  intk->subProblem().prepareFaceShapes(intk->variable().number(), _tid);
130  intk->subProblem().prepareNeighborShapes(intk->neighborVariable().number(), _tid);
131  intk->computeJacobian();
132  }
133 }
134 
135 void
137 {
139 
140  // Update variable Dependencies
141  std::set<MooseVariableFEBase *> needed_moose_vars;
147 
148  // Update material dependencies
149  std::set<unsigned int> needed_mat_props;
155 
157  _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
159 
160  // If users pass a empty vector or a full size of vector,
161  // we take all kernels
162  if (!_tags.size() || _tags.size() == _fe_problem.numMatrixTags())
163  {
164  _warehouse = &_kernels;
169  }
170  // If we have one tag only,
171  // We call tag based storage
172  else if (_tags.size() == 1)
173  {
179  }
180  // This one may be expensive, and hopefully we do not use it so often
181  else
182  {
188  }
189 }
190 
191 void
193 {
194  _fe_problem.prepare(elem, _tid);
195 
196  _fe_problem.reinitElem(elem, _tid);
197 
198  // Set up Sentinel class so that, even if reinitMaterials() throws, we
199  // still remember to swap back during stack unwinding.
202 
203  if (_nl.getScalarVariables(_tid).size() > 0)
205 
206  computeJacobian();
207 }
208 
209 void
210 ComputeJacobianThread::onBoundary(const Elem * elem, unsigned int side, BoundaryID bnd_id)
211 {
213  {
214  _fe_problem.reinitElemFace(elem, side, bnd_id, _tid);
215 
216  // Set up Sentinel class so that, even if reinitMaterials() throws, we
217  // still remember to swap back during stack unwinding.
219 
220  _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
222 
223  computeFaceJacobian(bnd_id);
224  }
225 }
226 
227 void
228 ComputeJacobianThread::onInternalSide(const Elem * elem, unsigned int side)
229 {
231  {
232  // Pointer to the neighbor we are currently working on.
233  const Elem * neighbor = elem->neighbor_ptr(side);
234 
235  // Get the global id of the element and the neighbor
236  const dof_id_type elem_id = elem->id(), neighbor_id = neighbor->id();
237 
238  if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) ||
239  (neighbor->level() < elem->level()))
240  {
241  _fe_problem.reinitNeighbor(elem, side, _tid);
242 
243  // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
244  // still remember to swap back during stack unwinding.
246  _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
247 
249  _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
250 
251  computeInternalFaceJacobian(neighbor);
252 
253  {
254  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
256  }
257  }
258  }
259 }
260 
261 void
262 ComputeJacobianThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
263 {
265  {
266  // Pointer to the neighbor we are currently working on.
267  const Elem * neighbor = elem->neighbor_ptr(side);
268 
269  if (neighbor->active())
270  {
271  _fe_problem.reinitNeighbor(elem, side, _tid);
272 
273  // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
274  // still remember to swap back during stack unwinding.
276  _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
278 
280  _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
281 
283 
284  {
285  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
287  }
288  }
289  }
290 }
291 
292 void
294 {
296  _num_cached++;
297 
298  if (_num_cached % 20 == 0)
299  {
300  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
302  }
303 }
304 
305 void
307 {
310 }
311 
312 void
314 {
315 }
virtual void post() override
Called after the element range loop.
void updateBoundaryVariableDependency(std::set< MooseVariableFEBase *> &needed_moose_vars, THREAD_ID tid=0) const
virtual void postElement(const Elem *) override
Called after the element assembly is done (including surface assembling)
virtual void prepare(const Elem *elem, THREAD_ID tid) override
Base class for assembly-like calculations.
const std::vector< MooseVariableScalar * > & getScalarVariables(THREAD_ID tid)
Definition: SystemBase.h:614
bool hasActiveBlockObjects(THREAD_ID tid=0) const
virtual void onElement(const Elem *elem) override
Assembly of the element (not including surface assembly)
void updateBlockMatPropDependency(SubdomainID id, std::set< unsigned int > &needed_mat_props, THREAD_ID tid=0) const
virtual void prepareMaterials(SubdomainID blk_id, THREAD_ID tid)
Add the MooseVariables that the current materials depend on to the dependency list.
virtual void setActiveElementalMooseVariables(const std::set< MooseVariableFEBase *> &moose_vars, THREAD_ID tid) override
Set the MOOSE variables to be reinited on each element.
const std::map< SubdomainID, std::vector< std::shared_ptr< T > > > & getActiveBlockObjects(THREAD_ID tid=0) const
MooseObjectTagWarehouse< KernelBase > & _kernels
NonlocalIntegratedBC is used for solving integral terms in integro-differential equations.
virtual void reinitMaterialsNeighbor(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
void updateBlockVariableDependency(SubdomainID id, std::set< MooseVariableFEBase *> &needed_moose_vars, THREAD_ID tid=0) const
virtual void addJacobianNeighbor(THREAD_ID tid) override
MooseObjectTagWarehouse< IntegratedBCBase > & _integrated_bcs
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
static PetscErrorCode Vec x
MooseObjectWarehouse< DGKernelBase > * _dg_warehouse
virtual void reinitElem(const Elem *elem, THREAD_ID tid) override
virtual void reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
virtual void reinitMaterialsBoundary(BoundaryID boundary_id, THREAD_ID tid, bool swap_stateful=true)
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
virtual void setActiveMaterialProperties(const std::set< unsigned int > &mat_prop_ids, THREAD_ID tid) override
Record and set the material properties required by the current computing thread.
virtual void subdomainSetup(SubdomainID subdomain, THREAD_ID tid)
MooseObjectWarehouse< InterfaceKernel > * _ik_warehouse
const std::set< TagID > & _tags
bool hasActiveBoundaryObjects(THREAD_ID tid=0) const
std::vector< std::string > split(const std::string &str, const std::string &delimiter)
Python like split function for strings.
Definition: MooseUtils.C:784
virtual void computeNonlocalJacobian() override
computeNonlocalJacobian and computeNonlocalOffDiagJacobian methods are introduced for providing the j...
virtual void addCachedJacobian(THREAD_ID tid) override
MooseObjectTagWarehouse< KernelBase > & _ad_jacobian_kernels
Reference to ADKernel<JACOBIAN> storage structure.
MooseObjectWarehouse< KernelBase > * _adjk_warehouse
Pointer to tags.
virtual void cacheJacobian(THREAD_ID tid) override
virtual void computeInternalInterFaceJacobian(BoundaryID bnd_id)
boundary_id_type BoundaryID
virtual void clearActiveMaterialProperties(THREAD_ID tid) override
Clear the active material properties.
MooseObjectTagWarehouse< DGKernelBase > & _dg_kernels
virtual bool checkNonlocalCouplingRequirement()
Definition: SubProblem.h:69
MooseObjectWarehouse< T > & getMatrixTagObjectWarehouse(TagID tag_id, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has the given matrix tag...
MooseObjectWarehouse< KernelBase > * _warehouse
virtual void swapBackMaterials(THREAD_ID tid)
MooseObjectWarehouse< IntegratedBCBase > * _ibc_warehouse
NonlinearSystemBase & _nl
virtual void swapBackMaterialsNeighbor(THREAD_ID tid)
const std::map< BoundaryID, std::vector< std::shared_ptr< T > > > & getActiveBoundaryObjects(THREAD_ID tid=0) const
virtual void subdomainChanged() override
Called every time the current subdomain changes (i.e.
virtual unsigned int numMatrixTags() const
The total number of tags.
Definition: SubProblem.h:157
ComputeJacobianThread(FEProblemBase &fe_problem, const std::set< TagID > &tags)
virtual void reinitElemFace(const Elem *elem, unsigned int side, BoundaryID bnd_id, THREAD_ID tid) override
virtual void onBoundary(const Elem *elem, unsigned int side, BoundaryID bnd_id) override
Called when doing boundary assembling.
virtual void reinitMaterialsFace(SubdomainID blk_id, THREAD_ID tid, bool swap_stateful=true)
MooseObjectWarehouse< T > & getMatrixTagsObjectWarehouse(const std::set< TagID > &tags, THREAD_ID tid)
Retrieve a moose object warehouse in which every moose object has one of the given matrix tags...
virtual void computeFaceJacobian(BoundaryID bnd_id)
void join(const ComputeJacobianThread &)
NonlocalKernel is used for solving integral terms in integro-differential equations.
virtual void reinitNeighbor(const Elem *elem, unsigned int side, THREAD_ID tid) override
virtual void computeInternalFaceJacobian(const Elem *neighbor)
SubdomainID _subdomain
The subdomain for the current element.
virtual void reinitOffDiagScalars(THREAD_ID tid) override
virtual void computeNonlocalJacobian() override
computeNonlocalJacobian and computeNonlocalOffDiagJacobian methods are introduced for providing the j...
virtual void onInterface(const Elem *elem, unsigned int side, BoundaryID bnd_id) override
Called when doing interface assembling.
MooseObjectTagWarehouse< InterfaceKernel > & _interface_kernels
virtual void swapBackMaterialsFace(THREAD_ID tid)
virtual void onInternalSide(const Elem *elem, unsigned int side) override
Called when doing internal edge assembling.
virtual void clearActiveElementalMooseVariables(THREAD_ID tid) override
Clear the active elemental MooseVariableFEBase.
void updateBoundaryMatPropDependency(std::set< unsigned int > &needed_mat_props, THREAD_ID tid=0) const
The "SwapBackSentinel" class&#39;s destructor guarantees that FEProblemBase::swapBackMaterials{Face,Neighbor}() is called even when an exception is thrown from FEProblemBase::reinitMaterials{Face,Neighbor}.