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