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 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"); 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 PetscFunctionBeginUser; 188 mx = user->mx; 189 my = user->my; 190 lambda = user->param; 191 192 hx = one / (PetscReal)(mx - 1); 193 hy = one / (PetscReal)(my - 1); 194 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) --------------------- */ 212 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; 221 222 PetscFunctionBeginUser; 223 mx = user->mx; 224 my = user->my; 225 lambda = user->param; 226 227 hx = one / (PetscReal)(mx - 1); 228 hy = one / (PetscReal)(my - 1); 229 sc = hx * hy; 230 hxdhy = hx / hy; 231 hydhx = hy / hx; 232 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) -------------------- */ 257 258 /* 259 Calculate the Jacobian matrix J(X,t). 260 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; 273 274 PetscFunctionBeginUser; 275 mx = user->mx; 276 my = user->my; 277 lambda = user->param; 278 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; 284 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 } 315 316 /*TEST 317 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 320 321 test: 322 suffix: 2 323 args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5 324 325 TEST*/ 326