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 3701 : AdvectiveFluxCalculatorBase::validParams()
19 : {
20 3701 : InputParameters params = ElementUserObject::validParams();
21 3701 : 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 7402 : MooseEnum flux_limiter_type("MinMod VanLeer MC superbee None", "VanLeer");
26 7402 : 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 11103 : params.addRangeCheckedParam<Real>(
31 : "allowable_MB_wastage",
32 7402 : 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 7402 : params.addRelationshipManager("ElementPointNeighborLayers",
39 : Moose::RelationshipManagerType::GEOMETRIC |
40 : Moose::RelationshipManagerType::ALGEBRAIC |
41 : Moose::RelationshipManagerType::COUPLING,
42 3675 : [](const InputParameters &, InputParameters & rm_params)
43 3675 : { rm_params.set<unsigned short>("layers") = 2; });
44 :
45 11103 : params.set<ExecFlagEnum>("execute_on", true) = {EXEC_LINEAR};
46 3701 : return params;
47 7402 : }
48 :
49 1664 : AdvectiveFluxCalculatorBase::AdvectiveFluxCalculatorBase(const InputParameters & parameters)
50 : : ElementUserObject(parameters),
51 1664 : _resizing_needed(true),
52 1664 : _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 1664 : _connections(),
61 1664 : _number_of_nodes(0),
62 1664 : _my_pid(processor_id()),
63 : _nodes_to_receive(),
64 : _nodes_to_send(),
65 : _pairs_to_receive(),
66 : _pairs_to_send(),
67 4992 : _allowable_MB_wastage(getParam<Real>("allowable_MB_wastage"))
68 : {
69 1664 : 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 1662 : }
77 :
78 : void
79 58342 : AdvectiveFluxCalculatorBase::timestepSetup()
80 : {
81 : // If needed, size and initialize quantities appropriately, and compute _valence
82 58342 : 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 2802 : _connections.clear();
95 52679 : for (const auto & elem : _fe_problem.getNonlinearEvaluableElementRange())
96 49877 : if (this->hasBlocks(elem->subdomain_id()))
97 172745 : for (unsigned i = 0; i < elem->n_nodes(); ++i)
98 126564 : _connections.addGlobalNode(elem->node_id(i));
99 2802 : _connections.finalizeAddingGlobalNodes();
100 52679 : for (const auto & elem : _fe_problem.getNonlinearEvaluableElementRange())
101 49877 : if (this->hasBlocks(elem->subdomain_id()))
102 172745 : for (unsigned i = 0; i < elem->n_nodes(); ++i)
103 532604 : for (unsigned j = 0; j < elem->n_nodes(); ++j)
104 406040 : _connections.addConnection(elem->node_id(i), elem->node_id(j));
105 2802 : _connections.finalizeAddingConnections();
106 :
107 2802 : _number_of_nodes = _connections.numNodes();
108 :
109 2802 : Real mb_wasted = (_connections.sizeSequential() - _number_of_nodes) * 4.0 / 1048576;
110 : gatherMax(mb_wasted);
111 2802 : 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 2802 : _kij.resize(_number_of_nodes);
123 59095 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
124 : _kij[sequential_i].assign(
125 112586 : _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 2802 : _valence.assign(_number_of_nodes, 0);
139 52679 : for (const auto & elem : _fe_problem.getNonlinearEvaluableElementRange())
140 49877 : if (this->hasBlocks(elem->subdomain_id()))
141 172745 : for (unsigned i = 0; i < elem->n_nodes(); ++i)
142 : {
143 126564 : const dof_id_type node_i = elem->node_id(i);
144 126564 : const dof_id_type sequential_i = _connections.sequentialID(node_i);
145 126564 : _valence[sequential_i] += 1;
146 : }
147 :
148 2802 : _u_nodal.assign(_number_of_nodes, 0.0);
149 2802 : _u_nodal_computed_by_thread.assign(_number_of_nodes, false);
150 2802 : _flux_out.assign(_number_of_nodes, 0.0);
151 2802 : _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
152 59095 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
153 56293 : zeroedConnection(_dflux_out_du[sequential_i], _connections.globalID(sequential_i));
154 2802 : _dflux_out_dKjk.resize(_number_of_nodes);
155 59095 : 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 56293 : _connections.sequentialConnectionsToSequentialID(sequential_i);
159 : const std::size_t num_con_i = con_i.size();
160 56293 : _dflux_out_dKjk[sequential_i].resize(num_con_i);
161 333018 : for (std::size_t j = 0; j < con_i.size(); ++j)
162 : {
163 276725 : const dof_id_type sequential_j = con_i[j];
164 : const std::size_t num_con_j =
165 276725 : _connections.sequentialConnectionsToSequentialID(sequential_j).size();
166 276725 : _dflux_out_dKjk[sequential_i][j].assign(num_con_j, 0.0);
167 : }
168 56293 : }
169 :
170 2802 : if (_app.n_processors() > 1)
171 2028 : buildCommLists();
172 :
173 2802 : _resizing_needed = false;
174 :
175 : // Clear and size all member vectors
176 2802 : _dij.assign(_number_of_nodes, {});
177 2802 : _dDij_dKij.assign(_number_of_nodes, {});
178 2802 : _dDij_dKji.assign(_number_of_nodes, {});
179 2802 : _dDii_dKij.assign(_number_of_nodes, {});
180 2802 : _dDii_dKji.assign(_number_of_nodes, {});
181 2802 : _lij.assign(_number_of_nodes, {});
182 2802 : _rP.assign(_number_of_nodes, 0.0);
183 2802 : _rM.assign(_number_of_nodes, 0.0);
184 2802 : _drP.assign(_number_of_nodes, {});
185 2802 : _drM.assign(_number_of_nodes, {});
186 2802 : _drP_dk.assign(_number_of_nodes, {});
187 2802 : _drM_dk.assign(_number_of_nodes, {});
188 2802 : _fa.assign(_number_of_nodes, {});
189 2802 : _dfa.assign(_number_of_nodes, {});
190 2802 : _dFij_dKik.assign(_number_of_nodes, {});
191 2802 : _dFij_dKjk.assign(_number_of_nodes, {});
192 : }
193 58342 : }
194 :
195 : void
196 1329 : AdvectiveFluxCalculatorBase::meshChanged()
197 : {
198 : ElementUserObject::meshChanged();
199 :
200 : // Signal that _kij, _valence, etc need to be rebuilt
201 1329 : _resizing_needed = true;
202 1329 : }
203 :
204 : void
205 195715 : AdvectiveFluxCalculatorBase::initialize()
206 : {
207 : // Zero _kij and falsify _u_nodal_computed_by_thread ready for building in execute() and
208 : // finalize()
209 195715 : _u_nodal_computed_by_thread.assign(_number_of_nodes, false);
210 5251309 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
211 5055594 : _kij[sequential_i].assign(_connections.sequentialConnectionsToSequentialID(sequential_i).size(),
212 10111188 : 0.0);
213 195715 : }
214 :
215 : void
216 2250490 : AdvectiveFluxCalculatorBase::execute()
217 : {
218 : // compute _kij contributions from this element that is local to this processor
219 : // and record _u_nodal
220 7820830 : for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
221 : {
222 5570340 : const dof_id_type node_i = _current_elem->node_id(i);
223 5570340 : const dof_id_type sequential_i = _connections.sequentialID(node_i);
224 5570340 : if (!_u_nodal_computed_by_thread[sequential_i])
225 : {
226 2802092 : _u_nodal[sequential_i] = computeU(i);
227 : _u_nodal_computed_by_thread[sequential_i] = true;
228 : }
229 21199660 : for (unsigned j = 0; j < _current_elem->n_nodes(); ++j)
230 : {
231 15629320 : const dof_id_type node_j = _current_elem->node_id(j);
232 66532120 : for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
233 50902800 : executeOnElement(node_i, node_j, i, j, qp);
234 : }
235 : }
236 2250490 : }
237 :
238 : void
239 50902800 : 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 50902800 : const dof_id_type sequential_i = _connections.sequentialID(global_i);
244 50902800 : const unsigned index_i_to_j = _connections.indexOfGlobalConnection(global_i, global_j);
245 50902800 : _kij[sequential_i][index_i_to_j] += _JxW[qp] * _coord[qp] * computeVelocity(local_i, local_j, qp);
246 50902800 : }
247 :
248 : void
249 88051 : 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 2239195 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
254 : {
255 : const std::size_t num_con_i =
256 2151144 : _connections.sequentialConnectionsToSequentialID(sequential_i).size();
257 12076502 : for (std::size_t j = 0; j < num_con_i; ++j)
258 9925358 : _kij[sequential_i][j] += afc._kij[sequential_i][j];
259 : }
260 :
261 : // gather the values of _u_nodal computed by different threads
262 2239195 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
263 2151144 : if (!_u_nodal_computed_by_thread[sequential_i] && afc._u_nodal_computed_by_thread[sequential_i])
264 724315 : _u_nodal[sequential_i] = afc._u_nodal[sequential_i];
265 88051 : }
266 :
267 : void
268 107664 : 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 107664 : if (_app.n_processors() > 1)
275 77980 : 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 3012114 : 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 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
286 : const std::size_t num_con_i = con_i.size();
287 2904450 : _dij[sequential_i].assign(num_con_i, 0.0);
288 2904450 : _dDij_dKij[sequential_i].assign(num_con_i, 0.0);
289 2904450 : _dDij_dKji[sequential_i].assign(num_con_i, 0.0);
290 2904450 : _dDii_dKij[sequential_i].assign(num_con_i, 0.0);
291 2904450 : _dDii_dKji[sequential_i].assign(num_con_i, 0.0);
292 : const unsigned index_i_to_i =
293 2904450 : _connections.indexOfSequentialConnection(sequential_i, sequential_i);
294 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
295 : {
296 12726050 : const dof_id_type sequential_j = con_i[j];
297 12726050 : if (sequential_i == sequential_j)
298 2904450 : continue;
299 : const unsigned index_j_to_i =
300 9821600 : _connections.indexOfSequentialConnection(sequential_j, sequential_i);
301 9821600 : const Real kij = _kij[sequential_i][j];
302 9821600 : const Real kji = _kij[sequential_j][index_j_to_i];
303 9821600 : if ((kij <= kji) && (kij < 0.0))
304 : {
305 4770415 : _dij[sequential_i][j] = -kij;
306 4770415 : _dDij_dKij[sequential_i][j] = -1.0;
307 4770415 : _dDii_dKij[sequential_i][j] += 1.0;
308 : }
309 5051185 : else if ((kji <= kij) && (kji < 0.0))
310 : {
311 4411621 : _dij[sequential_i][j] = -kji;
312 4411621 : _dDij_dKji[sequential_i][j] = -1.0;
313 4411621 : _dDii_dKji[sequential_i][j] += 1.0;
314 : }
315 : else
316 639564 : _dij[sequential_i][j] = 0.0;
317 9821600 : _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 3012114 : 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 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
327 : const std::size_t num_con_i = con_i.size();
328 2904450 : _lij[sequential_i].assign(num_con_i, 0.0);
329 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
330 12726050 : _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 3012114 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
336 : {
337 2904450 : _rP[sequential_i] = rPlus(sequential_i, _drP[sequential_i], _drP_dk[sequential_i]);
338 2904450 : _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 3012114 : 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 2904450 : _connections.globalConnectionsToSequentialID(sequential_i);
348 : const unsigned num_con_i = con_i.size();
349 2904450 : _fa[sequential_i].assign(num_con_i, 0.0);
350 2904450 : _dfa[sequential_i].resize(num_con_i);
351 2904450 : _dFij_dKik[sequential_i].resize(num_con_i);
352 2904450 : _dFij_dKjk[sequential_i].resize(num_con_i);
353 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
354 : {
355 85720400 : for (const auto & global_k : con_i)
356 72994350 : _dfa[sequential_i][j][global_k] = 0;
357 12726050 : const dof_id_type global_j = con_i[j];
358 12726050 : const std::vector<dof_id_type> & con_j = _connections.globalConnectionsToGlobalID(global_j);
359 : const unsigned num_con_j = con_j.size();
360 85720400 : for (const auto & global_k : con_j)
361 72994350 : _dfa[sequential_i][j][global_k] = 0;
362 12726050 : _dFij_dKik[sequential_i][j].assign(num_con_i, 0.0);
363 12726050 : _dFij_dKjk[sequential_i][j].assign(num_con_j, 0.0);
364 : }
365 : }
366 :
367 3012114 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
368 : {
369 2904450 : const dof_id_type global_i = _connections.globalID(sequential_i);
370 2904450 : const Real u_i = _u_nodal[sequential_i];
371 : const std::vector<dof_id_type> & con_i =
372 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
373 : const std::size_t num_con_i = con_i.size();
374 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
375 : {
376 12726050 : const dof_id_type sequential_j = con_i[j];
377 12726050 : if (sequential_i == sequential_j)
378 2904450 : continue;
379 : const unsigned index_j_to_i =
380 9821600 : _connections.indexOfSequentialConnection(sequential_j, sequential_i);
381 9821600 : const Real Lij = _lij[sequential_i][j];
382 9821600 : const Real Lji = _lij[sequential_j][index_j_to_i];
383 9821600 : if (Lji >= Lij) // node i is upwind of node j.
384 : {
385 5337009 : const Real Dij = _dij[sequential_i][j];
386 5337009 : const dof_id_type global_j = _connections.globalID(sequential_j);
387 5337009 : const Real u_j = _u_nodal[sequential_j];
388 : Real prefactor = 0.0;
389 : std::vector<Real> dprefactor_du(num_con_i,
390 5337009 : 0.0); // dprefactor_du[j] = d(prefactor)/du[sequential_j]
391 : std::vector<Real> dprefactor_dKij(
392 5337009 : 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 5337009 : if (u_i >= u_j)
395 : {
396 3197621 : if (Lji <= _rP[sequential_i] * Dij)
397 : {
398 : prefactor = Lji;
399 717146 : dprefactor_dKij[j] += _dDij_dKji[sequential_j][index_j_to_i];
400 717146 : dprefactor_dKji += 1.0 + _dDij_dKij[sequential_j][index_j_to_i];
401 : }
402 : else
403 : {
404 : prefactor = _rP[sequential_i] * Dij;
405 16119351 : for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
406 13638876 : dprefactor_du[ind_j] = _drP[sequential_i][ind_j] * Dij;
407 2480475 : dprefactor_dKij[j] += _rP[sequential_i] * _dDij_dKij[sequential_i][j];
408 2480475 : dprefactor_dKji += _rP[sequential_i] * _dDij_dKji[sequential_i][j];
409 16119351 : for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
410 13638876 : dprefactor_dKij[ind_j] += _drP_dk[sequential_i][ind_j] * Dij;
411 : }
412 : }
413 : else
414 : {
415 2139388 : if (Lji <= _rM[sequential_i] * Dij)
416 : {
417 : prefactor = Lji;
418 483920 : dprefactor_dKij[j] += _dDij_dKji[sequential_j][index_j_to_i];
419 483920 : dprefactor_dKji += 1.0 + _dDij_dKij[sequential_j][index_j_to_i];
420 : }
421 : else
422 : {
423 : prefactor = _rM[sequential_i] * Dij;
424 12634227 : for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
425 10978759 : dprefactor_du[ind_j] = _drM[sequential_i][ind_j] * Dij;
426 1655468 : dprefactor_dKij[j] += _rM[sequential_i] * _dDij_dKij[sequential_i][j];
427 1655468 : dprefactor_dKji += _rM[sequential_i] * _dDij_dKji[sequential_i][j];
428 12634227 : for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
429 10978759 : dprefactor_dKij[ind_j] += _drM_dk[sequential_i][ind_j] * Dij;
430 : }
431 : }
432 5337009 : _fa[sequential_i][j] = prefactor * (u_i - u_j);
433 5337009 : _dfa[sequential_i][j][global_i] = prefactor;
434 5337009 : _dfa[sequential_i][j][global_j] = -prefactor;
435 38516515 : for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
436 : {
437 33179506 : _dfa[sequential_i][j][_connections.globalID(con_i[ind_j])] +=
438 33179506 : dprefactor_du[ind_j] * (u_i - u_j);
439 33179506 : _dFij_dKik[sequential_i][j][ind_j] += dprefactor_dKij[ind_j] * (u_i - u_j);
440 : }
441 5337009 : _dFij_dKjk[sequential_i][j][index_j_to_i] += dprefactor_dKji * (u_i - u_j);
442 5337009 : }
443 : }
444 : }
445 :
446 3012114 : 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 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
450 : const std::size_t num_con_i = con_i.size();
451 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
452 : {
453 12726050 : const dof_id_type sequential_j = con_i[j];
454 12726050 : if (sequential_i == sequential_j)
455 2904450 : continue;
456 : const unsigned index_j_to_i =
457 9821600 : _connections.indexOfSequentialConnection(sequential_j, sequential_i);
458 9821600 : if (_lij[sequential_j][index_j_to_i] < _lij[sequential_i][j]) // node_i is downwind of node_j.
459 : {
460 4484591 : _fa[sequential_i][j] = -_fa[sequential_j][index_j_to_i];
461 41564125 : for (const auto & dof_deriv : _dfa[sequential_j][index_j_to_i])
462 37079534 : _dfa[sequential_i][j][dof_deriv.first] = -dof_deriv.second;
463 31573385 : for (std::size_t k = 0; k < num_con_i; ++k)
464 27088794 : _dFij_dKik[sequential_i][j][k] = -_dFij_dKjk[sequential_j][index_j_to_i][k];
465 : const std::size_t num_con_j =
466 4484591 : _connections.sequentialConnectionsToSequentialID(sequential_j).size();
467 31565137 : for (std::size_t k = 0; k < num_con_j; ++k)
468 27080546 : _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 107664 : _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 107664 : _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
480 3012114 : for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
481 15630500 : for (const auto & j : _connections.globalConnectionsToSequentialID(sequential_i))
482 : {
483 12726050 : _dflux_out_du[sequential_i][j] = 0.0;
484 85720400 : for (const auto & neighbors_j : _connections.globalConnectionsToGlobalID(j))
485 72994350 : _dflux_out_du[sequential_i][neighbors_j] = 0.0;
486 : }
487 107664 : _dflux_out_dKjk.resize(_number_of_nodes);
488 3012114 : 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 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
492 : const std::size_t num_con_i = con_i.size();
493 2904450 : _dflux_out_dKjk[sequential_i].resize(num_con_i);
494 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
495 : {
496 12726050 : const dof_id_type sequential_j = con_i[j];
497 : const std::vector<dof_id_type> & con_j =
498 12726050 : _connections.sequentialConnectionsToSequentialID(sequential_j);
499 12726050 : _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 3012114 : 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 2904450 : _connections.sequentialConnectionsToSequentialID(sequential_i);
509 : const size_t num_con_i = con_i.size();
510 : const dof_id_type index_i_to_i =
511 2904450 : _connections.indexOfSequentialConnection(sequential_i, sequential_i);
512 15630500 : for (std::size_t j = 0; j < num_con_i; ++j)
513 : {
514 12726050 : const dof_id_type sequential_j = con_i[j];
515 12726050 : const dof_id_type global_j = _connections.globalID(sequential_j);
516 12726050 : const Real u_j = _u_nodal[sequential_j];
517 :
518 : // negative sign because residual = -Lu (KT equation (19))
519 12726050 : _flux_out[sequential_i] -= _lij[sequential_i][j] * u_j + _fa[sequential_i][j];
520 :
521 12726050 : _dflux_out_du[sequential_i][global_j] -= _lij[sequential_i][j];
522 107116316 : for (const auto & dof_deriv : _dfa[sequential_i][j])
523 94390266 : _dflux_out_du[sequential_i][dof_deriv.first] -= dof_deriv.second;
524 :
525 12726050 : _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= 1.0 * u_j; // from the K in L = K + D
526 :
527 12726050 : if (sequential_j == sequential_i)
528 15630500 : for (dof_id_type k = 0; k < num_con_i; ++k)
529 12726050 : _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= _dDii_dKij[sequential_i][k] * u_j;
530 : else
531 9821600 : _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= _dDij_dKij[sequential_i][j] * u_j;
532 85720400 : for (dof_id_type k = 0; k < num_con_i; ++k)
533 72994350 : _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= _dFij_dKik[sequential_i][j][k];
534 :
535 12726050 : if (sequential_j == sequential_i)
536 15630500 : for (unsigned k = 0; k < con_i.size(); ++k)
537 : {
538 : const unsigned index_k_to_i =
539 12726050 : _connections.indexOfSequentialConnection(con_i[k], sequential_i);
540 12726050 : _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 9821600 : _connections.indexOfSequentialConnection(sequential_j, sequential_i);
546 9821600 : _dflux_out_dKjk[sequential_i][j][index_j_to_i] -= _dDij_dKji[sequential_i][j] * u_j;
547 : }
548 72994350 : for (unsigned k = 0;
549 85720400 : k < _connections.sequentialConnectionsToSequentialID(sequential_j).size();
550 : ++k)
551 72994350 : _dflux_out_dKjk[sequential_i][j][k] -= _dFij_dKjk[sequential_i][j][k];
552 : }
553 : }
554 107664 : }
555 :
556 : Real
557 2904450 : AdvectiveFluxCalculatorBase::rPlus(dof_id_type sequential_i,
558 : std::vector<Real> & dlimited_du,
559 : std::vector<Real> & dlimited_dk) const
560 : {
561 2904450 : const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
562 2904450 : dlimited_du.assign(num_con, 0.0);
563 2904450 : dlimited_dk.assign(num_con, 0.0);
564 2904450 : if (_flux_limiter_type == FluxLimiterTypeEnum::None)
565 : return 0.0;
566 : std::vector<Real> dp_du;
567 : std::vector<Real> dp_dk;
568 2865600 : const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PPlus, dp_du, dp_dk);
569 2865600 : 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 1830940 : const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QPlus, dq_du, dq_dk);
575 :
576 1830940 : const Real r = q / p;
577 : Real limited;
578 : Real dlimited_dr;
579 1830940 : limitFlux(1.0, r, limited, dlimited_dr);
580 :
581 : const Real p2 = std::pow(p, 2);
582 9787434 : for (std::size_t j = 0; j < num_con; ++j)
583 : {
584 7956494 : const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
585 7956494 : const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
586 7956494 : dlimited_du[j] = dlimited_dr * dr_du;
587 7956494 : dlimited_dk[j] = dlimited_dr * dr_dk;
588 : }
589 1830940 : return limited;
590 2865600 : }
591 :
592 : Real
593 2904450 : AdvectiveFluxCalculatorBase::rMinus(dof_id_type sequential_i,
594 : std::vector<Real> & dlimited_du,
595 : std::vector<Real> & dlimited_dk) const
596 : {
597 2904450 : const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
598 2904450 : dlimited_du.assign(num_con, 0.0);
599 2904450 : dlimited_dk.assign(num_con, 0.0);
600 2904450 : if (_flux_limiter_type == FluxLimiterTypeEnum::None)
601 : return 0.0;
602 : std::vector<Real> dp_du;
603 : std::vector<Real> dp_dk;
604 2865600 : const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PMinus, dp_du, dp_dk);
605 2865600 : 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 1036491 : const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QMinus, dq_du, dq_dk);
611 :
612 1036491 : const Real r = q / p;
613 : Real limited;
614 : Real dlimited_dr;
615 1036491 : limitFlux(1.0, r, limited, dlimited_dr);
616 :
617 : const Real p2 = std::pow(p, 2);
618 6759276 : for (std::size_t j = 0; j < num_con; ++j)
619 : {
620 5722785 : const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
621 5722785 : const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
622 5722785 : dlimited_du[j] = dlimited_dr * dr_du;
623 5722785 : dlimited_dk[j] = dlimited_dr * dr_dk;
624 : }
625 1036491 : return limited;
626 2865600 : }
627 :
628 : void
629 2867431 : AdvectiveFluxCalculatorBase::limitFlux(Real a, Real b, Real & limited, Real & dlimited_db) const
630 : {
631 2867431 : limited = 0.0;
632 2867431 : dlimited_db = 0.0;
633 2867431 : if (_flux_limiter_type == FluxLimiterTypeEnum::None)
634 : return;
635 :
636 2867431 : if ((a >= 0.0 && b <= 0.0) || (a <= 0.0 && b >= 0.0))
637 : return;
638 2325389 : const Real s = (a > 0.0 ? 1.0 : -1.0);
639 :
640 2325389 : const Real lal = std::abs(a);
641 2325389 : const Real lbl = std::abs(b);
642 2325389 : const Real dlbl = (b >= 0.0 ? 1.0 : -1.0); // d(lbl)/db
643 2325389 : switch (_flux_limiter_type)
644 : {
645 15947 : case FluxLimiterTypeEnum::MinMod:
646 : {
647 15947 : if (lal <= lbl)
648 : {
649 8330 : limited = s * lal;
650 8330 : dlimited_db = 0.0;
651 : }
652 : else
653 : {
654 7617 : limited = s * lbl;
655 7617 : dlimited_db = s * dlbl;
656 : }
657 : return;
658 : }
659 23550 : case FluxLimiterTypeEnum::VanLeer:
660 : {
661 23550 : limited = s * 2 * lal * lbl / (lal + lbl);
662 23550 : dlimited_db = s * 2 * lal * (dlbl / (lal + lbl) - lbl * dlbl / std::pow(lal + lbl, 2));
663 23550 : return;
664 : }
665 11174 : case FluxLimiterTypeEnum::MC:
666 : {
667 11174 : const Real av = 0.5 * std::abs(a + b);
668 11174 : if (2 * lal <= av && lal <= lbl)
669 : {
670 : // 2 * lal is the smallest
671 4812 : limited = s * 2.0 * lal;
672 4812 : dlimited_db = 0.0;
673 : }
674 6362 : else if (2 * lbl <= av && lbl <= lal)
675 : {
676 : // 2 * lbl is the smallest
677 2912 : limited = s * 2.0 * lbl;
678 2912 : dlimited_db = s * 2.0 * dlbl;
679 : }
680 : else
681 : {
682 : // av is the smallest
683 3450 : 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 3450 : dlimited_db = s * 0.5 * dlbl;
688 : }
689 : return;
690 : }
691 2274718 : case FluxLimiterTypeEnum::superbee:
692 : {
693 2274718 : const Real term1 = std::min(2.0 * lal, lbl);
694 2274718 : const Real term2 = std::min(lal, 2.0 * lbl);
695 2274718 : if (term1 >= term2)
696 : {
697 759321 : if (2.0 * lal <= lbl)
698 : {
699 407235 : limited = s * 2 * lal;
700 407235 : dlimited_db = 0.0;
701 : }
702 : else
703 : {
704 352086 : limited = s * lbl;
705 352086 : dlimited_db = s * dlbl;
706 : }
707 : }
708 : else
709 : {
710 1515397 : if (lal <= 2.0 * lbl)
711 : {
712 947290 : limited = s * lal;
713 947290 : dlimited_db = 0.0;
714 : }
715 : else
716 : {
717 568107 : limited = s * 2.0 * lbl;
718 568107 : 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 3368952 : AdvectiveFluxCalculatorBase::getdFluxOutdu(dof_id_type node_i) const
730 : {
731 3368952 : return _dflux_out_du[_connections.sequentialID(node_i)];
732 : }
733 :
734 : const std::vector<std::vector<Real>> &
735 1409604 : AdvectiveFluxCalculatorBase::getdFluxOutdKjk(dof_id_type node_i) const
736 : {
737 1409604 : return _dflux_out_dKjk[_connections.sequentialID(node_i)];
738 : }
739 :
740 : Real
741 5570340 : AdvectiveFluxCalculatorBase::getFluxOut(dof_id_type node_i) const
742 : {
743 5570340 : return _flux_out[_connections.sequentialID(node_i)];
744 : }
745 :
746 : unsigned
747 9407532 : AdvectiveFluxCalculatorBase::getValence(dof_id_type node_i) const
748 : {
749 9407532 : return _valence[_connections.sequentialID(node_i)];
750 : }
751 :
752 : void
753 56293 : AdvectiveFluxCalculatorBase::zeroedConnection(std::map<dof_id_type, Real> & the_map,
754 : dof_id_type node_i) const
755 : {
756 : the_map.clear();
757 333018 : for (const auto & node_j : _connections.globalConnectionsToGlobalID(node_i))
758 276725 : the_map[node_j] = 0.0;
759 56293 : }
760 :
761 : Real
762 8598631 : 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 8598631 : const Real u_i = _u_nodal[sequential_i];
769 :
770 : // Connections to sequential_i
771 : const std::vector<dof_id_type> con_i =
772 8598631 : _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 8598631 : const unsigned i_index_i = _connections.indexOfSequentialConnection(sequential_i, sequential_i);
776 :
777 : // Initialize the results
778 : Real result = 0.0;
779 8598631 : derivs.assign(num_con, 0.0);
780 8598631 : dpqdk.assign(num_con, 0.0);
781 :
782 : // Sum over all nodes connected with node_i.
783 47218486 : for (std::size_t j = 0; j < num_con; ++j)
784 : {
785 38619855 : const dof_id_type sequential_j = con_i[j];
786 38619855 : if (sequential_j == sequential_i)
787 8598631 : continue;
788 30021224 : const Real kentry = _kij[sequential_i][j];
789 :
790 : // Find the value of u at node_j
791 30021224 : const Real u_j = _u_nodal[sequential_j];
792 30021224 : const Real ujminusi = u_j - u_i;
793 :
794 : // Evaluate the i-j contribution to the result
795 30021224 : switch (pq_plus_minus)
796 : {
797 9604688 : case PQPlusMinusEnum::PPlus:
798 : {
799 9604688 : if (ujminusi < 0.0 && kentry < 0.0)
800 : {
801 2599157 : result += kentry * ujminusi;
802 2599157 : derivs[j] += kentry;
803 2599157 : derivs[i_index_i] -= kentry;
804 2599157 : dpqdk[j] += ujminusi;
805 : }
806 : break;
807 : }
808 9604688 : case PQPlusMinusEnum::PMinus:
809 : {
810 9604688 : if (ujminusi > 0.0 && kentry < 0.0)
811 : {
812 1910917 : result += kentry * ujminusi;
813 1910917 : derivs[j] += kentry;
814 1910917 : derivs[i_index_i] -= kentry;
815 1910917 : dpqdk[j] += ujminusi;
816 : }
817 : break;
818 : }
819 6125554 : case PQPlusMinusEnum::QPlus:
820 : {
821 6125554 : if (ujminusi > 0.0 && kentry > 0.0)
822 : {
823 2053203 : result += kentry * ujminusi;
824 2053203 : derivs[j] += kentry;
825 2053203 : derivs[i_index_i] -= kentry;
826 2053203 : dpqdk[j] += ujminusi;
827 : }
828 : break;
829 : }
830 4686294 : case PQPlusMinusEnum::QMinus:
831 : {
832 4686294 : if (ujminusi < 0.0 && kentry > 0.0)
833 : {
834 1380360 : result += kentry * ujminusi;
835 1380360 : derivs[j] += kentry;
836 1380360 : derivs[i_index_i] -= kentry;
837 1380360 : dpqdk[j] += ujminusi;
838 : }
839 : break;
840 : }
841 : }
842 : }
843 :
844 8598631 : return result;
845 8598631 : }
846 :
847 : void
848 2028 : 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 33142 : for (const auto & elem : _fe_problem.getNonlinearEvaluableElementRange())
865 31114 : if (this->hasBlocks(elem->subdomain_id()))
866 : {
867 28678 : const processor_id_type elem_pid = elem->processor_id();
868 28678 : if (elem_pid != _my_pid)
869 : {
870 6292 : if (_nodes_to_receive.find(elem_pid) == _nodes_to_receive.end())
871 3992 : _nodes_to_receive[elem_pid] = std::vector<dof_id_type>();
872 26068 : for (unsigned i = 0; i < elem->n_nodes(); ++i)
873 19776 : if (std::find(_nodes_to_receive[elem_pid].begin(),
874 19776 : _nodes_to_receive[elem_pid].end(),
875 39552 : elem->node_id(i)) == _nodes_to_receive[elem_pid].end())
876 11238 : _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 1996 : auto nodes_action_functor = [this](processor_id_type pid, const std::vector<dof_id_type> & nts)
883 1996 : { _nodes_to_send[pid] = nts; };
884 2028 : 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 4024 : for (auto & kv : _nodes_to_send)
892 : {
893 1996 : const processor_id_type pid = kv.first;
894 : const std::size_t num_nodes = kv.second.size();
895 13234 : for (unsigned i = 0; i < num_nodes; ++i)
896 11238 : _nodes_to_send[pid][i] = _connections.sequentialID(_nodes_to_send[pid][i]);
897 : }
898 4024 : for (auto & kv : _nodes_to_receive)
899 : {
900 1996 : const processor_id_type pid = kv.first;
901 : const std::size_t num_nodes = kv.second.size();
902 13234 : for (unsigned i = 0; i < num_nodes; ++i)
903 11238 : _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 33142 : for (const auto & elem : _fe_problem.getNonlinearEvaluableElementRange())
909 31114 : if (this->hasBlocks(elem->subdomain_id()))
910 : {
911 28678 : const processor_id_type elem_pid = elem->processor_id();
912 28678 : if (elem_pid != _my_pid)
913 : {
914 6292 : if (_pairs_to_receive.find(elem_pid) == _pairs_to_receive.end())
915 3992 : _pairs_to_receive[elem_pid] = std::vector<std::pair<dof_id_type, dof_id_type>>();
916 26068 : for (unsigned i = 0; i < elem->n_nodes(); ++i)
917 93760 : for (unsigned j = 0; j < elem->n_nodes(); ++j)
918 : {
919 73984 : std::pair<dof_id_type, dof_id_type> the_pair(elem->node_id(i), elem->node_id(j));
920 73984 : if (std::find(_pairs_to_receive[elem_pid].begin(),
921 73984 : _pairs_to_receive[elem_pid].end(),
922 73984 : the_pair) == _pairs_to_receive[elem_pid].end())
923 55690 : _pairs_to_receive[elem_pid].push_back(the_pair);
924 : }
925 : }
926 : }
927 :
928 : _pairs_to_send.clear();
929 : auto pairs_action_functor =
930 1996 : [this](processor_id_type pid, const std::vector<std::pair<dof_id_type, dof_id_type>> & pts)
931 1996 : { _pairs_to_send[pid] = pts; };
932 2028 : 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 4024 : for (auto & kv : _pairs_to_send)
940 : {
941 1996 : const processor_id_type pid = kv.first;
942 : const std::size_t num_pairs = kv.second.size();
943 57686 : for (unsigned i = 0; i < num_pairs; ++i)
944 : {
945 55690 : _pairs_to_send[pid][i].second = _connections.indexOfGlobalConnection(
946 55690 : _pairs_to_send[pid][i].first, _pairs_to_send[pid][i].second);
947 55690 : _pairs_to_send[pid][i].first = _connections.sequentialID(_pairs_to_send[pid][i].first);
948 : }
949 : }
950 4024 : for (auto & kv : _pairs_to_receive)
951 : {
952 1996 : const processor_id_type pid = kv.first;
953 : const std::size_t num_pairs = kv.second.size();
954 57686 : for (unsigned i = 0; i < num_pairs; ++i)
955 : {
956 55690 : _pairs_to_receive[pid][i].second = _connections.indexOfGlobalConnection(
957 55690 : _pairs_to_receive[pid][i].first, _pairs_to_receive[pid][i].second);
958 55690 : _pairs_to_receive[pid][i].first = _connections.sequentialID(_pairs_to_receive[pid][i].first);
959 : }
960 : }
961 2028 : }
962 :
963 : void
964 77980 : 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 155460 : for (const auto & kv : _nodes_to_send)
969 : {
970 77480 : const processor_id_type pid = kv.first;
971 154960 : unodal_to_send[pid] = std::vector<Real>();
972 504836 : for (const auto & nd : kv.second)
973 427356 : unodal_to_send[pid].push_back(_u_nodal[nd]);
974 : }
975 :
976 : auto unodal_action_functor =
977 77480 : [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 504836 : for (unsigned i = 0; i < msg_size; ++i)
985 427356 : _u_nodal[_nodes_to_receive[pid][i]] = unodal_received[i];
986 77980 : };
987 77980 : 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 155460 : for (const auto & kv : _pairs_to_send)
992 : {
993 77480 : const processor_id_type pid = kv.first;
994 154960 : kij_to_send[pid] = std::vector<Real>();
995 2089804 : for (const auto & pr : kv.second)
996 2012324 : kij_to_send[pid].push_back(_kij[pr.first][pr.second]);
997 : }
998 :
999 77480 : 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 2089804 : for (unsigned i = 0; i < msg_size; ++i)
1007 2012324 : _kij[_pairs_to_receive[pid][i].first][_pairs_to_receive[pid][i].second] += kij_received[i];
1008 77980 : };
1009 77980 : Parallel::push_parallel_vector_data(this->comm(), kij_to_send, kij_action_functor);
1010 77980 : }
|
