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 petscdm
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: TS ts
42: Vec X
43: Mat J
44: PetscInt mx
45: PetscBool OptionSaveToDisk
46: PetscErrorCode ierr
47: DM da
48: PetscReal ftime,dt
49: PetscReal one,pone
50: PetscInt im11,i2
51: PetscBool flg
53: PetscInt xs,xe,gxs,gxe,dof,gdof
54: PetscScalar shell_shift
55: Mat A
57: im11 = 11
58: i2 = 2
59: one = 1.0
60: pone = one / 10
62: PetscCallA(PetscInitialize(ierr))
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: ! Create distributed array (DMDA) to manage parallel grid and vectors
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER_ARRAY,da,ierr))
68: PetscCallA(DMSetFromOptions(da,ierr))
69: PetscCallA(DMSetUp(da,ierr))
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: ! Extract global vectors from DMDA
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: PetscCallA(DMCreateGlobalVector(da,X,ierr))
76: ! Initialize user application context
77: ! Use zero-based indexing for command line parameters to match ex22.c
78: user(user_a+1) = 1.0
79: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr))
80: user(user_a+2) = 0.0
81: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr))
82: user(user_k+1) = 1000000.0
83: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0', user(user_k+1),flg,ierr))
84: user(user_k+2) = 2*user(user_k+1)
85: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr))
86: user(user_s+1) = 0.0
87: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr))
88: user(user_s+2) = 1.0
89: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr))
91: OptionSaveToDisk=.FALSE.
92: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sdisk',OptionSaveToDisk,flg,ierr))
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: ! Create timestepping solver context
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
97: tscontext=ts
98: PetscCallA(TSSetDM(ts,da,ierr))
99: PetscCallA(TSSetType(ts,TSARKIMEX,ierr))
100: PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
102: ! - - - - - - - - -- - - - -
103: ! Matrix free setup
104: PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
105: dof=i2*(xe-xs+1)
106: gdof=i2*(gxe-gxs+1)
107: PetscCallA(MatCreateShell(PETSC_COMM_WORLD,dof,dof,gdof,gdof,shell_shift,A,ierr))
109: PetscCallA(MatShellSetOperation(A,MATOP_MULT,MyMult,ierr))
110: ! - - - - - - - - - - - -
112: PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr))
113: PetscCallA(DMSetMatType(da,MATAIJ,ierr))
114: PetscCallA(DMCreateMatrix(da,J,ierr))
116: Jmat=J
118: PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr))
119: PetscCallA(TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr))
121: ftime = 1.0
122: PetscCallA(TSSetMaxTime(ts,ftime,ierr))
123: PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Set initial conditions
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: PetscCallA(FormInitialSolution(ts,X,user,ierr))
129: PetscCallA(TSSetSolution(ts,X,ierr))
130: PetscCallA(VecGetSize(X,mx,ierr))
131: ! Advective CFL, I don't know why it needs so much safety factor.
132: dt = pone * max(user(user_a+1),user(user_a+2)) / mx
133: PetscCallA(TSSetTimeStep(ts,dt,ierr))
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! Set runtime options
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: PetscCallA(TSSetFromOptions(ts,ierr))
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: ! Solve nonlinear system
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: PetscCallA(TSSolve(ts,X,ierr))
145: if (OptionSaveToDisk) then
146: PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
147: dof=i2*(xe-xs+1)
148: gdof=i2*(gxe-gxs+1)
149: call SaveSolutionToDisk(da,X,gdof,xs,xe)
150: end if
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! Free work space.
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: PetscCallA(MatDestroy(A,ierr))
156: PetscCallA(MatDestroy(J,ierr))
157: PetscCallA(VecDestroy(X,ierr))
158: PetscCallA(TSDestroy(ts,ierr))
159: PetscCallA(DMDestroy(da,ierr))
160: PetscCallA(PetscFinalize(ierr))
161: contains
163: ! Small helper to extract the layout, result uses 1-based indexing.
164: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
165: use petscdm
166: implicit none
168: DM da
169: PetscInt mx,xs,xe,gxs,gxe
170: PetscErrorCode ierr
171: PetscInt xm,gxm
172: 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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,ierr))
173: PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
174: PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
175: xs = xs + 1
176: gxs = gxs + 1
177: xe = xs + xm - 1
178: gxe = gxs + gxm - 1
179: end subroutine GetLayout
181: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
182: implicit none
183: PetscInt mx,xs,xe,gxs,gxe
184: PetscScalar x(2,xs:xe)
185: PetscScalar xdot(2,xs:xe)
186: PetscScalar f(2,xs:xe)
187: PetscReal a(2),k(2),s(2)
188: PetscErrorCode ierr
189: PetscInt i
190: do i = xs,xe
191: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
192: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
193: end do
194: end subroutine FormIFunctionLocal
196: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
197: use petscdm
198: use petscts
199: implicit none
201: TS ts
202: PetscReal t
203: Vec X,Xdot,F
204: PetscReal user(6)
205: PetscErrorCode ierr
206: integer user_a,user_k,user_s
207: parameter (user_a = 1,user_k = 3,user_s = 5)
209: DM da
210: PetscInt mx,xs,xe,gxs,gxe
211: PetscScalar,pointer :: xx(:),xxdot(:),ff(:)
213: PetscCall(TSGetDM(ts,da,ierr))
214: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
216: ! Get access to vector data
217: PetscCall(VecGetArrayRead(X,xx,ierr))
218: PetscCall(VecGetArrayRead(Xdot,xxdot,ierr))
219: PetscCall(VecGetArray(F,ff,ierr))
221: PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx,xxdot,ff,user(user_a),user(user_k),user(user_s),ierr))
223: PetscCall(VecRestoreArrayRead(X,xx,ierr))
224: PetscCall(VecRestoreArrayRead(Xdot,xxdot,ierr))
225: PetscCall(VecRestoreArray(F,ff,ierr))
226: end subroutine FormIFunction
228: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
229: implicit none
230: PetscInt mx,xs,xe,gxs,gxe
231: PetscReal t
232: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
233: PetscReal a(2),k(2),s(2)
234: PetscErrorCode ierr
235: PetscInt i,j
236: PetscReal hx,u0t(2)
237: PetscReal one,two,three,four,six,twelve
238: PetscReal half,third,twothird,sixth
239: PetscReal twelfth
241: one = 1.0
242: two = 2.0
243: three = 3.0
244: four = 4.0
245: six = 6.0
246: twelve = 12.0
247: hx = one / mx
248: u0t(1) = one - sin(twelve*t)**four
249: u0t(2) = 0.0
250: half = one/two
251: third = one / three
252: twothird = two / three
253: sixth = one / six
254: twelfth = one / twelve
255: do i = xs,xe
256: do j = 1,2
257: if (i .eq. 1) then
258: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) + sixth*x(j,i+2))
259: else if (i .eq. 2) then
260: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2))
261: else if (i .eq. mx-1) then
262: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) - half*x(j,i) -third*x(j,i+1))
263: else if (i .eq. mx) then
264: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
265: else
266: 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))
267: end if
268: end do
269: end do
271: #ifdef EXPLICIT_INTEGRATOR22
272: do i = xs,xe
273: f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1))
274: f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2))
275: end do
276: #endif
278: end subroutine FormRHSFunctionLocal
280: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
281: use petscts
282: use petscdm
283: implicit none
285: TS ts
286: PetscReal t
287: Vec X,F
288: PetscReal user(6)
289: PetscErrorCode ierr
290: integer user_a,user_k,user_s
291: parameter (user_a = 1,user_k = 3,user_s = 5)
292: DM da
293: Vec Xloc
294: PetscInt mx,xs,xe,gxs,gxe
295: PetscScalar,pointer :: xx(:), ff(:)
297: PetscCall(TSGetDM(ts,da,ierr))
298: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
300: ! Scatter ghost points to local vector,using the 2-step process
301: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
302: ! By placing code between these two statements, computations can be
303: ! done while messages are in transition.
304: PetscCall(DMGetLocalVector(da,Xloc,ierr))
305: PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
306: PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))
308: ! Get access to vector data
309: PetscCall(VecGetArrayRead(Xloc,xx,ierr))
310: PetscCall(VecGetArray(F,ff,ierr))
312: PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx,ff,user(user_a),user(user_k),user(user_s),ierr))
314: PetscCall(VecRestoreArrayRead(Xloc,xx,ierr))
315: PetscCall(VecRestoreArray(F,ff,ierr))
316: PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
317: end subroutine FormRHSFunction
319: ! ---------------------------------------------------------------------
320: !
321: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
322: !
323: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
324: use petscts
325: use petscdm
326: implicit none
328: TS ts
329: PetscReal t,shift
330: Vec X,Xdot
331: Mat J,Jpre
332: PetscReal user(6)
333: PetscErrorCode ierr
334: integer user_a,user_k,user_s
335: parameter (user_a = 0,user_k = 2,user_s = 4)
337: DM da
338: PetscInt mx,xs,xe,gxs,gxe
339: PetscInt i,i1,row,col
340: PetscReal k1,k2
341: PetscScalar val(4)
343: PetscCall(TSGetDM(ts,da,ierr))
344: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
346: i1 = 1
347: k1 = user(user_k+1)
348: k2 = user(user_k+2)
349: do i = xs,xe
350: row = i-gxs
351: col = i-gxs
352: val(1) = shift + k1
353: val(2) = -k2
354: val(3) = -k1
355: val(4) = shift + k2
356: PetscCall(MatSetValuesBlockedLocal(Jpre,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
357: end do
358: PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
359: PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
360: if (J /= Jpre) then
361: PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
362: PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
363: end if
364: end subroutine FormIJacobian
366: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
367: implicit none
368: PetscInt mx,xs,xe,gxs,gxe
369: PetscScalar x(2,xs:xe)
370: PetscReal a(2),k(2),s(2)
371: PetscErrorCode ierr
373: PetscInt i
374: PetscReal one,hx,r,ik
375: one = 1.0
376: hx = one / mx
377: do i=xs,xe
378: r = i*hx
379: if (k(2) .ne. 0.0) then
380: ik = one/k(2)
381: else
382: ik = one
383: end if
384: x(1,i) = one + s(2)*r
385: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
386: end do
387: end subroutine FormInitialSolutionLocal
389: subroutine FormInitialSolution(ts,X,user,ierr)
390: use petscts
391: use petscdm
392: implicit none
394: TS ts
395: Vec X
396: PetscReal user(6)
397: PetscErrorCode ierr
398: integer user_a,user_k,user_s
399: parameter (user_a = 1,user_k = 3,user_s = 5)
401: DM da
402: PetscInt mx,xs,xe,gxs,gxe
403: PetscScalar,pointer :: xx(:)
405: PetscCall(TSGetDM(ts,da,ierr))
406: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
408: ! Get access to vector data
409: PetscCall(VecGetArray(X,xx,ierr))
411: PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx,user(user_a),user(user_k),user(user_s),ierr))
413: PetscCall(VecRestoreArray(X,xx,ierr))
414: end subroutine FormInitialSolution
416: ! ---------------------------------------------------------------------
417: !
418: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
419: !
420: subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
421: use ex22f_mfmodule
422: implicit none
423: TS ts
424: PetscReal t,shift
425: Vec X,Xdot
426: Mat J,Jpre
427: PetscReal user(6)
428: PetscErrorCode ierr
430: PETSC_SHIFT=shift
431: MFuser=user
433: end subroutine FormIJacobianMF
435: ! -------------------------------------------------------------------
436: !
437: ! MyMult - user provided matrix multiply
438: !
439: ! Input Parameters:
440: !. X - input vector
441: !
442: ! Output Parameter:
443: !. F - function vector
444: !
445: subroutine MyMult(A,X,F,ierr)
446: use ex22f_mfmodule
447: implicit none
449: Mat A
450: Vec X,F
452: PetscErrorCode ierr
453: PetscScalar shift
455: ! Mat J,Jpre
457: PetscReal user(6)
459: integer user_a,user_k,user_s
460: parameter (user_a = 0,user_k = 2,user_s = 4)
462: DM da
463: PetscInt mx,xs,xe,gxs,gxe
464: PetscInt i,i1,row,col
465: PetscReal k1,k2
466: PetscScalar val(4)
468: shift=PETSC_SHIFT
469: user=MFuser
471: PetscCall(TSGetDM(tscontext,da,ierr))
472: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
474: i1 = 1
475: k1 = user(user_k+1)
476: k2 = user(user_k+2)
478: do i = xs,xe
479: row = i-gxs
480: col = i-gxs
481: val(1) = shift + k1
482: val(2) = -k2
483: val(3) = -k1
484: val(4) = shift + k2
485: PetscCall(MatSetValuesBlockedLocal(Jmat,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
486: end do
488: ! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
489: ! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
490: ! if (J /= Jpre) then
491: PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr))
492: PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr))
493: ! end if
495: PetscCall(MatMult(Jmat,X,F,ierr))
497: end subroutine MyMult
499: !
500: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
501: use petscdm
502: implicit none
504: Vec X
505: DM da
506: PetscInt xs,xe,two
507: PetscInt gdof,i
508: PetscErrorCode ierr
509: PetscScalar data2(2,xs:xe),data(gdof)
510: PetscScalar,pointer :: xx(:)
512: PetscCall(VecGetArrayRead(X,xx,ierr))
514: two = 2
515: data2=reshape(xx(gdof:gdof),(/two,xe-xs+1/))
516: data=reshape(data2,(/gdof/))
517: open(1020,file='solution_out_ex22f_mf.txt')
518: do i=1,gdof
519: write(1020,'(e24.16,1x)') data(i)
520: end do
521: close(1020)
523: PetscCall(VecRestoreArrayRead(X,xx,ierr))
524: end subroutine SaveSolutionToDisk
525: end program main
527: !/*TEST
528: !
529: ! test:
530: ! args: -da_grid_x 200 -ts_arkimex_type 4
531: ! requires: !single
532: ! output_file: output/empty.out
533: !
534: !TEST*/