Actual source code: ex22f_mf.F90

petsc-3.7.3 2016-08-01
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: #define MF_EX22F_MF
 18: !#undef  MF_EX22F_MF


 21: #ifdef MF_EX22F_MF
 22:   module PETScShiftMod
 23: #include <petsc/finclude/petscsys.h>
 24: #include <petsc/finclude/petscts.h>
 25: #include <petsc/finclude/petscmat.h>
 26:     PetscScalar::PETSC_SHIFT
 27:     TS::tscontext
 28:     Mat::Jmat
 29:     PetscReal::MFuser(6)
 30:   end module PETScShiftMod
 31: #endif

 33: program main
 34: #ifdef MF_EX22F_MF
 35:   use PETScShiftMod, only :  tscontext,Jmat
 36: #endif
 37:   implicit none
 38: #include <petsc/finclude/petscsys.h>
 39: #include <petsc/finclude/petscvec.h>
 40: #include <petsc/finclude/petscmat.h>
 41: #include <petsc/finclude/petscsnes.h>
 42: #include <petsc/finclude/petscts.h>
 43: #include <petsc/finclude/petscdm.h>
 44: #include <petsc/finclude/petscdmda.h>
 45:   !
 46:   !     Create an application context to contain data needed by the
 47:   !     application-provided call-back routines, FormJacobian() and
 48:   !     FormFunction(). We use a double precision array with six
 49:   !     entries, two for each problem parameter a, k, s.
 50:   !
 51:   PetscReal user(6)
 52:   integer user_a,user_k,user_s
 53:   parameter (user_a = 0,user_k = 2,user_s = 4)

 55:   external FormRHSFunction,FormIFunction
 56:   external FormInitialSolution

 58: #ifndef MF_EX22F_MF
 59:   external FormIJacobian
 60: #endif

 62: #ifdef MF_EX22F_MF
 63:   external MyMult,FormIJacobianMF
 64: #endif

 66:   TS             ts
 67:   Vec            X
 68:   Mat            J
 69:   PetscInt       maxsteps,mx
 70:   PetscBool      OptionSaveToDisk
 71:   PetscErrorCode ierr
 72:   DM             da
 73:   PetscReal      ftime,dt
 74:   PetscReal      zero,one,pone
 75:   PetscInt       im11,i2
 76:   PetscBool      flg


 79:   PetscInt       xs,xe,gxs,gxe,dof,gdof
 80:   PetscScalar    shell_shift
 81:   Mat            A

 83:   im11 = -11
 84:   i2   = 2
 85:   zero = 0.0
 86:   one = 1.0
 87:   pone = one / 10

 89:   call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 91:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:   !  Create distributed array (DMDA) to manage parallel grid and vectors
 93:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:   call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2, &
 95:        PETSC_NULL_INTEGER,da,ierr)

 97:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:   !    Extract global vectors from DMDA;
 99:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:   call DMCreateGlobalVector(da,X,ierr)

102:   ! Initialize user application context
103:   ! Use zero-based indexing for command line parameters to match ex22.c
104:   user(user_a+1) = 1.0
105:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,           &
106:                                '-a0',user(user_a+1),flg,ierr)
107:   user(user_a+2) = 0.0
108:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,               &
109:                            '-a1',user(user_a+2),flg,ierr)
110:   user(user_k+1) = 1000000.0
111:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,       &
112:                            '-k0', user(user_k+1),flg,ierr)
113:   user(user_k+2) = 2*user(user_k+1)
114:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,          &
115:                            '-k1', user(user_k+2),flg,ierr)
116:   user(user_s+1) = 0.0
117:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,       &
118:                           '-s0',user(user_s+1),flg,ierr)
119:   user(user_s+2) = 1.0
120:   call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,      &
121:                            '-s1',user(user_s+2),flg,ierr)

123:   OptionSaveToDisk=.FALSE.
124:       call PetscOptionsGetBool(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,      &
125:      &                         '-sdisk',OptionSaveToDisk,flg,ierr)
126:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:   !    Create timestepping solver context
128:   !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:   call TSCreate(PETSC_COMM_WORLD,ts,ierr)
130: #ifdef MF_EX22F_MF
131:   tscontext=ts
132: #endif
133:   call TSSetDM(ts,da,ierr)
134:   call TSSetType(ts,TSARKIMEX,ierr)
135:   call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction,       &
136:        user,ierr)

138: #ifdef MF_EX22F_MF
139:   ! - - - - - - - - -- - - - -
140:   !   Matrix free setup
141:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
142:   dof=i2*(xe-xs+1)
143:   gdof=i2*(gxe-gxs+1)
144:   call MatCreateShell(PETSC_COMM_WORLD,                             &
145:        dof,dof,gdof,gdof,                                           &
146:        shell_shift,A,ierr);

148:   call MatShellSetOperation(A,MATOP_MULT,MyMult,ierr)
149:   ! - - - - - - - - - - - -
150: #endif

152:   call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
153:   call DMSetMatType(da,MATAIJ,ierr)
154:   call DMCreateMatrix(da,J,ierr)

156: #ifdef MF_EX22F_MF
157:   Jmat=J
158: #endif

160: #ifndef MF_EX22F_MF
161:   call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)
162: #endif
163: #ifdef MF_EX22F_MF
164:   call TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr)
165: #endif

167:   ftime = 1.0
168:   maxsteps = 10000
169:   call TSSetDuration(ts,maxsteps,ftime,ierr)
170:   call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)

172:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173:   !  Set initial conditions
174:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175:   call FormInitialSolution(ts,X,user,ierr)
176:   call TSSetSolution(ts,X,ierr)
177:   call VecGetSize(X,mx,ierr)
178:   !  Advective CFL, I don't know why it needs so much safety factor.
179:   dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
180:   call TSSetInitialTimeStep(ts,zero,dt,ierr)

182:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:   !   Set runtime options
184:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185:   call TSSetFromOptions(ts,ierr)

187:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:   !  Solve nonlinear system
189:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:   call TSSolve(ts,X,ierr)

192:   if (OptionSaveToDisk) then
193:      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
194:      dof=i2*(xe-xs+1)
195:      gdof=i2*(gxe-gxs+1)
196:      call SaveSolutionToDisk(da,X,gdof,xs,xe)
197:   end if

199:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200:   !  Free work space.
201: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: #ifdef MF_EX22F_MF
203:   call MatDestroy(A,ierr)
204: #endif
205:   call MatDestroy(J,ierr)
206:   call VecDestroy(X,ierr)
207:   call TSDestroy(ts,ierr)
208:   call DMDestroy(da,ierr)
209:   call PetscFinalize(ierr)
210: end program main

212: ! Small helper to extract the layout, result uses 1-based indexing.
213: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
214:   implicit none
215: #include <petsc/finclude/petscsys.h>
216: #include <petsc/finclude/petscdm.h>
217: #include <petsc/finclude/petscdmda.h>
218:   DM da
219:   PetscInt mx,xs,xe,gxs,gxe
220:   PetscErrorCode ierr
221:   PetscInt xm,gxm
222:   call DMDAGetInfo(da,PETSC_NULL_INTEGER,                     &
223:        mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                    &
224:        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
225:        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                       &
226:        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
227:        PETSC_NULL_INTEGER,ierr)
228:   call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,  &
229:        xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
230:   call DMDAGetGhostCorners(da,                                      &
231:        gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                   &
232:        gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
233:   xs = xs + 1
234:   gxs = gxs + 1
235:   xe = xs + xm - 1
236:   gxe = gxs + gxm - 1
237: end subroutine GetLayout

239: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,          &
240:      a,k,s,ierr)
241:   implicit none
242:   PetscInt mx,xs,xe,gxs,gxe
243:   PetscScalar x(2,xs:xe)
244:   PetscScalar xdot(2,xs:xe)
245:   PetscScalar f(2,xs:xe)
246:   PetscReal a(2),k(2),s(2)
247:   PetscErrorCode ierr
248:   PetscInt i
249:   do  i = xs,xe
250:      f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
251:      f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
252:   end do
253: end subroutine FormIFunctionLocal

255: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
256:   implicit none
257: #include <petsc/finclude/petscsys.h>
258: #include <petsc/finclude/petscvec.h>
259: #include <petsc/finclude/petscmat.h>
260: #include <petsc/finclude/petscsnes.h>
261: #include <petsc/finclude/petscts.h>
262: #include <petsc/finclude/petscdm.h>
263: #include <petsc/finclude/petscdmda.h>
264:   TS ts
265:   PetscReal t
266:   Vec X,Xdot,F
267:   PetscReal user(6)
268:   PetscErrorCode ierr
269:   integer user_a,user_k,user_s
270:   parameter (user_a = 1,user_k = 3,user_s = 5)

272:   DM             da
273:   PetscInt       mx,xs,xe,gxs,gxe
274:   PetscOffset    ixx,ixxdot,iff
275:   PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

277:   call TSGetDM(ts,da,ierr)
278:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

280:   ! Get access to vector data
281:   call VecGetArrayRead(X,xx,ixx,ierr)
282:   call VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)
283:   call VecGetArray(F,ff,iff,ierr)

285:   call FormIFunctionLocal(mx,xs,xe,gxs,gxe,                         &
286:        xx(ixx),xxdot(ixxdot),ff(iff),                               &
287:        user(user_a),user(user_k),user(user_s),ierr)

289:   call VecRestoreArrayRead(X,xx,ixx,ierr)
290:   call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)
291:   call VecRestoreArray(F,ff,iff,ierr)
292: end subroutine FormIFunction

294: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,           &
295:      a,k,s,ierr)
296:   implicit none
297:   PetscInt mx,xs,xe,gxs,gxe
298:   PetscReal t
299:   PetscScalar x(2,gxs:gxe),f(2,xs:xe)
300:   PetscReal a(2),k(2),s(2)
301:   PetscErrorCode ierr
302:   PetscInt i,j
303:   PetscReal hx,u0t(2)
304:   PetscReal one,two,three,four,six,twelve
305:   PetscReal half,third,twothird,sixth
306:   PetscReal twelfth

308:   one = 1.0
309:   two = 2.0
310:   three = 3.0
311:   four = 4.0
312:   six = 6.0
313:   twelve = 12.0
314:   hx = one / mx
315:   u0t(1) = one - sin(twelve*t)**four
316:   u0t(2) = 0.0
317:   half = one/two
318:   third = one / three
319:   twothird = two / three
320:   sixth = one / six
321:   twelfth = one / twelve
322:   do  i = xs,xe
323:      do  j = 1,2
324:         if (i .eq. 1) then
325:            f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
326:                 + sixth*x(j,i+2))
327:         else if (i .eq. 2) then
328:            f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
329:                 - twothird*x(j,i+1) + twelfth*x(j,i+2))
330:         else if (i .eq. mx-1) then
331:            f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1)             &
332:                 - half*x(j,i) -third*x(j,i+1))
333:         else if (i .eq. mx) then
334:            f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
335:         else
336:            f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
337:                 + twothird*x(j,i-1)                                 &
338:                 - twothird*x(j,i+1) + twelfth*x(j,i+2))
339:         end if
340:      end do
341:   end do

343: #ifdef EXPLICIT_INTEGRATOR22
344:   do  i = xs,xe
345:      f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1))
346:      f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2))
347:   end do
348: #endif

350: end subroutine FormRHSFunctionLocal

352: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
353:   implicit none
354: #include <petsc/finclude/petscsys.h>
355: #include <petsc/finclude/petscvec.h>
356: #include <petsc/finclude/petscmat.h>
357: #include <petsc/finclude/petscsnes.h>
358: #include <petsc/finclude/petscts.h>
359: #include <petsc/finclude/petscdm.h>
360: #include <petsc/finclude/petscdmda.h>
361:   TS ts
362:   PetscReal t
363:   Vec X,F
364:   PetscReal user(6)
365:   PetscErrorCode ierr
366:   integer user_a,user_k,user_s
367:   parameter (user_a = 1,user_k = 3,user_s = 5)
368:   DM             da
369:   Vec            Xloc
370:   PetscInt       mx,xs,xe,gxs,gxe
371:   PetscOffset    ixx,iff
372:   PetscScalar    xx(0:1),ff(0:1)

374:   call TSGetDM(ts,da,ierr)
375:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

377:   !     Scatter ghost points to local vector,using the 2-step process
378:   !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
379:   !     By placing code between these two statements, computations can be
380:   !     done while messages are in transition.
381:   call DMGetLocalVector(da,Xloc,ierr)
382:   call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
383:   call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)

385:   ! Get access to vector data
386:   call VecGetArrayRead(Xloc,xx,ixx,ierr)
387:   call VecGetArray(F,ff,iff,ierr)

389:   call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,                       &
390:        t,xx(ixx),ff(iff),                                           &
391:        user(user_a),user(user_k),user(user_s),ierr)

393:   call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
394:   call VecRestoreArray(F,ff,iff,ierr)
395:   call DMRestoreLocalVector(da,Xloc,ierr)
396: end subroutine FormRHSFunction

398: #ifndef MF_EX22F_MF
399: ! ---------------------------------------------------------------------
400: !
401: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
402: !
403: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
404:   implicit none
405: #include <petsc/finclude/petscsys.h>
406: #include <petsc/finclude/petscvec.h>
407: #include <petsc/finclude/petscmat.h>
408: #include <petsc/finclude/petscsnes.h>
409: #include <petsc/finclude/petscts.h>
410: #include <petsc/finclude/petscdm.h>
411: #include <petsc/finclude/petscdmda.h>
412:   TS ts
413:   PetscReal t,shift
414:   Vec X,Xdot
415:   Mat J,Jpre
416:   PetscReal user(6)
417:   PetscErrorCode ierr
418:   integer user_a,user_k,user_s
419:   parameter (user_a = 0,user_k = 2,user_s = 4)

421:   DM             da
422:   PetscInt       mx,xs,xe,gxs,gxe
423:   PetscInt       i,i1,row,col
424:   PetscReal      k1,k2;
425:   PetscScalar    val(4)

427:   call TSGetDM(ts,da,ierr)
428:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

430:   i1 = 1
431:   k1 = user(user_k+1)
432:   k2 = user(user_k+2)
433:   do i = xs,xe
434:      row = i-gxs
435:      col = i-gxs
436:      val(1) = shift + k1
437:      val(2) = -k2
438:      val(3) = -k1
439:      val(4) = shift + k2
440:      call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,          &
441:           INSERT_VALUES,ierr)
442:   end do
443:   call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
444:   call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
445:   if (J /= Jpre) then
446:      call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
447:      call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
448:   end if
449: end subroutine FormIJacobian
450: #endif

452: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
453:   implicit none
454:   PetscInt mx,xs,xe,gxs,gxe
455:   PetscScalar x(2,xs:xe)
456:   PetscReal a(2),k(2),s(2)
457:   PetscErrorCode ierr

459:   PetscInt i
460:   PetscReal one,hx,r,ik
461:   one = 1.0
462:   hx = one / mx
463:   do i=xs,xe
464:      r = i*hx
465:      if (k(2) .ne. 0.0) then
466:         ik = one/k(2)
467:      else
468:         ik = one
469:      end if
470:      x(1,i) = one + s(2)*r
471:      x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
472:   end do
473: end subroutine FormInitialSolutionLocal

475: subroutine FormInitialSolution(ts,X,user,ierr)
476:   implicit none
477: #include <petsc/finclude/petscsys.h>
478: #include <petsc/finclude/petscvec.h>
479: #include <petsc/finclude/petscmat.h>
480: #include <petsc/finclude/petscsnes.h>
481: #include <petsc/finclude/petscts.h>
482: #include <petsc/finclude/petscdm.h>
483: #include <petsc/finclude/petscdmda.h>
484:   TS ts
485:   Vec X
486:   PetscReal user(6)
487:   PetscErrorCode ierr
488:   integer user_a,user_k,user_s
489:   parameter (user_a = 1,user_k = 3,user_s = 5)

491:   DM             da
492:   PetscInt       mx,xs,xe,gxs,gxe
493:   PetscOffset    ixx
494:   PetscScalar    xx(0:1)

496:   call TSGetDM(ts,da,ierr)
497:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

499:   ! Get access to vector data
500:   call VecGetArray(X,xx,ixx,ierr)

502:   call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,                       &
503:        xx(ixx),user(user_a),user(user_k),user(user_s),ierr)

505:   call VecRestoreArray(X,xx,ixx,ierr)
506: end subroutine FormInitialSolution


509: #ifdef MF_EX22F_MF
510: ! ---------------------------------------------------------------------
511: !
512: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
513: !
514: subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
515:   use PETScShiftMod, only :  PETSC_SHIFT,MFuser
516:   implicit none
517: #include <petsc/finclude/petscsys.h>
518: #include <petsc/finclude/petscvec.h>
519: #include <petsc/finclude/petscmat.h>
520: #include <petsc/finclude/petscsnes.h>
521: #include <petsc/finclude/petscts.h>
522: #include <petsc/finclude/petscdm.h>
523: #include <petsc/finclude/petscdmda.h>
524:   TS ts
525:   PetscReal t,shift
526:   Vec X,Xdot
527:   Mat J,Jpre
528:   PetscReal user(6)
529:   PetscErrorCode ierr

531:   !  call MatShellSetContext(J,shift,ierr)
532:   PETSC_SHIFT=shift
533:   MFuser=user

535: end subroutine FormIJacobianMF

537: ! -------------------------------------------------------------------
538: !
539: !   MyMult - user provided matrix multiply
540: !
541: !   Input Parameters:
542: !.  X - input vector
543: !
544: !   Output Parameter:
545: !.  F - function vector
546: !
547: subroutine  MyMult(A,X,F,ierr)
548:   use PETScShiftMod, only :  PETSC_SHIFT,tscontext,Jmat,MFuser
549:   implicit none

551: #include <petsc/finclude/petscsys.h>
552: #include <petsc/finclude/petscvec.h>
553: #include <petsc/finclude/petscmat.h>
554: #include <petsc/finclude/petscpc.h>
555: #include <petsc/finclude/petscts.h>
556: #include <petsc/finclude/petscvec.h90>

558:   Mat     A
559:   Vec     X,F

561:   PetscErrorCode ierr
562:   PetscScalar shift

564: !  Mat J,Jpre

566:   PetscReal user(6)

568:   integer user_a,user_k,user_s
569:   parameter (user_a = 0,user_k = 2,user_s = 4)

571:   DM             da
572:   PetscInt       mx,xs,xe,gxs,gxe
573:   PetscInt       i,i1,row,col
574:   PetscReal      k1,k2;
575:   PetscScalar    val(4)

577:   !call MatShellGetContext(A,shift,ierr)
578:   shift=PETSC_SHIFT
579:   user=MFuser

581:   call TSGetDM(tscontext,da,ierr)
582:   call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

584:   i1 = 1
585:   k1 = user(user_k+1)
586:   k2 = user(user_k+2)

588:   do i = xs,xe
589:      row = i-gxs
590:      col = i-gxs
591:      val(1) = shift + k1
592:      val(2) = -k2
593:      val(3) = -k1
594:      val(4) = shift + k2
595:      call MatSetValuesBlockedLocal(Jmat,i1,row,i1,col,val,          &
596:           INSERT_VALUES,ierr)
597:   end do

599: !  call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
600: !  call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
601: !  if (J /= Jpre) then
602:      call MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr)
603:      call MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr)
604: !  end if

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

608:   return
609: end subroutine MyMult
610: #endif


613: !
614: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
615:   implicit none
616: #include <petsc/finclude/petscsys.h>
617: #include <petsc/finclude/petscvec.h>
618: #include <petsc/finclude/petscmat.h>
619: #include <petsc/finclude/petscsnes.h>
620: #include <petsc/finclude/petscts.h>
621: #include <petsc/finclude/petscdm.h>
622: #include <petsc/finclude/petscdmda.h>

624:   Vec X
625:   DM             da
626:   PetscInt xs,xe,two
627:   PetscInt gdof,i
628:   PetscErrorCode ierr
629:   PetscOffset    ixx
630:   PetscScalar data2(2,xs:xe),data(gdof)
631:   PetscScalar    xx(0:1)


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

636:   two = 2
637:   data2=reshape(xx(ixx:ixx+gdof),(/two,xe-xs+1/))
638:   data=reshape(data2,(/gdof/))
639:   open(1020,file='solution_out_ex22f_mf.txt')
640:   do i=1,gdof
641:      write(1020,'(e24.16,1x)') data(i)
642:   end do
643:   close(1020)

645:   call VecRestoreArrayRead(X,xx,ixx,ierr)
646: end subroutine SaveSolutionToDisk