Actual source code: ex22f_mf.F90

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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 Section 1.5 Writing Application Codes with PETSc context to contain data needed by the
 34:   !     Section 1.5 Writing Application Codes with PETSc-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 Section 1.5 Writing Application Codes with PETSc 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*/