xref: /petsc/src/snes/tutorials/ex19.c (revision 607e733f3db3ee7f6f605a13295c517df8dbb9c9)
1 /* Portions of this code are under:
2    Copyright (C) 2022 Advanced Micro Devices, Inc. All rights reserved.
3 */
4 static char help[] = "Nonlinear driven cavity with multigrid in 2d.\n \
5   \n\
6 The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
7 The flow can be driven with the lid or with buoyancy or both:\n\
8   -lidvelocity &ltlid&gt, where &ltlid&gt = dimensionless velocity of lid\n\
9   -grashof &ltgr&gt, where &ltgr&gt = dimensionless temperature gradent\n\
10   -prandtl &ltpr&gt, where &ltpr&gt = dimensionless thermal/momentum diffusity ratio\n\
11  -contours : draw contour plots of solution\n\n";
12 /* in HTML, '&lt' = '<' and '&gt' = '>' */
13 
14 /*
15       See src/ksp/ksp/tutorials/ex45.c
16 */
17 
18 /*F-----------------------------------------------------------------------
19 
20     We thank David E. Keyes for contributing the driven cavity discretization within this example code.
21 
22     This problem is modeled by the partial differential equation system
23 
24 \begin{eqnarray}
25         - \triangle U - \nabla_y \Omega & = & 0  \\
26         - \triangle V + \nabla_x\Omega & = & 0  \\
27         - \triangle \Omega + \nabla \cdot ([U*\Omega,V*\Omega]) - GR* \nabla_x T & = & 0  \\
28         - \triangle T + PR* \nabla \cdot ([U*T,V*T]) & = & 0
29 \end{eqnarray}
30 
31     in the unit square, which is uniformly discretized in each of x and y in this simple encoding.
32 
33     No-slip, rigid-wall Dirichlet conditions are used for $ [U,V]$.
34     Dirichlet conditions are used for Omega, based on the definition of
35     vorticity: $ \Omega = - \nabla_y U + \nabla_x V$, where along each
36     constant coordinate boundary, the tangential derivative is zero.
37     Dirichlet conditions are used for T on the left and right walls,
38     and insulation homogeneous Neumann conditions are used for T on
39     the top and bottom walls.
40 
41     A finite difference approximation with the usual 5-point stencil
42     is used to discretize the boundary value problem to obtain a
43     nonlinear system of equations.  Upwinding is used for the divergence
44     (convective) terms and central for the gradient (source) terms.
45 
46     The Jacobian can be either
47       * formed via finite differencing using coloring (the default), or
48       * applied matrix-free via the option -snes_mf
49         (for larger grid problems this variant may not converge
50         without a preconditioner due to ill-conditioning).
51 
52   ------------------------------------------------------------------------F*/
53 
54 /*
55    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
56    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
57    file automatically includes:
58      petscsys.h       - base PETSc routines   petscvec.h - vectors
59      petscmat.h - matrices
60      petscis.h     - index sets            petscksp.h - Krylov subspace methods
61      petscviewer.h - viewers               petscpc.h  - preconditioners
62      petscksp.h   - linear solvers
63 */
64 #if defined(PETSC_APPLE_FRAMEWORK)
65   #import <PETSc/petscsnes.h>
66   #import <PETSc/petscdmda.h>
67 #else
68   #include <petscsnes.h>
69   #include <petscdm.h>
70   #include <petscdmda.h>
71 #endif
72 
73 /*
74    User-defined routines and data structures
75 */
76 typedef struct {
77   PetscScalar u, v, omega, temp;
78 } Field;
79 
80 PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field **, Field **, void *);
81 
82 typedef struct {
83   PetscReal lidvelocity, prandtl, grashof; /* physical parameters */
84   PetscBool draw_contours;                 /* flag - 1 indicates drawing contours */
85 } AppCtx;
86 
87 extern PetscErrorCode FormInitialGuess(AppCtx *, DM, Vec);
88 extern PetscErrorCode NonlinearGS(SNES, Vec, Vec, void *);
89 
90 int main(int argc, char **argv)
91 {
92   AppCtx   user; /* user-defined work context */
93   PetscInt mx, my, its;
94   MPI_Comm comm;
95   SNES     snes;
96   DM       da;
97   Vec      x;
98 
99   PetscFunctionBeginUser;
100   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
101   comm = PETSC_COMM_WORLD;
102   PetscCall(SNESCreate(comm, &snes));
103 
104   /*
105       Create distributed array object to manage parallel grid and vectors
106       for principal unknowns (x) and governing residuals (f)
107   */
108   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da));
109   PetscCall(DMSetFromOptions(da));
110   PetscCall(DMSetUp(da));
111   PetscCall(SNESSetDM(snes, da));
112   PetscCall(SNESSetNGS(snes, NonlinearGS, (void *)&user));
113 
114   PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
115   /*
116      Problem parameters (velocity of lid, prandtl, and grashof numbers)
117   */
118   user.lidvelocity = 1.0 / (mx * my);
119   user.prandtl     = 1.0;
120   user.grashof     = 1.0;
121 
122   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lidvelocity", &user.lidvelocity, NULL));
123   PetscCall(PetscOptionsGetReal(NULL, NULL, "-prandtl", &user.prandtl, NULL));
124   PetscCall(PetscOptionsGetReal(NULL, NULL, "-grashof", &user.grashof, NULL));
125   PetscCall(PetscOptionsHasName(NULL, NULL, "-contours", &user.draw_contours));
126 
127   PetscCall(DMDASetFieldName(da, 0, "x_velocity"));
128   PetscCall(DMDASetFieldName(da, 1, "y_velocity"));
129   PetscCall(DMDASetFieldName(da, 2, "Omega"));
130   PetscCall(DMDASetFieldName(da, 3, "temperature"));
131 
132   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133      Create user context, set problem data, create vector data structures.
134      Also, compute the initial guess.
135      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 
137   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138      Create nonlinear solver context
139 
140      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141   PetscCall(DMSetApplicationContext(da, &user));
142   PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
143   PetscCall(SNESSetFromOptions(snes));
144   PetscCall(PetscPrintf(comm, "lid velocity = %g, prandtl # = %g, grashof # = %g\n", (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof));
145 
146   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147      Solve the nonlinear system
148      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149   PetscCall(DMCreateGlobalVector(da, &x));
150   PetscCall(FormInitialGuess(&user, da, x));
151 
152   PetscCall(SNESSolve(snes, NULL, x));
153 
154   PetscCall(SNESGetIterationNumber(snes, &its));
155   PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
156 
157   /*
158      Visualize solution
159   */
160   if (user.draw_contours) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
161 
162   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163      Free work space.  All PETSc objects should be destroyed when they
164      are no longer needed.
165      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166   PetscCall(VecDestroy(&x));
167   PetscCall(DMDestroy(&da));
168   PetscCall(SNESDestroy(&snes));
169   PetscCall(PetscFinalize());
170   return 0;
171 }
172 
173 /* ------------------------------------------------------------------- */
174 
175 /*
176    FormInitialGuess - Forms initial approximation.
177 
178    Input Parameters:
179    user - user-defined application context
180    X - vector
181 
182    Output Parameter:
183    X - vector
184 */
185 PetscErrorCode FormInitialGuess(AppCtx *user, DM da, Vec X)
186 {
187   PetscInt  i, j, mx, xs, ys, xm, ym;
188   PetscReal grashof, dx;
189   Field   **x;
190 
191   PetscFunctionBeginUser;
192   grashof = user->grashof;
193 
194   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
195   dx = 1.0 / (mx - 1);
196 
197   /*
198      Get local grid boundaries (for 2-dimensional DMDA):
199        xs, ys   - starting grid indices (no ghost points)
200        xm, ym   - widths of local grid (no ghost points)
201   */
202   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
203 
204   /*
205      Get a pointer to vector data.
206        - For default PETSc vectors, VecGetArray() returns a pointer to
207          the data array.  Otherwise, the routine is implementation dependent.
208        - You MUST call VecRestoreArray() when you no longer need access to
209          the array.
210   */
211   PetscCall(DMDAVecGetArrayWrite(da, X, &x));
212 
213   /*
214      Compute initial guess over the locally owned part of the grid
215      Initial condition is motionless fluid and equilibrium temperature
216   */
217   for (j = ys; j < ys + ym; j++) {
218     for (i = xs; i < xs + xm; i++) {
219       x[j][i].u     = 0.0;
220       x[j][i].v     = 0.0;
221       x[j][i].omega = 0.0;
222       x[j][i].temp  = (grashof > 0) * i * dx;
223     }
224   }
225 
226   /*
227      Restore vector
228   */
229   PetscCall(DMDAVecRestoreArrayWrite(da, X, &x));
230   PetscFunctionReturn(PETSC_SUCCESS);
231 }
232 
233 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr)
234 {
235   AppCtx     *user = (AppCtx *)ptr;
236   PetscInt    xints, xinte, yints, yinte, i, j;
237   PetscReal   hx, hy, dhx, dhy, hxdhy, hydhx;
238   PetscReal   grashof, prandtl, lid;
239   PetscScalar u, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;
240 
241   PetscFunctionBeginUser;
242   grashof = user->grashof;
243   prandtl = user->prandtl;
244   lid     = user->lidvelocity;
245 
246   /*
247      Define mesh intervals ratios for uniform grid.
248 
249      Note: FD formulae below are normalized by multiplying through by
250      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
251 
252   */
253   dhx   = (PetscReal)(info->mx - 1);
254   dhy   = (PetscReal)(info->my - 1);
255   hx    = 1.0 / dhx;
256   hy    = 1.0 / dhy;
257   hxdhy = hx * dhy;
258   hydhx = hy * dhx;
259 
260   xints = info->xs;
261   xinte = info->xs + info->xm;
262   yints = info->ys;
263   yinte = info->ys + info->ym;
264 
265   /* Test whether we are on the bottom edge of the global array */
266   if (yints == 0) {
267     j     = 0;
268     yints = yints + 1;
269     /* bottom edge */
270     for (i = info->xs; i < info->xs + info->xm; i++) {
271       f[j][i].u     = x[j][i].u;
272       f[j][i].v     = x[j][i].v;
273       f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy;
274       f[j][i].temp  = x[j][i].temp - x[j + 1][i].temp;
275     }
276   }
277 
278   /* Test whether we are on the top edge of the global array */
279   if (yinte == info->my) {
280     j     = info->my - 1;
281     yinte = yinte - 1;
282     /* top edge */
283     for (i = info->xs; i < info->xs + info->xm; i++) {
284       f[j][i].u     = x[j][i].u - lid;
285       f[j][i].v     = x[j][i].v;
286       f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy;
287       f[j][i].temp  = x[j][i].temp - x[j - 1][i].temp;
288     }
289   }
290 
291   /* Test whether we are on the left edge of the global array */
292   if (xints == 0) {
293     i     = 0;
294     xints = xints + 1;
295     /* left edge */
296     for (j = info->ys; j < info->ys + info->ym; j++) {
297       f[j][i].u     = x[j][i].u;
298       f[j][i].v     = x[j][i].v;
299       f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx;
300       f[j][i].temp  = x[j][i].temp;
301     }
302   }
303 
304   /* Test whether we are on the right edge of the global array */
305   if (xinte == info->mx) {
306     i     = info->mx - 1;
307     xinte = xinte - 1;
308     /* right edge */
309     for (j = info->ys; j < info->ys + info->ym; j++) {
310       f[j][i].u     = x[j][i].u;
311       f[j][i].v     = x[j][i].v;
312       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx;
313       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof > 0);
314     }
315   }
316 
317   /* Compute over the interior points */
318   for (j = yints; j < yinte; j++) {
319     for (i = xints; i < xinte; i++) {
320       /*
321        convective coefficients for upwinding
322       */
323       vx  = x[j][i].u;
324       avx = PetscAbsScalar(vx);
325       vxp = .5 * (vx + avx);
326       vxm = .5 * (vx - avx);
327       vy  = x[j][i].v;
328       avy = PetscAbsScalar(vy);
329       vyp = .5 * (vy + avy);
330       vym = .5 * (vy - avy);
331 
332       /* U velocity */
333       u         = x[j][i].u;
334       uxx       = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
335       uyy       = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
336       f[j][i].u = uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;
337 
338       /* V velocity */
339       u         = x[j][i].v;
340       uxx       = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
341       uyy       = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
342       f[j][i].v = uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;
343 
344       /* Omega */
345       u             = x[j][i].omega;
346       uxx           = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
347       uyy           = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
348       f[j][i].omega = uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy;
349 
350       /* Temperature */
351       u            = x[j][i].temp;
352       uxx          = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
353       uyy          = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
354       f[j][i].temp = uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx);
355     }
356   }
357 
358   /*
359      Flop count (multiply-adds are counted as 2 operations)
360   */
361   PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));
362   PetscFunctionReturn(PETSC_SUCCESS);
363 }
364 
365 /*
366     Performs sweeps of point block nonlinear Gauss-Seidel on all the local grid points
367 */
368 PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, PetscCtx ctx)
369 {
370   DMDALocalInfo info;
371   Field       **x, **b;
372   Vec           localX, localB;
373   DM            da;
374   PetscInt      xints, xinte, yints, yinte, i, j, k, l;
375   PetscInt      max_its, tot_its;
376   PetscInt      sweeps;
377   PetscReal     rtol, atol, stol;
378   PetscReal     hx, hy, dhx, dhy, hxdhy, hydhx;
379   PetscReal     grashof, prandtl, lid;
380   PetscScalar   u, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;
381   PetscScalar   fu, fv, fomega, ftemp;
382   PetscScalar   dfudu;
383   PetscScalar   dfvdv;
384   PetscScalar   dfodu, dfodv, dfodo;
385   PetscScalar   dftdu, dftdv, dftdt;
386   PetscScalar   yu = 0, yv = 0, yo = 0, yt = 0;
387   PetscScalar   bjiu, bjiv, bjiomega, bjitemp;
388   PetscBool     ptconverged;
389   PetscReal     pfnorm, pfnorm0, pynorm, pxnorm;
390   AppCtx       *user = (AppCtx *)ctx;
391 
392   PetscFunctionBeginUser;
393   grashof = user->grashof;
394   prandtl = user->prandtl;
395   lid     = user->lidvelocity;
396   tot_its = 0;
397   PetscCall(SNESNGSGetTolerances(snes, &rtol, &atol, &stol, &max_its));
398   PetscCall(SNESNGSGetSweeps(snes, &sweeps));
399   PetscCall(SNESGetDM(snes, &da));
400   PetscCall(DMGetLocalVector(da, &localX));
401   if (B) PetscCall(DMGetLocalVector(da, &localB));
402   /*
403      Scatter ghost points to local vector, using the 2-step process
404         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
405   */
406   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
407   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
408   if (B) {
409     PetscCall(DMGlobalToLocalBegin(da, B, INSERT_VALUES, localB));
410     PetscCall(DMGlobalToLocalEnd(da, B, INSERT_VALUES, localB));
411   }
412   PetscCall(DMDAGetLocalInfo(da, &info));
413   PetscCall(DMDAVecGetArrayWrite(da, localX, &x));
414   if (B) PetscCall(DMDAVecGetArrayRead(da, localB, &b));
415   /* looks like a combination of the formfunction / formjacobian routines */
416   dhx   = (PetscReal)(info.mx - 1);
417   dhy   = (PetscReal)(info.my - 1);
418   hx    = 1.0 / dhx;
419   hy    = 1.0 / dhy;
420   hxdhy = hx * dhy;
421   hydhx = hy * dhx;
422 
423   xints = info.xs;
424   xinte = info.xs + info.xm;
425   yints = info.ys;
426   yinte = info.ys + info.ym;
427 
428   /* Set the boundary conditions on the momentum equations */
429   /* Test whether we are on the bottom edge of the global array */
430   if (yints == 0) {
431     j = 0;
432     /* bottom edge */
433     for (i = info.xs; i < info.xs + info.xm; i++) {
434       if (B) {
435         bjiu = b[j][i].u;
436         bjiv = b[j][i].v;
437       } else {
438         bjiu = 0.0;
439         bjiv = 0.0;
440       }
441       x[j][i].u = 0.0 + bjiu;
442       x[j][i].v = 0.0 + bjiv;
443     }
444   }
445 
446   /* Test whether we are on the top edge of the global array */
447   if (yinte == info.my) {
448     j = info.my - 1;
449     /* top edge */
450     for (i = info.xs; i < info.xs + info.xm; i++) {
451       if (B) {
452         bjiu = b[j][i].u;
453         bjiv = b[j][i].v;
454       } else {
455         bjiu = 0.0;
456         bjiv = 0.0;
457       }
458       x[j][i].u = lid + bjiu;
459       x[j][i].v = bjiv;
460     }
461   }
462 
463   /* Test whether we are on the left edge of the global array */
464   if (xints == 0) {
465     i = 0;
466     /* left edge */
467     for (j = info.ys; j < info.ys + info.ym; j++) {
468       if (B) {
469         bjiu = b[j][i].u;
470         bjiv = b[j][i].v;
471       } else {
472         bjiu = 0.0;
473         bjiv = 0.0;
474       }
475       x[j][i].u = 0.0 + bjiu;
476       x[j][i].v = 0.0 + bjiv;
477     }
478   }
479 
480   /* Test whether we are on the right edge of the global array */
481   if (xinte == info.mx) {
482     i = info.mx - 1;
483     /* right edge */
484     for (j = info.ys; j < info.ys + info.ym; j++) {
485       if (B) {
486         bjiu = b[j][i].u;
487         bjiv = b[j][i].v;
488       } else {
489         bjiu = 0.0;
490         bjiv = 0.0;
491       }
492       x[j][i].u = 0.0 + bjiu;
493       x[j][i].v = 0.0 + bjiv;
494     }
495   }
496 
497   for (k = 0; k < sweeps; k++) {
498     for (j = info.ys; j < info.ys + info.ym; j++) {
499       for (i = info.xs; i < info.xs + info.xm; i++) {
500         ptconverged = PETSC_FALSE;
501         pfnorm0     = 0.0;
502         fu          = 0.0;
503         fv          = 0.0;
504         fomega      = 0.0;
505         ftemp       = 0.0;
506         /*  Run Newton's method on a single grid point */
507         for (l = 0; l < max_its && !ptconverged; l++) {
508           if (B) {
509             bjiu     = b[j][i].u;
510             bjiv     = b[j][i].v;
511             bjiomega = b[j][i].omega;
512             bjitemp  = b[j][i].temp;
513           } else {
514             bjiu     = 0.0;
515             bjiv     = 0.0;
516             bjiomega = 0.0;
517             bjitemp  = 0.0;
518           }
519 
520           if (i != 0 && i != info.mx - 1 && j != 0 && j != info.my - 1) {
521             /* U velocity */
522             u     = x[j][i].u;
523             uxx   = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
524             uyy   = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
525             fu    = uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx - bjiu;
526             dfudu = 2.0 * (hydhx + hxdhy);
527             /* V velocity */
528             u     = x[j][i].v;
529             uxx   = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
530             uyy   = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
531             fv    = uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy - bjiv;
532             dfvdv = 2.0 * (hydhx + hxdhy);
533             /*
534              convective coefficients for upwinding
535              */
536             vx  = x[j][i].u;
537             avx = PetscAbsScalar(vx);
538             vxp = .5 * (vx + avx);
539             vxm = .5 * (vx - avx);
540             vy  = x[j][i].v;
541             avy = PetscAbsScalar(vy);
542             vyp = .5 * (vy + avy);
543             vym = .5 * (vy - avy);
544             /* Omega */
545             u      = x[j][i].omega;
546             uxx    = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
547             uyy    = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
548             fomega = uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy - bjiomega;
549             /* convective coefficient derivatives */
550             dfodo = 2.0 * (hydhx + hxdhy) + ((vxp - vxm) * hy + (vyp - vym) * hx);
551             if (PetscRealPart(vx) > 0.0) dfodu = (u - x[j][i - 1].omega) * hy;
552             else dfodu = (x[j][i + 1].omega - u) * hy;
553 
554             if (PetscRealPart(vy) > 0.0) dfodv = (u - x[j - 1][i].omega) * hx;
555             else dfodv = (x[j + 1][i].omega - u) * hx;
556 
557             /* Temperature */
558             u     = x[j][i].temp;
559             uxx   = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
560             uyy   = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
561             ftemp = uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx) - bjitemp;
562             dftdt = 2.0 * (hydhx + hxdhy) + prandtl * ((vxp - vxm) * hy + (vyp - vym) * hx);
563             if (PetscRealPart(vx) > 0.0) dftdu = prandtl * (u - x[j][i - 1].temp) * hy;
564             else dftdu = prandtl * (x[j][i + 1].temp - u) * hy;
565 
566             if (PetscRealPart(vy) > 0.0) dftdv = prandtl * (u - x[j - 1][i].temp) * hx;
567             else dftdv = prandtl * (x[j + 1][i].temp - u) * hx;
568 
569             /* invert the system:
570              [ dfu / du     0        0        0    ][yu] = [fu]
571              [     0    dfv / dv     0        0    ][yv]   [fv]
572              [ dfo / du dfo / dv dfo / do     0    ][yo]   [fo]
573              [ dft / du dft / dv     0    dft / dt ][yt]   [ft]
574              by simple back-substitution
575            */
576             yu = fu / dfudu;
577             yv = fv / dfvdv;
578             yo = (fomega - (dfodu * yu + dfodv * yv)) / dfodo;
579             yt = (ftemp - (dftdu * yu + dftdv * yv)) / dftdt;
580 
581             x[j][i].u     = x[j][i].u - yu;
582             x[j][i].v     = x[j][i].v - yv;
583             x[j][i].temp  = x[j][i].temp - yt;
584             x[j][i].omega = x[j][i].omega - yo;
585           }
586           if (i == 0) {
587             fomega        = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx - bjiomega;
588             ftemp         = x[j][i].temp - bjitemp;
589             yo            = fomega;
590             yt            = ftemp;
591             x[j][i].omega = x[j][i].omega - fomega;
592             x[j][i].temp  = x[j][i].temp - ftemp;
593           }
594           if (i == info.mx - 1) {
595             fomega        = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx - bjiomega;
596             ftemp         = x[j][i].temp - (PetscReal)(grashof > 0) - bjitemp;
597             yo            = fomega;
598             yt            = ftemp;
599             x[j][i].omega = x[j][i].omega - fomega;
600             x[j][i].temp  = x[j][i].temp - ftemp;
601           }
602           if (j == 0) {
603             fomega        = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy - bjiomega;
604             ftemp         = x[j][i].temp - x[j + 1][i].temp - bjitemp;
605             yo            = fomega;
606             yt            = ftemp;
607             x[j][i].omega = x[j][i].omega - fomega;
608             x[j][i].temp  = x[j][i].temp - ftemp;
609           }
610           if (j == info.my - 1) {
611             fomega        = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy - bjiomega;
612             ftemp         = x[j][i].temp - x[j - 1][i].temp - bjitemp;
613             yo            = fomega;
614             yt            = ftemp;
615             x[j][i].omega = x[j][i].omega - fomega;
616             x[j][i].temp  = x[j][i].temp - ftemp;
617           }
618           tot_its++;
619           pfnorm = PetscRealPart(fu * fu + fv * fv + fomega * fomega + ftemp * ftemp);
620           pfnorm = PetscSqrtReal(pfnorm);
621           pynorm = PetscRealPart(yu * yu + yv * yv + yo * yo + yt * yt);
622           pynorm = PetscSqrtReal(pynorm);
623           pxnorm = PetscRealPart(x[j][i].u * x[j][i].u + x[j][i].v * x[j][i].v + x[j][i].omega * x[j][i].omega + x[j][i].temp * x[j][i].temp);
624           pxnorm = PetscSqrtReal(pxnorm);
625           if (l == 0) pfnorm0 = pfnorm;
626           if (rtol * pfnorm0 > pfnorm || atol > pfnorm || pxnorm * stol > pynorm) ptconverged = PETSC_TRUE;
627         }
628       }
629     }
630   }
631   PetscCall(DMDAVecRestoreArrayWrite(da, localX, &x));
632   if (B) PetscCall(DMDAVecRestoreArrayRead(da, localB, &b));
633   PetscCall(DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X));
634   PetscCall(DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X));
635   PetscCall(PetscLogFlops(tot_its * (84.0 + 41.0 + 26.0)));
636   PetscCall(DMRestoreLocalVector(da, &localX));
637   if (B) PetscCall(DMRestoreLocalVector(da, &localB));
638   PetscFunctionReturn(PETSC_SUCCESS);
639 }
640 
641 /*TEST
642 
643    test:
644       nsize: 2
645       args: -da_refine 3 -snes_monitor_short -pc_type mg -ksp_type fgmres -pc_mg_type full
646       requires: !single
647 
648    test:
649       suffix: 10
650       nsize: 3
651       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -pc_fieldsplit_type symmetric_multiplicative -snes_view -da_refine 1 -ksp_type fgmres
652       requires: !single
653 
654    test:
655       suffix: 11
656       nsize: 4
657       requires: pastix
658       args: -snes_monitor_short -pc_type redundant -dm_mat_type mpiaij -redundant_pc_factor_mat_solver_type pastix -mat_pastix_thread_nbr 1 -pc_redundant_number 2 -da_refine 4 -ksp_type fgmres
659 
660    test:
661       suffix: 12
662       nsize: 12
663       requires: pastix
664       args: -snes_monitor_short -pc_type redundant -dm_mat_type mpiaij -redundant_pc_factor_mat_solver_type pastix -mat_pastix_thread_nbr 1 -pc_redundant_number 5 -da_refine 4 -ksp_type fgmres
665 
666    test:
667       suffix: 13
668       nsize: 3
669       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -pc_fieldsplit_type multiplicative -snes_view -da_refine 1 -ksp_type fgmres -snes_mf_operator
670       requires: !single
671 
672    test:
673       suffix: 14
674       nsize: 4
675       args: -snes_monitor_short -pc_type mg -dm_mat_type baij -mg_coarse_pc_type bjacobi -da_refine 3 -ksp_type fgmres
676       requires: !single
677 
678    test:
679       suffix: 14_ds
680       nsize: 4
681       args: -snes_converged_reason -pc_type mg -dm_mat_type baij -mg_coarse_pc_type bjacobi -da_refine 3 -ksp_type fgmres -mat_fd_type ds
682       output_file: output/ex19_2.out
683       requires: !single
684 
685    test:
686       suffix: 17
687       args: -snes_monitor_short -ksp_pc_side right
688       requires: !single
689 
690    test:
691       suffix: 18
692       args: -snes_monitor_ksp draw::draw_lg -ksp_pc_side right
693       requires: x !single
694 
695    test:
696       suffix: 19
697       nsize: 2
698       args: -da_refine 3 -snes_monitor_short -pc_type mg -ksp_type fgmres -pc_mg_type full -snes_type newtontrdc
699       requires: !single
700 
701    test:
702       suffix: 20
703       nsize: 2
704       args: -da_refine 3 -snes_monitor_short -pc_type mg -ksp_type fgmres -pc_mg_type full -snes_type newtontrdc -snes_trdc_use_cauchy false
705       requires: !single
706 
707    test:
708       suffix: 2
709       nsize: 4
710       args: -da_refine 3 -snes_converged_reason -pc_type mg -mat_fd_type ds
711       requires: !single
712 
713    test:
714       suffix: 2_bcols1
715       nsize: 4
716       args: -da_refine 3 -snes_converged_reason -pc_type mg -mat_fd_type ds -mat_fd_coloring_bcols
717       output_file: output/ex19_2.out
718       requires: !single
719 
720    test:
721       suffix: 3
722       nsize: 4
723       requires: mumps
724       args: -da_refine 3 -snes_monitor_short -pc_type redundant -dm_mat_type mpiaij -redundant_ksp_type preonly -redundant_pc_factor_mat_solver_type mumps -pc_redundant_number 2
725 
726    test:
727       suffix: 4
728       nsize: 12
729       requires: mumps
730       args: -da_refine 3 -snes_monitor_short -pc_type redundant -dm_mat_type mpiaij -redundant_ksp_type preonly -redundant_pc_factor_mat_solver_type mumps -pc_redundant_number 5
731       output_file: output/ex19_3.out
732 
733    test:
734       suffix: 6
735       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -snes_view -ksp_type fgmres -da_refine 1
736       requires: !single
737 
738    test:
739       suffix: 7
740       nsize: 3
741       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -snes_view -da_refine 1 -ksp_type fgmres
742 
743       requires: !single
744    test:
745       suffix: 8
746       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -pc_fieldsplit_block_size 2 -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 0,1 -pc_fieldsplit_type multiplicative -snes_view -fieldsplit_pc_type lu -da_refine 1 -ksp_type fgmres
747       requires: !single
748 
749    test:
750       suffix: 9
751       nsize: 3
752       args: -snes_monitor_short -ksp_monitor_short -pc_type fieldsplit -pc_fieldsplit_type multiplicative -snes_view -da_refine 1 -ksp_type fgmres
753       requires: !single
754 
755    test:
756       suffix: aspin
757       nsize: 4
758       args: -da_refine 3 -da_overlap 2 -snes_monitor_short -snes_type aspin -grashof 4e4 -lidvelocity 100 -ksp_monitor_short -npc_sub_ksp_type preonly -npc_sub_pc_type lu
759       requires: !single
760 
761    test:
762       suffix: bcgsl
763       nsize: 2
764       args: -ksp_type bcgsl -ksp_monitor_short -da_refine 2 -ksp_bcgsl_ell 3 -snes_view
765       requires: !single
766 
767    test:
768       suffix: bcols1
769       nsize: 2
770       args: -da_refine 3 -snes_monitor_short -pc_type mg -ksp_type fgmres -pc_mg_type full -mat_fd_coloring_bcols 1
771       output_file: output/ex19_1.out
772       requires: !single
773 
774    test:
775       suffix: bjacobi
776       nsize: 4
777       args: -da_refine 4 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 2 -sub_ksp_type gmres -sub_ksp_max_it 2 -sub_pc_type bjacobi -sub_sub_ksp_type preonly -sub_sub_pc_type ilu -snes_monitor_short
778       requires: !single
779 
780    test:
781       suffix: cgne
782       args: -da_refine 2 -pc_type lu -ksp_type cgne -ksp_monitor_short -ksp_converged_reason -ksp_view -ksp_norm_type unpreconditioned
783       filter: grep -v HERMITIAN
784       requires: !single
785 
786    test:
787       suffix: cgs
788       args: -da_refine 1 -ksp_monitor_short -ksp_type cgs
789       requires: !single
790 
791    test:
792       suffix: composite_fieldsplit
793       args: -ksp_type fgmres -pc_type composite -pc_composite_type MULTIPLICATIVE -pc_composite_pcs fieldsplit,none -sub_0_pc_fieldsplit_block_size 4 -sub_0_pc_fieldsplit_type additive -sub_0_pc_fieldsplit_0_fields 0,1,2 -sub_0_pc_fieldsplit_1_fields 3 -snes_monitor_short -ksp_monitor_short
794       requires: !single
795 
796    test:
797       suffix: composite_fieldsplit_bjacobi
798       args: -ksp_type fgmres -pc_type composite -pc_composite_type MULTIPLICATIVE -pc_composite_pcs fieldsplit,bjacobi -sub_0_pc_fieldsplit_block_size 4 -sub_0_pc_fieldsplit_type additive -sub_0_pc_fieldsplit_0_fields 0,1,2 -sub_0_pc_fieldsplit_1_fields 3 -sub_1_pc_bjacobi_blocks 16 -sub_1_sub_pc_type lu -snes_monitor_short -ksp_monitor_short
799       requires: !single
800 
801    test:
802       suffix: composite_fieldsplit_bjacobi_2
803       nsize: 4
804       args: -ksp_type fgmres -pc_type composite -pc_composite_type MULTIPLICATIVE -pc_composite_pcs fieldsplit,bjacobi -sub_0_pc_fieldsplit_block_size 4 -sub_0_pc_fieldsplit_type additive -sub_0_pc_fieldsplit_0_fields 0,1,2 -sub_0_pc_fieldsplit_1_fields 3 -sub_1_pc_bjacobi_blocks 16 -sub_1_sub_pc_type lu -snes_monitor_short -ksp_monitor_short
805       requires: !single
806 
807    test:
808       suffix: composite_gs_newton
809       nsize: 2
810       args: -da_refine 3 -grashof 4e4 -lidvelocity 100 -snes_monitor_short -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses ngs,newtonls -sub_0_snes_max_it 20 -sub_1_pc_type mg
811       requires: !single
812 
813    test:
814       suffix: cuda
815       requires: cuda !single
816       args: -dm_vec_type cuda -dm_mat_type aijcusparse -pc_type none -ksp_type fgmres -snes_monitor_short -snes_rtol 1.e-5
817 
818    test:
819       suffix: hip
820       requires: hip !single
821       args: -dm_vec_type hip -dm_mat_type aijhipsparse -pc_type none -ksp_type fgmres -snes_monitor_short -snes_rtol 1.e-5
822 
823    test:
824       suffix: draw
825       args: -pc_type fieldsplit -snes_view draw -fieldsplit_x_velocity_pc_type mg -fieldsplit_x_velocity_pc_mg_galerkin pmat -fieldsplit_x_velocity_pc_mg_levels 2 -da_refine 1 -fieldsplit_x_velocity_mg_coarse_pc_type svd
826       requires: x !single
827 
828    test:
829       suffix: drawports
830       args: -snes_monitor_solution draw::draw_ports -da_refine 1
831       output_file: output/ex19_draw.out
832       requires: x !single
833 
834    test:
835       suffix: fas
836       args: -da_refine 4 -snes_monitor_short -snes_type fas -fas_levels_snes_type ngs -fas_levels_snes_ngs_sweeps 3 -fas_levels_snes_ngs_atol 0.0 -fas_levels_snes_ngs_stol 0.0 -grashof 4e4 -snes_fas_smoothup 6 -snes_fas_smoothdown 6 -lidvelocity 100
837       requires: !single
838 
839    test:
840       suffix: fas_full
841       args: -da_refine 4 -snes_monitor_short -snes_type fas -snes_fas_type full -snes_fas_full_downsweep -fas_levels_snes_type ngs -fas_levels_snes_ngs_sweeps 3 -fas_levels_snes_ngs_atol 0.0 -fas_levels_snes_ngs_stol 0.0 -grashof 4e4 -snes_fas_smoothup 6 -snes_fas_smoothdown 6 -lidvelocity 100
842       requires: !single
843 
844    test:
845       suffix: fdcoloring_ds
846       args: -da_refine 3 -snes_converged_reason -pc_type mg -mat_fd_type ds
847       output_file: output/ex19_2.out
848       requires: !single
849 
850    test:
851       suffix: fdcoloring_ds_baij
852       args: -da_refine 3 -snes_converged_reason -pc_type mg -mat_fd_type ds -dm_mat_type baij
853       output_file: output/ex19_2.out
854       requires: !single
855 
856    test:
857       suffix: fdcoloring_ds_bcols1
858       args: -da_refine 3 -snes_converged_reason -pc_type mg -mat_fd_type ds -mat_fd_coloring_bcols 1
859       output_file: output/ex19_2.out
860       requires: !single
861 
862    test:
863       suffix: fdcoloring_wp
864       args: -da_refine 3 -snes_monitor_short -pc_type mg
865       requires: !single
866 
867    test:
868       suffix: fdcoloring_wp_baij
869       args: -da_refine 3 -snes_monitor_short -pc_type mg -dm_mat_type baij
870       output_file: output/ex19_fdcoloring_wp.out
871       requires: !single
872 
873    test:
874       suffix: fdcoloring_wp_bcols1
875       args: -da_refine 3 -snes_monitor_short -pc_type mg -mat_fd_coloring_bcols 1
876       output_file: output/ex19_fdcoloring_wp.out
877       requires: !single
878 
879    test:
880       suffix: fieldsplit_2
881       args: -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_block_size 4 -pc_fieldsplit_type additive -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -snes_monitor_short -ksp_monitor_short
882       requires: !single
883 
884    test:
885       suffix: fieldsplit_3
886       args: -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_block_size 4 -pc_fieldsplit_type additive -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_1_pc_type lu -snes_monitor_short -ksp_monitor_short
887       requires: !single
888 
889    test:
890       suffix: fieldsplit_4
891       args: -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_block_size 4 -pc_fieldsplit_type SCHUR -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_1_pc_type lu -snes_monitor_short -ksp_monitor_short
892       requires: !single
893 
894    # HYPRE PtAP broken with complex numbers
895    test:
896       suffix: fieldsplit_hypre
897       nsize: 2
898       requires: hypre mumps !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
899       args: -pc_type fieldsplit -pc_fieldsplit_block_size 4 -pc_fieldsplit_type SCHUR -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_0_pc_factor_mat_solver_type mumps -fieldsplit_1_pc_type hypre -fieldsplit_1_pc_hypre_type boomeramg -snes_monitor_short -ksp_monitor_short
900 
901    test:
902       suffix: fieldsplit_mumps
903       nsize: 2
904       requires: mumps
905       args: -pc_type fieldsplit -pc_fieldsplit_block_size 4 -pc_fieldsplit_type SCHUR -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_1_pc_type lu -snes_monitor_short -ksp_monitor_short -fieldsplit_0_pc_factor_mat_solver_type mumps -fieldsplit_1_pc_factor_mat_solver_type mumps
906       output_file: output/ex19_fieldsplit_5.out
907 
908    test:
909       suffix: greedy_coloring
910       nsize: 2
911       args: -da_refine 3 -snes_monitor_short -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_weight_type lf -mat_coloring_view> ex19_greedy_coloring.tmp 2>&1
912       requires: !single
913 
914    # HYPRE PtAP broken with complex numbers
915    test:
916       suffix: hypre
917       nsize: 2
918       requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
919       args: -da_refine 3 -snes_monitor_short -pc_type hypre -ksp_norm_type unpreconditioned
920 
921    # ibcgs is broken when using device vectors
922    test:
923       suffix: ibcgs
924       nsize: 2
925       args: -ksp_type ibcgs -ksp_monitor_short -da_refine 2 -snes_view
926       requires: !complex !single
927 
928    test:
929       suffix: kaczmarz
930       nsize: 2
931       args: -pc_type kaczmarz -ksp_monitor_short -snes_monitor_short -snes_view
932       requires: !single
933 
934    test:
935       suffix: klu
936       requires: suitesparse
937       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type klu
938       output_file: output/ex19_superlu.out
939 
940    test:
941       suffix: klu_2
942       requires: suitesparse
943       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type klu -pc_factor_mat_ordering_type nd
944       output_file: output/ex19_superlu.out
945 
946    test:
947       suffix: klu_3
948       requires: suitesparse
949       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type klu -mat_klu_use_btf 0
950       output_file: output/ex19_superlu.out
951 
952    test:
953       suffix: ml
954       nsize: 2
955       requires: ml
956       args: -da_refine 3 -snes_monitor_short -pc_type ml
957 
958    test:
959       suffix: ngmres_fas
960       args: -da_refine 4 -snes_monitor_short -snes_type ngmres -npc_fas_levels_snes_type ngs -npc_fas_levels_snes_ngs_sweeps 3 -npc_fas_levels_snes_ngs_atol 0.0 -npc_fas_levels_snes_ngs_stol 0.0 -npc_snes_type fas -npc_fas_levels_snes_type ngs -npc_snes_max_it 1 -npc_snes_fas_smoothup 6 -npc_snes_fas_smoothdown 6 -lidvelocity 100 -grashof 4e4
961       requires: !single
962 
963    test:
964       suffix: ngmres_fas_gssecant
965       args: -da_refine 3 -snes_monitor_short -snes_type ngmres -npc_snes_type fas -npc_fas_levels_snes_type ngs -npc_fas_levels_snes_max_it 6 -npc_fas_levels_snes_ngs_secant -npc_fas_levels_snes_ngs_max_it 1 -npc_fas_coarse_snes_max_it 1 -lidvelocity 100 -grashof 4e4
966       requires: !single
967 
968    test:
969       suffix: ngmres_fas_ms
970       nsize: 2
971       args: -snes_grid_sequence 2 -lidvelocity 200 -grashof 1e4 -snes_monitor_short -snes_view -snes_converged_reason -snes_type ngmres -npc_snes_type fas -npc_fas_coarse_snes_type newtonls -npc_fas_coarse_ksp_type preonly -npc_snes_max_it 1
972       requires: !single
973 
974    test:
975       suffix: ngmres_nasm
976       nsize: 4
977       args: -da_refine 4 -da_overlap 2 -snes_monitor_short -snes_type ngmres -snes_max_it 10 -npc_snes_type nasm -npc_snes_nasm_type basic -grashof 4e4 -lidvelocity 100
978       requires: !single
979 
980    test:
981       suffix: ngs
982       args: -snes_type ngs -snes_view -snes_monitor -snes_rtol 1e-4
983       requires: !single
984 
985    test:
986       suffix: ngs_fd
987       args: -snes_type ngs -snes_ngs_secant -snes_view -snes_monitor -snes_rtol 1e-4
988       requires: !single
989 
990    test:
991       suffix: parms
992       nsize: 2
993       requires: parms
994       args: -pc_type parms -ksp_monitor_short -snes_view
995 
996    test:
997       suffix: superlu
998       requires: superlu
999       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu
1000 
1001    test:
1002       suffix: superlu_sell
1003       requires: superlu
1004       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu -dm_mat_type sell -pc_factor_mat_ordering_type natural
1005       output_file: output/ex19_superlu.out
1006 
1007    test:
1008       suffix: superlu_dist
1009       requires: superlu_dist
1010       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu_dist
1011       output_file: output/ex19_superlu.out
1012 
1013    test:
1014       suffix: superlu_dist_2
1015       nsize: 2
1016       requires: superlu_dist
1017       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu_dist
1018       output_file: output/ex19_superlu.out
1019 
1020    test:
1021       suffix: superlu_dist_3d
1022       nsize: 4
1023       requires: superlu_dist !defined(PETSCTEST_VALGRIND)
1024       filter: grep -v iam | grep -v openMP
1025       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_3d -mat_superlu_dist_d 2 -snes_view -snes_monitor -ksp_monitor
1026 
1027    test:
1028       suffix: superlu_dist_2s
1029       nsize: 2
1030       requires: superlu_dist defined(PETSC_HAVE_SUPERLU_DIST_SINGLE)
1031       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu_dist -pc_precision single
1032       output_file: output/ex19_superlu.out
1033 
1034    test:
1035       suffix: mumps_mixed
1036       nsize: 2
1037       requires: mumps defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
1038       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type mumps -pc_precision {{single double}}
1039       output_file: output/ex19_superlu.out
1040 
1041    test:
1042       suffix: superlu_dist_3ds
1043       nsize: 4
1044       requires: superlu_dist !defined(PETSCTEST_VALGRIND) defined(PETSC_HAVE_SUPERLU_DIST_SINGLE)
1045       filter: grep -v iam | grep -v openMP
1046       args: -da_grid_x 20 -da_grid_y 20 -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_3d -mat_superlu_dist_d 2 -snes_view -snes_monitor -ksp_monitor -pc_precision single
1047 
1048    test:
1049       suffix: superlu_equil
1050       requires: superlu
1051       args: -da_grid_x 20 -da_grid_y 20 -{snes,ksp}_monitor_short -pc_type lu -pc_factor_mat_solver_type superlu -mat_superlu_equil
1052 
1053    test:
1054       suffix: superlu_equil_sell
1055       requires: superlu
1056       args: -da_grid_x 20 -da_grid_y 20 -{snes,ksp}_monitor_short -pc_type lu -pc_factor_mat_solver_type superlu -mat_superlu_equil -dm_mat_type sell -pc_factor_mat_ordering_type natural
1057       output_file: output/ex19_superlu_equil.out
1058 
1059    test:
1060       suffix: tcqmr
1061       args: -da_refine 1 -ksp_monitor_short -ksp_type tcqmr
1062       requires: !single
1063 
1064    test:
1065       suffix: tfqmr
1066       args: -da_refine 1 -ksp_monitor_short -ksp_type tfqmr
1067       requires: !single
1068 
1069    test:
1070       suffix: umfpack
1071       requires: suitesparse
1072       args: -da_refine 2 -pc_type lu -pc_factor_mat_solver_type umfpack -snes_view -snes_monitor_short -ksp_monitor_short -pc_factor_mat_ordering_type external
1073 
1074    test:
1075       suffix: tut_1
1076       nsize: 4
1077       requires: !single
1078       args: -da_refine 5 -snes_monitor -ksp_monitor -snes_view
1079 
1080    test:
1081       suffix: tut_2
1082       nsize: 4
1083       requires: !single
1084       args: -da_refine 5 -snes_monitor -ksp_monitor -snes_view -pc_type mg
1085 
1086    # HYPRE PtAP broken with complex numbers
1087    test:
1088       suffix: tut_3
1089       nsize: 4
1090       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
1091       args: -da_refine 5 -snes_monitor -snes_converged_reason -pc_type hypre -dm_mat_type {{aij baij}}
1092 
1093    test:
1094       suffix: tut_3_seq
1095       nsize: 1
1096       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
1097       args: -da_refine 1 -snes_monitor -snes_converged_reason -pc_type hypre -dm_mat_type {{seqaij mpiaij seqbaij mpibaij}}
1098 
1099    test:
1100       suffix: tut_8
1101       nsize: 4
1102       requires: ml !single
1103       args: -da_refine 5 -snes_monitor -ksp_monitor -snes_view -pc_type ml
1104 
1105    test:
1106       suffix: tut_4
1107       nsize: 1
1108       requires: !single
1109       args: -da_refine 5 -log_view
1110       filter: head -n 2
1111       filter_output: head -n 2
1112 
1113    test:
1114       suffix: tut_5
1115       nsize: 1
1116       requires: !single
1117       args: -da_refine 5 -log_view -pc_type mg
1118       filter: head -n 2
1119       filter_output: head -n 2
1120 
1121    test:
1122       suffix: tut_6
1123       nsize: 4
1124       requires: !single
1125       args: -da_refine 5 -log_view
1126       filter: head -n 2
1127       filter_output: head -n 2
1128 
1129    test:
1130       suffix: tut_7
1131       nsize: 4
1132       requires: !single
1133       args: -da_refine 5 -log_view -pc_type mg
1134       filter: head -n 2
1135       filter_output: head -n 2
1136 
1137    test:
1138       suffix: cuda_1
1139       nsize: 1
1140       requires: cuda
1141       args: -snes_monitor -dm_mat_type seqaijcusparse -dm_vec_type seqcuda -pc_type gamg -ksp_monitor -mg_levels_ksp_max_it 1
1142 
1143    test:
1144       suffix: cuda_2
1145       nsize: 3
1146       requires: cuda !single
1147       args: -snes_monitor -dm_mat_type mpiaijcusparse -dm_vec_type mpicuda -pc_type gamg -ksp_monitor -mg_levels_ksp_max_it 1
1148 
1149    test:
1150       suffix: cuda_dm_bind_below
1151       nsize: 2
1152       requires: cuda defined(PETSC_USE_LOG)
1153       args: -dm_mat_type aijcusparse -dm_vec_type cuda -da_refine 3 -pc_type mg -mg_levels_ksp_type chebyshev -mg_levels_pc_type jacobi -log_view -pc_mg_log -dm_bind_below 10000
1154       filter: awk "/Level/ {print \$NF}"
1155 
1156    test:
1157       suffix: hip_1
1158       nsize: 1
1159       requires: hip
1160       args: -snes_monitor -dm_mat_type mpiaijhipsparse -dm_vec_type hip -pc_type gamg -ksp_monitor -mg_levels_ksp_max_it 1
1161 
1162    test:
1163       suffix: hip_2
1164       nsize: 3
1165       requires: hip !single
1166       args: -snes_monitor -dm_mat_type mpiaijhipsparse -dm_vec_type mpihip -pc_type gamg -ksp_monitor -mg_levels_ksp_max_it 1
1167 
1168    test:
1169       suffix: hip_dm_bind_below
1170       nsize: 2
1171       requires: hip
1172       args: -dm_mat_type aijhipsparse -dm_vec_type hip -da_refine 3 -pc_type mg -mg_levels_ksp_type chebyshev -mg_levels_pc_type jacobi -log_view -pc_mg_log -dm_bind_below 10000
1173       filter: awk "/Level/ {print \$NF}"
1174 
1175    test:
1176       suffix: viennacl_dm_bind_below
1177       nsize: 2
1178       requires: viennacl
1179       args: -dm_mat_type aijviennacl -dm_vec_type viennacl -da_refine 3 -pc_type mg -mg_levels_ksp_type chebyshev -mg_levels_pc_type jacobi -log_view -pc_mg_log -dm_bind_below 10000
1180       filter: awk "/Level/ {print \$NF}"
1181 
1182    test:
1183       suffix: seqbaijmkl
1184       nsize: 1
1185       requires: defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
1186       args: -dm_mat_type baij -snes_monitor -ksp_monitor -snes_view
1187 
1188    test:
1189       suffix: mpibaijmkl
1190       nsize: 2
1191       requires: defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
1192       args: -dm_mat_type baij -snes_monitor -ksp_monitor -snes_view
1193 
1194    test:
1195      suffix: cpardiso
1196      nsize: 4
1197      requires: mkl_cpardiso
1198      args: -pc_type lu -pc_factor_mat_solver_type mkl_cpardiso -ksp_monitor
1199 
1200    test:
1201      suffix: logviewmemory
1202      requires: defined(PETSC_USE_LOG) !defined(PETSC_HAVE_THREADSAFETY)
1203      args: -log_view -log_view_memory -da_refine 4
1204      filter: grep MatFDColorSetUp | wc -w | xargs -I % sh -c "expr % \> 21"
1205 
1206    test:
1207      suffix: fs
1208      requires: !single
1209      args: -pc_type fieldsplit -da_refine 3 -all_ksp_monitor -fieldsplit_y_velocity_pc_type lu -fieldsplit_temperature_pc_type lu -fieldsplit_x_velocity_pc_type lu -snes_view
1210 
1211    test:
1212      suffix: asm_matconvert
1213      args: -mat_type aij -pc_type asm -pc_asm_sub_mat_type dense -snes_view
1214 
1215    test:
1216       suffix: euclid
1217       nsize: 2
1218       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_MIXEDINT) !defined(PETSC_HAVE_HYPRE_DEVICE)
1219       args: -da_refine 2 -ksp_monitor -snes_monitor -snes_view -pc_type hypre -pc_hypre_type euclid
1220 
1221    test:
1222       suffix: euclid_bj
1223       nsize: 2
1224       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_MIXEDINT) !defined(PETSC_HAVE_HYPRE_DEVICE)
1225       args: -da_refine 2 -ksp_monitor -snes_monitor -snes_view -pc_type hypre -pc_hypre_type euclid -pc_hypre_euclid_bj
1226 
1227    test:
1228       suffix: euclid_droptolerance
1229       nsize: 1
1230       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_MIXEDINT) !defined(PETSC_HAVE_HYPRE_DEVICE)
1231       args: -da_refine 2 -ksp_monitor -snes_monitor -snes_view -pc_type hypre -pc_hypre_type euclid -pc_hypre_euclid_droptolerance .1
1232 
1233    test:
1234       suffix: failure_size
1235       nsize: 1
1236       requires: !defined(PETSC_USE_64BIT_INDICES) !defined(PETSCTEST_VALGRIND) !defined(PETSC_HAVE_SANITIZER)
1237       args: -da_refine 100 -petsc_ci_portable_error_output -error_output_stdout
1238       filter: grep -E -v "(memory block|leaked context|not freed before MPI_Finalize|Could be the program crashed)"
1239 
1240    testset:
1241       requires: hpddm cuda
1242       args: -snes_monitor -ksp_converged_reason -ksp_type hpddm -pc_type jacobi -dm_mat_type aijcusparse -dm_vec_type cuda
1243       test:
1244         suffix: hpddm_cuda
1245         filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 15/Linear solve converged due to CONVERGED_RTOL iterations 14/g"
1246         args: -ksp_hpddm_type {{gmres gcrodr}separate output} -ksp_hpddm_precision {{single double}shared output}
1247       test:
1248         suffix: hpddm_cuda_right
1249         filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 15/Linear solve converged due to CONVERGED_RTOL iterations 14/g"
1250         args: -ksp_hpddm_type gcrodr -ksp_pc_side right
1251         output_file: output/ex19_hpddm_cuda_ksp_hpddm_type-gcrodr.out
1252 
1253 TEST*/
1254