Actual source code: ex8.c
petsc-3.3-p7 2013-05-11
1: /* Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options] */
3: static char help[] = "Nonlinear PDE in 2d.\n\
4: We solve the Bratu equation in a 2D rectangular\n\
5: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";
7: /*T
8: Concepts: SNES^parallel Bratu example
9: Concepts: DMDA^using distributed arrays;
10: Processors: n
11: T*/
13: /* ------------------------------------------------------------------------
15: The Bratu equation is given by the partial differential equation
16:
17: -alpha*Laplacian u + lambda*e^u = f, 0 < x,y < 1,
18:
19: with boundary conditions
20:
21: u = 0 for x = 0, x = 1, y = 0, y = 1.
22:
23: A linear 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: static PetscScalar Kref[36] = { 0.5, 0.5, -0.5, 0, 0, -0.5,
55: 0.5, 0.5, -0.5, 0, 0, -0.5,
56: -0.5, -0.5, 0.5, 0, 0, 0.5,
57: 0, 0, 0, 0, 0, 0,
58: 0, 0, 0, 0, 0, 0,
59: -0.5, -0.5, 0.5, 0, 0, 0.5};
61: /* These are */
62: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
63: 0.07503111, 0.64494897,
64: 0.66639025, 0.15505103,
65: 0.28001992, 0.64494897};
66: static PetscScalar quadWeights[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
68: /*
69: User-defined routines
70: */
71: extern PetscErrorCode FormInitialGuess(DM,Vec);
72: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
73: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*);
74: extern PetscErrorCode L_2Error(DM, Vec, double *, AppCtx *);
78: int main(int argc,char **argv)
79: {
80: DM da;
81: SNES snes; /* nonlinear solver */
82: AppCtx *user; /* user-defined work context */
83: PetscBag bag;
84: PetscInt its; /* iterations for convergence */
85: SNESConvergedReason reason;
86: PetscErrorCode ierr;
87: PetscReal lambda_max = 6.81, lambda_min = 0.0, error;
88: Vec x;
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Initialize program
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscInitialize(&argc,&argv,(char *)0,help);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Initialize problem parameters
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
99: PetscBagGetData(bag, (void **) &user);
100: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
101: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
102: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
103: PetscBagSetFromOptions(bag);
104: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
105: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
106: if (user->lambda > lambda_max || user->lambda < lambda_min) {
107: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
108: }
110: SNESCreate(PETSC_COMM_WORLD,&snes);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create distributed array (DA) to manage parallel grid and vectors
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
116: 1,1,PETSC_NULL,PETSC_NULL,&da);
117: DMDASetFieldName(da, 0, "ooblek");
118: DMSetApplicationContext(da,&user);
119: SNESSetDM(snes, (DM) da);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Set the discretization functions
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
125: DMDASetLocalJacobian(da,(DMDALocalFunction1)FormJacobianLocal);
126: DMSetInitialGuess(da,FormInitialGuess);
127: SNESSetFromOptions(snes);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Solve nonlinear system
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: SNESSolve(snes,PETSC_NULL,PETSC_NULL);
133: SNESGetIterationNumber(snes,&its);
134: SNESGetConvergedReason(snes, &reason);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
138: DMDestroy(&da);
139: SNESGetDM(snes,&da);
140: SNESGetSolution(snes,&x);
141: L_2Error(da, x, &error, user);
142: PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Free work space. All PETSc objects should be destroyed when they
146: are no longer needed.
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: SNESDestroy(&snes);
149: PetscBagDestroy(&bag);
150: PetscFinalize();
151: return(0);
152: }
156: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
157: {
159: *u = x*x;
160: return(0);
161: }
165: PetscErrorCode FormInitialGuess(DM da,Vec X)
166: {
167: AppCtx *user;
168: PetscInt i,j,Mx,My,xs,ys,xm,ym;
170: PetscReal lambda,temp1,temp,hx,hy;
171: PetscScalar **x;
174: DMGetApplicationContext(da,&user);
175: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
176: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
178: lambda = user->lambda;
179: hx = 1.0/(PetscReal)(Mx-1);
180: hy = 1.0/(PetscReal)(My-1);
181: if (lambda == 0.0) {
182: temp1 = 0.0;
183: } else {
184: temp1 = lambda/(lambda + 1.0);
185: }
187: /*
188: Get a pointer to vector data.
189: - For default PETSc vectors, VecGetArray() returns a pointer to
190: the data array. Otherwise, the routine is implementation dependent.
191: - You MUST call VecRestoreArray() when you no longer need access to
192: the array.
193: */
194: DMDAVecGetArray(da,X,&x);
196: /*
197: Get local grid boundaries (for 2-dimensional DMDA):
198: xs, ys - starting grid indices (no ghost points)
199: xm, ym - widths of local grid (no ghost points)
201: */
202: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
204: /*
205: Compute initial guess over the locally owned part of the grid
206: */
207: for (j=ys; j<ys+ym; j++) {
208: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
209: for (i=xs; i<xs+xm; i++) {
211: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
212: /* boundary conditions are all zero Dirichlet */
213: x[j][i] = 0.0;
214: } else {
215: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
216: }
217: }
218: }
220: DMDAVecRestoreArray(da,X,&x);
221: return(0);
222: }
226: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
227: {
228: PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
229: PetscScalar phi[3] = {0.0, 0.0, 0.0};
230: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
231: PetscScalar res;
232: PetscInt q, k;
235: for(q = 0; q < 4; q++) {
236: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
237: phi[1] = quadPoints[q*2];
238: phi[2] = quadPoints[q*2+1];
239: /* These are currently wrong */
240: x = xI + quadPoints[q*2]*hx;
241: y = yI + quadPoints[q*2+1]*hy;
242: res = quadWeights[q]*(2.0);
243: for(k = 0; k < 3; k++) {
244: rLocal[k] += phi[k]*res;
245: }
246: }
247: for(k = 0; k < 3; k++) {
248: printf(" constLocal[%d] = %g\n", k, lambda*hxhy*rLocal[k]);
249: r[k] += lambda*hxhy*rLocal[k];
250: }
251: return(0);
252: }
256: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[])
257: {
258: PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
259: PetscScalar phi[3] = {0.0, 0.0, 0.0};
260: PetscScalar res;
261: PetscInt q;
264: for(q = 0; q < 4; q++) {
265: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
266: phi[1] = quadPoints[q*2];
267: phi[2] = quadPoints[q*2+1];
268: res = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
269: rLocal[0] += phi[0]*res;
270: rLocal[1] += phi[1]*res;
271: rLocal[2] += phi[2]*res;
272: }
273: r[0] += lambda*rLocal[0];
274: r[1] += lambda*rLocal[1];
275: r[2] += lambda*rLocal[2];
276: return(0);
277: }
281: /*
282: FormFunctionLocal - Evaluates nonlinear function, F(x).
284: Process adiC(36): FormFunctionLocal
286: */
287: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
288: {
289: PetscScalar uLocal[3];
290: PetscScalar rLocal[3];
291: PetscScalar G[4];
292: PetscScalar uExact;
293: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
294: PetscInt i,j,k,l;
299: /* Naive Jacobian calculation:
301: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
302: \ 0 1/hy / \ 0 hy /
303: */
304: alpha = user->alpha;
305: lambda = user->lambda;
306: hx = 1.0/(PetscReal)(info->mx-1);
307: hy = 1.0/(PetscReal)(info->my-1);
308: sc = hx*hy*lambda;
309: hxhy = hx*hy;
310: detJInv = hxhy;
311: G[0] = (1.0/(hx*hx)) * detJInv;
312: G[1] = 0.0;
313: G[2] = G[1];
314: G[3] = (1.0/(hy*hy)) * detJInv;
315: for(k = 0; k < 4; k++) {
316: printf("G[%d] = %g\n", k, G[k]);
317: }
319: /* Zero the vector */
320: PetscMemzero((void *) &(f[info->ys][info->xs]), info->xm*info->ym*sizeof(PetscScalar));
321: /* Compute function over the locally owned part of the grid. For each
322: vertex (i,j), we consider the element below:
324: 2 (1) (0)
325: i,j+1 --- i+1,j+1
326: | \ |
327: | \ |
328: | \ |
329: | \ |
330: | \ |
331: i,j --- i+1,j
332: 0 1 (2)
334: and therefore we do not loop over the last vertex in each dimension.
335: */
336: for(j = info->ys; j < info->ys+info->ym-1; j++) {
337: for(i = info->xs; i < info->xs+info->xm-1; i++) {
338: /* Lower element */
339: uLocal[0] = x[j][i];
340: uLocal[1] = x[j][i+1];
341: uLocal[2] = x[j+1][i];
342: printf("Solution ElementVector for (%d, %d)\n", i, j);
343: for(k = 0; k < 3; k++) {
344: printf(" uLocal[%d] = %g\n", k, uLocal[k]);
345: }
346: for(k = 0; k < 3; k++) {
347: rLocal[k] = 0.0;
348: for(l = 0; l < 3; l++) {
349: rLocal[k] += (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];
350: }
351: rLocal[k] *= alpha;
352: }
353: printf("Laplacian ElementVector for (%d, %d)\n", i, j);
354: for(k = 0; k < 3; k++) {
355: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
356: }
357: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
358: printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
359: for(k = 0; k < 3; k++) {
360: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
361: }
362: nonlinearResidual(0.0*sc, uLocal, rLocal);
363: printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
364: for(k = 0; k < 3; k++) {
365: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
366: }
367: f[j][i] += rLocal[0];
368: f[j][i+1] += rLocal[1];
369: f[j+1][i] += rLocal[2];
370: /* Upper element */
371: uLocal[0] = x[j+1][i+1];
372: uLocal[1] = x[j+1][i];
373: uLocal[2] = x[j][i+1];
374: printf("Solution ElementVector for (%d, %d)\n", i, j);
375: for(k = 0; k < 3; k++) {
376: printf(" uLocal[%d] = %g\n", k, uLocal[k]);
377: }
378: for(k = 0; k < 3; k++) {
379: rLocal[k] = 0.0;
380: for(l = 0; l < 3; l++) {
381: rLocal[k] += (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];
382: }
383: rLocal[k] *= alpha;
384: }
385: printf("Laplacian ElementVector for (%d, %d)\n", i, j);
386: for(k = 0; k < 3; k++) {
387: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
388: }
389: constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
390: printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
391: for(k = 0; k < 3; k++) {
392: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
393: }
394: nonlinearResidual(0.0*sc, uLocal, rLocal);
395: printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
396: for(k = 0; k < 3; k++) {
397: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
398: }
399: f[j+1][i+1] += rLocal[0];
400: f[j+1][i] += rLocal[1];
401: f[j][i+1] += rLocal[2];
402: /* Boundary conditions */
403: if (i == 0 || j == 0) {
404: ExactSolution(i*hx, j*hy, &uExact);
405: f[j][i] = x[j][i] - uExact;
406: }
407: if ((i == info->mx-2) || (j == 0)) {
408: ExactSolution((i+1)*hx, j*hy, &uExact);
409: f[j][i+1] = x[j][i+1] - uExact;
410: }
411: if ((i == info->mx-2) || (j == info->my-2)) {
412: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
413: f[j+1][i+1] = x[j+1][i+1] - uExact;
414: }
415: if ((i == 0) || (j == info->my-2)) {
416: ExactSolution(i*hx, (j+1)*hy, &uExact);
417: f[j+1][i] = x[j+1][i] - uExact;
418: }
419: }
420: }
422: for(j = info->ys+info->ym-1; j >= info->ys; j--) {
423: for(i = info->xs; i < info->xs+info->xm; i++) {
424: printf("f[%d][%d] = %g ", j, i, f[j][i]);
425: }
426: printf("\n");
427: }
428: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
429: return(0);
430: }
434: PetscErrorCode nonlinearJacobian(PetscReal lambda, PetscScalar u[], PetscScalar J[]) {
436: return(0);
437: }
441: /*
442: FormJacobianLocal - Evaluates Jacobian matrix.
443: */
444: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
445: {
446: PetscScalar JLocal[16], uLocal[4];
447: MatStencil rows[4], cols[4], ident;
448: PetscInt lowerRow[3] = {0, 1, 3};
449: PetscInt upperRow[3] = {2, 3, 1};
450: PetscInt hasLower[3], hasUpper[3], localRows[4];
451: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
452: PetscInt i,j,k,l,numRows;
456: alpha = user->alpha;
457: lambda = user->lambda;
458: hx = 1.0/(PetscReal)(info->mx-1);
459: hy = 1.0/(PetscReal)(info->my-1);
460: sc = hx*hy*lambda;
461: hxhy = hx*hy;
462: detJInv = hxhy;
463: G[0] = (1.0/(hx*hx)) * detJInv;
464: G[1] = 0.0;
465: G[2] = G[1];
466: G[3] = (1.0/(hy*hy)) * detJInv;
467: for(k = 0; k < 4; k++) {
468: printf("G[%d] = %g\n", k, G[k]);
469: }
471: MatZeroEntries(jac);
472: /*
473: Compute entries for the locally owned part of the Jacobian.
474: - Currently, all PETSc parallel matrix formats are partitioned by
475: contiguous chunks of rows across the processors.
476: - Each processor needs to insert only elements that it owns
477: locally (but any non-local elements will be sent to the
478: appropriate processor during matrix assembly).
479: - Here, we set all entries for a particular row at once.
480: - We can set matrix entries either using either
481: MatSetValuesLocal() or MatSetValues(), as discussed above.
482: */
483: for (j=info->ys; j<info->ys+info->ym-1; j++) {
484: for (i=info->xs; i<info->xs+info->xm-1; i++) {
485: PetscMemzero(JLocal, 16 * sizeof(PetscScalar));
486: numRows = 0;
487: /* Lower element */
488: uLocal[0] = x[j][i];
489: uLocal[1] = x[j][i+1];
490: uLocal[2] = x[j+1][i+1];
491: uLocal[3] = x[j+1][i];
492: /* i,j */
493: if (i == 0 || j == 0) {
494: hasLower[0] = 0;
495: ident.i = i; ident.j = j;
496: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
497: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
498: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
499: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
500: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
501: } else {
502: hasLower[0] = 1;
503: localRows[0] = numRows;
504: rows[numRows].i = i; rows[numRows].j = j;
505: numRows++;
506: }
507: cols[0].i = i; cols[0].j = j;
508: /* i+1,j */
509: if ((i == info->mx-2) || (j == 0)) {
510: hasLower[1] = 0;
511: hasUpper[2] = 0;
512: ident.i = i+1; ident.j = j;
513: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
514: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
515: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
516: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
517: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
518: } else {
519: localRows[1] = numRows;
520: hasLower[1] = 1;
521: hasUpper[2] = 1;
522: localRows[1] = numRows;
523: rows[numRows].i = i+1; rows[numRows].j = j;
524: numRows++;
525: }
526: cols[1].i = i+1; cols[1].j = j;
527: /* i+1,j+1 */
528: if ((i == info->mx-2) || (j == info->my-2)) {
529: hasUpper[0] = 0;
530: ident.i = i+1; ident.j = j+1;
531: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
532: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
533: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
534: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
535: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
536: } else {
537: hasUpper[0] = 1;
538: localRows[2] = numRows;
539: rows[numRows].i = i+1; rows[numRows].j = j+1;
540: numRows++;
541: }
542: cols[2].i = i+1; cols[2].j = j+1;
543: /* i,j+1 */
544: if ((i == 0) || (j == info->my-2)) {
545: hasLower[2] = 0;
546: hasUpper[1] = 0;
547: ident.i = i; ident.j = j+1;
548: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
549: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
550: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
551: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
552: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
553: } else {
554: hasLower[2] = 1;
555: hasUpper[1] = 1;
556: localRows[3] = numRows;
557: rows[numRows].i = i; rows[numRows].j = j+1;
558: numRows++;
559: }
560: cols[3].i = i; cols[3].j = j+1;
562: /* Lower Element */
563: for(k = 0; k < 3; k++) {
564: if (!hasLower[k]) continue;
565: for(l = 0; l < 3; l++) {
566: JLocal[localRows[lowerRow[k]]*4 + lowerRow[l]] += 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]);
567: }
568: }
569: /* Upper Element */
570: for(k = 0; k < 3; k++) {
571: if (!hasUpper[k]) continue;
572: for(l = 0; l < 3; l++) {
573: JLocal[localRows[upperRow[k]]*4 + upperRow[l]] += 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]);
574: }
575: }
577: nonlinearJacobian(-1.0*sc, uLocal, JLocal);
578: printf("Element matrix for (%d, %d)\n", i, j);
579: printf(" col ");
580: for(l = 0; l < 4; l++) {
581: printf("(%d, %d) ", cols[l].i, cols[l].j);
582: }
583: printf("\n");
584: for(k = 0; k < numRows; k++) {
585: printf("row (%d, %d): ", rows[k].i, rows[k].j);
586: for(l = 0; l < 4; l++) {
587: printf("%8.6g ", JLocal[k*4 + l]);
588: }
589: printf("\n");
590: }
591: MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
592: }
593: }
595: /*
596: Assemble matrix, using the 2-step process:
597: MatAssemblyBegin(), MatAssemblyEnd().
598: */
599: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
600: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
601: /*
602: Tell the matrix we will never add a new nonzero location to the
603: matrix. If we do, it will generate an error.
604: */
605: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
606: return(0);
607: }
611: /*
612: L_2Error - Integrate the L_2 error of our solution over each face
613: */
614: PetscErrorCode L_2Error(DM da, Vec fVec, double *error, AppCtx *user)
615: {
616: DMDALocalInfo info;
617: Vec fLocalVec;
618: PetscScalar **f;
619: PetscScalar u, uExact, uLocal[4];
620: PetscScalar hx, hy, hxhy, x, y, phi[3];
621: PetscInt i, j, q;
625: DMDAGetLocalInfo(da, &info);
626: DMGetLocalVector(da, &fLocalVec);
627: DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
628: DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
629: DMDAVecGetArray(da, fLocalVec, &f);
631: *error = 0.0;
632: hx = 1.0/(PetscReal)(info.mx-1);
633: hy = 1.0/(PetscReal)(info.my-1);
634: hxhy = hx*hy;
635: for(j = info.ys; j < info.ys+info.ym-1; j++) {
636: for(i = info.xs; i < info.xs+info.xm-1; i++) {
637: uLocal[0] = f[j][i];
638: uLocal[1] = f[j][i+1];
639: uLocal[2] = f[j+1][i+1];
640: uLocal[3] = f[j+1][i];
641: /* Lower element */
642: for(q = 0; q < 4; q++) {
643: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
644: phi[1] = quadPoints[q*2];
645: phi[2] = quadPoints[q*2+1];
646: u = uLocal[0]*phi[0] + uLocal[1]*phi[1] + uLocal[3]*phi[2];
647: x = (quadPoints[q*2] + i)*hx;
648: y = (quadPoints[q*2+1] + j)*hy;
649: ExactSolution(x, y, &uExact);
650: *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
651: }
652: /* Upper element */
653: /*
654: The affine map from the lower to the upper is
656: / x_U \ = / -1 0 \ / x_L \ + / hx \
657: \ y_U / \ 0 -1 / \ y_L / \ hy /
658: */
659: for(q = 0; q < 4; q++) {
660: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
661: phi[1] = quadPoints[q*2];
662: phi[2] = quadPoints[q*2+1];
663: u = uLocal[2]*phi[0] + uLocal[3]*phi[1] + uLocal[1]*phi[2];
664: x = (1.0 - quadPoints[q*2] + i)*hx;
665: y = (1.0 - quadPoints[q*2+1] + j)*hy;
666: ExactSolution(x, y, &uExact);
667: *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
668: }
669: }
670: }
672: DMDAVecRestoreArray(da, fLocalVec, &f);
673: /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
674: /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
675: DMRestoreLocalVector(da, &fLocalVec);
676: return(0);
677: }