Actual source code: ex1f.F90

      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

     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

     51:   i5 = 5
     52:   param_max = 6.81
     53:   param_min = 0

     55:   PetscCallA(PetscInitialize(ierr))
     56:   user(lmx) = 4
     57:   user(lmy) = 4
     58:   user(param) = 6.0

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

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

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

     94: !
     95: !     Create timestepper context
     96: !

     98:   PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
     99:   PetscCallA(TSSetProblemType(ts, TS_NONLINEAR, ierr))

    101: !
    102: !     Tell the timestepper context where to compute solutions
    103: !

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

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

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

    116: !
    117: !     Set the Jacobian matrix and the function used to compute
    118: !     Jacobians.
    119: !

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

    123: !
    124: !       For the initial guess for the problem
    125: !

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

    129: !
    130: !       This indicates that we are using pseudo timestepping to
    131: !     find a steady state solution to the nonlinear problem.
    132: !

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

    136: !
    137: !       Set the initial time to start at (this is arbitrary for
    138: !     steady state problems and the initial timestep given above
    139: !

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

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

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

    159:   PetscCallA(TSSetFromOptions(ts, ierr))

    161:   PetscCallA(TSSetUp(ts, ierr))

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

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

    174: !
    175: !     Free the data structures constructed above
    176: !

    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

    185: !
    186: !  --------------------  Form initial approximation -----------------
    187: !
    188: subroutine FormInitialGuess(X, user, ierr)
    189:   use petscts
    190:   implicit none

    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)

    202:   one = 1.0

    204:   mx = int(user(lmx))
    205:   my = int(user(lmy))
    206:   lambda = user(param)

    208:   hy = one/(my - 1)
    209:   hx = one/(mx - 1)

    211:   PetscCall(VecGetArray(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 == 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: 20    continue
    223: 10    continue
    224:       PetscCall(VecRestoreArray(X, xx, ierr))
    225:     end
    226: !
    227: !  --------------------  Evaluate Function F(x) ---------------------
    228: !
    229:     subroutine FormFunction(ts, t, X, F, user, ierr)
    230:       use petscts
    231:       implicit none

    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)

    247:       two = 2.0

    249:       mx = int(user(lmx))
    250:       my = int(user(lmy))
    251:       lambda = user(param)

    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

    259:       PetscCall(VecGetArrayRead(X, xx, ierr))
    260:       PetscCall(VecGetArray(F, ff, ierr))
    261:       do 10 j = 1, my
    262:         do 20 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: 20        continue
    277: 10        continue

    279:           PetscCall(VecRestoreArrayRead(X, xx, ierr))
    280:           PetscCall(VecRestoreArray(F, ff, ierr))
    281:         end
    282: !
    283: !  --------------------  Evaluate Jacobian of F(x) --------------------
    284: !
    285:         subroutine FormJacobian(ts, ctime, X, JJ, B, user, ierr)
    286:           use petscts
    287:           implicit none

    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

    302:           PetscInt param, lmx, lmy
    303:           parameter(param=1, lmx=2, lmy=3)

    305:           i1 = 1
    306:           i5 = 5
    307:           jac = B
    308:           two = 2.0
    309:           one = 1.0

    311:           mx = int(user(lmx))
    312:           my = int(user(lmy))
    313:           lambda = user(param)

    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

    321:           PetscCall(VecGetArrayRead(X, xx, ierr))
    322:           do 10 j = 1, my
    323:             do 20 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: 20            continue
    344: 10            continue
    345:               PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
    346:               PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
    347:               PetscCall(VecRestoreArray(X, xx, ierr))
    348:             end

    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*/