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

Computes Advective fluxes for a constant velocity. More...

#include <AdvectiveFluxCalculatorConstantVelocity.h>

Inheritance diagram for AdvectiveFluxCalculatorConstantVelocity:
[legend]

Public Member Functions

 AdvectiveFluxCalculatorConstantVelocity (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 override
 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 override
 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

RealVectorValue _velocity
 advection velocity More...
 
const VariableValue & _u_nodal
 the nodal values of u More...
 
const VariablePhiValue & _phi
 Kuzmin-Turek shape function. More...
 
const VariablePhiGradient & _grad_phi
 grad(Kuzmin-Turek shape function) More...
 
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, 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

Computes Advective fluxes for a constant velocity.

Definition at line 23 of file AdvectiveFluxCalculatorConstantVelocity.h.

Member Enumeration Documentation

◆ FluxLimiterTypeEnum

enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum
strongprotectedinherited

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

enum AdvectiveFluxCalculatorBase::PQPlusMinusEnum
strongprotectedinherited

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

◆ AdvectiveFluxCalculatorConstantVelocity()

AdvectiveFluxCalculatorConstantVelocity::AdvectiveFluxCalculatorConstantVelocity ( const InputParameters &  parameters)

Definition at line 26 of file AdvectiveFluxCalculatorConstantVelocity.C.

28  : AdvectiveFluxCalculatorBase(parameters),
29  _velocity(getParam<RealVectorValue>("velocity")),
30  _u_nodal(coupledNodalValue("u")),
31  _phi(_assembly.fePhi<Real>(getVar("u", 0)->feType())),
32  _grad_phi(_assembly.feGradPhi<Real>(getVar("u", 0)->feType()))
33 {
34 }
const VariablePhiValue & _phi
Kuzmin-Turek shape function.
const VariablePhiGradient & _grad_phi
grad(Kuzmin-Turek shape function)
AdvectiveFluxCalculatorBase(const InputParameters &parameters)
const VariableValue & _u_nodal
the nodal values of u

Member Function Documentation

◆ computeU()

Real AdvectiveFluxCalculatorConstantVelocity::computeU ( unsigned  i) const
overrideprotectedvirtual

Computes the value of u at the local node id of the current element (_current_elem)

Parameters
ilocal node id of the current element

Implements AdvectiveFluxCalculatorBase.

Definition at line 43 of file AdvectiveFluxCalculatorConstantVelocity.C.

44 {
45  return _u_nodal[i];
46 }
const VariableValue & _u_nodal
the nodal values of u

◆ computeVelocity()

Real AdvectiveFluxCalculatorConstantVelocity::computeVelocity ( unsigned  i,
unsigned  j,
unsigned  qp 
) const
overrideprotectedvirtual

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

Implements AdvectiveFluxCalculatorBase.

Definition at line 37 of file AdvectiveFluxCalculatorConstantVelocity.C.

38 {
39  return (_grad_phi[i][qp] * _velocity) * _phi[j][qp];
40 }
const VariablePhiValue & _phi
Kuzmin-Turek shape function.
const VariablePhiGradient & _grad_phi
grad(Kuzmin-Turek shape function)

◆ execute()

void AdvectiveFluxCalculatorBase::execute ( )
overridevirtualinherited

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 171 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::execute(), and AdvectiveFluxCalculatorBase::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 
)
virtualinherited

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 AdvectiveFluxCalculatorBase::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 ( )
overridevirtualinherited

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
inherited

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
inherited

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
inherited

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
protectedinherited

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
protectedinherited

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 AdvectiveFluxCalculatorBase::finalize(), and AdvectiveFluxCalculatorBase::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
protectedinherited

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 AdvectiveFluxCalculatorBase::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
inherited

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 ( )
overridevirtualinherited

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
protectedinherited

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 AdvectiveFluxCalculatorBase::rMinus(), and AdvectiveFluxCalculatorBase::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 ( )
overridevirtualinherited

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
protectedinherited

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 AdvectiveFluxCalculatorBase::rMinus(), and AdvectiveFluxCalculatorBase::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
protectedinherited

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 AdvectiveFluxCalculatorBase::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
protectedinherited

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 AdvectiveFluxCalculatorBase::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)
overridevirtualinherited

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 ( )
overridevirtualinherited

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
protectedinherited

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 AdvectiveFluxCalculatorBase::finalize(), AdvectiveFluxCalculatorBase::PQPlusMinus(), AdvectiveFluxCalculatorBase::rMinus(), AdvectiveFluxCalculatorBase::rPlus(), and AdvectiveFluxCalculatorBase::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
protectedinherited

_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 AdvectiveFluxCalculatorBase::finalize(), AdvectiveFluxCalculatorBase::getdFluxOutdKjk(), and AdvectiveFluxCalculatorBase::timestepSetup().

◆ _dflux_out_du

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

_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 AdvectiveFluxCalculatorBase::finalize(), AdvectiveFluxCalculatorBase::getdFluxOutdu(), and AdvectiveFluxCalculatorBase::timestepSetup().

◆ _flux_limiter_type

enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum AdvectiveFluxCalculatorBase::_flux_limiter_type
protectedinherited

◆ _flux_out

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

◆ _grad_phi

const VariablePhiGradient& AdvectiveFluxCalculatorConstantVelocity::_grad_phi
protected

grad(Kuzmin-Turek shape function)

Definition at line 43 of file AdvectiveFluxCalculatorConstantVelocity.h.

Referenced by computeVelocity().

◆ _kij

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

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 AdvectiveFluxCalculatorBase::executeOnElement(), AdvectiveFluxCalculatorBase::finalize(), AdvectiveFluxCalculatorBase::getKij(), AdvectiveFluxCalculatorBase::initialize(), PorousFlowAdvectiveFluxCalculatorBase::initialize(), AdvectiveFluxCalculatorBase::PQPlusMinus(), AdvectiveFluxCalculatorBase::threadJoin(), AdvectiveFluxCalculatorBase::timestepSetup(), PorousFlowAdvectiveFluxCalculatorBase::timestepSetup(), and AdvectiveFluxCalculatorBase::zeroedConnection().

◆ _phi

const VariablePhiValue& AdvectiveFluxCalculatorConstantVelocity::_phi
protected

Kuzmin-Turek shape function.

Definition at line 40 of file AdvectiveFluxCalculatorConstantVelocity.h.

Referenced by computeVelocity().

◆ _resizing_needed

bool AdvectiveFluxCalculatorBase::_resizing_needed
protectedinherited

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 AdvectiveFluxCalculatorBase::meshChanged(), AdvectiveFluxCalculatorBase::timestepSetup(), and PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

◆ _u_nodal

const VariableValue& AdvectiveFluxCalculatorConstantVelocity::_u_nodal
protected

the nodal values of u

Definition at line 37 of file AdvectiveFluxCalculatorConstantVelocity.h.

Referenced by computeU().

◆ _u_nodal_computed_by_thread

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

◆ _valence

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

_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 AdvectiveFluxCalculatorBase::getValence(), and AdvectiveFluxCalculatorBase::timestepSetup().

◆ _velocity

RealVectorValue AdvectiveFluxCalculatorConstantVelocity::_velocity
protected

advection velocity

Definition at line 34 of file AdvectiveFluxCalculatorConstantVelocity.h.

Referenced by computeVelocity().


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