Actual source code: ex13.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "This program is a replica of ex6.c except that it does 2 solves to avoid paging.\n\
3: This program demonstrates use of the SNES package to solve systems of\n\
4: nonlinear equations in parallel, using 2-dimensional distributed arrays.\n\
5: The 2-dim Bratu (SFI - solid fuel ignition) test problem is used, where\n\
6: analytic formation of the Jacobian is the default. The command line\n\
7: options are:\n\
8: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
9: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
10: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
11: -my <yg>, where <yg> = number of grid points in the y-direction\n\
12: -Nx <npx>, where <npx> = number of processors in the x-direction\n\
13: -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";
15: /*
16: 1) Solid Fuel Ignition (SFI) problem. This problem is modeled by
17: the partial differential equation
18:
19: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
20:
21: with boundary conditions
22:
23: u = 0 for x = 0, x = 1, y = 0, y = 1.
24:
25: A finite difference approximation with the usual 5-point stencil
26: is used to discretize the boundary value problem to obtain a nonlinear
27: system of equations.
28: */
30: #include <petscsnes.h>
31: #include <petscdmda.h>
33: /* User-defined application context */
34: typedef struct {
35: PetscReal param; /* test problem parameter */
36: PetscInt mx,my; /* discretization in x, y directions */
37: Vec localX,localF; /* ghosted local vector */
38: DM da; /* distributed array data structure */
39: } AppCtx;
41: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*),FormInitialGuess1(AppCtx*,Vec);
42: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
46: int main(int argc,char **argv)
47: {
48: SNES snes; /* nonlinear solver */
49: const SNESType type = SNESLS; /* nonlinear solution method */
50: Vec x,r; /* solution, residual vectors */
51: Mat J; /* Jacobian matrix */
52: AppCtx user; /* user-defined work context */
53: PetscInt i,its,N,Nx = PETSC_DECIDE,Ny = PETSC_DECIDE;
55: PetscBool matrix_free = PETSC_FALSE;
56: PetscMPIInt size;
57: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
58: #if defined(PETSC_USE_LOG)
59: PetscLogStage stages[2];
60: #endif
62: PetscInitialize(&argc,&argv,(char *)0,help);
64: PetscLogStageRegister("stage 1",&stages[0]);
65: PetscLogStageRegister("stage 2",&stages[1]);
66: for (i=0; i<2; i++) {
67: PetscLogStagePush(stages[i]);
68: user.mx = 4; user.my = 4; user.param = 6.0;
69:
70: if (i) {
71: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
72: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
73: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
74: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
75: }
76: N = user.mx*user.my;
78: MPI_Comm_size(PETSC_COMM_WORLD,&size);
79: PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
80: PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
81: if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE)) SETERRQ(PETSC_COMM_SELF,1,"Incompatible number of processors: Nx * Ny != size");
82:
83: /* Set up distributed array */
84: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,
85: PETSC_NULL,PETSC_NULL,&user.da);
86: DMCreateGlobalVector(user.da,&x);
87: VecDuplicate(x,&r);
88: DMCreateLocalVector(user.da,&user.localX);
89: VecDuplicate(user.localX,&user.localF);
91: /* Create nonlinear solver and set function evaluation routine */
92: SNESCreate(PETSC_COMM_WORLD,&snes);
93: SNESSetType(snes,type);
94: SNESSetFunction(snes,r,FormFunction1,&user);
96: /* Set default Jacobian evaluation routine. User can override with:
97: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
98: (unless user explicitly sets preconditioner)
99: -snes_fd : default finite differencing approximation of Jacobian
100: */
101: PetscOptionsGetBool(PETSC_NULL,"-snes_mf",&matrix_free,PETSC_NULL);
102: if (!matrix_free) {
103: PetscBool matrix_free_operator = PETSC_FALSE;
104: PetscOptionsGetBool(PETSC_NULL,"-snes_mf_operator",&matrix_free_operator,PETSC_NULL);
105: if (matrix_free_operator) matrix_free = PETSC_FALSE;
106: }
107: if (!matrix_free) {
108: MatCreate(PETSC_COMM_WORLD,&J);
109: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
110: MatSetFromOptions(J);
111: SNESSetJacobian(snes,J,J,FormJacobian1,&user);
112: }
114: /* Set PetscOptions, then solve nonlinear system */
115: SNESSetFromOptions(snes);
116: FormInitialGuess1(&user,x);
117: SNESSolve(snes,PETSC_NULL,x);
118: SNESGetIterationNumber(snes,&its);
119: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
121: /* Free data structures */
122: if (!matrix_free) {
123: MatDestroy(&J);
124: }
125: VecDestroy(&x);
126: VecDestroy(&r);
127: VecDestroy(&user.localX);
128: VecDestroy(&user.localF);
129: SNESDestroy(&snes);
130: DMDestroy(&user.da);
131: }
132: PetscFinalize();
134: return 0;
135: }/* -------------------- Form initial approximation ----------------- */
138: PetscErrorCode FormInitialGuess1(AppCtx *user,Vec X)
139: {
140: PetscInt i,j,row,mx,my,xs,ys,xm,ym,Xm,Ym,Xs,Ys;
142: PetscReal one = 1.0,lambda,temp1,temp,hx,hy;
143: PetscScalar *x;
144: Vec localX = user->localX;
146: mx = user->mx; my = user->my; lambda = user->param;
147: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
149: /* Get ghost points */
150: VecGetArray(localX,&x);
151: temp1 = lambda/(lambda + one);
152: DMDAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
153: DMDAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);
155: /* Compute initial guess */
156: for (j=ys; j<ys+ym; j++) {
157: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
158: for (i=xs; i<xs+xm; i++) {
159: row = i - Xs + (j - Ys)*Xm;
160: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
161: x[row] = 0.0;
162: continue;
163: }
164: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
165: }
166: }
167: VecRestoreArray(localX,&x);
169: /* Insert values into global vector */
170: DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);
171: DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);
172: return 0;
173: } /* -------------------- Evaluate Function F(x) --------------------- */
176: PetscErrorCode FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
177: {
178: AppCtx *user = (AppCtx*)ptr;
180: PetscInt i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym;
181: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
182: PetscScalar u,uxx,uyy,*x,*f;
183: Vec localX = user->localX,localF = user->localF;
185: mx = user->mx; my = user->my; lambda = user->param;
186: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
187: sc = hx*hy*lambda; hxdhy = hx/hy; hydhx = hy/hx;
189: /* Get ghost points */
190: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
191: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
192: VecGetArray(localX,&x);
193: VecGetArray(localF,&f);
194: DMDAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
195: DMDAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);
197: /* Evaluate function */
198: for (j=ys; j<ys+ym; j++) {
199: row = (j - Ys)*Xm + xs - Xs - 1;
200: for (i=xs; i<xs+xm; i++) {
201: row++;
202: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
203: f[row] = x[row];
204: continue;
205: }
206: u = x[row];
207: uxx = (two*u - x[row-1] - x[row+1])*hydhx;
208: uyy = (two*u - x[row-Xm] - x[row+Xm])*hxdhy;
209: f[row] = uxx + uyy - sc*PetscExpScalar(u);
210: }
211: }
212: VecRestoreArray(localX,&x);
213: VecRestoreArray(localF,&f);
215: /* Insert values into global vector */
216: DMLocalToGlobalBegin(user->da,localF,INSERT_VALUES,F);
217: DMLocalToGlobalEnd(user->da,localF,INSERT_VALUES,F);
218: PetscLogFlops(11.0*ym*xm);
219: return 0;
220: } /* -------------------- Evaluate Jacobian F'(x) --------------------- */
223: PetscErrorCode FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
224: {
225: AppCtx *user = (AppCtx*)ptr;
226: Mat jac = *J;
228: PetscInt i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym,col[5];
229: PetscInt nloc,*ltog,grow;
230: PetscScalar two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
231: Vec localX = user->localX;
233: mx = user->mx; my = user->my; lambda = user->param;
234: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
235: sc = hx*hy; hxdhy = hx/hy; hydhx = hy/hx;
237: /* Get ghost points */
238: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
239: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
240: VecGetArray(localX,&x);
241: DMDAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
242: DMDAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);
243: DMDAGetGlobalIndices(user->da,&nloc,<og);
245: /* Evaluate function */
246: for (j=ys; j<ys+ym; j++) {
247: row = (j - Ys)*Xm + xs - Xs - 1;
248: for (i=xs; i<xs+xm; i++) {
249: row++;
250: grow = ltog[row];
251: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
252: MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
253: continue;
254: }
255: v[0] = -hxdhy; col[0] = ltog[row - Xm];
256: v[1] = -hydhx; col[1] = ltog[row - 1];
257: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
258: v[3] = -hydhx; col[3] = ltog[row + 1];
259: v[4] = -hxdhy; col[4] = ltog[row + Xm];
260: MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
261: }
262: }
263: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
264: VecRestoreArray(X,&x);
265: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
266: *flag = SAME_NONZERO_PATTERN;
267: return 0;
268: }