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 365 : setAbsoluteTol(double tol)
41 : {
42 365 : abs_tol = tol;
43 365 : }
44 :
45 : void
46 365 : setRelativeTol(double tol)
47 : {
48 365 : rel_tol = tol;
49 365 : }
50 :
51 : void
52 836 : setNekSetupTime(const double & time)
53 : {
54 836 : setup_time = time;
55 836 : }
56 :
57 : double
58 825 : getNekSetupTime()
59 : {
60 825 : return setup_time;
61 : }
62 :
63 : void
64 766 : setStartTime(const double & start)
65 : {
66 1532 : platform->options.setArgs("START TIME", to_string_f(start));
67 766 : }
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 838 : buildOnly(int buildOnly)
102 : {
103 838 : build_only = buildOnly;
104 838 : }
105 :
106 : int
107 104131 : buildOnly()
108 : {
109 104131 : return build_only;
110 : }
111 :
112 : bool
113 0 : hasCHT()
114 : {
115 0 : return entireMesh()->cht;
116 : }
117 :
118 : bool
119 36347 : hasMovingMesh()
120 : {
121 36347 : return platform->options.compareArgs("MOVING MESH", "TRUE");
122 : }
123 :
124 : bool
125 32834 : hasVariableDt()
126 : {
127 32834 : 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 1572 : endControlElapsedTime()
144 : {
145 3144 : return !platform->options.getArgs("STOP AT ELAPSED TIME").empty();
146 : }
147 :
148 : bool
149 1572 : endControlTime()
150 : {
151 1572 : return endTime() > 0;
152 : }
153 :
154 : bool
155 786 : endControlNumSteps()
156 : {
157 786 : return !endControlElapsedTime() && !endControlTime();
158 : }
159 :
160 : bool
161 262393902 : hasTemperatureVariable()
162 : {
163 262393902 : nrs_t * nrs = (nrs_t *)nrsPtr();
164 786731650 : return nrs->Nscalar ? platform->options.compareArgs("SCALAR00 IS TEMPERATURE", "TRUE") : false;
165 : }
166 :
167 : bool
168 679 : hasTemperatureSolve()
169 : {
170 679 : nrs_t * nrs = (nrs_t *)nrsPtr();
171 679 : return hasTemperatureVariable() ? nrs->cds->compute[0] : false;
172 : }
173 :
174 : bool
175 1207 : hasScalarVariable(int scalarId)
176 : {
177 1207 : nrs_t * nrs = (nrs_t *)nrsPtr();
178 1207 : 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 834 : isInitialized()
189 : {
190 834 : nrs_t * nrs = (nrs_t *)nrsPtr();
191 834 : return nrs;
192 : }
193 :
194 : int
195 143813 : scalarFieldOffset()
196 : {
197 143813 : nrs_t * nrs = (nrs_t *)nrsPtr();
198 143813 : return nrs->cds->fieldOffset[0];
199 : }
200 :
201 : int
202 4358465 : velocityFieldOffset()
203 : {
204 4358465 : nrs_t * nrs = (nrs_t *)nrsPtr();
205 4358465 : return nrs->fieldOffset;
206 : }
207 :
208 : int
209 20582 : fieldOffset()
210 : {
211 20582 : if (hasTemperatureVariable())
212 20565 : return scalarFieldOffset();
213 : else
214 17 : return velocityFieldOffset();
215 : }
216 :
217 : mesh_t *
218 262284079 : entireMesh()
219 : {
220 262284079 : if (hasTemperatureVariable())
221 172403728 : return temperatureMesh();
222 : else
223 89880351 : return flowMesh();
224 : }
225 :
226 : mesh_t *
227 90306116 : flowMesh()
228 : {
229 90306116 : nrs_t * nrs = (nrs_t *)nrsPtr();
230 90306116 : return nrs->meshV;
231 : }
232 :
233 : mesh_t *
234 173240717 : temperatureMesh()
235 : {
236 173240717 : nrs_t * nrs = (nrs_t *)nrsPtr();
237 173240717 : return nrs->cds->mesh[0];
238 : }
239 :
240 : mesh_t *
241 408728 : getMesh(const nek_mesh::NekMeshEnum pp_mesh)
242 : {
243 408728 : if (pp_mesh == nek_mesh::fluid)
244 1328 : return flowMesh();
245 407400 : else if (pp_mesh == nek_mesh::all)
246 407399 : 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 13154032 : commRank()
254 : {
255 13154032 : return platform->comm.mpiRank;
256 : }
257 :
258 : int
259 1072797 : commSize()
260 : {
261 1072797 : return platform->comm.mpiCommSize;
262 : }
263 :
264 : bool
265 836 : scratchAvailable()
266 : {
267 836 : 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 836 : if (nrs->usrwrk)
277 1 : return false;
278 :
279 : return true;
280 : }
281 :
282 : void
283 835 : initializeScratch(const unsigned int & n_slots)
284 : {
285 835 : if (n_slots == 0)
286 : return;
287 :
288 439 : nrs_t * nrs = (nrs_t *)nrsPtr();
289 439 : mesh_t * mesh = entireMesh();
290 :
291 : // clear them just to be sure
292 439 : 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 439 : nrs->usrwrk = (double *)calloc(n_slots * fieldOffset(), sizeof(double));
298 439 : nrs->o_usrwrk = platform->device.malloc(n_slots * fieldOffset() * sizeof(double), nrs->usrwrk);
299 :
300 439 : n_usrwrk_slots = n_slots;
301 : }
302 :
303 : void
304 1197 : freeScratch()
305 : {
306 1197 : nrs_t * nrs = (nrs_t *)nrsPtr();
307 :
308 1197 : freePointer(nrs->usrwrk);
309 1197 : nrs->o_usrwrk.free();
310 1197 : }
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 1646 : interpolationMatrix(double * I, int starting_points, int ending_points)
342 : {
343 1646 : DegreeRaiseMatrix1D(starting_points - 1, ending_points - 1, I);
344 1646 : }
345 :
346 : void
347 2990262 : interpolateVolumeHex3D(const double * I, double * x, int N, double * Ix, int M)
348 : {
349 2990262 : double * Ix1 = (dfloat *)calloc(N * N * M, sizeof(double));
350 2990262 : double * Ix2 = (dfloat *)calloc(N * M * M, sizeof(double));
351 :
352 17190894 : for (int k = 0; k < N; ++k)
353 99730324 : for (int j = 0; j < N; ++j)
354 352395804 : for (int i = 0; i < M; ++i)
355 : {
356 : dfloat tmp = 0;
357 2077270560 : for (int n = 0; n < N; ++n)
358 1810404448 : tmp += I[i * N + n] * x[k * N * N + j * N + n];
359 266866112 : Ix1[k * N * M + j * M + i] = tmp;
360 : }
361 :
362 17190894 : for (int k = 0; k < N; ++k)
363 61027096 : for (int j = 0; j < M; ++j)
364 217922840 : for (int i = 0; i < M; ++i)
365 : {
366 : dfloat tmp = 0;
367 1045391880 : for (int n = 0; n < N; ++n)
368 874295504 : tmp += I[j * N + n] * Ix1[k * N * M + n * M + i];
369 171096376 : Ix2[k * M * M + j * M + i] = tmp;
370 : }
371 :
372 13768682 : for (int k = 0; k < M; ++k)
373 56344972 : for (int j = 0; j < M; ++j)
374 274845208 : for (int i = 0; i < M; ++i)
375 : {
376 : dfloat tmp = 0;
377 957360984 : for (int n = 0; n < N; ++n)
378 728082328 : tmp += I[k * N + n] * Ix2[n * M * M + j * M + i];
379 229278656 : Ix[k * M * M + j * M + i] = tmp;
380 : }
381 :
382 : freePointer(Ix1);
383 : freePointer(Ix2);
384 2990262 : }
385 :
386 : void
387 493962 : interpolateSurfaceFaceHex3D(
388 : double * scratch, const double * I, double * x, int N, double * Ix, int M)
389 : {
390 1955102 : for (int j = 0; j < N; ++j)
391 6689820 : for (int i = 0; i < M; ++i)
392 : {
393 : double tmp = 0;
394 24319128 : for (int n = 0; n < N; ++n)
395 : {
396 19090448 : tmp += I[i * N + n] * x[j * N + n];
397 : }
398 5228680 : scratch[j * M + i] = tmp;
399 : }
400 :
401 2398478 : for (int j = 0; j < M; ++j)
402 9681636 : for (int i = 0; i < M; ++i)
403 : {
404 : double tmp = 0;
405 27590304 : for (int n = 0; n < N; ++n)
406 : {
407 19813184 : tmp += I[j * N + n] * scratch[n * M + i];
408 : }
409 7777120 : Ix[j * M + i] = tmp;
410 : }
411 493962 : }
412 :
413 : void
414 95713 : displacementAndCounts(const std::vector<int> & base_counts,
415 : int * counts,
416 : int * displacement,
417 : const int multiplier = 1.0)
418 : {
419 487634 : for (int i = 0; i < commSize(); ++i)
420 391921 : counts[i] = base_counts[i] * multiplier;
421 :
422 95713 : displacement[0] = 0;
423 391921 : for (int i = 1; i < commSize(); i++)
424 296208 : displacement[i] = displacement[i - 1] + counts[i - 1];
425 95713 : }
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 833 : initializeHostMeshParameters()
551 : {
552 833 : mesh_t * mesh = entireMesh();
553 833 : sgeo = (dfloat *)calloc(mesh->o_sgeo.size(), sizeof(dfloat));
554 833 : vgeo = (dfloat *)calloc(mesh->o_vgeo.size(), sizeof(dfloat));
555 833 : }
556 :
557 : void
558 2337 : updateHostMeshParameters()
559 : {
560 2337 : mesh_t * mesh = entireMesh();
561 2337 : mesh->o_sgeo.copyTo(sgeo);
562 2337 : mesh->o_vgeo.copyTo(vgeo);
563 2337 : }
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, 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 6383351 : value = std::max(value, f(mesh->vmapM[offset + v], 0 /* unused */));
601 : else
602 6369628 : value = std::min(value, f(mesh->vmapM[offset + v], 0 /* unused */));
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, 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 662078764 : value = std::max(value, f(i * mesh->Np + j, 0 /* unused */));
655 : else
656 496912612 : value = std::min(value, f(i * mesh->Np + j, 0 /* unused */));
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 32192 : 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 32192 : integral *= nondimensionalDivisor(integrand);
818 :
819 : // scale the boundary integral
820 32192 : 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 32192 : auto add = nondimensionalAdditive(integrand);
825 32192 : if (std::abs(add) > 1e-8)
826 360 : integral += add * area(boundary_id, pp_mesh);
827 32192 : }
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, 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, 0 /* unused */) * 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 13740 : area(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
859 : {
860 13740 : mesh_t * mesh = getMesh(pp_mesh);
861 :
862 13740 : double integral = 0.0;
863 :
864 2569002 : for (int i = 0; i < mesh->Nelements; ++i)
865 : {
866 17886834 : for (int j = 0; j < mesh->Nfaces; ++j)
867 : {
868 15331572 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
869 :
870 15331572 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
871 : {
872 416204 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
873 8771664 : for (int v = 0; v < mesh->Nfp; ++v)
874 : {
875 8355460 : integral += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
876 : }
877 : }
878 : }
879 : }
880 :
881 : // sum across all processes
882 : double total_integral;
883 13740 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
884 :
885 13740 : dimensionalizeSideIntegral(field::unity, boundary_id, total_integral, pp_mesh);
886 :
887 13740 : return total_integral;
888 : }
889 :
890 : double
891 18417 : sideIntegral(const std::vector<int> & boundary_id, const field::NekFieldEnum & integrand,
892 : const nek_mesh::NekMeshEnum pp_mesh)
893 : {
894 18417 : mesh_t * mesh = getMesh(pp_mesh);
895 :
896 18416 : double integral = 0.0;
897 :
898 : double (*f)(int, int);
899 18416 : f = solutionPointer(integrand);
900 :
901 3062782 : for (int i = 0; i < mesh->Nelements; ++i)
902 : {
903 21310562 : for (int j = 0; j < mesh->Nfaces; ++j)
904 : {
905 18266196 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
906 :
907 18266196 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
908 : {
909 556758 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
910 11475402 : for (int v = 0; v < mesh->Nfp; ++v)
911 : {
912 10918644 : integral +=
913 10918644 : f(mesh->vmapM[offset + v], 0 /* unused */) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
914 : }
915 : }
916 : }
917 : }
918 :
919 : // sum across all processes
920 : double total_integral;
921 18416 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
922 :
923 18416 : dimensionalizeSideIntegral(integrand, boundary_id, total_integral, pp_mesh);
924 :
925 18416 : return total_integral;
926 : }
927 :
928 : double
929 868 : massFlowrate(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
930 : {
931 868 : mesh_t * mesh = getMesh(pp_mesh);
932 868 : nrs_t * nrs = (nrs_t *)nrsPtr();
933 :
934 : // TODO: This function only works correctly if the density is constant, because
935 : // otherwise we need to copy the density from device to host
936 : double rho;
937 868 : platform->options.getArgs("DENSITY", rho);
938 :
939 868 : double integral = 0.0;
940 :
941 412012 : for (int i = 0; i < mesh->Nelements; ++i)
942 : {
943 2878008 : for (int j = 0; j < mesh->Nfaces; ++j)
944 : {
945 2466864 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
946 :
947 2466864 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
948 : {
949 21112 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
950 586564 : for (int v = 0; v < mesh->Nfp; ++v)
951 : {
952 565452 : int vol_id = mesh->vmapM[offset + v];
953 565452 : int surf_offset = mesh->Nsgeo * (offset + v);
954 :
955 : double normal_velocity =
956 565452 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
957 565452 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
958 565452 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
959 :
960 565452 : integral += rho * normal_velocity * sgeo[surf_offset + WSJID];
961 : }
962 : }
963 : }
964 : }
965 :
966 : // sum across all processes
967 : double total_integral;
968 868 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
969 :
970 : // dimensionalize the mass flux and area
971 868 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
972 :
973 868 : return total_integral;
974 : }
975 :
976 : double
977 1280 : sideMassFluxWeightedIntegral(const std::vector<int> & boundary_id,
978 : const field::NekFieldEnum & integrand,
979 : const nek_mesh::NekMeshEnum pp_mesh)
980 : {
981 1280 : mesh_t * mesh = getMesh(pp_mesh);
982 1280 : nrs_t * nrs = (nrs_t *)nrsPtr();
983 :
984 : // TODO: This function only works correctly if the density is constant, because
985 : // otherwise we need to copy the density from device to host
986 : double rho;
987 1280 : platform->options.getArgs("DENSITY", rho);
988 :
989 1280 : double integral = 0.0;
990 :
991 : double (*f)(int, int);
992 1280 : f = solutionPointer(integrand);
993 :
994 605936 : for (int i = 0; i < mesh->Nelements; ++i)
995 : {
996 4232592 : for (int j = 0; j < mesh->Nfaces; ++j)
997 : {
998 3627936 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
999 :
1000 3627936 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1001 : {
1002 27498 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1003 914862 : for (int v = 0; v < mesh->Nfp; ++v)
1004 : {
1005 887364 : int vol_id = mesh->vmapM[offset + v];
1006 887364 : int surf_offset = mesh->Nsgeo * (offset + v);
1007 : double normal_velocity =
1008 887364 : nrs->U[vol_id + 0 * velocityFieldOffset()] * sgeo[surf_offset + NXID] +
1009 887364 : nrs->U[vol_id + 1 * velocityFieldOffset()] * sgeo[surf_offset + NYID] +
1010 887364 : nrs->U[vol_id + 2 * velocityFieldOffset()] * sgeo[surf_offset + NZID];
1011 887364 : integral += f(vol_id, 0 /* unused */) * rho * normal_velocity * sgeo[surf_offset + WSJID];
1012 : }
1013 : }
1014 : }
1015 : }
1016 :
1017 : // sum across all processes
1018 : double total_integral;
1019 1280 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1020 :
1021 : // dimensionalize the field if needed
1022 1280 : total_integral *= nondimensionalDivisor(integrand);
1023 :
1024 : // dimensionalize the mass flux and area
1025 1280 : total_integral *= scales.rho_ref * scales.U_ref * scales.A_ref;
1026 :
1027 : // for quantities with a relative scaling, we need to add back the reference
1028 : // contribution to the mass flux integral; we need this form here to avoid an infinite
1029 : // recursive loop
1030 1280 : auto add = nondimensionalAdditive(integrand);
1031 1280 : if (std::abs(add) > 1e-8)
1032 376 : total_integral += add * massFlowrate(boundary_id, pp_mesh);
1033 :
1034 1280 : return total_integral;
1035 : }
1036 :
1037 : double
1038 36 : pressureSurfaceForce(const std::vector<int> & boundary_id, const Point & direction, const nek_mesh::NekMeshEnum pp_mesh)
1039 : {
1040 36 : mesh_t * mesh = getMesh(pp_mesh);
1041 36 : nrs_t * nrs = (nrs_t *)nrsPtr();
1042 :
1043 36 : double integral = 0.0;
1044 :
1045 20232 : for (int i = 0; i < mesh->Nelements; ++i)
1046 : {
1047 141372 : for (int j = 0; j < mesh->Nfaces; ++j)
1048 : {
1049 121176 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1050 :
1051 121176 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1052 : {
1053 12780 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1054 472860 : for (int v = 0; v < mesh->Nfp; ++v)
1055 : {
1056 460080 : int vol_id = mesh->vmapM[offset + v];
1057 460080 : int surf_offset = mesh->Nsgeo * (offset + v);
1058 :
1059 460080 : double p_normal = nrs->P[vol_id] * (sgeo[surf_offset + NXID] * direction(0) +
1060 460080 : sgeo[surf_offset + NYID] * direction(1) +
1061 460080 : sgeo[surf_offset + NZID] * direction(2));
1062 :
1063 460080 : integral += p_normal * sgeo[surf_offset + WSJID];
1064 : }
1065 : }
1066 : }
1067 : }
1068 :
1069 : // sum across all processes
1070 : double total_integral;
1071 36 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1072 :
1073 36 : dimensionalizeSideIntegral(field::pressure, boundary_id, total_integral, pp_mesh);
1074 :
1075 36 : return total_integral;
1076 : }
1077 :
1078 : double
1079 11640 : heatFluxIntegral(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh)
1080 : {
1081 11640 : mesh_t * mesh = getMesh(pp_mesh);
1082 11640 : nrs_t * nrs = (nrs_t *)nrsPtr();
1083 :
1084 : // TODO: This function only works correctly if the conductivity is constant, because
1085 : // otherwise we need to copy the conductivity from device to host
1086 : double k;
1087 11640 : platform->options.getArgs("SCALAR00 DIFFUSIVITY", k);
1088 :
1089 11640 : double integral = 0.0;
1090 11640 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1091 :
1092 2259536 : for (int i = 0; i < mesh->Nelements; ++i)
1093 : {
1094 15735272 : for (int j = 0; j < mesh->Nfaces; ++j)
1095 : {
1096 13487376 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
1097 :
1098 13487376 : if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
1099 : {
1100 : // some inefficiency if an element has more than one face on the sideset of interest,
1101 : // because we will recompute the gradient in the element more than one time - but this
1102 : // is of little practical interest because this will be a minority of cases.
1103 224316 : gradient(mesh->Np, i, nrs->cds->S, grad_T, pp_mesh);
1104 :
1105 224316 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
1106 5751276 : for (int v = 0; v < mesh->Nfp; ++v)
1107 : {
1108 5526960 : int vol_id = mesh->vmapM[offset + v] - i * mesh->Np;
1109 5526960 : int surf_offset = mesh->Nsgeo * (offset + v);
1110 :
1111 5526960 : double normal_grad_T = grad_T[vol_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1112 5526960 : grad_T[vol_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1113 5526960 : grad_T[vol_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1114 :
1115 5526960 : integral += -k * normal_grad_T * sgeo[surf_offset + WSJID];
1116 : }
1117 : }
1118 : }
1119 : }
1120 :
1121 : freePointer(grad_T);
1122 :
1123 : // sum across all processes
1124 : double total_integral;
1125 11640 : MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
1126 :
1127 : // multiply by the reference heat flux and an area factor to dimensionalize
1128 11640 : total_integral *= scales.flux_ref * scales.A_ref;
1129 :
1130 11640 : return total_integral;
1131 : }
1132 :
1133 : void
1134 228348 : gradient(const int offset,
1135 : const int e,
1136 : const double * f,
1137 : double * grad_f,
1138 : const nek_mesh::NekMeshEnum pp_mesh)
1139 : {
1140 228348 : mesh_t * mesh = getMesh(pp_mesh);
1141 :
1142 1313484 : for (int k = 0; k < mesh->Nq; ++k)
1143 : {
1144 6630144 : for (int j = 0; j < mesh->Nq; ++j)
1145 : {
1146 36508224 : for (int i = 0; i < mesh->Nq; ++i)
1147 : {
1148 30963216 : const int gid = e * mesh->Np * mesh->Nvgeo + k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1149 30963216 : const double drdx = vgeo[gid + RXID * mesh->Np];
1150 30963216 : const double drdy = vgeo[gid + RYID * mesh->Np];
1151 30963216 : const double drdz = vgeo[gid + RZID * mesh->Np];
1152 30963216 : const double dsdx = vgeo[gid + SXID * mesh->Np];
1153 30963216 : const double dsdy = vgeo[gid + SYID * mesh->Np];
1154 30963216 : const double dsdz = vgeo[gid + SZID * mesh->Np];
1155 30963216 : const double dtdx = vgeo[gid + TXID * mesh->Np];
1156 30963216 : const double dtdy = vgeo[gid + TYID * mesh->Np];
1157 30963216 : const double dtdz = vgeo[gid + TZID * mesh->Np];
1158 :
1159 : // compute 'r' and 's' derivatives of (q_m) at node n
1160 : double dpdr = 0.f, dpds = 0.f, dpdt = 0.f;
1161 :
1162 222737184 : for (int n = 0; n < mesh->Nq; ++n)
1163 : {
1164 191773968 : const double Dr = mesh->D[i * mesh->Nq + n];
1165 191773968 : const double Ds = mesh->D[j * mesh->Nq + n];
1166 191773968 : const double Dt = mesh->D[k * mesh->Nq + n];
1167 :
1168 191773968 : dpdr += Dr * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + j * mesh->Nq + n];
1169 191773968 : dpds += Ds * f[e * mesh->Np + k * mesh->Nq * mesh->Nq + n * mesh->Nq + i];
1170 191773968 : dpdt += Dt * f[e * mesh->Np + n * mesh->Nq * mesh->Nq + j * mesh->Nq + i];
1171 : }
1172 :
1173 30963216 : const int id = k * mesh->Nq * mesh->Nq + j * mesh->Nq + i;
1174 30963216 : grad_f[id + 0 * offset] = drdx * dpdr + dsdx * dpds + dtdx * dpdt;
1175 30963216 : grad_f[id + 1 * offset] = drdy * dpdr + dsdy * dpds + dtdy * dpdt;
1176 30963216 : grad_f[id + 2 * offset] = drdz * dpdr + dsdz * dpds + dtdz * dpdt;
1177 : }
1178 : }
1179 : }
1180 228348 : }
1181 :
1182 : bool
1183 219 : isHeatFluxBoundary(const int boundary)
1184 : {
1185 : // the heat flux boundary is now named 'codedFixedGradient', but 'fixedGradient'
1186 : // will be present for backwards compatibility
1187 876 : return (bcMap::text(boundary, "scalar00") == "fixedGradient") ||
1188 1095 : (bcMap::text(boundary, "scalar00") == "codedFixedGradient");
1189 : }
1190 :
1191 : bool
1192 17 : isMovingMeshBoundary(const int boundary)
1193 : {
1194 34 : return bcMap::text(boundary, "mesh") == "codedFixedValue";
1195 : }
1196 :
1197 : bool
1198 0 : isTemperatureBoundary(const int boundary)
1199 : {
1200 0 : return bcMap::text(boundary, "scalar00") == "fixedValue";
1201 : }
1202 :
1203 : const std::string
1204 1 : temperatureBoundaryType(const int boundary)
1205 : {
1206 2 : return bcMap::text(boundary, "scalar00");
1207 : }
1208 :
1209 : int
1210 1655 : polynomialOrder()
1211 : {
1212 1655 : return entireMesh()->N;
1213 : }
1214 :
1215 : int
1216 832 : Nelements()
1217 : {
1218 832 : int n_local = entireMesh()->Nelements;
1219 : int n_global;
1220 832 : MPI_Allreduce(&n_local, &n_global, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm);
1221 832 : return n_global;
1222 : }
1223 :
1224 : int
1225 11875745 : Nfaces()
1226 : {
1227 11875745 : return entireMesh()->Nfaces;
1228 : }
1229 :
1230 : int
1231 833 : dim()
1232 : {
1233 833 : return entireMesh()->dim;
1234 : }
1235 :
1236 : int
1237 0 : NfaceVertices()
1238 : {
1239 0 : return entireMesh()->NfaceVertices;
1240 : }
1241 :
1242 : int
1243 832 : NboundaryFaces()
1244 : {
1245 832 : return entireMesh()->NboundaryFaces;
1246 : }
1247 :
1248 : int
1249 7471 : NboundaryID()
1250 : {
1251 7471 : if (entireMesh()->cht)
1252 182 : return nekData.NboundaryIDt;
1253 : else
1254 7289 : return nekData.NboundaryID;
1255 : }
1256 :
1257 : bool
1258 2603 : validBoundaryIDs(const std::vector<int> & boundary_id, int & first_invalid_id, int & n_boundaries)
1259 : {
1260 2603 : n_boundaries = NboundaryID();
1261 :
1262 : bool valid_boundary_ids = true;
1263 6033 : for (const auto & b : boundary_id)
1264 : {
1265 3430 : if ((b > n_boundaries) || (b <= 0))
1266 : {
1267 3 : first_invalid_id = b;
1268 : valid_boundary_ids = false;
1269 : }
1270 : }
1271 :
1272 2603 : return valid_boundary_ids;
1273 : }
1274 :
1275 : double
1276 35440 : get_scalar01(const int id, const int surf_offset = 0)
1277 : {
1278 35440 : nrs_t * nrs = (nrs_t *)nrsPtr();
1279 35440 : return nrs->cds->S[id + 1 * scalarFieldOffset()];
1280 : }
1281 :
1282 : double
1283 61552 : get_scalar02(const int id, const int surf_offset = 0)
1284 : {
1285 61552 : nrs_t * nrs = (nrs_t *)nrsPtr();
1286 61552 : return nrs->cds->S[id + 2 * scalarFieldOffset()];
1287 : }
1288 :
1289 : double
1290 26224 : get_scalar03(const int id, const int surf_offset = 0)
1291 : {
1292 26224 : nrs_t * nrs = (nrs_t *)nrsPtr();
1293 26224 : return nrs->cds->S[id + 3 * scalarFieldOffset()];
1294 : }
1295 :
1296 : double
1297 505472 : get_usrwrk00(const int id, const int surf_offset = 0)
1298 : {
1299 505472 : nrs_t * nrs = (nrs_t *)nrsPtr();
1300 505472 : return nrs->usrwrk[id];
1301 : }
1302 :
1303 : double
1304 20672 : get_usrwrk01(const int id, const int surf_offset = 0)
1305 : {
1306 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1307 20672 : return nrs->usrwrk[id + nrs->fieldOffset];
1308 : }
1309 :
1310 : double
1311 20672 : get_usrwrk02(const int id, const int surf_offset = 0)
1312 : {
1313 20672 : nrs_t * nrs = (nrs_t *)nrsPtr();
1314 20672 : return nrs->usrwrk[id + 2 * nrs->fieldOffset];
1315 : }
1316 :
1317 : double
1318 1328376256 : get_temperature(const int id, const int surf_offset)
1319 : {
1320 1328376256 : nrs_t * nrs = (nrs_t *)nrsPtr();
1321 1328376256 : return nrs->cds->S[id];
1322 : }
1323 :
1324 : double
1325 4032 : get_flux(const int id, const int surf_offset)
1326 : {
1327 : // TODO: this function does not support non-constant thermal conductivity
1328 : double k;
1329 4032 : platform->options.getArgs("SCALAR00 DIFFUSIVITY", k);
1330 :
1331 4032 : nrs_t * nrs = (nrs_t *)nrsPtr();
1332 :
1333 : // this call of nek_mesh::all should be fine because flux is not a 'field' which can be
1334 : // provided to the postprocessors which have the option to operate only on part of the mesh
1335 4032 : auto mesh = getMesh(nek_mesh::all);
1336 4032 : int elem_id = id / mesh->Np;
1337 4032 : int vertex_id = id % mesh->Np;
1338 :
1339 : // This function is slightly inefficient, because we compute grad(T) for all nodes in
1340 : // an element even though we only call this function for one node at a time
1341 4032 : double * grad_T = (double *)calloc(3 * mesh->Np, sizeof(double));
1342 4032 : gradient(mesh->Np, elem_id, nrs->cds->S, grad_T, nek_mesh::all);
1343 :
1344 4032 : double normal_grad_T = grad_T[vertex_id + 0 * mesh->Np] * sgeo[surf_offset + NXID] +
1345 4032 : grad_T[vertex_id + 1 * mesh->Np] * sgeo[surf_offset + NYID] +
1346 4032 : grad_T[vertex_id + 2 * mesh->Np] * sgeo[surf_offset + NZID];
1347 : freePointer(grad_T);
1348 :
1349 4032 : return -k * normal_grad_T;
1350 : }
1351 :
1352 : double
1353 223118976 : get_pressure(const int id, const int surf_offset)
1354 : {
1355 223118976 : nrs_t * nrs = (nrs_t *)nrsPtr();
1356 223118976 : return nrs->P[id];
1357 : }
1358 :
1359 : double
1360 86179604 : get_unity(const int /* id */, const int surf_offset)
1361 : {
1362 86179604 : return 1.0;
1363 : }
1364 :
1365 : double
1366 209818528 : get_velocity_x(const int id, const int surf_offset)
1367 : {
1368 209818528 : nrs_t * nrs = (nrs_t *)nrsPtr();
1369 209818528 : return nrs->U[id + 0 * nrs->fieldOffset];
1370 : }
1371 :
1372 : double
1373 147149728 : get_velocity_y(const int id, const int surf_offset)
1374 : {
1375 147149728 : nrs_t * nrs = (nrs_t *)nrsPtr();
1376 147149728 : return nrs->U[id + 1 * nrs->fieldOffset];
1377 : }
1378 :
1379 : double
1380 208312272 : get_velocity_z(const int id, const int surf_offset)
1381 : {
1382 208312272 : nrs_t * nrs = (nrs_t *)nrsPtr();
1383 208312272 : return nrs->U[id + 2 * nrs->fieldOffset];
1384 : }
1385 :
1386 : double
1387 2027844 : get_velocity(const int id, const int surf_offset)
1388 : {
1389 2027844 : nrs_t * nrs = (nrs_t *)nrsPtr();
1390 2027844 : int offset = nrs->fieldOffset;
1391 :
1392 2027844 : return std::sqrt(nrs->U[id + 0 * offset] * nrs->U[id + 0 * offset] +
1393 2027844 : nrs->U[id + 1 * offset] * nrs->U[id + 1 * offset] +
1394 2027844 : nrs->U[id + 2 * offset] * nrs->U[id + 2 * offset]);
1395 : }
1396 :
1397 : double
1398 47744 : get_velocity_x_squared(const int id, const int surf_offset)
1399 : {
1400 47744 : return std::pow(get_velocity_x(id, surf_offset), 2);
1401 : }
1402 :
1403 : double
1404 47744 : get_velocity_y_squared(const int id, const int surf_offset)
1405 : {
1406 47744 : return std::pow(get_velocity_y(id, surf_offset), 2);
1407 : }
1408 :
1409 : double
1410 52912 : get_velocity_z_squared(const int id, const int surf_offset)
1411 : {
1412 52912 : return std::pow(get_velocity_z(id, surf_offset), 2);
1413 : }
1414 :
1415 : void
1416 19640 : set_temperature(const int id, const dfloat value)
1417 : {
1418 19640 : nrs_t * nrs = (nrs_t *)nrsPtr();
1419 19640 : nrs->usrwrk[indices.temperature + id] = value;
1420 19640 : }
1421 :
1422 : void
1423 64817272 : set_flux(const int id, const dfloat value)
1424 : {
1425 64817272 : nrs_t * nrs = (nrs_t *)nrsPtr();
1426 64817272 : nrs->usrwrk[indices.flux + id] = value;
1427 64817272 : }
1428 :
1429 : void
1430 118273024 : set_heat_source(const int id, const dfloat value)
1431 : {
1432 118273024 : nrs_t * nrs = (nrs_t *)nrsPtr();
1433 118273024 : nrs->usrwrk[indices.heat_source + id] = value;
1434 118273024 : }
1435 :
1436 : void
1437 27648 : set_x_displacement(const int id, const dfloat value)
1438 : {
1439 27648 : mesh_t * mesh = entireMesh();
1440 27648 : mesh->x[id] = value;
1441 27648 : }
1442 :
1443 : void
1444 27648 : set_y_displacement(const int id, const dfloat value)
1445 : {
1446 27648 : mesh_t * mesh = entireMesh();
1447 27648 : mesh->y[id] = value;
1448 27648 : }
1449 :
1450 : void
1451 27648 : set_z_displacement(const int id, const dfloat value)
1452 : {
1453 27648 : mesh_t * mesh = entireMesh();
1454 27648 : mesh->z[id] = value;
1455 27648 : }
1456 :
1457 : void
1458 4032000 : set_mesh_velocity_x(const int id, const dfloat value)
1459 : {
1460 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1461 4032000 : nrs->usrwrk[indices.mesh_velocity_x + id] = value;
1462 4032000 : }
1463 :
1464 : void
1465 4032000 : set_mesh_velocity_y(const int id, const dfloat value)
1466 : {
1467 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1468 4032000 : nrs->usrwrk[indices.mesh_velocity_y + id] = value;
1469 4032000 : }
1470 :
1471 : void
1472 4032000 : set_mesh_velocity_z(const int id, const dfloat value)
1473 : {
1474 4032000 : nrs_t * nrs = (nrs_t *)nrsPtr();
1475 4032000 : nrs->usrwrk[indices.mesh_velocity_z + id] = value;
1476 4032000 : }
1477 :
1478 : void
1479 244 : checkFieldValidity(const field::NekWriteEnum & field)
1480 : {
1481 244 : switch (field)
1482 : {
1483 244 : case field::flux:
1484 244 : if (!hasTemperatureVariable())
1485 0 : mooseError("Cannot get NekRS heat flux "
1486 : "because your Nek case files do not have a temperature variable!");
1487 : break;
1488 0 : case field::heat_source:
1489 0 : if (!hasTemperatureVariable())
1490 0 : mooseError("Cannot get NekRS heat source "
1491 : "because your Nek case files do not have a temperature variable!");
1492 : break;
1493 : case field::x_displacement:
1494 : case field::y_displacement:
1495 : case field::z_displacement:
1496 : case field::mesh_velocity_x:
1497 : case field::mesh_velocity_y:
1498 : case field::mesh_velocity_z:
1499 : break;
1500 0 : default:
1501 0 : mooseError("Unhandled NekWriteEnum in checkFieldValidity!");
1502 : }
1503 244 : }
1504 :
1505 : void
1506 195308 : checkFieldValidity(const field::NekFieldEnum & field)
1507 : {
1508 : // by placing this check here, as opposed to inside the NekFieldInterface,
1509 : // we can also leverage this error checking for the 'outputs' of NekRSProblem,
1510 : // which does not inherit from NekFieldInterface but still accesses the solutionPointers.
1511 : // If this gets moved elsewhere, need to be sure to add dedicated testing for
1512 : // the 'outputs' on NekRSProblem.
1513 :
1514 : // TODO: would be nice for NekRSProblem to only access field information via the
1515 : // NekFieldInterface; refactor later
1516 :
1517 195308 : switch (field)
1518 : {
1519 87640 : case field::temperature:
1520 87640 : if (!hasTemperatureVariable())
1521 1 : mooseError("Cannot find 'temperature' "
1522 : "because your Nek case files do not have a temperature variable!");
1523 : break;
1524 85 : case field::scalar01:
1525 85 : if (!hasScalarVariable(1))
1526 1 : mooseError("Cannot find 'scalar01' "
1527 : "because your Nek case files do not have a scalar01 variable!");
1528 : break;
1529 201 : case field::scalar02:
1530 201 : if (!hasScalarVariable(2))
1531 1 : mooseError("Cannot find 'scalar02' "
1532 : "because your Nek case files do not have a scalar02 variable!");
1533 : break;
1534 57 : case field::scalar03:
1535 57 : if (!hasScalarVariable(3))
1536 1 : mooseError("Cannot find 'scalar03' "
1537 : "because your Nek case files do not have a scalar03 variable!");
1538 : break;
1539 454 : case field::usrwrk00:
1540 454 : if (n_usrwrk_slots < 1)
1541 2 : mooseError("Cannot find 'usrwrk00' because you have only allocated 'n_usrwrk_slots = " +
1542 1 : std::to_string(n_usrwrk_slots) + "'");
1543 : break;
1544 46 : case field::usrwrk01:
1545 46 : if (n_usrwrk_slots < 2)
1546 2 : mooseError("Cannot find 'usrwrk01' because you have only allocated 'n_usrwrk_slots = " +
1547 1 : std::to_string(n_usrwrk_slots) + "'");
1548 : break;
1549 46 : case field::usrwrk02:
1550 46 : if (n_usrwrk_slots < 3)
1551 2 : mooseError("Cannot find 'usrwrk02' because you have only allocated 'n_usrwrk_slots = " +
1552 1 : std::to_string(n_usrwrk_slots) + "'");
1553 : break;
1554 : }
1555 195301 : }
1556 :
1557 244 : double (*solutionPointer(const field::NekWriteEnum & field))(int, int)
1558 : {
1559 : double (*f)(int, int);
1560 :
1561 244 : checkFieldValidity(field);
1562 :
1563 244 : switch (field)
1564 : {
1565 244 : case field::flux:
1566 : f = &get_flux;
1567 : break;
1568 0 : default:
1569 0 : mooseError("Unhandled NekWriteEnum in solutionPointer!");
1570 : }
1571 :
1572 244 : return f;
1573 : }
1574 :
1575 191848 : double (*solutionPointer(const field::NekFieldEnum & field))(int, int)
1576 : {
1577 : // we include this here as well, in addition to within the NekFieldInterface, because
1578 : // the NekRSProblem accesses these methods without inheriting from NekFieldInterface
1579 191848 : checkFieldValidity(field);
1580 :
1581 : double (*f)(int, int);
1582 :
1583 191848 : switch (field)
1584 : {
1585 : case field::velocity_x:
1586 : f = &get_velocity_x;
1587 : break;
1588 17892 : case field::velocity_y:
1589 : f = &get_velocity_y;
1590 17892 : break;
1591 14984 : case field::velocity_z:
1592 : f = &get_velocity_z;
1593 14984 : break;
1594 196 : case field::velocity:
1595 : f = &get_velocity;
1596 196 : break;
1597 0 : case field::velocity_component:
1598 0 : mooseError("The 'velocity_component' field is not compatible with the solutionPointer "
1599 : "interface!");
1600 : break;
1601 48 : case field::velocity_x_squared:
1602 : f = &get_velocity_x_squared;
1603 48 : break;
1604 48 : case field::velocity_y_squared:
1605 : f = &get_velocity_y_squared;
1606 48 : break;
1607 52 : case field::velocity_z_squared:
1608 : f = &get_velocity_z_squared;
1609 52 : break;
1610 86376 : case field::temperature:
1611 : f = &get_temperature;
1612 86376 : break;
1613 31004 : case field::pressure:
1614 : f = &get_pressure;
1615 31004 : break;
1616 48 : case field::scalar01:
1617 : f = &get_scalar01;
1618 48 : break;
1619 176 : case field::scalar02:
1620 : f = &get_scalar02;
1621 176 : break;
1622 32 : case field::scalar03:
1623 : f = &get_scalar03;
1624 32 : break;
1625 9148 : case field::unity:
1626 : f = &get_unity;
1627 9148 : break;
1628 420 : case field::usrwrk00:
1629 : f = &get_usrwrk00;
1630 420 : break;
1631 16 : case field::usrwrk01:
1632 : f = &get_usrwrk01;
1633 16 : break;
1634 16 : case field::usrwrk02:
1635 : f = &get_usrwrk02;
1636 16 : break;
1637 0 : default:
1638 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1639 : }
1640 :
1641 191848 : return f;
1642 : }
1643 :
1644 1862524 : void (*solutionWritePointer(const field::NekWriteEnum & field))(int, dfloat)
1645 : {
1646 : void (*f)(int, dfloat);
1647 :
1648 1862524 : switch (field)
1649 : {
1650 : case field::flux:
1651 : f = &set_flux;
1652 : break;
1653 396748 : case field::heat_source:
1654 : f = &set_heat_source;
1655 396748 : break;
1656 1024 : case field::x_displacement:
1657 : f = &set_x_displacement;
1658 1024 : break;
1659 1024 : case field::y_displacement:
1660 : f = &set_y_displacement;
1661 1024 : break;
1662 1024 : case field::z_displacement:
1663 : f = &set_z_displacement;
1664 1024 : break;
1665 179200 : case field::mesh_velocity_x:
1666 : f = &set_mesh_velocity_x;
1667 179200 : break;
1668 179200 : case field::mesh_velocity_y:
1669 : f = &set_mesh_velocity_y;
1670 179200 : break;
1671 179200 : case field::mesh_velocity_z:
1672 : f = &set_mesh_velocity_z;
1673 179200 : break;
1674 0 : default:
1675 0 : throw std::runtime_error("Unhandled NekWriteEnum!");
1676 : }
1677 :
1678 1862524 : return f;
1679 : }
1680 :
1681 2972 : void (*solutionWritePointer(const field::NekFieldEnum & field))(int, dfloat)
1682 : {
1683 : void (*f)(int, dfloat);
1684 :
1685 2972 : switch (field)
1686 : {
1687 2972 : case field::temperature:
1688 : f = &set_temperature;
1689 : break;
1690 0 : default:
1691 0 : throw std::runtime_error("Unhandled NekFieldEnum in solutionWritePointer! Other write fields "
1692 0 : "have not been added to this interface yet.");
1693 : }
1694 :
1695 2972 : return f;
1696 : }
1697 :
1698 : void
1699 129 : initializeDimensionalScales(const double U,
1700 : const double T,
1701 : const double dT,
1702 : const double L,
1703 : const double rho,
1704 : const double Cp,
1705 : const double s01,
1706 : const double ds01,
1707 : const double s02,
1708 : const double ds02,
1709 : const double s03,
1710 : const double ds03)
1711 : {
1712 129 : scales.U_ref = U;
1713 129 : scales.T_ref = T;
1714 129 : scales.dT_ref = dT;
1715 129 : scales.L_ref = L;
1716 129 : scales.A_ref = L * L;
1717 129 : scales.V_ref = L * L * L;
1718 129 : scales.rho_ref = rho;
1719 129 : scales.Cp_ref = Cp;
1720 129 : scales.t_ref = L / U;
1721 129 : scales.P_ref = rho * U * U;
1722 :
1723 129 : scales.s01_ref = s01;
1724 129 : scales.ds01_ref = ds01;
1725 129 : scales.s02_ref = s02;
1726 129 : scales.ds02_ref = ds02;
1727 129 : scales.s03_ref = s03;
1728 129 : scales.ds03_ref = ds03;
1729 :
1730 129 : scales.flux_ref = rho * U * Cp * dT;
1731 129 : scales.source_ref = scales.flux_ref / L;
1732 129 : }
1733 :
1734 : double
1735 457 : referenceLength()
1736 : {
1737 457 : return scales.L_ref;
1738 : }
1739 :
1740 : double
1741 273522 : referenceTime()
1742 : {
1743 273522 : return scales.t_ref;
1744 : }
1745 :
1746 : double
1747 263 : referenceArea()
1748 : {
1749 263 : return scales.A_ref;
1750 : }
1751 :
1752 : double
1753 1132 : referenceVolume()
1754 : {
1755 1132 : return scales.V_ref;
1756 : }
1757 :
1758 : Real
1759 66784620 : nondimensionalAdditive(const field::NekFieldEnum & field)
1760 : {
1761 66784620 : switch (field)
1762 : {
1763 24766184 : case field::temperature:
1764 24766184 : return scales.T_ref;
1765 8660 : case field::scalar01:
1766 8660 : return scales.s01_ref;
1767 40904 : case field::scalar02:
1768 40904 : return scales.s02_ref;
1769 5576 : case field::scalar03:
1770 5576 : return scales.s03_ref;
1771 : default:
1772 : return 0;
1773 : }
1774 : }
1775 :
1776 : Real
1777 19382 : nondimensionalAdditive(const field::NekWriteEnum & field)
1778 : {
1779 19382 : switch (field)
1780 : {
1781 19382 : case field::flux:
1782 : case field::heat_source:
1783 : case field::x_displacement:
1784 : case field::y_displacement:
1785 : case field::z_displacement:
1786 : case field::mesh_velocity_x:
1787 : case field::mesh_velocity_y:
1788 : case field::mesh_velocity_z:
1789 19382 : return 0.0;
1790 0 : default:
1791 0 : mooseError("Unhandled NekWriteEnum in nondimensionalAdditive!");
1792 : }
1793 : }
1794 :
1795 : Real
1796 48825 : nondimensionalDivisor(const field::NekWriteEnum & field)
1797 : {
1798 48825 : switch (field)
1799 : {
1800 44402 : case field::flux:
1801 44402 : return scales.flux_ref;
1802 3519 : case field::heat_source:
1803 3519 : return scales.source_ref;
1804 904 : case field::x_displacement:
1805 : case field::y_displacement:
1806 : case field::z_displacement:
1807 904 : return scales.L_ref;
1808 0 : case field::mesh_velocity_x:
1809 : case field::mesh_velocity_y:
1810 : case field::mesh_velocity_z:
1811 0 : return scales.U_ref;
1812 0 : default:
1813 0 : mooseError("Unhandled NekWriteEnum in nondimensionalDivisor!");
1814 : }
1815 : }
1816 :
1817 : Real
1818 67322346 : nondimensionalDivisor(const field::NekFieldEnum & field)
1819 : {
1820 67322346 : switch (field)
1821 : {
1822 32127380 : case field::velocity_x:
1823 : case field::velocity_y:
1824 : case field::velocity_z:
1825 : case field::velocity:
1826 : case field::velocity_component:
1827 32127380 : return scales.U_ref;
1828 5336 : case field::velocity_x_squared:
1829 : case field::velocity_y_squared:
1830 : case field::velocity_z_squared:
1831 5336 : return scales.U_ref * scales.U_ref;
1832 24766184 : case field::temperature:
1833 24766184 : return scales.dT_ref;
1834 10342955 : case field::pressure:
1835 10342955 : return scales.P_ref;
1836 8660 : case field::scalar01:
1837 8660 : return scales.ds01_ref;
1838 40904 : case field::scalar02:
1839 40904 : return scales.ds02_ref;
1840 5576 : case field::scalar03:
1841 5576 : return scales.ds03_ref;
1842 : case field::unity:
1843 : // no dimensionalization needed
1844 : return 1.0;
1845 445 : case field::usrwrk00:
1846 445 : return scratchUnits(0);
1847 29 : case field::usrwrk01:
1848 29 : return scratchUnits(1);
1849 29 : case field::usrwrk02:
1850 29 : return scratchUnits(2);
1851 0 : default:
1852 0 : throw std::runtime_error("Unhandled 'NekFieldEnum'!");
1853 : }
1854 : }
1855 :
1856 : Real
1857 503 : scratchUnits(const int slot)
1858 : {
1859 503 : if (indices.flux != -1 && slot == indices.flux / nekrs::fieldOffset())
1860 424 : return scales.flux_ref;
1861 79 : else if (indices.heat_source != -1 && slot == indices.heat_source / nekrs::fieldOffset())
1862 4 : return scales.source_ref;
1863 75 : else if (is_nondimensional)
1864 : {
1865 : // TODO: we are lazy and did not include all the usrwrk indices
1866 66 : mooseDoOnce(mooseWarning(
1867 : "The units of 'usrwrk0" + std::to_string(slot) +
1868 : "' are unknown, so we cannot dimensionalize any objects using 'field = usrwrk0" +
1869 : std::to_string(slot) +
1870 : "'. The output for this quantity will be given in non-dimensional form.\n\nYou will need "
1871 : "to manipulate the data manually from Cardinal if you need to dimensionalize it."));
1872 : }
1873 :
1874 : return 1.0;
1875 : }
1876 :
1877 : void
1878 825 : nondimensional(const bool n)
1879 : {
1880 825 : is_nondimensional = n;
1881 825 : }
1882 :
1883 : template <>
1884 : MPI_Datatype
1885 183358 : resolveType<double>()
1886 : {
1887 183358 : return MPI_DOUBLE;
1888 : }
1889 :
1890 : template <>
1891 : MPI_Datatype
1892 7240 : resolveType<int>()
1893 : {
1894 7240 : return MPI_INT;
1895 : }
1896 :
1897 : } // end namespace nekrs
1898 :
1899 : #endif
|
