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;
55: PetscInitialize(&argc, &argv, NULL, help);
56: MPI_Comm_size(PETSC_COMM_WORLD, &size);
59: user.mx = 4;
60: user.my = 4;
61: user.param = 6.0;
63: /*
64: Allow user to set the grid dimensions and nonlinearity parameter at run-time
65: */
66: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL);
67: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL);
68: N = user.mx * user.my;
69: dt = .5 / PetscMax(user.mx, user.my);
70: PetscOptionsGetReal(NULL, NULL, "-param", &user.param, NULL);
73: /*
74: Create vectors to hold the solution and function value
75: */
76: VecCreateSeq(PETSC_COMM_SELF, N, &x);
77: VecDuplicate(x, &r);
79: /*
80: Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
81: in the sparse matrix. Note that this is not the optimal strategy; see
82: the Performance chapter of the users manual for information on
83: preallocating memory in sparse matrices.
84: */
85: MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5, 0, &J);
87: /*
88: Create timestepper context
89: */
90: TSCreate(PETSC_COMM_WORLD, &ts);
91: TSSetProblemType(ts, TS_NONLINEAR);
93: /*
94: Tell the timestepper context where to compute solutions
95: */
96: TSSetSolution(ts, x);
98: /*
99: Provide the call-back for the nonlinear function we are
100: evaluating. Thus whenever the timestepping routines need the
101: function they will call this routine. Note the final argument
102: is the application context used by the call-back functions.
103: */
104: TSSetRHSFunction(ts, NULL, FormFunction, &user);
106: /*
107: Set the Jacobian matrix and the function used to compute
108: Jacobians.
109: */
110: TSSetRHSJacobian(ts, J, J, FormJacobian, &user);
112: /*
113: Form the initial guess for the problem
114: */
115: FormInitialGuess(x, &user);
117: /*
118: This indicates that we are using pseudo timestepping to
119: find a steady state solution to the nonlinear problem.
120: */
121: TSSetType(ts, TSPSEUDO);
123: /*
124: Set the initial time to start at (this is arbitrary for
125: steady state problems); and the initial timestep given above
126: */
127: TSSetTimeStep(ts, dt);
129: /*
130: Set a large number of timesteps and final duration time
131: to insure convergence to steady state.
132: */
133: TSSetMaxSteps(ts, 10000);
134: TSSetMaxTime(ts, 1e12);
135: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
137: /*
138: Use the default strategy for increasing the timestep
139: */
140: TSPseudoSetTimeStep(ts, TSPseudoTimeStepDefault, 0);
142: /*
143: Set any additional options from the options database. This
144: includes all options for the nonlinear and linear solvers used
145: internally the timestepping routines.
146: */
147: TSSetFromOptions(ts);
149: TSSetUp(ts);
151: /*
152: Perform the solve. This is where the timestepping takes place.
153: */
154: TSSolve(ts, x);
155: TSGetSolveTime(ts, &ftime);
157: /*
158: Get the number of steps
159: */
160: TSGetStepNumber(ts, &its);
162: PetscPrintf(PETSC_COMM_WORLD, "Number of pseudo timesteps = %" PetscInt_FMT " final time %4.2e\n", its, (double)ftime);
164: /*
165: Free the data structures constructed above
166: */
167: VecDestroy(&x);
168: VecDestroy(&r);
169: MatDestroy(&J);
170: TSDestroy(&ts);
171: PetscFinalize();
172: return 0;
173: }
174: /* ------------------------------------------------------------------ */
175: /* Bratu (Solid Fuel Ignition) Test Problem */
176: /* ------------------------------------------------------------------ */
178: /* -------------------- Form initial approximation ----------------- */
180: PetscErrorCode FormInitialGuess(Vec X, AppCtx *user)
181: {
182: PetscInt i, j, row, mx, my;
183: PetscReal one = 1.0, lambda;
184: PetscReal temp1, temp, hx, hy;
185: PetscScalar *x;
187: mx = user->mx;
188: my = user->my;
189: lambda = user->param;
191: hx = one / (PetscReal)(mx - 1);
192: hy = one / (PetscReal)(my - 1);
194: VecGetArray(X, &x);
195: temp1 = lambda / (lambda + one);
196: for (j = 0; j < my; j++) {
197: temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
198: for (i = 0; i < mx; i++) {
199: row = i + j * mx;
200: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
201: x[row] = 0.0;
202: continue;
203: }
204: x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
205: }
206: }
207: VecRestoreArray(X, &x);
208: return 0;
209: }
210: /* -------------------- Evaluate Function F(x) --------------------- */
212: PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
213: {
214: AppCtx *user = (AppCtx *)ptr;
215: PetscInt i, j, row, mx, my;
216: PetscReal two = 2.0, one = 1.0, lambda;
217: PetscReal hx, hy, hxdhy, hydhx;
218: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f;
219: const PetscScalar *x;
221: mx = user->mx;
222: my = user->my;
223: lambda = user->param;
225: hx = one / (PetscReal)(mx - 1);
226: hy = one / (PetscReal)(my - 1);
227: sc = hx * hy;
228: hxdhy = hx / hy;
229: hydhx = hy / hx;
231: VecGetArrayRead(X, &x);
232: VecGetArray(F, &f);
233: for (j = 0; j < my; j++) {
234: for (i = 0; i < mx; i++) {
235: row = i + j * mx;
236: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
237: f[row] = x[row];
238: continue;
239: }
240: u = x[row];
241: ub = x[row - mx];
242: ul = x[row - 1];
243: ut = x[row + mx];
244: ur = x[row + 1];
245: uxx = (-ur + two * u - ul) * hydhx;
246: uyy = (-ut + two * u - ub) * hxdhy;
247: f[row] = -uxx + -uyy + sc * lambda * PetscExpScalar(u);
248: }
249: }
250: VecRestoreArrayRead(X, &x);
251: VecRestoreArray(F, &f);
252: return 0;
253: }
254: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
256: /*
257: Calculate the Jacobian matrix J(X,t).
259: Note: We put the Jacobian in the preconditioner storage B instead of J. This
260: way we can give the -snes_mf_operator flag to check our work. This replaces
261: J with a finite difference approximation, using our analytic Jacobian B for
262: the preconditioner.
263: */
264: PetscErrorCode FormJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
265: {
266: AppCtx *user = (AppCtx *)ptr;
267: PetscInt i, j, row, mx, my, col[5];
268: PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc;
269: const PetscScalar *x;
270: PetscReal hx, hy, hxdhy, hydhx;
272: mx = user->mx;
273: my = user->my;
274: lambda = user->param;
276: hx = 1.0 / (PetscReal)(mx - 1);
277: hy = 1.0 / (PetscReal)(my - 1);
278: sc = hx * hy;
279: hxdhy = hx / hy;
280: hydhx = hy / hx;
282: VecGetArrayRead(X, &x);
283: for (j = 0; j < my; j++) {
284: for (i = 0; i < mx; i++) {
285: row = i + j * mx;
286: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
287: MatSetValues(B, 1, &row, 1, &row, &one, INSERT_VALUES);
288: continue;
289: }
290: v[0] = hxdhy;
291: col[0] = row - mx;
292: v[1] = hydhx;
293: col[1] = row - 1;
294: v[2] = -two * (hydhx + hxdhy) + sc * lambda * PetscExpScalar(x[row]);
295: col[2] = row;
296: v[3] = hydhx;
297: col[3] = row + 1;
298: v[4] = hxdhy;
299: col[4] = row + mx;
300: MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES);
301: }
302: }
303: VecRestoreArrayRead(X, &x);
304: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
306: if (J != B) {
307: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
309: }
310: return 0;
311: }
313: /*TEST
315: test:
316: 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
318: test:
319: suffix: 2
320: args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5
322: TEST*/