Actual source code: ex1.c
petsc-3.7.5 2017-01-01
2: static char help[] ="Solves the time independent Bratu problem using pseudo-timestepping.";
4: /*
5: Concepts: TS^pseudo-timestepping
6: Concepts: 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/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 "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*);
53: int main(int argc,char **argv)
54: {
55: TS ts; /* timestepping context */
56: Vec x,r; /* solution, residual vectors */
57: Mat J; /* Jacobian matrix */
58: AppCtx user; /* user-defined work context */
59: PetscInt its,N; /* iterations for convergence */
61: PetscReal param_max = 6.81,param_min = 0.,dt;
62: PetscReal ftime;
63: PetscMPIInt size;
65: PetscInitialize(&argc,&argv,NULL,help);
66: MPI_Comm_size(PETSC_COMM_WORLD,&size);
67: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only");
69: user.mx = 4;
70: user.my = 4;
71: user.param = 6.0;
73: /*
74: Allow user to set the grid dimensions and nonlinearity parameter at run-time
75: */
76: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
77: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
78: N = user.mx*user.my;
79: dt = .5/PetscMax(user.mx,user.my);
80: PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL);
81: if (user.param >= param_max || user.param <= param_min) SETERRQ(PETSC_COMM_SELF,1,"Parameter is out of range");
83: /*
84: Create vectors to hold the solution and function value
85: */
86: VecCreateSeq(PETSC_COMM_SELF,N,&x);
87: VecDuplicate(x,&r);
89: /*
90: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
91: in the sparse matrix. Note that this is not the optimal strategy; see
92: the Performance chapter of the users manual for information on
93: preallocating memory in sparse matrices.
94: */
95: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);
97: /*
98: Create timestepper context
99: */
100: TSCreate(PETSC_COMM_WORLD,&ts);
101: TSSetProblemType(ts,TS_NONLINEAR);
103: /*
104: Tell the timestepper context where to compute solutions
105: */
106: TSSetSolution(ts,x);
108: /*
109: Provide the call-back for the nonlinear function we are
110: evaluating. Thus whenever the timestepping routines need the
111: function they will call this routine. Note the final argument
112: is the application context used by the call-back functions.
113: */
114: TSSetRHSFunction(ts,NULL,FormFunction,&user);
116: /*
117: Set the Jacobian matrix and the function used to compute
118: Jacobians.
119: */
120: TSSetRHSJacobian(ts,J,J,FormJacobian,&user);
122: /*
123: Form the initial guess for the problem
124: */
125: FormInitialGuess(x,&user);
127: /*
128: This indicates that we are using pseudo timestepping to
129: find a steady state solution to the nonlinear problem.
130: */
131: TSSetType(ts,TSPSEUDO);
133: /*
134: Set the initial time to start at (this is arbitrary for
135: steady state problems); and the initial timestep given above
136: */
137: TSSetInitialTimeStep(ts,0.0,dt);
139: /*
140: Set a large number of timesteps and final duration time
141: to insure convergence to steady state.
142: */
143: TSSetDuration(ts,1000,1.e12);
144: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
146: /*
147: Use the default strategy for increasing the timestep
148: */
149: TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,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 timestepping routines.
155: */
156: TSSetFromOptions(ts);
158: TSSetUp(ts);
160: /*
161: Perform the solve. This is where the timestepping takes place.
162: */
163: TSSolve(ts,x);
164: TSGetSolveTime(ts,&ftime);
166: /*
167: Get the number of steps
168: */
169: TSGetTimeStepNumber(ts,&its);
171: PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime);
173: /*
174: Free the data structures constructed above
175: */
176: VecDestroy(&x);
177: VecDestroy(&r);
178: MatDestroy(&J);
179: TSDestroy(&ts);
180: PetscFinalize();
182: return 0;
183: }
184: /* ------------------------------------------------------------------ */
185: /* Bratu (Solid Fuel Ignition) Test Problem */
186: /* ------------------------------------------------------------------ */
188: /* -------------------- Form initial approximation ----------------- */
192: PetscErrorCode FormInitialGuess(Vec X,AppCtx *user)
193: {
194: PetscInt i,j,row,mx,my;
196: PetscReal one = 1.0,lambda;
197: PetscReal temp1,temp,hx,hy;
198: PetscScalar *x;
200: mx = user->mx;
201: my = user->my;
202: lambda = user->param;
204: hx = one / (PetscReal)(mx-1);
205: hy = one / (PetscReal)(my-1);
207: VecGetArray(X,&x);
208: temp1 = lambda/(lambda + one);
209: for (j=0; j<my; j++) {
210: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
211: for (i=0; i<mx; i++) {
212: row = i + j*mx;
213: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
214: x[row] = 0.0;
215: continue;
216: }
217: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
218: }
219: }
220: VecRestoreArray(X,&x);
221: return 0;
222: }
223: /* -------------------- Evaluate Function F(x) --------------------- */
227: PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
228: {
229: AppCtx *user = (AppCtx*)ptr;
230: PetscErrorCode ierr;
231: PetscInt i,j,row,mx,my;
232: PetscReal two = 2.0,one = 1.0,lambda;
233: PetscReal hx,hy,hxdhy,hydhx;
234: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
235: const PetscScalar *x;
237: mx = user->mx;
238: my = user->my;
239: lambda = user->param;
241: hx = one / (PetscReal)(mx-1);
242: hy = one / (PetscReal)(my-1);
243: sc = hx*hy;
244: hxdhy = hx/hy;
245: hydhx = hy/hx;
247: VecGetArrayRead(X,&x);
248: VecGetArray(F,&f);
249: for (j=0; j<my; j++) {
250: for (i=0; i<mx; i++) {
251: row = i + j*mx;
252: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
253: f[row] = x[row];
254: continue;
255: }
256: u = x[row];
257: ub = x[row - mx];
258: ul = x[row - 1];
259: ut = x[row + mx];
260: ur = x[row + 1];
261: uxx = (-ur + two*u - ul)*hydhx;
262: uyy = (-ut + two*u - ub)*hxdhy;
263: f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u);
264: }
265: }
266: VecRestoreArrayRead(X,&x);
267: VecRestoreArray(F,&f);
268: return 0;
269: }
270: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
274: /*
275: Calculate the Jacobian matrix J(X,t).
277: Note: We put the Jacobian in the preconditioner storage B instead of J. This
278: way we can give the -snes_mf_operator flag to check our work. This replaces
279: J with a finite difference approximation, using our analytic Jacobian B for
280: the preconditioner.
281: */
282: PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr)
283: {
284: AppCtx *user = (AppCtx*)ptr;
285: PetscInt i,j,row,mx,my,col[5];
286: PetscErrorCode ierr;
287: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
288: const PetscScalar *x;
289: PetscReal hx,hy,hxdhy,hydhx;
292: mx = user->mx;
293: my = user->my;
294: lambda = user->param;
296: hx = 1.0 / (PetscReal)(mx-1);
297: hy = 1.0 / (PetscReal)(my-1);
298: sc = hx*hy;
299: hxdhy = hx/hy;
300: hydhx = hy/hx;
302: VecGetArrayRead(X,&x);
303: for (j=0; j<my; j++) {
304: for (i=0; i<mx; i++) {
305: row = i + j*mx;
306: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
307: MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES);
308: continue;
309: }
310: v[0] = hxdhy; col[0] = row - mx;
311: v[1] = hydhx; col[1] = row - 1;
312: v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row;
313: v[3] = hydhx; col[3] = row + 1;
314: v[4] = hxdhy; col[4] = row + mx;
315: MatSetValues(B,1,&row,5,col,v,INSERT_VALUES);
316: }
317: }
318: VecRestoreArrayRead(X,&x);
319: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
320: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
321: if (J != B) {
322: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
324: }
325: return 0;
326: }