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 2421 : NekRSProblem::validParams()
33 : {
34 2421 : InputParameters params = CardinalProblem::validParams();
35 4842 : 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 4842 : params.addParam<unsigned int>(
41 : "n_usrwrk_slots",
42 4842 : 7,
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 4842 : 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 4842 : 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 4842 : params.addParam<bool>(
59 : "write_fld_files",
60 4842 : 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 4842 : params.addParam<bool>(
65 4842 : "disable_fld_file_output", false, "Whether to turn off all NekRS field file output writing");
66 :
67 4842 : params.addParam<bool>("skip_final_field_file",
68 4842 : false,
69 : "By default, we write a NekRS field file "
70 : "on the last time step; set this to true to disable");
71 :
72 4842 : params.addParam<MooseEnum>(
73 : "synchronization_interval",
74 4842 : 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 4842 : params.addParam<unsigned int>("constant_interval",
78 4842 : 1,
79 : "Constant interval (in units of number of time steps) with which "
80 : "to synchronize the NekRS solution");
81 2421 : return params;
82 0 : }
83 :
84 799 : NekRSProblem::NekRSProblem(const InputParameters & params)
85 : : CardinalProblem(params),
86 799 : _serialized_solution(NumericVector<Number>::build(_communicator).release()),
87 1598 : _casename(getParam<std::string>("casename")),
88 1598 : _write_fld_files(getParam<bool>("write_fld_files")),
89 1598 : _disable_fld_file_output(getParam<bool>("disable_fld_file_output")),
90 1598 : _n_usrwrk_slots(getParam<unsigned int>("n_usrwrk_slots")),
91 1598 : _constant_interval(getParam<unsigned int>("constant_interval")),
92 1598 : _skip_final_field_file(getParam<bool>("skip_final_field_file")),
93 799 : _start_time(nekrs::startTime()),
94 799 : _elapsedStepSum(0.0),
95 1598 : _elapsedTime(nekrs::getNekSetupTime()),
96 799 : _tSolveStepMin(std::numeric_limits<double>::max()),
97 1598 : _tSolveStepMax(std::numeric_limits<double>::min())
98 : {
99 799 : const auto & actions = getMooseApp().actionWarehouse().getActions<DimensionalizeAction>();
100 799 : _nondimensional = actions.size();
101 799 : nekrs::nondimensional(_nondimensional);
102 :
103 1598 : _sync_interval = getParam<MooseEnum>("synchronization_interval")
104 : .getEnum<synchronization::SynchronizationEnum>();
105 799 : 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 799 : 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 799 : 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 798 : _nek_mesh = dynamic_cast<NekRSMesh *>(&mesh());
132 :
133 798 : 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 797 : 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 797 : _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 797 : if (!nekrs::buildOnly())
153 797 : _nek_mesh->printMeshInfo();
154 :
155 : // boundary-specific data
156 797 : _n_surface_elems = _nek_mesh->numSurfaceElems();
157 797 : _n_vertices_per_surface = _nek_mesh->numVerticesPerSurface();
158 :
159 : // volume-specific data
160 797 : _n_vertices_per_volume = _nek_mesh->numVerticesPerVolume();
161 :
162 797 : if (_nek_mesh->volume())
163 481 : _n_points =
164 481 : _nek_mesh->numVolumeElems() * _n_vertices_per_volume * _nek_mesh->nBuildPerVolumeElem();
165 : else
166 316 : _n_points = _n_surface_elems * _n_vertices_per_surface * _nek_mesh->nBuildPerSurfaceElem();
167 :
168 797 : 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 797 : _needs_interpolation = _nek_mesh->numQuadraturePoints1D() > 2;
179 :
180 3188 : checkJointParams(
181 : params, {"usrwrk_output", "usrwrk_output_prefix"}, "outputting usrwrk slots to field files");
182 :
183 1594 : 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 : "!");
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 2389 : }
199 :
200 1468 : 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 734 : 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 734 : if (nekrs::runTimeStatFreq())
219 734 : if (_t_step % nekrs::runTimeStatFreq())
220 734 : nekrs::printRuntimeStatistics(_t_step);
221 :
222 734 : freePointer(_interpolation_outgoing);
223 734 : freePointer(_interpolation_incoming);
224 734 : nekrs::freeScratch();
225 :
226 734 : nekrs::finalize();
227 1468 : }
228 :
229 : void
230 2197 : NekRSProblem::writeFieldFile(const Real & step_end_time, const int & step) const
231 : {
232 2197 : if (_disable_fld_file_output)
233 : return;
234 :
235 2197 : Real t = _timestepper->nondimensionalDT(step_end_time);
236 :
237 2197 : 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 2141 : nekrs::outfld(t, step);
258 : }
259 :
260 : void
261 797 : NekRSProblem::initializeInterpolationMatrices()
262 : {
263 797 : mesh_t * mesh = nekrs::entireMesh();
264 :
265 : // determine the interpolation matrix for the outgoing transfer
266 797 : int starting_points = mesh->Nq;
267 797 : int ending_points = _nek_mesh->numQuadraturePoints1D();
268 797 : _interpolation_outgoing = (double *)calloc(starting_points * ending_points, sizeof(double));
269 797 : 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 797 : _interpolation_incoming = (double *)calloc(starting_points * ending_points, sizeof(double));
274 797 : nekrs::interpolationMatrix(_interpolation_incoming, starting_points, ending_points);
275 797 : }
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 744 : NekRSProblem::initialSetup()
290 : {
291 744 : CardinalProblem::initialSetup();
292 :
293 744 : auto executioner = _app.getExecutioner();
294 744 : _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 744 : 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 743 : _timestepper = dynamic_cast<NekTimeStepper *>(stepper);
307 743 : 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 742 : const auto moose_start_time = _transient_executioner->getStartTime();
315 742 : nekrs::setStartTime(_timestepper->nondimensionalDT(moose_start_time));
316 742 : _start_time = moose_start_time;
317 :
318 742 : if (_sync_interval == synchronization::parent_app)
319 47 : _transfer_in = &getPostprocessorValueByName("transfer_in");
320 :
321 : // Find all of the data transfer objects
322 742 : TheWarehouse::Query query = theWarehouse().query().condition<AttribSystem>("FieldTransfer");
323 742 : query.queryInto(_field_transfers);
324 :
325 : // Find all of the scalar data transfer objects
326 742 : TheWarehouse::Query uo_query = theWarehouse().query().condition<AttribSystem>("ScalarTransfer");
327 742 : uo_query.queryInto(_scalar_transfers);
328 :
329 : // We require a NekMeshDeformation object to exist if the NekRS model has a moving mesh
330 742 : 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 741 : 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 3705 : {"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 1191 : for (const auto & field : field_usrwrk_map)
357 450 : _usrwrk_slots.insert(field.first);
358 790 : 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 5443 : for (int i = 0; i < _n_usrwrk_slots; ++i)
364 : {
365 4702 : std::string oudf = "bc->usrwrk[" + std::to_string(i) + " * bc->fieldOffset+bc->idM]";
366 4702 : std::string udf = "nrs->usrwrk[" + std::to_string(i) + " * nrs->fieldOffset + n]";
367 :
368 4702 : if (field_usrwrk_map.find(i) != field_usrwrk_map.end())
369 : {
370 : // a field transfer owns it
371 900 : 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 4461 : for (const auto & uo : _scalar_transfers)
378 : {
379 209 : 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 4252 : if (!owned_by_uo)
392 8422 : vt.addRow(i, "unused", oudf, udf);
393 : }
394 : }
395 :
396 741 : if (_n_usrwrk_slots == 0)
397 : {
398 20 : _console << "Skipping allocation of NekRS scratch space because 'n_usrwrk_slots' is 0\n"
399 20 : << std::endl;
400 : }
401 : else
402 : {
403 : _console
404 721 : << "\n ===================> MAPPING FROM MOOSE TO NEKRS <===================\n"
405 721 : << std::endl;
406 721 : _console << " Slot: slice in scratch space holding the data" << std::endl;
407 721 : _console << " MOOSE quantity: name of the AuxVariable or Postprocessor that gets written into"
408 721 : << std::endl;
409 721 : _console << " this slot. If 'unused', this means that the space has been"
410 721 : << std::endl;
411 721 : _console << " allocated, but Cardinal is not otherwise using it for coupling"
412 721 : << std::endl;
413 721 : _console << " How to Access: C++ code to use in NekRS files; for the .udf instructions,"
414 721 : << std::endl;
415 721 : _console << " 'n' indicates a loop variable over GLL points\n" << std::endl;
416 721 : vt.print(_console);
417 721 : _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 741 : if (_app.isUltimateMaster())
425 407 : isLastStep = !((nekrs::endTime() > nekrs::startTime() || nekrs::numSteps() > _t_step));
426 741 : nekrs::lastStep(isLastStep);
427 :
428 741 : nekrs::udfExecuteStep(
429 741 : _timestepper->nondimensionalDT(_start_time), _t_step, false /* not an output step */);
430 1482 : nekrs::resetTimer("udfExecuteStep");
431 2223 : }
432 :
433 : void
434 33358 : NekRSProblem::externalSolve()
435 : {
436 33358 : if (nekrs::buildOnly())
437 0 : return;
438 :
439 33358 : 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 33358 : 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 33358 : double step_start_time = _time - _dt;
457 33358 : double step_end_time = _time;
458 :
459 33358 : _is_output_step = isOutputStep();
460 :
461 : // tell NekRS what the value of nrs->isOutputStep should be
462 33358 : nekrs::outputStep(_is_output_step);
463 :
464 : // NekRS prints out verbose info for the first 1000 time steps
465 33358 : if (_t_step <= 1000)
466 24982 : nekrs::verboseInfo(true);
467 :
468 : // Tell NekRS what the time step size is
469 33358 : nekrs::initStep(_timestepper->nondimensionalDT(step_start_time),
470 33358 : _timestepper->nondimensionalDT(_dt),
471 33358 : _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 33358 : converged = nekrs::runStep(corrector++);
480 33358 : } while (!converged);
481 :
482 : // TODO: time is somehow corrected here
483 33358 : nekrs::finishStep();
484 :
485 : // copy-pasta from Nek's main() for calling timers and printing
486 33358 : if (nekrs::updateFileCheckFreq())
487 33358 : if (_t_step % nekrs::updateFileCheckFreq())
488 31847 : 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 33358 : nek::ocopyToNek(_timestepper->nondimensionalDT(step_end_time), _t_step);
498 :
499 33358 : if (nekrs::printInfoFreq())
500 33358 : if (_t_step % nekrs::printInfoFreq() == 0)
501 33358 : nekrs::printInfo(_timestepper->nondimensionalDT(_time), _t_step, false, true);
502 :
503 33358 : if (_is_output_step)
504 : {
505 2007 : 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 2007 : 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 33358 : MPI_Barrier(comm().get());
529 33358 : const double elapsedStep = MPI_Wtime() - timeStartStep;
530 33358 : _tSolveStepMin = std::min(elapsedStep, _tSolveStepMin);
531 33358 : _tSolveStepMax = std::max(elapsedStep, _tSolveStepMax);
532 33358 : nekrs::updateTimer("minSolveStep", _tSolveStepMin);
533 33358 : nekrs::updateTimer("maxSolveStep", _tSolveStepMax);
534 :
535 33358 : _elapsedStepSum += elapsedStep;
536 33358 : _elapsedTime += elapsedStep;
537 33358 : nekrs::updateTimer("elapsedStep", elapsedStep);
538 33358 : nekrs::updateTimer("elapsedStepSum", _elapsedStepSum);
539 33358 : nekrs::updateTimer("elapsed", _elapsedTime);
540 :
541 33358 : if (nekrs::printInfoFreq())
542 33358 : if (_t_step % nekrs::printInfoFreq() == 0)
543 33358 : nekrs::printInfo(_timestepper->nondimensionalDT(_time), _t_step, true, false);
544 :
545 33358 : if (nekrs::runTimeStatFreq())
546 33358 : if (_t_step % nekrs::runTimeStatFreq() == 0)
547 32 : nekrs::printRuntimeStatistics(_t_step);
548 :
549 33358 : _time += _dt;
550 : }
551 :
552 : bool
553 66719 : NekRSProblem::isDataTransferHappening(ExternalProblem::Direction direction)
554 : {
555 66719 : if (nekrs::buildOnly())
556 : return false;
557 :
558 66719 : switch (direction)
559 : {
560 33361 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
561 33361 : return synchronizeIn();
562 33358 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
563 33358 : return synchronizeOut();
564 0 : default:
565 0 : mooseError("Unhandled DirectionEnum in NekRSProblem!");
566 : }
567 : }
568 :
569 : void
570 66719 : NekRSProblem::syncSolutions(ExternalProblem::Direction direction)
571 : {
572 66719 : auto & solution = _aux->solution();
573 :
574 66719 : if (!isDataTransferHappening(direction))
575 : return;
576 :
577 64746 : switch (direction)
578 : {
579 32374 : case ExternalProblem::Direction::TO_EXTERNAL_APP:
580 : {
581 32374 : if (_first)
582 : {
583 740 : _serialized_solution->init(_aux->sys().n_dofs(), false, SERIAL);
584 740 : _first = false;
585 : }
586 :
587 32374 : solution.localize(*_serialized_solution);
588 :
589 : // execute all incoming field transfers
590 136154 : for (const auto & t : _field_transfers)
591 103782 : if (t->direction() == "to_nek")
592 15350 : t->sendDataToNek();
593 :
594 : // execute all incoming scalar transfers
595 33022 : for (const auto & t : _scalar_transfers)
596 650 : if (t->direction() == "to_nek")
597 650 : t->sendDataToNek();
598 :
599 32372 : if (udf.properties)
600 : {
601 16141 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
602 16141 : evaluateProperties(nrs, _timestepper->nondimensionalDT(_time));
603 : }
604 :
605 32372 : copyScratchToDevice();
606 :
607 32372 : break;
608 :
609 : return;
610 : }
611 32372 : case ExternalProblem::Direction::FROM_EXTERNAL_APP:
612 : {
613 : // execute all outgoing field transfers
614 136152 : for (const auto & t : _field_transfers)
615 103780 : if (t->direction() == "from_nek")
616 88432 : t->readDataFromNek();
617 :
618 : // execute all outgoing scalar transfers
619 33022 : 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 64744 : solution.close();
630 64744 : _aux->system().update();
631 : }
632 :
633 : bool
634 33361 : NekRSProblem::synchronizeIn()
635 : {
636 : bool synchronize = true;
637 :
638 33361 : 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 32384 : case synchronization::constant:
659 : {
660 32384 : synchronize = timeStep() % _constant_interval == 0;
661 32384 : break;
662 : }
663 0 : default:
664 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
665 : }
666 :
667 33360 : return synchronize;
668 : }
669 :
670 : bool
671 33358 : NekRSProblem::synchronizeOut()
672 : {
673 : bool synchronize = true;
674 :
675 33358 : 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 32382 : case synchronization::constant:
685 : {
686 32382 : synchronize = timeStep() % _constant_interval == 0;
687 32382 : break;
688 : }
689 0 : default:
690 0 : mooseError("Unhandled SynchronizationEnum in NekRSProblem!");
691 : }
692 :
693 33358 : return synchronize;
694 : }
695 :
696 : bool
697 33358 : NekRSProblem::isOutputStep() const
698 : {
699 33358 : if (_app.isUltimateMaster())
700 : {
701 17714 : bool last_step = nekrs::lastStep(
702 17714 : _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 17714 : 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 32928 : return nekrs::outputStep(_timestepper->nondimensionalDT(_time), _t_step);
716 : }
717 :
718 : void
719 774 : NekRSProblem::addExternalVariables()
720 : {
721 : // Creation of variables for data transfers is handled by the FieldTransferBase objects
722 :
723 774 : 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 774 : }
732 :
733 : void
734 68022 : NekRSProblem::volumeSolution(const field::NekFieldEnum & field, double * T)
735 : {
736 68022 : mesh_t * mesh = nekrs::entireMesh();
737 68022 : auto vc = _nek_mesh->volumeCoupling();
738 :
739 : double (*f)(int);
740 68022 : f = nekrs::solutionPointer(field);
741 :
742 68022 : int start_1d = mesh->Nq;
743 68022 : int end_1d = _moose_Nq;
744 68022 : int start_3d = start_1d * start_1d * start_1d;
745 68022 : int end_3d = end_1d * end_1d * end_1d;
746 :
747 68022 : int n_to_write = vc.n_elems * end_3d * _nek_mesh->nBuildPerVolumeElem();
748 :
749 : // allocate temporary space:
750 : // - Ttmp: results of the search for each process
751 : // - Telem: scratch space for volume interpolation to avoid reallocating a bunch (only used if
752 : // interpolating)
753 68022 : double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
754 68022 : double * Telem = (double *)calloc(start_3d, sizeof(double));
755 :
756 : auto indices = _nek_mesh->cornerIndices();
757 :
758 : int c = 0;
759 3365896 : for (int k = 0; k < mesh->Nelements; ++k)
760 : {
761 3297874 : int offset = k * start_3d;
762 :
763 7212788 : for (int build = 0; build < _nek_mesh->nBuildPerVolumeElem(); ++build)
764 : {
765 3914914 : if (_needs_interpolation)
766 : {
767 : // get the solution on the element
768 574516016 : for (int v = 0; v < start_3d; ++v)
769 572841440 : Telem[v] = f(offset + v);
770 :
771 : // and then interpolate it
772 1674576 : nekrs::interpolateVolumeHex3D(_interpolation_outgoing, Telem, start_1d, &(Ttmp[c]), end_1d);
773 1674576 : c += end_3d;
774 : }
775 : else
776 : {
777 : // get the solution on the element - no need to interpolate
778 20163042 : for (int v = 0; v < end_3d; ++v, ++c)
779 17922704 : Ttmp[c] = f(offset + indices[build][v]);
780 : }
781 : }
782 : }
783 :
784 : // dimensionalize the solution if needed
785 63204278 : for (int v = 0; v < n_to_write; ++v)
786 63136256 : Ttmp[v] = Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);
787 :
788 68022 : nekrs::allgatherv(vc.mirror_counts, Ttmp, T, end_3d);
789 :
790 : freePointer(Ttmp);
791 : freePointer(Telem);
792 68022 : }
793 :
794 : void
795 20410 : NekRSProblem::boundarySolution(const field::NekFieldEnum & field, double * T)
796 : {
797 20410 : mesh_t * mesh = nekrs::entireMesh();
798 20410 : auto bc = _nek_mesh->boundaryCoupling();
799 :
800 : double (*f)(int);
801 20410 : f = nekrs::solutionPointer(field);
802 :
803 20410 : int start_1d = mesh->Nq;
804 20410 : int end_1d = _moose_Nq;
805 20410 : int start_2d = start_1d * start_1d;
806 20410 : int end_2d = end_1d * end_1d;
807 :
808 20410 : int n_to_write = bc.n_faces * end_2d * _nek_mesh->nBuildPerSurfaceElem();
809 :
810 : // allocate temporary space:
811 : // - Ttmp: results of the search for each process
812 : // - Tface: scratch space for face solution to avoid reallocating a bunch (only used if
813 : // interpolating)
814 : // - scratch: scratch for the interpolatino process to avoid reallocating a bunch (only used if
815 : // interpolating0
816 20410 : double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
817 20410 : double * Tface = (double *)calloc(start_2d, sizeof(double));
818 20410 : double * scratch = (double *)calloc(start_1d * end_1d, sizeof(double));
819 :
820 : auto indices = _nek_mesh->cornerIndices();
821 :
822 : int c = 0;
823 2173114 : for (int k = 0; k < bc.total_n_faces; ++k)
824 : {
825 2152704 : if (bc.process[k] == nekrs::commRank())
826 : {
827 382640 : int i = bc.element[k];
828 382640 : int j = bc.face[k];
829 382640 : int offset = i * mesh->Nfaces * start_2d + j * start_2d;
830 :
831 1133568 : for (int build = 0; build < _nek_mesh->nBuildPerSurfaceElem(); ++build)
832 : {
833 750928 : if (_needs_interpolation)
834 : {
835 : // get the solution on the face
836 3623648 : for (int v = 0; v < start_2d; ++v)
837 : {
838 3529344 : int id = mesh->vmapM[offset + v];
839 3529344 : Tface[v] = f(id);
840 : }
841 :
842 : // and then interpolate it
843 94304 : nekrs::interpolateSurfaceFaceHex3D(
844 94304 : scratch, _interpolation_outgoing, Tface, start_1d, &(Ttmp[c]), end_1d);
845 94304 : c += end_2d;
846 : }
847 : else
848 : {
849 : // get the solution on the face - no need to interpolate
850 3283120 : for (int v = 0; v < end_2d; ++v, ++c)
851 : {
852 2626496 : int id = mesh->vmapM[offset + indices[build][v]];
853 2626496 : Ttmp[c] = f(id);
854 : }
855 : }
856 : }
857 : }
858 : }
859 :
860 : // dimensionalize the solution if needed
861 3495642 : for (int v = 0; v < n_to_write; ++v)
862 3475232 : Ttmp[v] = Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);
863 :
864 20410 : nekrs::allgatherv(bc.mirror_counts, Ttmp, T, end_2d);
865 :
866 : freePointer(Ttmp);
867 : freePointer(Tface);
868 : freePointer(scratch);
869 20410 : }
870 :
871 : void
872 1439356 : NekRSProblem::writeVolumeSolution(const int elem_id,
873 : const field::NekWriteEnum & field,
874 : double * T,
875 : const std::vector<double> * add)
876 : {
877 1439356 : mesh_t * mesh = nekrs::entireMesh();
878 : void (*write_solution)(int, dfloat);
879 1439356 : write_solution = nekrs::solutionPointer(field);
880 :
881 1439356 : auto vc = _nek_mesh->volumeCoupling();
882 1439356 : int id = vc.element[elem_id] * mesh->Np;
883 :
884 1439356 : if (_nek_mesh->exactMirror())
885 : {
886 : // can write directly into the NekRS solution
887 3807400 : for (int v = 0; v < mesh->Np; ++v)
888 : {
889 3681792 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
890 3681792 : write_solution(id + v, T[v] + extra);
891 : }
892 : }
893 : else
894 : {
895 : // need to interpolate onto the higher-order Nek mesh
896 1313748 : double * tmp = (double *)calloc(mesh->Np, sizeof(double));
897 :
898 1313748 : interpolateVolumeSolutionToNek(elem_id, T, tmp);
899 :
900 185363348 : for (int v = 0; v < mesh->Np; ++v)
901 : {
902 184049600 : double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
903 184049600 : write_solution(id + v, tmp[v] + extra);
904 : }
905 :
906 : freePointer(tmp);
907 : }
908 1439356 : }
909 :
910 : void
911 423168 : NekRSProblem::writeBoundarySolution(const int elem_id,
912 : const field::NekWriteEnum & field,
913 : double * T)
914 : {
915 423168 : mesh_t * mesh = nekrs::temperatureMesh();
916 : void (*write_solution)(int, dfloat);
917 423168 : write_solution = nekrs::solutionPointer(field);
918 :
919 423168 : const auto & bc = _nek_mesh->boundaryCoupling();
920 423168 : int offset = bc.element[elem_id] * mesh->Nfaces * mesh->Nfp + bc.face[elem_id] * mesh->Nfp;
921 :
922 423168 : if (_nek_mesh->exactMirror())
923 : {
924 : // can write directly into the NekRS solution
925 638144 : for (int i = 0; i < mesh->Nfp; ++i)
926 613600 : write_solution(mesh->vmapM[offset + i], T[i]);
927 : }
928 : else
929 : {
930 : // need to interpolate onto the higher-order Nek mesh
931 398624 : double * tmp = (double *)calloc(mesh->Nfp, sizeof(double));
932 398624 : interpolateBoundarySolutionToNek(T, tmp);
933 :
934 7322872 : for (int i = 0; i < mesh->Nfp; ++i)
935 6924248 : write_solution(mesh->vmapM[offset + i], tmp[i]);
936 :
937 : freePointer(tmp);
938 : }
939 423168 : }
940 :
941 : void
942 1313748 : NekRSProblem::interpolateVolumeSolutionToNek(const int elem_id,
943 : double * incoming_moose_value,
944 : double * outgoing_nek_value)
945 : {
946 1313748 : mesh_t * mesh = nekrs::entireMesh();
947 :
948 1313748 : nekrs::interpolateVolumeHex3D(
949 1313748 : _interpolation_incoming, incoming_moose_value, _moose_Nq, outgoing_nek_value, mesh->Nq);
950 1313748 : }
951 :
952 : void
953 398624 : NekRSProblem::interpolateBoundarySolutionToNek(double * incoming_moose_value,
954 : double * outgoing_nek_value)
955 : {
956 398624 : mesh_t * mesh = nekrs::temperatureMesh();
957 :
958 398624 : double * scratch = (double *)calloc(_moose_Nq * mesh->Nq, sizeof(double));
959 :
960 398624 : nekrs::interpolateSurfaceFaceHex3D(scratch,
961 398624 : _interpolation_incoming,
962 : incoming_moose_value,
963 : _moose_Nq,
964 : outgoing_nek_value,
965 : mesh->Nq);
966 :
967 : freePointer(scratch);
968 398624 : }
969 :
970 : void
971 32372 : NekRSProblem::copyScratchToDevice()
972 : {
973 50126 : for (const auto & slot : _usrwrk_slots)
974 17754 : copyIndividualScratchSlot(slot);
975 :
976 32372 : if (nekrs::hasMovingMesh())
977 1504 : nekrs::copyDeformationToDevice();
978 32372 : }
979 :
980 : bool
981 35 : NekRSProblem::isUsrWrkSlotReservedForCoupling(const unsigned int & slot) const
982 : {
983 35 : return std::find(_usrwrk_slots.begin(), _usrwrk_slots.end(), slot) != _usrwrk_slots.end();
984 : }
985 :
986 : void
987 17781 : NekRSProblem::copyIndividualScratchSlot(const unsigned int & slot) const
988 : {
989 17781 : auto n = nekrs::fieldOffset();
990 17781 : auto nbytes = n * sizeof(dfloat);
991 :
992 17781 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
993 17781 : nrs->o_usrwrk.copyFrom(nrs->usrwrk + slot * n, nbytes, slot * nbytes);
994 17781 : }
995 :
996 : void
997 423168 : NekRSProblem::mapFaceDataToNekFace(const unsigned int & e,
998 : const unsigned int & var_num,
999 : const Real & divisor_scale,
1000 : const Real & additive_scale,
1001 : double ** outgoing_data)
1002 : {
1003 423168 : auto sys_number = _aux->number();
1004 423168 : auto & mesh = _nek_mesh->getMesh();
1005 423168 : auto indices = _nek_mesh->cornerIndices();
1006 :
1007 1214496 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
1008 : {
1009 791328 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
1010 :
1011 : // Only work on elements we can find on our local chunk of a
1012 : // distributed mesh
1013 791328 : if (!elem_ptr)
1014 : {
1015 : libmesh_assert(!mesh.is_serial());
1016 0 : continue;
1017 : }
1018 :
1019 4628640 : for (unsigned int n = 0; n < _n_vertices_per_surface; n++)
1020 : {
1021 : auto node_ptr = elem_ptr->node_ptr(n);
1022 :
1023 : // convert libMesh node index into the ordering used by NekRS
1024 3837312 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->boundaryNodeIndex(n)]
1025 2266496 : : _nek_mesh->boundaryNodeIndex(n);
1026 :
1027 3837312 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
1028 3837312 : (*outgoing_data)[node_index] =
1029 3837312 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
1030 : }
1031 : }
1032 423168 : }
1033 :
1034 : void
1035 1039536 : NekRSProblem::mapFaceDataToNekVolume(const unsigned int & e,
1036 : const unsigned int & var_num,
1037 : const Real & divisor_scale,
1038 : const Real & additive_scale,
1039 : double ** outgoing_data)
1040 : {
1041 1039536 : auto sys_number = _aux->number();
1042 1039536 : auto & mesh = _nek_mesh->getMesh();
1043 1039536 : auto indices = _nek_mesh->cornerIndices();
1044 :
1045 2491232 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
1046 : {
1047 1451696 : int n_faces_on_boundary = _nek_mesh->facesOnBoundary(e);
1048 :
1049 : // the only meaningful values are on the coupling boundaries, so we can skip this
1050 : // interpolation if this volume element isn't on a coupling boundary
1051 1451696 : if (n_faces_on_boundary > 0)
1052 : {
1053 288568 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
1054 :
1055 : // Only work on elements we can find on our local chunk of a
1056 : // distributed mesh
1057 288568 : if (!elem_ptr)
1058 : {
1059 : libmesh_assert(!mesh.is_serial());
1060 0 : continue;
1061 : }
1062 :
1063 5606712 : for (unsigned int n = 0; n < _n_vertices_per_volume; ++n)
1064 : {
1065 : auto node_ptr = elem_ptr->node_ptr(n);
1066 :
1067 : // convert libMesh node index into the ordering used by NekRS
1068 5318144 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
1069 5154304 : : _nek_mesh->volumeNodeIndex(n);
1070 :
1071 5318144 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
1072 5318144 : (*outgoing_data)[node_index] =
1073 5318144 : ((*_serialized_solution)(dof_idx)-additive_scale) / divisor_scale;
1074 : }
1075 : }
1076 : }
1077 1039536 : }
1078 :
1079 : void
1080 399820 : NekRSProblem::mapVolumeDataToNekVolume(const unsigned int & e,
1081 : const unsigned int & var_num,
1082 : const Real & divisor,
1083 : const Real & additive,
1084 : double ** outgoing_data)
1085 : {
1086 399820 : auto sys_number = _aux->number();
1087 399820 : auto & mesh = _nek_mesh->getMesh();
1088 399820 : auto indices = _nek_mesh->cornerIndices();
1089 :
1090 1415848 : for (int build = 0; build < _nek_mesh->nMoosePerNek(); ++build)
1091 : {
1092 1016028 : auto elem_ptr = mesh.query_elem_ptr(e * _nek_mesh->nMoosePerNek() + build);
1093 :
1094 : // Only work on elements we can find on our local chunk of a
1095 : // distributed mesh
1096 1016028 : if (!elem_ptr)
1097 : {
1098 : libmesh_assert(!mesh.is_serial());
1099 0 : continue;
1100 : }
1101 10625416 : for (unsigned int n = 0; n < _n_vertices_per_volume; n++)
1102 : {
1103 : auto node_ptr = elem_ptr->node_ptr(n);
1104 :
1105 : // convert libMesh node index into the ordering used by NekRS
1106 9609388 : int node_index = _nek_mesh->exactMirror() ? indices[build][_nek_mesh->volumeNodeIndex(n)]
1107 4145900 : : _nek_mesh->volumeNodeIndex(n);
1108 :
1109 9609388 : auto dof_idx = node_ptr->dof_number(sys_number, var_num, 0);
1110 9609388 : (*outgoing_data)[node_index] = ((*_serialized_solution)(dof_idx)-additive) / divisor;
1111 : }
1112 : }
1113 399820 : }
1114 :
1115 : #endif
|
