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