1 static char help[] = "Solves the time independent Bratu problem using pseudo-timestepping.";
2
3 /* ------------------------------------------------------------------------
4
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.
9
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.
13
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.
17
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>
25
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;
36
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 *);
41
main(int argc,char ** argv)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;
52
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");
57
58 user.mx = 4;
59 user.my = 4;
60 user.param = 6.0;
61
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");
71
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));
77
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));
85
86 /*
87 Create timestepper context
88 */
89 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
90 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
91
92 /*
93 Tell the timestepper context where to compute solutions
94 */
95 PetscCall(TSSetSolution(ts, x));
96
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));
104
105 /*
106 Set the Jacobian matrix and the function used to compute
107 Jacobians.
108 */
109 PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &user));
110
111 /*
112 Form the initial guess for the problem
113 */
114 PetscCall(FormInitialGuess(x, &user));
115
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));
121
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));
127
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));
135
136 /*
137 Use the default strategy for increasing the timestep
138 */
139 PetscCall(TSPseudoSetTimeStep(ts, TSPseudoTimeStepDefault, 0));
140
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));
147
148 PetscCall(TSSetUp(ts));
149
150 /*
151 Perform the solve. This is where the timestepping takes place.
152 */
153 PetscCall(TSSolve(ts, x));
154 PetscCall(TSGetSolveTime(ts, &ftime));
155
156 /*
157 Get the number of steps
158 */
159 PetscCall(TSGetStepNumber(ts, &its));
160
161 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of pseudo timesteps = %" PetscInt_FMT " final time %4.2e\n", its, (double)ftime));
162
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 /* ------------------------------------------------------------------ */
176
177 /* -------------------- Form initial approximation ----------------- */
178
FormInitialGuess(Vec X,AppCtx * user)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;
185
186 PetscFunctionBeginUser;
187 mx = user->mx;
188 my = user->my;
189 lambda = user->param;
190
191 hx = one / (PetscReal)(mx - 1);
192 hy = one / (PetscReal)(my - 1);
193
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) --------------------- */
211
FormFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
220
221 PetscFunctionBeginUser;
222 mx = user->mx;
223 my = user->my;
224 lambda = user->param;
225
226 hx = one / (PetscReal)(mx - 1);
227 hy = one / (PetscReal)(my - 1);
228 sc = hx * hy;
229 hxdhy = hx / hy;
230 hydhx = hy / hx;
231
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) -------------------- */
256
257 /*
258 Calculate the Jacobian matrix J(X,t).
259
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 */
FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void * ptr)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;
272
273 PetscFunctionBeginUser;
274 mx = user->mx;
275 my = user->my;
276 lambda = user->param;
277
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;
283
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 }
314
315 /*TEST
316
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
319
320 test:
321 suffix: 2
322 args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5
323
324 TEST*/
325