Actual source code: ex9.c
petsc-3.3-p6 2013-02-11
2: static char help[] = "This program demonstrates use of the SNES package. Solve systems of\n\
3: nonlinear equations in parallel. This example uses matrix-free Krylov\n\
4: Newton methods with no preconditioner to solve the Bratu (SFI - solid fuel\n\
5: ignition) test problem. The command line options are:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
8: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
9: -my <yg>, where <yg> = number of grid points in the y-direction\n\
10: -mz <zg>, where <zg> = number of grid points in the z-direction\n\n";
12: /*
13: 1) Solid Fuel Ignition (SFI) problem. This problem is modeled by
14: the partial differential equation
15:
16: -Laplacian u - lambda*exp(u) = 0, 0 < x,y,z < 1,
17:
18: with boundary conditions
19:
20: u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
21:
22: A finite difference approximation with the usual 7-point stencil
23: is used to discretize the boundary value problem to obtain a nonlinear
24: system of equations.
25: */
27: #include <petscsnes.h>
28: #include <petscdmda.h>
30: typedef struct {
31: PetscReal param; /* test problem nonlinearity parameter */
32: PetscInt mx,my,mz; /* discretization in x,y,z-directions */
33: Vec localX,localF; /* ghosted local vectors */
34: DM da; /* distributed array datastructure */
35: } AppCtx;
37: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*),FormInitialGuess1(AppCtx*,Vec);
41: int main(int argc,char **argv)
42: {
43: SNES snes; /* nonlinear solver */
44: KSP ksp; /* linear solver */
45: PC pc; /* preconditioner */
46: Mat J; /* Jacobian matrix */
47: AppCtx user; /* user-defined application context */
48: Vec x,r; /* vectors */
49: DMDAStencilType stencil = DMDA_STENCIL_BOX;
51: PetscBool flg;
52: PetscInt Nx = PETSC_DECIDE,Ny = PETSC_DECIDE,Nz = PETSC_DECIDE,its;
53: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
55: PetscInitialize(&argc,&argv,(char *)0,help);
56: PetscOptionsHasName(PETSC_NULL,"-star",&flg);
57: if (flg) stencil = DMDA_STENCIL_STAR;
59: user.mx = 4;
60: user.my = 4;
61: user.mz = 4;
62: user.param = 6.0;
63: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
64: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
65: PetscOptionsGetInt(PETSC_NULL,"-mz",&user.mz,PETSC_NULL);
66: PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
67: PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
68: PetscOptionsGetInt(PETSC_NULL,"-Nz",&Nz,PETSC_NULL);
69: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
70: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
71:
72: /* Set up distributed array */
73: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,stencil,user.mx,user.my,user.mz,
74: Nx,Ny,Nz,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da);
75: DMCreateGlobalVector(user.da,&x);
76: VecDuplicate(x,&r);
77: DMCreateLocalVector(user.da,&user.localX);
78: VecDuplicate(user.localX,&user.localF);
80: /* Create nonlinear solver */
81: SNESCreate(PETSC_COMM_WORLD,&snes);
82: /* Set various routines and options */
83: SNESSetFunction(snes,r,FormFunction1,(void*)&user);
84: MatCreateSNESMF(snes,&J);
85: SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,&user);
86: SNESSetFromOptions(snes);
88: /* Force no preconditioning to be used. Note that this overrides whatever
89: choices may have been specified in the options database. */
90: SNESGetKSP(snes,&ksp);
91: KSPGetPC(ksp,&pc);
92: PCSetType(pc,PCNONE);
94: /* Solve nonlinear system */
95: FormInitialGuess1(&user,x);
96: SNESSolve(snes,PETSC_NULL,x);
97: SNESGetIterationNumber(snes,&its);
98: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
100: /* Free data structures */
101: VecDestroy(&user.localX);
102: VecDestroy(&user.localF);
103: DMDestroy(&user.da);
104: VecDestroy(&x); VecDestroy(&r);
105: MatDestroy(&J); SNESDestroy(&snes);
107: PetscFinalize();
108: return 0;
109: }/* -------------------- Form initial approximation ----------------- */
112: PetscErrorCode FormInitialGuess1(AppCtx *user,Vec X)
113: {
114: PetscInt i,j,k,loc,mx,my,mz,xs,ys,zs,xm,ym,zm,Xm,Ym,Zm,Xs,Ys,Zs,base1;
116: PetscReal one = 1.0,lambda,temp1,temp,Hx,Hy;
117: PetscScalar *x;
118: Vec localX = user->localX;
120: mx = user->mx; my = user->my; mz = user->mz; lambda = user->param;
121: Hx = one / (PetscReal)(mx-1); Hy = one / (PetscReal)(my-1);
123: VecGetArray(localX,&x);
124: temp1 = lambda/(lambda + one);
125: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
126: DMDAGetGhostCorners(user->da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm);
127:
128: for (k=zs; k<zs+zm; k++) {
129: base1 = (Xm*Ym)*(k-Zs);
130: for (j=ys; j<ys+ym; j++) {
131: temp = (PetscReal)(PetscMin(j,my-j-1))*Hy;
132: for (i=xs; i<xs+xm; i++) {
133: loc = base1 + i-Xs + (j-Ys)*Xm;
134: if (i == 0 || j == 0 || k == 0 || i==mx-1 || j==my-1 || k==mz-1) {
135: x[loc] = 0.0;
136: continue;
137: }
138: x[loc] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*Hx,temp));
139: }
140: }
141: }
143: VecRestoreArray(localX,&x);
144: /* stick values into global vector */
145: DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);
146: DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);
147: return 0;
148: }/* -------------------- Evaluate Function F(x) --------------------- */
151: PetscErrorCode FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
152: {
153: AppCtx *user = (AppCtx*)ptr;
155: PetscInt i,j,k,loc,mx,my,mz,xs,ys,zs,xm,ym,zm,Xs,Ys,Zs,Xm,Ym,Zm,base1,base2;
156: PetscReal two = 2.0,one = 1.0,lambda,Hx,Hy,Hz,HxHzdHy,HyHzdHx,HxHydHz;
157: PetscScalar u,uxx,uyy,sc,*x,*f,uzz;
158: Vec localX = user->localX,localF = user->localF;
160: mx = user->mx; my = user->my; mz = user->mz; lambda = user->param;
161: Hx = one / (PetscReal)(mx-1);
162: Hy = one / (PetscReal)(my-1);
163: Hz = one / (PetscReal)(mz-1);
164: sc = Hx*Hy*Hz*lambda; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
165: HxHydHz = Hx*Hy/Hz;
167: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
168: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
169: VecGetArray(localX,&x);
170: VecGetArray(localF,&f);
172: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
173: DMDAGetGhostCorners(user->da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm);
175: for (k=zs; k<zs+zm; k++) {
176: base1 = (Xm*Ym)*(k-Zs);
177: for (j=ys; j<ys+ym; j++) {
178: base2 = base1 + (j-Ys)*Xm;
179: for (i=xs; i<xs+xm; i++) {
180: loc = base2 + (i-Xs);
181: if (i == 0 || j == 0 || k== 0 || i == mx-1 || j == my-1 || k == mz-1) {
182: f[loc] = x[loc];
183: }
184: else {
185: u = x[loc];
186: uxx = (two*u - x[loc-1] - x[loc+1])*HyHzdHx;
187: uyy = (two*u - x[loc-Xm] - x[loc+Xm])*HxHzdHy;
188: uzz = (two*u - x[loc-Xm*Ym] - x[loc+Xm*Ym])*HxHydHz;
189: f[loc] = uxx + uyy + uzz - sc*PetscExpScalar(u);
190: }
191: }
192: }
193: }
194: VecRestoreArray(localX,&x);
195: VecRestoreArray(localF,&f);
196: /* stick values into global vector */
197: DMLocalToGlobalBegin(user->da,localF,INSERT_VALUES,F);
198: DMLocalToGlobalEnd(user->da,localF,INSERT_VALUES,F);
199: PetscLogFlops(11.0*ym*xm*zm);
200: return 0;
201: }
202:
207: