1 2 static char help[] ="Solves the time independent Bratu problem using pseudo-timestepping."; 3 4 /* 5 Concepts: TS^pseudo-timestepping 6 Concepts: TS^pseudo-timestepping 7 Concepts: TS^nonlinear problems 8 Processors: 1 9 10 */ 11 12 /* ------------------------------------------------------------------------ 13 14 This code demonstrates how one may solve a nonlinear problem 15 with pseudo-timestepping. In this simple example, the pseudo-timestep 16 is the same for all grid points, i.e., this is equivalent to using 17 the backward Euler method with a variable timestep. 18 19 Note: This example does not require pseudo-timestepping since it 20 is an easy nonlinear problem, but it is included to demonstrate how 21 the pseudo-timestepping may be done. 22 23 See snes/tutorials/ex4.c[ex4f.F] and 24 snes/tutorials/ex5.c[ex5f.F] where the problem is described 25 and solved using Newton's method alone. 26 27 ----------------------------------------------------------------------------- */ 28 /* 29 Include "petscts.h" to use the PETSc timestepping routines. Note that 30 this file automatically includes "petscsys.h" and other lower-level 31 PETSc include files. 32 */ 33 #include <petscts.h> 34 35 /* 36 Create an application context to contain data needed by the 37 application-provided call-back routines, FormJacobian() and 38 FormFunction(). 39 */ 40 typedef struct { 41 PetscReal param; /* test problem parameter */ 42 PetscInt mx; /* Discretization in x-direction */ 43 PetscInt my; /* Discretization in y-direction */ 44 } AppCtx; 45 46 /* 47 User-defined routines 48 */ 49 extern PetscErrorCode FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*), FormFunction(TS,PetscReal,Vec,Vec,void*), FormInitialGuess(Vec,AppCtx*); 50 51 int main(int argc,char **argv) 52 { 53 TS ts; /* timestepping context */ 54 Vec x,r; /* solution, residual vectors */ 55 Mat J; /* Jacobian matrix */ 56 AppCtx user; /* user-defined work context */ 57 PetscInt its,N; /* iterations for convergence */ 58 PetscErrorCode ierr; 59 PetscReal param_max = 6.81,param_min = 0.,dt; 60 PetscReal ftime; 61 PetscMPIInt size; 62 63 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 64 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 65 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only"); 66 67 user.mx = 4; 68 user.my = 4; 69 user.param = 6.0; 70 71 /* 72 Allow user to set the grid dimensions and nonlinearity parameter at run-time 73 */ 74 PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL); 75 PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL); 76 N = user.mx*user.my; 77 dt = .5/PetscMax(user.mx,user.my); 78 PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL); 79 if (user.param >= param_max || user.param <= param_min) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Parameter is out of range"); 80 81 /* 82 Create vectors to hold the solution and function value 83 */ 84 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr); 85 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 86 87 /* 88 Create matrix to hold Jacobian. Preallocate 5 nonzeros per row 89 in the sparse matrix. Note that this is not the optimal strategy; see 90 the Performance chapter of the users manual for information on 91 preallocating memory in sparse matrices. 92 */ 93 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J);CHKERRQ(ierr); 94 95 /* 96 Create timestepper context 97 */ 98 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 99 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 100 101 /* 102 Tell the timestepper context where to compute solutions 103 */ 104 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 105 106 /* 107 Provide the call-back for the nonlinear function we are 108 evaluating. Thus whenever the timestepping routines need the 109 function they will call this routine. Note the final argument 110 is the application context used by the call-back functions. 111 */ 112 ierr = TSSetRHSFunction(ts,NULL,FormFunction,&user);CHKERRQ(ierr); 113 114 /* 115 Set the Jacobian matrix and the function used to compute 116 Jacobians. 117 */ 118 ierr = TSSetRHSJacobian(ts,J,J,FormJacobian,&user);CHKERRQ(ierr); 119 120 /* 121 Form the initial guess for the problem 122 */ 123 ierr = FormInitialGuess(x,&user);CHKERRQ(ierr); 124 125 /* 126 This indicates that we are using pseudo timestepping to 127 find a steady state solution to the nonlinear problem. 128 */ 129 ierr = TSSetType(ts,TSPSEUDO);CHKERRQ(ierr); 130 131 /* 132 Set the initial time to start at (this is arbitrary for 133 steady state problems); and the initial timestep given above 134 */ 135 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 136 137 /* 138 Set a large number of timesteps and final duration time 139 to insure convergence to steady state. 140 */ 141 ierr = TSSetMaxSteps(ts,10000);CHKERRQ(ierr); 142 ierr = TSSetMaxTime(ts,1e12);CHKERRQ(ierr); 143 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 144 145 /* 146 Use the default strategy for increasing the timestep 147 */ 148 ierr = TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0);CHKERRQ(ierr); 149 150 /* 151 Set any additional options from the options database. This 152 includes all options for the nonlinear and linear solvers used 153 internally the timestepping routines. 154 */ 155 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 156 157 ierr = TSSetUp(ts);CHKERRQ(ierr); 158 159 /* 160 Perform the solve. This is where the timestepping takes place. 161 */ 162 ierr = TSSolve(ts,x);CHKERRQ(ierr); 163 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 164 165 /* 166 Get the number of steps 167 */ 168 ierr = TSGetStepNumber(ts,&its);CHKERRQ(ierr); 169 170 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime);CHKERRQ(ierr); 171 172 /* 173 Free the data structures constructed above 174 */ 175 ierr = VecDestroy(&x);CHKERRQ(ierr); 176 ierr = VecDestroy(&r);CHKERRQ(ierr); 177 ierr = MatDestroy(&J);CHKERRQ(ierr); 178 ierr = TSDestroy(&ts);CHKERRQ(ierr); 179 ierr = PetscFinalize(); 180 return ierr; 181 } 182 /* ------------------------------------------------------------------ */ 183 /* Bratu (Solid Fuel Ignition) Test Problem */ 184 /* ------------------------------------------------------------------ */ 185 186 /* -------------------- Form initial approximation ----------------- */ 187 188 PetscErrorCode FormInitialGuess(Vec X,AppCtx *user) 189 { 190 PetscInt i,j,row,mx,my; 191 PetscErrorCode ierr; 192 PetscReal one = 1.0,lambda; 193 PetscReal temp1,temp,hx,hy; 194 PetscScalar *x; 195 196 mx = user->mx; 197 my = user->my; 198 lambda = user->param; 199 200 hx = one / (PetscReal)(mx-1); 201 hy = one / (PetscReal)(my-1); 202 203 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 204 temp1 = lambda/(lambda + one); 205 for (j=0; j<my; j++) { 206 temp = (PetscReal)(PetscMin(j,my-j-1))*hy; 207 for (i=0; i<mx; i++) { 208 row = i + j*mx; 209 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 210 x[row] = 0.0; 211 continue; 212 } 213 x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp)); 214 } 215 } 216 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 217 return 0; 218 } 219 /* -------------------- Evaluate Function F(x) --------------------- */ 220 221 PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 222 { 223 AppCtx *user = (AppCtx*)ptr; 224 PetscErrorCode ierr; 225 PetscInt i,j,row,mx,my; 226 PetscReal two = 2.0,one = 1.0,lambda; 227 PetscReal hx,hy,hxdhy,hydhx; 228 PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f; 229 const PetscScalar *x; 230 231 mx = user->mx; 232 my = user->my; 233 lambda = user->param; 234 235 hx = one / (PetscReal)(mx-1); 236 hy = one / (PetscReal)(my-1); 237 sc = hx*hy; 238 hxdhy = hx/hy; 239 hydhx = hy/hx; 240 241 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 242 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 243 for (j=0; j<my; j++) { 244 for (i=0; i<mx; i++) { 245 row = i + j*mx; 246 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 247 f[row] = x[row]; 248 continue; 249 } 250 u = x[row]; 251 ub = x[row - mx]; 252 ul = x[row - 1]; 253 ut = x[row + mx]; 254 ur = x[row + 1]; 255 uxx = (-ur + two*u - ul)*hydhx; 256 uyy = (-ut + two*u - ub)*hxdhy; 257 f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u); 258 } 259 } 260 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 261 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 262 return 0; 263 } 264 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 265 266 /* 267 Calculate the Jacobian matrix J(X,t). 268 269 Note: We put the Jacobian in the preconditioner storage B instead of J. This 270 way we can give the -snes_mf_operator flag to check our work. This replaces 271 J with a finite difference approximation, using our analytic Jacobian B for 272 the preconditioner. 273 */ 274 PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr) 275 { 276 AppCtx *user = (AppCtx*)ptr; 277 PetscInt i,j,row,mx,my,col[5]; 278 PetscErrorCode ierr; 279 PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc; 280 const PetscScalar *x; 281 PetscReal hx,hy,hxdhy,hydhx; 282 283 mx = user->mx; 284 my = user->my; 285 lambda = user->param; 286 287 hx = 1.0 / (PetscReal)(mx-1); 288 hy = 1.0 / (PetscReal)(my-1); 289 sc = hx*hy; 290 hxdhy = hx/hy; 291 hydhx = hy/hx; 292 293 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 294 for (j=0; j<my; j++) { 295 for (i=0; i<mx; i++) { 296 row = i + j*mx; 297 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 298 ierr = MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES);CHKERRQ(ierr); 299 continue; 300 } 301 v[0] = hxdhy; col[0] = row - mx; 302 v[1] = hydhx; col[1] = row - 1; 303 v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row; 304 v[3] = hydhx; col[3] = row + 1; 305 v[4] = hxdhy; col[4] = row + mx; 306 ierr = MatSetValues(B,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr); 307 } 308 } 309 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 310 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 311 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 312 if (J != B) { 313 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 314 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 315 } 316 return 0; 317 } 318 319 /*TEST 320 321 test: 322 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 323 324 test: 325 suffix: 2 326 args: -ts_monitor_pseudo -ts_pseudo_frtol 1.e-5 327 328 TEST*/ 329