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>
 20: module ex1fmodule
 21:   use petscts
 22:   implicit none
 23: contains
 24: !
 25: !  --------------------  Evaluate Function F(x) ---------------------
 26: !
 27:   subroutine FormFunction(ts, t, X, F, user, ierr)

 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

 42:     two = 2.0

 44:     mx = int(user(lmx))
 45:     my = int(user(lmy))
 46:     lambda = user(param)

 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

 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

 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)

 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

 95:     PetscInt, parameter :: param = 1, lmx = 2, lmy = 3

 97:     i1 = 1
 98:     i5 = 5
 99:     jac = B
100:     two = 2.0
101:     one = 1.0

103:     mx = int(user(lmx))
104:     my = int(user(lmy))
105:     lambda = user(param)

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

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)

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

155:     one = 1.0

157:     mx = int(user(lmx))
158:     my = int(user(lmy))
159:     lambda = user(param)

161:     hy = one/(my - 1)
162:     hx = one/(mx - 1)

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
179: end module ex1fmodule
180: program 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

207:   i5 = 5
208:   param_max = 6.81
209:   param_min = 0

211:   PetscCallA(PetscInitialize(ierr))
212:   user(lmx) = 4
213:   user(lmy) = 4
214:   user(param) = 6.0

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))

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))

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))

250: !
251: !     Create timestepper context
252: !

254:   PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
255:   PetscCallA(TSSetProblemType(ts, TS_NONLINEAR, ierr))

257: !
258: !     Tell the timestepper context where to compute solutions
259: !

261:   PetscCallA(TSSetSolution(ts, x, ierr))

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: !

270:   PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormFunction, user, ierr))

272: !
273: !     Set the Jacobian matrix and the function used to compute
274: !     Jacobians.
275: !

277:   PetscCallA(TSSetRHSJacobian(ts, J, J, FormJacobian, user, ierr))

279: !
280: !       For the initial guess for the problem
281: !

283:   PetscCallA(FormInitialGuess(x, user, ierr))

285: !
286: !       This indicates that we are using pseudo timestepping to
287: !     find a steady state solution to the nonlinear problem.
288: !

290:   PetscCallA(TSSetType(ts, TSPSEUDO, ierr))

292: !
293: !       Set the initial time to start at (this is arbitrary for
294: !     steady state problems and the initial timestep given above
295: !

297:   PetscCallA(TSSetTimeStep(ts, dt, ierr))

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))

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: !

315:   PetscCallA(TSSetFromOptions(ts, ierr))

317:   PetscCallA(TSSetUp(ts, ierr))

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))

327:   write (6, 100) its, ftime, reason
328: 100 format('Number of pseudo time-steps ', i5, ' final time ', 1pe9.2, ' reason ', i3)

330: !
331: !     Free the data structures constructed above
332: !

334:   PetscCallA(VecDestroy(x, ierr))
335:   PetscCallA(VecDestroy(r, ierr))
336:   PetscCallA(MatDestroy(J, ierr))
337:   PetscCallA(TSDestroy(ts, ierr))
338:   PetscCallA(PetscFinalize(ierr))
339: end
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*/