Actual source code: ex22f_mf.F90
petsc-3.14.6 2021-03-30
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*/