Actual source code: ex8.c
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 <petscdmmg.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: static PetscScalar Kref[36] = { 0.5, 0.5, -0.5, 0, 0, -0.5,
56: 0.5, 0.5, -0.5, 0, 0, -0.5,
57: -0.5, -0.5, 0.5, 0, 0, 0.5,
58: 0, 0, 0, 0, 0, 0,
59: 0, 0, 0, 0, 0, 0,
60: -0.5, -0.5, 0.5, 0, 0, 0.5};
62: /* These are */
63: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
64: 0.07503111, 0.64494897,
65: 0.66639025, 0.15505103,
66: 0.28001992, 0.64494897};
67: static PetscScalar quadWeights[4] = {0.15902069, 0.09097931, 0.15902069, 0.09097931};
69: /*
70: User-defined routines
71: */
80: int main(int argc,char **argv)
81: {
82: DMMG *dmmg; /* hierarchy manager */
83: DM da;
84: SNES snes; /* nonlinear solver */
85: AppCtx *user; /* user-defined work context */
86: PetscBag bag;
87: PetscInt its; /* iterations for convergence */
88: SNESConvergedReason reason;
89: PetscBool drawContours; /* flag for drawing contours */
90: PetscErrorCode ierr;
91: PetscReal lambda_max = 6.81, lambda_min = 0.0, error;
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Initialize program
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: PetscInitialize(&argc,&argv,(char *)0,help);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Initialize problem parameters
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
102: PetscBagGetData(bag, (void **) &user);
103: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
104: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
105: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
106: PetscBagSetFromOptions(bag);
107: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
108: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
109: if (user->lambda > lambda_max || user->lambda < lambda_min) {
110: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
111: }
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create multilevel DM data structure (DMMG) to manage hierarchical solvers
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: DMMGCreate(PETSC_COMM_WORLD,1,user,&dmmg);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Create distributed array (DA) to manage parallel grid and vectors
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
122: 1,1,PETSC_NULL,PETSC_NULL,&da);
123: DMDASetFieldName(da, 0, "ooblek");
124: DMMGSetDM(dmmg, (DM) da);
125: DMDestroy(&da);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Set the discretization functions
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);
131: DMMGSetFromOptions(dmmg);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Evaluate initial guess
135: Note: The user should initialize the vector, x, with the initial guess
136: for the nonlinear solver prior to calling SNESSolve(). In particular,
137: to employ an initial guess of zero, the user should explicitly set
138: this vector to zero by calling VecSet().
139: */
140: DMMGSetInitialGuess(dmmg, FormInitialGuess);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve nonlinear system
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: DMMGSolve(dmmg);
146: snes = DMMGGetSNES(dmmg);
147: SNESGetIterationNumber(snes,&its);
148: SNESGetConvergedReason(snes, &reason);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
152: L_2Error(DMMGGetDM(dmmg), DMMGGetx(dmmg), &error, user);
153: PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
154: PrintVector(dmmg[0], DMMGGetx(dmmg));
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Free work space. All PETSc objects should be destroyed when they
158: are no longer needed.
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: DMMGDestroy(dmmg);
161: PetscBagDestroy(&bag);
162: PetscFinalize();
163: return(0);
164: }
168: PetscErrorCode PrintVector(DMMG dmmg, Vec U)
169: {
170: DM da = dmmg->dm;
171: PetscScalar **u;
172: PetscInt i,j,xs,ys,xm,ym;
176: DMDAVecGetArray(da,U,&u);
177: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
178: for(j = ys+ym-1; j >= ys; j--) {
179: for(i = xs; i < xs+xm; i++) {
180: printf("u[%d][%d] = %G ", j, i, u[j][i]);
181: }
182: printf("\n");
183: }
184: DMDAVecRestoreArray(da,U,&u);
185: return(0);
186: }
190: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
191: {
193: *u = x*x;
194: return(0);
195: }
199: /*
200: FormInitialGuess - Forms initial approximation.
202: Input Parameters:
203: dmmg - The DMMG context
204: X - vector
206: Output Parameter:
207: X - vector
208: */
209: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
210: {
211: AppCtx *user = (AppCtx *) dmmg->user;
212: DM da = dmmg->dm;
213: PetscInt i,j,Mx,My,xs,ys,xm,ym;
215: PetscReal lambda,temp1,temp,hx,hy;
216: PetscScalar **x;
219: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
220: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
222: lambda = user->lambda;
223: hx = 1.0/(PetscReal)(Mx-1);
224: hy = 1.0/(PetscReal)(My-1);
225: if (lambda == 0.0) {
226: temp1 = 0.0;
227: } else {
228: temp1 = lambda/(lambda + 1.0);
229: }
231: /*
232: Get a pointer to vector data.
233: - For default PETSc vectors, VecGetArray() returns a pointer to
234: the data array. Otherwise, the routine is implementation dependent.
235: - You MUST call VecRestoreArray() when you no longer need access to
236: the array.
237: */
238: DMDAVecGetArray(da,X,&x);
240: /*
241: Get local grid boundaries (for 2-dimensional DMDA):
242: xs, ys - starting grid indices (no ghost points)
243: xm, ym - widths of local grid (no ghost points)
245: */
246: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
248: /*
249: Compute initial guess over the locally owned part of the grid
250: */
251: for (j=ys; j<ys+ym; j++) {
252: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
253: for (i=xs; i<xs+xm; i++) {
255: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
256: /* boundary conditions are all zero Dirichlet */
257: x[j][i] = 0.0;
258: } else {
259: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
260: }
261: }
262: }
264: DMDAVecRestoreArray(da,X,&x);
265: PrintVector(dmmg, X);
266: return(0);
267: }
271: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
272: {
273: PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
274: PetscScalar phi[3] = {0.0, 0.0, 0.0};
275: PetscReal xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
276: PetscScalar res;
277: PetscInt q, k;
280: for(q = 0; q < 4; q++) {
281: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
282: phi[1] = quadPoints[q*2];
283: phi[2] = quadPoints[q*2+1];
284: /* These are currently wrong */
285: x = xI + quadPoints[q*2]*hx;
286: y = yI + quadPoints[q*2+1]*hy;
287: res = quadWeights[q]*(2.0);
288: for(k = 0; k < 3; k++) {
289: rLocal[k] += phi[k]*res;
290: }
291: }
292: for(k = 0; k < 3; k++) {
293: printf(" constLocal[%d] = %g\n", k, lambda*hxhy*rLocal[k]);
294: r[k] += lambda*hxhy*rLocal[k];
295: }
296: return(0);
297: }
299: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[]) {
300: PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
301: PetscScalar phi[3] = {0.0, 0.0, 0.0};
302: PetscScalar res;
303: PetscInt q;
306: for(q = 0; q < 4; q++) {
307: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
308: phi[1] = quadPoints[q*2];
309: phi[2] = quadPoints[q*2+1];
310: res = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
311: rLocal[0] += phi[0]*res;
312: rLocal[1] += phi[1]*res;
313: rLocal[2] += phi[2]*res;
314: }
315: r[0] += lambda*rLocal[0];
316: r[1] += lambda*rLocal[1];
317: r[2] += lambda*rLocal[2];
318: return(0);
319: }
323: /*
324: FormFunctionLocal - Evaluates nonlinear function, F(x).
326: Process adiC(36): FormFunctionLocal
328: */
329: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
330: {
331: PetscScalar uLocal[3];
332: PetscScalar rLocal[3];
333: PetscScalar G[4];
334: PetscScalar uExact;
335: PetscReal alpha,lambda,hx,hy,hxhy,sc,detJInv;
336: PetscInt i,j,k,l;
341: /* Naive Jacobian calculation:
343: J = / 1/hx 0 \ J^{-1} = / hx 0 \ 1/|J| = hx*hy = |J^{-1}|
344: \ 0 1/hy / \ 0 hy /
345: */
346: alpha = user->alpha;
347: lambda = user->lambda;
348: hx = 1.0/(PetscReal)(info->mx-1);
349: hy = 1.0/(PetscReal)(info->my-1);
350: sc = hx*hy*lambda;
351: hxhy = hx*hy;
352: detJInv = hxhy;
353: G[0] = (1.0/(hx*hx)) * detJInv;
354: G[1] = 0.0;
355: G[2] = G[1];
356: G[3] = (1.0/(hy*hy)) * detJInv;
357: for(k = 0; k < 4; k++) {
358: printf("G[%d] = %g\n", k, G[k]);
359: }
361: /* Zero the vector */
362: PetscMemzero((void *) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(PetscScalar));
363: /* Compute function over the locally owned part of the grid. For each
364: vertex (i,j), we consider the element below:
366: 2 (1) (0)
367: i,j+1 --- i+1,j+1
368: | \ |
369: | \ |
370: | \ |
371: | \ |
372: | \ |
373: i,j --- i+1,j
374: 0 1 (2)
376: and therefore we do not loop over the last vertex in each dimension.
377: */
378: for(j = info->ys; j < info->ys+info->ym-1; j++) {
379: for(i = info->xs; i < info->xs+info->xm-1; i++) {
380: /* Lower element */
381: uLocal[0] = x[j][i];
382: uLocal[1] = x[j][i+1];
383: uLocal[2] = x[j+1][i];
384: printf("Solution ElementVector for (%d, %d)\n", i, j);
385: for(k = 0; k < 3; k++) {
386: printf(" uLocal[%d] = %g\n", k, uLocal[k]);
387: }
388: for(k = 0; k < 3; k++) {
389: rLocal[k] = 0.0;
390: for(l = 0; l < 3; l++) {
391: 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];
392: }
393: rLocal[k] *= alpha;
394: }
395: printf("Laplacian 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: constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
400: printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
401: for(k = 0; k < 3; k++) {
402: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
403: }
404: nonlinearResidual(0.0*sc, uLocal, rLocal);
405: printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
406: for(k = 0; k < 3; k++) {
407: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
408: }
409: f[j][i] += rLocal[0];
410: f[j][i+1] += rLocal[1];
411: f[j+1][i] += rLocal[2];
412: /* Upper element */
413: uLocal[0] = x[j+1][i+1];
414: uLocal[1] = x[j+1][i];
415: uLocal[2] = x[j][i+1];
416: printf("Solution ElementVector for (%d, %d)\n", i, j);
417: for(k = 0; k < 3; k++) {
418: printf(" uLocal[%d] = %g\n", k, uLocal[k]);
419: }
420: for(k = 0; k < 3; k++) {
421: rLocal[k] = 0.0;
422: for(l = 0; l < 3; l++) {
423: 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];
424: }
425: rLocal[k] *= alpha;
426: }
427: printf("Laplacian ElementVector for (%d, %d)\n", i, j);
428: for(k = 0; k < 3; k++) {
429: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
430: }
431: constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
432: printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
433: for(k = 0; k < 3; k++) {
434: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
435: }
436: nonlinearResidual(0.0*sc, uLocal, rLocal);
437: printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
438: for(k = 0; k < 3; k++) {
439: printf(" rLocal[%d] = %g\n", k, rLocal[k]);
440: }
441: f[j+1][i+1] += rLocal[0];
442: f[j+1][i] += rLocal[1];
443: f[j][i+1] += rLocal[2];
444: /* Boundary conditions */
445: if (i == 0 || j == 0) {
446: ExactSolution(i*hx, j*hy, &uExact);
447: f[j][i] = x[j][i] - uExact;
448: }
449: if ((i == info->mx-2) || (j == 0)) {
450: ExactSolution((i+1)*hx, j*hy, &uExact);
451: f[j][i+1] = x[j][i+1] - uExact;
452: }
453: if ((i == info->mx-2) || (j == info->my-2)) {
454: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
455: f[j+1][i+1] = x[j+1][i+1] - uExact;
456: }
457: if ((i == 0) || (j == info->my-2)) {
458: ExactSolution(i*hx, (j+1)*hy, &uExact);
459: f[j+1][i] = x[j+1][i] - uExact;
460: }
461: }
462: }
464: for(j = info->ys+info->ym-1; j >= info->ys; j--) {
465: for(i = info->xs; i < info->xs+info->xm; i++) {
466: printf("f[%d][%d] = %g ", j, i, f[j][i]);
467: }
468: printf("\n");
469: }
470: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
471: return(0);
472: }
474: PetscErrorCode nonlinearJacobian(PetscReal lambda, PetscScalar u[], PetscScalar J[]) {
476: return(0);
477: }
481: /*
482: FormJacobianLocal - Evaluates Jacobian matrix.
483: */
484: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
485: {
486: PetscScalar JLocal[16], uLocal[4];
487: MatStencil rows[4], cols[4], ident;
488: PetscInt lowerRow[3] = {0, 1, 3};
489: PetscInt upperRow[3] = {2, 3, 1};
490: PetscInt hasLower[3], hasUpper[3], localRows[4];
491: PetscScalar alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
492: PetscInt i,j,k,l,numRows;
496: alpha = user->alpha;
497: lambda = user->lambda;
498: hx = 1.0/(PetscReal)(info->mx-1);
499: hy = 1.0/(PetscReal)(info->my-1);
500: sc = hx*hy*lambda;
501: hxhy = hx*hy;
502: detJInv = hxhy;
503: G[0] = (1.0/(hx*hx)) * detJInv;
504: G[1] = 0.0;
505: G[2] = G[1];
506: G[3] = (1.0/(hy*hy)) * detJInv;
507: for(k = 0; k < 4; k++) {
508: printf("G[%d] = %g\n", k, G[k]);
509: }
511: MatZeroEntries(jac);
512: /*
513: Compute entries for the locally owned part of the Jacobian.
514: - Currently, all PETSc parallel matrix formats are partitioned by
515: contiguous chunks of rows across the processors.
516: - Each processor needs to insert only elements that it owns
517: locally (but any non-local elements will be sent to the
518: appropriate processor during matrix assembly).
519: - Here, we set all entries for a particular row at once.
520: - We can set matrix entries either using either
521: MatSetValuesLocal() or MatSetValues(), as discussed above.
522: */
523: for (j=info->ys; j<info->ys+info->ym-1; j++) {
524: for (i=info->xs; i<info->xs+info->xm-1; i++) {
525: PetscMemzero(JLocal, 16 * sizeof(PetscScalar));
526: numRows = 0;
527: /* Lower element */
528: uLocal[0] = x[j][i];
529: uLocal[1] = x[j][i+1];
530: uLocal[2] = x[j+1][i+1];
531: uLocal[3] = x[j+1][i];
532: /* i,j */
533: if (i == 0 || j == 0) {
534: hasLower[0] = 0;
535: ident.i = i; ident.j = j;
536: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
537: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
538: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
539: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
540: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
541: } else {
542: hasLower[0] = 1;
543: localRows[0] = numRows;
544: rows[numRows].i = i; rows[numRows].j = j;
545: numRows++;
546: }
547: cols[0].i = i; cols[0].j = j;
548: /* i+1,j */
549: if ((i == info->mx-2) || (j == 0)) {
550: hasLower[1] = 0;
551: hasUpper[2] = 0;
552: ident.i = i+1; ident.j = j;
553: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
554: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
555: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
556: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
557: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
558: } else {
559: localRows[1] = numRows;
560: hasLower[1] = 1;
561: hasUpper[2] = 1;
562: localRows[1] = numRows;
563: rows[numRows].i = i+1; rows[numRows].j = j;
564: numRows++;
565: }
566: cols[1].i = i+1; cols[1].j = j;
567: /* i+1,j+1 */
568: if ((i == info->mx-2) || (j == info->my-2)) {
569: hasUpper[0] = 0;
570: ident.i = i+1; ident.j = j+1;
571: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
572: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
573: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
574: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
575: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
576: } else {
577: hasUpper[0] = 1;
578: localRows[2] = numRows;
579: rows[numRows].i = i+1; rows[numRows].j = j+1;
580: numRows++;
581: }
582: cols[2].i = i+1; cols[2].j = j+1;
583: /* i,j+1 */
584: if ((i == 0) || (j == info->my-2)) {
585: hasLower[2] = 0;
586: hasUpper[1] = 0;
587: ident.i = i; ident.j = j+1;
588: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
589: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
590: MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
591: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
592: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
593: } else {
594: hasLower[2] = 1;
595: hasUpper[1] = 1;
596: localRows[3] = numRows;
597: rows[numRows].i = i; rows[numRows].j = j+1;
598: numRows++;
599: }
600: cols[3].i = i; cols[3].j = j+1;
602: /* Lower Element */
603: for(k = 0; k < 3; k++) {
604: if (!hasLower[k]) continue;
605: for(l = 0; l < 3; l++) {
606: 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]);
607: }
608: }
609: /* Upper Element */
610: for(k = 0; k < 3; k++) {
611: if (!hasUpper[k]) continue;
612: for(l = 0; l < 3; l++) {
613: 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]);
614: }
615: }
617: nonlinearJacobian(-1.0*sc, uLocal, JLocal);
618: printf("Element matrix for (%d, %d)\n", i, j);
619: printf(" col ");
620: for(l = 0; l < 4; l++) {
621: printf("(%d, %d) ", cols[l].i, cols[l].j);
622: }
623: printf("\n");
624: for(k = 0; k < numRows; k++) {
625: printf("row (%d, %d): ", rows[k].i, rows[k].j);
626: for(l = 0; l < 4; l++) {
627: printf("%8.6g ", JLocal[k*4 + l]);
628: }
629: printf("\n");
630: }
631: MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
632: }
633: }
635: /*
636: Assemble matrix, using the 2-step process:
637: MatAssemblyBegin(), MatAssemblyEnd().
638: */
639: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
640: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
641: /*
642: Tell the matrix we will never add a new nonzero location to the
643: matrix. If we do, it will generate an error.
644: */
645: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
646: return(0);
647: }
651: /*
652: L_2Error - Integrate the L_2 error of our solution over each face
653: */
654: PetscErrorCode L_2Error(DM da, Vec fVec, double *error, AppCtx *user)
655: {
656: DMDALocalInfo info;
657: Vec fLocalVec;
658: PetscScalar **f;
659: PetscScalar u, uExact, uLocal[4];
660: PetscScalar hx, hy, hxhy, x, y, phi[3];
661: PetscInt i, j, q;
665: DMDAGetLocalInfo(da, &info);
666: DMGetLocalVector(da, &fLocalVec);
667: DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
668: DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
669: DMDAVecGetArray(da, fLocalVec, &f);
671: *error = 0.0;
672: hx = 1.0/(PetscReal)(info.mx-1);
673: hy = 1.0/(PetscReal)(info.my-1);
674: hxhy = hx*hy;
675: for(j = info.ys; j < info.ys+info.ym-1; j++) {
676: for(i = info.xs; i < info.xs+info.xm-1; i++) {
677: uLocal[0] = f[j][i];
678: uLocal[1] = f[j][i+1];
679: uLocal[2] = f[j+1][i+1];
680: uLocal[3] = f[j+1][i];
681: /* Lower element */
682: for(q = 0; q < 4; q++) {
683: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
684: phi[1] = quadPoints[q*2];
685: phi[2] = quadPoints[q*2+1];
686: u = uLocal[0]*phi[0] + uLocal[1]*phi[1] + uLocal[3]*phi[2];
687: x = (quadPoints[q*2] + i)*hx;
688: y = (quadPoints[q*2+1] + j)*hy;
689: ExactSolution(x, y, &uExact);
690: *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
691: }
692: /* Upper element */
693: /*
694: The affine map from the lower to the upper is
696: / x_U \ = / -1 0 \ / x_L \ + / hx \
697: \ y_U / \ 0 -1 / \ y_L / \ hy /
698: */
699: for(q = 0; q < 4; q++) {
700: phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
701: phi[1] = quadPoints[q*2];
702: phi[2] = quadPoints[q*2+1];
703: u = uLocal[2]*phi[0] + uLocal[3]*phi[1] + uLocal[1]*phi[2];
704: x = (1.0 - quadPoints[q*2] + i)*hx;
705: y = (1.0 - quadPoints[q*2+1] + j)*hy;
706: ExactSolution(x, y, &uExact);
707: *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
708: }
709: }
710: }
712: DMDAVecRestoreArray(da, fLocalVec, &f);
713: /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
714: /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
715: DMRestoreLocalVector(da, &fLocalVec);
716: return(0);
717: }