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:       external FormRHSFunction,FormIFunction,FormIJacobian
     34:       external FormInitialSolution

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

     49:       im11 = 11
     50:       i2   = 2
     51:       one = 1.0
     52:       pone = one / 10

     54:       PetscCallA(PetscInitialize(ierr))

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

     63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     64: !    Extract global vectors from DMDA;
     65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:       PetscCallA(DMCreateGlobalVector(da,X,ierr))

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

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

     95:       PetscCallA(TSGetSNES(ts,snes,ierr))
     96:       PetscCallA(SNESGetLineSearch(snes,linesearch,ierr))
     97:       PetscCallA(SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr))

     99:       ftime = 1.0
    100:       PetscCallA(TSSetMaxTime(ts,ftime,ierr))
    101:       PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))

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

    113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114: !   Set runtime options
    115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:       PetscCallA(TSSetFromOptions(ts,ierr))

    118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    119: !  Solve nonlinear system
    120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121:       PetscCallA(TSSolve(ts,X,ierr))

    123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    124: !  Free work space.
    125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    126:       PetscCallA(MatDestroy(J,ierr))
    127:       PetscCallA(VecDestroy(X,ierr))
    128:       PetscCallA(TSDestroy(ts,ierr))
    129:       PetscCallA(DMDestroy(da,ierr))
    130:       PetscCallA(PetscFinalize(ierr))
    131:       end program

    133: ! Small helper to extract the layout, result uses 1-based indexing.
    134:       subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
    135:       use petscdmda
    136:       implicit none

    138:       DM da
    139:       PetscInt mx,xs,xe,gxs,gxe
    140:       PetscErrorCode ierr
    141:       PetscInt xm,gxm
    142:       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))
    143:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
    144:       PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
    145:       xs = xs + 1
    146:       gxs = gxs + 1
    147:       xe = xs + xm - 1
    148:       gxe = gxs + gxm - 1
    149:       end subroutine

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

    166:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
    167:       use petscts
    168:       implicit none

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

    178:       DM             da
    179:       PetscInt       mx,xs,xe,gxs,gxe
    180:       PetscScalar,pointer :: xx(:),xxdot(:),ff(:)

    182:       PetscCall(TSGetDM(ts,da,ierr))
    183:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    185: ! Get access to vector data
    186:       PetscCall(VecGetArrayReadF90(X,xx,ierr))
    187:       PetscCall(VecGetArrayReadF90(Xdot,xxdot,ierr))
    188:       PetscCall(VecGetArrayF90(F,ff,ierr))

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

    192:       PetscCall(VecRestoreArrayReadF90(X,xx,ierr))
    193:       PetscCall(VecRestoreArrayReadF90(Xdot,xxdot,ierr))
    194:       PetscCall(VecRestoreArrayF90(F,ff,ierr))
    195:       end subroutine

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

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

    247:       subroutine FormRHSFunction(ts,t,X,F,user,ierr)
    248:       use petscts
    249:       implicit none

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

    263:       PetscCall(TSGetDM(ts,da,ierr))
    264:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

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

    274: ! Get access to vector data
    275:       PetscCall(VecGetArrayReadF90(Xloc,xx,ierr))
    276:       PetscCall(VecGetArrayF90(F,ff,ierr))

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

    280:       PetscCall(VecRestoreArrayReadF90(Xloc,xx,ierr))
    281:       PetscCall(VecRestoreArrayF90(F,ff,ierr))
    282:       PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
    283:       end subroutine

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

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

    302:       DM             da
    303:       PetscInt       mx,xs,xe,gxs,gxe
    304:       PetscInt       i,i1,row,col
    305:       PetscReal      k1,k2;
    306:       PetscScalar    val(4)

    308:       PetscCall(TSGetDM(ts,da,ierr))
    309:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

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

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

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

    354:       subroutine FormInitialSolution(ts,X,user,ierr)
    355:       use petscts
    356:       implicit none

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

    365:       DM             da
    366:       PetscInt       mx,xs,xe,gxs,gxe
    367:       PetscScalar,pointer :: xx(:)

    369:       PetscCall(TSGetDM(ts,da,ierr))
    370:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    372: ! Get access to vector data
    373:       PetscCall(VecGetArrayF90(X,xx,ierr))

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

    377:       PetscCall(VecRestoreArrayF90(X,xx,ierr))
    378:       end subroutine

    380: !/*TEST
    381: !
    382: !    test:
    383: !      args: -da_grid_x 200 -ts_arkimex_type 4
    384: !      requires: !single
    385: !
    386: !TEST*/