Actual source code: ex7.c
petsc-3.3-p7 2013-05-11
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 <petscsnes.h>
46: /*
47: User-defined application context - contains data needed by the
48: application-provided call-back routines, FormJacobianLocal() and
49: FormFunctionLocal().
50: */
51: typedef struct {
52: PetscReal alpha; /* parameter controlling linearity */
53: PetscReal lambda; /* parameter controlling nonlinearity */
54: } AppCtx;
56: typedef struct {
57: PetscScalar u;
58: PetscScalar v;
59: PetscScalar p;
60: } Field;
62: static PetscScalar Kref[36] = { 0.5, 0.5, -0.5, 0, 0, -0.5,
63: 0.5, 0.5, -0.5, 0, 0, -0.5,
64: -0.5, -0.5, 0.5, 0, 0, 0.5,
65: 0, 0, 0, 0, 0, 0,
66: 0, 0, 0, 0, 0, 0,
67: -0.5, -0.5, 0.5, 0, 0, 0.5};
69: static PetscScalar Gradient[18] = {-0.1666667, -0.1666667, -0.1666667,
70: -0.1666667, -0.1666667, -0.1666667,
71: 0.1666667, 0.1666667, 0.1666667,
72: 0.0, 0.0, 0.0,
73: 0.0, 0.0, 0.0,
74: 0.1666667, 0.1666667, 0.1666667};
76: static PetscScalar Divergence[18] = {-0.1666667, 0.1666667, 0.0,
77: -0.1666667, 0.0, 0.1666667,
78: -0.1666667, 0.1666667, 0.0,
79: -0.1666667, 0.0, 0.1666667,
80: -0.1666667, 0.1666667, 0.0,
81: -0.1666667, 0.0, 0.1666667};
83: /* These are */
84: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
85: 0.07503111, 0.64494897,
86: 0.66639025, 0.15505103,
87: 0.28001992, 0.64494897};
88: static PetscScalar quadWeights[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
90: /*
91: User-defined routines
92: */
93: extern PetscErrorCode CreateNullSpace(DM, Vec*);
94: extern PetscErrorCode FormInitialGuess(DM,Vec);
95: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,AppCtx*);
96: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,Field**,Mat,AppCtx*);
97: extern PetscErrorCode L_2Error(DM, Vec, double *, AppCtx *);
101: int main(int argc,char **argv)
102: {
103: DM da;
104: SNES snes; /* nonlinear solver */
105: AppCtx *user; /* user-defined work context */
106: PetscBag bag;
107: PetscInt its; /* iterations for convergence */
108: SNESConvergedReason reason;
109: PetscErrorCode ierr;
110: PetscReal lambda_max = 6.81, lambda_min = 0.0, error;
111: Vec x;
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Initialize program
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: PetscInitialize(&argc,&argv,(char *)0,help);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Initialize problem parameters
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
122: PetscBagGetData(bag, (void **) &user);
123: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
124: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
125: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
126: PetscBagSetFromOptions(bag);
127: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
128: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
129: if (user->lambda > lambda_max || user->lambda < lambda_min) {
130: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
131: }
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Create multilevel DM data structure (SNES) to manage hierarchical solvers
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: SNESCreate(PETSC_COMM_WORLD,&snes);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Create distributed array (DMDA) to manage parallel grid and vectors
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
142: 3,1,PETSC_NULL,PETSC_NULL,&da);
143: DMDASetFieldName(da, 0, "ooblek");
144: DMSetApplicationContext(da,&user);
145: SNESSetDM(snes, (DM) da);
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Set the discretization functions
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
151: DMDASetLocalJacobian(da, (DMDALocalFunction1)FormJacobianLocal);
152: SNESSetFromOptions(snes);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Evaluate initial guess
156: Note: The user should initialize the vector, x, with the initial guess
157: for the nonlinear solver prior to calling SNESSolve(). In particular,
158: to employ an initial guess of zero, the user should explicitly set
159: this vector to zero by calling VecSet().
160: */
161: DMSetInitialGuess(da, FormInitialGuess);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Solve nonlinear system
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: SNESSolve(snes,PETSC_NULL,PETSC_NULL);
167: SNESGetIterationNumber(snes,&its);
168: SNESGetConvergedReason(snes, &reason);
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
172: DMDestroy(&da);
173: SNESGetDM(snes,&da);
174: SNESGetSolution(snes,&x);
175: L_2Error(da, x, &error, user);
176: PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Free work space. All PETSc objects should be destroyed when they
180: are no longer needed.
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: SNESDestroy(&snes);
183: PetscBagDestroy(&bag);
184: PetscFinalize();
185: return(0);
186: }
190: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, Field *u)
191: {
193: (*u).u = x;
194: (*u).v = -y;
195: (*u).p = 0.0;
196: return(0);
197: }
201: PetscErrorCode CreateNullSpace(DM da, Vec *N)
202: {
203: Field **x;
204: PetscInt xs,ys,xm,ym,i,j;
208: DMDAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL);
209: DMCreateGlobalVector(da,N);
210: DMDAVecGetArray(da, *N, &x);
211: for(j = ys; j < ys+ym; j++) {
212: for(i = xs; i < xs+xm; i++) {
213: x[j][i].u = 0.0;
214: x[j][i].v = 0.0;
215: x[j][i].p = 1.0;
216: }
217: }
218: DMDAVecRestoreArray(da, *N, &x);
219: return(0);
220: }
224: /*
225: FormInitialGuess - Forms initial approximation.
227: Input Parameters:
228: dm - The DM context
229: X - vector
231: Output Parameter:
232: X - vector
233: */
234: PetscErrorCode FormInitialGuess(DM da,Vec X)
235: {
236: AppCtx *user;
237: PetscInt i,j,Mx,My,xs,ys,xm,ym;
239: PetscReal lambda,temp1,temp,hx,hy;
240: Field **x;
243: DMGetApplicationContext(da,&user);
244: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
245: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
247: lambda = user->lambda;
248: hx = 1.0/(PetscReal)(Mx-1);
249: hy = 1.0/(PetscReal)(My-1);
250: if (lambda == 0.0) {
251: temp1 = 0.0;
252: } else {
253: temp1 = lambda/(lambda + 1.0);
254: }
256: /*
257: Get a pointer to vector data.
258: - For default PETSc vectors, VecGetArray() returns a pointer to
259: the data array. Otherwise, the routine is implementation dependent.
260: - You MUST call VecRestoreArray() when you no longer need access to
261: the array.
262: */
263: DMDAVecGetArray(da,X,&x);
265: /*
266: Get local grid boundaries (for 2-dimensional DMDA):
267: xs, ys - starting grid indices (no ghost points)
268: xm, ym - widths of local grid (no ghost points)
270: */
271: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
273: /*
274: Compute initial guess over the locally owned part of the grid
275: */
276: for (j=ys; j<ys+ym; j++) {
277: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
278: for (i=xs; i<xs+xm; i++) {
279: #define CHECK_SOLUTION
280: #ifdef CHECK_SOLUTION
281: ExactSolution(i*hx, j*hy, &x[j][i]);
282: #else
283: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
284: /* Boundary conditions are usually zero Dirichlet */
285: ExactSolution(i*hx, j*hy, &x[j][i]);
286: } else {
287: x[j][i].u = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
288: x[j][i].v = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
289: x[j][i].p = 1.0;
290: }
291: #endif
292: }
293: }
295: /*
296: Restore vector
297: */
298: DMDAVecRestoreArray(da,X,&x);
299: return(0);
300: }
304: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, Field r[])
305: {
306: Field rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
307: PetscScalar phi[3] = {0.0, 0.0, 0.0};
308: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
309: Field res;
310: PetscInt q, k;
313: for(q = 0; q < 4; q++) {
314: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
315: phi[1] = quadPoints[q*2];
316: phi[2] = quadPoints[q*2+1];
317: if (isLower) {
318: x = xI + quadPoints[q*2]*hx;
319: y = yI + quadPoints[q*2+1]*hy;
320: } else {
321: x = xI + 1.0 - quadPoints[q*2]*hx;
322: y = yI + 1.0 - quadPoints[q*2+1]*hy;
323: }
324: res.u = quadWeights[q]*(0.0);
325: res.v = quadWeights[q]*(0.0);
326: res.p = quadWeights[q]*(0.0);
327: for(k = 0; k < 3; k++) {
328: rLocal[k].u += phi[k]*res.u;
329: rLocal[k].v += phi[k]*res.v;
330: rLocal[k].p += phi[k]*res.p;
331: }
332: }
333: for(k = 0; k < 3; k++) {
334: /*printf(" constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p); */
335: r[k].u += lambda*hxhy*rLocal[k].u;
336: r[k].v += lambda*hxhy*rLocal[k].v;
337: r[k].p += hxhy*rLocal[k].p;
338: }
339: return(0);
340: }
344: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[]) {
345: Field rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
346: PetscScalar phi[3] = {0.0, 0.0, 0.0};
347: Field res;
348: PetscInt q;
351: for(q = 0; q < 4; q++) {
352: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
353: phi[1] = quadPoints[q*2];
354: phi[2] = quadPoints[q*2+1];
355: res.u = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
356: res.v = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);
357: rLocal[0].u += phi[0]*res.u;
358: rLocal[0].v += phi[0]*res.v;
359: rLocal[1].u += phi[1]*res.u;
360: rLocal[1].v += phi[1]*res.v;
361: rLocal[2].u += phi[2]*res.u;
362: rLocal[2].v += phi[2]*res.v;
363: }
364: r[0].u += lambda*rLocal[0].u;
365: r[0].v += lambda*rLocal[0].v;
366: r[1].u += lambda*rLocal[1].u;
367: r[1].v += lambda*rLocal[1].v;
368: r[2].u += lambda*rLocal[2].u;
369: r[2].v += lambda*rLocal[2].v;
370: return(0);
371: }
375: /*
376: FormFunctionLocal - Evaluates nonlinear function, F(x).
378: Process adiC(36): FormFunctionLocal
380: */
381: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
382: {
383: Field uLocal[3];
384: Field rLocal[3];
385: PetscScalar G[4];
386: Field uExact;
387: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
388: PetscInt i,j,k,l;
393: /* Naive Jacobian calculation:
395: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
396: \ 0 1/hy / \ 0 hy /
397: */
398: alpha = user->alpha;
399: lambda = user->lambda;
400: hx = 1.0/(PetscReal)(info->mx-1);
401: hy = 1.0/(PetscReal)(info->my-1);
402: sc = hx*hy*lambda;
403: hxhy = hx*hy;
404: detJInv = hxhy;
405: G[0] = (1.0/(hx*hx)) * detJInv;
406: G[1] = 0.0;
407: G[2] = G[1];
408: G[3] = (1.0/(hy*hy)) * detJInv;
409: for(k = 0; k < 4; k++) {
410: /* printf("G[%d] = %g\n", k, G[k]);*/
411: }
413: /* Zero the vector */
414: PetscMemzero((void *) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(Field));
415: /* Compute function over the locally owned part of the grid. For each
416: vertex (i,j), we consider the element below:
418: 2 (1) (0)
419: i,j+1 --- i+1,j+1
420: | \ |
421: | \ |
422: | \ |
423: | \ |
424: | \ |
425: i,j --- i+1,j
426: 0 1 (2)
428: and therefore we do not loop over the last vertex in each dimension.
429: */
430: for(j = info->ys; j < info->ys+info->ym-1; j++) {
431: for(i = info->xs; i < info->xs+info->xm-1; i++) {
432: /* Lower element */
433: uLocal[0] = x[j][i];
434: uLocal[1] = x[j][i+1];
435: uLocal[2] = x[j+1][i];
436: /* printf("Lower Solution ElementVector for (%d, %d)\n", i, j);*/
437: for(k = 0; k < 3; k++) {
438: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
439: }
440: for(k = 0; k < 3; k++) {
441: rLocal[k].u = 0.0;
442: rLocal[k].v = 0.0;
443: rLocal[k].p = 0.0;
444: for(l = 0; l < 3; l++) {
445: 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;
446: 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;
447: /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
448: /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
449: /* Gradient */
450: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
451: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
452: /* Divergence */
453: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
454: }
455: }
456: /* printf("Lower Laplacian ElementVector for (%d, %d)\n", i, j);*/
457: for(k = 0; k < 3; k++) {
458: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
459: }
460: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
461: /* printf("Lower Laplacian+Constant ElementVector for (%d, %d)\n", i, j);*/
462: for(k = 0; k < 3; k++) {
463: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
464: }
465: nonlinearResidual(0.0*sc, uLocal, rLocal);
466: /* printf("Lower Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
467: for(k = 0; k < 3; k++) {
468: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
469: }
470: f[j][i].u += rLocal[0].u;
471: f[j][i].v += rLocal[0].v;
472: f[j][i].p += rLocal[0].p;
473: f[j][i+1].u += rLocal[1].u;
474: f[j][i+1].v += rLocal[1].v;
475: f[j][i+1].p += rLocal[1].p;
476: f[j+1][i].u += rLocal[2].u;
477: f[j+1][i].v += rLocal[2].v;
478: f[j+1][i].p += rLocal[2].p;
479: /* Upper element */
480: uLocal[0] = x[j+1][i+1];
481: uLocal[1] = x[j+1][i];
482: uLocal[2] = x[j][i+1];
483: /* printf("Upper Solution ElementVector for (%d, %d)\n", i, j);*/
484: for(k = 0; k < 3; k++) {
485: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
486: }
487: for(k = 0; k < 3; k++) {
488: rLocal[k].u = 0.0;
489: rLocal[k].v = 0.0;
490: rLocal[k].p = 0.0;
491: for(l = 0; l < 3; l++) {
492: 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;
493: 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;
494: /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
495: /* Gradient */
496: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
497: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
498: /* Divergence */
499: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
500: }
501: }
502: /* printf("Upper Laplacian ElementVector for (%d, %d)\n", i, j);*/
503: for(k = 0; k < 3; k++) {
504: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
505: }
506: constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
507: /* printf("Upper Laplacian+Constant ElementVector for (%d, %d)\n", i, j);*/
508: for(k = 0; k < 3; k++) {
509: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
510: }
511: nonlinearResidual(0.0*sc, uLocal, rLocal);
512: /* printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
513: for(k = 0; k < 3; k++) {
514: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
515: }
516: f[j+1][i+1].u += rLocal[0].u;
517: f[j+1][i+1].v += rLocal[0].v;
518: f[j+1][i+1].p += rLocal[0].p;
519: f[j+1][i].u += rLocal[1].u;
520: f[j+1][i].v += rLocal[1].v;
521: f[j+1][i].p += rLocal[1].p;
522: f[j][i+1].u += rLocal[2].u;
523: f[j][i+1].v += rLocal[2].v;
524: f[j][i+1].p += rLocal[2].p;
525: /* Boundary conditions */
526: if (i == 0 || j == 0) {
527: ExactSolution(i*hx, j*hy, &uExact);
528: f[j][i].u = x[j][i].u - uExact.u;
529: f[j][i].v = x[j][i].v - uExact.v;
530: f[j][i].p = x[j][i].p - uExact.p;
531: }
532: if ((i == info->mx-2) || (j == 0)) {
533: ExactSolution((i+1)*hx, j*hy, &uExact);
534: f[j][i+1].u = x[j][i+1].u - uExact.u;
535: f[j][i+1].v = x[j][i+1].v - uExact.v;
536: f[j][i+1].p = x[j][i+1].p - uExact.p;
537: }
538: if ((i == info->mx-2) || (j == info->my-2)) {
539: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
540: f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
541: f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
542: f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
543: }
544: if ((i == 0) || (j == info->my-2)) {
545: ExactSolution(i*hx, (j+1)*hy, &uExact);
546: f[j+1][i].u = x[j+1][i].u - uExact.u;
547: f[j+1][i].v = x[j+1][i].v - uExact.v;
548: f[j+1][i].p = x[j+1][i].p - uExact.p;
549: }
550: }
551: }
553: for(j = info->ys+info->ym-1; j >= info->ys; j--) {
554: for(i = info->xs; i < info->xs+info->xm; i++) {
555: /* printf("f[%d][%d] = (%g, %g, %g) ", j, i, f[j][i].u, f[j][i].v, f[j][i].p);*/
556: }
557: /*printf("\n");*/
558: }
559: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
560: return(0);
561: }
565: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[]) {
567: return(0);
568: }
572: /*
573: FormJacobianLocal - Evaluates Jacobian matrix.
574: */
575: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat jac, AppCtx *user)
576: {
577: Field uLocal[4];
578: PetscScalar JLocal[144];
579: MatStencil rows[4*3], cols[4*3], ident;
580: PetscInt lowerRow[3] = {0, 1, 3};
581: PetscInt upperRow[3] = {2, 3, 1};
582: PetscInt hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
583: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
584: PetscInt i,j,k,l,numRows,dof = info->dof;
586: MatNullSpace nullspace;
587: Vec N;
590: alpha = user->alpha;
591: lambda = user->lambda;
592: hx = 1.0/(PetscReal)(info->mx-1);
593: hy = 1.0/(PetscReal)(info->my-1);
594: sc = hx*hy*lambda;
595: hxhy = hx*hy;
596: detJInv = hxhy;
597: G[0] = (1.0/(hx*hx)) * detJInv;
598: G[1] = 0.0;
599: G[2] = G[1];
600: G[3] = (1.0/(hy*hy)) * detJInv;
601: for(k = 0; k < 4; k++) {
602: /* printf("G[%d] = %g\n", k, G[k]);*/
603: }
605: MatZeroEntries(jac);
606: /*
607: Compute entries for the locally owned part of the Jacobian.
608: - Currently, all PETSc parallel matrix formats are partitioned by
609: contiguous chunks of rows across the processors.
610: - Each processor needs to insert only elements that it owns
611: locally (but any non-local elements will be sent to the
612: appropriate processor during matrix assembly).
613: - Here, we set all entries for a particular row at once.
614: - We can set matrix entries either using either
615: MatSetValuesLocal() or MatSetValues(), as discussed above.
616: */
617: #define NOT_PRES_BC 1
618: for (j=info->ys; j<info->ys+info->ym-1; j++) {
619: for (i=info->xs; i<info->xs+info->xm-1; i++) {
620: PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
621: numRows = 0;
622: /* Lower element */
623: uLocal[0] = x[j][i];
624: uLocal[1] = x[j][i+1];
625: uLocal[2] = x[j+1][i+1];
626: uLocal[3] = x[j+1][i];
627: /* i,j */
628: if (i == 0 || j == 0) {
629: hasLower[0] = 0;
630: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
631: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
632: ident.i = i; ident.j = j; ident.c = 0;
633: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
634: ident.i = i; ident.j = j; ident.c = 1;
635: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
636: #ifdef PRES_BC
637: ident.i = i; ident.j = j; ident.c = 2;
638: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
639: #endif
640: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
641: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
642: } else {
643: hasLower[0] = 1;
644: velocityRows[0] = numRows;
645: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
646: numRows++;
647: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
648: numRows++;
649: #ifdef PRES_BC
650: pressureRows[0] = numRows;
651: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
652: numRows++;
653: #endif
654: }
655: #ifndef PRES_BC
656: pressureRows[0] = numRows;
657: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
658: numRows++;
659: #endif
660: cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
661: cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
662: cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
663: /* i+1,j */
664: if ((i == info->mx-2) || (j == 0)) {
665: hasLower[1] = 0;
666: hasUpper[2] = 0;
667: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
668: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
669: ident.i = i+1; ident.j = j; ident.c = 0;
670: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
671: ident.i = i+1; ident.j = j; ident.c = 1;
672: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
673: #ifdef PRES_BC
674: ident.i = i+1; ident.j = j; ident.c = 2;
675: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
676: #endif
677: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
678: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
679: } else {
680: hasLower[1] = 1;
681: hasUpper[2] = 1;
682: velocityRows[1] = numRows;
683: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
684: numRows++;
685: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
686: numRows++;
687: #ifdef PRES_BC
688: pressureRows[1] = numRows;
689: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
690: numRows++;
691: #endif
692: }
693: #ifndef PRES_BC
694: pressureRows[1] = numRows;
695: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
696: numRows++;
697: #endif
698: cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
699: cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
700: cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
701: /* i+1,j+1 */
702: if ((i == info->mx-2) || (j == info->my-2)) {
703: hasUpper[0] = 0;
704: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
705: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
706: ident.i = i+1; ident.j = j+1; ident.c = 0;
707: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
708: ident.i = i+1; ident.j = j+1; ident.c = 1;
709: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
710: #ifdef PRES_BC
711: ident.i = i+1; ident.j = j+1; ident.c = 2;
712: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
713: #endif
714: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
715: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
716: } else {
717: hasUpper[0] = 1;
718: velocityRows[2] = numRows;
719: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
720: numRows++;
721: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
722: numRows++;
723: #ifdef PRES_BC
724: pressureRows[2] = numRows;
725: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
726: numRows++;
727: #endif
728: }
729: #ifndef PRES_BC
730: pressureRows[2] = numRows;
731: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
732: numRows++;
733: #endif
734: cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
735: cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
736: cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
737: /* i,j+1 */
738: if ((i == 0) || (j == info->my-2)) {
739: hasLower[2] = 0;
740: hasUpper[1] = 0;
741: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
742: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
743: ident.i = i; ident.j = j+1; ident.c = 0;
744: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
745: ident.i = i; ident.j = j+1; ident.c = 1;
746: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
747: #ifdef PRES_BC
748: ident.i = i; ident.j = j+1; ident.c = 2;
749: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
750: #endif
751: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
752: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
753: } else {
754: hasLower[2] = 1;
755: hasUpper[1] = 1;
756: velocityRows[3] = numRows;
757: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
758: numRows++;
759: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
760: numRows++;
761: #ifdef PRES_BC
762: pressureRows[3] = numRows;
763: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
764: numRows++;
765: #endif
766: }
767: #ifndef PRES_BC
768: pressureRows[3] = numRows;
769: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
770: numRows++;
771: #endif
772: cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
773: cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
774: cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;
776: /* Lower Element */
777: for(k = 0; k < 3; k++) {
778: #ifdef PRES_BC
779: if (!hasLower[k]) continue;
780: #endif
781: for(l = 0; l < 3; l++) {
782: /* Divergence */
783: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
784: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
785: /* JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += Identity[k*3 + l]; */
786: }
787: if (!hasLower[k]) continue;
788: for(l = 0; l < 3; l++) {
789: /* Laplacian */
790: 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]);
791: 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]);
792: /* JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += Identity[k*3 + l]; */
793: /* JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += Identity[k*3 + l]; */
794: /* Gradient */
795: JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
796: JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
797: }
798: }
799: /* Upper Element */
800: for(k = 0; k < 3; k++) {
801: #ifdef PRES_BC
802: if (!hasUpper[k]) continue;
803: #endif
804: for(l = 0; l < 3; l++) {
805: /* Divergence */
806: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
807: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
808: /* JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += Identity[k*3 + l]; */
809: }
810: if (!hasUpper[k]) continue;
811: for(l = 0; l < 3; l++) {
812: /* Laplacian */
813: 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]);
814: 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]);
815: /* Gradient */
816: JLocal[velocityRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
817: JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
818: }
819: }
821: nonlinearJacobian(-1.0*sc, uLocal, JLocal);
822: /* printf("Element matrix for (%d, %d)\n", i, j);*/
823: /* printf(" col ");*/
824: for(l = 0; l < 4*3; l++) {
825: /* printf("(%d, %d, %d) ", cols[l].i, cols[l].j, cols[l].c);*/
826: }
827: /* printf("\n");*/
828: for(k = 0; k < numRows; k++) {
829: /* printf("row (%d, %d, %d): ", rows[k].i, rows[k].j, rows[k].c);*/
830: for(l = 0; l < 4; l++) {
831: /* 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]);*/
832: }
833: /* printf("\n");*/
834: }
835: MatSetValuesStencil(jac,numRows,rows,4*dof,cols, JLocal,ADD_VALUES);
836: }
837: }
839: /*
840: Assemble matrix, using the 2-step process:
841: MatAssemblyBegin(), MatAssemblyEnd().
842: */
843: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
844: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
845: /*
846: Tell the matrix we will never add a new nonzero location to the
847: matrix. If we do, it will generate an error.
848: */
849: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
851: CreateNullSpace(info->da,&N);
852: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&N,&nullspace);
853: VecDestroy(&N);
854: MatSetNullSpace(jac,nullspace);
855: MatNullSpaceDestroy(&nullspace);
856: return(0);
857: }
861: /*
862: L_2Error - Integrate the L_2 error of our solution over each face
863: */
864: PetscErrorCode L_2Error(DM da, Vec fVec, double *error, AppCtx *user)
865: {
866: DMDALocalInfo info;
867: Vec fLocalVec;
868: Field **f;
869: Field u, uExact, uLocal[4];
870: PetscScalar hx, hy, hxhy, x, y, phi[3];
871: PetscInt i, j, q;
875: DMDAGetLocalInfo(da, &info);
876: DMGetLocalVector(da, &fLocalVec);
877: DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
878: DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
879: DMDAVecGetArray(da, fLocalVec, &f);
881: *error = 0.0;
882: hx = 1.0/(PetscReal)(info.mx-1);
883: hy = 1.0/(PetscReal)(info.my-1);
884: hxhy = hx*hy;
885: for(j = info.ys; j < info.ys+info.ym-1; j++) {
886: for(i = info.xs; i < info.xs+info.xm-1; i++) {
887: uLocal[0] = f[j][i];
888: uLocal[1] = f[j][i+1];
889: uLocal[2] = f[j+1][i+1];
890: uLocal[3] = f[j+1][i];
891: /* Lower element */
892: for(q = 0; q < 4; q++) {
893: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
894: phi[1] = quadPoints[q*2];
895: phi[2] = quadPoints[q*2+1];
896: u.u = uLocal[0].u*phi[0]+ uLocal[1].u*phi[1] + uLocal[3].u*phi[2];
897: u.v = uLocal[0].v*phi[0]+ uLocal[1].v*phi[1] + uLocal[3].v*phi[2];
898: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
899: x = (quadPoints[q*2] + i)*hx;
900: y = (quadPoints[q*2+1] + j)*hy;
901: ExactSolution(x, y, &uExact);
902: *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));
903: }
904: /* Upper element */
905: /*
906: The affine map from the lower to the upper is
908: / x_U \ = / -1 0 \ / x_L \ + / hx \
909: \ y_U / \ 0 -1 / \ y_L / \ hy /
910: */
911: for(q = 0; q < 4; q++) {
912: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
913: phi[1] = quadPoints[q*2];
914: phi[2] = quadPoints[q*2+1];
915: u.u = uLocal[2].u*phi[0]+ uLocal[3].u*phi[1] + uLocal[1].u*phi[2];
916: u.v = uLocal[2].v*phi[0]+ uLocal[3].v*phi[1] + uLocal[1].v*phi[2];
917: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
918: x = (1.0 - quadPoints[q*2] + i)*hx;
919: y = (1.0 - quadPoints[q*2+1] + j)*hy;
920: ExactSolution(x, y, &uExact);
921: *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));
922: }
923: }
924: }
926: DMDAVecRestoreArray(da, fLocalVec, &f);
927: /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
928: /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
929: DMRestoreLocalVector(da, &fLocalVec);
930: return(0);
931: }