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 PetscReal one,lambda 197 PetscReal temp1,temp,hx,hy 198 PetscScalar,pointer :: xx(:) 199 PetscInt param,lmx,lmy 200 parameter (param = 1,lmx = 2,lmy = 3) 201 202 one = 1.0 203 204 mx = int(user(lmx)) 205 my = int(user(lmy)) 206 lambda = user(param) 207 208 hy = one / (my-1) 209 hx = one / (mx-1) 210 211 PetscCall(VecGetArrayF90(X,xx,ierr)) 212 temp1 = lambda/(lambda + one) 213 do 10, j=1,my 214 temp = min(j-1,my-j)*hy 215 do 20 i=1,mx 216 row = i + (j-1)*mx 217 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 218 xx(row) = 0.0 219 else 220 xx(row) = temp1*sqrt(min(min(i-1,mx-i)*hx,temp)) 221 endif 222 20 continue 223 10 continue 224 PetscCall(VecRestoreArrayF90(X,xx,ierr)) 225 return 226 end 227! 228! -------------------- Evaluate Function F(x) --------------------- 229! 230 subroutine FormFunction(ts,t,X,F,user,ierr) 231 use petscts 232 implicit none 233 234 TS ts 235 PetscReal t 236 Vec X,F 237 PetscReal user(3) 238 PetscErrorCode ierr 239 PetscInt i,j,row,mx,my 240 PetscReal two,lambda 241 PetscReal hx,hy,hxdhy,hydhx 242 PetscScalar ut,ub,ul,ur,u 243 PetscScalar uxx,uyy,sc 244 PetscScalar,pointer :: xx(:), ff(:) 245 PetscInt param,lmx,lmy 246 parameter (param = 1,lmx = 2,lmy = 3) 247 248 two = 2.0 249 250 mx = int(user(lmx)) 251 my = int(user(lmy)) 252 lambda = user(param) 253 254 hx = 1.0 / real(mx-1) 255 hy = 1.0 / real(my-1) 256 sc = hx*hy 257 hxdhy = hx/hy 258 hydhx = hy/hx 259 260 PetscCall(VecGetArrayReadF90(X,xx,ierr)) 261 PetscCall(VecGetArrayF90(F,ff,ierr)) 262 do 10 j=1,my 263 do 20 i=1,mx 264 row = i + (j-1)*mx 265 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 266 ff(row) = xx(row) 267 else 268 u = xx(row) 269 ub = xx(row - mx) 270 ul = xx(row - 1) 271 ut = xx(row + mx) 272 ur = xx(row + 1) 273 uxx = (-ur + two*u - ul)*hydhx 274 uyy = (-ut + two*u - ub)*hxdhy 275 ff(row) = -uxx - uyy + sc*lambda*exp(u) 276 endif 277 20 continue 278 10 continue 279 280 PetscCall(VecRestoreArrayReadF90(X,xx,ierr)) 281 PetscCall(VecRestoreArrayF90(F,ff,ierr)) 282 return 283 end 284! 285! -------------------- Evaluate Jacobian of F(x) -------------------- 286! 287 subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr) 288 use petscts 289 implicit none 290 291 TS ts 292 Vec X 293 Mat JJ,B 294 PetscReal user(3),ctime 295 PetscErrorCode ierr 296 Mat jac 297 PetscInt i,j,row(1),mx,my 298 PetscInt col(5),i1,i5 299 PetscScalar two,one,lambda 300 PetscScalar v(5),sc 301 PetscScalar,pointer :: xx(:) 302 PetscReal hx,hy,hxdhy,hydhx 303 304 PetscInt param,lmx,lmy 305 parameter (param = 1,lmx = 2,lmy = 3) 306 307 i1 = 1 308 i5 = 5 309 jac = B 310 two = 2.0 311 one = 1.0 312 313 mx = int(user(lmx)) 314 my = int(user(lmy)) 315 lambda = user(param) 316 317 hx = 1.0 / real(mx-1) 318 hy = 1.0 / real(my-1) 319 sc = hx*hy 320 hxdhy = hx/hy 321 hydhx = hy/hx 322 323 PetscCall(VecGetArrayReadF90(X,xx,ierr)) 324 do 10 j=1,my 325 do 20 i=1,mx 326! 327! When inserting into PETSc matrices, indices start at 0 328! 329 row(1) = i - 1 + (j-1)*mx 330 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 331 PetscCall(MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)) 332 else 333 v(1) = hxdhy 334 col(1) = row(1) - mx 335 v(2) = hydhx 336 col(2) = row(1) - 1 337 v(3) = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1))) 338 col(3) = row(1) 339 v(4) = hydhx 340 col(4) = row(1) + 1 341 v(5) = hxdhy 342 col(5) = row(1) + mx 343 PetscCall(MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)) 344 endif 345 20 continue 346 10 continue 347 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 348 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 349 PetscCall(VecRestoreArrayF90(X,xx,ierr)) 350 return 351 end 352 353!/*TEST 354! 355! test: 356! TODO: broken 357! 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 358! 359!TEST*/ 360