Actual source code: ex4.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 Lane-Emden 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 Lane-Emden example
9: Concepts: DMDA^using distributed arrays;
10: Processors: n
11: T*/
13: /* ------------------------------------------------------------------------
15: The Lane-Emden equation is given by the partial differential equation
16:
17: -alpha*Laplacian u - lambda*u^3 = 0, 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 bilinear finite element approximation is used to discretize the boundary
24: value problem to obtain a nonlinear system of equations.
26: ------------------------------------------------------------------------- */
28: /*
29: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
30: Include "petscsnes.h" so that we can use SNES solvers. Note that this
31: file automatically includes:
32: petscsys.h - base PETSc routines petscvec.h - vectors
33: petscmat.h - matrices
34: petscis.h - index sets petscksp.h - Krylov subspace methods
35: petscviewer.h - viewers petscpc.h - preconditioners
36: petscksp.h - linear solvers
37: */
38: #include <petscsys.h>
39: #include <petscbag.h>
40: #include <petscdmda.h>
41: #include <petscsnes.h>
43: /*
44: User-defined application context - contains data needed by the
45: application-provided call-back routines, FormJacobianLocal() and
46: FormFunctionLocal().
47: */
48: typedef struct {
49: PetscReal alpha; /* parameter controlling linearity */
50: PetscReal lambda; /* parameter controlling nonlinearity */
51: } AppCtx;
53: static PetscScalar Kref[16] = { 0.666667, -0.166667, -0.333333, -0.166667,
54: -0.166667, 0.666667, -0.166667, -0.333333,
55: -0.333333, -0.166667, 0.666667, -0.166667,
56: -0.166667, -0.333333, -0.166667, 0.666667};
58: /* These are */
59: static PetscScalar quadPoints[8] = {0.211325, 0.211325,
60: 0.788675, 0.211325,
61: 0.788675, 0.788675,
62: 0.211325, 0.788675};
63: static PetscScalar quadWeights[4] = {0.25, 0.25, 0.25, 0.25};
65: /*
66: User-defined routines
67: */
68: extern PetscErrorCode FormInitialGuess(DM,Vec);
69: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
70: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*);
71: extern PetscErrorCode PrintVector(DM, Vec);
75: int main(int argc,char **argv)
76: {
77: DM da;
78: SNES snes; /* nonlinear solver */
79: AppCtx *user; /* user-defined work context */
80: PetscBag bag;
81: PetscInt its; /* iterations for convergence */
82: SNESConvergedReason reason;
83: PetscErrorCode ierr;
84: PetscReal lambda_max = 6.81, lambda_min = 0.0;
85: Vec x;
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Initialize program
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: PetscInitialize(&argc,&argv,(char *)0,help);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Initialize problem parameters
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
96: PetscBagGetData(bag, (void **) &user);
97: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
98: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
99: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
100: PetscBagSetFromOptions(bag);
101: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
102: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
103: if (user->lambda > lambda_max || user->lambda < lambda_min) {
104: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
105: }
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create SNES to manage hierarchical solvers
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: SNESCreate(PETSC_COMM_WORLD,&snes);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create distributed array (DMDA) 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,1,1,PETSC_NULL,PETSC_NULL,&da);
116: DMDASetFieldName(da, 0, "ooblek");
117: DMSetApplicationContext(da,user);
118: SNESSetDM(snes, (DM) da);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Set the discretization functions
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
124: DMDASetLocalJacobian(da,(DMDALocalFunction1)FormJacobianLocal);
125: SNESSetFromOptions(snes);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Evaluate initial guess
129: Note: The user should initialize the vector, x, with the initial guess
130: for the nonlinear solver prior to calling SNESSolve(). In particular,
131: to employ an initial guess of zero, the user should explicitly set
132: this vector to zero by calling VecSet().
133: */
134: DMSetInitialGuess(da, FormInitialGuess);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Solve nonlinear system
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: SNESSolve(snes,PETSC_NULL,PETSC_NULL);
140: SNESGetIterationNumber(snes,&its);
141: SNESGetConvergedReason(snes, &reason);
142: DMDestroy(&da);
143: SNESGetDM(snes,&da);
144: SNESGetSolution(snes,&x);
145: PrintVector(da, x);
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Free work space. All PETSc objects should be destroyed when they
149: are no longer needed.
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: SNESDestroy(&snes);
152: PetscBagDestroy(&bag);
153: PetscFinalize();
154: return(0);
155: }
159: PetscErrorCode PrintVector(DM da, Vec U)
160: {
161: PetscScalar **u;
162: PetscInt i,j,xs,ys,xm,ym;
166: DMDAVecGetArray(da,U,&u);
167: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
168: for(j = ys+ym-1; j >= ys; j--) {
169: for(i = xs; i < xs+xm; i++) {
170: PetscPrintf(PETSC_COMM_SELF,"u[%d][%d] = %G ", j, i, PetscRealPart(u[j][i]));
171: }
172: PetscPrintf(PETSC_COMM_SELF,"\n");
173: }
174: DMDAVecRestoreArray(da,U,&u);
175: return(0);
176: }
180: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
181: {
183: *u = x*x;
184: return(0);
185: }
189: /*
190: FormInitialGuess - Forms initial approximation.
192: Input Parameters:
193: X - vector
195: Output Parameter:
196: X - vector
197: */
198: PetscErrorCode FormInitialGuess(DM da,Vec X)
199: {
200: AppCtx *user;
201: PetscInt i,j,Mx,My,xs,ys,xm,ym;
203: PetscReal lambda,temp1,temp,hx,hy;
204: PetscScalar **x;
207: DMGetApplicationContext(da,&user);
208: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
209: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
211: lambda = user->lambda;
212: hx = 1.0/(PetscReal)(Mx-1);
213: hy = 1.0/(PetscReal)(My-1);
214: if (lambda == 0.0) {
215: temp1 = 1.0;
216: } else {
217: temp1 = lambda/(lambda + 1.0);
218: }
220: /*
221: Get a pointer to vector data.
222: - For default PETSc vectors, VecGetArray() returns a pointer to
223: the data array. Otherwise, the routine is implementation dependent.
224: - You MUST call VecRestoreArray() when you no longer need access to
225: the array.
226: */
227: DMDAVecGetArray(da,X,&x);
229: /*
230: Get local grid boundaries (for 2-dimensional DMDA):
231: xs, ys - starting grid indices (no ghost points)
232: xm, ym - widths of local grid (no ghost points)
234: */
235: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
237: /*
238: Compute initial guess over the locally owned part of the grid
239: */
240: for (j=ys; j<ys+ym; j++) {
241: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
242: for (i=xs; i<xs+xm; i++) {
244: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
245: /* boundary conditions are all zero Dirichlet */
246: x[j][i] = 0.0;
247: } else {
248: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
249: }
250: }
251: }
253: /*
254: Restore vector
255: */
256: DMDAVecRestoreArray(da,X,&x);
258: return(0);
259: }
263: PetscErrorCode constantResidual(PetscReal lambda, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
264: {
265: PetscScalar rLocal[4] = {0.0, 0.0, 0.0};
266: PetscScalar phi[4] = {0.0, 0.0, 0.0, 0.0};
267: /* PetscReal xI = i*hx, yI = j*hy, x, y; */
268: PetscReal hxhy = hx*hy;
269: PetscScalar res;
270: PetscInt q, k;
273: for(q = 0; q < 4; q++) {
274: phi[0] = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
275: phi[1] = quadPoints[q*2] *(1.0 - quadPoints[q*2+1]);
276: phi[2] = quadPoints[q*2] * quadPoints[q*2+1];
277: phi[3] = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
278: /*
279: x = xI + quadPoints[q*2]*hx;
280: y = yI + quadPoints[q*2+1]*hy;
281: */
282: res = quadWeights[q]*(2.0);
283: for(k = 0; k < 4; k++) {
284: rLocal[k] += phi[k]*res;
285: }
286: }
287: for(k = 0; k < 4; k++) {
288: r[k] += lambda*hxhy*rLocal[k];
289: }
290: return(0);
291: }
295: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[]) {
297: r[0] += lambda*(48.0*u[0]*u[0]*u[0] + 12.0*u[1]*u[1]*u[1] + 9.0*u[0]*u[0]*(4.0*u[1] + u[2] + 4.0*u[3]) + u[1]*u[1]*(9.0*u[2] + 6.0*u[3]) + u[1]*(6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3])
298: + 3.0*(u[2]*u[2]*u[2] + 2.0*u[2]*u[2]*u[3] + 3.0*u[2]*u[3]*u[3] + 4.0*u[3]*u[3]*u[3])
299: + 2.0*u[0]*(12.0*u[1]*u[1] + u[1]*(6.0*u[2] + 9.0*u[3]) + 2.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3])))/1200.0;
300: r[1] += lambda*(12.0*u[0]*u[0]*u[0] + 48.0*u[1]*u[1]*u[1] + 9.0*u[1]*u[1]*(4.0*u[2] + u[3]) + 3.0*u[0]*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3])
301: + 4.0*u[1]*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]) + 3.0*(4.0*u[2]*u[2]*u[2] + 3.0*u[2]*u[2]*u[3] + 2.0*u[2]*u[3]*u[3] + u[3]*u[3]*u[3])
302: + 2.0*u[0]*((18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3])))/1200.0;
303: r[2] += lambda*(3.0*u[0]*u[0]*u[0] + u[0]*u[0]*(6.0*u[1] + 4.0*u[2] + 6.0*u[3]) + u[0]*(9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]))
304: + 6.0*(2.0*u[1]*u[1]*u[1] + u[1]*u[1]*(4.0*u[2] + u[3]) + u[1]*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]) + 2.0*(4.0*u[2]*u[2]*u[2] + 3.0*u[2]*u[2]*u[3] + 2.0*u[2]*u[3]*u[3] + u[3]*u[3]*u[3])))/1200.0;
305: r[3] += lambda*(12.0*u[0]*u[0]*u[0] + 3.0*u[1]*u[1]*u[1] + u[1]*u[1]*(6.0*u[2] + 4.0*u[3]) + 3.0*u[0]*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3])
306: + 3.0*u[1]*(3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3]) + 12.0*(u[2]*u[2]*u[2] + 2.0*u[2]*u[2]*u[3] + 3.0*u[2]*u[3]*u[3] + 4.0*u[3]*u[3]*u[3])
307: + 2.0*u[0]*(3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3])))/1200.0;
308: return(0);
309: }
313: PetscErrorCode nonlinearResidualBratu(PetscReal lambda, PetscScalar u[], PetscScalar r[]) {
314: PetscScalar rLocal[4] = {0.0, 0.0, 0.0, 0.0};
315: PetscScalar phi[4] = {0.0, 0.0, 0.0, 0.0};
316: PetscScalar res;
317: PetscInt q;
320: for(q = 0; q < 4; q++) {
321: phi[0] = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
322: phi[1] = quadPoints[q*2] *(1.0 - quadPoints[q*2+1]);
323: phi[2] = quadPoints[q*2] * quadPoints[q*2+1];
324: phi[3] = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
325: res = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
326: rLocal[0] += phi[0]*res;
327: rLocal[1] += phi[1]*res;
328: rLocal[2] += phi[2]*res;
329: rLocal[3] += phi[3]*res;
330: }
331: r[0] += lambda*rLocal[0];
332: r[1] += lambda*rLocal[1];
333: r[2] += lambda*rLocal[2];
334: r[3] += lambda*rLocal[3];
335: return(0);
336: }
340: PetscErrorCode nonlinearJacobian(PetscScalar lambda, PetscScalar u[], PetscScalar J[]) {
342: J[0] = lambda*(72.0*u[0]*u[0] + 12.0*u[1]*u[1] + 9.0*u[0]*(4.0*u[1] + u[2] + 4.0*u[3]) + u[1]*(6.0*u[2] + 9.0*u[3]) + 2.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;
343: J[1] = lambda*(18.0*u[0]*u[0] + 18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3] + 3.0*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3]))/600.0;
344: J[2] = lambda*( 9.0*u[0]*u[0] + 9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
345: J[3] = lambda*(18.0*u[0]*u[0] + 3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;
347: J[4] = lambda*(18.0*u[0]*u[0] + 18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3] + 3.0*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3]))/600.0;
348: J[5] = lambda*(12.0*u[0]*u[0] + 72.0*u[1]*u[1] + 9.0*u[1]*(4.0*u[2] + u[3]) + u[0]*(36.0*u[1] + 9.0*u[2] + 6.0*u[3]) + 2.0*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]))/600.0;
349: J[6] = lambda*( 3.0*u[0]*u[0] + u[0]*(9.0*u[1] + 6.0*u[2] + 4.0*u[3]) + 3.0*(6.0*u[1]*u[1] + 6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3] + 2.0*u[1]*(4.0*u[2] + u[3])))/600.0;
350: J[7] = lambda*( 9.0*u[0]*u[0] + 9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
352: J[8] = lambda*( 9.0*u[0]*u[0] + 9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
353: J[9] = lambda*( 3.0*u[0]*u[0] + u[0]*(9.0*u[1] + 6.0*u[2] + 4.0*u[3]) + 3.0*(6.0*u[1]*u[1] + 6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3] + 2.0*u[1]*(4.0*u[2] + u[3])))/600.0;
354: J[10] = lambda*( 2.0*u[0]*u[0] + u[0]*(6.0*u[1] + 9.0*u[2] + 6.0*u[3]) + 3.0*(4.0*u[1]*u[1] + 3.0*u[1]*(4.0*u[2] + u[3]) + 4.0*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3])))/600.0;
355: J[11] = lambda*( 3.0*u[0]*u[0] + u[0]*(4.0*u[1] + 6.0*u[2] + 9.0*u[3]) + 3.0*(u[1]*u[1] + 6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3] + u[1]*(3.0*u[2] + 2.0*u[3])))/600.0;
357: J[12] = lambda*(18.0*u[0]*u[0] + 3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;
358: J[13] = lambda*( 9.0*u[0]*u[0] + 9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
359: J[14] = lambda*( 3.0*u[0]*u[0] + u[0]*(4.0*u[1] + 6.0*u[2] + 9.0*u[3]) + 3.0*(u[1]*u[1] + 6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3] + u[1]*(3.0*u[2] + 2.0*u[3])))/600.0;
360: J[15] = lambda*(12.0*u[0]*u[0] + 2.0*u[1]*u[1] + u[1]*(6.0*u[2] + 9.0*u[3]) + 12.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]) + u[0]*(6.0*u[1] + 9.0*(u[2] + 4.0*u[3])))/600.0;
361: return(0);
362: }
366: /*
367: FormFunctionLocal - Evaluates nonlinear function, F(x).
369: Process adiC(36): FormFunctionLocal
371: */
372: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
373: {
374: PetscScalar uLocal[4];
375: PetscScalar rLocal[4];
376: PetscScalar uExact;
377: PetscReal alpha,lambda,hx,hy,hxhy,sc;
378: PetscInt i,j,k,l;
382: alpha = user->alpha;
383: lambda = user->lambda;
384: hx = 1.0/(PetscReal)(info->mx-1);
385: hy = 1.0/(PetscReal)(info->my-1);
386: sc = hx*hy*lambda;
387: hxhy = hx*hy;
389: /* Zero the vector */
390: PetscMemzero((void *) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(PetscScalar));
391: /* Compute function over the locally owned part of the grid. For each
392: vertex (i,j), we consider the element below:
394: 3 2
395: i,j+1 --- i+1,j+1
396: | |
397: | |
398: i,j --- i+1,j
399: 0 1
401: and therefore we do not loop over the last vertex in each dimension.
402: */
403: for(j = info->ys; j < info->ys+info->ym-1; j++) {
404: for(i = info->xs; i < info->xs+info->xm-1; i++) {
405: uLocal[0] = x[j][i];
406: uLocal[1] = x[j][i+1];
407: uLocal[2] = x[j+1][i+1];
408: uLocal[3] = x[j+1][i];
409: for(k = 0; k < 4; k++) {
410: rLocal[k] = 0.0;
411: for(l = 0; l < 4; l++) {
412: rLocal[k] += Kref[k*4 + l]*uLocal[l];
413: }
414: rLocal[k] *= hxhy*alpha;
415: }
416: constantResidual(1.0, i, j, hx, hy, rLocal);
417: nonlinearResidual(-1.0*sc, uLocal, rLocal);
418: f[j][i] += rLocal[0];
419: f[j][i+1] += rLocal[1];
420: f[j+1][i+1] += rLocal[2];
421: f[j+1][i] += rLocal[3];
422: if (i == 0 || j == 0) {
423: ExactSolution(i*hx, j*hy, &uExact);
424: f[j][i] = x[j][i] - uExact;
425: }
426: if ((i == info->mx-2) || (j == 0)) {
427: ExactSolution((i+1)*hx, j*hy, &uExact);
428: f[j][i+1] = x[j][i+1] - uExact;
429: }
430: if ((i == info->mx-2) || (j == info->my-2)) {
431: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
432: f[j+1][i+1] = x[j+1][i+1] - uExact;
433: }
434: if ((i == 0) || (j == info->my-2)) {
435: ExactSolution(i*hx, (j+1)*hy, &uExact);
436: f[j+1][i] = x[j+1][i] - uExact;
437: }
438: }
439: }
441: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
442: return(0);
443: }
447: /*
448: FormJacobianLocal - Evaluates Jacobian matrix.
449: */
450: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
451: {
452: PetscScalar JLocal[16], ELocal[16], uLocal[4];
453: MatStencil rows[4], cols[4], ident;
454: PetscInt localRows[4];
455: PetscScalar alpha,lambda,hx,hy,hxhy,sc;
456: PetscInt i,j,k,l,numRows;
460: alpha = user->alpha;
461: lambda = user->lambda;
462: hx = 1.0/(PetscReal)(info->mx-1);
463: hy = 1.0/(PetscReal)(info->my-1);
464: sc = hx*hy*lambda;
465: hxhy = hx*hy;
467: MatZeroEntries(jac);
468: /*
469: Compute entries for the locally owned part of the Jacobian.
470: - Currently, all PETSc parallel matrix formats are partitioned by
471: contiguous chunks of rows across the processors.
472: - Each processor needs to insert only elements that it owns
473: locally (but any non-local elements will be sent to the
474: appropriate processor during matrix assembly).
475: - Here, we set all entries for a particular row at once.
476: - We can set matrix entries either using either
477: MatSetValuesLocal() or MatSetValues(), as discussed above.
478: */
479: for (j=info->ys; j<info->ys+info->ym-1; j++) {
480: for (i=info->xs; i<info->xs+info->xm-1; i++) {
481: numRows = 0;
482: uLocal[0] = x[j][i];
483: uLocal[1] = x[j][i+1];
484: uLocal[2] = x[j+1][i+1];
485: uLocal[3] = x[j+1][i];
486: /* i,j */
487: if (i == 0 || j == 0) {
488: ident.i = i; ident.j = j;
489: JLocal[0] = 1.0;
490: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
491: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
492: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
493: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
494: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
495: } else {
496: localRows[numRows] = 0;
497: rows[numRows].i = i; rows[numRows].j = j;
498: numRows++;
499: }
500: cols[0].i = i; cols[0].j = j;
501: /* i+1,j */
502: if ((i == info->mx-2) || (j == 0)) {
503: ident.i = i+1; ident.j = j;
504: JLocal[0] = 1.0;
505: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
506: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
507: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
508: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
509: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
510: } else {
511: localRows[numRows] = 1;
512: rows[numRows].i = i+1; rows[numRows].j = j;
513: numRows++;
514: }
515: cols[1].i = i+1; cols[1].j = j;
516: /* i+1,j+1 */
517: if ((i == info->mx-2) || (j == info->my-2)) {
518: ident.i = i+1; ident.j = j+1;
519: JLocal[0] = 1.0;
520: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
521: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
522: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
523: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
524: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
525: } else {
526: localRows[numRows] = 2;
527: rows[numRows].i = i+1; rows[numRows].j = j+1;
528: numRows++;
529: }
530: cols[2].i = i+1; cols[2].j = j+1;
531: /* i,j+1 */
532: if ((i == 0) || (j == info->my-2)) {
533: ident.i = i; ident.j = j+1;
534: JLocal[0] = 1.0;
535: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
536: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
537: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
538: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
539: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
540: } else {
541: localRows[numRows] = 3;
542: rows[numRows].i = i; rows[numRows].j = j+1;
543: numRows++;
544: }
545: cols[3].i = i; cols[3].j = j+1;
546: nonlinearJacobian(-1.0*sc, uLocal, ELocal);
547: for(k = 0; k < numRows; k++) {
548: for(l = 0; l < 4; l++) {
549: JLocal[k*4 + l] = (hxhy*alpha*Kref[localRows[k]*4 + l] + ELocal[localRows[k]*4 + l]);
550: }
551: }
552: MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
553: }
554: }
556: /*
557: Assemble matrix, using the 2-step process:
558: MatAssemblyBegin(), MatAssemblyEnd().
559: */
560: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
561: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
562: /*
563: Tell the matrix we will never add a new nonzero location to the
564: matrix. If we do, it will generate an error.
565: */
566: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
567: return(0);
568: }