Actual source code: ex7.c
petsc-3.7.5 2017-01-01
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,NULL,"-alpha",&user->alpha,NULL);
130: PetscOptionsGetReal(NULL,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;
230: PetscErrorCode ierr;
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: {
298: Field rLocal[3];
299: PetscScalar phi[3] = {0.0, 0.0, 0.0};
300: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy;
301: Field res;
302: PetscInt q, k;
305: PetscMemzero(&rLocal,sizeof(rLocal));
306: for (q = 0; q < 4; q++) {
307: PETSC_UNUSED PetscReal x, y;
308: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
309: phi[1] = quadPoints[q*2];
310: phi[2] = quadPoints[q*2+1];
311: if (isLower) {
312: x = xI + PetscAbsScalar(quadPoints[q*2])*hx;
313: y = yI + PetscAbsScalar(quadPoints[q*2+1])*hy;
314: } else {
315: x = xI + 1.0 - PetscAbsScalar(quadPoints[q*2])*hx;
316: y = yI + 1.0 - PetscAbsScalar(quadPoints[q*2+1])*hy;
317: }
318: res.u = quadWeights[q]*(0.0);
319: res.v = quadWeights[q]*(0.0);
320: res.p = quadWeights[q]*(0.0);
321: for (k = 0; k < 3; k++) {
322: rLocal[k].u += phi[k]*res.u;
323: rLocal[k].v += phi[k]*res.v;
324: rLocal[k].p += phi[k]*res.p;
325: }
326: }
327: for (k = 0; k < 3; k++) {
328: /*printf(" constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p); */
329: r[k].u += lambda*hxhy*rLocal[k].u;
330: r[k].v += lambda*hxhy*rLocal[k].v;
331: r[k].p += hxhy*rLocal[k].p;
332: }
333: return(0);
334: }
338: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[])
339: {
341: Field rLocal[3];
342: PetscScalar phi[3] = {0.0, 0.0, 0.0};
343: Field res;
344: PetscInt q;
347: PetscMemzero(&rLocal,sizeof(rLocal));
348: for (q = 0; q < 4; q++) {
349: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
350: phi[1] = quadPoints[q*2];
351: phi[2] = quadPoints[q*2+1];
352: res.u = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
353: res.v = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);
355: rLocal[0].u += phi[0]*res.u;
356: rLocal[0].v += phi[0]*res.v;
357: rLocal[1].u += phi[1]*res.u;
358: rLocal[1].v += phi[1]*res.v;
359: rLocal[2].u += phi[2]*res.u;
360: rLocal[2].v += phi[2]*res.v;
361: }
362: r[0].u += lambda*rLocal[0].u;
363: r[0].v += lambda*rLocal[0].v;
364: r[1].u += lambda*rLocal[1].u;
365: r[1].v += lambda*rLocal[1].v;
366: r[2].u += lambda*rLocal[2].u;
367: r[2].v += lambda*rLocal[2].v;
368: return(0);
369: }
373: /*
374: FormFunctionLocal - Evaluates nonlinear function, F(x).
376: */
377: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
378: {
379: Field uLocal[3];
380: Field rLocal[3];
381: PetscScalar G[4];
382: Field uExact;
383: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
384: PetscInt i,j,k,l;
388: /* Naive Jacobian calculation:
390: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
391: \ 0 1/hy / \ 0 hy /
392: */
393: alpha = user->alpha;
394: lambda = user->lambda;
395: hx = 1.0/(PetscReal)(info->mx-1);
396: hy = 1.0/(PetscReal)(info->my-1);
397: sc = hx*hy*lambda;
398: hxhy = hx*hy;
399: detJInv = hxhy;
401: G[0] = (1.0/(hx*hx)) * detJInv;
402: G[1] = 0.0;
403: G[2] = G[1];
404: G[3] = (1.0/(hy*hy)) * detJInv;
405: for (k = 0; k < 4; k++) {
406: /* printf("G[%d] = %g\n", k, G[k]);*/
407: }
409: /* Zero the vector */
410: PetscMemzero((void*) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(Field));
411: /* Compute function over the locally owned part of the grid. For each
412: vertex (i,j), we consider the element below:
414: 2 (1) (0)
415: i,j+1 --- i+1,j+1
416: | \ |
417: | \ |
418: | \ |
419: | \ |
420: | \ |
421: i,j --- i+1,j
422: 0 1 (2)
424: and therefore we do not loop over the last vertex in each dimension.
425: */
426: for (j = info->ys; j < info->ys+info->ym-1; j++) {
427: for (i = info->xs; i < info->xs+info->xm-1; i++) {
428: /* Lower element */
429: uLocal[0] = x[j][i];
430: uLocal[1] = x[j][i+1];
431: uLocal[2] = x[j+1][i];
432: /* printf("Lower Solution ElementVector for (%d, %d)\n", i, j);*/
433: for (k = 0; k < 3; k++) {
434: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
435: }
436: for (k = 0; k < 3; k++) {
437: rLocal[k].u = 0.0;
438: rLocal[k].v = 0.0;
439: rLocal[k].p = 0.0;
440: for (l = 0; l < 3; l++) {
441: 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;
442: 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;
443: /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
444: /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
445: /* Gradient */
446: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
447: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
448: /* Divergence */
449: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
450: }
451: }
452: /* printf("Lower Laplacian 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: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
457: /* printf("Lower Laplacian+Constant 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: nonlinearResidual(0.0*sc, uLocal, rLocal);
462: /* printf("Lower Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
463: for (k = 0; k < 3; k++) {
464: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
465: }
466: f[j][i].u += rLocal[0].u;
467: f[j][i].v += rLocal[0].v;
468: f[j][i].p += rLocal[0].p;
469: f[j][i+1].u += rLocal[1].u;
470: f[j][i+1].v += rLocal[1].v;
471: f[j][i+1].p += rLocal[1].p;
472: f[j+1][i].u += rLocal[2].u;
473: f[j+1][i].v += rLocal[2].v;
474: f[j+1][i].p += rLocal[2].p;
475: /* Upper element */
476: uLocal[0] = x[j+1][i+1];
477: uLocal[1] = x[j+1][i];
478: uLocal[2] = x[j][i+1];
479: /* printf("Upper Solution ElementVector for (%d, %d)\n", i, j);*/
480: for (k = 0; k < 3; k++) {
481: /* printf(" uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
482: }
483: for (k = 0; k < 3; k++) {
484: rLocal[k].u = 0.0;
485: rLocal[k].v = 0.0;
486: rLocal[k].p = 0.0;
487: for (l = 0; l < 3; l++) {
488: 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;
489: 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;
490: /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
491: /* Gradient */
492: rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
493: rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
494: /* Divergence */
495: rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
496: }
497: }
498: /* printf("Upper Laplacian 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: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
503: /* printf("Upper Laplacian+Constant 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: nonlinearResidual(0.0*sc, uLocal, rLocal);
508: /* printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
509: for (k = 0; k < 3; k++) {
510: /* printf(" rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
511: }
512: f[j+1][i+1].u += rLocal[0].u;
513: f[j+1][i+1].v += rLocal[0].v;
514: f[j+1][i+1].p += rLocal[0].p;
515: f[j+1][i].u += rLocal[1].u;
516: f[j+1][i].v += rLocal[1].v;
517: f[j+1][i].p += rLocal[1].p;
518: f[j][i+1].u += rLocal[2].u;
519: f[j][i+1].v += rLocal[2].v;
520: f[j][i+1].p += rLocal[2].p;
521: /* Boundary conditions */
522: if (i == 0 || j == 0) {
523: ExactSolution(i*hx, j*hy, &uExact);
525: f[j][i].u = x[j][i].u - uExact.u;
526: f[j][i].v = x[j][i].v - uExact.v;
527: f[j][i].p = x[j][i].p - uExact.p;
528: }
529: if ((i == info->mx-2) || (j == 0)) {
530: ExactSolution((i+1)*hx, j*hy, &uExact);
532: f[j][i+1].u = x[j][i+1].u - uExact.u;
533: f[j][i+1].v = x[j][i+1].v - uExact.v;
534: f[j][i+1].p = x[j][i+1].p - uExact.p;
535: }
536: if ((i == info->mx-2) || (j == info->my-2)) {
537: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
539: f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
540: f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
541: f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
542: }
543: if ((i == 0) || (j == info->my-2)) {
544: 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[])
566: {
568: return(0);
569: }
573: /*
574: FormJacobianLocal - Evaluates Jacobian matrix.
575: */
576: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat A,Mat jac, AppCtx *user)
577: {
578: Field uLocal[4];
579: PetscScalar JLocal[144];
580: MatStencil rows[4*3], cols[4*3], ident;
581: PetscInt lowerRow[3] = {0, 1, 3};
582: PetscInt upperRow[3] = {2, 3, 1};
583: PetscInt hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
584: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
585: PetscInt i,j,k,l,numRows,dof = info->dof;
587: MatNullSpace nullspace;
588: Vec N;
591: alpha = user->alpha;
592: lambda = user->lambda;
593: hx = 1.0/(PetscReal)(info->mx-1);
594: hy = 1.0/(PetscReal)(info->my-1);
595: sc = hx*hy*lambda;
596: hxhy = hx*hy;
597: detJInv = hxhy;
599: G[0] = (1.0/(hx*hx)) * detJInv;
600: G[1] = 0.0;
601: G[2] = G[1];
602: G[3] = (1.0/(hy*hy)) * detJInv;
603: for (k = 0; k < 4; k++) {
604: /* printf("G[%d] = %g\n", k, G[k]);*/
605: }
607: MatZeroEntries(jac);
608: /*
609: Compute entries for the locally owned part of the Jacobian.
610: - Currently, all PETSc parallel matrix formats are partitioned by
611: contiguous chunks of rows across the processors.
612: - Each processor needs to insert only elements that it owns
613: locally (but any non-local elements will be sent to the
614: appropriate processor during matrix assembly).
615: - Here, we set all entries for a particular row at once.
616: - We can set matrix entries either using either
617: MatSetValuesLocal() or MatSetValues(), as discussed above.
618: */
619: #define NOT_PRES_BC 1
620: for (j=info->ys; j<info->ys+info->ym-1; j++) {
621: for (i=info->xs; i<info->xs+info->xm-1; i++) {
622: PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
623: numRows = 0;
624: /* Lower element */
625: uLocal[0] = x[j][i];
626: uLocal[1] = x[j][i+1];
627: uLocal[2] = x[j+1][i+1];
628: uLocal[3] = x[j+1][i];
629: /* i,j */
630: if (i == 0 || j == 0) {
631: hasLower[0] = 0;
633: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
634: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
636: ident.i = i; ident.j = j; ident.c = 0;
638: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
640: ident.i = i; ident.j = j; ident.c = 1;
642: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
643: #if defined(PRES_BC)
644: ident.i = i; ident.j = j; ident.c = 2;
646: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
647: #endif
648: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
649: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
650: } else {
651: hasLower[0] = 1;
652: velocityRows[0] = numRows;
653: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
654: numRows++;
655: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
656: numRows++;
657: #if defined(PRES_BC)
658: pressureRows[0] = numRows;
659: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
660: numRows++;
661: #endif
662: }
663: #if !defined(PRES_BC)
664: pressureRows[0] = numRows;
665: rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
666: numRows++;
667: #endif
668: cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
669: cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
670: cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
671: /* i+1,j */
672: if ((i == info->mx-2) || (j == 0)) {
673: hasLower[1] = 0;
674: hasUpper[2] = 0;
676: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
677: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
679: ident.i = i+1; ident.j = j; ident.c = 0;
681: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
683: ident.i = i+1; ident.j = j; ident.c = 1;
685: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
686: #if defined(PRES_BC)
687: ident.i = i+1; ident.j = j; ident.c = 2;
689: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
690: #endif
691: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
692: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
693: } else {
694: hasLower[1] = 1;
695: hasUpper[2] = 1;
696: velocityRows[1] = numRows;
697: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
698: numRows++;
699: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
700: numRows++;
701: #if defined(PRES_BC)
702: pressureRows[1] = numRows;
703: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
704: numRows++;
705: #endif
706: }
707: #if !defined(PRES_BC)
708: pressureRows[1] = numRows;
709: rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
710: numRows++;
711: #endif
712: cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
713: cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
714: cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
715: /* i+1,j+1 */
716: if ((i == info->mx-2) || (j == info->my-2)) {
717: hasUpper[0] = 0;
718: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
719: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
721: ident.i = i+1; ident.j = j+1; ident.c = 0;
723: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
725: ident.i = i+1; ident.j = j+1; ident.c = 1;
727: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
728: #if defined(PRES_BC)
729: ident.i = i+1; ident.j = j+1; ident.c = 2;
731: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
732: #endif
733: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
734: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
735: } else {
736: hasUpper[0] = 1;
737: velocityRows[2] = numRows;
738: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
739: numRows++;
740: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
741: numRows++;
742: #if defined(PRES_BC)
743: pressureRows[2] = numRows;
744: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
745: numRows++;
746: #endif
747: }
748: #if !defined(PRES_BC)
749: pressureRows[2] = numRows;
750: rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
751: numRows++;
752: #endif
753: cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
754: cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
755: cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
756: /* i,j+1 */
757: if ((i == 0) || (j == info->my-2)) {
758: hasLower[2] = 0;
759: hasUpper[1] = 0;
760: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
761: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
763: ident.i = i; ident.j = j+1; ident.c = 0;
765: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
767: ident.i = i; ident.j = j+1; ident.c = 1;
769: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
770: #if defined(PRES_BC)
771: ident.i = i; ident.j = j+1; ident.c = 2;
773: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
774: #endif
775: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
776: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
777: } else {
778: hasLower[2] = 1;
779: hasUpper[1] = 1;
780: velocityRows[3] = numRows;
781: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
782: numRows++;
783: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
784: numRows++;
785: #if defined(PRES_BC)
786: pressureRows[3] = numRows;
787: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
788: numRows++;
789: #endif
790: }
791: #if !defined(PRES_BC)
792: pressureRows[3] = numRows;
793: rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
794: numRows++;
795: #endif
796: cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
797: cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
798: cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;
800: /* Lower Element */
801: for (k = 0; k < 3; k++) {
802: #if defined(PRES_BC)
803: if (!hasLower[k]) continue;
804: #endif
805: for (l = 0; l < 3; l++) {
806: /* Divergence */
807: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
808: JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
809: /* JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += Identity[k*3 + l]; */
810: }
811: if (!hasLower[k]) continue;
812: for (l = 0; l < 3; l++) {
813: /* Laplacian */
814: 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]);
815: 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]);
816: /* JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += Identity[k*3 + l]; */
817: /* JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += Identity[k*3 + l]; */
818: /* Gradient */
819: JLocal[velocityRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
820: JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
821: }
822: }
823: /* Upper Element */
824: for (k = 0; k < 3; k++) {
825: #if defined(PRES_BC)
826: if (!hasUpper[k]) continue;
827: #endif
828: for (l = 0; l < 3; l++) {
829: /* Divergence */
830: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
831: JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
832: /* JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += Identity[k*3 + l]; */
833: }
834: if (!hasUpper[k]) continue;
835: for (l = 0; l < 3; l++) {
836: /* Laplacian */
837: 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]);
838: 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]);
839: /* Gradient */
840: JLocal[velocityRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
841: JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
842: }
843: }
845: nonlinearJacobian(-1.0*PetscAbsScalar(sc), uLocal, JLocal);
846: /* printf("Element matrix for (%d, %d)\n", i, j);*/
847: /* printf(" col ");*/
848: for (l = 0; l < 4*3; l++) {
849: /* printf("(%d, %d, %d) ", cols[l].i, cols[l].j, cols[l].c);*/
850: }
851: /* printf("\n");*/
852: for (k = 0; k < numRows; k++) {
853: /* printf("row (%d, %d, %d): ", rows[k].i, rows[k].j, rows[k].c);*/
854: for (l = 0; l < 4; l++) {
855: /* 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]);*/
856: }
857: /* printf("\n");*/
858: }
859: MatSetValuesStencil(jac,numRows,rows,4*dof,cols, JLocal,ADD_VALUES);
860: }
861: }
863: /*
864: Assemble matrix, using the 2-step process:
865: MatAssemblyBegin(), MatAssemblyEnd().
866: */
867: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
868: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
869: if (A != jac) {
870: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
871: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
872: }
874: /*
875: Tell the matrix we will never add a new nonzero location to the
876: matrix. If we do, it will generate an error.
877: */
878: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
880: CreateNullSpace(info->da,&N);
881: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&N,&nullspace);
882: VecDestroy(&N);
883: MatSetNullSpace(A,nullspace);
884: MatNullSpaceDestroy(&nullspace);
885: return(0);
886: }
890: /*
891: L_2Error - Integrate the L_2 error of our solution over each face
892: */
893: PetscErrorCode L_2Error(DM da, Vec fVec, PetscReal *error, AppCtx *user)
894: {
895: DMDALocalInfo info;
896: Vec fLocalVec;
897: Field **f;
898: Field u, uExact, uLocal[4];
899: PetscScalar hx, hy, hxhy, x, y, phi[3];
900: PetscInt i, j, q;
904: DMDAGetLocalInfo(da, &info);
905: DMGetLocalVector(da, &fLocalVec);
906: DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
907: DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
908: DMDAVecGetArray(da, fLocalVec, &f);
910: *error = 0.0;
911: hx = 1.0/(PetscReal)(info.mx-1);
912: hy = 1.0/(PetscReal)(info.my-1);
913: hxhy = hx*hy;
914: for (j = info.ys; j < info.ys+info.ym-1; j++) {
915: for (i = info.xs; i < info.xs+info.xm-1; i++) {
916: uLocal[0] = f[j][i];
917: uLocal[1] = f[j][i+1];
918: uLocal[2] = f[j+1][i+1];
919: uLocal[3] = f[j+1][i];
920: /* Lower element */
921: for (q = 0; q < 4; q++) {
922: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
923: phi[1] = quadPoints[q*2];
924: phi[2] = quadPoints[q*2+1];
925: u.u = uLocal[0].u*phi[0]+ uLocal[1].u*phi[1] + uLocal[3].u*phi[2];
926: u.v = uLocal[0].v*phi[0]+ uLocal[1].v*phi[1] + uLocal[3].v*phi[2];
927: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
928: x = (quadPoints[q*2] + (PetscReal)i)*hx;
929: y = (quadPoints[q*2+1] + (PetscReal)j)*hy;
930: ExactSolution(PetscAbsScalar(x), PetscAbsScalar(y), &uExact);
931: *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)));
932: }
933: /* Upper element */
934: /*
935: The affine map from the lower to the upper is
937: / x_U \ = / -1 0 \ / x_L \ + / hx \
938: \ y_U / \ 0 -1 / \ y_L / \ hy /
939: */
940: for (q = 0; q < 4; q++) {
941: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
942: phi[1] = quadPoints[q*2];
943: phi[2] = quadPoints[q*2+1];
944: u.u = uLocal[2].u*phi[0]+ uLocal[3].u*phi[1] + uLocal[1].u*phi[2];
945: u.v = uLocal[2].v*phi[0]+ uLocal[3].v*phi[1] + uLocal[1].v*phi[2];
946: u.p = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
947: x = (1.0 - quadPoints[q*2] + (PetscReal)i)*hx;
948: y = (1.0 - quadPoints[q*2+1] + (PetscReal)j)*hy;
949: ExactSolution(PetscAbsScalar(x), PetscAbsScalar(y), &uExact);
950: *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)));
951: }
952: }
953: }
955: DMDAVecRestoreArray(da, fLocalVec, &f);
956: DMRestoreLocalVector(da, &fLocalVec);
957: return(0);
958: }