PetIBM 0.5.4
Toolbox and applications of the immersed-boundary method for distributed-memory architectures
cartesianmesh2d_dirichlet.cpp
Go to the documentation of this file.
1
8#include <petsc.h>
9
10#include <gtest/gtest.h>
11#include <yaml-cpp/yaml.h>
12
13#include <petibm/mesh.h>
14
15class CartesianMeshTest2D_AllDirichlet : public ::testing::Test
16{
17protected:
18 CartesianMeshTest2D_AllDirichlet(){};
19
20 virtual ~CartesianMeshTest2D_AllDirichlet(){};
21
22 static void SetUpTestCase()
23 {
24 using namespace YAML;
25
26 Node config;
27
28 config["mesh"].push_back(Node(NodeType::Map));
29
30 config["mesh"][0]["direction"] = "x";
31 config["mesh"][0]["start"] = "0.1";
32 config["mesh"][0]["subDomains"].push_back(Node(NodeType::Map));
33 config["mesh"][0]["subDomains"][0]["end"] = 1.6;
34 config["mesh"][0]["subDomains"][0]["cells"] = 4;
35 config["mesh"][0]["subDomains"][0]["stretchRatio"] = 0.5;
36 config["mesh"][0]["subDomains"][1]["end"] = 1.9;
37 config["mesh"][0]["subDomains"][1]["cells"] = 3;
38 config["mesh"][0]["subDomains"][1]["stretchRatio"] = 1;
39 config["mesh"][0]["subDomains"][2]["end"] = 5.0;
40 config["mesh"][0]["subDomains"][2]["cells"] = 5;
41 config["mesh"][0]["subDomains"][2]["stretchRatio"] = 2.0;
42
43 config["mesh"][1]["direction"] = "y";
44 config["mesh"][1]["start"] = "0.05";
45 config["mesh"][1]["subDomains"].push_back(Node(NodeType::Map));
46 config["mesh"][1]["subDomains"][0]["end"] = 0.8625;
47 config["mesh"][1]["subDomains"][0]["cells"] = 4;
48 config["mesh"][1]["subDomains"][0]["stretchRatio"] = 0.666666666666666;
49 config["mesh"][1]["subDomains"][1]["end"] = 1.1625;
50 config["mesh"][1]["subDomains"][1]["cells"] = 3;
51 config["mesh"][1]["subDomains"][1]["stretchRatio"] = 1.0;
52 config["mesh"][1]["subDomains"][2]["end"] = 1.975;
53 config["mesh"][1]["subDomains"][2]["cells"] = 4;
54 config["mesh"][1]["subDomains"][2]["stretchRatio"] = 1.5;
55
56 config["flow"] = YAML::Node(NodeType::Map);
57 config["flow"]["boundaryConditions"].push_back(Node(NodeType::Map));
58 config["flow"]["boundaryConditions"][0]["location"] = "xMinus";
59 config["flow"]["boundaryConditions"][1]["location"] = "xPlus";
60 config["flow"]["boundaryConditions"][2]["location"] = "yMinus";
61 config["flow"]["boundaryConditions"][3]["location"] = "yPlus";
62
63 for (unsigned int i = 0; i < 4; ++i)
64 {
65 config["flow"]["boundaryConditions"][i]["u"][0] = "DIRICHLET";
66 config["flow"]["boundaryConditions"][i]["u"][1] = 0.0;
67 config["flow"]["boundaryConditions"][i]["v"][0] = "DIRICHLET";
68 config["flow"]["boundaryConditions"][i]["v"][1] = 0.0;
69 }
70 config["flow"]["boundaryConditions"][3]["u"][1] = 1.0;
71
72 PetscErrorCode ierr;
73 PetscInt nProc[3];
74
75 ierr = petibm::mesh::createMesh(PETSC_COMM_WORLD, config, mesh);
76 ASSERT_FALSE(ierr);
77
78 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED,
79 DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, 12, 11,
80 PETSC_DECIDE, PETSC_DECIDE, 1, 1, nullptr, nullptr,
81 &da[3]);
82 ASSERT_FALSE(ierr);
83 ierr = DMSetUp(da[3]);
84 ASSERT_FALSE(ierr);
85
86 ierr = DMDAGetInfo(da[3], nullptr, nullptr, nullptr, nullptr, &nProc[0],
87 &nProc[1], &nProc[2], nullptr, nullptr, nullptr,
88 nullptr, nullptr, nullptr);
89 ASSERT_FALSE(ierr);
90
91 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED,
92 DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, 11, 11,
93 nProc[0], nProc[1], 1, 1, nullptr, nullptr, &da[0]);
94 ASSERT_FALSE(ierr);
95 ierr = DMSetUp(da[0]);
96 ASSERT_FALSE(ierr);
97
98 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED,
99 DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, 12, 10,
100 nProc[0], nProc[1], 1, 1, nullptr, nullptr, &da[1]);
101 ASSERT_FALSE(ierr);
102 ierr = DMSetUp(da[1]);
103 ASSERT_FALSE(ierr);
104 };
105
106 virtual void SetUp(){};
107
108 virtual void TearDown(){};
109
110 static void TearDownTestCase()
111 {
112 PetscErrorCode ierr;
113 ierr = DMDestroy(&da[0]);
114 ASSERT_FALSE(ierr);
115 ierr = DMDestroy(&da[1]);
116 ASSERT_FALSE(ierr);
117 ierr = DMDestroy(&da[3]);
118 ASSERT_FALSE(ierr);
119 };
120
121 static petibm::type::Mesh mesh;
122 static DM da[5];
123}; // CartesianMeshTest
124
125petibm::type::Mesh CartesianMeshTest2D_AllDirichlet::mesh = nullptr;
126DM CartesianMeshTest2D_AllDirichlet::da[5]{nullptr, nullptr, nullptr, nullptr,
127 nullptr};
128
129// test dimension
130TEST_F(CartesianMeshTest2D_AllDirichlet, check_dim) { ASSERT_EQ(2, mesh->dim); }
131
132// check min
133TEST_F(CartesianMeshTest2D_AllDirichlet, check_min)
134{
135 ASSERT_DOUBLE_EQ(0.1, mesh->min[0]);
136 ASSERT_DOUBLE_EQ(0.05, mesh->min[1]);
137 ASSERT_DOUBLE_EQ(0.0, mesh->min[2]);
138}
139
140// check max
141TEST_F(CartesianMeshTest2D_AllDirichlet, check_max)
142{
143 ASSERT_DOUBLE_EQ(5.0, mesh->max[0]);
144 ASSERT_DOUBLE_EQ(1.975, mesh->max[1]);
145 ASSERT_DOUBLE_EQ(1.0, mesh->max[2]);
146}
147
148// check n
149TEST_F(CartesianMeshTest2D_AllDirichlet, check_n)
150{
151 petibm::type::RealVec2D n{
152 {11, 11, 1}, {12, 10, 1}, {1, 1, 1}, {12, 11, 1}, {13, 12, 1}};
153
154 for (unsigned int f = 0; f < 5; ++f)
155 for (unsigned int d = 0; d < 3; ++d) ASSERT_EQ(n[f][d], mesh->n[f][d]);
156}
157
158// check periodic
159TEST_F(CartesianMeshTest2D_AllDirichlet, check_periodic)
160{
161 petibm::type::BoolVec2D periodic{{PETSC_FALSE, PETSC_FALSE, PETSC_FALSE},
162 {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE},
163 {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE}};
164
165 for (unsigned int f = 0; f < 3; ++f)
166 for (unsigned int d = 0; d < 3; ++d)
167 ASSERT_EQ(periodic[f][d], mesh->periodic[f][d]);
168}
169
170// check coord
171TEST_F(CartesianMeshTest2D_AllDirichlet, check_coord)
172{
173 petibm::type::RealVec3D coordTrue{
174 {{0.1, 0.9, 1.3, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.6, 3.4, 5.0},
175 {-0.11875, 0.21875, 0.5, 0.6875, 0.8125, 0.9125, 1.0125, 1.1125,
176 1.2125, 1.3375, 1.525, 1.80625, 2.14375},
177 {0.0}},
178 {{-0.3, 0.5, 1.1, 1.4, 1.55, 1.65, 1.75, 1.85, 1.95, 2.1, 2.4, 3.0, 4.2,
179 5.8},
180 {0.05, 0.3875, 0.6125, 0.7625, 0.8625, 0.9625, 1.0625, 1.1625, 1.2625,
181 1.4125, 1.6375, 1.975},
182 {0.0}},
183 {{0.0}, {0.0}, {0.0}},
184 {{0.5, 1.1, 1.4, 1.55, 1.65, 1.75, 1.85, 1.95, 2.1, 2.4, 3.0, 4.2},
185 {0.21875, 0.5, 0.6875, 0.8125, 0.9125, 1.0125, 1.1125, 1.2125, 1.3375,
186 1.525, 1.80625},
187 {0.0}},
188 {{0.1, 0.9, 1.3, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.6, 3.4, 5.0},
189 {0.05, 0.3875, 0.6125, 0.7625, 0.8625, 0.9625, 1.0625, 1.1625, 1.2625,
190 1.4125, 1.6375, 1.975},
191 {0.0}}};
192
193 petibm::type::GhostedVec3D coord{
194 {&coordTrue[0][0][1], &coordTrue[0][1][1], &coordTrue[0][2][0]},
195 {&coordTrue[1][0][1], &coordTrue[1][1][1], &coordTrue[1][2][0]},
196 {&coordTrue[2][0][0], &coordTrue[2][1][0], &coordTrue[2][2][0]},
197 {&coordTrue[3][0][0], &coordTrue[3][1][0], &coordTrue[3][2][0]},
198 {&coordTrue[4][0][0], &coordTrue[4][1][0], &coordTrue[4][2][0]},
199 };
200
201 for (unsigned int f = 0; f < 2; ++f)
202 {
203 for (unsigned int d = 0; d < 2; ++d)
204 for (int i = -1; i < mesh->n[f][d] + 1; ++i)
205 ASSERT_NEAR(coord[f][d][i], mesh->coord[f][d][i], 1e-12)
206 << "dL incorrect at field " << f << ", direction " << d
207 << ", index " << i << std::endl;
208
209 ASSERT_NEAR(coord[f][2][0], mesh->coord[f][2][0], 1e-12)
210 << "dL incorrect at field " << f << ", direction " << 2
211 << ", index " << 0 << std::endl;
212 }
213
214 for (unsigned int f = 2; f < 5; ++f)
215 for (unsigned int d = 0; d < 3; ++d)
216 for (int i = 0; i < mesh->n[f][d]; ++i)
217 ASSERT_NEAR(coord[f][d][i], mesh->coord[f][d][i], 1e-12)
218 << "dL incorrect at field " << f << ", direction " << d
219 << ", index " << i << std::endl;
220}
221
222// check dL
223TEST_F(CartesianMeshTest2D_AllDirichlet, check_dL)
224{
225 petibm::type::RealVec3D dLTrue{
226 {{0.8, 0.6, 0.3, 0.15, 0.1, 0.1, 0.1, 0.1, 0.15, 0.3, 0.6, 1.2, 1.6},
227 {0.3375, 0.3375, 0.225, 0.15, 0.1, 0.1, 0.1, 0.1, 0.1, 0.15, 0.225,
228 0.3375, 0.3375},
229 {1.0}},
230 {{0.8, 0.8, 0.4, 0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.4, 0.8, 1.6, 1.6},
231 {0.3375, 0.28125, 0.1875, 0.125, 0.1, 0.1, 0.1, 0.1, 0.125, 0.1875,
232 0.28125, 0.3375},
233 {1.0}},
234 {{1.0}, {1.0}, {1.0}},
235 {{0.8, 0.4, 0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.4, 0.8, 1.6},
236 {0.3375, 0.225, 0.15, 0.1, 0.1, 0.1, 0.1, 0.1, 0.15, 0.225, 0.3375},
237 {1.0}},
238 {{1.0}, {1.0}, {1.0}}};
239
240 petibm::type::GhostedVec3D dL{
241 {&dLTrue[0][0][1], &dLTrue[0][1][1], &dLTrue[0][2][0]},
242 {&dLTrue[1][0][1], &dLTrue[1][1][1], &dLTrue[1][2][0]},
243 {&dLTrue[2][0][0], &dLTrue[2][1][0], &dLTrue[2][2][0]},
244 {&dLTrue[3][0][0], &dLTrue[3][1][0], &dLTrue[3][2][0]},
245 {nullptr, nullptr, nullptr}};
246
247 for (unsigned int f = 0; f < 2; ++f)
248 {
249 for (unsigned int d = 0; d < 2; ++d)
250 for (int i = -1; i < mesh->n[f][d] + 1; ++i)
251 ASSERT_NEAR(dL[f][d][i], mesh->dL[f][d][i], 1e-12)
252 << "dL incorrect at field " << f << ", direction " << d
253 << ", index " << i << std::endl;
254
255 ASSERT_NEAR(dL[f][2][0], mesh->dL[f][2][0], 1e-12)
256 << "dL incorrect at field " << f << ", direction " << 2
257 << ", index " << 0 << std::endl;
258 }
259
260 for (unsigned int f = 2; f < 4; ++f)
261 for (unsigned int d = 0; d < 3; ++d)
262 for (int i = 0; i < mesh->n[f][d]; ++i)
263 ASSERT_NEAR(dL[f][d][i], mesh->dL[f][d][i], 1e-12)
264 << "dL incorrect at field " << f << ", direction " << d
265 << ", index " << i << std::endl;
266
267 for (unsigned int d = 0; d < 3; ++d)
268 ASSERT_FALSE(mesh->dL[4][d])
269 << "The dL of the vertex mesh is not nullptr." << std::endl;
270}
271
272// check UN
273TEST_F(CartesianMeshTest2D_AllDirichlet, check_UN)
274{
275 ASSERT_EQ(241, mesh->UN)
276 << "The total number of velocity points is not correct." << std::endl;
277}
278
279// check pN
280TEST_F(CartesianMeshTest2D_AllDirichlet, check_pN)
281{
282 ASSERT_EQ(132, mesh->pN)
283 << "The total number of pressure points is not correct." << std::endl;
284}
285
286// check da
287TEST_F(CartesianMeshTest2D_AllDirichlet, check_da)
288{
289 PetscErrorCode ierr;
290 DMType type;
291
292 ierr = DMGetType(mesh->da[0], &type);
293 ASSERT_FALSE(ierr);
294 ASSERT_STREQ(type, DMDA)
295 << "u-velocity mesh is " << type << ", not expected da" << std::endl;
296
297 ierr = DMGetType(mesh->da[1], &type);
298 ASSERT_FALSE(ierr);
299 ASSERT_STREQ(type, DMDA)
300 << "v-velocity mesh is " << type << ", not expected da" << std::endl;
301
302 ierr = DMGetType(mesh->da[3], &type);
303 ASSERT_FALSE(ierr);
304 ASSERT_STREQ(type, DMDA)
305 << "v-velocity mesh is " << type << ", not expected da" << std::endl;
306
307 ASSERT_FALSE(mesh->da[2])
308 << "w-velocity mesh should not be initialized." << std::endl;
309 ASSERT_FALSE(mesh->da[4])
310 << "vertex mesh should not be initialized." << std::endl;
311}
312
313// check nProc
314TEST_F(CartesianMeshTest2D_AllDirichlet, check_nProc)
315{
316 PetscErrorCode ierr;
317
318 PetscInt nProc[3];
319
320 ierr = DMDAGetInfo(da[3], nullptr, nullptr, nullptr, nullptr, &nProc[0],
321 &nProc[1], &nProc[2], nullptr, nullptr, nullptr, nullptr,
322 nullptr, nullptr);
323 ASSERT_FALSE(ierr);
324
325 ASSERT_EQ(nProc[0], mesh->nProc[0]);
326 ASSERT_EQ(nProc[1], mesh->nProc[1]);
327 ASSERT_EQ(nProc[2], mesh->nProc[2]);
328}
329
330// check bg
331TEST_F(CartesianMeshTest2D_AllDirichlet, check_bg)
332{
333 PetscErrorCode ierr;
334 for (unsigned int f = 0; f < 5; f++)
335 {
336 if ((f == 2) || (f == 4)) continue;
337 DMDALocalInfo info;
338 ierr = DMDAGetLocalInfo(da[f], &info);
339 ASSERT_FALSE(ierr);
340 ASSERT_EQ(info.xs, mesh->bg[f][0]) << "f = " << f << std::endl;
341 ASSERT_EQ(info.ys, mesh->bg[f][1]) << "f = " << f << std::endl;
342 ASSERT_EQ(info.zs, mesh->bg[f][2]) << "f = " << f << std::endl;
343 }
344
345 ASSERT_EQ(0, mesh->bg[2][0]) << "f = " << 2 << std::endl;
346 ASSERT_EQ(0, mesh->bg[2][1]) << "f = " << 2 << std::endl;
347 ASSERT_EQ(0, mesh->bg[2][2]) << "f = " << 2 << std::endl;
348}
349
350// check ed
351TEST_F(CartesianMeshTest2D_AllDirichlet, check_ed)
352{
353 PetscErrorCode ierr;
354 for (unsigned int f = 0; f < 5; f++)
355 {
356 if ((f == 2) || (f == 4)) continue;
357 DMDALocalInfo info;
358 ierr = DMDAGetLocalInfo(da[f], &info);
359 ASSERT_FALSE(ierr);
360 ASSERT_EQ(info.xs + info.xm, mesh->ed[f][0])
361 << "f = " << f << std::endl;
362 ASSERT_EQ(info.ys + info.ym, mesh->ed[f][1])
363 << "f = " << f << std::endl;
364 ASSERT_EQ(info.zs + info.zm, mesh->ed[f][2])
365 << "f = " << f << std::endl;
366 }
367
368 ASSERT_EQ(1, mesh->ed[2][0]) << "f = " << 2 << std::endl;
369 ASSERT_EQ(1, mesh->ed[2][1]) << "f = " << 2 << std::endl;
370 ASSERT_EQ(1, mesh->ed[2][2]) << "f = " << 2 << std::endl;
371}
372
373// check m
374TEST_F(CartesianMeshTest2D_AllDirichlet, check_m)
375{
376 PetscErrorCode ierr;
377 for (unsigned int f = 0; f < 5; f++)
378 {
379 if ((f == 2) || (f == 4)) continue;
380 DMDALocalInfo info;
381 ierr = DMDAGetLocalInfo(da[f], &info);
382 ASSERT_FALSE(ierr);
383 ASSERT_EQ(info.xm, mesh->m[f][0]) << "f = " << f << std::endl;
384 ASSERT_EQ(info.ym, mesh->m[f][1]) << "f = " << f << std::endl;
385 ASSERT_EQ(info.zm, mesh->m[f][2]) << "f = " << f << std::endl;
386 }
387
388 ASSERT_EQ(0, mesh->m[2][0]) << "f = " << 2 << std::endl;
389 ASSERT_EQ(0, mesh->m[2][1]) << "f = " << 2 << std::endl;
390 ASSERT_EQ(0, mesh->m[2][2]) << "f = " << 2 << std::endl;
391}
392
393// check UNLocal
394TEST_F(CartesianMeshTest2D_AllDirichlet, check_UNLocal)
395{
396 int N = mesh->m[0][0] * mesh->m[0][1] + mesh->m[1][0] * mesh->m[1][1];
397 ASSERT_EQ(N, mesh->UNLocal);
398}
399
400// check pNLocal
401TEST_F(CartesianMeshTest2D_AllDirichlet, check_pNLocal)
402{
403 int N = mesh->m[3][0] * mesh->m[3][1];
404 ASSERT_EQ(N, mesh->pNLocal);
405}
406
407// check UPack
408TEST_F(CartesianMeshTest2D_AllDirichlet, check_UPack)
409{
410 PetscErrorCode ierr;
411 DMType type;
412
413 ierr = DMGetType(mesh->UPack, &type);
414 ASSERT_FALSE(ierr);
415 ASSERT_STREQ(type, DMCOMPOSITE);
416
417 DM dms[3]{nullptr, nullptr, nullptr};
418 ierr = DMCompositeGetEntriesArray(mesh->UPack, dms);
419 ASSERT_FALSE(ierr);
420 ASSERT_EQ(dms[0], mesh->da[0]);
421 ASSERT_EQ(dms[1], mesh->da[1]);
422 ASSERT_EQ(dms[2], mesh->da[2]);
423
424 ierr = DMDestroy(&dms[0]);
425 ASSERT_FALSE(ierr);
426 ierr = DMDestroy(&dms[1]);
427 ASSERT_FALSE(ierr);
428}
429
430// check mpi
431TEST_F(CartesianMeshTest2D_AllDirichlet, check_mpi)
432{
433 PetscErrorCode ierr;
434 PetscMPIInt rank, size;
435
436 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);
437 ASSERT_FALSE(ierr);
438 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
439 ASSERT_FALSE(ierr);
440
441 ASSERT_EQ(PETSC_COMM_WORLD, mesh->comm);
442 ASSERT_EQ(size, mesh->mpiSize);
443 ASSERT_EQ(rank, mesh->mpiRank);
444}
445
446// check getLocalIndex with i, j, k
447TEST_F(CartesianMeshTest2D_AllDirichlet, check_getLocalIndex_ijk)
448{
449 for (int f = 0; f < 5; ++f)
450 {
451 switch (f)
452 {
453 case 2:
454 case 4:
455 break;
456 default:
457 for (int k = mesh->bg[f][2]; k < mesh->ed[f][2]; k++)
458 for (int j = mesh->bg[f][1] - 1; j < mesh->ed[f][1] + 1;
459 j++)
460 for (int i = mesh->bg[f][0] - 1; i < mesh->ed[f][0] + 1;
461 i++)
462 {
463 PetscInt idx;
464 mesh->getLocalIndex(f, i, j, k, idx);
465
466 PetscInt i2, j2, k2, idx2;
467
468 i2 = i - mesh->bg[f][0] + 1;
469 j2 = j - mesh->bg[f][1] + 1;
470 k2 = k - mesh->bg[f][2];
471
472 idx2 = i2 + j2 * (mesh->m[f][0] + 2) +
473 k2 * ((mesh->m[f][0] + 2) *
474 (mesh->m[f][1] + 2));
475
476 ASSERT_EQ(idx2, idx)
477 << "getLocalIndex wrong at field " << f
478 << " and (" << i << ", " << j << ", " << k
479 << ") " << std::endl;
480 }
481 }
482 }
483}
484
485// check getNaturalIndex with i, j, k
486TEST_F(CartesianMeshTest2D_AllDirichlet, check_getNaturalIndex_ijk)
487{
488 for (int f = 0; f < 5; ++f)
489 {
490 // internal points
491 for (int k = 0; k < mesh->n[f][2]; k++)
492 for (int j = 0; j < mesh->n[f][1]; j++)
493 for (int i = 0; i < mesh->n[f][0]; i++)
494 {
495 PetscInt idx, idx2;
496 mesh->getNaturalIndex(f, i, j, k, idx);
497
498 idx2 = i + j * mesh->n[f][0] +
499 k * (mesh->n[f][0] * mesh->n[f][1]);
500
501 ASSERT_EQ(idx2, idx)
502 << "getNatrualIndex wrong at field " << f << " and ("
503 << i << ", " << j << ", " << k << ") " << std::endl;
504 }
505
506 // ghost points at top and bottom
507 for (int k = 0; k < mesh->n[f][2]; k++)
508 for (int i = 0; i < mesh->n[f][0]; i++)
509 {
510 PetscInt idx, j;
511
512 j = mesh->n[f][1];
513 mesh->getNaturalIndex(f, i, j, k, idx);
514 ASSERT_EQ(-1, idx)
515 << "getNatrualIndex wrong at field " << f << " and (" << i
516 << ", " << j << ", " << k << ") " << std::endl;
517
518 j = -1;
519 mesh->getNaturalIndex(f, i, j, k, idx);
520 ASSERT_EQ(-1, idx)
521 << "getNatrualIndex wrong at field " << f << " and (" << i
522 << ", " << j << ", " << k << ") " << std::endl;
523 }
524
525 // ghost points at right and left
526 for (int k = 0; k < mesh->n[f][2]; k++)
527 for (int j = 0; j < mesh->n[f][1]; j++)
528 {
529 PetscInt idx, i;
530
531 i = mesh->n[f][0];
532 mesh->getNaturalIndex(f, i, j, k, idx);
533 ASSERT_EQ(-1, idx)
534 << "getNatrualIndex wrong at field " << f << " and (" << i
535 << ", " << j << ", " << k << ") " << std::endl;
536
537 i = -1;
538 mesh->getNaturalIndex(f, i, j, k, idx);
539 ASSERT_EQ(-1, idx)
540 << "getNatrualIndex wrong at field " << f << " and (" << i
541 << ", " << j << ", " << k << ") " << std::endl;
542 }
543
544 // corners
545 for (int k = 0; k < mesh->n[f][2]; k++)
546 {
547 PetscInt idx;
548
549 mesh->getNaturalIndex(f, -1, -1, 0, idx);
550 ASSERT_EQ(-1, idx)
551 << "getNatrualIndex wrong at field " << f << " and (" << -1
552 << ", " << -1 << ", " << k << ") " << std::endl;
553
554 mesh->getNaturalIndex(f, -1, mesh->n[f][1], 0, idx);
555 ASSERT_EQ(-1, idx)
556 << "getNatrualIndex wrong at field " << f << " and (" << -1
557 << ", " << mesh->n[f][1] << ", " << k << ") " << std::endl;
558
559 mesh->getNaturalIndex(f, mesh->n[f][0], -1, 0, idx);
560 ASSERT_EQ(-1, idx) << "getNatrualIndex wrong at field " << f
561 << " and (" << mesh->n[f][0] << ", " << -1
562 << ", " << k << ") " << std::endl;
563
564 mesh->getNaturalIndex(f, mesh->n[f][0], mesh->n[f][1], 0, idx);
565 ASSERT_EQ(-1, idx)
566 << "getNatrualIndex wrong at field " << f << " and ("
567 << mesh->n[f][0] << ", " << mesh->n[f][1] << ", " << k << ") "
568 << std::endl;
569 }
570 }
571}
572
573// check getGlobalIndex with i, j, k of internal points
574TEST_F(CartesianMeshTest2D_AllDirichlet, check_getGlobalIndex_ijk_internal)
575{
576 for (int f = 0; f < 5; ++f)
577 {
578 if ((f == 2) || (f == 4)) continue;
579
580 PetscErrorCode ierr;
581
582 ISLocalToGlobalMapping mapping;
583 ierr = DMGetLocalToGlobalMapping(mesh->da[f], &mapping);
584 ASSERT_FALSE(ierr);
585
586 Vec vini, vinj, vink;
587 ierr = DMGetGlobalVector(mesh->da[f], &vini);
588 ASSERT_FALSE(ierr);
589 ierr = DMGetGlobalVector(mesh->da[f], &vinj);
590 ASSERT_FALSE(ierr);
591 ierr = DMGetGlobalVector(mesh->da[f], &vink);
592 ASSERT_FALSE(ierr);
593
594 for (int k = mesh->bg[f][2]; k < mesh->ed[f][2]; k++)
595 for (int j = mesh->bg[f][1]; j < mesh->ed[f][1]; j++)
596 for (int i = mesh->bg[f][0]; i < mesh->ed[f][0]; i++)
597 {
598 PetscInt idx;
599
600 ierr = DMDAConvertToCell(mesh->da[f], {k, j, i, 0}, &idx);
601 ASSERT_FALSE(ierr);
602 ierr = ISLocalToGlobalMappingApply(mapping, 1, &idx, &idx);
603 ASSERT_FALSE(ierr);
604
605 ierr = VecSetValue(vini, idx, i, INSERT_VALUES);
606 ASSERT_FALSE(ierr);
607 ierr = VecSetValue(vinj, idx, j, INSERT_VALUES);
608 ASSERT_FALSE(ierr);
609 ierr = VecSetValue(vink, idx, k, INSERT_VALUES);
610 ASSERT_FALSE(ierr);
611 }
612 ierr = VecAssemblyBegin(vini);
613 ASSERT_FALSE(ierr);
614 ierr = VecAssemblyBegin(vinj);
615 ASSERT_FALSE(ierr);
616 ierr = VecAssemblyBegin(vink);
617 ASSERT_FALSE(ierr);
618 ierr = VecAssemblyEnd(vini);
619 ASSERT_FALSE(ierr);
620 ierr = VecAssemblyEnd(vinj);
621 ASSERT_FALSE(ierr);
622 ierr = VecAssemblyEnd(vink);
623 ASSERT_FALSE(ierr);
624
625 Vec vouti, voutj, voutk;
626 VecScatter ctx;
627 ierr = VecScatterCreateToAll(vini, &ctx, &vouti);
628 ASSERT_FALSE(ierr);
629 ierr = VecDuplicate(vouti, &voutj);
630 ASSERT_FALSE(ierr);
631 ierr = VecDuplicate(vouti, &voutk);
632 ASSERT_FALSE(ierr);
633 ierr =
634 VecScatterBegin(ctx, vini, vouti, INSERT_VALUES, SCATTER_FORWARD);
635 ASSERT_FALSE(ierr);
636 ierr = VecScatterEnd(ctx, vini, vouti, INSERT_VALUES, SCATTER_FORWARD);
637 ASSERT_FALSE(ierr);
638 ierr =
639 VecScatterBegin(ctx, vinj, voutj, INSERT_VALUES, SCATTER_FORWARD);
640 ASSERT_FALSE(ierr);
641 ierr = VecScatterEnd(ctx, vinj, voutj, INSERT_VALUES, SCATTER_FORWARD);
642 ASSERT_FALSE(ierr);
643 ierr =
644 VecScatterBegin(ctx, vink, voutk, INSERT_VALUES, SCATTER_FORWARD);
645 ASSERT_FALSE(ierr);
646 ierr = VecScatterEnd(ctx, vink, voutk, INSERT_VALUES, SCATTER_FORWARD);
647 ASSERT_FALSE(ierr);
648
649 // internal points
650 for (int k = 0; k < mesh->n[f][2]; k++)
651 for (int j = 0; j < mesh->n[f][1]; j++)
652 for (int i = 0; i < mesh->n[f][0]; i++)
653 {
654 PetscInt idx;
655 ierr = mesh->getGlobalIndex(f, i, j, k, idx);
656 ASSERT_FALSE(ierr);
657
658 PetscReal vi, vj, vk;
659 ierr = VecGetValues(vouti, 1, &idx, &vi);
660 ASSERT_FALSE(ierr);
661 ierr = VecGetValues(voutj, 1, &idx, &vj);
662 ASSERT_FALSE(ierr);
663 ierr = VecGetValues(voutk, 1, &idx, &vk);
664 ASSERT_FALSE(ierr);
665
666 ASSERT_EQ(vi, PetscReal(i));
667 ASSERT_EQ(vj, PetscReal(j));
668 ASSERT_EQ(vk, PetscReal(k));
669 }
670
671 ierr = VecScatterDestroy(&ctx);
672 ASSERT_FALSE(ierr);
673 ierr = VecDestroy(&voutk);
674 ASSERT_FALSE(ierr);
675 ierr = VecDestroy(&voutj);
676 ASSERT_FALSE(ierr);
677 ierr = VecDestroy(&vouti);
678 ASSERT_FALSE(ierr);
679 ierr = VecDestroy(&vink);
680 ASSERT_FALSE(ierr);
681 ierr = VecDestroy(&vinj);
682 ASSERT_FALSE(ierr);
683 ierr = VecDestroy(&vini);
684 ASSERT_FALSE(ierr);
685
686 // ierr = ISLocalToGlobalMappingDestroy(&mapping); ASSERT_FALSE(ierr);
687 // Seem there's a bug in ISLocalToGlobalMappingDestroy. It should just
688 // decrease the reference number by one. Instead, it now destroys the
689 // underlying object completely....
690 mapping = PETSC_NULL;
691 }
692}
693
694// check getGlobalIndex with i, j, k of ghost points
695TEST_F(CartesianMeshTest2D_AllDirichlet, check_getGlobalIndex_ijk_ghosts)
696{
697 for (int f = 0; f < 5; ++f)
698 {
699 if ((f == 2) || (f == 4)) continue;
700
701 // ghost points at top and bottom
702 for (int k = 0; k < mesh->n[f][2]; k++)
703 for (int i = 0; i < mesh->n[f][0]; i++)
704 {
705 PetscInt idx, j;
706
707 j = mesh->n[f][1];
708 mesh->getGlobalIndex(f, i, j, k, idx);
709 ASSERT_EQ(-1, idx)
710 << "getGlobalIndex wrong at field " << f << " and (" << i
711 << ", " << j << ", " << k << ") " << std::endl;
712
713 j = -1;
714 mesh->getGlobalIndex(f, i, j, k, idx);
715 ASSERT_EQ(-1, idx)
716 << "getGlobalIndex wrong at field " << f << " and (" << i
717 << ", " << j << ", " << k << ") " << std::endl;
718 }
719
720 // ghost points at right and left
721 for (int k = 0; k < mesh->n[f][2]; k++)
722 for (int j = 0; j < mesh->n[f][1]; j++)
723 {
724 PetscInt idx, i;
725
726 i = mesh->n[f][0];
727 mesh->getGlobalIndex(f, i, j, k, idx);
728 ASSERT_EQ(-1, idx)
729 << "getGlobalIndex wrong at field " << f << " and (" << i
730 << ", " << j << ", " << k << ") " << std::endl;
731
732 i = -1;
733 mesh->getGlobalIndex(f, i, j, k, idx);
734 ASSERT_EQ(-1, idx)
735 << "getGlobalIndex wrong at field " << f << " and (" << i
736 << ", " << j << ", " << k << ") " << std::endl;
737 }
738
739 // corners
740 for (int k = 0; k < mesh->n[f][2]; k++)
741 {
742 PetscInt idx;
743
744 mesh->getGlobalIndex(f, -1, -1, 0, idx);
745 ASSERT_EQ(-1, idx)
746 << "getGlobalIndex wrong at field " << f << " and (" << -1
747 << ", " << -1 << ", " << k << ") " << std::endl;
748
749 mesh->getGlobalIndex(f, -1, mesh->n[f][1], 0, idx);
750 ASSERT_EQ(-1, idx)
751 << "getGlobalIndex wrong at field " << f << " and (" << -1
752 << ", " << mesh->n[f][1] << ", " << k << ") " << std::endl;
753
754 mesh->getGlobalIndex(f, mesh->n[f][0], -1, 0, idx);
755 ASSERT_EQ(-1, idx) << "getGlobalIndex wrong at field " << f
756 << " and (" << mesh->n[f][0] << ", " << -1
757 << ", " << k << ") " << std::endl;
758
759 mesh->getGlobalIndex(f, mesh->n[f][0], mesh->n[f][1], 0, idx);
760 ASSERT_EQ(-1, idx)
761 << "getGlobalIndex wrong at field " << f << " and ("
762 << mesh->n[f][0] << ", " << mesh->n[f][1] << ", " << k << ") "
763 << std::endl;
764 }
765 }
766}
767
768// check getPackedGlobalIndex with i, j, k of internal points
769TEST_F(CartesianMeshTest2D_AllDirichlet,
770 check_getPackedGlobalIndex_ijk_internal)
771{
772 PetscErrorCode ierr;
773
774 DM UPack;
775
776 ierr = DMCompositeCreate(PETSC_COMM_WORLD, &UPack);
777 ASSERT_FALSE(ierr);
778 ierr = DMCompositeAddDM(UPack, da[0]);
779 ASSERT_FALSE(ierr);
780 ierr = DMCompositeAddDM(UPack, da[1]);
781 ASSERT_FALSE(ierr);
782 ierr = DMSetUp(UPack);
783 ASSERT_FALSE(ierr);
784
785 Vec fp, ip, jp, kp;
786 ierr = DMCreateGlobalVector(mesh->UPack, &fp);
787 ASSERT_FALSE(ierr);
788 ierr = DMCreateGlobalVector(mesh->UPack, &ip);
789 ASSERT_FALSE(ierr);
790 ierr = DMCreateGlobalVector(mesh->UPack, &jp);
791 ASSERT_FALSE(ierr);
792 ierr = DMCreateGlobalVector(mesh->UPack, &kp);
793 ASSERT_FALSE(ierr);
794
795 Vec fv[2], iv[2], jv[2], kv[2];
796 ierr = DMCompositeGetAccessArray(mesh->UPack, fp, 2, nullptr, fv);
797 ASSERT_FALSE(ierr);
798 ierr = DMCompositeGetAccessArray(mesh->UPack, ip, 2, nullptr, iv);
799 ASSERT_FALSE(ierr);
800 ierr = DMCompositeGetAccessArray(mesh->UPack, jp, 2, nullptr, jv);
801 ASSERT_FALSE(ierr);
802 ierr = DMCompositeGetAccessArray(mesh->UPack, kp, 2, nullptr, kv);
803 ASSERT_FALSE(ierr);
804
805 for (int f = 0; f < 2; ++f)
806 {
807 ISLocalToGlobalMapping mapping;
808 ierr = DMGetLocalToGlobalMapping(mesh->da[f], &mapping);
809 ASSERT_FALSE(ierr);
810
811 for (int k = mesh->bg[f][2]; k < mesh->ed[f][2]; ++k)
812 for (int j = mesh->bg[f][1]; j < mesh->ed[f][1]; ++j)
813 for (int i = mesh->bg[f][0]; i < mesh->ed[f][0]; ++i)
814 {
815 PetscInt idx;
816 ierr = DMDAConvertToCell(mesh->da[f], {k, j, i, 0}, &idx);
817 ASSERT_FALSE(ierr);
818 ierr = ISLocalToGlobalMappingApply(mapping, 1, &idx, &idx);
819 ASSERT_FALSE(ierr);
820
821 ierr = VecSetValue(fv[f], idx, (PetscReal)f, INSERT_VALUES);
822 ASSERT_FALSE(ierr);
823 ierr = VecSetValue(iv[f], idx, (PetscReal)i, INSERT_VALUES);
824 ASSERT_FALSE(ierr);
825 ierr = VecSetValue(jv[f], idx, (PetscReal)j, INSERT_VALUES);
826 ASSERT_FALSE(ierr);
827 ierr = VecSetValue(kv[f], idx, (PetscReal)k, INSERT_VALUES);
828 ASSERT_FALSE(ierr);
829 }
830
831 ierr = VecAssemblyBegin(fv[f]);
832 ASSERT_FALSE(ierr);
833 ierr = VecAssemblyBegin(iv[f]);
834 ASSERT_FALSE(ierr);
835 ierr = VecAssemblyBegin(jv[f]);
836 ASSERT_FALSE(ierr);
837 ierr = VecAssemblyBegin(kv[f]);
838 ASSERT_FALSE(ierr);
839 ierr = VecAssemblyEnd(fv[f]);
840 ASSERT_FALSE(ierr);
841 ierr = VecAssemblyEnd(iv[f]);
842 ASSERT_FALSE(ierr);
843 ierr = VecAssemblyEnd(jv[f]);
844 ASSERT_FALSE(ierr);
845 ierr = VecAssemblyEnd(kv[f]);
846 ASSERT_FALSE(ierr);
847
848 ierr = ISLocalToGlobalMappingDestroy(&mapping);
849 ASSERT_FALSE(ierr);
850 }
851
852 ierr = DMCompositeRestoreAccessArray(mesh->UPack, fp, 2, nullptr, fv);
853 ASSERT_FALSE(ierr);
854 ierr = DMCompositeRestoreAccessArray(mesh->UPack, ip, 2, nullptr, iv);
855 ASSERT_FALSE(ierr);
856 ierr = DMCompositeRestoreAccessArray(mesh->UPack, jp, 2, nullptr, jv);
857 ASSERT_FALSE(ierr);
858 ierr = DMCompositeRestoreAccessArray(mesh->UPack, kp, 2, nullptr, kv);
859 ASSERT_FALSE(ierr);
860
861 VecScatter ctx;
862 Vec fps, ips, jps, kps;
863
864 ierr = VecScatterCreateToAll(fp, &ctx, &fps);
865 ASSERT_FALSE(ierr);
866 ierr = VecDuplicate(fps, &ips);
867 ASSERT_FALSE(ierr);
868 ierr = VecDuplicate(fps, &jps);
869 ASSERT_FALSE(ierr);
870 ierr = VecDuplicate(fps, &kps);
871 ASSERT_FALSE(ierr);
872
873 ierr = VecScatterBegin(ctx, fp, fps, INSERT_VALUES, SCATTER_FORWARD);
874 ASSERT_FALSE(ierr);
875 ierr = VecScatterEnd(ctx, fp, fps, INSERT_VALUES, SCATTER_FORWARD);
876 ASSERT_FALSE(ierr);
877 ierr = VecScatterBegin(ctx, ip, ips, INSERT_VALUES, SCATTER_FORWARD);
878 ASSERT_FALSE(ierr);
879 ierr = VecScatterEnd(ctx, ip, ips, INSERT_VALUES, SCATTER_FORWARD);
880 ASSERT_FALSE(ierr);
881 ierr = VecScatterBegin(ctx, jp, jps, INSERT_VALUES, SCATTER_FORWARD);
882 ASSERT_FALSE(ierr);
883 ierr = VecScatterEnd(ctx, jp, jps, INSERT_VALUES, SCATTER_FORWARD);
884 ASSERT_FALSE(ierr);
885 ierr = VecScatterBegin(ctx, kp, kps, INSERT_VALUES, SCATTER_FORWARD);
886 ASSERT_FALSE(ierr);
887 ierr = VecScatterEnd(ctx, kp, kps, INSERT_VALUES, SCATTER_FORWARD);
888 ASSERT_FALSE(ierr);
889
890 ierr = VecScatterDestroy(&ctx);
891 ASSERT_FALSE(ierr);
892 ierr = VecDestroy(&kp);
893 ASSERT_FALSE(ierr);
894 ierr = VecDestroy(&jp);
895 ASSERT_FALSE(ierr);
896 ierr = VecDestroy(&ip);
897 ASSERT_FALSE(ierr);
898 ierr = VecDestroy(&fp);
899 ASSERT_FALSE(ierr);
900
901 for (int f = 0; f < 2; ++f)
902 {
903 // internal points
904 for (int k = 0; k < mesh->n[f][2]; k++)
905 for (int j = 0; j < mesh->n[f][1]; j++)
906 for (int i = 0; i < mesh->n[f][0]; i++)
907 {
908 PetscInt idx;
909 ierr = mesh->getPackedGlobalIndex(f, i, j, k, idx);
910 ASSERT_FALSE(ierr);
911
912 PetscReal vf, vi, vj, vk;
913 ierr = VecGetValues(fps, 1, &idx, &vf);
914 ASSERT_FALSE(ierr);
915 ierr = VecGetValues(ips, 1, &idx, &vi);
916 ASSERT_FALSE(ierr);
917 ierr = VecGetValues(jps, 1, &idx, &vj);
918 ASSERT_FALSE(ierr);
919 ierr = VecGetValues(kps, 1, &idx, &vk);
920 ASSERT_FALSE(ierr);
921
922 ASSERT_EQ(vf, f) << f << ".(" << i << ", " << j << ", " << k
923 << ")" << std::endl;
924 ASSERT_EQ(vi, i) << f << ".(" << i << ", " << j << ", " << k
925 << ")" << std::endl;
926 ASSERT_EQ(vj, j) << f << ".(" << i << ", " << j << ", " << k
927 << ")" << std::endl;
928 ASSERT_EQ(vk, k) << f << ".(" << i << ", " << j << ", " << k
929 << ")" << std::endl;
930 }
931 }
932
933 ierr = VecDestroy(&kps);
934 ASSERT_FALSE(ierr);
935 ierr = VecDestroy(&jps);
936 ASSERT_FALSE(ierr);
937 ierr = VecDestroy(&ips);
938 ASSERT_FALSE(ierr);
939 ierr = VecDestroy(&fps);
940 ASSERT_FALSE(ierr);
941}
942
943// check getPackedGlobalIndex with i, j, k of ghosts
944TEST_F(CartesianMeshTest2D_AllDirichlet, check_getPackedGlobalIndex_ijk_ghost)
945{
946 for (int f = 0; f < 2; ++f)
947 {
948 // ghost points at top and bottom
949 for (int k = 0; k < mesh->n[f][2]; k++)
950 for (int i = 0; i < mesh->n[f][0]; i++)
951 {
952 PetscInt idx, j;
953
954 j = mesh->n[f][1];
955 mesh->getPackedGlobalIndex(f, i, j, k, idx);
956 ASSERT_EQ(-1, idx)
957 << "getPackedGlobalIndex wrong at field " << f << " and ("
958 << i << ", " << j << ", " << k << ") " << std::endl;
959
960 j = -1;
961 mesh->getPackedGlobalIndex(f, i, j, k, idx);
962 ASSERT_EQ(-1, idx)
963 << "getPackedGlobalIndex wrong at field " << f << " and ("
964 << i << ", " << j << ", " << k << ") " << std::endl;
965 }
966
967 // ghost points at right and left
968 for (int k = 0; k < mesh->n[f][2]; k++)
969 for (int j = 0; j < mesh->n[f][1]; j++)
970 {
971 PetscInt idx, i;
972
973 i = mesh->n[f][0];
974 mesh->getPackedGlobalIndex(f, i, j, k, idx);
975 ASSERT_EQ(-1, idx)
976 << "getPackedGlobalIndex wrong at field " << f << " and ("
977 << i << ", " << j << ", " << k << ") " << std::endl;
978
979 i = -1;
980 mesh->getPackedGlobalIndex(f, i, j, k, idx);
981 ASSERT_EQ(-1, idx)
982 << "getPackedGlobalIndex wrong at field " << f << " and ("
983 << i << ", " << j << ", " << k << ") " << std::endl;
984 }
985
986 // corners
987 for (int k = 0; k < mesh->n[f][2]; k++)
988 {
989 PetscInt idx;
990
991 mesh->getPackedGlobalIndex(f, -1, -1, 0, idx);
992 ASSERT_EQ(-1, idx)
993 << "getPackedGlobalIndex wrong at field " << f << " and (" << -1
994 << ", " << -1 << ", " << k << ") " << std::endl;
995
996 mesh->getPackedGlobalIndex(f, -1, mesh->n[f][1], 0, idx);
997 ASSERT_EQ(-1, idx)
998 << "getPackedGlobalIndex wrong at field " << f << " and (" << -1
999 << ", " << mesh->n[f][1] << ", " << k << ") " << std::endl;
1000
1001 mesh->getPackedGlobalIndex(f, mesh->n[f][0], -1, 0, idx);
1002 ASSERT_EQ(-1, idx) << "getPackedGlobalIndex wrong at field " << f
1003 << " and (" << mesh->n[f][0] << ", " << -1
1004 << ", " << k << ") " << std::endl;
1005
1006 mesh->getPackedGlobalIndex(f, mesh->n[f][0], mesh->n[f][1], 0, idx);
1007 ASSERT_EQ(-1, idx)
1008 << "getPackedGlobalIndex wrong at field " << f << " and ("
1009 << mesh->n[f][0] << ", " << mesh->n[f][1] << ", " << k << ") "
1010 << std::endl;
1011 }
1012 }
1013}
TEST_F
TEST_F(CartesianMeshTest2D_AllDirichlet, check_dim)
Definition: cartesianmesh2d_dirichlet.cpp:130
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::type::BoolVec2D
std::vector< BoolVec1D > BoolVec2D
2D std::vector holding PetscBool.
Definition: type.h:149
petibm::type::RealVec2D
std::vector< RealVec1D > RealVec2D
2D std::vector holding PetscReal.
Definition: type.h:142
petibm::type::RealVec3D
std::vector< RealVec2D > RealVec3D
3D std::vector holding PetscReal.
Definition: type.h:144
mesh.h
Prototype of mesh::MeshBase, type::Mesh, and factory function.
YAML
Definition: parser.cpp:45
petibm::type::GhostedVec3D
std::vector< GhostedVec2D > GhostedVec3D
a vector of vector pointers to mimic ghosted 2D vectors. type
Definition: type.h:157