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#include <petsc/finclude/petscts.h> 20program main 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) >= param_max .or. user(param) <= param_min) then 69 print *, 'Parameter is out of range' 70 end if 71 if (user(lmx) > user(lmy)) then 72 dt = .5/user(lmx) 73 else 74 dt = .5/user(lmy) 75 end if 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_ARRAY, 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 172100 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)) 183end 184 185! 186! -------------------- Form initial approximation ----------------- 187! 188subroutine 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(VecGetArray(X, xx, ierr)) 212 temp1 = lambda/(lambda + one) 213 do j = 1, my 214 temp = min(j - 1, my - j)*hy 215 do i = 1, mx 216 row = i + (j - 1)*mx 217 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 218 xx(row) = 0.0 219 else 220 xx(row) = temp1*sqrt(min(min(i - 1, mx - i)*hx, temp)) 221 end if 222 end do 223 end do 224 PetscCall(VecRestoreArray(X, xx, ierr)) 225end 226! 227! -------------------- Evaluate Function F(x) --------------------- 228! 229subroutine FormFunction(ts, t, X, F, user, ierr) 230 use petscts 231 implicit none 232 233 TS ts 234 PetscReal t 235 Vec X, F 236 PetscReal user(3) 237 PetscErrorCode ierr 238 PetscInt i, j, row, mx, my 239 PetscReal two, lambda 240 PetscReal hx, hy, hxdhy, hydhx 241 PetscScalar ut, ub, ul, ur, u 242 PetscScalar uxx, uyy, sc 243 PetscScalar, pointer :: xx(:), ff(:) 244 PetscInt param, lmx, lmy 245 parameter(param=1, lmx=2, lmy=3) 246 247 two = 2.0 248 249 mx = int(user(lmx)) 250 my = int(user(lmy)) 251 lambda = user(param) 252 253 hx = 1.0/real(mx - 1) 254 hy = 1.0/real(my - 1) 255 sc = hx*hy 256 hxdhy = hx/hy 257 hydhx = hy/hx 258 259 PetscCall(VecGetArrayRead(X, xx, ierr)) 260 PetscCall(VecGetArray(F, ff, ierr)) 261 do j = 1, my 262 do i = 1, mx 263 row = i + (j - 1)*mx 264 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 265 ff(row) = xx(row) 266 else 267 u = xx(row) 268 ub = xx(row - mx) 269 ul = xx(row - 1) 270 ut = xx(row + mx) 271 ur = xx(row + 1) 272 uxx = (-ur + two*u - ul)*hydhx 273 uyy = (-ut + two*u - ub)*hxdhy 274 ff(row) = -uxx - uyy + sc*lambda*exp(u) 275 end if 276 end do 277 end do 278 279 PetscCall(VecRestoreArrayRead(X, xx, ierr)) 280 PetscCall(VecRestoreArray(F, ff, ierr)) 281end 282! 283! -------------------- Evaluate Jacobian of F(x) -------------------- 284! 285subroutine FormJacobian(ts, ctime, X, JJ, B, user, ierr) 286 use petscts 287 implicit none 288 289 TS ts 290 Vec X 291 Mat JJ, B 292 PetscReal user(3), ctime 293 PetscErrorCode ierr 294 Mat jac 295 PetscInt i, j, row(1), mx, my 296 PetscInt col(5), i1, i5 297 PetscScalar two, one, lambda 298 PetscScalar v(5), sc 299 PetscScalar, pointer :: xx(:) 300 PetscReal hx, hy, hxdhy, hydhx 301 302 PetscInt param, lmx, lmy 303 parameter(param=1, lmx=2, lmy=3) 304 305 i1 = 1 306 i5 = 5 307 jac = B 308 two = 2.0 309 one = 1.0 310 311 mx = int(user(lmx)) 312 my = int(user(lmy)) 313 lambda = user(param) 314 315 hx = 1.0/real(mx - 1) 316 hy = 1.0/real(my - 1) 317 sc = hx*hy 318 hxdhy = hx/hy 319 hydhx = hy/hx 320 321 PetscCall(VecGetArrayRead(X, xx, ierr)) 322 do j = 1, my 323 do i = 1, mx 324! 325! When inserting into PETSc matrices, indices start at 0 326! 327 row(1) = i - 1 + (j - 1)*mx 328 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 329 PetscCall(MatSetValues(jac, i1, [row], i1, [row], [one], INSERT_VALUES, ierr)) 330 else 331 v(1) = hxdhy 332 col(1) = row(1) - mx 333 v(2) = hydhx 334 col(2) = row(1) - 1 335 v(3) = -two*(hydhx + hxdhy) + sc*lambda*exp(xx(row(1))) 336 col(3) = row(1) 337 v(4) = hydhx 338 col(4) = row(1) + 1 339 v(5) = hxdhy 340 col(5) = row(1) + mx 341 PetscCall(MatSetValues(jac, i1, [row], i5, col, v, INSERT_VALUES, ierr)) 342 end if 343 end do 344 end do 345 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 346 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 347 PetscCall(VecRestoreArray(X, xx, ierr)) 348end 349 350!/*TEST 351! 352! test: 353! TODO: broken 354! 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 355! 356!TEST*/ 357