Actual source code: ex1.c
2: static char help[] ="Solves the time independent Bratu problem using pseudo-timestepping.";
4: /*
5: Concepts: TS^pseudo-timestepping
6: Concepts: TS^pseudo-timestepping
7: Concepts: TS^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/tutorials/ex4.c[ex4f.F] and
24: snes/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 "petscsys.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: */
49: extern PetscErrorCode FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*), FormFunction(TS,PetscReal,Vec,Vec,void*), FormInitialGuess(Vec,AppCtx*);
51: int main(int argc,char **argv)
52: {
53: TS ts; /* timestepping context */
54: Vec x,r; /* solution, residual vectors */
55: Mat J; /* Jacobian matrix */
56: AppCtx user; /* user-defined work context */
57: PetscInt its,N; /* iterations for convergence */
58: PetscReal param_max = 6.81,param_min = 0.,dt;
59: PetscReal ftime;
60: PetscMPIInt size;
62: PetscInitialize(&argc,&argv,NULL,help);
63: MPI_Comm_size(PETSC_COMM_WORLD,&size);
66: user.mx = 4;
67: user.my = 4;
68: user.param = 6.0;
70: /*
71: Allow user to set the grid dimensions and nonlinearity parameter at run-time
72: */
73: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
74: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
75: N = user.mx*user.my;
76: dt = .5/PetscMax(user.mx,user.my);
77: PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL);
80: /*
81: Create vectors to hold the solution and function value
82: */
83: VecCreateSeq(PETSC_COMM_SELF,N,&x);
84: VecDuplicate(x,&r);
86: /*
87: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
88: in the sparse matrix. Note that this is not the optimal strategy; see
89: the Performance chapter of the users manual for information on
90: preallocating memory in sparse matrices.
91: */
92: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);
94: /*
95: Create timestepper context
96: */
97: TSCreate(PETSC_COMM_WORLD,&ts);
98: TSSetProblemType(ts,TS_NONLINEAR);
100: /*
101: Tell the timestepper context where to compute solutions
102: */
103: TSSetSolution(ts,x);
105: /*
106: Provide the call-back for the nonlinear function we are
107: evaluating. Thus whenever the timestepping routines need the
108: function they will call this routine. Note the final argument
109: is the application context used by the call-back functions.
110: */
111: TSSetRHSFunction(ts,NULL,FormFunction,&user);
113: /*
114: Set the Jacobian matrix and the function used to compute
115: Jacobians.
116: */
117: TSSetRHSJacobian(ts,J,J,FormJacobian,&user);
119: /*
120: Form the initial guess for the problem
121: */
122: FormInitialGuess(x,&user);
124: /*
125: This indicates that we are using pseudo timestepping to
126: find a steady state solution to the nonlinear problem.
127: */
128: TSSetType(ts,TSPSEUDO);
130: /*
131: Set the initial time to start at (this is arbitrary for
132: steady state problems); and the initial timestep given above
133: */
134: TSSetTimeStep(ts,dt);
136: /*
137: Set a large number of timesteps and final duration time
138: to insure convergence to steady state.
139: */
140: TSSetMaxSteps(ts,10000);
141: TSSetMaxTime(ts,1e12);
142: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
144: /*
145: Use the default strategy for increasing the timestep
146: */
147: TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0);
149: /*
150: Set any additional options from the options database. This
151: includes all options for the nonlinear and linear solvers used
152: internally the timestepping routines.
153: */
154: TSSetFromOptions(ts);
156: TSSetUp(ts);
158: /*
159: Perform the solve. This is where the timestepping takes place.
160: */
161: TSSolve(ts,x);
162: TSGetSolveTime(ts,&ftime);
164: /*
165: Get the number of steps
166: */
167: TSGetStepNumber(ts,&its);
169: PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime);
171: /*
172: Free the data structures constructed above
173: */
174: VecDestroy(&x);
175: VecDestroy(&r);
176: MatDestroy(&J);
177: TSDestroy(&ts);
178: PetscFinalize();
179: return 0;
180: }
181: /* ------------------------------------------------------------------ */
182: /* Bratu (Solid Fuel Ignition) Test Problem */
183: /* ------------------------------------------------------------------ */
185: /* -------------------- Form initial approximation ----------------- */
187: PetscErrorCode FormInitialGuess(Vec X,AppCtx *user)
188: {
189: 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*PetscSqrtReal(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: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
220: {
221: AppCtx *user = (AppCtx*)ptr;
222: PetscInt i,j,row,mx,my;
223: PetscReal two = 2.0,one = 1.0,lambda;
224: PetscReal hx,hy,hxdhy,hydhx;
225: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
226: const PetscScalar *x;
228: mx = user->mx;
229: my = user->my;
230: lambda = user->param;
232: hx = one / (PetscReal)(mx-1);
233: hy = one / (PetscReal)(my-1);
234: sc = hx*hy;
235: hxdhy = hx/hy;
236: hydhx = hy/hx;
238: VecGetArrayRead(X,&x);
239: VecGetArray(F,&f);
240: for (j=0; j<my; j++) {
241: for (i=0; i<mx; i++) {
242: row = i + j*mx;
243: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
244: f[row] = x[row];
245: continue;
246: }
247: u = x[row];
248: ub = x[row - mx];
249: ul = x[row - 1];
250: ut = x[row + mx];
251: ur = x[row + 1];
252: uxx = (-ur + two*u - ul)*hydhx;
253: uyy = (-ut + two*u - ub)*hxdhy;
254: f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
255: }
256: }
257: VecRestoreArrayRead(X,&x);
258: VecRestoreArray(F,&f);
259: return 0;
260: }
261: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
263: /*
264: Calculate the Jacobian matrix J(X,t).
266: Note: We put the Jacobian in the preconditioner storage B instead of J. This
267: way we can give the -snes_mf_operator flag to check our work. This replaces
268: J with a finite difference approximation, using our analytic Jacobian B for
269: the preconditioner.
270: */
271: PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr)
272: {
273: AppCtx *user = (AppCtx*)ptr;
274: PetscInt i,j,row,mx,my,col[5];
275: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
276: const PetscScalar *x;
277: PetscReal hx,hy,hxdhy,hydhx;
279: mx = user->mx;
280: my = user->my;
281: lambda = user->param;
283: hx = 1.0 / (PetscReal)(mx-1);
284: hy = 1.0 / (PetscReal)(my-1);
285: sc = hx*hy;
286: hxdhy = hx/hy;
287: hydhx = hy/hx;
289: VecGetArrayRead(X,&x);
290: for (j=0; j<my; j++) {
291: for (i=0; i<mx; i++) {
292: row = i + j*mx;
293: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
294: MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES);
295: continue;
296: }
297: v[0] = hxdhy; col[0] = row - mx;
298: v[1] = hydhx; col[1] = row - 1;
299: v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
300: v[3] = hydhx; col[3] = row + 1;
301: v[4] = hxdhy; col[4] = row + mx;
302: MatSetValues(B,1,&row,5,col,v,INSERT_VALUES);
303: }
304: }
305: VecRestoreArrayRead(X,&x);
306: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
307: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
308: if (J != B) {
309: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
310: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
311: }
312: return 0;
313: }
315: /*TEST
317: test:
318: args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo -snes_atol 1.e-7 -ts_pseudo_frtol 1.e-5 -ts_view draw:tikz:fig.tex
320: test:
321: suffix: 2
322: args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5
324: TEST*/