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 "NekRSMesh.h"
22 : #include "libmesh/face_quad4.h"
23 : #include "libmesh/face_quad9.h"
24 : #include "libmesh/cell_hex8.h"
25 : #include "libmesh/cell_hex27.h"
26 : #include "inipp.hpp"
27 : #include "nekrs.hpp"
28 : #include "CardinalUtils.h"
29 : #include "VariadicTable.h"
30 :
31 : registerMooseObject("CardinalApp", NekRSMesh);
32 :
33 : InputParameters
34 1663 : NekRSMesh::validParams()
35 : {
36 1663 : InputParameters params = MooseMesh::validParams();
37 3326 : params.addParam<std::vector<int>>("boundary",
38 : "Boundary ID(s) through which nekRS will be coupled to MOOSE");
39 3326 : params.addParam<bool>("volume", false, "Whether the nekRS volume will be coupled to MOOSE");
40 3326 : params.addParam<bool>("exact", false, "Whether the mesh mirror is an exact replica of the NekRS mesh");
41 3326 : params.addParam<MooseEnum>(
42 3326 : "order", getNekOrderEnum(), "Order of the mesh interpolation between nekRS and MOOSE");
43 4989 : params.addRangeCheckedParam<Real>(
44 3326 : "scaling", 1.0, "scaling > 0.0", "Scaling factor to apply to the mesh");
45 3326 : params.addParam<unsigned int>("fluid_block_id", 0, "Subdomain ID to use for the fluid mesh mirror");
46 3326 : params.addParam<unsigned int>("solid_block_id", 1, "Subdomain ID to use for the solid mesh mirror");
47 1663 : params.addClassDescription(
48 : "Construct a mirror of the NekRS mesh in boundary and/or volume format");
49 1663 : return params;
50 0 : }
51 :
52 809 : NekRSMesh::NekRSMesh(const InputParameters & parameters)
53 : : MooseMesh(parameters),
54 809 : _volume(getParam<bool>("volume")),
55 2406 : _boundary(isParamValid("boundary") ? &getParam<std::vector<int>>("boundary") : nullptr),
56 1618 : _order(getParam<MooseEnum>("order").getEnum<order::NekOrderEnum>()),
57 1618 : _exact(getParam<bool>("exact")),
58 1618 : _scaling(getParam<Real>("scaling")),
59 1618 : _fluid_block_id(getParam<unsigned int>("fluid_block_id")),
60 1618 : _solid_block_id(getParam<unsigned int>("solid_block_id")),
61 809 : _n_surface_elems(0),
62 1618 : _n_volume_elems(0)
63 : {
64 809 : if (_exact && _order != order::first)
65 1 : mooseError("When building an exact mesh mirror, the 'order' must be FIRST!");
66 :
67 808 : if (!_boundary && !_volume)
68 0 : mooseError("This mesh requires at least 'volume = true' or a list of IDs in 'boundary'!");
69 :
70 808 : if (_boundary && _boundary->empty())
71 0 : paramError("boundary", "The length of 'boundary' must be greater than zero!");
72 :
73 : // see if NekRS's mesh even exists
74 808 : if (!nekrs::isInitialized())
75 1 : mooseError("This mesh can only be used with wrapped Nek cases! "
76 : "You need to change the problem type to NekRSProblem.");
77 :
78 807 : _nek_internal_mesh = nekrs::entireMesh();
79 :
80 807 : nekrs::initializeHostMeshParameters();
81 807 : nekrs::updateHostMeshParameters();
82 :
83 : // nekRS will only ever support 3-D meshes. Just to be sure that this remains
84 : // the case for future Cardinal developers, throw an error if the mesh isn't 3-D
85 : // (since this would affect how we construct the mesh here).
86 807 : int dimension = nekrs::dim();
87 807 : if (dimension != 3)
88 0 : mooseError("This mesh assumes that the nekRS mesh dimension is 3!\n\nYour mesh is "
89 0 : "dimension " +
90 0 : std::to_string(dimension) + ".");
91 :
92 : // if doing a JIT build, the boundary information does not exist yet
93 807 : if (!nekrs::buildOnly() && _boundary)
94 : {
95 : int first_invalid_id, n_boundaries;
96 393 : bool valid_ids = nekrs::validBoundaryIDs(*_boundary, first_invalid_id, n_boundaries);
97 :
98 393 : if (!valid_ids)
99 1 : mooseError("Invalid 'boundary' entry: ",
100 : first_invalid_id,
101 : "\n\n"
102 : "nekRS assumes the boundary IDs are ordered contiguously beginning at 1. "
103 : "For this problem, nekRS has ",
104 : n_boundaries,
105 : " boundaries. "
106 : "Did you enter a valid 'boundary'?");
107 : }
108 :
109 1612 : _corner_indices = nekrs::cornerGLLIndices(nekrs::entireMesh()->N, _exact);
110 806 : }
111 :
112 : void
113 25 : NekRSMesh::saveInitialVolMesh()
114 : {
115 : // save the initial mesh structure in case we are applying displacements
116 : // (which are additive to the initial mesh structure)
117 :
118 25 : long ngllpts = _nek_internal_mesh->Nelements * _nek_internal_mesh->Np;
119 :
120 25 : _initial_x.resize(ngllpts,0.0);
121 25 : _initial_y.resize(ngllpts,0.0);
122 25 : _initial_z.resize(ngllpts,0.0);
123 :
124 25 : memcpy(_initial_x.data(), _nek_internal_mesh->x, ngllpts * sizeof(double));
125 25 : memcpy(_initial_y.data(), _nek_internal_mesh->y, ngllpts * sizeof(double));
126 25 : memcpy(_initial_z.data(), _nek_internal_mesh->z, ngllpts * sizeof(double));
127 25 : }
128 :
129 : void
130 16 : NekRSMesh::initializePreviousDisplacements()
131 : {
132 16 : long disp_length = _volume ? _n_vertices_per_volume * _n_volume_elems
133 8 : : _n_vertices_per_surface * _n_surface_elems;
134 16 : _prev_disp_x.resize(disp_length,0.0);
135 16 : _prev_disp_y.resize(disp_length,0.0);
136 16 : _prev_disp_z.resize(disp_length,0.0);
137 16 : }
138 :
139 : void
140 797 : NekRSMesh::printMeshInfo() const
141 : {
142 797 : _console << "\nNekRS mesh mapping to MOOSE:" << std::endl;
143 : VariadicTable<std::string, int, std::string, int, int> vt(
144 4782 : {"", "Order", "Boundaries", "# Side Elems", "# Volume Elems"});
145 :
146 : std::vector<int> nek_bids;
147 4668 : for (int i = 1; i <= nekrs::NboundaryID(); ++i)
148 3871 : nek_bids.push_back(i);
149 :
150 1594 : vt.addRow("NekRS mesh",
151 : nekrs::polynomialOrder(),
152 1594 : Moose::stringify(nek_bids),
153 797 : _nek_n_surface_elems,
154 797 : _nek_n_volume_elems);
155 :
156 797 : std::string boundaries = "";
157 797 : if (_boundary)
158 774 : boundaries = Moose::stringify(*_boundary);
159 : else
160 820 : boundaries = Moose::stringify(nek_bids);
161 1594 : vt.addRow("NekRS mirror", _order + 1, boundaries, _n_surface_elems * _n_build_per_surface_elem,
162 797 : _n_volume_elems * _n_build_per_volume_elem);
163 :
164 797 : vt.print(_console);
165 797 : _console << std::endl;
166 2391 : }
167 :
168 : std::unique_ptr<MooseMesh>
169 44 : NekRSMesh::safeClone() const
170 : {
171 44 : return _app.getFactory().copyConstruct(*this);
172 : }
173 :
174 : int
175 1594 : NekRSMesh::numQuadraturePoints1D() const
176 : {
177 1594 : return _order + 2;
178 : }
179 :
180 : int
181 200 : NekRSMesh::nekNumQuadraturePoints1D() const
182 : {
183 200 : return _nek_polynomial_order + 1;
184 : }
185 :
186 : void
187 806 : NekRSMesh::initializeMeshParams()
188 : {
189 806 : _nek_polynomial_order = nekrs::polynomialOrder();
190 806 : _n_build_per_surface_elem = _exact ? _nek_polynomial_order * _nek_polynomial_order : 1;
191 806 : _n_build_per_volume_elem = _exact ? std::pow(_nek_polynomial_order, 3) : 1;
192 :
193 : /**
194 : * The libMesh face numbering for a 3-D hexagonal element is
195 : * 0
196 : * ^ 3
197 : * | /
198 : * o--------o
199 : * /: | / /|
200 : * / : |/ / |
201 : * / : / |
202 : * o--------o -|-> 2
203 : * 4<-|- o....|...o
204 : * | . | /
205 : * | . /| | /
206 : * |. / | |/
207 : * o--------o
208 : * / |
209 : * 1 5
210 : *
211 : * but for nekRS it is
212 : * 3
213 : * ^ 5
214 : * | /
215 : * o--------o
216 : * /: | / /|
217 : * / : |/ / |
218 : * / : / |
219 : * o--------o -|-> 2
220 : * 4<-|- o....|...o
221 : * | . | /
222 : * | . /| | /
223 : * |. / | |/
224 : * o--------o
225 : * / |
226 : * 0 1
227 :
228 : */
229 806 : _side_index = {1, 5, 2, 0, 4, 3};
230 :
231 806 : switch (_order)
232 : {
233 675 : case order::first:
234 675 : _n_vertices_per_surface = 4;
235 675 : _n_vertices_per_volume = 8;
236 :
237 : /** The libMesh node numbering for Quad 4 is
238 : *
239 : * 3 -- 2
240 : * | |
241 : * 0 -- 1
242 : *
243 : * but for nekRS it is
244 : *
245 : * 2 -- 3
246 : * | |
247 : * 0 -- 1
248 : **/
249 675 : _bnd_node_index = {0, 1, 3, 2};
250 :
251 : /** The libMesh node number for Hex 8 is
252 : * 3 2
253 : * o--------o
254 : * /: /|
255 : * / : / |
256 : * 0 / : 1 / |
257 : * o--------o |
258 : * | o....|...o 6
259 : * | .7 | /
260 : * | . | /
261 : * |. |/
262 : * o--------o
263 : * 4 5
264 : *
265 : * but for nekRS it is
266 : *
267 : * 6 7
268 : * o--------o
269 : * /: /|
270 : * / : / |
271 : * 2 / : 3 / |
272 : * o--------o |
273 : * | o....|...o 5
274 : * | .4 | /
275 : * | . | /
276 : * |. |/
277 : * o--------o
278 : * 0 1
279 : */
280 675 : _vol_node_index = {2, 3, 7, 6, 0, 1, 5, 4};
281 :
282 675 : break;
283 131 : case order::second:
284 131 : _n_vertices_per_surface = 9;
285 131 : _n_vertices_per_volume = 27;
286 :
287 : /** The libMesh node numbering for Quad 9 is
288 : *
289 : * 3 - 6 - 2
290 : * | |
291 : * 7 8 5
292 : * | |
293 : * 0 - 4 - 1
294 : *
295 : * but for nekRS it is
296 : *
297 : * 6 - 7 - 8
298 : * | |
299 : * 3 4 5
300 : * | |
301 : * 0 - 1 - 2
302 : **/
303 131 : _bnd_node_index = {0, 2, 8, 6, 1, 5, 7, 3, 4};
304 :
305 : /** The libMesh node numbering for Hex 27 is
306 : *
307 : *
308 : * 3 10 2
309 : * o--------------o--------------o
310 : * /: / /|
311 : * / : / / |
312 : * / : / / |
313 : * 11/ : 20/ 9/ |
314 : * o--------------o--------------o |
315 : * / : / /| |
316 : * / 15o / 23o / | 14o
317 : * / : / / | /|
318 : * 0/ : 8/ 1/ | / |
319 : * o--------------o--------------o | / |
320 : * | : | 26 | |/ |
321 : * | 24o : | o | 22o |
322 : * | : | 18 | /| |
323 : * | 7o....|.........o....|../.|....o
324 : * | . | | / | / 6
325 : * | . 21| 13|/ | /
326 : * 12 o--------------o--------------o | /
327 : * | . | | |/
328 : * | 19o | 25o | o
329 : * | . | | / 17
330 : * | . | | /
331 : * | . | | /
332 : * |. | |/
333 : * o--------------o--------------o
334 : * 4 16 5
335 : *
336 : * but for nekRS it is
337 : *
338 : * 24 25 26
339 : * o--------------o--------------o
340 : * /: / /|
341 : * / : / / |
342 : * / : / / |
343 : * 15/ : 16/ 17/ |
344 : * o--------------o--------------o |
345 : * / : / /| |
346 : * / 21o / 22o / | 23o
347 : * / : / / | /|
348 : * 6/ : 7/ 8/ | / |
349 : * o--------------o--------------o | / |
350 : * | : | 13 | |/ |
351 : * | 12o : | o | 14o |
352 : * | : | 19 | /| |
353 : * | 18o....|.........o....|../.|....o
354 : * | . | | / | / 20
355 : * | . 4| 5|/ | /
356 : * 3 o--------------o--------------o | /
357 : * | . | | |/
358 : * | 9o | 10o | o
359 : * | . | | / 11
360 : * | . | | /
361 : * | . | | /
362 : * |. | |/
363 : * o--------------o--------------o
364 : * 0 1 2
365 : */
366 262 : _vol_node_index = {6, 8, 26, 24, 0, 2, 20, 18, 7, 17, 25, 15, 3, 5,
367 131 : 23, 21, 1, 11, 19, 9, 16, 4, 14, 22, 12, 10, 13};
368 :
369 131 : break;
370 0 : default:
371 0 : mooseError("Unhandled 'NekOrderEnum' in 'NekRSMesh'!");
372 : }
373 806 : }
374 :
375 : void
376 0 : NekRSMesh::buildDummyMesh()
377 : {
378 : int e = 1;
379 0 : auto elem = new Quad4;
380 0 : elem->set_id() = e;
381 0 : elem->processor_id() = 0;
382 0 : _mesh->add_elem(elem);
383 :
384 : Point pt1(0.0, 0.0, 0.0);
385 : Point pt2(1.0, 0.0, 0.0);
386 : Point pt3(1.0, 1.0, 0.0);
387 : Point pt4(0.0, 1.0, 0.0);
388 :
389 0 : elem->set_node(0) = _mesh->add_point(pt1);
390 0 : elem->set_node(1) = _mesh->add_point(pt2);
391 0 : elem->set_node(2) = _mesh->add_point(pt3);
392 0 : elem->set_node(3) = _mesh->add_point(pt4);
393 :
394 0 : _mesh->prepare_for_use();
395 0 : }
396 :
397 : void
398 392 : NekRSMesh::storeBoundaryCoupling()
399 : {
400 392 : int rank = nekrs::commRank();
401 392 : int max_possible_surfaces = _nek_internal_mesh->NboundaryFaces;
402 :
403 392 : int * etmp = (int *)malloc(max_possible_surfaces * sizeof(int));
404 392 : int * ftmp = (int *)malloc(max_possible_surfaces * sizeof(int));
405 392 : int * ptmp = (int *)malloc(max_possible_surfaces * sizeof(int));
406 392 : int * btmp = (int *)malloc(max_possible_surfaces * sizeof(int));
407 392 : int * element = (int *)malloc(max_possible_surfaces * sizeof(int));
408 392 : int * face = (int *)malloc(max_possible_surfaces * sizeof(int));
409 392 : int * process = (int *)malloc(max_possible_surfaces * sizeof(int));
410 392 : int * boundary_id = (int *)malloc(max_possible_surfaces * sizeof(int));
411 :
412 : // number of faces on boundary of interest for this process
413 392 : int Nfaces = 0;
414 :
415 : int d = 0;
416 183743 : for (int i = 0; i < _nek_internal_mesh->Nelements; ++i)
417 : {
418 1283457 : for (int j = 0; j < _nek_internal_mesh->Nfaces; ++j)
419 : {
420 1100106 : int face_id = _nek_internal_mesh->EToB[i * _nek_internal_mesh->Nfaces + j];
421 :
422 1100106 : if (std::find(_boundary->begin(), _boundary->end(), face_id) != _boundary->end())
423 : {
424 39901 : Nfaces += 1;
425 :
426 39901 : etmp[d] = i;
427 39901 : ftmp[d] = j;
428 39901 : ptmp[d] = rank;
429 39901 : btmp[d] = face_id;
430 39901 : d++;
431 : }
432 : }
433 : }
434 :
435 : // gather all the boundary face counters and make available in N
436 392 : MPI_Allreduce(&Nfaces, &_n_surface_elems, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm);
437 392 : _boundary_coupling.n_faces = Nfaces;
438 392 : _boundary_coupling.total_n_faces = _n_surface_elems;
439 :
440 : // make available to all processes the number of faces owned by each process
441 392 : _boundary_coupling.counts.resize(nekrs::commSize());
442 392 : MPI_Allgather(
443 : &Nfaces, 1, MPI_INT, &_boundary_coupling.counts[0], 1, MPI_INT, platform->comm.mpiComm);
444 :
445 392 : int N_mirror_faces = Nfaces * _n_build_per_surface_elem;
446 392 : _boundary_coupling.mirror_counts.resize(nekrs::commSize());
447 392 : MPI_Allgather(
448 : &N_mirror_faces, 1, MPI_INT, &_boundary_coupling.mirror_counts[0], 1, MPI_INT, platform->comm.mpiComm);
449 :
450 : // compute the counts and displacements for face-based data exchange
451 392 : int * recvCounts = (int *)calloc(nekrs::commSize(), sizeof(int));
452 392 : int * displacement = (int *)calloc(nekrs::commSize(), sizeof(int));
453 392 : nekrs::displacementAndCounts(_boundary_coupling.counts, recvCounts, displacement, 1);
454 :
455 392 : _boundary_coupling.offset = displacement[rank];
456 :
457 392 : nekrs::allgatherv(_boundary_coupling.counts, etmp, element);
458 392 : nekrs::allgatherv(_boundary_coupling.counts, ftmp, face);
459 392 : nekrs::allgatherv(_boundary_coupling.counts, ptmp, process);
460 392 : nekrs::allgatherv(_boundary_coupling.counts, btmp, boundary_id);
461 :
462 467914 : for (int i = 0; i < max_possible_surfaces; ++i)
463 : {
464 467522 : _boundary_coupling.element.push_back(element[i]);
465 467522 : _boundary_coupling.face.push_back(face[i]);
466 467522 : _boundary_coupling.process.push_back(process[i]);
467 467522 : _boundary_coupling.boundary_id.push_back(boundary_id[i]);
468 : }
469 :
470 : freePointer(recvCounts);
471 : freePointer(displacement);
472 : freePointer(etmp);
473 : freePointer(ftmp);
474 : freePointer(ptmp);
475 : freePointer(btmp);
476 : freePointer(element);
477 : freePointer(face);
478 : freePointer(process);
479 : freePointer(boundary_id);
480 392 : }
481 :
482 : void
483 486 : NekRSMesh::storeVolumeCoupling()
484 : {
485 486 : int rank = nekrs::commRank();
486 :
487 486 : _volume_coupling.n_elems = _nek_internal_mesh->Nelements;
488 486 : MPI_Allreduce(
489 : &_volume_coupling.n_elems, &_n_volume_elems, 1, MPI_INT, MPI_SUM, platform->comm.mpiComm);
490 486 : _volume_coupling.total_n_elems = _n_volume_elems;
491 :
492 486 : _volume_coupling.counts.resize(nekrs::commSize());
493 486 : MPI_Allgather(&_volume_coupling.n_elems,
494 : 1,
495 : MPI_INT,
496 : &_volume_coupling.counts[0],
497 : 1,
498 : MPI_INT,
499 : platform->comm.mpiComm);
500 :
501 486 : _volume_coupling.mirror_counts.resize(nekrs::commSize());
502 486 : int N_mirror_elems = _volume_coupling.n_elems * _n_build_per_volume_elem;
503 486 : MPI_Allgather(&N_mirror_elems,
504 : 1,
505 : MPI_INT,
506 : &_volume_coupling.mirror_counts[0],
507 : 1,
508 : MPI_INT,
509 : platform->comm.mpiComm);
510 :
511 : // Save information regarding the volume mesh coupling in terms of the process-local
512 : // element IDs and process ownership; the 'tmp' arrays hold the rank-local data,
513 : // while the other arrays hold the result of the allgatherv
514 486 : int * etmp = (int *)malloc(_n_volume_elems * sizeof(int));
515 486 : int * ptmp = (int *)malloc(_n_volume_elems * sizeof(int));
516 486 : int * btmp = (int *)malloc(_n_volume_elems * _nek_internal_mesh->Nfaces * sizeof(int));
517 486 : int * element = (int *)malloc(_n_volume_elems * sizeof(int));
518 486 : int * process = (int *)malloc(_n_volume_elems * sizeof(int));
519 486 : int * boundary = (int *)malloc(_n_volume_elems * _nek_internal_mesh->Nfaces * sizeof(int));
520 :
521 255181 : for (int i = 0; i < _nek_internal_mesh->Nelements; ++i)
522 : {
523 254695 : etmp[i] = i;
524 254695 : ptmp[i] = rank;
525 :
526 1782865 : for (int j = 0; j < _nek_internal_mesh->Nfaces; ++j)
527 : {
528 1528170 : int id = i * _nek_internal_mesh->Nfaces + j;
529 1528170 : btmp[id] = _nek_internal_mesh->EToB[id];
530 : }
531 : }
532 :
533 486 : nekrs::allgatherv(_volume_coupling.counts, etmp, element, 1);
534 486 : nekrs::allgatherv(_volume_coupling.counts, ptmp, process, 1);
535 :
536 486 : int * ftmp = (int *)calloc(_n_volume_elems, sizeof(int));
537 486 : int * n_faces_on_boundary = (int *)calloc(_n_volume_elems, sizeof(int));
538 :
539 486 : int b_start = _boundary_coupling.offset;
540 10156 : for (int i = 0; i < _boundary_coupling.n_faces; ++i)
541 : {
542 9670 : int e = _boundary_coupling.element[b_start + i];
543 9670 : ftmp[e] += 1;
544 : }
545 :
546 486 : nekrs::allgatherv(_volume_coupling.counts, ftmp, n_faces_on_boundary, 1);
547 486 : nekrs::allgatherv(_volume_coupling.counts, btmp, boundary, _nek_internal_mesh->Nfaces);
548 :
549 4112136 : for (int i = 0; i < _n_volume_elems * _nek_internal_mesh->Nfaces; ++i)
550 4111650 : _volume_coupling.boundary.push_back(boundary[i]);
551 :
552 685761 : for (int i = 0; i < _n_volume_elems; ++i)
553 : {
554 685275 : _volume_coupling.element.push_back(element[i]);
555 685275 : _volume_coupling.process.push_back(process[i]);
556 685275 : _volume_coupling.n_faces_on_boundary.push_back(n_faces_on_boundary[i]);
557 : }
558 :
559 : freePointer(etmp);
560 : freePointer(ptmp);
561 : freePointer(ftmp);
562 : freePointer(btmp);
563 : freePointer(element);
564 : freePointer(process);
565 : freePointer(boundary);
566 : freePointer(n_faces_on_boundary);
567 486 : }
568 :
569 : void
570 806 : NekRSMesh::buildMesh()
571 : {
572 806 : if (nekrs::buildOnly())
573 : {
574 0 : buildDummyMesh();
575 0 : return;
576 : }
577 :
578 806 : _nek_n_surface_elems = nekrs::NboundaryFaces();
579 806 : _nek_n_volume_elems = nekrs::Nelements();
580 :
581 : // initialize the mesh mapping parameters that depend on order
582 806 : initializeMeshParams();
583 :
584 : // Loop through the mesh to establish a data structure (_boundary_coupling)
585 : // that holds the rank-local element ID, element-local face ID, and owning rank.
586 : // This data structure is used internally by nekRS during the transfer portion.
587 : // We must call this before the volume portion so that we can map the boundary
588 : // coupling to the volume coupling.
589 806 : if (_boundary)
590 392 : storeBoundaryCoupling();
591 :
592 : // Loop through the mesh to establish a data structure (_volume_coupling)
593 : // that holds the rank-local element ID and owning rank.
594 : // This data structure is used internally by nekRS during the transfer portion.
595 806 : if (_volume)
596 486 : storeVolumeCoupling();
597 :
598 806 : if (_boundary && !_volume)
599 320 : extractSurfaceMesh();
600 :
601 806 : if (_volume)
602 486 : extractVolumeMesh();
603 :
604 806 : addElems();
605 :
606 : // We're looking up the elements by id, so we can't let the ids get
607 : // renumbered.
608 : _mesh->allow_renumbering(false);
609 :
610 : // If we have a DistributedMesh then:
611 806 : if (!_mesh->is_replicated())
612 : {
613 : // we've already partitioned the elements to match the nekrs
614 : // mesh, and libMesh shouldn't try to improve on that. We won't
615 : // ever be doing any element deletion or coarsening, so we don't
616 : // even need libMesh's "critical" partitioning.
617 : _mesh->skip_partitioning(true);
618 :
619 : // But that means we have to update the partitioning metadata
620 : // ourselves
621 404 : _mesh->recalculate_n_partitions();
622 :
623 : // But, we haven't yet partitioned nodes, and if we tell libMesh
624 : // not to do that automatically then we need to do it manually
625 404 : libMesh::Partitioner::set_node_processor_ids(*_mesh);
626 : }
627 :
628 806 : _mesh->prepare_for_use();
629 : }
630 :
631 : void
632 806 : NekRSMesh::addElems()
633 : {
634 : BoundaryInfo & boundary_info = _mesh->get_boundary_info();
635 806 : auto nested_elems_on_face = nekrs::nestedElementsOnFace(_nek_polynomial_order);
636 :
637 788051 : for (int e = 0; e < _n_elems; e++)
638 : {
639 2746354 : for (int build = 0; build < _n_moose_per_nek; ++build)
640 : {
641 1959109 : auto elem = (this->*_new_elem)();
642 1959109 : elem->set_id() = e * _n_moose_per_nek + build;
643 1959109 : elem->processor_id() = (this->*_elem_processor_id)(e);
644 1959109 : _mesh->add_elem(elem);
645 :
646 : // add one point for each vertex of the face element
647 19202285 : for (int n = 0; n < _n_vertices_per_elem; n++)
648 : {
649 17243176 : int node = (*_node_index)[n];
650 :
651 17243176 : auto node_offset = (e * _n_moose_per_nek + build) * _n_vertices_per_elem + node;
652 17243176 : Point p(_x[node_offset], _y[node_offset], _z[node_offset]);
653 17243176 : p *= _scaling;
654 :
655 17243176 : auto node_ptr = _mesh->add_point(p);
656 17243176 : elem->set_node(n) = node_ptr;
657 : }
658 :
659 : // add sideset IDs to the mesh if we have volume coupling (this only adds the
660 : // sidesets associated with the coupling)
661 1959109 : if (_volume)
662 : {
663 11866645 : for (int f = 0; f < nekrs::Nfaces(); ++f)
664 : {
665 10171410 : int b_id = boundary_id(e, f);
666 10171410 : if (b_id != -1 /* NekRS's setting to indicate not on a sideset */)
667 : {
668 1054526 : if (_exact)
669 : {
670 605576 : auto faces = nested_elems_on_face[f];
671 605576 : if (!std::count(faces.begin(), faces.end(), build))
672 : continue;
673 : }
674 :
675 652110 : boundary_info.add_side(elem, _side_index[f], b_id);
676 : }
677 : }
678 :
679 1695235 : if (_phase[e * _n_moose_per_nek + build])
680 1248 : elem->subdomain_id() = _solid_block_id;
681 : else
682 1693987 : elem->subdomain_id() = _fluid_block_id;
683 : }
684 : }
685 : }
686 806 : }
687 :
688 : void
689 320 : NekRSMesh::faceVertices()
690 : {
691 320 : int n_vertices_in_mirror = _n_build_per_surface_elem * _n_surface_elems * _n_vertices_per_surface;
692 320 : double * x = (double *) malloc(n_vertices_in_mirror * sizeof(double));
693 320 : double * y = (double *) malloc(n_vertices_in_mirror * sizeof(double));
694 320 : double * z = (double *) malloc(n_vertices_in_mirror * sizeof(double));
695 :
696 320 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
697 320 : int rank = nekrs::commRank();
698 :
699 : mesh_t * mesh;
700 : int Nfp_mirror;
701 :
702 320 : if (_order == 0)
703 : {
704 : // For a first-order mesh mirror, we can take a shortcut and instead just fetch the
705 : // corner nodes. In this case, 'mesh' is no longer a custom-build mesh copy, but the
706 : // actual mesh for computation
707 298 : mesh = _nek_internal_mesh;
708 : Nfp_mirror = 4;
709 : }
710 : else
711 : {
712 : // Create a duplicate of the solution mesh, but with the desired order of the mesh
713 : // interpolation. Then we can just read the coordinates of the GLL points to find the libMesh
714 : // node positions. We only need to do this if the NekRS mesh is not already order 2.
715 22 : if (_nek_internal_mesh->N == 2)
716 : mesh = _nek_internal_mesh;
717 : else
718 : mesh =
719 13 : createMesh(platform->comm.mpiComm, _order + 1, 0, nrs->cht, *(nrs->kernelInfo));
720 :
721 22 : Nfp_mirror = mesh->Nfp;
722 : }
723 :
724 : // Allocate space for the coordinates that are on this rank
725 320 : int n_vertices_on_rank = _n_build_per_surface_elem * _boundary_coupling.n_faces * Nfp_mirror;
726 320 : double * xtmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
727 320 : double * ytmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
728 320 : double * ztmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
729 :
730 : int c = 0;
731 102290 : for (int k = 0; k < _boundary_coupling.total_n_faces; ++k)
732 : {
733 101970 : if (_boundary_coupling.process[k] == rank)
734 : {
735 30231 : int i = _boundary_coupling.element[k];
736 30231 : int j = _boundary_coupling.face[k];
737 30231 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
738 :
739 83094 : for (int build = 0; build < _n_build_per_surface_elem; ++build)
740 : {
741 274275 : for (int v = 0; v < Nfp_mirror; ++v, ++c)
742 : {
743 221412 : int vertex_offset = _order == 0 ? _corner_indices[build][v] : v;
744 221412 : int id = mesh->vmapM[offset + vertex_offset];
745 :
746 221412 : xtmp[c] = mesh->x[id];
747 221412 : ytmp[c] = mesh->y[id];
748 221412 : ztmp[c] = mesh->z[id];
749 : }
750 : }
751 : }
752 : }
753 :
754 320 : nekrs::allgatherv(_boundary_coupling.mirror_counts, xtmp, x, Nfp_mirror);
755 320 : nekrs::allgatherv(_boundary_coupling.mirror_counts, ytmp, y, Nfp_mirror);
756 320 : nekrs::allgatherv(_boundary_coupling.mirror_counts, ztmp, z, Nfp_mirror);
757 :
758 1084696 : for (int i = 0; i < n_vertices_in_mirror; ++i)
759 : {
760 1084376 : _x.push_back(x[i]);
761 1084376 : _y.push_back(y[i]);
762 1084376 : _z.push_back(z[i]);
763 : }
764 :
765 : freePointer(x);
766 : freePointer(y);
767 : freePointer(z);
768 : freePointer(xtmp);
769 : freePointer(ytmp);
770 : freePointer(ztmp);
771 320 : }
772 :
773 : void
774 486 : NekRSMesh::volumeVertices()
775 : {
776 : // nekRS has already performed a global operation such that all processes know the
777 : // toal number of volume elements and their phase
778 486 : int n_vertices_in_mirror = _n_build_per_volume_elem * _n_volume_elems * _n_vertices_per_volume;
779 486 : double * x = (double *) malloc(n_vertices_in_mirror * sizeof(double));
780 486 : double * y = (double *) malloc(n_vertices_in_mirror * sizeof(double));
781 486 : double * z = (double *) malloc(n_vertices_in_mirror * sizeof(double));
782 486 : double * p = (double *) malloc(_n_build_per_volume_elem * _n_volume_elems * sizeof(double));
783 :
784 486 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
785 486 : int rank = nekrs::commRank();
786 :
787 : mesh_t * mesh;
788 : int Np_mirror;
789 :
790 486 : if (_order == 0)
791 : {
792 : // For a first-order mesh mirror, we can take a shortcut and instead just fetch the
793 : // corner nodes. In this case, 'mesh' is no longer a custom-build mesh copy, but the
794 : // actual mesh for computation
795 377 : mesh = _nek_internal_mesh;
796 : Np_mirror = 8;
797 : }
798 : else
799 : {
800 : // Create a duplicate of the solution mesh, but with the desired order of the mesh interpolation.
801 : // Then we can just read the coordinates of the GLL points to find the libMesh node positions.
802 : // We only need to do this if the mesh is not already N = 2.
803 109 : if (_nek_internal_mesh->N == 2)
804 : mesh = _nek_internal_mesh;
805 : else
806 75 : mesh = createMesh(platform->comm.mpiComm, _order + 1, 0, nrs->cht, *(nrs->kernelInfo));
807 109 : Np_mirror = mesh->Np;
808 : }
809 :
810 : // Allocate space for the coordinates and phase that are on this rank
811 486 : int n_vertices_on_rank = _n_build_per_volume_elem * _volume_coupling.n_elems * Np_mirror;
812 486 : double * xtmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
813 486 : double * ytmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
814 486 : double * ztmp = (double *) malloc(n_vertices_on_rank * sizeof(double));
815 486 : double * ptmp = (double *) malloc(_n_build_per_volume_elem * _volume_coupling.n_elems * sizeof(double));
816 :
817 : int c = 0;
818 : int d = 0;
819 685761 : for (int k = 0; k < _volume_coupling.total_n_elems; ++k)
820 : {
821 685275 : if (_volume_coupling.process[k] == rank)
822 : {
823 254695 : int i = _volume_coupling.element[k];
824 254695 : int offset = i * mesh->Np;
825 :
826 678466 : for (int build = 0; build < _n_build_per_volume_elem; ++build)
827 : {
828 423771 : ptmp[d++] = i >= nekrs::flowMesh()->Nelements;
829 4353691 : for (int v = 0; v < Np_mirror; ++v, ++c)
830 : {
831 3929920 : int vertex_offset = _order == 0 ? _corner_indices[build][v] : v;
832 3929920 : int id = offset + vertex_offset;
833 :
834 3929920 : xtmp[c] = mesh->x[id];
835 3929920 : ytmp[c] = mesh->y[id];
836 3929920 : ztmp[c] = mesh->z[id];
837 : }
838 : }
839 : }
840 : }
841 :
842 486 : nekrs::allgatherv(_volume_coupling.mirror_counts, xtmp, x, Np_mirror);
843 486 : nekrs::allgatherv(_volume_coupling.mirror_counts, ytmp, y, Np_mirror);
844 486 : nekrs::allgatherv(_volume_coupling.mirror_counts, ztmp, z, Np_mirror);
845 486 : nekrs::allgatherv(_volume_coupling.mirror_counts, ptmp, p);
846 :
847 1695721 : for (int i = 0; i < _n_build_per_volume_elem * _n_volume_elems; ++i)
848 1695235 : _phase.push_back(p[i]);
849 :
850 16159286 : for (int i = 0; i < n_vertices_in_mirror; ++i)
851 : {
852 16158800 : _x.push_back(x[i]);
853 16158800 : _y.push_back(y[i]);
854 16158800 : _z.push_back(z[i]);
855 : }
856 :
857 : freePointer(x);
858 : freePointer(y);
859 : freePointer(z);
860 : freePointer(p);
861 : freePointer(xtmp);
862 : freePointer(ytmp);
863 : freePointer(ztmp);
864 : freePointer(ptmp);
865 486 : }
866 :
867 : void
868 320 : NekRSMesh::extractSurfaceMesh()
869 : {
870 : // Find the global vertex IDs that are on the _boundary. Note that nekRS performs a
871 : // global communciation here such that each nekRS process has knowledge of all the
872 : // boundary information.
873 320 : faceVertices();
874 :
875 320 : _new_elem = &NekRSMesh::boundaryElem;
876 320 : _n_elems = _n_surface_elems;
877 320 : _n_vertices_per_elem = _n_vertices_per_surface;
878 320 : _n_moose_per_nek = _n_build_per_surface_elem;
879 320 : _node_index = &_bnd_node_index;
880 320 : _elem_processor_id = &NekRSMesh::boundaryElemProcessorID;
881 320 : }
882 :
883 : void
884 486 : NekRSMesh::extractVolumeMesh()
885 : {
886 : // Find the global vertex IDs in the volume. Note that nekRS performs a
887 : // global communciation here such that each nekRS process has knowledge of all the
888 : // volume information.
889 486 : volumeVertices();
890 :
891 486 : _new_elem = &NekRSMesh::volumeElem;
892 486 : _n_elems = _n_volume_elems;
893 486 : _n_vertices_per_elem = _n_vertices_per_volume;
894 486 : _n_moose_per_nek = _n_build_per_volume_elem;
895 486 : _node_index = &_vol_node_index;
896 486 : _elem_processor_id = &NekRSMesh::volumeElemProcessorID;
897 486 : }
898 :
899 : Elem *
900 263874 : NekRSMesh::boundaryElem() const
901 : {
902 263874 : switch (_order)
903 : {
904 258098 : case order::first:
905 258098 : return new Quad4;
906 : break;
907 5776 : case order::second:
908 5776 : return new Quad9;
909 : break;
910 0 : default:
911 0 : mooseError("Unhandled 'NekOrderEnum' in 'NekRSMesh'!");
912 : }
913 : }
914 :
915 : Elem *
916 1695235 : NekRSMesh::volumeElem() const
917 : {
918 1695235 : switch (_order)
919 : {
920 1558555 : case order::first:
921 1558555 : return new Hex8;
922 : break;
923 136680 : case order::second:
924 136680 : return new Hex27;
925 : break;
926 0 : default:
927 0 : mooseError("Unhandled 'NekOrderEnum' in 'NekRSMesh'!");
928 : }
929 : }
930 :
931 : int
932 263874 : NekRSMesh::boundaryElemProcessorID(const int elem_id)
933 : {
934 263874 : return _boundary_coupling.processor_id(elem_id);
935 : }
936 :
937 : int
938 1695235 : NekRSMesh::volumeElemProcessorID(const int elem_id)
939 : {
940 1695235 : return _volume_coupling.processor_id(elem_id);
941 : }
942 :
943 : int
944 10171410 : NekRSMesh::boundary_id(const int elem_id, const int face_id)
945 : {
946 10171410 : return _volume_coupling.boundary[elem_id * _nek_internal_mesh->Nfaces + face_id];
947 : }
948 :
949 : int
950 1451696 : NekRSMesh::facesOnBoundary(const int elem_id) const
951 : {
952 1451696 : return _volume_coupling.n_faces_on_boundary[elem_id];
953 : }
954 :
955 : void
956 537600 : NekRSMesh::updateDisplacement(const int e, const double *src, const field::NekWriteEnum field)
957 : {
958 537600 : int nsrc = _volume? _n_vertices_per_volume : _n_vertices_per_surface;
959 537600 : int offset = e * nsrc;
960 :
961 537600 : switch (field)
962 : {
963 179200 : case field::x_displacement:
964 179200 : memcpy(&_prev_disp_x[offset], src, nsrc * sizeof(double));
965 179200 : break;
966 179200 : case field::y_displacement:
967 179200 : memcpy(&_prev_disp_y[offset], src, nsrc * sizeof(double));
968 179200 : break;
969 179200 : case field::z_displacement:
970 179200 : memcpy(&_prev_disp_z[offset], src, nsrc * sizeof(double));
971 179200 : break;
972 0 : default:
973 0 : throw std::runtime_error("Unhandled NekWriteEnum in NekRSMesh::copyToDisplacement!\n");
974 : }
975 537600 : }
976 : #endif
|
