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::vector< std::vector< Real > > & getdFluxOutdKjk (dof_id_type node_i) const
 Returns r where r[j][k] = d(flux out of global node i)/dK[connected node j][connected node k] used in Jacobian computations. More...
 
unsigned getValence (dof_id_type node_i) const
 Returns the valence of the global node i Valence is the number of times the node 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) 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 sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
 Returns the value of R_{i}^{+}, Eqn (49) of KT. More...
 
Real rMinus (dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
 Returns the value of R_{i}^{-}, Eqn (49) of KT. More...
 
Real PQPlusMinus (dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< 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...
 

Protected Attributes

RealVectorValue _velocity
 advection velocity More...
 
const VariableValue & _u_at_nodes
 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::vector< std::vector< Real > > _kij
 Kuzmin-Turek K_ij matrix. More...
 
std::vector< Real > _flux_out
 _flux_out[i] = flux of "heat" from sequential node i More...
 
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
 _dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). More...
 
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
 _dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]). More...
 
std::vector< unsigned > _valence
 _valence[i] = 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 sequential node i is encountered More...
 
std::vector< Real > _u_nodal
 _u_nodal[i] = value of _u at sequential node number i More...
 
std::vector< 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...
 
PorousFlowConnectedNodes _connections
 Holds the sequential and global nodal IDs, and info regarding mesh connections between them. More...
 
std::size_t _number_of_nodes
 Number of nodes held by the _connections object. 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 154 of file AdvectiveFluxCalculatorBase.h.

154 { 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 195 of file AdvectiveFluxCalculatorBase.h.

196  {
197  PPlus,
198  PMinus,
199  QPlus,
200  QMinus
201  };

Constructor & Destructor Documentation

◆ AdvectiveFluxCalculatorConstantVelocity()

AdvectiveFluxCalculatorConstantVelocity::AdvectiveFluxCalculatorConstantVelocity ( const InputParameters &  parameters)

Definition at line 27 of file AdvectiveFluxCalculatorConstantVelocity.C.

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

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 44 of file AdvectiveFluxCalculatorConstantVelocity.C.

45 {
46  return _u_at_nodes[i];
47 }
const VariableValue & _u_at_nodes
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 38 of file AdvectiveFluxCalculatorConstantVelocity.C.

39 {
40  return (_grad_phi[i][qp] * _velocity) * _phi[j][qp];
41 }
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 163 of file AdvectiveFluxCalculatorBase.C.

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

164 {
165  // compute _kij contributions from this element that is local to this processor
166  // and record _u_nodal
167  for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
168  {
169  const dof_id_type node_i = _current_elem->node_id(i);
170  const dof_id_type sequential_i = _connections.sequentialID(node_i);
171  if (!_u_nodal_computed_by_thread[sequential_i])
172  {
173  _u_nodal[sequential_i] = computeU(i);
174  _u_nodal_computed_by_thread[sequential_i] = true;
175  }
176  for (unsigned j = 0; j < _current_elem->n_nodes(); ++j)
177  {
178  const dof_id_type node_j = _current_elem->node_id(j);
179  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
180  executeOnElement(node_i, node_j, i, j, qp);
181  }
182  }
183 }
std::vector< 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 ...
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
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...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
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 186 of file AdvectiveFluxCalculatorBase.C.

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

188 {
189  // KT Eqn (20)
190  const dof_id_type sequential_i = _connections.sequentialID(global_i);
191  const unsigned index_i_to_j = _connections.indexOfGlobalConnection(global_i, global_j);
192  _kij[sequential_i][index_i_to_j] += _JxW[qp] * _coord[qp] * computeVelocity(local_i, local_j, qp);
193 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
unsigned indexOfGlobalConnection(dof_id_type global_node_ID_from, dof_id_type global_node_ID_to) const
Return the index of global_node_ID_to in the globalConnectionsToGlobalID(global_node_ID_from) vector...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
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::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ finalize()

void AdvectiveFluxCalculatorBase::finalize ( )
overridevirtualinherited

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 215 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::finalize().

216 {
217  // Overall: ensure _kij is fully built, then compute Kuzmin-Turek D, L, f^a and
218  // relevant Jacobian information, and then the relevant quantities into _flux_out and
219  // _dflux_out_du, _dflux_out_dKjk
220 
221  /*
222  * MULTIPROC NOTE: execute() only computed contributions from the elements
223  * local to this processor. Now do the same for the 2 layers of ghosted elements
224  */
225  for (const auto & elem_id : _fe_problem.ghostedElems())
226  {
227  _current_elem = _mesh.elem(elem_id);
228  if (this->hasBlocks(_current_elem->subdomain_id()))
229  {
230  _fe_problem.prepare(_current_elem, _tid);
231  _fe_problem.reinitElem(_current_elem, _tid);
232  execute();
233  }
234  }
235 
236  // Calculate KuzminTurek D matrix
237  // See Eqn (32)
238  // This adds artificial diffusion, which eliminates any spurious oscillations
239  // The idea is that D will remove all negative off-diagonal elements when it is added to K
240  // This is identical to full upwinding
241  std::vector<std::vector<Real>> dij(_number_of_nodes);
242  std::vector<std::vector<Real>> dDij_dKij(
243  _number_of_nodes); // dDij_dKij[i][j] = d(D[i][j])/d(K[i][j]) for i!=j
244  std::vector<std::vector<Real>> dDij_dKji(
245  _number_of_nodes); // dDij_dKji[i][j] = d(D[i][j])/d(K[j][i]) for i!=j
246  std::vector<std::vector<Real>> dDii_dKij(
247  _number_of_nodes); // dDii_dKij[i][j] = d(D[i][i])/d(K[i][j])
248  std::vector<std::vector<Real>> dDii_dKji(
249  _number_of_nodes); // dDii_dKji[i][j] = d(D[i][i])/d(K[j][i])
250  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
251  {
252  const std::vector<dof_id_type> con_i =
254  const std::size_t num_con_i = con_i.size();
255  dij[sequential_i].assign(num_con_i, 0.0);
256  dDij_dKij[sequential_i].assign(num_con_i, 0.0);
257  dDij_dKji[sequential_i].assign(num_con_i, 0.0);
258  dDii_dKij[sequential_i].assign(num_con_i, 0.0);
259  dDii_dKji[sequential_i].assign(num_con_i, 0.0);
260  const unsigned index_i_to_i =
261  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
262  for (std::size_t j = 0; j < num_con_i; ++j)
263  {
264  const dof_id_type sequential_j = con_i[j];
265  if (sequential_i == sequential_j)
266  continue;
267  const unsigned index_j_to_i =
268  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
269  const Real kij = _kij[sequential_i][j];
270  const Real kji = _kij[sequential_j][index_j_to_i];
271  if ((kij <= kji) && (kij < 0.0))
272  {
273  dij[sequential_i][j] = -kij;
274  dDij_dKij[sequential_i][j] = -1.0;
275  dDii_dKij[sequential_i][j] += 1.0;
276  }
277  else if ((kji <= kij) && (kji < 0.0))
278  {
279  dij[sequential_i][j] = -kji;
280  dDij_dKji[sequential_i][j] = -1.0;
281  dDii_dKji[sequential_i][j] += 1.0;
282  }
283  else
284  dij[sequential_i][j] = 0.0;
285  dij[sequential_i][index_i_to_i] -= dij[sequential_i][j];
286  }
287  }
288 
289  // Calculate KuzminTurek L matrix
290  // See Fig 2: L = K + D
291  std::vector<std::vector<Real>> lij(_number_of_nodes);
292  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
293  {
294  const std::vector<dof_id_type> con_i =
296  const std::size_t num_con_i = con_i.size();
297  lij[sequential_i].assign(num_con_i, 0.0);
298  for (std::size_t j = 0; j < num_con_i; ++j)
299  lij[sequential_i][j] = _kij[sequential_i][j] + dij[sequential_i][j];
300  }
301 
302  // Compute KuzminTurek R matrices
303  // See Eqns (49) and (12)
304  std::vector<Real> rP(_number_of_nodes);
305  std::vector<Real> rM(_number_of_nodes);
306  std::vector<std::vector<Real>> drP(_number_of_nodes); // drP[i][j] = d(rP[i])/d(u[j]). Here j
307  // indexes the j^th node connected to i
308  std::vector<std::vector<Real>> drM(_number_of_nodes);
309  std::vector<std::vector<Real>> drP_dk(
310  _number_of_nodes); // drP_dk[i][j] = d(rP[i])/d(K[i][j]). Here j indexes the j^th node
311  // connected to i
312  std::vector<std::vector<Real>> drM_dk(_number_of_nodes);
313  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
314  {
315  rP[sequential_i] = rPlus(sequential_i, drP[sequential_i], drP_dk[sequential_i]);
316  rM[sequential_i] = rMinus(sequential_i, drM[sequential_i], drM_dk[sequential_i]);
317  }
318 
319  // Calculate KuzminTurek f^{a} matrix
320  // This is the antidiffusive flux
321  // See Eqn (50)
322  std::vector<std::vector<Real>> fa(
323  _number_of_nodes); // fa[sequential_i][j] sequential_j is the j^th connection to sequential_i
324  // The derivatives are a bit complicated.
325  // If i is upwind of j then fa[i][j] depends on all nodes connected to i.
326  // But if i is downwind of j then fa[i][j] depends on all nodes connected to j.
327  std::vector<std::vector<std::map<dof_id_type, Real>>> dfa(
328  _number_of_nodes); // dfa[sequential_i][j][global_k] = d(fa[sequential_i][j])/du[global_k].
329  // Here global_k can be a neighbor to sequential_i or a neighbour to
330  // sequential_j (sequential_j is the j^th connection to sequential_i)
331  std::vector<std::vector<std::vector<Real>>> dFij_dKik(
332  _number_of_nodes); // dFij_dKik[sequential_i][j][k] =
333  // d(fa[sequential_i][j])/d(K[sequential_i][k]). Here j denotes the j^th
334  // connection to sequential_i, while k denotes the k^th connection to
335  // sequential_i
336  std::vector<std::vector<std::vector<Real>>> dFij_dKjk(
337  _number_of_nodes); // dFij_dKjk[sequential_i][j][k] =
338  // d(fa[sequential_i][j])/d(K[sequential_j][k]). Here sequential_j is the
339  // j^th connection to sequential_i, and k denotes the k^th connection to
340  // sequential_j (this will include sequential_i itself)
341  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
342  {
343  const std::vector<dof_id_type> con_i =
345  const unsigned num_con_i = con_i.size();
346  fa[sequential_i].assign(num_con_i, 0.0);
347  dfa[sequential_i].resize(num_con_i);
348  dFij_dKik[sequential_i].resize(num_con_i);
349  dFij_dKjk[sequential_i].resize(num_con_i);
350  for (std::size_t j = 0; j < num_con_i; ++j)
351  {
352  for (const auto & global_k : con_i)
353  dfa[sequential_i][j][global_k] = 0;
354  const dof_id_type global_j = con_i[j];
355  const std::vector<dof_id_type> con_j = _connections.globalConnectionsToGlobalID(global_j);
356  const unsigned num_con_j = con_j.size();
357  for (const auto & global_k : con_j)
358  dfa[sequential_i][j][global_k] = 0;
359  dFij_dKik[sequential_i][j].assign(num_con_i, 0.0);
360  dFij_dKjk[sequential_i][j].assign(num_con_j, 0.0);
361  }
362  }
363 
364  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
365  {
366  const dof_id_type global_i = _connections.globalID(sequential_i);
367  const Real u_i = _u_nodal[sequential_i];
368  const std::vector<dof_id_type> con_i =
370  const std::size_t num_con_i = con_i.size();
371  for (std::size_t j = 0; j < num_con_i; ++j)
372  {
373  const dof_id_type sequential_j = con_i[j];
374  if (sequential_i == sequential_j)
375  continue;
376  const unsigned index_j_to_i =
377  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
378  const Real Lij = lij[sequential_i][j];
379  const Real Lji = lij[sequential_j][index_j_to_i];
380  if (Lji >= Lij) // node i is upwind of node j.
381  {
382  const Real Dij = dij[sequential_i][j];
383  const dof_id_type global_j = _connections.globalID(sequential_j);
384  const Real u_j = _u_nodal[sequential_j];
385  Real prefactor = 0.0;
386  std::vector<Real> dprefactor_du(num_con_i,
387  0.0); // dprefactor_du[j] = d(prefactor)/du[sequential_j]
388  std::vector<Real> dprefactor_dKij(
389  num_con_i, 0.0); // dprefactor_dKij[j] = d(prefactor)/dKij[sequential_i][j].
390  Real dprefactor_dKji = 0; // dprefactor_dKji = d(prefactor)/dKij[sequential_j][index_j_to_i]
391  if (u_i >= u_j)
392  {
393  if (Lji <= rP[sequential_i] * Dij)
394  {
395  prefactor = Lji;
396  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
397  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
398  }
399  else
400  {
401  prefactor = rP[sequential_i] * Dij;
402  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
403  dprefactor_du[ind_j] = drP[sequential_i][ind_j] * Dij;
404  dprefactor_dKij[j] += rP[sequential_i] * dDij_dKij[sequential_i][j];
405  dprefactor_dKji += rP[sequential_i] * dDij_dKji[sequential_i][j];
406  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
407  dprefactor_dKij[ind_j] += drP_dk[sequential_i][ind_j] * Dij;
408  }
409  }
410  else
411  {
412  if (Lji <= rM[sequential_i] * Dij)
413  {
414  prefactor = Lji;
415  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
416  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
417  }
418  else
419  {
420  prefactor = rM[sequential_i] * Dij;
421  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
422  dprefactor_du[ind_j] = drM[sequential_i][ind_j] * Dij;
423  dprefactor_dKij[j] += rM[sequential_i] * dDij_dKij[sequential_i][j];
424  dprefactor_dKji += rM[sequential_i] * dDij_dKji[sequential_i][j];
425  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
426  dprefactor_dKij[ind_j] += drM_dk[sequential_i][ind_j] * Dij;
427  }
428  }
429  fa[sequential_i][j] = prefactor * (u_i - u_j);
430  dfa[sequential_i][j][global_i] = prefactor;
431  dfa[sequential_i][j][global_j] = -prefactor;
432  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
433  {
434  dfa[sequential_i][j][_connections.globalID(con_i[ind_j])] +=
435  dprefactor_du[ind_j] * (u_i - u_j);
436  dFij_dKik[sequential_i][j][ind_j] += dprefactor_dKij[ind_j] * (u_i - u_j);
437  }
438  dFij_dKjk[sequential_i][j][index_j_to_i] += dprefactor_dKji * (u_i - u_j);
439  }
440  }
441  }
442 
443  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
444  {
445  const std::vector<dof_id_type> con_i =
447  const std::size_t num_con_i = con_i.size();
448  for (std::size_t j = 0; j < num_con_i; ++j)
449  {
450  const dof_id_type sequential_j = con_i[j];
451  if (sequential_i == sequential_j)
452  continue;
453  const unsigned index_j_to_i =
454  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
455  if (lij[sequential_j][index_j_to_i] < lij[sequential_i][j]) // node_i is downwind of node_j.
456  {
457  fa[sequential_i][j] = -fa[sequential_j][index_j_to_i];
458  for (const auto & dof_deriv : dfa[sequential_j][index_j_to_i])
459  dfa[sequential_i][j][dof_deriv.first] = -dof_deriv.second;
460  for (std::size_t k = 0; k < num_con_i; ++k)
461  dFij_dKik[sequential_i][j][k] = -dFij_dKjk[sequential_j][index_j_to_i][k];
462  const std::size_t num_con_j =
464  for (std::size_t k = 0; k < num_con_j; ++k)
465  dFij_dKjk[sequential_i][j][k] = -dFij_dKik[sequential_j][index_j_to_i][k];
466  }
467  }
468  }
469 
470  // zero _flux_out and its derivatives
471  _flux_out.assign(_number_of_nodes, 0.0);
472  // The derivatives are a bit complicated.
473  // If i is upwind of a node "j" then _flux_out[i] depends on all nodes connected to i.
474  // But if i is downwind of a node "j" then _flux_out depends on all nodes connected with node
475  // j.
476  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
477  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
478  for (const auto & j : _connections.globalConnectionsToSequentialID(sequential_i))
479  {
480  _dflux_out_du[sequential_i][j] = 0.0;
481  for (const auto & neighbors_j : _connections.globalConnectionsToGlobalID(j))
482  _dflux_out_du[sequential_i][neighbors_j] = 0.0;
483  }
485  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
486  {
487  const std::vector<dof_id_type> con_i =
489  const std::size_t num_con_i = con_i.size();
490  _dflux_out_dKjk[sequential_i].resize(num_con_i);
491  for (std::size_t j = 0; j < num_con_i; ++j)
492  {
493  const dof_id_type sequential_j = con_i[j];
494  std::vector<dof_id_type> con_j =
496  _dflux_out_dKjk[sequential_i][j].assign(con_j.size(), 0.0);
497  }
498  }
499 
500  // Add everything together
501  // See step 3 in Fig 2, noting Eqn (36)
502  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
503  {
504  std::vector<dof_id_type> con_i = _connections.sequentialConnectionsToSequentialID(sequential_i);
505  const size_t num_con_i = con_i.size();
506  const dof_id_type index_i_to_i =
507  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
508  for (std::size_t j = 0; j < num_con_i; ++j)
509  {
510  const dof_id_type sequential_j = con_i[j];
511  const dof_id_type global_j = _connections.globalID(sequential_j);
512  const Real u_j = _u_nodal[sequential_j];
513 
514  // negative sign because residual = -Lu (KT equation (19))
515  _flux_out[sequential_i] -= lij[sequential_i][j] * u_j + fa[sequential_i][j];
516 
517  _dflux_out_du[sequential_i][global_j] -= lij[sequential_i][j];
518  for (const auto & dof_deriv : dfa[sequential_i][j])
519  _dflux_out_du[sequential_i][dof_deriv.first] -= dof_deriv.second;
520 
521  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= 1.0 * u_j; // from the K in L = K + D
522 
523  if (sequential_j == sequential_i)
524  for (dof_id_type k = 0; k < num_con_i; ++k)
525  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dDii_dKij[sequential_i][k] * u_j;
526  else
527  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= dDij_dKij[sequential_i][j] * u_j;
528  for (dof_id_type k = 0; k < num_con_i; ++k)
529  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dFij_dKik[sequential_i][j][k];
530 
531  if (sequential_j == sequential_i)
532  for (unsigned k = 0; k < con_i.size(); ++k)
533  {
534  const unsigned index_k_to_i =
535  _connections.indexOfSequentialConnection(con_i[k], sequential_i);
536  _dflux_out_dKjk[sequential_i][k][index_k_to_i] -= dDii_dKji[sequential_i][k] * u_j;
537  }
538  else
539  {
540  const unsigned index_j_to_i =
541  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
542  _dflux_out_dKjk[sequential_i][j][index_j_to_i] -= dDij_dKji[sequential_i][j] * u_j;
543  }
544  for (unsigned k = 0;
545  k < _connections.sequentialConnectionsToSequentialID(sequential_j).size();
546  ++k)
547  _dflux_out_dKjk[sequential_i][j][k] -= dFij_dKjk[sequential_i][j][k];
548  }
549  }
550 }
dof_id_type globalID(dof_id_type sequential_node_ID) const
Return the global node ID (node number in the mesh) corresponding to the provided sequential node ID...
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real rMinus(dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
Returns the value of R_{i}^{-}, Eqn (49) of KT.
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
const std::vector< dof_id_type > & globalConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (global node IDs) connected to the given sequential node ID. ...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
Real rPlus(dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
Returns the value of R_{i}^{+}, Eqn (49) of KT.
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
const std::vector< dof_id_type > & globalConnectionsToGlobalID(dof_id_type global_node_ID) const
Return all the nodes (global node IDs) connected to the given global node ID.
unsigned indexOfSequentialConnection(dof_id_type sequential_node_ID_from, dof_id_type sequential_node_ID_to) const
Return the index of sequential_node_ID_to in the sequentialConnectionsToSequentialID(sequential_node_...

◆ getdFluxOutdKjk()

const std::vector< std::vector< Real > > & AdvectiveFluxCalculatorBase::getdFluxOutdKjk ( dof_id_type  node_i) const
inherited

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

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

Definition at line 732 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::finalize().

733 {
734  return _dflux_out_dKjk[_connections.sequentialID(node_i)];
735 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ 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 726 of file AdvectiveFluxCalculatorBase.C.

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

727 {
728  return _dflux_out_du[_connections.sequentialID(node_i)];
729 }
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ 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 738 of file AdvectiveFluxCalculatorBase.C.

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

739 {
740  return _flux_out[_connections.sequentialID(node_i)];
741 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ getValence()

unsigned AdvectiveFluxCalculatorBase::getValence ( dof_id_type  node_i) const
inherited

Returns the valence of the global node i Valence is the number of times the node 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_igloal id of i^th node
Returns
valence of the node

Definition at line 744 of file AdvectiveFluxCalculatorBase.C.

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

745 {
746  return _valence[_connections.sequentialID(node_i)];
747 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< unsigned > _valence
_valence[i] = 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 sequential node i is encountered

◆ initialize()

void AdvectiveFluxCalculatorBase::initialize ( )
overridevirtualinherited

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 152 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::initialize().

153 {
154  // Zero _kij and falsify _u_nodal_computed_by_thread ready for building in execute() and
155  // finalize()
157  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
158  _kij[sequential_i].assign(_connections.sequentialConnectionsToSequentialID(sequential_i).size(),
159  0.0);
160 }
std::vector< 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 ...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< std::vector< 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 625 of file AdvectiveFluxCalculatorBase.C.

Referenced by AdvectiveFluxCalculatorBase::rMinus(), and AdvectiveFluxCalculatorBase::rPlus().

626 {
627  limited = 0.0;
628  dlimited_db = 0.0;
630  return;
631 
632  if ((a >= 0.0 && b <= 0.0) || (a <= 0.0 && b >= 0.0))
633  return;
634  const Real s = (a > 0.0 ? 1.0 : -1.0);
635 
636  const Real lal = std::abs(a);
637  const Real lbl = std::abs(b);
638  const Real dlbl = (b >= 0.0 ? 1.0 : -1.0); // d(lbl)/db
639  switch (_flux_limiter_type)
640  {
642  {
643  if (lal <= lbl)
644  {
645  limited = s * lal;
646  dlimited_db = 0.0;
647  }
648  else
649  {
650  limited = s * lbl;
651  dlimited_db = s * dlbl;
652  }
653  return;
654  }
656  {
657  limited = s * 2 * lal * lbl / (lal + lbl);
658  dlimited_db = s * 2 * lal * (dlbl / (lal + lbl) - lbl * dlbl / std::pow(lal + lbl, 2));
659  return;
660  }
662  {
663  const Real av = 0.5 * std::abs(a + b);
664  if (2 * lal <= av && lal <= lbl)
665  {
666  // 2 * lal is the smallest
667  limited = s * 2.0 * lal;
668  dlimited_db = 0.0;
669  }
670  else if (2 * lbl <= av && lbl <= lal)
671  {
672  // 2 * lbl is the smallest
673  limited = s * 2.0 * lbl;
674  dlimited_db = s * 2.0 * dlbl;
675  }
676  else
677  {
678  // av is the smallest
679  limited = s * av;
680  // if (a>0 and b>0) then d(av)/db = 0.5 = 0.5 * dlbl
681  // if (a<0 and b<0) then d(av)/db = -0.5 = 0.5 * dlbl
682  // if a and b have different sign then limited=0, above
683  dlimited_db = s * 0.5 * dlbl;
684  }
685  return;
686  }
688  {
689  const Real term1 = std::min(2.0 * lal, lbl);
690  const Real term2 = std::min(lal, 2.0 * lbl);
691  if (term1 >= term2)
692  {
693  if (2.0 * lal <= lbl)
694  {
695  limited = s * 2 * lal;
696  dlimited_db = 0.0;
697  }
698  else
699  {
700  limited = s * lbl;
701  dlimited_db = s * dlbl;
702  }
703  }
704  else
705  {
706  if (lal <= 2.0 * lbl)
707  {
708  limited = s * lal;
709  dlimited_db = 0.0;
710  }
711  else
712  {
713  limited = s * 2.0 * lbl;
714  dlimited_db = s * 2.0 * dlbl;
715  }
716  }
717  return;
718  }
719  default:
720  return;
721  }
722 }
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ meshChanged()

void AdvectiveFluxCalculatorBase::meshChanged ( )
overridevirtualinherited

Definition at line 143 of file AdvectiveFluxCalculatorBase.C.

144 {
145  ElementUserObject::meshChanged();
146 
147  // Signal that _kij, _valence, etc need to be rebuilt
148  _resizing_needed = true;
149 }
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  sequential_i,
const PQPlusMinusEnum  pq_plus_minus,
std::vector< Real > &  derivs,
std::vector< 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
sequential_isequential 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[sequential_j]). Here sequential_j is the j^th connection to sequential_i
dpq_dk[out]dpq_dk[j] = d(result)/d(K[node_i][j]). Here j indexes a connection to sequential_i. Recall that d(result)/d(K[l][m]) are zero unless l=sequential_i

Definition at line 759 of file AdvectiveFluxCalculatorBase.C.

Referenced by AdvectiveFluxCalculatorBase::rMinus(), and AdvectiveFluxCalculatorBase::rPlus().

763 {
764  // Find the value of u at sequential_i
765  const Real u_i = _u_nodal[sequential_i];
766 
767  // Connections to sequential_i
768  const std::vector<dof_id_type> con_i =
770  const std::size_t num_con = con_i.size();
771  // The neighbor number of sequential_i to sequential_i
772  const unsigned i_index_i = _connections.indexOfSequentialConnection(sequential_i, sequential_i);
773 
774  // Initialize the results
775  Real result = 0.0;
776  derivs.assign(num_con, 0.0);
777  dpqdk.assign(num_con, 0.0);
778 
779  // Sum over all nodes connected with node_i.
780  for (std::size_t j = 0; j < num_con; ++j)
781  {
782  const dof_id_type sequential_j = con_i[j];
783  if (sequential_j == sequential_i)
784  continue;
785  const Real kentry = _kij[sequential_i][j];
786 
787  // Find the value of u at node_j
788  const Real u_j = _u_nodal[sequential_j];
789  const Real ujminusi = u_j - u_i;
790 
791  // Evaluate the i-j contribution to the result
792  switch (pq_plus_minus)
793  {
795  {
796  if (ujminusi < 0.0 && kentry < 0.0)
797  {
798  result += kentry * ujminusi;
799  derivs[j] += kentry;
800  derivs[i_index_i] -= kentry;
801  dpqdk[j] += ujminusi;
802  }
803  break;
804  }
806  {
807  if (ujminusi > 0.0 && kentry < 0.0)
808  {
809  result += kentry * ujminusi;
810  derivs[j] += kentry;
811  derivs[i_index_i] -= kentry;
812  dpqdk[j] += ujminusi;
813  }
814  break;
815  }
817  {
818  if (ujminusi > 0.0 && kentry > 0.0)
819  {
820  result += kentry * ujminusi;
821  derivs[j] += kentry;
822  derivs[i_index_i] -= kentry;
823  dpqdk[j] += ujminusi;
824  }
825  break;
826  }
828  {
829  if (ujminusi < 0.0 && kentry > 0.0)
830  {
831  result += kentry * ujminusi;
832  derivs[j] += kentry;
833  derivs[i_index_i] -= kentry;
834  dpqdk[j] += ujminusi;
835  }
836  break;
837  }
838  }
839  }
840 
841  return result;
842 }
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
unsigned indexOfSequentialConnection(dof_id_type sequential_node_ID_from, dof_id_type sequential_node_ID_to) const
Return the index of sequential_node_ID_to in the sequentialConnectionsToSequentialID(sequential_node_...

◆ rMinus()

Real AdvectiveFluxCalculatorBase::rMinus ( dof_id_type  sequential_i,
std::vector< Real > &  dlimited_du,
std::vector< Real > &  dlimited_dk 
) const
protectedinherited

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

Parameters
sequential_iSequential nodal ID
dlimited_du[out]dlimited_du[j] = d(R_{sequential_i}^{-})/du[sequential_j]. Here sequential_j is the j^th connection to sequential_i
dlimited_dk[out]dlimited_dk[j] = d(R_{sequential_i}^{-})/d(K[sequential_i][j]). Note Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 589 of file AdvectiveFluxCalculatorBase.C.

Referenced by AdvectiveFluxCalculatorBase::finalize().

592 {
593  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
594  dlimited_du.assign(num_con, 0.0);
595  dlimited_dk.assign(num_con, 0.0);
597  return 0.0;
598  std::vector<Real> dp_du;
599  std::vector<Real> dp_dk;
600  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PMinus, dp_du, dp_dk);
601  if (p == 0.0)
602  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
603  return 1.0;
604  std::vector<Real> dq_du;
605  std::vector<Real> dq_dk;
606  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QMinus, dq_du, dq_dk);
607 
608  const Real r = q / p;
609  Real limited;
610  Real dlimited_dr;
611  limitFlux(1.0, r, limited, dlimited_dr);
612 
613  const Real p2 = std::pow(p, 2);
614  for (std::size_t j = 0; j < num_con; ++j)
615  {
616  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
617  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
618  dlimited_du[j] = dlimited_dr * dr_du;
619  dlimited_dk[j] = dlimited_dr * dr_dk;
620  }
621  return limited;
622 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real PQPlusMinus(dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< 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
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ rPlus()

Real AdvectiveFluxCalculatorBase::rPlus ( dof_id_type  sequential_i,
std::vector< Real > &  dlimited_du,
std::vector< Real > &  dlimited_dk 
) const
protectedinherited

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

Parameters
node_inodal id
dlimited_du[out]dlimited_du[j] = d(R_{sequential_i}^{+})/du[sequential_j]. Here sequential_j is the j^th connection to sequential_i
dlimited_dk[out]dlimited_dk[j] = d(R_{sequential_i}^{+})/d(K[sequential_i][j]). Note Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 553 of file AdvectiveFluxCalculatorBase.C.

Referenced by AdvectiveFluxCalculatorBase::finalize().

556 {
557  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
558  dlimited_du.assign(num_con, 0.0);
559  dlimited_dk.assign(num_con, 0.0);
561  return 0.0;
562  std::vector<Real> dp_du;
563  std::vector<Real> dp_dk;
564  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PPlus, dp_du, dp_dk);
565  if (p == 0.0)
566  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
567  return 1.0;
568  std::vector<Real> dq_du;
569  std::vector<Real> dq_dk;
570  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QPlus, dq_du, dq_dk);
571 
572  const Real r = q / p;
573  Real limited;
574  Real dlimited_dr;
575  limitFlux(1.0, r, limited, dlimited_dr);
576 
577  const Real p2 = std::pow(p, 2);
578  for (std::size_t j = 0; j < num_con; ++j)
579  {
580  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
581  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
582  dlimited_du[j] = dlimited_dr * dr_du;
583  dlimited_dk[j] = dlimited_dr * dr_dk;
584  }
585  return limited;
586 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real PQPlusMinus(dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< 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
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ threadJoin()

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

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 196 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::threadJoin().

197 {
198  const AdvectiveFluxCalculatorBase & afc = static_cast<const AdvectiveFluxCalculatorBase &>(uo);
199  // add the values of _kij computed by different threads
200  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
201  {
202  const std::size_t num_con_i =
204  for (std::size_t j = 0; j < num_con_i; ++j)
205  _kij[sequential_i][j] += afc._kij[sequential_i][j];
206  }
207 
208  // gather the values of _u_nodal computed by different threads
209  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
210  if (!_u_nodal_computed_by_thread[sequential_i] && afc._u_nodal_computed_by_thread[sequential_i])
211  _u_nodal[sequential_i] = afc._u_nodal[sequential_i];
212 }
std::vector< 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 ...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
Base class to compute Advective fluxes.
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ timestepSetup()

void AdvectiveFluxCalculatorBase::timestepSetup ( )
overridevirtualinherited

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 57 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

58 {
59  // If needed, size and initialize quantities appropriately, and compute _valence
60  if (_resizing_needed)
61  {
63  // QUERY: does the following properly account for multiple processors, threading, and other
64  // things i haven't thought about?
65 
66  /*
67  * Populate _connections for all nodes that can be seen by this processor and on relevant blocks
68  *
69  * MULTIPROC NOTE: this must loop over local elements and 2 layers of ghosted elements.
70  * The Kernel will only loop over local elements, so will only use _kij, etc, for
71  * linked node-node pairs that appear in the local elements. Nevertheless, we
72  * need to build _kij, etc, for the nodes in the ghosted elements in order to simplify
73  * Jacobian computations
74  */
75  for (const auto & elem : _subproblem.mesh().getMesh().active_element_ptr_range())
76  if (this->hasBlocks(elem->subdomain_id()))
77  for (unsigned i = 0; i < elem->n_nodes(); ++i)
78  _connections.addGlobalNode(elem->node_id(i));
80  for (const auto & elem : _subproblem.mesh().getMesh().active_element_ptr_range())
81  if (this->hasBlocks(elem->subdomain_id()))
82  for (unsigned i = 0; i < elem->n_nodes(); ++i)
83  for (unsigned j = 0; j < elem->n_nodes(); ++j)
84  _connections.addConnection(elem->node_id(i), elem->node_id(j));
86 
88 
89  // initialize _kij
90  _kij.resize(_number_of_nodes);
91  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
92  _kij[sequential_i].assign(
93  _connections.sequentialConnectionsToSequentialID(sequential_i).size(), 0.0);
94 
95  /*
96  * Build _valence[i], which is the number of times the sequential node i is encountered when
97  * looping over the entire mesh (and on relevant blocks)
98  *
99  * MULTIPROC NOTE: this must loop over local elements and >=1 layer of ghosted elements.
100  * The Kernel will only loop over local elements, so will only use _valence for
101  * linked node-node pairs that appear in the local elements. But other processors will
102  * loop over neighboring elements, so avoid multiple counting of the residual and Jacobian
103  * this processor must record how many times each node-node link of its local elements
104  * appears in the local elements and >=1 layer, and pass that info to the Kernel
105  */
106  _valence.assign(_number_of_nodes, 0);
107  for (const auto & elem : _fe_problem.getEvaluableElementRange())
108  if (this->hasBlocks(elem->subdomain_id()))
109  for (unsigned i = 0; i < elem->n_nodes(); ++i)
110  {
111  const dof_id_type node_i = elem->node_id(i);
112  const dof_id_type sequential_i = _connections.sequentialID(node_i);
113  _valence[sequential_i] += 1;
114  }
115 
116  _u_nodal.assign(_number_of_nodes, 0.0);
118  _flux_out.assign(_number_of_nodes, 0.0);
119  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
120  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
121  zeroedConnection(_dflux_out_du[sequential_i], _connections.globalID(sequential_i));
123  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
124  {
125  const std::vector<dof_id_type> con_i =
127  const std::size_t num_con_i = con_i.size();
128  _dflux_out_dKjk[sequential_i].resize(num_con_i);
129  for (std::size_t j = 0; j < con_i.size(); ++j)
130  {
131  const dof_id_type sequential_j = con_i[j];
132  const std::size_t num_con_j =
134  _dflux_out_dKjk[sequential_i][j].assign(num_con_j, 0.0);
135  }
136  }
137 
138  _resizing_needed = false;
139  }
140 }
void clear()
clear all data in readiness for adding global nodes and connections
std::vector< 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 ...
dof_id_type globalID(dof_id_type sequential_node_ID) const
Return the global node ID (node number in the mesh) corresponding to the provided sequential node ID...
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
void addConnection(dof_id_type global_node_from, dof_id_type global_node_to)
Specifies that global_node_to is connected to global_node_from.
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
void finalizeAddingConnections()
Signal that all global node IDs have been added to the internal data structures.
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::size_t numNodes() const
number of nodes known by this class
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
void finalizeAddingGlobalNodes()
Signal that all global node IDs have been added to the internal data structures.
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::vector< unsigned > _valence
_valence[i] = 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 sequential node i is encountered
void addGlobalNode(dof_id_type global_node_ID)
Add the given global_node_ID to the internal data structures If the global node ID has already been a...

◆ 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 750 of file AdvectiveFluxCalculatorBase.C.

Referenced by AdvectiveFluxCalculatorBase::timestepSetup().

752 {
753  the_map.clear();
754  for (const auto & node_j : _connections.globalConnectionsToGlobalID(node_i))
755  the_map[node_j] = 0.0;
756 }
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
const std::vector< dof_id_type > & globalConnectionsToGlobalID(dof_id_type global_node_ID) const
Return all the nodes (global node IDs) connected to the given global node ID.

Member Data Documentation

◆ _connections

PorousFlowConnectedNodes AdvectiveFluxCalculatorBase::_connections
protectedinherited

◆ _dflux_out_dKjk

std::vector<std::vector<std::vector<Real> > > AdvectiveFluxCalculatorBase::_dflux_out_dKjk
protectedinherited

_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).

Here sequential_i is a sequential node number according to the _connections object, and j represents the j^th node connected to i according to the _connections object, and k represents the k^th node connected to j according to the _connections object. Here j must be connected to i (this does include (the sequential version of j) == i), and k must be connected to j (this does include (the sequential version of k) = i and (the sequential version of k) == (sequential version of j)))

Definition at line 175 of file AdvectiveFluxCalculatorBase.h.

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

◆ _dflux_out_du

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

_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).

Here i is a sequential node number according to the _connections object, and j (global ID) must be connected to i, or to a node that is connected to i.

Definition at line 168 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::vector<Real> AdvectiveFluxCalculatorBase::_flux_out
protectedinherited

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

Definition at line 164 of file AdvectiveFluxCalculatorBase.h.

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

◆ _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::vector<std::vector<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. Here i is a sequential node numbers according to the _connections object, and j represents the j^th node connected to i according to the _connections object.

Definition at line 161 of file AdvectiveFluxCalculatorBase.h.

Referenced by AdvectiveFluxCalculatorBase::executeOnElement(), AdvectiveFluxCalculatorBase::finalize(), AdvectiveFluxCalculatorBase::initialize(), AdvectiveFluxCalculatorBase::PQPlusMinus(), AdvectiveFluxCalculatorBase::threadJoin(), and AdvectiveFluxCalculatorBase::timestepSetup().

◆ _number_of_nodes

std::size_t AdvectiveFluxCalculatorBase::_number_of_nodes
protectedinherited

◆ _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 115 of file AdvectiveFluxCalculatorBase.h.

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

◆ _u_at_nodes

const VariableValue& AdvectiveFluxCalculatorConstantVelocity::_u_at_nodes
protected

the nodal values of u

Definition at line 37 of file AdvectiveFluxCalculatorConstantVelocity.h.

Referenced by computeU().

◆ _u_nodal

std::vector<Real> AdvectiveFluxCalculatorBase::_u_nodal
protectedinherited

◆ _u_nodal_computed_by_thread

std::vector<bool> AdvectiveFluxCalculatorBase::_u_nodal_computed_by_thread
protectedinherited

_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 186 of file AdvectiveFluxCalculatorBase.h.

Referenced by AdvectiveFluxCalculatorBase::execute(), AdvectiveFluxCalculatorBase::initialize(), AdvectiveFluxCalculatorBase::threadJoin(), and AdvectiveFluxCalculatorBase::timestepSetup().

◆ _valence

std::vector<unsigned> AdvectiveFluxCalculatorBase::_valence
protectedinherited

_valence[i] = 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 sequential node i is encountered

Definition at line 180 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: