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