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 2499 : NekRSProblem::validParams()
33 : {
34 2499 : InputParameters params = CardinalProblem::validParams();
35 4998 : 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 4998 : params.addParam<unsigned int>(
41 : "n_usrwrk_slots",
42 4998 : 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 4998 : 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 4998 : 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 4998 : params.addParam<bool>(
59 : "write_fld_files",
60 4998 : 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 4998 : params.addParam<bool>(
65 4998 : "disable_fld_file_output", false, "Whether to turn off all NekRS field file output writing");
66 :
67 4998 : params.addParam<bool>("skip_final_field_file",
68 4998 : false,
69 : "By default, we write a NekRS field file "
70 : "on the last time step; set this to true to disable");
71 :
72 4998 : params.addParam<MooseEnum>(
73 : "synchronization_interval",
74 4998 : 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 4998 : params.addParam<unsigned int>("constant_interval",
78 4998 : 1,
79 : "Constant interval (in units of number of time steps) with which "
80 : "to synchronize the NekRS solution");
81 2499 : return params;
82 0 : }
83 :
84 825 : NekRSProblem::NekRSProblem(const InputParameters & params)
85 : : CardinalProblem(params),
86 825 : _serialized_solution(NumericVector<Number>::build(_communicator).release()),
87 1650 : _casename(getParam<std::string>("casename")),
88 1650 : _write_fld_files(getParam<bool>("write_fld_files")),
89 1650 : _disable_fld_file_output(getParam<bool>("disable_fld_file_output")),
90 1650 : _n_usrwrk_slots(getParam<unsigned int>("n_usrwrk_slots")),
91 1650 : _constant_interval(getParam<unsigned int>("constant_interval")),
92 1650 : _skip_final_field_file(getParam<bool>("skip_final_field_file")),
93 825 : _start_time(nekrs::startTime()),
94 825 : _elapsedStepSum(0.0),
95 1650 : _elapsedTime(nekrs::getNekSetupTime()),
96 825 : _tSolveStepMin(std::numeric_limits<double>::max()),
97 1650 : _tSolveStepMax(std::numeric_limits<double>::min())
98 : {
99 825 : const auto & actions = getMooseApp().actionWarehouse().getActions<DimensionalizeAction>();
100 825 : _nondimensional = actions.size();
101 825 : nekrs::nondimensional(_nondimensional);
102 :
103 1650 : _sync_interval = getParam<MooseEnum>("synchronization_interval")
104 : .getEnum<synchronization::SynchronizationEnum>();
105 825 : 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 47 : 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 94 : checkUnusedParam(params, "constant_interval", "synchronizing based on the 'parent_app'");
121 : }
122 :
123 825 : 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 825 : 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 824 : _nek_mesh = dynamic_cast<NekRSMesh *>(&mesh());
132 :
133 824 : 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 823 : 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 823 : _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 823 : if (!nekrs::buildOnly())
153 823 : _nek_mesh->printMeshInfo();
154 :
155 : // boundary-specific data
156 823 : _n_surface_elems = _nek_mesh->numSurfaceElems();
157 823 : _n_vertices_per_surface = _nek_mesh->numVerticesPerSurface();
158 :
159 : // volume-specific data
160 823 : _n_vertices_per_volume = _nek_mesh->numVerticesPerVolume();
161 :
162 823 : if (_nek_mesh->volume())
163 486 : _n_points =
164 486 : _nek_mesh->numVolumeElems() * _n_vertices_per_volume * _nek_mesh->nBuildPerVolumeElem();
165 : else
166 337 : _n_points = _n_surface_elems * _n_vertices_per_surface * _nek_mesh->nBuildPerSurfaceElem();
167 :
168 823 : 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 823 : _needs_interpolation = _nek_mesh->numQuadraturePoints1D() > 2;
179 :
180 3292 : checkJointParams(
181 : params, {"usrwrk_output", "usrwrk_output_prefix"}, "outputting usrwrk slots to field files");
182 :
183 1646 : if (isParamValid("usrwrk_output"))
184 : {
185 52 : _usrwrk_output = &getParam<std::vector<unsigned int>>("usrwrk_output");
186 52 : _usrwrk_output_prefix = &getParam<std::vector<std::string>>("usrwrk_output_prefix");
187 :
188 64 : for (const auto & s : *_usrwrk_output)
189 39 : 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 25 : 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 2467 : }
199 :
200 1516 : 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 758 : if (!_is_output_step && !_skip_final_field_file)
205 : {
206 190 : 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 190 : writeFieldFile(_time, _t_step);
216 : }
217 :
218 758 : if (nekrs::runTimeStatFreq())
219 758 : if (_t_step % nekrs::runTimeStatFreq())
220 758 : nekrs::printRuntimeStatistics(_t_step);
221 :
222 758 : freePointer(_interpolation_outgoing);
223 758 : freePointer(_interpolation_incoming);
224 758 : nekrs::freeScratch();
225 :
226 758 : nekrs::finalize();
227 1516 : }
228 :
229 : void
230 2241 : NekRSProblem::writeFieldFile(const Real & step_end_time, const int & step) const
231 : {
232 2241 : if (_disable_fld_file_output)
233 : return;
234 :
235 2241 : Real t = _timestepper->nondimensionalDT(step_end_time);
236 :
237 2241 : 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 2185 : nekrs::outfld(t, step);
258 : }
259 :
260 : void
261 823 : NekRSProblem::initializeInterpolationMatrices()
262 : {
263 823 : mesh_t * mesh = nekrs::entireMesh();
264 :
265 : // determine the interpolation matrix for the outgoing transfer
266 823 : int starting_points = mesh->Nq;
267 823 : int ending_points = _nek_mesh->numQuadraturePoints1D();
268 823 : _interpolation_outgoing = (double *)calloc(starting_points * ending_points, sizeof(double));
269 823 : 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 823 : _interpolation_incoming = (double *)calloc(starting_points * ending_points, sizeof(double));
274 823 : nekrs::interpolationMatrix(_interpolation_incoming, starting_points, ending_points);
275 823 : }
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 768 : NekRSProblem::initialSetup()
290 : {
291 768 : CardinalProblem::initialSetup();
292 :
293 768 : auto executioner = _app.getExecutioner();
294 768 : _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 768 : 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 767 : _timestepper = dynamic_cast<NekTimeStepper *>(stepper);
307 767 : 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 766 : const auto moose_start_time = _transient_executioner->getStartTime();
315 766 : nekrs::setStartTime(_timestepper->nondimensionalDT(moose_start_time));
316 766 : _start_time = moose_start_time;
317 :
318 766 : if (_sync_interval == synchronization::parent_app)
319 47 : _transfer_in = &getPostprocessorValueByName("transfer_in");
320 :
321 : // Find all of the data transfer objects
322 766 : TheWarehouse::Query query = theWarehouse().query().condition<AttribSystem>("FieldTransfer");
323 766 : query.queryInto(_field_transfers);
324 :
325 : // Find all of the scalar data transfer objects
326 766 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("ScalarTransfer");
327 766 : uo_query.queryInto(_scalar_transfers);
328 :
329 : // We require a NekMeshDeformation object to exist if the NekRS model has a moving mesh
330 766 : 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 765 : 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 3825 : {"Slot", "MOOSE quantity", "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 1235 : for (const auto & field : field_usrwrk_map)
357 470 : _usrwrk_slots.insert(field.first);
358 814 : for (const auto & uo : _scalar_transfers)
359 49 : _usrwrk_slots.insert(uo->usrwrkSlot());
360 :
361 : // fill out table, being careful to only write information if owned by a field transfer,
362 : // a user object, or neither
363 1494 : for (int i = 0; i < _n_usrwrk_slots; ++i)
364 : {
365 729 : std::string oudf = "bc->usrwrk[" + std::to_string(i) + " * bc->fieldOffset+bc->idM]";
366 729 : std::string udf = "nrs->usrwrk[" + std::to_string(i) + " * nrs->fieldOffset + n]";
367 :
368 729 : if (field_usrwrk_map.find(i) != field_usrwrk_map.end())
369 : {
370 : // a field transfer owns it
371 940 : vt.addRow(i, field_usrwrk_map[i], oudf, udf);
372 : }
373 : else
374 : {
375 : // a user object might own it, or it could be unused
376 : bool owned_by_uo = false;
377 350 : for (const auto & uo : _scalar_transfers)
378 : {
379 91 : if (uo->usrwrkSlot() == i)
380 : {
381 : owned_by_uo = true;
382 49 : auto slot = std::to_string(uo->usrwrkSlot());
383 49 : auto count = std::to_string(uo->offset());
384 98 : vt.addRow(i,
385 49 : uo->name(),
386 147 : "bc->usrwrk[" + slot + " * bc->fieldOffset + " + count + "]",
387 147 : "nrs->usrwrk[" + slot + " * nrs->fieldOffset + " + count + "]");
388 : }
389 : }
390 :
391 259 : if (!owned_by_uo)
392 436 : vt.addRow(i, "unused", oudf, udf);
393 : }
394 : }
395 :
396 765 : if (_n_usrwrk_slots == 0)
397 : {
398 353 : _console << "Skipping allocation of NekRS scratch space because 'n_usrwrk_slots' is 0\n"
399 353 : << std::endl;
400 : }
401 : else
402 : {
403 : _console
404 412 : << "\n ===================> MAPPING FROM MOOSE TO NEKRS <===================\n"
405 412 : << std::endl;
406 412 : _console << " Slot: slice in scratch space holding the data" << std::endl;
407 412 : _console << " MOOSE quantity: name of the AuxVariable or Postprocessor that gets written into"
408 412 : << std::endl;
409 412 : _console << " this slot. If 'unused', this means that the space has been"
410 412 : << std::endl;
411 412 : _console << " allocated, but Cardinal is not otherwise using it for coupling"
412 412 : << std::endl;
413 412 : _console << " How to Access: C++ code to use in NekRS files; for the .udf instructions,"
414 412 : << std::endl;
415 412 : _console << " 'n' indicates a loop variable over GLL points\n" << std::endl;
416 412 : vt.print(_console);
417 412 : _console << std::endl;
418 : }
419 :
420 : // nekRS calls UDF_ExecuteStep once before the time stepping begins; the isLastStep stuff is
421 : // copy-pasta from NekRS main(), except that if Nek is a sub-app, we give full control of
422 : // time stepping to the main app
423 : bool isLastStep = false;
424 765 : if (_app.isUltimateMaster())
425 419 : isLastStep = !((nekrs::endTime() > nekrs::startTime() || nekrs::numSteps() > _t_step));
426 765 : nekrs::lastStep(isLastStep);
427 :
428 765 : nekrs::udfExecuteStep(
429 765 : _timestepper->nondimensionalDT(_start_time), _t_step, false /* not an output step */);
430 1530 : nekrs::resetTimer("udfExecuteStep");
431 2295 : }
432 :
433 : void
434 33618 : NekRSProblem::externalSolve()
435 : {
436 33618 : if (nekrs::buildOnly())
437 0 : return;
438 :
439 33618 : const double timeStartStep = MPI_Wtime();
440 :
441 : // _dt reflects the time step that MOOSE wants Nek to
442 : // take. For instance, if Nek is controlled by a master app and subcycling is used,
443 : // Nek must advance to the time interval taken by the master app. If the time step
444 : // that MOOSE wants nekRS to take (i.e. _dt) is smaller than we'd like nekRS to take, error.
445 33618 : if (_dt < _timestepper->minDT())
446 0 : mooseError("Requested time step of " + std::to_string(_dt) +
447 : " is smaller than the minimum "
448 0 : "time step of " +
449 0 : Moose::stringify(_timestepper->minDT()) +
450 : " allowed in NekRS!\n\n"
451 : "You can control this behavior with the 'min_dt' parameter on 'NekTimeStepper'.");
452 :
453 : // _time represents the time that we're simulating _to_, but we need to pass sometimes slightly
454 : // different times into the nekRS routines, which assume that the "time" passed into their
455 : // routines is sometimes a different interpretation.
456 33618 : double step_start_time = _time - _dt;
457 33618 : double step_end_time = _time;
458 :
459 33618 : _is_output_step = isOutputStep();
460 :
461 : // tell NekRS what the value of nrs->isOutputStep should be
462 33618 : nekrs::outputStep(_is_output_step);
463 :
464 : // NekRS prints out verbose info for the first 1000 time steps
465 33618 : if (_t_step <= 1000)
466 25242 : nekrs::verboseInfo(true);
467 :
468 : // Tell NekRS what the time step size is
469 33618 : nekrs::initStep(_timestepper->nondimensionalDT(step_start_time),
470 33618 : _timestepper->nondimensionalDT(_dt),
471 33618 : _t_step);
472 :
473 : // Run a nekRS time step. After the time step, this also calls UDF_ExecuteStep,
474 : // evaluated at (step_end_time, _t_step) == (nek_step_start_time + nek_dt, t_step)
475 : int corrector = 1;
476 : bool converged = false;
477 : do
478 : {
479 33618 : converged = nekrs::runStep(corrector++);
480 33618 : } while (!converged);
481 :
482 : // TODO: time is somehow corrected here
483 33618 : nekrs::finishStep();
484 :
485 : // copy-pasta from Nek's main() for calling timers and printing
486 33618 : if (nekrs::updateFileCheckFreq())
487 33618 : if (_t_step % nekrs::updateFileCheckFreq())
488 32095 : nekrs::processUpdFile();
489 :
490 : // Note: here, we copy to both the nrs solution arrays and to the Nek5000 backend arrays,
491 : // because it is possible that users may interact using the legacy usr-file approach.
492 : // If we move away from the Nek5000 backend entirely, we could replace this line with
493 : // direct OCCA memcpy calls. But we do definitely need some type of copy here for _every_
494 : // time step, even if we're not technically passing data to another app, because we have
495 : // postprocessors that touch the `nrs` arrays that can be called in an arbitrary fashion
496 : // by the user.
497 33618 : nek::ocopyToNek(_timestepper->nondimensionalDT(step_end_time), _t_step);
498 :
499 33618 : if (nekrs::printInfoFreq())
500 33618 : if (_t_step % nekrs::printInfoFreq() == 0)
501 33618 : nekrs::printInfo(_timestepper->nondimensionalDT(_time), _t_step, false, true);
502 :
503 33618 : if (_is_output_step)
504 : {
505 2051 : writeFieldFile(step_end_time, _t_step);
506 :
507 : // TODO: I could not figure out why this can't be called from the destructor, to
508 : // add another field file on Cardinal's last time step. Revisit in the future.
509 2051 : if (_usrwrk_output)
510 : {
511 32 : static std::vector<bool> first_fld(_usrwrk_output->size(), true);
512 :
513 84 : for (unsigned int i = 0; i < _usrwrk_output->size(); ++i)
514 : {
515 52 : bool write_coords = first_fld[i] ? true : false;
516 :
517 52 : nekrs::write_usrwrk_field_file((*_usrwrk_output)[i],
518 52 : (*_usrwrk_output_prefix)[i],
519 52 : _timestepper->nondimensionalDT(step_end_time),
520 52 : _t_step,
521 : write_coords);
522 :
523 : first_fld[i] = false;
524 : }
525 : }
526 : }
527 :
528 33618 : MPI_Barrier(comm().get());
529 33618 : const double elapsedStep = MPI_Wtime() - timeStartStep;
530 33618 : _tSolveStepMin = std::min(elapsedStep, _tSolveStepMin);
531 33618 : _tSolveStepMax = std::max(elapsedStep, _tSolveStepMax);
532 33618 : nekrs::updateTimer("minSolveStep", _tSolveStepMin);
533 33618 : nekrs::updateTimer("maxSolveStep", _tSolveStepMax);
534 :
535 33618 : _elapsedStepSum += elapsedStep;
536 33618 : _elapsedTime += elapsedStep;
537 33618 : nekrs::updateTimer("elapsedStep", elapsedStep);
538 33618 : nekrs::updateTimer("elapsedStepSum", _elapsedStepSum);
539 33618 : nekrs::updateTimer("elapsed", _elapsedTime);
540 :
541 33618 : if (nekrs::printInfoFreq())
542 33618 : if (_t_step % nekrs::printInfoFreq() == 0)
543 33618 : nekrs::printInfo(_timestepper->nondimensionalDT(_time), _t_step, true, false);
544 :
545 33618 : if (nekrs::runTimeStatFreq())
546 33618 : if (_t_step % nekrs::runTimeStatFreq() == 0)
547 32 : nekrs::printRuntimeStatistics(_t_step);
548 :
549 33618 : _time += _dt;
550 : }
551 :
552 : bool
553 67239 : NekRSProblem::isDataTransferHappening(ExternalProblem::Direction direction)
554 : {
555 67239 : if (nekrs::buildOnly())
556 : return false;
557 :
558 67239 : switch (direction)
559 : {
560 33621 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
561 33621 : return synchronizeIn();
562 33618 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
563 33618 : return synchronizeOut();
564 0 : default:
565 0 : mooseError("Unhandled DirectionEnum in NekRSProblem!");
566 : }
567 : }
568 :
569 : void
570 67239 : NekRSProblem::syncSolutions(ExternalProblem::Direction direction)
571 : {
572 67239 : auto & solution = _aux->solution();
573 :
574 67239 : if (!isDataTransferHappening(direction))
575 : return;
576 :
577 65266 : switch (direction)
578 : {
579 32634 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
580 : {
581 32634 : if (_first)
582 : {
583 764 : _serialized_solution->init(_aux->sys().n_dofs(), false, SERIAL);
584 764 : _first = false;
585 : }
586 :
587 32634 : solution.localize(*_serialized_solution);
588 :
589 : // execute all incoming field transfers
590 136930 : for (const auto & t : _field_transfers)
591 104298 : if (t->direction() == "to_nek")
592 15606 : t->sendDataToNek();
593 :
594 : // execute all incoming scalar transfers
595 33282 : for (const auto & t : _scalar_transfers)
596 650 : if (t->direction() == "to_nek")
597 650 : t->sendDataToNek();
598 :
599 32632 : if (udf.properties)
600 : {
601 16141 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
602 16141 : evaluateProperties(nrs, _timestepper->nondimensionalDT(_time));
603 : }
604 :
605 32632 : copyScratchToDevice();
606 :
607 32632 : break;
608 :
609 : return;
610 : }
611 32632 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
612 : {
613 : // execute all outgoing field transfers
614 136928 : for (const auto & t : _field_transfers)
615 104296 : if (t->direction() == "from_nek")
616 88692 : t->readDataFromNek();
617 :
618 : // execute all outgoing scalar transfers
619 33282 : for (const auto & t : _scalar_transfers)
620 650 : if (t->direction() == "from_nek")
621 0 : t->sendDataToNek();
622 :
623 : break;
624 : }
625 0 : default:
626 0 : mooseError("Unhandled Transfer::DIRECTION enum!");
627 : }
628 :
629 65264 : solution.close();
630 65264 : _aux->system().update();
631 : }
632 :
633 : bool
634 33621 : NekRSProblem::synchronizeIn()
635 : {
636 : bool synchronize = true;
637 :
638 33621 : switch (_sync_interval)
639 : {
640 977 : case synchronization::parent_app:
641 : {
642 : // For the minimized incoming synchronization to work correctly, the value
643 : // of the incoming postprocessor must not be zero. We only need to check this for the very
644 : // first time we evaluate this function. This ensures that you don't accidentally set a
645 : // zero value as a default in the master application's postprocessor.
646 977 : if (_first && *_transfer_in == false)
647 1 : mooseError("The default value for the 'transfer_in' postprocessor received by nekRS "
648 : "must not be false! Make sure that the master application's "
649 : "postprocessor is not zero.");
650 :
651 976 : if (*_transfer_in == false)
652 : synchronize = false;
653 : else
654 828 : setPostprocessorValueByName("transfer_in", false, 0);
655 :
656 : break;
657 : }
658 32644 : case synchronization::constant:
659 : {
660 32644 : synchronize = timeStep() % _constant_interval == 0;
661 32644 : break;
662 : }
663 0 : default:
664 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
665 : }
666 :
667 33620 : return synchronize;
668 : }
669 :
670 : bool
671 33618 : NekRSProblem::synchronizeOut()
672 : {
673 : bool synchronize = true;
674 :
675 33618 : switch (_sync_interval)
676 : {
677 976 : case synchronization::parent_app:
678 : {
679 976 : if (std::abs(_time - _dt - _transient_executioner->getTargetTime()) >
680 976 : _transient_executioner->timestepTol())
681 : synchronize = false;
682 : break;
683 : }
684 32642 : case synchronization::constant:
685 : {
686 32642 : synchronize = timeStep() % _constant_interval == 0;
687 32642 : break;
688 : }
689 0 : default:
690 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
691 : }
692 :
693 33618 : return synchronize;
694 : }
695 :
696 : bool
697 33618 : NekRSProblem::isOutputStep() const
698 : {
699 33618 : if (_app.isUltimateMaster())
700 : {
701 17734 : bool last_step = nekrs::lastStep(
702 17734 : _timestepper->nondimensionalDT(_time), _t_step, 0.0 /* dummy elapsed time */);
703 :
704 : // if Nek is controlled by a master application, then the last time step
705 : // is controlled by that master application, in which case we don't want to
706 : // write at what nekRS thinks is the last step (since it may or may not be
707 : // the actual end step), especially because we already ensure that we write on the
708 : // last time step from MOOSE's perspective in NekRSProblem's destructor.
709 17734 : if (last_step)
710 : return true;
711 : }
712 :
713 : // this routine does not check if we are on the last step - just whether we have
714 : // met the requested runtime or time step interval
715 33176 : return nekrs::outputStep(_timestepper->nondimensionalDT(_time), _t_step);
716 : }
717 :
718 : void
719 798 : NekRSProblem::addExternalVariables()
720 : {
721 : // Creation of variables for data transfers is handled by the FieldTransferBase objects
722 :
723 798 : if (_sync_interval == synchronization::parent_app)
724 : {
725 94 : auto pp_params = _factory.getValidParams("Receiver");
726 141 : pp_params.set<std::vector<OutputName>>("outputs") = {"none"};
727 :
728 : // we do not need to check for duplicate names because MOOSE already handles it
729 94 : addPostprocessor("Receiver", "transfer_in", pp_params);
730 47 : }
731 798 : }
732 :
733 : void
734 1315686 : NekRSProblem::interpolateVolumeSolutionToNek(const int elem_id,
735 : double * incoming_moose_value,
736 : double * outgoing_nek_value)
737 : {
738 1315686 : mesh_t * mesh = nekrs::entireMesh();
739 :
740 1315686 : nekrs::interpolateVolumeHex3D(
741 1315686 : _interpolation_incoming, incoming_moose_value, _moose_Nq, outgoing_nek_value, mesh->Nq);
742 1315686 : }
743 :
744 : void
745 399658 : NekRSProblem::interpolateBoundarySolutionToNek(double * incoming_moose_value,
746 : double * outgoing_nek_value)
747 : {
748 399658 : mesh_t * mesh = nekrs::temperatureMesh();
749 :
750 399658 : double * scratch = (double *)calloc(_moose_Nq * mesh->Nq, sizeof(double));
751 :
752 399658 : nekrs::interpolateSurfaceFaceHex3D(scratch,
753 399658 : _interpolation_incoming,
754 : incoming_moose_value,
755 : _moose_Nq,
756 : outgoing_nek_value,
757 : mesh->Nq);
758 :
759 : freePointer(scratch);
760 399658 : }
761 :
762 : void
763 32632 : NekRSProblem::copyScratchToDevice()
764 : {
765 50642 : for (const auto & slot : _usrwrk_slots)
766 18010 : copyIndividualScratchSlot(slot);
767 :
768 32632 : if (nekrs::hasMovingMesh())
769 1504 : nekrs::copyDeformationToDevice();
770 32632 : }
771 :
772 : bool
773 35 : NekRSProblem::isUsrWrkSlotReservedForCoupling(const unsigned int & slot) const
774 : {
775 35 : return std::find(_usrwrk_slots.begin(), _usrwrk_slots.end(), slot) != _usrwrk_slots.end();
776 : }
777 :
778 : void
779 18037 : NekRSProblem::copyIndividualScratchSlot(const unsigned int & slot) const
780 : {
781 18037 : auto n = nekrs::fieldOffset();
782 18037 : auto nbytes = n * sizeof(dfloat);
783 :
784 18037 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
785 18037 : nrs->o_usrwrk.copyFrom(nrs->usrwrk + slot * n, nbytes, slot * nbytes);
786 18037 : }
787 :
788 : void
789 424202 : NekRSProblem::mapFaceDataToNekFace(const unsigned int & e,
790 : const unsigned int & var_num,
791 : const Real & divisor_scale,
792 : const Real & additive_scale,
793 : double ** outgoing_data)
794 : {
795 424202 : auto sys_number = _aux->number();
796 424202 : auto & mesh = _nek_mesh->getMesh();
797 424202 : auto indices = _nek_mesh->cornerIndices();
798 :
799 1216564 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
800 : {
801 792362 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
802 :
803 : // Only work on elements we can find on our local chunk of a
804 : // distributed mesh
805 792362 : if (!elem_ptr)
806 : {
807 : libmesh_assert(!mesh.is_serial());
808 0 : continue;
809 : }
810 :
811 4633810 : for (unsigned int n = 0; n < _n_vertices_per_surface; n++)
812 : {
813 : auto node_ptr = elem_ptr->node_ptr(n);
814 :
815 : // convert libMesh node index into the ordering used by NekRS
816 3841448 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->boundaryNodeIndex(n)]
817 2270632 : : _nek_mesh->boundaryNodeIndex(n);
818 :
819 3841448 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
820 3841448 : (*outgoing_data)[node_index] =
821 3841448 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
822 : }
823 : }
824 424202 : }
825 :
826 : void
827 1039536 : NekRSProblem::mapFaceDataToNekVolume(const unsigned int & e,
828 : const unsigned int & var_num,
829 : const Real & divisor_scale,
830 : const Real & additive_scale,
831 : double ** outgoing_data)
832 : {
833 1039536 : auto sys_number = _aux->number();
834 1039536 : auto & mesh = _nek_mesh->getMesh();
835 1039536 : auto indices = _nek_mesh->cornerIndices();
836 :
837 2491232 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
838 : {
839 1451696 : int n_faces_on_boundary = _nek_mesh->facesOnBoundary(e);
840 :
841 : // the only meaningful values are on the coupling boundaries, so we can skip this
842 : // interpolation if this volume element isn't on a coupling boundary
843 1451696 : if (n_faces_on_boundary > 0)
844 : {
845 288568 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
846 :
847 : // Only work on elements we can find on our local chunk of a
848 : // distributed mesh
849 288568 : if (!elem_ptr)
850 : {
851 : libmesh_assert(!mesh.is_serial());
852 0 : continue;
853 : }
854 :
855 5606712 : for (unsigned int n = 0; n < _n_vertices_per_volume; ++n)
856 : {
857 : auto node_ptr = elem_ptr->node_ptr(n);
858 :
859 : // convert libMesh node index into the ordering used by NekRS
860 5318144 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
861 5154304 : : _nek_mesh->volumeNodeIndex(n);
862 :
863 5318144 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
864 5318144 : (*outgoing_data)[node_index] =
865 5318144 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
866 : }
867 : }
868 : }
869 1039536 : }
870 :
871 : void
872 401758 : NekRSProblem::mapVolumeDataToNekVolume(const unsigned int & e,
873 : const unsigned int & var_num,
874 : const Real & divisor,
875 : const Real & additive,
876 : double ** outgoing_data)
877 : {
878 401758 : auto sys_number = _aux->number();
879 401758 : auto & mesh = _nek_mesh->getMesh();
880 401758 : auto indices = _nek_mesh->cornerIndices();
881 :
882 1419724 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
883 : {
884 1017966 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
885 :
886 : // Only work on elements we can find on our local chunk of a
887 : // distributed mesh
888 1017966 : if (!elem_ptr)
889 : {
890 : libmesh_assert(!mesh.is_serial());
891 0 : continue;
892 : }
893 10642858 : for (unsigned int n = 0; n < _n_vertices_per_volume; n++)
894 : {
895 : auto node_ptr = elem_ptr->node_ptr(n);
896 :
897 : // convert libMesh node index into the ordering used by NekRS
898 9624892 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
899 4161404 : : _nek_mesh->volumeNodeIndex(n);
900 :
901 9624892 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
902 9624892 : (*outgoing_data)[node_index] = ((*_serialized_solution)(dof_idx)-additive) / divisor;
903 : }
904 : }
905 401758 : }
906 :
907 : #endif
|
