1 static char help[] = "Poiseuille flow problem. Viscous, laminar flow in a 2D channel with parabolic velocity\n\
2 profile and linear pressure drop, exact solution of the 2D Stokes\n";
3
4 /*
5 M A R I T I M E R E S E A R C H I N S T I T U T E N E T H E R L A N D S
6 author : Christiaan M. Klaij
7
8 Poiseuille flow problem.
9
10 Viscous, laminar flow in a 2D channel with parabolic velocity
11 profile and linear pressure drop, exact solution of the 2D Stokes
12 equations.
13
14 Discretized with the cell-centered finite-volume method on a
15 Cartesian grid with co-located variables. Variables ordered as
16 [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with
17 PCFIELDSPLIT. Lower factorization is used to mimic the Semi-Implicit
18 Method for Pressure Linked Equations (SIMPLE) used as preconditioner
19 instead of solver.
20
21 Disclaimer: does not contain the pressure-weighed interpolation
22 method needed to suppress spurious pressure modes in real-life
23 problems.
24
25 Usage:
26 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none
27
28 Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur
29 complement because A11 is zero. FGMRES is needed because
30 PCFIELDSPLIT is a variable preconditioner.
31
32 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
33
34 Same as above but with user defined PC for the true Schur
35 complement. PC based on the SIMPLE-type approximation (inverse of
36 A00 approximated by inverse of its diagonal).
37
38 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_ksp
39
40 Replace the true Schur complement with a user defined Schur
41 complement based on the SIMPLE-type approximation. Same matrix is
42 used as PC.
43
44 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi
45
46 Out-of-the-box SIMPLE-type preconditioning. The major advantage
47 is that the user neither needs to provide the approximation of
48 the Schur complement, nor the corresponding preconditioner.
49 */
50
51 #include <petscksp.h>
52
53 typedef struct {
54 PetscBool userPC, userKSP, matsymmetric; /* user defined preconditioner and matrix for the Schur complement */
55 PetscInt nx, ny; /* nb of cells in x- and y-direction */
56 PetscReal hx, hy; /* mesh size in x- and y-direction */
57 Mat A; /* block matrix */
58 Mat subA[4]; /* the four blocks */
59 Mat myS; /* the approximation of the Schur complement */
60 Vec x, b, y; /* solution, rhs and temporary vector */
61 IS isg[2]; /* index sets of split "0" and "1" */
62 } Stokes;
63
64 PetscErrorCode StokesSetupMatBlock00(Stokes *); /* setup the block Q */
65 PetscErrorCode StokesSetupMatBlock01(Stokes *); /* setup the block G */
66 PetscErrorCode StokesSetupMatBlock10(Stokes *); /* setup the block D (equal to the transpose of G) */
67 PetscErrorCode StokesSetupMatBlock11(Stokes *); /* setup the block C (equal to zero) */
68
69 PetscErrorCode StokesGetPosition(Stokes *, PetscInt, PetscInt *, PetscInt *); /* row number j*nx+i corresponds to position (i,j) in grid */
70
71 PetscErrorCode StokesStencilLaplacian(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Laplacian operator */
72 PetscErrorCode StokesStencilGradientX(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (x-component) */
73 PetscErrorCode StokesStencilGradientY(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (y-component) */
74
75 PetscErrorCode StokesRhs(Stokes *); /* rhs vector */
76 PetscErrorCode StokesRhsMomX(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (x-component) */
77 PetscErrorCode StokesRhsMomY(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (y-component) */
78 PetscErrorCode StokesRhsMass(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of pressure */
79
80 PetscErrorCode StokesSetupApproxSchur(Stokes *); /* approximation of the Schur complement */
81
82 PetscErrorCode StokesExactSolution(Stokes *); /* exact solution vector */
83 PetscErrorCode StokesWriteSolution(Stokes *); /* write solution to file */
84
85 /* exact solution for the velocity (x-component, y-component is zero) */
StokesExactVelocityX(const PetscScalar y)86 PetscScalar StokesExactVelocityX(const PetscScalar y)
87 {
88 return 4.0 * y * (1.0 - y);
89 }
90
91 /* exact solution for the pressure */
StokesExactPressure(const PetscScalar x)92 PetscScalar StokesExactPressure(const PetscScalar x)
93 {
94 return 8.0 * (2.0 - x);
95 }
96
StokesSetupPC(Stokes * s,KSP ksp)97 PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp)
98 {
99 KSP *subksp;
100 PC pc;
101 PetscInt n = 1;
102
103 PetscFunctionBeginUser;
104 PetscCall(KSPGetPC(ksp, &pc));
105 PetscCall(PCFieldSplitSetIS(pc, "0", s->isg[0]));
106 PetscCall(PCFieldSplitSetIS(pc, "1", s->isg[1]));
107 if (s->userPC) PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS));
108 if (s->userKSP) {
109 PetscCall(PCSetUp(pc));
110 PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
111 PetscCall(KSPSetOperators(subksp[1], s->myS, s->myS));
112 PetscCall(PetscFree(subksp));
113 }
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
StokesWriteSolution(Stokes * s)117 PetscErrorCode StokesWriteSolution(Stokes *s)
118 {
119 PetscMPIInt size;
120 PetscInt n, i, j;
121 const PetscScalar *array;
122
123 PetscFunctionBeginUser;
124 /* write data (*warning* only works sequential) */
125 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
126 if (size == 1) {
127 PetscViewer viewer;
128 PetscCall(VecGetArrayRead(s->x, &array));
129 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer));
130 PetscCall(PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n"));
131 for (j = 0; j < s->ny; j++) {
132 for (i = 0; i < s->nx; i++) {
133 n = j * s->nx + i;
134 PetscCall(PetscViewerASCIIPrintf(viewer, "%.12g %.12g %.12g %.12g %.12g\n", (double)(i * s->hx + s->hx / 2), (double)(j * s->hy + s->hy / 2), (double)PetscRealPart(array[n]), (double)PetscRealPart(array[n + s->nx * s->ny]), (double)PetscRealPart(array[n + 2 * s->nx * s->ny])));
135 }
136 }
137 PetscCall(VecRestoreArrayRead(s->x, &array));
138 PetscCall(PetscViewerDestroy(&viewer));
139 }
140 PetscFunctionReturn(PETSC_SUCCESS);
141 }
142
StokesSetupIndexSets(Stokes * s)143 PetscErrorCode StokesSetupIndexSets(Stokes *s)
144 {
145 PetscFunctionBeginUser;
146 /* the two index sets */
147 PetscCall(MatNestGetISs(s->A, s->isg, NULL));
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
StokesSetupVectors(Stokes * s)151 PetscErrorCode StokesSetupVectors(Stokes *s)
152 {
153 PetscFunctionBeginUser;
154 /* solution vector x */
155 PetscCall(VecCreate(PETSC_COMM_WORLD, &s->x));
156 PetscCall(VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny));
157 PetscCall(VecSetType(s->x, VECMPI));
158
159 /* exact solution y */
160 PetscCall(VecDuplicate(s->x, &s->y));
161 PetscCall(StokesExactSolution(s));
162
163 /* rhs vector b */
164 PetscCall(VecDuplicate(s->x, &s->b));
165 PetscCall(StokesRhs(s));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
StokesGetPosition(Stokes * s,PetscInt row,PetscInt * i,PetscInt * j)169 PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
170 {
171 PetscInt n;
172
173 PetscFunctionBeginUser;
174 /* cell number n=j*nx+i has position (i,j) in grid */
175 n = row % (s->nx * s->ny);
176 *i = n % s->nx;
177 *j = (n - (*i)) / s->nx;
178 PetscFunctionReturn(PETSC_SUCCESS);
179 }
180
StokesExactSolution(Stokes * s)181 PetscErrorCode StokesExactSolution(Stokes *s)
182 {
183 PetscInt row, start, end, i, j;
184 PetscScalar val;
185 Vec y0, y1;
186
187 PetscFunctionBeginUser;
188 /* velocity part */
189 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
190 PetscCall(VecGetOwnershipRange(y0, &start, &end));
191 for (row = start; row < end; row++) {
192 PetscCall(StokesGetPosition(s, row, &i, &j));
193 if (row < s->nx * s->ny) {
194 val = StokesExactVelocityX(j * s->hy + s->hy / 2);
195 } else {
196 val = 0;
197 }
198 PetscCall(VecSetValue(y0, row, val, INSERT_VALUES));
199 }
200 PetscCall(VecAssemblyBegin(y0));
201 PetscCall(VecAssemblyEnd(y0));
202 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
203
204 /* pressure part */
205 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
206 PetscCall(VecGetOwnershipRange(y1, &start, &end));
207 for (row = start; row < end; row++) {
208 PetscCall(StokesGetPosition(s, row, &i, &j));
209 val = StokesExactPressure(i * s->hx + s->hx / 2);
210 PetscCall(VecSetValue(y1, row, val, INSERT_VALUES));
211 }
212 PetscCall(VecAssemblyBegin(y1));
213 PetscCall(VecAssemblyEnd(y1));
214 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217
StokesRhs(Stokes * s)218 PetscErrorCode StokesRhs(Stokes *s)
219 {
220 PetscInt row, start, end, i, j;
221 PetscScalar val;
222 Vec b0, b1;
223
224 PetscFunctionBeginUser;
225 /* velocity part */
226 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
227 PetscCall(VecGetOwnershipRange(b0, &start, &end));
228 for (row = start; row < end; row++) {
229 PetscCall(StokesGetPosition(s, row, &i, &j));
230 if (row < s->nx * s->ny) {
231 PetscCall(StokesRhsMomX(s, i, j, &val));
232 } else {
233 PetscCall(StokesRhsMomY(s, i, j, &val));
234 }
235 PetscCall(VecSetValue(b0, row, val, INSERT_VALUES));
236 }
237 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
238
239 /* pressure part */
240 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
241 PetscCall(VecGetOwnershipRange(b1, &start, &end));
242 for (row = start; row < end; row++) {
243 PetscCall(StokesGetPosition(s, row, &i, &j));
244 PetscCall(StokesRhsMass(s, i, j, &val));
245 if (s->matsymmetric) val = -1.0 * val;
246 PetscCall(VecSetValue(b1, row, val, INSERT_VALUES));
247 }
248 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
249 PetscFunctionReturn(PETSC_SUCCESS);
250 }
251
StokesSetupMatBlock00(Stokes * s)252 PetscErrorCode StokesSetupMatBlock00(Stokes *s)
253 {
254 PetscInt row, start, end, sz, i, j;
255 PetscInt cols[5];
256 PetscScalar vals[5];
257
258 PetscFunctionBeginUser;
259 /* A[0] is 2N-by-2N */
260 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0]));
261 PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_"));
262 PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny));
263 PetscCall(MatSetType(s->subA[0], MATMPIAIJ));
264 PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL));
265 PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end));
266
267 for (row = start; row < end; row++) {
268 PetscCall(StokesGetPosition(s, row, &i, &j));
269 /* first part: rows 0 to (nx*ny-1) */
270 PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals));
271 /* second part: rows (nx*ny) to (2*nx*ny-1) */
272 if (row >= s->nx * s->ny) {
273 for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny;
274 }
275 for (i = 0; i < sz; i++) vals[i] = -1.0 * vals[i]; /* dynamic viscosity coef mu=-1 */
276 PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES));
277 }
278 PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY));
279 PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY));
280 PetscFunctionReturn(PETSC_SUCCESS);
281 }
282
StokesSetupMatBlock01(Stokes * s)283 PetscErrorCode StokesSetupMatBlock01(Stokes *s)
284 {
285 PetscInt row, start, end, sz, i, j;
286 PetscInt cols[5];
287 PetscScalar vals[5];
288
289 PetscFunctionBeginUser;
290 /* A[1] is 2N-by-N */
291 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1]));
292 PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_"));
293 PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny));
294 PetscCall(MatSetType(s->subA[1], MATMPIAIJ));
295 PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL));
296 PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end));
297
298 PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
299
300 for (row = start; row < end; row++) {
301 PetscCall(StokesGetPosition(s, row, &i, &j));
302 /* first part: rows 0 to (nx*ny-1) */
303 if (row < s->nx * s->ny) {
304 PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals));
305 } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
306 PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals));
307 }
308 PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES));
309 }
310 PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY));
311 PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY));
312 PetscFunctionReturn(PETSC_SUCCESS);
313 }
314
StokesSetupMatBlock10(Stokes * s)315 PetscErrorCode StokesSetupMatBlock10(Stokes *s)
316 {
317 PetscFunctionBeginUser;
318 /* A[2] is minus transpose of A[1] */
319 PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]));
320 if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0));
321 PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_"));
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
StokesSetupMatBlock11(Stokes * s)325 PetscErrorCode StokesSetupMatBlock11(Stokes *s)
326 {
327 PetscFunctionBeginUser;
328 /* A[3] is N-by-N null matrix */
329 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3]));
330 PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_"));
331 PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny));
332 PetscCall(MatSetType(s->subA[3], MATMPIAIJ));
333 PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL));
334 PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY));
335 PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY));
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338
StokesSetupApproxSchur(Stokes * s)339 PetscErrorCode StokesSetupApproxSchur(Stokes *s)
340 {
341 Vec diag;
342
343 PetscFunctionBeginUser;
344 /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
345 /* note: A11 is zero */
346 /* note: in real life this matrix would be build directly, */
347 /* i.e. without MatMatMult */
348
349 /* inverse of diagonal of A00 */
350 PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
351 PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny));
352 PetscCall(VecSetType(diag, VECMPI));
353 PetscCall(MatGetDiagonal(s->subA[0], diag));
354 PetscCall(VecReciprocal(diag));
355
356 /* compute: - A10 inv(DIAGFORM(A00)) A01 */
357 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */
358 PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &s->myS));
359 PetscCall(MatScale(s->myS, -1.0));
360
361 /* restore A10 */
362 PetscCall(MatGetDiagonal(s->subA[0], diag));
363 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL));
364 PetscCall(VecDestroy(&diag));
365 PetscFunctionReturn(PETSC_SUCCESS);
366 }
367
StokesSetupMatrix(Stokes * s)368 PetscErrorCode StokesSetupMatrix(Stokes *s)
369 {
370 PetscFunctionBeginUser;
371 PetscCall(StokesSetupMatBlock00(s));
372 PetscCall(StokesSetupMatBlock01(s));
373 PetscCall(StokesSetupMatBlock10(s));
374 PetscCall(StokesSetupMatBlock11(s));
375 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A));
376 PetscCall(StokesSetupApproxSchur(s));
377 PetscFunctionReturn(PETSC_SUCCESS);
378 }
379
StokesStencilLaplacian(Stokes * s,PetscInt i,PetscInt j,PetscInt * sz,PetscInt * cols,PetscScalar * vals)380 PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
381 {
382 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx;
383 PetscScalar ae = s->hy / s->hx, aeb = 0;
384 PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2);
385 PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2);
386 PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2);
387
388 PetscFunctionBeginUser;
389 if (i == 0 && j == 0) { /* south-west corner */
390 *sz = 3;
391 cols[0] = p;
392 vals[0] = -(ae + awb + asb + an);
393 cols[1] = e;
394 vals[1] = ae;
395 cols[2] = n;
396 vals[2] = an;
397 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
398 *sz = 3;
399 cols[0] = s2;
400 vals[0] = as;
401 cols[1] = p;
402 vals[1] = -(ae + awb + as + anb);
403 cols[2] = e;
404 vals[2] = ae;
405 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
406 *sz = 3;
407 cols[0] = w;
408 vals[0] = aw;
409 cols[1] = p;
410 vals[1] = -(aeb + aw + asb + an);
411 cols[2] = n;
412 vals[2] = an;
413 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
414 *sz = 3;
415 cols[0] = s2;
416 vals[0] = as;
417 cols[1] = w;
418 vals[1] = aw;
419 cols[2] = p;
420 vals[2] = -(aeb + aw + as + anb);
421 } else if (i == 0) { /* west boundary */
422 *sz = 4;
423 cols[0] = s2;
424 vals[0] = as;
425 cols[1] = p;
426 vals[1] = -(ae + awb + as + an);
427 cols[2] = e;
428 vals[2] = ae;
429 cols[3] = n;
430 vals[3] = an;
431 } else if (i == s->nx - 1) { /* east boundary */
432 *sz = 4;
433 cols[0] = s2;
434 vals[0] = as;
435 cols[1] = w;
436 vals[1] = aw;
437 cols[2] = p;
438 vals[2] = -(aeb + aw + as + an);
439 cols[3] = n;
440 vals[3] = an;
441 } else if (j == 0) { /* south boundary */
442 *sz = 4;
443 cols[0] = w;
444 vals[0] = aw;
445 cols[1] = p;
446 vals[1] = -(ae + aw + asb + an);
447 cols[2] = e;
448 vals[2] = ae;
449 cols[3] = n;
450 vals[3] = an;
451 } else if (j == s->ny - 1) { /* north boundary */
452 *sz = 4;
453 cols[0] = s2;
454 vals[0] = as;
455 cols[1] = w;
456 vals[1] = aw;
457 cols[2] = p;
458 vals[2] = -(ae + aw + as + anb);
459 cols[3] = e;
460 vals[3] = ae;
461 } else { /* interior */
462 *sz = 5;
463 cols[0] = s2;
464 vals[0] = as;
465 cols[1] = w;
466 vals[1] = aw;
467 cols[2] = p;
468 vals[2] = -(ae + aw + as + an);
469 cols[3] = e;
470 vals[3] = ae;
471 cols[4] = n;
472 vals[4] = an;
473 }
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
StokesStencilGradientX(Stokes * s,PetscInt i,PetscInt j,PetscInt * sz,PetscInt * cols,PetscScalar * vals)477 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
478 {
479 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1;
480 PetscScalar ae = s->hy / 2, aeb = s->hy;
481 PetscScalar aw = -s->hy / 2, awb = 0;
482
483 PetscFunctionBeginUser;
484 if (i == 0 && j == 0) { /* south-west corner */
485 *sz = 2;
486 cols[0] = p;
487 vals[0] = -(ae + awb);
488 cols[1] = e;
489 vals[1] = ae;
490 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
491 *sz = 2;
492 cols[0] = p;
493 vals[0] = -(ae + awb);
494 cols[1] = e;
495 vals[1] = ae;
496 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
497 *sz = 2;
498 cols[0] = w;
499 vals[0] = aw;
500 cols[1] = p;
501 vals[1] = -(aeb + aw);
502 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
503 *sz = 2;
504 cols[0] = w;
505 vals[0] = aw;
506 cols[1] = p;
507 vals[1] = -(aeb + aw);
508 } else if (i == 0) { /* west boundary */
509 *sz = 2;
510 cols[0] = p;
511 vals[0] = -(ae + awb);
512 cols[1] = e;
513 vals[1] = ae;
514 } else if (i == s->nx - 1) { /* east boundary */
515 *sz = 2;
516 cols[0] = w;
517 vals[0] = aw;
518 cols[1] = p;
519 vals[1] = -(aeb + aw);
520 } else if (j == 0) { /* south boundary */
521 *sz = 3;
522 cols[0] = w;
523 vals[0] = aw;
524 cols[1] = p;
525 vals[1] = -(ae + aw);
526 cols[2] = e;
527 vals[2] = ae;
528 } else if (j == s->ny - 1) { /* north boundary */
529 *sz = 3;
530 cols[0] = w;
531 vals[0] = aw;
532 cols[1] = p;
533 vals[1] = -(ae + aw);
534 cols[2] = e;
535 vals[2] = ae;
536 } else { /* interior */
537 *sz = 3;
538 cols[0] = w;
539 vals[0] = aw;
540 cols[1] = p;
541 vals[1] = -(ae + aw);
542 cols[2] = e;
543 vals[2] = ae;
544 }
545 PetscFunctionReturn(PETSC_SUCCESS);
546 }
547
StokesStencilGradientY(Stokes * s,PetscInt i,PetscInt j,PetscInt * sz,PetscInt * cols,PetscScalar * vals)548 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
549 {
550 PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx;
551 PetscScalar as = -s->hx / 2, asb = 0;
552 PetscScalar an = s->hx / 2, anb = 0;
553
554 PetscFunctionBeginUser;
555 if (i == 0 && j == 0) { /* south-west corner */
556 *sz = 2;
557 cols[0] = p;
558 vals[0] = -(asb + an);
559 cols[1] = n;
560 vals[1] = an;
561 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
562 *sz = 2;
563 cols[0] = s2;
564 vals[0] = as;
565 cols[1] = p;
566 vals[1] = -(as + anb);
567 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
568 *sz = 2;
569 cols[0] = p;
570 vals[0] = -(asb + an);
571 cols[1] = n;
572 vals[1] = an;
573 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
574 *sz = 2;
575 cols[0] = s2;
576 vals[0] = as;
577 cols[1] = p;
578 vals[1] = -(as + anb);
579 } else if (i == 0) { /* west boundary */
580 *sz = 3;
581 cols[0] = s2;
582 vals[0] = as;
583 cols[1] = p;
584 vals[1] = -(as + an);
585 cols[2] = n;
586 vals[2] = an;
587 } else if (i == s->nx - 1) { /* east boundary */
588 *sz = 3;
589 cols[0] = s2;
590 vals[0] = as;
591 cols[1] = p;
592 vals[1] = -(as + an);
593 cols[2] = n;
594 vals[2] = an;
595 } else if (j == 0) { /* south boundary */
596 *sz = 2;
597 cols[0] = p;
598 vals[0] = -(asb + an);
599 cols[1] = n;
600 vals[1] = an;
601 } else if (j == s->ny - 1) { /* north boundary */
602 *sz = 2;
603 cols[0] = s2;
604 vals[0] = as;
605 cols[1] = p;
606 vals[1] = -(as + anb);
607 } else { /* interior */
608 *sz = 3;
609 cols[0] = s2;
610 vals[0] = as;
611 cols[1] = p;
612 vals[1] = -(as + an);
613 cols[2] = n;
614 vals[2] = an;
615 }
616 PetscFunctionReturn(PETSC_SUCCESS);
617 }
618
StokesRhsMomX(Stokes * s,PetscInt i,PetscInt j,PetscScalar * val)619 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
620 {
621 PetscScalar y = j * s->hy + s->hy / 2;
622 PetscScalar awb = s->hy / (s->hx / 2);
623
624 PetscFunctionBeginUser;
625 if (i == 0) { /* west boundary */
626 *val = awb * StokesExactVelocityX(y);
627 } else {
628 *val = 0.0;
629 }
630 PetscFunctionReturn(PETSC_SUCCESS);
631 }
632
StokesRhsMomY(Stokes * s,PetscInt i,PetscInt j,PetscScalar * val)633 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
634 {
635 PetscFunctionBeginUser;
636 *val = 0.0;
637 PetscFunctionReturn(PETSC_SUCCESS);
638 }
639
StokesRhsMass(Stokes * s,PetscInt i,PetscInt j,PetscScalar * val)640 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
641 {
642 PetscScalar y = j * s->hy + s->hy / 2;
643 PetscScalar aeb = s->hy;
644
645 PetscFunctionBeginUser;
646 if (i == 0) { /* west boundary */
647 *val = aeb * StokesExactVelocityX(y);
648 } else {
649 *val = 0.0;
650 }
651 PetscFunctionReturn(PETSC_SUCCESS);
652 }
653
StokesCalcResidual(Stokes * s)654 PetscErrorCode StokesCalcResidual(Stokes *s)
655 {
656 PetscReal val;
657 Vec b0, b1;
658
659 PetscFunctionBeginUser;
660 /* residual Ax-b (*warning* overwrites b) */
661 PetscCall(VecScale(s->b, -1.0));
662 PetscCall(MatMultAdd(s->A, s->x, s->b, s->b));
663
664 /* residual velocity */
665 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
666 PetscCall(VecNorm(b0, NORM_2, &val));
667 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val));
668 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
669
670 /* residual pressure */
671 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
672 PetscCall(VecNorm(b1, NORM_2, &val));
673 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val));
674 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
675
676 /* total residual */
677 PetscCall(VecNorm(s->b, NORM_2, &val));
678 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val));
679 PetscFunctionReturn(PETSC_SUCCESS);
680 }
681
StokesCalcError(Stokes * s)682 PetscErrorCode StokesCalcError(Stokes *s)
683 {
684 PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny);
685 PetscReal val;
686 Vec y0, y1;
687
688 PetscFunctionBeginUser;
689 /* error y-x */
690 PetscCall(VecAXPY(s->y, -1.0, s->x));
691
692 /* error in velocity */
693 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
694 PetscCall(VecNorm(y0, NORM_2, &val));
695 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale))));
696 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
697
698 /* error in pressure */
699 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
700 PetscCall(VecNorm(y1, NORM_2, &val));
701 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale))));
702 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
703
704 /* total error */
705 PetscCall(VecNorm(s->y, NORM_2, &val));
706 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart(val / scale)));
707 PetscFunctionReturn(PETSC_SUCCESS);
708 }
709
main(int argc,char ** argv)710 int main(int argc, char **argv)
711 {
712 Stokes s;
713 KSP ksp;
714
715 PetscFunctionBeginUser;
716 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
717 s.nx = 4;
718 s.ny = 6;
719 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL));
720 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL));
721 s.hx = 2.0 / s.nx;
722 s.hy = 1.0 / s.ny;
723 s.matsymmetric = PETSC_FALSE;
724 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL));
725 s.userPC = s.userKSP = PETSC_FALSE;
726 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC));
727 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP));
728
729 PetscCall(StokesSetupMatrix(&s));
730 PetscCall(StokesSetupIndexSets(&s));
731 PetscCall(StokesSetupVectors(&s));
732
733 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
734 PetscCall(KSPSetOperators(ksp, s.A, s.A));
735 PetscCall(KSPSetFromOptions(ksp));
736 PetscCall(StokesSetupPC(&s, ksp));
737 PetscCall(KSPSolve(ksp, s.b, s.x));
738
739 /* don't trust, verify! */
740 PetscCall(StokesCalcResidual(&s));
741 PetscCall(StokesCalcError(&s));
742 PetscCall(StokesWriteSolution(&s));
743
744 PetscCall(KSPDestroy(&ksp));
745 PetscCall(MatDestroy(&s.subA[0]));
746 PetscCall(MatDestroy(&s.subA[1]));
747 PetscCall(MatDestroy(&s.subA[2]));
748 PetscCall(MatDestroy(&s.subA[3]));
749 PetscCall(MatDestroy(&s.A));
750 PetscCall(VecDestroy(&s.x));
751 PetscCall(VecDestroy(&s.b));
752 PetscCall(VecDestroy(&s.y));
753 PetscCall(MatDestroy(&s.myS));
754 PetscCall(PetscFinalize());
755 return 0;
756 }
757
758 /*TEST
759
760 test:
761 nsize: 2
762 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none
763
764 test:
765 suffix: 2
766 nsize: 2
767 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
768
769 test:
770 suffix: 3
771 nsize: 2
772 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
773
774 test:
775 suffix: 4
776 nsize: 2
777 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi
778
779 test:
780 suffix: 4_pcksp
781 nsize: 2
782 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi
783
784 test:
785 suffix: 5
786 nsize: 2
787 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10
788
789 test:
790 suffix: 6
791 nsize: 2
792 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10
793
794 test:
795 suffix: 7
796 nsize: 2
797 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -pc_fieldsplit_gkb_tol 1e-4 -pc_fieldsplit_gkb_nu 5 -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-6
798
799 TEST*/
800