1 2 static char help[] = "Solves the time independent Bratu problem using pseudo-timestepping."; 3 4 /* ------------------------------------------------------------------------ 5 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. 10 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. 14 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. 18 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> 26 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; 37 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 *); 42 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; 53 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"); 58 59 user.mx = 4; 60 user.my = 4; 61 user.param = 6.0; 62 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); 71 PetscCheck(user.param < param_max && user.param >= param_min, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Parameter is out of range"); 72 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)); 78 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)); 86 87 /* 88 Create timestepper context 89 */ 90 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 91 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 92 93 /* 94 Tell the timestepper context where to compute solutions 95 */ 96 PetscCall(TSSetSolution(ts, x)); 97 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)); 105 106 /* 107 Set the Jacobian matrix and the function used to compute 108 Jacobians. 109 */ 110 PetscCall(TSSetRHSJacobian(ts, J, J, FormJacobian, &user)); 111 112 /* 113 Form the initial guess for the problem 114 */ 115 PetscCall(FormInitialGuess(x, &user)); 116 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)); 122 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)); 128 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)); 136 137 /* 138 Use the default strategy for increasing the timestep 139 */ 140 PetscCall(TSPseudoSetTimeStep(ts, TSPseudoTimeStepDefault, 0)); 141 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)); 148 149 PetscCall(TSSetUp(ts)); 150 151 /* 152 Perform the solve. This is where the timestepping takes place. 153 */ 154 PetscCall(TSSolve(ts, x)); 155 PetscCall(TSGetSolveTime(ts, &ftime)); 156 157 /* 158 Get the number of steps 159 */ 160 PetscCall(TSGetStepNumber(ts, &its)); 161 162 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of pseudo timesteps = %" PetscInt_FMT " final time %4.2e\n", its, (double)ftime)); 163 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 /* ------------------------------------------------------------------ */ 177 178 /* -------------------- Form initial approximation ----------------- */ 179 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; 186 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 return 0; 209 } 210 /* -------------------- Evaluate Function F(x) --------------------- */ 211 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 mx = user->mx; 222 my = user->my; 223 lambda = user->param; 224 225 hx = one / (PetscReal)(mx - 1); 226 hy = one / (PetscReal)(my - 1); 227 sc = hx * hy; 228 hxdhy = hx / hy; 229 hydhx = hy / hx; 230 231 PetscCall(VecGetArrayRead(X, &x)); 232 PetscCall(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 PetscCall(VecRestoreArrayRead(X, &x)); 251 PetscCall(VecRestoreArray(F, &f)); 252 return 0; 253 } 254 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 255 256 /* 257 Calculate the Jacobian matrix J(X,t). 258 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; 271 272 mx = user->mx; 273 my = user->my; 274 lambda = user->param; 275 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; 281 282 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES)); 301 } 302 } 303 PetscCall(VecRestoreArrayRead(X, &x)); 304 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 305 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 306 if (J != B) { 307 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 308 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 309 } 310 return 0; 311 } 312 313 /*TEST 314 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 317 318 test: 319 suffix: 2 320 args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5 321 322 TEST*/ 323