Actual source code: ex46.c
1: static char help[] = "Surface processes in geophysics.\n\n";
3: /*T
4: Concepts: SNES^parallel Surface process example
5: Concepts: DMDA^using distributed arrays;
6: Concepts: IS coloirng types;
7: Processors: n
8: T*/
10: /* ------------------------------------------------------------------------
12: Solid Fuel Ignition (SFI) problem. This problem is modeled by
13: the partial differential equation
14:
15: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
16:
17: with boundary conditions
18:
19: u = 0 for x = 0, x = 1, y = 0, y = 1.
20:
21: A finite difference approximation with the usual 5-point stencil
22: is used to discretize the boundary value problem to obtain a nonlinear
23: system of equations.
25: Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options]
26: e.g.,
27: ./ex5 -fd_jacobian -mat_fd_coloring_view_draw -draw_pause -1
28: mpiexec -n 2 ./ex5 -fd_jacobian_ghosted -log_summary
30: ------------------------------------------------------------------------- */
32: /*
33: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
34: Include "petscsnes.h" so that we can use SNES solvers. Note that this
35: file automatically includes:
36: petscsys.h - base PETSc routines petscvec.h - vectors
37: petscmat.h - matrices
38: petscis.h - index sets petscksp.h - Krylov subspace methods
39: petscviewer.h - viewers petscpc.h - preconditioners
40: petscksp.h - linear solvers
41: */
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: DM da; /* distributed array data structure */
52: PassiveReal D; /* The diffusion coefficient */
53: PassiveReal K; /* The advection coefficient */
54: PetscInt m; /* Exponent for A */
55: } AppCtx;
57: /*
58: User-defined routines
59: */
66: int main(int argc,char **argv)
67: {
68: DMMG *dmmg;
69: SNES snes; /* nonlinear solver */
70: AppCtx user; /* user-defined work context */
71: PetscInt its; /* iterations for convergence */
72: PetscErrorCode ierr;
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Initialize program
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: PetscInitialize(&argc,&argv,(char *)0,help);
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Initialize problem parameters
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Surface Process Problem Options", "DMMG");
84: user.D = 1.0;
85: PetscOptionsReal("-D", "The diffusion coefficient D", __FILE__, user.D, &user.D, PETSC_NULL);
86: user.K = 1.0;
87: PetscOptionsReal("-K", "The advection coefficient K", __FILE__, user.K, &user.K, PETSC_NULL);
88: user.m = 1;
89: PetscOptionsInt("-m", "The exponent for A", __FILE__, user.m, &user.m, PETSC_NULL);
90: PetscOptionsEnd();
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create distributed array (DMDA) to manage parallel grid and vectors
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
96: 1,1,PETSC_NULL,PETSC_NULL,&user.da);
97: DMDASetUniformCoordinates(user.da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
98: DMMGCreate(PETSC_COMM_WORLD, 1, &user, &dmmg);
99: DMMGSetDM(dmmg, (DM) user.da);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Set local function evaluation routine
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Customize solver; set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: DMMGSetFromOptions(dmmg);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Form initial guess
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: //FormInitialGuess(&user,DMMGGetx(dmmg));
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Solve nonlinear system
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: DMMGSolve(dmmg);
120: snes = DMMGGetSNES(dmmg);
121: SNESGetIterationNumber(snes,&its);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
127: VecAssemblyBegin(DMMGGetx(dmmg));
128: VecAssemblyEnd(DMMGGetx(dmmg));
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Free work space. All PETSc objects should be destroyed when they
131: are no longer needed.
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: DMMGDestroy(dmmg);
135: DMDestroy(&user.da);
137: PetscFinalize();
138: return(0);
139: }
140: /* ------------------------------------------------------------------- */
143: /*
144: FormInitialGuess - Forms initial approximation.
146: Input Parameters:
147: user - user-defined application context
148: X - vector
150: Output Parameter:
151: X - vector
152: */
153: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
154: {
155: PetscInt i,j,Mx,My,xs,ys,xm,ym;
157: PetscReal D,temp1,temp,hx,hy;
158: PetscScalar **x;
161: DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
162: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
164: D = user->D;
165: hx = 1.0/(PetscReal)(Mx-1);
166: hy = 1.0/(PetscReal)(My-1);
167: temp1 = D/(D + 1.0);
169: /*
170: Get a pointer to vector data.
171: - For default PETSc vectors, VecGetArray() returns a pointer to
172: the data array. Otherwise, the routine is implementation dependent.
173: - You MUST call VecRestoreArray() when you no longer need access to
174: the array.
175: */
176: DMDAVecGetArray(user->da,X,&x);
178: /*
179: Get local grid boundaries (for 2-dimensional DMDA):
180: xs, ys - starting grid indices (no ghost points)
181: xm, ym - widths of local grid (no ghost points)
183: */
184: DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
186: /*
187: Compute initial guess over the locally owned part of the grid
188: */
189: for (j=ys; j<ys+ym; j++) {
190: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
191: for (i=xs; i<xs+xm; i++) {
192: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
193: /* boundary conditions are all zero Dirichlet */
194: x[j][i] = 0.0;
195: } else {
196: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
197: }
198: }
199: }
201: /*
202: Restore vector
203: */
204: DMDAVecRestoreArray(user->da,X,&x);
206: return(0);
207: }
211: PetscScalar funcU(DMDACoor2d *coords)
212: {
213: return coords->x + coords->y;
214: }
218: PetscScalar funcA(PetscScalar z, AppCtx *user)
219: {
220: PetscScalar v = 1.0;
221: PetscInt i;
223: for(i = 0; i < user->m; ++i) {
224: v *= z;
225: }
226: return v;
227: }
231: PetscScalar funcADer(PetscScalar z, AppCtx *user)
232: {
233: PetscScalar v = 1.0;
234: PetscInt i;
236: for(i = 0; i < user->m-1; ++i) {
237: v *= z;
238: }
239: return user->m*v;
240: }
244: /*
245: FormFunctionLocal - Evaluates nonlinear function, F(x).
246: */
247: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
248: {
249: DM coordDA;
250: Vec coordinates;
251: DMDACoor2d **coords;
252: PetscScalar u, ux, uy, uxx, uyy;
253: PetscReal D, K, hx, hy, hxdhy, hydhx;
254: PetscInt i,j;
259: D = user->D;
260: K = user->K;
261: hx = 1.0/(PetscReal)(info->mx-1);
262: hy = 1.0/(PetscReal)(info->my-1);
263: hxdhy = hx/hy;
264: hydhx = hy/hx;
265: /*
266: Compute function over the locally owned part of the grid
267: */
268: DMDAGetCoordinateDA(user->da, &coordDA);
269: DMDAGetCoordinates(user->da, &coordinates);
270: DMDAVecGetArray(coordDA, coordinates, &coords);
271: for (j=info->ys; j<info->ys+info->ym; j++) {
272: for (i=info->xs; i<info->xs+info->xm; i++) {
273: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
274: f[j][i] = x[j][i];
275: } else {
276: u = x[j][i];
277: ux = (x[j][i+1] - x[j][i])/hx;
278: uy = (x[j+1][i] - x[j][i])/hy;
279: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
280: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
281: f[j][i] = D*(uxx + uyy) - (K*funcA(x[j][i], user)*sqrt(ux*ux + uy*uy) + funcU(&coords[j][i]))*hx*hy;
282: if (PetscIsInfOrNanScalar(f[j][i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(f[j][i]));
283: }
284: }
285: }
286: DMDAVecRestoreArray(coordDA, coordinates, &coords);
287: PetscLogFlops(11*info->ym*info->xm);
288: return(0);
289: }
293: /*
294: FormJacobianLocal - Evaluates Jacobian matrix.
295: */
296: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
297: {
298: MatStencil col[5], row;
299: PetscScalar D, K, A, v[5], hx, hy, hxdhy, hydhx, ux, uy, normGradZ;
300: PetscInt i, j,k;
304: D = user->D;
305: K = user->K;
306: hx = 1.0/(PetscReal)(info->mx-1);
307: hy = 1.0/(PetscReal)(info->my-1);
308: hxdhy = hx/hy;
309: hydhx = hy/hx;
311: /*
312: Compute entries for the locally owned part of the Jacobian.
313: - Currently, all PETSc parallel matrix formats are partitioned by
314: contiguous chunks of rows across the processors.
315: - Each processor needs to insert only elements that it owns
316: locally (but any non-local elements will be sent to the
317: appropriate processor during matrix assembly).
318: - Here, we set all entries for a particular row at once.
319: - We can set matrix entries either using either
320: MatSetValuesLocal() or MatSetValues(), as discussed above.
321: */
322: for (j=info->ys; j<info->ys+info->ym; j++) {
323: for (i=info->xs; i<info->xs+info->xm; i++) {
324: row.j = j; row.i = i;
325: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
326: /* boundary points */
327: v[0] = 1.0;
328: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
329: } else {
330: /* interior grid points */
331: ux = (x[j][i+1] - x[j][i])/hx;
332: uy = (x[j+1][i] - x[j][i])/hy;
333: normGradZ = sqrt(ux*ux + uy*uy);
334: //PetscPrintf(PETSC_COMM_SELF, "i: %d j: %d normGradZ: %g\n", i, j, normGradZ);
335: if (normGradZ < 1.0e-8) {
336: normGradZ = 1.0e-8;
337: }
338: A = funcA(x[j][i], user);
340: v[0] = -D*hxdhy; col[0].j = j - 1; col[0].i = i;
341: v[1] = -D*hydhx; col[1].j = j; col[1].i = i-1;
342: v[2] = D*2.0*(hydhx + hxdhy) + K*(funcADer(x[j][i], user)*normGradZ - A/normGradZ)*hx*hy; col[2].j = row.j; col[2].i = row.i;
343: v[3] = -D*hydhx + K*A*hx*hy/(2.0*normGradZ); col[3].j = j; col[3].i = i+1;
344: v[4] = -D*hxdhy + K*A*hx*hy/(2.0*normGradZ); col[4].j = j + 1; col[4].i = i;
345: for(k = 0; k < 5; ++k) {
346: if (PetscIsInfOrNanScalar(v[k])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(v[k]));
347: }
348: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
349: }
350: }
351: }
353: /*
354: Assemble matrix, using the 2-step process:
355: MatAssemblyBegin(), MatAssemblyEnd().
356: */
357: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
358: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
359: /*
360: Tell the matrix we will never add a new nonzero location to the
361: matrix. If we do, it will generate an error.
362: */
363: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
364: return(0);
365: }