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*/