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