Actual source code: ex7.c
1: /* Program usage: mpiexec -n <procs> ex7 [-help] [all PETSc options] */
3: static char help[] = "Nonlinear PDE in 2d.\n\
4: We solve the Stokes equation in a 2D rectangular\n\
5: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";
7: /*T
8: Concepts: SNES^parallel Stokes example
9: Concepts: DMDA^using distributed arrays;
10: Processors: n
11: T*/
13: /* ------------------------------------------------------------------------
15: The Stokes equation is given by the partial differential equation
16:
17: -alpha*Laplacian u - \nabla p = f, 0 < x,y < 1,
19: \nabla \cdot u = 0
20:
21: with boundary conditions
22:
23: u = 0 for x = 0, x = 1, y = 0, y = 1.
24:
25: A P2/P1 finite element approximation is used to discretize the boundary
26: value problem on the two triangles which make up each rectangle in the DMDA
27: to obtain a nonlinear system of equations.
29: ------------------------------------------------------------------------- */
31: /*
32: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
33: Include "petscsnes.h" so that we can use SNES solvers. Note that this
34: file automatically includes:
35: petscsys.h - base PETSc routines petscvec.h - vectors
36: petscmat.h - matrices
37: petscis.h - index sets petscksp.h - Krylov subspace methods
38: petscviewer.h - viewers petscpc.h - preconditioners
39: petscksp.h - linear solvers
40: */
41: #include <petscsys.h>
42: #include <petscbag.h>
43: #include <petscdmda.h>
44: #include <petscdmmg.h>
45: #include <petscsnes.h>
47: /*
48: User-defined application context - contains data needed by the
49: application-provided call-back routines, FormJacobianLocal() and
50: FormFunctionLocal().
51: */
52: typedef struct {
53: PetscReal alpha; /* parameter controlling linearity */
54: PetscReal lambda; /* parameter controlling nonlinearity */
55: } AppCtx;
57: typedef struct {
58: PetscScalar u;
59: PetscScalar v;
60: PetscScalar p;
61: } Field;
63: static PetscScalar Kref[36] = { 0.5, 0.5, -0.5, 0, 0, -0.5,
64: 0.5, 0.5, -0.5, 0, 0, -0.5,
65: -0.5, -0.5, 0.5, 0, 0, 0.5,
66: 0, 0, 0, 0, 0, 0,
67: 0, 0, 0, 0, 0, 0,
68: -0.5, -0.5, 0.5, 0, 0, 0.5};
70: static PetscScalar Identity[9] = {0.0833333, 0.0416667, 0.0416667,
71: 0.0416667, 0.0833333, 0.0416667,
72: 0.0416667, 0.0416667, 0.0833333};
74: static PetscScalar Gradient[18] = {-0.1666667, -0.1666667, -0.1666667,
75: -0.1666667, -0.1666667, -0.1666667,
76: 0.1666667, 0.1666667, 0.1666667,
77: 0.0, 0.0, 0.0,
78: 0.0, 0.0, 0.0,
79: 0.1666667, 0.1666667, 0.1666667};
81: static PetscScalar Divergence[18] = {-0.1666667, 0.1666667, 0.0,
82: -0.1666667, 0.0, 0.1666667,
83: -0.1666667, 0.1666667, 0.0,
84: -0.1666667, 0.0, 0.1666667,
85: -0.1666667, 0.1666667, 0.0,
86: -0.1666667, 0.0, 0.1666667};
88: /* These are */
89: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
90: 0.07503111, 0.64494897,
91: 0.66639025, 0.15505103,
92: 0.28001992, 0.64494897};
93: static PetscScalar quadWeights[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
95: /*
96: User-defined routines
97: */
107: int main(int argc,char **argv)
108: {
109: DMMG *dmmg; /* hierarchy manager */
110: DM da;
111: SNES snes; /* nonlinear solver */
112: AppCtx *user; /* user-defined work context */
113: PetscBag bag;
114: PetscInt its; /* iterations for convergence */
115: SNESConvergedReason reason;
116: PetscBool drawContours; /* flag for drawing contours */
117: PetscErrorCode ierr;
118: PetscReal lambda_max = 6.81, lambda_min = 0.0, error;
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Initialize program
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: PetscInitialize(&argc,&argv,(char *)0,help);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Initialize problem parameters
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
129: PetscBagGetData(bag, (void **) &user);
130: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
131: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
132: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
133: PetscBagSetFromOptions(bag);
134: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
135: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
136: if (user->lambda > lambda_max || user->lambda < lambda_min) {
137: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
138: }
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create multilevel DM data structure (DMMG) to manage hierarchical solvers
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: DMMGCreate(PETSC_COMM_WORLD,1,user,&dmmg);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Create distributed array (DMDA) to manage parallel grid and vectors
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
149: 3,1,PETSC_NULL,PETSC_NULL,&da);
150: DMDASetFieldName(da, 0, "ooblek");
151: DMMGSetDM(dmmg, (DM) da);
152: DMDestroy(&da);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Set the discretization functions
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);
158: DMMGSetFromOptions(dmmg);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Set the problem null space
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: DMMGSetNullSpace(dmmg, PETSC_FALSE, 1, CreateNullSpace);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Evaluate initial guess
167: Note: The user should initialize the vector, x, with the initial guess
168: for the nonlinear solver prior to calling SNESSolve(). In particular,
169: to employ an initial guess of zero, the user should explicitly set
170: this vector to zero by calling VecSet().
171: */
172: DMMGSetInitialGuess(dmmg, FormInitialGuess);
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Solve nonlinear system
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: DMMGSolve(dmmg);
178: snes = DMMGGetSNES(dmmg);
179: SNESGetIterationNumber(snes,&its);
180: SNESGetConvergedReason(snes, &reason);
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
184: L_2Error(DMMGGetDM(dmmg), DMMGGetx(dmmg), &error, user);
185: PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
186: PrintVector(dmmg[0], DMMGGetx(dmmg));
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Free work space. All PETSc objects should be destroyed when they
190: are no longer needed.
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192: DMMGDestroy(dmmg);
193: PetscBagDestroy(&bag);
194: PetscFinalize();
195: return(0);
196: }
200: PetscErrorCode PrintVector(DMMG dmmg, Vec U)
201: {
202: DM da = dmmg->dm;
203: Field **u;
204: PetscInt i,j,xs,ys,xm,ym;
208: DMDAVecGetArray(da,U,&u);
209: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
210: for(j = ys+ym-1; j >= ys; j--) {
211: for(i = xs; i < xs+xm; i++) {
212: printf("u[%d][%d] = (%g, %g, %g) ", j, i, u[j][i].u, u[j][i].v, u[j][i].p);
213: }
214: printf("\n");
215: }
216: DMDAVecRestoreArray(da,U,&u);
217: return(0);
218: }
222: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, Field *u)
223: {
225: (*u).u = x;
226: (*u).v = -y;
227: (*u).p = 0.0;
228: return(0);
229: }
233: PetscErrorCode CreateNullSpace(DMMG dmmg, Vec *nulls)
234: {
235: DM da = dmmg->dm;
236: Vec X = nulls[0];
237: Field **x;
238: PetscInt xs,ys,xm,ym,i,j;
242: DMDAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL);
243: DMDAVecGetArray(da, X, &x);
244: for(j = ys; j < ys+ym; j++) {
245: for(i = xs; i < xs+xm; i++) {
246: x[j][i].u = 0.0;
247: x[j][i].v = 0.0;
248: x[j][i].p = 1.0;
249: }
250: }
251: DMDAVecRestoreArray(da, X, &x);
252: PrintVector(dmmg, X);
253: return(0);
254: }
258: /*
259: FormInitialGuess - Forms initial approximation.
261: Input Parameters:
262: dmmg - The DMMG context
263: X - vector
265: Output Parameter:
266: X - vector
267: */
268: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
269: {
270: AppCtx *user = (AppCtx *) dmmg->user;
271: DM da = dmmg->dm;
272: PetscInt i,j,Mx,My,xs,ys,xm,ym;
274: PetscReal lambda,temp1,temp,hx,hy;
275: Field **x;
278: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
279: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
281: lambda = user->lambda;
282: hx = 1.0/(PetscReal)(Mx-1);
283: hy = 1.0/(PetscReal)(My-1);
284: if (lambda == 0.0) {
285: temp1 = 0.0;
286: } else {
287: temp1 = lambda/(lambda + 1.0);
288: }
290: /*
291: Get a pointer to vector data.
292: - For default PETSc vectors, VecGetArray() returns a pointer to
293: the data array. Otherwise, the routine is implementation dependent.
294: - You MUST call VecRestoreArray() when you no longer need access to
295: the array.
296: */
297: DMDAVecGetArray(da,X,&x);
299: /*
300: Get local grid boundaries (for 2-dimensional DMDA):
301: xs, ys - starting grid indices (no ghost points)
302: xm, ym - widths of local grid (no ghost points)
304: */
305: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
307: /*
308: Compute initial guess over the locally owned part of the grid
309: */
310: for (j=ys; j<ys+ym; j++) {
311: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
312: for (i=xs; i<xs+xm; i++) {
313: #define CHECK_SOLUTION
314: #ifdef CHECK_SOLUTION
315: ExactSolution(i*hx, j*hy, &x[j][i]);
316: #else
317: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
318: /* Boundary conditions are usually zero Dirichlet */
319: ExactSolution(i*hx, j*hy, &x[j][i]);
320: } else {
321: x[j][i].u = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
322: x[j][i].v = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
323: x[j][i].p = 1.0;
324: }
325: #endif
326: }
327: }
329: /*
330: Restore vector
331: */
332: DMDAVecRestoreArray(da,X,&x);
333: PrintVector(dmmg, X);
334: return(0);
335: }
339: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, Field r[])
340: {
341: Field rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
342: PetscScalar phi[3] = {0.0, 0.0, 0.0};
343: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
344: Field res;
345: PetscInt q, k;
348: for(q = 0; q < 4; q++) {
349: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
350: phi[1] = quadPoints[q*2];
351: phi[2] = quadPoints[q*2+1];
352: if (isLower) {
353: x = xI + quadPoints[q*2]*hx;
354: y = yI + quadPoints[q*2+1]*hy;
355: } else {
356: x = xI + 1.0 - quadPoints[q*2]*hx;
357: y = yI + 1.0 - quadPoints[q*2+1]*hy;
358: }
359: res.u = quadWeights[q]*(0.0);
360: res.v = quadWeights[q]*(0.0);
361: res.p = quadWeights[q]*(0.0);
362: for(k = 0; k < 3; k++) {
363: rLocal[k].u += phi[k]*res.u;
364: rLocal[k].v += phi[k]*res.v;
365: rLocal[k].p += phi[k]*res.p;
366: }
367: }
368: for(k = 0; k < 3; k++) {
369: printf(" constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p);
370: r[k].u += lambda*hxhy*rLocal[k].u;
371: r[k].v += lambda*hxhy*rLocal[k].v;
372: r[k].p += hxhy*rLocal[k].p;
373: }
374: return(0);
375: }
377: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[]) {
378: Field rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
379: PetscScalar phi[3] = {0.0, 0.0, 0.0};
380: Field res;
381: PetscInt q;
384: for(q = 0; q < 4; q++) {
385: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
386: phi[1] = quadPoints[q*2];
387: phi[2] = quadPoints[q*2+1];
388: res.u = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
389: res.v = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);
390: rLocal[0].u += phi[0]*res.u;
391: rLocal[0].v += phi[0]*res.v;
392: rLocal[1].u += phi[1]*res.u;
393: rLocal[1].v += phi[1]*res.v;
394: rLocal[2].u += phi[2]*res.u;
395: rLocal[2].v += phi[2]*res.v;
396: }
397: r[0].u += lambda*rLocal[0].u;
398: r[0].v += lambda*rLocal[0].v;
399: r[1].u += lambda*rLocal[1].u;
400: r[1].v += lambda*rLocal[1].v;
401: r[2].u += lambda*rLocal[2].u;
402: r[2].v += lambda*rLocal[2].v;
403: return(0);
404: }
408: /*
409: FormFunctionLocal - Evaluates nonlinear function, F(x).
411: Process adiC(36): FormFunctionLocal
413: */
414: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
415: {
416: Field uLocal[3];
417: Field rLocal[3];
418: PetscScalar G[4];
419: Field uExact;
420: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
421: PetscInt i,j,k,l;
426: /* Naive Jacobian calculation:
428: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
429: \ 0 1/hy / \ 0 hy /
430: */
431: alpha = user->alpha;
432: lambda = user->lambda;
433: hx = 1.0/(PetscReal)(info->mx-1);
434: hy = 1.0/(PetscReal)(info->my-1);
435: sc = hx*hy*lambda;
436: hxhy = hx*hy;
437: detJInv = hxhy;
438: G[0] = (1.0/(hx*hx)) * detJInv;
439: G[1] = 0.0;
440: G[2] = G[1];
441: G[3] = (1.0/(hy*hy)) * detJInv;
442: for(k = 0; k < 4; k++) {
443: printf("G[%d] = %g\n", k, G[k]);
444: }
446: /* Zero the vector */
447: PetscMemzero((void *) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(Field));
448: /* Compute function over the locally owned part of the grid. For each
449: vertex (i,j), we consider the element below:
451: 2 (1) (0)
452: i,j+1 --- i+1,j+1
453: | \ |
454: | \ |
455: | \ |
456: | \ |
457: | \ |
458: i,j --- i+1,j
459: 0 1 (2)
461: and therefore we do not loop over the last vertex in each dimension.
462: */
463: for(j = info->ys; j < info->ys+info->ym-1; j++) {
464: for(i = info->xs; i < info->xs+info->xm-1; i++) {
465: /* Lower element */
466: uLocal[0] = x[j][i];
467: uLocal[1] = x[j][i+1];
468: uLocal[2] = x[j+1][i];
469: printf("Lower Solution ElementVector for (%d, %d)\n", i, j);
470: for(k = 0; k < 3; k++) {
471: printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);
472: }
473: for(k = 0; k < 3; k++) {
474: rLocal[k].u = 0.0;
475: rLocal[k].v = 0.0;
476: rLocal[k].p = 0.0;
477: for(l = 0; l < 3; l++) {
478: rLocal[k].u += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].u;
479: rLocal[k].v += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].v;
480: /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
481: /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
482: /* Gradient */
483: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
484: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
485: /* Divergence */
486: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
487: }
488: }
489: printf("Lower Laplacian ElementVector for (%d, %d)\n", i, j);
490: for(k = 0; k < 3; k++) {
491: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
492: }
493: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
494: printf("Lower Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
495: for(k = 0; k < 3; k++) {
496: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
497: }
498: nonlinearResidual(0.0*sc, uLocal, rLocal);
499: printf("Lower Full nonlinear ElementVector for (%d, %d)\n", i, j);
500: for(k = 0; k < 3; k++) {
501: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
502: }
503: f[j][i].u += rLocal[0].u;
504: f[j][i].v += rLocal[0].v;
505: f[j][i].p += rLocal[0].p;
506: f[j][i+1].u += rLocal[1].u;
507: f[j][i+1].v += rLocal[1].v;
508: f[j][i+1].p += rLocal[1].p;
509: f[j+1][i].u += rLocal[2].u;
510: f[j+1][i].v += rLocal[2].v;
511: f[j+1][i].p += rLocal[2].p;
512: /* Upper element */
513: uLocal[0] = x[j+1][i+1];
514: uLocal[1] = x[j+1][i];
515: uLocal[2] = x[j][i+1];
516: printf("Upper Solution ElementVector for (%d, %d)\n", i, j);
517: for(k = 0; k < 3; k++) {
518: printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);
519: }
520: for(k = 0; k < 3; k++) {
521: rLocal[k].u = 0.0;
522: rLocal[k].v = 0.0;
523: rLocal[k].p = 0.0;
524: for(l = 0; l < 3; l++) {
525: rLocal[k].u += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].u;
526: rLocal[k].v += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].v;
527: /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
528: /* Gradient */
529: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
530: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
531: /* Divergence */
532: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
533: }
534: }
535: printf("Upper Laplacian ElementVector for (%d, %d)\n", i, j);
536: for(k = 0; k < 3; k++) {
537: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
538: }
539: constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
540: printf("Upper Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
541: for(k = 0; k < 3; k++) {
542: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
543: }
544: nonlinearResidual(0.0*sc, uLocal, rLocal);
545: printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);
546: for(k = 0; k < 3; k++) {
547: printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
548: }
549: f[j+1][i+1].u += rLocal[0].u;
550: f[j+1][i+1].v += rLocal[0].v;
551: f[j+1][i+1].p += rLocal[0].p;
552: f[j+1][i].u += rLocal[1].u;
553: f[j+1][i].v += rLocal[1].v;
554: f[j+1][i].p += rLocal[1].p;
555: f[j][i+1].u += rLocal[2].u;
556: f[j][i+1].v += rLocal[2].v;
557: f[j][i+1].p += rLocal[2].p;
558: /* Boundary conditions */
559: if (i == 0 || j == 0) {
560: ExactSolution(i*hx, j*hy, &uExact);
561: f[j][i].u = x[j][i].u - uExact.u;
562: f[j][i].v = x[j][i].v - uExact.v;
563: f[j][i].p = x[j][i].p - uExact.p;
564: }
565: if ((i == info->mx-2) || (j == 0)) {
566: ExactSolution((i+1)*hx, j*hy, &uExact);
567: f[j][i+1].u = x[j][i+1].u - uExact.u;
568: f[j][i+1].v = x[j][i+1].v - uExact.v;
569: f[j][i+1].p = x[j][i+1].p - uExact.p;
570: }
571: if ((i == info->mx-2) || (j == info->my-2)) {
572: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
573: f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
574: f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
575: f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
576: }
577: if ((i == 0) || (j == info->my-2)) {
578: ExactSolution(i*hx, (j+1)*hy, &uExact);
579: f[j+1][i].u = x[j+1][i].u - uExact.u;
580: f[j+1][i].v = x[j+1][i].v - uExact.v;
581: f[j+1][i].p = x[j+1][i].p - uExact.p;
582: }
583: }
584: }
586: for(j = info->ys+info->ym-1; j >= info->ys; j--) {
587: for(i = info->xs; i < info->xs+info->xm; i++) {
588: printf("f[%d][%d] = (%g, %g, %g) ", j, i, f[j][i].u, f[j][i].v, f[j][i].p);
589: }
590: printf("\n");
591: }
592: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
593: return(0);
594: }
596: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[]) {
598: return(0);
599: }
603: /*
604: FormJacobianLocal - Evaluates Jacobian matrix.
605: */
606: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat jac, AppCtx *user)
607: {
608: Field uLocal[4];
609: PetscScalar JLocal[144];
610: MatStencil rows[4*3], cols[4*3], ident;
611: PetscInt lowerRow[3] = {0, 1, 3};
612: PetscInt upperRow[3] = {2, 3, 1};
613: PetscInt hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
614: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
615: PetscInt i,j,k,l,numRows,dof = info->dof;
619: alpha = user->alpha;
620: lambda = user->lambda;
621: hx = 1.0/(PetscReal)(info->mx-1);
622: hy = 1.0/(PetscReal)(info->my-1);
623: sc = hx*hy*lambda;
624: hxhy = hx*hy;
625: detJInv = hxhy;
626: G[0] = (1.0/(hx*hx)) * detJInv;
627: G[1] = 0.0;
628: G[2] = G[1];
629: G[3] = (1.0/(hy*hy)) * detJInv;
630: for(k = 0; k < 4; k++) {
631: printf("G[%d] = %g\n", k, G[k]);
632: }
634: MatZeroEntries(jac);
635: /*
636: Compute entries for the locally owned part of the Jacobian.
637: - Currently, all PETSc parallel matrix formats are partitioned by
638: contiguous chunks of rows across the processors.
639: - Each processor needs to insert only elements that it owns
640: locally (but any non-local elements will be sent to the
641: appropriate processor during matrix assembly).
642: - Here, we set all entries for a particular row at once.
643: - We can set matrix entries either using either
644: MatSetValuesLocal() or MatSetValues(), as discussed above.
645: */
646: #define NOT_PRES_BC 1
647: for (j=info->ys; j<info->ys+info->ym-1; j++) {
648: for (i=info->xs; i<info->xs+info->xm-1; i++) {
649: PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
650: numRows = 0;
651: /* Lower element */
652: uLocal[0] = x[j][i];
653: uLocal[1] = x[j][i+1];
654: uLocal[2] = x[j+1][i+1];
655: uLocal[3] = x[j+1][i];
656: /* i,j */
657: if (i == 0 || j == 0) {
658: hasLower[0] = 0;
659: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
660: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
661: ident.i = i; ident.j = j; ident.c = 0;
662: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
663: ident.i = i; ident.j = j; ident.c = 1;
664: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
665: #ifdef PRES_BC
666: ident.i = i; ident.j = j; ident.c = 2;
667: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
668: #endif
669: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
670: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
671: } else {
672: hasLower[0] = 1;
673: velocityRows[0] = numRows;
674: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
675: numRows++;
676: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
677: numRows++;
678: #ifdef PRES_BC
679: pressureRows[0] = numRows;
680: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
681: numRows++;
682: #endif
683: }
684: #ifndef PRES_BC
685: pressureRows[0] = numRows;
686: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
687: numRows++;
688: #endif
689: cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
690: cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
691: cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
692: /* i+1,j */
693: if ((i == info->mx-2) || (j == 0)) {
694: hasLower[1] = 0;
695: hasUpper[2] = 0;
696: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
697: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
698: ident.i = i+1; ident.j = j; ident.c = 0;
699: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
700: ident.i = i+1; ident.j = j; ident.c = 1;
701: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
702: #ifdef PRES_BC
703: ident.i = i+1; ident.j = j; ident.c = 2;
704: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
705: #endif
706: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
707: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
708: } else {
709: hasLower[1] = 1;
710: hasUpper[2] = 1;
711: velocityRows[1] = numRows;
712: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
713: numRows++;
714: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
715: numRows++;
716: #ifdef PRES_BC
717: pressureRows[1] = numRows;
718: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
719: numRows++;
720: #endif
721: }
722: #ifndef PRES_BC
723: pressureRows[1] = numRows;
724: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
725: numRows++;
726: #endif
727: cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
728: cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
729: cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
730: /* i+1,j+1 */
731: if ((i == info->mx-2) || (j == info->my-2)) {
732: hasUpper[0] = 0;
733: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
734: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
735: ident.i = i+1; ident.j = j+1; ident.c = 0;
736: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
737: ident.i = i+1; ident.j = j+1; ident.c = 1;
738: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
739: #ifdef PRES_BC
740: ident.i = i+1; ident.j = j+1; ident.c = 2;
741: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
742: #endif
743: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
744: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
745: } else {
746: hasUpper[0] = 1;
747: velocityRows[2] = numRows;
748: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
749: numRows++;
750: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
751: numRows++;
752: #ifdef PRES_BC
753: pressureRows[2] = numRows;
754: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
755: numRows++;
756: #endif
757: }
758: #ifndef PRES_BC
759: pressureRows[2] = numRows;
760: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
761: numRows++;
762: #endif
763: cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
764: cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
765: cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
766: /* i,j+1 */
767: if ((i == 0) || (j == info->my-2)) {
768: hasLower[2] = 0;
769: hasUpper[1] = 0;
770: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
771: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
772: ident.i = i; ident.j = j+1; ident.c = 0;
773: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
774: ident.i = i; ident.j = j+1; ident.c = 1;
775: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
776: #ifdef PRES_BC
777: ident.i = i; ident.j = j+1; ident.c = 2;
778: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
779: #endif
780: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
781: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
782: } else {
783: hasLower[2] = 1;
784: hasUpper[1] = 1;
785: velocityRows[3] = numRows;
786: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
787: numRows++;
788: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
789: numRows++;
790: #ifdef PRES_BC
791: pressureRows[3] = numRows;
792: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
793: numRows++;
794: #endif
795: }
796: #ifndef PRES_BC
797: pressureRows[3] = numRows;
798: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
799: numRows++;
800: #endif
801: cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
802: cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
803: cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;
805: /* Lower Element */
806: for(k = 0; k < 3; k++) {
807: #ifdef PRES_BC
808: if (!hasLower[k]) continue;
809: #endif
810: for(l = 0; l < 3; l++) {
811: /* Divergence */
812: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
813: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
814: /* JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += Identity[k*3 + l]; */
815: }
816: if (!hasLower[k]) continue;
817: for(l = 0; l < 3; l++) {
818: /* Laplacian */
819: JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
820: JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
821: /* JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += Identity[k*3 + l]; */
822: /* JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += Identity[k*3 + l]; */
823: /* Gradient */
824: JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
825: JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
826: }
827: }
828: /* Upper Element */
829: for(k = 0; k < 3; k++) {
830: #ifdef PRES_BC
831: if (!hasUpper[k]) continue;
832: #endif
833: for(l = 0; l < 3; l++) {
834: /* Divergence */
835: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
836: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
837: /* JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += Identity[k*3 + l]; */
838: }
839: if (!hasUpper[k]) continue;
840: for(l = 0; l < 3; l++) {
841: /* Laplacian */
842: JLocal[velocityRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
843: JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+1] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
844: /* Gradient */
845: JLocal[velocityRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
846: JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
847: }
848: }
850: nonlinearJacobian(-1.0*sc, uLocal, JLocal);
851: printf("Element matrix for (%d, %d)\n", i, j);
852: printf(" col ");
853: for(l = 0; l < 4*3; l++) {
854: printf("(%d, %d, %d) ", cols[l].i, cols[l].j, cols[l].c);
855: }
856: printf("\n");
857: for(k = 0; k < numRows; k++) {
858: printf("row (%d, %d, %d): ", rows[k].i, rows[k].j, rows[k].c);
859: for(l = 0; l < 4; l++) {
860: printf("%9.6g %9.6g %9.6g ", JLocal[k*dof*4 + l*dof+0], JLocal[k*dof*4 + l*dof+1], JLocal[k*dof*4 + l*dof+2]);
861: }
862: printf("\n");
863: }
864: MatSetValuesStencil(jac,numRows,rows,4*dof,cols, JLocal,ADD_VALUES);
865: }
866: }
868: /*
869: Assemble matrix, using the 2-step process:
870: MatAssemblyBegin(), MatAssemblyEnd().
871: */
872: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
873: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
874: /*
875: Tell the matrix we will never add a new nonzero location to the
876: matrix. If we do, it will generate an error.
877: */
878: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
879: return(0);
880: }
884: /*
885: L_2Error - Integrate the L_2 error of our solution over each face
886: */
887: PetscErrorCode L_2Error(DM da, Vec fVec, double *error, AppCtx *user)
888: {
889: DMDALocalInfo info;
890: Vec fLocalVec;
891: Field **f;
892: Field u, uExact, uLocal[4];
893: PetscScalar hx, hy, hxhy, x, y, phi[3];
894: PetscInt i, j, q;
898: DMDAGetLocalInfo(da, &info);
899: DMGetLocalVector(da, &fLocalVec);
900: DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
901: DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
902: DMDAVecGetArray(da, fLocalVec, &f);
904: *error = 0.0;
905: hx = 1.0/(PetscReal)(info.mx-1);
906: hy = 1.0/(PetscReal)(info.my-1);
907: hxhy = hx*hy;
908: for(j = info.ys; j < info.ys+info.ym-1; j++) {
909: for(i = info.xs; i < info.xs+info.xm-1; i++) {
910: uLocal[0] = f[j][i];
911: uLocal[1] = f[j][i+1];
912: uLocal[2] = f[j+1][i+1];
913: uLocal[3] = f[j+1][i];
914: /* Lower element */
915: for(q = 0; q < 4; q++) {
916: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
917: phi[1] = quadPoints[q*2];
918: phi[2] = quadPoints[q*2+1];
919: u.u = uLocal[0].u*phi[0]+ uLocal[1].u*phi[1] + uLocal[3].u*phi[2];
920: u.v = uLocal[0].v*phi[0]+ uLocal[1].v*phi[1] + uLocal[3].v*phi[2];
921: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
922: x = (quadPoints[q*2] + i)*hx;
923: y = (quadPoints[q*2+1] + j)*hy;
924: ExactSolution(x, y, &uExact);
925: *error += hxhy*quadWeights[q]*((u.u - uExact.u)*(u.u - uExact.u) + (u.v - uExact.v)*(u.v - uExact.v) + (u.p - uExact.p)*(u.p - uExact.p));
926: }
927: /* Upper element */
928: /*
929: The affine map from the lower to the upper is
931: / x_U \ = / -1 0 \ / x_L \ + / hx \
932: \ y_U / \ 0 -1 / \ y_L / \ hy /
933: */
934: for(q = 0; q < 4; q++) {
935: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
936: phi[1] = quadPoints[q*2];
937: phi[2] = quadPoints[q*2+1];
938: u.u = uLocal[2].u*phi[0]+ uLocal[3].u*phi[1] + uLocal[1].u*phi[2];
939: u.v = uLocal[2].v*phi[0]+ uLocal[3].v*phi[1] + uLocal[1].v*phi[2];
940: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
941: x = (1.0 - quadPoints[q*2] + i)*hx;
942: y = (1.0 - quadPoints[q*2+1] + j)*hy;
943: ExactSolution(x, y, &uExact);
944: *error += hxhy*quadWeights[q]*((u.u - uExact.u)*(u.u - uExact.u) + (u.v - uExact.v)*(u.v - uExact.v) + (u.p - uExact.p)*(u.p - uExact.p));
945: }
946: }
947: }
949: DMDAVecRestoreArray(da, fLocalVec, &f);
950: /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
951: /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
952: DMRestoreLocalVector(da, &fLocalVec);
953: return(0);
954: }