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 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 PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL); 66 PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL); 67 N = user.mx * user.my; 68 dt = .5 / PetscMax(user.mx, user.my); 69 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 PetscInt i, j, row, mx, my; 181 PetscReal one = 1.0, lambda; 182 PetscReal temp1, temp, hx, hy; 183 PetscScalar *x; 184 185 mx = user->mx; 186 my = user->my; 187 lambda = user->param; 188 189 hx = one / (PetscReal)(mx - 1); 190 hy = one / (PetscReal)(my - 1); 191 192 PetscCall(VecGetArray(X, &x)); 193 temp1 = lambda / (lambda + one); 194 for (j = 0; j < my; j++) { 195 temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy; 196 for (i = 0; i < mx; i++) { 197 row = i + j * mx; 198 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 199 x[row] = 0.0; 200 continue; 201 } 202 x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp)); 203 } 204 } 205 PetscCall(VecRestoreArray(X, &x)); 206 return 0; 207 } 208 /* -------------------- Evaluate Function F(x) --------------------- */ 209 210 PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 211 AppCtx *user = (AppCtx *)ptr; 212 PetscInt i, j, row, mx, my; 213 PetscReal two = 2.0, one = 1.0, lambda; 214 PetscReal hx, hy, hxdhy, hydhx; 215 PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f; 216 const PetscScalar *x; 217 218 mx = user->mx; 219 my = user->my; 220 lambda = user->param; 221 222 hx = one / (PetscReal)(mx - 1); 223 hy = one / (PetscReal)(my - 1); 224 sc = hx * hy; 225 hxdhy = hx / hy; 226 hydhx = hy / hx; 227 228 PetscCall(VecGetArrayRead(X, &x)); 229 PetscCall(VecGetArray(F, &f)); 230 for (j = 0; j < my; j++) { 231 for (i = 0; i < mx; i++) { 232 row = i + j * mx; 233 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 234 f[row] = x[row]; 235 continue; 236 } 237 u = x[row]; 238 ub = x[row - mx]; 239 ul = x[row - 1]; 240 ut = x[row + mx]; 241 ur = x[row + 1]; 242 uxx = (-ur + two * u - ul) * hydhx; 243 uyy = (-ut + two * u - ub) * hxdhy; 244 f[row] = -uxx + -uyy + sc * lambda * PetscExpScalar(u); 245 } 246 } 247 PetscCall(VecRestoreArrayRead(X, &x)); 248 PetscCall(VecRestoreArray(F, &f)); 249 return 0; 250 } 251 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 252 253 /* 254 Calculate the Jacobian matrix J(X,t). 255 256 Note: We put the Jacobian in the preconditioner storage B instead of J. This 257 way we can give the -snes_mf_operator flag to check our work. This replaces 258 J with a finite difference approximation, using our analytic Jacobian B for 259 the preconditioner. 260 */ 261 PetscErrorCode FormJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat B, void *ptr) { 262 AppCtx *user = (AppCtx *)ptr; 263 PetscInt i, j, row, mx, my, col[5]; 264 PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc; 265 const PetscScalar *x; 266 PetscReal hx, hy, hxdhy, hydhx; 267 268 mx = user->mx; 269 my = user->my; 270 lambda = user->param; 271 272 hx = 1.0 / (PetscReal)(mx - 1); 273 hy = 1.0 / (PetscReal)(my - 1); 274 sc = hx * hy; 275 hxdhy = hx / hy; 276 hydhx = hy / hx; 277 278 PetscCall(VecGetArrayRead(X, &x)); 279 for (j = 0; j < my; j++) { 280 for (i = 0; i < mx; i++) { 281 row = i + j * mx; 282 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 283 PetscCall(MatSetValues(B, 1, &row, 1, &row, &one, INSERT_VALUES)); 284 continue; 285 } 286 v[0] = hxdhy; 287 col[0] = row - mx; 288 v[1] = hydhx; 289 col[1] = row - 1; 290 v[2] = -two * (hydhx + hxdhy) + sc * lambda * PetscExpScalar(x[row]); 291 col[2] = row; 292 v[3] = hydhx; 293 col[3] = row + 1; 294 v[4] = hxdhy; 295 col[4] = row + mx; 296 PetscCall(MatSetValues(B, 1, &row, 5, col, v, INSERT_VALUES)); 297 } 298 } 299 PetscCall(VecRestoreArrayRead(X, &x)); 300 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 301 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 302 if (J != B) { 303 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 304 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 305 } 306 return 0; 307 } 308 309 /*TEST 310 311 test: 312 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 313 314 test: 315 suffix: 2 316 args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5 317 318 TEST*/ 319