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