Actual source code: ex1.c
1: static char help[] = "Solves the time independent Bratu problem using pseudo-timestepping.";
3: /* ------------------------------------------------------------------------
5: This code demonstrates how one may solve a nonlinear problem
6: with pseudo-timestepping. In this simple example, the pseudo-timestep
7: is the same for all grid points, i.e., this is equivalent to using
8: the backward Euler method with a variable timestep.
10: Note: This example does not require pseudo-timestepping since it
11: is an easy nonlinear problem, but it is included to demonstrate how
12: the pseudo-timestepping may be done.
14: See snes/tutorials/ex4.c[ex4f.F] and
15: snes/tutorials/ex5.c[ex5f.F] where the problem is described
16: and solved using Newton's method alone.
18: ----------------------------------------------------------------------------- */
19: /*
20: Include "petscts.h" to use the PETSc timestepping routines. Note that
21: this file automatically includes "petscsys.h" and other lower-level
22: PETSc include files.
23: */
24: #include <petscts.h>
26: /*
27: Create an application context to contain data needed by the
28: application-provided call-back routines, FormJacobian() and
29: FormFunction().
30: */
31: typedef struct {
32: PetscReal param; /* test problem parameter */
33: PetscInt mx; /* Discretization in x-direction */
34: PetscInt my; /* Discretization in y-direction */
35: } AppCtx;
37: /*
38: User-defined routines
39: */
40: extern PetscErrorCode FormJacobian(TS, PetscReal, Vec, Mat, Mat, void *), FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialGuess(Vec, AppCtx *);
42: int main(int argc, char **argv)
43: {
44: TS ts; /* timestepping context */
45: Vec x, r; /* solution, residual vectors */
46: Mat J; /* Jacobian matrix */
47: AppCtx user; /* user-defined work context */
48: PetscInt its, N; /* iterations for convergence */
49: PetscReal param_max = 6.81, param_min = 0., dt;
50: PetscReal ftime;
51: PetscMPIInt size;
53: PetscFunctionBeginUser;
54: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
55: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
56: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
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: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
66: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
67: N = user.mx * user.my;
68: dt = .5 / PetscMax(user.mx, user.my);
69: PetscCall(PetscOptionsGetReal(NULL, NULL, "-param", &user.param, NULL));
70: PetscCheck(user.param < param_max && user.param >= param_min, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Parameter is out of range");
72: /*
73: Create vectors to hold the solution and function value
74: */
75: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
76: PetscCall(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: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5, 0, &J));
86: /*
87: Create timestepper context
88: */
89: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
90: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
92: /*
93: Tell the timestepper context where to compute solutions
94: */
95: PetscCall(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: PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user));
105: /*
106: Set the Jacobian matrix and the function used to compute
107: Jacobians.
108: */
109: PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &user));
111: /*
112: Form the initial guess for the problem
113: */
114: PetscCall(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: PetscCall(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: PetscCall(TSSetTimeStep(ts, dt));
128: /*
129: Set a large number of timesteps and final duration time
130: to insure convergence to steady state.
131: */
132: PetscCall(TSSetMaxSteps(ts, 10000));
133: PetscCall(TSSetMaxTime(ts, 1e12));
134: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
136: /*
137: Use the default strategy for increasing the timestep
138: */
139: PetscCall(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: PetscCall(TSSetFromOptions(ts));
148: PetscCall(TSSetUp(ts));
150: /*
151: Perform the solve. This is where the timestepping takes place.
152: */
153: PetscCall(TSSolve(ts, x));
154: PetscCall(TSGetSolveTime(ts, &ftime));
156: /*
157: Get the number of steps
158: */
159: PetscCall(TSGetStepNumber(ts, &its));
161: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of pseudo timesteps = %" PetscInt_FMT " final time %4.2e\n", its, (double)ftime));
163: /*
164: Free the data structures constructed above
165: */
166: PetscCall(VecDestroy(&x));
167: PetscCall(VecDestroy(&r));
168: PetscCall(MatDestroy(&J));
169: PetscCall(TSDestroy(&ts));
170: PetscCall(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: PetscFunctionBeginUser;
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: PetscCall(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: PetscCall(VecRestoreArray(X, &x));
208: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
222: mx = user->mx;
223: my = user->my;
224: lambda = user->param;
226: hx = one / (PetscReal)(mx - 1);
227: hy = one / (PetscReal)(my - 1);
228: sc = hx * hy;
229: hxdhy = hx / hy;
230: hydhx = hy / hx;
232: PetscCall(VecGetArrayRead(X, &x));
233: PetscCall(VecGetArray(F, &f));
234: for (j = 0; j < my; j++) {
235: for (i = 0; i < mx; i++) {
236: row = i + j * mx;
237: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
238: f[row] = x[row];
239: continue;
240: }
241: u = x[row];
242: ub = x[row - mx];
243: ul = x[row - 1];
244: ut = x[row + mx];
245: ur = x[row + 1];
246: uxx = (-ur + two * u - ul) * hydhx;
247: uyy = (-ut + two * u - ub) * hxdhy;
248: f[row] = -uxx + -uyy + sc * lambda * PetscExpScalar(u);
249: }
250: }
251: PetscCall(VecRestoreArrayRead(X, &x));
252: PetscCall(VecRestoreArray(F, &f));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
255: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
257: /*
258: Calculate the Jacobian matrix J(X,t).
260: Note: We put the Jacobian in the preconditioner storage B instead of J. This
261: way we can give the -snes_mf_operator flag to check our work. This replaces
262: J with a finite difference approximation, using our analytic Jacobian B for
263: the preconditioner.
264: */
265: PetscErrorCode FormJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
266: {
267: AppCtx *user = (AppCtx *)ptr;
268: PetscInt i, j, row, mx, my, col[5];
269: PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc;
270: const PetscScalar *x;
271: PetscReal hx, hy, hxdhy, hydhx;
273: PetscFunctionBeginUser;
274: mx = user->mx;
275: my = user->my;
276: lambda = user->param;
278: hx = 1.0 / (PetscReal)(mx - 1);
279: hy = 1.0 / (PetscReal)(my - 1);
280: sc = hx * hy;
281: hxdhy = hx / hy;
282: hydhx = hy / hx;
284: PetscCall(VecGetArrayRead(X, &x));
285: for (j = 0; j < my; j++) {
286: for (i = 0; i < mx; i++) {
287: row = i + j * mx;
288: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
289: PetscCall(MatSetValues(B, 1, &row, 1, &row, &one, INSERT_VALUES));
290: continue;
291: }
292: v[0] = hxdhy;
293: col[0] = row - mx;
294: v[1] = hydhx;
295: col[1] = row - 1;
296: v[2] = -two * (hydhx + hxdhy) + sc * lambda * PetscExpScalar(x[row]);
297: col[2] = row;
298: v[3] = hydhx;
299: col[3] = row + 1;
300: v[4] = hxdhy;
301: col[4] = row + mx;
302: PetscCall(MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES));
303: }
304: }
305: PetscCall(VecRestoreArrayRead(X, &x));
306: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
307: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
308: if (J != B) {
309: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
310: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
311: }
312: PetscFunctionReturn(PETSC_SUCCESS);
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*/