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