Actual source code: ex22f.F90

      1: ! Time-dependent advection-reaction PDE in 1d. Demonstrates IMEX methods
      2: !
      3: !   u_t + a1*u_x = -k1*u + k2*v + s1
      4: !   v_t + a2*v_x = k1*u - k2*v + s2
      5: !   0 < x < 1
      6: !   a1 = 1, k1 = 10^6, s1 = 0,
      7: !   a2 = 0, k2 = 2*k1, s2 = 1
      8: !
      9: !   Initial conditions:
     10: !   u(x,0) = 1 + s2*x
     11: !   v(x,0) = k0/k1*u(x,0) + s1/k1
     12: !
     13: !   Upstream boundary conditions:
     14: !   u(0,t) = 1-sin(12*t)^4
     15: !

     17: program main
     18: #include <petsc/finclude/petscts.h>
     19: #include <petsc/finclude/petscdmda.h>
     20:   use petscts
     21:   implicit none

     23: !
     24: !  Create an application context to contain data needed by the
     25: !  application-provided call-back routines, FormJacobian() and
     26: !  FormFunction(). We use a double precision array with six
     27: !  entries, two for each problem parameter a, k, s.
     28: !
     29:   PetscReal user(6)
     30:   integer user_a, user_k, user_s
     31:   parameter(user_a=0, user_k=2, user_s=4)

     33:   TS ts
     34:   SNES snes
     35:   SNESLineSearch linesearch
     36:   Vec X
     37:   Mat J
     38:   PetscInt mx
     39:   PetscErrorCode ierr
     40:   DM da
     41:   PetscReal ftime, dt
     42:   PetscReal one, pone
     43:   PetscInt im11, i2
     44:   PetscBool flg

     46:   im11 = 11
     47:   i2 = 2
     48:   one = 1.0
     49:   pone = one/10

     51:   PetscCallA(PetscInitialize(ierr))

     53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     54: !  Create distributed array (DMDA) to manage parallel grid and vectors
     55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     56:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, im11, i2, i2, PETSC_NULL_INTEGER, da, ierr))
     57:   PetscCallA(DMSetFromOptions(da, ierr))
     58:   PetscCallA(DMSetUp(da, ierr))

     60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     61: !    Extract global vectors from DMDA
     62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     63:   PetscCallA(DMCreateGlobalVector(da, X, ierr))

     65: ! Initialize user application context
     66: ! Use zero-based indexing for command line parameters to match ex22.c
     67:   user(user_a + 1) = 1.0
     68:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', user(user_a + 1), flg, ierr))
     69:   user(user_a + 2) = 0.0
     70:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', user(user_a + 2), flg, ierr))
     71:   user(user_k + 1) = 1000000.0
     72:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', user(user_k + 1), flg, ierr))
     73:   user(user_k + 2) = 2*user(user_k + 1)
     74:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', user(user_k + 2), flg, ierr))
     75:   user(user_s + 1) = 0.0
     76:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', user(user_s + 1), flg, ierr))
     77:   user(user_s + 2) = 1.0
     78:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', user(user_s + 2), flg, ierr))

     80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     81: !    Create timestepping solver context
     82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83:   PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
     84:   PetscCallA(TSSetDM(ts, da, ierr))
     85:   PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
     86:   PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, user, ierr))
     87:   PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, user, ierr))
     88:   PetscCallA(DMSetMatType(da, MATAIJ, ierr))
     89:   PetscCallA(DMCreateMatrix(da, J, ierr))
     90:   PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, user, ierr))

     92:   PetscCallA(TSGetSNES(ts, snes, ierr))
     93:   PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
     94:   PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))

     96:   ftime = 1.0
     97:   PetscCallA(TSSetMaxTime(ts, ftime, ierr))
     98:   PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))

    100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101: !  Set initial conditions
    102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    103:   PetscCallA(FormInitialSolution(ts, X, user, ierr))
    104:   PetscCallA(TSSetSolution(ts, X, ierr))
    105:   PetscCallA(VecGetSize(X, mx, ierr))
    106: !  Advective CFL, I don't know why it needs so much safety factor.
    107:   dt = pone*max(user(user_a + 1), user(user_a + 2))/mx
    108:   PetscCallA(TSSetTimeStep(ts, dt, ierr))

    110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    111: !   Set runtime options
    112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    113:   PetscCallA(TSSetFromOptions(ts, ierr))

    115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116: !  Solve nonlinear system
    117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118:   PetscCallA(TSSolve(ts, X, ierr))

    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121: !  Free work space.
    122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:   PetscCallA(MatDestroy(J, ierr))
    124:   PetscCallA(VecDestroy(X, ierr))
    125:   PetscCallA(TSDestroy(ts, ierr))
    126:   PetscCallA(DMDestroy(da, ierr))
    127:   PetscCallA(PetscFinalize(ierr))
    128: contains

    130: ! Small helper to extract the layout, result uses 1-based indexing.
    131:   subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
    132:     use petscdm
    133:     implicit none

    135:     DM da
    136:     PetscInt mx, xs, xe, gxs, gxe
    137:     PetscErrorCode ierr
    138:     PetscInt xm, gxm
    139:     PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_ENUM, PETSC_NULL_ENUM, PETSC_NULL_ENUM, PETSC_NULL_ENUM, ierr))
    140:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    141:     PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
    142:     xs = xs + 1
    143:     gxs = gxs + 1
    144:     xe = xs + xm - 1
    145:     gxe = gxs + gxm - 1
    146:   end subroutine

    148:   subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
    149:     implicit none
    150:     PetscInt mx, xs, xe, gxs, gxe
    151:     PetscScalar x(2, xs:xe)
    152:     PetscScalar xdot(2, xs:xe)
    153:     PetscScalar f(2, xs:xe)
    154:     PetscReal a(2), k(2), s(2)
    155:     PetscErrorCode ierr
    156:     PetscInt i
    157:     do 10 i = xs, xe
    158:       f(1, i) = xdot(1, i) + k(1)*x(1, i) - k(2)*x(2, i) - s(1)
    159:       f(2, i) = xdot(2, i) - k(1)*x(1, i) + k(2)*x(2, i) - s(2)
    160: 10    continue
    161:       end subroutine

    163:       subroutine FormIFunction(ts, t, X, Xdot, F, user, ierr)
    164:         use petscts
    165:         implicit none

    167:         TS ts
    168:         PetscReal t
    169:         Vec X, Xdot, F
    170:         PetscReal user(6)
    171:         PetscErrorCode ierr
    172:         integer user_a, user_k, user_s
    173:         parameter(user_a=1, user_k=3, user_s=5)

    175:         DM da
    176:         PetscInt mx, xs, xe, gxs, gxe
    177:         PetscScalar, pointer :: xx(:), xxdot(:), ff(:)

    179:         PetscCall(TSGetDM(ts, da, ierr))
    180:         PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    182: ! Get access to vector data
    183:         PetscCall(VecGetArrayRead(X, xx, ierr))
    184:         PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
    185:         PetscCall(VecGetArray(F, ff, ierr))

    187:         PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, user(user_a), user(user_k), user(user_s), ierr))

    189:         PetscCall(VecRestoreArrayRead(X, xx, ierr))
    190:         PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
    191:         PetscCall(VecRestoreArray(F, ff, ierr))
    192:       end subroutine

    194:       subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
    195:         implicit none
    196:         PetscInt mx, xs, xe, gxs, gxe
    197:         PetscReal t
    198:         PetscScalar x(2, gxs:gxe), f(2, xs:xe)
    199:         PetscReal a(2), k(2), s(2)
    200:         PetscErrorCode ierr
    201:         PetscInt i, j
    202:         PetscReal hx, u0t(2)
    203:         PetscReal one, two, three, four, six, twelve
    204:         PetscReal half, third, twothird, sixth
    205:         PetscReal twelfth

    207:         one = 1.0
    208:         two = 2.0
    209:         three = 3.0
    210:         four = 4.0
    211:         six = 6.0
    212:         twelve = 12.0
    213:         hx = one/mx
    214: !     The Fortran standard only allows positive base for power functions; Nag compiler fails on this
    215:         u0t(1) = one - abs(sin(twelve*t))**four
    216:         u0t(2) = 0.0
    217:         half = one/two
    218:         third = one/three
    219:         twothird = two/three
    220:         sixth = one/six
    221:         twelfth = one/twelve
    222:         do 20 i = xs, xe
    223:           do 10 j = 1, 2
    224:             if (i == 1) then
    225:               f(j, i) = a(j)/hx*(third*u0t(j) + half*x(j, i) - x(j, i + 1)  &
    226:     &              + sixth*x(j, i + 2))
    227:             else if (i == 2) then
    228:               f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1)    &
    229:     &              - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
    230:             else if (i == mx - 1) then
    231:               f(j, i) = a(j)/hx*(-sixth*x(j, i - 2) + x(j, i - 1)             &
    232:     &         - half*x(j, i) - third*x(j, i + 1))
    233:             else if (i == mx) then
    234:               f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
    235:             else
    236:               f(j, i) = a(j)/hx*(-twelfth*x(j, i - 2)                      &
    237:     &              + twothird*x(j, i - 1)                                 &
    238:     &              - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
    239:             end if
    240: 10          continue
    241: 20          continue
    242:             end subroutine

    244:             subroutine FormRHSFunction(ts, t, X, F, user, ierr)
    245:               use petscts
    246:               implicit none

    248:               TS ts
    249:               PetscReal t
    250:               Vec X, F
    251:               PetscReal user(6)
    252:               PetscErrorCode ierr
    253:               integer user_a, user_k, user_s
    254:               parameter(user_a=1, user_k=3, user_s=5)
    255:               DM da
    256:               Vec Xloc
    257:               PetscInt mx, xs, xe, gxs, gxe
    258:               PetscScalar, pointer :: xx(:), ff(:)

    260:               PetscCall(TSGetDM(ts, da, ierr))
    261:               PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    263: !     Scatter ghost points to local vector,using the 2-step process
    264: !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    265: !     By placing code between these two statements, computations can be
    266: !     done while messages are in transition.
    267:               PetscCall(DMGetLocalVector(da, Xloc, ierr))
    268:               PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
    269:               PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))

    271: ! Get access to vector data
    272:               PetscCall(VecGetArrayRead(Xloc, xx, ierr))
    273:               PetscCall(VecGetArray(F, ff, ierr))

    275:               PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, user(user_a), user(user_k), user(user_s), ierr))

    277:               PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
    278:               PetscCall(VecRestoreArray(F, ff, ierr))
    279:               PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
    280:             end subroutine

    282: ! ---------------------------------------------------------------------
    283: !
    284: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    285: !
    286:             subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, user, ierr)
    287:               use petscts
    288:               implicit none

    290:               TS ts
    291:               PetscReal t, shift
    292:               Vec X, Xdot
    293:               Mat J, Jpre
    294:               PetscReal user(6)
    295:               PetscErrorCode ierr
    296:               integer user_a, user_k, user_s
    297:               parameter(user_a=0, user_k=2, user_s=4)

    299:               DM da
    300:               PetscInt mx, xs, xe, gxs, gxe
    301:               PetscInt i, i1, row, col
    302:               PetscReal k1, k2
    303:               PetscScalar val(4)

    305:               PetscCall(TSGetDM(ts, da, ierr))
    306:               PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    308:               i1 = 1
    309:               k1 = user(user_k + 1)
    310:               k2 = user(user_k + 2)
    311:               do 10 i = xs, xe
    312:                 row = i - gxs
    313:                 col = i - gxs
    314:                 val(1) = shift + k1
    315:                 val(2) = -k2
    316:                 val(3) = -k1
    317:                 val(4) = shift + k2
    318:                 PetscCall(MatSetValuesBlockedLocal(Jpre, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
    319: 10              continue
    320:                 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
    321:                 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
    322:                 if (J /= Jpre) then
    323:                   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
    324:                   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
    325:                 end if
    326:                 end subroutine

    328:                 subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
    329:                   implicit none
    330:                   PetscInt mx, xs, xe, gxs, gxe
    331:                   PetscScalar x(2, xs:xe)
    332:                   PetscReal a(2), k(2), s(2)
    333:                   PetscErrorCode ierr

    335:                   PetscInt i
    336:                   PetscReal one, hx, r, ik
    337:                   one = 1.0
    338:                   hx = one/mx
    339:                   do 10 i = xs, xe
    340:                     r = i*hx
    341:                     if (k(2) /= 0.0) then
    342:                       ik = one/k(2)
    343:                     else
    344:                       ik = one
    345:                     end if
    346:                     x(1, i) = one + s(2)*r
    347:                     x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
    348: 10                  continue
    349:                     end subroutine

    351:                     subroutine FormInitialSolution(ts, X, user, ierr)
    352:                       use petscts
    353:                       implicit none

    355:                       TS ts
    356:                       Vec X
    357:                       PetscReal user(6)
    358:                       PetscErrorCode ierr
    359:                       integer user_a, user_k, user_s
    360:                       parameter(user_a=1, user_k=3, user_s=5)

    362:                       DM da
    363:                       PetscInt mx, xs, xe, gxs, gxe
    364:                       PetscScalar, pointer :: xx(:)

    366:                       PetscCall(TSGetDM(ts, da, ierr))
    367:                       PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    369: ! Get access to vector data
    370:                       PetscCall(VecGetArray(X, xx, ierr))

    372:                       PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, user(user_a), user(user_k), user(user_s), ierr))

    374:                       PetscCall(VecRestoreArray(X, xx, ierr))
    375:                     end subroutine
    376:                     end program
    377: !/*TEST
    378: !
    379: !    test:
    380: !      args: -da_grid_x 200 -ts_arkimex_type 4
    381: !      requires: !single
    382: !      output_file: output/empty.out
    383: !
    384: !TEST*/