1! 2! Solves the time dependent Bratu problem using pseudo-timestepping 3! 4! This code demonstrates how one may solve a nonlinear problem 5! with pseudo-timestepping. In this simple example, the pseudo-timestep 6! is the same for all grid points, i.e., this is equivalent to using 7! the backward Euler method with a variable timestep. 8! 9! Note: This example does not require pseudo-timestepping since it 10! is an easy nonlinear problem, but it is included to demonstrate how 11! the pseudo-timestepping may be done. 12! 13! See snes/tutorials/ex4.c[ex4f.F] and 14! snes/tutorials/ex5.c[ex5f.F] where the problem is described 15! and solved using the method of Newton alone. 16! 17! 18!23456789012345678901234567890123456789012345678901234567890123456789012 19 program main 20#include <petsc/finclude/petscts.h> 21 use petscts 22 implicit none 23 24! 25! Create an application context to contain data needed by the 26! application-provided call-back routines, FormJacobian() and 27! FormFunction(). We use a double precision array with three 28! entries indexed by param, lmx, lmy. 29! 30 PetscReal user(3) 31 PetscInt param,lmx,lmy,i5 32 parameter (param = 1,lmx = 2,lmy = 3) 33! 34! User-defined routines 35! 36 external FormJacobian,FormFunction 37! 38! Data for problem 39! 40 TS ts 41 Vec x,r 42 Mat J 43 PetscInt its,N,i1000,itmp 44 PetscBool flg 45 PetscErrorCode ierr 46 PetscReal param_max,param_min,dt 47 PetscReal tmax 48 PetscReal ftime 49 TSConvergedReason reason 50 51 i5 = 5 52 param_max = 6.81 53 param_min = 0 54 55 PetscCallA(PetscInitialize(ierr)) 56 user(lmx) = 4 57 user(lmy) = 4 58 user(param) = 6.0 59 60! 61! Allow user to set the grid dimensions and nonlinearity parameter at run-time 62! 63 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',user(lmx),flg,ierr)) 64 itmp = 4 65 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',itmp,flg,ierr)) 66 user(lmy) = itmp 67 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-param',user(param),flg,ierr)) 68 if (user(param) .ge. param_max .or. user(param) .le. param_min) then 69 print*,'Parameter is out of range' 70 endif 71 if (user(lmx) .gt. user(lmy)) then 72 dt = .5/user(lmx) 73 else 74 dt = .5/user(lmy) 75 endif 76 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)) 77 N = int(user(lmx)*user(lmy)) 78 79! 80! Create vectors to hold the solution and function value 81! 82 PetscCallA(VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)) 83 PetscCallA(VecDuplicate(x,r,ierr)) 84 85! 86! Create matrix to hold Jacobian. Preallocate 5 nonzeros per row 87! in the sparse matrix. Note that this is not the optimal strategy see 88! the Performance chapter of the users manual for information on 89! preallocating memory in sparse matrices. 90! 91 i5 = 5 92 PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER,J,ierr)) 93 94! 95! Create timestepper context 96! 97 98 PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr)) 99 PetscCallA(TSSetProblemType(ts,TS_NONLINEAR,ierr)) 100 101! 102! Tell the timestepper context where to compute solutions 103! 104 105 PetscCallA(TSSetSolution(ts,x,ierr)) 106 107! 108! Provide the call-back for the nonlinear function we are 109! evaluating. Thus whenever the timestepping routines need the 110! function they will call this routine. Note the final argument 111! is the application context used by the call-back functions. 112! 113 114 PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormFunction,user,ierr)) 115 116! 117! Set the Jacobian matrix and the function used to compute 118! Jacobians. 119! 120 121 PetscCallA(TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr)) 122 123! 124! For the initial guess for the problem 125! 126 127 PetscCallA(FormInitialGuess(x,user,ierr)) 128 129! 130! This indicates that we are using pseudo timestepping to 131! find a steady state solution to the nonlinear problem. 132! 133 134 PetscCallA(TSSetType(ts,TSPSEUDO,ierr)) 135 136! 137! Set the initial time to start at (this is arbitrary for 138! steady state problems and the initial timestep given above 139! 140 141 PetscCallA(TSSetTimeStep(ts,dt,ierr)) 142 143! 144! Set a large number of timesteps and final duration time 145! to insure convergence to steady state. 146! 147 i1000 = 1000 148 tmax = 1.e12 149 PetscCallA(TSSetMaxSteps(ts,i1000,ierr)) 150 PetscCallA(TSSetMaxTime(ts,tmax,ierr)) 151 PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)) 152 153! 154! Set any additional options from the options database. This 155! includes all options for the nonlinear and linear solvers used 156! internally the timestepping routines. 157! 158 159 PetscCallA(TSSetFromOptions(ts,ierr)) 160 161 PetscCallA(TSSetUp(ts,ierr)) 162 163! 164! Perform the solve. This is where the timestepping takes place. 165! 166 PetscCallA(TSSolve(ts,x,ierr)) 167 PetscCallA(TSGetSolveTime(ts,ftime,ierr)) 168 PetscCallA(TSGetStepNumber(ts,its,ierr)) 169 PetscCallA(TSGetConvergedReason(ts,reason,ierr)) 170 171 write(6,100) its,ftime,reason 172 100 format('Number of pseudo time-steps ',i5,' final time ',1pe9.2,' reason ',i3) 173 174! 175! Free the data structures constructed above 176! 177 178 PetscCallA(VecDestroy(x,ierr)) 179 PetscCallA(VecDestroy(r,ierr)) 180 PetscCallA(MatDestroy(J,ierr)) 181 PetscCallA(TSDestroy(ts,ierr)) 182 PetscCallA(PetscFinalize(ierr)) 183 end 184 185! 186! -------------------- Form initial approximation ----------------- 187! 188 subroutine FormInitialGuess(X,user,ierr) 189 use petscts 190 implicit none 191 192 Vec X 193 PetscReal user(3) 194 PetscInt i,j,row,mx,my 195 PetscErrorCode ierr 196 PetscOffset xidx 197 PetscReal one,lambda 198 PetscReal temp1,temp,hx,hy 199 PetscScalar xx(2) 200 PetscInt param,lmx,lmy 201 parameter (param = 1,lmx = 2,lmy = 3) 202 203 one = 1.0 204 205 mx = int(user(lmx)) 206 my = int(user(lmy)) 207 lambda = user(param) 208 209 hy = one / (my-1) 210 hx = one / (mx-1) 211 212 PetscCall(VecGetArray(X,xx,xidx,ierr)) 213 temp1 = lambda/(lambda + one) 214 do 10, j=1,my 215 temp = min(j-1,my-j)*hy 216 do 20 i=1,mx 217 row = i + (j-1)*mx 218 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 219 xx(row+xidx) = 0.0 220 else 221 xx(row+xidx) = temp1*sqrt(min(min(i-1,mx-i)*hx,temp)) 222 endif 223 20 continue 224 10 continue 225 PetscCall(VecRestoreArray(X,xx,xidx,ierr)) 226 return 227 end 228! 229! -------------------- Evaluate Function F(x) --------------------- 230! 231 subroutine FormFunction(ts,t,X,F,user,ierr) 232 use petscts 233 implicit none 234 235 TS ts 236 PetscReal t 237 Vec X,F 238 PetscReal user(3) 239 PetscErrorCode ierr 240 PetscInt i,j,row,mx,my 241 PetscOffset xidx,fidx 242 PetscReal two,lambda 243 PetscReal hx,hy,hxdhy,hydhx 244 PetscScalar ut,ub,ul,ur,u 245 PetscScalar uxx,uyy,sc 246 PetscScalar xx(2),ff(2) 247 PetscInt param,lmx,lmy 248 parameter (param = 1,lmx = 2,lmy = 3) 249 250 two = 2.0 251 252 mx = int(user(lmx)) 253 my = int(user(lmy)) 254 lambda = user(param) 255 256 hx = 1.0 / real(mx-1) 257 hy = 1.0 / real(my-1) 258 sc = hx*hy 259 hxdhy = hx/hy 260 hydhx = hy/hx 261 262 PetscCall(VecGetArrayRead(X,xx,xidx,ierr)) 263 PetscCall(VecGetArray(F,ff,fidx,ierr)) 264 do 10 j=1,my 265 do 20 i=1,mx 266 row = i + (j-1)*mx 267 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 268 ff(row+fidx) = xx(row+xidx) 269 else 270 u = xx(row + xidx) 271 ub = xx(row - mx + xidx) 272 ul = xx(row - 1 + xidx) 273 ut = xx(row + mx + xidx) 274 ur = xx(row + 1 + xidx) 275 uxx = (-ur + two*u - ul)*hydhx 276 uyy = (-ut + two*u - ub)*hxdhy 277 ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u) 278 u = -uxx - uyy + sc*lambda*exp(u) 279 endif 280 20 continue 281 10 continue 282 283 PetscCall(VecRestoreArrayRead(X,xx,xidx,ierr)) 284 PetscCall(VecRestoreArray(F,ff,fidx,ierr)) 285 return 286 end 287! 288! -------------------- Evaluate Jacobian of F(x) -------------------- 289! 290 subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr) 291 use petscts 292 implicit none 293 294 TS ts 295 Vec X 296 Mat JJ,B 297 PetscReal user(3),ctime 298 PetscErrorCode ierr 299 Mat jac 300 PetscOffset xidx 301 PetscInt i,j,row(1),mx,my 302 PetscInt col(5),i1,i5 303 PetscScalar two,one,lambda 304 PetscScalar v(5),sc,xx(2) 305 PetscReal hx,hy,hxdhy,hydhx 306 307 PetscInt param,lmx,lmy 308 parameter (param = 1,lmx = 2,lmy = 3) 309 310 i1 = 1 311 i5 = 5 312 jac = B 313 two = 2.0 314 one = 1.0 315 316 mx = int(user(lmx)) 317 my = int(user(lmy)) 318 lambda = user(param) 319 320 hx = 1.0 / real(mx-1) 321 hy = 1.0 / real(my-1) 322 sc = hx*hy 323 hxdhy = hx/hy 324 hydhx = hy/hx 325 326 PetscCall(VecGetArrayRead(X,xx,xidx,ierr)) 327 do 10 j=1,my 328 do 20 i=1,mx 329! 330! When inserting into PETSc matrices, indices start at 0 331! 332 row(1) = i - 1 + (j-1)*mx 333 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 334 PetscCall(MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)) 335 else 336 v(1) = hxdhy 337 col(1) = row(1) - mx 338 v(2) = hydhx 339 col(2) = row(1) - 1 340 v(3) = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx)) 341 col(3) = row(1) 342 v(4) = hydhx 343 col(4) = row(1) + 1 344 v(5) = hxdhy 345 col(5) = row(1) + mx 346 PetscCall(MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)) 347 endif 348 20 continue 349 10 continue 350 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 351 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 352 PetscCall(VecRestoreArray(X,xx,xidx,ierr)) 353 return 354 end 355 356!/*TEST 357! 358! test: 359! TODO: broken 360! args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo -ts_max_snes_failures 3 -ts_pseudo_frtol 1.e-5 -snes_stol 1e-5 361! 362!TEST*/ 363 364