Actual source code: ex1.c
1: /*$Id: ex1.c,v 1.31 2001/08/07 21:31:30 bsmith Exp $*/
3: static char help[] ="Solves the time dependent Bratu problem using pseudo-timestepping.";
5: /*
6: Concepts: TS^pseudo-timestepping
7: Concepts: pseudo-timestepping
8: Concepts: nonlinear problems
9: Processors: 1
11: */
13: /* ------------------------------------------------------------------------
15: This code demonstrates how one may solve a nonlinear problem
16: with pseudo-timestepping. In this simple example, the pseudo-timestep
17: is the same for all grid points, i.e., this is equivalent to using
18: the backward Euler method with a variable timestep.
20: Note: This example does not require pseudo-timestepping since it
21: is an easy nonlinear problem, but it is included to demonstrate how
22: the pseudo-timestepping may be done.
24: See snes/examples/tutorials/ex4.c[ex4f.F] and
25: snes/examples/tutorials/ex5.c[ex5f.F] where the problem is described
26: and solved using Newton's method alone.
28: ----------------------------------------------------------------------------- */
29: /*
30: Include "petscts.h" to use the PETSc timestepping routines. Note that
31: this file automatically includes "petsc.h" and other lower-level
32: PETSc include files.
33: */
34: #include petscts.h
36: /*
37: Create an application context to contain data needed by the
38: application-provided call-back routines, FormJacobian() and
39: FormFunction().
40: */
41: typedef struct {
42: PetscReal param; /* test problem parameter */
43: int mx; /* Discretization in x-direction */
44: int my; /* Discretization in y-direction */
45: } AppCtx;
47: /*
48: User-defined routines
49: */
50: extern int FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*),
51: FormFunction(TS,PetscReal,Vec,Vec,void*),
52: FormInitialGuess(Vec,AppCtx*);
54: #undef __FUNCT__
56: int main(int argc,char **argv)
57: {
58: TS ts; /* timestepping context */
59: Vec x,r; /* solution, residual vectors */
60: Mat J; /* Jacobian matrix */
61: AppCtx user; /* user-defined work context */
62: int its; /* iterations for convergence */
63: int ierr,N;
64: PetscReal param_max = 6.81,param_min = 0.,dt;
65: PetscReal ftime;
67: PetscInitialize(&argc,&argv,PETSC_NULL,help);
68: user.mx = 4;
69: user.my = 4;
70: user.param = 6.0;
71:
72: /*
73: Allow user to set the grid dimensions and nonlinearity parameter at run-time
74: */
75: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
76: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
77: PetscOptionsGetReal(PETSC_NULL,"-param",&user.param,PETSC_NULL);
78: if (user.param >= param_max || user.param <= param_min) {
79: SETERRQ(1,"Parameter is out of range");
80: }
81: dt = .5/PetscMax(user.mx,user.my);
82: PetscOptionsGetReal(PETSC_NULL,"-dt",&dt,PETSC_NULL);
83: N = user.mx*user.my;
84:
85: /*
86: Create vectors to hold the solution and function value
87: */
88: VecCreateSeq(PETSC_COMM_SELF,N,&x);
89: VecDuplicate(x,&r);
91: /*
92: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
93: in the sparse matrix. Note that this is not the optimal strategy; see
94: the Performance chapter of the users manual for information on
95: preallocating memory in sparse matrices.
96: */
97: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);
99: /*
100: Create timestepper context
101: */
102: TSCreate(PETSC_COMM_WORLD,&ts);
103: TSSetProblemType(ts,TS_NONLINEAR);
105: /*
106: Tell the timestepper context where to compute solutions
107: */
108: TSSetSolution(ts,x);
110: /*
111: Provide the call-back for the nonlinear function we are
112: evaluating. Thus whenever the timestepping routines need the
113: function they will call this routine. Note the final argument
114: is the application context used by the call-back functions.
115: */
116: TSSetRHSFunction(ts,FormFunction,&user);
118: /*
119: Set the Jacobian matrix and the function used to compute
120: Jacobians.
121: */
122: TSSetRHSJacobian(ts,J,J,FormJacobian,&user);
124: /*
125: For the initial guess for the problem
126: */
127: FormInitialGuess(x,&user);
129: /*
130: This indicates that we are using pseudo timestepping to
131: find a steady state solution to the nonlinear problem.
132: */
133: TSSetType(ts,TS_PSEUDO);
135: /*
136: Set the initial time to start at (this is arbitrary for
137: steady state problems; and the initial timestep given above
138: */
139: TSSetInitialTimeStep(ts,0.0,dt);
141: /*
142: Set a large number of timesteps and final duration time
143: to insure convergence to steady state.
144: */
145: TSSetDuration(ts,1000,1.e12);
147: /*
148: Use the default strategy for increasing the timestep
149: */
150: TSPseudoSetTimeStep(ts,TSPseudoDefaultTimeStep,0);
152: /*
153: Set any additional options from the options database. This
154: includes all options for the nonlinear and linear solvers used
155: internally the the timestepping routines.
156: */
157: TSSetFromOptions(ts);
159: TSSetUp(ts);
161: /*
162: Perform the solve. This is where the timestepping takes place.
163: */
164: TSStep(ts,&its,&ftime);
165:
166: printf("Number of pseudo timesteps = %d final time %4.2en",its,ftime);
168: /*
169: Free the data structures constructed above
170: */
171: VecDestroy(x);
172: VecDestroy(r);
173: MatDestroy(J);
174: TSDestroy(ts);
175: PetscFinalize();
177: return 0;
178: }
179: /* ------------------------------------------------------------------ */
180: /* Bratu (Solid Fuel Ignition) Test Problem */
181: /* ------------------------------------------------------------------ */
183: /* -------------------- Form initial approximation ----------------- */
185: #undef __FUNCT__
187: int FormInitialGuess(Vec X,AppCtx *user)
188: {
189: int i,j,row,mx,my,ierr;
190: PetscReal one = 1.0,lambda;
191: PetscReal temp1,temp,hx,hy;
192: PetscScalar *x;
194: mx = user->mx;
195: my = user->my;
196: lambda = user->param;
198: hx = one / (PetscReal)(mx-1);
199: hy = one / (PetscReal)(my-1);
201: VecGetArray(X,&x);
202: temp1 = lambda/(lambda + one);
203: for (j=0; j<my; j++) {
204: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
205: for (i=0; i<mx; i++) {
206: row = i + j*mx;
207: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
208: x[row] = 0.0;
209: continue;
210: }
211: x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
212: }
213: }
214: VecRestoreArray(X,&x);
215: return 0;
216: }
217: /* -------------------- Evaluate Function F(x) --------------------- */
219: #undef __FUNCT__
221: int FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
222: {
223: AppCtx *user = (AppCtx*)ptr;
224: int ierr,i,j,row,mx,my;
225: PetscReal two = 2.0,one = 1.0,lambda;
226: PetscReal hx,hy,hxdhy,hydhx;
227: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*x,*f;
229: mx = user->mx;
230: my = user->my;
231: lambda = user->param;
233: hx = one / (PetscReal)(mx-1);
234: hy = one / (PetscReal)(my-1);
235: sc = hx*hy;
236: hxdhy = hx/hy;
237: hydhx = hy/hx;
239: VecGetArray(X,&x);
240: VecGetArray(F,&f);
241: for (j=0; j<my; j++) {
242: for (i=0; i<mx; i++) {
243: row = i + j*mx;
244: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
245: f[row] = x[row];
246: continue;
247: }
248: u = x[row];
249: ub = x[row - mx];
250: ul = x[row - 1];
251: ut = x[row + mx];
252: ur = x[row + 1];
253: uxx = (-ur + two*u - ul)*hydhx;
254: uyy = (-ut + two*u - ub)*hxdhy;
255: f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
256: }
257: }
258: VecRestoreArray(X,&x);
259: VecRestoreArray(F,&f);
260: return 0;
261: }
262: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
264: #undef __FUNCT__
266: int FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
267: {
268: AppCtx *user = (AppCtx*)ptr;
269: Mat jac = *B;
270: int i,j,row,mx,my,col[5],ierr;
271: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc,*x;
272: PetscReal hx,hy,hxdhy,hydhx;
275: mx = user->mx;
276: my = user->my;
277: lambda = user->param;
279: hx = 1.0 / (PetscReal)(mx-1);
280: hy = 1.0 / (PetscReal)(my-1);
281: sc = hx*hy;
282: hxdhy = hx/hy;
283: hydhx = hy/hx;
285: VecGetArray(X,&x);
286: for (j=0; j<my; j++) {
287: for (i=0; i<mx; i++) {
288: row = i + j*mx;
289: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
290: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
291: continue;
292: }
293: v[0] = hxdhy; col[0] = row - mx;
294: v[1] = hydhx; col[1] = row - 1;
295: v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
296: v[3] = hydhx; col[3] = row + 1;
297: v[4] = hxdhy; col[4] = row + mx;
298: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
299: }
300: }
301: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
302: VecRestoreArray(X,&x);
303: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
304: *flag = SAME_NONZERO_PATTERN;
305: return 0;
306: }