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 ex22f_mfmodule
     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 ex22f_mfmodule

     26: program main
     27:   use ex22f_mfmodule
     28:   use petscdmda
     29:   implicit none

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

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

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

     58:   PetscInt       xs,xe,gxs,gxe,dof,gdof
     59:   PetscScalar    shell_shift
     60:   Mat            A

     62:   im11 = 11
     63:   i2   = 2
     64:   one = 1.0
     65:   pone = one / 10

     67:   PetscCallA(PetscInitialize(ierr))

     69:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     70:   !  Create distributed array (DMDA) to manage parallel grid and vectors
     71:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     72:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER_ARRAY,da,ierr))
     73:   PetscCallA(DMSetFromOptions(da,ierr))
     74:   PetscCallA(DMSetUp(da,ierr))

     76:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     77:   !    Extract global vectors from DMDA;
     78:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     79:   PetscCallA(DMCreateGlobalVector(da,X,ierr))

     81:   ! Initialize user application context
     82:   ! Use zero-based indexing for command line parameters to match ex22.c
     83:   user(user_a+1) = 1.0
     84:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr))
     85:   user(user_a+2) = 0.0
     86:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr))
     87:   user(user_k+1) = 1000000.0
     88:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0', user(user_k+1),flg,ierr))
     89:   user(user_k+2) = 2*user(user_k+1)
     90:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr))
     91:   user(user_s+1) = 0.0
     92:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr))
     93:   user(user_s+2) = 1.0
     94:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr))

     96:   OptionSaveToDisk=.FALSE.
     97:       PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sdisk',OptionSaveToDisk,flg,ierr))
     98:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     99:   !    Create timestepping solver context
    100:   !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101:   PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
    102:   tscontext=ts
    103:   PetscCallA(TSSetDM(ts,da,ierr))
    104:   PetscCallA(TSSetType(ts,TSARKIMEX,ierr))
    105:   PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))

    107:   ! - - - - - - - - -- - - - -
    108:   !   Matrix free setup
    109:   PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
    110:   dof=i2*(xe-xs+1)
    111:   gdof=i2*(gxe-gxs+1)
    112:   PetscCallA(MatCreateShell(PETSC_COMM_WORLD,dof,dof,gdof,gdof,shell_shift,A,ierr))

    114:   PetscCallA(MatShellSetOperation(A,MATOP_MULT,MyMult,ierr))
    115:   ! - - - - - - - - - - - -

    117:   PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr))
    118:   PetscCallA(DMSetMatType(da,MATAIJ,ierr))
    119:   PetscCallA(DMCreateMatrix(da,J,ierr))

    121:   Jmat=J

    123:   PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr))
    124:   PetscCallA(TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr))

    126:   ftime = 1.0
    127:   PetscCallA(TSSetMaxTime(ts,ftime,ierr))
    128:   PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))

    130:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131:   !  Set initial conditions
    132:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    133:   PetscCallA(FormInitialSolution(ts,X,user,ierr))
    134:   PetscCallA(TSSetSolution(ts,X,ierr))
    135:   PetscCallA(VecGetSize(X,mx,ierr))
    136:   !  Advective CFL, I don't know why it needs so much safety factor.
    137:   dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
    138:   PetscCallA(TSSetTimeStep(ts,dt,ierr))

    140:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141:   !   Set runtime options
    142:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    143:   PetscCallA(TSSetFromOptions(ts,ierr))

    145:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    146:   !  Solve nonlinear system
    147:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:   PetscCallA(TSSolve(ts,X,ierr))

    150:   if (OptionSaveToDisk) then
    151:      PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
    152:      dof=i2*(xe-xs+1)
    153:      gdof=i2*(gxe-gxs+1)
    154:      call SaveSolutionToDisk(da,X,gdof,xs,xe)
    155:   end if

    157:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158:   !  Free work space.
    159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    160:   PetscCallA(MatDestroy(A,ierr))
    161:   PetscCallA(MatDestroy(J,ierr))
    162:   PetscCallA(VecDestroy(X,ierr))
    163:   PetscCallA(TSDestroy(ts,ierr))
    164:   PetscCallA(DMDestroy(da,ierr))
    165:   PetscCallA(PetscFinalize(ierr))
    166: end program main

    168: ! Small helper to extract the layout, result uses 1-based indexing.
    169:   subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
    170:   use petscdmda
    171:   implicit none

    173:   DM da
    174:   PetscInt mx,xs,xe,gxs,gxe
    175:   PetscErrorCode ierr
    176:   PetscInt xm,gxm
    177:   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))
    178:   PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
    179:   PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
    180:   xs = xs + 1
    181:   gxs = gxs + 1
    182:   xe = xs + xm - 1
    183:   gxe = gxs + gxm - 1
    184: end subroutine GetLayout

    186: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
    187:   implicit none
    188:   PetscInt mx,xs,xe,gxs,gxe
    189:   PetscScalar x(2,xs:xe)
    190:   PetscScalar xdot(2,xs:xe)
    191:   PetscScalar f(2,xs:xe)
    192:   PetscReal a(2),k(2),s(2)
    193:   PetscErrorCode ierr
    194:   PetscInt i
    195:   do  i = xs,xe
    196:      f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
    197:      f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
    198:   end do
    199: end subroutine FormIFunctionLocal

    201: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
    202:   use petscdmda
    203:   use petscts
    204:   implicit none

    206:   TS ts
    207:   PetscReal t
    208:   Vec X,Xdot,F
    209:   PetscReal user(6)
    210:   PetscErrorCode ierr
    211:   integer user_a,user_k,user_s
    212:   parameter (user_a = 1,user_k = 3,user_s = 5)

    214:   DM             da
    215:   PetscInt       mx,xs,xe,gxs,gxe
    216:   PetscScalar,pointer :: xx(:),xxdot(:),ff(:)

    218:   PetscCall(TSGetDM(ts,da,ierr))
    219:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    221:   ! Get access to vector data
    222:   PetscCall(VecGetArrayReadF90(X,xx,ierr))
    223:   PetscCall(VecGetArrayReadF90(Xdot,xxdot,ierr))
    224:   PetscCall(VecGetArrayF90(F,ff,ierr))

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

    228:   PetscCall(VecRestoreArrayReadF90(X,xx,ierr))
    229:   PetscCall(VecRestoreArrayReadF90(Xdot,xxdot,ierr))
    230:   PetscCall(VecRestoreArrayF90(F,ff,ierr))
    231: end subroutine FormIFunction

    233: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
    234:   implicit none
    235:   PetscInt mx,xs,xe,gxs,gxe
    236:   PetscReal t
    237:   PetscScalar x(2,gxs:gxe),f(2,xs:xe)
    238:   PetscReal a(2),k(2),s(2)
    239:   PetscErrorCode ierr
    240:   PetscInt i,j
    241:   PetscReal hx,u0t(2)
    242:   PetscReal one,two,three,four,six,twelve
    243:   PetscReal half,third,twothird,sixth
    244:   PetscReal twelfth

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

    276: #ifdef EXPLICIT_INTEGRATOR22
    277:   do  i = xs,xe
    278:      f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1))
    279:      f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2))
    280:   end do
    281: #endif

    283: end subroutine FormRHSFunctionLocal

    285: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
    286:   use petscts
    287:   use petscdmda
    288:   implicit none

    290:   TS ts
    291:   PetscReal t
    292:   Vec X,F
    293:   PetscReal user(6)
    294:   PetscErrorCode ierr
    295:   integer user_a,user_k,user_s
    296:   parameter (user_a = 1,user_k = 3,user_s = 5)
    297:   DM             da
    298:   Vec            Xloc
    299:   PetscInt       mx,xs,xe,gxs,gxe
    300:   PetscScalar,pointer :: xx(:), ff(:)

    302:   PetscCall(TSGetDM(ts,da,ierr))
    303:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    305:   !     Scatter ghost points to local vector,using the 2-step process
    306:   !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    307:   !     By placing code between these two statements, computations can be
    308:   !     done while messages are in transition.
    309:   PetscCall(DMGetLocalVector(da,Xloc,ierr))
    310:   PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
    311:   PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))

    313:   ! Get access to vector data
    314:   PetscCall(VecGetArrayReadF90(Xloc,xx,ierr))
    315:   PetscCall(VecGetArrayF90(F,ff,ierr))

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

    319:   PetscCall(VecRestoreArrayReadF90(Xloc,xx,ierr))
    320:   PetscCall(VecRestoreArrayF90(F,ff,ierr))
    321:   PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
    322: end subroutine FormRHSFunction

    324: ! ---------------------------------------------------------------------
    325: !
    326: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    327: !
    328: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
    329:   use petscts
    330:   use petscdmda
    331:   implicit none

    333:   TS ts
    334:   PetscReal t,shift
    335:   Vec X,Xdot
    336:   Mat J,Jpre
    337:   PetscReal user(6)
    338:   PetscErrorCode ierr
    339:   integer user_a,user_k,user_s
    340:   parameter (user_a = 0,user_k = 2,user_s = 4)

    342:   DM             da
    343:   PetscInt       mx,xs,xe,gxs,gxe
    344:   PetscInt       i,i1,row,col
    345:   PetscReal      k1,k2;
    346:   PetscScalar    val(4)

    348:   PetscCall(TSGetDM(ts,da,ierr))
    349:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    351:   i1 = 1
    352:   k1 = user(user_k+1)
    353:   k2 = user(user_k+2)
    354:   do i = xs,xe
    355:      row = i-gxs
    356:      col = i-gxs
    357:      val(1) = shift + k1
    358:      val(2) = -k2
    359:      val(3) = -k1
    360:      val(4) = shift + k2
    361:      PetscCall(MatSetValuesBlockedLocal(Jpre,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
    362:   end do
    363:   PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
    364:   PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
    365:   if (J /= Jpre) then
    366:      PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
    367:      PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
    368:   end if
    369: end subroutine FormIJacobian

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

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

    394: subroutine FormInitialSolution(ts,X,user,ierr)
    395:   use petscts
    396:   use petscdmda
    397:   implicit none

    399:   TS ts
    400:   Vec X
    401:   PetscReal user(6)
    402:   PetscErrorCode ierr
    403:   integer user_a,user_k,user_s
    404:   parameter (user_a = 1,user_k = 3,user_s = 5)

    406:   DM             da
    407:   PetscInt       mx,xs,xe,gxs,gxe
    408:   PetscScalar,pointer :: xx(:)

    410:   PetscCall(TSGetDM(ts,da,ierr))
    411:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    413:   ! Get access to vector data
    414:   PetscCall(VecGetArrayF90(X,xx,ierr))

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

    418:   PetscCall(VecRestoreArrayF90(X,xx,ierr))
    419: end subroutine FormInitialSolution

    421: ! ---------------------------------------------------------------------
    422: !
    423: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
    424: !
    425: subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
    426:   use ex22f_mfmodule
    427:   implicit none
    428:   TS ts
    429:   PetscReal t,shift
    430:   Vec X,Xdot
    431:   Mat J,Jpre
    432:   PetscReal user(6)
    433:   PetscErrorCode ierr

    435:   PETSC_SHIFT=shift
    436:   MFuser=user

    438: end subroutine FormIJacobianMF

    440: ! -------------------------------------------------------------------
    441: !
    442: !   MyMult - user provided matrix multiply
    443: !
    444: !   Input Parameters:
    445: !.  X - input vector
    446: !
    447: !   Output Parameter:
    448: !.  F - function vector
    449: !
    450: subroutine  MyMult(A,X,F,ierr)
    451:   use ex22f_mfmodule
    452:   implicit none

    454:   Mat     A
    455:   Vec     X,F

    457:   PetscErrorCode ierr
    458:   PetscScalar shift

    460: !  Mat J,Jpre

    462:   PetscReal user(6)

    464:   integer user_a,user_k,user_s
    465:   parameter (user_a = 0,user_k = 2,user_s = 4)

    467:   DM             da
    468:   PetscInt       mx,xs,xe,gxs,gxe
    469:   PetscInt       i,i1,row,col
    470:   PetscReal      k1,k2;
    471:   PetscScalar    val(4)

    473:   shift=PETSC_SHIFT
    474:   user=MFuser

    476:   PetscCall(TSGetDM(tscontext,da,ierr))
    477:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

    479:   i1 = 1
    480:   k1 = user(user_k+1)
    481:   k2 = user(user_k+2)

    483:   do i = xs,xe
    484:      row = i-gxs
    485:      col = i-gxs
    486:      val(1) = shift + k1
    487:      val(2) = -k2
    488:      val(3) = -k1
    489:      val(4) = shift + k2
    490:      PetscCall(MatSetValuesBlockedLocal(Jmat,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
    491:   end do

    493: !  PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
    494: !  PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
    495: !  if (J /= Jpre) then
    496:      PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr))
    497:      PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr))
    498: !  end if

    500:   PetscCall(MatMult(Jmat,X,F,ierr))

    502: end subroutine MyMult

    504: !
    505: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
    506:   use petscdmda
    507:   implicit none

    509:   Vec X
    510:   DM             da
    511:   PetscInt xs,xe,two
    512:   PetscInt gdof,i
    513:   PetscErrorCode ierr
    514:   PetscScalar data2(2,xs:xe),data(gdof)
    515:   PetscScalar,pointer :: xx(:)

    517:   PetscCall(VecGetArrayReadF90(X,xx,ierr))

    519:   two = 2
    520:   data2=reshape(xx(gdof:gdof),(/two,xe-xs+1/))
    521:   data=reshape(data2,(/gdof/))
    522:   open(1020,file='solution_out_ex22f_mf.txt')
    523:   do i=1,gdof
    524:      write(1020,'(e24.16,1x)') data(i)
    525:   end do
    526:   close(1020)

    528:   PetscCall(VecRestoreArrayReadF90(X,xx,ierr))
    529: end subroutine SaveSolutionToDisk

    531: !/*TEST
    532: !
    533: !    test:
    534: !      args: -da_grid_x 200 -ts_arkimex_type 4
    535: !      requires: !single
    536: !      output_file: output/ex22f_mf_1.out
    537: !
    538: !TEST*/