xref: /petsc/src/snes/tutorials/ex71.c (revision 8ebe3e4e9e00d86ece2e9fcd0cc84910b0ad437c)
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";
4 
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$.
15 
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.
18 
19 For visualization, use
20 
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*/
23 
24 #include <petscdmplex.h>
25 #include <petscsnes.h>
26 #include <petscds.h>
27 #include <petscbag.h>
28 
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;
35 
36 typedef struct {
37   PetscBag bag;    /* Holds problem parameters */
38 } AppCtx;
39 
40 /*
41   In 2D, plane Poiseuille flow has exact solution:
42 
43     u = \Delta/(2 \nu) y (1 - y) + u_0
44     v = 0
45     p = -\Delta x
46     f = 0
47 
48   so that
49 
50     -\nu \Delta u + \nabla p + f = <\Delta, 0> + <-\Delta, 0> + <0, 0> = 0
51     \nabla \cdot u               = 0 + 0                               = 0
52 
53   In 3D we use exact solution:
54 
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
60 
61   so that
62 
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
65 
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;
76 
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 0;
81 }
82 
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;
87 
88   p[0] = -Delta * X[0];
89   return 0;
90 }
91 
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;
97 
98   u[0] = u_0;
99   for (d = 1; d < dim; ++d) u[d] = 0.0;
100   return 0;
101 }
102 
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,
106           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
107           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
108           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
109 {
110   const PetscReal nu = PetscRealPart(constants[1]);
111   const PetscInt  Nc = dim;
112   PetscInt        c, d;
113 
114   for (c = 0; c < Nc; ++c) {
115     for (d = 0; d < dim; ++d) {
116       /* f1[c*dim+d] = 0.5*nu*(u_x[c*dim+d] + u_x[d*dim+c]); */
117       f1[c*dim+d] = nu*u_x[c*dim+d];
118     }
119     f1[c*dim+c] -= u[uOff[1]];
120   }
121 }
122 
123 /* 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} */
124 void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
125           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
126           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
127           PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
128 {
129   PetscInt d;
130   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d];
131 }
132 
133 /* Residual functions are in reference coordinates */
134 static void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
135                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
136                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
137                     PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
138 {
139   const PetscReal Delta = PetscRealPart(constants[0]);
140   PetscReal       alpha = PetscRealPart(constants[3]);
141   PetscReal       X     = PetscCosReal(alpha)*x[0] + PetscSinReal(alpha)*x[1];
142   PetscInt        d;
143 
144   for (d = 0; d < dim; ++d) {
145     f0[d] = -Delta * X * n[d];
146   }
147 }
148 
149 /* < q, \nabla\cdot u >
150    NcompI = 1, NcompJ = dim */
151 void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
152            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
153            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
154            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
155 {
156   PetscInt d;
157   for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
158 }
159 
160 /* -< \nabla\cdot v, p >
161     NcompI = dim, NcompJ = 1 */
162 void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux,
163            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
164            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
165            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
166 {
167   PetscInt d;
168   for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
169 }
170 
171 /* < \nabla v, \nabla u + {\nabla u}^T >
172    This just gives \nabla u, give the perdiagonal for the transpose */
173 void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
174            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
175            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
176            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
177 {
178   const PetscReal nu = PetscRealPart(constants[1]);
179   const PetscInt  Nc = dim;
180   PetscInt        c, d;
181 
182   for (c = 0; c < Nc; ++c) {
183     for (d = 0; d < dim; ++d) {
184       g3[((c*Nc+c)*dim+d)*dim+d] = nu;
185     }
186   }
187 }
188 
189 static PetscErrorCode SetupParameters(AppCtx *user)
190 {
191   PetscBag       bag;
192   Parameter     *p;
193   PetscErrorCode ierr;
194 
195   PetscFunctionBeginUser;
196   /* setup PETSc parameter bag */
197   ierr = PetscBagGetData(user->bag, (void **) &p);CHKERRQ(ierr);
198   ierr = PetscBagSetName(user->bag, "par", "Poiseuille flow parameters");CHKERRQ(ierr);
199   bag  = user->bag;
200   ierr = PetscBagRegisterReal(bag, &p->Delta, 1.0, "Delta", "Pressure drop per unit length");CHKERRQ(ierr);
201   ierr = PetscBagRegisterReal(bag, &p->nu,    1.0, "nu",    "Kinematic viscosity");CHKERRQ(ierr);
202   ierr = PetscBagRegisterReal(bag, &p->u_0,   0.0, "u_0",   "Tangential velocity at the wall");CHKERRQ(ierr);
203   ierr = PetscBagRegisterReal(bag, &p->alpha, 0.0, "alpha", "Angle of pipe wall to x-axis");CHKERRQ(ierr);
204   PetscFunctionReturn(0);
205 }
206 
207 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
208 {
209   PetscErrorCode ierr;
210 
211   PetscFunctionBeginUser;
212   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
213   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
214   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
215   {
216     Parameter   *param;
217     Vec          coordinates;
218     PetscScalar *coords;
219     PetscReal    alpha;
220     PetscInt     cdim, N, bs, i;
221 
222     ierr = DMGetCoordinateDim(*dm, &cdim);CHKERRQ(ierr);
223     ierr = DMGetCoordinates(*dm, &coordinates);CHKERRQ(ierr);
224     ierr = VecGetLocalSize(coordinates, &N);CHKERRQ(ierr);
225     ierr = VecGetBlockSize(coordinates, &bs);CHKERRQ(ierr);
226     if (bs != cdim) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %D != embedding dimension %D", bs, cdim);
227     ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
228     ierr = PetscBagGetData(user->bag, (void **) &param);CHKERRQ(ierr);
229     alpha = param->alpha;
230     for (i = 0; i < N; i += cdim) {
231       PetscScalar x = coords[i+0];
232       PetscScalar y = coords[i+1];
233 
234       coords[i+0] = PetscCosReal(alpha)*x - PetscSinReal(alpha)*y;
235       coords[i+1] = PetscSinReal(alpha)*x + PetscCosReal(alpha)*y;
236     }
237     ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
238     ierr = DMSetCoordinates(*dm, coordinates);CHKERRQ(ierr);
239   }
240   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
241   PetscFunctionReturn(0);
242 }
243 
244 PetscErrorCode SetupProblem(DM dm, AppCtx *user)
245 {
246   PetscDS        ds;
247   PetscWeakForm  wf;
248   DMLabel        label;
249   Parameter     *ctx;
250   PetscInt       id, bd;
251   PetscErrorCode ierr;
252 
253   PetscFunctionBeginUser;
254   ierr = PetscBagGetData(user->bag, (void **) &ctx);CHKERRQ(ierr);
255   ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
256   ierr = PetscDSSetResidual(ds, 0, NULL, f1_u);CHKERRQ(ierr);
257   ierr = PetscDSSetResidual(ds, 1, f0_p, NULL);CHKERRQ(ierr);
258   ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL,  NULL,  g3_uu);CHKERRQ(ierr);
259   ierr = PetscDSSetJacobian(ds, 0, 1, NULL, NULL,  g2_up, NULL);CHKERRQ(ierr);
260   ierr = PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL,  NULL);CHKERRQ(ierr);
261 
262   id   = 2;
263   ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr);
264   ierr = DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd);CHKERRQ(ierr);
265   ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
266   ierr = PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_bd_u, 0, NULL);CHKERRQ(ierr);
267   /* Setup constants */
268   {
269     Parameter  *param;
270     PetscScalar constants[4];
271 
272     ierr = PetscBagGetData(user->bag, (void **) &param);CHKERRQ(ierr);
273 
274     constants[0] = param->Delta;
275     constants[1] = param->nu;
276     constants[2] = param->u_0;
277     constants[3] = param->alpha;
278     ierr = PetscDSSetConstants(ds, 4, constants);CHKERRQ(ierr);
279   }
280   /* Setup Boundary Conditions */
281   id   = 3;
282   ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall",    label, 1, &id, 0, 0, NULL, (void (*)(void)) wall_velocity, NULL, ctx, NULL);CHKERRQ(ierr);
283   id   = 1;
284   ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) wall_velocity, NULL, ctx, NULL);CHKERRQ(ierr);
285   /* Setup exact solution */
286   ierr = PetscDSSetExactSolution(ds, 0, quadratic_u, ctx);CHKERRQ(ierr);
287   ierr = PetscDSSetExactSolution(ds, 1, linear_p, ctx);CHKERRQ(ierr);
288   PetscFunctionReturn(0);
289 }
290 
291 PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
292 {
293   DM             cdm = dm;
294   PetscFE        fe[2];
295   Parameter     *param;
296   PetscBool      simplex;
297   PetscInt       dim;
298   MPI_Comm       comm;
299   PetscErrorCode ierr;
300 
301   PetscFunctionBeginUser;
302   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
303   ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr);
304   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
305   ierr = PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr);
306   ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr);
307   ierr = PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr);
308   ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr);
309   ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr);
310   /* Set discretization and boundary conditions for each mesh */
311   ierr = DMSetField(dm, 0, NULL, (PetscObject) fe[0]);CHKERRQ(ierr);
312   ierr = DMSetField(dm, 1, NULL, (PetscObject) fe[1]);CHKERRQ(ierr);
313   ierr = DMCreateDS(dm);CHKERRQ(ierr);
314   ierr = SetupProblem(dm, user);CHKERRQ(ierr);
315   ierr = PetscBagGetData(user->bag, (void **) &param);CHKERRQ(ierr);
316   while (cdm) {
317     ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
318     ierr = DMPlexCreateBasisRotation(cdm, param->alpha, 0.0, 0.0);CHKERRQ(ierr);
319     ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
320   }
321   ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr);
322   ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr);
323   PetscFunctionReturn(0);
324 }
325 
326 int main(int argc, char **argv)
327 {
328   SNES           snes; /* nonlinear solver */
329   DM             dm;   /* problem definition */
330   Vec            u, r; /* solution and residual */
331   AppCtx         user; /* user-defined work context */
332   PetscErrorCode ierr;
333 
334   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
335   ierr = PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);CHKERRQ(ierr);
336   ierr = SetupParameters(&user);CHKERRQ(ierr);
337   ierr = PetscBagSetFromOptions(user.bag);CHKERRQ(ierr);
338   ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
339   ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
340   ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
341   ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);
342   /* Setup problem */
343   ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr);
344   ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr);
345 
346   ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
347   ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
348 
349   ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr);
350 
351   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
352 
353   {
354     PetscDS              ds;
355     PetscSimplePointFunc exactFuncs[2];
356     void                *ctxs[2];
357 
358     ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
359     ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);CHKERRQ(ierr);
360     ierr = PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);CHKERRQ(ierr);
361     ierr = DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
362     ierr = PetscObjectSetName((PetscObject) u, "Exact Solution");CHKERRQ(ierr);
363     ierr = VecViewFromOptions(u, NULL, "-exact_vec_view");CHKERRQ(ierr);
364   }
365   ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr);
366   ierr = VecSet(u, 0.0);CHKERRQ(ierr);
367   ierr = PetscObjectSetName((PetscObject) u, "Solution");CHKERRQ(ierr);
368   ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
369   ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr);
370 
371   ierr = VecDestroy(&u);CHKERRQ(ierr);
372   ierr = VecDestroy(&r);CHKERRQ(ierr);
373   ierr = DMDestroy(&dm);CHKERRQ(ierr);
374   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
375   ierr = PetscBagDestroy(&user.bag);CHKERRQ(ierr);
376   ierr = PetscFinalize();
377   return ierr;
378 }
379 
380 /*TEST
381 
382   # Convergence
383   test:
384     suffix: 2d_quad_q1_p0_conv
385     requires: !single
386     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 \
387       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
388       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
389       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
390       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
391         -fieldsplit_velocity_pc_type lu \
392         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
393   test:
394     suffix: 2d_quad_q1_p0_conv_u0
395     requires: !single
396     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 \
397       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
398       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
399       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
400       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
401         -fieldsplit_velocity_pc_type lu \
402         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
403   test:
404     suffix: 2d_quad_q1_p0_conv_u0_alpha
405     requires: !single
406     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 -alpha 0.3927 \
407       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
408       -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
409       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
410       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
411         -fieldsplit_velocity_pc_type lu \
412         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
413   test:
414     suffix: 2d_quad_q1_p0_conv_gmg_vanka
415     requires: !single long_runtime
416     args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
417       -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
418       -snes_convergence_estimate -convest_num_refine 1 -snes_error_if_not_converged \
419       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
420       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
421         -fieldsplit_velocity_pc_type mg \
422           -fieldsplit_velocity_mg_levels_pc_type patch -fieldsplit_velocity_mg_levels_pc_patch_exclude_subspaces 1 \
423           -fieldsplit_velocity_mg_levels_pc_patch_construct_codim 0 -fieldsplit_velocity_mg_levels_pc_patch_construct_type vanka \
424         -fieldsplit_pressure_ksp_rtol 1e-5 -fieldsplit_pressure_pc_type jacobi
425   test:
426     suffix: 2d_tri_p2_p1_conv
427     requires: triangle !single
428     args: -dm_plex_separate_marker -dm_refine 1 \
429       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
430       -dmsnes_check .001 -snes_error_if_not_converged \
431       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
432       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
433         -fieldsplit_velocity_pc_type lu \
434         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
435   test:
436     suffix: 2d_tri_p2_p1_conv_u0_alpha
437     requires: triangle !single
438     args: -dm_plex_separate_marker -dm_refine 0 -u_0 0.125 -alpha 0.3927 \
439       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
440       -dmsnes_check .001 -snes_error_if_not_converged \
441       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
442       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
443         -fieldsplit_velocity_pc_type lu \
444         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
445   test:
446     suffix: 2d_tri_p2_p1_conv_gmg_vcycle
447     requires: triangle !single
448     args: -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
449       -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
450       -dmsnes_check .001 -snes_error_if_not_converged \
451       -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
452       -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
453         -fieldsplit_velocity_pc_type mg \
454         -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
455 TEST*/
456