Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #pragma once
20 :
21 : #include "CardinalProblem.h"
22 : #include "NekInterface.h"
23 : #include "NekTimeStepper.h"
24 : #include "NekScalarValue.h"
25 : #include "NekRSMesh.h"
26 : #include "Transient.h"
27 :
28 : class FieldTransferBase;
29 :
30 : /**
31 : * \brief Solve nekRS wrapped as a MOOSE app.
32 : *
33 : * This object controls all of the execution of and data transfers to/from nekRS.
34 : * Any number of data transfers can be added using the [FieldTransfers] block.
35 : */
36 : class NekRSProblem : public CardinalProblem
37 : {
38 : public:
39 : NekRSProblem(const InputParameters & params);
40 :
41 : static InputParameters validParams();
42 :
43 : ~NekRSProblem();
44 :
45 : /**
46 : * The number of points (DOFs) on the mesh mirror for this rank; this is used to define the
47 : * size of various allocated terms
48 : * @return number of points on this rank
49 : */
50 1087 : int nPoints() const { return _n_points; }
51 :
52 : /**
53 : * Whether a data transfer to/from Nek is occurring
54 : * @param[in] direction direction of data transfer
55 : * @return whether a data transfer to Nek is about to occur
56 : */
57 : bool isDataTransferHappening(ExternalProblem::Direction direction);
58 :
59 537600 : Transient * transientExecutioner() const { return _transient_executioner; }
60 :
61 : /**
62 : * Whether a given slot space is reserved for coupling
63 : * @param[in] slot slot in usrwrk array
64 : * @return whether a usrwrk slot space is reserved by Cardinal
65 : */
66 : bool isUsrWrkSlotReservedForCoupling(const unsigned int & slot) const;
67 :
68 : /**
69 : * Copy an individual slice in the usrwrk array from host to device
70 : * @param[in] slot slot in usrwrk array
71 : */
72 : void copyIndividualScratchSlot(const unsigned int & slot) const;
73 :
74 : /**
75 : * Write NekRS solution field file
76 : * @param[in] time solution time in NekRS (if NekRS is non-dimensional, this will be
77 : * non-dimensional)
78 : * @param[in] step time step index
79 : */
80 : void writeFieldFile(const Real & time, const int & step) const;
81 :
82 : /**
83 : * \brief Whether nekRS should write an output file for the current time step
84 : *
85 : * A nekRS output file (suffix .f000xx) is written if the time step is an integer
86 : * multiple of the output writing interval or if the time step is the last time step.
87 : * \return whether to write a nekRS output file
88 : **/
89 : virtual bool isOutputStep() const;
90 :
91 : virtual void initialSetup() override;
92 :
93 : virtual void externalSolve() override;
94 :
95 33972 : virtual bool converged(unsigned int) override { return true; }
96 :
97 : virtual void addExternalVariables() override;
98 :
99 : virtual void syncSolutions(ExternalProblem::Direction direction) override;
100 :
101 : /**
102 : * Whether the solve is in nondimensional form
103 : * @return whether solve is in nondimensional form
104 : */
105 297 : virtual bool nondimensional() const { return _nondimensional; }
106 :
107 : /**
108 : * Whether the mesh is moving
109 : * @return whether the mesh is moving
110 : */
111 0 : virtual const bool hasMovingNekMesh() const { return false; }
112 :
113 : /**
114 : * Whether data should be synchronized in to nekRS
115 : * \return whether inward data synchronization should occur
116 : */
117 : virtual bool synchronizeIn();
118 :
119 : /**
120 : * Whether data should be synchronized out of nekRS
121 : * \return whether outward data synchronization should occur
122 : */
123 : virtual bool synchronizeOut();
124 :
125 : /**
126 : * Get the number of usrwrk slots allocated
127 : * @return number of allocated usrwrk slots
128 : */
129 606 : unsigned int nUsrWrkSlots() const { return _n_usrwrk_slots; }
130 :
131 : /**
132 : * Interpolate the NekRS volume solution onto the volume MOOSE mesh mirror (re2 -> mirror)
133 : * @param[in] f field to interpolate
134 : * @param[out] s interpolated volume value
135 : */
136 : template <typename T>
137 66712 : void volumeSolution(const T & field, double * s)
138 : {
139 66712 : mesh_t * mesh = nekrs::entireMesh();
140 66712 : auto vc = _nek_mesh->volumeCoupling();
141 :
142 : double (*f)(int, int);
143 66712 : f = nekrs::solutionPointer(field);
144 :
145 66712 : int start_1d = mesh->Nq;
146 66712 : int end_1d = _moose_Nq;
147 66712 : int start_3d = start_1d * start_1d * start_1d;
148 66712 : int end_3d = end_1d * end_1d * end_1d;
149 66712 : int n_to_write = vc.n_elems * end_3d * _nek_mesh->nBuildPerVolumeElem();
150 :
151 : // allocate temporary space:
152 : // - Ttmp: results of the search for each process
153 : // - Telem: scratch space for volume interpolation to avoid reallocating a bunch (only used if
154 : // interpolating)
155 66712 : double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
156 66712 : double * Telem = (double *)calloc(start_3d, sizeof(double));
157 : auto indices = _nek_mesh->cornerIndices();
158 :
159 : int c = 0;
160 2673034 : for (int k = 0; k < mesh->Nelements; ++k)
161 : {
162 2606322 : int offset = k * start_3d;
163 :
164 6139820 : for (int build = 0; build < _nek_mesh->nBuildPerVolumeElem(); ++build)
165 : {
166 3533498 : if (_needs_interpolation)
167 : {
168 : // get the solution on the element
169 240687826 : for (int v = 0; v < start_3d; ++v)
170 239667680 : Telem[v] = f(offset + v, 0 /* unused for volumes */);
171 :
172 : // and then interpolate it
173 1020146 : nekrs::interpolateVolumeHex3D(
174 1020146 : _interpolation_outgoing, Telem, start_1d, &(Ttmp[c]), end_1d);
175 1020146 : c += end_3d;
176 : }
177 : else
178 : {
179 : // get the solution on the element - no need to interpolate
180 22620168 : for (int v = 0; v < end_3d; ++v, ++c)
181 20106816 : Ttmp[c] = f(offset + indices[build][v], 0 /* unused for volumes */);
182 : }
183 : }
184 : }
185 :
186 : // dimensionalize the solution if needed
187 47717470 : for (int v = 0; v < n_to_write; ++v)
188 47650758 : Ttmp[v] =
189 47650758 : Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);
190 :
191 66712 : nekrs::allgatherv(vc.mirror_counts, Ttmp, s, end_3d);
192 :
193 : freePointer(Ttmp);
194 : freePointer(Telem);
195 66712 : }
196 :
197 : /**
198 : * Interpolate the NekRS boundary solution onto the boundary MOOSE mesh mirror (re2 -> mirror)
199 : * @param[in] f field to interpolate
200 : * @param[out] s interpolated boundary value
201 : */
202 : template <typename T>
203 20038 : void boundarySolution(const T & field, double * s)
204 : {
205 20038 : mesh_t * mesh = nekrs::entireMesh();
206 20038 : auto bc = _nek_mesh->boundaryCoupling();
207 :
208 : double (*f)(int, int);
209 20038 : f = nekrs::solutionPointer(field);
210 :
211 20038 : int start_1d = mesh->Nq;
212 20038 : int end_1d = _moose_Nq;
213 20038 : int start_2d = start_1d * start_1d;
214 20038 : int end_2d = end_1d * end_1d;
215 :
216 20038 : int n_to_write = bc.n_faces * end_2d * _nek_mesh->nBuildPerSurfaceElem();
217 :
218 : // allocate temporary space:
219 : // - Ttmp: results of the search for each process
220 : // - Tface: scratch space for face solution to avoid reallocating a bunch (only used if
221 : // interpolating)
222 : // - scratch: scratch for the interpolatino process to avoid reallocating a bunch (only used if
223 : // interpolating0
224 20038 : double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
225 20038 : double * Tface = (double *)calloc(start_2d, sizeof(double));
226 20038 : double * scratch = (double *)calloc(start_1d * end_1d, sizeof(double));
227 :
228 : auto indices = _nek_mesh->cornerIndices();
229 :
230 : int c = 0;
231 2151306 : for (int k = 0; k < bc.total_n_faces; ++k)
232 : {
233 2131268 : if (bc.process[k] == nekrs::commRank())
234 : {
235 368002 : int i = bc.element[k];
236 368002 : int j = bc.face[k];
237 368002 : int offset = i * mesh->Nfaces * start_2d + j * start_2d;
238 :
239 1104292 : for (int build = 0; build < _nek_mesh->nBuildPerSurfaceElem(); ++build)
240 : {
241 736290 : if (_needs_interpolation)
242 : {
243 : // get the solution on the face
244 3623648 : for (int v = 0; v < start_2d; ++v)
245 : {
246 3529344 : int id = mesh->vmapM[offset + v];
247 3529344 : Tface[v] = f(id, mesh->Nsgeo * (offset + v));
248 : }
249 :
250 : // and then interpolate it
251 94304 : nekrs::interpolateSurfaceFaceHex3D(
252 94304 : scratch, _interpolation_outgoing, Tface, start_1d, &(Ttmp[c]), end_1d);
253 94304 : c += end_2d;
254 : }
255 : else
256 : {
257 : // get the solution on the face - no need to interpolate
258 3209930 : for (int v = 0; v < end_2d; ++v, ++c)
259 : {
260 2567944 : int id = mesh->vmapM[offset + indices[build][v]];
261 2567944 : Ttmp[c] = f(id, mesh->Nsgeo * (offset + v));
262 : }
263 : }
264 : }
265 : }
266 : }
267 :
268 : // dimensionalize the solution if needed
269 3436718 : for (int v = 0; v < n_to_write; ++v)
270 3416680 : Ttmp[v] =
271 3416680 : Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);
272 :
273 20038 : nekrs::allgatherv(bc.mirror_counts, Ttmp, s, end_2d);
274 :
275 : freePointer(Ttmp);
276 : freePointer(Tface);
277 : freePointer(scratch);
278 20038 : }
279 :
280 : /**
281 : * Write a mesh displacement into NekRS
282 : * @param[in] elem_id element ID
283 : * @param[in] s solution values to write for the field for the given element
284 : * @param[in] f which component of the displacement to write
285 : * @param[in] add optional vector of values to add to each value set on the NekRS end
286 : */
287 : void writeVolumeDisplacement(const int elem_id,
288 : double * s,
289 : const field::NekWriteEnum f,
290 : const std::vector<double> * add = nullptr);
291 : /**
292 : * Write into the NekRS solution space for coupling volumes; for setting a mesh position in terms
293 : * of a displacement, we need to add the displacement to the initial mesh coordinates. For this,
294 : * the 'add' parameter lets you pass in a vector of values (in NekRS's mesh order, i.e. the re2
295 : * order) to add.
296 : * @param[in] slot the slot in the usrwrk array to populate
297 : * @param[in] elem_id element ID
298 : * @param[in] s solution values to write for the field for the given element
299 : * @param[in] add optional vector of values to add to each value set on the NekRS end
300 : */
301 : void writeVolumeSolution(const int slot,
302 : const int elem_id,
303 : double * s,
304 : const std::vector<double> * add = nullptr);
305 :
306 : /**
307 : * Write into the NekRS solution space for coupling boundaries; for setting a mesh position in
308 : * terms of a displacement, we need to add the displacement to the initial mesh coordinates.
309 : * @param[in] slot the slot in the usrwrk array to populate
310 : * @param[in] elem_id element ID
311 : * @param[in] s solution values to write for the field for the given element
312 : */
313 : void writeBoundarySolution(const int slot, const int elem_id, double * s);
314 :
315 : /**
316 : * The casename (prefix) of the NekRS files
317 : * @return casename
318 : */
319 26 : std::string casename() const { return _casename; }
320 :
321 : /**
322 : * Map nodal points on a MOOSE face element to the GLL points on a Nek face element.
323 : * @param[in] e MOOSE element ID
324 : * @param[in] var_num variable index to fetch MOOSE data from
325 : * @param[in] divisor number to divide MOOSE data by before sending to Nek (to non-dimensionalize
326 : * it)
327 : * @param[in] additive number to subtract from MOOSE data, before dividing by divisor and sending
328 : * to Nek (to non-dimensionalize)
329 : * @param[out] outgoing_data data represented on Nek's GLL points, ready to be applied in Nek
330 : */
331 : void mapFaceDataToNekFace(const unsigned int & e,
332 : const unsigned int & var_num,
333 : const Real & divisor,
334 : const Real & additive,
335 : double ** outgoing_data);
336 :
337 : /**
338 : * Map nodal points on a MOOSE volume element to the GLL points on a Nek volume element.
339 : * @param[in] e MOOSE element ID
340 : * @param[in] var_num variable index to fetch MOOSE data from
341 : * @param[in] divisor number to divide MOOSE data by before sending to Nek (to non-dimensionalize
342 : * it)
343 : * @param[in] additive number to subtract from MOOSE data, before dividing by divisor and sending
344 : * to Nek (to non-dimensionalize)
345 : * @param[out] outgoing_data data represented on Nek's GLL points, ready to be applied in Nek
346 : */
347 : void mapVolumeDataToNekVolume(const unsigned int & e,
348 : const unsigned int & var_num,
349 : const Real & divisor,
350 : const Real & additive,
351 : double ** outgoing_data);
352 :
353 : /**
354 : * \brief Map nodal points on a MOOSE face element to the GLL points on a Nek volume element.
355 : *
356 : * This function is to be used when MOOSE variables are defined over the entire volume
357 : * (maybe the MOOSE transfer only sent meaningful values to the coupling boundaries), so we
358 : * need to do a volume interpolation of the incoming MOOSE data into nrs->usrwrk, rather
359 : * than a face interpolation. This could be optimized in the future to truly only just write
360 : * the boundary values into the nekRS scratch space rather than the volume values, but it
361 : * looks right now that our biggest expense occurs in the MOOSE transfer system, not these
362 : * transfers internally to nekRS.
363 : *
364 : * @param[in] e MOOSE element ID
365 : * @param[in] var_num variable index to fetch MOOSE data from
366 : * @param[in] divisor number to divide MOOSE data by before sending to Nek (to non-dimensionalize
367 : * it)
368 : * @param[in] additive number to subtract from MOOSE data, before dividing by divisor and sending
369 : * to Nek (to non-dimensionalize)
370 : * @param[out] outgoing_data data represented on Nek's GLL points, ready to be applied in Nek
371 : */
372 : void mapFaceDataToNekVolume(const unsigned int & e,
373 : const unsigned int & var_num,
374 : const Real & divisor,
375 : const Real & additive,
376 : double ** outgoing_data);
377 :
378 : protected:
379 : /// Copy the data sent from MOOSE->Nek from host to device.
380 : void copyScratchToDevice();
381 :
382 : /**
383 : * Interpolate the MOOSE mesh mirror solution onto the NekRS boundary mesh (mirror -> re2)
384 : * @param[in] incoming_moose_value MOOSE face values
385 : * @param[out] outgoing_nek_value interpolated MOOSE face values onto the NekRS boundary mesh
386 : */
387 : void interpolateBoundarySolutionToNek(double * incoming_moose_value, double * outgoing_nek_value);
388 :
389 : /**
390 : * Interpolate the MOOSE mesh mirror solution onto the NekRS volume mesh (mirror -> re2)
391 : * @param[in] elem_id element ID
392 : * @param[in] incoming_moose_value MOOSE face values
393 : * @param[out] outgoing_nek_value interpolated MOOSE face values onto the NekRS volume mesh
394 : */
395 : void interpolateVolumeSolutionToNek(const int elem_id,
396 : double * incoming_moose_value,
397 : double * outgoing_nek_value);
398 :
399 : /// Initialize interpolation matrices for transfers in/out of nekRS
400 : void initializeInterpolationMatrices();
401 :
402 : std::unique_ptr<NumericVector<Number>> _serialized_solution;
403 :
404 : /**
405 : * Get a three-character prefix for use in writing output files for repeated
406 : * Nek sibling apps.
407 : * @param[in] number multi-app number
408 : */
409 : std::string fieldFilePrefix(const int & number) const;
410 :
411 : /// NekRS casename
412 : const std::string & _casename;
413 :
414 : /**
415 : * Whether to disable output file writing by NekRS and replace it by output
416 : * file writing in Cardinal. Suppose the case name is 'channel'. If this parameter
417 : * is false, then NekRS will write output files as usual, with names like
418 : * channel0.f00001 for write step 1, channel0.f00002 for write step 2, and so on.
419 : * If true, then NekRS itself does not output any files like this, and instead
420 : * output files are written with names a01channel0.f00001, a01channel0.f00002 (for
421 : * first Nek app), a02channel0.f00001, a02channel0.f00002 (for second Nek app),
422 : * and so on. This feature should only be used when running repeated Nek sub
423 : * apps so that the output from each app is retained. Otherwise, if running N
424 : * repeated Nek sub apps, only a single output file is obtained because each app
425 : * overwrites the output files of the other apps in the order that the apps
426 : * reach the nekrs::outfld function.
427 : */
428 : const bool & _write_fld_files;
429 :
430 : /// Whether to turn off all field file writing
431 : const bool & _disable_fld_file_output;
432 :
433 : /// Whether a dimensionalization action has been added
434 : bool _nondimensional;
435 :
436 : /**
437 : * Number of slices/slots to allocate in nrs->usrwrk to hold fields
438 : * for coupling (i.e. data going into NekRS, written by Cardinal), or
439 : * used for custom user actions, but not for coupling. By default, we just
440 : * allocate 7 slots (no inherent reason, just a fairly big amount). For
441 : * memory-limited cases, you can reduce this number to just the bare
442 : * minimum necessary for your use case.
443 : */
444 : unsigned int _n_usrwrk_slots;
445 :
446 : /// For constant synchronization intervals, the desired frequency (in units of Nek time steps)
447 : const unsigned int & _constant_interval;
448 :
449 : /// Whether to skip writing a field file on NekRS's last time steo
450 : const bool & _skip_final_field_file;
451 :
452 : /// Number of surface elements in the data transfer mesh, across all processes
453 : int _n_surface_elems;
454 :
455 : /// Number of vertices per surface element of the transfer mesh
456 : int _n_vertices_per_surface;
457 :
458 : /// Number of volume elements in the data transfer mesh, across all processes
459 : int _n_volume_elems;
460 :
461 : /// Number of vertices per volume element of the transfer mesh
462 : int _n_vertices_per_volume;
463 :
464 : /// Start time of the simulation based on NekRS's .par file
465 : double _start_time;
466 :
467 : /// Whether the most recent time step was an output file writing step
468 : bool _is_output_step;
469 :
470 : /**
471 : * Underlying mesh object on which NekRS exchanges fields with MOOSE
472 : * or extracts NekRS's solution for I/O features
473 : */
474 : NekRSMesh * _nek_mesh;
475 :
476 : /// The time stepper used for selection of time step size
477 : NekTimeStepper * _timestepper = nullptr;
478 :
479 : /// Underlying executioner
480 : Transient * _transient_executioner = nullptr;
481 :
482 : /**
483 : * Whether an interpolation needs to be performed on the nekRS temperature solution, or
484 : * if we can just grab the solution at specified points
485 : */
486 : bool _needs_interpolation;
487 :
488 : /// Number of points for interpolated fields on the MOOSE mesh
489 : int _n_points;
490 :
491 : /// Postprocessor containing the signal of when a synchronization has occurred
492 : const PostprocessorValue * _transfer_in = nullptr;
493 :
494 : /// Vandermonde interpolation matrix (for outgoing transfers)
495 : double * _interpolation_outgoing = nullptr;
496 :
497 : /// Vandermonde interpolation matrix (for incoming transfers)
498 : double * _interpolation_incoming = nullptr;
499 :
500 : /// For the MOOSE mesh, the number of quadrature points in each coordinate direction
501 : int _moose_Nq;
502 :
503 : /// Slots in the nrs->o_usrwrk array to write to a field file
504 : const std::vector<unsigned int> * _usrwrk_output = nullptr;
505 :
506 : /// Filename prefix to use for naming the field files containing the nrs->o_usrwrk array slots
507 : const std::vector<std::string> * _usrwrk_output_prefix = nullptr;
508 :
509 : /// Sum of the elapsed time in NekRS solves
510 : double _elapsedStepSum;
511 :
512 : /// Sum of the total elapsed time in NekRS solves
513 : double _elapsedTime;
514 :
515 : /// Minimum step solve time
516 : double _tSolveStepMin;
517 :
518 : /// Maximum step solve time
519 : double _tSolveStepMax;
520 :
521 : /**
522 : * \brief When to synchronize the NekRS solution with the mesh mirror
523 : *
524 : * This parameter determines when to synchronize the NekRS solution with the mesh
525 : * mirror - this entails:
526 : *
527 : * - Mapping from the NekRS spectral element mesh to the finite element mesh mirror,
528 : * to extract information from NekRS and make it available to MOOSE
529 : * - Mapping from the finite element mesh mirror into the NekRS spectral element mesh,
530 : * to send information from MOOSE into NekRS
531 : *
532 : * Several options are available:
533 : * - 'constant' will simply keep the NekRS solution and the mesh mirror entirely
534 : * consistent with one another on a given constant frequency of time steps. By
535 : * default, the 'constant_interval' is 1, so that NekRS and MOOSE communicate
536 : * with each other on every single time step
537 : *
538 : * - 'parent_app' will only send data between NekRS and a parent application
539 : * when (1) the main application has just sent "new" information to NekRS, and
540 : * when (2) the main application is just about to run a new time step (with
541 : * updated BCs/source terms from NekRS).
542 : *
543 : * nekRS is often subcycled relative to the application controlling it -
544 : * that is, nekRS may be run with a time step 10x smaller than a conduction MOOSE app.
545 : * If 'interpolate_transfers = false'
546 : * in the master application, then the data going into nekRS is fixed for each
547 : * of the subcycled time steps it takes, so these extra data transfers are
548 : * completely unnecssary. This flag indicates that the information sent from MOOSE
549 : * to NekRS should only be updated if the data from MOOSE is "new", and likewise
550 : * whether the NekRS solution should only be interpolated to the mesh mirror once
551 : * MOOSE is actually "ready" to solve a time step using it.
552 : *
553 : * NOTE: if 'interpolate_transfers = true' in the master application, then the data
554 : * coming into nekRS is _unique_ on each subcycled time step, so setting this to
555 : * true will in effect override `interpolate_transfers` to be false. For the best
556 : * performance, you should set `interpolate_transfers` to false so that you don't
557 : * even bother computing the interpolated data, since it's not used if this parameter
558 : * is set to true.
559 : */
560 : synchronization::SynchronizationEnum _sync_interval;
561 :
562 : /// flag to indicate whether this is the first pass to serialize the solution
563 : static bool _first;
564 :
565 : /// All of the FieldTransfer objects which pass data in/out of NekRS
566 : std::vector<FieldTransferBase *> _field_transfers;
567 :
568 : /// All of the ScalarTransfer objecst which pass data in/out of NekRS
569 : std::vector<ScalarTransferBase *> _scalar_transfers;
570 :
571 : /// Usrwrk slots managed by Cardinal
572 : std::set<unsigned int> _usrwrk_slots;
573 : };
|
