Actual source code: ex43.c

```  1: static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 2d on the unit domain \n\
2: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
3: The models defined utilise free slip boundary conditions on all sides. \n\
4: Options: \n"
5:                      "\
6:      -mx : Number of elements in the x-direction \n\
7:      -my : Number of elements in the y-direction \n\
8:      -o : Specify output filename for solution (will be petsc binary format or paraview format if the extension is .vts) \n\
9:      -gnuplot : Output Gauss point coordinates, coefficients and u,p solution in gnuplot format \n\
10:      -glvis : Visualizes coefficients and u,p solution through GLVIs (use -viewer_glvis_dmda_bs 2,1 to visualize velocity as a vector)\n\
11:      -c_str : Indicates the structure of the coefficients to use \n"
12:                      "\
13:           -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
14:                       This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
15:                          Parameters: \n\
16:                               -solcx_eta0  : Viscosity to the left of the interface \n\
17:                               -solcx_eta1  : Viscosity to the right of the interface \n\
18:                               -solcx_xc    : Location of the interface \n\
19:                               -solcx_nz    : Wavenumber in the y direction \n"
20:                      "\
21:           -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
22:                          Parameters: \n\
23:                               -sinker_eta0 : Viscosity of the background fluid \n\
24:                               -sinker_eta1 : Viscosity of the blob \n\
25:                               -sinker_dx   : Width of the blob \n\
26:                               -sinker_dy   : Height of the blob \n"
27:                      "\
28:           -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
29:                          Parameters: \n\
30:                               -sinker_eta0 : Viscosity of the background fluid \n\
31:                               -sinker_eta1 : Viscosity of the blob \n\
32:                               -sinker_r    : Radius of the blob \n"
33:                      "\
34:           -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
35:                               -sinker_eta0 : Viscosity of the background fluid \n\
36:                               -sinker_eta1 : Viscosity of the two inclusions \n\
37:                               -sinker_r    : Radius of the circular inclusion \n\
38:                               -sinker_c0x  : Origin (x-coord) of the circular inclusion \n\
39:                               -sinker_c0y  : Origin (y-coord) of the circular inclusion \n\
40:                               -sinker_dx   : Width of the rectangular inclusion \n\
41:                               -sinker_dy   : Height of the rectangular inclusion \n\
42:                               -sinker_phi  : Rotation angle of the rectangular inclusion \n"
43:                      "\
44:           -c_str 4 => Coefficient definition for checkerboard jumps aligned with the domain decomposition \n\
45:                               -jump_eta0      : Viscosity for black subdomains \n\
46:                               -jump_magnitude : Magnitude of jumps. White subdomains will have eta = eta0*10^magnitude \n\
47:                               -jump_nz        : Wavenumber in the y direction for rhs \n"
48:                      "\
49:      -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
50:                       By default, the viscosity and force term are evaluated at the element center and applied as a constant over the entire element \n";

52: /* Contributed by Dave May */

54: #include <petscksp.h>
55: #include <petscdm.h>
56: #include <petscdmda.h>

58: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
59: #include "ex43-solcx.h"

61: static PetscErrorCode DMDABCApplyFreeSlip(DM, Mat, Vec);

63: #define NSD          2 /* number of spatial dimensions */
64: #define NODES_PER_EL 4 /* nodes per element */
65: #define U_DOFS       2 /* degrees of freedom per velocity node */
66: #define P_DOFS       1 /* degrees of freedom per pressure node */
67: #define GAUSS_POINTS 4

69: /* Gauss point based evaluation 8+4+4+4 = 20 */
70: typedef struct {
71:   PetscScalar gp_coords[2 * GAUSS_POINTS];
72:   PetscScalar eta[GAUSS_POINTS];
73:   PetscScalar fx[GAUSS_POINTS];
74:   PetscScalar fy[GAUSS_POINTS];
75: } GaussPointCoefficients;

77: typedef struct {
78:   PetscScalar u_dof;
79:   PetscScalar v_dof;
80:   PetscScalar p_dof;
81: } StokesDOF;

83: static PetscErrorCode glvis_extract_eta(PetscObject oV, PetscInt nf, PetscObject oVf[], void *ctx)
84: {
85:   DM                       properties_da = (DM)(ctx), stokes_da;
86:   Vec                      V = (Vec)oV, *Vf = (Vec *)oVf;
87:   GaussPointCoefficients **props;
88:   PetscInt                 sex, sey, mx, my;
89:   PetscInt                 ei, ej, p, cum;
90:   PetscScalar             *array;

92:   PetscFunctionBeginUser;
93:   PetscCall(VecGetDM(Vf[0], &stokes_da));
95:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
96:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
97:   PetscCall(VecGetArray(Vf[0], &array));
98:   cum = 0;
99:   for (ej = sey; ej < sey + my; ej++) {
100:     for (ei = sex; ei < sex + mx; ei++) {
101:       for (p = 0; p < GAUSS_POINTS; p++) array[cum++] = props[ej][ei].eta[p];
102:     }
103:   }
104:   PetscCall(VecRestoreArray(Vf[0], &array));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: /*
110:  Element: Local basis function ordering
111:  1-----2
112:  |     |
113:  |     |
114:  0-----3
115:  */
116: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
117: {
118:   PetscScalar xi  = _xi[0];
119:   PetscScalar eta = _xi[1];

121:   Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
122:   Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
123:   Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
124:   Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
125: }

127: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
128: {
129:   PetscScalar xi  = _xi[0];
130:   PetscScalar eta = _xi[1];

132:   GNi[0][0] = -0.25 * (1.0 - eta);
133:   GNi[0][1] = -0.25 * (1.0 + eta);
134:   GNi[0][2] = 0.25 * (1.0 + eta);
135:   GNi[0][3] = 0.25 * (1.0 - eta);

137:   GNi[1][0] = -0.25 * (1.0 - xi);
138:   GNi[1][1] = 0.25 * (1.0 - xi);
139:   GNi[1][2] = 0.25 * (1.0 + xi);
140:   GNi[1][3] = -0.25 * (1.0 + xi);
141: }

143: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
144: {
145:   PetscScalar J00, J01, J10, J11, J;
146:   PetscScalar iJ00, iJ01, iJ10, iJ11;
147:   PetscInt    i;

149:   J00 = J01 = J10 = J11 = 0.0;
150:   for (i = 0; i < NODES_PER_EL; i++) {
151:     PetscScalar cx = coords[2 * i];
152:     PetscScalar cy = coords[2 * i + 1];

154:     J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
155:     J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
156:     J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
157:     J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
158:   }
159:   J = (J00 * J11) - (J01 * J10);

161:   iJ00 = J11 / J;
162:   iJ01 = -J01 / J;
163:   iJ10 = -J10 / J;
164:   iJ11 = J00 / J;

166:   for (i = 0; i < NODES_PER_EL; i++) {
167:     GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
168:     GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
169:   }

171:   if (det_J) *det_J = J;
172: }

174: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
175: {
176:   *ngp         = 4;
177:   gp_xi[0][0]  = -0.57735026919;
178:   gp_xi[0][1]  = -0.57735026919;
179:   gp_xi[1][0]  = -0.57735026919;
180:   gp_xi[1][1]  = 0.57735026919;
181:   gp_xi[2][0]  = 0.57735026919;
182:   gp_xi[2][1]  = 0.57735026919;
183:   gp_xi[3][0]  = 0.57735026919;
184:   gp_xi[3][1]  = -0.57735026919;
185:   gp_weight[0] = 1.0;
186:   gp_weight[1] = 1.0;
187:   gp_weight[2] = 1.0;
188:   gp_weight[3] = 1.0;
189: }

191: /*
192: i,j are the element indices
193: The unknown is a vector quantity.
194: The s[].c is used to indicate the degree of freedom.
195: */
196: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j)
197: {
198:   PetscFunctionBeginUser;
199:   /* velocity */
200:   /* node 0 */
201:   s_u[0].i = i;
202:   s_u[0].j = j;
203:   s_u[0].c = 0; /* Vx0 */
204:   s_u[1].i = i;
205:   s_u[1].j = j;
206:   s_u[1].c = 1; /* Vy0 */

208:   /* node 1 */
209:   s_u[2].i = i;
210:   s_u[2].j = j + 1;
211:   s_u[2].c = 0; /* Vx1 */
212:   s_u[3].i = i;
213:   s_u[3].j = j + 1;
214:   s_u[3].c = 1; /* Vy1 */

216:   /* node 2 */
217:   s_u[4].i = i + 1;
218:   s_u[4].j = j + 1;
219:   s_u[4].c = 0; /* Vx2 */
220:   s_u[5].i = i + 1;
221:   s_u[5].j = j + 1;
222:   s_u[5].c = 1; /* Vy2 */

224:   /* node 3 */
225:   s_u[6].i = i + 1;
226:   s_u[6].j = j;
227:   s_u[6].c = 0; /* Vx3 */
228:   s_u[7].i = i + 1;
229:   s_u[7].j = j;
230:   s_u[7].c = 1; /* Vy3 */

232:   /* pressure */
233:   s_p[0].i = i;
234:   s_p[0].j = j;
235:   s_p[0].c = 2; /* P0 */
236:   s_p[1].i = i;
237:   s_p[1].j = j + 1;
238:   s_p[1].c = 2; /* P1 */
239:   s_p[2].i = i + 1;
240:   s_p[2].j = j + 1;
241:   s_p[2].c = 2; /* P2 */
242:   s_p[3].i = i + 1;
243:   s_p[3].j = j;
244:   s_p[3].c = 2; /* P3 */
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
249: {
250:   PetscMPIInt  rank;
251:   PetscInt     proc_I, proc_J;
252:   PetscInt     cpu_x, cpu_y;
253:   PetscInt     local_mx, local_my;
254:   Vec          vlx, vly;
255:   PetscInt    *LX, *LY, i;
256:   PetscScalar *_a;
257:   Vec          V_SEQ;
258:   VecScatter   ctx;

260:   PetscFunctionBeginUser;
261:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

263:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));

265:   proc_J = rank / cpu_x;
266:   proc_I = rank - cpu_x * proc_J;

268:   PetscCall(PetscMalloc1(cpu_x, &LX));
269:   PetscCall(PetscMalloc1(cpu_y, &LY));

271:   PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
272:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
273:   PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
274:   PetscCall(VecSetFromOptions(vlx));

276:   PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
277:   PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
278:   PetscCall(VecSetFromOptions(vly));

280:   PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
281:   PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
282:   PetscCall(VecAssemblyBegin(vlx));
283:   PetscCall(VecAssemblyEnd(vlx));
284:   PetscCall(VecAssemblyBegin(vly));
285:   PetscCall(VecAssemblyEnd(vly));

287:   PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
288:   PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
289:   PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
290:   PetscCall(VecGetArray(V_SEQ, &_a));
291:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
292:   PetscCall(VecRestoreArray(V_SEQ, &_a));
293:   PetscCall(VecScatterDestroy(&ctx));
294:   PetscCall(VecDestroy(&V_SEQ));

296:   PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
297:   PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
298:   PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
299:   PetscCall(VecGetArray(V_SEQ, &_a));
300:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
301:   PetscCall(VecRestoreArray(V_SEQ, &_a));
302:   PetscCall(VecScatterDestroy(&ctx));
303:   PetscCall(VecDestroy(&V_SEQ));

305:   *_lx = LX;
306:   *_ly = LY;

308:   PetscCall(VecDestroy(&vlx));
309:   PetscCall(VecDestroy(&vly));
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
314: {
315:   DM           cda;
316:   Vec          coords;
317:   DMDACoor2d **_coords;
318:   PetscInt     si, sj, nx, ny, i, j;
319:   FILE        *fp;
320:   char         fname[PETSC_MAX_PATH_LEN];
321:   PetscMPIInt  rank;

323:   PetscFunctionBeginUser;
324:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
325:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
326:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
327:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

329:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));

331:   PetscCall(DMGetCoordinateDM(da, &cda));
332:   PetscCall(DMGetCoordinatesLocal(da, &coords));
334:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
335:   for (j = sj; j < sj + ny - 1; j++) {
336:     for (i = si; i < si + nx - 1; i++) {
337:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
338:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
339:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y)));
340:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
341:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
342:     }
343:   }

346:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
351: {
352:   DM           cda;
353:   Vec          coords, local_fields;
354:   DMDACoor2d **_coords;
355:   FILE        *fp;
356:   char         fname[PETSC_MAX_PATH_LEN];
357:   PetscMPIInt  rank;
358:   PetscInt     si, sj, nx, ny, i, j;
359:   PetscInt     n_dofs, d;
360:   PetscScalar *_fields;

362:   PetscFunctionBeginUser;
363:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
364:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
365:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
366:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

368:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
369:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
370:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
371:   for (d = 0; d < n_dofs; d++) {
372:     const char *field_name;
373:     PetscCall(DMDAGetFieldName(da, d, &field_name));
374:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
375:   }
376:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

378:   PetscCall(DMGetCoordinateDM(da, &cda));
379:   PetscCall(DMGetCoordinatesLocal(da, &coords));
380:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
381:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

383:   PetscCall(DMCreateLocalVector(da, &local_fields));
384:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
385:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
386:   PetscCall(VecGetArray(local_fields, &_fields));

388:   for (j = sj; j < sj + ny; j++) {
389:     for (i = si; i < si + nx; i++) {
390:       PetscScalar coord_x, coord_y;
391:       PetscScalar field_d;

393:       coord_x = _coords[j][i].x;
394:       coord_y = _coords[j][i].y;

396:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
397:       for (d = 0; d < n_dofs; d++) {
398:         field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
399:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
400:       }
401:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
402:     }
403:   }
404:   PetscCall(VecRestoreArray(local_fields, &_fields));
405:   PetscCall(VecDestroy(&local_fields));

407:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

409:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
414: {
415:   DM                       cda;
416:   Vec                      local_fields;
417:   FILE                    *fp;
418:   char                     fname[PETSC_MAX_PATH_LEN];
419:   PetscMPIInt              rank;
420:   PetscInt                 si, sj, nx, ny, i, j, p;
421:   PetscInt                 n_dofs, d;
422:   GaussPointCoefficients **_coefficients;

424:   PetscFunctionBeginUser;
425:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
426:   PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
427:   PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
428:   PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");

430:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
431:   PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
432:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
433:   for (d = 0; d < n_dofs; d++) {
434:     const char *field_name;
435:     PetscCall(DMDAGetFieldName(da, d, &field_name));
436:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
437:   }
438:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));

440:   PetscCall(DMGetCoordinateDM(da, &cda));
441:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));

443:   PetscCall(DMCreateLocalVector(da, &local_fields));
444:   PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
445:   PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
446:   PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));

448:   for (j = sj; j < sj + ny; j++) {
449:     for (i = si; i < si + nx; i++) {
450:       PetscScalar coord_x, coord_y;

452:       for (p = 0; p < GAUSS_POINTS; p++) {
453:         coord_x = _coefficients[j][i].gp_coords[2 * p];
454:         coord_y = _coefficients[j][i].gp_coords[2 * p + 1];

456:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));

458:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e", (double)PetscRealPart(_coefficients[j][i].eta[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]), (double)PetscRealPart(_coefficients[j][i].fy[p])));
459:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
460:       }
461:     }
462:   }
463:   PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
464:   PetscCall(VecDestroy(&local_fields));

466:   PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi, PetscInt wd, PetscInt w_NPE, PetscInt w_dof, PetscInt ui, PetscInt ud, PetscInt u_NPE, PetscInt u_dof)
471: {
472:   PetscInt ij;
473:   PetscInt r, c, nc;

475:   nc = u_NPE * u_dof;
476:   r  = w_dof * wi + wd;
477:   c  = u_dof * ui + ud;
478:   ij = r * nc + c;
479:   return ij;
480: }

482: /*
483:  D = [ 2.eta   0   0   ]
484:      [   0   2.eta 0   ]
485:      [   0     0   eta ]

487:  B = [ d_dx   0   ]
488:      [  0    d_dy ]
489:      [ d_dy  d_dx ]
490: */
491: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
492: {
493:   PetscInt    ngp;
494:   PetscScalar gp_xi[GAUSS_POINTS][2];
495:   PetscScalar gp_weight[GAUSS_POINTS];
496:   PetscInt    p, i, j, k;
497:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
498:   PetscScalar J_p, tildeD[3];
499:   PetscScalar B[3][U_DOFS * NODES_PER_EL];

501:   /* define quadrature rule */

504:   /* evaluate integral */
505:   for (p = 0; p < ngp; p++) {
506:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
507:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);

509:     for (i = 0; i < NODES_PER_EL; i++) {
510:       PetscScalar d_dx_i = GNx_p[0][i];
511:       PetscScalar d_dy_i = GNx_p[1][i];

513:       B[0][2 * i]     = d_dx_i;
514:       B[0][2 * i + 1] = 0.0;
515:       B[1][2 * i]     = 0.0;
516:       B[1][2 * i + 1] = d_dy_i;
517:       B[2][2 * i]     = d_dy_i;
518:       B[2][2 * i + 1] = d_dx_i;
519:     }

521:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
522:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
523:     tildeD[2] = gp_weight[p] * J_p * eta[p];

525:     /* form Bt tildeD B */
526:     /*
527:     Ke_ij = Bt_ik . D_kl . B_lj
528:           = B_ki . D_kl . B_lj
529:           = B_ki . D_kk . B_kj
530:     */
531:     for (i = 0; i < 8; i++) {
532:       for (j = 0; j < 8; j++) {
533:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
534:           Ke[i + 8 * j] = Ke[i + 8 * j] + B[k][i] * tildeD[k] * B[k][j];
535:         }
536:       }
537:     }
538:   }
539: }

541: static void FormGradientOperatorQ1(PetscScalar Ke[], PetscScalar coords[])
542: {
543:   PetscInt    ngp;
544:   PetscScalar gp_xi[GAUSS_POINTS][2];
545:   PetscScalar gp_weight[GAUSS_POINTS];
546:   PetscInt    p, i, j, di;
547:   PetscScalar Ni_p[NODES_PER_EL];
548:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
549:   PetscScalar J_p, fac;

551:   /* define quadrature rule */

554:   /* evaluate integral */
555:   for (p = 0; p < ngp; p++) {
556:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
557:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
558:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
559:     fac = gp_weight[p] * J_p;

561:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
562:       for (di = 0; di < NSD; di++) {     /* u dofs */
563:         for (j = 0; j < 4; j++) {        /* p nodes, p dofs = 1 (ie no loop) */
564:           PetscInt IJ;
565:           IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);

567:           Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
568:         }
569:       }
570:     }
571:   }
572: }

574: static void FormDivergenceOperatorQ1(PetscScalar De[], PetscScalar coords[])
575: {
576:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
577:   PetscInt    i, j;
578:   PetscInt    nr_g, nc_g;

580:   PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));

583:   nr_g = U_DOFS * NODES_PER_EL;
584:   nc_g = P_DOFS * NODES_PER_EL;

586:   for (i = 0; i < nr_g; i++) {
587:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
588:   }
589: }

591: static void FormStabilisationOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
592: {
593:   PetscInt    ngp;
594:   PetscScalar gp_xi[GAUSS_POINTS][2];
595:   PetscScalar gp_weight[GAUSS_POINTS];
596:   PetscInt    p, i, j;
597:   PetscScalar Ni_p[NODES_PER_EL];
598:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
599:   PetscScalar J_p, fac, eta_avg;

601:   /* define quadrature rule */

604:   /* evaluate integral */
605:   for (p = 0; p < ngp; p++) {
606:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
607:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
608:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
609:     fac = gp_weight[p] * J_p;

611:     for (i = 0; i < NODES_PER_EL; i++) {
612:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * (Ni_p[i] * Ni_p[j] - 0.0625);
613:     }
614:   }

616:   /* scale */
617:   eta_avg = 0.0;
618:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
619:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
620:   fac     = 1.0 / eta_avg;
621:   for (i = 0; i < NODES_PER_EL; i++) {
622:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
623:   }
624: }

626: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
627: {
628:   PetscInt    ngp;
629:   PetscScalar gp_xi[GAUSS_POINTS][2];
630:   PetscScalar gp_weight[GAUSS_POINTS];
631:   PetscInt    p, i, j;
632:   PetscScalar Ni_p[NODES_PER_EL];
633:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
634:   PetscScalar J_p, fac, eta_avg;

636:   /* define quadrature rule */

639:   /* evaluate integral */
640:   for (p = 0; p < ngp; p++) {
641:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
642:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
643:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
644:     fac = gp_weight[p] * J_p;

646:     for (i = 0; i < NODES_PER_EL; i++) {
647:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * Ni_p[i] * Ni_p[j];
648:     }
649:   }

651:   /* scale */
652:   eta_avg = 0.0;
653:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
654:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
655:   fac     = 1.0 / eta_avg;
656:   for (i = 0; i < NODES_PER_EL; i++) {
657:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
658:   }
659: }

661: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
662: {
663:   PetscInt    ngp;
664:   PetscScalar gp_xi[GAUSS_POINTS][2];
665:   PetscScalar gp_weight[GAUSS_POINTS];
666:   PetscInt    p, i;
667:   PetscScalar Ni_p[NODES_PER_EL];
668:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
669:   PetscScalar J_p, fac;

671:   /* define quadrature rule */

674:   /* evaluate integral */
675:   for (p = 0; p < ngp; p++) {
676:     ConstructQ12D_Ni(gp_xi[p], Ni_p);
677:     ConstructQ12D_GNi(gp_xi[p], GNi_p);
678:     ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
679:     fac = gp_weight[p] * J_p;

681:     for (i = 0; i < NODES_PER_EL; i++) {
682:       Fe[NSD * i] += fac * Ni_p[i] * fx[p];
683:       Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
684:     }
685:   }
686: }

688: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
689: {
690:   PetscFunctionBeginUser;
691:   /* get coords for the element */
692:   el_coords[NSD * 0]     = _coords[ej][ei].x;
693:   el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
694:   el_coords[NSD * 1]     = _coords[ej + 1][ei].x;
695:   el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
696:   el_coords[NSD * 2]     = _coords[ej + 1][ei + 1].x;
697:   el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
698:   el_coords[NSD * 3]     = _coords[ej][ei + 1].x;
699:   el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
700:   PetscFunctionReturn(PETSC_SUCCESS);
701: }

703: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
704: {
705:   DM                       cda;
706:   Vec                      coords;
707:   DMDACoor2d             **_coords;
708:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
709:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
710:   PetscInt                 sex, sey, mx, my;
711:   PetscInt                 ei, ej;
712:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
713:   PetscScalar              Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
714:   PetscScalar              De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
715:   PetscScalar              Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
716:   PetscScalar              el_coords[NODES_PER_EL * NSD];
717:   Vec                      local_properties;
718:   GaussPointCoefficients **props;
719:   PetscScalar             *prop_eta;

721:   PetscFunctionBeginUser;
722:   /* setup for coords */
723:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
724:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
725:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

727:   /* setup for coefficients */
728:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
729:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
730:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
731:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

733:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
734:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
735:   for (ej = sey; ej < sey + my; ej++) {
736:     for (ei = sex; ei < sex + mx; ei++) {
737:       /* get coords for the element */
738:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

740:       /* get coefficients for the element */
741:       prop_eta = props[ej][ei].eta;

743:       /* initialise element stiffness matrix */
744:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));
745:       PetscCall(PetscMemzero(Ge, sizeof(Ge)));
746:       PetscCall(PetscMemzero(De, sizeof(De)));
747:       PetscCall(PetscMemzero(Ce, sizeof(Ce)));

749:       /* form element stiffness matrix */
750:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
752:       FormDivergenceOperatorQ1(De, el_coords);
753:       FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);

755:       /* insert element matrix into global matrix */
756:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
757:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
758:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
759:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
760:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
761:     }
762:   }
763:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
764:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

766:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

768:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
769:   PetscCall(VecDestroy(&local_properties));
770:   PetscFunctionReturn(PETSC_SUCCESS);
771: }

773: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
774: {
775:   DM                       cda;
776:   Vec                      coords;
777:   DMDACoor2d             **_coords;
778:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
779:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
780:   PetscInt                 sex, sey, mx, my;
781:   PetscInt                 ei, ej;
782:   PetscScalar              Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
783:   PetscScalar              Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
784:   PetscScalar              De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
785:   PetscScalar              Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
786:   PetscScalar              el_coords[NODES_PER_EL * NSD];
787:   Vec                      local_properties;
788:   GaussPointCoefficients **props;
789:   PetscScalar             *prop_eta;

791:   PetscFunctionBeginUser;
792:   /* setup for coords */
793:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
794:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
795:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

797:   /* setup for coefficients */
798:   PetscCall(DMCreateLocalVector(properties_da, &local_properties));
799:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
800:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
801:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

803:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
804:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
805:   for (ej = sey; ej < sey + my; ej++) {
806:     for (ei = sex; ei < sex + mx; ei++) {
807:       /* get coords for the element */
808:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

810:       /* get coefficients for the element */
811:       prop_eta = props[ej][ei].eta;

813:       /* initialise element stiffness matrix */
814:       PetscCall(PetscMemzero(Ae, sizeof(Ae)));
815:       PetscCall(PetscMemzero(Ge, sizeof(Ge)));
816:       PetscCall(PetscMemzero(De, sizeof(De)));
817:       PetscCall(PetscMemzero(Ce, sizeof(Ce)));

819:       /* form element stiffness matrix */
820:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
822:       FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);

824:       /* insert element matrix into global matrix */
825:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
826:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
827:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
828:       PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
829:     }
830:   }
831:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
832:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

834:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

836:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
837:   PetscCall(VecDestroy(&local_properties));
838:   PetscFunctionReturn(PETSC_SUCCESS);
839: }

841: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
842: {
843:   PetscInt n;

845:   PetscFunctionBeginUser;
846:   for (n = 0; n < 4; n++) {
847:     fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof         = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof + Fe_u[2 * n];
848:     fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof + Fe_u[2 * n + 1];
849:     fields_F[p_eqn[n].j][p_eqn[n].i].p_dof                 = fields_F[p_eqn[n].j][p_eqn[n].i].p_dof + Fe_p[n];
850:   }
851:   PetscFunctionReturn(PETSC_SUCCESS);
852: }

854: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, DM properties_da, Vec properties)
855: {
856:   DM                       cda;
857:   Vec                      coords;
858:   DMDACoor2d             **_coords;
859:   MatStencil               u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
860:   MatStencil               p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
861:   PetscInt                 sex, sey, mx, my;
862:   PetscInt                 ei, ej;
863:   PetscScalar              Fe[NODES_PER_EL * U_DOFS];
864:   PetscScalar              He[NODES_PER_EL * P_DOFS];
865:   PetscScalar              el_coords[NODES_PER_EL * NSD];
866:   Vec                      local_properties;
867:   GaussPointCoefficients **props;
868:   PetscScalar             *prop_fx, *prop_fy;
869:   Vec                      local_F;
870:   StokesDOF              **ff;

872:   PetscFunctionBeginUser;
873:   /* setup for coords */
874:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
875:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
876:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

878:   /* setup for coefficients */
879:   PetscCall(DMGetLocalVector(properties_da, &local_properties));
880:   PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
881:   PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
882:   PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));

885:   PetscCall(DMGetLocalVector(stokes_da, &local_F));
886:   PetscCall(VecZeroEntries(local_F));
887:   PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));

889:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
890:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
891:   for (ej = sey; ej < sey + my; ej++) {
892:     for (ei = sex; ei < sex + mx; ei++) {
893:       /* get coords for the element */
894:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));

896:       /* get coefficients for the element */
897:       prop_fx = props[ej][ei].fx;
898:       prop_fy = props[ej][ei].fy;

900:       /* initialise element stiffness matrix */
901:       PetscCall(PetscMemzero(Fe, sizeof(Fe)));
902:       PetscCall(PetscMemzero(He, sizeof(He)));

904:       /* form element stiffness matrix */
905:       FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);

907:       /* insert element matrix into global matrix */
908:       PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));

910:       PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
911:     }
912:   }

914:   PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
917:   PetscCall(DMRestoreLocalVector(stokes_da, &local_F));

919:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

921:   PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
922:   PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
923:   PetscFunctionReturn(PETSC_SUCCESS);
924: }

926: static PetscErrorCode DMDACreateSolCx(PetscReal eta0, PetscReal eta1, PetscReal xc, PetscInt nz, PetscInt mx, PetscInt my, DM *_da, Vec *_X)
927: {
928:   DM           da, cda;
929:   Vec          X;
930:   StokesDOF  **_stokes;
931:   Vec          coords;
932:   DMDACoor2d **_coords;
933:   PetscInt     si, sj, ei, ej, i, j;

935:   PetscFunctionBeginUser;
936:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, &da));
937:   PetscCall(DMSetFromOptions(da));
938:   PetscCall(DMSetUp(da));
939:   PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
940:   PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
941:   PetscCall(DMDASetFieldName(da, 2, "analytic_P"));

943:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.));

945:   PetscCall(DMGetCoordinatesLocal(da, &coords));
946:   PetscCall(DMGetCoordinateDM(da, &cda));
947:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

949:   PetscCall(DMCreateGlobalVector(da, &X));
950:   PetscCall(DMDAVecGetArray(da, X, &_stokes));

952:   PetscCall(DMDAGetCorners(da, &si, &sj, 0, &ei, &ej, 0));
953:   for (j = sj; j < sj + ej; j++) {
954:     for (i = si; i < si + ei; i++) {
955:       PetscReal pos[2], pressure, vel[2], total_stress[3], strain_rate[3];

957:       pos[0] = PetscRealPart(_coords[j][i].x);
958:       pos[1] = PetscRealPart(_coords[j][i].y);

960:       evaluate_solCx(pos, eta0, eta1, xc, nz, vel, &pressure, total_stress, strain_rate);

962:       _stokes[j][i].u_dof = vel[0];
963:       _stokes[j][i].v_dof = vel[1];
964:       _stokes[j][i].p_dof = pressure;
965:     }
966:   }
967:   PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
968:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

970:   *_da = da;
971:   *_X  = X;
972:   PetscFunctionReturn(PETSC_SUCCESS);
973: }

975: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
976: {
977:   PetscFunctionBeginUser;
978:   /* get the nodal fields */
979:   nodal_fields[0].u_dof = fields[ej][ei].u_dof;
980:   nodal_fields[0].v_dof = fields[ej][ei].v_dof;
981:   nodal_fields[0].p_dof = fields[ej][ei].p_dof;
982:   nodal_fields[1].u_dof = fields[ej + 1][ei].u_dof;
983:   nodal_fields[1].v_dof = fields[ej + 1][ei].v_dof;
984:   nodal_fields[1].p_dof = fields[ej + 1][ei].p_dof;
985:   nodal_fields[2].u_dof = fields[ej + 1][ei + 1].u_dof;
986:   nodal_fields[2].v_dof = fields[ej + 1][ei + 1].v_dof;
987:   nodal_fields[2].p_dof = fields[ej + 1][ei + 1].p_dof;
988:   nodal_fields[3].u_dof = fields[ej][ei + 1].u_dof;
989:   nodal_fields[3].v_dof = fields[ej][ei + 1].v_dof;
990:   nodal_fields[3].p_dof = fields[ej][ei + 1].p_dof;
991:   PetscFunctionReturn(PETSC_SUCCESS);
992: }

994: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da, Vec X, Vec X_analytic)
995: {
996:   DM           cda;
997:   Vec          coords, X_analytic_local, X_local;
998:   DMDACoor2d **_coords;
999:   PetscInt     sex, sey, mx, my;
1000:   PetscInt     ei, ej;
1001:   PetscScalar  el_coords[NODES_PER_EL * NSD];
1002:   StokesDOF  **stokes_analytic, **stokes;
1003:   StokesDOF    stokes_analytic_e[4], stokes_e[4];

1005:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1006:   PetscScalar Ni_p[NODES_PER_EL];
1007:   PetscInt    ngp;
1008:   PetscScalar gp_xi[GAUSS_POINTS][2];
1009:   PetscScalar gp_weight[GAUSS_POINTS];
1010:   PetscInt    p, i;
1011:   PetscScalar J_p, fac;
1012:   PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1013:   PetscInt    M;
1014:   PetscReal   xymin[2], xymax[2];

1016:   PetscFunctionBeginUser;
1017:   /* define quadrature rule */

1020:   /* setup for coords */
1021:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1022:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1023:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));

1025:   /* setup for analytic */
1026:   PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
1027:   PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1028:   PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1029:   PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));

1031:   /* setup for solution */
1032:   PetscCall(DMCreateLocalVector(stokes_da, &X_local));
1033:   PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
1034:   PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
1035:   PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));

1037:   PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1038:   PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));

1040:   h = (xymax[0] - xymin[0]) / ((PetscReal)M);

1042:   tp_L2 = tu_L2 = tu_H1 = 0.0;

1044:   PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
1045:   PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
1046:   for (ej = sey; ej < sey + my; ej++) {
1047:     for (ei = sex; ei < sex + mx; ei++) {
1048:       /* get coords for the element */
1049:       PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
1050:       PetscCall(StokesDAGetNodalFields(stokes, ei, ej, stokes_e));
1051:       PetscCall(StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e));

1053:       /* evaluate integral */
1054:       p_e_L2 = 0.0;
1055:       u_e_L2 = 0.0;
1056:       u_e_H1 = 0.0;
1057:       for (p = 0; p < ngp; p++) {
1058:         ConstructQ12D_Ni(gp_xi[p], Ni_p);
1059:         ConstructQ12D_GNi(gp_xi[p], GNi_p);
1060:         ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
1061:         fac = gp_weight[p] * J_p;

1063:         for (i = 0; i < NODES_PER_EL; i++) {
1064:           PetscScalar u_error, v_error;

1066:           p_e_L2 = p_e_L2 + fac * Ni_p[i] * (stokes_e[i].p_dof - stokes_analytic_e[i].p_dof) * (stokes_e[i].p_dof - stokes_analytic_e[i].p_dof);

1068:           u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1069:           v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1070:           u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);

1072:           u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error     /* du/dx */
1073:                                    + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error   /* du/dy */
1074:                                    + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error   /* dv/dx */
1075:                                    + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
1076:         }
1077:       }

1079:       tp_L2 += p_e_L2;
1080:       tu_L2 += u_e_L2;
1081:       tu_H1 += u_e_H1;
1082:     }
1083:   }
1084:   PetscCall(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1085:   PetscCall(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1086:   PetscCall(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1087:   p_L2 = PetscSqrtScalar(p_L2);
1088:   u_L2 = PetscSqrtScalar(u_L2);
1089:   u_H1 = PetscSqrtScalar(u_H1);

1091:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%1.4e   %1.4e   %1.4e   %1.4e\n", (double)PetscRealPart(h), (double)PetscRealPart(p_L2), (double)PetscRealPart(u_L2), (double)PetscRealPart(u_H1)));

1093:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));

1095:   PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1096:   PetscCall(VecDestroy(&X_analytic_local));
1097:   PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1098:   PetscCall(VecDestroy(&X_local));
1099:   PetscFunctionReturn(PETSC_SUCCESS);
1100: }

1102: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
1103: {
1104:   DM                       da_Stokes, da_prop;
1105:   PetscInt                 u_dof, p_dof, dof, stencil_width;
1106:   Mat                      A, B;
1107:   DM                       prop_cda, vel_cda;
1108:   Vec                      prop_coords, vel_coords;
1109:   PetscInt                 si, sj, nx, ny, i, j, p;
1110:   Vec                      f, X;
1111:   PetscInt                 prop_dof, prop_stencil_width;
1112:   Vec                      properties, l_properties;
1113:   PetscReal                dx, dy;
1114:   PetscInt                 M, N;
1115:   DMDACoor2d             **_prop_coords, **_vel_coords;
1116:   GaussPointCoefficients **element_props;
1117:   PetscInt                 its;
1118:   KSP                      ksp_S;
1119:   PetscInt                 coefficient_structure = 0;
1120:   PetscInt                 cpu_x, cpu_y, *lx = NULL, *ly = NULL;
1121:   PetscBool                use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
1122:   char                     filename[PETSC_MAX_PATH_LEN];

1124:   PetscFunctionBeginUser;
1125:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL));
1126:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL));
1127:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL));

1129:   /* Generate the da for velocity and pressure */
1130:   /*
1131:   We use Q1 elements for the temperature.
1132:   FEM has a 9-point stencil (BOX) or connectivity pattern
1133:   Num nodes in each direction is mx+1, my+1
1134:   */
1135:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1136:   p_dof         = P_DOFS; /* p - pressure */
1137:   dof           = u_dof + p_dof;
1138:   stencil_width = 1;
1139:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &da_Stokes));

1141:   PetscCall(DMSetMatType(da_Stokes, MATAIJ));
1142:   PetscCall(DMSetFromOptions(da_Stokes));
1143:   PetscCall(DMSetUp(da_Stokes));
1144:   PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
1145:   PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
1146:   PetscCall(DMDASetFieldName(da_Stokes, 2, "P"));

1148:   /* unit box [0,1] x [0,1] */
1149:   PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.));

1151:   /* Generate element properties, we will assume all material properties are constant over the element */
1152:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!!  */
1153:   PetscCall(DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
1154:   PetscCall(DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly));

1156:   prop_dof           = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1157:   prop_stencil_width = 0;
1158:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, cpu_x, cpu_y, prop_dof, prop_stencil_width, lx, ly, &da_prop));
1159:   PetscCall(DMSetFromOptions(da_prop));
1160:   PetscCall(DMSetUp(da_prop));
1161:   PetscCall(PetscFree(lx));
1162:   PetscCall(PetscFree(ly));

1164:   /* define centroid positions */
1165:   PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1166:   dx = 1.0 / ((PetscReal)(M));
1167:   dy = 1.0 / ((PetscReal)(N));

1169:   PetscCall(DMDASetUniformCoordinates(da_prop, 0.0 + 0.5 * dx, 1.0 - 0.5 * dx, 0.0 + 0.5 * dy, 1.0 - 0.5 * dy, 0., 0));

1171:   /* define coefficients */
1172:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));

1174:   PetscCall(DMCreateGlobalVector(da_prop, &properties));
1175:   PetscCall(DMCreateLocalVector(da_prop, &l_properties));
1176:   PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));

1178:   PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
1179:   PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
1180:   PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));

1182:   PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));

1184:   PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
1185:   PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
1186:   PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));

1188:   /* interpolate the coordinates */
1189:   for (j = sj; j < sj + ny; j++) {
1190:     for (i = si; i < si + nx; i++) {
1191:       PetscInt    ngp;
1192:       PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
1193:       PetscScalar el_coords[8];

1195:       PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));

1198:       for (p = 0; p < GAUSS_POINTS; p++) {
1199:         PetscScalar gp_x, gp_y;
1200:         PetscInt    n;
1201:         PetscScalar xi_p[2], Ni_p[4];

1203:         xi_p[0] = gp_xi[p][0];
1204:         xi_p[1] = gp_xi[p][1];
1205:         ConstructQ12D_Ni(xi_p, Ni_p);

1207:         gp_x = 0.0;
1208:         gp_y = 0.0;
1209:         for (n = 0; n < NODES_PER_EL; n++) {
1210:           gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
1211:           gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
1212:         }
1213:         element_props[j][i].gp_coords[2 * p]     = gp_x;
1214:         element_props[j][i].gp_coords[2 * p + 1] = gp_y;
1215:       }
1216:     }
1217:   }

1219:   /* define the coefficients */
1220:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL));

1222:   for (j = sj; j < sj + ny; j++) {
1223:     for (i = si; i < si + nx; i++) {
1224:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1225:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1226:       PetscReal coord_x, coord_y;

1228:       if (coefficient_structure == 0) {
1229:         PetscReal opts_eta0, opts_eta1, opts_xc;
1230:         PetscInt  opts_nz;

1232:         opts_eta0 = 1.0;
1233:         opts_eta1 = 1.0;
1234:         opts_xc   = 0.5;
1235:         opts_nz   = 1;

1237:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1238:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1239:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1240:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));

1242:         for (p = 0; p < GAUSS_POINTS; p++) {
1243:           coord_x = centroid_x;
1244:           coord_y = centroid_y;
1245:           if (use_gp_coords) {
1246:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1247:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1248:           }

1250:           element_props[j][i].eta[p] = opts_eta0;
1251:           if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;

1253:           element_props[j][i].fx[p] = 0.0;
1254:           element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1255:         }
1256:       } else if (coefficient_structure == 1) { /* square sinker */
1257:         PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;

1259:         opts_eta0 = 1.0;
1260:         opts_eta1 = 1.0;
1261:         opts_dx   = 0.50;
1262:         opts_dy   = 0.50;

1264:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1265:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1266:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1267:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));

1269:         for (p = 0; p < GAUSS_POINTS; p++) {
1270:           coord_x = centroid_x;
1271:           coord_y = centroid_y;
1272:           if (use_gp_coords) {
1273:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1274:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1275:           }

1277:           element_props[j][i].eta[p] = opts_eta0;
1278:           element_props[j][i].fx[p]  = 0.0;
1279:           element_props[j][i].fy[p]  = 0.0;

1281:           if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
1282:             if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
1283:               element_props[j][i].eta[p] = opts_eta1;
1284:               element_props[j][i].fx[p]  = 0.0;
1285:               element_props[j][i].fy[p]  = -1.0;
1286:             }
1287:           }
1288:         }
1289:       } else if (coefficient_structure == 2) { /* circular sinker */
1290:         PetscReal opts_eta0, opts_eta1, opts_r, radius2;

1292:         opts_eta0 = 1.0;
1293:         opts_eta1 = 1.0;
1294:         opts_r    = 0.25;

1296:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1297:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1298:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));

1300:         for (p = 0; p < GAUSS_POINTS; p++) {
1301:           coord_x = centroid_x;
1302:           coord_y = centroid_y;
1303:           if (use_gp_coords) {
1304:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1305:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1306:           }

1308:           element_props[j][i].eta[p] = opts_eta0;
1309:           element_props[j][i].fx[p]  = 0.0;
1310:           element_props[j][i].fy[p]  = 0.0;

1312:           radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
1313:           if (radius2 < opts_r * opts_r) {
1314:             element_props[j][i].eta[p] = opts_eta1;
1315:             element_props[j][i].fx[p]  = 0.0;
1316:             element_props[j][i].fy[p]  = -1.0;
1317:           }
1318:         }
1319:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1320:         PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;

1322:         opts_eta0 = 1.0;
1323:         opts_eta1 = 1.0;
1324:         opts_r    = 0.25;
1325:         opts_c0x  = 0.35; /* circle center */
1326:         opts_c0y  = 0.35;
1327:         opts_s0x  = 0.7; /* square center */
1328:         opts_s0y  = 0.7;
1329:         opts_dx   = 0.25;
1330:         opts_dy   = 0.25;
1331:         opts_phi  = 25;

1333:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1334:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1335:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1336:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL));
1337:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL));
1338:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL));
1339:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL));
1340:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1341:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1342:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL));
1343:         opts_phi *= PETSC_PI / 180;

1345:         for (p = 0; p < GAUSS_POINTS; p++) {
1346:           coord_x = centroid_x;
1347:           coord_y = centroid_y;
1348:           if (use_gp_coords) {
1349:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1350:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1351:           }

1353:           element_props[j][i].eta[p] = opts_eta0;
1354:           element_props[j][i].fx[p]  = 0.0;
1355:           element_props[j][i].fy[p]  = -0.2;

1357:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1358:           if (radius2 < opts_r * opts_r || (PetscAbs(+(coord_x - opts_s0x) * PetscCosReal(opts_phi) + (coord_y - opts_s0y) * PetscSinReal(opts_phi)) < opts_dx / 2 && PetscAbs(-(coord_x - opts_s0x) * PetscSinReal(opts_phi) + (coord_y - opts_s0y) * PetscCosReal(opts_phi)) < opts_dy / 2)) {
1359:             element_props[j][i].eta[p] = opts_eta1;
1360:             element_props[j][i].fx[p]  = 0.0;
1361:             element_props[j][i].fy[p]  = -1.0;
1362:           }
1363:         }
1364:       } else if (coefficient_structure == 4) { /* subdomain jump */
1365:         PetscReal opts_mag, opts_eta0;
1366:         PetscInt  opts_nz, px, py;
1367:         PetscBool jump;

1369:         opts_mag  = 1.0;
1370:         opts_eta0 = 1.0;
1371:         opts_nz   = 1;

1373:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
1374:         PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
1375:         PetscCall(PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL));
1376:         PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1377:         if (px % 2) {
1378:           jump = (PetscBool)(PetscGlobalRank % 2);
1379:         } else {
1380:           jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
1381:         }
1382:         for (p = 0; p < GAUSS_POINTS; p++) {
1383:           coord_x = centroid_x;
1384:           coord_y = centroid_y;
1385:           if (use_gp_coords) {
1386:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1387:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1388:           }

1390:           element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1391:           element_props[j][i].fx[p]  = 0.0;
1392:           element_props[j][i].fy[p]  = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1393:         }
1394:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
1395:     }
1396:   }
1397:   PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));

1399:   PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));

1401:   PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));

1405:   if (output_gnuplot) {
1406:     PetscCall(DMDACoordViewGnuplot2d(da_Stokes, "mesh"));
1407:     PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties"));
1408:   }

1410:   if (glvis) {
1411:     Vec         glv_prop, etaf;
1412:     PetscViewer view;
1413:     PetscInt    dim;
1414:     const char *fec = {"FiniteElementCollection: L2_2D_P1"};

1416:     PetscCall(DMGetDimension(da_Stokes, &dim));
1417:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf));
1418:     PetscCall(PetscObjectSetName((PetscObject)etaf, "viscosity"));
1419:     PetscCall(PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view));
1420:     PetscCall(PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL));
1421:     PetscCall(DMGetLocalVector(da_prop, &glv_prop));
1422:     PetscCall(DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop));
1423:     PetscCall(DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop));
1424:     PetscCall(VecSetDM(etaf, da_Stokes));
1425:     PetscCall(VecView(glv_prop, view));
1426:     PetscCall(PetscViewerDestroy(&view));
1427:     PetscCall(DMRestoreLocalVector(da_prop, &glv_prop));
1428:     PetscCall(VecDestroy(&etaf));
1429:   }

1431:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1432:   PetscCall(DMCreateMatrix(da_Stokes, &A));
1433:   PetscCall(DMCreateMatrix(da_Stokes, &B));
1434:   PetscCall(DMCreateGlobalVector(da_Stokes, &f));
1435:   PetscCall(DMCreateGlobalVector(da_Stokes, &X));

1437:   /* assemble A11 */
1438:   PetscCall(MatZeroEntries(A));
1439:   PetscCall(MatZeroEntries(B));
1440:   PetscCall(VecZeroEntries(f));

1442:   PetscCall(AssembleA_Stokes(A, da_Stokes, da_prop, properties));
1443:   PetscCall(AssembleA_PCStokes(B, da_Stokes, da_prop, properties));
1444:   /* build force vector */
1445:   PetscCall(AssembleF_Stokes(f, da_Stokes, da_prop, properties));

1447:   PetscCall(DMDABCApplyFreeSlip(da_Stokes, A, f));
1448:   PetscCall(DMDABCApplyFreeSlip(da_Stokes, B, NULL));

1450:   /* SOLVE */
1451:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
1452:   PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
1453:   PetscCall(KSPSetDM(ksp_S, da_Stokes));
1454:   PetscCall(KSPSetDMActive(ksp_S, PETSC_FALSE));
1455:   PetscCall(KSPSetOperators(ksp_S, A, B));
1456:   PetscCall(KSPSetFromOptions(ksp_S));
1457:   {
1458:     PC             pc;
1459:     const PetscInt ufields[] = {0, 1}, pfields[1] = {2};

1461:     PetscCall(KSPGetPC(ksp_S, &pc));
1462:     PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1463:     PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1464:     PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1465:   }

1467:   {
1468:     PC        pc;
1469:     PetscBool same = PETSC_FALSE;
1470:     PetscCall(KSPGetPC(ksp_S, &pc));
1471:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
1472:     if (same) {
1473:       PetscBool usedivmat = PETSC_FALSE;
1474:       PetscCall(KSPSetOperators(ksp_S, A, A));

1476:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL));
1477:       if (usedivmat) {
1478:         IS      *fields, vel;
1479:         PetscInt i, nf;

1481:         PetscCall(DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL));
1482:         PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel));
1483:         PetscCall(MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL)); /* we put 1.0 on the diagonal to pick the pressure average too */
1484:         PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
1485:         PetscCall(MatZeroRowsIS(B, vel, 0.0, NULL, NULL));
1486:         PetscCall(ISDestroy(&vel));
1487:         PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
1488:         for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
1489:         PetscCall(PetscFree(fields));
1490:       }
1491:     }
1492:   }

1494:   {
1495:     PC        pc_S;
1496:     KSP      *sub_ksp, ksp_U;
1497:     PetscInt  nsplits;
1498:     DM        da_U;
1499:     PetscBool is_pcfs;

1501:     PetscCall(KSPSetUp(ksp_S));
1502:     PetscCall(KSPGetPC(ksp_S, &pc_S));

1504:     is_pcfs = PETSC_FALSE;
1505:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs));

1507:     if (is_pcfs) {
1508:       PetscCall(PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp));
1509:       ksp_U = sub_ksp[0];
1510:       PetscCall(PetscFree(sub_ksp));

1512:       if (nsplits == 2) {
1513:         PetscCall(DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U));

1515:         PetscCall(KSPSetDM(ksp_U, da_U));
1516:         PetscCall(KSPSetDMActive(ksp_U, PETSC_FALSE));
1517:         PetscCall(DMDestroy(&da_U));
1518:         if (change) {
1519:           const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1520:           PetscCall(PetscOptionsInsertString(NULL, opt));
1521:           PetscCall(KSPSetFromOptions(ksp_U));
1522:         }
1523:         PetscCall(KSPSetUp(ksp_U));
1524:       }
1525:     }
1526:   }

1528:   PetscCall(KSPSolve(ksp_S, f, X));

1530:   PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set));
1531:   if (set) {
1532:     char       *ext;
1533:     PetscViewer viewer;
1534:     PetscBool   flg;
1535:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1536:     PetscCall(PetscStrrchr(filename, '.', &ext));
1537:     PetscCall(PetscStrcmp("vts", ext, &flg));
1538:     if (flg) {
1539:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1540:     } else {
1541:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
1542:     }
1543:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1544:     PetscCall(PetscViewerFileSetName(viewer, filename));
1545:     PetscCall(VecView(X, viewer));
1546:     PetscCall(PetscViewerDestroy(&viewer));
1547:   }
1548:   if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X"));

1550:   if (glvis) {
1551:     PetscViewer view;

1553:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &view));
1554:     PetscCall(PetscViewerSetType(view, PETSCVIEWERGLVIS));
1555:     PetscCall(VecView(X, view));
1556:     PetscCall(PetscViewerDestroy(&view));
1557:   }

1559:   PetscCall(KSPGetIterationNumber(ksp_S, &its));

1561:   if (coefficient_structure == 0) {
1562:     PetscReal opts_eta0, opts_eta1, opts_xc;
1563:     PetscInt  opts_nz, N;
1564:     DM        da_Stokes_analytic;
1565:     Vec       X_analytic;
1566:     PetscReal nrm1[3], nrm2[3], nrmI[3];

1568:     opts_eta0 = 1.0;
1569:     opts_eta1 = 1.0;
1570:     opts_xc   = 0.5;
1571:     opts_nz   = 1;

1573:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1574:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1575:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1576:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));

1578:     PetscCall(DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic));
1579:     if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic"));
1580:     PetscCall(DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic));

1582:     PetscCall(VecAXPY(X_analytic, -1.0, X));
1583:     PetscCall(VecGetSize(X_analytic, &N));
1584:     N = N / 3;

1586:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]));
1587:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]));
1588:     PetscCall(VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]));

1590:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]));
1591:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]));
1592:     PetscCall(VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]));

1594:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]));
1595:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]));
1596:     PetscCall(VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]));

1598:     PetscCall(DMDestroy(&da_Stokes_analytic));
1599:     PetscCall(VecDestroy(&X_analytic));
1600:   }

1602:   PetscCall(KSPDestroy(&ksp_S));
1603:   PetscCall(VecDestroy(&X));
1604:   PetscCall(VecDestroy(&f));
1605:   PetscCall(MatDestroy(&A));
1606:   PetscCall(MatDestroy(&B));

1608:   PetscCall(DMDestroy(&da_Stokes));
1609:   PetscCall(DMDestroy(&da_prop));

1611:   PetscCall(VecDestroy(&properties));
1612:   PetscCall(VecDestroy(&l_properties));
1613:   PetscFunctionReturn(PETSC_SUCCESS);
1614: }

1616: int main(int argc, char **args)
1617: {
1618:   PetscInt mx, my;

1620:   PetscFunctionBeginUser;
1621:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1622:   mx = my = 10;
1623:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1624:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1625:   PetscCall(solve_stokes_2d_coupled(mx, my));
1626:   PetscCall(PetscFinalize());
1627:   return 0;
1628: }

1630: /* -------------------------- helpers for boundary conditions -------------------------------- */
1631: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1632: {
1633:   DM                     cda;
1634:   Vec                    coords;
1635:   PetscInt               si, sj, nx, ny, i, j;
1636:   PetscInt               M, N;
1637:   DMDACoor2d           **_coords;
1638:   const PetscInt        *g_idx;
1639:   PetscInt              *bc_global_ids;
1640:   PetscScalar           *bc_vals;
1641:   PetscInt               nbcs;
1642:   PetscInt               n_dofs;
1643:   ISLocalToGlobalMapping ltogm;

1645:   PetscFunctionBeginUser;
1646:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1647:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1649:   PetscCall(DMGetCoordinateDM(da, &cda));
1650:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1651:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1652:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1653:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1655:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1656:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1658:   /* init the entries to -1 so VecSetValues will ignore them */
1659:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1661:   i = nx - 1;
1662:   for (j = 0; j < ny; j++) {
1663:     PetscInt local_id;

1665:     local_id = i + j * nx;

1667:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1669:     bc_vals[j] = 0.0;
1670:   }
1671:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1672:   nbcs = 0;
1673:   if ((si + nx) == (M)) nbcs = ny;

1675:   if (b) {
1676:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1677:     PetscCall(VecAssemblyBegin(b));
1678:     PetscCall(VecAssemblyEnd(b));
1679:   }
1680:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1682:   PetscCall(PetscFree(bc_vals));
1683:   PetscCall(PetscFree(bc_global_ids));

1685:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1686:   PetscFunctionReturn(PETSC_SUCCESS);
1687: }

1689: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1690: {
1691:   DM                     cda;
1692:   Vec                    coords;
1693:   PetscInt               si, sj, nx, ny, i, j;
1694:   PetscInt               M, N;
1695:   DMDACoor2d           **_coords;
1696:   const PetscInt        *g_idx;
1697:   PetscInt              *bc_global_ids;
1698:   PetscScalar           *bc_vals;
1699:   PetscInt               nbcs;
1700:   PetscInt               n_dofs;
1701:   ISLocalToGlobalMapping ltogm;

1703:   PetscFunctionBeginUser;
1704:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1705:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1707:   PetscCall(DMGetCoordinateDM(da, &cda));
1708:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1709:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1710:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1711:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1713:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1714:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1716:   /* init the entries to -1 so VecSetValues will ignore them */
1717:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1719:   i = 0;
1720:   for (j = 0; j < ny; j++) {
1721:     PetscInt local_id;

1723:     local_id = i + j * nx;

1725:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1727:     bc_vals[j] = 0.0;
1728:   }
1729:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1730:   nbcs = 0;
1731:   if (si == 0) nbcs = ny;

1733:   if (b) {
1734:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1735:     PetscCall(VecAssemblyBegin(b));
1736:     PetscCall(VecAssemblyEnd(b));
1737:   }

1739:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1741:   PetscCall(PetscFree(bc_vals));
1742:   PetscCall(PetscFree(bc_global_ids));

1744:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1745:   PetscFunctionReturn(PETSC_SUCCESS);
1746: }

1748: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1749: {
1750:   DM                     cda;
1751:   Vec                    coords;
1752:   PetscInt               si, sj, nx, ny, i, j;
1753:   PetscInt               M, N;
1754:   DMDACoor2d           **_coords;
1755:   const PetscInt        *g_idx;
1756:   PetscInt              *bc_global_ids;
1757:   PetscScalar           *bc_vals;
1758:   PetscInt               nbcs;
1759:   PetscInt               n_dofs;
1760:   ISLocalToGlobalMapping ltogm;

1762:   PetscFunctionBeginUser;
1763:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1764:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1766:   PetscCall(DMGetCoordinateDM(da, &cda));
1767:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1768:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1769:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1770:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1772:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
1773:   PetscCall(PetscMalloc1(nx, &bc_vals));

1775:   /* init the entries to -1 so VecSetValues will ignore them */
1776:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1778:   j = ny - 1;
1779:   for (i = 0; i < nx; i++) {
1780:     PetscInt local_id;

1782:     local_id = i + j * nx;

1784:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

1786:     bc_vals[i] = 0.0;
1787:   }
1788:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1789:   nbcs = 0;
1790:   if ((sj + ny) == (N)) nbcs = nx;

1792:   if (b) {
1793:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1794:     PetscCall(VecAssemblyBegin(b));
1795:     PetscCall(VecAssemblyEnd(b));
1796:   }
1797:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1799:   PetscCall(PetscFree(bc_vals));
1800:   PetscCall(PetscFree(bc_global_ids));

1802:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1803:   PetscFunctionReturn(PETSC_SUCCESS);
1804: }

1806: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1807: {
1808:   DM                     cda;
1809:   Vec                    coords;
1810:   PetscInt               si, sj, nx, ny, i, j;
1811:   PetscInt               M, N;
1812:   DMDACoor2d           **_coords;
1813:   const PetscInt        *g_idx;
1814:   PetscInt              *bc_global_ids;
1815:   PetscScalar           *bc_vals;
1816:   PetscInt               nbcs;
1817:   PetscInt               n_dofs;
1818:   ISLocalToGlobalMapping ltogm;

1820:   PetscFunctionBeginUser;
1821:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1822:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1824:   PetscCall(DMGetCoordinateDM(da, &cda));
1825:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1826:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1827:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1828:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1830:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
1831:   PetscCall(PetscMalloc1(nx, &bc_vals));

1833:   /* init the entries to -1 so VecSetValues will ignore them */
1834:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1836:   j = 0;
1837:   for (i = 0; i < nx; i++) {
1838:     PetscInt local_id;

1840:     local_id = i + j * nx;

1842:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

1844:     bc_vals[i] = 0.0;
1845:   }
1846:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1847:   nbcs = 0;
1848:   if (sj == 0) nbcs = nx;

1850:   if (b) {
1851:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1852:     PetscCall(VecAssemblyBegin(b));
1853:     PetscCall(VecAssemblyEnd(b));
1854:   }
1855:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1857:   PetscCall(PetscFree(bc_vals));
1858:   PetscCall(PetscFree(bc_global_ids));

1860:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1861:   PetscFunctionReturn(PETSC_SUCCESS);
1862: }

1864: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1865: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
1866: {
1867:   PetscFunctionBeginUser;
1868:   PetscCall(BCApplyZero_NORTH(da_Stokes, 1, A, f));
1869:   PetscCall(BCApplyZero_EAST(da_Stokes, 0, A, f));
1870:   PetscCall(BCApplyZero_SOUTH(da_Stokes, 1, A, f));
1871:   PetscCall(BCApplyZero_WEST(da_Stokes, 0, A, f));
1872:   PetscFunctionReturn(PETSC_SUCCESS);
1873: }

1875: /*TEST

1877:    build:
1878:       requires: !complex !single

1880:    test:
1881:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short

1883:    testset:
1884:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_fieldsplit_u_ksp_type preonly -stokes_fieldsplit_u_pc_type lu -stokes_fieldsplit_p_ksp_type preonly -stokes_fieldsplit_p_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1885:       test:
1886:          suffix: 2
1887:          args:
1888:          output_file: output/ex43_1.out
1889:       test:
1890:          requires: mumps
1891:          suffix: 2_mumps
1892:          args: -change true -stokes_ksp_view
1893:          output_file: output/ex43_2_mumps.out
1894:          # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1895:          filter:  grep -E -v "(INFOG\([^7]|INFO\(|\[0\])" | grep -v "            total: nonzeros="

1897:    test:
1898:       suffix: 3
1899:       nsize: 4
1900:       args: -stokes_pc_use_amat -stokes_ksp_type gcr -stokes_ksp_gcr_restart 60 -stokes_ksp_norm_type unpreconditioned -stokes_ksp_rtol 1.e-2 -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

1902:    test:
1903:       suffix: 4
1904:       nsize: 4
1905:       args: -stokes_ksp_type pipegcr -stokes_ksp_pipegcr_mmax 60 -stokes_ksp_pipegcr_unroll_w 1 -stokes_ksp_norm_type natural -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

1907:    test:
1908:       suffix: 5
1909:       nsize: 4
1910:       args: -stokes_pc_fieldsplit_off_diag_use_amat -stokes_ksp_type pipegcr -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type bjacobi -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type bjacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short -stokes_ksp_view

1912:    test:
1913:       suffix: 6
1914:       nsize: 8
1915:       args: -stokes_ksp_view -stokes_pc_type mg -stokes_pc_mg_levels 2 -stokes_mg_coarse_pc_type telescope -stokes_mg_coarse_pc_telescope_reduction_factor 2 -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_coarse_pc_telescope_subcomm_type contiguous

1917:    test:
1918:       suffix: bjacobi
1919:       nsize: 4
1920:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason

1922:    test:
1923:       suffix: bjacobi_baij
1924:       nsize: 4
1925:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1926:       output_file: output/ex43_bjacobi.out

1928:    test:
1929:       suffix: nested_gmg
1930:       nsize: 4
1931:       args: -stokes_pc_fieldsplit_off_diag_use_amat -mx 16 -my 16 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_fieldsplit_u_pc_type mg -stokes_fieldsplit_u_pc_mg_levels 5 -stokes_fieldsplit_u_pc_mg_galerkin pmat -stokes_fieldsplit_u_ksp_type cg -stokes_fieldsplit_u_ksp_rtol 1.0e-4 -stokes_fieldsplit_u_mg_levels_pc_type jacobi -solcx_eta0 1.0e4 -stokes_fieldsplit_u_ksp_converged_reason -stokes_ksp_converged_reason -stokes_fieldsplit_p_sub_pc_factor_zeropivot 1.e-8

1933:    test:
1934:       suffix: fetidp
1935:       nsize: 8
1936:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_converged_reason -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1

1938:    test:
1939:       suffix: fetidp_unsym
1940:       nsize: 8
1941:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_monitor_true_residual -stokes_ksp_converged_reason -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd

1943:    test:
1944:       suffix: bddc_stokes_deluxe
1945:       nsize: 8
1946:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

1948:    test:
1949:       suffix: bddc_stokes_subdomainjump_deluxe
1950:       nsize: 9
1951:       args: -c_str 4 -jump_magnitude 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

1953: TEST*/
```