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