Actual source code: ex1.c
2: static char help[] ="Solves the time dependent Bratu problem using pseudo-timestepping.";
4: /*
5: Concepts: TS^pseudo-timestepping
6: Concepts: pseudo-timestepping
7: Concepts: nonlinear problems
8: Processors: 1
10: */
12: /* ------------------------------------------------------------------------
14: This code demonstrates how one may solve a nonlinear problem
15: with pseudo-timestepping. In this simple example, the pseudo-timestep
16: is the same for all grid points, i.e., this is equivalent to using
17: the backward Euler method with a variable timestep.
19: Note: This example does not require pseudo-timestepping since it
20: is an easy nonlinear problem, but it is included to demonstrate how
21: the pseudo-timestepping may be done.
23: See snes/examples/tutorials/ex4.c[ex4f.F] and
24: snes/examples/tutorials/ex5.c[ex5f.F] where the problem is described
25: and solved using Newton's method alone.
27: ----------------------------------------------------------------------------- */
28: /*
29: Include "petscts.h" to use the PETSc timestepping routines. Note that
30: this file automatically includes "petsc.h" and other lower-level
31: PETSc include files.
32: */
33: #include petscts.h
35: /*
36: Create an application context to contain data needed by the
37: application-provided call-back routines, FormJacobian() and
38: FormFunction().
39: */
40: typedef struct {
41: PetscReal param; /* test problem parameter */
42: PetscInt mx; /* Discretization in x-direction */
43: PetscInt my; /* Discretization in y-direction */
44: } AppCtx;
46: /*
47: User-defined routines
48: */
50: FormFunction(TS,PetscReal,Vec,Vec,void*),
51: FormInitialGuess(Vec,AppCtx*);
55: int main(int argc,char **argv)
56: {
57: TS ts; /* timestepping context */
58: Vec x,r; /* solution, residual vectors */
59: Mat J; /* Jacobian matrix */
60: AppCtx user; /* user-defined work context */
61: PetscInt its,N; /* iterations for convergence */
63: PetscReal param_max = 6.81,param_min = 0.,dt;
64: PetscReal ftime;
66: PetscInitialize(&argc,&argv,PETSC_NULL,help);
67: user.mx = 4;
68: user.my = 4;
69: user.param = 6.0;
70:
71: /*
72: Allow user to set the grid dimensions and nonlinearity parameter at run-time
73: */
74: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
75: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
76: PetscOptionsGetReal(PETSC_NULL,"-param",&user.param,PETSC_NULL);
77: if (user.param >= param_max || user.param <= param_min) {
78: SETERRQ(1,"Parameter is out of range");
79: }
80: dt = .5/PetscMax(user.mx,user.my);
81: PetscOptionsGetReal(PETSC_NULL,"-dt",&dt,PETSC_NULL);
82: N = user.mx*user.my;
83:
84: /*
85: Create vectors to hold the solution and function value
86: */
87: VecCreateSeq(PETSC_COMM_SELF,N,&x);
88: VecDuplicate(x,&r);
90: /*
91: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
92: in the sparse matrix. Note that this is not the optimal strategy; see
93: the Performance chapter of the users manual for information on
94: preallocating memory in sparse matrices.
95: */
96: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);
98: /*
99: Create timestepper context
100: */
101: TSCreate(PETSC_COMM_WORLD,&ts);
102: TSSetProblemType(ts,TS_NONLINEAR);
104: /*
105: Tell the timestepper context where to compute solutions
106: */
107: TSSetSolution(ts,x);
109: /*
110: Provide the call-back for the nonlinear function we are
111: evaluating. Thus whenever the timestepping routines need the
112: function they will call this routine. Note the final argument
113: is the application context used by the call-back functions.
114: */
115: TSSetRHSFunction(ts,FormFunction,&user);
117: /*
118: Set the Jacobian matrix and the function used to compute
119: Jacobians.
120: */
121: TSSetRHSJacobian(ts,J,J,FormJacobian,&user);
123: /*
124: For the initial guess for the problem
125: */
126: FormInitialGuess(x,&user);
128: /*
129: This indicates that we are using pseudo timestepping to
130: find a steady state solution to the nonlinear problem.
131: */
132: TSSetType(ts,TS_PSEUDO);
134: /*
135: Set the initial time to start at (this is arbitrary for
136: steady state problems; and the initial timestep given above
137: */
138: TSSetInitialTimeStep(ts,0.0,dt);
140: /*
141: Set a large number of timesteps and final duration time
142: to insure convergence to steady state.
143: */
144: TSSetDuration(ts,1000,1.e12);
146: /*
147: Use the default strategy for increasing the timestep
148: */
149: TSPseudoSetTimeStep(ts,TSPseudoDefaultTimeStep,0);
151: /*
152: Set any additional options from the options database. This
153: includes all options for the nonlinear and linear solvers used
154: internally the the timestepping routines.
155: */
156: TSSetFromOptions(ts);
158: TSSetUp(ts);
160: /*
161: Perform the solve. This is where the timestepping takes place.
162: */
163: TSStep(ts,&its,&ftime);
164:
165: printf("Number of pseudo timesteps = %d final time %4.2e\n",(int)its,ftime);
167: /*
168: Free the data structures constructed above
169: */
170: VecDestroy(x);
171: VecDestroy(r);
172: MatDestroy(J);
173: TSDestroy(ts);
174: PetscFinalize();
176: return 0;
177: }
178: /* ------------------------------------------------------------------ */
179: /* Bratu (Solid Fuel Ignition) Test Problem */
180: /* ------------------------------------------------------------------ */
182: /* -------------------- Form initial approximation ----------------- */
186: PetscErrorCode FormInitialGuess(Vec X,AppCtx *user)
187: {
188: PetscInt i,j,row,mx,my;
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) --------------------- */
221: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
222: {
223: AppCtx *user = (AppCtx*)ptr;
225: PetscInt i,j,row,mx,my;
226: PetscReal two = 2.0,one = 1.0,lambda;
227: PetscReal hx,hy,hxdhy,hydhx;
228: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*x,*f;
230: mx = user->mx;
231: my = user->my;
232: lambda = user->param;
234: hx = one / (PetscReal)(mx-1);
235: hy = one / (PetscReal)(my-1);
236: sc = hx*hy;
237: hxdhy = hx/hy;
238: hydhx = hy/hx;
240: VecGetArray(X,&x);
241: VecGetArray(F,&f);
242: for (j=0; j<my; j++) {
243: for (i=0; i<mx; i++) {
244: row = i + j*mx;
245: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
246: f[row] = x[row];
247: continue;
248: }
249: u = x[row];
250: ub = x[row - mx];
251: ul = x[row - 1];
252: ut = x[row + mx];
253: ur = x[row + 1];
254: uxx = (-ur + two*u - ul)*hydhx;
255: uyy = (-ut + two*u - ub)*hxdhy;
256: f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
257: }
258: }
259: VecRestoreArray(X,&x);
260: VecRestoreArray(F,&f);
261: return 0;
262: }
263: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
267: /*
268: Calculate the Jacobian matrix J(X,t).
270: Note: We put the Jacobian in the preconditioner storage B instead of J. This
271: way we can give the -snes_mf_operator flag to check our work. This replaces
272: J with a finite difference approximation, using our analytic Jacobian B for
273: the preconditioner.
274: */
275: PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
276: {
277: AppCtx *user = (AppCtx*)ptr;
278: Mat jac = *B;
279: PetscInt i,j,row,mx,my,col[5];
281: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc,*x;
282: PetscReal hx,hy,hxdhy,hydhx;
285: mx = user->mx;
286: my = user->my;
287: lambda = user->param;
289: hx = 1.0 / (PetscReal)(mx-1);
290: hy = 1.0 / (PetscReal)(my-1);
291: sc = hx*hy;
292: hxdhy = hx/hy;
293: hydhx = hy/hx;
295: VecGetArray(X,&x);
296: for (j=0; j<my; j++) {
297: for (i=0; i<mx; i++) {
298: row = i + j*mx;
299: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
300: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
301: continue;
302: }
303: v[0] = hxdhy; col[0] = row - mx;
304: v[1] = hydhx; col[1] = row - 1;
305: v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
306: v[3] = hydhx; col[3] = row + 1;
307: v[4] = hxdhy; col[4] = row + mx;
308: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
309: }
310: }
311: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
312: VecRestoreArray(X,&x);
313: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
314: *flag = SAME_NONZERO_PATTERN;
315: return 0;
316: }