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 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 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 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 */ 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