PetIBM 0.5.4
Toolbox and applications of the immersed-boundary method for distributed-memory architectures
createconvection.cpp
Go to the documentation of this file.
1
8// TODO: investigate if exact interpolation is necessary
9
10// STL
11#include <memory>
12
13// PETSc
14#include <petscmat.h>
15
16// PetIBM
17#include <petibm/boundary.h>
18#include <petibm/mesh.h>
19
20namespace // anonymous namespace for internal linkage only
21{
22// a private struct used in MatShell
23struct NonLinearCtx
24{
25 const petibm::type::Mesh mesh;
26 const petibm::type::Boundary bc;
27 std::vector<Vec> qLocal;
28
29 NonLinearCtx(const petibm::type::Mesh &_mesh,
30 const petibm::type::Boundary &_bc)
31 : mesh(_mesh), bc(_bc), qLocal(_mesh->dim)
32 {
33 // create necessary local vectors
34 for (PetscInt f = 0; f < mesh->dim; ++f)
35 DMCreateLocalVector(mesh->da[f], &qLocal[f]);
36 };
37};
38
39// a private kernel for the convection at a u-velocity point in 2D.
40inline PetscReal kernelU(NonLinearCtx const *const &ctx,
41 const std::vector<PetscReal **> &flux,
42 const PetscInt &i, const PetscInt &j)
43{
44 PetscReal uSelf;
45 PetscReal uS, uN, uW, uE;
46 PetscReal vS, vN;
47
48 // prepare self
49 uSelf = flux[0][j][i];
50
51 // prepare u
52 uW = (uSelf + flux[0][j][i - 1]) / 2.0;
53 uE = (uSelf + flux[0][j][i + 1]) / 2.0;
54 uS = (uSelf + flux[0][j - 1][i]) / 2.0;
55 uN = (uSelf + flux[0][j + 1][i]) / 2.0;
56
57 // prepare v
58 vS = (flux[1][j - 1][i] + flux[1][j - 1][i + 1]) / 2.0;
59 vN = (flux[1][j][i] + flux[1][j][i + 1]) / 2.0;
60
61 return (uE * uE - uW * uW) / ctx->mesh->dL[0][0][i] +
62 (vN * uN - vS * uS) / ctx->mesh->dL[0][1][j];
63} // kernelU
64
65// a private kernel for the convection at a v-velocity point in 2D
66inline PetscReal kernelV(NonLinearCtx const *const &ctx,
67 const std::vector<PetscReal **> &flux,
68 const PetscInt &i, const PetscInt &j)
69{
70 PetscReal vSelf;
71 PetscReal uW, uE;
72 PetscReal vS, vN, vW, vE;
73
74 // prepare self
75 vSelf = flux[1][j][i];
76
77 // prepare u
78 uW = (flux[0][j][i - 1] + flux[0][j + 1][i - 1]) / 2.0;
79 uE = (flux[0][j][i] + flux[0][j + 1][i]) / 2.0;
80
81 // prepare v
82 vW = (vSelf + flux[1][j][i - 1]) / 2.0;
83 vE = (vSelf + flux[1][j][i + 1]) / 2.0;
84 vS = (vSelf + flux[1][j - 1][i]) / 2.0;
85 vN = (vSelf + flux[1][j + 1][i]) / 2.0;
86
87 return (uE * vE - uW * vW) / ctx->mesh->dL[1][0][i] +
88 (vN * vN - vS * vS) / ctx->mesh->dL[1][1][j];
89} // kernelV
90
91// a private kernel for the convection at a u-velocity point in 3D
92inline PetscReal kernelU(NonLinearCtx const *const &ctx,
93 const std::vector<PetscReal ***> &flux,
94 const PetscInt &i, const PetscInt &j,
95 const PetscInt &k)
96{
97 PetscReal uSelf;
98 PetscReal uS, uN, uW, uE, uB, uF;
99 PetscReal vS, vN;
100 PetscReal wB, wF;
101
102 // prepare self
103 uSelf = flux[0][k][j][i];
104
105 // prepare u
106 uW = (uSelf + flux[0][k][j][i - 1]) / 2.0;
107 uE = (uSelf + flux[0][k][j][i + 1]) / 2.0;
108 uS = (uSelf + flux[0][k][j - 1][i]) / 2.0;
109 uN = (uSelf + flux[0][k][j + 1][i]) / 2.0;
110 uB = (uSelf + flux[0][k - 1][j][i]) / 2.0;
111 uF = (uSelf + flux[0][k + 1][j][i]) / 2.0;
112
113 // prepare v
114 vS = (flux[1][k][j - 1][i] + flux[1][k][j - 1][i + 1]) / 2.0;
115 vN = (flux[1][k][j][i] + flux[1][k][j][i + 1]) / 2.0;
116
117 // prepare w
118 wB = (flux[2][k - 1][j][i] + flux[2][k - 1][j][i + 1]) / 2.0;
119 wF = (flux[2][k][j][i] + flux[2][k][j][i + 1]) / 2.0;
120
121 return (uE * uE - uW * uW) / ctx->mesh->dL[0][0][i] +
122 (vN * uN - vS * uS) / ctx->mesh->dL[0][1][j] +
123 (wF * uF - wB * uB) / ctx->mesh->dL[0][2][k];
124} // kernelU
125
126// a private kernel for the convection at a v-velocity point in 3D
127inline PetscReal kernelV(NonLinearCtx const *const &ctx,
128 const std::vector<PetscReal ***> &flux,
129 const PetscInt &i, const PetscInt &j,
130 const PetscInt &k)
131{
132 PetscReal vSelf;
133 PetscReal uW, uE;
134 PetscReal vS, vN, vW, vE, vB, vF;
135 PetscReal wB, wF;
136
137 // prepare self
138 vSelf = flux[1][k][j][i];
139
140 // prepare u
141 uW = (flux[0][k][j][i - 1] + flux[0][k][j + 1][i - 1]) / 2.0;
142 uE = (flux[0][k][j][i] + flux[0][k][j + 1][i]) / 2.0;
143
144 // prepare v
145 vW = (vSelf + flux[1][k][j][i - 1]) / 2.0;
146 vE = (vSelf + flux[1][k][j][i + 1]) / 2.0;
147 vS = (vSelf + flux[1][k][j - 1][i]) / 2.0;
148 vN = (vSelf + flux[1][k][j + 1][i]) / 2.0;
149 vB = (vSelf + flux[1][k - 1][j][i]) / 2.0;
150 vF = (vSelf + flux[1][k + 1][j][i]) / 2.0;
151
152 // prepare w
153 wB = (flux[2][k - 1][j][i] + flux[2][k - 1][j + 1][i]) / 2.0;
154 wF = (flux[2][k][j][i] + flux[2][k][j + 1][i]) / 2.0;
155
156 return (uE * vE - uW * vW) / ctx->mesh->dL[1][0][i] +
157 (vN * vN - vS * vS) / ctx->mesh->dL[1][1][j] +
158 (wF * vF - wB * vB) / ctx->mesh->dL[1][2][k];
159
160} // kernelV
161
162// a private kernel for the convection at a w-velocity point in 3D
163inline PetscReal kernelW(NonLinearCtx const *const &ctx,
164 const std::vector<PetscReal ***> &flux,
165 const PetscInt &i, const PetscInt &j,
166 const PetscInt &k)
167{
168 PetscReal wSelf;
169 PetscReal uW, uE;
170 PetscReal vS, vN;
171 PetscReal wS, wN, wW, wE, wB, wF;
172
173 // prepare self
174 wSelf = flux[2][k][j][i];
175
176 // prepare u
177 uW = (flux[0][k][j][i - 1] + flux[0][k + 1][j][i - 1]) / 2.0;
178 uE = (flux[0][k][j][i] + flux[0][k + 1][j][i]) / 2.0;
179
180 // prepare v
181 vS = (flux[1][k][j - 1][i] + flux[1][k + 1][j - 1][i]) / 2.0;
182 vN = (flux[1][k][j][i] + flux[1][k + 1][j][i]) / 2.0;
183
184 // prepare w
185 wW = (wSelf + flux[2][k][j][i - 1]) / 2.0;
186 wE = (wSelf + flux[2][k][j][i + 1]) / 2.0;
187 wS = (wSelf + flux[2][k][j - 1][i]) / 2.0;
188 wN = (wSelf + flux[2][k][j + 1][i]) / 2.0;
189 wB = (wSelf + flux[2][k - 1][j][i]) / 2.0;
190 wF = (wSelf + flux[2][k + 1][j][i]) / 2.0;
191
192 return (uE * wE - uW * wW) / ctx->mesh->dL[2][0][i] +
193 (vN * wN - vS * wS) / ctx->mesh->dL[2][1][j] +
194 (wF * wF - wB * wB) / ctx->mesh->dL[2][2][k];
195} // kernelW
196
197// a private function for convection operator's MatMult in 2D.
198// For user-defined PETSc Mat operations, please refer to PETSc manual.
199PetscErrorCode ConvectionMult2D(Mat mat, Vec x, Vec y)
200{
201 PetscFunctionBeginUser;
202
203 PetscErrorCode ierr;
204
205 NonLinearCtx *ctx;
206
207 std::vector<Vec> unPacked(2);
208
209 std::vector<PetscReal **> xArry(2);
210
211 std::vector<PetscReal **> yArry(2);
212
213 // get the context
214 ierr = MatShellGetContext(mat, (void *)&ctx); CHKERRQ(ierr);
215
216 // get local (including overlapped points) values of x
217 ierr = DMCompositeScatterArray(ctx->mesh->UPack, x, ctx->qLocal.data());
218 CHKERRQ(ierr);
219
220 // set the values of ghost points in local vectors
221 ierr = ctx->bc->copyValues2LocalVecs(ctx->qLocal); CHKERRQ(ierr);
222
223 // get unPacked vectors of y
224 ierr = DMCompositeGetAccessArray(ctx->mesh->UPack, y, ctx->mesh->dim,
225 nullptr, unPacked.data()); CHKERRQ(ierr);
226
227 // get underlying data of Vecs
228 for (PetscInt f = 0; f < ctx->mesh->dim; ++f)
229 {
230 ierr = DMDAVecGetArrayRead(ctx->mesh->da[f], ctx->qLocal[f], &xArry[f]);
231 CHKERRQ(ierr);
232
233 ierr = DMDAVecGetArray(ctx->mesh->da[f], unPacked[f], &yArry[f]);
234 CHKERRQ(ierr);
235 }
236
237 // u-velocity field
238 for (PetscInt j = ctx->mesh->bg[0][1]; j < ctx->mesh->ed[0][1]; ++j)
239 for (PetscInt i = ctx->mesh->bg[0][0]; i < ctx->mesh->ed[0][0]; ++i)
240 yArry[0][j][i] = kernelU(ctx, xArry, i, j);
241
242 // v-velocity field
243 for (PetscInt j = ctx->mesh->bg[1][1]; j < ctx->mesh->ed[1][1]; ++j)
244 for (PetscInt i = ctx->mesh->bg[1][0]; i < ctx->mesh->ed[1][0]; ++i)
245 yArry[1][j][i] = kernelV(ctx, xArry, i, j);
246
247 // return underlying arrays of Vecs
248 for (PetscInt f = 0; f < ctx->mesh->dim; ++f)
249 {
250 ierr = DMDAVecRestoreArrayRead(ctx->mesh->da[f], ctx->qLocal[f],
251 &xArry[f]); CHKERRQ(ierr);
252
253 ierr = DMDAVecRestoreArray(ctx->mesh->da[f], unPacked[f], &yArry[f]);
254 CHKERRQ(ierr);
255 }
256
257 // return the ownership of unpacked vectors back to packed vectors
258 ierr = DMCompositeRestoreAccessArray(ctx->mesh->UPack, y, ctx->mesh->dim,
259 nullptr, unPacked.data());
260 CHKERRQ(ierr);
261
262 PetscFunctionReturn(0);
263} // ConvectionMult2D
264
265// a private function for convection operator's MatMult in 3D.
266// For user-defined PETSc Mat operations, please refer to PETSc manual.
267PetscErrorCode ConvectionMult3D(Mat mat, Vec x, Vec y)
268{
269 PetscFunctionBeginUser;
270
271 PetscErrorCode ierr;
272
273 NonLinearCtx *ctx;
274
275 std::vector<Vec> unPacked(3);
276
277 std::vector<PetscReal ***> xArry(3);
278
279 std::vector<PetscReal ***> yArry(3);
280
281 // get the context
282 ierr = MatShellGetContext(mat, (void *)&ctx); CHKERRQ(ierr);
283
284 // get local (including overlapped points) vectors of x
285 ierr = DMCompositeScatterArray(ctx->mesh->UPack, x, ctx->qLocal.data());
286 CHKERRQ(ierr);
287
288 // set the values of ghost points in local vectors
289 ierr = ctx->bc->copyValues2LocalVecs(ctx->qLocal); CHKERRQ(ierr);
290
291 // get unPacked vectors of y
292 ierr = DMCompositeGetAccessArray(ctx->mesh->UPack, y, ctx->mesh->dim,
293 nullptr, unPacked.data()); CHKERRQ(ierr);
294
295 // get underlying data of Vecs
296 for (PetscInt f = 0; f < ctx->mesh->dim; ++f)
297 {
298 ierr = DMDAVecGetArrayRead(ctx->mesh->da[f], ctx->qLocal[f], &xArry[f]);
299 CHKERRQ(ierr);
300
301 ierr = DMDAVecGetArray(ctx->mesh->da[f], unPacked[f], &yArry[f]);
302 CHKERRQ(ierr);
303 }
304
305 // u-velocity field
306 for (PetscInt k = ctx->mesh->bg[0][2]; k < ctx->mesh->ed[0][2]; ++k)
307 for (PetscInt j = ctx->mesh->bg[0][1]; j < ctx->mesh->ed[0][1]; ++j)
308 for (PetscInt i = ctx->mesh->bg[0][0]; i < ctx->mesh->ed[0][0]; ++i)
309 yArry[0][k][j][i] = kernelU(ctx, xArry, i, j, k);
310
311 // v-velocity field
312 for (PetscInt k = ctx->mesh->bg[1][2]; k < ctx->mesh->ed[1][2]; ++k)
313 for (PetscInt j = ctx->mesh->bg[1][1]; j < ctx->mesh->ed[1][1]; ++j)
314 for (PetscInt i = ctx->mesh->bg[1][0]; i < ctx->mesh->ed[1][0]; ++i)
315 yArry[1][k][j][i] = kernelV(ctx, xArry, i, j, k);
316
317 // w-velocity field
318 for (PetscInt k = ctx->mesh->bg[2][2]; k < ctx->mesh->ed[2][2]; ++k)
319 for (PetscInt j = ctx->mesh->bg[2][1]; j < ctx->mesh->ed[2][1]; ++j)
320 for (PetscInt i = ctx->mesh->bg[2][0]; i < ctx->mesh->ed[2][0]; ++i)
321 yArry[2][k][j][i] = kernelW(ctx, xArry, i, j, k);
322
323 // return underlying arrays of Vecs
324 for (PetscInt f = 0; f < ctx->mesh->dim; ++f)
325 {
326 ierr = DMDAVecRestoreArrayRead(ctx->mesh->da[f], ctx->qLocal[f],
327 &xArry[f]); CHKERRQ(ierr);
328
329 ierr = DMDAVecRestoreArray(ctx->mesh->da[f], unPacked[f], &yArry[f]);
330 CHKERRQ(ierr);
331 }
332
333 // return the ownership of unpacked vectors back to packed vectors
334 ierr = DMCompositeRestoreAccessArray(ctx->mesh->UPack, y, ctx->mesh->dim,
335 nullptr, unPacked.data());
336 CHKERRQ(ierr);
337
338 PetscFunctionReturn(0);
339} // ConvectionMult3D
340
341// a private function for convection operator's destroying.
342// For user-defined PETSc Mat operations, please refer to PETSc manual.
343PetscErrorCode ConvectionDestroy(Mat mat)
344{
345 PetscFunctionBeginUser;
346
347 PetscErrorCode ierr;
348
349 NonLinearCtx *ctx;
350
351 // get the context
352 ierr = MatShellGetContext(mat, (void *)&ctx); CHKERRQ(ierr);
353
354 // destroy qLocal
355 for (PetscInt f = 0; f < ctx->mesh->dim; f++)
356 {
357 ierr = VecDestroy(&ctx->qLocal[f]); CHKERRQ(ierr);
358 }
359
360 // deallocate the memory space pointed by ctx
361 delete ctx;
362
363 PetscFunctionReturn(0);
364} // ConvectionDestroy
365} // end of anonymous namespace
366
367namespace petibm
368{
369namespace operators
370{
371// implementation of petibm::operators::createConvection
372PetscErrorCode createConvection(const type::Mesh &mesh,
373 const type::Boundary &bd, Mat &H)
374{
375 PetscFunctionBeginUser;
376
377 PetscErrorCode ierr;
378
379 NonLinearCtx *ctx;
380
381 // global dimension
382 PetscInt N =
383 mesh->n[0][0] * mesh->n[0][1] * mesh->n[0][2] +
384 mesh->n[1][0] * mesh->n[1][1] * mesh->n[1][2] +
385 ((mesh->dim == 3) ? mesh->n[2][0] * mesh->n[2][1] * mesh->n[2][2] : 0);
386
387 // allocate space for ctx
388 ctx = new NonLinearCtx(mesh, bd);
389
390 // create a matrix-free operator
391 ierr = MatCreateShell(mesh->comm, mesh->UNLocal, mesh->UNLocal, N, N,
392 (void *)ctx, &H); CHKERRQ(ierr);
393
394 // bind MatMult
395 if (mesh->dim == 2)
396 {
397 ierr = MatShellSetOperation(H, MATOP_MULT,
398 (void (*)(void))ConvectionMult2D);
399 CHKERRQ(ierr);
400 }
401 else // assume the dim is either 2 or 3
402 {
403 ierr = MatShellSetOperation(H, MATOP_MULT,
404 (void (*)(void))ConvectionMult3D);
405 CHKERRQ(ierr);
406 }
407
408 // bind MatDestroy
409 ierr = MatShellSetOperation(H, MATOP_DESTROY,
410 (void (*)(void))ConvectionDestroy);
411 CHKERRQ(ierr);
412
413 PetscFunctionReturn(0);
414} // createConvection
415
416} // end of namespace operators
417} // end of namespace petibm
boundary.h
boundary::BoundaryBase, type::Boundary, and factory function.
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::operators::createConvection
PetscErrorCode createConvection(const type::Mesh &mesh, const type::Boundary &bc, Mat &H)
Create a matrix-free Mat for convection operator, .
Definition: createconvection.cpp:372
petibm::type::x
@ x
Definition: type.h:69
petibm::type::y
@ y
Definition: type.h:70
mesh.h
Prototype of mesh::MeshBase, type::Mesh, and factory function.
petibm
A toolbox for building flow solvers.
Definition: bodypack.h:52