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 : #ifdef ENABLE_NEK_COUPLING
20 :
21 : #include "NekRSProblem.h"
22 : #include "CardinalUtils.h"
23 : #include "DimensionalizeAction.h"
24 : #include "FieldTransferBase.h"
25 : #include "NekMeshDeformation.h"
26 :
27 : registerMooseObject("CardinalApp", NekRSProblem);
28 :
29 : bool NekRSProblem::_first = true;
30 :
31 : InputParameters
32 2279 : NekRSProblem::validParams()
33 : {
34 2279 : InputParameters params = CardinalProblem::validParams();
35 4558 : params.addRequiredParam<std::string>(
36 : "casename",
37 : "Case name for the NekRS input files; "
38 : "this is <case> in <case>.par, <case>.udf, <case>.oudf, and <case>.re2.");
39 :
40 4558 : params.addParam<unsigned int>(
41 : "n_usrwrk_slots",
42 4558 : 0,
43 : "Number of slots to allocate in platform->app->bc->o_usrwrk to hold fields either related to "
44 : "coupling "
45 : "(which will be populated by Cardinal), or other custom usages, such as a distance-to-wall "
46 : "calculation (which will be populated by the user from the case files)");
47 :
48 4558 : params.addParam<std::vector<unsigned int>>(
49 : "usrwrk_output",
50 : "Usrwrk slot(s) to output to NekRS field files; this can be used for viewing the quantities "
51 : "passed from MOOSE to NekRS after interpolation to the CFD mesh. Can also be used for any "
52 : "slots "
53 : "in usrwrk that are written by the user, but unused for coupling.");
54 4558 : params.addParam<std::vector<std::string>>(
55 : "usrwrk_output_prefix",
56 : "String prefix to use for naming the field file(s); "
57 : "only the first three characters are used in the name based on limitations in NekRS");
58 :
59 4558 : params.addParam<bool>(
60 : "write_fld_files",
61 4558 : false,
62 : "Whether to write NekRS field file output "
63 : "from Cardinal. If true, this will disable any output writing by NekRS itself, and "
64 : "instead produce output files with names a01...a99pin, b01...b99pin, etc.");
65 4558 : params.addParam<bool>(
66 4558 : "disable_fld_file_output", false, "Whether to turn off all NekRS field file output writing");
67 :
68 4558 : params.addParam<bool>("skip_final_field_file",
69 4558 : false,
70 : "By default, we write a NekRS field file "
71 : "on the last time step; set this to true to disable");
72 :
73 4558 : params.addParam<MooseEnum>(
74 : "synchronization_interval",
75 4558 : getSynchronizationEnum(),
76 : "When to synchronize the NekRS solution with the mesh mirror. By default, the NekRS solution "
77 : "is mapped to/receives data from the mesh mirror for every time step.");
78 4558 : params.addParam<unsigned int>("constant_interval",
79 4558 : 1,
80 : "Constant interval (in units of number of time steps) with which "
81 : "to synchronize the NekRS solution");
82 2279 : return params;
83 0 : }
84 :
85 749 : NekRSProblem::NekRSProblem(const InputParameters & params)
86 : : CardinalProblem(params),
87 749 : _serialized_solution(NumericVector<Number>::build(_communicator).release()),
88 1498 : _casename(getParam<std::string>("casename")),
89 1498 : _write_fld_files(getParam<bool>("write_fld_files")),
90 1498 : _disable_fld_file_output(getParam<bool>("disable_fld_file_output")),
91 1498 : _n_usrwrk_slots(getParam<unsigned int>("n_usrwrk_slots")),
92 1498 : _constant_interval(getParam<unsigned int>("constant_interval")),
93 1498 : _skip_final_field_file(getParam<bool>("skip_final_field_file")),
94 749 : _start_time(nekrs::startTime()),
95 749 : _elapsedStepSum(0.0),
96 1498 : _elapsedTime(nekrs::getNekSetupTime()),
97 749 : _tSolveStepMin(std::numeric_limits<double>::max()),
98 1498 : _tSolveStepMax(std::numeric_limits<double>::min())
99 : {
100 : // fetch data from device to host; this needs to be here in case there are any Cardinal objects
101 : // executing on 'initial', so that this will execute even before NekRS has run any time steps
102 749 : nekrs::copySolutionToHost();
103 :
104 749 : const auto & actions = getMooseApp().actionWarehouse().getActions<DimensionalizeAction>();
105 749 : _nondimensional = actions.size();
106 749 : nekrs::nondimensional(_nondimensional);
107 :
108 1498 : _sync_interval = getParam<MooseEnum>("synchronization_interval")
109 : .getEnum<synchronization::SynchronizationEnum>();
110 749 : if (_sync_interval == synchronization::parent_app)
111 : {
112 : // the way the data transfers are detected depend on nekRS being a sub-application,
113 : // so these settings are not invalid if nekRS is the master app (though you could
114 : // relax this in the future by reversing the synchronization step identification
115 : // from the nekRS-subapp case to the nekRS-master app case - it's just not implemented yet).
116 21 : if (_app.isUltimateMaster())
117 : {
118 0 : mooseWarning("The 'synchronization_interval = parent_app' capability "
119 : "requires that nekRS is receiving and sending data to a parent application, but "
120 : "in your case nekRS is the main application.\n\n"
121 : "We are reverting synchronization_interval to 'constant'.");
122 0 : _sync_interval = synchronization::constant;
123 : }
124 :
125 42 : checkUnusedParam(params, "constant_interval", "synchronizing based on the 'parent_app'");
126 : }
127 :
128 749 : if (_disable_fld_file_output && _write_fld_files)
129 0 : mooseError("Cannot both disable all field file output and write custom field files! "
130 : "'write_fld_files' and 'disable_fld_file_output' cannot both be true!");
131 :
132 749 : if (_app.isUltimateMaster() && _write_fld_files)
133 1 : mooseError("The 'write_fld_files' setting should only be true when multiple Nek simulations "
134 : "are run as sub-apps on a master app. Your input has Nek as the master app.");
135 :
136 748 : _nek_mesh = dynamic_cast<NekRSMesh *>(&mesh());
137 :
138 748 : if (!_nek_mesh)
139 2 : mooseError("The mesh for NekRSProblem must be of type 'NekRSMesh', but you have specified a '" +
140 0 : mesh().type() + "'!");
141 :
142 : // The mesh movement error checks are triggered based on whether the NekRS input files
143 : // have a moving mesh. From there, we impose the necessary checks on the [Mesh] block
144 : // and the existence of the NekMeshDeformation object.
145 747 : if (nekrs::hasMovingMesh())
146 : {
147 4 : if (!_nek_mesh->getMesh().is_replicated())
148 0 : mooseError("Distributed mesh features are not yet implemented for moving mesh cases!");
149 : }
150 :
151 747 : _moose_Nq = _nek_mesh->order() + 2;
152 :
153 : // the Problem constructor is called right after building the mesh. In order
154 : // to have pretty screen output without conflicting with the timed print messages,
155 : // print diagnostic info related to the mesh here. If running in JIT mode, this
156 : // diagnostic info was never set, so the numbers that would be printed are garbage.
157 747 : if (!nekrs::buildOnly())
158 747 : _nek_mesh->printMeshInfo();
159 :
160 : // boundary-specific data
161 747 : _n_surface_elems = _nek_mesh->numSurfaceElems();
162 747 : _n_vertices_per_surface = _nek_mesh->numVerticesPerSurface();
163 :
164 : // volume-specific data
165 747 : _n_vertices_per_volume = _nek_mesh->numVerticesPerVolume();
166 :
167 747 : if (_nek_mesh->volume())
168 446 : _n_points =
169 446 : _nek_mesh->numVolumeElems() * _n_vertices_per_volume * _nek_mesh->nBuildPerVolumeElem();
170 : else
171 301 : _n_points = _n_surface_elems * _n_vertices_per_surface * _nek_mesh->nBuildPerSurfaceElem();
172 :
173 747 : initializeInterpolationMatrices();
174 :
175 : // we can save some effort for the low-order situations where the interpolation
176 : // matrix is the identity matrix (i.e. for which equi-spaced libMesh nodes are an
177 : // exact subset of the nekRS GLL points). This will happen for any first-order mesh,
178 : // and if a second-order mesh is used with a polynomial order of 2 in nekRS. Because
179 : // we pretty much always use a polynomial order greater than 2 in nekRS, let's just
180 : // check the first case because this will simplify our code in the nekrs::boundarySolution
181 : // function. If you change this line, you MUST change the innermost if/else statement
182 : // in nekrs::boundarySolution!
183 747 : _needs_interpolation = _nek_mesh->numQuadraturePoints1D() > 2;
184 :
185 747 : checkJointParams(
186 : params, {"usrwrk_output", "usrwrk_output_prefix"}, "outputting usrwrk slots to field files");
187 :
188 1494 : if (isParamValid("usrwrk_output"))
189 : {
190 40 : _usrwrk_output = &getParam<std::vector<unsigned int>>("usrwrk_output");
191 40 : _usrwrk_output_prefix = &getParam<std::vector<std::string>>("usrwrk_output_prefix");
192 :
193 52 : for (const auto & s : *_usrwrk_output)
194 33 : if (s >= _n_usrwrk_slots)
195 1 : mooseError("Cannot write field file for usrwrk slot greater than the total number of "
196 : "allocated slots: ",
197 1 : _n_usrwrk_slots,
198 : "! Please increase 'n_usrwrk_slots'.");
199 :
200 19 : if (_usrwrk_output->size() != _usrwrk_output_prefix->size())
201 1 : mooseError("The length of 'usrwrk_output' must match the length of 'usrwrk_output_prefix'!");
202 : }
203 745 : }
204 :
205 678 : NekRSProblem::~NekRSProblem()
206 : {
207 : // write nekRS solution to output if not already written for this step; nekRS does this
208 : // behavior, so we duplicate it
209 678 : if (!_is_output_step && !_skip_final_field_file)
210 : {
211 134 : if (_write_fld_files)
212 0 : mooseWarning(
213 : "When 'write_fld_files' is enabled, we skip Nek field file writing on end time!\n"
214 : "Depending on how many ranks you used, MOOSE may use the same object to run multiple\n"
215 : "sub-applications. By the time we get to the last time step, we've collapsed back to\n"
216 : "this singular state and don't have access to the individual Nek solves, so we cannot\n"
217 : "write the last time step solution to field files.\n\n"
218 : "To hide this warning, set 'skip_final_field_file = true'.");
219 : else
220 134 : writeFieldFile(_time, _t_step);
221 : }
222 :
223 678 : if (nekrs::runTimeStatFreq())
224 678 : if (_t_step % nekrs::runTimeStatFreq())
225 678 : nekrs::printRuntimeStatistics(_t_step);
226 :
227 678 : freePointer(_interpolation_outgoing);
228 678 : freePointer(_interpolation_incoming);
229 678 : nekrs::freeScratch();
230 :
231 678 : nekrs::finalize();
232 678 : }
233 :
234 : void
235 1945 : NekRSProblem::writeFieldFile(const Real & step_end_time, const int & step) const
236 : {
237 1945 : if (_disable_fld_file_output)
238 : return;
239 :
240 1945 : Real t = _timestepper->nondimensionalDT(step_end_time);
241 :
242 1945 : if (_write_fld_files)
243 : {
244 : // this is the app number, but a single app may run Nek multiple times
245 56 : auto app_number = std::to_string(_app.multiAppNumber());
246 :
247 : // apps may also have numbers in their names, so we first need to get the actual raw app name;
248 : // we strip out the app_number from the end of the app name
249 56 : if (!stringHasEnding(_app.name(), app_number))
250 0 : mooseError("Internal error: app name '" + _app.name() +
251 0 : "' does not end with app number: " + app_number);
252 :
253 56 : auto name = _app.name().substr(0, _app.name().size() - app_number.size());
254 56 : auto full_path = _app.getOutputFileBase();
255 56 : std::string last_element(full_path.substr(full_path.rfind(name) + name.size()));
256 :
257 112 : auto prefix = fieldFilePrefix(std::stoi(last_element)) + casename();
258 :
259 56 : nekrs::write_field_file(prefix, t, step);
260 : }
261 : else
262 1889 : nekrs::writeCheckpoint(t);
263 : }
264 :
265 : void
266 747 : NekRSProblem::initializeInterpolationMatrices()
267 : {
268 747 : mesh_t * mesh = nekrs::entireMesh();
269 :
270 : // determine the interpolation matrix for the outgoing transfer
271 747 : int starting_points = mesh->Nq;
272 747 : int ending_points = _nek_mesh->numQuadraturePoints1D();
273 747 : _interpolation_outgoing = (double *)calloc(starting_points * ending_points, sizeof(double));
274 747 : nekrs::interpolationMatrix(_interpolation_outgoing, starting_points, ending_points);
275 :
276 : // determine the interpolation matrix for the incoming transfer
277 : std::swap(starting_points, ending_points);
278 747 : _interpolation_incoming = (double *)calloc(starting_points * ending_points, sizeof(double));
279 747 : nekrs::interpolationMatrix(_interpolation_incoming, starting_points, ending_points);
280 747 : }
281 :
282 : std::string
283 56 : NekRSProblem::fieldFilePrefix(const int & number) const
284 : {
285 56 : const std::string alphabet = "abcdefghijklmnopqrstuvwxyz";
286 56 : int letter = number / 26;
287 56 : int remainder = number % 100;
288 56 : std::string s = remainder < 10 ? "0" : "";
289 :
290 168 : return alphabet[letter] + s + std::to_string(remainder);
291 : }
292 :
293 : void
294 688 : NekRSProblem::initialSetup()
295 : {
296 688 : CardinalProblem::initialSetup();
297 :
298 688 : auto executioner = _app.getExecutioner();
299 688 : _transient_executioner = dynamic_cast<Transient *>(executioner);
300 :
301 : // NekRS only supports transient simulations - therefore, it does not make
302 : // sense to use anything except a Transient-derived executioner
303 688 : if (!_transient_executioner)
304 2 : mooseError(
305 1 : "A 'Transient' executioner must be used with NekRSProblem, but you have specified the '" +
306 : executioner->type() + "' executioner!");
307 :
308 : // To get the correct time stepping information on the MOOSE side, we also
309 : // must use the NekTimeStepper
310 : TimeStepper * stepper = _transient_executioner->getTimeStepper();
311 687 : _timestepper = dynamic_cast<NekTimeStepper *>(stepper);
312 687 : if (!_timestepper)
313 2 : mooseError("The 'NekTimeStepper' stepper must be used with NekRSProblem, but you have "
314 1 : "specified the '" +
315 : stepper->type() + "' time stepper!");
316 :
317 : // Set the NekRS start time to whatever is set on Executioner/start_time; print
318 : // a message if those times don't match the .par file
319 686 : const auto moose_start_time = _transient_executioner->getStartTime();
320 686 : nekrs::setStartTime(_timestepper->nondimensionalDT(moose_start_time));
321 686 : _start_time = moose_start_time;
322 :
323 686 : if (_sync_interval == synchronization::parent_app)
324 21 : _transfer_in = &getPostprocessorValueByName("transfer_in");
325 :
326 : // Find all of the data transfer objects
327 686 : TheWarehouse::Query query = theWarehouse().query().condition<AttribSystem>("FieldTransfer");
328 686 : query.queryInto(_field_transfers);
329 :
330 : // Find all of the scalar data transfer objects
331 686 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("ScalarTransfer");
332 686 : uo_query.queryInto(_scalar_transfers);
333 :
334 : // We require a NekMeshDeformation object to exist if the NekRS model has a moving mesh
335 686 : if (nekrs::hasMovingMesh())
336 : {
337 : bool has_deformation = false;
338 3 : for (const auto & t : _field_transfers)
339 : {
340 2 : NekMeshDeformation * deform = dynamic_cast<NekMeshDeformation *>(t);
341 2 : if (deform)
342 : has_deformation = true;
343 : }
344 :
345 1 : if (has_deformation && !_app.actionWarehouse().displacedMesh())
346 1 : mooseError("Moving mesh problems require 'displacements' in the [Mesh] block! The names of "
347 : "the 'displacements' variables must match the variables created by a "
348 : "NekMeshDeformation object.");
349 : }
350 :
351 : // save initial mesh for moving mesh problems to match deformation in exodus output files
352 685 : if (nekrs::hasMovingMesh() && !_disable_fld_file_output)
353 0 : nekrs::writeCheckpoint(_timestepper->nondimensionalDT(_time));
354 :
355 : VariadicTable<int, std::string, std::string, std::string> vt(
356 685 : {"Slot", "Data", "How to Access (.oudf)", "How to Access (.udf)"});
357 :
358 : // fill a set with all of the slots managed by Cardinal, coming from either field transfers
359 : // or userobjects
360 : auto field_usrwrk_map = FieldTransferBase::usrwrkMap();
361 : auto field_usrwrk_scales = FieldTransferBase::usrwrkScales();
362 1053 : for (const auto & field : field_usrwrk_map)
363 368 : _usrwrk_slots.insert(field.first);
364 728 : for (const auto & uo : _scalar_transfers)
365 43 : _usrwrk_slots.insert(uo->usrwrkSlot());
366 :
367 : // fill out table, being careful to only write information if owned by a field transfer,
368 : // a user object, or neither
369 1190 : for (int i = 0; i < _n_usrwrk_slots; ++i)
370 : {
371 505 : std::string oudf = "bc->usrwrk[" + std::to_string(i) + "*bc->fieldOffset+bc->idxVol]";
372 505 : std::string udf = "platform->app->bc->o_usrwrk[" + std::to_string(i) + "*nrs->fieldOffset+n]";
373 :
374 505 : if (field_usrwrk_map.find(i) != field_usrwrk_map.end())
375 : {
376 : // a field transfer owns it
377 368 : auto scales = field_usrwrk_scales[field_usrwrk_map[i]];
378 :
379 : std::string top;
380 368 : if (MooseUtils::absoluteFuzzyEqual(scales.first, 0.0))
381 720 : top = field_usrwrk_map[i];
382 : else
383 16 : top = field_usrwrk_map[i] + "-" + std::to_string(scales.first);
384 :
385 368 : if (!MooseUtils::absoluteFuzzyEqual(scales.second, 1.0))
386 : {
387 28 : if (MooseUtils::absoluteFuzzyEqual(scales.first, 0.0))
388 60 : top = "(" + top + ")/" + std::to_string(scales.second);
389 : else
390 16 : top = top + "/" + std::to_string(scales.second);
391 : }
392 :
393 736 : vt.addRow(i, top, oudf, udf);
394 : }
395 : else
396 : {
397 : // a user object might own it, or it could be unused
398 : bool owned_by_uo = false;
399 222 : for (const auto & uo : _scalar_transfers)
400 : {
401 85 : if (uo->usrwrkSlot() == i)
402 : {
403 : owned_by_uo = true;
404 43 : auto slot = std::to_string(uo->usrwrkSlot());
405 43 : auto count = std::to_string(uo->offset());
406 :
407 43 : std::string top = uo->name();
408 43 : if (!MooseUtils::absoluteFuzzyEqual(uo->scaling(), 1.0))
409 0 : top += top + "*" + std::to_string(uo->scaling());
410 :
411 86 : vt.addRow(i,
412 : top,
413 129 : "bc->usrwrk[" + slot + "*bc->fieldOffset+" + count + "]",
414 129 : "platform->app->bc->o_usrwrk[" + slot + "*nrs->fieldOffset+" + count + "]");
415 : }
416 : }
417 :
418 137 : if (!owned_by_uo)
419 204 : vt.addRow(i, "unused", oudf, udf);
420 : }
421 : }
422 :
423 685 : if (_n_usrwrk_slots > 0)
424 : {
425 : _console
426 326 : << "\n ===================> MAPPING FROM MOOSE TO NEKRS <===================\n"
427 326 : << std::endl;
428 326 : _console << " Slot: slice in scratch space holding the data\n" << std::endl;
429 326 : _console << " Data: data that gets written into this slot. This data is shown"
430 326 : << std::endl;
431 326 : _console << " in the form actually written into NekRS (which will be"
432 326 : << std::endl;
433 326 : _console << " non-dimensional quantities if using the [Dimensionalize]"
434 326 : << std::endl;
435 326 : _console << " block). Words refer to MOOSE AuxVariables/Postprocessors."
436 326 : << std::endl;
437 326 : _console << " If 'unused', this means that the space has been allocated,"
438 326 : << std::endl;
439 326 : _console << " but Cardinal is not otherwise using it for coupling.\n"
440 326 : << std::endl;
441 326 : _console << " How to Access: C++ code to use in NekRS files; for the .udf instructions,"
442 326 : << std::endl;
443 326 : _console << " 'n' indicates a loop variable over GLL points\n" << std::endl;
444 326 : vt.print(_console);
445 326 : _console << std::endl;
446 : }
447 :
448 : // nekRS calls UDF_ExecuteStep once before the time stepping begins; the isLastStep stuff is
449 : // copy-pasta from NekRS main(), except that if Nek is a sub-app, we give full control of
450 : // time stepping to the main app
451 : bool isLastStep = false;
452 685 : if (_app.isUltimateMaster())
453 425 : isLastStep = !((nekrs::endTime() > nekrs::startTime() || nekrs::numSteps() > _t_step));
454 685 : nekrs::lastStep(isLastStep);
455 :
456 685 : nekrs::udfExecuteStep(
457 685 : _timestepper->nondimensionalDT(_start_time), _t_step, false /* not an output step */);
458 1370 : nekrs::resetTimer("udfExecuteStep");
459 685 : }
460 :
461 : void
462 21084 : NekRSProblem::externalSolve()
463 : {
464 21084 : if (nekrs::buildOnly())
465 0 : return;
466 :
467 21084 : const double timeStartStep = MPI_Wtime();
468 :
469 : // _dt reflects the time step that MOOSE wants Nek to
470 : // take. For instance, if Nek is controlled by a master app and subcycling is used,
471 : // Nek must advance to the time interval taken by the master app. If the time step
472 : // that MOOSE wants nekRS to take (i.e. _dt) is smaller than we'd like nekRS to take, error.
473 21084 : if (_dt < _timestepper->minDT())
474 0 : mooseError("Requested time step of " + std::to_string(_dt) +
475 : " is smaller than the minimum "
476 0 : "time step of " +
477 0 : Moose::stringify(_timestepper->minDT()) +
478 : " allowed in NekRS!\n\n"
479 : "You can control this behavior with the 'min_dt' parameter on 'NekTimeStepper'.");
480 :
481 : // _time represents the time that we're simulating _to_, but we need to pass sometimes slightly
482 : // different times into the nekRS routines, which assume that the "time" passed into their
483 : // routines is sometimes a different interpretation.
484 21084 : double step_start_time = _time - _dt;
485 21084 : double step_end_time = _time;
486 :
487 21084 : _is_output_step = isOutputStep();
488 :
489 : // tell NekRS what the value of nrs->isOutputStep should be
490 21084 : nekrs::checkpointStep(_is_output_step);
491 :
492 : // Tell NekRS what the time step size is
493 21084 : nekrs::initStep(_timestepper->nondimensionalDT(step_start_time),
494 21084 : _timestepper->nondimensionalDT(_dt),
495 21084 : _t_step);
496 :
497 : // Run a nekRS time step. After the time step, this also calls UDF_ExecuteStep,
498 : // evaluated at (step_end_time, _t_step) == (nek_step_start_time + nek_dt, t_step)
499 : int corrector = 1;
500 : bool converged = false;
501 : do
502 : {
503 21084 : converged = nekrs::runStep(corrector++);
504 21084 : } while (!converged);
505 :
506 : // TODO: time is somehow corrected here
507 21084 : nekrs::finishStep();
508 :
509 : // copy-pasta from Nek's main() for calling timers and printing
510 21084 : if (nekrs::updateFileCheckFreq())
511 21084 : if (_t_step % nekrs::updateFileCheckFreq())
512 20161 : nekrs::processUpdFile();
513 :
514 : // Note: here, we copy to both the nrs solution arrays and to the Nek5000 backend arrays,
515 : // because it is possible that users may interact using the legacy usr-file approach.
516 : // If we move away from the Nek5000 backend entirely, we could replace this line with
517 : // direct OCCA memcpy calls. But we do definitely need some type of copy here for _every_
518 : // time step, even if we're not technically passing data to another app, because we have
519 : // postprocessors that touch the `nrs` arrays that can be called in an arbitrary fashion
520 : // by the user.
521 21084 : auto nrs = nekrs::nrsPtr();
522 21084 : nrs->copyToNek(_timestepper->nondimensionalDT(step_end_time), _t_step);
523 :
524 21084 : if (nekrs::printStepInfoFreq())
525 21084 : if (_t_step % nekrs::printStepInfoFreq() == 0)
526 21084 : nekrs::printStepInfo(_timestepper->nondimensionalDT(_time), _t_step, false, true);
527 :
528 21084 : if (_is_output_step)
529 : {
530 1811 : writeFieldFile(step_end_time, _t_step);
531 :
532 : // TODO: I could not figure out why this can't be called from the destructor, to
533 : // add another field file on Cardinal's last time step. Revisit in the future.
534 1811 : if (_usrwrk_output)
535 : {
536 34 : static std::vector<bool> first_fld(_usrwrk_output->size(), true);
537 :
538 96 : for (unsigned int i = 0; i < _usrwrk_output->size(); ++i)
539 : {
540 62 : bool write_coords = first_fld[i] ? true : false;
541 :
542 124 : nekrs::write_usrwrk_field_file(_usrwrk_output->size(),
543 124 : i,
544 62 : (*_usrwrk_output)[i],
545 62 : (*_usrwrk_output_prefix)[i],
546 62 : _timestepper->nondimensionalDT(step_end_time),
547 62 : _t_step,
548 : write_coords);
549 :
550 : first_fld[i] = false;
551 : }
552 : }
553 : }
554 :
555 21084 : MPI_Barrier(comm().get());
556 21084 : const double elapsedStep = MPI_Wtime() - timeStartStep;
557 21084 : _tSolveStepMin = std::min(elapsedStep, _tSolveStepMin);
558 21084 : _tSolveStepMax = std::max(elapsedStep, _tSolveStepMax);
559 21084 : nekrs::updateTimer("minSolveStep", _tSolveStepMin);
560 21084 : nekrs::updateTimer("maxSolveStep", _tSolveStepMax);
561 :
562 21084 : _elapsedStepSum += elapsedStep;
563 21084 : _elapsedTime += elapsedStep;
564 21084 : nekrs::updateTimer("elapsedStep", elapsedStep);
565 21084 : nekrs::updateTimer("elapsedStepSum", _elapsedStepSum);
566 21084 : nekrs::updateTimer("elapsed", _elapsedTime);
567 :
568 21084 : if (nekrs::printStepInfoFreq())
569 21084 : if (_t_step % nekrs::printStepInfoFreq() == 0)
570 21084 : nekrs::printStepInfo(_timestepper->nondimensionalDT(_time), _t_step, true, false);
571 :
572 21084 : if (nekrs::runTimeStatFreq())
573 21084 : if (_t_step % nekrs::runTimeStatFreq() == 0)
574 12 : nekrs::printRuntimeStatistics(_t_step);
575 :
576 21084 : _time += _dt;
577 : }
578 :
579 : bool
580 42171 : NekRSProblem::isDataTransferHappening(ExternalProblem::Direction direction)
581 : {
582 42171 : if (nekrs::buildOnly())
583 : return false;
584 :
585 42171 : switch (direction)
586 : {
587 21087 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
588 21087 : return synchronizeIn();
589 21084 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
590 21084 : return synchronizeOut();
591 0 : default:
592 0 : mooseError("Unhandled DirectionEnum in NekRSProblem!");
593 : }
594 : }
595 :
596 : void
597 42171 : NekRSProblem::syncSolutions(ExternalProblem::Direction direction)
598 : {
599 42171 : auto & solution = _aux->solution();
600 :
601 42171 : if (!isDataTransferHappening(direction))
602 : return;
603 :
604 40682 : switch (direction)
605 : {
606 20342 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
607 : {
608 20342 : if (_first)
609 : {
610 684 : _serialized_solution->init(_aux->sys().n_dofs(), false, SERIAL);
611 684 : _first = false;
612 : }
613 :
614 20342 : solution.localize(*_serialized_solution);
615 :
616 : // execute all incoming field transfers
617 64090 : for (const auto & t : _field_transfers)
618 43750 : if (t->direction() == "to_nek")
619 13354 : t->sendDataToNek();
620 :
621 : // execute all incoming scalar transfers
622 20590 : for (const auto & t : _scalar_transfers)
623 250 : if (t->direction() == "to_nek")
624 250 : t->sendDataToNek();
625 :
626 : // update any user-defined properties (could be used for UQ)
627 20340 : auto nrs = nekrs::nrsPtr();
628 20340 : if (nrs->userProperties)
629 4541 : nrs->evaluateProperties(_timestepper->nondimensionalDT(_time));
630 :
631 : // copy host-side arrays which were filled to the device
632 20340 : copyHostToDevice();
633 :
634 20340 : break;
635 :
636 : return;
637 : }
638 20340 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
639 : {
640 : // fetch data from device to host
641 20340 : nekrs::copySolutionToHost();
642 :
643 : // execute all outgoing field transfers
644 64088 : for (const auto & t : _field_transfers)
645 43748 : if (t->direction() == "from_nek")
646 30396 : t->readDataFromNek();
647 :
648 : // execute all outgoing scalar transfers
649 : // TODO: is this the wrong trransfer?
650 20590 : for (const auto & t : _scalar_transfers)
651 250 : if (t->direction() == "from_nek")
652 0 : t->sendDataToNek();
653 :
654 : break;
655 : }
656 0 : default:
657 0 : mooseError("Unhandled Transfer::DIRECTION enum!");
658 : }
659 :
660 40680 : solution.close();
661 40680 : _aux->system().update();
662 : }
663 :
664 : bool
665 21087 : NekRSProblem::synchronizeIn()
666 : {
667 : bool synchronize = true;
668 :
669 21087 : switch (_sync_interval)
670 : {
671 681 : case synchronization::parent_app:
672 : {
673 : // For the minimized incoming synchronization to work correctly, the value
674 : // of the incoming postprocessor must not be zero. We only need to check this for the very
675 : // first time we evaluate this function. This ensures that you don't accidentally set a
676 : // zero value as a default in the master application's postprocessor.
677 681 : if (_first && *_transfer_in == false)
678 1 : mooseError("The default value for the 'transfer_in' postprocessor received by nekRS "
679 : "must not be false! Make sure that the master application's "
680 : "postprocessor is not zero.");
681 :
682 680 : if (*_transfer_in == false)
683 : synchronize = false;
684 : else
685 720 : setPostprocessorValueByName("transfer_in", false, 0);
686 :
687 : break;
688 : }
689 20406 : case synchronization::constant:
690 : {
691 20406 : synchronize = timeStep() % _constant_interval == 0;
692 20406 : break;
693 : }
694 0 : default:
695 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
696 : }
697 :
698 21086 : return synchronize;
699 : }
700 :
701 : bool
702 21084 : NekRSProblem::synchronizeOut()
703 : {
704 : bool synchronize = true;
705 :
706 21084 : switch (_sync_interval)
707 : {
708 680 : case synchronization::parent_app:
709 : {
710 680 : if (std::abs(_time - _dt - _transient_executioner->getTargetTime()) >
711 680 : _transient_executioner->timestepTol())
712 : synchronize = false;
713 : break;
714 : }
715 20404 : case synchronization::constant:
716 : {
717 20404 : synchronize = timeStep() % _constant_interval == 0;
718 20404 : break;
719 : }
720 0 : default:
721 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
722 : }
723 :
724 21084 : return synchronize;
725 : }
726 :
727 : bool
728 21084 : NekRSProblem::isOutputStep() const
729 : {
730 21084 : if (_app.isUltimateMaster())
731 : {
732 7670 : bool last_step = nekrs::lastStep(
733 7670 : _timestepper->nondimensionalDT(_time), _t_step, 0.0 /* dummy elapsed time */);
734 :
735 : // if Nek is controlled by a master application, then the last time step
736 : // is controlled by that master application, in which case we don't want to
737 : // write at what nekRS thinks is the last step (since it may or may not be
738 : // the actual end step), especially because we already ensure that we write on the
739 : // last time step from MOOSE's perspective in NekRSProblem's destructor.
740 7670 : if (last_step)
741 : return true;
742 : }
743 :
744 : // this routine does not check if we are on the last step - just whether we have
745 : // met the requested runtime or time step interval
746 :
747 20660 : return nekrs::checkpointStep(_timestepper->nondimensionalDT(_time), _t_step);
748 : }
749 :
750 : void
751 722 : NekRSProblem::addExternalVariables()
752 : {
753 : // Creation of variables for data transfers is handled by the FieldTransferBase objects
754 :
755 722 : if (_sync_interval == synchronization::parent_app)
756 : {
757 42 : auto pp_params = _factory.getValidParams("Receiver");
758 63 : pp_params.set<std::vector<OutputName>>("outputs") = {"none"};
759 :
760 : // we do not need to check for duplicate names because MOOSE already handles it
761 42 : addPostprocessor("Receiver", "transfer_in", pp_params);
762 21 : }
763 722 : }
764 :
765 : void
766 661428 : NekRSProblem::interpolateVolumeSolutionToNek(const int elem_id,
767 : double * incoming_moose_value,
768 : double * outgoing_nek_value)
769 : {
770 661428 : mesh_t * mesh = nekrs::entireMesh();
771 :
772 661428 : nekrs::interpolateVolumeHex3D(
773 661428 : _interpolation_incoming, incoming_moose_value, _moose_Nq, outgoing_nek_value, mesh->Nq);
774 661428 : }
775 :
776 : void
777 237046 : NekRSProblem::interpolateBoundarySolutionToNek(double * incoming_moose_value,
778 : double * outgoing_nek_value)
779 : {
780 237046 : mesh_t * mesh = nekrs::temperatureMesh();
781 :
782 237046 : double * scratch = (double *)calloc(_moose_Nq * mesh->Nq, sizeof(double));
783 :
784 237046 : nekrs::interpolateSurfaceFaceHex3D(scratch,
785 237046 : _interpolation_incoming,
786 : incoming_moose_value,
787 : _moose_Nq,
788 : outgoing_nek_value,
789 : mesh->Nq);
790 :
791 : freePointer(scratch);
792 237046 : }
793 :
794 : void
795 20340 : NekRSProblem::copyHostToDevice()
796 : {
797 33890 : for (const auto & slot : _usrwrk_slots)
798 13550 : copyIndividualScratchSlot(slot);
799 :
800 20340 : if (nekrs::hasMovingMesh())
801 0 : nekrs::copyDeformationToDevice();
802 20340 : }
803 :
804 : bool
805 27 : NekRSProblem::isUsrWrkSlotReservedForCoupling(const unsigned int & slot) const
806 : {
807 27 : return std::find(_usrwrk_slots.begin(), _usrwrk_slots.end(), slot) != _usrwrk_slots.end();
808 : }
809 :
810 : void
811 13553 : NekRSProblem::copyIndividualScratchSlot(const unsigned int & slot) const
812 : {
813 13553 : auto n = nekrs::fieldOffset();
814 13553 : auto nbytes = n * sizeof(dfloat);
815 :
816 13553 : auto nrs = nekrs::nrsPtr();
817 13553 : platform->app->bc->o_usrwrk.copyFrom(nekrs::host_wrk() + slot * n, n, slot * n);
818 13553 : }
819 :
820 : void
821 255590 : NekRSProblem::mapFaceDataToNekFace(const unsigned int & e,
822 : const unsigned int & var_num,
823 : const Real & divisor_scale,
824 : const Real & additive_scale,
825 : double ** outgoing_data)
826 : {
827 255590 : auto sys_number = _aux->number();
828 255590 : auto & mesh = _nek_mesh->getMesh();
829 255590 : auto indices = _nek_mesh->cornerIndices();
830 :
831 747340 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
832 : {
833 491750 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
834 :
835 : // Only work on elements we can find on our local chunk of a
836 : // distributed mesh
837 491750 : if (!elem_ptr)
838 : {
839 : libmesh_assert(!mesh.is_serial());
840 0 : continue;
841 : }
842 :
843 2458750 : for (unsigned int n = 0; n < _n_vertices_per_surface; n++)
844 : {
845 : auto node_ptr = elem_ptr->node_ptr(n);
846 :
847 : // convert libMesh node index into the ordering used by NekRS
848 1967000 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->boundaryNodeIndex(n)]
849 948184 : : _nek_mesh->boundaryNodeIndex(n);
850 :
851 1967000 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
852 1967000 : (*outgoing_data)[node_index] =
853 1967000 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
854 : }
855 : }
856 255590 : }
857 :
858 : void
859 515120 : NekRSProblem::mapFaceDataToNekVolume(const unsigned int & e,
860 : const unsigned int & var_num,
861 : const Real & divisor_scale,
862 : const Real & additive_scale,
863 : double ** outgoing_data)
864 : {
865 515120 : auto sys_number = _aux->number();
866 515120 : auto & mesh = _nek_mesh->getMesh();
867 515120 : auto indices = _nek_mesh->cornerIndices();
868 :
869 1854560 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
870 : {
871 1339440 : int n_faces_on_boundary = _nek_mesh->facesOnBoundary(e);
872 :
873 : // the only meaningful values are on the coupling boundaries, so we can skip this
874 : // interpolation if this volume element isn't on a coupling boundary
875 1339440 : if (n_faces_on_boundary > 0)
876 : {
877 148080 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
878 :
879 : // Only work on elements we can find on our local chunk of a
880 : // distributed mesh
881 148080 : if (!elem_ptr)
882 : {
883 : libmesh_assert(!mesh.is_serial());
884 0 : continue;
885 : }
886 :
887 1355520 : for (unsigned int n = 0; n < _n_vertices_per_volume; ++n)
888 : {
889 : auto node_ptr = elem_ptr->node_ptr(n);
890 :
891 : // convert libMesh node index into the ordering used by NekRS
892 1207440 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
893 879760 : : _nek_mesh->volumeNodeIndex(n);
894 :
895 1207440 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
896 1207440 : (*outgoing_data)[node_index] =
897 1207440 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
898 : }
899 : }
900 : }
901 515120 : }
902 :
903 : void
904 403188 : NekRSProblem::mapVolumeDataToNekVolume(const unsigned int & e,
905 : const unsigned int & var_num,
906 : const Real & divisor,
907 : const Real & additive,
908 : double ** outgoing_data)
909 : {
910 403188 : auto sys_number = _aux->number();
911 403188 : auto & mesh = _nek_mesh->getMesh();
912 403188 : auto indices = _nek_mesh->cornerIndices();
913 :
914 3404752 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
915 : {
916 3001564 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
917 :
918 : // Only work on elements we can find on our local chunk of a
919 : // distributed mesh
920 3001564 : if (!elem_ptr)
921 : {
922 : libmesh_assert(!mesh.is_serial());
923 0 : continue;
924 : }
925 28143246 : for (unsigned int n = 0; n < _n_vertices_per_volume; n++)
926 : {
927 : auto node_ptr = elem_ptr->node_ptr(n);
928 :
929 : // convert libMesh node index into the ordering used by NekRS
930 25141682 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
931 3241714 : : _nek_mesh->volumeNodeIndex(n);
932 :
933 25141682 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
934 25141682 : (*outgoing_data)[node_index] = ((*_serialized_solution)(dof_idx)-additive) / divisor;
935 : }
936 : }
937 403188 : }
938 :
939 : void
940 0 : NekRSProblem::writeVolumeDisplacement(const int elem_id,
941 : double * s,
942 : const field::NekWriteEnum f,
943 : const std::vector<double> * add)
944 : {
945 0 : mesh_t * mesh = nekrs::entireMesh();
946 0 : auto nrs = nekrs::nrsPtr();
947 :
948 0 : auto vc = _nek_mesh->volumeCoupling();
949 0 : int id = vc.element[elem_id] * mesh->Np;
950 :
951 0 : auto [x, y, z] = nekrs::host_xyz();
952 :
953 0 : if (_nek_mesh->exactMirror())
954 : {
955 : // can write directly into the NekRS solution
956 0 : for (int v = 0; v < mesh->Np; ++v)
957 : {
958 0 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
959 :
960 0 : if (f == field::x_displacement)
961 0 : x[id + v] = s[v] + extra;
962 0 : else if (f == field::y_displacement)
963 0 : y[id + v] = s[v] + extra;
964 0 : else if (f == field::z_displacement)
965 0 : z[id + v] = s[v] + extra;
966 : else
967 0 : mooseError("Unhandled NekWriteEnum in writeVolumeDisplacement!");
968 : }
969 : }
970 : else
971 : {
972 : // need to interpolate onto the higher-order Nek mesh
973 0 : double * tmp = (double *)calloc(mesh->Np, sizeof(double));
974 :
975 0 : interpolateVolumeSolutionToNek(elem_id, s, tmp);
976 :
977 0 : for (int v = 0; v < mesh->Np; ++v)
978 : {
979 0 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
980 0 : if (f == field::x_displacement)
981 0 : x[id + v] = s[v] + extra;
982 0 : else if (f == field::y_displacement)
983 0 : y[id + v] = s[v] + extra;
984 0 : else if (f == field::z_displacement)
985 0 : z[id + v] = s[v] + extra;
986 : else
987 0 : mooseError("Unhandled NekWriteEnum in writeVolumeDisplacement!");
988 : }
989 :
990 : freePointer(tmp);
991 : }
992 0 : }
993 :
994 : void
995 918308 : NekRSProblem::writeVolumeSolution(const int slot,
996 : const int elem_id,
997 : double * s,
998 : const std::vector<double> * add)
999 : {
1000 918308 : mesh_t * mesh = nekrs::entireMesh();
1001 918308 : auto nrs = nekrs::nrsPtr();
1002 :
1003 918308 : auto vc = _nek_mesh->volumeCoupling();
1004 918308 : int id = vc.element[elem_id] * mesh->Np;
1005 :
1006 918308 : auto usrwrk = nekrs::host_wrk();
1007 :
1008 918308 : if (_nek_mesh->exactMirror())
1009 : {
1010 : // can write directly into the NekRS solution
1011 9812432 : for (int v = 0; v < mesh->Np; ++v)
1012 : {
1013 9555552 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
1014 9555552 : usrwrk[slot + id + v] = s[v] + extra;
1015 : }
1016 : }
1017 : else
1018 : {
1019 : // need to interpolate onto the higher-order Nek mesh
1020 661428 : double * tmp = (double *)calloc(mesh->Np, sizeof(double));
1021 :
1022 661428 : interpolateVolumeSolutionToNek(elem_id, s, tmp);
1023 :
1024 125149652 : for (int v = 0; v < mesh->Np; ++v)
1025 : {
1026 124488224 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
1027 124488224 : usrwrk[slot + id + v] = tmp[v] + extra;
1028 : }
1029 :
1030 : freePointer(tmp);
1031 : }
1032 918308 : }
1033 :
1034 : void
1035 255590 : NekRSProblem::writeBoundarySolution(const int slot, const int elem_id, double * s)
1036 : {
1037 255590 : mesh_t * mesh = nekrs::temperatureMesh();
1038 255590 : auto nrs = nekrs::nrsPtr();
1039 :
1040 255590 : const auto & bc = _nek_mesh->boundaryCoupling();
1041 255590 : int offset = bc.element[elem_id] * mesh->Nfaces * mesh->Nfp + bc.face[elem_id] * mesh->Nfp;
1042 :
1043 255590 : auto usrwrk = nekrs::host_wrk();
1044 :
1045 255590 : if (_nek_mesh->exactMirror())
1046 : {
1047 : // can write directly into the NekRS solution
1048 428144 : for (int i = 0; i < mesh->Nfp; ++i)
1049 409600 : usrwrk[slot + mesh->vmapM[offset + i]] = s[i];
1050 : }
1051 : else
1052 : {
1053 : // need to interpolate onto the higher-order Nek mesh
1054 237046 : double * tmp = (double *)calloc(mesh->Nfp, sizeof(double));
1055 237046 : interpolateBoundarySolutionToNek(s, tmp);
1056 :
1057 5233822 : for (int i = 0; i < mesh->Nfp; ++i)
1058 4996776 : usrwrk[slot + mesh->vmapM[offset + i]] = tmp[i];
1059 :
1060 : freePointer(tmp);
1061 : }
1062 255590 : }
1063 : #endif
|
