Actual source code: ex7.c
petsc-3.5.4 2015-05-23
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 <petscdm.h>
42: #include <petscdmda.h>
43: #include <petscsnes.h>
45: /*
46: User-defined application context - contains data needed by the
47: application-provided call-back routines, FormJacobianLocal() and
48: FormFunctionLocal().
49: */
50: typedef struct {
51: PetscReal alpha; /* parameter controlling linearity */
52: PetscReal lambda; /* parameter controlling nonlinearity */
53: } AppCtx;
55: typedef struct {
56: PetscScalar u;
57: PetscScalar v;
58: PetscScalar p;
59: } Field;
61: static PetscScalar Kref[36] = { 0.5, 0.5, -0.5, 0, 0, -0.5,
62: 0.5, 0.5, -0.5, 0, 0, -0.5,
63: -0.5, -0.5, 0.5, 0, 0, 0.5,
64: 0, 0, 0, 0, 0, 0,
65: 0, 0, 0, 0, 0, 0,
66: -0.5, -0.5, 0.5, 0, 0, 0.5};
68: static PetscScalar Gradient[18] = {-0.1666667, -0.1666667, -0.1666667,
69: -0.1666667, -0.1666667, -0.1666667,
70: 0.1666667, 0.1666667, 0.1666667,
71: 0.0, 0.0, 0.0,
72: 0.0, 0.0, 0.0,
73: 0.1666667, 0.1666667, 0.1666667};
75: static PetscScalar Divergence[18] = {-0.1666667, 0.1666667, 0.0,
76: -0.1666667, 0.0, 0.1666667,
77: -0.1666667, 0.1666667, 0.0,
78: -0.1666667, 0.0, 0.1666667,
79: -0.1666667, 0.1666667, 0.0,
80: -0.1666667, 0.0, 0.1666667};
82: /* These are */
83: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
84: 0.07503111, 0.64494897,
85: 0.66639025, 0.15505103,
86: 0.28001992, 0.64494897};
87: static PetscScalar quadWeights[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
89: /*
90: User-defined routines
91: */
92: extern PetscErrorCode CreateNullSpace(DM, Vec*);
93: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
94: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,AppCtx*);
95: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,Field**,Mat,Mat,AppCtx*);
96: extern PetscErrorCode L_2Error(DM, Vec, PetscReal*, AppCtx*);
100: int main(int argc,char **argv)
101: {
102: DM da;
103: SNES snes; /* nonlinear solver */
104: AppCtx *user; /* user-defined work context */
105: PetscBag bag;
106: PetscInt its; /* iterations for convergence */
107: PetscMPIInt size;
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);
117: MPI_Comm_size(PETSC_COMM_WORLD,&size);
118: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example only works for one process.");
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Initialize problem parameters
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
124: PetscBagGetData(bag, (void**) &user);
125: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
126: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
127: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
128: PetscBagSetFromOptions(bag);
129: PetscOptionsGetReal(NULL,"-alpha",&user->alpha,NULL);
130: PetscOptionsGetReal(NULL,"-lambda",&user->lambda,NULL);
131: if (user->lambda > lambda_max || user->lambda < lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda %g is out of range [%g, %g]", (double)user->lambda, (double)lambda_min, (double)lambda_max);
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, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
142: DMDASetFieldName(da, 0, "ooblek");
143: DMSetApplicationContext(da,user);
144: SNESSetDM(snes, (DM) da);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Set the discretization functions
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,user);
150: DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,user);
151: SNESSetFromOptions(snes);
153: SNESSetComputeInitialGuess(snes, FormInitialGuess,NULL);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Solve nonlinear system
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: SNESSolve(snes,NULL,NULL);
159: SNESGetIterationNumber(snes,&its);
160: SNESGetConvergedReason(snes, &reason);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
164: DMDestroy(&da);
165: SNESGetDM(snes,&da);
166: SNESGetSolution(snes,&x);
167: L_2Error(da, x, &error, user);
168: PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %g\n", (double)error);
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Free work space. All PETSc objects should be destroyed when they
172: are no longer needed.
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: SNESDestroy(&snes);
175: PetscBagDestroy(&bag);
176: PetscFinalize();
177: return(0);
178: }
182: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, Field *u)
183: {
185: (*u).u = x;
186: (*u).v = -y;
187: (*u).p = 0.0;
188: return(0);
189: }
193: PetscErrorCode CreateNullSpace(DM da, Vec *N)
194: {
195: Field **x;
196: PetscInt xs,ys,xm,ym,i,j;
200: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
201: DMCreateGlobalVector(da,N);
202: DMDAVecGetArray(da, *N, &x);
203: for (j = ys; j < ys+ym; j++) {
204: for (i = xs; i < xs+xm; i++) {
205: x[j][i].u = 0.0;
206: x[j][i].v = 0.0;
207: x[j][i].p = 1.0;
208: }
209: }
210: DMDAVecRestoreArray(da, *N, &x);
211: return(0);
212: }
216: /*
217: FormInitialGuess - Forms initial approximation.
219: Input Parameters:
220: dm - The DM context
221: X - vector
223: Output Parameter:
224: X - vector
225: */
226: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
227: {
228: AppCtx *user;
229: PetscInt i,j,Mx,My,xs,ys,xm,ym;
231: PetscReal lambda,hx,hy;
232: PETSC_UNUSED PetscReal temp1;
233: Field **x;
234: DM da;
237: SNESGetDM(snes,&da);
238: DMGetApplicationContext(da,&user);
239: 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);
241: lambda = user->lambda;
242: hx = 1.0/(PetscReal)(Mx-1);
243: hy = 1.0/(PetscReal)(My-1);
244: if (lambda == 0.0) temp1 = 0.0;
245: else temp1 = lambda/(lambda + 1.0);
247: /*
248: Get a pointer to vector data.
249: - For default PETSc vectors, VecGetArray() returns a pointer to
250: the data array. Otherwise, the routine is implementation dependent.
251: - You MUST call VecRestoreArray() when you no longer need access to
252: the array.
253: */
254: DMDAVecGetArray(da,X,&x);
256: /*
257: Get local grid boundaries (for 2-dimensional DMDA):
258: xs, ys - starting grid indices (no ghost points)
259: xm, ym - widths of local grid (no ghost points)
261: */
262: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
264: /*
265: Compute initial guess over the locally owned part of the grid
266: */
267: for (j=ys; j<ys+ym; j++) {
268: for (i=xs; i<xs+xm; i++) {
269: #define CHECK_SOLUTION
270: #if defined(CHECK_SOLUTION)
271: ExactSolution(i*hx, j*hy, &x[j][i]);
272: #else
273: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
274: /* Boundary conditions are usually zero Dirichlet */
275: ExactSolution(i*hx, j*hy, &x[j][i]);
276: } else {
277: PetscReal temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
278: x[j][i].u = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
279: x[j][i].v = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
280: x[j][i].p = 1.0;
281: }
282: #endif
283: }
284: }
286: /*
287: Restore vector
288: */
289: DMDAVecRestoreArray(da,X,&x);
290: return(0);
291: }
295: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, Field r[])
296: {
297: Field rLocal[3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
298: PetscScalar phi[3] = {0.0, 0.0, 0.0};
299: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy;
300: Field res;
301: PetscInt q, k;
304: for (q = 0; q < 4; q++) {
305: PETSC_UNUSED PetscReal x, y;
306: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
307: phi[1] = quadPoints[q*2];
308: phi[2] = quadPoints[q*2+1];
309: if (isLower) {
310: x = xI + PetscAbsScalar(quadPoints[q*2])*hx;
311: y = yI + PetscAbsScalar(quadPoints[q*2+1])*hy;
312: } else {
313: x = xI + 1.0 - PetscAbsScalar(quadPoints[q*2])*hx;
314: y = yI + 1.0 - PetscAbsScalar(quadPoints[q*2+1])*hy;
315: }
316: res.u = quadWeights[q]*(0.0);
317: res.v = quadWeights[q]*(0.0);
318: res.p = quadWeights[q]*(0.0);
319: for (k = 0; k < 3; k++) {
320: rLocal[k].u += phi[k]*res.u;
321: rLocal[k].v += phi[k]*res.v;
322: rLocal[k].p += phi[k]*res.p;
323: }
324: }
325: for (k = 0; k < 3; k++) {
326: /*printf(" constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p); */
327: r[k].u += lambda*hxhy*rLocal[k].u;
328: r[k].v += lambda*hxhy*rLocal[k].v;
329: r[k].p += hxhy*rLocal[k].p;
330: }
331: return(0);
332: }
336: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[])
337: {
338: Field rLocal[3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
339: PetscScalar phi[3] = {0.0, 0.0, 0.0};
340: Field res;
341: PetscInt q;
344: for (q = 0; q < 4; q++) {
345: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
346: phi[1] = quadPoints[q*2];
347: phi[2] = quadPoints[q*2+1];
348: res.u = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
349: res.v = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);
351: rLocal[0].u += phi[0]*res.u;
352: rLocal[0].v += phi[0]*res.v;
353: rLocal[1].u += phi[1]*res.u;
354: rLocal[1].v += phi[1]*res.v;
355: rLocal[2].u += phi[2]*res.u;
356: rLocal[2].v += phi[2]*res.v;
357: }
358: r[0].u += lambda*rLocal[0].u;
359: r[0].v += lambda*rLocal[0].v;
360: r[1].u += lambda*rLocal[1].u;
361: r[1].v += lambda*rLocal[1].v;
362: r[2].u += lambda*rLocal[2].u;
363: r[2].v += lambda*rLocal[2].v;
364: return(0);
365: }
369: /*
370: FormFunctionLocal - Evaluates nonlinear function, F(x).
372: */
373: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
374: {
375: Field uLocal[3];
376: Field rLocal[3];
377: PetscScalar G[4];
378: Field uExact;
379: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
380: PetscInt i,j,k,l;
384: /* Naive Jacobian calculation:
386: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
387: \ 0 1/hy / \ 0 hy /
388: */
389: alpha = user->alpha;
390: lambda = user->lambda;
391: hx = 1.0/(PetscReal)(info->mx-1);
392: hy = 1.0/(PetscReal)(info->my-1);
393: sc = hx*hy*lambda;
394: hxhy = hx*hy;
395: detJInv = hxhy;
397: G[0] = (1.0/(hx*hx)) * detJInv;
398: G[1] = 0.0;
399: G[2] = G[1];
400: G[3] = (1.0/(hy*hy)) * detJInv;
401: for (k = 0; k < 4; k++) {
402: /* printf("G[%d] = %g\n", k, G[k]);*/
403: }
405: /* Zero the vector */
406: PetscMemzero((void*) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(Field));
407: /* Compute function over the locally owned part of the grid. For each
408: vertex (i,j), we consider the element below:
410: 2 (1) (0)
411: i,j+1 --- i+1,j+1
412: | \ |
413: | \ |
414: | \ |
415: | \ |
416: | \ |
417: i,j --- i+1,j
418: 0 1 (2)
420: and therefore we do not loop over the last vertex in each dimension.
421: */
422: for (j = info->ys; j < info->ys+info->ym-1; j++) {
423: for (i = info->xs; i < info->xs+info->xm-1; i++) {
424: /* Lower element */
425: uLocal[0] = x[j][i];
426: uLocal[1] = x[j][i+1];
427: uLocal[2] = x[j+1][i];
428: /* printf("Lower Solution ElementVector for (%d, %d)\n", i, j);*/
429: for (k = 0; k < 3; k++) {
430: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
431: }
432: for (k = 0; k < 3; k++) {
433: rLocal[k].u = 0.0;
434: rLocal[k].v = 0.0;
435: rLocal[k].p = 0.0;
436: for (l = 0; l < 3; l++) {
437: 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;
438: 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;
439: /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
440: /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
441: /* Gradient */
442: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
443: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
444: /* Divergence */
445: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
446: }
447: }
448: /* printf("Lower Laplacian ElementVector for (%d, %d)\n", i, j);*/
449: for (k = 0; k < 3; k++) {
450: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
451: }
452: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
453: /* printf("Lower Laplacian+Constant ElementVector for (%d, %d)\n", i, j);*/
454: for (k = 0; k < 3; k++) {
455: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
456: }
457: nonlinearResidual(0.0*sc, uLocal, rLocal);
458: /* printf("Lower Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
459: for (k = 0; k < 3; k++) {
460: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
461: }
462: f[j][i].u += rLocal[0].u;
463: f[j][i].v += rLocal[0].v;
464: f[j][i].p += rLocal[0].p;
465: f[j][i+1].u += rLocal[1].u;
466: f[j][i+1].v += rLocal[1].v;
467: f[j][i+1].p += rLocal[1].p;
468: f[j+1][i].u += rLocal[2].u;
469: f[j+1][i].v += rLocal[2].v;
470: f[j+1][i].p += rLocal[2].p;
471: /* Upper element */
472: uLocal[0] = x[j+1][i+1];
473: uLocal[1] = x[j+1][i];
474: uLocal[2] = x[j][i+1];
475: /* printf("Upper Solution ElementVector for (%d, %d)\n", i, j);*/
476: for (k = 0; k < 3; k++) {
477: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
478: }
479: for (k = 0; k < 3; k++) {
480: rLocal[k].u = 0.0;
481: rLocal[k].v = 0.0;
482: rLocal[k].p = 0.0;
483: for (l = 0; l < 3; l++) {
484: 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;
485: 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;
486: /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
487: /* Gradient */
488: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
489: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
490: /* Divergence */
491: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
492: }
493: }
494: /* printf("Upper Laplacian 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: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
499: /* printf("Upper Laplacian+Constant 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: nonlinearResidual(0.0*sc, uLocal, rLocal);
504: /* printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
505: for (k = 0; k < 3; k++) {
506: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
507: }
508: f[j+1][i+1].u += rLocal[0].u;
509: f[j+1][i+1].v += rLocal[0].v;
510: f[j+1][i+1].p += rLocal[0].p;
511: f[j+1][i].u += rLocal[1].u;
512: f[j+1][i].v += rLocal[1].v;
513: f[j+1][i].p += rLocal[1].p;
514: f[j][i+1].u += rLocal[2].u;
515: f[j][i+1].v += rLocal[2].v;
516: f[j][i+1].p += rLocal[2].p;
517: /* Boundary conditions */
518: if (i == 0 || j == 0) {
519: ExactSolution(i*hx, j*hy, &uExact);
521: f[j][i].u = x[j][i].u - uExact.u;
522: f[j][i].v = x[j][i].v - uExact.v;
523: f[j][i].p = x[j][i].p - uExact.p;
524: }
525: if ((i == info->mx-2) || (j == 0)) {
526: ExactSolution((i+1)*hx, j*hy, &uExact);
528: f[j][i+1].u = x[j][i+1].u - uExact.u;
529: f[j][i+1].v = x[j][i+1].v - uExact.v;
530: f[j][i+1].p = x[j][i+1].p - uExact.p;
531: }
532: if ((i == info->mx-2) || (j == info->my-2)) {
533: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
535: f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
536: f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
537: f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
538: }
539: if ((i == 0) || (j == info->my-2)) {
540: ExactSolution(i*hx, (j+1)*hy, &uExact);
542: f[j+1][i].u = x[j+1][i].u - uExact.u;
543: f[j+1][i].v = x[j+1][i].v - uExact.v;
544: f[j+1][i].p = x[j+1][i].p - uExact.p;
545: }
546: }
547: }
549: for (j = info->ys+info->ym-1; j >= info->ys; j--) {
550: for (i = info->xs; i < info->xs+info->xm; i++) {
551: /* printf("f[%d][%d] = (%g, %g, %g) ", j, i, f[j][i].u, f[j][i].v, f[j][i].p);*/
552: }
553: /*printf("\n");*/
554: }
555: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
556: return(0);
557: }
561: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[])
562: {
564: return(0);
565: }
569: /*
570: FormJacobianLocal - Evaluates Jacobian matrix.
571: */
572: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat A,Mat jac, AppCtx *user)
573: {
574: Field uLocal[4];
575: PetscScalar JLocal[144];
576: MatStencil rows[4*3], cols[4*3], ident;
577: PetscInt lowerRow[3] = {0, 1, 3};
578: PetscInt upperRow[3] = {2, 3, 1};
579: PetscInt hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
580: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
581: PetscInt i,j,k,l,numRows,dof = info->dof;
583: MatNullSpace nullspace;
584: Vec N;
587: alpha = user->alpha;
588: lambda = user->lambda;
589: hx = 1.0/(PetscReal)(info->mx-1);
590: hy = 1.0/(PetscReal)(info->my-1);
591: sc = hx*hy*lambda;
592: hxhy = hx*hy;
593: detJInv = hxhy;
595: G[0] = (1.0/(hx*hx)) * detJInv;
596: G[1] = 0.0;
597: G[2] = G[1];
598: G[3] = (1.0/(hy*hy)) * detJInv;
599: for (k = 0; k < 4; k++) {
600: /* printf("G[%d] = %g\n", k, G[k]);*/
601: }
603: MatZeroEntries(jac);
604: /*
605: Compute entries for the locally owned part of the Jacobian.
606: - Currently, all PETSc parallel matrix formats are partitioned by
607: contiguous chunks of rows across the processors.
608: - Each processor needs to insert only elements that it owns
609: locally (but any non-local elements will be sent to the
610: appropriate processor during matrix assembly).
611: - Here, we set all entries for a particular row at once.
612: - We can set matrix entries either using either
613: MatSetValuesLocal() or MatSetValues(), as discussed above.
614: */
615: #define NOT_PRES_BC 1
616: for (j=info->ys; j<info->ys+info->ym-1; j++) {
617: for (i=info->xs; i<info->xs+info->xm-1; i++) {
618: PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
619: numRows = 0;
620: /* Lower element */
621: uLocal[0] = x[j][i];
622: uLocal[1] = x[j][i+1];
623: uLocal[2] = x[j+1][i+1];
624: uLocal[3] = x[j+1][i];
625: /* i,j */
626: if (i == 0 || j == 0) {
627: hasLower[0] = 0;
629: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
630: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
632: ident.i = i; ident.j = j; ident.c = 0;
634: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
636: ident.i = i; ident.j = j; ident.c = 1;
638: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
639: #if defined(PRES_BC)
640: ident.i = i; ident.j = j; ident.c = 2;
642: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
643: #endif
644: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
645: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
646: } else {
647: hasLower[0] = 1;
648: velocityRows[0] = numRows;
649: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
650: numRows++;
651: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
652: numRows++;
653: #if defined(PRES_BC)
654: pressureRows[0] = numRows;
655: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
656: numRows++;
657: #endif
658: }
659: #if !defined(PRES_BC)
660: pressureRows[0] = numRows;
661: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
662: numRows++;
663: #endif
664: cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
665: cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
666: cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
667: /* i+1,j */
668: if ((i == info->mx-2) || (j == 0)) {
669: hasLower[1] = 0;
670: hasUpper[2] = 0;
672: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
673: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
675: ident.i = i+1; ident.j = j; ident.c = 0;
677: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
679: ident.i = i+1; ident.j = j; ident.c = 1;
681: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
682: #if defined(PRES_BC)
683: ident.i = i+1; ident.j = j; ident.c = 2;
685: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
686: #endif
687: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
688: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
689: } else {
690: hasLower[1] = 1;
691: hasUpper[2] = 1;
692: velocityRows[1] = numRows;
693: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
694: numRows++;
695: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
696: numRows++;
697: #if defined(PRES_BC)
698: pressureRows[1] = numRows;
699: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
700: numRows++;
701: #endif
702: }
703: #if !defined(PRES_BC)
704: pressureRows[1] = numRows;
705: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
706: numRows++;
707: #endif
708: cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
709: cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
710: cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
711: /* i+1,j+1 */
712: if ((i == info->mx-2) || (j == info->my-2)) {
713: hasUpper[0] = 0;
714: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
715: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
717: ident.i = i+1; ident.j = j+1; ident.c = 0;
719: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
721: ident.i = i+1; ident.j = j+1; ident.c = 1;
723: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
724: #if defined(PRES_BC)
725: ident.i = i+1; ident.j = j+1; ident.c = 2;
727: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
728: #endif
729: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
730: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
731: } else {
732: hasUpper[0] = 1;
733: velocityRows[2] = numRows;
734: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
735: numRows++;
736: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
737: numRows++;
738: #if defined(PRES_BC)
739: pressureRows[2] = numRows;
740: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
741: numRows++;
742: #endif
743: }
744: #if !defined(PRES_BC)
745: pressureRows[2] = numRows;
746: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
747: numRows++;
748: #endif
749: cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
750: cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
751: cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
752: /* i,j+1 */
753: if ((i == 0) || (j == info->my-2)) {
754: hasLower[2] = 0;
755: hasUpper[1] = 0;
756: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
757: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
759: ident.i = i; ident.j = j+1; ident.c = 0;
761: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
763: ident.i = i; ident.j = j+1; ident.c = 1;
765: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
766: #if defined(PRES_BC)
767: ident.i = i; ident.j = j+1; ident.c = 2;
769: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
770: #endif
771: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
772: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
773: } else {
774: hasLower[2] = 1;
775: hasUpper[1] = 1;
776: velocityRows[3] = numRows;
777: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
778: numRows++;
779: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
780: numRows++;
781: #if defined(PRES_BC)
782: pressureRows[3] = numRows;
783: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
784: numRows++;
785: #endif
786: }
787: #if !defined(PRES_BC)
788: pressureRows[3] = numRows;
789: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
790: numRows++;
791: #endif
792: cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
793: cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
794: cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;
796: /* Lower Element */
797: for (k = 0; k < 3; k++) {
798: #if defined(PRES_BC)
799: if (!hasLower[k]) continue;
800: #endif
801: for (l = 0; l < 3; l++) {
802: /* Divergence */
803: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
804: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
805: /* JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += Identity[k*3 + l]; */
806: }
807: if (!hasLower[k]) continue;
808: for (l = 0; l < 3; l++) {
809: /* Laplacian */
810: 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]);
811: 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]);
812: /* JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += Identity[k*3 + l]; */
813: /* JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += Identity[k*3 + l]; */
814: /* Gradient */
815: JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
816: JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
817: }
818: }
819: /* Upper Element */
820: for (k = 0; k < 3; k++) {
821: #if defined(PRES_BC)
822: if (!hasUpper[k]) continue;
823: #endif
824: for (l = 0; l < 3; l++) {
825: /* Divergence */
826: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
827: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
828: /* JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += Identity[k*3 + l]; */
829: }
830: if (!hasUpper[k]) continue;
831: for (l = 0; l < 3; l++) {
832: /* Laplacian */
833: 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]);
834: 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]);
835: /* Gradient */
836: JLocal[velocityRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
837: JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
838: }
839: }
841: nonlinearJacobian(-1.0*PetscAbsScalar(sc), uLocal, JLocal);
842: /* printf("Element matrix for (%d, %d)\n", i, j);*/
843: /* printf(" col ");*/
844: for (l = 0; l < 4*3; l++) {
845: /* printf("(%d, %d, %d) ", cols[l].i, cols[l].j, cols[l].c);*/
846: }
847: /* printf("\n");*/
848: for (k = 0; k < numRows; k++) {
849: /* printf("row (%d, %d, %d): ", rows[k].i, rows[k].j, rows[k].c);*/
850: for (l = 0; l < 4; l++) {
851: /* 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]);*/
852: }
853: /* printf("\n");*/
854: }
855: MatSetValuesStencil(jac,numRows,rows,4*dof,cols, JLocal,ADD_VALUES);
856: }
857: }
859: /*
860: Assemble matrix, using the 2-step process:
861: MatAssemblyBegin(), MatAssemblyEnd().
862: */
863: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
864: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
865: if (A != jac) {
866: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
867: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
868: }
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: }