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 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 65 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"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 PetscCheck(user.param < param_max && user.param >= param_min,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 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,N,&x)); 85 CHKERRQ(VecDuplicate(x,&r)); 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 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J)); 94 95 /* 96 Create timestepper context 97 */ 98 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 99 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 100 101 /* 102 Tell the timestepper context where to compute solutions 103 */ 104 CHKERRQ(TSSetSolution(ts,x)); 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 CHKERRQ(TSSetRHSFunction(ts,NULL,FormFunction,&user)); 113 114 /* 115 Set the Jacobian matrix and the function used to compute 116 Jacobians. 117 */ 118 CHKERRQ(TSSetRHSJacobian(ts,J,J,FormJacobian,&user)); 119 120 /* 121 Form the initial guess for the problem 122 */ 123 CHKERRQ(FormInitialGuess(x,&user)); 124 125 /* 126 This indicates that we are using pseudo timestepping to 127 find a steady state solution to the nonlinear problem. 128 */ 129 CHKERRQ(TSSetType(ts,TSPSEUDO)); 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 CHKERRQ(TSSetTimeStep(ts,dt)); 136 137 /* 138 Set a large number of timesteps and final duration time 139 to insure convergence to steady state. 140 */ 141 CHKERRQ(TSSetMaxSteps(ts,10000)); 142 CHKERRQ(TSSetMaxTime(ts,1e12)); 143 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 144 145 /* 146 Use the default strategy for increasing the timestep 147 */ 148 CHKERRQ(TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0)); 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 CHKERRQ(TSSetFromOptions(ts)); 156 157 CHKERRQ(TSSetUp(ts)); 158 159 /* 160 Perform the solve. This is where the timestepping takes place. 161 */ 162 CHKERRQ(TSSolve(ts,x)); 163 CHKERRQ(TSGetSolveTime(ts,&ftime)); 164 165 /* 166 Get the number of steps 167 */ 168 CHKERRQ(TSGetStepNumber(ts,&its)); 169 170 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime)); 171 172 /* 173 Free the data structures constructed above 174 */ 175 CHKERRQ(VecDestroy(&x)); 176 CHKERRQ(VecDestroy(&r)); 177 CHKERRQ(MatDestroy(&J)); 178 CHKERRQ(TSDestroy(&ts)); 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 PetscReal one = 1.0,lambda; 192 PetscReal temp1,temp,hx,hy; 193 PetscScalar *x; 194 195 mx = user->mx; 196 my = user->my; 197 lambda = user->param; 198 199 hx = one / (PetscReal)(mx-1); 200 hy = one / (PetscReal)(my-1); 201 202 CHKERRQ(VecGetArray(X,&x)); 203 temp1 = lambda/(lambda + one); 204 for (j=0; j<my; j++) { 205 temp = (PetscReal)(PetscMin(j,my-j-1))*hy; 206 for (i=0; i<mx; i++) { 207 row = i + j*mx; 208 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 209 x[row] = 0.0; 210 continue; 211 } 212 x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp)); 213 } 214 } 215 CHKERRQ(VecRestoreArray(X,&x)); 216 return 0; 217 } 218 /* -------------------- Evaluate Function F(x) --------------------- */ 219 220 PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 221 { 222 AppCtx *user = (AppCtx*)ptr; 223 PetscInt i,j,row,mx,my; 224 PetscReal two = 2.0,one = 1.0,lambda; 225 PetscReal hx,hy,hxdhy,hydhx; 226 PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f; 227 const PetscScalar *x; 228 229 mx = user->mx; 230 my = user->my; 231 lambda = user->param; 232 233 hx = one / (PetscReal)(mx-1); 234 hy = one / (PetscReal)(my-1); 235 sc = hx*hy; 236 hxdhy = hx/hy; 237 hydhx = hy/hx; 238 239 CHKERRQ(VecGetArrayRead(X,&x)); 240 CHKERRQ(VecGetArray(F,&f)); 241 for (j=0; j<my; j++) { 242 for (i=0; i<mx; i++) { 243 row = i + j*mx; 244 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 245 f[row] = x[row]; 246 continue; 247 } 248 u = x[row]; 249 ub = x[row - mx]; 250 ul = x[row - 1]; 251 ut = x[row + mx]; 252 ur = x[row + 1]; 253 uxx = (-ur + two*u - ul)*hydhx; 254 uyy = (-ut + two*u - ub)*hxdhy; 255 f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u); 256 } 257 } 258 CHKERRQ(VecRestoreArrayRead(X,&x)); 259 CHKERRQ(VecRestoreArray(F,&f)); 260 return 0; 261 } 262 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 263 264 /* 265 Calculate the Jacobian matrix J(X,t). 266 267 Note: We put the Jacobian in the preconditioner storage B instead of J. This 268 way we can give the -snes_mf_operator flag to check our work. This replaces 269 J with a finite difference approximation, using our analytic Jacobian B for 270 the preconditioner. 271 */ 272 PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr) 273 { 274 AppCtx *user = (AppCtx*)ptr; 275 PetscInt i,j,row,mx,my,col[5]; 276 PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc; 277 const PetscScalar *x; 278 PetscReal hx,hy,hxdhy,hydhx; 279 280 mx = user->mx; 281 my = user->my; 282 lambda = user->param; 283 284 hx = 1.0 / (PetscReal)(mx-1); 285 hy = 1.0 / (PetscReal)(my-1); 286 sc = hx*hy; 287 hxdhy = hx/hy; 288 hydhx = hy/hx; 289 290 CHKERRQ(VecGetArrayRead(X,&x)); 291 for (j=0; j<my; j++) { 292 for (i=0; i<mx; i++) { 293 row = i + j*mx; 294 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 295 CHKERRQ(MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES)); 296 continue; 297 } 298 v[0] = hxdhy; col[0] = row - mx; 299 v[1] = hydhx; col[1] = row - 1; 300 v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row; 301 v[3] = hydhx; col[3] = row + 1; 302 v[4] = hxdhy; col[4] = row + mx; 303 CHKERRQ(MatSetValues(B,1,&row,5,col,v,INSERT_VALUES)); 304 } 305 } 306 CHKERRQ(VecRestoreArrayRead(X,&x)); 307 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 308 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 309 if (J != B) { 310 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 311 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 312 } 313 return 0; 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