Actual source code: ex22f_mf.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:   module PETScShiftMod
     18: #include <petsc/finclude/petscts.h>
     19:     use petscts
     20:     PetscScalar::PETSC_SHIFT
     21:     TS::tscontext
     22:     Mat::Jmat
     23:     PetscReal::MFuser(6)
     24:   end module PETScShiftMod

     27: program main
     28:   use PETScShiftMod
     29:   use petscdmda
     30:   implicit none

     32:   !
     33:   !     Create an application context to contain data needed by the
     34:   !     application-provided call-back routines, FormJacobian() and
     35:   !     FormFunction(). We use a double precision array with six
     36:   !     entries, two for each problem parameter a, k, s.
     37:   !
     38:   PetscReal user(6)
     39:   integer user_a,user_k,user_s
     40:   parameter (user_a = 0,user_k = 2,user_s = 4)

     42:   external FormRHSFunction,FormIFunction
     43:   external FormInitialSolution
     44:   external FormIJacobian
     45:   external MyMult,FormIJacobianMF

     47:   TS             ts
     48:   Vec            X
     49:   Mat            J
     50:   PetscInt       mx
     51:   PetscBool      OptionSaveToDisk
     52:   PetscErrorCode ierr
     53:   DM             da
     54:   PetscReal      ftime,dt
     55:   PetscReal      one,pone
     56:   PetscInt       im11,i2
     57:   PetscBool      flg

     60:   PetscInt       xs,xe,gxs,gxe,dof,gdof
     61:   PetscScalar    shell_shift
     62:   Mat            A

     64:   im11 = 11
     65:   i2   = 2
     66:   one = 1.0
     67:   pone = one / 10

     69:   call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     70:   if (ierr .ne. 0) then
     71:     print*,'PetscInitialize failed'
     72:     stop
     73:   endif

     75:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:   !  Create distributed array (DMDA) to manage parallel grid and vectors
     77:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     78:   call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr)
     79:   call DMSetFromOptions(da,ierr);CHKERRA(ierr)
     80:   call DMSetUp(da,ierr);CHKERRA(ierr)

     82:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83:   !    Extract global vectors from DMDA;
     84:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     85:   call DMCreateGlobalVector(da,X,ierr);CHKERRA(ierr)

     87:   ! Initialize user application context
     88:   ! Use zero-based indexing for command line parameters to match ex22.c
     89:   user(user_a+1) = 1.0
     90:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr);CHKERRA(ierr)
     91:   user(user_a+2) = 0.0
     92:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr);CHKERRA(ierr)
     93:   user(user_k+1) = 1000000.0
     94:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0', user(user_k+1),flg,ierr);CHKERRA(ierr)
     95:   user(user_k+2) = 2*user(user_k+1)
     96:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr);CHKERRA(ierr)
     97:   user(user_s+1) = 0.0
     98:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr);CHKERRA(ierr)
     99:   user(user_s+2) = 1.0
    100:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr);CHKERRA(ierr)

    102:   OptionSaveToDisk=.FALSE.
    103:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sdisk',OptionSaveToDisk,flg,ierr);CHKERRA(ierr)
    104:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    105:   !    Create timestepping solver context
    106:   !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    107:   call TSCreate(PETSC_COMM_WORLD,ts,ierr);CHKERRA(ierr)
    108:   tscontext=ts
    109:   call TSSetDM(ts,da,ierr);CHKERRA(ierr)
    110:   call TSSetType(ts,TSARKIMEX,ierr);CHKERRA(ierr)
    111:   call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr);CHKERRA(ierr)

    113:   ! - - - - - - - - -- - - - -
    114:   !   Matrix free setup
    115:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRA(ierr)
    116:   dof=i2*(xe-xs+1)
    117:   gdof=i2*(gxe-gxs+1)
    118:   call MatCreateShell(PETSC_COMM_WORLD,dof,dof,gdof,gdof,shell_shift,A,ierr);CHKERRA(ierr)

    120:   call MatShellSetOperation(A,MATOP_MULT,MyMult,ierr);CHKERRA(ierr)
    121:   ! - - - - - - - - - - - -

    123:   call TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr);CHKERRA(ierr)
    124:   call DMSetMatType(da,MATAIJ,ierr);CHKERRA(ierr)
    125:   call DMCreateMatrix(da,J,ierr);CHKERRA(ierr)

    127:   Jmat=J

    129:   call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr);CHKERRA(ierr)
    130:   call TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr);CHKERRA(ierr)

    133:   ftime = 1.0
    134:   call TSSetMaxTime(ts,ftime,ierr);CHKERRA(ierr)
    135:   call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr);CHKERRA(ierr)

    137:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    138:   !  Set initial conditions
    139:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    140:   call FormInitialSolution(ts,X,user,ierr);CHKERRA(ierr)
    141:   call TSSetSolution(ts,X,ierr);CHKERRA(ierr)
    142:   call VecGetSize(X,mx,ierr);CHKERRA(ierr)
    143:   !  Advective CFL, I don't know why it needs so much safety factor.
    144:   dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
    145:   call TSSetTimeStep(ts,dt,ierr);CHKERRA(ierr)

    147:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:   !   Set runtime options
    149:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    150:   call TSSetFromOptions(ts,ierr);CHKERRA(ierr)

    152:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    153:   !  Solve nonlinear system
    154:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    155:   call TSSolve(ts,X,ierr);CHKERRA(ierr)

    157:   if (OptionSaveToDisk) then
    158:      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRA(ierr)
    159:      dof=i2*(xe-xs+1)
    160:      gdof=i2*(gxe-gxs+1)
    161:      call SaveSolutionToDisk(da,X,gdof,xs,xe);CHKERRA(ierr)
    162:   end if

    164:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    165:   !  Free work space.
    166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    167:   call MatDestroy(A,ierr);CHKERRA(ierr)
    168:   call MatDestroy(J,ierr);CHKERRA(ierr)
    169:   call VecDestroy(X,ierr);CHKERRA(ierr)
    170:   call TSDestroy(ts,ierr);CHKERRA(ierr)
    171:   call DMDestroy(da,ierr);CHKERRA(ierr)
    172:   call PetscFinalize(ierr)
    173: end program main

    175: ! Small helper to extract the layout, result uses 1-based indexing.
    176:   subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
    177:   use petscdmda
    178:   implicit none

    180:   DM da
    181:   PetscInt mx,xs,xe,gxs,gxe
    182:   PetscErrorCode ierr
    183:   PetscInt xm,gxm
    184:   call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
    185:        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
    186:   call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
    187:   call DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
    188:   xs = xs + 1
    189:   gxs = gxs + 1
    190:   xe = xs + xm - 1
    191:   gxe = gxs + gxm - 1
    192: end subroutine GetLayout

    194: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
    195:   implicit none
    196:   PetscInt mx,xs,xe,gxs,gxe
    197:   PetscScalar x(2,xs:xe)
    198:   PetscScalar xdot(2,xs:xe)
    199:   PetscScalar f(2,xs:xe)
    200:   PetscReal a(2),k(2),s(2)
    201:   PetscErrorCode ierr
    202:   PetscInt i
    203:   do  i = xs,xe
    204:      f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
    205:      f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
    206:   end do
    207: end subroutine FormIFunctionLocal

    209: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
    210:   use petscdmda
    211:   use petscts
    212:   implicit none

    214:   TS ts
    215:   PetscReal t
    216:   Vec X,Xdot,F
    217:   PetscReal user(6)
    218:   PetscErrorCode ierr
    219:   integer user_a,user_k,user_s
    220:   parameter (user_a = 1,user_k = 3,user_s = 5)

    222:   DM             da
    223:   PetscInt       mx,xs,xe,gxs,gxe
    224:   PetscOffset    ixx,ixxdot,iff
    225:   PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

    227:   call TSGetDM(ts,da,ierr);CHKERRQ(ierr)
    228:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRQ(ierr)

    230:   ! Get access to vector data
    231:   call VecGetArrayRead(X,xx,ixx,ierr);CHKERRQ(ierr)
    232:   call VecGetArrayRead(Xdot,xxdot,ixxdot,ierr);CHKERRQ(ierr)
    233:   call VecGetArray(F,ff,iff,ierr);CHKERRQ(ierr)

    235:   call FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr);CHKERRQ(ierr)

    237:   call VecRestoreArrayRead(X,xx,ixx,ierr);CHKERRQ(ierr)
    238:   call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr);CHKERRQ(ierr)
    239:   call VecRestoreArray(F,ff,iff,ierr);CHKERRQ(ierr)
    240: end subroutine FormIFunction

    242: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
    243:   implicit none
    244:   PetscInt mx,xs,xe,gxs,gxe
    245:   PetscReal t
    246:   PetscScalar x(2,gxs:gxe),f(2,xs:xe)
    247:   PetscReal a(2),k(2),s(2)
    248:   PetscErrorCode ierr
    249:   PetscInt i,j
    250:   PetscReal hx,u0t(2)
    251:   PetscReal one,two,three,four,six,twelve
    252:   PetscReal half,third,twothird,sixth
    253:   PetscReal twelfth

    255:   one = 1.0
    256:   two = 2.0
    257:   three = 3.0
    258:   four = 4.0
    259:   six = 6.0
    260:   twelve = 12.0
    261:   hx = one / mx
    262:   u0t(1) = one - sin(twelve*t)**four
    263:   u0t(2) = 0.0
    264:   half = one/two
    265:   third = one / three
    266:   twothird = two / three
    267:   sixth = one / six
    268:   twelfth = one / twelve
    269:   do  i = xs,xe
    270:      do  j = 1,2
    271:         if (i .eq. 1) then
    272:            f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) + sixth*x(j,i+2))
    273:         else if (i .eq. 2) then
    274:            f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2))
    275:         else if (i .eq. mx-1) then
    276:            f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) - half*x(j,i) -third*x(j,i+1))
    277:         else if (i .eq. mx) then
    278:            f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
    279:         else
    280:            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))
    281:         end if
    282:      end do
    283:   end do

    285: #ifdef EXPLICIT_INTEGRATOR22
    286:   do  i = xs,xe
    287:      f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1))
    288:      f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2))
    289:   end do
    290: #endif

    292: end subroutine FormRHSFunctionLocal

    294: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
    295:   use petscts
    296:   use petscdmda
    297:   implicit none

    299:   TS ts
    300:   PetscReal t
    301:   Vec X,F
    302:   PetscReal user(6)
    303:   PetscErrorCode ierr
    304:   integer user_a,user_k,user_s
    305:   parameter (user_a = 1,user_k = 3,user_s = 5)
    306:   DM             da
    307:   Vec            Xloc
    308:   PetscInt       mx,xs,xe,gxs,gxe
    309:   PetscOffset    ixx,iff
    310:   PetscScalar    xx(0:1),ff(0:1)

    312:   call TSGetDM(ts,da,ierr);CHKERRQ(ierr)
    313:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRQ(ierr)

    315:   !     Scatter ghost points to local vector,using the 2-step process
    316:   !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    317:   !     By placing code between these two statements, computations can be
    318:   !     done while messages are in transition.
    319:   call DMGetLocalVector(da,Xloc,ierr);CHKERRQ(ierr)
    320:   call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr);CHKERRQ(ierr)
    321:   call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr);CHKERRQ(ierr)

    323:   ! Get access to vector data
    324:   call VecGetArrayRead(Xloc,xx,ixx,ierr);CHKERRQ(ierr)
    325:   call VecGetArray(F,ff,iff,ierr);CHKERRQ(ierr)

    327:   call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx(ixx),ff(iff),user(user_a),user(user_k),user(user_s),ierr);CHKERRQ(ierr)

    329:   call VecRestoreArrayRead(Xloc,xx,ixx,ierr);CHKERRQ(ierr)
    330:   call VecRestoreArray(F,ff,iff,ierr);CHKERRQ(ierr)
    331:   call DMRestoreLocalVector(da,Xloc,ierr);CHKERRQ(ierr)
    332: end subroutine FormRHSFunction

    334: ! ---------------------------------------------------------------------
    335: !
    336: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    337: !
    338: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
    339:   use petscts
    340:   use petscdmda
    341:   implicit none

    343:   TS ts
    344:   PetscReal t,shift
    345:   Vec X,Xdot
    346:   Mat J,Jpre
    347:   PetscReal user(6)
    348:   PetscErrorCode ierr
    349:   integer user_a,user_k,user_s
    350:   parameter (user_a = 0,user_k = 2,user_s = 4)

    352:   DM             da
    353:   PetscInt       mx,xs,xe,gxs,gxe
    354:   PetscInt       i,i1,row,col
    355:   PetscReal      k1,k2;
    356:   PetscScalar    val(4)

    358:   call TSGetDM(ts,da,ierr);CHKERRQ(ierr)
    359:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRQ(ierr)

    361:   i1 = 1
    362:   k1 = user(user_k+1)
    363:   k2 = user(user_k+2)
    364:   do i = xs,xe
    365:      row = i-gxs
    366:      col = i-gxs
    367:      val(1) = shift + k1
    368:      val(2) = -k2
    369:      val(3) = -k1
    370:      val(4) = shift + k2
    371:      call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,INSERT_VALUES,ierr);CHKERRQ(ierr)
    372:   end do
    373:   call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    374:   call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    375:   if (J /= Jpre) then
    376:      call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    377:      call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
    378:   end if
    379: end subroutine FormIJacobian

    381: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
    382:   implicit none
    383:   PetscInt mx,xs,xe,gxs,gxe
    384:   PetscScalar x(2,xs:xe)
    385:   PetscReal a(2),k(2),s(2)
    386:   PetscErrorCode ierr

    388:   PetscInt i
    389:   PetscReal one,hx,r,ik
    390:   one = 1.0
    391:   hx = one / mx
    392:   do i=xs,xe
    393:      r = i*hx
    394:      if (k(2) .ne. 0.0) then
    395:         ik = one/k(2)
    396:      else
    397:         ik = one
    398:      end if
    399:      x(1,i) = one + s(2)*r
    400:      x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
    401:   end do
    402: end subroutine FormInitialSolutionLocal

    404: subroutine FormInitialSolution(ts,X,user,ierr)
    405:   use petscts
    406:   use petscdmda
    407:   implicit none

    409:   TS ts
    410:   Vec X
    411:   PetscReal user(6)
    412:   PetscErrorCode ierr
    413:   integer user_a,user_k,user_s
    414:   parameter (user_a = 1,user_k = 3,user_s = 5)

    416:   DM             da
    417:   PetscInt       mx,xs,xe,gxs,gxe
    418:   PetscOffset    ixx
    419:   PetscScalar    xx(0:1)

    421:   call TSGetDM(ts,da,ierr);CHKERRQ(ierr)
    422:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr);CHKERRQ(ierr)

    424:   ! Get access to vector data
    425:   call VecGetArray(X,xx,ixx,ierr);CHKERRQ(ierr)

    427:   call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx(ixx),user(user_a),user(user_k),user(user_s),ierr);CHKERRQ(ierr)

    429:   call VecRestoreArray(X,xx,ixx,ierr);CHKERRQ(ierr)
    430: end subroutine FormInitialSolution

    433: ! ---------------------------------------------------------------------
    434: !
    435: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    436: !
    437: subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
    438:   use PETScShiftMod
    439:   implicit none
    440:   TS ts
    441:   PetscReal t,shift
    442:   Vec X,Xdot
    443:   Mat J,Jpre
    444:   PetscReal user(6)
    445:   PetscErrorCode ierr

    447:   !  call MatShellSetContext(J,shift,ierr)
    448:   PETSC_SHIFT=shift
    449:   MFuser=user

    451: end subroutine FormIJacobianMF

    453: ! -------------------------------------------------------------------
    454: !
    455: !   MyMult - user provided matrix multiply
    456: !
    457: !   Input Parameters:
    458: !.  X - input vector
    459: !
    460: !   Output Parameter:
    461: !.  F - function vector
    462: !
    463: subroutine  MyMult(A,X,F,ierr)
    464:   use PETScShiftMod
    465:   implicit none

    467:   Mat     A
    468:   Vec     X,F

    470:   PetscErrorCode ierr
    471:   PetscScalar shift

    473: !  Mat J,Jpre

    475:   PetscReal user(6)

    477:   integer user_a,user_k,user_s
    478:   parameter (user_a = 0,user_k = 2,user_s = 4)

    480:   DM             da
    481:   PetscInt       mx,xs,xe,gxs,gxe
    482:   PetscInt       i,i1,row,col
    483:   PetscReal      k1,k2;
    484:   PetscScalar    val(4)

    486:   !call MatShellGetContext(A,shift,ierr)
    487:   shift=PETSC_SHIFT
    488:   user=MFuser

    490:   call TSGetDM(tscontext,da,ierr)
    491:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

    493:   i1 = 1
    494:   k1 = user(user_k+1)
    495:   k2 = user(user_k+2)

    497:   do i = xs,xe
    498:      row = i-gxs
    499:      col = i-gxs
    500:      val(1) = shift + k1
    501:      val(2) = -k2
    502:      val(3) = -k1
    503:      val(4) = shift + k2
    504:      call MatSetValuesBlockedLocal(Jmat,i1,row,i1,col,val,INSERT_VALUES,ierr)
    505:   end do

    507: !  call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
    508: !  call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
    509: !  if (J /= Jpre) then
    510:      call MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr)
    511:      call MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr)
    512: !  end if

    514:   call MatMult(Jmat,X,F,ierr)

    516:   return
    517: end subroutine MyMult

    520: !
    521: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
    522:   use petscdmda
    523:   implicit none

    525:   Vec X
    526:   DM             da
    527:   PetscInt xs,xe,two
    528:   PetscInt gdof,i
    529:   PetscErrorCode ierr
    530:   PetscOffset    ixx
    531:   PetscScalar data2(2,xs:xe),data(gdof)
    532:   PetscScalar    xx(0:1)

    534:   call VecGetArrayRead(X,xx,ixx,ierr)

    536:   two = 2
    537:   data2=reshape(xx(ixx:ixx+gdof),(/two,xe-xs+1/))
    538:   data=reshape(data2,(/gdof/))
    539:   open(1020,file='solution_out_ex22f_mf.txt')
    540:   do i=1,gdof
    541:      write(1020,'(e24.16,1x)') data(i)
    542:   end do
    543:   close(1020)

    545:   call VecRestoreArrayRead(X,xx,ixx,ierr)
    546: end subroutine SaveSolutionToDisk

    548: !/*TEST
    549: !
    550: !    test:
    551: !      args: -da_grid_x 200 -ts_arkimex_type 4
    552: !      requires: !single
    553: !      output_file: output/ex22f_mf_1.out
    554: !
    555: !TEST*/