PetIBM 0.5.4
Toolbox and applications of the immersed-boundary method for distributed-memory architectures
main.cpp
Go to the documentation of this file.
1
10#include <iomanip>
11#include <sstream>
12
13#include <petsc.h>
14#include <petscviewerhdf5.h>
15#include <yaml-cpp/yaml.h>
16
17#include <petibm/boundary.h>
18#include <petibm/io.h>
19#include <petibm/mesh.h>
20#include <petibm/parser.h>
21#include <petibm/solution.h>
22#include <petibm/type.h>
23
43PetscErrorCode initVorticityMesh(const petibm::type::Mesh &mesh,
44 petibm::type::IntVec2D &n,
45 petibm::type::RealVec3D &coord,
46 std::vector<DM> &wDMs,
47 std::vector<std::string> &names,
48 std::vector<Vec> &vecs);
49
50PetscErrorCode calculateVorticity2D(const petibm::type::Mesh &mesh,
51 const petibm::type::Boundary &bds,
52 const petibm::type::Solution &soln,
53 const std::vector<DM> &wDm,
54 std::vector<Vec> &w);
55
56PetscErrorCode calculateVorticity3D(const petibm::type::Mesh &mesh,
57 const petibm::type::Boundary &bds,
58 const petibm::type::Solution &soln,
59 const std::vector<DM> &wDm,
60 std::vector<Vec> &w);
61
62int main(int argc, char **argv)
63{
64 PetscErrorCode ierr;
65
66 ierr = PetscInitialize(&argc, &argv, nullptr, nullptr); CHKERRQ(ierr);
67
68 YAML::Node setting;
69 petibm::type::Mesh mesh;
70 petibm::type::Boundary bc;
71 petibm::type::Solution soln;
72
73 PetscInt bg;
74 PetscInt ed;
75 PetscInt step;
76
77 petibm::type::IntVec2D wN;
78 petibm::type::RealVec3D wCoord;
79 std::vector<DM> wDMs;
80 std::vector<Vec> w;
81 std::vector<std::string> wNames;
82
83 PetscBool isSet;
84
85 // read the config.yaml
86 ierr = petibm::parser::getSettings(setting); CHKERRQ(ierr);
87
88 // initialize mesh, bc, soln
89 ierr = petibm::mesh::createMesh(PETSC_COMM_WORLD, setting, mesh);
90 CHKERRQ(ierr);
91 ierr = petibm::boundary::createBoundary(mesh, setting, bc); CHKERRQ(ierr);
92 ierr = petibm::solution::createSolution(mesh, soln); CHKERRQ(ierr);
93
94 // initialize vorticity mesh
95 ierr = initVorticityMesh(mesh, wN, wCoord, wDMs, wNames, w); CHKERRQ(ierr);
96
97 // output grid coordinates to existing gird.h5
98 if (mesh->mpiRank == 0)
99 {
100 for (unsigned int i = 0; i < wNames.size(); ++i)
101 {
102 ierr = petibm::io::writeHDF5Vecs(
103 PETSC_COMM_SELF,
104 setting["output"].as<std::string>() + "/grid.h5", wNames[i],
105 {"x", "y", "z"}, wCoord[i], FILE_MODE_APPEND); CHKERRQ(ierr);
106 }
107 }
108 ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
109
110 // get the range of solutions that are going to be calculated
111 ierr = PetscOptionsGetInt(nullptr, nullptr, "-bg", &bg, &isSet);
112 CHKERRQ(ierr);
113 if (!isSet) bg = setting["parameters"]["startStep"].as<PetscInt>(0);
114 ierr = PetscOptionsGetInt(nullptr, nullptr, "-ed", &ed, &isSet);
115 CHKERRQ(ierr);
116 if (!isSet) ed = bg + setting["parameters"]["nt"].as<PetscInt>();
117 ierr = PetscOptionsGetInt(nullptr, nullptr, "-step", &step, &isSet);
118 CHKERRQ(ierr);
119 if (!isSet) step = setting["parameters"]["nsave"].as<PetscInt>();
120
121 // start calculating vorticity
122 ierr = PetscPrintf(PETSC_COMM_WORLD,
123 "====================================\n"
124 "Start to calculate vorticity fields!\n"
125 "====================================\n\n");
126 CHKERRQ(ierr);
127
128 ierr = PetscPrintf(PETSC_COMM_WORLD,
129 "Beginning: %d\n"
130 "End: %d\n"
131 "Step: %d\n\n",
132 bg, ed, step); CHKERRQ(ierr);
133
134 for (int i = bg; i <= ed; i += step)
135 {
136 ierr = PetscPrintf(PETSC_COMM_WORLD,
137 "Calculating vorticity fields for time step %d ... ",
138 i); CHKERRQ(ierr);
139
140 // read solution
141 std::stringstream ss;
142 ss << setting["output"].as<std::string>() << "/" << std::setfill('0')
143 << std::setw(7) << i << ".h5";
144 ierr = soln->read(ss.str()); CHKERRQ(ierr);
145
146 // calculate values at ghost points based on the current solution
147 ierr = bc->setGhostICs(soln); CHKERRQ(ierr);
148
149 // calculate vorticity
150 if (mesh->dim == 2)
151 {
152 ierr = calculateVorticity2D(mesh, bc, soln, wDMs, w);
153 CHKERRQ(ierr);
154 }
155 else
156 {
157 ierr = calculateVorticity3D(mesh, bc, soln, wDMs, w);
158 CHKERRQ(ierr);
159 }
160
161 // append to existing HDF5 file
162 ierr = petibm::io::writeHDF5Vecs(mesh->comm, ss.str(), "/", wNames, w,
163 FILE_MODE_APPEND); CHKERRQ(ierr);
164
165 ierr = PetscPrintf(PETSC_COMM_WORLD, "done.\n"); CHKERRQ(ierr);
166 }
167
168 // destroy everything from PETSc
169 for (unsigned int f = 0; f < w.size(); ++f)
170 {
171 ierr = VecDestroy(&w[f]); CHKERRQ(ierr);
172 ierr = DMDestroy(&wDMs[f]); CHKERRQ(ierr);
173 }
174
175 // manually destroy PETSc objects inside the PetIBM objects
176 ierr = soln->destroy(); CHKERRQ(ierr);
177 ierr = bc->destroy(); CHKERRQ(ierr);
178 ierr = mesh->destroy(); CHKERRQ(ierr);
179
180 ierr = PetscFinalize(); CHKERRQ(ierr);
181
182 return 0;
183} // main
184
185PetscErrorCode calculateVorticity2D(const petibm::type::Mesh &mesh,
186 const petibm::type::Boundary &bds,
187 const petibm::type::Solution &soln,
188 const std::vector<DM> &wDm,
189 std::vector<Vec> &w)
190{
191 PetscErrorCode ierr;
192
193 PetscFunctionBeginUser;
194
195 std::vector<Vec> uLocals;
196 std::vector<PetscReal **> uLocalArrys;
197 PetscReal **wArray;
198 DMDALocalInfo info;
199
200 // get local info for vorticity meshes
201 ierr = DMDAGetLocalInfo(wDm[0], &info); CHKERRQ(ierr);
202
203 // resize the STL vectors
204 uLocals.resize(mesh->dim);
205 uLocalArrys.resize(mesh->dim);
206
207 // create local Vecs
208 for (PetscInt f = 0; f < mesh->dim; ++f)
209 {
210 ierr = DMCreateLocalVector(mesh->da[f], &uLocals[f]); CHKERRQ(ierr);
211 }
212
213 // scatter from global Vecs to local Vecs
214 ierr = DMCompositeScatterArray(mesh->UPack, soln->UGlobal, uLocals.data());
215 CHKERRQ(ierr);
216
217 // copy the values at ghost points from Boundary instance to local Vecs
218 ierr = bds->copyValues2LocalVecs(uLocals); CHKERRQ(ierr);
219
220 // get read-only raw array from local Vecs
221 for (PetscInt f = 0; f < mesh->dim; ++f)
222 {
223 ierr = DMDAVecGetArrayRead(mesh->da[f], uLocals[f], &uLocalArrys[f]);
224 CHKERRQ(ierr);
225 }
226
227 // get raw array from vorticity Vecs
228 ierr = DMDAVecGetArray(wDm[0], w[0], &wArray); CHKERRQ(ierr);
229
230 // calculate z-vorticity
231 for (PetscInt j = info.ys; j < info.ys + info.ym; ++j)
232 {
233 for (PetscInt i = info.xs; i < info.xs + info.xm; ++i)
234 {
235 wArray[j][i] =
236 (uLocalArrys[1][j - 1][i] - uLocalArrys[1][j - 1][i - 1]) /
237 (mesh->coord[1][0][i] - mesh->coord[1][0][i - 1]) -
238 (uLocalArrys[0][j][i - 1] - uLocalArrys[0][j - 1][i - 1]) /
239 (mesh->coord[0][1][j] - mesh->coord[0][1][j - 1]);
240 }
241 }
242
243 // return the ownership of vorticity raw array
244 ierr = DMDAVecRestoreArray(wDm[0], w[0], &wArray); CHKERRQ(ierr);
245
246 // return the ownership of local velocity array and destroy Vecs
247 for (PetscInt f = 0; f < mesh->dim; ++f)
248 {
249 ierr =
250 DMDAVecRestoreArrayRead(mesh->da[f], uLocals[f], &uLocalArrys[f]);
251 CHKERRQ(ierr);
252
253 ierr = VecDestroy(&uLocals[f]); CHKERRQ(ierr);
254 }
255
256 PetscFunctionReturn(0);
257} // calculateVorticity2D
258
259PetscErrorCode calculateVorticity3D(const petibm::type::Mesh &mesh,
260 const petibm::type::Boundary &bds,
261 const petibm::type::Solution &soln,
262 const std::vector<DM> &wDm,
263 std::vector<Vec> &w)
264{
265 PetscErrorCode ierr;
266
267 PetscFunctionBeginUser;
268
269 std::vector<Vec> uLocals;
270 std::vector<PetscReal ***> uLocalArrys;
271 std::vector<DMDALocalInfo> info;
272 PetscReal ***wArray;
273
274 // resize the STL vectors
275 info.resize(mesh->dim);
276 uLocals.resize(mesh->dim);
277 uLocalArrys.resize(mesh->dim);
278
279 // create local Vecs and local infos
280 for (PetscInt f = 0; f < mesh->dim; ++f)
281 {
282 ierr = DMDAGetLocalInfo(wDm[f], &info[f]); CHKERRQ(ierr);
283 ierr = DMCreateLocalVector(mesh->da[f], &uLocals[f]); CHKERRQ(ierr);
284 }
285
286 // scatter from global Vecs to local Vecs
287 ierr = DMCompositeScatterArray(mesh->UPack, soln->UGlobal, uLocals.data());
288 CHKERRQ(ierr);
289
290 // copy the values at ghost points from Boundary instance to local Vecs
291 ierr = bds->copyValues2LocalVecs(uLocals); CHKERRQ(ierr);
292
293 // get read-only raw array from local Vecs
294 for (PetscInt f = 0; f < mesh->dim; ++f)
295 {
296 ierr = DMDAVecGetArrayRead(mesh->da[f], uLocals[f], &uLocalArrys[f]);
297 CHKERRQ(ierr);
298 }
299
300 // compute the x component of the vorticity field
301 ierr = DMDAVecGetArray(wDm[0], w[0], &wArray); CHKERRQ(ierr);
302 for (PetscInt k = info[0].zs; k < info[0].zs + info[0].zm; ++k)
303 {
304 for (PetscInt j = info[0].ys; j < info[0].ys + info[0].ym; ++j)
305 {
306 for (PetscInt i = info[0].xs; i < info[0].xs + info[0].xm; ++i)
307 {
308 wArray[k][j][i] =
309 (uLocalArrys[2][k - 1][j][i - 1] -
310 uLocalArrys[2][k - 1][j - 1][i - 1]) /
311 (mesh->coord[2][1][j] - mesh->coord[2][1][j - 1]) -
312 (uLocalArrys[1][k][j - 1][i - 1] -
313 uLocalArrys[1][k - 1][j - 1][i - 1]) /
314 (mesh->coord[1][2][k] - mesh->coord[1][2][k - 1]);
315 }
316 }
317 }
318 ierr = DMDAVecRestoreArray(wDm[0], w[0], &wArray); CHKERRQ(ierr);
319
320 // compute the y component of the vorticity field
321 ierr = DMDAVecGetArray(wDm[1], w[1], &wArray); CHKERRQ(ierr);
322 for (PetscInt k = info[1].zs; k < info[1].zs + info[1].zm; ++k)
323 {
324 for (PetscInt j = info[1].ys; j < info[1].ys + info[1].ym; ++j)
325 {
326 for (PetscInt i = info[1].xs; i < info[1].xs + info[1].xm; ++i)
327 {
328 wArray[k][j][i] =
329 (uLocalArrys[0][k][j - 1][i - 1] -
330 uLocalArrys[0][k - 1][j - 1][i - 1]) /
331 (mesh->coord[0][2][k] - mesh->coord[0][2][k - 1]) -
332 (uLocalArrys[2][k - 1][j - 1][i] -
333 uLocalArrys[2][k - 1][j - 1][i - 1]) /
334 (mesh->coord[2][0][i] - mesh->coord[2][0][i - 1]);
335 }
336 }
337 }
338 ierr = DMDAVecRestoreArray(wDm[1], w[1], &wArray); CHKERRQ(ierr);
339
340 // compute the z component of the vorticity field
341 ierr = DMDAVecGetArray(wDm[2], w[2], &wArray); CHKERRQ(ierr);
342 for (PetscInt k = info[2].zs; k < info[2].zs + info[2].zm; ++k)
343 {
344 for (PetscInt j = info[2].ys; j < info[2].ys + info[2].ym; ++j)
345 {
346 for (PetscInt i = info[2].xs; i < info[2].xs + info[2].xm; ++i)
347 {
348 wArray[k][j][i] =
349 (uLocalArrys[1][k][j - 1][i] -
350 uLocalArrys[1][k][j - 1][i - 1]) /
351 (mesh->coord[1][0][i] - mesh->coord[1][0][i - 1]) -
352 (uLocalArrys[0][k][j][i - 1] -
353 uLocalArrys[0][k][j - 1][i - 1]) /
354 (mesh->coord[0][1][j] - mesh->coord[0][1][j - 1]);
355 }
356 }
357 }
358 ierr = DMDAVecRestoreArray(wDm[2], w[2], &wArray); CHKERRQ(ierr);
359
360 // return the ownership of local velocity array and destroy local Vecs
361 for (PetscInt f = 0; f < mesh->dim; ++f)
362 {
363 ierr =
364 DMDAVecRestoreArrayRead(mesh->da[f], uLocals[f], &uLocalArrys[f]);
365 CHKERRQ(ierr);
366
367 ierr = VecDestroy(&uLocals[f]); CHKERRQ(ierr);
368 }
369
370 PetscFunctionReturn(0);
371} // calculateVorticity3D
372
373PetscErrorCode initVorticityMesh(const petibm::type::Mesh &mesh,
374 petibm::type::IntVec2D &n,
375 petibm::type::RealVec3D &coord,
376 std::vector<DM> &wDMs,
377 std::vector<std::string> &names,
378 std::vector<Vec> &vecs)
379{
380 PetscErrorCode ierr;
381
382 PetscFunctionBeginUser;
383
384 if (mesh->dim == 2)
385 {
386 // initialize size of meshes
387 n = petibm::type::IntVec2D(1, petibm::type::IntVec1D(3));
388
389 // set values of n for wz
390 n[0][0] = mesh->n[4][0];
391 n[0][1] = mesh->n[4][1];
392 n[0][2] = mesh->n[3][2]; // this value should be one
393
394 // initialize vectors holding coordinates
395 coord = petibm::type::RealVec3D(1, petibm::type::RealVec2D(3));
396
397 // set coordinates for wz
398 coord[0][0].assign(mesh->coord[4][0], mesh->coord[4][0] + n[0][0]);
399 coord[0][1].assign(mesh->coord[4][1], mesh->coord[4][1] + n[0][1]);
400 coord[0][2].assign(mesh->coord[3][2], mesh->coord[3][2] + n[0][2]);
401
402 // name of vorticity fields
403 names = std::vector<std::string>(1);
404
405 // set name for wz
406 names[0] = "wz";
407
408 // initialize vectors holding DMs and Vecs
409 wDMs = std::vector<DM>(1);
410 vecs = std::vector<Vec>(1);
411
412 // create DMs
413 ierr = DMDACreate2d(mesh->comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
414 DMDA_STENCIL_BOX, n[0][0], n[0][1], mesh->nProc[0],
415 mesh->nProc[1], 1, 1, nullptr, nullptr, &wDMs[0]);
416 CHKERRQ(ierr);
417
418 ierr = DMSetUp(wDMs[0]); CHKERRQ(ierr);
419
420 // create Vecs
421 ierr = DMCreateGlobalVector(wDMs[0], &vecs[0]); CHKERRQ(ierr);
422
423 // set names
424 ierr = PetscObjectSetName(PetscObject(vecs[0]), "wz"); CHKERRQ(ierr);
425 }
426 else
427 {
428 // initialize size of meshes
429 n = petibm::type::IntVec2D(3, petibm::type::IntVec1D(3));
430
431 // set values of n for wx
432 n[0][0] = mesh->n[3][0];
433 n[0][1] = mesh->n[4][1];
434 n[0][2] = mesh->n[4][2];
435 // set values of n for wy
436 n[1][0] = mesh->n[4][0];
437 n[1][1] = mesh->n[3][1];
438 n[1][2] = mesh->n[4][2];
439 // set values of n for wz
440 n[2][0] = mesh->n[4][0];
441 n[2][1] = mesh->n[4][1];
442 n[2][2] = mesh->n[3][2];
443
444 // initialize vectors holding coordinates
445 coord = petibm::type::RealVec3D(3, petibm::type::RealVec2D(3));
446
447 // set coordinates for wx
448 coord[0][0].assign(mesh->coord[3][0], mesh->coord[3][0] + n[0][0]);
449 coord[0][1].assign(mesh->coord[4][1], mesh->coord[4][1] + n[0][1]);
450 coord[0][2].assign(mesh->coord[4][2], mesh->coord[4][2] + n[0][2]);
451 // set coordinates for wy
452 coord[1][0].assign(mesh->coord[4][0], mesh->coord[4][0] + n[1][0]);
453 coord[1][1].assign(mesh->coord[3][1], mesh->coord[3][1] + n[1][1]);
454 coord[1][2].assign(mesh->coord[4][2], mesh->coord[4][2] + n[1][2]);
455 // set coordinates for wz
456 coord[2][0].assign(mesh->coord[4][0], mesh->coord[4][0] + n[2][0]);
457 coord[2][1].assign(mesh->coord[4][1], mesh->coord[4][1] + n[2][1]);
458 coord[2][2].assign(mesh->coord[3][2], mesh->coord[3][2] + n[2][2]);
459
460 // name of vorticity fields
461 names = std::vector<std::string>(3);
462
463 // set name for wx
464 names[0] = "wx";
465 // set name for wy
466 names[1] = "wy";
467 // set name for wz
468 names[2] = "wz";
469
470 wDMs = std::vector<DM>(3);
471 vecs = std::vector<Vec>(3);
472
473 for (PetscInt f = 0; f < 3; ++f)
474 {
475 // create DMs
476 ierr = DMDACreate3d(mesh->comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
477 DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, n[f][0],
478 n[f][1], n[f][2], mesh->nProc[0],
479 mesh->nProc[1], mesh->nProc[2], 1, 1, nullptr,
480 nullptr, nullptr, &wDMs[f]); CHKERRQ(ierr);
481
482 ierr = DMSetUp(wDMs[f]); CHKERRQ(ierr);
483
484 // create Vecs
485 ierr = DMCreateGlobalVector(wDMs[f], &vecs[f]); CHKERRQ(ierr);
486
487 // set name
488 ierr = PetscObjectSetName(PetscObject(vecs[f]), names[f].c_str());
489 CHKERRQ(ierr);
490 }
491 }
492
493 PetscFunctionReturn(0);
494} // initVorticityMesh
main
int main(int argc, char **argv)
Definition: main.cpp:48
calculateVorticity2D
PetscErrorCode calculateVorticity2D(const petibm::type::Mesh &mesh, const petibm::type::Boundary &bds, const petibm::type::Solution &soln, const std::vector< DM > &wDm, std::vector< Vec > &w)
Definition: main.cpp:185
calculateVorticity3D
PetscErrorCode calculateVorticity3D(const petibm::type::Mesh &mesh, const petibm::type::Boundary &bds, const petibm::type::Solution &soln, const std::vector< DM > &wDm, std::vector< Vec > &w)
Definition: main.cpp:259
initVorticityMesh
PetscErrorCode initVorticityMesh(const petibm::type::Mesh &mesh, petibm::type::IntVec2D &n, petibm::type::RealVec3D &coord, std::vector< DM > &wDMs, std::vector< std::string > &names, std::vector< Vec > &vecs)
Definition: main.cpp:373
boundary.h
boundary::BoundaryBase, type::Boundary, and factory function.
petibm::boundary::createBoundary
PetscErrorCode createBoundary(const type::Mesh &mesh, const YAML::Node &node, type::Boundary &boundary)
Create a Boundary object.
Definition: boundary.cpp:41
petibm::type::Boundary
std::shared_ptr< boundary::BoundaryBase > Boundary
Type definition of petibm::type::Boundary.
Definition: boundary.h:175
petibm::type::Mesh
std::shared_ptr< mesh::MeshBase > Mesh
Type definition of Mesh.
Definition: mesh.h:348
petibm::mesh::createMesh
PetscErrorCode createMesh(const MPI_Comm &comm, const YAML::Node &node, type::Mesh &mesh)
Factory function for creating a Mesh object.
Definition: mesh.cpp:86
petibm::parser::getSettings
PetscErrorCode getSettings(YAML::Node &node)
Get configuration settings.
Definition: parser.cpp:175
petibm::io::writeHDF5Vecs
PetscErrorCode writeHDF5Vecs(const MPI_Comm comm, const std::string &filePath, const std::string &loc, const std::vector< std::string > &names, const std::vector< Vec > &vecs, const PetscFileMode mode=FILE_MODE_WRITE)
Write a vector of Vec objects to a HDF5 file.
Definition: io.cpp:137
petibm::type::Solution
std::shared_ptr< solution::SolutionBase > Solution
Type definition of solution object.
Definition: solution.h:210
petibm::solution::createSolution
PetscErrorCode createSolution(const type::Mesh &mesh, type::Solution &solution)
Factory function to create a petibm::solution::Solution object.
Definition: solution.cpp:65
petibm::type::IntVec2D
std::vector< IntVec1D > IntVec2D
2D std::vector holding PetscInt.
Definition: type.h:135
petibm::type::RealVec2D
std::vector< RealVec1D > RealVec2D
2D std::vector holding PetscReal.
Definition: type.h:142
petibm::type::IntVec1D
std::vector< PetscInt > IntVec1D
1D std::vector holding PetscInt.
Definition: type.h:133
petibm::type::RealVec3D
std::vector< RealVec2D > RealVec3D
3D std::vector holding PetscReal.
Definition: type.h:144
petibm::type::w
@ w
Definition: type.h:83
io.h
Prototypes of I/O functions.
mesh.h
Prototype of mesh::MeshBase, type::Mesh, and factory function.
parser.h
Prototypes of parser functions.
solution.h
Definition of the class petibm::solution::SolutionBase, the type definition petibm::type::Solution,...
type.h
Definition of user-defined types for convenience.