Actual source code: ex22f.F
petsc-3.4.2 2013-07-02

      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:       implicit none
     19: #include <finclude/petscsys.h>
     20: #include <finclude/petscvec.h>
     21: #include <finclude/petscmat.h>
     22: #include <finclude/petscsnes.h>
     23: #include <finclude/petscts.h>
     24: #include <finclude/petscdmda.h>
     25: !
     26: !  Create an application context to contain data needed by the
     27: !  application-provided call-back routines, FormJacobian() and
     28: !  FormFunction(). We use a double precision array with six
     29: !  entries, two for each problem parameter a, k, s.
     30: !
     31:       PetscReal user(6)
     32:       integer user_a,user_k,user_s
     33:       parameter (user_a = 0,user_k = 2,user_s = 4)

     35:       external FormRHSFunction,FormIFunction,FormIJacobian
     36:       external FormInitialSolution

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

     51:       im11 = -11
     52:       i2   = 2
     53:       zero = 0.0
     54:       one = 1.0
     55:       pone = one / 10

     57:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

     59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60: !  Create distributed array (DMDA) to manage parallel grid and vectors
     61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     62:       call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,im11,i2,i2, &
     63:      &     PETSC_NULL_INTEGER,da,ierr)

     65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66: !    Extract global vectors from DMDA;
     67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     68:       call DMCreateGlobalVector(da,X,ierr)

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

     91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     92: !    Create timestepping solver context
     93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     94:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
     95:       call TSSetDM(ts,da,ierr)
     96:       call TSSetType(ts,TSARKIMEX,ierr)
     97:       call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction,       &
     98:      &     user,ierr)
     99:       call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
    100:       call DMCreateMatrix(da,MATAIJ,J,ierr)
    101:       call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)

    103:       call TSGetSNES(ts,snes,ierr)
    104:       call SNESGetLineSearch(snes,linesearch,ierr)
    105:       call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)

    107:       ftime = 1.0
    108:       maxsteps = 10000
    109:       call TSSetDuration(ts,maxsteps,ftime,ierr)

    111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    112: !  Set initial conditions
    113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114:       call FormInitialSolution(ts,X,user,ierr)
    115:       call TSSetSolution(ts,X,ierr)
    116:       call VecGetSize(X,mx,ierr)
    117: !  Advective CFL, I don't know why it needs so much safety factor.
    118:       dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
    119:       call TSSetInitialTimeStep(ts,zero,dt,ierr)

    121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    122: !   Set runtime options
    123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    124:       call TSSetFromOptions(ts,ierr)

    126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127: !  Solve nonlinear system
    128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    129:       call TSSolve(ts,X,ierr)

    131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    132: !  Free work space.
    133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    134:       call MatDestroy(J,ierr)
    135:       call VecDestroy(X,ierr)
    136:       call TSDestroy(ts,ierr)
    137:       call DMDestroy(da,ierr)
    138:       call PetscFinalize(ierr)
    139:       end program

    141: ! Small helper to extract the layout, result uses 1-based indexing.
    142:       subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
    143:       implicit none
    144: #include <finclude/petscsys.h>
    145: #include <finclude/petscdmda.h>
    146:       DM da
    147:       PetscInt mx,xs,xe,gxs,gxe
    148:       PetscErrorCode ierr
    149:       PetscInt xm,gxm
    150:             call DMDAGetInfo(da,PETSC_NULL_INTEGER,                     &
    151:      &     mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                    &
    152:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
    153:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                       &
    154:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
    155:      &     PETSC_NULL_INTEGER,ierr)
    156:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,  &
    157:      &     xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
    158:       call DMDAGetGhostCorners(da,                                      &
    159:      &     gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                   &
    160:      &     gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
    161:       xs = xs + 1
    162:       gxs = gxs + 1
    163:       xe = xs + xm - 1
    164:       gxe = gxs + gxm - 1
    165:       end subroutine

    167:       subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,          &
    168:      &     a,k,s,ierr)
    169:       implicit none
    170:       PetscInt mx,xs,xe,gxs,gxe
    171:       PetscScalar x(2,xs:xe)
    172:       PetscScalar xdot(2,xs:xe)
    173:       PetscScalar f(2,xs:xe)
    174:       PetscReal a(2),k(2),s(2)
    175:       PetscErrorCode ierr
    176:       PetscInt i
    177:       do 10 i = xs,xe
    178:          f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
    179:          f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
    180:  10   continue
    181:       end subroutine

    183:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
    184:       implicit none
    185: #include <finclude/petscsys.h>
    186: #include <finclude/petscvec.h>
    187: #include <finclude/petscmat.h>
    188: #include <finclude/petscsnes.h>
    189: #include <finclude/petscts.h>
    190: #include <finclude/petscdmda.h>
    191:       TS ts
    192:       PetscReal t
    193:       Vec X,Xdot,F
    194:       PetscReal user(6)
    195:       PetscErrorCode ierr
    196:       integer user_a,user_k,user_s
    197:       parameter (user_a = 1,user_k = 3,user_s = 5)

    199:       DM             da
    200:       PetscInt       mx,xs,xe,gxs,gxe
    201:       PetscOffset    ixx,ixxdot,iff
    202:       PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

    204:       call TSGetDM(ts,da,ierr)
    205:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

    207: ! Get access to vector data
    208:       call VecGetArray(X,xx,ixx,ierr)
    209:       call VecGetArray(Xdot,xxdot,ixxdot,ierr)
    210:       call VecGetArray(F,ff,iff,ierr)

    212:       call FormIFunctionLocal(mx,xs,xe,gxs,gxe,                         &
    213:      &     xx(ixx),xxdot(ixxdot),ff(iff),                               &
    214:      &     user(user_a),user(user_k),user(user_s),ierr)

    216:       call VecRestoreArray(X,xx,ixx,ierr)
    217:       call VecRestoreArray(Xdot,xxdot,ixxdot,ierr)
    218:       call VecRestoreArray(F,ff,iff,ierr)
    219:       end subroutine

    221:       subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,           &
    222:      &     a,k,s,ierr)
    223:       implicit none
    224:       PetscInt mx,xs,xe,gxs,gxe
    225:       PetscReal t
    226:       PetscScalar x(2,gxs:gxe),f(2,xs:xe)
    227:       PetscReal a(2),k(2),s(2)
    228:       PetscErrorCode ierr
    229:       PetscInt i,j
    230:       PetscReal hx,u0t(2)
    231:       PetscReal one,two,three,four,six,twelve
    232:       PetscReal half,third,twothird,sixth
    233:       PetscReal twelfth

    235:       one = 1.0
    236:       two = 2.0
    237:       three = 3.0
    238:       four = 4.0
    239:       six = 6.0
    240:       twelve = 12.0
    241:       hx = one / mx
    242:       u0t(1) = one - sin(twelve*t)**four
    243:       u0t(2) = 0.0
    244:       half = one/two
    245:       third = one / three
    246:       twothird = two / three
    247:       sixth = one / six
    248:       twelfth = one / twelve
    249:       do 20 i = xs,xe
    250:          do 10 j = 1,2
    251:             if (i .eq. 1) then
    252:                f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
    253:      &              + sixth*x(j,i+2))
    254:             else if (i .eq. 2) then
    255:                f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
    256:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
    257:             else if (i .eq. mx-1) then
    258:                f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1)             &
    259:      &         - half*x(j,i) -third*x(j,i+1))
    260:             else if (i .eq. mx) then
    261:                f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
    262:             else
    263:                f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
    264:      &              + twothird*x(j,i-1)                                 &
    265:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
    266:             end if
    267:  10      continue
    268:  20   continue
    269:       end subroutine

    271:       subroutine FormRHSFunction(ts,t,X,F,user,ierr)
    272:       implicit none
    273: #include <finclude/petscsys.h>
    274: #include <finclude/petscvec.h>
    275: #include <finclude/petscmat.h>
    276: #include <finclude/petscsnes.h>
    277: #include <finclude/petscts.h>
    278: #include <finclude/petscdmda.h>
    279:       TS ts
    280:       PetscReal t
    281:       Vec X,F
    282:       PetscReal user(6)
    283:       PetscErrorCode ierr
    284:       integer user_a,user_k,user_s
    285:       parameter (user_a = 1,user_k = 3,user_s = 5)
    286:       DM             da
    287:       Vec            Xloc
    288:       PetscInt       mx,xs,xe,gxs,gxe
    289:       PetscOffset    ixx,iff
    290:       PetscScalar    xx(0:1),ff(0:1)

    292:       call TSGetDM(ts,da,ierr)
    293:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

    295: !     Scatter ghost points to local vector,using the 2-step process
    296: !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    297: !     By placing code between these two statements, computations can be
    298: !     done while messages are in transition.
    299:       call DMGetLocalVector(da,Xloc,ierr)
    300:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
    301:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)

    303: ! Get access to vector data
    304:       call VecGetArray(Xloc,xx,ixx,ierr)
    305:       call VecGetArray(F,ff,iff,ierr)

    307:       call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,                       &
    308:      &     t,xx(ixx),ff(iff),                                           &
    309:      &     user(user_a),user(user_k),user(user_s),ierr)

    311:       call VecRestoreArray(Xloc,xx,ixx,ierr)
    312:       call VecRestoreArray(F,ff,iff,ierr)
    313:       call DMRestoreLocalVector(da,Xloc,ierr)
    314:       end subroutine

    316: ! ---------------------------------------------------------------------
    317: !
    318: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    319: !
    320:       subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,mstr,user,ierr)
    321:       implicit none
    322: #include <finclude/petscsys.h>
    323: #include <finclude/petscvec.h>
    324: #include <finclude/petscmat.h>
    325: #include <finclude/petscsnes.h>
    326: #include <finclude/petscts.h>
    327: #include <finclude/petscdmda.h>
    328:       TS ts
    329:       PetscReal t,shift
    330:       Vec X,Xdot
    331:       Mat J,Jpre
    332:       MatStructure mstr
    333:       PetscReal user(6)
    334:       PetscErrorCode ierr
    335:       integer user_a,user_k,user_s
    336:       parameter (user_a = 0,user_k = 2,user_s = 4)

    338:       DM             da
    339:       PetscInt       mx,xs,xe,gxs,gxe
    340:       PetscInt       i,i1,row,col
    341:       PetscReal      k1,k2;
    342:       PetscScalar    val(4)

    344:       call TSGetDM(ts,da,ierr)
    345:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

    347:       i1 = 1
    348:       k1 = user(user_k+1)
    349:       k2 = user(user_k+2)
    350:       do 10 i = xs,xe
    351:          row = i-gxs
    352:          col = i-gxs
    353:          val(1) = shift + k1
    354:          val(2) = -k2
    355:          val(3) = -k1
    356:          val(4) = shift + k2
    357:          call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,          &
    358:      &        INSERT_VALUES,ierr)
    359:  10   continue
    360:       call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
    361:       call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
    362:       if (J /= Jpre) then
    363:          call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
    364:          call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
    365:       end if
    366:       mstr = SAME_NONZERO_PATTERN
    367:       end subroutine

    369:       subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
    370:       implicit none
    371:       PetscInt mx,xs,xe,gxs,gxe
    372:       PetscScalar x(2,xs:xe)
    373:       PetscReal a(2),k(2),s(2)
    374:       PetscErrorCode ierr

    376:       PetscInt i
    377:       PetscReal one,hx,r,ik
    378:       one = 1.0
    379:       hx = one / mx
    380:       do 10 i=xs,xe
    381:          r = i*hx
    382:          if (k(2) .ne. 0.0) then
    383:             ik = one/k(2)
    384:          else
    385:             ik = one
    386:          end if
    387:          x(1,i) = one + s(2)*r
    388:          x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
    389:  10   continue
    390:       end subroutine

    392:       subroutine FormInitialSolution(ts,X,user,ierr)
    393:       implicit none
    394: #include <finclude/petscsys.h>
    395: #include <finclude/petscvec.h>
    396: #include <finclude/petscmat.h>
    397: #include <finclude/petscsnes.h>
    398: #include <finclude/petscts.h>
    399: #include <finclude/petscdmda.h>
    400:       TS ts
    401:       Vec X
    402:       PetscReal user(6)
    403:       PetscErrorCode ierr
    404:       integer user_a,user_k,user_s
    405:       parameter (user_a = 1,user_k = 3,user_s = 5)

    407:       DM             da
    408:       PetscInt       mx,xs,xe,gxs,gxe
    409:       PetscOffset    ixx
    410:       PetscScalar    xx(0:1)

    412:       call TSGetDM(ts,da,ierr)
    413:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

    415: ! Get access to vector data
    416:       call VecGetArray(X,xx,ixx,ierr)

    418:       call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,                       &
    419:      &     xx(ixx),user(user_a),user(user_k),user(user_s),ierr)

    421:       call VecRestoreArray(X,xx,ixx,ierr)
    422:       end subroutine