www.mooseframework.org
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
AdvectiveFluxCalculatorBase Class Referenceabstract

Base class to compute Advective fluxes. More...

#include <AdvectiveFluxCalculatorBase.h>

Inheritance diagram for AdvectiveFluxCalculatorBase:
[legend]

Public Member Functions

 AdvectiveFluxCalculatorBase (const InputParameters &parameters)
 
virtual void timestepSetup () override
 
virtual void meshChanged () override
 
virtual void initialize () override
 
virtual void threadJoin (const UserObject &uo) override
 
virtual void finalize () override
 
virtual void execute () override
 
virtual void executeOnElement (dof_id_type global_i, dof_id_type global_j, unsigned local_i, unsigned local_j, unsigned qp)
 This is called by multiple times in execute() in a double loop over _current_elem's nodes (local_i and local_j) nested in a loop over each of _current_elem's quadpoints (qp). More...
 
Real getFluxOut (dof_id_type node_i) const
 Returns the flux out of lobal node id. More...
 
const std::map< dof_id_type, Real > & getdFluxOutdu (dof_id_type node_i) const
 Returns r where r[j] = d(flux out of global node i)/du(global node j) used in Jacobian computations. More...
 
const std::map< dof_id_type, std::map< dof_id_type, Real > > & getdFluxOutdKjk (dof_id_type node_i) const
 Returns r where r[j] = d(flux out of global node i)/dK[global node j][global node k] used in Jacobian computations. More...
 
unsigned getValence (dof_id_type node_i, dof_id_type node_j) const
 Returns the valence of the i-j edge. More...
 

Protected Types

enum  FluxLimiterTypeEnum {
  FluxLimiterTypeEnum::MinMod, FluxLimiterTypeEnum::VanLeer, FluxLimiterTypeEnum::MC, FluxLimiterTypeEnum::superbee,
  FluxLimiterTypeEnum::None
}
 Determines Flux Limiter type (Page 135 of Kuzmin and Turek) "None" means that limitFlux=0 always, which implies zero antidiffusion will be added. More...
 
enum  PQPlusMinusEnum { PQPlusMinusEnum::PPlus, PQPlusMinusEnum::PMinus, PQPlusMinusEnum::QPlus, PQPlusMinusEnum::QMinus }
 Signals to the PQPlusMinus method what should be computed. More...
 

Protected Member Functions

virtual Real computeVelocity (unsigned i, unsigned j, unsigned qp) const =0
 Computes the transfer velocity between current node i and current node j at the current qp in the current element (_current_elem). More...
 
virtual Real computeU (unsigned i) const =0
 Computes the value of u at the local node id of the current element (_current_elem) More...
 
void limitFlux (Real a, Real b, Real &limited, Real &dlimited_db) const
 flux limiter, L, on Page 135 of Kuzmin and Turek More...
 
Real rPlus (dof_id_type node_i, std::map< dof_id_type, Real > &dlimited_du, std::map< dof_id_type, Real > &dlimited_dk) const
 Returns the value of R_{i}^{+}, Eqn (49) of KT. More...
 
Real rMinus (dof_id_type node_i, std::map< dof_id_type, Real > &dlimited_du, std::map< dof_id_type, Real > &dlimited_dk) const
 Returns the value of R_{i}^{-}, Eqn (49) of KT. More...
 
Real getKij (dof_id_type node_i, dof_id_type node_j) const
 Returns the value of k_ij as computed by KT Eqns (18)-(20). More...
 
Real PQPlusMinus (dof_id_type node_i, const PQPlusMinusEnum pq_plus_minus, std::map< dof_id_type, Real > &derivs, std::map< dof_id_type, Real > &dpq_dk) const
 Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which are defined in Eqns (47) and (48) of KT. More...
 
void zeroedConnection (std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
 Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i. More...
 
Real getUnodal (dof_id_type node_i) const
 Returns the value of u at the global id node_i. More...
 
bool getUnodalComputedByThread (dof_id_type node_i) const
 Returns the value of _u_nodal_computed_by_thread at the global id node_i. More...
 

Protected Attributes

bool _resizing_needed
 whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timestep More...
 
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
 
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
 Kuzmin-Turek K_ij matrix. More...
 
std::map< dof_id_type, Real > _flux_out
 _flux_out[i] = flux of "heat" from node i More...
 
std::map< dof_id_type, std::map< dof_id_type, Real > > _dflux_out_du
 _dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is connected to i. More...
 
std::map< dof_id_type, std::map< dof_id_type, std::map< dof_id_type, Real > > > _dflux_out_dKjk
 _dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does include j = i), and k must be connected to j (this does include k = i and k = j) More...
 
std::map< std::pair< dof_id_type, dof_id_type >, unsigned > _valence
 _valence[(i, j)] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the i-j edge is encountered More...
 
std::map< dof_id_type, Real > _u_nodal
 _u_nodal[i] = value of _u at global node number i More...
 
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
 _u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on this processor More...
 

Detailed Description

Base class to compute Advective fluxes.

Specifically, computes K_ij, D_ij, L_ij, R+, R-, f^a_ij detailed in D Kuzmin and S Turek "High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter" Journal of Computational Physics 198 (2004) 131-158 Then combines the results to produce the residual and Jacobian contributions that may be used by Kernels

K_ij is a measure of flux from node i to node j D_ij is a diffusion matrix that stabilizes K_ij (K+D has the LED property) L = K + D R^+_i and R^-_i and f^a_{ij} quantify how much antidiffusion to allow around node i

Definition at line 33 of file AdvectiveFluxCalculatorBase.h.

Member Enumeration Documentation

◆ FluxLimiterTypeEnum

Determines Flux Limiter type (Page 135 of Kuzmin and Turek) "None" means that limitFlux=0 always, which implies zero antidiffusion will be added.

Enumerator
MinMod 
VanLeer 
MC 
superbee 
None 

Definition at line 162 of file AdvectiveFluxCalculatorBase.h.

162 { MinMod, VanLeer, MC, superbee, None } _flux_limiter_type;
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type

◆ PQPlusMinusEnum

Signals to the PQPlusMinus method what should be computed.

Enumerator
PPlus 
PMinus 
QPlus 
QMinus 

Definition at line 193 of file AdvectiveFluxCalculatorBase.h.

194  {
195  PPlus,
196  PMinus,
197  QPlus,
198  QMinus
199  };

Constructor & Destructor Documentation

◆ AdvectiveFluxCalculatorBase()

AdvectiveFluxCalculatorBase::AdvectiveFluxCalculatorBase ( const InputParameters &  parameters)

Definition at line 33 of file AdvectiveFluxCalculatorBase.C.

34  : ElementUserObject(parameters),
35  _resizing_needed(true),
36  _flux_limiter_type(getParam<MooseEnum>("flux_limiter_type").getEnum<FluxLimiterTypeEnum>()),
37  _kij({}),
38  _flux_out({}),
39  _dflux_out_du({}),
40  _dflux_out_dKjk({}),
41  _valence({}),
42  _u_nodal({}),
44 {
45  if (!_execute_enum.contains(EXEC_LINEAR))
46  paramError(
47  "execute_on",
48  "The AdvectiveFluxCalculator UserObject " + name() +
49  " execute_on parameter must include, at least, 'linear'. This is to ensure that "
50  "this UserObject computes all necessary quantities just before the Kernels evaluate "
51  "their Residuals");
52 }
std::map< dof_id_type, Real > _u_nodal
_u_nodal[i] = value of _u at global node number i
std::map< dof_id_type, std::map< dof_id_type, std::map< dof_id_type, Real > > > _dflux_out_dKjk
_dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does includ...
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...
std::map< std::pair< dof_id_type, dof_id_type >, unsigned > _valence
_valence[(i, j)] = number of times, in a loop over elements seen by this processor (viz...
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.
std::map< dof_id_type, Real > _flux_out
_flux_out[i] = flux of "heat" from node i
std::map< dof_id_type, std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is con...

Member Function Documentation

◆ computeU()

virtual Real AdvectiveFluxCalculatorBase::computeU ( unsigned  i) const
protectedpure virtual

◆ computeVelocity()

virtual Real AdvectiveFluxCalculatorBase::computeVelocity ( unsigned  i,
unsigned  j,
unsigned  qp 
) const
protectedpure virtual

Computes the transfer velocity between current node i and current node j at the current qp in the current element (_current_elem).

For instance, (_grad_phi[i][qp] * _velocity) * _phi[j][qp];

Parameters
inode number in the current element
jnode number in the current element
qpquadpoint number in the current element

Implemented in PorousFlowAdvectiveFluxCalculatorBase, and AdvectiveFluxCalculatorConstantVelocity.

Referenced by executeOnElement().

◆ execute()

void AdvectiveFluxCalculatorBase::execute ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 171 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::execute(), and finalize().

172 {
173  // compute _kij contributions from this element that is local to this processor
174  // and record _u_nodal
175  for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
176  {
177  const dof_id_type node_i = _current_elem->node_id(i);
178  if (!_u_nodal_computed_by_thread[node_i])
179  {
180  _u_nodal[node_i] = computeU(i);
181  _u_nodal_computed_by_thread[node_i] = true;
182  }
183  for (unsigned j = 0; j < _current_elem->n_nodes(); ++j)
184  {
185  const dof_id_type node_j = _current_elem->node_id(j);
186  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
187  executeOnElement(node_i, node_j, i, j, qp);
188  }
189  }
190 }
virtual void executeOnElement(dof_id_type global_i, dof_id_type global_j, unsigned local_i, unsigned local_j, unsigned qp)
This is called by multiple times in execute() in a double loop over _current_elem&#39;s nodes (local_i an...
std::map< dof_id_type, Real > _u_nodal
_u_nodal[i] = value of _u at global node number i
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...
virtual Real computeU(unsigned i) const =0
Computes the value of u at the local node id of the current element (_current_elem) ...

◆ executeOnElement()

void AdvectiveFluxCalculatorBase::executeOnElement ( dof_id_type  global_i,
dof_id_type  global_j,
unsigned  local_i,
unsigned  local_j,
unsigned  qp 
)
virtual

This is called by multiple times in execute() in a double loop over _current_elem's nodes (local_i and local_j) nested in a loop over each of _current_elem's quadpoints (qp).

It is used to compute _kij and its derivatives

Parameters
global_iglobal node id corresponding to the local node local_i
global_jglobal node id corresponding to the local node local_j
local_ilocal node number of the _current_elem
local_jlocal node number of the _current_elem
qpquadpoint number of the _current_elem

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 193 of file AdvectiveFluxCalculatorBase.C.

Referenced by execute(), and PorousFlowAdvectiveFluxCalculatorBase::executeOnElement().

195 {
196  // KT Eqn (20)
197  _kij[global_i][global_j] += _JxW[qp] * _coord[qp] * computeVelocity(local_i, local_j, qp);
198 }
virtual Real computeVelocity(unsigned i, unsigned j, unsigned qp) const =0
Computes the transfer velocity between current node i and current node j at the current qp in the cur...
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ finalize()

void AdvectiveFluxCalculatorBase::finalize ( )
overridevirtual

Definition at line 225 of file AdvectiveFluxCalculatorBase.C.

226 {
227  // Overall: ensure _kij is fully built, then compute Kuzmin-Turek D, L, f^a and
228  // relevant Jacobian information, and then the relevant quantities into _flux_out and
229  // _dflux_out_du, _dflux_out_dKjk
230 
231  /*
232  * MULTIPROC NOTE: execute() only computed contributions from the elements
233  * local to this processor. Now do the same for the 2 layers of ghosted elements
234  */
235  for (const auto & elem_id : _fe_problem.ghostedElems())
236  {
237  _current_elem = _mesh.elem(elem_id);
238  if (this->hasBlocks(_current_elem->subdomain_id()))
239  {
240  _fe_problem.prepare(_current_elem, _tid);
241  _fe_problem.reinitElem(_current_elem, _tid);
242  execute();
243  }
244  }
245 
246  // Calculate KuzminTurek D matrix
247  // See Eqn (32)
248  // This adds artificial diffusion, which eliminates any spurious oscillations
249  // The idea is that D will remove all negative off-diagonal elements when it is added to K
250  // This is identical to full upwinding
251  std::map<dof_id_type, std::map<dof_id_type, Real>> dij;
252  std::map<dof_id_type, std::map<dof_id_type, Real>>
253  dDij_dKij; // dDij_dKij[i][j] = d(D[i][j])/d(K[i][j]) for i!=j
254  std::map<dof_id_type, std::map<dof_id_type, Real>>
255  dDij_dKji; // dDij_dKji[i][j] = d(D[i][j])/d(K[j][i]) for i!=j
256  std::map<dof_id_type, std::map<dof_id_type, Real>>
257  dDii_dKij; // dDii_dKij[i][j] = d(D[i][i])/d(K[i][j])
258  std::map<dof_id_type, std::map<dof_id_type, Real>>
259  dDii_dKji; // dDii_dKji[i][j] = d(D[i][i])/d(K[j][i])
260  for (auto & nodes : _kij)
261  {
262  const dof_id_type i = nodes.first;
263  dij[i] = {};
264  dij[i][i] = 0.0;
265  dDij_dKij[i] = {};
266  zeroedConnection(dDij_dKij[i], i);
267  dDij_dKji[i] = {};
268  zeroedConnection(dDij_dKji[i], i);
269  dDii_dKij[i] = {};
270  zeroedConnection(dDii_dKij[i], i);
271  dDii_dKji[i] = {};
272  zeroedConnection(dDii_dKji[i], i);
273  for (auto & neighbors : nodes.second)
274  {
275  const dof_id_type j = neighbors.first;
276  if (i == j)
277  continue;
278  if ((_kij[i][j] <= _kij[j][i]) && (_kij[i][j] < 0))
279  {
280  dij[i][j] = -_kij[i][j];
281  dDij_dKij[i][j] = -1.0;
282  dDii_dKij[i][j] += 1.0;
283  }
284  else if ((_kij[j][i] <= _kij[i][j]) && (_kij[j][i] < 0))
285  {
286  dij[i][j] = -_kij[j][i];
287  dDij_dKji[i][j] = -1.0;
288  dDii_dKji[i][j] += 1.0;
289  }
290  else
291  dij[i][j] = 0.0;
292  dij[i][i] -= dij[i][j];
293  }
294  }
295 
296  // Calculate KuzminTurek L matrix
297  // See Fig 2: L = K + D
298  std::map<dof_id_type, std::map<dof_id_type, Real>> lij;
299  for (auto & nodes : _kij)
300  {
301  const dof_id_type i = nodes.first;
302  lij[i] = {};
303  for (auto & neighbors : nodes.second)
304  {
305  const dof_id_type j = neighbors.first;
306  lij[i][j] = _kij[i][j] + dij[i][j];
307  }
308  }
309 
310  // Compute KuzminTurek R matrices
311  // See Eqns (49) and (12)
312  std::map<dof_id_type, Real> rP;
313  std::map<dof_id_type, Real> rM;
314  std::map<dof_id_type, std::map<dof_id_type, Real>>
315  drP; // drP[i][j] = d(rP[i])/d(u[j]). Here j must be connected to i (ie _kij[i][j] exists)
316  // for there to be a nonzero derivative
317  std::map<dof_id_type, std::map<dof_id_type, Real>> drM;
318  std::map<dof_id_type, std::map<dof_id_type, Real>>
319  drP_dk; // drP_dk[i][j] = d(rP[i])/d(K[i][j]). Here j must be connected to i (ie _kij[i][j]
320  // exists) for there to be a nonzero derivative
321  std::map<dof_id_type, std::map<dof_id_type, Real>> drM_dk;
322  for (const auto & nodes : _kij)
323  {
324  const dof_id_type i = nodes.first;
325  drP[i] = {};
326  drP_dk[i] = {};
327  rP[i] = rPlus(i, drP[i], drP_dk[i]);
328  drM[i] = {};
329  drM_dk[i] = {};
330  rM[i] = rMinus(i, drM[i], drM_dk[i]);
331  }
332 
333  // Calculate KuzminTurek f^{a} matrix
334  // This is the antidiffusive flux
335  // See Eqn (50)
336  std::map<dof_id_type, std::map<dof_id_type, Real>> fa;
337  std::map<dof_id_type, std::map<dof_id_type, std::map<dof_id_type, Real>>>
338  dfa; // dfa[i][j][k] = d(fa[i][j])/du[k]. Here k can be a neighbour to i or a neighbour to j
339  std::map<dof_id_type, std::map<dof_id_type, std::map<dof_id_type, Real>>>
340  dFij_dKik; // dFij_dKik[i][j][k] = d(fa[i][j])/d(K[i][k]). Here j must be connected to i, and
341  // k connected to i
342  std::map<dof_id_type, std::map<dof_id_type, std::map<dof_id_type, Real>>>
343  dFij_dKjk; // dFij_dKjk[i][j][k] = d(fa[i][j])/d(K[j][k]). Here j must be connected to i, and
344  // k connected to j (including i itself)
345  for (const auto & nodes : _kij)
346  {
347  const dof_id_type i = nodes.first;
348  fa[i] = {};
349  dfa[i] = {};
350  dFij_dKik[i] = {};
351  dFij_dKjk[i] = {};
352  for (const auto & neighbors : nodes.second)
353  {
354  const dof_id_type j = neighbors.first;
355  fa[i][j] = 0.0;
356  dfa[i][j] = {};
357  dFij_dKik[i][j] = {};
358  dFij_dKjk[i][j] = {};
359  // The derivatives are a bit complicated.
360  // If i is upwind of j then fa[i][j] depends on all nodes connected to i.
361  // But if i is downwind of j then fa[i][j] depends on all nodes connected to j.
362  for (const auto & neighbor_to_i : _kij[i])
363  {
364  dfa[i][j][neighbor_to_i.first] = 0.0;
365  dFij_dKik[i][j][neighbor_to_i.first] = 0.0;
366  }
367  for (const auto & neighbor_to_j : _kij[j])
368  {
369  dfa[i][j][neighbor_to_j.first] = 0.0;
370  dFij_dKjk[i][j][neighbor_to_j.first] = 0.0;
371  }
372  }
373  }
374  for (const auto & nodes : _kij)
375  {
376  const dof_id_type i = nodes.first;
377  const Real u_i = getUnodal(i);
378  for (const auto & neighbors : nodes.second)
379  {
380  const dof_id_type j = neighbors.first;
381  const Real u_j = getUnodal(j);
382  if (i == j)
383  continue;
384  if (lij[j][i] >= lij[i][j]) // node i is upwind of node j.
385  {
386  Real prefactor = 0.0;
387  std::map<dof_id_type, Real>
388  dprefactor_du; // dprefactor_du[global_id] = d(prefactor)/du[global_id];
389  for (const auto & dof_deriv : drP[i])
390  dprefactor_du[dof_deriv.first] = 0.0;
391  std::map<dof_id_type, Real>
392  dprefactor_dKij; // dprefactor_dKij[global_id] = d(prefactor)/dKij[i][global_id]. Here
393  // global_id must be connected to i for a nonzero derivative
394  zeroedConnection(dprefactor_dKij, i);
395  Real dprefactor_dKji = 0; // dprefactor_dKji = d(prefactor)/dKij[j][i]
396  if (u_i >= u_j)
397  {
398  if (lij[j][i] <= rP[i] * dij[i][j])
399  {
400  prefactor = lij[j][i];
401  dprefactor_dKij[j] += dDij_dKji[j][i];
402  dprefactor_dKji += 1.0 + dDij_dKij[j][i];
403  }
404  else
405  {
406  prefactor = rP[i] * dij[i][j];
407  for (const auto & dof_deriv : drP[i])
408  dprefactor_du[dof_deriv.first] = dof_deriv.second * dij[i][j];
409  dprefactor_dKij[j] += rP[i] * dDij_dKij[i][j];
410  dprefactor_dKji += rP[i] * dDij_dKji[i][j];
411  for (const auto & dof_deriv : drP_dk[i])
412  dprefactor_dKij[dof_deriv.first] += dof_deriv.second * dij[i][j];
413  }
414  }
415  else
416  {
417  if (lij[j][i] <= rM[i] * dij[i][j])
418  {
419  prefactor = lij[j][i];
420  dprefactor_dKij[j] += dDij_dKji[j][i];
421  dprefactor_dKji += 1.0 + dDij_dKij[j][i];
422  }
423  else
424  {
425  prefactor = rM[i] * dij[i][j];
426  for (const auto & dof_deriv : drM[i])
427  dprefactor_du[dof_deriv.first] = dof_deriv.second * dij[i][j];
428  dprefactor_dKij[j] += rM[i] * dDij_dKij[i][j];
429  dprefactor_dKji += rM[i] * dDij_dKji[i][j];
430  for (const auto & dof_deriv : drM_dk[i])
431  dprefactor_dKij[dof_deriv.first] += dof_deriv.second * dij[i][j];
432  }
433  }
434  fa[i][j] = prefactor * (u_i - u_j);
435  dfa[i][j][i] = prefactor;
436  dfa[i][j][j] = -prefactor;
437  for (const auto & dof_deriv : dprefactor_du)
438  dfa[i][j][dof_deriv.first] += dof_deriv.second * (u_i - u_j);
439  for (const auto & dof_deriv : dprefactor_dKij)
440  dFij_dKik[i][j][dof_deriv.first] += dof_deriv.second * (u_i - u_j);
441  dFij_dKjk[i][j][i] += dprefactor_dKji * (u_i - u_j);
442  }
443  }
444  }
445  for (const auto & nodes : _kij)
446  {
447  const dof_id_type i = nodes.first;
448  for (const auto & neighbors : nodes.second)
449  {
450  const dof_id_type j = neighbors.first;
451  if (i == j)
452  continue;
453  if (lij[j][i] < lij[i][j]) // node_i is downwind of node_j.
454  {
455  fa[i][j] = -fa[j][i];
456  for (const auto & dof_deriv : dfa[j][i])
457  dfa[i][j][dof_deriv.first] = -dof_deriv.second;
458  for (const auto & dof_deriv : dFij_dKjk[j][i])
459  dFij_dKik[i][j][dof_deriv.first] = -dof_deriv.second;
460  for (const auto & dof_deriv : dFij_dKik[j][i])
461  dFij_dKjk[i][j][dof_deriv.first] = -dof_deriv.second;
462  }
463  }
464  }
465 
466  // zero _flux_out and its derivatives
467  for (const auto & nodes : _kij)
468  {
469  const dof_id_type i = nodes.first;
470  _flux_out[i] = 0.0;
471  // The derivatives are a bit complicated.
472  // If i is upwind of a node "j" then _flux_out[i] depends on all nodes connected to i.
473  // But if i is downwind of a node "j" then _flux_out depends on all nodes connected with node
474  // j.
475  for (const auto & neighbors_i : nodes.second)
476  {
477  const dof_id_type j = neighbors_i.first;
478  _dflux_out_du[i][j] = 0.0;
479  for (const auto & neighbors_j : _kij[j])
480  {
481  _dflux_out_du[i][neighbors_j.first] = 0.0;
482  _dflux_out_dKjk[i][j][neighbors_j.first] = 0.0;
483  }
484  }
485  }
486 
487  // Add everything together
488  // See step 3 in Fig 2, noting Eqn (36)
489  for (auto & nodes : _kij)
490  {
491  const dof_id_type i = nodes.first;
492  for (const auto & neighbors : nodes.second)
493  {
494  const dof_id_type j = neighbors.first;
495  const Real u_j = getUnodal(j);
496 
497  // negative sign because residual = -Lu (KT equation (19))
498  _flux_out[i] -= lij[i][j] * u_j + fa[i][j];
499 
500  _dflux_out_du[i][j] -= lij[i][j];
501  for (const auto & dof_deriv : dfa[i][j])
502  _dflux_out_du[i][dof_deriv.first] -= dof_deriv.second;
503 
504  _dflux_out_dKjk[i][i][j] -= 1.0 * u_j; // from the K in L = K + D
505 
506  if (j == i)
507  for (const auto & dof_deriv : dDii_dKij[i])
508  _dflux_out_dKjk[i][i][dof_deriv.first] -= dof_deriv.second * u_j;
509  else
510  _dflux_out_dKjk[i][i][j] -= dDij_dKij[i][j] * u_j;
511  for (const auto & dof_deriv : dFij_dKik[i][j])
512  _dflux_out_dKjk[i][i][dof_deriv.first] -= dof_deriv.second;
513 
514  if (j == i)
515  for (const auto & dof_deriv : dDii_dKji[i])
516  _dflux_out_dKjk[i][dof_deriv.first][i] -= dof_deriv.second * u_j;
517  else
518  _dflux_out_dKjk[i][j][i] -= dDij_dKji[i][j] * u_j;
519  for (const auto & dof_deriv : dFij_dKjk[i][j])
520  _dflux_out_dKjk[i][j][dof_deriv.first] -= dof_deriv.second;
521  }
522  }
523 }
std::map< dof_id_type, std::map< dof_id_type, std::map< dof_id_type, Real > > > _dflux_out_dKjk
_dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does includ...
Real rPlus(dof_id_type node_i, std::map< dof_id_type, Real > &dlimited_du, std::map< dof_id_type, Real > &dlimited_dk) const
Returns the value of R_{i}^{+}, Eqn (49) of KT.
Real getUnodal(dof_id_type node_i) const
Returns the value of u at the global id node_i.
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.
std::map< dof_id_type, Real > _flux_out
_flux_out[i] = flux of "heat" from node i
Real rMinus(dof_id_type node_i, std::map< dof_id_type, Real > &dlimited_du, std::map< dof_id_type, Real > &dlimited_dk) const
Returns the value of R_{i}^{-}, Eqn (49) of KT.
std::map< dof_id_type, std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is con...

◆ getdFluxOutdKjk()

const std::map< dof_id_type, std::map< dof_id_type, Real > > & AdvectiveFluxCalculatorBase::getdFluxOutdKjk ( dof_id_type  node_i) const

Returns r where r[j] = d(flux out of global node i)/dK[global node j][global node k] used in Jacobian computations.

Parameters
node_iglobal id of node
Returns
the derivatives (after applying the KT procedure)

Definition at line 723 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::getdFluxOut_dvars().

724 {
725  const auto & row_find = _dflux_out_dKjk.find(node_i);
726  if (row_find == _dflux_out_dKjk.end())
727  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
728  " _dflux_out_dKjk does not contain node " + Moose::stringify(node_i));
729  return row_find->second;
730 }
std::map< dof_id_type, std::map< dof_id_type, std::map< dof_id_type, Real > > > _dflux_out_dKjk
_dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does includ...
const std::string name
Definition: Setup.h:22

◆ getdFluxOutdu()

const std::map< dof_id_type, Real > & AdvectiveFluxCalculatorBase::getdFluxOutdu ( dof_id_type  node_i) const

Returns r where r[j] = d(flux out of global node i)/du(global node j) used in Jacobian computations.

Parameters
node_iglobal id of node
Returns
the derivatives (after applying the KT procedure)

Definition at line 713 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeJacobian(), and PorousFlowAdvectiveFluxCalculatorBase::getdFluxOut_dvars().

714 {
715  const auto & row_find = _dflux_out_du.find(node_i);
716  if (row_find == _dflux_out_du.end())
717  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
718  " _dflux_out_du does not contain node " + Moose::stringify(node_i));
719  return row_find->second;
720 }
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is con...

◆ getFluxOut()

Real AdvectiveFluxCalculatorBase::getFluxOut ( dof_id_type  node_i) const

Returns the flux out of lobal node id.

Parameters
node_iid of node
Returns
advective flux out of node after applying the KT procedure

Definition at line 733 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeResidual(), and PorousFlowFluxLimitedTVDAdvection::computeResidual().

734 {
735  const auto & entry_find = _flux_out.find(node_i);
736  if (entry_find == _flux_out.end())
737  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
738  " _flux_out does not contain node " + Moose::stringify(node_i));
739  return entry_find->second;
740 }
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, Real > _flux_out
_flux_out[i] = flux of "heat" from node i

◆ getKij()

Real AdvectiveFluxCalculatorBase::getKij ( dof_id_type  node_i,
dof_id_type  node_j 
) const
protected

Returns the value of k_ij as computed by KT Eqns (18)-(20).

Parameters
node_iid of i^th node
node_jid of j^th node
Returns
k_ij of KT

Definition at line 696 of file AdvectiveFluxCalculatorBase.C.

697 {
698 
699  const auto & row_find = _kij.find(node_i);
700  if (row_find == _kij.end())
701  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() + " Kij does not contain node " +
702  Moose::stringify(node_i));
703  const std::map<dof_id_type, Real> & kij_row = row_find->second;
704  const auto & entry_find = kij_row.find(node_j);
705  if (entry_find == kij_row.end())
706  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() + " Kij on row " +
707  Moose::stringify(node_i) + " does not contain node " + Moose::stringify(node_j));
708 
709  return entry_find->second;
710 }
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ getUnodal()

Real AdvectiveFluxCalculatorBase::getUnodal ( dof_id_type  node_i) const
protected

Returns the value of u at the global id node_i.

This is _u_nodal[node_i] that has been computed by computeU during execute()

Parameters
node_iglobal node id

Definition at line 853 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize(), and PQPlusMinus().

854 {
855  const auto & node_u = _u_nodal.find(node_i);
856  if (node_u == _u_nodal.end())
857  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
858  " _u_nodal does not contain node " + Moose::stringify(node_i));
859  return node_u->second;
860 }
std::map< dof_id_type, Real > _u_nodal
_u_nodal[i] = value of _u at global node number i
const std::string name
Definition: Setup.h:22

◆ getUnodalComputedByThread()

bool AdvectiveFluxCalculatorBase::getUnodalComputedByThread ( dof_id_type  node_i) const
protected

Returns the value of _u_nodal_computed_by_thread at the global id node_i.

This will have been initialized to false in initialize() and potentially set true during execute()

Parameters
node_iglobal node id

Definition at line 863 of file AdvectiveFluxCalculatorBase.C.

Referenced by threadJoin().

864 {
865  const auto & node_u = _u_nodal_computed_by_thread.find(node_i);
866  if (node_u == _u_nodal_computed_by_thread.end())
867  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
868  " _u_nodal_computed_by_thread does not contain node " + Moose::stringify(node_i));
869  return node_u->second;
870 }
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...

◆ getValence()

unsigned AdvectiveFluxCalculatorBase::getValence ( dof_id_type  node_i,
dof_id_type  node_j 
) const

Returns the valence of the i-j edge.

Valence is the number of times the edge is encountered in a loop over elements (that have appropriate subdomain_id, if the user has employed the "blocks=" parameter) seen by this processor (including ghosted elements)

Parameters
node_iid of i^th node
node_jid of j^th node
Returns
valence of the i-j edge

Definition at line 743 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeJacobian(), PorousFlowFluxLimitedTVDAdvection::computeJacobian(), FluxLimitedTVDAdvection::computeResidual(), and PorousFlowFluxLimitedTVDAdvection::computeResidual().

744 {
745  const std::pair<dof_id_type, dof_id_type> i_j(node_i, node_j);
746  const auto & entry_find = _valence.find(i_j);
747  if (entry_find == _valence.end())
748  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() +
749  " Valence does not contain node-pair " + Moose::stringify(node_i) + " to " +
750  Moose::stringify(node_j));
751  return entry_find->second;
752 }
const std::string name
Definition: Setup.h:22
std::map< std::pair< dof_id_type, dof_id_type >, unsigned > _valence
_valence[(i, j)] = number of times, in a loop over elements seen by this processor (viz...

◆ initialize()

void AdvectiveFluxCalculatorBase::initialize ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 158 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::initialize().

159 {
160  // Zero _kij and falsify _u_nodal_computed_by_thread ready for building in execute() and
161  // finalize()
162  for (auto & nodes : _kij)
163  {
164  _u_nodal_computed_by_thread[nodes.first] = false;
165  for (auto & neighbors : nodes.second)
166  neighbors.second = 0.0;
167  }
168 }
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ limitFlux()

void AdvectiveFluxCalculatorBase::limitFlux ( Real  a,
Real  b,
Real &  limited,
Real &  dlimited_db 
) const
protected

flux limiter, L, on Page 135 of Kuzmin and Turek

Parameters
aKT's "a" parameter
bKT's "b" parameter
limited[out]The value of the flux limiter, L
dlimited_db[out]The derivative dL/db

Definition at line 596 of file AdvectiveFluxCalculatorBase.C.

Referenced by rMinus(), and rPlus().

597 {
598  limited = 0.0;
599  dlimited_db = 0.0;
601  return;
602 
603  if ((a >= 0.0 && b <= 0.0) || (a <= 0.0 && b >= 0.0))
604  return;
605  const Real s = (a > 0.0 ? 1.0 : -1.0);
606 
607  const Real lal = std::abs(a);
608  const Real lbl = std::abs(b);
609  const Real dlbl = (b >= 0.0 ? 1.0 : -1.0); // d(lbl)/db
610  switch (_flux_limiter_type)
611  {
613  {
614  if (lal <= lbl)
615  {
616  limited = s * lal;
617  dlimited_db = 0.0;
618  }
619  else
620  {
621  limited = s * lbl;
622  dlimited_db = s * dlbl;
623  }
624  return;
625  }
627  {
628  limited = s * 2 * lal * lbl / (lal + lbl);
629  dlimited_db = s * 2 * lal * (dlbl / (lal + lbl) - lbl * dlbl / std::pow(lal + lbl, 2));
630  return;
631  }
633  {
634  const Real av = 0.5 * std::abs(a + b);
635  if (2 * lal <= av && lal <= lbl)
636  {
637  // 2 * lal is the smallest
638  limited = s * 2.0 * lal;
639  dlimited_db = 0.0;
640  }
641  else if (2 * lbl <= av && lbl <= lal)
642  {
643  // 2 * lbl is the smallest
644  limited = s * 2.0 * lbl;
645  dlimited_db = s * 2.0 * dlbl;
646  }
647  else
648  {
649  // av is the smallest
650  limited = s * av;
651  // if (a>0 and b>0) then d(av)/db = 0.5 = 0.5 * dlbl
652  // if (a<0 and b<0) then d(av)/db = -0.5 = 0.5 * dlbl
653  // if a and b have different sign then limited=0, above
654  dlimited_db = s * 0.5 * dlbl;
655  }
656  return;
657  }
659  {
660  const Real term1 = std::min(2.0 * lal, lbl);
661  const Real term2 = std::min(lal, 2.0 * lbl);
662  if (term1 >= term2)
663  {
664  if (2.0 * lal <= lbl)
665  {
666  limited = s * 2 * lal;
667  dlimited_db = 0.0;
668  }
669  else
670  {
671  limited = s * lbl;
672  dlimited_db = s * dlbl;
673  }
674  }
675  else
676  {
677  if (lal <= 2.0 * lbl)
678  {
679  limited = s * lal;
680  dlimited_db = 0.0;
681  }
682  else
683  {
684  limited = s * 2.0 * lbl;
685  dlimited_db = s * 2.0 * dlbl;
686  }
687  }
688  return;
689  }
690  default:
691  return;
692  }
693 }
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ meshChanged()

void AdvectiveFluxCalculatorBase::meshChanged ( )
overridevirtual

Definition at line 149 of file AdvectiveFluxCalculatorBase.C.

150 {
151  ElementUserObject::meshChanged();
152 
153  // Signal that _kij, _valence, etc need to be rebuilt
154  _resizing_needed = true;
155 }
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...

◆ PQPlusMinus()

Real AdvectiveFluxCalculatorBase::PQPlusMinus ( dof_id_type  node_i,
const PQPlusMinusEnum  pq_plus_minus,
std::map< dof_id_type, Real > &  derivs,
std::map< dof_id_type, Real > &  dpq_dk 
) const
protected

Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which are defined in Eqns (47) and (48) of KT.

Parameters
node_iglobal nodal id
pq_plus_minusindicates whether P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} should be returned
derivs[out]derivs[j] = d(result)/d(u[global node j])
dpq_dk[out]dpq_dk[j] = d(result)/d(K[node_i][global node j]). Recall that d(result)/d(K[l][m]) are zero unless l=node_i

Definition at line 768 of file AdvectiveFluxCalculatorBase.C.

Referenced by rMinus(), and rPlus().

772 {
773  // Find the value of u at node_i
774  const Real u_i = getUnodal(node_i);
775 
776  // We're going to sum over all nodes connected with node_i
777  // These nodes are found in _kij[node_i]
778  const auto & row_find = _kij.find(node_i);
779  if (row_find == _kij.end())
780  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() + " Kij does not contain node " +
781  Moose::stringify(node_i));
782  const std::map<dof_id_type, Real> nodej_and_kij = row_find->second;
783 
784  // Initialize the results
785  Real result = 0.0;
786  zeroedConnection(derivs, node_i);
787  zeroedConnection(dpqdk, node_i);
788 
789  // Sum over all nodes connected with node_i.
790  for (const auto & nk : nodej_and_kij)
791  {
792  const dof_id_type node_j = nk.first;
793  const Real kentry = nk.second;
794 
795  // Find the value of u at node_j
796  const Real u_j = getUnodal(node_j);
797  const Real ujminusi = u_j - u_i;
798 
799  // Evaluate the i-j contribution to the result
800  switch (pq_plus_minus)
801  {
803  {
804  if (ujminusi < 0.0 && kentry < 0.0)
805  {
806  result += kentry * ujminusi;
807  derivs[node_j] += kentry;
808  derivs[node_i] -= kentry;
809  dpqdk[node_j] += ujminusi;
810  }
811  break;
812  }
814  {
815  if (ujminusi > 0.0 && kentry < 0.0)
816  {
817  result += kentry * ujminusi;
818  derivs[node_j] += kentry;
819  derivs[node_i] -= kentry;
820  dpqdk[node_j] += ujminusi;
821  }
822  break;
823  }
825  {
826  if (ujminusi > 0.0 && kentry > 0.0)
827  {
828  result += kentry * ujminusi;
829  derivs[node_j] += kentry;
830  derivs[node_i] -= kentry;
831  dpqdk[node_j] += ujminusi;
832  }
833  break;
834  }
836  {
837  if (ujminusi < 0.0 && kentry > 0.0)
838  {
839  result += kentry * ujminusi;
840  derivs[node_j] += kentry;
841  derivs[node_i] -= kentry;
842  dpqdk[node_j] += ujminusi;
843  }
844  break;
845  }
846  }
847  }
848 
849  return result;
850 }
const std::string name
Definition: Setup.h:22
Real getUnodal(dof_id_type node_i) const
Returns the value of u at the global id node_i.
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.

◆ rMinus()

Real AdvectiveFluxCalculatorBase::rMinus ( dof_id_type  node_i,
std::map< dof_id_type, Real > &  dlimited_du,
std::map< dof_id_type, Real > &  dlimited_dk 
) const
protected

Returns the value of R_{i}^{-}, Eqn (49) of KT.

Parameters
node_inodal id
dlimited_du[out]dlimited_du[node_id] = d(R_{i}^{-})/du[node_id]
dlimited_dk[out]dlimited_dk[node_id] = d(R_{i}^{-})/d(K[i][node_id]). Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 561 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize().

564 {
565  zeroedConnection(dlimited_du, node_i);
566  zeroedConnection(dlimited_dk, node_i);
568  return 0.0;
569  std::map<dof_id_type, Real> dp_du;
570  std::map<dof_id_type, Real> dp_dk;
571  const Real p = PQPlusMinus(node_i, PQPlusMinusEnum::PMinus, dp_du, dp_dk);
572  if (p == 0.0)
573  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
574  return 1.0;
575  std::map<dof_id_type, Real> dq_du;
576  std::map<dof_id_type, Real> dq_dk;
577  const Real q = PQPlusMinus(node_i, PQPlusMinusEnum::QMinus, dq_du, dq_dk);
578  const Real r = q / p;
579  std::map<dof_id_type, Real> dr_du;
580  for (const auto & node_deriv : dp_du)
581  dr_du[node_deriv.first] = dq_du[node_deriv.first] / p - q * node_deriv.second / std::pow(p, 2);
582  std::map<dof_id_type, Real> dr_dk;
583  for (const auto & node_deriv : dp_dk)
584  dr_dk[node_deriv.first] = dq_dk[node_deriv.first] / p - q * node_deriv.second / std::pow(p, 2);
585  Real limited;
586  Real dlimited_dr;
587  limitFlux(1.0, r, limited, dlimited_dr);
588  for (const auto & node_deriv : dr_du)
589  dlimited_du[node_deriv.first] = dlimited_dr * node_deriv.second;
590  for (const auto & node_deriv : dr_dk)
591  dlimited_dk[node_deriv.first] = dlimited_dr * node_deriv.second;
592  return limited;
593 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
Real PQPlusMinus(dof_id_type node_i, const PQPlusMinusEnum pq_plus_minus, std::map< dof_id_type, Real > &derivs, std::map< dof_id_type, Real > &dpq_dk) const
Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which ...
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.

◆ rPlus()

Real AdvectiveFluxCalculatorBase::rPlus ( dof_id_type  node_i,
std::map< dof_id_type, Real > &  dlimited_du,
std::map< dof_id_type, Real > &  dlimited_dk 
) const
protected

Returns the value of R_{i}^{+}, Eqn (49) of KT.

Parameters
node_inodal id
dlimited_du[out]dlimited_du[node_id] = d(R_{i}^{+})/du[node_id] (Here node_id is a global node id)
dlimited_dk[out]dlimited_dk[node_id] = d(R_{i}^{+})/d(K[i][node_id]). Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 526 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize().

529 {
530  zeroedConnection(dlimited_du, node_i);
531  zeroedConnection(dlimited_dk, node_i);
533  return 0.0;
534  std::map<dof_id_type, Real> dp_du;
535  std::map<dof_id_type, Real> dp_dk;
536  const Real p = PQPlusMinus(node_i, PQPlusMinusEnum::PPlus, dp_du, dp_dk);
537  if (p == 0.0)
538  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
539  return 1.0;
540  std::map<dof_id_type, Real> dq_du;
541  std::map<dof_id_type, Real> dq_dk;
542  const Real q = PQPlusMinus(node_i, PQPlusMinusEnum::QPlus, dq_du, dq_dk);
543  const Real r = q / p;
544  std::map<dof_id_type, Real> dr_du;
545  for (const auto & node_deriv : dp_du)
546  dr_du[node_deriv.first] = dq_du[node_deriv.first] / p - q * node_deriv.second / std::pow(p, 2);
547  std::map<dof_id_type, Real> dr_dk;
548  for (const auto & node_deriv : dp_dk)
549  dr_dk[node_deriv.first] = dq_dk[node_deriv.first] / p - q * node_deriv.second / std::pow(p, 2);
550  Real limited;
551  Real dlimited_dr;
552  limitFlux(1.0, r, limited, dlimited_dr);
553  for (const auto & node_deriv : dr_du)
554  dlimited_du[node_deriv.first] = dlimited_dr * node_deriv.second;
555  for (const auto & node_deriv : dr_dk)
556  dlimited_dk[node_deriv.first] = dlimited_dr * node_deriv.second;
557  return limited;
558 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
Real PQPlusMinus(dof_id_type node_i, const PQPlusMinusEnum pq_plus_minus, std::map< dof_id_type, Real > &derivs, std::map< dof_id_type, Real > &dpq_dk) const
Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which ...
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.

◆ threadJoin()

void AdvectiveFluxCalculatorBase::threadJoin ( const UserObject &  uo)
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 201 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::threadJoin().

202 {
203  const AdvectiveFluxCalculatorBase & afc = static_cast<const AdvectiveFluxCalculatorBase &>(uo);
204  // add the values of _kij computed by different threads
205  for (const auto & nodes : afc._kij)
206  {
207  const dof_id_type i = nodes.first;
208  for (const auto & neighbors : nodes.second)
209  {
210  const dof_id_type j = neighbors.first;
211  const Real afc_val = neighbors.second;
212  _kij[i][j] += afc_val;
213  }
214  }
215  // gather the values of _u_nodal computed by different threads
216  for (const auto & node_u : afc._u_nodal)
217  {
218  const dof_id_type i = node_u.first;
220  _u_nodal[i] = node_u.second;
221  }
222 }
std::map< dof_id_type, Real > _u_nodal
_u_nodal[i] = value of _u at global node number i
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...
Base class to compute Advective fluxes.
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.
bool getUnodalComputedByThread(dof_id_type node_i) const
Returns the value of _u_nodal_computed_by_thread at the global id node_i.

◆ timestepSetup()

void AdvectiveFluxCalculatorBase::timestepSetup ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 55 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

56 {
57  // If needed, size and initialize quantities appropriately, and compute _valence
58  if (_resizing_needed)
59  {
60  _kij.clear();
61 
62  // QUERY: does this properly account for multiple processors, threading, and other things i
63  // haven't thought about?
64 
65  /*
66  * Initialize _kij for all nodes that can be seen by this processor and on relevant blocks
67  *
68  * MULTIPROC NOTE: this must loop over local elements and 2 layers of ghosted elements.
69  * The Kernel will only loop over local elements, so will only use _kij for
70  * linked node-node pairs that appear in the local elements. Nevertheless, we
71  * need to build _kij for the nodes in the ghosted elements in order to simplify
72  * Jacobian computations
73  */
74  for (const auto & elem : _subproblem.mesh().getMesh().active_element_ptr_range())
75  {
76  if (this->hasBlocks(elem->subdomain_id()))
77  {
78  for (unsigned i = 0; i < elem->n_nodes(); ++i)
79  {
80  const dof_id_type node_i = elem->node_id(i);
81  if (_kij.find(node_i) == _kij.end())
82  _kij[node_i] = {};
83  for (unsigned j = 0; j < elem->n_nodes(); ++j)
84  {
85  const dof_id_type node_j = elem->node_id(j);
86  _kij[node_i][node_j] = 0.0;
87  }
88  }
89  }
90  }
91 
92  /*
93  * Build _valence[i_j], which is the number of times the i_j pair is encountered when looping
94  * over the entire mesh (and on relevant blocks)
95  *
96  * MULTIPROC NOTE: this must loop over local elements and >=1 layer of ghosted elements.
97  * The Kernel will only loop over local elements, so will only use _valence for
98  * linked node-node pairs that appear in the local elements. But other processors will
99  * loop over neighboring elements, so avoid multiple counting of the residual and Jacobian
100  * this processor must record how many times each node-node link of its local elements
101  * appears in the *global* mesh and pass that info to the Kernel
102  */
103  _valence.clear();
104  for (const auto & elem : _fe_problem.getEvaluableElementRange())
105  {
106  if (this->hasBlocks(elem->subdomain_id()))
107  {
108  for (unsigned i = 0; i < elem->n_nodes(); ++i)
109  {
110  const dof_id_type node_i = elem->node_id(i);
111  for (unsigned j = 0; j < elem->n_nodes(); ++j)
112  {
113  const dof_id_type node_j = elem->node_id(j);
114  const std::pair<dof_id_type, dof_id_type> i_j(node_i, node_j);
115  if (_valence.find(i_j) == _valence.end())
116  _valence[i_j] = 0;
117  _valence[i_j] += 1;
118  }
119  }
120  }
121  }
122 
123  _u_nodal.clear();
125  _flux_out.clear();
126  _dflux_out_du.clear();
127  _dflux_out_dKjk.clear();
128  for (auto & nodes : _kij)
129  {
130  const dof_id_type node_i = nodes.first;
131  _u_nodal[node_i] = 0.0;
132  _u_nodal_computed_by_thread[node_i] = false;
133  _flux_out[node_i] = 0.0;
134  _dflux_out_du[node_i] = {};
135  zeroedConnection(_dflux_out_du[node_i], node_i);
136  _dflux_out_dKjk[node_i] = {};
137  for (const auto & neighbour_to_i : nodes.second)
138  {
139  const dof_id_type node_j = neighbour_to_i.first;
140  _dflux_out_dKjk[node_i][node_j] = {};
141  zeroedConnection(_dflux_out_dKjk[node_i][node_j], node_j);
142  }
143  }
144  _resizing_needed = false;
145  }
146 }
std::map< dof_id_type, Real > _u_nodal
_u_nodal[i] = value of _u at global node number i
std::map< dof_id_type, std::map< dof_id_type, std::map< dof_id_type, Real > > > _dflux_out_dKjk
_dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does includ...
std::map< dof_id_type, bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on ...
std::map< std::pair< dof_id_type, dof_id_type >, unsigned > _valence
_valence[(i, j)] = number of times, in a loop over elements seen by this processor (viz...
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.
std::map< dof_id_type, Real > _flux_out
_flux_out[i] = flux of "heat" from node i
std::map< dof_id_type, std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is con...

◆ zeroedConnection()

void AdvectiveFluxCalculatorBase::zeroedConnection ( std::map< dof_id_type, Real > &  the_map,
dof_id_type  node_i 
) const
protected

Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.

Parameters
[out]the_mapthe map to be zeroed appropriately
[in]node_inodal id

Definition at line 755 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize(), PQPlusMinus(), rMinus(), rPlus(), and timestepSetup().

757 {
758  the_map.clear();
759  const auto & row_find = _kij.find(node_i);
760  if (row_find == _kij.end())
761  mooseError("AdvectiveFluxCalculatorBase UserObject " + name() + " Kij does not contain node " +
762  Moose::stringify(node_i));
763  for (const auto & nk : row_find->second)
764  the_map[nk.first] = 0.0;
765 }
const std::string name
Definition: Setup.h:22
std::map< dof_id_type, std::map< dof_id_type, Real > > _kij
Kuzmin-Turek K_ij matrix.

Member Data Documentation

◆ _dflux_out_dKjk

std::map<dof_id_type, std::map<dof_id_type, std::map<dof_id_type, Real> > > AdvectiveFluxCalculatorBase::_dflux_out_dKjk
protected

_dflux_out_dKjk[i][j][k] = d(flux_out[i])/d(K[j][k]). Here j must be connected to i (this does include j = i), and k must be connected to j (this does include k = i and k = j)

Definition at line 179 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getdFluxOutdKjk(), and timestepSetup().

◆ _dflux_out_du

std::map<dof_id_type, std::map<dof_id_type, Real> > AdvectiveFluxCalculatorBase::_dflux_out_du
protected

_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). Here j must be connected to i, or to a node that is connected to i.

Definition at line 176 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getdFluxOutdu(), and timestepSetup().

◆ _flux_limiter_type

enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum AdvectiveFluxCalculatorBase::_flux_limiter_type
protected

Referenced by limitFlux(), rMinus(), and rPlus().

◆ _flux_out

std::map<dof_id_type, Real> AdvectiveFluxCalculatorBase::_flux_out
protected

_flux_out[i] = flux of "heat" from node i

Definition at line 173 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getFluxOut(), and timestepSetup().

◆ _kij

std::map<dof_id_type, std::map<dof_id_type, Real> > AdvectiveFluxCalculatorBase::_kij
protected

Kuzmin-Turek K_ij matrix.

Along with R+ and R-, this is the key quantity computed by this UserObject. _kij[i][j] = k_ij corresponding to the i-j node pair. Because it explicitly holds info which nodes are paired with which other nodes, it is also used to perform: given a node_id, loop over all neighbouring nodes. This occurs in the computation of R+ and R-.

Definition at line 170 of file AdvectiveFluxCalculatorBase.h.

Referenced by executeOnElement(), finalize(), getKij(), initialize(), PorousFlowAdvectiveFluxCalculatorBase::initialize(), PQPlusMinus(), threadJoin(), timestepSetup(), PorousFlowAdvectiveFluxCalculatorBase::timestepSetup(), and zeroedConnection().

◆ _resizing_needed

bool AdvectiveFluxCalculatorBase::_resizing_needed
protected

whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timestep

Definition at line 116 of file AdvectiveFluxCalculatorBase.h.

Referenced by meshChanged(), timestepSetup(), and PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

◆ _u_nodal

std::map<dof_id_type, Real> AdvectiveFluxCalculatorBase::_u_nodal
protected

_u_nodal[i] = value of _u at global node number i

Definition at line 187 of file AdvectiveFluxCalculatorBase.h.

Referenced by execute(), getUnodal(), threadJoin(), and timestepSetup().

◆ _u_nodal_computed_by_thread

std::map<dof_id_type, bool> AdvectiveFluxCalculatorBase::_u_nodal_computed_by_thread
protected

_u_nodal_computed_by_thread[i] = true if _u_nodal[i] has been computed in execute() by the thread on this processor

Definition at line 190 of file AdvectiveFluxCalculatorBase.h.

Referenced by execute(), getUnodalComputedByThread(), initialize(), threadJoin(), and timestepSetup().

◆ _valence

std::map<std::pair<dof_id_type, dof_id_type>, unsigned> AdvectiveFluxCalculatorBase::_valence
protected

_valence[(i, j)] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the i-j edge is encountered

Definition at line 184 of file AdvectiveFluxCalculatorBase.h.

Referenced by getValence(), and timestepSetup().


The documentation for this class was generated from the following files: