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 PetscReal param_max = 6.81,param_min = 0.,dt; 59 PetscReal ftime; 60 PetscMPIInt size; 61 62 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 63 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 64 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only"); 65 66 user.mx = 4; 67 user.my = 4; 68 user.param = 6.0; 69 70 /* 71 Allow user to set the grid dimensions and nonlinearity parameter at run-time 72 */ 73 PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL); 74 PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL); 75 N = user.mx*user.my; 76 dt = .5/PetscMax(user.mx,user.my); 77 PetscOptionsGetReal(NULL,NULL,"-param",&user.param,NULL); 78 PetscCheck(user.param < param_max && user.param >= param_min,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Parameter is out of range"); 79 80 /* 81 Create vectors to hold the solution and function value 82 */ 83 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,N,&x)); 84 CHKERRQ(VecDuplicate(x,&r)); 85 86 /* 87 Create matrix to hold Jacobian. Preallocate 5 nonzeros per row 88 in the sparse matrix. Note that this is not the optimal strategy; see 89 the Performance chapter of the users manual for information on 90 preallocating memory in sparse matrices. 91 */ 92 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&J)); 93 94 /* 95 Create timestepper context 96 */ 97 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 98 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 99 100 /* 101 Tell the timestepper context where to compute solutions 102 */ 103 CHKERRQ(TSSetSolution(ts,x)); 104 105 /* 106 Provide the call-back for the nonlinear function we are 107 evaluating. Thus whenever the timestepping routines need the 108 function they will call this routine. Note the final argument 109 is the application context used by the call-back functions. 110 */ 111 CHKERRQ(TSSetRHSFunction(ts,NULL,FormFunction,&user)); 112 113 /* 114 Set the Jacobian matrix and the function used to compute 115 Jacobians. 116 */ 117 CHKERRQ(TSSetRHSJacobian(ts,J,J,FormJacobian,&user)); 118 119 /* 120 Form the initial guess for the problem 121 */ 122 CHKERRQ(FormInitialGuess(x,&user)); 123 124 /* 125 This indicates that we are using pseudo timestepping to 126 find a steady state solution to the nonlinear problem. 127 */ 128 CHKERRQ(TSSetType(ts,TSPSEUDO)); 129 130 /* 131 Set the initial time to start at (this is arbitrary for 132 steady state problems); and the initial timestep given above 133 */ 134 CHKERRQ(TSSetTimeStep(ts,dt)); 135 136 /* 137 Set a large number of timesteps and final duration time 138 to insure convergence to steady state. 139 */ 140 CHKERRQ(TSSetMaxSteps(ts,10000)); 141 CHKERRQ(TSSetMaxTime(ts,1e12)); 142 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 143 144 /* 145 Use the default strategy for increasing the timestep 146 */ 147 CHKERRQ(TSPseudoSetTimeStep(ts,TSPseudoTimeStepDefault,0)); 148 149 /* 150 Set any additional options from the options database. This 151 includes all options for the nonlinear and linear solvers used 152 internally the timestepping routines. 153 */ 154 CHKERRQ(TSSetFromOptions(ts)); 155 156 CHKERRQ(TSSetUp(ts)); 157 158 /* 159 Perform the solve. This is where the timestepping takes place. 160 */ 161 CHKERRQ(TSSolve(ts,x)); 162 CHKERRQ(TSGetSolveTime(ts,&ftime)); 163 164 /* 165 Get the number of steps 166 */ 167 CHKERRQ(TSGetStepNumber(ts,&its)); 168 169 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Number of pseudo timesteps = %D final time %4.2e\n",its,(double)ftime)); 170 171 /* 172 Free the data structures constructed above 173 */ 174 CHKERRQ(VecDestroy(&x)); 175 CHKERRQ(VecDestroy(&r)); 176 CHKERRQ(MatDestroy(&J)); 177 CHKERRQ(TSDestroy(&ts)); 178 CHKERRQ(PetscFinalize()); 179 return 0; 180 } 181 /* ------------------------------------------------------------------ */ 182 /* Bratu (Solid Fuel Ignition) Test Problem */ 183 /* ------------------------------------------------------------------ */ 184 185 /* -------------------- Form initial approximation ----------------- */ 186 187 PetscErrorCode FormInitialGuess(Vec X,AppCtx *user) 188 { 189 PetscInt i,j,row,mx,my; 190 PetscReal one = 1.0,lambda; 191 PetscReal temp1,temp,hx,hy; 192 PetscScalar *x; 193 194 mx = user->mx; 195 my = user->my; 196 lambda = user->param; 197 198 hx = one / (PetscReal)(mx-1); 199 hy = one / (PetscReal)(my-1); 200 201 CHKERRQ(VecGetArray(X,&x)); 202 temp1 = lambda/(lambda + one); 203 for (j=0; j<my; j++) { 204 temp = (PetscReal)(PetscMin(j,my-j-1))*hy; 205 for (i=0; i<mx; i++) { 206 row = i + j*mx; 207 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 208 x[row] = 0.0; 209 continue; 210 } 211 x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp)); 212 } 213 } 214 CHKERRQ(VecRestoreArray(X,&x)); 215 return 0; 216 } 217 /* -------------------- Evaluate Function F(x) --------------------- */ 218 219 PetscErrorCode FormFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 220 { 221 AppCtx *user = (AppCtx*)ptr; 222 PetscInt i,j,row,mx,my; 223 PetscReal two = 2.0,one = 1.0,lambda; 224 PetscReal hx,hy,hxdhy,hydhx; 225 PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f; 226 const PetscScalar *x; 227 228 mx = user->mx; 229 my = user->my; 230 lambda = user->param; 231 232 hx = one / (PetscReal)(mx-1); 233 hy = one / (PetscReal)(my-1); 234 sc = hx*hy; 235 hxdhy = hx/hy; 236 hydhx = hy/hx; 237 238 CHKERRQ(VecGetArrayRead(X,&x)); 239 CHKERRQ(VecGetArray(F,&f)); 240 for (j=0; j<my; j++) { 241 for (i=0; i<mx; i++) { 242 row = i + j*mx; 243 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 244 f[row] = x[row]; 245 continue; 246 } 247 u = x[row]; 248 ub = x[row - mx]; 249 ul = x[row - 1]; 250 ut = x[row + mx]; 251 ur = x[row + 1]; 252 uxx = (-ur + two*u - ul)*hydhx; 253 uyy = (-ut + two*u - ub)*hxdhy; 254 f[row] = -uxx + -uyy + sc*lambda*PetscExpScalar(u); 255 } 256 } 257 CHKERRQ(VecRestoreArrayRead(X,&x)); 258 CHKERRQ(VecRestoreArray(F,&f)); 259 return 0; 260 } 261 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 262 263 /* 264 Calculate the Jacobian matrix J(X,t). 265 266 Note: We put the Jacobian in the preconditioner storage B instead of J. This 267 way we can give the -snes_mf_operator flag to check our work. This replaces 268 J with a finite difference approximation, using our analytic Jacobian B for 269 the preconditioner. 270 */ 271 PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat J,Mat B,void *ptr) 272 { 273 AppCtx *user = (AppCtx*)ptr; 274 PetscInt i,j,row,mx,my,col[5]; 275 PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc; 276 const PetscScalar *x; 277 PetscReal hx,hy,hxdhy,hydhx; 278 279 mx = user->mx; 280 my = user->my; 281 lambda = user->param; 282 283 hx = 1.0 / (PetscReal)(mx-1); 284 hy = 1.0 / (PetscReal)(my-1); 285 sc = hx*hy; 286 hxdhy = hx/hy; 287 hydhx = hy/hx; 288 289 CHKERRQ(VecGetArrayRead(X,&x)); 290 for (j=0; j<my; j++) { 291 for (i=0; i<mx; i++) { 292 row = i + j*mx; 293 if (i == 0 || j == 0 || i == mx-1 || j == my-1) { 294 CHKERRQ(MatSetValues(B,1,&row,1,&row,&one,INSERT_VALUES)); 295 continue; 296 } 297 v[0] = hxdhy; col[0] = row - mx; 298 v[1] = hydhx; col[1] = row - 1; 299 v[2] = -two*(hydhx + hxdhy) + sc*lambda*PetscExpScalar(x[row]); col[2] = row; 300 v[3] = hydhx; col[3] = row + 1; 301 v[4] = hxdhy; col[4] = row + mx; 302 CHKERRQ(MatSetValues(B,1,&row,5,col,v,INSERT_VALUES)); 303 } 304 } 305 CHKERRQ(VecRestoreArrayRead(X,&x)); 306 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 307 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 308 if (J != B) { 309 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 310 CHKERRQ(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 311 } 312 return 0; 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