Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PETScDMDAMesh.h"
11 :
12 : #include "libmesh/mesh_generation.h"
13 : #include "libmesh/string_to_enum.h"
14 : #include "libmesh/periodic_boundaries.h"
15 : #include "libmesh/periodic_boundary_base.h"
16 : #include "libmesh/unstructured_mesh.h"
17 : #include "libmesh/partitioner.h"
18 : #include "libmesh/metis_csr_graph.h"
19 : #include "libmesh/edge_edge2.h"
20 : #include "libmesh/edge_edge3.h"
21 : #include "libmesh/edge_edge4.h"
22 : #include "libmesh/mesh_communication.h"
23 : #include "libmesh/remote_elem.h"
24 : #include "libmesh/face_quad4.h"
25 : #include "libmesh/cell_hex8.h"
26 : #include "libmesh/petsc_solver_exception.h"
27 : #include "ExternalPetscSolverApp.h"
28 : // C++ includes
29 : #include <cmath> // provides round, not std::round (see http://www.cplusplus.com/reference/cmath/round/)
30 :
31 : registerMooseObject("ExternalPetscSolverApp", PETScDMDAMesh);
32 :
33 : InputParameters
34 740 : PETScDMDAMesh::validParams()
35 : {
36 740 : InputParameters params = MooseMesh::validParams();
37 :
38 740 : params.addClassDescription("Create a square mesh from PETSc DMDA.");
39 :
40 : // This mesh is always distributed
41 1480 : params.set<MooseEnum>("parallel_type") = "DISTRIBUTED";
42 :
43 740 : return params;
44 0 : }
45 :
46 370 : PETScDMDAMesh::PETScDMDAMesh(const InputParameters & parameters) : MooseMesh(parameters)
47 : {
48 : // Get TS from ExternalPetscSolverApp
49 370 : ExternalPetscSolverApp * petsc_app = dynamic_cast<ExternalPetscSolverApp *>(&_app);
50 :
51 370 : if (petsc_app && petsc_app->getPetscTS())
52 : // Retrieve mesh from TS
53 370 : LibmeshPetscCall(TSGetDM(petsc_app->getPetscTS(), &_dmda));
54 : else
55 0 : mooseError(" PETSc external solver TS does not exist or this is not a petsc external solver");
56 370 : }
57 :
58 : std::unique_ptr<MooseMesh>
59 0 : PETScDMDAMesh::safeClone() const
60 : {
61 0 : return _app.getFactory().copyConstruct(*this);
62 : }
63 :
64 : inline dof_id_type
65 : node_id_Quad4(const dof_id_type nx,
66 : const dof_id_type /*ny*/,
67 : const dof_id_type i,
68 : const dof_id_type j,
69 : const dof_id_type /*k*/)
70 :
71 : {
72 : // Transform a grid coordinate (i, j) to its global node ID
73 : // This match what PETSc does
74 61200 : return i + j * (nx + 1);
75 : }
76 :
77 : void
78 15300 : add_element_Quad4(DM da,
79 : const dof_id_type nx,
80 : const dof_id_type ny,
81 : const dof_id_type i,
82 : const dof_id_type j,
83 : const dof_id_type elem_id,
84 : const processor_id_type pid,
85 : MeshBase & mesh)
86 : {
87 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
88 :
89 : // Mx: number of grid points in x direction for all processors
90 : // My: number of grid points in y direction for all processors
91 : // xp: number of processors in x direction
92 : // yp: number of processors in y direction
93 : PetscInt Mx, My, xp, yp;
94 15300 : LibmeshPetscCallA(mesh.comm().get(),
95 : DMDAGetInfo(da,
96 : PETSC_IGNORE,
97 : &Mx,
98 : &My,
99 : PETSC_IGNORE,
100 : &xp,
101 : &yp,
102 : PETSC_IGNORE,
103 : PETSC_IGNORE,
104 : PETSC_IGNORE,
105 : PETSC_IGNORE,
106 : PETSC_IGNORE,
107 : PETSC_IGNORE,
108 : PETSC_IGNORE));
109 :
110 : const PetscInt *lx, *ly;
111 : PetscInt *lxo, *lyo;
112 : // PETSc-3.8.x or older use PetscDataType
113 : #if PETSC_VERSION_LESS_THAN(3, 9, 0)
114 : LibmeshPetscCallA(mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, PETSC_INT, &lxo));
115 : #else
116 : // PETSc-3.9.x or newer use MPI_DataType
117 15300 : LibmeshPetscCallA(mesh.comm().get(), DMGetWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
118 : #endif
119 :
120 : // Gets the ranges of indices in the x, y and z direction that are owned by each process
121 : // Ranges here are different from what we have in Mat and Vec.
122 : // It means how many points each processor holds
123 15300 : LibmeshPetscCallA(mesh.comm().get(), DMDAGetOwnershipRanges(da, &lx, &ly, NULL));
124 15300 : lxo[0] = 0;
125 33000 : for (PetscInt i = 0; i < xp; i++)
126 17700 : lxo[i + 1] = lxo[i] + lx[i];
127 :
128 15300 : lyo = lxo + xp + 1;
129 15300 : lyo[0] = 0;
130 43800 : for (PetscInt i = 0; i < yp; i++)
131 28500 : lyo[i + 1] = lyo[i] + ly[i];
132 :
133 : // Try to calculate processor-grid coordinate (xpid, ypid)
134 : PetscInt xpid, ypid, xpidplus, ypidplus;
135 : // Finds integer in a sorted array of integers
136 : // Loc: the location if found, otherwise -(slot+1)
137 : // where slot is the place the value would go
138 15300 : LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(i, xp + 1, lxo, &xpid));
139 :
140 15300 : xpid = xpid < 0 ? -xpid - 1 - 1 : xpid;
141 :
142 15300 : LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(i + 1, xp + 1, lxo, &xpidplus));
143 :
144 15300 : xpidplus = xpidplus < 0 ? -xpidplus - 1 - 1 : xpidplus;
145 :
146 15300 : LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(j, yp + 1, lyo, &ypid));
147 :
148 15300 : ypid = ypid < 0 ? -ypid - 1 - 1 : ypid;
149 :
150 15300 : LibmeshPetscCallA(mesh.comm().get(), PetscFindInt(j + 1, yp + 1, lyo, &ypidplus));
151 :
152 15300 : ypidplus = ypidplus < 0 ? -ypidplus - 1 - 1 : ypidplus;
153 : #if PETSC_VERSION_LESS_THAN(3, 9, 0)
154 : LibmeshPetscCallA(mesh.comm().get(), DMRestoreWorkArray(da, xp + yp + 2, PETSC_INT, &lxo));
155 : #else
156 15300 : LibmeshPetscCallA(mesh.comm().get(), DMRestoreWorkArray(da, xp + yp + 2, MPIU_INT, &lxo));
157 : #endif
158 :
159 : // Bottom Left
160 15300 : auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
161 : node_id_Quad4(nx, 0, i, j, 0));
162 : node0_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j, 0));
163 15300 : node0_ptr->set_id() = node0_ptr->unique_id();
164 : // xpid + ypid * xp is the global processor ID
165 15300 : node0_ptr->processor_id() = xpid + ypid * xp;
166 :
167 : // Bottom Right
168 : auto node1_ptr =
169 15300 : mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j) / ny, 0),
170 : node_id_Quad4(nx, 0, i + 1, j, 0));
171 : node1_ptr->set_unique_id(node_id_Quad4(nx, 0, i + 1, j, 0));
172 15300 : node1_ptr->set_id() = node1_ptr->unique_id();
173 15300 : node1_ptr->processor_id() = xpidplus + ypid * xp;
174 :
175 : // Top Right
176 : auto node2_ptr =
177 15300 : mesh.add_point(Point(static_cast<Real>(i + 1) / nx, static_cast<Real>(j + 1) / ny, 0),
178 : node_id_Quad4(nx, 0, i + 1, j + 1, 0));
179 : node2_ptr->set_unique_id(node_id_Quad4(nx, 0, i + 1, j + 1, 0));
180 15300 : node2_ptr->set_id() = node2_ptr->unique_id();
181 15300 : node2_ptr->processor_id() = xpidplus + ypidplus * xp;
182 :
183 : // Top Left
184 : auto node3_ptr =
185 15300 : mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j + 1) / ny, 0),
186 : node_id_Quad4(nx, 0, i, j + 1, 0));
187 : node3_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j + 1, 0));
188 15300 : node3_ptr->set_id() = node3_ptr->unique_id();
189 15300 : node3_ptr->processor_id() = xpid + ypidplus * xp;
190 :
191 : // New an element and attach four nodes to it
192 15300 : Elem * elem = new Quad4;
193 : elem->set_id(elem_id);
194 15300 : elem->processor_id() = pid;
195 : // Make sure our unique_id doesn't overlap any nodes'
196 15300 : elem->set_unique_id(elem_id + (nx + 1) * (ny + 1));
197 15300 : elem = mesh.add_elem(elem);
198 15300 : elem->set_node(0, node0_ptr);
199 15300 : elem->set_node(1, node1_ptr);
200 15300 : elem->set_node(2, node2_ptr);
201 15300 : elem->set_node(3, node3_ptr);
202 :
203 : // Bottom
204 15300 : if (j == 0)
205 1530 : boundary_info.add_side(elem, 0, 0);
206 :
207 : // Right
208 15300 : if (i == nx - 1)
209 1530 : boundary_info.add_side(elem, 1, 1);
210 :
211 : // Top
212 15300 : if (j == ny - 1)
213 1530 : boundary_info.add_side(elem, 2, 2);
214 :
215 : // Left
216 15300 : if (i == 0)
217 1530 : boundary_info.add_side(elem, 3, 3);
218 15300 : }
219 :
220 : void
221 345 : set_boundary_names_Quad4(BoundaryInfo & boundary_info)
222 : {
223 345 : boundary_info.sideset_name(0) = "bottom";
224 345 : boundary_info.sideset_name(1) = "right";
225 345 : boundary_info.sideset_name(2) = "top";
226 345 : boundary_info.sideset_name(3) = "left";
227 345 : }
228 :
229 : void
230 0 : get_indices_Quad4(const dof_id_type nx,
231 : const dof_id_type /*ny*/,
232 : const dof_id_type elem_id,
233 : dof_id_type & i,
234 : dof_id_type & j)
235 : {
236 0 : i = elem_id % nx;
237 0 : j = (elem_id - i) / nx;
238 0 : }
239 :
240 : dof_id_type
241 0 : elem_id_Quad4(const dof_id_type nx,
242 : const dof_id_type /*nx*/,
243 : const dof_id_type i,
244 : const dof_id_type j,
245 : const dof_id_type /*k*/)
246 : {
247 0 : return (j * nx) + i;
248 : }
249 :
250 : inline void
251 0 : get_neighbors_Quad4(const dof_id_type nx,
252 : const dof_id_type ny,
253 : const dof_id_type i,
254 : const dof_id_type j,
255 : std::vector<dof_id_type> & neighbors)
256 : {
257 : std::fill(neighbors.begin(), neighbors.end(), Elem::invalid_id);
258 :
259 : // Bottom
260 0 : if (j != 0)
261 0 : neighbors[0] = elem_id_Quad4(nx, 0, i, j - 1, 0);
262 :
263 : // Right
264 0 : if (i != nx - 1)
265 0 : neighbors[1] = elem_id_Quad4(nx, 0, i + 1, j, 0);
266 :
267 : // Top
268 0 : if (j != ny - 1)
269 0 : neighbors[2] = elem_id_Quad4(nx, 0, i, j + 1, 0);
270 :
271 : // Left
272 0 : if (i != 0)
273 0 : neighbors[3] = elem_id_Quad4(nx, 0, i - 1, j, 0);
274 0 : }
275 :
276 : void
277 0 : get_ghost_neighbors_Quad4(const dof_id_type nx,
278 : const dof_id_type ny,
279 : const MeshBase & mesh,
280 : std::set<dof_id_type> & ghost_elems)
281 : {
282 : auto & boundary_info = mesh.get_boundary_info();
283 :
284 : dof_id_type i, j;
285 :
286 0 : std::vector<dof_id_type> neighbors(4);
287 :
288 0 : for (auto elem_ptr : mesh.element_ptr_range())
289 : {
290 0 : for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
291 : {
292 : // No current neighbor
293 0 : if (!elem_ptr->neighbor_ptr(s))
294 : {
295 : // Not on a boundary
296 0 : if (!boundary_info.n_boundary_ids(elem_ptr, s))
297 : {
298 : auto elem_id = elem_ptr->id();
299 :
300 0 : get_indices_Quad4(nx, 0, elem_id, i, j);
301 :
302 0 : get_neighbors_Quad4(nx, ny, i, j, neighbors);
303 :
304 0 : ghost_elems.insert(neighbors[s]);
305 : }
306 : }
307 : }
308 0 : }
309 0 : }
310 :
311 : void
312 0 : add_node_Qua4(dof_id_type nx,
313 : dof_id_type ny,
314 : dof_id_type i,
315 : dof_id_type j,
316 : processor_id_type pid,
317 : MeshBase & mesh)
318 : {
319 : // Bottom Left
320 0 : auto node0_ptr = mesh.add_point(Point(static_cast<Real>(i) / nx, static_cast<Real>(j) / ny, 0),
321 : node_id_Quad4(nx, 0, i, j, 0));
322 : node0_ptr->set_unique_id(node_id_Quad4(nx, 0, i, j, 0));
323 0 : node0_ptr->processor_id() = pid;
324 0 : }
325 :
326 : void
327 345 : build_cube_Quad4(UnstructuredMesh & mesh, DM da)
328 : {
329 : const auto pid = mesh.comm().rank();
330 :
331 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
332 : // xs: start grid point (not element) index on local in x direction
333 : // ys: start grid point index on local in y direction
334 : // xm: number of grid points owned by the local processor in x direction
335 : // ym: number of grid points owned by the local processor in y direction
336 : // Mx: number of grid points on all processors in x direction
337 : // My: number of grid points on all processors in y direction
338 : // xp: number of processor cores in x direction
339 : // yp: number of processor cores in y direction
340 : PetscInt xs, ys, xm, ym, Mx, My, xp, yp;
341 :
342 : /* Get local grid boundaries */
343 345 : LibmeshPetscCallA(mesh.comm().get(),
344 : DMDAGetCorners(da, &xs, &ys, PETSC_IGNORE, &xm, &ym, PETSC_IGNORE));
345 345 : LibmeshPetscCallA(mesh.comm().get(),
346 : DMDAGetInfo(da,
347 : PETSC_IGNORE,
348 : &Mx,
349 : &My,
350 : PETSC_IGNORE,
351 : &xp,
352 : &yp,
353 : PETSC_IGNORE,
354 : PETSC_IGNORE,
355 : PETSC_IGNORE,
356 : PETSC_IGNORE,
357 : PETSC_IGNORE,
358 : PETSC_IGNORE,
359 : PETSC_IGNORE));
360 :
361 2292 : for (PetscInt j = ys; j < ys + ym; j++)
362 20460 : for (PetscInt i = xs; i < xs + xm; i++)
363 : {
364 : // We loop over grid points, but we are
365 : // here building elements. So that we just
366 : // simply skip the first x and y points since the
367 : // number of grid ponts is one more than
368 : // the number of grid elements
369 18513 : if (!i || !j)
370 3213 : continue;
371 :
372 15300 : dof_id_type ele_id = (i - 1) + (j - 1) * (Mx - 1);
373 :
374 15300 : add_element_Quad4(da, Mx - 1, My - 1, i - 1, j - 1, ele_id, pid, mesh);
375 : }
376 :
377 : // If there is no element at the given processor
378 : // We need to manually add all mesh nodes
379 345 : if ((ys == 0 && ym == 1) || (xs == 0 && xm == 1))
380 0 : for (PetscInt j = ys; j < ys + ym; j++)
381 0 : for (PetscInt i = xs; i < xs + xm; i++)
382 0 : add_node_Qua4(Mx, My, i, j, pid, mesh);
383 :
384 : // Need to link up the local elements before we can know what's missing
385 345 : mesh.find_neighbors();
386 :
387 345 : mesh.find_neighbors(true);
388 :
389 : // Set RemoteElem neighbors
390 31290 : for (auto & elem_ptr : mesh.element_ptr_range())
391 76500 : for (unsigned int s = 0; s < elem_ptr->n_sides(); s++)
392 61200 : if (!elem_ptr->neighbor_ptr(s) && !boundary_info.n_boundary_ids(elem_ptr, s))
393 3465 : elem_ptr->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
394 :
395 345 : set_boundary_names_Quad4(boundary_info);
396 : // Already partitioned!
397 : mesh.skip_partitioning(true);
398 :
399 : // We've already set our own unique_id values; now make sure that
400 : // future mesh modifications use subsequent values.
401 345 : mesh.set_next_unique_id(Mx * My + (Mx - 1) * (My - 1));
402 :
403 : // No need to renumber or find neighbors - done did it.
404 : mesh.allow_renumbering(false);
405 : mesh.allow_find_neighbors(false);
406 345 : mesh.prepare_for_use();
407 : mesh.allow_find_neighbors(true);
408 345 : }
409 :
410 : void
411 345 : PETScDMDAMesh::buildMesh()
412 : {
413 : // Reference to the libmesh mesh
414 345 : MeshBase & mesh = getMesh();
415 :
416 345 : build_cube_Quad4(dynamic_cast<UnstructuredMesh &>(mesh), _dmda);
417 345 : }
|
