Actual source code: ex1.c
2: static char help[] ="Solves the time independent Bratu problem using pseudo-timestepping.";
4: /* ------------------------------------------------------------------------
6: This code demonstrates how one may solve a nonlinear problem
7: with pseudo-timestepping. In this simple example, the pseudo-timestep
8: is the same for all grid points, i.e., this is equivalent to using
9: the backward Euler method with a variable timestep.
11: Note: This example does not require pseudo-timestepping since it
12: is an easy nonlinear problem, but it is included to demonstrate how
13: the pseudo-timestepping may be done.
15: See snes/tutorials/ex4.c[ex4f.F] and
16: snes/tutorials/ex5.c[ex5f.F] where the problem is described
17: and solved using Newton's method alone.
19: ----------------------------------------------------------------------------- */
20: /*
21: Include "petscts.h" to use the PETSc timestepping routines. Note that
22: this file automatically includes "petscsys.h" and other lower-level
23: PETSc include files.
24: */
25: #include <petscts.h>
27: /*
28: Create an application context to contain data needed by the
29: application-provided call-back routines, FormJacobian() and
30: FormFunction().
31: */
32: typedef struct {
33: PetscReal param; /* test problem parameter */
34: PetscInt mx; /* Discretization in x-direction */
35: PetscInt my; /* Discretization in y-direction */
36: } AppCtx;
38: /*
39: User-defined routines
40: */
41: extern PetscErrorCode FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*), FormFunction(TS,PetscReal,Vec,Vec,void*), FormInitialGuess(Vec,AppCtx*);
43: int main(int argc,char **argv)
44: {
45: TS ts; /* timestepping context */
46: Vec x,r; /* solution, residual vectors */
47: Mat J; /* Jacobian matrix */
48: AppCtx user; /* user-defined work context */
49: PetscInt its,N; /* iterations for convergence */
50: PetscReal param_max = 6.81,param_min = 0.,dt;
51: PetscReal ftime;
52: PetscMPIInt size;
54: PetscInitialize(&argc,&argv,NULL,help);
55: MPI_Comm_size(PETSC_COMM_WORLD,&size);
58: user.mx = 4;
59: user.my = 4;
60: user.param = 6.0;
62: /*
63: Allow user to set the grid dimensions and nonlinearity parameter at run-time
64: */
65: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
66: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
67: N = user.mx*user.my;
68: dt = .5/PetscMax(user.mx,user.my);
69: PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL);
72: /*
73: Create vectors to hold the solution and function value
74: */
75: VecCreateSeq(PETSC_COMM_SELF,N,&x);
76: VecDuplicate(x,&r);
78: /*
79: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
80: in the sparse matrix. Note that this is not the optimal strategy; see
81: the Performance chapter of the users manual for information on
82: preallocating memory in sparse matrices.
83: */
84: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);
86: /*
87: Create timestepper context
88: */
89: TSCreate(PETSC_COMM_WORLD,&ts);
90: TSSetProblemType(ts,TS_NONLINEAR);
92: /*
93: Tell the timestepper context where to compute solutions
94: */
95: TSSetSolution(ts,x);
97: /*
98: Provide the call-back for the nonlinear function we are
99: evaluating. Thus whenever the timestepping routines need the
100: function they will call this routine. Note the final argument
101: is the application context used by the call-back functions.
102: */
103: TSSetRHSFunction(ts,NULL,FormFunction,&user);
105: /*
106: Set the Jacobian matrix and the function used to compute
107: Jacobians.
108: */
109: TSSetRHSJacobian(ts,J,J,FormJacobian,&user);
111: /*
112: Form the initial guess for the problem
113: */
114: FormInitialGuess(x,&user);
116: /*
117: This indicates that we are using pseudo timestepping to
118: find a steady state solution to the nonlinear problem.
119: */
120: TSSetType(ts,TSPSEUDO);
122: /*
123: Set the initial time to start at (this is arbitrary for
124: steady state problems); and the initial timestep given above
125: */
126: TSSetTimeStep(ts,dt);
128: /*
129: Set a large number of timesteps and final duration time
130: to insure convergence to steady state.
131: */
132: TSSetMaxSteps(ts,10000);
133: TSSetMaxTime(ts,1e12);
134: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
136: /*
137: Use the default strategy for increasing the timestep
138: */
139: TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0);
141: /*
142: Set any additional options from the options database. This
143: includes all options for the nonlinear and linear solvers used
144: internally the timestepping routines.
145: */
146: TSSetFromOptions(ts);
148: TSSetUp(ts);
150: /*
151: Perform the solve. This is where the timestepping takes place.
152: */
153: TSSolve(ts,x);
154: TSGetSolveTime(ts,&ftime);
156: /*
157: Get the number of steps
158: */
159: TSGetStepNumber(ts,&its);
161: PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime);
163: /*
164: Free the data structures constructed above
165: */
166: VecDestroy(&x);
167: VecDestroy(&r);
168: MatDestroy(&J);
169: TSDestroy(&ts);
170: PetscFinalize();
171: return 0;
172: }
173: /* ------------------------------------------------------------------ */
174: /* Bratu (Solid Fuel Ignition) Test Problem */
175: /* ------------------------------------------------------------------ */
177: /* -------------------- Form initial approximation ----------------- */
179: PetscErrorCode FormInitialGuess(Vec X,AppCtx *user)
180: {
181: PetscInt i,j,row,mx,my;
182: PetscReal one = 1.0,lambda;
183: PetscReal temp1,temp,hx,hy;
184: PetscScalar *x;
186: mx = user->mx;
187: my = user->my;
188: lambda = user->param;
190: hx = one / (PetscReal)(mx-1);
191: hy = one / (PetscReal)(my-1);
193: VecGetArray(X,&x);
194: temp1 = lambda/(lambda + one);
195: for (j=0; j<my; j++) {
196: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
197: for (i=0; i<mx; i++) {
198: row = i + j*mx;
199: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
200: x[row] = 0.0;
201: continue;
202: }
203: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
204: }
205: }
206: VecRestoreArray(X,&x);
207: return 0;
208: }
209: /* -------------------- Evaluate Function F(x) --------------------- */
211: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
212: {
213: AppCtx *user = (AppCtx*)ptr;
214: PetscInt i,j,row,mx,my;
215: PetscReal two = 2.0,one = 1.0,lambda;
216: PetscReal hx,hy,hxdhy,hydhx;
217: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
218: const PetscScalar *x;
220: mx = user->mx;
221: my = user->my;
222: lambda = user->param;
224: hx = one / (PetscReal)(mx-1);
225: hy = one / (PetscReal)(my-1);
226: sc = hx*hy;
227: hxdhy = hx/hy;
228: hydhx = hy/hx;
230: VecGetArrayRead(X,&x);
231: VecGetArray(F,&f);
232: for (j=0; j<my; j++) {
233: for (i=0; i<mx; i++) {
234: row = i + j*mx;
235: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
236: f[row] = x[row];
237: continue;
238: }
239: u = x[row];
240: ub = x[row - mx];
241: ul = x[row - 1];
242: ut = x[row + mx];
243: ur = x[row + 1];
244: uxx = (-ur + two*u - ul)*hydhx;
245: uyy = (-ut + two*u - ub)*hxdhy;
246: f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
247: }
248: }
249: VecRestoreArrayRead(X,&x);
250: VecRestoreArray(F,&f);
251: return 0;
252: }
253: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
255: /*
256: Calculate the Jacobian matrix J(X,t).
258: Note: We put the Jacobian in the preconditioner storage B instead of J. This
259: way we can give the -snes_mf_operator flag to check our work. This replaces
260: J with a finite difference approximation, using our analytic Jacobian B for
261: the preconditioner.
262: */
263: PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr)
264: {
265: AppCtx *user = (AppCtx*)ptr;
266: PetscInt i,j,row,mx,my,col[5];
267: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
268: const PetscScalar *x;
269: PetscReal hx,hy,hxdhy,hydhx;
271: mx = user->mx;
272: my = user->my;
273: lambda = user->param;
275: hx = 1.0 / (PetscReal)(mx-1);
276: hy = 1.0 / (PetscReal)(my-1);
277: sc = hx*hy;
278: hxdhy = hx/hy;
279: hydhx = hy/hx;
281: VecGetArrayRead(X,&x);
282: for (j=0; j<my; j++) {
283: for (i=0; i<mx; i++) {
284: row = i + j*mx;
285: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
286: MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES);
287: continue;
288: }
289: v[0] = hxdhy; col[0] = row - mx;
290: v[1] = hydhx; col[1] = row - 1;
291: v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
292: v[3] = hydhx; col[3] = row + 1;
293: v[4] = hxdhy; col[4] = row + mx;
294: MatSetValues(B,1,&row,5,col,v,INSERT_VALUES);
295: }
296: }
297: VecRestoreArrayRead(X,&x);
298: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
299: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
300: if (J != B) {
301: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
302: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
303: }
304: return 0;
305: }
307: /*TEST
309: test:
310: 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
312: test:
313: suffix: 2
314: args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5
316: TEST*/