Actual source code: ex22f_mf.F90
petsc-3.7.5 2017-01-01
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