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: PetscFunctionBeginUser;
55: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
56: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
57: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
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: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
67: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
68: N = user.mx * user.my;
69: dt = .5 / PetscMax(user.mx, user.my);
70: PetscCall(PetscOptionsGetReal(NULL, NULL, "-param", &user.param, NULL));
71: PetscCheck(user.param < param_max && user.param >= param_min, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Parameter is out of range");
73: /*
74: Create vectors to hold the solution and function value
75: */
76: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
77: PetscCall(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: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5, 0, &J));
87: /*
88: Create timestepper context
89: */
90: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
91: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
93: /*
94: Tell the timestepper context where to compute solutions
95: */
96: PetscCall(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: PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user));
106: /*
107: Set the Jacobian matrix and the function used to compute
108: Jacobians.
109: */
110: PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &user));
112: /*
113: Form the initial guess for the problem
114: */
115: PetscCall(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: PetscCall(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: PetscCall(TSSetTimeStep(ts, dt));
129: /*
130: Set a large number of timesteps and final duration time
131: to insure convergence to steady state.
132: */
133: PetscCall(TSSetMaxSteps(ts, 10000));
134: PetscCall(TSSetMaxTime(ts, 1e12));
135: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
137: /*
138: Use the default strategy for increasing the timestep
139: */
140: PetscCall(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: PetscCall(TSSetFromOptions(ts));
149: PetscCall(TSSetUp(ts));
151: /*
152: Perform the solve. This is where the timestepping takes place.
153: */
154: PetscCall(TSSolve(ts, x));
155: PetscCall(TSGetSolveTime(ts, &ftime));
157: /*
158: Get the number of steps
159: */
160: PetscCall(TSGetStepNumber(ts, &its));
162: PetscCall(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: PetscCall(VecDestroy(&x));
168: PetscCall(VecDestroy(&r));
169: PetscCall(MatDestroy(&J));
170: PetscCall(TSDestroy(&ts));
171: PetscCall(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: PetscFunctionBeginUser;
188: mx = user->mx;
189: my = user->my;
190: lambda = user->param;
192: hx = one / (PetscReal)(mx - 1);
193: hy = one / (PetscReal)(my - 1);
195: PetscCall(VecGetArray(X, &x));
196: temp1 = lambda / (lambda + one);
197: for (j = 0; j < my; j++) {
198: temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
199: for (i = 0; i < mx; i++) {
200: row = i + j * mx;
201: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
202: x[row] = 0.0;
203: continue;
204: }
205: x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
206: }
207: }
208: PetscCall(VecRestoreArray(X, &x));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
211: /* -------------------- Evaluate Function F(x) --------------------- */
213: PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
214: {
215: AppCtx *user = (AppCtx *)ptr;
216: PetscInt i, j, row, mx, my;
217: PetscReal two = 2.0, one = 1.0, lambda;
218: PetscReal hx, hy, hxdhy, hydhx;
219: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f;
220: const PetscScalar *x;
222: PetscFunctionBeginUser;
223: mx = user->mx;
224: my = user->my;
225: lambda = user->param;
227: hx = one / (PetscReal)(mx - 1);
228: hy = one / (PetscReal)(my - 1);
229: sc = hx * hy;
230: hxdhy = hx / hy;
231: hydhx = hy / hx;
233: PetscCall(VecGetArrayRead(X, &x));
234: PetscCall(VecGetArray(F, &f));
235: for (j = 0; j < my; j++) {
236: for (i = 0; i < mx; i++) {
237: row = i + j * mx;
238: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
239: f[row] = x[row];
240: continue;
241: }
242: u = x[row];
243: ub = x[row - mx];
244: ul = x[row - 1];
245: ut = x[row + mx];
246: ur = x[row + 1];
247: uxx = (-ur + two * u - ul) * hydhx;
248: uyy = (-ut + two * u - ub) * hxdhy;
249: f[row] = -uxx + -uyy + sc * lambda * PetscExpScalar(u);
250: }
251: }
252: PetscCall(VecRestoreArrayRead(X, &x));
253: PetscCall(VecRestoreArray(F, &f));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
256: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
258: /*
259: Calculate the Jacobian matrix J(X,t).
261: Note: We put the Jacobian in the preconditioner storage B instead of J. This
262: way we can give the -snes_mf_operator flag to check our work. This replaces
263: J with a finite difference approximation, using our analytic Jacobian B for
264: the preconditioner.
265: */
266: PetscErrorCode FormJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr)
267: {
268: AppCtx *user = (AppCtx *)ptr;
269: PetscInt i, j, row, mx, my, col[5];
270: PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc;
271: const PetscScalar *x;
272: PetscReal hx, hy, hxdhy, hydhx;
274: PetscFunctionBeginUser;
275: mx = user->mx;
276: my = user->my;
277: lambda = user->param;
279: hx = 1.0 / (PetscReal)(mx - 1);
280: hy = 1.0 / (PetscReal)(my - 1);
281: sc = hx * hy;
282: hxdhy = hx / hy;
283: hydhx = hy / hx;
285: PetscCall(VecGetArrayRead(X, &x));
286: for (j = 0; j < my; j++) {
287: for (i = 0; i < mx; i++) {
288: row = i + j * mx;
289: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
290: PetscCall(MatSetValues(B, 1, &row, 1, &row, &one, INSERT_VALUES));
291: continue;
292: }
293: v[0] = hxdhy;
294: col[0] = row - mx;
295: v[1] = hydhx;
296: col[1] = row - 1;
297: v[2] = -two * (hydhx + hxdhy) + sc * lambda * PetscExpScalar(x[row]);
298: col[2] = row;
299: v[3] = hydhx;
300: col[3] = row + 1;
301: v[4] = hxdhy;
302: col[4] = row + mx;
303: PetscCall(MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES));
304: }
305: }
306: PetscCall(VecRestoreArrayRead(X, &x));
307: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
308: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
309: if (J != B) {
310: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
311: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
312: }
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*TEST
318: test:
319: 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
321: test:
322: suffix: 2
323: args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5
325: TEST*/