Actual source code: ex4.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Solves the Lane-Emden equation in a 2D rectangular\n\
2: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";
4: /*T
5: Concepts: SNES^parallel Lane-Emden example
6: Concepts: DMDA^using distributed arrays;
7: Processors: n
8: T*/
10: /* ------------------------------------------------------------------------
12: The Lane-Emden equation is given by the partial differential equation
14: -alpha*Laplacian u - lambda*u^3 = 0, 0 < x,y < 1,
16: with boundary conditions
18: u = 0 for x = 0, x = 1, y = 0, y = 1.
20: A bilinear finite element approximation is used to discretize the boundary
21: value problem to obtain a nonlinear system of equations.
23: ------------------------------------------------------------------------- */
25: /*
26: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
27: Include "petscsnes.h" so that we can use SNES solvers. Note that this
28: file automatically includes:
29: petscsys.h - base PETSc routines petscvec.h - vectors
30: petscmat.h - matrices
31: petscis.h - index sets petscksp.h - Krylov subspace methods
32: petscviewer.h - viewers petscpc.h - preconditioners
33: petscksp.h - linear solvers
34: */
35: #include <petscsys.h>
36: #include <petscbag.h>
37: #include <petscdmda.h>
38: #include <petscsnes.h>
40: /*
41: User-defined application context - contains data needed by the
42: application-provided call-back routines, FormJacobianLocal() and
43: FormFunctionLocal().
44: */
45: typedef struct {
46: PetscReal alpha; /* parameter controlling linearity */
47: PetscReal lambda; /* parameter controlling nonlinearity */
48: } AppCtx;
50: static PetscScalar Kref[16] = { 0.666667, -0.166667, -0.333333, -0.166667,
51: -0.166667, 0.666667, -0.166667, -0.333333,
52: -0.333333, -0.166667, 0.666667, -0.166667,
53: -0.166667, -0.333333, -0.166667, 0.666667};
55: /* These are */
56: static PetscScalar quadPoints[8] = {0.211325, 0.211325,
57: 0.788675, 0.211325,
58: 0.788675, 0.788675,
59: 0.211325, 0.788675};
60: static PetscScalar quadWeights[4] = {0.25, 0.25, 0.25, 0.25};
62: /*
63: User-defined routines
64: */
65: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
66: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
67: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,MatStructure*,AppCtx*);
68: extern PetscErrorCode PrintVector(DM, Vec);
72: int main(int argc,char **argv)
73: {
74: DM da;
75: SNES snes; /* nonlinear solver */
76: AppCtx *user; /* user-defined work context */
77: PetscBag bag;
78: PetscInt its; /* iterations for convergence */
79: SNESConvergedReason reason;
80: PetscErrorCode ierr;
81: PetscReal lambda_max = 6.81, lambda_min = 0.0;
82: Vec x;
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Initialize program
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: PetscInitialize(&argc,&argv,(char*)0,help);
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Initialize problem parameters
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
93: PetscBagGetData(bag, (void**) &user);
94: PetscBagSetName(bag, "params", "Parameters for SNES example 4");
95: PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
96: PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
97: PetscBagSetFromOptions(bag);
98: PetscOptionsGetReal(NULL,"-alpha",&user->alpha,NULL);
99: PetscOptionsGetReal(NULL,"-lambda",&user->lambda,NULL);
100: if (user->lambda > lambda_max || user->lambda < lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create SNES to manage hierarchical solvers
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: SNESCreate(PETSC_COMM_WORLD,&snes);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create distributed array (DMDA) to manage parallel grid and vectors
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
111: DMDASetFieldName(da, 0, "ooblek");
112: DMSetApplicationContext(da,user);
113: SNESSetDM(snes, (DM) da);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Set the discretization functions
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,user);
119: DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*))FormJacobianLocal,user);
120: SNESSetFromOptions(snes);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Evaluate initial guess
124: Note: The user should initialize the vector, x, with the initial guess
125: for the nonlinear solver prior to calling SNESSolve(). In particular,
126: to employ an initial guess of zero, the user should explicitly set
127: this vector to zero by calling VecSet().
128: */
129: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Solve nonlinear system
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: SNESSolve(snes,NULL,NULL);
135: SNESGetIterationNumber(snes,&its);
136: SNESGetConvergedReason(snes, &reason);
137: DMDestroy(&da);
138: SNESGetDM(snes,&da);
139: SNESGetSolution(snes,&x);
140: PrintVector(da, x);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Free work space. All PETSc objects should be destroyed when they
144: are no longer needed.
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: SNESDestroy(&snes);
147: PetscBagDestroy(&bag);
148: PetscFinalize();
149: return(0);
150: }
154: PetscErrorCode PrintVector(DM da, Vec U)
155: {
156: PetscScalar **u;
157: PetscInt i,j,xs,ys,xm,ym;
161: DMDAVecGetArray(da,U,&u);
162: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
163: for (j = ys+ym-1; j >= ys; j--) {
164: for (i = xs; i < xs+xm; i++) {
165: PetscPrintf(PETSC_COMM_SELF,"u[%d][%d] = %G ", j, i, PetscRealPart(u[j][i]));
166: }
167: PetscPrintf(PETSC_COMM_SELF,"\n");
168: }
169: DMDAVecRestoreArray(da,U,&u);
170: return(0);
171: }
175: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
176: {
178: *u = x*x;
179: return(0);
180: }
184: /*
185: FormInitialGuess - Forms initial approximation.
187: Input Parameters:
188: X - vector
190: Output Parameter:
191: X - vector
192: */
193: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
194: {
195: AppCtx *user;
196: PetscInt i,j,Mx,My,xs,ys,xm,ym;
198: PetscReal lambda,temp1,temp,hx,hy;
199: PetscScalar **x;
200: DM da;
203: SNESGetDM(snes,&da);
204: DMGetApplicationContext(da,&user);
205: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
206: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
208: lambda = user->lambda;
209: hx = 1.0/(PetscReal)(Mx-1);
210: hy = 1.0/(PetscReal)(My-1);
211: if (lambda == 0.0) temp1 = 1.0;
212: else temp1 = lambda/(lambda + 1.0);
214: /*
215: Get a pointer to vector data.
216: - For default PETSc vectors, VecGetArray() returns a pointer to
217: the data array. Otherwise, the routine is implementation dependent.
218: - You MUST call VecRestoreArray() when you no longer need access to
219: the array.
220: */
221: DMDAVecGetArray(da,X,&x);
223: /*
224: Get local grid boundaries (for 2-dimensional DMDA):
225: xs, ys - starting grid indices (no ghost points)
226: xm, ym - widths of local grid (no ghost points)
228: */
229: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
231: /*
232: Compute initial guess over the locally owned part of the grid
233: */
234: for (j=ys; j<ys+ym; j++) {
235: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
236: for (i=xs; i<xs+xm; i++) {
238: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) x[j][i] = 0.0; /* boundary conditions are all zero Dirichlet */
239: else x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
240: }
241: }
243: /*
244: Restore vector
245: */
246: DMDAVecRestoreArray(da,X,&x);
247: return(0);
248: }
252: PetscErrorCode constantResidual(PetscReal lambda, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
253: {
254: PetscScalar rLocal[4] = {0.0, 0.0, 0.0};
255: PetscScalar phi[4] = {0.0, 0.0, 0.0, 0.0};
256: /* PetscReal xI = i*hx, yI = j*hy, x, y; */
257: PetscReal hxhy = hx*hy;
258: PetscScalar res;
259: PetscInt q, k;
262: for (q = 0; q < 4; q++) {
263: phi[0] = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
264: phi[1] = quadPoints[q*2] *(1.0 - quadPoints[q*2+1]);
265: phi[2] = quadPoints[q*2] * quadPoints[q*2+1];
266: phi[3] = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
267: /*
268: x = xI + quadPoints[q*2]*hx;
269: y = yI + quadPoints[q*2+1]*hy;
270: */
271: res = quadWeights[q]*(2.0);
272: for (k = 0; k < 4; k++) rLocal[k] += phi[k]*res;
273: }
274: for (k = 0; k < 4; k++) r[k] += lambda*hxhy*rLocal[k];
275: return(0);
276: }
280: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[])
281: {
283: 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])
284: + 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])
285: + 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;
286: 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])
287: + 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])
288: + 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;
289: 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]))
290: + 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;
291: 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])
292: + 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])
293: + 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;
294: return(0);
295: }
299: PetscErrorCode nonlinearResidualBratu(PetscReal lambda, PetscScalar u[], PetscScalar r[])
300: {
301: PetscScalar rLocal[4] = {0.0, 0.0, 0.0, 0.0};
302: PetscScalar phi[4] = {0.0, 0.0, 0.0, 0.0};
303: PetscScalar res;
304: PetscInt q;
307: for (q = 0; q < 4; q++) {
308: phi[0] = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
309: phi[1] = quadPoints[q*2] *(1.0 - quadPoints[q*2+1]);
310: phi[2] = quadPoints[q*2] * quadPoints[q*2+1];
311: phi[3] = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
312: res = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
313: rLocal[0] += phi[0]*res;
314: rLocal[1] += phi[1]*res;
315: rLocal[2] += phi[2]*res;
316: rLocal[3] += phi[3]*res;
317: }
318: r[0] += lambda*rLocal[0];
319: r[1] += lambda*rLocal[1];
320: r[2] += lambda*rLocal[2];
321: r[3] += lambda*rLocal[3];
322: return(0);
323: }
327: PetscErrorCode nonlinearJacobian(PetscScalar lambda, PetscScalar u[], PetscScalar J[])
328: {
330: 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;
331: 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;
332: 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;
333: 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;
335: 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;
336: 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;
337: 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;
338: 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;
340: 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;
341: 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;
342: 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;
343: 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;
345: 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;
346: 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;
347: 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;
348: 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;
349: return(0);
350: }
354: /*
355: FormFunctionLocal - Evaluates nonlinear function, F(x).
357: */
358: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
359: {
360: PetscScalar uLocal[4];
361: PetscScalar rLocal[4];
362: PetscScalar uExact;
363: PetscReal alpha,lambda,hx,hy,hxhy,sc;
364: PetscInt i,j,k,l;
368: alpha = user->alpha;
369: lambda = user->lambda;
370: hx = 1.0/(PetscReal)(info->mx-1);
371: hy = 1.0/(PetscReal)(info->my-1);
372: sc = hx*hy*lambda;
373: hxhy = hx*hy;
375: /* Zero the vector */
376: PetscMemzero((void*) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(PetscScalar));
377: /* Compute function over the locally owned part of the grid. For each
378: vertex (i,j), we consider the element below:
380: 3 2
381: i,j+1 --- i+1,j+1
382: | |
383: | |
384: i,j --- i+1,j
385: 0 1
387: and therefore we do not loop over the last vertex in each dimension.
388: */
389: for (j = info->ys; j < info->ys+info->ym-1; j++) {
390: for (i = info->xs; i < info->xs+info->xm-1; i++) {
391: uLocal[0] = x[j][i];
392: uLocal[1] = x[j][i+1];
393: uLocal[2] = x[j+1][i+1];
394: uLocal[3] = x[j+1][i];
395: for (k = 0; k < 4; k++) {
396: rLocal[k] = 0.0;
397: for (l = 0; l < 4; l++) rLocal[k] += Kref[k*4 + l]*uLocal[l];
398: rLocal[k] *= hxhy*alpha;
399: }
400: constantResidual(1.0, i, j, hx, hy, rLocal);
401: nonlinearResidual(-1.0*sc, uLocal, rLocal);
403: f[j][i] += rLocal[0];
404: f[j][i+1] += rLocal[1];
405: f[j+1][i+1] += rLocal[2];
406: f[j+1][i] += rLocal[3];
408: if (i == 0 || j == 0) {
409: ExactSolution(i*hx, j*hy, &uExact);
411: f[j][i] = x[j][i] - uExact;
412: }
413: if ((i == info->mx-2) || (j == 0)) {
414: ExactSolution((i+1)*hx, j*hy, &uExact);
416: f[j][i+1] = x[j][i+1] - uExact;
417: }
418: if ((i == info->mx-2) || (j == info->my-2)) {
419: ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
421: f[j+1][i+1] = x[j+1][i+1] - uExact;
422: }
423: if ((i == 0) || (j == info->my-2)) {
424: ExactSolution(i*hx, (j+1)*hy, &uExact);
426: f[j+1][i] = x[j+1][i] - uExact;
427: }
428: }
429: }
431: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
432: return(0);
433: }
437: /*
438: FormJacobianLocal - Evaluates Jacobian matrix.
439: */
440: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat A,Mat jac,MatStructure *str,AppCtx *user)
441: {
442: PetscScalar JLocal[16], ELocal[16], uLocal[4];
443: MatStencil rows[4], cols[4], ident;
444: PetscInt localRows[4];
445: PetscScalar alpha,lambda,hx,hy,hxhy,sc;
446: PetscInt i,j,k,l,numRows;
450: alpha = user->alpha;
451: lambda = user->lambda;
452: hx = 1.0/(PetscReal)(info->mx-1);
453: hy = 1.0/(PetscReal)(info->my-1);
454: sc = hx*hy*lambda;
455: hxhy = hx*hy;
457: MatZeroEntries(jac);
458: /*
459: Compute entries for the locally owned part of the Jacobian.
460: - Currently, all PETSc parallel matrix formats are partitioned by
461: contiguous chunks of rows across the processors.
462: - Each processor needs to insert only elements that it owns
463: locally (but any non-local elements will be sent to the
464: appropriate processor during matrix assembly).
465: - Here, we set all entries for a particular row at once.
466: - We can set matrix entries either using either
467: MatSetValuesLocal() or MatSetValues(), as discussed above.
468: */
469: for (j=info->ys; j<info->ys+info->ym-1; j++) {
470: for (i=info->xs; i<info->xs+info->xm-1; i++) {
471: numRows = 0;
472: uLocal[0] = x[j][i];
473: uLocal[1] = x[j][i+1];
474: uLocal[2] = x[j+1][i+1];
475: uLocal[3] = x[j+1][i];
476: /* i,j */
477: if (i == 0 || j == 0) {
478: ident.i = i; ident.j = j;
479: JLocal[0] = 1.0;
481: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
482: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
483: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
484: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
485: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
486: } else {
487: localRows[numRows] = 0;
488: rows[numRows].i = i; rows[numRows].j = j;
489: numRows++;
490: }
491: cols[0].i = i; cols[0].j = j;
492: /* i+1,j */
493: if ((i == info->mx-2) || (j == 0)) {
494: ident.i = i+1; ident.j = j;
495: JLocal[0] = 1.0;
497: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
498: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
499: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
500: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
501: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
502: } else {
503: localRows[numRows] = 1;
504: rows[numRows].i = i+1; rows[numRows].j = j;
505: numRows++;
506: }
507: cols[1].i = i+1; cols[1].j = j;
508: /* i+1,j+1 */
509: if ((i == info->mx-2) || (j == info->my-2)) {
510: ident.i = i+1; ident.j = j+1;
511: JLocal[0] = 1.0;
513: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
514: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
515: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
516: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
517: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
518: } else {
519: localRows[numRows] = 2;
520: rows[numRows].i = i+1; rows[numRows].j = j+1;
521: numRows++;
522: }
523: cols[2].i = i+1; cols[2].j = j+1;
524: /* i,j+1 */
525: if ((i == 0) || (j == info->my-2)) {
526: ident.i = i; ident.j = j+1;
527: JLocal[0] = 1.0;
529: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
530: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
531: MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
532: MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
533: MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
534: } else {
535: localRows[numRows] = 3;
536: rows[numRows].i = i; rows[numRows].j = j+1;
537: numRows++;
538: }
539: cols[3].i = i; cols[3].j = j+1;
540: nonlinearJacobian(-1.0*sc, uLocal, ELocal);
541: for (k = 0; k < numRows; k++) {
542: for (l = 0; l < 4; l++) {
543: JLocal[k*4 + l] = (hxhy*alpha*Kref[localRows[k]*4 + l] + ELocal[localRows[k]*4 + l]);
544: }
545: }
546: MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
547: }
548: }
550: /*
551: Assemble matrix, using the 2-step process:
552: MatAssemblyBegin(), MatAssemblyEnd().
553: */
554: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
555: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
556: if (A != jac) {
557: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
558: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
559: }
560: *str = SAME_NONZERO_PATTERN;
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: }