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 : #pragma once
20 :
21 : #include "CardinalEnums.h"
22 : #include "MooseTypes.h"
23 : #include "NekBoundaryCoupling.h"
24 : #include "NekVolumeCoupling.h"
25 : #include "Function.h"
26 :
27 : #include "inipp.hpp"
28 : #include "nekrs.hpp"
29 : #include "nrs.hpp"
30 : #include "udf.hpp"
31 : #include "inipp.hpp"
32 : #include "mesh.h"
33 :
34 : #include "libmesh/point.h"
35 :
36 : #include <string>
37 : #include <vector>
38 :
39 : /**
40 : * \brief Cardinal-specific nekRS API
41 : *
42 : * nekRS ships with a rudimentary API in their nekrs namespace, but we need additional
43 : * functionality from within Cardinal. Many of these functions are quite basic and could
44 : * eventually be ported back into nekRS itself.
45 : */
46 : namespace nekrs
47 : {
48 :
49 : static int build_only;
50 :
51 : /**
52 : * Number of passive scalars
53 : * @return number of passive scalars
54 : */
55 : int Nscalar();
56 :
57 : /// Allocate memory for the host mesh parameters
58 : void initializeHostMeshParameters();
59 :
60 : /// Update the mesh parameters on host
61 : void updateHostMeshParameters();
62 :
63 : dfloat * host_x();
64 : dfloat * host_y();
65 : dfloat * host_z();
66 :
67 : nrs_t * nrsPtr();
68 :
69 : dfloat * getSgeo();
70 : dfloat * getVgeo();
71 :
72 : /**
73 : * Check that the field specified can be accessed, e.g., if a user is requesting
74 : * to access temperature, the problem must have a temperature variable
75 : * @param[in] field field to check
76 : */
77 : void checkFieldValidity(const field::NekFieldEnum & field);
78 : void checkFieldValidity(const field::NekWriteEnum & field);
79 :
80 : void copyDeviceToHost();
81 :
82 : /**
83 : * Compute the minimum distance to a wall as a mesh field
84 : * @param[in] boundary_id boundary(s) on which to compute the distance to the wall
85 : */
86 : void computeWallDistance(const std::vector<int> & boundary_id);
87 :
88 : /**
89 : * Compute y+ on the NekRS mesh
90 : * @param[in] boundary_id boundary(s) on which to compute y+
91 : * @return max, min, average y+
92 : */
93 : std::vector<dfloat> yPlus(const std::vector<int> & boundary_id);
94 :
95 : /**
96 : * Compute the three components of the viscous drag along a given boundary
97 : * @param[in] boundary boundary IDs for drag computation
98 : * @return components of drag (force fluid exerts on walls)
99 : */
100 : std::vector<dfloat> viscousDrag(const std::vector<int> & boundary);
101 :
102 : /**
103 : * Set the absolute tolerance for checking energy conservation in data transfers to Nek
104 : * @param[in] tol tolerance
105 : */
106 : void setAbsoluteTol(double tol);
107 :
108 : /**
109 : * Return the reference units for a usrwrk slot
110 : * @param[in] slot usrwrk slot
111 : * @return value by which to multiply the usrwrk slot to go from non-dimensional form into
112 : * dimensional form
113 : */
114 : Real scratchUnits(const int slot);
115 :
116 : /**
117 : * Inform backend if dimensionalization should be performed
118 : * @param[in] n if dimensionalize should be performed
119 : */
120 : void nondimensional(const bool n);
121 :
122 : /**
123 : * Set the relative tolerance for checking energy conservation in data transfers to Nek
124 : * @param[in] tol tolerance
125 : */
126 : void setRelativeTol(double tol);
127 :
128 : /**
129 : * Nek's runtime statistics are formed by collecting a timer of both the initialization
130 : * and accumulated run time. We unfortunately have to split this across multiple classes,
131 : * so if we want correct times we need to have NekInitAction save the value of the time
132 : * spent on initialization.
133 : * @param[in] time time spent on initialization
134 : */
135 : void setNekSetupTime(const double & time);
136 :
137 : /**
138 : * Get time spent on initialization
139 : * @return time spent on initialization
140 : */
141 : double getNekSetupTime();
142 :
143 : /**
144 : * Set the start time used by NekRS
145 : * @param[in] start start time
146 : */
147 : void setStartTime(const double & start);
148 :
149 : /**
150 : * Whether NekRS itself has been initialized yet
151 : * @return whether NekRS is initialized
152 : */
153 : bool isInitialized();
154 :
155 : /**
156 : * Write a field file containing a specific slot of the nrs->usrwrk scratch space;
157 : * this will write the field to the 'temperature' slot in a field file.
158 : * @param[in] slot index in the nrs->usrwrk array to write
159 : * @param[in] prefix prefix for file name
160 : * @param[in] time simulation time to write file for
161 : * @param[in] step time step index
162 : * @param[in] write_coords whether to write the mesh coordinates
163 : */
164 : void write_usrwrk_field_file(const int & size,
165 : const int & index,
166 : const int & slot,
167 : const std::string & prefix,
168 : const dfloat & time,
169 : const int & step,
170 : const bool & write_coords);
171 :
172 : /**
173 : * Write a field file containing pressure, velocity, and scalars with given prefix
174 : * @param[in] prefix three-character prefix
175 : * @param[in] time time
176 : * @param[in] step time step index
177 : */
178 : void write_field_file(const std::string & prefix, const dfloat time, const int & step);
179 :
180 : /**
181 : * Indicate whether NekRS was run in build-only mode (this doesn't actually
182 : * cause NekRS to run in build-only mode, but only provides an interface to
183 : * this information elsewhere).
184 : * @param[in] buildOnly whether NekRS is to be run in build-only mode
185 : */
186 : void buildOnly(int buildOnly);
187 :
188 : /**
189 : * Whether NekRS was run in JIT build-only mode
190 : * @return whether NekRS was run in build-only mode
191 : */
192 : int buildOnly();
193 :
194 : /**
195 : * Interpolate a volume between NekRS's GLL points and a given-order receiving/sending mesh
196 : */
197 : void interpolateVolumeHex3D(const double * I, double * x, int N, double * Ix, int M);
198 :
199 : /**
200 : * Whether nekRS's input file has CHT
201 : * @return whether nekRS input files model CHT
202 : */
203 : bool hasCHT();
204 :
205 : /**
206 : * Whether nekRS's input file indicates a moving mesh
207 : * @return whether nekRS's input file indicates a moving mesh
208 : */
209 : bool hasMovingMesh();
210 :
211 : /**
212 : * Whether nekRS's input file indicates a variable time stepping scheme
213 : * @return whether nekRS's input file indicates a variable time stepping
214 : */
215 : bool hasVariableDt();
216 :
217 : /**
218 : * Whether nekRS's input file has the blending mesh solver
219 : * @return whether nekRS's input file has a non-user [MESH] solver
220 : */
221 : bool hasBlendingSolver();
222 :
223 : /**
224 : * Whether nekRS's input file has the user mesh solver
225 : * @return whether nekRS's input file has [MESH] solver = user
226 : */
227 : bool hasUserMeshSolver();
228 :
229 : /**
230 : * Whether nekRS's input file intends to terminate the simulation based on a wall time
231 : * @return whether a wall time is used in nekRS to end the simulation
232 : */
233 : bool endControlElapsedTime();
234 :
235 : /**
236 : * Whether nekRS's input file intends to terminate the simulation based on an end time
237 : * @return whether an end time is used in nekRS to end the simulation
238 : */
239 : bool endControlTime();
240 :
241 : /**
242 : * Whether nekRS's input file intends to terminate the simulation based on a number of steps
243 : * @return whether a time step interval is used in nekRS to end the simulation
244 : */
245 : bool endControlNumSteps();
246 :
247 : /**
248 : * Offset increment for indexing into multi-volume arrays for the scalar fields.
249 : * This assumes that all scalars are the same length as the temperature scalar.
250 : * TODO: evaluate whether this works if nekRS uses CHT
251 : * @return scalar field offset
252 : */
253 : int scalarFieldOffset();
254 :
255 : /**
256 : * Offset increment for indexing into the velocity array
257 : * @return velocity field offset
258 : */
259 : int velocityFieldOffset();
260 :
261 : /**
262 : * Offset increment to use for generic slice indexing
263 : * @return field offset
264 : */
265 : int fieldOffset();
266 :
267 : /**
268 : * Get the "entire" NekRS mesh. For cases with a temperature scalar, this returns
269 : * nrs->meshT, which will cover both the fluid and solid regions if CHT is present.
270 : * For flow-only cases, this will return the flow mesh.
271 : * @return entire NekRS mesh
272 : */
273 : mesh_t * entireMesh();
274 :
275 : /**
276 : * Get the mesh for the flow solve
277 : * @return flow mesh
278 : */
279 : mesh_t * flowMesh();
280 :
281 : /**
282 : * Get the mesh for the temperature scalar
283 : * @return temperature mesh
284 : */
285 : mesh_t * temperatureMesh();
286 :
287 : /**
288 : * Get the mesh to act on
289 : * @param[in] pp_mesh which NekRS mesh to operate on
290 : * @return mesh to act on
291 : */
292 : mesh_t * getMesh(const nek_mesh::NekMeshEnum pp_mesh);
293 :
294 : /**
295 : * Get the process rank
296 : * @return process rank
297 : */
298 : int commRank();
299 :
300 : /**
301 : * Get the communicator size
302 : * @return communicator size
303 : */
304 : int commSize();
305 :
306 : /**
307 : * Whether nekRS's input file indicates that the problem has a temperature variable
308 : * @return whether the nekRS problem includes a temperature variable
309 : */
310 : bool hasTemperatureVariable();
311 :
312 : /**
313 : * Whether nekRS actually solves for temperature (as opposed to setting its solver to 'none')
314 : * @return whether nekRS will solve for temperature
315 : */
316 : bool hasTemperatureSolve();
317 :
318 : /**
319 : * Whether nekRS's input file indicates that the problem has a scalar0(scalarId) variable
320 : * @param[in] scalarId scalar number, i.e. for scalar03 scalarId=3
321 : * @return whether the nekRS problem includes the scalar0(scalarId) variable
322 : */
323 : bool hasScalarVariable(int scalarId);
324 :
325 : /**
326 : * Whether nekRS contains an OCCA kernel to apply a source to the passive scalar equations
327 : * @return whether nekRS has an OCCA kernel for apply a passive scalar source
328 : */
329 : bool hasHeatSourceKernel();
330 :
331 : /**
332 : * Whether the scratch space has already been allocated by the user
333 : * @return whether scratch space is already allocated
334 : */
335 : bool scratchAvailable();
336 :
337 : /**
338 : * Initialize scratch space for data to get sent into NekRS
339 : * @param[in] n_slots number of slots (for volume arrays) to allocate
340 : */
341 : void initializeScratch(const unsigned int & n_slots);
342 :
343 : /// Free the scratch space
344 : void freeScratch();
345 :
346 : /**
347 : * Get the viscosity used in the definition of the Reynolds number; note that
348 : * for dimensional cases, this is only guaranteed to be correct if the viscosity is constant.
349 : * @return constant dynamic viscosity
350 : */
351 : double viscosity();
352 :
353 : /**
354 : * Get the Prandtl number; note that for dimensional cases, this is only guaranteed
355 : * to be correct if the density, viscosity, heat capacity, and conductivity are constant.
356 : * @return constant Prandtl number
357 : */
358 : double Pr();
359 :
360 : /// Copy the deformation from host to device
361 : void copyDeformationToDevice();
362 :
363 : template <typename T>
364 : void allgatherv(const std::vector<int> & base_counts,
365 : const T * input,
366 : T * output,
367 : const int multiplier = 1);
368 :
369 : /**
370 : * Determine the receiving counts and displacements for all gather routines
371 : * @param[in] base_counts unit-wise receiving counts for each process
372 : * @param[out] counts receiving counts from each process
373 : * @param[out] displacement displacement for each process's counts
374 : * @param[in] multiplier optional multiplier on the face-based data
375 : */
376 : void displacementAndCounts(const std::vector<int> & base_counts,
377 : int * counts,
378 : int * displacement,
379 : const int multiplier);
380 :
381 : /**
382 : * Form the 2-D interpolation matrix from a starting GLL quadrature rule to an ending
383 : * GLL quadrature rule.
384 : * @param[out] I interpolation matrix
385 : * @param[in] starting_points number of points in the source quadrature rule
386 : * @param[in] ending_points number of points in the end quadrature rule
387 : */
388 : void interpolationMatrix(double * I, int starting_points, int ending_points);
389 :
390 : /**
391 : * Interpolate face data onto a new set of points
392 : * @param[in] scratch available scratch space for the calculation
393 : * @param[in] I interpolation matrix
394 : * @param[in] x face data to be interpolated
395 : * @param[in] N number of points in 1-D to be interpolated
396 : * @param[out] Ix interpolated data
397 : * @param[in] M resulting number of interpolated points in 1-D
398 : */
399 : void interpolateSurfaceFaceHex3D(
400 : double * scratch, const double * I, double * x, int N, double * Ix, int M);
401 :
402 : /**
403 : * Compute the face centroid given a local element ID and face ID (NOTE: returns in dimensional
404 : * form)
405 : * @param[in] local_elem_id local element ID on this rank
406 : * @param[in] local_face_id local face ID on the element
407 : * @return centroid
408 : */
409 : Point centroidFace(int local_elem_id, int local_face_id);
410 :
411 : /**
412 : * Compute the centroid given a local element ID (NOTE: returns in dimensional form)
413 : * @param[in] local_elem_id local element ID on this rank
414 : * @return centroid
415 : */
416 : Point centroid(int local_elem_id);
417 :
418 : /**
419 : * Get the coordinate given a local element ID and local node ID (NOTE: returns in dimensional form)
420 : * @param[in] local_elem_id local element ID on this rank
421 : * @param[in] local_node_id local node ID on this element
422 : * @return point
423 : */
424 : Point gllPoint(int local_elem_id, int local_node_id);
425 :
426 : /**
427 : * Get the coordinate given a local element ID, a local face ID, and local node ID (NOTE: returns in
428 : * dimensional form)
429 : * @param[in] local_elem_id local element ID on this rank
430 : * @param[in] local_face_id local face ID on this element
431 : * @param[in] local_node_id local node ID on this element
432 : * @return point
433 : */
434 : Point gllPointFace(int local_elem_id, int local_face_id, int local_node_id);
435 :
436 : /**
437 : * Integrate the scratch space over boundaries
438 : * @param[in] slot slot in scratch space
439 : * @param[in] boundary boundaries over which to integrate the scratch space
440 : * @param[in] pp_mesh portion of NekRS mesh to integrate over
441 : * @return boundary integrated scratch space, with one value per sideset
442 : */
443 : std::vector<double> usrwrkSideIntegral(const unsigned int & slot,
444 : const std::vector<int> & boundary,
445 : const nek_mesh::NekMeshEnum pp_mesh);
446 :
447 : /**
448 : * Volume integrate the scratch space
449 : * @param[in] slot slot in scratch space to i ntegrat
450 : * @param[in] pp_mesh NekRS mesh to integrate over
451 : * @return volume integrated scratch space
452 : */
453 : double usrwrkVolumeIntegral(const unsigned int & slot, const nek_mesh::NekMeshEnum pp_mesh);
454 :
455 : /**
456 : * Scale a slot in the usrwrk by a fixed value (multiplication)
457 : * @param[in] slot slot in usrwrk to modify
458 : * @param[in] value value to multiply on scratch slot
459 : */
460 : void scaleUsrwrk(const unsigned int & slot, const dfloat & value);
461 :
462 : /**
463 : * Compute the area of a set of boundary IDs
464 : * @param[in] boundary_id nekRS boundary IDs for which to perform the integral
465 : * @param[in] pp_mesh which NekRS mesh to operate on
466 : * @return area integral
467 : */
468 : double area(const std::vector<int> & boundary_id, const nek_mesh::NekMeshEnum pp_mesh);
469 :
470 : /**
471 : * Compute the area integral of a given integrand over a set of boundary IDs
472 : * @param[in] boundary_id nekRS boundary IDs for which to perform the integral
473 : * @param[in] integrand field to integrate
474 : * @param[in] pp_mesh which NekRS mesh to operate on
475 : * @return area integral of a field
476 : */
477 : double sideIntegral(const std::vector<int> & boundary_id, const field::NekFieldEnum & integrand,
478 : const nek_mesh::NekMeshEnum pp_mesh);
479 :
480 : /**
481 : * Compute the volume over the entire scalar mesh
482 : * @param[in] pp_mesh which NekRS mesh to operate on
483 : * @return volume integral
484 : */
485 : double volume(const nek_mesh::NekMeshEnum pp_mesh);
486 :
487 : /**
488 : * Dimensionalize a volume
489 : * @param[in] integral integral to dimensionalize
490 : */
491 : void dimensionalizeVolume(double & integral);
492 :
493 : /**
494 : * Dimensionalize an area
495 : * @param[in] integral integral to dimensionalize
496 : */
497 : void dimensionalizeArea(double & integral);
498 :
499 : /**
500 : * Dimensionalize a given integral of f over volume, i.e. fdV
501 : * @param[in] integrand field to dimensionalize
502 : * @param[in] volume volume of the domain (only used for dimensionalizing temperature)
503 : * @param[in] integral integral to dimensionalize
504 : */
505 : void dimensionalizeVolumeIntegral(const field::NekFieldEnum & integrand,
506 : const Real & volume,
507 : double & integral);
508 :
509 : /**
510 : * Dimensionalize a given integral of f over a side, i.e. fdS
511 : * @param[in] integrand field to dimensionalize
512 : * @param[in] area area of the boundary
513 : * @param[in] integral integral to dimensionalize
514 : */
515 : void dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
516 : const Real & area,
517 : double & integral);
518 :
519 : /**
520 : * Dimensionalize a given integral of f over a side, i.e. fdS
521 : * @param[in] integrand field to dimensionalize
522 : * @param[in] boundary_id boundary IDs for the integral
523 : * @param[in] integral integral to dimensionalize
524 : * @param[in] pp_mesh which NekRS mesh to operate on
525 : */
526 : void dimensionalizeSideIntegral(const field::NekFieldEnum & integrand,
527 : const std::vector<int> & boundary_id,
528 : double & integral,
529 : const nek_mesh::NekMeshEnum pp_mesh);
530 :
531 : /**
532 : * Compute the volume integral of a given integrand over the entire scalar mesh
533 : * @param[in] integrand field to integrate
534 : * @param[in] volume volume of the domain (only used for dimensionalizing temperature)
535 : * @param[in] pp_mesh which NekRS mesh to operate on
536 : * @return volume integral of a field
537 : */
538 : double volumeIntegral(const field::NekFieldEnum & integrand,
539 : const double & volume,
540 : const nek_mesh::NekMeshEnum pp_mesh);
541 :
542 : /**
543 : * Transform the point and time passed into a function into dimensional form, because
544 : * the functions on the MOOSE side are defined in dimensional form.
545 : * @param[in] f function to query
546 : * @param[in] time time
547 : * @param[in] id element ID
548 : */
549 : double evaluateFunctionOnMesh(const Function * f, const Real time, const int id);
550 :
551 : /**
552 : * Compute the L^N norm of a given integrand over the mesh
553 : * @param[in] integrand field to integrate
554 : * @param[in] pp_mesh which NekRS mesh to operate on
555 : * @param[in] function MOOSE function to use to shift the field
556 : * @param[in] time time to evaluate function at
557 : * @param[in] N order of the norm
558 : * @return integrated L^N norm of the NekRS field, relative to a function
559 : */
560 : double volumeNorm(const field::NekFieldEnum & integrand,
561 : const nek_mesh::NekMeshEnum pp_mesh,
562 : const Function * function,
563 : const Real & time,
564 : const Real & N);
565 :
566 : /**
567 : * Compute the mass flowrate over a set of boundary IDs
568 : * @param[in] boundary_id nekRS boundary IDs for which to compute the mass flowrate
569 : * @param[in] pp_mesh which NekRS mesh to operate on
570 : * @return mass flowrate
571 : */
572 : double massFlowrate(const std::vector<int> & boundary_id,
573 : const nek_mesh::NekMeshEnum pp_mesh);
574 :
575 : /**
576 : * Compute the mass flux weighted integral of a given integrand over a set of boundary IDs
577 : * @param[in] boundary_id nekRS boundary IDs for which to perform the integral
578 : * @param[in] integrand field to integrate and weight by mass flux
579 : * @param[in] pp_mesh which NekRS mesh to operate on
580 : * @return mass flux weighted area average of a field
581 : */
582 : double sideMassFluxWeightedIntegral(const std::vector<int> & boundary_id,
583 : const field::NekFieldEnum & integrand,
584 : const nek_mesh::NekMeshEnum pp_mesh);
585 :
586 : /**
587 : * Compute the integral of pressure on a surface, multiplied by the unit normal
588 : * of the surface with a specified direction vector. This represents the force
589 : * that the fluid exerts ON the boundary.
590 : * @param[in] boundary_id NekRS boundary IDs for which to perform the integral
591 : * @param[in] direction unit vector to dot with the boundary surface normal
592 : * @param[in] pp_mesh which NekRS mesh to operate on
593 : * @return pressure surface force, along a particular direction
594 : */
595 : double pressureSurfaceForce(const std::vector<int> & boundary_id, const Point & direction, const nek_mesh::NekMeshEnum pp_mesh);
596 :
597 : /**
598 : * Compute the heat flux over a set of boundary IDs
599 : * @param[in] boundary_id nekRS boundary IDs for which to perform the integral
600 : * @param[in] pp_mesh which NekRS mesh to operate on
601 : * @return heat flux area integral
602 : */
603 : double heatFluxIntegral(const std::vector<int> & boundary_id,
604 : const nek_mesh::NekMeshEnum pp_mesh);
605 :
606 : /**
607 : * Limit the temperature in nekRS to within the range of [min_T, max_T]
608 : * @param[in] min_T minimum temperature allowable in nekRS
609 : * @param[in] max_T maximum temperature allowable in nekRS
610 : */
611 : void limitTemperature(const double * min_T, const double * max_T);
612 :
613 : /**
614 : * Compute the gradient of a volume field
615 : * @param[in] offset in the gradient field for each component (grad_x, grad_y, or grad_z)
616 : * @param[in] e element ID to compute gradient
617 : * @param[in] f field to compute the gradient of
618 : * @param[in] pp_mesh which NekRS mesh to operate on
619 : * @param[out] grad_f gradient of field
620 : */
621 : void gradient(const int offset,
622 : const int e,
623 : const double * f,
624 : double * grad_f,
625 : const nek_mesh::NekMeshEnum pp_mesh);
626 :
627 : /**
628 : * Find the extreme value of a given field over the entire nekRS domain
629 : * @param[in] field field to find the minimum value of
630 : * @param[in] pp_mesh which NekRS mesh to operate on
631 : * @param[in] max whether to take the maximum (or if false, the minimum)
632 : * @return max or min value of field in volume
633 : */
634 : double volumeExtremeValue(const field::NekFieldEnum & field,
635 : const nek_mesh::NekMeshEnum pp_mesh,
636 : const bool max);
637 :
638 : /**
639 : * Find the extreme of a given field over a set of boundary IDs
640 : * @param[in] boundary_id nekRS boundary IDs for which to find the extreme value
641 : * @param[in] field field to find the maximum value of
642 : * @param[in] pp_mesh which NekRS mesh to operate on
643 : * @param[in] max whether to take the maximum (or if false, the minimum)
644 : * @return max or min value of field on boundary
645 : */
646 : double sideExtremeValue(const std::vector<int> & boundary_id, const field::NekFieldEnum & field,
647 : const nek_mesh::NekMeshEnum pp_mesh, const bool max);
648 :
649 : /**
650 : * Number of faces per element; because NekRS only supports HEX20, this should be 6
651 : * @return number of faces per mesh element
652 : */
653 : int Nfaces();
654 :
655 : /**
656 : * Whether the specific boundary is a flux boundary
657 : * @param[in] boundary boundary ID
658 : * @return whether boundary is a flux boundary
659 : */
660 : bool isHeatFluxBoundary(const int boundary);
661 :
662 : /**
663 : * Whether the specific boundary is a moving mesh boundary
664 : * @param[in] boundary boundary ID
665 : * @return whether boundary is a moving mesh boundary
666 : */
667 : bool isMovingMeshBoundary(const int boundary);
668 :
669 : /**
670 : * Whether the specific boundary is a specified temperature boundary
671 : * @param[in] boundary boundary ID
672 : * @return whether boundary is a temperature boundary
673 : */
674 : bool isTemperatureBoundary(const int boundary);
675 :
676 : /**
677 : * String name indicating the temperature boundary condition type on a given boundary
678 : * @param[in] boundary boundary ID
679 : * @return string name of boundary condition type
680 : */
681 : const std::string temperatureBoundaryType(const int boundary);
682 :
683 : /**
684 : * Polynomial order used in nekRS solution
685 : * @return polynomial order
686 : */
687 : int polynomialOrder();
688 :
689 : /**
690 : * Total number of volume elements in nekRS mesh summed over all processes
691 : * @return number of volume elements
692 : */
693 : int Nelements();
694 :
695 : /**
696 : * Mesh dimension
697 : * @return mesh dimension
698 : */
699 : int dim();
700 :
701 : /**
702 : * \brief Number of vertices required to define an element face
703 : * Vertices refer to the points required to place the "corners" of an element face,
704 : * and _not_ the quadrature points. For instance, for hexahedral elements, the number of vertices
705 : * per face is 4 regardless of the polynomial order.
706 : * @return Number of vertices per element face
707 : */
708 : int NfaceVertices();
709 :
710 : /**
711 : * Total number of element faces on a boundary of the nekRS mesh summed over all processes
712 : * @return number of boundary element faces
713 : */
714 : int NboundaryFaces();
715 :
716 : /**
717 : * Number of boundary IDs in the nekRS mesh
718 : * @return number of boundary IDs
719 : */
720 : int NboundaryID();
721 :
722 : /**
723 : * Whether the provided boundary IDs are all valid in the nekRS mesh
724 : * @param[in] boundary_id vector of boundary IDs to check
725 : * @param[out] first_invalid_id first invalid ID encountered for printing an error on the MOOSE side
726 : * @param[out] n_boundaries maximum valid boundary ID for printing an error on the MOOSE side
727 : * @return whether all boundaries are valid
728 : */
729 : bool
730 : validBoundaryIDs(const std::vector<int> & boundary_id, int & first_invalid_id, int & n_boundaries);
731 :
732 : /**
733 : * Store the rank-local element, element-local face, and rank ownership for boundary coupling
734 : * @param[in] boundary_id boundaries through which nekRS will be coupled
735 : * @param[out] N total number of surface elements
736 : */
737 : void storeBoundaryCoupling(const std::vector<int> & boundary_id, int & N);
738 :
739 : /**
740 : * Characteristic scales assumed in nekRS if using a non-dimensional solution; initial values
741 : * are applied, which will be overridden by the DimensionalizeAction in Cardinal.
742 : */
743 : struct characteristicScales
744 : {
745 : double U_ref = 1;
746 : double T_ref = 0;
747 : double dT_ref = 1;
748 : double P_ref = 1;
749 : double L_ref = 1;
750 : double A_ref = 1;
751 : double V_ref = 1;
752 : double rho_ref = 1;
753 : double Cp_ref = 1;
754 : double flux_ref = 1;
755 : double source_ref = 1;
756 : double t_ref = 1;
757 : double s01_ref = 0;
758 : double ds01_ref = 1;
759 : double s02_ref = 0;
760 : double ds02_ref = 1;
761 : double s03_ref = 0;
762 : double ds03_ref = 1;
763 : };
764 :
765 : /**
766 : * Get pointer to various solution functions (for reading only) based on enumeration
767 : * @param[in] field field to return a pointer to
768 : * @return function pointer to method that returns said field as a function of GLL index
769 : */
770 : double (*solutionPointer(const field::NekFieldEnum & field))(int, int);
771 : double (*solutionPointer(const field::NekWriteEnum & field))(int, int);
772 :
773 : /**
774 : * Determine the index in the scalar array for the non-temperature scalars
775 : * @param[in] id id of the scalar; 1 corresponds to first non-temperature scalar, etc.
776 : * @return slot in scalar array that holds this variable
777 : */
778 : int scalarSlot(const int id);
779 :
780 : /**
781 : * Get the scalar01 solution at given GLL index
782 : * @param[in] id GLL index
783 : * @return scalar01 value at index
784 : */
785 : double get_scalar01(const int id, const int surf_offset);
786 :
787 : /**
788 : * Get the scalar02 solution at given GLL index
789 : * @param[in] id GLL index
790 : * @return scalar02 value at index
791 : */
792 : double get_scalar02(const int id, const int surf_offset);
793 :
794 : /**
795 : * Get the scalar03 solution at given GLL index
796 : * @param[in] id GLL index
797 : * @return scalar03 value at index
798 : */
799 : double get_scalar03(const int id, const int surf_offset);
800 :
801 : /**
802 : * Get the usrwrk zeroth slice at given GLL index
803 : * @param[in] id GLL index
804 : * @return zeroth slice of usrwrk value at index
805 : */
806 : double get_usrwrk00(const int id, const int surf_offset);
807 :
808 : /**
809 : * Get the usrwrk first slice at given GLL index
810 : * @param[in] id GLL index
811 : * @return first slice of usrwrk value at index
812 : */
813 : double get_usrwrk01(const int id, const int surf_offset);
814 :
815 : /**
816 : * Get the usrwrk second slice at given GLL index
817 : * @param[in] id GLL index
818 : * @return second slice of usrwrk value at index
819 : */
820 : double get_usrwrk02(const int id, const int surf_offset);
821 :
822 : /**
823 : * Get the temperature solution at given GLL index
824 : * @param[in] id GLL index
825 : * @return temperature value at index
826 : */
827 : double get_temperature(const int id, const int surf_offset);
828 :
829 : /**
830 : * Get the pressure solution at given GLL index
831 : * @param[in] id GLL index
832 : * @return pressure value at index
833 : */
834 : double get_pressure(const int id, const int surf_offset);
835 :
836 : /**
837 : * Return unity, for cases where the integrand or operator we are generalizing acts on 1
838 : * @param[in] id GLL index
839 : * @return unity
840 : */
841 : double get_unity(const int id, const int surf_offset);
842 :
843 : /**
844 : * Get the x-velocity at given GLL index
845 : * @param[in] id GLL index
846 : * @return x-velocity at index
847 : */
848 : double get_velocity_x(const int id, const int surf_offset);
849 :
850 : /**
851 : * Get the y-velocity at given GLL index
852 : * @param[in] id GLL index
853 : * @return y-velocity at index
854 : */
855 : double get_velocity_y(const int id, const int surf_offset);
856 :
857 : /**
858 : * Get the z-velocity at given GLL index
859 : * @param[in] id GLL index
860 : * @return z-velocity at index
861 : */
862 : double get_velocity_z(const int id, const int surf_offset);
863 :
864 : /**
865 : * Get the magnitude of the velocity solution at given GLL index
866 : * @param[in] id GLL index
867 : * @return velocity magnitude at index
868 : */
869 : double get_velocity(const int id, const int surf_offset);
870 :
871 : /**
872 : * Get the x-velocity squared at given GLL index
873 : * @param[in] id GLL index
874 : * @return square of x-velocity at index
875 : */
876 : double get_velocity_x_squared(const int id, const int surf_offset);
877 :
878 : /**
879 : * Get the y-velocity squared at given GLL index
880 : * @param[in] id GLL index
881 : * @return square of y-velocity at index
882 : */
883 : double get_velocity_y_squared(const int id, const int surf_offset);
884 :
885 : /**
886 : * Get the z-velocity squared at given GLL index
887 : * @param[in] id GLL index
888 : * @return square of z-velocity at index
889 : */
890 : double get_velocity_z_squared(const int id, const int surf_offset);
891 :
892 : /**
893 : * Initialize the characteristic scales for a nondimesional solution
894 : * @param[in] U reference velocity
895 : * @param[in] T reference temperature
896 : * @param[in] dT reference temperature range
897 : * @param[in] L reference length scale
898 : * @param[in] rho reference density
899 : * @param[in] Cp reference heat capacity
900 : * @param[in] s01 reference scalar01
901 : * @param[in] ds01 reference s01 range
902 : * @param[in] s02 reference scalar02
903 : * @param[in] ds02 reference s02 range
904 : * @param[in] s03 reference scalar03
905 : * @param[in] ds03 reference s03 range
906 : */
907 : void initializeDimensionalScales(const double U,
908 : const double T,
909 : const double dT,
910 : const double L,
911 : const double rho,
912 : const double Cp,
913 : const double s01,
914 : const double ds01,
915 : const double s02,
916 : const double ds02,
917 : const double s03,
918 : const double ds03);
919 :
920 : /**
921 : * \brief Return the reference divisor scale that defines the non-dimensional field
922 : *
923 : * All fields in NekRS are assumed non-dimensionalized according to the general form
924 : * f (non-dimensional) = (f - f_ref) / df
925 : *
926 : * so that to define a nondimensionalization requires two scales: the divisor scale
927 : * (df) and the additive scale (f_ref).
928 : *
929 : * @param[in] field physical interpretation of value
930 : * @param[out] value nondimensional divisor scale (df)
931 : */
932 : Real nondimensionalDivisor(const field::NekFieldEnum & field);
933 : Real nondimensionalDivisor(const field::NekWriteEnum & field);
934 :
935 : /**
936 : * All fields in NekRS are assumed non-dimensionalized according to the general form
937 : * f (non-dimensional) = (f - f_ref) / df
938 : *
939 : * so that to define a nondimensionalization requires two scales: the divisor scale
940 : * (df) and the additive scale (f_ref).
941 : *
942 : * @param[in] field physical interpretation of value
943 : * @param[out] value nondimensional additive scale (f_ref)
944 : *
945 : */
946 : Real nondimensionalAdditive(const field::NekFieldEnum & field);
947 : Real nondimensionalAdditive(const field::NekWriteEnum & field);
948 :
949 : /**
950 : * Get the reference length scale
951 : * @return reference length scale
952 : */
953 : double referenceLength();
954 :
955 : /**
956 : * Get the reference time scale
957 : * @return reference time scale
958 : */
959 : double referenceTime();
960 :
961 : /**
962 : * Get the reference area scale
963 : * @return reference area scale
964 : */
965 : double referenceArea();
966 :
967 : /**
968 : * Get the reference volume scale
969 : * @return reference volume scale
970 : */
971 : double referenceVolume();
972 :
973 : // useful concept from Stack Overflow for templating MPI calls
974 : template <typename T>
975 : MPI_Datatype resolveType();
976 :
977 : /**
978 : * Helper function for MPI_Allgatherv of results in NekRS
979 : * @param[in] base_counts once multiplied by 'multiplier', the number of counts on each rank
980 : * @param[in] input rank-local data
981 : * @param[out] output collected result
982 : * @param[in] multiplier constant multiplier to set on each count indicator
983 : */
984 : template <typename T>
985 : void
986 36379 : allgatherv(const std::vector<int> & base_counts, const T * input, T * output, const int multiplier)
987 : {
988 36379 : int * recvCounts = (int *)calloc(commSize(), sizeof(int));
989 36379 : int * displacement = (int *)calloc(commSize(), sizeof(int));
990 36379 : displacementAndCounts(base_counts, recvCounts, displacement, multiplier);
991 :
992 36379 : MPI_Allgatherv(input,
993 : recvCounts[commRank()],
994 : resolveType<T>(),
995 : output,
996 : (const int *)recvCounts,
997 : (const int *)displacement,
998 : resolveType<T>(),
999 : platform->comm.mpiComm());
1000 :
1001 36379 : free(recvCounts);
1002 36379 : free(displacement);
1003 36379 : }
1004 :
1005 : void initializeNekHostArrays();
1006 :
1007 : // Accessors for NekRS host arrays
1008 : std::tuple<dfloat *, dfloat *, dfloat *> host_xyz();
1009 : std::vector<dfloat> & host_U();
1010 : std::vector<dfloat> & host_P();
1011 : std::vector<dfloat> & host_S();
1012 : dfloat * host_wrk();
1013 :
1014 : mesh_t * createMesh2(mesh_t * _mesh, int Nc);
1015 :
1016 : } // end namespace nekrs
|
