www.mooseframework.org
AdvectiveFluxCalculatorBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "Conversion.h" // for stringify
12 
13 #include "libmesh/string_to_enum.h"
14 #include "libmesh/mesh_tools.h"
15 #include "libmesh/parallel_sync.h"
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<ElementUserObject>();
22  params.addClassDescription(
23  "Base class to compute K_ij (a measure of advective flux from node i to node j) "
24  "and R+ and R- (which quantify amount of antidiffusion to add) in the "
25  "Kuzmin-Turek FEM-TVD multidimensional scheme");
26  MooseEnum flux_limiter_type("MinMod VanLeer MC superbee None", "VanLeer");
27  params.addParam<MooseEnum>("flux_limiter_type",
28  flux_limiter_type,
29  "Type of flux limiter to use. 'None' means that no antidiffusion "
30  "will be added in the Kuzmin-Turek scheme");
31  params.addRangeCheckedParam<Real>(
32  "allowable_MB_wastage",
33  5.0,
34  "allowable_MB_wastage > 0.0",
35  "This object will issue a memory warning if the internal node-numbering data structure "
36  "wastes more than allowable_MB_wastage megabytes. This data structure uses sequential "
37  "node-numbering which is optimized for speed rather than memory efficiency");
38 
39  params.addRelationshipManager("ElementPointNeighborLayers",
40  Moose::RelationshipManagerType::GEOMETRIC |
41  Moose::RelationshipManagerType::ALGEBRAIC,
42  [](const InputParameters &, InputParameters & rm_params) {
43  rm_params.set<unsigned short>("layers") = 2;
44  });
45 
46  params.set<ExecFlagEnum>("execute_on", true) = {EXEC_LINEAR};
47  return params;
48 }
49 
51  : ElementUserObject(parameters),
52  _resizing_needed(true),
53  _flux_limiter_type(getParam<MooseEnum>("flux_limiter_type").getEnum<FluxLimiterTypeEnum>()),
54  _kij(),
55  _flux_out(),
56  _dflux_out_du(),
57  _dflux_out_dKjk(),
58  _valence(),
59  _u_nodal(),
60  _u_nodal_computed_by_thread(),
61  _connections(),
62  _number_of_nodes(0),
63  _my_pid(processor_id()),
64  _nodes_to_receive(),
65  _nodes_to_send(),
66  _pairs_to_receive(),
67  _pairs_to_send(),
68  _allowable_MB_wastage(getParam<Real>("allowable_MB_wastage"))
69 {
70  if (!_execute_enum.contains(EXEC_LINEAR))
71  paramError(
72  "execute_on",
73  "The AdvectiveFluxCalculator UserObject " + name() +
74  " execute_on parameter must include, at least, 'linear'. This is to ensure that "
75  "this UserObject computes all necessary quantities just before the Kernels evaluate "
76  "their Residuals");
77 }
78 
79 void
81 {
82  // If needed, size and initialize quantities appropriately, and compute _valence
83  if (_resizing_needed)
84  {
85  /*
86  * Populate _connections for all nodes that can be seen by this processor and on relevant
87  * blocks
88  *
89  * MULTIPROC NOTE: this must loop over local elements and 2 layers of ghosted elements.
90  * The Kernel will only loop over local elements, so will only use _kij, etc, for
91  * linked node-node pairs that appear in the local elements. Nevertheless, we
92  * need to build _kij, etc, for the nodes in the ghosted elements in order to simplify
93  * Jacobian computations
94  */
96  for (const auto & elem : _fe_problem.getEvaluableElementRange())
97  if (this->hasBlocks(elem->subdomain_id()))
98  for (unsigned i = 0; i < elem->n_nodes(); ++i)
99  _connections.addGlobalNode(elem->node_id(i));
101  for (const auto & elem : _fe_problem.getEvaluableElementRange())
102  if (this->hasBlocks(elem->subdomain_id()))
103  for (unsigned i = 0; i < elem->n_nodes(); ++i)
104  for (unsigned j = 0; j < elem->n_nodes(); ++j)
105  _connections.addConnection(elem->node_id(i), elem->node_id(j));
107 
109 
110  Real mb_wasted = (_connections.sizeSequential() - _number_of_nodes) * 4.0 / 1048576;
111  gatherMax(mb_wasted);
112  if (mb_wasted > _allowable_MB_wastage)
113  mooseWarning(
114  "In at least one processor, the sequential node-numbering internal data structure used "
115  "in " +
116  name() + " is using memory inefficiently.\nThe memory wasted is " +
117  Moose::stringify(mb_wasted) +
118  " megabytes.\n The node-numbering data structure has been optimized for speed at the "
119  "expense of memory, but that may not be an appropriate optimization for your case, "
120  "because the node numbering is not close to sequential in your case.\n");
121 
122  // initialize _kij
123  _kij.resize(_number_of_nodes);
124  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
125  _kij[sequential_i].assign(
126  _connections.sequentialConnectionsToSequentialID(sequential_i).size(), 0.0);
127 
128  /*
129  * Build _valence[i], which is the number of times the sequential node i is encountered when
130  * looping over the entire mesh (and on relevant blocks)
131  *
132  * MULTIPROC NOTE: this must loop over local elements and >=1 layer of ghosted elements.
133  * The Kernel will only loop over local elements, so will only use _valence for
134  * linked node-node pairs that appear in the local elements. But other processors will
135  * loop over neighboring elements, so avoid multiple counting of the residual and Jacobian
136  * this processor must record how many times each node-node link of its local elements
137  * appears in the local elements and >=1 layer, and pass that info to the Kernel
138  */
139  _valence.assign(_number_of_nodes, 0);
140  for (const auto & elem : _fe_problem.getEvaluableElementRange())
141  if (this->hasBlocks(elem->subdomain_id()))
142  for (unsigned i = 0; i < elem->n_nodes(); ++i)
143  {
144  const dof_id_type node_i = elem->node_id(i);
145  const dof_id_type sequential_i = _connections.sequentialID(node_i);
146  _valence[sequential_i] += 1;
147  }
148 
149  _u_nodal.assign(_number_of_nodes, 0.0);
151  _flux_out.assign(_number_of_nodes, 0.0);
152  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
153  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
154  zeroedConnection(_dflux_out_du[sequential_i], _connections.globalID(sequential_i));
156  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
157  {
158  const std::vector<dof_id_type> con_i =
160  const std::size_t num_con_i = con_i.size();
161  _dflux_out_dKjk[sequential_i].resize(num_con_i);
162  for (std::size_t j = 0; j < con_i.size(); ++j)
163  {
164  const dof_id_type sequential_j = con_i[j];
165  const std::size_t num_con_j =
167  _dflux_out_dKjk[sequential_i][j].assign(num_con_j, 0.0);
168  }
169  }
170 
171  if (_app.n_processors() > 1)
172  buildCommLists();
173 
174  _resizing_needed = false;
175  }
176 }
177 
178 void
180 {
181  ElementUserObject::meshChanged();
182 
183  // Signal that _kij, _valence, etc need to be rebuilt
184  _resizing_needed = true;
185 }
186 
187 void
189 {
190  // Zero _kij and falsify _u_nodal_computed_by_thread ready for building in execute() and
191  // finalize()
193  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
194  _kij[sequential_i].assign(_connections.sequentialConnectionsToSequentialID(sequential_i).size(),
195  0.0);
196 }
197 
198 void
200 {
201  // compute _kij contributions from this element that is local to this processor
202  // and record _u_nodal
203  for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
204  {
205  const dof_id_type node_i = _current_elem->node_id(i);
206  const dof_id_type sequential_i = _connections.sequentialID(node_i);
207  if (!_u_nodal_computed_by_thread[sequential_i])
208  {
209  _u_nodal[sequential_i] = computeU(i);
210  _u_nodal_computed_by_thread[sequential_i] = true;
211  }
212  for (unsigned j = 0; j < _current_elem->n_nodes(); ++j)
213  {
214  const dof_id_type node_j = _current_elem->node_id(j);
215  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
216  executeOnElement(node_i, node_j, i, j, qp);
217  }
218  }
219 }
220 
221 void
223  dof_id_type global_i, dof_id_type global_j, unsigned local_i, unsigned local_j, unsigned qp)
224 {
225  // KT Eqn (20)
226  const dof_id_type sequential_i = _connections.sequentialID(global_i);
227  const unsigned index_i_to_j = _connections.indexOfGlobalConnection(global_i, global_j);
228  _kij[sequential_i][index_i_to_j] += _JxW[qp] * _coord[qp] * computeVelocity(local_i, local_j, qp);
229 }
230 
231 void
233 {
234  const AdvectiveFluxCalculatorBase & afc = static_cast<const AdvectiveFluxCalculatorBase &>(uo);
235  // add the values of _kij computed by different threads
236  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
237  {
238  const std::size_t num_con_i =
240  for (std::size_t j = 0; j < num_con_i; ++j)
241  _kij[sequential_i][j] += afc._kij[sequential_i][j];
242  }
243 
244  // gather the values of _u_nodal computed by different threads
245  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
246  if (!_u_nodal_computed_by_thread[sequential_i] && afc._u_nodal_computed_by_thread[sequential_i])
247  _u_nodal[sequential_i] = afc._u_nodal[sequential_i];
248 }
249 
250 void
252 {
253  // Overall: ensure _kij is fully built, then compute Kuzmin-Turek D, L, f^a and
254  // relevant Jacobian information, and then the relevant quantities into _flux_out and
255  // _dflux_out_du, _dflux_out_dKjk
256 
257  if (_app.n_processors() > 1)
259 
260  // Calculate KuzminTurek D matrix
261  // See Eqn (32)
262  // This adds artificial diffusion, which eliminates any spurious oscillations
263  // The idea is that D will remove all negative off-diagonal elements when it is added to K
264  // This is identical to full upwinding
265  std::vector<std::vector<Real>> dij(_number_of_nodes);
266  std::vector<std::vector<Real>> dDij_dKij(
267  _number_of_nodes); // dDij_dKij[i][j] = d(D[i][j])/d(K[i][j]) for i!=j
268  std::vector<std::vector<Real>> dDij_dKji(
269  _number_of_nodes); // dDij_dKji[i][j] = d(D[i][j])/d(K[j][i]) for i!=j
270  std::vector<std::vector<Real>> dDii_dKij(
271  _number_of_nodes); // dDii_dKij[i][j] = d(D[i][i])/d(K[i][j])
272  std::vector<std::vector<Real>> dDii_dKji(
273  _number_of_nodes); // dDii_dKji[i][j] = d(D[i][i])/d(K[j][i])
274  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
275  {
276  const std::vector<dof_id_type> con_i =
278  const std::size_t num_con_i = con_i.size();
279  dij[sequential_i].assign(num_con_i, 0.0);
280  dDij_dKij[sequential_i].assign(num_con_i, 0.0);
281  dDij_dKji[sequential_i].assign(num_con_i, 0.0);
282  dDii_dKij[sequential_i].assign(num_con_i, 0.0);
283  dDii_dKji[sequential_i].assign(num_con_i, 0.0);
284  const unsigned index_i_to_i =
285  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
286  for (std::size_t j = 0; j < num_con_i; ++j)
287  {
288  const dof_id_type sequential_j = con_i[j];
289  if (sequential_i == sequential_j)
290  continue;
291  const unsigned index_j_to_i =
292  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
293  const Real kij = _kij[sequential_i][j];
294  const Real kji = _kij[sequential_j][index_j_to_i];
295  if ((kij <= kji) && (kij < 0.0))
296  {
297  dij[sequential_i][j] = -kij;
298  dDij_dKij[sequential_i][j] = -1.0;
299  dDii_dKij[sequential_i][j] += 1.0;
300  }
301  else if ((kji <= kij) && (kji < 0.0))
302  {
303  dij[sequential_i][j] = -kji;
304  dDij_dKji[sequential_i][j] = -1.0;
305  dDii_dKji[sequential_i][j] += 1.0;
306  }
307  else
308  dij[sequential_i][j] = 0.0;
309  dij[sequential_i][index_i_to_i] -= dij[sequential_i][j];
310  }
311  }
312 
313  // Calculate KuzminTurek L matrix
314  // See Fig 2: L = K + D
315  std::vector<std::vector<Real>> lij(_number_of_nodes);
316  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
317  {
318  const std::vector<dof_id_type> con_i =
320  const std::size_t num_con_i = con_i.size();
321  lij[sequential_i].assign(num_con_i, 0.0);
322  for (std::size_t j = 0; j < num_con_i; ++j)
323  lij[sequential_i][j] = _kij[sequential_i][j] + dij[sequential_i][j];
324  }
325 
326  // Compute KuzminTurek R matrices
327  // See Eqns (49) and (12)
328  std::vector<Real> rP(_number_of_nodes);
329  std::vector<Real> rM(_number_of_nodes);
330  std::vector<std::vector<Real>> drP(_number_of_nodes); // drP[i][j] = d(rP[i])/d(u[j]). Here j
331  // indexes the j^th node connected to i
332  std::vector<std::vector<Real>> drM(_number_of_nodes);
333  std::vector<std::vector<Real>> drP_dk(
334  _number_of_nodes); // drP_dk[i][j] = d(rP[i])/d(K[i][j]). Here j indexes the j^th node
335  // connected to i
336  std::vector<std::vector<Real>> drM_dk(_number_of_nodes);
337  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
338  {
339  rP[sequential_i] = rPlus(sequential_i, drP[sequential_i], drP_dk[sequential_i]);
340  rM[sequential_i] = rMinus(sequential_i, drM[sequential_i], drM_dk[sequential_i]);
341  }
342 
343  // Calculate KuzminTurek f^{a} matrix
344  // This is the antidiffusive flux
345  // See Eqn (50)
346  std::vector<std::vector<Real>> fa(_number_of_nodes); // fa[sequential_i][j] sequential_j is the
347  // j^th connection to sequential_i
348  // The derivatives are a bit complicated.
349  // If i is upwind of j then fa[i][j] depends on all nodes connected to i.
350  // But if i is downwind of j then fa[i][j] depends on all nodes connected to j.
351  std::vector<std::vector<std::map<dof_id_type, Real>>> dfa(
352  _number_of_nodes); // dfa[sequential_i][j][global_k] = d(fa[sequential_i][j])/du[global_k].
353  // Here global_k can be a neighbor to sequential_i or a neighbour to
354  // sequential_j (sequential_j is the j^th connection to sequential_i)
355  std::vector<std::vector<std::vector<Real>>> dFij_dKik(
356  _number_of_nodes); // dFij_dKik[sequential_i][j][k] =
357  // d(fa[sequential_i][j])/d(K[sequential_i][k]). Here j denotes the j^th
358  // connection to sequential_i, while k denotes the k^th connection to
359  // sequential_i
360  std::vector<std::vector<std::vector<Real>>> dFij_dKjk(
361  _number_of_nodes); // dFij_dKjk[sequential_i][j][k] =
362  // d(fa[sequential_i][j])/d(K[sequential_j][k]). Here sequential_j is
363  // the j^th connection to sequential_i, and k denotes the k^th connection
364  // to sequential_j (this will include sequential_i itself)
365  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
366  {
367  const std::vector<dof_id_type> con_i =
369  const unsigned num_con_i = con_i.size();
370  fa[sequential_i].assign(num_con_i, 0.0);
371  dfa[sequential_i].resize(num_con_i);
372  dFij_dKik[sequential_i].resize(num_con_i);
373  dFij_dKjk[sequential_i].resize(num_con_i);
374  for (std::size_t j = 0; j < num_con_i; ++j)
375  {
376  for (const auto & global_k : con_i)
377  dfa[sequential_i][j][global_k] = 0;
378  const dof_id_type global_j = con_i[j];
379  const std::vector<dof_id_type> con_j = _connections.globalConnectionsToGlobalID(global_j);
380  const unsigned num_con_j = con_j.size();
381  for (const auto & global_k : con_j)
382  dfa[sequential_i][j][global_k] = 0;
383  dFij_dKik[sequential_i][j].assign(num_con_i, 0.0);
384  dFij_dKjk[sequential_i][j].assign(num_con_j, 0.0);
385  }
386  }
387 
388  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
389  {
390  const dof_id_type global_i = _connections.globalID(sequential_i);
391  const Real u_i = _u_nodal[sequential_i];
392  const std::vector<dof_id_type> con_i =
394  const std::size_t num_con_i = con_i.size();
395  for (std::size_t j = 0; j < num_con_i; ++j)
396  {
397  const dof_id_type sequential_j = con_i[j];
398  if (sequential_i == sequential_j)
399  continue;
400  const unsigned index_j_to_i =
401  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
402  const Real Lij = lij[sequential_i][j];
403  const Real Lji = lij[sequential_j][index_j_to_i];
404  if (Lji >= Lij) // node i is upwind of node j.
405  {
406  const Real Dij = dij[sequential_i][j];
407  const dof_id_type global_j = _connections.globalID(sequential_j);
408  const Real u_j = _u_nodal[sequential_j];
409  Real prefactor = 0.0;
410  std::vector<Real> dprefactor_du(num_con_i,
411  0.0); // dprefactor_du[j] = d(prefactor)/du[sequential_j]
412  std::vector<Real> dprefactor_dKij(
413  num_con_i, 0.0); // dprefactor_dKij[j] = d(prefactor)/dKij[sequential_i][j].
414  Real dprefactor_dKji = 0; // dprefactor_dKji = d(prefactor)/dKij[sequential_j][index_j_to_i]
415  if (u_i >= u_j)
416  {
417  if (Lji <= rP[sequential_i] * Dij)
418  {
419  prefactor = Lji;
420  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
421  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
422  }
423  else
424  {
425  prefactor = rP[sequential_i] * Dij;
426  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
427  dprefactor_du[ind_j] = drP[sequential_i][ind_j] * Dij;
428  dprefactor_dKij[j] += rP[sequential_i] * dDij_dKij[sequential_i][j];
429  dprefactor_dKji += rP[sequential_i] * dDij_dKji[sequential_i][j];
430  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
431  dprefactor_dKij[ind_j] += drP_dk[sequential_i][ind_j] * Dij;
432  }
433  }
434  else
435  {
436  if (Lji <= rM[sequential_i] * Dij)
437  {
438  prefactor = Lji;
439  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
440  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
441  }
442  else
443  {
444  prefactor = rM[sequential_i] * Dij;
445  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
446  dprefactor_du[ind_j] = drM[sequential_i][ind_j] * Dij;
447  dprefactor_dKij[j] += rM[sequential_i] * dDij_dKij[sequential_i][j];
448  dprefactor_dKji += rM[sequential_i] * dDij_dKji[sequential_i][j];
449  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
450  dprefactor_dKij[ind_j] += drM_dk[sequential_i][ind_j] * Dij;
451  }
452  }
453  fa[sequential_i][j] = prefactor * (u_i - u_j);
454  dfa[sequential_i][j][global_i] = prefactor;
455  dfa[sequential_i][j][global_j] = -prefactor;
456  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
457  {
458  dfa[sequential_i][j][_connections.globalID(con_i[ind_j])] +=
459  dprefactor_du[ind_j] * (u_i - u_j);
460  dFij_dKik[sequential_i][j][ind_j] += dprefactor_dKij[ind_j] * (u_i - u_j);
461  }
462  dFij_dKjk[sequential_i][j][index_j_to_i] += dprefactor_dKji * (u_i - u_j);
463  }
464  }
465  }
466 
467  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
468  {
469  const std::vector<dof_id_type> con_i =
471  const std::size_t num_con_i = con_i.size();
472  for (std::size_t j = 0; j < num_con_i; ++j)
473  {
474  const dof_id_type sequential_j = con_i[j];
475  if (sequential_i == sequential_j)
476  continue;
477  const unsigned index_j_to_i =
478  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
479  if (lij[sequential_j][index_j_to_i] < lij[sequential_i][j]) // node_i is downwind of node_j.
480  {
481  fa[sequential_i][j] = -fa[sequential_j][index_j_to_i];
482  for (const auto & dof_deriv : dfa[sequential_j][index_j_to_i])
483  dfa[sequential_i][j][dof_deriv.first] = -dof_deriv.second;
484  for (std::size_t k = 0; k < num_con_i; ++k)
485  dFij_dKik[sequential_i][j][k] = -dFij_dKjk[sequential_j][index_j_to_i][k];
486  const std::size_t num_con_j =
488  for (std::size_t k = 0; k < num_con_j; ++k)
489  dFij_dKjk[sequential_i][j][k] = -dFij_dKik[sequential_j][index_j_to_i][k];
490  }
491  }
492  }
493 
494  // zero _flux_out and its derivatives
495  _flux_out.assign(_number_of_nodes, 0.0);
496  // The derivatives are a bit complicated.
497  // If i is upwind of a node "j" then _flux_out[i] depends on all nodes connected to i.
498  // But if i is downwind of a node "j" then _flux_out depends on all nodes connected with node
499  // j.
500  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
501  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
502  for (const auto & j : _connections.globalConnectionsToSequentialID(sequential_i))
503  {
504  _dflux_out_du[sequential_i][j] = 0.0;
505  for (const auto & neighbors_j : _connections.globalConnectionsToGlobalID(j))
506  _dflux_out_du[sequential_i][neighbors_j] = 0.0;
507  }
509  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
510  {
511  const std::vector<dof_id_type> con_i =
513  const std::size_t num_con_i = con_i.size();
514  _dflux_out_dKjk[sequential_i].resize(num_con_i);
515  for (std::size_t j = 0; j < num_con_i; ++j)
516  {
517  const dof_id_type sequential_j = con_i[j];
518  std::vector<dof_id_type> con_j =
520  _dflux_out_dKjk[sequential_i][j].assign(con_j.size(), 0.0);
521  }
522  }
523 
524  // Add everything together
525  // See step 3 in Fig 2, noting Eqn (36)
526  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
527  {
528  std::vector<dof_id_type> con_i = _connections.sequentialConnectionsToSequentialID(sequential_i);
529  const size_t num_con_i = con_i.size();
530  const dof_id_type index_i_to_i =
531  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
532  for (std::size_t j = 0; j < num_con_i; ++j)
533  {
534  const dof_id_type sequential_j = con_i[j];
535  const dof_id_type global_j = _connections.globalID(sequential_j);
536  const Real u_j = _u_nodal[sequential_j];
537 
538  // negative sign because residual = -Lu (KT equation (19))
539  _flux_out[sequential_i] -= lij[sequential_i][j] * u_j + fa[sequential_i][j];
540 
541  _dflux_out_du[sequential_i][global_j] -= lij[sequential_i][j];
542  for (const auto & dof_deriv : dfa[sequential_i][j])
543  _dflux_out_du[sequential_i][dof_deriv.first] -= dof_deriv.second;
544 
545  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= 1.0 * u_j; // from the K in L = K + D
546 
547  if (sequential_j == sequential_i)
548  for (dof_id_type k = 0; k < num_con_i; ++k)
549  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dDii_dKij[sequential_i][k] * u_j;
550  else
551  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= dDij_dKij[sequential_i][j] * u_j;
552  for (dof_id_type k = 0; k < num_con_i; ++k)
553  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dFij_dKik[sequential_i][j][k];
554 
555  if (sequential_j == sequential_i)
556  for (unsigned k = 0; k < con_i.size(); ++k)
557  {
558  const unsigned index_k_to_i =
559  _connections.indexOfSequentialConnection(con_i[k], sequential_i);
560  _dflux_out_dKjk[sequential_i][k][index_k_to_i] -= dDii_dKji[sequential_i][k] * u_j;
561  }
562  else
563  {
564  const unsigned index_j_to_i =
565  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
566  _dflux_out_dKjk[sequential_i][j][index_j_to_i] -= dDij_dKji[sequential_i][j] * u_j;
567  }
568  for (unsigned k = 0;
569  k < _connections.sequentialConnectionsToSequentialID(sequential_j).size();
570  ++k)
571  _dflux_out_dKjk[sequential_i][j][k] -= dFij_dKjk[sequential_i][j][k];
572  }
573  }
574 }
575 
576 Real
577 AdvectiveFluxCalculatorBase::rPlus(dof_id_type sequential_i,
578  std::vector<Real> & dlimited_du,
579  std::vector<Real> & dlimited_dk) const
580 {
581  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
582  dlimited_du.assign(num_con, 0.0);
583  dlimited_dk.assign(num_con, 0.0);
585  return 0.0;
586  std::vector<Real> dp_du;
587  std::vector<Real> dp_dk;
588  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PPlus, dp_du, dp_dk);
589  if (p == 0.0)
590  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
591  return 1.0;
592  std::vector<Real> dq_du;
593  std::vector<Real> dq_dk;
594  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QPlus, dq_du, dq_dk);
595 
596  const Real r = q / p;
597  Real limited;
598  Real dlimited_dr;
599  limitFlux(1.0, r, limited, dlimited_dr);
600 
601  const Real p2 = std::pow(p, 2);
602  for (std::size_t j = 0; j < num_con; ++j)
603  {
604  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
605  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
606  dlimited_du[j] = dlimited_dr * dr_du;
607  dlimited_dk[j] = dlimited_dr * dr_dk;
608  }
609  return limited;
610 }
611 
612 Real
613 AdvectiveFluxCalculatorBase::rMinus(dof_id_type sequential_i,
614  std::vector<Real> & dlimited_du,
615  std::vector<Real> & dlimited_dk) const
616 {
617  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
618  dlimited_du.assign(num_con, 0.0);
619  dlimited_dk.assign(num_con, 0.0);
621  return 0.0;
622  std::vector<Real> dp_du;
623  std::vector<Real> dp_dk;
624  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PMinus, dp_du, dp_dk);
625  if (p == 0.0)
626  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
627  return 1.0;
628  std::vector<Real> dq_du;
629  std::vector<Real> dq_dk;
630  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QMinus, dq_du, dq_dk);
631 
632  const Real r = q / p;
633  Real limited;
634  Real dlimited_dr;
635  limitFlux(1.0, r, limited, dlimited_dr);
636 
637  const Real p2 = std::pow(p, 2);
638  for (std::size_t j = 0; j < num_con; ++j)
639  {
640  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
641  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
642  dlimited_du[j] = dlimited_dr * dr_du;
643  dlimited_dk[j] = dlimited_dr * dr_dk;
644  }
645  return limited;
646 }
647 
648 void
649 AdvectiveFluxCalculatorBase::limitFlux(Real a, Real b, Real & limited, Real & dlimited_db) const
650 {
651  limited = 0.0;
652  dlimited_db = 0.0;
654  return;
655 
656  if ((a >= 0.0 && b <= 0.0) || (a <= 0.0 && b >= 0.0))
657  return;
658  const Real s = (a > 0.0 ? 1.0 : -1.0);
659 
660  const Real lal = std::abs(a);
661  const Real lbl = std::abs(b);
662  const Real dlbl = (b >= 0.0 ? 1.0 : -1.0); // d(lbl)/db
663  switch (_flux_limiter_type)
664  {
666  {
667  if (lal <= lbl)
668  {
669  limited = s * lal;
670  dlimited_db = 0.0;
671  }
672  else
673  {
674  limited = s * lbl;
675  dlimited_db = s * dlbl;
676  }
677  return;
678  }
680  {
681  limited = s * 2 * lal * lbl / (lal + lbl);
682  dlimited_db = s * 2 * lal * (dlbl / (lal + lbl) - lbl * dlbl / std::pow(lal + lbl, 2));
683  return;
684  }
686  {
687  const Real av = 0.5 * std::abs(a + b);
688  if (2 * lal <= av && lal <= lbl)
689  {
690  // 2 * lal is the smallest
691  limited = s * 2.0 * lal;
692  dlimited_db = 0.0;
693  }
694  else if (2 * lbl <= av && lbl <= lal)
695  {
696  // 2 * lbl is the smallest
697  limited = s * 2.0 * lbl;
698  dlimited_db = s * 2.0 * dlbl;
699  }
700  else
701  {
702  // av is the smallest
703  limited = s * av;
704  // if (a>0 and b>0) then d(av)/db = 0.5 = 0.5 * dlbl
705  // if (a<0 and b<0) then d(av)/db = -0.5 = 0.5 * dlbl
706  // if a and b have different sign then limited=0, above
707  dlimited_db = s * 0.5 * dlbl;
708  }
709  return;
710  }
712  {
713  const Real term1 = std::min(2.0 * lal, lbl);
714  const Real term2 = std::min(lal, 2.0 * lbl);
715  if (term1 >= term2)
716  {
717  if (2.0 * lal <= lbl)
718  {
719  limited = s * 2 * lal;
720  dlimited_db = 0.0;
721  }
722  else
723  {
724  limited = s * lbl;
725  dlimited_db = s * dlbl;
726  }
727  }
728  else
729  {
730  if (lal <= 2.0 * lbl)
731  {
732  limited = s * lal;
733  dlimited_db = 0.0;
734  }
735  else
736  {
737  limited = s * 2.0 * lbl;
738  dlimited_db = s * 2.0 * dlbl;
739  }
740  }
741  return;
742  }
743  default:
744  return;
745  }
746 }
747 
748 const std::map<dof_id_type, Real> &
750 {
751  return _dflux_out_du[_connections.sequentialID(node_i)];
752 }
753 
754 const std::vector<std::vector<Real>> &
756 {
757  return _dflux_out_dKjk[_connections.sequentialID(node_i)];
758 }
759 
760 Real
762 {
763  return _flux_out[_connections.sequentialID(node_i)];
764 }
765 
766 unsigned
768 {
769  return _valence[_connections.sequentialID(node_i)];
770 }
771 
772 void
773 AdvectiveFluxCalculatorBase::zeroedConnection(std::map<dof_id_type, Real> & the_map,
774  dof_id_type node_i) const
775 {
776  the_map.clear();
777  for (const auto & node_j : _connections.globalConnectionsToGlobalID(node_i))
778  the_map[node_j] = 0.0;
779 }
780 
781 Real
783  const PQPlusMinusEnum pq_plus_minus,
784  std::vector<Real> & derivs,
785  std::vector<Real> & dpqdk) const
786 {
787  // Find the value of u at sequential_i
788  const Real u_i = _u_nodal[sequential_i];
789 
790  // Connections to sequential_i
791  const std::vector<dof_id_type> con_i =
793  const std::size_t num_con = con_i.size();
794  // The neighbor number of sequential_i to sequential_i
795  const unsigned i_index_i = _connections.indexOfSequentialConnection(sequential_i, sequential_i);
796 
797  // Initialize the results
798  Real result = 0.0;
799  derivs.assign(num_con, 0.0);
800  dpqdk.assign(num_con, 0.0);
801 
802  // Sum over all nodes connected with node_i.
803  for (std::size_t j = 0; j < num_con; ++j)
804  {
805  const dof_id_type sequential_j = con_i[j];
806  if (sequential_j == sequential_i)
807  continue;
808  const Real kentry = _kij[sequential_i][j];
809 
810  // Find the value of u at node_j
811  const Real u_j = _u_nodal[sequential_j];
812  const Real ujminusi = u_j - u_i;
813 
814  // Evaluate the i-j contribution to the result
815  switch (pq_plus_minus)
816  {
818  {
819  if (ujminusi < 0.0 && kentry < 0.0)
820  {
821  result += kentry * ujminusi;
822  derivs[j] += kentry;
823  derivs[i_index_i] -= kentry;
824  dpqdk[j] += ujminusi;
825  }
826  break;
827  }
829  {
830  if (ujminusi > 0.0 && kentry < 0.0)
831  {
832  result += kentry * ujminusi;
833  derivs[j] += kentry;
834  derivs[i_index_i] -= kentry;
835  dpqdk[j] += ujminusi;
836  }
837  break;
838  }
840  {
841  if (ujminusi > 0.0 && kentry > 0.0)
842  {
843  result += kentry * ujminusi;
844  derivs[j] += kentry;
845  derivs[i_index_i] -= kentry;
846  dpqdk[j] += ujminusi;
847  }
848  break;
849  }
851  {
852  if (ujminusi < 0.0 && kentry > 0.0)
853  {
854  result += kentry * ujminusi;
855  derivs[j] += kentry;
856  derivs[i_index_i] -= kentry;
857  dpqdk[j] += ujminusi;
858  }
859  break;
860  }
861  }
862  }
863 
864  return result;
865 }
866 
867 void
869 {
883  _nodes_to_receive.clear();
884  for (const auto & elem : _fe_problem.getEvaluableElementRange())
885  if (this->hasBlocks(elem->subdomain_id()))
886  {
887  const processor_id_type elem_pid = elem->processor_id();
888  if (elem_pid != _my_pid)
889  {
890  if (_nodes_to_receive.find(elem_pid) == _nodes_to_receive.end())
891  _nodes_to_receive[elem_pid] = std::vector<dof_id_type>();
892  for (unsigned i = 0; i < elem->n_nodes(); ++i)
893  if (std::find(_nodes_to_receive[elem_pid].begin(),
894  _nodes_to_receive[elem_pid].end(),
895  elem->node_id(i)) == _nodes_to_receive[elem_pid].end())
896  _nodes_to_receive[elem_pid].push_back(elem->node_id(i));
897  }
898  }
899 
900  // exchange this info with other processors, building _nodes_to_send at the same time
901  _nodes_to_send.clear();
902  auto nodes_action_functor = [this](processor_id_type pid, const std::vector<dof_id_type> & nts) {
903  _nodes_to_send[pid] = nts;
904  };
905  Parallel::push_parallel_vector_data(this->comm(), _nodes_to_receive, nodes_action_functor);
906 
907  // At the moment, _nodes_to_send and _nodes_to_receive contain global node numbers
908  // It is slightly more efficient to convert to sequential node numbers
909  // so that we don't have to keep doing things like
910  // _u_nodal[_connections.sequentialID(_nodes_to_send[pid][i])
911  // every time we send/receive
912  for (auto & kv : _nodes_to_send)
913  {
914  const processor_id_type pid = kv.first;
915  const std::size_t num_nodes = kv.second.size();
916  for (unsigned i = 0; i < num_nodes; ++i)
918  }
919  for (auto & kv : _nodes_to_receive)
920  {
921  const processor_id_type pid = kv.first;
922  const std::size_t num_nodes = kv.second.size();
923  for (unsigned i = 0; i < num_nodes; ++i)
925  }
926 
927  // Build pairs_to_receive
928  _pairs_to_receive.clear();
929  for (const auto & elem : _fe_problem.getEvaluableElementRange())
930  if (this->hasBlocks(elem->subdomain_id()))
931  {
932  const processor_id_type elem_pid = elem->processor_id();
933  if (elem_pid != _my_pid)
934  {
935  if (_pairs_to_receive.find(elem_pid) == _pairs_to_receive.end())
936  _pairs_to_receive[elem_pid] = std::vector<std::pair<dof_id_type, dof_id_type>>();
937  for (unsigned i = 0; i < elem->n_nodes(); ++i)
938  for (unsigned j = 0; j < elem->n_nodes(); ++j)
939  {
940  std::pair<dof_id_type, dof_id_type> the_pair(elem->node_id(i), elem->node_id(j));
941  if (std::find(_pairs_to_receive[elem_pid].begin(),
942  _pairs_to_receive[elem_pid].end(),
943  the_pair) == _pairs_to_receive[elem_pid].end())
944  _pairs_to_receive[elem_pid].push_back(the_pair);
945  }
946  }
947  }
948 
949  _pairs_to_send.clear();
950  auto pairs_action_functor = [this](processor_id_type pid,
951  const std::vector<std::pair<dof_id_type, dof_id_type>> & pts) {
952  _pairs_to_send[pid] = pts;
953  };
954  Parallel::push_parallel_vector_data(this->comm(), _pairs_to_receive, pairs_action_functor);
955 
956  // _pairs_to_send and _pairs_to_receive have been built using global node IDs
957  // since all processors know about that. However, using global IDs means
958  // that every time we send/receive, we keep having to do things like
959  // _kij[_connections.sequentialID(_pairs_to_send[pid][i].first)][_connections.indexOfGlobalConnection(_pairs_to_send[pid][i].first,
960  // _pairs_to_send[pid][i].second)] which is quite inefficient. So:
961  for (auto & kv : _pairs_to_send)
962  {
963  const processor_id_type pid = kv.first;
964  const std::size_t num_pairs = kv.second.size();
965  for (unsigned i = 0; i < num_pairs; ++i)
966  {
968  _pairs_to_send[pid][i].first, _pairs_to_send[pid][i].second);
969  _pairs_to_send[pid][i].first = _connections.sequentialID(_pairs_to_send[pid][i].first);
970  }
971  }
972  for (auto & kv : _pairs_to_receive)
973  {
974  const processor_id_type pid = kv.first;
975  const std::size_t num_pairs = kv.second.size();
976  for (unsigned i = 0; i < num_pairs; ++i)
977  {
979  _pairs_to_receive[pid][i].first, _pairs_to_receive[pid][i].second);
980  _pairs_to_receive[pid][i].first = _connections.sequentialID(_pairs_to_receive[pid][i].first);
981  }
982  }
983 }
984 
985 void
987 {
988  // Exchange ghosted _u_nodal information with other processors
989  std::map<processor_id_type, std::vector<Real>> unodal_to_send;
990  for (const auto & kv : _nodes_to_send)
991  {
992  const processor_id_type pid = kv.first;
993  unodal_to_send[pid] = std::vector<Real>();
994  for (const auto & nd : kv.second)
995  unodal_to_send[pid].push_back(_u_nodal[nd]);
996  }
997 
998  auto unodal_action_functor = [this](processor_id_type pid,
999  const std::vector<Real> & unodal_received) {
1000  const std::size_t msg_size = unodal_received.size();
1001  mooseAssert(msg_size == _nodes_to_receive[pid].size(),
1002  "Message size, " << msg_size
1003  << ", incompatible with nodes_to_receive, which has size "
1004  << _nodes_to_receive[pid].size());
1005  for (unsigned i = 0; i < msg_size; ++i)
1006  _u_nodal[_nodes_to_receive[pid][i]] = unodal_received[i];
1007  };
1008  Parallel::push_parallel_vector_data(this->comm(), unodal_to_send, unodal_action_functor);
1009 
1010  // Exchange _kij information with other processors
1011  std::map<processor_id_type, std::vector<Real>> kij_to_send;
1012  for (const auto & kv : _pairs_to_send)
1013  {
1014  const processor_id_type pid = kv.first;
1015  kij_to_send[pid] = std::vector<Real>();
1016  for (const auto & pr : kv.second)
1017  kij_to_send[pid].push_back(_kij[pr.first][pr.second]);
1018  }
1019 
1020  auto kij_action_functor = [this](processor_id_type pid, const std::vector<Real> & kij_received) {
1021  const std::size_t msg_size = kij_received.size();
1022  mooseAssert(msg_size == _pairs_to_receive[pid].size(),
1023  "Message size, " << msg_size
1024  << ", incompatible with pairs_to_receive, which has size "
1025  << _pairs_to_receive[pid].size());
1026  for (unsigned i = 0; i < msg_size; ++i)
1027  _kij[_pairs_to_receive[pid][i].first][_pairs_to_receive[pid][i].second] += kij_received[i];
1028  };
1029  Parallel::push_parallel_vector_data(this->comm(), kij_to_send, kij_action_functor);
1030 }
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 ...
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 ...
FluxLimiterTypeEnum
Determines Flux Limiter type (Page 135 of Kuzmin and Turek) "None" means that limitFlux=0 always...
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...
Real getFluxOut(dof_id_type node_i) const
Returns the flux out of lobal node id.
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
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 ...
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.
void finalizeAddingConnections()
Signal that all global node IDs have been added to the internal data structures.
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...
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::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]).
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 ...
InputParameters validParams< AdvectiveFluxCalculatorBase >()
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. ...
std::size_t sizeSequential() const
Return the size of _sequential_id, for checking memory efficiency.
const std::string name
Definition: Setup.h:21
Base class to compute Advective fluxes.
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
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::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_send
_nodes_to_send[proc_id] = list of sequential nodal IDs.
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...
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...
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::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_send
_pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is fi...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
processor_id_type _my_pid
processor ID of this object
virtual void exchangeGhostedInfo()
Sends and receives multi-processor information regarding u_nodal and k_ij.
const Real _allowable_MB_wastage
A mooseWarning is issued if mb_wasted = (_connections.sizeSequential() - _connections.numNodes()) * 4 / 1048576 > _allowable_MB_wastage.
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...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual Real computeU(unsigned i) const =0
Computes the value of u at the local node id of the current element (_current_elem) ...
std::size_t numNodes() const
number of nodes known by this class
AdvectiveFluxCalculatorBase(const InputParameters &parameters)
virtual void buildCommLists()
When using multiple processors, other processors will compute:
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_receive
_nodes_to_receive[proc_id] = list of sequential nodal IDs.
PQPlusMinusEnum
Signals to the PQPlusMinus method what should be computed.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_receive
_pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to...
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.
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
virtual void threadJoin(const UserObject &uo) override
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...
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_...