Actual source code: ex71.c

  1: static char help[] = "Poiseuille Flow in 2d and 3d channels with finite elements.\n\
  2: We solve the Poiseuille flow problem in a rectangular\n\
  3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";

  5: /*F
  6: A Poiseuille flow is a steady-state isoviscous Stokes flow in a pipe of constant cross-section. We discretize using the
  7: finite element method on an unstructured mesh. The weak form equations are
  8: \begin{align*}
  9:   < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > + < v, \Delta \hat n >_{\Gamma_o} = 0
 10:   < q, \nabla\cdot u >                                                                                 = 0
 11: \end{align*}
 12: where $\nu$ is the kinematic viscosity, $\Delta$ is the pressure drop per unit length, assuming that pressure is 0 on
 13: the left edge, and $\Gamma_o$ is the outlet boundary at the right edge of the pipe. The normal velocity will be zero at
 14: the wall, but we will allow a fixed tangential velocity $u_0$.

 16: In order to test our global to local basis transformation, we will allow the pipe to be at an angle $\alpha$ to the
 17: coordinate axes.

 19: For visualization, use

 21:   -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append
 22: F*/

 24: #include <petscdmplex.h>
 25: #include <petscsnes.h>
 26: #include <petscds.h>
 27: #include <petscbag.h>

 29: typedef struct {
 30:   PetscReal Delta; /* Pressure drop per unit length */
 31:   PetscReal nu;    /* Kinematic viscosity */
 32:   PetscReal u_0;   /* Tangential velocity at the wall */
 33:   PetscReal alpha; /* Angle of pipe wall to x-axis */
 34: } Parameter;

 36: typedef struct {
 37:   PetscBag bag; /* Holds problem parameters */
 38: } AppCtx;

 40: /*
 41:   In 2D, plane Poiseuille flow has exact solution:

 43:     u = \Delta/(2 \nu) y (1 - y) + u_0
 44:     v = 0
 45:     p = -\Delta x
 46:     f = 0

 48:   so that

 50:     -\nu \Delta u + \nabla p + f = <\Delta, 0> + <-\Delta, 0> + <0, 0> = 0
 51:     \nabla \cdot u               = 0 + 0                               = 0

 53:   In 3D we use exact solution:

 55:     u = \Delta/(4 \nu) (y (1 - y) + z (1 - z)) + u_0
 56:     v = 0
 57:     w = 0
 58:     p = -\Delta x
 59:     f = 0

 61:   so that

 63:     -\nu \Delta u + \nabla p + f = <Delta, 0, 0> + <-Delta, 0, 0> + <0, 0, 0> = 0
 64:     \nabla \cdot u               = 0 + 0 + 0                                  = 0

 66:   Note that these functions use coordinates X in the global (rotated) frame
 67: */
 68: PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
 69: {
 70:   Parameter *param = (Parameter *)ctx;
 71:   PetscReal  Delta = param->Delta;
 72:   PetscReal  nu    = param->nu;
 73:   PetscReal  u_0   = param->u_0;
 74:   PetscReal  fac   = (PetscReal)(dim - 1);

 76:   u[0] = u_0;
 77:   for (PetscInt d = 1; d < dim; ++d) u[0] += Delta / (fac * 2.0 * nu) * X[d] * (1.0 - X[d]);
 78:   for (PetscInt d = 1; d < dim; ++d) u[d] = 0.0;
 79:   return PETSC_SUCCESS;
 80: }

 82: PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
 83: {
 84:   Parameter *param = (Parameter *)ctx;
 85:   PetscReal  Delta = param->Delta;

 87:   p[0] = -Delta * X[0];
 88:   return PETSC_SUCCESS;
 89: }

 91: PetscErrorCode wall_velocity(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
 92: {
 93:   Parameter *param = (Parameter *)ctx;
 94:   PetscReal  u_0   = param->u_0;

 96:   u[0] = u_0;
 97:   for (PetscInt d = 1; d < dim; ++d) u[d] = 0.0;
 98:   return PETSC_SUCCESS;
 99: }

101: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z}
102:    u[Ncomp]          = {p} */
103: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
104: {
105:   const PetscReal nu = PetscRealPart(constants[1]);
106:   const PetscInt  Nc = dim;

108:   for (PetscInt c = 0; c < Nc; ++c) {
109:     for (PetscInt d = 0; d < dim; ++d) {
110:       /* f1[c*dim+d] = 0.5*nu*(u_x[c*dim+d] + u_x[d*dim+c]); */
111:       f1[c * dim + d] = nu * u_x[c * dim + d];
112:     }
113:     f1[c * dim + c] -= u[uOff[1]];
114:   }
115: }

117: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z} */
118: void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
119: {
120:   f0[0] = 0.0;
121:   for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
122: }

124: /* Residual functions are in reference coordinates */
125: static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
126: {
127:   const PetscReal Delta = PetscRealPart(constants[0]);
128:   PetscReal       alpha = PetscRealPart(constants[3]);
129:   PetscReal       X     = PetscCosReal(alpha) * x[0] + PetscSinReal(alpha) * x[1];

131:   for (PetscInt d = 0; d < dim; ++d) f0[d] = -Delta * X * n[d];
132: }

134: /* < q, \nabla\cdot u >
135:    NcompI = 1, NcompJ = dim */
136: void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
137: {
138:   for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
139: }

141: /* -< \nabla\cdot v, p >
142:     NcompI = dim, NcompJ = 1 */
143: void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
144: {
145:   for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
146: }

148: /* < \nabla v, \nabla u + {\nabla u}^T >
149:    This just gives \nabla u, give the perdiagonal for the transpose */
150: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
151: {
152:   const PetscReal nu = PetscRealPart(constants[1]);
153:   const PetscInt  Nc = dim;

155:   for (PetscInt c = 0; c < Nc; ++c) {
156:     for (PetscInt d = 0; d < dim; ++d) g3[((c * Nc + c) * dim + d) * dim + d] = nu;
157:   }
158: }

160: static PetscErrorCode SetupParameters(AppCtx *user)
161: {
162:   PetscBag   bag;
163:   Parameter *p;

165:   PetscFunctionBeginUser;
166:   /* setup PETSc parameter bag */
167:   PetscCall(PetscBagGetData(user->bag, &p));
168:   PetscCall(PetscBagSetName(user->bag, "par", "Poiseuille flow parameters"));
169:   bag = user->bag;
170:   PetscCall(PetscBagRegisterReal(bag, &p->Delta, 1.0, "Delta", "Pressure drop per unit length"));
171:   PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
172:   PetscCall(PetscBagRegisterReal(bag, &p->u_0, 0.0, "u_0", "Tangential velocity at the wall"));
173:   PetscCall(PetscBagRegisterReal(bag, &p->alpha, 0.0, "alpha", "Angle of pipe wall to x-axis"));
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
178: {
179:   PetscFunctionBeginUser;
180:   PetscCall(DMCreate(comm, dm));
181:   PetscCall(DMSetType(*dm, DMPLEX));
182:   PetscCall(DMSetFromOptions(*dm));
183:   {
184:     Parameter   *param;
185:     Vec          coordinates;
186:     PetscScalar *coords;
187:     PetscReal    alpha;
188:     PetscInt     cdim, N, bs, i;

190:     PetscCall(DMGetCoordinateDim(*dm, &cdim));
191:     PetscCall(DMGetCoordinates(*dm, &coordinates));
192:     PetscCall(VecGetLocalSize(coordinates, &N));
193:     PetscCall(VecGetBlockSize(coordinates, &bs));
194:     PetscCheck(bs == cdim, comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %" PetscInt_FMT " != embedding dimension %" PetscInt_FMT, bs, cdim);
195:     PetscCall(VecGetArray(coordinates, &coords));
196:     PetscCall(PetscBagGetData(user->bag, &param));
197:     alpha = param->alpha;
198:     for (i = 0; i < N; i += cdim) {
199:       PetscScalar x = coords[i + 0];
200:       PetscScalar y = coords[i + 1];

202:       coords[i + 0] = PetscCosReal(alpha) * x - PetscSinReal(alpha) * y;
203:       coords[i + 1] = PetscSinReal(alpha) * x + PetscCosReal(alpha) * y;
204:     }
205:     PetscCall(VecRestoreArray(coordinates, &coords));
206:     PetscCall(DMSetCoordinates(*dm, coordinates));
207:   }
208:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
213: {
214:   PetscDS       ds;
215:   PetscWeakForm wf;
216:   DMLabel       label;
217:   Parameter    *ctx;
218:   PetscInt      id, bd;

220:   PetscFunctionBeginUser;
221:   PetscCall(PetscBagGetData(user->bag, &ctx));
222:   PetscCall(DMGetDS(dm, &ds));
223:   PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
224:   PetscCall(PetscDSSetResidual(ds, 1, f0_p, NULL));
225:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
226:   PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL));
227:   PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL));

229:   id = 2;
230:   PetscCall(DMGetLabel(dm, "marker", &label));
231:   PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
232:   PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
233:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_bd_u, 0, NULL));
234:   /* Setup constants */
235:   {
236:     Parameter  *param;
237:     PetscScalar constants[4];

239:     PetscCall(PetscBagGetData(user->bag, &param));

241:     constants[0] = param->Delta;
242:     constants[1] = param->nu;
243:     constants[2] = param->u_0;
244:     constants[3] = param->alpha;
245:     PetscCall(PetscDSSetConstants(ds, 4, constants));
246:   }
247:   /* Setup Boundary Conditions */
248:   id = 3;
249:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)wall_velocity, NULL, ctx, NULL));
250:   id = 1;
251:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)wall_velocity, NULL, ctx, NULL));
252:   /* Setup exact solution */
253:   PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, ctx));
254:   PetscCall(PetscDSSetExactSolution(ds, 1, linear_p, ctx));
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }

258: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
259: {
260:   DM         cdm = dm;
261:   PetscFE    fe[2];
262:   Parameter *param;
263:   PetscBool  simplex;
264:   PetscInt   dim;
265:   MPI_Comm   comm;

267:   PetscFunctionBeginUser;
268:   PetscCall(DMGetDimension(dm, &dim));
269:   PetscCall(DMPlexIsSimplex(dm, &simplex));
270:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
271:   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
272:   PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
273:   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
274:   PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
275:   PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
276:   /* Set discretization and boundary conditions for each mesh */
277:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
278:   PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
279:   PetscCall(DMCreateDS(dm));
280:   PetscCall(SetupProblem(dm, user));
281:   PetscCall(PetscBagGetData(user->bag, &param));
282:   while (cdm) {
283:     PetscCall(DMCopyDisc(dm, cdm));
284:     PetscCall(DMPlexCreateBasisRotation(cdm, param->alpha, 0.0, 0.0));
285:     PetscCall(DMGetCoarseDM(cdm, &cdm));
286:   }
287:   PetscCall(PetscFEDestroy(&fe[0]));
288:   PetscCall(PetscFEDestroy(&fe[1]));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: int main(int argc, char **argv)
293: {
294:   SNES   snes; /* nonlinear solver */
295:   DM     dm;   /* problem definition */
296:   Vec    u, r; /* solution and residual */
297:   AppCtx user; /* user-defined work context */

299:   PetscFunctionBeginUser;
300:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
301:   PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
302:   PetscCall(SetupParameters(&user));
303:   PetscCall(PetscBagSetFromOptions(user.bag));
304:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
305:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
306:   PetscCall(SNESSetDM(snes, dm));
307:   PetscCall(DMSetApplicationContext(dm, &user));
308:   /* Setup problem */
309:   PetscCall(SetupDiscretization(dm, &user));
310:   PetscCall(DMPlexCreateClosureIndex(dm, NULL));

312:   PetscCall(DMCreateGlobalVector(dm, &u));
313:   PetscCall(VecDuplicate(u, &r));

315:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));

317:   PetscCall(SNESSetFromOptions(snes));

319:   {
320:     PetscDS             ds;
321:     PetscSimplePointFn *exactFuncs[2];
322:     void               *ctxs[2];

324:     PetscCall(DMGetDS(dm, &ds));
325:     PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
326:     PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
327:     PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u));
328:     PetscCall(PetscObjectSetName((PetscObject)u, "Exact Solution"));
329:     PetscCall(VecViewFromOptions(u, NULL, "-exact_vec_view"));
330:   }
331:   PetscCall(DMSNESCheckFromOptions(snes, u));
332:   PetscCall(VecSet(u, 0.0));
333:   PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
334:   PetscCall(SNESSolve(snes, NULL, u));
335:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

337:   PetscCall(VecDestroy(&u));
338:   PetscCall(VecDestroy(&r));
339:   PetscCall(DMDestroy(&dm));
340:   PetscCall(SNESDestroy(&snes));
341:   PetscCall(PetscBagDestroy(&user.bag));
342:   PetscCall(PetscFinalize());
343:   return 0;
344: }

346: /*TEST

348:   # Convergence
349:   test:
350:     suffix: 2d_quad_q1_p0_conv
351:     requires: !single
352:     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 \
353:       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
354:       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
355:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
356:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
357:         -fieldsplit_velocity_pc_type lu \
358:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
359:   test:
360:     suffix: 2d_quad_q1_p0_conv_u0
361:     requires: !single
362:     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 \
363:       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
364:       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
365:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
366:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
367:         -fieldsplit_velocity_pc_type lu \
368:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
369:   test:
370:     suffix: 2d_quad_q1_p0_conv_u0_alpha
371:     requires: !single
372:     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 -alpha 0.3927 \
373:       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
374:       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
375:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
376:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
377:         -fieldsplit_velocity_pc_type lu \
378:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
379:   test:
380:     suffix: 2d_quad_q1_p0_conv_gmg_vanka
381:     requires: !single long_runtime
382:     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
383:       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
384:       -snes_convergence_estimate -convest_num_refine 1 -snes_error_if_not_converged \
385:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
386:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
387:         -fieldsplit_velocity_pc_type mg \
388:           -fieldsplit_velocity_mg_levels_pc_type patch -fieldsplit_velocity_mg_levels_pc_patch_exclude_subspaces 1 \
389:           -fieldsplit_velocity_mg_levels_pc_patch_construct_codim 0 -fieldsplit_velocity_mg_levels_pc_patch_construct_type vanka \
390:         -fieldsplit_pressure_ksp_rtol 1e-5 -fieldsplit_pressure_pc_type jacobi
391:   test:
392:     suffix: 2d_tri_p2_p1_conv
393:     requires: triangle !single
394:     args: -dm_plex_separate_marker -dm_refine 1 \
395:       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
396:       -dmsnes_check .001 -snes_error_if_not_converged \
397:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
398:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
399:         -fieldsplit_velocity_pc_type lu \
400:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
401:   test:
402:     suffix: 2d_tri_p2_p1_conv_u0_alpha
403:     requires: triangle !single
404:     args: -dm_plex_separate_marker -dm_refine 0 -u_0 0.125 -alpha 0.3927 \
405:       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
406:       -dmsnes_check .001 -snes_error_if_not_converged \
407:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
408:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
409:         -fieldsplit_velocity_pc_type lu \
410:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
411:   test:
412:     suffix: 2d_tri_p2_p1_conv_gmg_vcycle
413:     TODO: broken (requires subDMs hooks)
414:     requires: triangle !single
415:     args: -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
416:       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
417:       -dmsnes_check .001 -snes_error_if_not_converged \
418:       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
419:       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
420:         -fieldsplit_velocity_pc_type mg \
421:         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
422: TEST*/