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,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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,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:   PetscOffset    ixx,ixxdot,iff
217:   PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

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

222:   ! Get access to vector data
223:   PetscCall(VecGetArrayRead(X,xx,ixx,ierr))
224:   PetscCall(VecGetArrayRead(Xdot,xxdot,ixxdot,ierr))
225:   PetscCall(VecGetArray(F,ff,iff,ierr))

227:   PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr))

229:   PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr))
230:   PetscCall(VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr))
231:   PetscCall(VecRestoreArray(F,ff,iff,ierr))
232: end subroutine FormIFunction

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

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

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

284: end subroutine FormRHSFunctionLocal

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

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

304:   PetscCall(TSGetDM(ts,da,ierr))
305:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

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

315:   ! Get access to vector data
316:   PetscCall(VecGetArrayRead(Xloc,xx,ixx,ierr))
317:   PetscCall(VecGetArray(F,ff,iff,ierr))

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

321:   PetscCall(VecRestoreArrayRead(Xloc,xx,ixx,ierr))
322:   PetscCall(VecRestoreArray(F,ff,iff,ierr))
323:   PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
324: end subroutine FormRHSFunction

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

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

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

350:   PetscCall(TSGetDM(ts,da,ierr))
351:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

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

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

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

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

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

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

413:   PetscCall(TSGetDM(ts,da,ierr))
414:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

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

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

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

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

438:   PETSC_SHIFT=shift
439:   MFuser=user

441: end subroutine FormIJacobianMF

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

457:   Mat     A
458:   Vec     X,F

460:   PetscErrorCode ierr
461:   PetscScalar shift

463: !  Mat J,Jpre

465:   PetscReal user(6)

467:   integer user_a,user_k,user_s
468:   parameter (user_a = 0,user_k = 2,user_s = 4)

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

476:   shift=PETSC_SHIFT
477:   user=MFuser

479:   PetscCall(TSGetDM(tscontext,da,ierr))
480:   PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

482:   i1 = 1
483:   k1 = user(user_k+1)
484:   k2 = user(user_k+2)

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

496: !  PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
497: !  PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
498: !  if (J /= Jpre) then
499:      PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr))
500:      PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr))
501: !  end if

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

505:   return
506: end subroutine MyMult

508: !
509: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
510:   use petscdmda
511:   implicit none

513:   Vec X
514:   DM             da
515:   PetscInt xs,xe,two
516:   PetscInt gdof,i
517:   PetscErrorCode ierr
518:   PetscOffset    ixx
519:   PetscScalar data2(2,xs:xe),data(gdof)
520:   PetscScalar    xx(0:1)

522:   PetscCall(VecGetArrayRead(X,xx,ixx,ierr))

524:   two = 2
525:   data2=reshape(xx(ixx:ixx+gdof),(/two,xe-xs+1/))
526:   data=reshape(data2,(/gdof/))
527:   open(1020,file='solution_out_ex22f_mf.txt')
528:   do i=1,gdof
529:      write(1020,'(e24.16,1x)') data(i)
530:   end do
531:   close(1020)

533:   PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr))
534: end subroutine SaveSolutionToDisk

536: !/*TEST
537: !
538: !    test:
539: !      args: -da_grid_x 200 -ts_arkimex_type 4
540: !      requires: !single
541: !      output_file: output/ex22f_mf_1.out
542: !
543: !TEST*/