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 dfloat * sgeo;
26 : static dfloat * vgeo;
27 : static unsigned int n_usrwrk_slots;
28 : static bool is_nondimensional;
29 :
30 : namespace nekrs
31 : {
32 : static double setup_time;
33 :
34 : // various constants for controlling tolerances
35 : static double abs_tol;
36 : static double rel_tol;
37 :
38 : void
39 339 : setAbsoluteTol(double tol)
40 : {
41 339 : abs_tol = tol;
42 339 : }
43 :
44 : void
45 339 : setRelativeTol(double tol)
46 : {
47 339 : rel_tol = tol;
48 339 : }
49 :
50 : void
51 811 : setNekSetupTime(const double & time)
52 : {
53 811 : setup_time = time;
54 811 : }
55 :
56 : double
57 800 : getNekSetupTime()
58 : {
59 800 : return setup_time;
60 : }
61 :
62 : void
63 740 : setStartTime(const double & start)
64 : {
65 1480 : platform->options.setArgs("START TIME", to_string_f(start));
66 740 : }
67 :
68 : void
69 62 : write_usrwrk_field_file(const int & slot, const std::string & prefix, const dfloat & time, const int & step, const bool & write_coords)
70 : {
71 62 : int num_bytes = fieldOffset() * sizeof(dfloat);
72 :
73 62 : nrs_t * nrs = (nrs_t *)nrsPtr();
74 62 : occa::memory o_write = platform->device.malloc(num_bytes);
75 124 : o_write.copyFrom(nrs->o_usrwrk, num_bytes /* length we are copying */,
76 62 : 0 /* where to place data */, num_bytes * slot /* where to source data */);
77 :
78 62 : occa::memory o_null;
79 62 : writeFld(prefix.c_str(), time, step, write_coords, 1 /* FP64 */, o_null, o_null, o_write, 1);
80 62 : }
81 :
82 : void
83 56 : write_field_file(const std::string & prefix, const dfloat time, const int & step)
84 : {
85 56 : nrs_t * nrs = (nrs_t *)nrsPtr();
86 :
87 : int Nscalar = 0;
88 56 : occa::memory o_s;
89 56 : if (nrs->Nscalar)
90 : {
91 56 : o_s = nrs->cds->o_S;
92 56 : Nscalar = nrs->Nscalar;
93 : }
94 :
95 56 : writeFld(
96 56 : prefix.c_str(), time, step, 1 /* coords */, 1 /* FP64 */, nrs->o_U, nrs->o_P, o_s, Nscalar);
97 56 : }
98 :
99 : void
100 813 : buildOnly(int buildOnly)
101 : {
102 813 : build_only = buildOnly;
103 813 : }
104 :
105 : int
106 105093 : buildOnly()
107 : {
108 105093 : return build_only;
109 : }
110 :
111 : bool
112 0 : hasCHT()
113 : {
114 0 : return entireMesh()->cht;
115 : }
116 :
117 : bool
118 36604 : hasMovingMesh()
119 : {
120 36604 : return platform->options.compareArgs("MOVING MESH", "TRUE");
121 : }
122 :
123 : bool
124 33218 : hasVariableDt()
125 : {
126 33218 : return platform->options.compareArgs("VARIABLE DT", "TRUE");
127 : }
128 :
129 : bool
130 953 : hasBlendingSolver()
131 : {
132 1906 : return !platform->options.compareArgs("MESH SOLVER", "NONE") && hasMovingMesh();
133 : }
134 :
135 : bool
136 954 : hasUserMeshSolver()
137 : {
138 1908 : return platform->options.compareArgs("MESH SOLVER", "NONE") && hasMovingMesh();
139 : }
140 :
141 : bool
142 1522 : endControlElapsedTime()
143 : {
144 3044 : return !platform->options.getArgs("STOP AT ELAPSED TIME").empty();
145 : }
146 :
147 : bool
148 1522 : endControlTime()
149 : {
150 1522 : return endTime() > 0;
151 : }
152 :
153 : bool
154 761 : endControlNumSteps()
155 : {
156 761 : return !endControlElapsedTime() && !endControlTime();
157 : }
158 :
159 : bool
160 148986037 : hasTemperatureVariable()
161 : {
162 148986037 : nrs_t * nrs = (nrs_t *)nrsPtr();
163 446485487 : return nrs->Nscalar ? platform->options.compareArgs("SCALAR00 IS TEMPERATURE", "TRUE") : false;
164 : }
165 :
166 : bool
167 627 : hasTemperatureSolve()
168 : {
169 627 : nrs_t * nrs = (nrs_t *)nrsPtr();
170 627 : return hasTemperatureVariable() ? nrs->cds->compute[0] : false;
171 : }
172 :
173 : bool
174 1207 : hasScalarVariable(int scalarId)
175 : {
176 1207 : nrs_t * nrs = (nrs_t *)nrsPtr();
177 1207 : return (scalarId < nrs->Nscalar);
178 : }
179 :
180 : bool
181 93 : hasHeatSourceKernel()
182 : {
183 93 : return static_cast<bool>(udf.sEqnSource);
184 : }
185 :
186 : bool
187 809 : isInitialized()
188 : {
189 809 : nrs_t * nrs = (nrs_t *)nrsPtr();
190 809 : return nrs;
191 : }
192 :
193 : int
194 9125063 : scalarFieldOffset()
195 : {
196 9125063 : nrs_t * nrs = (nrs_t *)nrsPtr();
197 9125063 : return nrs->cds->fieldOffset[0];
198 : }
199 :
200 : int
201 1831210 : velocityFieldOffset()
202 : {
203 1831210 : nrs_t * nrs = (nrs_t *)nrsPtr();
204 1831210 : return nrs->fieldOffset;
205 : }
206 :
207 : int
208 9001831 : fieldOffset()
209 : {
210 9001831 : if (hasTemperatureVariable())
211 9001815 : return scalarFieldOffset();
212 : else
213 16 : return velocityFieldOffset();
214 : }
215 :
216 : mesh_t *
217 139898829 : entireMesh()
218 : {
219 139898829 : if (hasTemperatureVariable())
220 139662537 : return temperatureMesh();
221 : else
222 236292 : return flowMesh();
223 : }
224 :
225 : mesh_t *
226 584145 : flowMesh()
227 : {
228 584145 : nrs_t * nrs = (nrs_t *)nrsPtr();
229 584145 : return nrs->meshV;
230 : }
231 :
232 : mesh_t *
233 140467810 : temperatureMesh()
234 : {
235 140467810 : nrs_t * nrs = (nrs_t *)nrsPtr();
236 140467810 : return nrs->cds->mesh[0];
237 : }
238 :
239 : mesh_t *
240 396170 : getMesh(const nek_mesh::NekMeshEnum pp_mesh)
241 : {
242 396170 : if (pp_mesh == nek_mesh::fluid)
243 2188 : return flowMesh();
244 393982 : else if (pp_mesh == nek_mesh::all)
245 393981 : return entireMesh();
246 : else
247 1 : mooseError("This object does not support operations on the solid part of the NekRS mesh!\n"
248 : "Valid options for 'mesh' are 'fluid' or 'all'.");
249 : }
250 :
251 : int
252 12765306 : commRank()
253 : {
254 12765306 : return platform->comm.mpiRank;
255 : }
256 :
257 : int
258 1048657 : commSize()
259 : {
260 1048657 : return platform->comm.mpiCommSize;
261 : }
262 :
263 : bool
264 811 : scratchAvailable()
265 : {
266 811 : nrs_t * nrs = (nrs_t *)nrsPtr();
267 :
268 : // Because these scratch spaces are available for whatever the user sees fit, it is
269 : // possible that the user wants to use these arrays for a _different_ purpose aside from
270 : // transferring in MOOSE values. In nekrs::setup, we call the UDF_Setup0, UDF_Setup,
271 : // and UDF_ExecuteStep routines. These scratch space arrays aren't initialized anywhere
272 : // else in the core base, so we will make sure to throw an error from MOOSE if these
273 : // arrays are already in use, because otherwise our MOOSE transfer might get overwritten
274 : // by whatever other operation the user is trying to do.
275 811 : if (nrs->usrwrk)
276 1 : return false;
277 :
278 : return true;
279 : }
280 :
281 : void
282 810 : initializeScratch(const unsigned int & n_slots)
283 : {
284 810 : if (n_slots == 0)
285 : return;
286 :
287 413 : nrs_t * nrs = (nrs_t *)nrsPtr();
288 413 : mesh_t * mesh = entireMesh();
289 :
290 : // clear them just to be sure
291 413 : freeScratch();
292 :
293 : // In order to make indexing simpler in the device user functions (which is where the
294 : // boundary conditions are then actually applied), we define these scratch arrays
295 : // as volume arrays.
296 413 : nrs->usrwrk = (double *)calloc(n_slots * fieldOffset(), sizeof(double));
297 413 : nrs->o_usrwrk = platform->device.malloc(n_slots * fieldOffset() * sizeof(double), nrs->usrwrk);
298 :
299 413 : n_usrwrk_slots = n_slots;
300 : }
301 :
302 : void
303 1145 : freeScratch()
304 : {
305 1145 : nrs_t * nrs = (nrs_t *)nrsPtr();
306 :
307 1145 : freePointer(nrs->usrwrk);
308 1145 : nrs->o_usrwrk.free();
309 1145 : }
310 :
311 : double
312 276 : viscosity()
313 : {
314 : dfloat mu;
315 276 : setupAide & options = platform->options;
316 276 : options.getArgs("VISCOSITY", mu);
317 :
318 : // because we set rho_ref, U_ref, and L_ref all equal to 1 if our input is dimensional,
319 : // we don't need to have separate treatments for dimensional vs. nondimensional cases
320 276 : dfloat Re = 1.0 / mu;
321 276 : return scales.rho_ref * scales.U_ref * scales.L_ref / Re;
322 : }
323 :
324 : double
325 92 : Pr()
326 : {
327 : dfloat rho, rho_cp, k;
328 92 : setupAide & options = platform->options;
329 92 : options.getArgs("DENSITY", rho);
330 92 : options.getArgs("SCALAR00 DENSITY", rho_cp);
331 92 : options.getArgs("SCALAR00 DIFFUSIVITY", k);
332 :
333 92 : dfloat Pe = 1.0 / k;
334 92 : dfloat conductivity = scales.rho_ref * scales.U_ref * scales.Cp_ref * scales.L_ref / Pe;
335 92 : dfloat Cp = rho_cp / rho * scales.Cp_ref;
336 92 : return viscosity() * Cp / conductivity;
337 : }
338 :
339 : void
340 1596 : interpolationMatrix(double * I, int starting_points, int ending_points)
341 : {
342 1596 : DegreeRaiseMatrix1D(starting_points - 1, ending_points - 1, I);
343 1596 : }
344 :
345 : void
346 2198994 : interpolateVolumeHex3D(const double * I, double * x, int N, double * Ix, int M)
347 : {
348 2198994 : double * Ix1 = (dfloat *)calloc(N * N * M, sizeof(double));
349 2198994 : double * Ix2 = (dfloat *)calloc(N * M * M, sizeof(double));
350 :
351 10888880 : for (int k = 0; k < N; ++k)
352 51895756 : for (int j = 0; j < N; ++j)
353 180941182 : for (int i = 0; i < M; ++i)
354 : {
355 : dfloat tmp = 0;
356 939543840 : for (int n = 0; n < N; ++n)
357 801808528 : tmp += I[i * N + n] * x[k * N * N + j * N + n];
358 137735312 : Ix1[k * N * M + j * M + i] = tmp;
359 : }
360 :
361 10888880 : for (int k = 0; k < N; ++k)
362 38012670 : for (int j = 0; j < M; ++j)
363 141667368 : for (int i = 0; i < M; ++i)
364 : {
365 : dfloat tmp = 0;
366 585034584 : for (int n = 0; n < N; ++n)
367 472690000 : tmp += I[j * N + n] * Ix1[k * N * M + n * M + i];
368 112344584 : Ix2[k * M * M + j * M + i] = tmp;
369 : }
370 :
371 10153964 : for (int k = 0; k < M; ++k)
372 42219396 : for (int j = 0; j < M; ++j)
373 211512576 : for (int i = 0; i < M; ++i)
374 : {
375 : dfloat tmp = 0;
376 688209246 : for (int n = 0; n < N; ++n)
377 510961096 : tmp += I[k * N + n] * Ix2[n * M * M + j * M + i];
378 177248150 : Ix[k * M * M + j * M + i] = tmp;
379 : }
380 :
381 : freePointer(Ix1);
382 : freePointer(Ix2);
383 2198994 : }
384 :
385 : void
386 478350 : interpolateSurfaceFaceHex3D(
387 : double * scratch, const double * I, double * x, int N, double * Ix, int M)
388 : {
389 1908266 : for (int j = 0; j < N; ++j)
390 6512948 : for (int i = 0; i < M; ++i)
391 : {
392 : double tmp = 0;
393 23882184 : for (int n = 0; n < N; ++n)
394 : {
395 18799152 : tmp += I[i * N + n] * x[j * N + n];
396 : }
397 5083032 : scratch[j * M + i] = tmp;
398 : }
399 :
400 2310042 : for (int j = 0; j < M; ++j)
401 9255804 : for (int i = 0; i < M; ++i)
402 : {
403 : double tmp = 0;
404 26531280 : for (int n = 0; n < N; ++n)
405 : {
406 19107168 : tmp += I[j * N + n] * scratch[n * M + i];
407 : }
408 7424112 : Ix[j * M + i] = tmp;
409 : }
410 478350 : }
411 :
412 : void
413 93511 : displacementAndCounts(const std::vector<int> & base_counts,
414 : int * counts,
415 : int * displacement,
416 : const int multiplier = 1.0)
417 : {
418 476702 : for (int i = 0; i < commSize(); ++i)
419 383191 : counts[i] = base_counts[i] * multiplier;
420 :
421 93511 : displacement[0] = 0;
422 383191 : for (int i = 1; i < commSize(); i++)
423 289680 : displacement[i] = displacement[i - 1] + counts[i - 1];
424 93511 : }
425 :
426 : double
427 1940 : usrwrkVolumeIntegral(const unsigned int & slot, const nek_mesh::NekMeshEnum pp_mesh)
428 : {
429 1940 : nrs_t * nrs = (nrs_t *)nrsPtr();
430 1940 : const auto & mesh = getMesh(pp_mesh);
431 :
432 1940 : double integral = 0.0;
433 :
434 763608 : for (int k = 0; k < mesh->Nelements; ++k)
435 : {
436 761668 : int offset = k * mesh->Np;
437 :
438 199804868 : for (int v = 0; v < mesh->Np; ++v)
439 199043200 : integral += nrs->usrwrk[slot + offset + v] * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
440 : }
441 :
442 : // sum across all processes
443 : double total_integral;
444 1940 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
445 :
446 1940 : return total_integral;
447 : }
448 :
449 : void
450 958 : scaleUsrwrk(const unsigned int & slot, const dfloat & value)
451 : {
452 958 : nrs_t * nrs = (nrs_t *)nrsPtr();
453 958 : mesh_t * mesh = getMesh(nek_mesh::all);
454 :
455 381664 : for (int k = 0; k < mesh->Nelements; ++k)
456 : {
457 380706 : int id = k * mesh->Np;
458 :
459 99898850 : for (int v = 0; v < mesh->Np; ++v)
460 99518144 : nrs->usrwrk[slot + id + v] *= value;
461 : }
462 958 : }
463 :
464 : std::vector<double>
465 25359 : usrwrkSideIntegral(const unsigned int & slot,
466 : const std::vector<int> & boundary,
467 : const nek_mesh::NekMeshEnum pp_mesh)
468 : {
469 25359 : nrs_t * nrs = (nrs_t *)nrsPtr();
470 25359 : const auto & mesh = getMesh(pp_mesh);
471 :
472 25359 : std::vector<double> integral(boundary.size(), 0.0);
473 :
474 7860674 : for (int i = 0; i < mesh->Nelements; ++i)
475 : {
476 54847205 : for (int j = 0; j < mesh->Nfaces; ++j)
477 : {
478 47011890 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
479 :
480 47011890 : if (std::find(boundary.begin(), boundary.end(), face_id) != boundary.end())
481 : {
482 814188 : auto it = std::find(boundary.begin(), boundary.end(), face_id);
483 : auto b_index = it - boundary.begin();
484 :
485 814188 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
486 :
487 15160204 : for (int v = 0; v < mesh->Nfp; ++v)
488 14346016 : integral[b_index] += nrs->usrwrk[slot + mesh->vmapM[offset + v]] *
489 14346016 : sgeo[mesh->Nsgeo * (offset + v) + WSJID];
490 : }
491 : }
492 : }
493 :
494 : // sum across all processes; this can probably be done more efficiently
495 25359 : std::vector<double> total_integral(boundary.size(), 0.0);
496 50839 : for (std::size_t i = 0; i < boundary.size(); ++i)
497 25480 : MPI_Allreduce(&integral[i], &total_integral[i], 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
498 :
499 25359 : return total_integral;
500 25359 : }
501 :
502 : void
503 0 : limitTemperature(const double * min_T, const double * max_T)
504 : {
505 : // if no limiters are provided, simply return
506 0 : if (!min_T && !max_T)
507 : return;
508 :
509 0 : double minimum = min_T ? *min_T : std::numeric_limits<double>::min();
510 0 : double maximum = max_T ? *max_T : std::numeric_limits<double>::max();
511 :
512 : // nondimensionalize if necessary
513 0 : minimum = (minimum - scales.T_ref) / scales.dT_ref;
514 0 : maximum = (maximum - scales.T_ref) / scales.dT_ref;
515 :
516 0 : nrs_t * nrs = (nrs_t *)nrsPtr();
517 0 : mesh_t * mesh = temperatureMesh();
518 :
519 0 : for (int i = 0; i < mesh->Nelements; ++i)
520 : {
521 0 : for (int j = 0; j < mesh->Np; ++j)
522 : {
523 0 : int id = i * mesh->Np + j;
524 :
525 0 : if (nrs->cds->S[id] < minimum)
526 0 : nrs->cds->S[id] = minimum;
527 0 : if (nrs->cds->S[id] > maximum)
528 0 : nrs->cds->S[id] = maximum;
529 : }
530 : }
531 :
532 : // when complete, copy to device
533 0 : nrs->cds->o_S.copyFrom(nrs->cds->S);
534 : }
535 :
536 : void
537 1504 : copyDeformationToDevice()
538 : {
539 1504 : mesh_t * mesh = entireMesh();
540 1504 : mesh->o_x.copyFrom(mesh->x);
541 1504 : mesh->o_y.copyFrom(mesh->y);
542 1504 : mesh->o_z.copyFrom(mesh->z);
543 1504 : mesh->update();
544 :
545 1504 : updateHostMeshParameters();
546 1504 : }
547 :
548 : void
549 808 : initializeHostMeshParameters()
550 : {
551 808 : mesh_t * mesh = entireMesh();
552 808 : sgeo = (dfloat *)calloc(mesh->o_sgeo.size(), sizeof(dfloat));
553 808 : vgeo = (dfloat *)calloc(mesh->o_vgeo.size(), sizeof(dfloat));
554 808 : }
555 :
556 : void
557 2312 : updateHostMeshParameters()
558 : {
559 2312 : mesh_t * mesh = entireMesh();
560 2312 : mesh->o_sgeo.copyTo(sgeo);
561 2312 : mesh->o_vgeo.copyTo(vgeo);
562 2312 : }
563 :
564 : dfloat *
565 104 : getSgeo()
566 : {
567 104 : return sgeo;
568 : }
569 :
570 : dfloat *
571 355 : getVgeo()
572 : {
573 355 : return vgeo;
574 : }
575 :
576 : double
577 24332 : sideExtremeValue(const std::vector<int> & boundary_id, const field::NekFieldEnum & field,
578 : const nek_mesh::NekMeshEnum pp_mesh, const bool max)
579 : {
580 24332 : mesh_t * mesh = getMesh(pp_mesh);
581 :
582 24332 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
583 :
584 : double (*f)(int, int);
585 24332 : f = solutionPointer(field);
586 :
587 2694568 : for (int i = 0; i < mesh->Nelements; ++i)
588 : {
589 18691652 : for (int j = 0; j < mesh->Nfaces; ++j)
590 : {
591 16021416 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
592 :
593 16021416 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
594 : {
595 270524 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
596 11791388 : for (int v = 0; v < mesh->Nfp; ++v)
597 : {
598 11520864 : if (max)
599 5835328 : value = std::max(value, f(mesh->vmapM[offset + v], 0 /* unused */));
600 : else
601 5809910 : value = std::min(value, f(mesh->vmapM[offset + v], 0 /* unused */));
602 : }
603 : }
604 : }
605 : }
606 :
607 : // find extreme value across all processes
608 : double reduced_value;
609 24332 : auto op = max ? MPI_MAX : MPI_MIN;
610 24332 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm);
611 :
612 : // dimensionalize the field if needed
613 24332 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
614 :
615 24332 : return reduced_value;
616 : }
617 :
618 : double
619 46788 : volumeExtremeValue(const field::NekFieldEnum & field, const nek_mesh::NekMeshEnum pp_mesh, const bool max)
620 : {
621 46788 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
622 :
623 : double (*f)(int, int);
624 46788 : f = solutionPointer(field);
625 :
626 : mesh_t * mesh;
627 : int start_id;
628 :
629 46788 : switch (pp_mesh)
630 : {
631 46776 : case nek_mesh::fluid:
632 : case nek_mesh::all:
633 : {
634 46776 : mesh = getMesh(pp_mesh);
635 : start_id = 0;
636 : break;
637 : }
638 12 : case nek_mesh::solid:
639 : {
640 12 : mesh = entireMesh();
641 12 : start_id = flowMesh()->Nelements;
642 12 : break;
643 : }
644 0 : default:
645 0 : mooseError("Unhandled NekMeshEnum in volumeExtremeValue");
646 : }
647 :
648 7914462 : for (int i = start_id; i < mesh->Nelements; ++i)
649 : {
650 1093968674 : for (int j = 0; j < mesh->Np; ++j)
651 : {
652 1086101000 : if (max)
653 614096211 : value = std::max(value, f(i * mesh->Np + j, 0 /* unused */));
654 : else
655 474404620 : value = std::min(value, f(i * mesh->Np + j, 0 /* unused */));
656 : }
657 : }
658 :
659 : // find extreme value across all processes
660 : double reduced_value;
661 46788 : auto op = max ? MPI_MAX : MPI_MIN;
662 46788 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm);
663 :
664 : // dimensionalize the field if needed
665 46788 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
666 :
667 46788 : return reduced_value;
668 : }
669 :
670 : Point
671 99589536 : gllPoint(int local_elem_id, int local_node_id)
672 : {
673 99589536 : mesh_t * mesh = entireMesh();
674 :
675 99589536 : int id = local_elem_id * mesh->Np + local_node_id;
676 99589536 : Point p(mesh->x[id], mesh->y[id], mesh->z[id]);
677 : p *= scales.L_ref;
678 99589536 : return p;
679 : }
680 :
681 : Point
682 933120 : gllPointFace(int local_elem_id, int local_face_id, int local_node_id)
683 : {
684 933120 : mesh_t * mesh = entireMesh();
685 933120 : int face_id = mesh->EToB[local_elem_id * mesh->Nfaces + local_face_id];
686 933120 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
687 933120 : int id = mesh->vmapM[offset + local_node_id];
688 933120 : Point p(mesh->x[id], mesh->y[id], mesh->z[id]);
689 : p *= scales.L_ref;
690 933120 : return p;
691 : }
692 :
693 : Point
694 234240 : centroidFace(int local_elem_id, int local_face_id)
695 : {
696 234240 : mesh_t * mesh = entireMesh();
697 :
698 : double x_c = 0.0;
699 : double y_c = 0.0;
700 : double z_c = 0.0;
701 : double mass = 0.0;
702 :
703 234240 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
704 14150400 : for (int v = 0; v < mesh->Np; ++v)
705 : {
706 13916160 : int id = mesh->vmapM[offset + v];
707 13916160 : double mass_matrix = sgeo[mesh->Nsgeo * (offset + v) + WSJID];
708 13916160 : x_c += mesh->x[id] * mass_matrix;
709 13916160 : y_c += mesh->y[id] * mass_matrix;
710 13916160 : z_c += mesh->z[id] * mass_matrix;
711 13916160 : mass += mass_matrix;
712 : }
713 :
714 : Point c(x_c, y_c, z_c);
715 234240 : return c / mass * scales.L_ref;
716 : }
717 :
718 : Point
719 26997504 : centroid(int local_elem_id)
720 : {
721 26997504 : mesh_t * mesh = entireMesh();
722 :
723 : double x_c = 0.0;
724 : double y_c = 0.0;
725 : double z_c = 0.0;
726 : double mass = 0.0;
727 :
728 6192016128 : for (int v = 0; v < mesh->Np; ++v)
729 : {
730 6165018624 : int id = local_elem_id * mesh->Np + v;
731 6165018624 : double mass_matrix = vgeo[local_elem_id * mesh->Np * mesh->Nvgeo + JWID * mesh->Np + v];
732 6165018624 : x_c += mesh->x[id] * mass_matrix;
733 6165018624 : y_c += mesh->y[id] * mass_matrix;
734 6165018624 : z_c += mesh->z[id] * mass_matrix;
735 6165018624 : mass += mass_matrix;
736 : }
737 :
738 : Point c(x_c, y_c, z_c);
739 26997504 : return c / mass * scales.L_ref;
740 : }
741 :
742 : double
743 16250 : volume(const nek_mesh::NekMeshEnum pp_mesh)
744 : {
745 16250 : mesh_t * mesh = getMesh(pp_mesh);
746 16250 : double integral = 0.0;
747 :
748 3261198 : for (int k = 0; k < mesh->Nelements; ++k)
749 : {
750 3244948 : int offset = k * mesh->Np;
751 :
752 428948804 : for (int v = 0; v < mesh->Np; ++v)
753 425703856 : integral += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
754 : }
755 :
756 : // sum across all processes
757 : double total_integral;
758 16250 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
759 :
760 16250 : total_integral *= scales.V_ref;
761 :
762 16250 : return total_integral;
763 : }
764 :
765 : void
766 44275 : dimensionalizeVolume(double & integral)
767 : {
768 44275 : integral *= scales.V_ref;
769 44275 : }
770 :
771 : void
772 196 : dimensionalizeArea(double & integral)
773 : {
774 196 : integral *= scales.A_ref;
775 196 : }
776 :
777 : void
778 24400 : dimensionalizeVolumeIntegral(const field::NekFieldEnum & integrand,
779 : const Real & volume,
780 : double & integral)
781 : {
782 : // dimensionalize the field if needed
783 24400 : integral *= nondimensionalDivisor(integrand);
784 :
785 : // scale the volume integral
786 24400 : integral *= scales.V_ref;
787 :
788 : // for quantities with a relative scaling, we need to add back the reference
789 : // contribution to the volume integral
790 24400 : integral += nondimensionalAdditive(integrand) * volume;
791 24400 : }
792 :
793 : void
794 196 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
795 : const Real & area,
796 : double & integral)
797 : {
798 : // dimensionalize the field if needed
799 196 : integral *= nondimensionalDivisor(integrand);
800 :
801 : // scale the boundary integral
802 196 : integral *= scales.A_ref;
803 :
804 : // for quantities with a relative scaling, we need to add back the reference
805 : // contribution to the side integral
806 196 : integral += nondimensionalAdditive(integrand) * area;
807 196 : }
808 :
809 : void
810 30864 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
811 : const std::vector<int> & boundary_id,
812 : double & integral,
813 : const nek_mesh::NekMeshEnum pp_mesh)
814 : {
815 : // dimensionalize the field if needed
816 30864 : integral *= nondimensionalDivisor(integrand);
817 :
818 : // scale the boundary integral
819 30864 : integral *= scales.A_ref;
820 :
821 : // for quantities with a relative scaling, we need to add back the reference
822 : // contribution to the side integral; we need this form here to avoid a recursive loop
823 30864 : auto add = nondimensionalAdditive(integrand);
824 30864 : if (std::abs(add) > 1e-8)
825 120 : integral += add * area(boundary_id, pp_mesh);
826 30864 : }
827 :
828 : double
829 9888 : volumeIntegral(const field::NekFieldEnum & integrand, const Real & volume,
830 : const nek_mesh::NekMeshEnum pp_mesh)
831 : {
832 9888 : mesh_t * mesh = getMesh(pp_mesh);
833 :
834 9888 : double integral = 0.0;
835 :
836 : double (*f)(int, int);
837 9888 : f = solutionPointer(integrand);
838 :
839 1795120 : for (int k = 0; k < mesh->Nelements; ++k)
840 : {
841 1785232 : int offset = k * mesh->Np;
842 :
843 263259328 : for (int v = 0; v < mesh->Np; ++v)
844 261474096 : integral += f(offset + v, 0 /* unused */) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
845 : }
846 :
847 : // sum across all processes
848 : double total_integral;
849 9888 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
850 :
851 9888 : dimensionalizeVolumeIntegral(integrand, volume, total_integral);
852 :
853 9888 : return total_integral;
854 : }
855 :
856 : double
857 13070 : area(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
858 : {
859 13070 : mesh_t * mesh = getMesh(pp_mesh);
860 :
861 13070 : double integral = 0.0;
862 :
863 2026532 : for (int i = 0; i < mesh->Nelements; ++i)
864 : {
865 14094234 : for (int j = 0; j < mesh->Nfaces; ++j)
866 : {
867 12080772 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
868 :
869 12080772 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
870 : {
871 364424 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
872 7277964 : for (int v = 0; v < mesh->Nfp; ++v)
873 : {
874 6913540 : integral += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
875 : }
876 : }
877 : }
878 : }
879 :
880 : // sum across all processes
881 : double total_integral;
882 13070 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
883 :
884 13070 : dimensionalizeSideIntegral(field::unity, boundary_id, total_integral, pp_mesh);
885 :
886 13070 : return total_integral;
887 : }
888 :
889 : std::vector<dfloat>
890 12 : yPlus(const std::vector<int> & boundary_id, const unsigned int & index)
891 : {
892 12 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
893 :
894 : // compute the rate of strain tensor on device
895 12 : auto o_Sij = platform->o_memPool.reserve<dfloat>(2 * nrs->NVfields * nrs->fieldOffset);
896 12 : postProcessing::strainRate(nrs, true, nrs->o_U, o_Sij);
897 :
898 : // copy to host for evaluating wall-parallel stress
899 12 : dfloat * Sij = (dfloat *)calloc(o_Sij.size(), sizeof(dfloat));
900 12 : o_Sij.copyTo(Sij);
901 12 : o_Sij.free();
902 :
903 12 : double * wall_distance = (double *)nek::scPtr(index);
904 :
905 : // integrate over the boundaries in the mesh; each rank will compute contributions to the
906 : // x, y, and z components
907 12 : mesh_t * mesh = getMesh(nek_mesh::fluid);
908 :
909 : // TODO: This function only works correctly if the viscosity and rho is constant, because
910 : // otherwise we need to copy the viscosity and rho from device to host
911 : double mu;
912 12 : platform->options.getArgs("VISCOSITY", mu);
913 : double rho;
914 12 : platform->options.getArgs("DENSITY", rho);
915 12 : double nu = mu / rho;
916 :
917 12 : dfloat max_yp = -std::numeric_limits<float>::max();
918 12 : dfloat min_yp = std::numeric_limits<float>::max();
919 12 : dfloat avg_yp = 0.0;
920 12 : dfloat denom_yp = 0.0;
921 :
922 : std::vector<int> istride = {
923 12 : mesh->Nq * mesh->Nq, mesh->Nq, -1, -mesh->Nq, 1, -mesh->Nq * mesh->Nq};
924 :
925 1932 : for (int i = 0; i < mesh->Nelements; ++i)
926 : {
927 13440 : for (int j = 0; j < mesh->Nfaces; ++j)
928 : {
929 11520 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
930 :
931 11520 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
932 : {
933 384 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
934 :
935 6528 : for (int v = 0; v < mesh->Nfp; ++v)
936 : {
937 6144 : int surf_offset = mesh->Nsgeo * (offset + v);
938 6144 : int vol_id = mesh->vmapM[offset + v];
939 6144 : dfloat sWJ = sgeo[surf_offset + WSJID];
940 6144 : dfloat scale = 2 * mu;
941 :
942 6144 : dfloat n1 = sgeo[surf_offset + NXID];
943 6144 : dfloat n2 = sgeo[surf_offset + NYID];
944 6144 : dfloat n3 = sgeo[surf_offset + NZID];
945 :
946 6144 : dfloat s11 = Sij[vol_id + 0 * nrs->fieldOffset];
947 6144 : dfloat s12 = Sij[vol_id + 3 * nrs->fieldOffset];
948 6144 : dfloat s13 = Sij[vol_id + 5 * nrs->fieldOffset];
949 : dfloat s21 = s12;
950 6144 : dfloat s22 = Sij[vol_id + 1 * nrs->fieldOffset];
951 6144 : dfloat s23 = Sij[vol_id + 4 * nrs->fieldOffset];
952 : dfloat s31 = s13;
953 : dfloat s32 = s23;
954 6144 : dfloat s33 = Sij[vol_id + 2 * nrs->fieldOffset];
955 :
956 : // tau_{ij}n_j - (tau_{jk} n_k n_j) * n_i
957 6144 : dfloat f1 = scale * (s11 * n1 + s12 * n2 + s13 * n3);
958 6144 : dfloat f2 = scale * (s21 * n1 + s22 * n2 + s23 * n3);
959 6144 : dfloat f3 = scale * (s31 * n1 + s32 * n2 + s33 * n3);
960 :
961 6144 : f1 -= f1 * n1 * n1;
962 6144 : f2 -= f2 * n2 * n2;
963 6144 : f3 -= f3 * n3 * n3;
964 :
965 6144 : dfloat tauw = sqrt(f1 * f1 + f2 * f2 + f3 * f3);
966 6144 : dfloat utau = sqrt(tauw / rho);
967 :
968 : // need to shift when evaluating the wall distance, because we want the wall distance
969 : // at the nearest node away from the face, not precisely on the face
970 6144 : dfloat wd = wall_distance[vol_id + istride[j]];
971 :
972 : // check that we are not at a corner
973 6144 : if (wd > 1e-8)
974 : {
975 6144 : dfloat yplus = wd * utau / nu;
976 6144 : max_yp = std::max(yplus, max_yp);
977 6144 : min_yp = std::min(yplus, min_yp);
978 6144 : avg_yp += yplus * sWJ;
979 6144 : denom_yp += sWJ;
980 : }
981 : }
982 : }
983 : }
984 : }
985 :
986 : // max across all processes
987 : double total_max_yp;
988 12 : MPI_Allreduce(&max_yp, &total_max_yp, 1, MPI_DOUBLE, MPI_MAX, platform->comm.mpiComm);
989 :
990 : // min across all processes
991 : double total_min_yp;
992 12 : MPI_Allreduce(&min_yp, &total_min_yp, 1, MPI_DOUBLE, MPI_MIN, platform->comm.mpiComm);
993 :
994 : // sum across all processes
995 : double total_avg_yp;
996 12 : MPI_Allreduce(&avg_yp, &total_avg_yp, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
997 :
998 12 : if (total_avg_yp == 0)
999 0 : mooseError("Failed to find any eligible points on boundaries for computing y+!");
1000 :
1001 : double total_denom_yp;
1002 12 : MPI_Allreduce(&denom_yp, &total_denom_yp, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1003 :
1004 : // TODO: dimensionalize
1005 : // dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
1006 :
1007 : freePointer(Sij);
1008 24 : return {total_max_yp, total_min_yp, total_avg_yp / total_denom_yp};
1009 12 : }
1010 :
1011 : std::vector<dfloat>
1012 848 : viscousDrag(const std::vector<int> & boundary_id)
1013 : {
1014 848 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
1015 :
1016 : // compute the rate of strain tensor on device
1017 848 : auto o_Sij = platform->o_memPool.reserve<dfloat>(2 * nrs->NVfields * nrs->fieldOffset);
1018 848 : postProcessing::strainRate(nrs, true, nrs->o_U, o_Sij);
1019 :
1020 : // copy to host for evaluating drag
1021 848 : dfloat * Sij = (dfloat *)calloc(o_Sij.size(), sizeof(dfloat));
1022 848 : o_Sij.copyTo(Sij);
1023 848 : o_Sij.free();
1024 :
1025 : // integrate over the boundaries in the mesh; each rank will compute contributions to the
1026 : // x, y, and z components
1027 848 : mesh_t * mesh = getMesh(nek_mesh::fluid);
1028 :
1029 848 : double integral_x = 0.0;
1030 848 : double integral_y = 0.0;
1031 848 : double integral_z = 0.0;
1032 :
1033 : // TODO: This function only works correctly if the viscosity is constant, because
1034 : // otherwise we need to copy the viscosity from device to host
1035 : double mu;
1036 848 : platform->options.getArgs("VISCOSITY", mu);
1037 :
1038 51776 : for (int i = 0; i < mesh->Nelements; ++i)
1039 : {
1040 356496 : for (int j = 0; j < mesh->Nfaces; ++j)
1041 : {
1042 305568 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1043 :
1044 305568 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1045 : {
1046 4840 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1047 235080 : for (int v = 0; v < mesh->Nfp; ++v)
1048 : {
1049 230240 : int surf_offset = mesh->Nsgeo * (offset + v);
1050 230240 : int vol_id = mesh->vmapM[offset + v];
1051 230240 : dfloat sWJ = sgeo[surf_offset + WSJID];
1052 230240 : dfloat scale = -2 * mu * sWJ;
1053 :
1054 230240 : dfloat n1 = sgeo[surf_offset + NXID];
1055 230240 : dfloat n2 = sgeo[surf_offset + NYID];
1056 230240 : dfloat n3 = sgeo[surf_offset + NZID];
1057 :
1058 230240 : dfloat s11 = Sij[vol_id + 0 * nrs->fieldOffset];
1059 230240 : dfloat s12 = Sij[vol_id + 3 * nrs->fieldOffset];
1060 230240 : dfloat s13 = Sij[vol_id + 5 * nrs->fieldOffset];
1061 : dfloat s21 = s12;
1062 230240 : dfloat s22 = Sij[vol_id + 1 * nrs->fieldOffset];
1063 230240 : dfloat s23 = Sij[vol_id + 4 * nrs->fieldOffset];
1064 : dfloat s31 = s13;
1065 : dfloat s32 = s23;
1066 230240 : dfloat s33 = Sij[vol_id + 2 * nrs->fieldOffset];
1067 :
1068 230240 : dfloat dragx = scale * (s11 * n1 + s12 * n2 + s13 * n3);
1069 230240 : dfloat dragy = scale * (s21 * n1 + s22 * n2 + s23 * n3);
1070 230240 : dfloat dragz = scale * (s31 * n1 + s32 * n2 + s33 * n3);
1071 :
1072 230240 : integral_x += dragx;
1073 230240 : integral_y += dragy;
1074 230240 : integral_z += dragz;
1075 : }
1076 : }
1077 : }
1078 : }
1079 :
1080 : // sum across all processes
1081 : double total_integral_x;
1082 848 : MPI_Allreduce(&integral_x, &total_integral_x, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1083 :
1084 : double total_integral_y;
1085 848 : MPI_Allreduce(&integral_y, &total_integral_y, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1086 :
1087 : double total_integral_z;
1088 848 : MPI_Allreduce(&integral_z, &total_integral_z, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1089 :
1090 : // TODO: dimensionalize
1091 : // dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
1092 :
1093 : freePointer(Sij);
1094 1696 : return {total_integral_x, total_integral_y, total_integral_z};
1095 848 : }
1096 :
1097 : double
1098 17759 : sideIntegral(const std::vector<int> & boundary_id, const field::NekFieldEnum & integrand,
1099 : const nek_mesh::NekMeshEnum pp_mesh)
1100 : {
1101 17759 : mesh_t * mesh = getMesh(pp_mesh);
1102 :
1103 17758 : double integral = 0.0;
1104 :
1105 : double (*f)(int, int);
1106 17758 : f = solutionPointer(integrand);
1107 :
1108 2522262 : for (int i = 0; i < mesh->Nelements; ++i)
1109 : {
1110 17531528 : for (int j = 0; j < mesh->Nfaces; ++j)
1111 : {
1112 15027024 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1113 :
1114 15027024 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1115 : {
1116 505080 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1117 9982212 : for (int v = 0; v < mesh->Nfp; ++v)
1118 : {
1119 9477132 : integral +=
1120 9477132 : f(mesh->vmapM[offset + v], 0 /* unused */) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
1121 : }
1122 : }
1123 : }
1124 : }
1125 :
1126 : // sum across all processes
1127 : double total_integral;
1128 17758 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1129 :
1130 17758 : dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
1131 :
1132 17758 : return total_integral;
1133 : }
1134 :
1135 : double
1136 350 : massFlowrate(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1137 : {
1138 350 : mesh_t * mesh = getMesh(pp_mesh);
1139 350 : nrs_t * nrs = (nrs_t *)nrsPtr();
1140 :
1141 : // TODO: This function only works correctly if the density is constant, because
1142 : // otherwise we need to copy the density from device to host
1143 : double rho;
1144 350 : platform->options.getArgs("DENSITY", rho);
1145 :
1146 350 : double integral = 0.0;
1147 :
1148 139694 : for (int i = 0; i < mesh->Nelements; ++i)
1149 : {
1150 975408 : for (int j = 0; j < mesh->Nfaces; ++j)
1151 : {
1152 836064 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1153 :
1154 836064 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1155 : {
1156 7810 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1157 228040 : for (int v = 0; v < mesh->Nfp; ++v)
1158 : {
1159 220230 : int vol_id = mesh->vmapM[offset + v];
1160 220230 : int surf_offset = mesh->Nsgeo * (offset + v);
1161 :
1162 : double normal_velocity =
1163 220230 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1164 220230 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1165 220230 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1166 :
1167 220230 : integral += rho * normal_velocity * sgeo[surf_offset + WSJID];
1168 : }
1169 : }
1170 : }
1171 : }
1172 :
1173 : // sum across all processes
1174 : double total_integral;
1175 350 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1176 :
1177 : // dimensionalize the mass flux and area
1178 350 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1179 :
1180 350 : return total_integral;
1181 : }
1182 :
1183 : double
1184 692 : sideMassFluxWeightedIntegral(const std::vector<int> & boundary_id,
1185 : const field::NekFieldEnum & integrand,
1186 : const nek_mesh::NekMeshEnum pp_mesh)
1187 : {
1188 692 : mesh_t * mesh = getMesh(pp_mesh);
1189 692 : nrs_t * nrs = (nrs_t *)nrsPtr();
1190 :
1191 : // TODO: This function only works correctly if the density is constant, because
1192 : // otherwise we need to copy the density from device to host
1193 : double rho;
1194 692 : platform->options.getArgs("DENSITY", rho);
1195 :
1196 692 : double integral = 0.0;
1197 :
1198 : double (*f)(int, int);
1199 692 : f = solutionPointer(integrand);
1200 :
1201 186548 : for (int i = 0; i < mesh->Nelements; ++i)
1202 : {
1203 1300992 : for (int j = 0; j < mesh->Nfaces; ++j)
1204 : {
1205 1115136 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1206 :
1207 1115136 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1208 : {
1209 11982 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1210 402150 : for (int v = 0; v < mesh->Nfp; ++v)
1211 : {
1212 390168 : int vol_id = mesh->vmapM[offset + v];
1213 390168 : int surf_offset = mesh->Nsgeo * (offset + v);
1214 : double normal_velocity =
1215 390168 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1216 390168 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1217 390168 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1218 390168 : integral += f(vol_id, 0 /* unused */) * rho * normal_velocity * sgeo[surf_offset + WSJID];
1219 : }
1220 : }
1221 : }
1222 : }
1223 :
1224 : // sum across all processes
1225 : double total_integral;
1226 692 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1227 :
1228 : // dimensionalize the field if needed
1229 692 : total_integral *= nondimensionalDivisor(integrand);
1230 :
1231 : // dimensionalize the mass flux and area
1232 692 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1233 :
1234 : // for quantities with a relative scaling, we need to add back the reference
1235 : // contribution to the mass flux integral; we need this form here to avoid an infinite
1236 : // recursive loop
1237 692 : auto add = nondimensionalAdditive(integrand);
1238 692 : if (std::abs(add) > 1e-8)
1239 120 : total_integral += add * massFlowrate(boundary_id, pp_mesh);
1240 :
1241 692 : return total_integral;
1242 : }
1243 :
1244 : double
1245 36 : pressureSurfaceForce(const std::vector<int> & boundary_id, const Point & direction, const nek_mesh::NekMeshEnum pp_mesh)
1246 : {
1247 36 : mesh_t * mesh = getMesh(pp_mesh);
1248 36 : nrs_t * nrs = (nrs_t *)nrsPtr();
1249 :
1250 36 : double integral = 0.0;
1251 :
1252 20232 : for (int i = 0; i < mesh->Nelements; ++i)
1253 : {
1254 141372 : for (int j = 0; j < mesh->Nfaces; ++j)
1255 : {
1256 121176 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1257 :
1258 121176 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1259 : {
1260 12780 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1261 472860 : for (int v = 0; v < mesh->Nfp; ++v)
1262 : {
1263 460080 : int vol_id = mesh->vmapM[offset + v];
1264 460080 : int surf_offset = mesh->Nsgeo * (offset + v);
1265 :
1266 460080 : double p_normal = nrs->P[vol_id] * (sgeo[surf_offset + NXID] * direction(0) +
1267 460080 : sgeo[surf_offset + NYID] * direction(1) +
1268 460080 : sgeo[surf_offset + NZID] * direction(2));
1269 :
1270 460080 : integral += p_normal * sgeo[surf_offset + WSJID];
1271 : }
1272 : }
1273 : }
1274 : }
1275 :
1276 : // sum across all processes
1277 : double total_integral;
1278 36 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1279 :
1280 36 : dimensionalizeSideIntegral(field::pressure, boundary_id, total_integral, pp_mesh);
1281 :
1282 36 : return total_integral;
1283 : }
1284 :
1285 : double
1286 13416 : heatFluxIntegral(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1287 : {
1288 13416 : mesh_t * mesh = getMesh(pp_mesh);
1289 13416 : nrs_t * nrs = (nrs_t *)nrsPtr();
1290 :
1291 : // TODO: This function only works correctly if the conductivity is constant, because
1292 : // otherwise we need to copy the conductivity from device to host
1293 : double k;
1294 13416 : platform->options.getArgs("SCALAR00 DIFFUSIVITY", k);
1295 :
1296 13416 : double integral = 0.0;
1297 13416 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1298 :
1299 2946344 : for (int i = 0; i < mesh->Nelements; ++i)
1300 : {
1301 20530496 : for (int j = 0; j < mesh->Nfaces; ++j)
1302 : {
1303 17597568 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1304 :
1305 17597568 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1306 : {
1307 : // some inefficiency if an element has more than one face on the sideset of interest,
1308 : // because we will recompute the gradient in the element more than one time - but this
1309 : // is of little practical interest because this will be a minority of cases.
1310 216420 : gradient(mesh->Np, i, nrs->cds->S, grad_T, pp_mesh);
1311 :
1312 216420 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1313 5007924 : for (int v = 0; v < mesh->Nfp; ++v)
1314 : {
1315 : // special use of vol_id only when calling gradient(...), since we have written this
1316 : // function internally here so that it computes the gradient in a given element (as
1317 : // opposed to doing so for all elements at once, like in NekRS's
1318 : // gradientVolumeHex3D kernel
1319 4791504 : int vol_id = mesh->vmapM[offset + v] - i * mesh->Np;
1320 4791504 : int surf_offset = mesh->Nsgeo * (offset + v);
1321 :
1322 4791504 : double normal_grad_T = grad_T[vol_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1323 4791504 : grad_T[vol_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1324 4791504 : grad_T[vol_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1325 :
1326 4791504 : integral += -k * normal_grad_T * sgeo[surf_offset + WSJID];
1327 : }
1328 : }
1329 : }
1330 : }
1331 :
1332 : freePointer(grad_T);
1333 :
1334 : // sum across all processes
1335 : double total_integral;
1336 13416 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1337 :
1338 : // multiply by the reference heat flux and an area factor to dimensionalize
1339 13416 : total_integral *= scales.flux_ref * scales.A_ref;
1340 :
1341 13416 : return total_integral;
1342 : }
1343 :
1344 : void
1345 220452 : gradient(const int offset,
1346 : const int e,
1347 : const double * f,
1348 : double * grad_f,
1349 : const nek_mesh::NekMeshEnum pp_mesh)
1350 : {
1351 220452 : mesh_t * mesh = getMesh(pp_mesh);
1352 :
1353 1218324 : for (int k = 0; k < mesh->Nq; ++k)
1354 : {
1355 5807424 : for (int j = 0; j < mesh->Nq; ++j)
1356 : {
1357 30019584 : for (int i = 0; i < mesh->Nq; ++i)
1358 : {
1359 25210032 : const int gid = e * mesh->Np * mesh->Nvgeo + k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1360 25210032 : const double drdx = vgeo[gid + RXID * mesh->Np];
1361 25210032 : const double drdy = vgeo[gid + RYID * mesh->Np];
1362 25210032 : const double drdz = vgeo[gid + RZID * mesh->Np];
1363 25210032 : const double dsdx = vgeo[gid + SXID * mesh->Np];
1364 25210032 : const double dsdy = vgeo[gid + SYID * mesh->Np];
1365 25210032 : const double dsdz = vgeo[gid + SZID * mesh->Np];
1366 25210032 : const double dtdx = vgeo[gid + TXID * mesh->Np];
1367 25210032 : const double dtdy = vgeo[gid + TYID * mesh->Np];
1368 25210032 : const double dtdz = vgeo[gid + TZID * mesh->Np];
1369 :
1370 : // compute 'r' and 's' derivatives of (q_m) at node n
1371 : double dpdr = 0.f, dpds = 0.f, dpdt = 0.f;
1372 :
1373 172864224 : for (int n = 0; n < mesh->Nq; ++n)
1374 : {
1375 147654192 : const double Dr = mesh->D[i * mesh->Nq + n];
1376 147654192 : const double Ds = mesh->D[j * mesh->Nq + n];
1377 147654192 : const double Dt = mesh->D[k * mesh->Nq + n];
1378 :
1379 147654192 : dpdr += Dr * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + j * mesh->Nq + n];
1380 147654192 : dpds += Ds * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + n * mesh->Nq + i];
1381 147654192 : dpdt += Dt * f[e * mesh->Np + n * mesh->Nq * mesh->Nq + j * mesh->Nq + i];
1382 : }
1383 :
1384 25210032 : const int id = k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1385 25210032 : grad_f[id + 0 * offset] = drdx * dpdr + dsdx * dpds + dtdx * dpdt;
1386 25210032 : grad_f[id + 1 * offset] = drdy * dpdr + dsdy * dpds + dtdy * dpdt;
1387 25210032 : grad_f[id + 2 * offset] = drdz * dpdr + dsdz * dpds + dtdz * dpdt;
1388 : }
1389 : }
1390 : }
1391 220452 : }
1392 :
1393 : bool
1394 201 : isHeatFluxBoundary(const int boundary)
1395 : {
1396 : // the heat flux boundary is now named 'codedFixedGradient', but 'fixedGradient'
1397 : // will be present for backwards compatibility
1398 402 : return (bcMap::text(boundary, "scalar00") == "fixedGradient") ||
1399 804 : (bcMap::text(boundary, "scalar00") == "codedFixedGradient");
1400 : }
1401 :
1402 : bool
1403 17 : isMovingMeshBoundary(const int boundary)
1404 : {
1405 34 : return bcMap::text(boundary, "mesh") == "codedFixedValue";
1406 : }
1407 :
1408 : bool
1409 0 : isTemperatureBoundary(const int boundary)
1410 : {
1411 0 : return bcMap::text(boundary, "scalar00") == "fixedValue";
1412 : }
1413 :
1414 : const std::string
1415 1 : temperatureBoundaryType(const int boundary)
1416 : {
1417 2 : return bcMap::text(boundary, "scalar00");
1418 : }
1419 :
1420 : int
1421 1605 : polynomialOrder()
1422 : {
1423 1605 : return entireMesh()->N;
1424 : }
1425 :
1426 : int
1427 807 : Nelements()
1428 : {
1429 807 : int n_local = entireMesh()->Nelements;
1430 : int n_global;
1431 807 : MPI_Allreduce(&n_local, &n_global, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm);
1432 807 : return n_global;
1433 : }
1434 :
1435 : int
1436 9047269 : Nfaces()
1437 : {
1438 9047269 : return entireMesh()->Nfaces;
1439 : }
1440 :
1441 : int
1442 808 : dim()
1443 : {
1444 808 : return entireMesh()->dim;
1445 : }
1446 :
1447 : int
1448 0 : NfaceVertices()
1449 : {
1450 0 : return entireMesh()->NfaceVertices;
1451 : }
1452 :
1453 : int
1454 807 : NboundaryFaces()
1455 : {
1456 807 : return entireMesh()->NboundaryFaces;
1457 : }
1458 :
1459 : int
1460 7341 : NboundaryID()
1461 : {
1462 7341 : if (entireMesh()->cht)
1463 182 : return nekData.NboundaryIDt;
1464 : else
1465 7159 : return nekData.NboundaryID;
1466 : }
1467 :
1468 : bool
1469 2553 : validBoundaryIDs(const std::vector<int> & boundary_id, int & first_invalid_id, int & n_boundaries)
1470 : {
1471 2553 : n_boundaries = NboundaryID();
1472 :
1473 : bool valid_boundary_ids = true;
1474 5893 : for (const auto & b : boundary_id)
1475 : {
1476 3340 : if ((b > n_boundaries) || (b <= 0))
1477 : {
1478 3 : first_invalid_id = b;
1479 : valid_boundary_ids = false;
1480 : }
1481 : }
1482 :
1483 2553 : return valid_boundary_ids;
1484 : }
1485 :
1486 : double
1487 35440 : get_scalar01(const int id, const int surf_offset = 0)
1488 : {
1489 35440 : nrs_t * nrs = (nrs_t *)nrsPtr();
1490 35440 : return nrs->cds->S[id + 1 * scalarFieldOffset()];
1491 : }
1492 :
1493 : double
1494 61552 : get_scalar02(const int id, const int surf_offset = 0)
1495 : {
1496 61552 : nrs_t * nrs = (nrs_t *)nrsPtr();
1497 61552 : return nrs->cds->S[id + 2 * scalarFieldOffset()];
1498 : }
1499 :
1500 : double
1501 26224 : get_scalar03(const int id, const int surf_offset = 0)
1502 : {
1503 26224 : nrs_t * nrs = (nrs_t *)nrsPtr();
1504 26224 : return nrs->cds->S[id + 3 * scalarFieldOffset()];
1505 : }
1506 :
1507 : double
1508 505472 : get_usrwrk00(const int id, const int surf_offset = 0)
1509 : {
1510 505472 : nrs_t * nrs = (nrs_t *)nrsPtr();
1511 505472 : return nrs->usrwrk[id];
1512 : }
1513 :
1514 : double
1515 20672 : get_usrwrk01(const int id, const int surf_offset = 0)
1516 : {
1517 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1518 20672 : return nrs->usrwrk[id + nrs->fieldOffset];
1519 : }
1520 :
1521 : double
1522 20672 : get_usrwrk02(const int id, const int surf_offset = 0)
1523 : {
1524 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1525 20672 : return nrs->usrwrk[id + 2 * nrs->fieldOffset];
1526 : }
1527 :
1528 : double
1529 1233561952 : get_temperature(const int id, const int surf_offset)
1530 : {
1531 1233561952 : nrs_t * nrs = (nrs_t *)nrsPtr();
1532 1233561952 : return nrs->cds->S[id];
1533 : }
1534 :
1535 : double
1536 4032 : get_flux(const int id, const int surf_offset)
1537 : {
1538 : // TODO: this function does not support non-constant thermal conductivity
1539 : double k;
1540 4032 : platform->options.getArgs("SCALAR00 DIFFUSIVITY", k);
1541 :
1542 4032 : nrs_t * nrs = (nrs_t *)nrsPtr();
1543 :
1544 : // this call of nek_mesh::all should be fine because flux is not a 'field' which can be
1545 : // provided to the postprocessors which have the option to operate only on part of the mesh
1546 4032 : auto mesh = getMesh(nek_mesh::all);
1547 4032 : int elem_id = id / mesh->Np;
1548 4032 : int vertex_id = id % mesh->Np;
1549 :
1550 : // This function is slightly inefficient, because we compute grad(T) for all nodes in
1551 : // an element even though we only call this function for one node at a time
1552 4032 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1553 4032 : gradient(mesh->Np, elem_id, nrs->cds->S, grad_T, nek_mesh::all);
1554 :
1555 4032 : double normal_grad_T = grad_T[vertex_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1556 4032 : grad_T[vertex_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1557 4032 : grad_T[vertex_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1558 : freePointer(grad_T);
1559 :
1560 4032 : return -k * normal_grad_T;
1561 : }
1562 :
1563 : double
1564 134030256 : get_pressure(const int id, const int surf_offset)
1565 : {
1566 134030256 : nrs_t * nrs = (nrs_t *)nrsPtr();
1567 134030256 : return nrs->P[id];
1568 : }
1569 :
1570 : double
1571 55253552 : get_unity(const int /* id */, const int surf_offset)
1572 : {
1573 55253552 : return 1.0;
1574 : }
1575 :
1576 : double
1577 128478784 : get_velocity_x(const int id, const int surf_offset)
1578 : {
1579 128478784 : nrs_t * nrs = (nrs_t *)nrsPtr();
1580 128478784 : return nrs->U[id + 0 * nrs->fieldOffset];
1581 : }
1582 :
1583 : double
1584 65809984 : get_velocity_y(const int id, const int surf_offset)
1585 : {
1586 65809984 : nrs_t * nrs = (nrs_t *)nrsPtr();
1587 65809984 : return nrs->U[id + 1 * nrs->fieldOffset];
1588 : }
1589 :
1590 : double
1591 45871728 : get_velocity_z(const int id, const int surf_offset)
1592 : {
1593 45871728 : nrs_t * nrs = (nrs_t *)nrsPtr();
1594 45871728 : return nrs->U[id + 2 * nrs->fieldOffset];
1595 : }
1596 :
1597 : double
1598 2027844 : get_velocity(const int id, const int surf_offset)
1599 : {
1600 2027844 : nrs_t * nrs = (nrs_t *)nrsPtr();
1601 2027844 : int offset = nrs->fieldOffset;
1602 :
1603 2027844 : return std::sqrt(nrs->U[id + 0 * offset] * nrs->U[id + 0 * offset] +
1604 2027844 : nrs->U[id + 1 * offset] * nrs->U[id + 1 * offset] +
1605 2027844 : nrs->U[id + 2 * offset] * nrs->U[id + 2 * offset]);
1606 : }
1607 :
1608 : double
1609 47744 : get_velocity_x_squared(const int id, const int surf_offset)
1610 : {
1611 47744 : return std::pow(get_velocity_x(id, surf_offset), 2);
1612 : }
1613 :
1614 : double
1615 47744 : get_velocity_y_squared(const int id, const int surf_offset)
1616 : {
1617 47744 : return std::pow(get_velocity_y(id, surf_offset), 2);
1618 : }
1619 :
1620 : double
1621 52912 : get_velocity_z_squared(const int id, const int surf_offset)
1622 : {
1623 52912 : return std::pow(get_velocity_z(id, surf_offset), 2);
1624 : }
1625 :
1626 : void
1627 244 : checkFieldValidity(const field::NekWriteEnum & field)
1628 : {
1629 244 : switch (field)
1630 : {
1631 244 : case field::flux:
1632 244 : if (!hasTemperatureVariable())
1633 0 : mooseError("Cannot get NekRS heat flux "
1634 : "because your Nek case files do not have a temperature variable!");
1635 : break;
1636 0 : case field::heat_source:
1637 0 : if (!hasTemperatureVariable())
1638 0 : mooseError("Cannot get NekRS heat source "
1639 : "because your Nek case files do not have a temperature variable!");
1640 : break;
1641 : case field::x_displacement:
1642 : case field::y_displacement:
1643 : case field::z_displacement:
1644 : case field::mesh_velocity_x:
1645 : case field::mesh_velocity_y:
1646 : case field::mesh_velocity_z:
1647 : break;
1648 0 : default:
1649 0 : mooseError("Unhandled NekWriteEnum in checkFieldValidity!");
1650 : }
1651 244 : }
1652 :
1653 : void
1654 189430 : checkFieldValidity(const field::NekFieldEnum & field)
1655 : {
1656 : // by placing this check here, as opposed to inside the NekFieldInterface,
1657 : // we can also leverage this error checking for the 'outputs' of NekRSProblem,
1658 : // which does not inherit from NekFieldInterface but still accesses the solutionPointers.
1659 : // If this gets moved elsewhere, need to be sure to add dedicated testing for
1660 : // the 'outputs' on NekRSProblem.
1661 :
1662 : // TODO: would be nice for NekRSProblem to only access field information via the
1663 : // NekFieldInterface; refactor later
1664 :
1665 189430 : switch (field)
1666 : {
1667 83860 : case field::temperature:
1668 83860 : if (!hasTemperatureVariable())
1669 1 : mooseError("Cannot find 'temperature' "
1670 : "because your Nek case files do not have a temperature variable!");
1671 : break;
1672 85 : case field::scalar01:
1673 85 : if (!hasScalarVariable(1))
1674 1 : mooseError("Cannot find 'scalar01' "
1675 : "because your Nek case files do not have a scalar01 variable!");
1676 : break;
1677 201 : case field::scalar02:
1678 201 : if (!hasScalarVariable(2))
1679 1 : mooseError("Cannot find 'scalar02' "
1680 : "because your Nek case files do not have a scalar02 variable!");
1681 : break;
1682 57 : case field::scalar03:
1683 57 : if (!hasScalarVariable(3))
1684 1 : mooseError("Cannot find 'scalar03' "
1685 : "because your Nek case files do not have a scalar03 variable!");
1686 : break;
1687 454 : case field::usrwrk00:
1688 454 : if (n_usrwrk_slots < 1)
1689 2 : mooseError("Cannot find 'usrwrk00' because you have only allocated 'n_usrwrk_slots = " +
1690 1 : std::to_string(n_usrwrk_slots) + "'");
1691 : break;
1692 42 : case field::usrwrk01:
1693 42 : if (n_usrwrk_slots < 2)
1694 2 : mooseError("Cannot find 'usrwrk01' because you have only allocated 'n_usrwrk_slots = " +
1695 1 : std::to_string(n_usrwrk_slots) + "'");
1696 : break;
1697 42 : case field::usrwrk02:
1698 42 : if (n_usrwrk_slots < 3)
1699 2 : mooseError("Cannot find 'usrwrk02' because you have only allocated 'n_usrwrk_slots = " +
1700 1 : std::to_string(n_usrwrk_slots) + "'");
1701 : break;
1702 : }
1703 189423 : }
1704 :
1705 244 : double (*solutionPointer(const field::NekWriteEnum & field))(int, int)
1706 : {
1707 : double (*f)(int, int);
1708 :
1709 244 : checkFieldValidity(field);
1710 :
1711 244 : switch (field)
1712 : {
1713 244 : case field::flux:
1714 : f = &get_flux;
1715 : break;
1716 0 : default:
1717 0 : mooseError("Unhandled NekWriteEnum in solutionPointer!");
1718 : }
1719 :
1720 244 : return f;
1721 : }
1722 :
1723 186212 : double (*solutionPointer(const field::NekFieldEnum & field))(int, int)
1724 : {
1725 : // we include this here as well, in addition to within the NekFieldInterface, because
1726 : // the NekRSProblem accesses these methods without inheriting from NekFieldInterface
1727 186212 : checkFieldValidity(field);
1728 :
1729 : double (*f)(int, int);
1730 :
1731 186212 : switch (field)
1732 : {
1733 : case field::velocity_x:
1734 : f = &get_velocity_x;
1735 : break;
1736 17796 : case field::velocity_y:
1737 : f = &get_velocity_y;
1738 17796 : break;
1739 14808 : case field::velocity_z:
1740 : f = &get_velocity_z;
1741 14808 : break;
1742 196 : case field::velocity:
1743 : f = &get_velocity;
1744 196 : break;
1745 0 : case field::velocity_component:
1746 0 : mooseError("The 'velocity_component' field is not compatible with the solutionPointer "
1747 : "interface!");
1748 : break;
1749 48 : case field::velocity_x_squared:
1750 : f = &get_velocity_x_squared;
1751 48 : break;
1752 48 : case field::velocity_y_squared:
1753 : f = &get_velocity_y_squared;
1754 48 : break;
1755 52 : case field::velocity_z_squared:
1756 : f = &get_velocity_z_squared;
1757 52 : break;
1758 82736 : case field::temperature:
1759 : f = &get_temperature;
1760 82736 : break;
1761 30416 : case field::pressure:
1762 : f = &get_pressure;
1763 30416 : break;
1764 48 : case field::scalar01:
1765 : f = &get_scalar01;
1766 48 : break;
1767 176 : case field::scalar02:
1768 : f = &get_scalar02;
1769 176 : break;
1770 32 : case field::scalar03:
1771 : f = &get_scalar03;
1772 32 : break;
1773 8108 : case field::unity:
1774 : f = &get_unity;
1775 8108 : break;
1776 420 : case field::usrwrk00:
1777 : f = &get_usrwrk00;
1778 420 : break;
1779 16 : case field::usrwrk01:
1780 : f = &get_usrwrk01;
1781 16 : break;
1782 16 : case field::usrwrk02:
1783 : f = &get_usrwrk02;
1784 16 : break;
1785 0 : default:
1786 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1787 : }
1788 :
1789 186212 : return f;
1790 : }
1791 :
1792 : void
1793 129 : initializeDimensionalScales(const double U,
1794 : const double T,
1795 : const double dT,
1796 : const double L,
1797 : const double rho,
1798 : const double Cp,
1799 : const double s01,
1800 : const double ds01,
1801 : const double s02,
1802 : const double ds02,
1803 : const double s03,
1804 : const double ds03)
1805 : {
1806 129 : scales.U_ref = U;
1807 129 : scales.T_ref = T;
1808 129 : scales.dT_ref = dT;
1809 129 : scales.L_ref = L;
1810 129 : scales.A_ref = L * L;
1811 129 : scales.V_ref = L * L * L;
1812 129 : scales.rho_ref = rho;
1813 129 : scales.Cp_ref = Cp;
1814 129 : scales.t_ref = L / U;
1815 129 : scales.P_ref = rho * U * U;
1816 :
1817 129 : scales.s01_ref = s01;
1818 129 : scales.ds01_ref = ds01;
1819 129 : scales.s02_ref = s02;
1820 129 : scales.ds02_ref = ds02;
1821 129 : scales.s03_ref = s03;
1822 129 : scales.ds03_ref = ds03;
1823 :
1824 129 : scales.flux_ref = rho * U * Cp * dT;
1825 129 : scales.source_ref = scales.flux_ref / L;
1826 129 : }
1827 :
1828 : double
1829 433 : referenceLength()
1830 : {
1831 433 : return scales.L_ref;
1832 : }
1833 :
1834 : double
1835 276488 : referenceTime()
1836 : {
1837 276488 : return scales.t_ref;
1838 : }
1839 :
1840 : double
1841 241 : referenceArea()
1842 : {
1843 241 : return scales.A_ref;
1844 : }
1845 :
1846 : double
1847 982 : referenceVolume()
1848 : {
1849 982 : return scales.V_ref;
1850 : }
1851 :
1852 : Real
1853 51191106 : nondimensionalAdditive(const field::NekFieldEnum & field)
1854 : {
1855 51191106 : switch (field)
1856 : {
1857 26472776 : case field::temperature:
1858 26472776 : return scales.T_ref;
1859 8660 : case field::scalar01:
1860 8660 : return scales.s01_ref;
1861 40904 : case field::scalar02:
1862 40904 : return scales.s02_ref;
1863 5576 : case field::scalar03:
1864 5576 : return scales.s03_ref;
1865 : default:
1866 : return 0;
1867 : }
1868 : }
1869 :
1870 : Real
1871 18964 : nondimensionalAdditive(const field::NekWriteEnum & field)
1872 : {
1873 18964 : switch (field)
1874 : {
1875 18964 : case field::flux:
1876 : case field::heat_source:
1877 : case field::x_displacement:
1878 : case field::y_displacement:
1879 : case field::z_displacement:
1880 : case field::mesh_velocity_x:
1881 : case field::mesh_velocity_y:
1882 : case field::mesh_velocity_z:
1883 18964 : return 0.0;
1884 0 : default:
1885 0 : mooseError("Unhandled NekWriteEnum in nondimensionalAdditive!");
1886 : }
1887 : }
1888 :
1889 : Real
1890 46821 : nondimensionalDivisor(const field::NekWriteEnum & field)
1891 : {
1892 46821 : switch (field)
1893 : {
1894 42705 : case field::flux:
1895 42705 : return scales.flux_ref;
1896 3167 : case field::heat_source:
1897 3167 : return scales.source_ref;
1898 904 : case field::x_displacement:
1899 : case field::y_displacement:
1900 : case field::z_displacement:
1901 904 : return scales.L_ref;
1902 45 : case field::mesh_velocity_x:
1903 : case field::mesh_velocity_y:
1904 : case field::mesh_velocity_z:
1905 45 : return scales.U_ref;
1906 0 : default:
1907 0 : mooseError("Unhandled NekWriteEnum in nondimensionalDivisor!");
1908 : }
1909 : }
1910 :
1911 : Real
1912 51728832 : nondimensionalDivisor(const field::NekFieldEnum & field)
1913 : {
1914 51728832 : switch (field)
1915 : {
1916 19106132 : case field::velocity_x:
1917 : case field::velocity_y:
1918 : case field::velocity_z:
1919 : case field::velocity:
1920 : case field::velocity_component:
1921 19106132 : return scales.U_ref;
1922 5336 : case field::velocity_x_squared:
1923 : case field::velocity_y_squared:
1924 : case field::velocity_z_squared:
1925 5336 : return scales.U_ref * scales.U_ref;
1926 26472776 : case field::temperature:
1927 26472776 : return scales.dT_ref;
1928 6065815 : case field::pressure:
1929 6065815 : return scales.P_ref;
1930 8660 : case field::scalar01:
1931 8660 : return scales.ds01_ref;
1932 40904 : case field::scalar02:
1933 40904 : return scales.ds02_ref;
1934 5576 : case field::scalar03:
1935 5576 : return scales.ds03_ref;
1936 : case field::unity:
1937 : // no dimensionalization needed
1938 : return 1.0;
1939 445 : case field::usrwrk00:
1940 445 : return scratchUnits(0);
1941 25 : case field::usrwrk01:
1942 25 : return scratchUnits(1);
1943 25 : case field::usrwrk02:
1944 25 : return scratchUnits(2);
1945 0 : default:
1946 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1947 : }
1948 : }
1949 :
1950 : Real
1951 495 : scratchUnits(const int slot)
1952 : {
1953 : // if (indices.flux != -1 && slot == indices.flux / nekrs::fieldOffset())
1954 : // return scales.flux_ref;
1955 : // else if (indices.heat_source != -1 && slot == indices.heat_source / nekrs::fieldOffset())
1956 : // return scales.source_ref;
1957 495 : if (is_nondimensional)
1958 : {
1959 66 : mooseDoOnce(mooseWarning(
1960 : "The units of 'usrwrk0" + std::to_string(slot) +
1961 : "' are unknown, so we cannot dimensionalize any objects using 'field = usrwrk0" +
1962 : std::to_string(slot) +
1963 : "'. The output for this quantity will be given in non-dimensional form.\n\nYou will need "
1964 : "to manipulate the data manually from Cardinal if you need to dimensionalize it."));
1965 : }
1966 :
1967 492 : return 1.0;
1968 : }
1969 :
1970 : void
1971 800 : nondimensional(const bool n)
1972 : {
1973 800 : is_nondimensional = n;
1974 800 : }
1975 :
1976 : template <>
1977 : MPI_Datatype
1978 179286 : resolveType<double>()
1979 : {
1980 179286 : return MPI_DOUBLE;
1981 : }
1982 :
1983 : template <>
1984 : MPI_Datatype
1985 6944 : resolveType<int>()
1986 : {
1987 6944 : return MPI_INT;
1988 : }
1989 :
1990 : } // end namespace nekrs
1991 :
1992 : #endif
|
