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, void *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);
75:   PetscInt   d;

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

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

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

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

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

103: /* 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}
104:    u[Ncomp]          = {p} */
105: 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[])
106: {
107:   const PetscReal nu = PetscRealPart(constants[1]);
108:   const PetscInt  Nc = dim;
109:   PetscInt        c, d;

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

120: /* 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} */
121: 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[])
122: {
123:   PetscInt d;
124:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
125: }

127: /* Residual functions are in reference coordinates */
128: 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[])
129: {
130:   const PetscReal Delta = PetscRealPart(constants[0]);
131:   PetscReal       alpha = PetscRealPart(constants[3]);
132:   PetscReal       X     = PetscCosReal(alpha) * x[0] + PetscSinReal(alpha) * x[1];
133:   PetscInt        d;

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

138: /* < q, \nabla\cdot u >
139:    NcompI = 1, NcompJ = dim */
140: 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[])
141: {
142:   PetscInt d;
143:   for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
144: }

146: /* -< \nabla\cdot v, p >
147:     NcompI = dim, NcompJ = 1 */
148: 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[])
149: {
150:   PetscInt d;
151:   for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
152: }

154: /* < \nabla v, \nabla u + {\nabla u}^T >
155:    This just gives \nabla u, give the perdiagonal for the transpose */
156: 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[])
157: {
158:   const PetscReal nu = PetscRealPart(constants[1]);
159:   const PetscInt  Nc = dim;
160:   PetscInt        c, d;

162:   for (c = 0; c < Nc; ++c) {
163:     for (d = 0; d < dim; ++d) g3[((c * Nc + c) * dim + d) * dim + d] = nu;
164:   }
165: }

167: static PetscErrorCode SetupParameters(AppCtx *user)
168: {
169:   PetscBag   bag;
170:   Parameter *p;

172:   PetscFunctionBeginUser;
173:   /* setup PETSc parameter bag */
174:   PetscCall(PetscBagGetData(user->bag, (void **)&p));
175:   PetscCall(PetscBagSetName(user->bag, "par", "Poiseuille flow parameters"));
176:   bag = user->bag;
177:   PetscCall(PetscBagRegisterReal(bag, &p->Delta, 1.0, "Delta", "Pressure drop per unit length"));
178:   PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
179:   PetscCall(PetscBagRegisterReal(bag, &p->u_0, 0.0, "u_0", "Tangential velocity at the wall"));
180:   PetscCall(PetscBagRegisterReal(bag, &p->alpha, 0.0, "alpha", "Angle of pipe wall to x-axis"));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
185: {
186:   PetscFunctionBeginUser;
187:   PetscCall(DMCreate(comm, dm));
188:   PetscCall(DMSetType(*dm, DMPLEX));
189:   PetscCall(DMSetFromOptions(*dm));
190:   {
191:     Parameter   *param;
192:     Vec          coordinates;
193:     PetscScalar *coords;
194:     PetscReal    alpha;
195:     PetscInt     cdim, N, bs, i;

197:     PetscCall(DMGetCoordinateDim(*dm, &cdim));
198:     PetscCall(DMGetCoordinates(*dm, &coordinates));
199:     PetscCall(VecGetLocalSize(coordinates, &N));
200:     PetscCall(VecGetBlockSize(coordinates, &bs));
201:     PetscCheck(bs == cdim, comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %" PetscInt_FMT " != embedding dimension %" PetscInt_FMT, bs, cdim);
202:     PetscCall(VecGetArray(coordinates, &coords));
203:     PetscCall(PetscBagGetData(user->bag, (void **)&param));
204:     alpha = param->alpha;
205:     for (i = 0; i < N; i += cdim) {
206:       PetscScalar x = coords[i + 0];
207:       PetscScalar y = coords[i + 1];

209:       coords[i + 0] = PetscCosReal(alpha) * x - PetscSinReal(alpha) * y;
210:       coords[i + 1] = PetscSinReal(alpha) * x + PetscCosReal(alpha) * y;
211:     }
212:     PetscCall(VecRestoreArray(coordinates, &coords));
213:     PetscCall(DMSetCoordinates(*dm, coordinates));
214:   }
215:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }

219: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
220: {
221:   PetscDS       ds;
222:   PetscWeakForm wf;
223:   DMLabel       label;
224:   Parameter    *ctx;
225:   PetscInt      id, bd;

227:   PetscFunctionBeginUser;
228:   PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
229:   PetscCall(DMGetDS(dm, &ds));
230:   PetscCall(PetscDSSetResidual(ds, 0, NULL, f1_u));
231:   PetscCall(PetscDSSetResidual(ds, 1, f0_p, NULL));
232:   PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
233:   PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL));
234:   PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL));

236:   id = 2;
237:   PetscCall(DMGetLabel(dm, "marker", &label));
238:   PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
239:   PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
240:   PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_bd_u, 0, NULL));
241:   /* Setup constants */
242:   {
243:     Parameter  *param;
244:     PetscScalar constants[4];

246:     PetscCall(PetscBagGetData(user->bag, (void **)&param));

248:     constants[0] = param->Delta;
249:     constants[1] = param->nu;
250:     constants[2] = param->u_0;
251:     constants[3] = param->alpha;
252:     PetscCall(PetscDSSetConstants(ds, 4, constants));
253:   }
254:   /* Setup Boundary Conditions */
255:   id = 3;
256:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall", label, 1, &id, 0, 0, NULL, (void (*)(void))wall_velocity, NULL, ctx, NULL));
257:   id = 1;
258:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall", label, 1, &id, 0, 0, NULL, (void (*)(void))wall_velocity, NULL, ctx, NULL));
259:   /* Setup exact solution */
261:   PetscCall(PetscDSSetExactSolution(ds, 1, linear_p, ctx));
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

265: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
266: {
267:   DM         cdm = dm;
268:   PetscFE    fe[2];
269:   Parameter *param;
270:   PetscBool  simplex;
271:   PetscInt   dim;
272:   MPI_Comm   comm;

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

299: int main(int argc, char **argv)
300: {
301:   SNES   snes; /* nonlinear solver */
302:   DM     dm;   /* problem definition */
303:   Vec    u, r; /* solution and residual */
304:   AppCtx user; /* user-defined work context */

306:   PetscFunctionBeginUser;
307:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
308:   PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
309:   PetscCall(SetupParameters(&user));
310:   PetscCall(PetscBagSetFromOptions(user.bag));
311:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
312:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
313:   PetscCall(SNESSetDM(snes, dm));
314:   PetscCall(DMSetApplicationContext(dm, &user));
315:   /* Setup problem */
316:   PetscCall(SetupDiscretization(dm, &user));
317:   PetscCall(DMPlexCreateClosureIndex(dm, NULL));

319:   PetscCall(DMCreateGlobalVector(dm, &u));
320:   PetscCall(VecDuplicate(u, &r));

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

324:   PetscCall(SNESSetFromOptions(snes));

326:   {
327:     PetscDS             ds;
328:     PetscSimplePointFn *exactFuncs[2];
329:     void               *ctxs[2];

331:     PetscCall(DMGetDS(dm, &ds));
332:     PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
333:     PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
334:     PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u));
335:     PetscCall(PetscObjectSetName((PetscObject)u, "Exact Solution"));
336:     PetscCall(VecViewFromOptions(u, NULL, "-exact_vec_view"));
337:   }
338:   PetscCall(DMSNESCheckFromOptions(snes, u));
339:   PetscCall(VecSet(u, 0.0));
340:   PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
341:   PetscCall(SNESSolve(snes, NULL, u));
342:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

344:   PetscCall(VecDestroy(&u));
345:   PetscCall(VecDestroy(&r));
346:   PetscCall(DMDestroy(&dm));
347:   PetscCall(SNESDestroy(&snes));
348:   PetscCall(PetscBagDestroy(&user.bag));
349:   PetscCall(PetscFinalize());
350:   return 0;
351: }

353: /*TEST

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