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;

 93:   VecGetDM(Vf[0], &stokes_da);
 94:   DMDAVecGetArrayRead(properties_da, V, &props);
 95:   DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL);
 96:   DMDAGetElementsSizes(stokes_da, &mx, &my, NULL);
 97:   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:   VecRestoreArray(Vf[0], &array);
105:   DMDAVecRestoreArrayRead(properties_da, V, &props);
106:   return 0;
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: {
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:   return 0;
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;

261:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

263:   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:   PetscMalloc1(cpu_x, &LX);
269:   PetscMalloc1(cpu_y, &LY);

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

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

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

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

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

303:   *_lx = LX;
304:   *_ly = LY;

306:   VecDestroy(&vlx);
307:   VecDestroy(&vly);
308:   return 0;
309: }

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

322:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
323:   PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
324:   PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);

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

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

344:   PetscFClose(PETSC_COMM_SELF, fp);
345:   return 0;
346: }

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

361:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
362:   PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
363:   PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);

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

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

381:   DMCreateLocalVector(da, &local_fields);
382:   DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
383:   DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
384:   VecGetArray(local_fields, &_fields);

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

391:       coord_x = _coords[j][i].x;
392:       coord_y = _coords[j][i].y;

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

405:   DMDAVecRestoreArray(cda, coords, &_coords);

407:   PetscFClose(PETSC_COMM_SELF, fp);
408:   return 0;
409: }

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

423:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
424:   PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
425:   PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);

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

438:   DMGetCoordinateDM(da, &cda);
439:   DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);

441:   DMCreateLocalVector(da, &local_fields);
442:   DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
443:   DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
444:   DMDAVecGetArray(da, local_fields, &_coefficients);

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

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

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

456:         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]));
457:         PetscFPrintf(PETSC_COMM_SELF, fp, "\n");
458:       }
459:     }
460:   }
461:   DMDAVecRestoreArray(da, local_fields, &_coefficients);
462:   VecDestroy(&local_fields);

464:   PetscFClose(PETSC_COMM_SELF, fp);
465:   return 0;
466: }

468: 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)
469: {
470:   PetscInt ij;
471:   PetscInt r, c, nc;

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

480: /*
481:  D = [ 2.eta   0   0   ]
482:      [   0   2.eta 0   ]
483:      [   0     0   eta ]

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

499:   /* define quadrature rule */
500:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

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

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

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

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

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

549:   /* define quadrature rule */
550:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

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

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

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

578:   PetscMemzero(Ge, sizeof(Ge));
579:   FormGradientOperatorQ1(Ge, coords);

581:   nr_g = U_DOFS * NODES_PER_EL;
582:   nc_g = P_DOFS * NODES_PER_EL;

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

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

599:   /* define quadrature rule */
600:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

609:     for (i = 0; i < NODES_PER_EL; i++) {
610:       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);
611:     }
612:   }

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

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

634:   /* define quadrature rule */
635:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

644:     for (i = 0; i < NODES_PER_EL; i++) {
645:       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];
646:     }
647:   }

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

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

669:   /* define quadrature rule */
670:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

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

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

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

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

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

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

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

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

747:       /* form element stiffness matrix */
748:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
749:       FormGradientOperatorQ1(Ge, el_coords);
750:       FormDivergenceOperatorQ1(De, el_coords);
751:       FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);

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

764:   DMDAVecRestoreArray(cda, coords, &_coords);

766:   DMDAVecRestoreArray(properties_da, local_properties, &props);
767:   VecDestroy(&local_properties);
768:   return 0;
769: }

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

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

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

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

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

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

817:       /* form element stiffness matrix */
818:       FormStressOperatorQ1(Ae, el_coords, prop_eta);
819:       FormGradientOperatorQ1(Ge, el_coords);
820:       FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);

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

832:   DMDAVecRestoreArray(cda, coords, &_coords);

834:   DMDAVecRestoreArray(properties_da, local_properties, &props);
835:   VecDestroy(&local_properties);
836:   return 0;
837: }

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

844:   for (n = 0; n < 4; n++) {
845:     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];
846:     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];
847:     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];
848:   }
849:   return 0;
850: }

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

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

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

882:   /* get access to the vector */
883:   DMGetLocalVector(stokes_da, &local_F);
884:   VecZeroEntries(local_F);
885:   DMDAVecGetArray(stokes_da, local_F, &ff);

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

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

898:       /* initialise element stiffness matrix */
899:       PetscMemzero(Fe, sizeof(Fe));
900:       PetscMemzero(He, sizeof(He));

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

905:       /* insert element matrix into global matrix */
906:       DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej);

908:       DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He);
909:     }
910:   }

912:   DMDAVecRestoreArray(stokes_da, local_F, &ff);
913:   DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F);
914:   DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F);
915:   DMRestoreLocalVector(stokes_da, &local_F);

917:   DMDAVecRestoreArray(cda, coords, &_coords);

919:   DMDAVecRestoreArray(properties_da, local_properties, &props);
920:   DMRestoreLocalVector(properties_da, &local_properties);
921:   return 0;
922: }

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

934:   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);
935:   DMSetFromOptions(da);
936:   DMSetUp(da);
937:   DMDASetFieldName(da, 0, "anlytic_Vx");
938:   DMDASetFieldName(da, 1, "anlytic_Vy");
939:   DMDASetFieldName(da, 2, "analytic_P");

941:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.);

943:   DMGetCoordinatesLocal(da, &coords);
944:   DMGetCoordinateDM(da, &cda);
945:   DMDAVecGetArray(cda, coords, &_coords);

947:   DMCreateGlobalVector(da, &X);
948:   DMDAVecGetArray(da, X, &_stokes);

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

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

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

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

968:   *_da = da;
969:   *_X  = X;
970:   return 0;
971: }

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

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

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

1015:   /* define quadrature rule */
1016:   ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

1018:   /* setup for coords */
1019:   DMGetCoordinateDM(stokes_da, &cda);
1020:   DMGetCoordinatesLocal(stokes_da, &coords);
1021:   DMDAVecGetArray(cda, coords, &_coords);

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

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

1035:   DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1036:   DMGetBoundingBox(stokes_da, xymin, xymax);

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

1040:   tp_L2 = tu_L2 = tu_H1 = 0.0;

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

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

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

1064:           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);

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

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

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

1089:   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));

1091:   DMDAVecRestoreArray(cda, coords, &_coords);

1093:   DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic);
1094:   VecDestroy(&X_analytic_local);
1095:   DMDAVecRestoreArray(stokes_da, X_local, &stokes);
1096:   VecDestroy(&X_local);
1097:   return 0;
1098: }

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


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

1128:   /* Generate the da for velocity and pressure */
1129:   /*
1130:   We use Q1 elements for the temperature.
1131:   FEM has a 9-point stencil (BOX) or connectivity pattern
1132:   Num nodes in each direction is mx+1, my+1
1133:   */
1134:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1135:   p_dof         = P_DOFS; /* p - pressure */
1136:   dof           = u_dof + p_dof;
1137:   stencil_width = 1;
1138:   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);

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

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

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

1155:   prop_dof           = (int)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1156:   prop_stencil_width = 0;
1157:   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);
1158:   DMSetFromOptions(da_prop);
1159:   DMSetUp(da_prop);
1160:   PetscFree(lx);
1161:   PetscFree(ly);

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

1168:   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);

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

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

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

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

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

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

1194:       GetElementCoords(_vel_coords, i, j, el_coords);
1195:       ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1356:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1357:           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)) {
1358:             element_props[j][i].eta[p] = opts_eta1;
1359:             element_props[j][i].fx[p]  = 0.0;
1360:             element_props[j][i].fy[p]  = -1.0;
1361:           }
1362:         }
1363:       } else if (coefficient_structure == 4) { /* subdomain jump */
1364:         PetscReal opts_mag, opts_eta0;
1365:         PetscInt  opts_nz, px, py;
1366:         PetscBool jump;

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

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

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

1398:   DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords);

1400:   DMDAVecRestoreArray(da_prop, l_properties, &element_props);
1401:   DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties);
1402:   DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties);

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

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

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

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

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

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

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

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

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

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

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

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

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

1500:     KSPSetUp(ksp_S);
1501:     KSPGetPC(ksp_S, &pc_S);

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

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

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

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

1527:   KSPSolve(ksp_S, f, X);

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

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

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

1558:   KSPGetIterationNumber(ksp_S, &its);

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

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

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

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

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

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

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

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

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

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

1607:   DMDestroy(&da_Stokes);
1608:   DMDestroy(&da_prop);

1610:   VecDestroy(&properties);
1611:   VecDestroy(&l_properties);
1612:   return 0;
1613: }

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

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

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

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

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

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

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

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

1664:     local_id = i + j * nx;

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

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

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

1681:   PetscFree(bc_vals);
1682:   PetscFree(bc_global_ids);

1684:   DMDAVecRestoreArray(cda, coords, &_coords);
1685:   return 0;
1686: }

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

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

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

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

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

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

1722:     local_id = i + j * nx;

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

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

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

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

1740:   PetscFree(bc_vals);
1741:   PetscFree(bc_global_ids);

1743:   DMDAVecRestoreArray(cda, coords, &_coords);
1744:   return 0;
1745: }

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

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

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

1771:   PetscMalloc1(nx, &bc_global_ids);
1772:   PetscMalloc1(nx, &bc_vals);

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

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

1781:     local_id = i + j * nx;

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

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

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

1798:   PetscFree(bc_vals);
1799:   PetscFree(bc_global_ids);

1801:   DMDAVecRestoreArray(cda, coords, &_coords);
1802:   return 0;
1803: }

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

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

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

1829:   PetscMalloc1(nx, &bc_global_ids);
1830:   PetscMalloc1(nx, &bc_vals);

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

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

1839:     local_id = i + j * nx;

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

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

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

1856:   PetscFree(bc_vals);
1857:   PetscFree(bc_global_ids);

1859:   DMDAVecRestoreArray(cda, coords, &_coords);
1860:   return 0;
1861: }

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

1874: /*TEST

1876:    build:
1877:       requires: !complex !single

1879:    test:
1880:       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

1882:    testset:
1883:       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
1884:       test:
1885:          suffix: 2
1886:          args:
1887:          output_file: output/ex43_1.out
1888:       test:
1889:          requires: mumps
1890:          suffix: 2_mumps
1891:          args: -change true -stokes_ksp_view
1892:          output_file: output/ex43_2_mumps.out
1893:          # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1894:          filter:  grep -E -v "(INFOG\([^7]|INFO\(|\[0\])"

1896:    test:
1897:       suffix: 3
1898:       nsize: 4
1899:       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

1901:    test:
1902:       suffix: 4
1903:       nsize: 4
1904:       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

1906:    test:
1907:       suffix: 5
1908:       nsize: 4
1909:       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

1911:    test:
1912:       suffix: 6
1913:       nsize: 8
1914:       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

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

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

1927:    test:
1928:       suffix: nested_gmg
1929:       nsize: 4
1930:       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

1932:    test:
1933:       suffix: fetidp
1934:       nsize: 8
1935:       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

1937:    test:
1938:       suffix: fetidp_unsym
1939:       nsize: 8
1940:       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

1942:    test:
1943:       suffix: bddc_stokes_deluxe
1944:       nsize: 8
1945:       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

1947:    test:
1948:       suffix: bddc_stokes_subdomainjump_deluxe
1949:       nsize: 9
1950:       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

1952: TEST*/