Actual source code: ex19.c
petsc-3.6.1 2015-08-06
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 <lid>, where <lid> = dimensionless velocity of lid\n\
7: -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
8: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
9: -contours : draw contour plots of solution\n\n" ;
10: /* in HTML, '<' = '<' and '>' = '>' */
12: /*
13: See src/ksp/ksp/examples/tutorials/ex45.c
14: */
16: /*T
17: Concepts: SNES ^solving a system of nonlinear equations (parallel multicomponent example);
18: Concepts: DMDA ^using distributed arrays;
19: Concepts: multicomponent
20: Processors: n
21: T*/
We thank David E. Keyes for contributing the driven cavity discretization within this example code.
This problem is modeled by the partial differential equation system
\begin{eqnarray}
- \triangle U - \nabla_y \Omega & = & 0 \\
- \triangle V + \nabla_x\Omega & = & 0 \\
- \triangle \Omega + \nabla \cdot ([U*\Omega,V*\Omega]) - GR* \nabla_x T & = & 0 \\
- \triangle T + PR* \nabla \cdot ([U*T,V*T]) & = & 0
\end{eqnarray}
in the unit square, which is uniformly discretized in each of x and y in this simple encoding.
No-slip, rigid-wall Dirichlet conditions are used for $ [U,V]$.
Dirichlet conditions are used for Omega, based on the definition of
vorticity: $ \Omega = - \nabla_y U + \nabla_x V$, where along each
constant coordinate boundary, the tangential derivative is zero.
Dirichlet conditions are used for T on the left and right walls,
and insulation homogeneous Neumann conditions are used for T on
the top and bottom walls.
A finite difference approximation with the usual 5-point stencil
is used to discretize the boundary value problem to obtain a
nonlinear system of equations. Upwinding is used for the divergence
(convective) terms and central for the gradient (source) terms.
The Jacobian can be either
* formed via finite differencing using coloring (the default), or
* applied matrix-free via the option -snes_mf
(for larger grid problems this variant may not converge
without a preconditioner due to ill-conditioning).
58: /*
59: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
60: Include "petscsnes.h" so that we can use SNES solvers. Note that this
61: file automatically includes:
62: petscsys.h - base PETSc routines petscvec.h - vectors
63: petscmat.h - matrices
64: petscis.h - index sets petscksp.h - Krylov subspace methods
65: petscviewer.h - viewers petscpc.h - preconditioners
66: petscksp.h - linear solvers
67: */
68: #if defined(PETSC_APPLE_FRAMEWORK)
69: #import <PETSc/petscsnes.h>
70: #import <PETSc/petscdmda.h>
71: #else
72: #include <petscsnes.h>
73: #include <petscdm.h>
74: #include <petscdmda.h>
75: #endif
77: /*
78: User-defined routines and data structures
79: */
80: typedef struct {
81: PetscScalar u,v,omega,temp;
82: } Field;
84: PetscErrorCode FormFunctionLocal(DMDALocalInfo *,Field**,Field**,void*) ;
86: typedef struct {
87: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
88: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
89: } AppCtx;
91: extern PetscErrorCode FormInitialGuess(AppCtx*,DM ,Vec ) ;
92: extern PetscErrorCode NonlinearGS(SNES ,Vec ,Vec ,void*) ;
96: int main(int argc,char **argv)
97: {
98: AppCtx user; /* user-defined work context */
99: PetscInt mx,my,its;
101: MPI_Comm comm;
102: SNES snes;
103: DM da;
104: Vec x;
106: PetscInitialize (&argc,&argv,(char*)0,help);if (ierr) return (1);
109: comm = PETSC_COMM_WORLD ;
110: SNESCreate (comm,&snes);
112: /*
113: Create distributed array object to manage parallel grid and vectors
114: for principal unknowns (x) and governing residuals (f)
115: */
116: DMDACreate2d (PETSC_COMM_WORLD ,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR ,-4,-4,PETSC_DECIDE ,PETSC_DECIDE ,4,1,0,0,&da);
117: SNESSetDM (snes,(DM )da);
118: SNESSetNGS (snes, NonlinearGS, (void*)&user);
120: DMDAGetInfo (da,0,&mx,&my,PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE ,
121: PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE ,PETSC_IGNORE );
122: /*
123: Problem parameters (velocity of lid, prandtl, and grashof numbers)
124: */
125: user.lidvelocity = 1.0/(mx*my);
126: user.prandtl = 1.0;
127: user.grashof = 1.0;
129: PetscOptionsGetReal (NULL,"-lidvelocity" ,&user.lidvelocity,NULL);
130: PetscOptionsGetReal (NULL,"-prandtl" ,&user.prandtl,NULL);
131: PetscOptionsGetReal (NULL,"-grashof" ,&user.grashof,NULL);
132: PetscOptionsHasName (NULL,"-contours" ,&user.draw_contours);
134: DMDASetFieldName (da,0,"x_velocity" );
135: DMDASetFieldName (da,1,"y_velocity" );
136: DMDASetFieldName (da,2,"Omega" );
137: DMDASetFieldName (da,3,"temperature" );
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Create user context, set problem data, create vector data structures.
141: Also, compute the initial guess.
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Create nonlinear solver context
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: DMSetApplicationContext (da,&user);
149: DMDASNESSetFunctionLocal (da,INSERT_VALUES ,(PetscErrorCode (*)(DMDALocalInfo *,void*,void*,void*))FormFunctionLocal,&user);
150: SNESSetFromOptions (snes);
151: PetscPrintf (comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n" ,(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Solve the nonlinear system
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: DMCreateGlobalVector (da,&x);
158: FormInitialGuess(&user,da,x);
160: SNESSolve (snes,NULL,x);
162: SNESGetIterationNumber (snes,&its);
163: PetscPrintf (comm,"Number of SNES iterations = %D\n" , its);
165: /*
166: Visualize solution
167: */
168: if (user.draw_contours) {
169: VecView (x,PETSC_VIEWER_DRAW_WORLD );
170: }
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Free work space. All PETSc objects should be destroyed when they
174: are no longer needed.
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: VecDestroy (&x);
177: DMDestroy (&da);
178: SNESDestroy (&snes);
179: PetscFinalize ();
180: return 0;
181: }
183: /* ------------------------------------------------------------------- */
188: /*
189: FormInitialGuess - Forms initial approximation.
191: Input Parameters:
192: user - user-defined application context
193: X - vector
195: Output Parameter:
196: X - vector
197: */
198: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
199: {
200: PetscInt i,j,mx,xs,ys,xm,ym;
202: PetscReal grashof,dx;
203: Field **x;
206: grashof = user->grashof;
208: DMDAGetInfo (da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
209: dx = 1.0/(mx-1);
211: /*
212: Get local grid boundaries (for 2-dimensional DMDA ):
213: xs, ys - starting grid indices (no ghost points)
214: xm, ym - widths of local grid (no ghost points)
215: */
216: DMDAGetCorners (da,&xs,&ys,NULL,&xm,&ym,NULL);
218: /*
219: Get a pointer to vector data.
220: - For default PETSc vectors, VecGetArray () returns a pointer to
221: the data array. Otherwise, the routine is implementation dependent.
222: - You MUST call VecRestoreArray () when you no longer need access to
223: the array.
224: */
225: DMDAVecGetArray (da,X,&x);
227: /*
228: Compute initial guess over the locally owned part of the grid
229: Initial condition is motionless fluid and equilibrium temperature
230: */
231: for (j=ys; j<ys+ym; j++) {
232: for (i=xs; i<xs+xm; i++) {
233: x[j][i].u = 0.0;
234: x[j][i].v = 0.0;
235: x[j][i].omega = 0.0;
236: x[j][i].temp = (grashof>0)*i*dx;
237: }
238: }
240: /*
241: Restore vector
242: */
243: DMDAVecRestoreArray (da,X,&x);
244: return (0);
245: }
249: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
250: {
251: AppCtx *user = (AppCtx*)ptr;
253: PetscInt xints,xinte,yints,yinte,i,j;
254: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
255: PetscReal grashof,prandtl,lid;
256: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
259: grashof = user->grashof;
260: prandtl = user->prandtl;
261: lid = user->lidvelocity;
263: /*
264: Define mesh intervals ratios for uniform grid.
266: Note: FD formulae below are normalized by multiplying through by
267: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
270: */
271: dhx = (PetscReal )(info->mx-1); dhy = (PetscReal )(info->my-1);
272: hx = 1.0/dhx; hy = 1.0/dhy;
273: hxdhy = hx*dhy; hydhx = hy*dhx;
275: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
277: /* Test whether we are on the bottom edge of the global array */
278: if (yints == 0) {
279: j = 0;
280: yints = yints + 1;
281: /* bottom edge */
282: for (i=info->xs; i<info->xs+info->xm; i++) {
283: f[j][i].u = x[j][i].u;
284: f[j][i].v = x[j][i].v;
285: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
286: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
287: }
288: }
290: /* Test whether we are on the top edge of the global array */
291: if (yinte == info->my) {
292: j = info->my - 1;
293: yinte = yinte - 1;
294: /* top edge */
295: for (i=info->xs; i<info->xs+info->xm; i++) {
296: f[j][i].u = x[j][i].u - lid;
297: f[j][i].v = x[j][i].v;
298: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
299: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
300: }
301: }
303: /* Test whether we are on the left edge of the global array */
304: if (xints == 0) {
305: i = 0;
306: xints = xints + 1;
307: /* left edge */
308: for (j=info->ys; j<info->ys+info->ym; j++) {
309: f[j][i].u = x[j][i].u;
310: f[j][i].v = x[j][i].v;
311: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
312: f[j][i].temp = x[j][i].temp;
313: }
314: }
316: /* Test whether we are on the right edge of the global array */
317: if (xinte == info->mx) {
318: i = info->mx - 1;
319: xinte = xinte - 1;
320: /* right edge */
321: for (j=info->ys; j<info->ys+info->ym; j++) {
322: f[j][i].u = x[j][i].u;
323: f[j][i].v = x[j][i].v;
324: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
325: f[j][i].temp = x[j][i].temp - (PetscReal )(grashof>0);
326: }
327: }
329: /* Compute over the interior points */
330: for (j=yints; j<yinte; j++) {
331: for (i=xints; i<xinte; i++) {
333: /*
334: convective coefficients for upwinding
335: */
336: vx = x[j][i].u; avx = PetscAbsScalar(vx);
337: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
338: vy = x[j][i].v; avy = PetscAbsScalar(vy);
339: vyp = .5*(vy+avy); vym = .5*(vy-avy);
341: /* U velocity */
342: u = x[j][i].u;
343: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
344: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
345: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
347: /* V velocity */
348: u = x[j][i].v;
349: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
350: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
351: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
353: /* Omega */
354: u = x[j][i].omega;
355: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
356: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
357: f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
358: (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
359: .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy;
361: /* Temperature */
362: u = x[j][i].temp;
363: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
364: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
365: f[j][i].temp = uxx + uyy + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
366: (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx);
367: }
368: }
370: /*
371: Flop count (multiply-adds are counted as 2 operations)
372: */
373: PetscLogFlops (84.0*info->ym*info->xm);
374: return (0);
375: }
380: PetscErrorCode NonlinearGS(SNES snes, Vec X, Vec B, void *ctx)
381: {
382: DMDALocalInfo info;
383: Field **x,**b;
385: Vec localX, localB;
386: DM da;
387: PetscInt xints,xinte,yints,yinte,i,j,k,l;
388: PetscInt max_its,tot_its;
389: PetscInt sweeps;
390: PetscReal rtol,atol,stol;
391: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
392: PetscReal grashof,prandtl,lid;
393: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
394: PetscScalar fu, fv, fomega, ftemp;
395: PetscScalar dfudu;
396: PetscScalar dfvdv;
397: PetscScalar dfodu, dfodv, dfodo;
398: PetscScalar dftdu, dftdv, dftdt;
399: PetscScalar yu=0, yv=0, yo=0, yt=0;
400: PetscScalar bjiu, bjiv, bjiomega, bjitemp;
401: PetscBool ptconverged;
402: PetscReal pfnorm,pfnorm0,pynorm,pxnorm;
403: AppCtx *user = (AppCtx*)ctx;
406: grashof = user->grashof;
407: prandtl = user->prandtl;
408: lid = user->lidvelocity;
409: tot_its = 0;
410: SNESNGSGetTolerances (snes,&rtol,&atol,&stol,&max_its);
411: SNESNGSGetSweeps (snes,&sweeps);
412: SNESGetDM (snes,(DM *)&da);
413: DMGetLocalVector (da,&localX);
414: if (B) {
415: DMGetLocalVector (da,&localB);
416: }
417: /*
418: Scatter ghost points to local vector, using the 2-step process
419: DMGlobalToLocalBegin (), DMGlobalToLocalEnd ().
420: */
421: DMGlobalToLocalBegin (da,X,INSERT_VALUES ,localX);
422: DMGlobalToLocalEnd (da,X,INSERT_VALUES ,localX);
423: if (B) {
424: DMGlobalToLocalBegin (da,B,INSERT_VALUES ,localB);
425: DMGlobalToLocalEnd (da,B,INSERT_VALUES ,localB);
426: }
427: DMDAGetLocalInfo (da,&info);
428: DMDAVecGetArray (da,localX,&x);
429: if (B) {
430: DMDAVecGetArrayRead (da,localB,&b);
431: }
432: /* looks like a combination of the formfunction / formjacobian routines */
433: dhx = (PetscReal )(info.mx-1);dhy = (PetscReal )(info.my-1);
434: hx = 1.0/dhx; hy = 1.0/dhy;
435: hxdhy = hx*dhy; hydhx = hy*dhx;
437: xints = info.xs; xinte = info.xs+info.xm; yints = info.ys; yinte = info.ys+info.ym;
439: /* Set the boundary conditions on the momentum equations */
440: /* Test whether we are on the bottom edge of the global array */
441: if (yints == 0) {
442: j = 0;
443: yints = yints + 1;
444: /* bottom edge */
445: for (i=info.xs; i<info.xs+info.xm; i++) {
447: if (B) {
448: bjiu = b[j][i].u;
449: bjiv = b[j][i].v;
450: } else {
451: bjiu = 0.0;
452: bjiv = 0.0;
453: }
454: x[j][i].u = 0.0 + bjiu;
455: x[j][i].v = 0.0 + bjiv;
456: }
457: }
459: /* Test whether we are on the top edge of the global array */
460: if (yinte == info.my) {
461: j = info.my - 1;
462: yinte = yinte - 1;
463: /* top edge */
464: for (i=info.xs; i<info.xs+info.xm; i++) {
465: if (B) {
466: bjiu = b[j][i].u;
467: bjiv = b[j][i].v;
468: } else {
469: bjiu = 0.0;
470: bjiv = 0.0;
471: }
472: x[j][i].u = lid + bjiu;
473: x[j][i].v = bjiv;
474: }
475: }
477: /* Test whether we are on the left edge of the global array */
478: if (xints == 0) {
479: i = 0;
480: xints = xints + 1;
481: /* left edge */
482: for (j=info.ys; j<info.ys+info.ym; j++) {
483: if (B) {
484: bjiu = b[j][i].u;
485: bjiv = b[j][i].v;
486: } else {
487: bjiu = 0.0;
488: bjiv = 0.0;
489: }
490: x[j][i].u = 0.0 + bjiu;
491: x[j][i].v = 0.0 + bjiv;
492: }
493: }
495: /* Test whether we are on the right edge of the global array */
496: if (xinte == info.mx) {
497: i = info.mx - 1;
498: xinte = xinte - 1;
499: /* right edge */
500: for (j=info.ys; j<info.ys+info.ym; j++) {
501: if (B) {
502: bjiu = b[j][i].u;
503: bjiv = b[j][i].v;
504: } else {
505: bjiu = 0.0;
506: bjiv = 0.0;
507: }
508: x[j][i].u = 0.0 + bjiu;
509: x[j][i].v = 0.0 + bjiv;
510: }
511: }
513: for (k=0; k < sweeps; k++) {
514: for (j=info.ys; j<info.ys + info.ym; j++) {
515: for (i=info.xs; i<info.xs + info.xm; i++) {
516: ptconverged = PETSC_FALSE ;
517: pfnorm0 = 0.0;
518: pfnorm = 0.0;
519: fu = 0.0;
520: fv = 0.0;
521: fomega = 0.0;
522: ftemp = 0.0;
523: for (l = 0; l < max_its && !ptconverged; l++) {
524: if (B) {
525: bjiu = b[j][i].u;
526: bjiv = b[j][i].v;
527: bjiomega = b[j][i].omega;
528: bjitemp = b[j][i].temp;
529: } else {
530: bjiu = 0.0;
531: bjiv = 0.0;
532: bjiomega = 0.0;
533: bjitemp = 0.0;
534: }
536: if (i != 0 && i != info.mx - 1 && j != 0 && j != info.my-1) {
537: /* U velocity */
538: u = x[j][i].u;
539: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
540: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
541: fu = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx - bjiu;
542: dfudu = 2.0*(hydhx + hxdhy);
543: /* V velocity */
544: u = x[j][i].v;
545: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
546: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
547: fv = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy - bjiv;
548: dfvdv = 2.0*(hydhx + hxdhy);
549: /*
550: convective coefficients for upwinding
551: */
552: vx = x[j][i].u; avx = PetscAbsScalar(vx);
553: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
554: vy = x[j][i].v; avy = PetscAbsScalar(vy);
555: vyp = .5*(vy+avy); vym = .5*(vy-avy);
556: /* Omega */
557: u = x[j][i].omega;
558: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
559: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
560: fomega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
561: (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
562: .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy - bjiomega;
563: /* convective coefficient derivatives */
564: dfodo = 2.0*(hydhx + hxdhy) + ((vxp - vxm)*hy + (vyp - vym)*hx);
565: if (PetscRealPart(vx) > 0.0) dfodu = (u - x[j][i-1].omega)*hy;
566: else dfodu = (x[j][i+1].omega - u)*hy;
568: if (PetscRealPart(vy) > 0.0) dfodv = (u - x[j-1][i].omega)*hx;
569: else dfodv = (x[j+1][i].omega - u)*hx;
571: /* Temperature */
572: u = x[j][i].temp;
573: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
574: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
575: ftemp = uxx + uyy + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
576: (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx) - bjitemp;
577: dftdt = 2.0*(hydhx + hxdhy) + prandtl*((vxp - vxm)*hy + (vyp - vym)*hx);
578: if (PetscRealPart(vx) > 0.0) dftdu = prandtl*(u - x[j][i-1].temp)*hy;
579: else dftdu = prandtl*(x[j][i+1].temp - u)*hy;
581: if (PetscRealPart(vy) > 0.0) dftdv = prandtl*(u - x[j-1][i].temp)*hx;
582: else dftdv = prandtl*(x[j+1][i].temp - u)*hx;
584: /* invert the system:
585: [ dfu / du 0 0 0 ][yu] = [fu]
586: [ 0 dfv / dv 0 0 ][yv] [fv]
587: [ dfo / du dfo / dv dfo / do 0 ][yo] [fo]
588: [ dft / du dft / dv 0 dft / dt ][yt] [ft]
589: by simple back-substitution
590: */
591: yu = fu / dfudu;
592: yv = fv / dfvdv;
593: yo = fomega / dfodo;
594: yt = ftemp / dftdt;
595: yo = (fomega - (dfodu*yu + dfodv*yv)) / dfodo;
596: yt = (ftemp - (dftdu*yu + dftdv*yv)) / dftdt;
598: x[j][i].u = x[j][i].u - yu;
599: x[j][i].v = x[j][i].v - yv;
600: x[j][i].temp = x[j][i].temp - yt;
601: x[j][i].omega = x[j][i].omega - yo;
602: }
603: if (i == 0) {
604: fomega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx - bjiomega;
605: ftemp = x[j][i].temp - bjitemp;
606: yo = fomega;
607: yt = ftemp;
608: x[j][i].omega = x[j][i].omega - fomega;
609: x[j][i].temp = x[j][i].temp - ftemp;
610: }
611: if (i == info.mx - 1) {
612: fomega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx - bjiomega;
613: ftemp = x[j][i].temp - (PetscReal )(grashof>0) - bjitemp;
614: yo = fomega;
615: yt = ftemp;
616: x[j][i].omega = x[j][i].omega - fomega;
617: x[j][i].temp = x[j][i].temp - ftemp;
618: }
619: if (j == 0) {
620: fomega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy - bjiomega;
621: ftemp = x[j][i].temp-x[j+1][i].temp - bjitemp;
622: yo = fomega;
623: yt = ftemp;
624: x[j][i].omega = x[j][i].omega - fomega;
625: x[j][i].temp = x[j][i].temp - ftemp;
626: }
627: if (j == info.my - 1) {
628: fomega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy - bjiomega;
629: ftemp = x[j][i].temp-x[j-1][i].temp - bjitemp;
630: yo = fomega;
631: yt = ftemp;
632: x[j][i].omega = x[j][i].omega - fomega;
633: x[j][i].temp = x[j][i].temp - ftemp;
634: }
635: tot_its++;
636: pfnorm = PetscRealPart(fu*fu + fv*fv + fomega*fomega + ftemp*ftemp);
637: pfnorm = PetscSqrtReal(pfnorm);
638: pynorm = PetscRealPart(yu*yu + yv*yv + yo*yo + yt*yt);
639: pfnorm = PetscSqrtReal(pynorm);
640: 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);
641: pxnorm = PetscSqrtReal(pxnorm);
642: if (l == 0) pfnorm0 = pfnorm;
643: if (rtol*pfnorm0 > pfnorm ||
644: atol > pfnorm ||
645: pxnorm*stol > pynorm) ptconverged = PETSC_TRUE ;
646: }
647: }
648: }
649: }
650: DMDAVecRestoreArray (da,localX,&x);
651: if (B) {
652: DMDAVecRestoreArrayRead (da,localB,&b);
653: }
654: DMLocalToGlobalBegin (da,localX,INSERT_VALUES ,X);
655: DMLocalToGlobalEnd (da,localX,INSERT_VALUES ,X);
656: PetscLogFlops (tot_its*(84.0 + 41.0 + 26.0));
657: DMRestoreLocalVector (da,&localX);
658: if (B) {
659: DMRestoreLocalVector (da,&localB);
660: }
661: return (0);
662: }