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 : nekrs::usrwrkIndices indices;
30 :
31 : namespace nekrs
32 : {
33 : static double setup_time;
34 :
35 : // various constants for controlling tolerances
36 : static double abs_tol;
37 : static double rel_tol;
38 :
39 : void
40 348 : setAbsoluteTol(double tol)
41 : {
42 348 : abs_tol = tol;
43 348 : }
44 :
45 : void
46 348 : setRelativeTol(double tol)
47 : {
48 348 : rel_tol = tol;
49 348 : }
50 :
51 : void
52 810 : setNekSetupTime(const double & time)
53 : {
54 810 : setup_time = time;
55 810 : }
56 :
57 : double
58 799 : getNekSetupTime()
59 : {
60 799 : return setup_time;
61 : }
62 :
63 : void
64 742 : setStartTime(const double & start)
65 : {
66 1484 : platform->options.setArgs("START TIME", to_string_f(start));
67 742 : }
68 :
69 : void
70 52 : write_usrwrk_field_file(const int & slot, const std::string & prefix, const dfloat & time, const int & step, const bool & write_coords)
71 : {
72 52 : int num_bytes = fieldOffset() * sizeof(dfloat);
73 :
74 52 : nrs_t * nrs = (nrs_t *)nrsPtr();
75 52 : occa::memory o_write = platform->device.malloc(num_bytes);
76 104 : o_write.copyFrom(nrs->o_usrwrk, num_bytes /* length we are copying */,
77 52 : 0 /* where to place data */, num_bytes * slot /* where to source data */);
78 :
79 52 : occa::memory o_null;
80 52 : writeFld(prefix.c_str(), time, step, write_coords, 1 /* FP64 */, o_null, o_null, o_write, 1);
81 52 : }
82 :
83 : void
84 56 : write_field_file(const std::string & prefix, const dfloat time, const int & step)
85 : {
86 56 : nrs_t * nrs = (nrs_t *)nrsPtr();
87 :
88 : int Nscalar = 0;
89 56 : occa::memory o_s;
90 56 : if (nrs->Nscalar)
91 : {
92 56 : o_s = nrs->cds->o_S;
93 56 : Nscalar = nrs->Nscalar;
94 : }
95 :
96 56 : writeFld(
97 56 : prefix.c_str(), time, step, 1 /* coords */, 1 /* FP64 */, nrs->o_U, nrs->o_P, o_s, Nscalar);
98 56 : }
99 :
100 : void
101 812 : buildOnly(int buildOnly)
102 : {
103 812 : build_only = buildOnly;
104 812 : }
105 :
106 : int
107 103249 : buildOnly()
108 : {
109 103249 : return build_only;
110 : }
111 :
112 : bool
113 0 : hasCHT()
114 : {
115 0 : return entireMesh()->cht;
116 : }
117 :
118 : bool
119 36013 : hasMovingMesh()
120 : {
121 36013 : return platform->options.compareArgs("MOVING MESH", "TRUE");
122 : }
123 :
124 : bool
125 32598 : hasVariableDt()
126 : {
127 32598 : return platform->options.compareArgs("VARIABLE DT", "TRUE");
128 : }
129 :
130 : bool
131 953 : hasBlendingSolver()
132 : {
133 1906 : return !platform->options.compareArgs("MESH SOLVER", "NONE") && hasMovingMesh();
134 : }
135 :
136 : bool
137 954 : hasUserMeshSolver()
138 : {
139 1908 : return platform->options.compareArgs("MESH SOLVER", "NONE") && hasMovingMesh();
140 : }
141 :
142 : bool
143 1524 : endControlElapsedTime()
144 : {
145 3048 : return !platform->options.getArgs("STOP AT ELAPSED TIME").empty();
146 : }
147 :
148 : bool
149 1524 : endControlTime()
150 : {
151 1524 : return endTime() > 0;
152 : }
153 :
154 : bool
155 762 : endControlNumSteps()
156 : {
157 762 : return !endControlElapsedTime() && !endControlTime();
158 : }
159 :
160 : bool
161 262370574 : hasTemperatureVariable()
162 : {
163 262370574 : nrs_t * nrs = (nrs_t *)nrsPtr();
164 786661596 : return nrs->Nscalar ? platform->options.compareArgs("SCALAR00 IS TEMPERATURE", "TRUE") : false;
165 : }
166 :
167 : bool
168 680 : hasTemperatureSolve()
169 : {
170 680 : nrs_t * nrs = (nrs_t *)nrsPtr();
171 680 : return hasTemperatureVariable() ? nrs->cds->compute[0] : false;
172 : }
173 :
174 : bool
175 1123 : hasScalarVariable(int scalarId)
176 : {
177 1123 : nrs_t * nrs = (nrs_t *)nrsPtr();
178 1123 : return (scalarId < nrs->Nscalar);
179 : }
180 :
181 : bool
182 93 : hasHeatSourceKernel()
183 : {
184 93 : return static_cast<bool>(udf.sEqnSource);
185 : }
186 :
187 : bool
188 808 : isInitialized()
189 : {
190 808 : nrs_t * nrs = (nrs_t *)nrsPtr();
191 808 : return nrs;
192 : }
193 :
194 : int
195 144214 : scalarFieldOffset()
196 : {
197 144214 : nrs_t * nrs = (nrs_t *)nrsPtr();
198 144214 : return nrs->cds->fieldOffset[0];
199 : }
200 :
201 : int
202 4358488 : velocityFieldOffset()
203 : {
204 4358488 : nrs_t * nrs = (nrs_t *)nrsPtr();
205 4358488 : return nrs->fieldOffset;
206 : }
207 :
208 : int
209 21006 : fieldOffset()
210 : {
211 21006 : if (hasTemperatureVariable())
212 20966 : return scalarFieldOffset();
213 : else
214 40 : return velocityFieldOffset();
215 : }
216 :
217 : mesh_t *
218 262260886 : entireMesh()
219 : {
220 262260886 : if (hasTemperatureVariable())
221 172380524 : return temperatureMesh();
222 : else
223 89880362 : return flowMesh();
224 : }
225 :
226 : mesh_t *
227 90305473 : flowMesh()
228 : {
229 90305473 : nrs_t * nrs = (nrs_t *)nrsPtr();
230 90305473 : return nrs->meshV;
231 : }
232 :
233 : mesh_t *
234 173215445 : temperatureMesh()
235 : {
236 173215445 : nrs_t * nrs = (nrs_t *)nrsPtr();
237 173215445 : return nrs->cds->mesh[0];
238 : }
239 :
240 : mesh_t *
241 398932 : getMesh(const nek_mesh::NekMeshEnum pp_mesh)
242 : {
243 398932 : if (pp_mesh == nek_mesh::fluid)
244 1328 : return flowMesh();
245 397604 : else if (pp_mesh == nek_mesh::all)
246 397603 : return entireMesh();
247 : else
248 1 : mooseError("This object does not support operations on the solid part of the NekRS mesh!\n"
249 : "Valid options for 'mesh' are 'fluid' or 'all'.");
250 : }
251 :
252 : int
253 13137500 : commRank()
254 : {
255 13137500 : return platform->comm.mpiRank;
256 : }
257 :
258 : int
259 1066670 : commSize()
260 : {
261 1066670 : return platform->comm.mpiCommSize;
262 : }
263 :
264 : bool
265 810 : scratchAvailable()
266 : {
267 810 : nrs_t * nrs = (nrs_t *)nrsPtr();
268 :
269 : // Because these scratch spaces are available for whatever the user sees fit, it is
270 : // possible that the user wants to use these arrays for a _different_ purpose aside from
271 : // transferring in MOOSE values. In nekrs::setup, we call the UDF_Setup0, UDF_Setup,
272 : // and UDF_ExecuteStep routines. These scratch space arrays aren't initialized anywhere
273 : // else in the core base, so we will make sure to throw an error from MOOSE if these
274 : // arrays are already in use, because otherwise our MOOSE transfer might get overwritten
275 : // by whatever other operation the user is trying to do.
276 810 : if (nrs->usrwrk)
277 1 : return false;
278 :
279 : return true;
280 : }
281 :
282 : void
283 809 : initializeScratch(const unsigned int & n_slots)
284 : {
285 809 : if (n_slots == 0)
286 : return;
287 :
288 788 : nrs_t * nrs = (nrs_t *)nrsPtr();
289 788 : mesh_t * mesh = entireMesh();
290 :
291 : // clear them just to be sure
292 788 : freeScratch();
293 :
294 : // In order to make indexing simpler in the device user functions (which is where the
295 : // boundary conditions are then actually applied), we define these scratch arrays
296 : // as volume arrays.
297 788 : nrs->usrwrk = (double *)calloc(n_slots * fieldOffset(), sizeof(double));
298 788 : nrs->o_usrwrk = platform->device.malloc(n_slots * fieldOffset() * sizeof(double), nrs->usrwrk);
299 :
300 788 : n_usrwrk_slots = n_slots;
301 : }
302 :
303 : void
304 1522 : freeScratch()
305 : {
306 1522 : nrs_t * nrs = (nrs_t *)nrsPtr();
307 :
308 1522 : freePointer(nrs->usrwrk);
309 1522 : nrs->o_usrwrk.free();
310 1522 : }
311 :
312 : double
313 276 : viscosity()
314 : {
315 : dfloat mu;
316 276 : setupAide & options = platform->options;
317 276 : options.getArgs("VISCOSITY", mu);
318 :
319 : // because we set rho_ref, U_ref, and L_ref all equal to 1 if our input is dimensional,
320 : // we don't need to have separate treatments for dimensional vs. nondimensional cases
321 276 : dfloat Re = 1.0 / mu;
322 276 : return scales.rho_ref * scales.U_ref * scales.L_ref / Re;
323 : }
324 :
325 : double
326 92 : Pr()
327 : {
328 : dfloat rho, rho_cp, k;
329 92 : setupAide & options = platform->options;
330 92 : options.getArgs("DENSITY", rho);
331 92 : options.getArgs("SCALAR00 DENSITY", rho_cp);
332 92 : options.getArgs("SCALAR00 DIFFUSIVITY", k);
333 :
334 92 : dfloat Pe = 1.0 / k;
335 92 : dfloat conductivity = scales.rho_ref * scales.U_ref * scales.Cp_ref * scales.L_ref / Pe;
336 92 : dfloat Cp = rho_cp / rho * scales.Cp_ref;
337 92 : return viscosity() * Cp / conductivity;
338 : }
339 :
340 : void
341 1594 : interpolationMatrix(double * I, int starting_points, int ending_points)
342 : {
343 1594 : DegreeRaiseMatrix1D(starting_points - 1, ending_points - 1, I);
344 1594 : }
345 :
346 : void
347 2988324 : interpolateVolumeHex3D(const double * I, double * x, int N, double * Ix, int M)
348 : {
349 2988324 : double * Ix1 = (dfloat *)calloc(N * N * M, sizeof(double));
350 2988324 : double * Ix2 = (dfloat *)calloc(N * M * M, sizeof(double));
351 :
352 17185080 : for (int k = 0; k < N; ++k)
353 99718696 : for (int j = 0; j < N; ++j)
354 352372548 : for (int i = 0; i < M; ++i)
355 : {
356 : dfloat tmp = 0;
357 2077224048 : for (int n = 0; n < N; ++n)
358 1810373440 : tmp += I[i * N + n] * x[k * N * N + j * N + n];
359 266850608 : Ix1[k * N * M + j * M + i] = tmp;
360 : }
361 :
362 17185080 : for (int k = 0; k < N; ++k)
363 61015468 : for (int j = 0; j < M; ++j)
364 217899584 : for (int i = 0; i < M; ++i)
365 : {
366 : dfloat tmp = 0;
367 1045345368 : for (int n = 0; n < N; ++n)
368 874264496 : tmp += I[j * N + n] * Ix1[k * N * M + n * M + i];
369 171080872 : Ix2[k * M * M + j * M + i] = tmp;
370 : }
371 :
372 13762868 : for (int k = 0; k < M; ++k)
373 56333344 : for (int j = 0; j < M; ++j)
374 274821952 : for (int i = 0; i < M; ++i)
375 : {
376 : dfloat tmp = 0;
377 957314472 : for (int n = 0; n < N; ++n)
378 728051320 : tmp += I[k * N + n] * Ix2[n * M * M + j * M + i];
379 229263152 : Ix[k * M * M + j * M + i] = tmp;
380 : }
381 :
382 : freePointer(Ix1);
383 : freePointer(Ix2);
384 2988324 : }
385 :
386 : void
387 492928 : interpolateSurfaceFaceHex3D(
388 : double * scratch, const double * I, double * x, int N, double * Ix, int M)
389 : {
390 1952000 : for (int j = 0; j < N; ++j)
391 6683616 : for (int i = 0; i < M; ++i)
392 : {
393 : double tmp = 0;
394 24306720 : for (int n = 0; n < N; ++n)
395 : {
396 19082176 : tmp += I[i * N + n] * x[j * N + n];
397 : }
398 5224544 : scratch[j * M + i] = tmp;
399 : }
400 :
401 2395376 : for (int j = 0; j < M; ++j)
402 9675432 : for (int i = 0; i < M; ++i)
403 : {
404 : double tmp = 0;
405 27577896 : for (int n = 0; n < N; ++n)
406 : {
407 19804912 : tmp += I[j * N + n] * scratch[n * M + i];
408 : }
409 7772984 : Ix[j * M + i] = tmp;
410 : }
411 492928 : }
412 :
413 : void
414 95240 : displacementAndCounts(const std::vector<int> & base_counts,
415 : int * counts,
416 : int * displacement,
417 : const int multiplier = 1.0)
418 : {
419 484834 : for (int i = 0; i < commSize(); ++i)
420 389594 : counts[i] = base_counts[i] * multiplier;
421 :
422 95240 : displacement[0] = 0;
423 389594 : for (int i = 1; i < commSize(); i++)
424 294354 : displacement[i] = displacement[i - 1] + counts[i - 1];
425 95240 : }
426 :
427 : double
428 2240 : usrwrkVolumeIntegral(const unsigned int & slot, const nek_mesh::NekMeshEnum pp_mesh)
429 : {
430 2240 : nrs_t * nrs = (nrs_t *)nrsPtr();
431 2240 : const auto & mesh = getMesh(pp_mesh);
432 :
433 2240 : double integral = 0.0;
434 :
435 795480 : for (int k = 0; k < mesh->Nelements; ++k)
436 : {
437 793240 : int offset = k * mesh->Np;
438 :
439 237332376 : for (int v = 0; v < mesh->Np; ++v)
440 236539136 : integral += nrs->usrwrk[slot + offset + v] * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
441 : }
442 :
443 : // sum across all processes
444 : double total_integral;
445 2240 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
446 :
447 2240 : return total_integral;
448 : }
449 :
450 : void
451 1108 : scaleUsrwrk(const unsigned int & slot, const dfloat & value)
452 : {
453 1108 : nrs_t * nrs = (nrs_t *)nrsPtr();
454 1108 : mesh_t * mesh = getMesh(nek_mesh::all);
455 :
456 397600 : for (int k = 0; k < mesh->Nelements; ++k)
457 : {
458 396492 : int id = k * mesh->Np;
459 :
460 118662604 : for (int v = 0; v < mesh->Np; ++v)
461 118266112 : nrs->usrwrk[slot + id + v] *= value;
462 : }
463 1108 : }
464 :
465 : std::vector<double>
466 26459 : usrwrkSideIntegral(const unsigned int & slot,
467 : const std::vector<int> & boundary,
468 : const nek_mesh::NekMeshEnum pp_mesh)
469 : {
470 26459 : nrs_t * nrs = (nrs_t *)nrsPtr();
471 26459 : const auto & mesh = getMesh(pp_mesh);
472 :
473 26459 : std::vector<double> integral(boundary.size(), 0.0);
474 :
475 8260330 : for (int i = 0; i < mesh->Nelements; ++i)
476 : {
477 57637097 : for (int j = 0; j < mesh->Nfaces; ++j)
478 : {
479 49403226 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
480 :
481 49403226 : if (std::find(boundary.begin(), boundary.end(), face_id) != boundary.end())
482 : {
483 841040 : auto it = std::find(boundary.begin(), boundary.end(), face_id);
484 : auto b_index = it - boundary.begin();
485 :
486 841040 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
487 :
488 16010176 : for (int v = 0; v < mesh->Nfp; ++v)
489 15169136 : integral[b_index] += nrs->usrwrk[slot + mesh->vmapM[offset + v]] *
490 15169136 : sgeo[mesh->Nsgeo * (offset + v) + WSJID];
491 : }
492 : }
493 : }
494 :
495 : // sum across all processes; this can probably be done more efficiently
496 26459 : std::vector<double> total_integral(boundary.size(), 0.0);
497 53039 : for (std::size_t i = 0; i < boundary.size(); ++i)
498 26580 : MPI_Allreduce(&integral[i], &total_integral[i], 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
499 :
500 26459 : return total_integral;
501 : }
502 :
503 : void
504 0 : limitTemperature(const double * min_T, const double * max_T)
505 : {
506 : // if no limiters are provided, simply return
507 0 : if (!min_T && !max_T)
508 : return;
509 :
510 0 : double minimum = min_T ? *min_T : std::numeric_limits<double>::min();
511 0 : double maximum = max_T ? *max_T : std::numeric_limits<double>::max();
512 :
513 : // nondimensionalize if necessary
514 0 : minimum = (minimum - scales.T_ref) / scales.dT_ref;
515 0 : maximum = (maximum - scales.T_ref) / scales.dT_ref;
516 :
517 0 : nrs_t * nrs = (nrs_t *)nrsPtr();
518 0 : mesh_t * mesh = temperatureMesh();
519 :
520 0 : for (int i = 0; i < mesh->Nelements; ++i)
521 : {
522 0 : for (int j = 0; j < mesh->Np; ++j)
523 : {
524 0 : int id = i * mesh->Np + j;
525 :
526 0 : if (nrs->cds->S[id] < minimum)
527 0 : nrs->cds->S[id] = minimum;
528 0 : if (nrs->cds->S[id] > maximum)
529 0 : nrs->cds->S[id] = maximum;
530 : }
531 : }
532 :
533 : // when complete, copy to device
534 0 : nrs->cds->o_S.copyFrom(nrs->cds->S);
535 : }
536 :
537 : void
538 1504 : copyDeformationToDevice()
539 : {
540 1504 : mesh_t * mesh = entireMesh();
541 1504 : mesh->o_x.copyFrom(mesh->x);
542 1504 : mesh->o_y.copyFrom(mesh->y);
543 1504 : mesh->o_z.copyFrom(mesh->z);
544 1504 : mesh->update();
545 :
546 1504 : updateHostMeshParameters();
547 1504 : }
548 :
549 : void
550 807 : initializeHostMeshParameters()
551 : {
552 807 : mesh_t * mesh = entireMesh();
553 807 : sgeo = (dfloat *)calloc(mesh->o_sgeo.size(), sizeof(dfloat));
554 807 : vgeo = (dfloat *)calloc(mesh->o_vgeo.size(), sizeof(dfloat));
555 807 : }
556 :
557 : void
558 2311 : updateHostMeshParameters()
559 : {
560 2311 : mesh_t * mesh = entireMesh();
561 2311 : mesh->o_sgeo.copyTo(sgeo);
562 2311 : mesh->o_vgeo.copyTo(vgeo);
563 2311 : }
564 :
565 : dfloat *
566 128 : getSgeo()
567 : {
568 128 : return sgeo;
569 : }
570 :
571 : dfloat *
572 507 : getVgeo()
573 : {
574 507 : return vgeo;
575 : }
576 :
577 : double
578 24812 : sideExtremeValue(const std::vector<int> & boundary_id, const field::NekFieldEnum & field,
579 : const nek_mesh::NekMeshEnum pp_mesh, const bool max)
580 : {
581 24812 : mesh_t * mesh = getMesh(pp_mesh);
582 :
583 24812 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
584 :
585 : double (*f)(int);
586 24812 : f = solutionPointer(field);
587 :
588 2935048 : for (int i = 0; i < mesh->Nelements; ++i)
589 : {
590 20371652 : for (int j = 0; j < mesh->Nfaces; ++j)
591 : {
592 17461416 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
593 :
594 17461416 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
595 : {
596 282524 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
597 12898268 : for (int v = 0; v < mesh->Nfp; ++v)
598 : {
599 12615744 : if (max)
600 6383348 : value = std::max(value, f(mesh->vmapM[offset + v]));
601 : else
602 6369614 : value = std::min(value, f(mesh->vmapM[offset + v]));
603 : }
604 : }
605 : }
606 : }
607 :
608 : // find extreme value across all processes
609 : double reduced_value;
610 24812 : auto op = max ? MPI_MAX : MPI_MIN;
611 24812 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm);
612 :
613 : // dimensionalize the field if needed
614 24812 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
615 :
616 24812 : return reduced_value;
617 : }
618 :
619 : double
620 47966 : volumeExtremeValue(const field::NekFieldEnum & field, const nek_mesh::NekMeshEnum pp_mesh, const bool max)
621 : {
622 47966 : double value = max ? -std::numeric_limits<double>::max() : std::numeric_limits<double>::max();
623 :
624 : double (*f)(int);
625 47966 : f = solutionPointer(field);
626 :
627 : mesh_t * mesh;
628 : int start_id;
629 :
630 47966 : switch (pp_mesh)
631 : {
632 47954 : case nek_mesh::fluid:
633 : case nek_mesh::all:
634 : {
635 47954 : mesh = getMesh(pp_mesh);
636 : start_id = 0;
637 : break;
638 : }
639 12 : case nek_mesh::solid:
640 : {
641 12 : mesh = entireMesh();
642 12 : start_id = flowMesh()->Nelements;
643 12 : break;
644 : }
645 0 : default:
646 0 : mooseError("Unhandled NekMeshEnum in volumeExtremeValue");
647 : }
648 :
649 8299594 : for (int i = start_id; i < mesh->Nelements; ++i)
650 : {
651 1164789844 : for (int j = 0; j < mesh->Np; ++j)
652 : {
653 1156538216 : if (max)
654 662079182 : value = std::max(value, f(i * mesh->Np + j));
655 : else
656 496912621 : value = std::min(value, f(i * mesh->Np + j));
657 : }
658 : }
659 :
660 : // find extreme value across all processes
661 : double reduced_value;
662 47966 : auto op = max ? MPI_MAX : MPI_MIN;
663 47966 : MPI_Allreduce(&value, &reduced_value, 1, MPI_DOUBLE, op, platform->comm.mpiComm);
664 :
665 : // dimensionalize the field if needed
666 47966 : reduced_value = reduced_value * nondimensionalDivisor(field) + nondimensionalAdditive(field);
667 :
668 47966 : return reduced_value;
669 : }
670 :
671 : Point
672 218869152 : gllPoint(int local_elem_id, int local_node_id)
673 : {
674 218869152 : mesh_t * mesh = entireMesh();
675 :
676 218869152 : int id = local_elem_id * mesh->Np + local_node_id;
677 218869152 : Point p(mesh->x[id], mesh->y[id], mesh->z[id]);
678 : p *= scales.L_ref;
679 218869152 : return p;
680 : }
681 :
682 : Point
683 951552 : gllPointFace(int local_elem_id, int local_face_id, int local_node_id)
684 : {
685 951552 : mesh_t * mesh = entireMesh();
686 951552 : int face_id = mesh->EToB[local_elem_id * mesh->Nfaces + local_face_id];
687 951552 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
688 951552 : int id = mesh->vmapM[offset + local_node_id];
689 951552 : Point p(mesh->x[id], mesh->y[id], mesh->z[id]);
690 : p *= scales.L_ref;
691 951552 : return p;
692 : }
693 :
694 : Point
695 234240 : centroidFace(int local_elem_id, int local_face_id)
696 : {
697 234240 : mesh_t * mesh = entireMesh();
698 :
699 : double x_c = 0.0;
700 : double y_c = 0.0;
701 : double z_c = 0.0;
702 : double mass = 0.0;
703 :
704 234240 : int offset = local_elem_id * mesh->Nfaces * mesh->Nfp + local_face_id * mesh->Nfp;
705 14150400 : for (int v = 0; v < mesh->Np; ++v)
706 : {
707 13916160 : int id = mesh->vmapM[offset + v];
708 13916160 : double mass_matrix = sgeo[mesh->Nsgeo * (offset + v) + WSJID];
709 13916160 : x_c += mesh->x[id] * mass_matrix;
710 13916160 : y_c += mesh->y[id] * mass_matrix;
711 13916160 : z_c += mesh->z[id] * mass_matrix;
712 13916160 : mass += mass_matrix;
713 : }
714 :
715 : Point c(x_c, y_c, z_c);
716 234240 : return c / mass * scales.L_ref;
717 : }
718 :
719 : Point
720 26997504 : centroid(int local_elem_id)
721 : {
722 26997504 : mesh_t * mesh = entireMesh();
723 :
724 : double x_c = 0.0;
725 : double y_c = 0.0;
726 : double z_c = 0.0;
727 : double mass = 0.0;
728 :
729 6192016128 : for (int v = 0; v < mesh->Np; ++v)
730 : {
731 6165018624 : int id = local_elem_id * mesh->Np + v;
732 6165018624 : double mass_matrix = vgeo[local_elem_id * mesh->Np * mesh->Nvgeo + JWID * mesh->Np + v];
733 6165018624 : x_c += mesh->x[id] * mass_matrix;
734 6165018624 : y_c += mesh->y[id] * mass_matrix;
735 6165018624 : z_c += mesh->z[id] * mass_matrix;
736 6165018624 : mass += mass_matrix;
737 : }
738 :
739 : Point c(x_c, y_c, z_c);
740 26997504 : return c / mass * scales.L_ref;
741 : }
742 :
743 : double
744 17264 : volume(const nek_mesh::NekMeshEnum pp_mesh)
745 : {
746 17264 : mesh_t * mesh = getMesh(pp_mesh);
747 17264 : double integral = 0.0;
748 :
749 3673920 : for (int k = 0; k < mesh->Nelements; ++k)
750 : {
751 3656656 : int offset = k * mesh->Np;
752 :
753 519575744 : for (int v = 0; v < mesh->Np; ++v)
754 515919088 : integral += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
755 : }
756 :
757 : // sum across all processes
758 : double total_integral;
759 17264 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
760 :
761 17264 : total_integral *= scales.V_ref;
762 :
763 17264 : return total_integral;
764 : }
765 :
766 : void
767 49731 : dimensionalizeVolume(double & integral)
768 : {
769 49731 : integral *= scales.V_ref;
770 49731 : }
771 :
772 : void
773 216 : dimensionalizeArea(double & integral)
774 : {
775 216 : integral *= scales.A_ref;
776 216 : }
777 :
778 : void
779 50546 : dimensionalizeVolumeIntegral(const field::NekFieldEnum & integrand,
780 : const Real & volume,
781 : double & integral)
782 : {
783 : // dimensionalize the field if needed
784 50546 : integral *= nondimensionalDivisor(integrand);
785 :
786 : // scale the volume integral
787 50546 : integral *= scales.V_ref;
788 :
789 : // for quantities with a relative scaling, we need to add back the reference
790 : // contribution to the volume integral
791 50546 : integral += nondimensionalAdditive(integrand) * volume;
792 50546 : }
793 :
794 : void
795 296 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
796 : const Real & area,
797 : double & integral)
798 : {
799 : // dimensionalize the field if needed
800 296 : integral *= nondimensionalDivisor(integrand);
801 :
802 : // scale the boundary integral
803 296 : integral *= scales.A_ref;
804 :
805 : // for quantities with a relative scaling, we need to add back the reference
806 : // contribution to the side integral
807 296 : integral += nondimensionalAdditive(integrand) * area;
808 296 : }
809 :
810 : void
811 31712 : dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
812 : const std::vector<int> & boundary_id,
813 : double & integral,
814 : const nek_mesh::NekMeshEnum pp_mesh)
815 : {
816 : // dimensionalize the field if needed
817 31712 : integral *= nondimensionalDivisor(integrand);
818 :
819 : // scale the boundary integral
820 31712 : integral *= scales.A_ref;
821 :
822 : // for quantities with a relative scaling, we need to add back the reference
823 : // contribution to the side integral; we need this form here to avoid a recursive loop
824 31712 : auto add = nondimensionalAdditive(integrand);
825 31712 : if (std::abs(add) > 1e-8)
826 360 : integral += add * area(boundary_id, pp_mesh);
827 31712 : }
828 :
829 : double
830 10530 : volumeIntegral(const field::NekFieldEnum & integrand, const Real & volume,
831 : const nek_mesh::NekMeshEnum pp_mesh)
832 : {
833 10530 : mesh_t * mesh = getMesh(pp_mesh);
834 :
835 10530 : double integral = 0.0;
836 :
837 : double (*f)(int);
838 10530 : f = solutionPointer(integrand);
839 :
840 2054116 : for (int k = 0; k < mesh->Nelements; ++k)
841 : {
842 2043586 : int offset = k * mesh->Np;
843 :
844 312285298 : for (int v = 0; v < mesh->Np; ++v)
845 310241712 : integral += f(offset + v) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
846 : }
847 :
848 : // sum across all processes
849 : double total_integral;
850 10530 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
851 :
852 10530 : dimensionalizeVolumeIntegral(integrand, volume, total_integral);
853 :
854 10530 : return total_integral;
855 : }
856 :
857 : double
858 13500 : area(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
859 : {
860 13500 : mesh_t * mesh = getMesh(pp_mesh);
861 :
862 13500 : double integral = 0.0;
863 :
864 2520762 : for (int i = 0; i < mesh->Nelements; ++i)
865 : {
866 17550834 : for (int j = 0; j < mesh->Nfaces; ++j)
867 : {
868 15043572 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
869 :
870 15043572 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
871 : {
872 415204 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
873 8766664 : for (int v = 0; v < mesh->Nfp; ++v)
874 : {
875 8351460 : integral += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
876 : }
877 : }
878 : }
879 : }
880 :
881 : // sum across all processes
882 : double total_integral;
883 13500 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
884 :
885 13500 : dimensionalizeSideIntegral(field::unity, boundary_id, total_integral, pp_mesh);
886 :
887 13500 : return total_integral;
888 : }
889 :
890 : double
891 18177 : sideIntegral(const std::vector<int> & boundary_id, const field::NekFieldEnum & integrand,
892 : const nek_mesh::NekMeshEnum pp_mesh)
893 : {
894 18177 : mesh_t * mesh = getMesh(pp_mesh);
895 :
896 18176 : double integral = 0.0;
897 :
898 : double (*f)(int);
899 18176 : f = solutionPointer(integrand);
900 :
901 3014542 : for (int i = 0; i < mesh->Nelements; ++i)
902 : {
903 20974562 : for (int j = 0; j < mesh->Nfaces; ++j)
904 : {
905 17978196 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
906 :
907 17978196 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
908 : {
909 555758 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
910 11470402 : for (int v = 0; v < mesh->Nfp; ++v)
911 : {
912 10914644 : integral += f(mesh->vmapM[offset + v]) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
913 : }
914 : }
915 : }
916 : }
917 :
918 : // sum across all processes
919 : double total_integral;
920 18176 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
921 :
922 18176 : dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
923 :
924 18176 : return total_integral;
925 : }
926 :
927 : double
928 868 : massFlowrate(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
929 : {
930 868 : mesh_t * mesh = getMesh(pp_mesh);
931 868 : nrs_t * nrs = (nrs_t *)nrsPtr();
932 :
933 : // TODO: This function only works correctly if the density is constant, because
934 : // otherwise we need to copy the density from device to host
935 : double rho;
936 868 : platform->options.getArgs("DENSITY", rho);
937 :
938 868 : double integral = 0.0;
939 :
940 412012 : for (int i = 0; i < mesh->Nelements; ++i)
941 : {
942 2878008 : for (int j = 0; j < mesh->Nfaces; ++j)
943 : {
944 2466864 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
945 :
946 2466864 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
947 : {
948 21112 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
949 586564 : for (int v = 0; v < mesh->Nfp; ++v)
950 : {
951 565452 : int vol_id = mesh->vmapM[offset + v];
952 565452 : int surf_offset = mesh->Nsgeo * (offset + v);
953 :
954 : double normal_velocity =
955 565452 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
956 565452 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
957 565452 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
958 :
959 565452 : integral += rho * normal_velocity * sgeo[surf_offset + WSJID];
960 : }
961 : }
962 : }
963 : }
964 :
965 : // sum across all processes
966 : double total_integral;
967 868 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
968 :
969 : // dimensionalize the mass flux and area
970 868 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
971 :
972 868 : return total_integral;
973 : }
974 :
975 : double
976 1280 : sideMassFluxWeightedIntegral(const std::vector<int> & boundary_id,
977 : const field::NekFieldEnum & integrand,
978 : const nek_mesh::NekMeshEnum pp_mesh)
979 : {
980 1280 : mesh_t * mesh = getMesh(pp_mesh);
981 1280 : nrs_t * nrs = (nrs_t *)nrsPtr();
982 :
983 : // TODO: This function only works correctly if the density is constant, because
984 : // otherwise we need to copy the density from device to host
985 : double rho;
986 1280 : platform->options.getArgs("DENSITY", rho);
987 :
988 1280 : double integral = 0.0;
989 :
990 : double (*f)(int);
991 1280 : f = solutionPointer(integrand);
992 :
993 605936 : for (int i = 0; i < mesh->Nelements; ++i)
994 : {
995 4232592 : for (int j = 0; j < mesh->Nfaces; ++j)
996 : {
997 3627936 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
998 :
999 3627936 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1000 : {
1001 27498 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1002 914862 : for (int v = 0; v < mesh->Nfp; ++v)
1003 : {
1004 887364 : int vol_id = mesh->vmapM[offset + v];
1005 887364 : int surf_offset = mesh->Nsgeo * (offset + v);
1006 : double normal_velocity =
1007 887364 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1008 887364 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1009 887364 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1010 887364 : integral += f(vol_id) * rho * normal_velocity * sgeo[surf_offset + WSJID];
1011 : }
1012 : }
1013 : }
1014 : }
1015 :
1016 : // sum across all processes
1017 : double total_integral;
1018 1280 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1019 :
1020 : // dimensionalize the field if needed
1021 1280 : total_integral *= nondimensionalDivisor(integrand);
1022 :
1023 : // dimensionalize the mass flux and area
1024 1280 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1025 :
1026 : // for quantities with a relative scaling, we need to add back the reference
1027 : // contribution to the mass flux integral; we need this form here to avoid an infinite
1028 : // recursive loop
1029 1280 : auto add = nondimensionalAdditive(integrand);
1030 1280 : if (std::abs(add) > 1e-8)
1031 376 : total_integral += add * massFlowrate(boundary_id, pp_mesh);
1032 :
1033 1280 : return total_integral;
1034 : }
1035 :
1036 : double
1037 36 : pressureSurfaceForce(const std::vector<int> & boundary_id, const Point & direction, const nek_mesh::NekMeshEnum pp_mesh)
1038 : {
1039 36 : mesh_t * mesh = getMesh(pp_mesh);
1040 36 : nrs_t * nrs = (nrs_t *)nrsPtr();
1041 :
1042 36 : double integral = 0.0;
1043 :
1044 20232 : for (int i = 0; i < mesh->Nelements; ++i)
1045 : {
1046 141372 : for (int j = 0; j < mesh->Nfaces; ++j)
1047 : {
1048 121176 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1049 :
1050 121176 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1051 : {
1052 12780 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1053 472860 : for (int v = 0; v < mesh->Nfp; ++v)
1054 : {
1055 460080 : int vol_id = mesh->vmapM[offset + v];
1056 460080 : int surf_offset = mesh->Nsgeo * (offset + v);
1057 :
1058 460080 : double p_normal = nrs->P[vol_id] * (sgeo[surf_offset + NXID] * direction(0) +
1059 460080 : sgeo[surf_offset + NYID] * direction(1) +
1060 460080 : sgeo[surf_offset + NZID] * direction(2));
1061 :
1062 460080 : integral += p_normal * sgeo[surf_offset + WSJID];
1063 : }
1064 : }
1065 : }
1066 : }
1067 :
1068 : // sum across all processes
1069 : double total_integral;
1070 36 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1071 :
1072 36 : dimensionalizeSideIntegral(field::pressure, boundary_id, total_integral, pp_mesh);
1073 :
1074 36 : return total_integral;
1075 : }
1076 :
1077 : double
1078 11396 : heatFluxIntegral(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1079 : {
1080 11396 : mesh_t * mesh = getMesh(pp_mesh);
1081 11396 : nrs_t * nrs = (nrs_t *)nrsPtr();
1082 :
1083 : // TODO: This function only works correctly if the conductivity is constant, because
1084 : // otherwise we need to copy the conductivity from device to host
1085 : double k;
1086 11396 : platform->options.getArgs("SCALAR00 DIFFUSIVITY", k);
1087 :
1088 11396 : double integral = 0.0;
1089 11396 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1090 :
1091 2211276 : for (int i = 0; i < mesh->Nelements; ++i)
1092 : {
1093 15399160 : for (int j = 0; j < mesh->Nfaces; ++j)
1094 : {
1095 13199280 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1096 :
1097 13199280 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1098 : {
1099 : // some inefficiency if an element has more than one face on the sideset of interest,
1100 : // because we will recompute the gradient in the element more than one time - but this
1101 : // is of little practical interest because this will be a minority of cases.
1102 223308 : gradient(mesh->Np, i, nrs->cds->S, grad_T, pp_mesh);
1103 :
1104 223308 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1105 5745756 : for (int v = 0; v < mesh->Nfp; ++v)
1106 : {
1107 5522448 : int vol_id = mesh->vmapM[offset + v] - i * mesh->Np;
1108 5522448 : int surf_offset = mesh->Nsgeo * (offset + v);
1109 :
1110 5522448 : double normal_grad_T = grad_T[vol_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1111 5522448 : grad_T[vol_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1112 5522448 : grad_T[vol_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1113 :
1114 5522448 : integral += -k * normal_grad_T * sgeo[surf_offset + WSJID];
1115 : }
1116 : }
1117 : }
1118 : }
1119 :
1120 : freePointer(grad_T);
1121 :
1122 : // sum across all processes
1123 : double total_integral;
1124 11396 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1125 :
1126 : // multiply by the reference heat flux and an area factor to dimensionalize
1127 11396 : total_integral *= scales.flux_ref * scales.A_ref;
1128 :
1129 11396 : return total_integral;
1130 : }
1131 :
1132 : void
1133 223308 : gradient(const int offset,
1134 : const int e,
1135 : const double * f,
1136 : double * grad_f,
1137 : const nek_mesh::NekMeshEnum pp_mesh)
1138 : {
1139 223308 : mesh_t * mesh = getMesh(pp_mesh);
1140 :
1141 1298124 : for (int k = 0; k < mesh->Nq; ++k)
1142 : {
1143 6597264 : for (int j = 0; j < mesh->Nq; ++j)
1144 : {
1145 36425184 : for (int i = 0; i < mesh->Nq; ++i)
1146 : {
1147 30902736 : const int gid = e * mesh->Np * mesh->Nvgeo + k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1148 30902736 : const double drdx = vgeo[gid + RXID * mesh->Np];
1149 30902736 : const double drdy = vgeo[gid + RYID * mesh->Np];
1150 30902736 : const double drdz = vgeo[gid + RZID * mesh->Np];
1151 30902736 : const double dsdx = vgeo[gid + SXID * mesh->Np];
1152 30902736 : const double dsdy = vgeo[gid + SYID * mesh->Np];
1153 30902736 : const double dsdz = vgeo[gid + SZID * mesh->Np];
1154 30902736 : const double dtdx = vgeo[gid + TXID * mesh->Np];
1155 30902736 : const double dtdy = vgeo[gid + TYID * mesh->Np];
1156 30902736 : const double dtdz = vgeo[gid + TZID * mesh->Np];
1157 :
1158 : // compute 'r' and 's' derivatives of (q_m) at node n
1159 : double dpdr = 0.f, dpds = 0.f, dpdt = 0.f;
1160 :
1161 222432864 : for (int n = 0; n < mesh->Nq; ++n)
1162 : {
1163 191530128 : const double Dr = mesh->D[i * mesh->Nq + n];
1164 191530128 : const double Ds = mesh->D[j * mesh->Nq + n];
1165 191530128 : const double Dt = mesh->D[k * mesh->Nq + n];
1166 :
1167 191530128 : dpdr += Dr * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + j * mesh->Nq + n];
1168 191530128 : dpds += Ds * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + n * mesh->Nq + i];
1169 191530128 : dpdt += Dt * f[e * mesh->Np + n * mesh->Nq * mesh->Nq + j * mesh->Nq + i];
1170 : }
1171 :
1172 30902736 : const int id = k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1173 30902736 : grad_f[id + 0 * offset] = drdx * dpdr + dsdx * dpds + dtdx * dpdt;
1174 30902736 : grad_f[id + 1 * offset] = drdy * dpdr + dsdy * dpds + dtdy * dpdt;
1175 30902736 : grad_f[id + 2 * offset] = drdz * dpdr + dsdz * dpds + dtdz * dpdt;
1176 : }
1177 : }
1178 : }
1179 223308 : }
1180 :
1181 : bool
1182 219 : isHeatFluxBoundary(const int boundary)
1183 : {
1184 : // the heat flux boundary is now named 'codedFixedGradient', but 'fixedGradient'
1185 : // will be present for backwards compatibility
1186 876 : return (bcMap::text(boundary, "scalar00") == "fixedGradient") ||
1187 1095 : (bcMap::text(boundary, "scalar00") == "codedFixedGradient");
1188 : }
1189 :
1190 : bool
1191 17 : isMovingMeshBoundary(const int boundary)
1192 : {
1193 34 : return bcMap::text(boundary, "mesh") == "codedFixedValue";
1194 : }
1195 :
1196 : bool
1197 0 : isTemperatureBoundary(const int boundary)
1198 : {
1199 0 : return bcMap::text(boundary, "scalar00") == "fixedValue";
1200 : }
1201 :
1202 : const std::string
1203 1 : temperatureBoundaryType(const int boundary)
1204 : {
1205 2 : return bcMap::text(boundary, "scalar00");
1206 : }
1207 :
1208 : int
1209 1603 : polynomialOrder()
1210 : {
1211 1603 : return entireMesh()->N;
1212 : }
1213 :
1214 : int
1215 806 : Nelements()
1216 : {
1217 806 : int n_local = entireMesh()->Nelements;
1218 : int n_global;
1219 806 : MPI_Allreduce(&n_local, &n_global, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm);
1220 806 : return n_global;
1221 : }
1222 :
1223 : int
1224 11866645 : Nfaces()
1225 : {
1226 11866645 : return entireMesh()->Nfaces;
1227 : }
1228 :
1229 : int
1230 807 : dim()
1231 : {
1232 807 : return entireMesh()->dim;
1233 : }
1234 :
1235 : int
1236 0 : NfaceVertices()
1237 : {
1238 0 : return entireMesh()->NfaceVertices;
1239 : }
1240 :
1241 : int
1242 806 : NboundaryFaces()
1243 : {
1244 806 : return entireMesh()->NboundaryFaces;
1245 : }
1246 :
1247 : int
1248 7221 : NboundaryID()
1249 : {
1250 7221 : if (entireMesh()->cht)
1251 182 : return nekData.NboundaryIDt;
1252 : else
1253 7039 : return nekData.NboundaryID;
1254 : }
1255 :
1256 : bool
1257 2553 : validBoundaryIDs(const std::vector<int> & boundary_id, int & first_invalid_id, int & n_boundaries)
1258 : {
1259 2553 : n_boundaries = NboundaryID();
1260 :
1261 : bool valid_boundary_ids = true;
1262 5933 : for (const auto & b : boundary_id)
1263 : {
1264 3380 : if ((b > n_boundaries) || (b <= 0))
1265 : {
1266 3 : first_invalid_id = b;
1267 : valid_boundary_ids = false;
1268 : }
1269 : }
1270 :
1271 2553 : return valid_boundary_ids;
1272 : }
1273 :
1274 : double
1275 35440 : scalar01(const int id)
1276 : {
1277 35440 : nrs_t * nrs = (nrs_t *)nrsPtr();
1278 35440 : return nrs->cds->S[id + 1 * scalarFieldOffset()];
1279 : }
1280 :
1281 : double
1282 61552 : scalar02(const int id)
1283 : {
1284 61552 : nrs_t * nrs = (nrs_t *)nrsPtr();
1285 61552 : return nrs->cds->S[id + 2 * scalarFieldOffset()];
1286 : }
1287 :
1288 : double
1289 26224 : scalar03(const int id)
1290 : {
1291 26224 : nrs_t * nrs = (nrs_t *)nrsPtr();
1292 26224 : return nrs->cds->S[id + 3 * scalarFieldOffset()];
1293 : }
1294 :
1295 : double
1296 505472 : usrwrk00(const int id)
1297 : {
1298 505472 : nrs_t * nrs = (nrs_t *)nrsPtr();
1299 505472 : return nrs->usrwrk[id];
1300 : }
1301 :
1302 : double
1303 20672 : usrwrk01(const int id)
1304 : {
1305 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1306 20672 : return nrs->usrwrk[id + nrs->fieldOffset];
1307 : }
1308 :
1309 : double
1310 20672 : usrwrk02(const int id)
1311 : {
1312 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1313 20672 : return nrs->usrwrk[id + 2 * nrs->fieldOffset];
1314 : }
1315 :
1316 : double
1317 1328356616 : temperature(const int id)
1318 : {
1319 1328356616 : nrs_t * nrs = (nrs_t *)nrsPtr();
1320 1328356616 : return nrs->cds->S[id];
1321 : }
1322 :
1323 : double
1324 223118976 : pressure(const int id)
1325 : {
1326 223118976 : nrs_t * nrs = (nrs_t *)nrsPtr();
1327 223118976 : return nrs->P[id];
1328 : }
1329 :
1330 : double
1331 86179604 : unity(const int /* id */)
1332 : {
1333 86179604 : return 1.0;
1334 : }
1335 :
1336 : double
1337 209818528 : velocity_x(const int id)
1338 : {
1339 209818528 : nrs_t * nrs = (nrs_t *)nrsPtr();
1340 209818528 : return nrs->U[id + 0 * nrs->fieldOffset];
1341 : }
1342 :
1343 : double
1344 147149728 : velocity_y(const int id)
1345 : {
1346 147149728 : nrs_t * nrs = (nrs_t *)nrsPtr();
1347 147149728 : return nrs->U[id + 1 * nrs->fieldOffset];
1348 : }
1349 :
1350 : double
1351 208312272 : velocity_z(const int id)
1352 : {
1353 208312272 : nrs_t * nrs = (nrs_t *)nrsPtr();
1354 208312272 : return nrs->U[id + 2 * nrs->fieldOffset];
1355 : }
1356 :
1357 : double
1358 2027844 : velocity(const int id)
1359 : {
1360 2027844 : nrs_t * nrs = (nrs_t *)nrsPtr();
1361 2027844 : int offset = nrs->fieldOffset;
1362 :
1363 2027844 : return std::sqrt(nrs->U[id + 0 * offset] * nrs->U[id + 0 * offset] +
1364 2027844 : nrs->U[id + 1 * offset] * nrs->U[id + 1 * offset] +
1365 2027844 : nrs->U[id + 2 * offset] * nrs->U[id + 2 * offset]);
1366 : }
1367 :
1368 : double
1369 47744 : velocity_x_squared(const int id)
1370 : {
1371 47744 : return std::pow(velocity_x(id), 2);
1372 : }
1373 :
1374 : double
1375 47744 : velocity_y_squared(const int id)
1376 : {
1377 47744 : return std::pow(velocity_y(id), 2);
1378 : }
1379 :
1380 : double
1381 52912 : velocity_z_squared(const int id)
1382 : {
1383 52912 : return std::pow(velocity_z(id), 2);
1384 : }
1385 :
1386 : void
1387 64817272 : flux(const int id, const dfloat value)
1388 : {
1389 64817272 : nrs_t * nrs = (nrs_t *)nrsPtr();
1390 64817272 : nrs->usrwrk[indices.flux + id] = value;
1391 64817272 : }
1392 :
1393 : void
1394 118273024 : heat_source(const int id, const dfloat value)
1395 : {
1396 118273024 : nrs_t * nrs = (nrs_t *)nrsPtr();
1397 118273024 : nrs->usrwrk[indices.heat_source + id] = value;
1398 118273024 : }
1399 :
1400 : void
1401 27648 : x_displacement(const int id, const dfloat value)
1402 : {
1403 27648 : mesh_t * mesh = entireMesh();
1404 27648 : mesh->x[id] = value;
1405 27648 : }
1406 :
1407 : void
1408 27648 : y_displacement(const int id, const dfloat value)
1409 : {
1410 27648 : mesh_t * mesh = entireMesh();
1411 27648 : mesh->y[id] = value;
1412 27648 : }
1413 :
1414 : void
1415 27648 : z_displacement(const int id, const dfloat value)
1416 : {
1417 27648 : mesh_t * mesh = entireMesh();
1418 27648 : mesh->z[id] = value;
1419 27648 : }
1420 :
1421 : void
1422 4032000 : mesh_velocity_x(const int id, const dfloat value)
1423 : {
1424 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1425 4032000 : nrs->usrwrk[indices.mesh_velocity_x + id] = value;
1426 4032000 : }
1427 :
1428 : void
1429 4032000 : mesh_velocity_y(const int id, const dfloat value)
1430 : {
1431 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1432 4032000 : nrs->usrwrk[indices.mesh_velocity_y + id] = value;
1433 4032000 : }
1434 :
1435 : void
1436 4032000 : mesh_velocity_z(const int id, const dfloat value)
1437 : {
1438 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1439 4032000 : nrs->usrwrk[indices.mesh_velocity_z + id] = value;
1440 4032000 : }
1441 :
1442 : void
1443 195040 : checkFieldValidity(const field::NekFieldEnum & field)
1444 : {
1445 : // by placing this check here, as opposed to inside the NekFieldInterface,
1446 : // we can also leverage this error checking for the 'outputs' of NekRSProblem,
1447 : // which does not inherit from NekFieldInterface but still accesses the solutionPointers.
1448 : // If this gets moved elsewhere, need to be sure to add dedicated testing for
1449 : // the 'outputs' on NekRSProblem.
1450 :
1451 : // TODO: would be nice for NekRSProblem to only access field information via the
1452 : // NekFieldInterface; refactor later
1453 :
1454 195040 : switch (field)
1455 : {
1456 87372 : case field::temperature:
1457 87372 : if (!hasTemperatureVariable())
1458 1 : mooseError("Cannot find 'temperature' "
1459 : "because your Nek case files do not have a temperature variable!");
1460 : break;
1461 85 : case field::scalar01:
1462 85 : if (!hasScalarVariable(1))
1463 1 : mooseError("Cannot find 'scalar01' "
1464 : "because your Nek case files do not have a scalar01 variable!");
1465 : break;
1466 201 : case field::scalar02:
1467 201 : if (!hasScalarVariable(2))
1468 1 : mooseError("Cannot find 'scalar02' "
1469 : "because your Nek case files do not have a scalar02 variable!");
1470 : break;
1471 57 : case field::scalar03:
1472 57 : if (!hasScalarVariable(3))
1473 1 : mooseError("Cannot find 'scalar03' "
1474 : "because your Nek case files do not have a scalar03 variable!");
1475 : break;
1476 454 : case field::usrwrk00:
1477 454 : if (n_usrwrk_slots < 1)
1478 2 : mooseError("Cannot find 'usrwrk00' because you have only allocated 'n_usrwrk_slots = " +
1479 1 : std::to_string(n_usrwrk_slots) + "'");
1480 : break;
1481 46 : case field::usrwrk01:
1482 46 : if (n_usrwrk_slots < 2)
1483 2 : mooseError("Cannot find 'usrwrk01' because you have only allocated 'n_usrwrk_slots = " +
1484 1 : std::to_string(n_usrwrk_slots) + "'");
1485 : break;
1486 46 : case field::usrwrk02:
1487 46 : if (n_usrwrk_slots < 3)
1488 2 : mooseError("Cannot find 'usrwrk02' because you have only allocated 'n_usrwrk_slots = " +
1489 1 : std::to_string(n_usrwrk_slots) + "'");
1490 : break;
1491 : }
1492 195033 : }
1493 :
1494 191592 : double (*solutionPointer(const field::NekFieldEnum & field))(int)
1495 : {
1496 : // we include this here as well, in addition to within the NekFieldInterface, because
1497 : // the NekRSProblem accesses these methods without inheriting from NekFieldInterface
1498 191592 : checkFieldValidity(field);
1499 :
1500 : double (*f)(int);
1501 :
1502 191592 : switch (field)
1503 : {
1504 : case field::velocity_x:
1505 : f = &velocity_x;
1506 : break;
1507 17892 : case field::velocity_y:
1508 : f = &velocity_y;
1509 17892 : break;
1510 14984 : case field::velocity_z:
1511 : f = &velocity_z;
1512 14984 : break;
1513 196 : case field::velocity:
1514 : f = &velocity;
1515 196 : break;
1516 0 : case field::velocity_component:
1517 0 : mooseError("The 'velocity_component' field is not compatible with the solutionPointer "
1518 : "interface!");
1519 : break;
1520 48 : case field::velocity_x_squared:
1521 : f = &velocity_x_squared;
1522 48 : break;
1523 48 : case field::velocity_y_squared:
1524 : f = &velocity_y_squared;
1525 48 : break;
1526 52 : case field::velocity_z_squared:
1527 : f = &velocity_z_squared;
1528 52 : break;
1529 86120 : case field::temperature:
1530 : f = &temperature;
1531 86120 : break;
1532 31004 : case field::pressure:
1533 : f = &pressure;
1534 31004 : break;
1535 48 : case field::scalar01:
1536 : f = &scalar01;
1537 48 : break;
1538 176 : case field::scalar02:
1539 : f = &scalar02;
1540 176 : break;
1541 32 : case field::scalar03:
1542 : f = &scalar03;
1543 32 : break;
1544 9148 : case field::unity:
1545 : f = &unity;
1546 9148 : break;
1547 420 : case field::usrwrk00:
1548 : f = &usrwrk00;
1549 420 : break;
1550 16 : case field::usrwrk01:
1551 : f = &usrwrk01;
1552 16 : break;
1553 16 : case field::usrwrk02:
1554 : f = &usrwrk02;
1555 16 : break;
1556 0 : default:
1557 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1558 : }
1559 :
1560 191592 : return f;
1561 : }
1562 :
1563 1862524 : void (*solutionPointer(const field::NekWriteEnum & field))(int, dfloat)
1564 : {
1565 : void (*f)(int, dfloat);
1566 :
1567 1862524 : switch (field)
1568 : {
1569 : case field::flux:
1570 : f = &flux;
1571 : break;
1572 396748 : case field::heat_source:
1573 : f = &heat_source;
1574 396748 : break;
1575 1024 : case field::x_displacement:
1576 : f = &x_displacement;
1577 1024 : break;
1578 1024 : case field::y_displacement:
1579 : f = &y_displacement;
1580 1024 : break;
1581 1024 : case field::z_displacement:
1582 : f = &z_displacement;
1583 1024 : break;
1584 179200 : case field::mesh_velocity_x:
1585 : f = &mesh_velocity_x;
1586 179200 : break;
1587 179200 : case field::mesh_velocity_y:
1588 : f = &mesh_velocity_y;
1589 179200 : break;
1590 179200 : case field::mesh_velocity_z:
1591 : f = &mesh_velocity_z;
1592 179200 : break;
1593 0 : default:
1594 0 : throw std::runtime_error("Unhandled NekWriteEnum!");
1595 : }
1596 :
1597 1862524 : return f;
1598 : }
1599 :
1600 : void
1601 115 : initializeDimensionalScales(const double U,
1602 : const double T,
1603 : const double dT,
1604 : const double L,
1605 : const double rho,
1606 : const double Cp,
1607 : const double s01,
1608 : const double ds01,
1609 : const double s02,
1610 : const double ds02,
1611 : const double s03,
1612 : const double ds03)
1613 : {
1614 115 : scales.U_ref = U;
1615 115 : scales.T_ref = T;
1616 115 : scales.dT_ref = dT;
1617 115 : scales.L_ref = L;
1618 115 : scales.A_ref = L * L;
1619 115 : scales.V_ref = L * L * L;
1620 115 : scales.rho_ref = rho;
1621 115 : scales.Cp_ref = Cp;
1622 115 : scales.t_ref = L / U;
1623 115 : scales.P_ref = rho * U * U;
1624 :
1625 115 : scales.s01_ref = s01;
1626 115 : scales.ds01_ref = ds01;
1627 115 : scales.s02_ref = s02;
1628 115 : scales.ds02_ref = ds02;
1629 115 : scales.s03_ref = s03;
1630 115 : scales.ds03_ref = ds03;
1631 :
1632 115 : scales.flux_ref = rho * U * Cp * dT;
1633 115 : scales.source_ref = scales.flux_ref / L;
1634 115 : }
1635 :
1636 : double
1637 457 : referenceLength()
1638 : {
1639 457 : return scales.L_ref;
1640 : }
1641 :
1642 : double
1643 271566 : referenceTime()
1644 : {
1645 271566 : return scales.t_ref;
1646 : }
1647 :
1648 : double
1649 246 : referenceArea()
1650 : {
1651 246 : return scales.A_ref;
1652 : }
1653 :
1654 : double
1655 1132 : referenceVolume()
1656 : {
1657 1132 : return scales.V_ref;
1658 : }
1659 :
1660 : Real
1661 66768244 : nondimensionalAdditive(const field::NekFieldEnum & field)
1662 : {
1663 66768244 : switch (field)
1664 : {
1665 24750048 : case field::temperature:
1666 24750048 : return scales.T_ref;
1667 8660 : case field::scalar01:
1668 8660 : return scales.s01_ref;
1669 40904 : case field::scalar02:
1670 40904 : return scales.s02_ref;
1671 5576 : case field::scalar03:
1672 5576 : return scales.s03_ref;
1673 : default:
1674 : return 0;
1675 : }
1676 : }
1677 :
1678 : Real
1679 15350 : nondimensionalAdditive(const field::NekWriteEnum & field)
1680 : {
1681 15350 : switch (field)
1682 : {
1683 15350 : case field::flux:
1684 : case field::heat_source:
1685 : case field::x_displacement:
1686 : case field::y_displacement:
1687 : case field::z_displacement:
1688 : case field::mesh_velocity_x:
1689 : case field::mesh_velocity_y:
1690 : case field::mesh_velocity_z:
1691 15350 : return 0.0;
1692 0 : default:
1693 0 : mooseError("Unhandled NekWriteEnum in nondimensionalAdditive!");
1694 : }
1695 : }
1696 :
1697 : Real
1698 44752 : nondimensionalDivisor(const field::NekWriteEnum & field)
1699 : {
1700 44752 : switch (field)
1701 : {
1702 40341 : case field::flux:
1703 40341 : return scales.flux_ref;
1704 3507 : case field::heat_source:
1705 3507 : return scales.source_ref;
1706 904 : case field::x_displacement:
1707 : case field::y_displacement:
1708 : case field::z_displacement:
1709 904 : return scales.L_ref;
1710 0 : case field::mesh_velocity_x:
1711 : case field::mesh_velocity_y:
1712 : case field::mesh_velocity_z:
1713 0 : return scales.U_ref;
1714 0 : default:
1715 0 : mooseError("Unhandled NekWriteEnum in nondimensionalDivisor!");
1716 : }
1717 : }
1718 :
1719 : Real
1720 67305958 : nondimensionalDivisor(const field::NekFieldEnum & field)
1721 : {
1722 67305958 : switch (field)
1723 : {
1724 32127380 : case field::velocity_x:
1725 : case field::velocity_y:
1726 : case field::velocity_z:
1727 : case field::velocity:
1728 : case field::velocity_component:
1729 32127380 : return scales.U_ref;
1730 5336 : case field::velocity_x_squared:
1731 : case field::velocity_y_squared:
1732 : case field::velocity_z_squared:
1733 5336 : return scales.U_ref * scales.U_ref;
1734 24750048 : case field::temperature:
1735 24750048 : return scales.dT_ref;
1736 10342943 : case field::pressure:
1737 10342943 : return scales.P_ref;
1738 8660 : case field::scalar01:
1739 8660 : return scales.ds01_ref;
1740 40904 : case field::scalar02:
1741 40904 : return scales.ds02_ref;
1742 5576 : case field::scalar03:
1743 5576 : return scales.ds03_ref;
1744 : case field::unity:
1745 : // no dimensionalization needed
1746 : return 1.0;
1747 445 : case field::usrwrk00:
1748 445 : return scratchUnits(0);
1749 29 : case field::usrwrk01:
1750 29 : return scratchUnits(1);
1751 29 : case field::usrwrk02:
1752 29 : return scratchUnits(2);
1753 0 : default:
1754 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1755 : }
1756 : }
1757 :
1758 : Real
1759 503 : scratchUnits(const int slot)
1760 : {
1761 503 : if (indices.flux != -1 && slot == indices.flux / nekrs::fieldOffset())
1762 424 : return scales.flux_ref;
1763 79 : else if (indices.heat_source != -1 && slot == indices.heat_source / nekrs::fieldOffset())
1764 4 : return scales.source_ref;
1765 75 : else if (is_nondimensional)
1766 : {
1767 : // TODO: we are lazy and did not include all the usrwrk indices
1768 66 : mooseDoOnce(mooseWarning(
1769 : "The units of 'usrwrk0" + std::to_string(slot) +
1770 : "' are unknown, so we cannot dimensionalize any objects using 'field = usrwrk0" +
1771 : std::to_string(slot) +
1772 : "'. The output for this quantity will be given in non-dimensional form.\n\nYou will need "
1773 : "to manipulate the data manually from Cardinal if you need to dimensionalize it."));
1774 : }
1775 :
1776 : return 1.0;
1777 : }
1778 :
1779 : void
1780 799 : nondimensional(const bool n)
1781 : {
1782 799 : is_nondimensional = n;
1783 799 : }
1784 :
1785 : template <>
1786 : MPI_Datatype
1787 182672 : resolveType<double>()
1788 : {
1789 182672 : return MPI_DOUBLE;
1790 : }
1791 :
1792 : template <>
1793 : MPI_Datatype
1794 7024 : resolveType<int>()
1795 : {
1796 7024 : return MPI_INT;
1797 : }
1798 :
1799 : } // end namespace nekrs
1800 :
1801 : #endif
|
