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 "NekInterface.h"
22 : #include "CardinalUtils.h"
23 :
24 : static nekrs::characteristicScales scales;
25 : static unsigned int n_usrwrk_slots;
26 : static bool is_nondimensional;
27 :
28 : /// Host arrays for essential NekRS fields
29 : static dfloat * sgeo;
30 : static dfloat * vgeo;
31 : static dfloat * x;
32 : static dfloat * y;
33 : static dfloat * z;
34 : static std::vector<dfloat> U;
35 : static std::vector<dfloat> P;
36 : static std::vector<dfloat> S;
37 : static dfloat * wall_distance;
38 :
39 : dfloat * usrwrk = nullptr;
40 :
41 : nrs_t * nrs;
42 :
43 : namespace nekrs
44 : {
45 : static double setup_time;
46 :
47 : // various constants for controlling tolerances
48 : static double abs_tol;
49 : static double rel_tol;
50 :
51 : void
52 295 : setAbsoluteTol(double tol)
53 : {
54 295 : abs_tol = tol;
55 295 : }
56 :
57 : void
58 295 : setRelativeTol(double tol)
59 : {
60 295 : rel_tol = tol;
61 295 : }
62 :
63 : void
64 764 : setNekSetupTime(const double & time)
65 : {
66 764 : setup_time = time;
67 764 : }
68 :
69 : double
70 749 : getNekSetupTime()
71 : {
72 749 : return setup_time;
73 : }
74 :
75 : void
76 686 : setStartTime(const double & start)
77 : {
78 1372 : platform->options.setArgs("START TIME", to_string_f(start));
79 686 : }
80 :
81 : void
82 62 : write_usrwrk_field_file(const int & usrWriterSize,
83 : const int & usrWriterIndex,
84 : const int & slot,
85 : const std::string & prefix,
86 : const dfloat & time,
87 : const int & step,
88 : const bool & write_coords)
89 : {
90 62 : static std::vector<std::unique_ptr<iofld>> usrWriterVector(usrWriterSize);
91 :
92 62 : auto & usrWriter = usrWriterVector[usrWriterIndex];
93 :
94 62 : if (!usrWriter)
95 : {
96 60 : usrWriter = iofldFactory::create();
97 30 : auto mesh = entireMesh();
98 :
99 30 : usrWriter->open(mesh, iofld::mode::write, prefix.c_str());
100 :
101 60 : if (platform->options.compareArgs("CHECKPOINT PRECISION", "FP32"))
102 : {
103 60 : usrWriter->writeAttribute("precision", "32");
104 : }
105 : else
106 : {
107 0 : usrWriter->writeAttribute("precision", "64");
108 : }
109 30 : usrWriter->addVariable("scalar00",
110 90 : std::vector<deviceMemory<dfloat>>{platform->app->bc->o_usrwrk.slice(
111 30 : slot * fieldOffset(), mesh->Nlocal)});
112 : }
113 :
114 156 : usrWriter->writeAttribute("outputmesh", write_coords ? "true" : "false");
115 124 : usrWriter->addVariable("time", const_cast<double &>(time));
116 62 : usrWriter->process();
117 92 : }
118 :
119 : void
120 56 : write_field_file(const std::string & prefix, const dfloat time, const int & step)
121 : {
122 56 : auto checkpointWriter = iofldFactory::create();
123 :
124 56 : auto visMesh = nrs->meshT;
125 56 : checkpointWriter->open(visMesh, iofld::mode::write, prefix.c_str());
126 :
127 112 : if (nrs->fluid &&
128 336 : platform->options.compareArgs(upperCase(nrs->fluid->name) + " CHECKPOINTING", "TRUE"))
129 : {
130 : std::vector<occa::memory> o_V;
131 224 : for (int i = 0; i < nrs->meshV->dim; i++)
132 336 : o_V.push_back(nrs->fluid->o_U.slice(i * nrs->fluid->fieldOffset, visMesh->Nlocal));
133 :
134 112 : checkpointWriter->addVariable("velocity", o_V);
135 :
136 112 : std::vector<occa::memory> o_p = {nrs->fluid->o_P.slice(0, visMesh->Nlocal)};
137 112 : checkpointWriter->addVariable("pressure", o_p);
138 56 : }
139 :
140 124 : for (int i = 0; i < nrs->Nscalar; i++)
141 : {
142 204 : if (platform->options.compareArgs("SCALAR" + scalarDigitStr(i) + " CHECKPOINTING", "TRUE"))
143 : {
144 68 : const auto temperatureExists = hasTemperatureVariable();
145 :
146 : std::vector<occa::memory> o_Si = {
147 136 : nrs->scalar->o_S.slice(nrs->scalar->fieldOffsetScan[i], visMesh->Nlocal)};
148 :
149 68 : if (i == 0 && temperatureExists)
150 112 : checkpointWriter->addVariable("temperature", o_Si);
151 : else
152 24 : checkpointWriter->addVariable("scalar" + scalarDigitStr(temperatureExists ? i - 1 : i),
153 : o_Si);
154 68 : }
155 : }
156 :
157 : int N;
158 112 : platform->options.getArgs("POLYNOMIAL DEGREE", N);
159 112 : checkpointWriter->writeAttribute("polynomialOrder", std::to_string(N));
160 :
161 168 : checkpointWriter->writeAttribute(
162 112 : "precision", platform->options.compareArgs("CHECKPOINT PRECISION", "FP64") ? "64" : "32");
163 :
164 168 : checkpointWriter->writeAttribute(
165 : "outputMesh",
166 112 : platform->options.compareArgs("CHECKPOINT OUTPUT MESH", "TRUE") ? "true" : "false");
167 :
168 112 : checkpointWriter->addVariable("time", const_cast<double &>(time));
169 56 : checkpointWriter->process();
170 56 : checkpointWriter->close();
171 180 : }
172 :
173 : void
174 766 : buildOnly(int buildOnly)
175 : {
176 766 : build_only = buildOnly;
177 766 : }
178 :
179 : int
180 66234 : buildOnly()
181 : {
182 66234 : return build_only;
183 : }
184 :
185 : bool
186 6968 : hasCHT()
187 : {
188 6968 : auto ns = Nscalar();
189 13884 : for (int is = 0; is < ns; is++)
190 : {
191 21294 : if (platform->options.compareArgs("SCALAR" + scalarDigitStr(is) + " MESH", "SOLID"))
192 : {
193 : return true;
194 : }
195 : }
196 : return false;
197 : }
198 :
199 : bool
200 22718 : hasMovingMesh()
201 : {
202 22718 : return platform->options.compareArgs("MOVING MESH", "TRUE");
203 : }
204 :
205 : bool
206 20384 : hasVariableDt()
207 : {
208 20384 : return platform->options.compareArgs("VARIABLE DT", "TRUE");
209 : }
210 :
211 : bool
212 9 : hasBlendingSolver()
213 : {
214 18 : return !platform->options.compareArgs("GEOM SOLVER", "NONE") && hasMovingMesh();
215 : }
216 :
217 : bool
218 2 : hasUserMeshSolver()
219 : {
220 4 : return platform->options.compareArgs("GEOM SOLVER", "NONE") && hasMovingMesh();
221 : }
222 :
223 : bool
224 1422 : endControlElapsedTime()
225 : {
226 2844 : return !platform->options.getArgs("STOP AT ELAPSED TIME").empty();
227 : }
228 :
229 : bool
230 1422 : endControlTime()
231 : {
232 1422 : return endTime() > 0;
233 : }
234 :
235 : bool
236 711 : endControlNumSteps()
237 : {
238 711 : return !endControlElapsedTime() && !endControlTime();
239 : }
240 :
241 : bool
242 164452221 : hasTemperatureVariable()
243 : {
244 483693317 : return Nscalar() ? platform->options.compareArgs("SCALAR00 NAME", "TEMPERATURE") : false;
245 : }
246 :
247 : bool
248 539 : hasTemperatureSolve()
249 : {
250 539 : return hasTemperatureVariable() ? nrs->scalar->compute[0] : false;
251 : }
252 :
253 : bool
254 1255 : hasScalarVariable(int scalarId)
255 : {
256 1255 : if (hasTemperatureVariable())
257 1148 : return scalarId < nrs->Nscalar;
258 : else
259 107 : return scalarId <= nrs->Nscalar;
260 : }
261 :
262 : bool
263 69 : hasHeatSourceKernel()
264 : {
265 69 : return static_cast<bool>(nrs->userSource);
266 : }
267 :
268 : bool
269 762 : isInitialized()
270 : {
271 762 : return nrs;
272 : }
273 :
274 : int
275 1128944493 : scalarFieldOffset()
276 : {
277 1128944493 : return nrs->scalar->fieldOffset(); // same for all scalars
278 : }
279 :
280 : int
281 910966 : velocityFieldOffset()
282 : {
283 910966 : return nrs->fluid->fieldOffset;
284 : }
285 :
286 : int
287 7204063 : fieldOffset()
288 : {
289 7204063 : if (hasTemperatureVariable())
290 7190541 : return scalarFieldOffset();
291 : else
292 13522 : return velocityFieldOffset();
293 : }
294 :
295 : mesh_t *
296 150619531 : entireMesh()
297 : {
298 150619531 : if (hasTemperatureVariable())
299 137930561 : return temperatureMesh();
300 : else
301 12688970 : return flowMesh();
302 : }
303 :
304 : mesh_t *
305 13879084 : flowMesh()
306 : {
307 13879084 : return nrs->meshV;
308 : }
309 :
310 : mesh_t *
311 138435268 : temperatureMesh()
312 : {
313 276870536 : return nrs->scalar->mesh("temperature");
314 : }
315 :
316 : mesh_t *
317 329754 : getMesh(const nek_mesh::NekMeshEnum pp_mesh)
318 : {
319 329754 : if (pp_mesh == nek_mesh::fluid)
320 1354 : return flowMesh();
321 328400 : else if (pp_mesh == nek_mesh::all)
322 328399 : return entireMesh();
323 : else
324 1 : mooseError("This object does not support operations on the solid part of the NekRS mesh!\n"
325 : "Valid options for 'mesh' are 'fluid' or 'all'.");
326 : }
327 :
328 : int
329 10412580 : commRank()
330 : {
331 10412580 : return platform->comm.mpiRank();
332 : }
333 :
334 : int
335 413537 : commSize()
336 : {
337 413537 : return platform->comm.mpiCommSize();
338 : }
339 :
340 : bool
341 764 : scratchAvailable()
342 : {
343 :
344 : // Because these scratch spaces are available for whatever the user sees fit, it is
345 : // possible that the user wants to use these arrays for a _different_ purpose aside from
346 : // transferring in MOOSE values. In nekrs::setup, we call the UDF_Setup0, UDF_Setup,
347 : // and UDF_ExecuteStep routines. These scratch space arrays aren't initialized anywhere
348 : // else in the core base, so we will make sure to throw an error from MOOSE if these
349 : // arrays are already in use, because otherwise our MOOSE transfer might get overwritten
350 : // by whatever other operation the user is trying to do.
351 764 : if (platform->app->bc->o_usrwrk.size() != 0)
352 1 : return false;
353 :
354 : return true;
355 : }
356 :
357 : void
358 763 : initializeScratch(const unsigned int & n_slots)
359 : {
360 763 : if (n_slots == 0)
361 : return;
362 :
363 : // In order to make indexing simpler in the device user functions (which is where the
364 : // boundary conditions are then actually applied), we define these scratch arrays
365 : // as volume arrays.
366 354 : usrwrk = (dfloat *)calloc(n_slots * fieldOffset(), sizeof(dfloat));
367 354 : platform->app->bc->o_usrwrk.resize(n_slots * fieldOffset());
368 :
369 354 : n_usrwrk_slots = n_slots;
370 : }
371 :
372 : void
373 678 : freeScratch()
374 : {
375 678 : freePointer(x);
376 678 : freePointer(y);
377 678 : freePointer(z);
378 678 : freePointer(sgeo);
379 678 : freePointer(vgeo);
380 678 : freePointer(usrwrk);
381 678 : freePointer(wall_distance);
382 678 : }
383 :
384 : double
385 276 : viscosity()
386 : {
387 : dfloat mu;
388 276 : setupAide & options = platform->options;
389 276 : options.getArgs("FLUID VISCOSITY", mu);
390 :
391 : // because we set rho_ref, U_ref, and L_ref all equal to 1 if our input is dimensional,
392 : // we don't need to have separate treatments for dimensional vs. nondimensional cases
393 276 : dfloat Re = 1.0 / mu;
394 276 : return scales.rho_ref * scales.U_ref * scales.L_ref / Re;
395 : }
396 :
397 : double
398 92 : Pr()
399 : {
400 : dfloat rho, rho_cp, k;
401 92 : setupAide & options = platform->options;
402 92 : options.getArgs("FLUID DENSITY", rho);
403 92 : options.getArgs("SCALAR00 TRANSPORTCOEFF", rho_cp);
404 92 : options.getArgs("SCALAR00 DIFFUSIONCOEFF", k);
405 :
406 92 : dfloat Pe = 1.0 / k;
407 92 : dfloat conductivity = scales.rho_ref * scales.U_ref * scales.Cp_ref * scales.L_ref / Pe;
408 92 : dfloat Cp = rho_cp / rho * scales.Cp_ref;
409 92 : return viscosity() * Cp / conductivity;
410 : }
411 :
412 : void
413 1494 : interpolationMatrix(double * I, int starting_points, int ending_points)
414 : {
415 1494 : DegreeRaiseMatrix1D(starting_points - 1, ending_points - 1, I);
416 1494 : }
417 :
418 : void
419 944150 : interpolateVolumeHex3D(const double * I, double * x, int N, double * Ix, int M)
420 : {
421 944150 : double * Ix1 = (dfloat *)calloc(N * N * M, sizeof(double));
422 944150 : double * Ix2 = (dfloat *)calloc(N * M * M, sizeof(double));
423 :
424 4501628 : for (int k = 0; k < N; ++k)
425 23418372 : for (int j = 0; j < N; ++j)
426 85493774 : for (int i = 0; i < M; ++i)
427 : {
428 : dfloat tmp = 0;
429 498403104 : for (int n = 0; n < N; ++n)
430 432770224 : tmp += I[i * N + n] * x[k * N * N + j * N + n];
431 65632880 : Ix1[k * N * M + j * M + i] = tmp;
432 : }
433 :
434 4501628 : for (int k = 0; k < N; ++k)
435 16791286 : for (int j = 0; j < M; ++j)
436 73875960 : for (int i = 0; i < M; ++i)
437 : {
438 : dfloat tmp = 0;
439 306733848 : for (int n = 0; n < N; ++n)
440 246091696 : tmp += I[j * N + n] * Ix1[k * N * M + n * M + i];
441 60642152 : Ix2[k * M * M + j * M + i] = tmp;
442 : }
443 :
444 4902712 : for (int k = 0; k < M; ++k)
445 25086012 : for (int j = 0; j < M; ++j)
446 153249168 : for (int i = 0; i < M; ++i)
447 : {
448 : dfloat tmp = 0;
449 470844510 : for (int n = 0; n < N; ++n)
450 338722792 : tmp += I[k * N + n] * Ix2[n * M * M + j * M + i];
451 132121718 : Ix[k * M * M + j * M + i] = tmp;
452 : }
453 :
454 : freePointer(Ix1);
455 : freePointer(Ix2);
456 944150 : }
457 :
458 : void
459 286550 : interpolateSurfaceFaceHex3D(
460 : double * scratch, const double * I, double * x, int N, double * Ix, int M)
461 : {
462 1153666 : for (int j = 0; j < N; ++j)
463 4199348 : for (int i = 0; i < M; ++i)
464 : {
465 : double tmp = 0;
466 17016984 : for (int n = 0; n < N; ++n)
467 : {
468 13684752 : tmp += I[i * N + n] * x[j * N + n];
469 : }
470 3332232 : scratch[j * M + i] = tmp;
471 : }
472 :
473 1511642 : for (int j = 0; j < M; ++j)
474 6667404 : for (int i = 0; i < M; ++i)
475 : {
476 : double tmp = 0;
477 18973080 : for (int n = 0; n < N; ++n)
478 : {
479 13530768 : tmp += I[j * N + n] * scratch[n * M + i];
480 : }
481 5442312 : Ix[j * M + i] = tmp;
482 : }
483 286550 : }
484 :
485 : void
486 36741 : displacementAndCounts(const std::vector<int> & base_counts,
487 : int * counts,
488 : int * displacement,
489 : const int multiplier = 1.0)
490 : {
491 187582 : for (int i = 0; i < commSize(); ++i)
492 150841 : counts[i] = base_counts[i] * multiplier;
493 :
494 36741 : displacement[0] = 0;
495 150841 : for (int i = 1; i < commSize(); i++)
496 114100 : displacement[i] = displacement[i - 1] + counts[i - 1];
497 36741 : }
498 :
499 : double
500 1756 : usrwrkVolumeIntegral(const unsigned int & slot, const nek_mesh::NekMeshEnum pp_mesh)
501 : {
502 1756 : const auto & mesh = getMesh(pp_mesh);
503 :
504 1756 : double integral = 0.0;
505 :
506 761632 : for (int k = 0; k < mesh->Nelements; ++k)
507 : {
508 759876 : int offset = k * mesh->Np;
509 :
510 199754692 : for (int v = 0; v < mesh->Np; ++v)
511 198994816 : integral += usrwrk[slot + offset + v] * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
512 : }
513 :
514 : // sum across all processes
515 : double total_integral;
516 1756 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
517 :
518 1756 : return total_integral;
519 : }
520 :
521 : void
522 878 : scaleUsrwrk(const unsigned int & slot, const dfloat & value)
523 : {
524 878 : mesh_t * mesh = getMesh(nek_mesh::all);
525 :
526 380816 : for (int k = 0; k < mesh->Nelements; ++k)
527 : {
528 379938 : int id = k * mesh->Np;
529 :
530 99877346 : for (int v = 0; v < mesh->Np; ++v)
531 99497408 : usrwrk[slot + id + v] *= value;
532 : }
533 878 : }
534 :
535 : std::vector<double>
536 24227 : usrwrkSideIntegral(const unsigned int & slot,
537 : const std::vector<int> & boundary,
538 : const nek_mesh::NekMeshEnum pp_mesh)
539 : {
540 24227 : const auto & mesh = getMesh(pp_mesh);
541 :
542 24227 : std::vector<double> integral(boundary.size(), 0.0);
543 :
544 7445542 : for (int i = 0; i < mesh->Nelements; ++i)
545 : {
546 51949205 : for (int j = 0; j < mesh->Nfaces; ++j)
547 : {
548 44527890 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
549 :
550 44527890 : if (std::find(boundary.begin(), boundary.end(), face_id) != boundary.end())
551 : {
552 731388 : auto it = std::find(boundary.begin(), boundary.end(), face_id);
553 : auto b_index = it - boundary.begin();
554 :
555 731388 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
556 :
557 12791404 : for (int v = 0; v < mesh->Nfp; ++v)
558 12060016 : integral[b_index] +=
559 12060016 : usrwrk[slot + mesh->vmapM[offset + v]] * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
560 : }
561 : }
562 : }
563 :
564 : // sum across all processes; this can probably be done more efficiently
565 24227 : std::vector<double> total_integral(boundary.size(), 0.0);
566 48575 : for (std::size_t i = 0; i < boundary.size(); ++i)
567 24348 : MPI_Allreduce(
568 : &integral[i], &total_integral[i], 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
569 :
570 24227 : return total_integral;
571 24227 : }
572 :
573 : void
574 0 : limitTemperature(const double * min_T, const double * max_T)
575 : {
576 : // if no limiters are provided, simply return
577 0 : if (!min_T && !max_T)
578 0 : return;
579 :
580 0 : double minimum = min_T ? *min_T : std::numeric_limits<double>::min();
581 0 : double maximum = max_T ? *max_T : std::numeric_limits<double>::max();
582 :
583 : // nondimensionalize if necessary
584 0 : minimum = (minimum - scales.T_ref) / scales.dT_ref;
585 0 : maximum = (maximum - scales.T_ref) / scales.dT_ref;
586 :
587 0 : mesh_t * mesh = temperatureMesh();
588 :
589 0 : const auto sid = nrs->scalar->nameToIndex.find("temperature")->second;
590 0 : const auto offset = scalarFieldOffset();
591 :
592 0 : for (int i = 0; i < mesh->Nelements; ++i)
593 : {
594 0 : for (int j = 0; j < mesh->Np; ++j)
595 : {
596 0 : int id = i * mesh->Np + j + sid * offset;
597 :
598 0 : if (S[id] < minimum)
599 0 : S[id] = minimum;
600 0 : if (S[id] > maximum)
601 0 : S[id] = maximum;
602 : }
603 : }
604 :
605 : // when complete, copy to device
606 0 : auto o_temperature = nrs->scalar->o_solution("temperature");
607 0 : o_temperature.copyFrom(S, offset, sid * offset, 0);
608 : }
609 :
610 : void
611 21089 : copySolutionToHost()
612 : {
613 21089 : nrs->fluid->o_U.copyTo(U.data(), U.size());
614 21089 : nrs->fluid->o_P.copyTo(P.data(), P.size());
615 :
616 21089 : if (Nscalar())
617 20164 : nrs->scalar->o_S.copyTo(S.data(), S.size());
618 21089 : }
619 :
620 : void
621 761 : copyMeshToHost()
622 : {
623 761 : mesh_t * mesh = entireMesh();
624 761 : mesh->o_x.copyTo(x);
625 761 : mesh->o_y.copyTo(y);
626 761 : mesh->o_z.copyTo(z);
627 761 : }
628 :
629 : void
630 0 : copyDeformationToDevice()
631 : {
632 0 : mesh_t * mesh = entireMesh();
633 :
634 0 : mesh->o_x.copyFrom(x);
635 0 : mesh->o_y.copyFrom(y);
636 0 : mesh->o_z.copyFrom(z);
637 0 : mesh->update();
638 :
639 0 : updateHostMeshParameters();
640 0 : }
641 :
642 : void
643 761 : initializeHostMeshParameters()
644 : {
645 761 : mesh_t * mesh = entireMesh();
646 761 : sgeo = (dfloat *)calloc(mesh->o_sgeo.size(), sizeof(dfloat));
647 761 : vgeo = (dfloat *)calloc(mesh->o_vgeo.size(), sizeof(dfloat));
648 761 : x = (dfloat *)calloc(mesh->o_x.size(), sizeof(dfloat));
649 761 : y = (dfloat *)calloc(mesh->o_y.size(), sizeof(dfloat));
650 761 : z = (dfloat *)calloc(mesh->o_z.size(), sizeof(dfloat));
651 :
652 761 : copyMeshToHost();
653 761 : }
654 :
655 : void
656 761 : updateHostMeshParameters()
657 : {
658 761 : mesh_t * mesh = entireMesh();
659 761 : mesh->o_sgeo.copyTo(sgeo);
660 761 : mesh->o_vgeo.copyTo(vgeo);
661 761 : }
662 :
663 : dfloat *
664 152 : getSgeo()
665 : {
666 152 : return sgeo;
667 : }
668 :
669 : dfloat *
670 355 : getVgeo()
671 : {
672 355 : return vgeo;
673 : }
674 :
675 : double
676 23972 : sideExtremeValue(const std::vector<int> & boundary_id, const field::NekFieldEnum & field,
677 : const nek_mesh::NekMeshEnum pp_mesh, const bool max)
678 : {
679 23972 : mesh_t * mesh = getMesh(pp_mesh);
680 :
681 23972 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
682 :
683 : double (*f)(int, int);
684 23972 : f = solutionPointer(field);
685 :
686 2514208 : for (int i = 0; i < mesh->Nelements; ++i)
687 : {
688 17431652 : for (int j = 0; j < mesh->Nfaces; ++j)
689 : {
690 14941416 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
691 :
692 14941416 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
693 : {
694 261524 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
695 11557388 : for (int v = 0; v < mesh->Nfp; ++v)
696 : {
697 11295864 : if (max)
698 5724844 : value = std::max(value, f(mesh->vmapM[offset + v], 0 /* unused */));
699 : else
700 5698298 : value = std::min(value, f(mesh->vmapM[offset + v], 0 /* unused */));
701 : }
702 : }
703 : }
704 : }
705 :
706 : // find extreme value across all processes
707 : double reduced_value;
708 23972 : auto op = max ? MPI_MAX : MPI_MIN;
709 23972 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm());
710 :
711 : // dimensionalize the field if needed
712 23972 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
713 :
714 23972 : return reduced_value;
715 : }
716 :
717 : double
718 36284 : volumeExtremeValue(const field::NekFieldEnum & field, const nek_mesh::NekMeshEnum pp_mesh, const bool max)
719 : {
720 36284 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
721 :
722 : double (*f)(int, int);
723 36284 : f = solutionPointer(field);
724 :
725 : mesh_t * mesh;
726 : int start_id;
727 :
728 36284 : switch (pp_mesh)
729 : {
730 36272 : case nek_mesh::fluid:
731 : case nek_mesh::all:
732 : {
733 36272 : mesh = getMesh(pp_mesh);
734 : start_id = 0;
735 : break;
736 : }
737 12 : case nek_mesh::solid:
738 : {
739 12 : mesh = entireMesh();
740 12 : start_id = flowMesh()->Nelements;
741 12 : break;
742 : }
743 0 : default:
744 0 : mooseError("Unhandled NekMeshEnum in volumeExtremeValue");
745 : }
746 :
747 7581654 : for (int i = start_id; i < mesh->Nelements; ++i)
748 : {
749 996250626 : for (int j = 0; j < mesh->Np; ++j)
750 : {
751 988705256 : if (max)
752 526435657 : value = std::max(value, f(i * mesh->Np + j, 0 /* unused */));
753 : else
754 464091746 : value = std::min(value, f(i * mesh->Np + j, 0 /* unused */));
755 : }
756 : }
757 :
758 : // find extreme value across all processes
759 : double reduced_value;
760 36284 : auto op = max ? MPI_MAX : MPI_MIN;
761 36284 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm());
762 :
763 : // dimensionalize the field if needed
764 36284 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
765 :
766 36284 : return reduced_value;
767 : }
768 :
769 : Point
770 99589536 : gllPoint(int local_elem_id, int local_node_id)
771 : {
772 99589536 : mesh_t * mesh = entireMesh();
773 :
774 99589536 : int id = local_elem_id * mesh->Np + local_node_id;
775 99589536 : Point p(x[id], y[id], z[id]);
776 : p *= scales.L_ref;
777 99589536 : return p;
778 : }
779 :
780 : Point
781 933120 : gllPointFace(int local_elem_id, int local_face_id, int local_node_id)
782 : {
783 933120 : mesh_t * mesh = entireMesh();
784 933120 : int face_id = mesh->EToB[local_elem_id * mesh->Nfaces + local_face_id];
785 933120 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
786 933120 : int id = mesh->vmapM[offset + local_node_id];
787 :
788 933120 : Point p(x[id], y[id], z[id]);
789 : p *= scales.L_ref;
790 933120 : return p;
791 : }
792 :
793 : Point
794 400128 : centroidFace(int local_elem_id, int local_face_id)
795 : {
796 400128 : mesh_t * mesh = entireMesh();
797 :
798 : double x_c = 0.0;
799 : double y_c = 0.0;
800 : double z_c = 0.0;
801 : double mass = 0.0;
802 :
803 400128 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
804 6571776 : for (int v = 0; v < mesh->Nfp; ++v)
805 : {
806 6171648 : int id = mesh->vmapM[offset + v];
807 6171648 : double mass_matrix = sgeo[mesh->Nsgeo * (offset + v) + WSJID];
808 6171648 : x_c += x[id] * mass_matrix;
809 6171648 : y_c += y[id] * mass_matrix;
810 6171648 : z_c += z[id] * mass_matrix;
811 6171648 : mass += mass_matrix;
812 : }
813 :
814 : Point c(x_c, y_c, z_c);
815 400128 : return c / mass * scales.L_ref;
816 : }
817 :
818 : Point
819 26945124 : centroid(int local_elem_id)
820 : {
821 26945124 : mesh_t * mesh = entireMesh();
822 :
823 : double x_c = 0.0;
824 : double y_c = 0.0;
825 : double z_c = 0.0;
826 : double mass = 0.0;
827 :
828 6163730928 : for (int v = 0; v < mesh->Np; ++v)
829 : {
830 6136785804 : int id = local_elem_id * mesh->Np + v;
831 6136785804 : double mass_matrix = vgeo[local_elem_id * mesh->Np * mesh->Nvgeo + JWID * mesh->Np + v];
832 6136785804 : x_c += x[id] * mass_matrix;
833 6136785804 : y_c += y[id] * mass_matrix;
834 6136785804 : z_c += z[id] * mass_matrix;
835 6136785804 : mass += mass_matrix;
836 : }
837 :
838 : Point c(x_c, y_c, z_c);
839 26945124 : return c / mass * scales.L_ref;
840 : }
841 :
842 : double
843 15130 : volume(const nek_mesh::NekMeshEnum pp_mesh)
844 : {
845 15130 : mesh_t * mesh = getMesh(pp_mesh);
846 15130 : double integral = 0.0;
847 :
848 2963308 : for (int k = 0; k < mesh->Nelements; ++k)
849 : {
850 2948178 : int offset = k * mesh->Np;
851 :
852 393276610 : for (int v = 0; v < mesh->Np; ++v)
853 390328432 : integral += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
854 : }
855 :
856 : // sum across all processes
857 : double total_integral;
858 15130 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
859 :
860 15130 : total_integral *= scales.V_ref;
861 :
862 15130 : return total_integral;
863 : }
864 :
865 : void
866 44275 : dimensionalizeVolume(double & integral)
867 : {
868 44275 : integral *= scales.V_ref;
869 44275 : }
870 :
871 : void
872 484 : dimensionalizeArea(double & integral)
873 : {
874 484 : integral *= scales.A_ref;
875 484 : }
876 :
877 : void
878 23460 : dimensionalizeVolumeIntegral(const field::NekFieldEnum & integrand,
879 : const Real & volume,
880 : double & integral)
881 : {
882 : // dimensionalize the field if needed
883 23460 : integral *= nondimensionalDivisor(integrand);
884 :
885 : // scale the volume integral
886 23460 : integral *= scales.V_ref;
887 :
888 : // for quantities with a relative scaling, we need to add back the reference
889 : // contribution to the volume integral
890 23460 : integral += nondimensionalAdditive(integrand) * volume;
891 23460 : }
892 :
893 : void
894 484 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
895 : const Real & area,
896 : double & integral)
897 : {
898 : // dimensionalize the field if needed
899 484 : integral *= nondimensionalDivisor(integrand);
900 :
901 : // scale the boundary integral
902 484 : integral *= scales.A_ref;
903 :
904 : // for quantities with a relative scaling, we need to add back the reference
905 : // contribution to the side integral
906 484 : integral += nondimensionalAdditive(integrand) * area;
907 484 : }
908 :
909 : void
910 28180 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
911 : const std::vector<int> & boundary_id,
912 : double & integral,
913 : const nek_mesh::NekMeshEnum pp_mesh)
914 : {
915 : // dimensionalize the field if needed
916 28180 : integral *= nondimensionalDivisor(integrand);
917 :
918 : // scale the boundary integral
919 28180 : integral *= scales.A_ref;
920 :
921 : // for quantities with a relative scaling, we need to add back the reference
922 : // contribution to the side integral; we need this form here to avoid a recursive loop
923 28180 : auto add = nondimensionalAdditive(integrand);
924 28180 : if (std::abs(add) > 1e-8)
925 0 : integral += add * area(boundary_id, pp_mesh);
926 28180 : }
927 :
928 : double
929 558144 : evaluateFunctionOnMesh(const Function * f, const Real time, const int id)
930 : {
931 : double shift = 0.0;
932 558144 : if (f)
933 : {
934 : // the function is given in dimensional form from MOOSE, so we need to
935 : // convert the x,y,z points we loop through on NekRS's mesh into the
936 : // dimensional form before passing them into the dimensional function
937 418608 : Point p(x[id], y[id], z[id]);
938 : p *= scales.L_ref;
939 418608 : auto t = time * scales.t_ref;
940 418608 : shift = f->value(t, p);
941 : }
942 :
943 558144 : return shift;
944 : }
945 :
946 : double
947 16 : volumeNorm(const field::NekFieldEnum & integrand,
948 : const nek_mesh::NekMeshEnum pp_mesh,
949 : const Function * function,
950 : const Real & time,
951 : const Real & N)
952 : {
953 16 : mesh_t * mesh = getMesh(pp_mesh);
954 :
955 : double (*f)(int, int);
956 16 : f = solutionPointer(integrand);
957 16 : double integral = 0.0;
958 16 : double total_integral = 0.0;
959 :
960 2600 : for (int k = 0; k < mesh->Nelements; ++k)
961 : {
962 2584 : const int offset = k * mesh->Np;
963 :
964 560728 : for (int v = 0; v < mesh->Np; ++v)
965 : {
966 558144 : const int n = offset + v;
967 558144 : const auto shift = evaluateFunctionOnMesh(function, time, n);
968 558144 : const double error = std::abs(f(n, 0 /* unused */) - shift);
969 :
970 558144 : if (std::isinf(N))
971 279072 : integral = std::max(integral, error);
972 : else
973 279072 : integral += std::pow(error, N) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
974 : }
975 : }
976 :
977 16 : auto reduction_type = std::isinf(N) ? MPI_MAX : MPI_SUM;
978 16 : MPI_Allreduce(
979 : &integral, &total_integral, 1, MPI_DOUBLE, reduction_type, platform->comm.mpiComm());
980 :
981 16 : return std::isinf(N) ? total_integral : std::pow(total_integral, 1.0 / N);
982 : }
983 :
984 : double
985 8948 : volumeIntegral(const field::NekFieldEnum & integrand, const Real & volume,
986 : const nek_mesh::NekMeshEnum pp_mesh)
987 : {
988 8948 : mesh_t * mesh = getMesh(pp_mesh);
989 :
990 8948 : double integral = 0.0;
991 :
992 : double (*f)(int, int);
993 8948 : f = solutionPointer(integrand);
994 :
995 1587410 : for (int k = 0; k < mesh->Nelements; ++k)
996 : {
997 1578462 : int offset = k * mesh->Np;
998 :
999 238927134 : for (int v = 0; v < mesh->Np; ++v)
1000 237348672 : integral += f(offset + v, 0 /* unused */) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
1001 : }
1002 :
1003 : // sum across all processes
1004 : double total_integral;
1005 8948 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1006 :
1007 8948 : dimensionalizeVolumeIntegral(integrand, volume, total_integral);
1008 :
1009 8948 : return total_integral;
1010 : }
1011 :
1012 : double
1013 12366 : area(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1014 : {
1015 12366 : mesh_t * mesh = getMesh(pp_mesh);
1016 :
1017 12366 : double integral = 0.0;
1018 :
1019 1695828 : for (int i = 0; i < mesh->Nelements; ++i)
1020 : {
1021 11784234 : for (int j = 0; j < mesh->Nfaces; ++j)
1022 : {
1023 10100772 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1024 :
1025 10100772 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1026 : {
1027 325424 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1028 6263964 : for (int v = 0; v < mesh->Nfp; ++v)
1029 : {
1030 5938540 : integral += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
1031 : }
1032 : }
1033 : }
1034 : }
1035 :
1036 : // sum across all processes
1037 : double total_integral;
1038 12366 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1039 :
1040 12366 : dimensionalizeSideIntegral(field::unity, boundary_id, total_integral, pp_mesh);
1041 :
1042 12366 : return total_integral;
1043 : }
1044 :
1045 : void
1046 14 : computeWallDistance(const std::vector<int> & boundary_id)
1047 : {
1048 14 : mesh_t * mesh = getMesh(nek_mesh::fluid);
1049 14 : freePointer(wall_distance);
1050 14 : auto o_wbID = platform->device.malloc<int>(boundary_id.size(), boundary_id.data());
1051 14 : auto o_ywd = mesh->minDistance(boundary_id.size(), o_wbID, "cheap_dist");
1052 14 : wall_distance = (dfloat *)calloc(o_ywd.size(), sizeof(dfloat));
1053 14 : o_ywd.copyTo(wall_distance);
1054 14 : o_ywd.free();
1055 14 : o_wbID.free();
1056 14 : }
1057 :
1058 : std::vector<dfloat>
1059 12 : yPlus(const std::vector<int> & boundary_id)
1060 : {
1061 12 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
1062 :
1063 : // integrate over the boundaries in the mesh; each rank will compute contributions to the
1064 : // x, y, and z components
1065 12 : mesh_t * mesh = getMesh(nek_mesh::fluid);
1066 :
1067 : // compute the rate of strain tensor on device
1068 12 : auto o_Sij = nrs->strainRate();
1069 :
1070 : // copy to host for evaluating wall-parallel stress
1071 12 : dfloat * Sij = (dfloat *)calloc(o_Sij.size(), sizeof(dfloat));
1072 12 : o_Sij.copyTo(Sij);
1073 12 : o_Sij.free();
1074 :
1075 : // TODO: This function only works correctly if the viscosity and rho is constant, because
1076 : // otherwise we need to copy the viscosity and rho from device to host
1077 : double mu;
1078 12 : platform->options.getArgs("FLUID VISCOSITY", mu);
1079 : double rho;
1080 12 : platform->options.getArgs("FLUID DENSITY", rho);
1081 12 : double nu = mu / rho;
1082 :
1083 12 : dfloat max_yp = -std::numeric_limits<float>::max();
1084 12 : dfloat min_yp = std::numeric_limits<float>::max();
1085 12 : dfloat avg_yp = 0.0;
1086 12 : dfloat denom_yp = 0.0;
1087 :
1088 : std::vector<int> istride = {
1089 12 : mesh->Nq * mesh->Nq, mesh->Nq, -1, -mesh->Nq, 1, -mesh->Nq * mesh->Nq};
1090 :
1091 1932 : for (int i = 0; i < mesh->Nelements; ++i)
1092 : {
1093 13440 : for (int j = 0; j < mesh->Nfaces; ++j)
1094 : {
1095 11520 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1096 :
1097 11520 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1098 : {
1099 384 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1100 :
1101 6528 : for (int v = 0; v < mesh->Nfp; ++v)
1102 : {
1103 6144 : int surf_offset = mesh->Nsgeo * (offset + v);
1104 6144 : int vol_id = mesh->vmapM[offset + v];
1105 6144 : dfloat sWJ = sgeo[surf_offset + WSJID];
1106 6144 : dfloat scale = 2 * mu;
1107 :
1108 6144 : dfloat n1 = sgeo[surf_offset + NXID];
1109 6144 : dfloat n2 = sgeo[surf_offset + NYID];
1110 6144 : dfloat n3 = sgeo[surf_offset + NZID];
1111 :
1112 6144 : dfloat s11 = Sij[vol_id + 0 * nrs->fieldOffset];
1113 6144 : dfloat s12 = Sij[vol_id + 1 * nrs->fieldOffset];
1114 6144 : dfloat s13 = Sij[vol_id + 2 * nrs->fieldOffset];
1115 : dfloat s21 = s12;
1116 6144 : dfloat s22 = Sij[vol_id + 3 * nrs->fieldOffset];
1117 6144 : dfloat s23 = Sij[vol_id + 4 * nrs->fieldOffset];
1118 : dfloat s31 = s13;
1119 : dfloat s32 = s23;
1120 6144 : dfloat s33 = Sij[vol_id + 5 * nrs->fieldOffset];
1121 :
1122 : // tau_{ij}n_j - (tau_{jk} n_k n_j) * n_i
1123 6144 : dfloat f1 = scale * (s11 * n1 + s12 * n2 + s13 * n3);
1124 6144 : dfloat f2 = scale * (s21 * n1 + s22 * n2 + s23 * n3);
1125 6144 : dfloat f3 = scale * (s31 * n1 + s32 * n2 + s33 * n3);
1126 :
1127 6144 : f1 -= f1 * n1 * n1;
1128 6144 : f2 -= f2 * n2 * n2;
1129 6144 : f3 -= f3 * n3 * n3;
1130 :
1131 6144 : dfloat tauw = sqrt(f1 * f1 + f2 * f2 + f3 * f3);
1132 6144 : dfloat utau = sqrt(tauw / rho);
1133 :
1134 : // need to shift when evaluating the wall distance, because we want the wall distance
1135 : // at the nearest node away from the face, not precisely on the face
1136 6144 : dfloat wd = wall_distance[vol_id + istride[j]];
1137 :
1138 : // check that we are not at a corner
1139 6144 : if (wd > 1e-8)
1140 : {
1141 6144 : dfloat yplus = wd * utau / nu;
1142 6144 : max_yp = std::max(yplus, max_yp);
1143 6144 : min_yp = std::min(yplus, min_yp);
1144 6144 : avg_yp += yplus * sWJ;
1145 6144 : denom_yp += sWJ;
1146 : }
1147 : }
1148 : }
1149 : }
1150 : }
1151 :
1152 : // max across all processes
1153 : double total_max_yp;
1154 12 : MPI_Allreduce(&max_yp, &total_max_yp, 1, MPI_DOUBLE, MPI_MAX, platform->comm.mpiComm());
1155 :
1156 : // min across all processes
1157 : double total_min_yp;
1158 12 : MPI_Allreduce(&min_yp, &total_min_yp, 1, MPI_DOUBLE, MPI_MIN, platform->comm.mpiComm());
1159 :
1160 : // sum across all processes
1161 : double total_avg_yp;
1162 12 : MPI_Allreduce(&avg_yp, &total_avg_yp, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1163 :
1164 12 : if (total_avg_yp == 0)
1165 0 : mooseError("Failed to find any eligible points on boundaries for computing y+!");
1166 :
1167 : double total_denom_yp;
1168 12 : MPI_Allreduce(&denom_yp, &total_denom_yp, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1169 :
1170 : // TODO: dimensionalize
1171 : // dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
1172 :
1173 : freePointer(Sij);
1174 24 : return {total_max_yp, total_min_yp, total_avg_yp / total_denom_yp};
1175 12 : }
1176 :
1177 : double
1178 15779 : sideIntegral(const std::vector<int> & boundary_id, const field::NekFieldEnum & integrand,
1179 : const nek_mesh::NekMeshEnum pp_mesh)
1180 : {
1181 15779 : mesh_t * mesh = getMesh(pp_mesh);
1182 :
1183 15778 : double integral = 0.0;
1184 :
1185 : double (*f)(int, int);
1186 15778 : f = solutionPointer(integrand);
1187 :
1188 1862106 : for (int i = 0; i < mesh->Nelements; ++i)
1189 : {
1190 12924296 : for (int j = 0; j < mesh->Nfaces; ++j)
1191 : {
1192 11077968 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1193 :
1194 11077968 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1195 : {
1196 359188 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1197 7760792 : for (int v = 0; v < mesh->Nfp; ++v)
1198 : {
1199 7401604 : integral +=
1200 7401604 : f(mesh->vmapM[offset + v], 0 /* unused */) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
1201 : }
1202 : }
1203 : }
1204 : }
1205 :
1206 : // sum across all processes
1207 : double total_integral;
1208 15778 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1209 :
1210 15778 : dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
1211 :
1212 15778 : return total_integral;
1213 : }
1214 :
1215 : double
1216 104 : massFlowrate(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1217 : {
1218 104 : mesh_t * mesh = getMesh(pp_mesh);
1219 :
1220 : // TODO: This function only works correctly if the density is constant, because
1221 : // otherwise we need to copy the density from device to host
1222 : double rho;
1223 104 : platform->options.getArgs("FLUID DENSITY", rho);
1224 :
1225 104 : double integral = 0.0;
1226 :
1227 16448 : for (int i = 0; i < mesh->Nelements; ++i)
1228 : {
1229 114408 : for (int j = 0; j < mesh->Nfaces; ++j)
1230 : {
1231 98064 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1232 :
1233 98064 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1234 : {
1235 1660 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1236 68140 : for (int v = 0; v < mesh->Nfp; ++v)
1237 : {
1238 66480 : int vol_id = mesh->vmapM[offset + v];
1239 66480 : int surf_offset = mesh->Nsgeo * (offset + v);
1240 :
1241 : double normal_velocity =
1242 66480 : U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1243 66480 : U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1244 66480 : U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1245 :
1246 66480 : integral += rho * normal_velocity * sgeo[surf_offset + WSJID];
1247 : }
1248 : }
1249 : }
1250 : }
1251 :
1252 : // sum across all processes
1253 : double total_integral;
1254 104 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1255 :
1256 : // dimensionalize the mass flux and area
1257 104 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1258 :
1259 104 : return total_integral;
1260 : }
1261 :
1262 : double
1263 436 : sideMassFluxWeightedIntegral(const std::vector<int> & boundary_id,
1264 : const field::NekFieldEnum & integrand,
1265 : const nek_mesh::NekMeshEnum pp_mesh)
1266 : {
1267 436 : mesh_t * mesh = getMesh(pp_mesh);
1268 :
1269 : // TODO: This function only works correctly if the density is constant, because
1270 : // otherwise we need to copy the density from device to host
1271 : double rho;
1272 436 : platform->options.getArgs("FLUID DENSITY", rho);
1273 :
1274 436 : double integral = 0.0;
1275 :
1276 : double (*f)(int, int);
1277 436 : f = solutionPointer(integrand);
1278 :
1279 60292 : for (int i = 0; i < mesh->Nelements; ++i)
1280 : {
1281 418992 : for (int j = 0; j < mesh->Nfaces; ++j)
1282 : {
1283 359136 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1284 :
1285 359136 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1286 : {
1287 5682 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1288 238350 : for (int v = 0; v < mesh->Nfp; ++v)
1289 : {
1290 232668 : int vol_id = mesh->vmapM[offset + v];
1291 232668 : int surf_offset = mesh->Nsgeo * (offset + v);
1292 : double normal_velocity =
1293 232668 : U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1294 232668 : U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1295 232668 : U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1296 232668 : integral += f(vol_id, 0 /* unused */) * rho * normal_velocity * sgeo[surf_offset + WSJID];
1297 : }
1298 : }
1299 : }
1300 : }
1301 :
1302 : // sum across all processes
1303 : double total_integral;
1304 436 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1305 :
1306 : // dimensionalize the field if needed
1307 436 : total_integral *= nondimensionalDivisor(integrand);
1308 :
1309 : // dimensionalize the mass flux and area
1310 436 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1311 :
1312 : // for quantities with a relative scaling, we need to add back the reference
1313 : // contribution to the mass flux integral; we need this form here to avoid an infinite
1314 : // recursive loop
1315 436 : auto add = nondimensionalAdditive(integrand);
1316 436 : if (std::abs(add) > 1e-8)
1317 0 : total_integral += add * massFlowrate(boundary_id, pp_mesh);
1318 :
1319 436 : return total_integral;
1320 : }
1321 :
1322 : double
1323 36 : pressureSurfaceForce(const std::vector<int> & boundary_id, const Point & direction, const nek_mesh::NekMeshEnum pp_mesh)
1324 : {
1325 36 : mesh_t * mesh = getMesh(pp_mesh);
1326 :
1327 36 : double integral = 0.0;
1328 :
1329 20232 : for (int i = 0; i < mesh->Nelements; ++i)
1330 : {
1331 141372 : for (int j = 0; j < mesh->Nfaces; ++j)
1332 : {
1333 121176 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1334 :
1335 121176 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1336 : {
1337 12780 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1338 472860 : for (int v = 0; v < mesh->Nfp; ++v)
1339 : {
1340 460080 : int vol_id = mesh->vmapM[offset + v];
1341 460080 : int surf_offset = mesh->Nsgeo * (offset + v);
1342 :
1343 460080 : double p_normal = P[vol_id] * (sgeo[surf_offset + NXID] * direction(0) +
1344 460080 : sgeo[surf_offset + NYID] * direction(1) +
1345 460080 : sgeo[surf_offset + NZID] * direction(2));
1346 :
1347 460080 : integral += p_normal * sgeo[surf_offset + WSJID];
1348 : }
1349 : }
1350 : }
1351 : }
1352 :
1353 : // sum across all processes
1354 : double total_integral;
1355 36 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1356 :
1357 36 : dimensionalizeSideIntegral(field::pressure, boundary_id, total_integral, pp_mesh);
1358 :
1359 36 : return total_integral;
1360 : }
1361 :
1362 : double
1363 12124 : heatFluxIntegral(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1364 : {
1365 12124 : mesh_t * mesh = getMesh(pp_mesh);
1366 :
1367 : // TODO: This function only works correctly if the conductivity is constant, because
1368 : // otherwise we need to copy the conductivity from device to host
1369 : double k;
1370 12124 : platform->options.getArgs("SCALAR00 DIFFUSIONCOEFF", k);
1371 :
1372 12124 : double integral = 0.0;
1373 12124 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1374 :
1375 12124 : const auto sid = nrs->scalar->nameToIndex.find("temperature")->second;
1376 12124 : const int offset = sid * scalarFieldOffset();
1377 :
1378 2728652 : for (int i = 0; i < mesh->Nelements; ++i)
1379 : {
1380 19015696 : for (int j = 0; j < mesh->Nfaces; ++j)
1381 : {
1382 16299168 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1383 :
1384 16299168 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1385 : {
1386 : // some inefficiency if an element has more than one face on the sideset of interest,
1387 : // because we will recompute the gradient in the element more than one time - but this
1388 : // is of little practical interest because this will be a minority of cases.
1389 169620 : gradient(mesh->Np, i, S.data() + offset, grad_T, pp_mesh);
1390 :
1391 169620 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1392 3431124 : for (int v = 0; v < mesh->Nfp; ++v)
1393 : {
1394 : // special use of vol_id only when calling gradient(...), since we have written this
1395 : // function internally here so that it computes the gradient in a given element (as
1396 : // opposed to doing so for all elements at once, like in NekRS's
1397 : // gradientVolumeHex3D kernel
1398 3261504 : int vol_id = mesh->vmapM[offset + v] - i * mesh->Np;
1399 3261504 : int surf_offset = mesh->Nsgeo * (offset + v);
1400 :
1401 3261504 : double normal_grad_T = grad_T[vol_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1402 3261504 : grad_T[vol_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1403 3261504 : grad_T[vol_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1404 :
1405 3261504 : integral += -k * normal_grad_T * sgeo[surf_offset + WSJID];
1406 : }
1407 : }
1408 : }
1409 : }
1410 :
1411 : freePointer(grad_T);
1412 :
1413 : // sum across all processes
1414 : double total_integral;
1415 12124 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm());
1416 :
1417 : // multiply by the reference heat flux and an area factor to dimensionalize
1418 12124 : total_integral *= scales.flux_ref * scales.A_ref;
1419 :
1420 12124 : return total_integral;
1421 : }
1422 :
1423 : void
1424 173652 : gradient(const int offset,
1425 : const int e,
1426 : const double * f,
1427 : double * grad_f,
1428 : const nek_mesh::NekMeshEnum pp_mesh)
1429 : {
1430 173652 : mesh_t * mesh = getMesh(pp_mesh);
1431 :
1432 913524 : for (int k = 0; k < mesh->Nq; ++k)
1433 : {
1434 4019424 : for (int j = 0; j < mesh->Nq; ++j)
1435 : {
1436 18439584 : for (int i = 0; i < mesh->Nq; ++i)
1437 : {
1438 15160032 : const int gid = e * mesh->Np * mesh->Nvgeo + k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1439 15160032 : const double drdx = vgeo[gid + RXID * mesh->Np];
1440 15160032 : const double drdy = vgeo[gid + RYID * mesh->Np];
1441 15160032 : const double drdz = vgeo[gid + RZID * mesh->Np];
1442 15160032 : const double dsdx = vgeo[gid + SXID * mesh->Np];
1443 15160032 : const double dsdy = vgeo[gid + SYID * mesh->Np];
1444 15160032 : const double dsdz = vgeo[gid + SZID * mesh->Np];
1445 15160032 : const double dtdx = vgeo[gid + TXID * mesh->Np];
1446 15160032 : const double dtdy = vgeo[gid + TYID * mesh->Np];
1447 15160032 : const double dtdz = vgeo[gid + TZID * mesh->Np];
1448 :
1449 : // compute 'r' and 's' derivatives of (q_m) at node n
1450 : double dpdr = 0.f, dpds = 0.f, dpdt = 0.f;
1451 :
1452 88564224 : for (int n = 0; n < mesh->Nq; ++n)
1453 : {
1454 73404192 : const double Dr = mesh->D[i * mesh->Nq + n];
1455 73404192 : const double Ds = mesh->D[j * mesh->Nq + n];
1456 73404192 : const double Dt = mesh->D[k * mesh->Nq + n];
1457 :
1458 73404192 : dpdr += Dr * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + j * mesh->Nq + n];
1459 73404192 : dpds += Ds * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + n * mesh->Nq + i];
1460 73404192 : dpdt += Dt * f[e * mesh->Np + n * mesh->Nq * mesh->Nq + j * mesh->Nq + i];
1461 : }
1462 :
1463 15160032 : const int id = k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1464 15160032 : grad_f[id + 0 * offset] = drdx * dpdr + dsdx * dpds + dtdx * dpdt;
1465 15160032 : grad_f[id + 1 * offset] = drdy * dpdr + dsdy * dpds + dtdy * dpdt;
1466 15160032 : grad_f[id + 2 * offset] = drdz * dpdr + dsdz * dpds + dtdz * dpdt;
1467 : }
1468 : }
1469 : }
1470 173652 : }
1471 :
1472 : bool
1473 181 : isHeatFluxBoundary(const int boundary)
1474 : {
1475 181 : auto sid = scalarDigitStr(nrs->scalar->nameToIndex.find("temperature")->second);
1476 181 : auto bcType = platform->app->bc->typeId(boundary, "scalar" + sid);
1477 :
1478 : // the purpose of this function is to check if the user has a non-zero flux boundary
1479 : // condition so that if MOOSE sends data to NekRS it will actually be used (even
1480 : // though zeroNeumann is technically a heat flux condition, it would not indicate
1481 : // the user has set up their model correctly)
1482 181 : return bcType == bdryBase::bcType_udfNeumann;
1483 : }
1484 :
1485 : bool
1486 1 : isMovingMeshBoundary(const int boundary)
1487 : {
1488 1 : auto bcType = platform->app->bc->typeId(boundary, "geom");
1489 1 : return bcType == bdryBase::bcType_udfDirichlet;
1490 : }
1491 :
1492 : bool
1493 0 : isTemperatureBoundary(const int boundary)
1494 : {
1495 0 : auto sid = scalarDigitStr(nrs->scalar->nameToIndex.find("temperature")->second);
1496 :
1497 0 : auto bcType = platform->app->bc->typeId(boundary, "scalar" + sid);
1498 0 : return bcType == bdryBase::bcType_udfDirichlet;
1499 : }
1500 :
1501 : const std::string
1502 1 : temperatureBoundaryType(const int boundary)
1503 : {
1504 1 : auto sid = scalarDigitStr(nrs->scalar->nameToIndex.find("temperature")->second);
1505 :
1506 1 : auto bcType = platform->app->bc->typeId(boundary, "scalar" + sid);
1507 2 : return platform->app->bc->sBcIDToText.at(bcType);
1508 : }
1509 :
1510 : int
1511 1507 : polynomialOrder()
1512 : {
1513 1507 : return entireMesh()->N;
1514 : }
1515 :
1516 : int
1517 760 : Nelements()
1518 : {
1519 760 : int n_local = entireMesh()->Nelements;
1520 : int n_global;
1521 760 : MPI_Allreduce(&n_local, &n_global, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm());
1522 760 : return n_global;
1523 : }
1524 :
1525 : int
1526 20803440 : Nfaces()
1527 : {
1528 20803440 : return entireMesh()->Nfaces;
1529 : }
1530 :
1531 : int
1532 761 : dim()
1533 : {
1534 761 : return entireMesh()->dim;
1535 : }
1536 :
1537 : int
1538 0 : NfaceVertices()
1539 : {
1540 0 : return entireMesh()->NfaceVertices;
1541 : }
1542 :
1543 : int
1544 760 : NboundaryFaces()
1545 : {
1546 760 : return entireMesh()->NboundaryFaces;
1547 : }
1548 :
1549 : int
1550 6968 : NboundaryID()
1551 : {
1552 6968 : if (hasCHT())
1553 182 : return nekData.NboundaryIDt;
1554 : else
1555 6786 : return nekData.NboundaryID;
1556 : }
1557 :
1558 : bool
1559 2359 : validBoundaryIDs(const std::vector<int> & boundary_id, int & first_invalid_id, int & n_boundaries)
1560 : {
1561 2359 : n_boundaries = NboundaryID();
1562 :
1563 : bool valid_boundary_ids = true;
1564 5489 : for (const auto & b : boundary_id)
1565 : {
1566 3130 : if ((b > n_boundaries) || (b <= 0))
1567 : {
1568 3 : first_invalid_id = b;
1569 : valid_boundary_ids = false;
1570 : }
1571 : }
1572 :
1573 2359 : return valid_boundary_ids;
1574 : }
1575 :
1576 : int
1577 6565572 : scalarSlot(const int id)
1578 : {
1579 6565572 : if (hasTemperatureVariable())
1580 : return id;
1581 : else
1582 6442332 : return id - 1;
1583 : }
1584 :
1585 : double
1586 6467424 : get_scalar01(const int id, const int surf_offset = 0)
1587 : {
1588 6467424 : return S[id + scalarSlot(1) * scalarFieldOffset()];
1589 : }
1590 :
1591 : double
1592 66720 : get_scalar02(const int id, const int surf_offset = 0)
1593 : {
1594 66720 : return S[id + scalarSlot(2) * scalarFieldOffset()];
1595 : }
1596 :
1597 : double
1598 31392 : get_scalar03(const int id, const int surf_offset = 0)
1599 : {
1600 31392 : return S[id + scalarSlot(3) * scalarFieldOffset()];
1601 : }
1602 :
1603 : double
1604 505472 : get_usrwrk00(const int id, const int surf_offset = 0)
1605 : {
1606 505472 : return usrwrk[id];
1607 : }
1608 :
1609 : double
1610 20672 : get_usrwrk01(const int id, const int surf_offset = 0)
1611 : {
1612 20672 : return usrwrk[id + nrs->fieldOffset];
1613 : }
1614 :
1615 : double
1616 20672 : get_usrwrk02(const int id, const int surf_offset = 0)
1617 : {
1618 20672 : return usrwrk[id + 2 * nrs->fieldOffset];
1619 : }
1620 :
1621 : double
1622 1115172216 : get_temperature(const int id, const int surf_offset)
1623 : {
1624 1115172216 : const auto sid = nrs->scalar->nameToIndex.find("temperature")->second;
1625 1115172216 : return S[id + sid * scalarFieldOffset()];
1626 : }
1627 :
1628 : double
1629 4032 : get_flux(const int id, const int surf_offset)
1630 : {
1631 : // TODO: this function does not support non-constant thermal conductivity
1632 : double k;
1633 4032 : platform->options.getArgs("SCALAR00 DIFFUSIONCOEFF", k);
1634 :
1635 : // this call of nek_mesh::all should be fine because flux is not a 'field' which can be
1636 : // provided to the postprocessors which have the option to operate only on part of the mesh
1637 4032 : auto mesh = getMesh(nek_mesh::all);
1638 4032 : int elem_id = id / mesh->Np;
1639 4032 : int vertex_id = id % mesh->Np;
1640 :
1641 4032 : const auto sid = nrs->scalar->nameToIndex.find("temperature")->second;
1642 4032 : const int offset = sid * scalarFieldOffset();
1643 : // This function is slightly inefficient, because we compute grad(T) for all nodes in
1644 : // an element even though we only call this function for one node at a time
1645 4032 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1646 4032 : gradient(mesh->Np, elem_id, S.data() + offset, grad_T, nek_mesh::all);
1647 :
1648 4032 : double normal_grad_T = grad_T[vertex_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1649 4032 : grad_T[vertex_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1650 4032 : grad_T[vertex_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1651 : freePointer(grad_T);
1652 :
1653 4032 : return -k * normal_grad_T;
1654 : }
1655 :
1656 : double
1657 96803984 : get_pressure(const int id, const int surf_offset)
1658 : {
1659 96803984 : return P[id];
1660 : }
1661 :
1662 : double
1663 37000524 : get_unity(const int /* id */, const int surf_offset)
1664 : {
1665 37000524 : return 1.0;
1666 : }
1667 :
1668 : double
1669 106814400 : get_velocity_x(const int id, const int surf_offset)
1670 : {
1671 106814400 : return U[id + 0 * nrs->fieldOffset];
1672 : }
1673 :
1674 : double
1675 44145600 : get_velocity_y(const int id, const int surf_offset)
1676 : {
1677 44145600 : return U[id + 1 * nrs->fieldOffset];
1678 : }
1679 :
1680 : double
1681 24067808 : get_velocity_z(const int id, const int surf_offset)
1682 : {
1683 24067808 : return U[id + 2 * nrs->fieldOffset];
1684 : }
1685 :
1686 : double
1687 2055492 : get_velocity(const int id, const int surf_offset)
1688 : {
1689 2055492 : int offset = nrs->fieldOffset;
1690 :
1691 2055492 : return std::sqrt(U[id + 0 * offset] * U[id + 0 * offset] +
1692 2055492 : U[id + 1 * offset] * U[id + 1 * offset] +
1693 2055492 : U[id + 2 * offset] * U[id + 2 * offset]);
1694 : }
1695 :
1696 : double
1697 47744 : get_velocity_x_squared(const int id, const int surf_offset)
1698 : {
1699 47744 : return std::pow(get_velocity_x(id, surf_offset), 2);
1700 : }
1701 :
1702 : double
1703 47744 : get_velocity_y_squared(const int id, const int surf_offset)
1704 : {
1705 47744 : return std::pow(get_velocity_y(id, surf_offset), 2);
1706 : }
1707 :
1708 : double
1709 52912 : get_velocity_z_squared(const int id, const int surf_offset)
1710 : {
1711 52912 : return std::pow(get_velocity_z(id, surf_offset), 2);
1712 : }
1713 :
1714 : void
1715 244 : checkFieldValidity(const field::NekWriteEnum & field)
1716 : {
1717 244 : switch (field)
1718 : {
1719 244 : case field::flux:
1720 244 : if (!hasTemperatureVariable())
1721 0 : mooseError("Cannot get NekRS heat flux "
1722 : "because your Nek case files do not have a temperature variable!");
1723 : break;
1724 0 : case field::heat_source:
1725 0 : if (!hasTemperatureVariable())
1726 0 : mooseError("Cannot get NekRS heat source "
1727 : "because your Nek case files do not have a temperature variable!");
1728 : break;
1729 : case field::x_displacement:
1730 : case field::y_displacement:
1731 : case field::z_displacement:
1732 : case field::mesh_velocity_x:
1733 : case field::mesh_velocity_y:
1734 : case field::mesh_velocity_z:
1735 : break;
1736 0 : default:
1737 0 : mooseError("Unhandled NekWriteEnum in checkFieldValidity!");
1738 : }
1739 244 : }
1740 :
1741 : void
1742 118903 : checkFieldValidity(const field::NekFieldEnum & field)
1743 : {
1744 : // by placing this check here, as opposed to inside the NekFieldInterface,
1745 : // we can also leverage this error checking for the 'outputs' of NekRSProblem,
1746 : // which does not inherit from NekFieldInterface but still accesses the solutionPointers.
1747 : // If this gets moved elsewhere, need to be sure to add dedicated testing for
1748 : // the 'outputs' on NekRSProblem.
1749 :
1750 : // TODO: would be nice for NekRSProblem to only access field information via the
1751 : // NekFieldInterface; refactor later
1752 :
1753 118903 : switch (field)
1754 : {
1755 60377 : case field::temperature:
1756 60377 : if (!hasTemperatureVariable())
1757 2 : mooseError("Cannot find 'temperature' "
1758 : "because your Nek case files do not have a temperature variable!");
1759 : break;
1760 137 : case field::scalar01:
1761 137 : if (!hasScalarVariable(1))
1762 1 : mooseError("Cannot find 'scalar01' "
1763 : "because your Nek case files do not have a scalar01 variable!");
1764 : break;
1765 213 : case field::scalar02:
1766 213 : if (!hasScalarVariable(2))
1767 1 : mooseError("Cannot find 'scalar02' "
1768 : "because your Nek case files do not have a scalar02 variable!");
1769 : break;
1770 69 : case field::scalar03:
1771 69 : if (!hasScalarVariable(3))
1772 1 : mooseError("Cannot find 'scalar03' "
1773 : "because your Nek case files do not have a scalar03 variable!");
1774 : break;
1775 454 : case field::usrwrk00:
1776 454 : if (n_usrwrk_slots < 1)
1777 2 : mooseError("Cannot find 'usrwrk00' because you have only allocated 'n_usrwrk_slots = " +
1778 1 : std::to_string(n_usrwrk_slots) + "'");
1779 : break;
1780 42 : case field::usrwrk01:
1781 42 : if (n_usrwrk_slots < 2)
1782 2 : mooseError("Cannot find 'usrwrk01' because you have only allocated 'n_usrwrk_slots = " +
1783 1 : std::to_string(n_usrwrk_slots) + "'");
1784 : break;
1785 42 : case field::usrwrk02:
1786 42 : if (n_usrwrk_slots < 3)
1787 2 : mooseError("Cannot find 'usrwrk02' because you have only allocated 'n_usrwrk_slots = " +
1788 1 : std::to_string(n_usrwrk_slots) + "'");
1789 : break;
1790 : }
1791 118895 : }
1792 :
1793 244 : double (*solutionPointer(const field::NekWriteEnum & field))(int, int)
1794 : {
1795 : double (*f)(int, int);
1796 :
1797 244 : checkFieldValidity(field);
1798 :
1799 244 : switch (field)
1800 : {
1801 244 : case field::flux:
1802 : f = &get_flux;
1803 : break;
1804 0 : default:
1805 0 : mooseError("Unhandled NekWriteEnum in solutionPointer!");
1806 : }
1807 :
1808 244 : return f;
1809 : }
1810 :
1811 115858 : double (*solutionPointer(const field::NekFieldEnum & field))(int, int)
1812 : {
1813 : // we include this here as well, in addition to within the NekFieldInterface, because
1814 : // the NekRSProblem accesses these methods without inheriting from NekFieldInterface
1815 115858 : checkFieldValidity(field);
1816 :
1817 : double (*f)(int, int);
1818 :
1819 115858 : switch (field)
1820 : {
1821 : case field::velocity_x:
1822 : f = &get_velocity_x;
1823 : break;
1824 6776 : case field::velocity_y:
1825 : f = &get_velocity_y;
1826 6776 : break;
1827 3784 : case field::velocity_z:
1828 : f = &get_velocity_z;
1829 3784 : break;
1830 204 : case field::velocity:
1831 : f = &get_velocity;
1832 204 : break;
1833 0 : case field::velocity_component:
1834 0 : mooseError("The 'velocity_component' field is not compatible with the solutionPointer "
1835 : "interface!");
1836 : break;
1837 48 : case field::velocity_x_squared:
1838 : f = &get_velocity_x_squared;
1839 48 : break;
1840 48 : case field::velocity_y_squared:
1841 : f = &get_velocity_y_squared;
1842 48 : break;
1843 52 : case field::velocity_z_squared:
1844 : f = &get_velocity_z_squared;
1845 52 : break;
1846 59330 : case field::temperature:
1847 : f = &get_temperature;
1848 59330 : break;
1849 18800 : case field::pressure:
1850 : f = &get_pressure;
1851 18800 : break;
1852 76 : case field::scalar01:
1853 : f = &get_scalar01;
1854 76 : break;
1855 180 : case field::scalar02:
1856 : f = &get_scalar02;
1857 180 : break;
1858 36 : case field::scalar03:
1859 : f = &get_scalar03;
1860 36 : break;
1861 5796 : case field::unity:
1862 : f = &get_unity;
1863 5796 : break;
1864 420 : case field::usrwrk00:
1865 : f = &get_usrwrk00;
1866 420 : break;
1867 16 : case field::usrwrk01:
1868 : f = &get_usrwrk01;
1869 16 : break;
1870 16 : case field::usrwrk02:
1871 : f = &get_usrwrk02;
1872 16 : break;
1873 0 : default:
1874 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1875 : }
1876 :
1877 115858 : return f;
1878 : }
1879 :
1880 : void
1881 129 : initializeDimensionalScales(const double U,
1882 : const double T,
1883 : const double dT,
1884 : const double L,
1885 : const double rho,
1886 : const double Cp,
1887 : const double s01,
1888 : const double ds01,
1889 : const double s02,
1890 : const double ds02,
1891 : const double s03,
1892 : const double ds03)
1893 : {
1894 129 : scales.U_ref = U;
1895 129 : scales.T_ref = T;
1896 129 : scales.dT_ref = dT;
1897 129 : scales.L_ref = L;
1898 129 : scales.A_ref = L * L;
1899 129 : scales.V_ref = L * L * L;
1900 129 : scales.rho_ref = rho;
1901 129 : scales.Cp_ref = Cp;
1902 129 : scales.t_ref = L / U;
1903 129 : scales.P_ref = rho * U * U;
1904 :
1905 129 : scales.s01_ref = s01;
1906 129 : scales.ds01_ref = ds01;
1907 129 : scales.s02_ref = s02;
1908 129 : scales.ds02_ref = ds02;
1909 129 : scales.s03_ref = s03;
1910 129 : scales.ds03_ref = ds03;
1911 :
1912 129 : scales.flux_ref = rho * U * Cp * dT;
1913 129 : scales.source_ref = scales.flux_ref / L;
1914 129 : }
1915 :
1916 : double
1917 469 : referenceLength()
1918 : {
1919 469 : return scales.L_ref;
1920 : }
1921 :
1922 : double
1923 163562 : referenceTime()
1924 : {
1925 163562 : return scales.t_ref;
1926 : }
1927 :
1928 : double
1929 221 : referenceArea()
1930 : {
1931 221 : return scales.A_ref;
1932 : }
1933 :
1934 : double
1935 878 : referenceVolume()
1936 : {
1937 878 : return scales.V_ref;
1938 : }
1939 :
1940 : Real
1941 28979014 : nondimensionalAdditive(const field::NekFieldEnum & field)
1942 : {
1943 28979014 : switch (field)
1944 : {
1945 17394940 : case field::temperature:
1946 17394940 : return scales.T_ref;
1947 4496684 : case field::scalar01:
1948 4496684 : return scales.s01_ref;
1949 40912 : case field::scalar02:
1950 40912 : return scales.s02_ref;
1951 5584 : case field::scalar03:
1952 5584 : return scales.s03_ref;
1953 : default:
1954 : return 0;
1955 : }
1956 : }
1957 :
1958 : Real
1959 17306 : nondimensionalAdditive(const field::NekWriteEnum & field)
1960 : {
1961 17306 : switch (field)
1962 : {
1963 17306 : case field::flux:
1964 : case field::heat_source:
1965 : case field::x_displacement:
1966 : case field::y_displacement:
1967 : case field::z_displacement:
1968 : case field::mesh_velocity_x:
1969 : case field::mesh_velocity_y:
1970 : case field::mesh_velocity_z:
1971 17306 : return 0.0;
1972 0 : default:
1973 0 : mooseError("Unhandled NekWriteEnum in nondimensionalAdditive!");
1974 : }
1975 : }
1976 :
1977 : Real
1978 43795 : nondimensionalDivisor(const field::NekWriteEnum & field)
1979 : {
1980 43795 : switch (field)
1981 : {
1982 40963 : case field::flux:
1983 40963 : return scales.flux_ref;
1984 2827 : case field::heat_source:
1985 2827 : return scales.source_ref;
1986 0 : case field::x_displacement:
1987 : case field::y_displacement:
1988 : case field::z_displacement:
1989 0 : return scales.L_ref;
1990 5 : case field::mesh_velocity_x:
1991 : case field::mesh_velocity_y:
1992 : case field::mesh_velocity_z:
1993 5 : return scales.U_ref;
1994 0 : default:
1995 0 : mooseError("Unhandled NekWriteEnum in nondimensionalDivisor!");
1996 : }
1997 : }
1998 :
1999 : Real
2000 28979136 : nondimensionalDivisor(const field::NekFieldEnum & field)
2001 : {
2002 28979136 : switch (field)
2003 : {
2004 5354476 : case field::velocity_x:
2005 : case field::velocity_y:
2006 : case field::velocity_z:
2007 : case field::velocity:
2008 : case field::velocity_component:
2009 5354476 : return scales.U_ref;
2010 5336 : case field::velocity_x_squared:
2011 : case field::velocity_y_squared:
2012 : case field::velocity_z_squared:
2013 5336 : return scales.U_ref * scales.U_ref;
2014 17394940 : case field::temperature:
2015 17394940 : return scales.dT_ref;
2016 1660587 : case field::pressure:
2017 1660587 : return scales.P_ref;
2018 4496684 : case field::scalar01:
2019 4496684 : return scales.ds01_ref;
2020 40912 : case field::scalar02:
2021 40912 : return scales.ds02_ref;
2022 5584 : case field::scalar03:
2023 5584 : return scales.ds03_ref;
2024 : case field::unity:
2025 : // no dimensionalization needed
2026 : return 1.0;
2027 445 : case field::usrwrk00:
2028 445 : return scratchUnits(0);
2029 25 : case field::usrwrk01:
2030 25 : return scratchUnits(1);
2031 25 : case field::usrwrk02:
2032 25 : return scratchUnits(2);
2033 0 : default:
2034 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
2035 : }
2036 : }
2037 :
2038 : Real
2039 495 : scratchUnits(const int slot)
2040 : {
2041 : // if (indices.flux != -1 && slot == indices.flux / nekrs::fieldOffset())
2042 : // return scales.flux_ref;
2043 : // else if (indices.heat_source != -1 && slot == indices.heat_source / nekrs::fieldOffset())
2044 : // return scales.source_ref;
2045 495 : if (is_nondimensional)
2046 : {
2047 66 : mooseDoOnce(mooseWarning(
2048 : "The units of 'usrwrk0" + std::to_string(slot) +
2049 : "' are unknown, so we cannot dimensionalize any objects using 'field = usrwrk0" +
2050 : std::to_string(slot) +
2051 : "'. The output for this quantity will be given in non-dimensional form.\n\nYou will need "
2052 : "to manipulate the data manually from Cardinal if you need to dimensionalize it."));
2053 : }
2054 :
2055 492 : return 1.0;
2056 : }
2057 :
2058 : void
2059 749 : nondimensional(const bool n)
2060 : {
2061 749 : is_nondimensional = n;
2062 749 : }
2063 :
2064 : template <>
2065 : MPI_Datatype
2066 66254 : resolveType<double>()
2067 : {
2068 66254 : return MPI_DOUBLE;
2069 : }
2070 :
2071 : template <>
2072 : MPI_Datatype
2073 6504 : resolveType<int>()
2074 : {
2075 6504 : return MPI_INT;
2076 : }
2077 :
2078 : int
2079 164481042 : Nscalar()
2080 : {
2081 : int Nscalar;
2082 164481042 : platform->options.getArgs("NUMBER OF SCALARS", Nscalar);
2083 164481042 : return Nscalar;
2084 : }
2085 :
2086 : void
2087 764 : initializeNekHostArrays()
2088 : {
2089 764 : nrs = dynamic_cast<nrs_t *>(platform->app);
2090 :
2091 764 : mesh_t * mesh = entireMesh();
2092 :
2093 764 : U.resize(mesh->dim * nrs->fluid->fieldOffset);
2094 764 : P.resize(mesh->Nlocal);
2095 :
2096 764 : if (Nscalar())
2097 633 : S.resize(nrs->scalar->NSfields * nrs->scalar->fieldOffset()); // offset is same for all scalars
2098 764 : }
2099 :
2100 : dfloat *
2101 1 : host_x()
2102 : {
2103 1 : return x;
2104 : }
2105 : dfloat *
2106 1 : host_y()
2107 : {
2108 1 : return y;
2109 : }
2110 : dfloat *
2111 1 : host_z()
2112 : {
2113 1 : return z;
2114 : }
2115 :
2116 : // std::tuple<std::vector<dfloat>&, std::vector<dfloat>&, std::vector<dfloat>&>
2117 : std::tuple<dfloat *, dfloat *, dfloat *>
2118 1 : host_xyz()
2119 : {
2120 1 : return {x, y, z};
2121 : }
2122 :
2123 : std::vector<dfloat> &
2124 0 : host_U()
2125 : {
2126 0 : return U;
2127 : }
2128 :
2129 : std::vector<dfloat> &
2130 0 : host_P()
2131 : {
2132 0 : return P;
2133 : }
2134 :
2135 : std::vector<dfloat> &
2136 0 : host_S()
2137 : {
2138 0 : return S;
2139 : }
2140 :
2141 : dfloat *
2142 1199772 : host_wrk()
2143 : {
2144 1199772 : return usrwrk;
2145 : }
2146 :
2147 : nrs_t *
2148 1230336 : nrsPtr()
2149 : {
2150 1230336 : return nrs;
2151 : }
2152 :
2153 : mesh_t *
2154 54 : createMesh2(mesh_t * _mesh, int Nc)
2155 : {
2156 54 : mesh_t * mesh = new mesh_t();
2157 54 : memcpy(mesh, _mesh, sizeof(mesh_t));
2158 :
2159 : const int cubN = 0;
2160 54 : meshLoadReferenceNodesHex3D(mesh, Nc, cubN);
2161 :
2162 54 : mesh->geometricFactorsKernel = nullptr;
2163 54 : mesh->surfaceGeometricFactorsKernel = nullptr;
2164 54 : mesh->cubatureGeometricFactorsKernel = nullptr;
2165 :
2166 54 : mesh->o_D = platform->device.malloc<dfloat>(mesh->Nq * mesh->Nq, mesh->D);
2167 :
2168 54 : dfloat * DT = (dfloat *)calloc(mesh->Nq * mesh->Nq, sizeof(dfloat));
2169 216 : for (int j = 0; j < mesh->Nq; j++)
2170 : {
2171 648 : for (int i = 0; i < mesh->Nq; i++)
2172 : {
2173 486 : DT[j * mesh->Nq + i] = mesh->D[i * mesh->Nq + j];
2174 : }
2175 : }
2176 54 : mesh->o_DT = platform->device.malloc<dfloat>(mesh->Nq * mesh->Nq, DT);
2177 54 : free(DT);
2178 :
2179 54 : meshPhysicalNodesHex3D(mesh);
2180 :
2181 54 : meshConnectFaceNodes3D(mesh);
2182 :
2183 54 : meshGlobalIds(mesh);
2184 :
2185 54 : meshParallelGatherScatterSetup(
2186 54 : mesh, mesh->Nlocal, mesh->globalIds, platform->comm.mpiComm(), OOGS_AUTO, 0);
2187 :
2188 54 : return mesh;
2189 : }
2190 :
2191 : } // end namespace nekrs
2192 :
2193 : #endif
|
