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: #include <petsc/finclude/petscts.h>
     18: #include <petsc/finclude/petscdmda.h>

     20: module ex22f_modctx
     21:   use petscts
     22:   use petscdm
     23:   implicit none
     24:   type AppCtx
     25:     PetscReal a(2), k(2), s(2)
     26:   end type AppCtx
     27: contains

     29: ! Small helper to extract the layout, result uses 1-based indexing.
     30:   subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
     31:     DM da
     32:     PetscInt mx, xs, xe, gxs, gxe
     33:     PetscErrorCode, intent(out) :: ierr
     34:     PetscInt xm, gxm
     35:     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))
     36:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
     37:     PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
     38:     xs = xs + 1
     39:     gxs = gxs + 1
     40:     xe = xs + xm - 1
     41:     gxe = gxs + gxm - 1
     42:   end subroutine

     44:   subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
     45:     PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
     46:     PetscScalar, dimension(2, xs:xe), intent(in) :: x, xdot
     47:     PetscScalar, dimension(2, xs:xe), intent(out) :: f
     48:     PetscReal, dimension(2) :: a, k, s
     49:     PetscErrorCode, intent(out) :: ierr

     51:     f(1, :) = xdot(1, :) + k(1)*x(1, :) - k(2)*x(2, :) - s(1)
     52:     f(2, :) = xdot(2, :) - k(1)*x(1, :) + k(2)*x(2, :) - s(2)
     53:     ierr = 0
     54:   end subroutine

     56:   subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr)
     57:     TS ts
     58:     type(AppCtx) ctx
     59:     PetscReal t
     60:     Vec X, Xdot, F
     61:     PetscErrorCode, intent(out) :: ierr

     63:     DM da
     64:     PetscInt mx, xs, xe, gxs, gxe
     65:     PetscScalar, pointer :: xx(:), xxdot(:), ff(:)

     67:     PetscCall(TSGetDM(ts, da, ierr))
     68:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

     70: ! Get access to vector data
     71:     PetscCall(VecGetArrayRead(X, xx, ierr))
     72:     PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
     73:     PetscCall(VecGetArray(F, ff, ierr))

     75:     PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr))

     77:     PetscCall(VecRestoreArrayRead(X, xx, ierr))
     78:     PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
     79:     PetscCall(VecRestoreArray(F, ff, ierr))
     80:   end subroutine

     82:   subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
     83:     PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
     84:     PetscReal, intent(in) :: t
     85:     PetscScalar x(2, gxs:gxe), f(2, xs:xe)
     86:     PetscReal, dimension(2) :: a, k, s, u0t
     87:     PetscErrorCode ierr
     88:     PetscInt i, j
     89:     PetscReal hx
     90:     PetscReal, parameter :: twelfth = 1._PETSC_REAL_KIND/12._PETSC_REAL_KIND, twothird = 2._PETSC_REAL_KIND/3._PETSC_REAL_KIND

     92:     hx = 1.0_PETSC_REAL_KIND/mx
     93:     ! The Fortran standard only allows positive base for power functions; Nag compiler fails on this
     94:     u0t = [1.0_PETSC_REAL_KIND - abs(sin(12._PETSC_REAL_KIND*t))**4, 0._PETSC_REAL_KIND]
     95:     do i = xs, xe
     96:       do j = 1, 2
     97:         if (i == 1) then
     98:           f(j, i) = a(j)/hx*(u0t(j)/3._PETSC_REAL_KIND + .5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1) + x(j, i + 2)/6._PETSC_REAL_KIND)
     99:         else if (i == 2) then
    100:           f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
    101:         else if (i == mx - 1) then
    102:           f(j, i) = a(j)/hx*(-x(j, i - 2)/6._PETSC_REAL_KIND + x(j, i - 1) - .5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1)/3._PETSC_REAL_KIND)
    103:         else if (i == mx) then
    104:           f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
    105:         else
    106:           f(j, i) = a(j)/hx*(-twelfth*x(j, i - 2) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
    107:         end if
    108:       end do
    109:     end do
    110:   end subroutine

    112:   subroutine FormRHSFunction(ts, t, X, F, ctx, ierr)
    113:     type(AppCtx) ctx
    114:     TS ts
    115:     PetscReal t
    116:     Vec X, F
    117:     PetscErrorCode, intent(out) :: ierr
    118:     DM da
    119:     Vec Xloc
    120:     PetscInt mx, xs, xe, gxs, gxe
    121:     PetscScalar, pointer :: xx(:), ff(:)

    123:     PetscCall(TSGetDM(ts, da, ierr))
    124:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    126: !     Scatter ghost points to local vector,using the 2-step process
    127: !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    128: !     By placing code between these two statements, computations can be
    129: !     done while messages are in transition.
    130:     PetscCall(DMGetLocalVector(da, Xloc, ierr))
    131:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
    132:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))

    134: ! Get access to vector data
    135:     PetscCall(VecGetArrayRead(Xloc, xx, ierr))
    136:     PetscCall(VecGetArray(F, ff, ierr))

    138:     PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr))

    140:     PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
    141:     PetscCall(VecRestoreArray(F, ff, ierr))
    142:     PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
    143:   end subroutine

    145: ! ---------------------------------------------------------------------
    146: !
    147: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    148: !
    149:   subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
    150:     type(AppCtx) ctx
    151:     TS ts
    152:     PetscReal t, shift
    153:     Vec X, Xdot
    154:     Mat J, Jpre
    155:     PetscErrorCode ierr

    157:     DM da
    158:     PetscInt mx, xs, xe, gxs, gxe
    159:     PetscInt i, row, col
    160:     PetscReal k1, k2
    161:     PetscScalar val(4)

    163:     PetscCall(TSGetDM(ts, da, ierr))
    164:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    166:     k1 = ctx%k(1)
    167:     k2 = ctx%k(2)
    168:     do i = xs, xe
    169:       row = i - gxs
    170:       col = i - gxs
    171:       val = [shift + k1, -k2, -k1, shift + k2]
    172:       PetscCall(MatSetValuesBlockedLocal(Jpre, 1_PETSC_INT_KIND, [row], 1_PETSC_INT_KIND, [col], val, INSERT_VALUES, ierr))
    173:     end do
    174:     PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
    175:     PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
    176:     if (J /= Jpre) then
    177:       PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
    178:       PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
    179:     end if
    180:   end subroutine

    182:   subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
    183:     PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
    184:     PetscScalar, intent(out) :: x(2, xs:xe)
    185:     PetscReal, dimension(2) :: a, k, s
    186:     PetscErrorCode, intent(out) :: ierr
    187:     PetscInt i
    188:     PetscReal hx, r, ik

    190:     hx = 1._PETSC_REAL_KIND/mx
    191:     do i = xs, xe
    192:       r = i*hx
    193:       if (k(2) /= 0.0) then
    194:         ik = 1._PETSC_REAL_KIND/k(2)
    195:       else
    196:         ik = 1._PETSC_REAL_KIND
    197:       end if
    198:       x(1, i) = 1._PETSC_REAL_KIND + s(2)*r
    199:       x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
    200:     end do
    201:     ierr = 0
    202:   end subroutine

    204:   subroutine FormInitialSolution(ts, X, ctx, ierr)
    205:     type(AppCtx) ctx
    206:     TS ts
    207:     Vec X
    208:     PetscErrorCode, intent(out) :: ierr
    209:     DM da
    210:     PetscInt mx, xs, xe, gxs, gxe
    211:     PetscScalar, pointer :: xx(:)

    213:     PetscCall(TSGetDM(ts, da, ierr))
    214:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

    216: ! Get access to vector data
    217:     PetscCall(VecGetArray(X, xx, ierr))

    219:     PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr))

    221:     PetscCall(VecRestoreArray(X, xx, ierr))
    222:   end subroutine

    224: end module ex22f_modctx
    225: program main
    226:   use ex22f_modctx
    227:   implicit none

    229: !
    230: !  Create an application context to contain data needed by the
    231: !  application-provided call-back routines, FormJacobian() and
    232: !  FormFunction(). We use a double precision array with six
    233: !  entries, two for each problem parameter a, k, s.
    234: !

    236:   TS ts
    237:   SNES snes
    238:   SNESLineSearch linesearch
    239:   Vec X
    240:   Mat J
    241:   PetscInt mx
    242:   PetscErrorCode ierr
    243:   DM da
    244:   PetscReal, parameter :: ftime = 1.0
    245:   PetscReal dt
    246:   PetscBool flg
    247:   type(AppCtx) ctx

    249:   PetscCallA(PetscInitialize(ierr))

    251: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    252: !  Create distributed array (DMDA) to manage parallel grid and vectors
    253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    254:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11_PETSC_INT_KIND, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER, da, ierr))
    255:   PetscCallA(DMSetFromOptions(da, ierr))
    256:   PetscCallA(DMSetUp(da, ierr))

    258: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    259: !    Extract global vectors from DMDA
    260: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    261:   PetscCallA(DMCreateGlobalVector(da, X, ierr))

    263: ! Initialize user application context
    264: ! Use zero-based indexing for command line parameters to match ex22.c
    265:   ctx%a(1) = 1.0
    266:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr))
    267:   ctx%a(2) = 0.0
    268:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr))
    269:   ctx%k(1) = 1000000.0
    270:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr))
    271:   ctx%k(2) = 2*ctx%k(1)
    272:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr))
    273:   ctx%s(1) = 0.0
    274:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr))
    275:   ctx%s(2) = 1.0
    276:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr))

    278: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    279: !    Create timestepping solver context
    280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    281:   PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
    282:   PetscCallA(TSSetDM(ts, da, ierr))
    283:   PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
    284:   PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr))
    285:   PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr))
    286:   PetscCallA(DMSetMatType(da, MATAIJ, ierr))
    287:   PetscCallA(DMCreateMatrix(da, J, ierr))
    288:   PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr))

    290:   PetscCallA(TSGetSNES(ts, snes, ierr))
    291:   PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
    292:   PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))

    294:   PetscCallA(TSSetMaxTime(ts, ftime, ierr))
    295:   PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))

    297: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    298: !  Set initial conditions
    299: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    300:   PetscCallA(FormInitialSolution(ts, X, ctx, ierr))
    301:   PetscCallA(TSSetSolution(ts, X, ierr))
    302:   PetscCallA(VecGetSize(X, mx, ierr))
    303: !  Advective CFL, I don't know why it needs so much safety factor.
    304:   dt = .1_PETSC_REAL_KIND*max(ctx%a(1), ctx%a(2))/mx
    305:   PetscCallA(TSSetTimeStep(ts, dt, ierr))

    307: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    308: !   Set runtime options
    309: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    310:   PetscCallA(TSSetFromOptions(ts, ierr))

    312: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    313: !  Solve nonlinear system
    314: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    315:   PetscCallA(TSSolve(ts, X, ierr))

    317: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    318: !  Free work space.
    319: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    320:   PetscCallA(MatDestroy(J, ierr))
    321:   PetscCallA(VecDestroy(X, ierr))
    322:   PetscCallA(TSDestroy(ts, ierr))
    323:   PetscCallA(DMDestroy(da, ierr))
    324:   PetscCallA(PetscFinalize(ierr))
    325: end program
    326: !/*TEST
    327: !
    328: !    test:
    329: !      args: -da_grid_x 200 -ts_arkimex_type 4
    330: !      requires: !single
    331: !      output_file: output/empty.out
    332: !
    333: !TEST*/