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_ARRAY,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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,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: PetscScalar,pointer :: xx(:),xxdot(:),ff(:)
218: PetscCall(TSGetDM(ts,da,ierr))
219: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
221: ! Get access to vector data
222: PetscCall(VecGetArrayReadF90(X,xx,ierr))
223: PetscCall(VecGetArrayReadF90(Xdot,xxdot,ierr))
224: PetscCall(VecGetArrayF90(F,ff,ierr))
226: PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx,xxdot,ff,user(user_a),user(user_k),user(user_s),ierr))
228: PetscCall(VecRestoreArrayReadF90(X,xx,ierr))
229: PetscCall(VecRestoreArrayReadF90(Xdot,xxdot,ierr))
230: PetscCall(VecRestoreArrayF90(F,ff,ierr))
231: end subroutine FormIFunction
233: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
234: implicit none
235: PetscInt mx,xs,xe,gxs,gxe
236: PetscReal t
237: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
238: PetscReal a(2),k(2),s(2)
239: PetscErrorCode ierr
240: PetscInt i,j
241: PetscReal hx,u0t(2)
242: PetscReal one,two,three,four,six,twelve
243: PetscReal half,third,twothird,sixth
244: PetscReal twelfth
246: one = 1.0
247: two = 2.0
248: three = 3.0
249: four = 4.0
250: six = 6.0
251: twelve = 12.0
252: hx = one / mx
253: u0t(1) = one - sin(twelve*t)**four
254: u0t(2) = 0.0
255: half = one/two
256: third = one / three
257: twothird = two / three
258: sixth = one / six
259: twelfth = one / twelve
260: do i = xs,xe
261: do j = 1,2
262: if (i .eq. 1) then
263: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) + sixth*x(j,i+2))
264: else if (i .eq. 2) then
265: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2))
266: else if (i .eq. mx-1) then
267: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) - half*x(j,i) -third*x(j,i+1))
268: else if (i .eq. mx) then
269: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
270: else
271: 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))
272: end if
273: end do
274: end do
276: #ifdef EXPLICIT_INTEGRATOR22
277: do i = xs,xe
278: f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1))
279: f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2))
280: end do
281: #endif
283: end subroutine FormRHSFunctionLocal
285: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
286: use petscts
287: use petscdmda
288: implicit none
290: TS ts
291: PetscReal t
292: Vec X,F
293: PetscReal user(6)
294: PetscErrorCode ierr
295: integer user_a,user_k,user_s
296: parameter (user_a = 1,user_k = 3,user_s = 5)
297: DM da
298: Vec Xloc
299: PetscInt mx,xs,xe,gxs,gxe
300: PetscScalar,pointer :: xx(:), ff(:)
302: PetscCall(TSGetDM(ts,da,ierr))
303: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
305: ! Scatter ghost points to local vector,using the 2-step process
306: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
307: ! By placing code between these two statements, computations can be
308: ! done while messages are in transition.
309: PetscCall(DMGetLocalVector(da,Xloc,ierr))
310: PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
311: PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))
313: ! Get access to vector data
314: PetscCall(VecGetArrayReadF90(Xloc,xx,ierr))
315: PetscCall(VecGetArrayF90(F,ff,ierr))
317: PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx,ff,user(user_a),user(user_k),user(user_s),ierr))
319: PetscCall(VecRestoreArrayReadF90(Xloc,xx,ierr))
320: PetscCall(VecRestoreArrayF90(F,ff,ierr))
321: PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
322: end subroutine FormRHSFunction
324: ! ---------------------------------------------------------------------
325: !
326: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
327: !
328: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
329: use petscts
330: use petscdmda
331: implicit none
333: TS ts
334: PetscReal t,shift
335: Vec X,Xdot
336: Mat J,Jpre
337: PetscReal user(6)
338: PetscErrorCode ierr
339: integer user_a,user_k,user_s
340: parameter (user_a = 0,user_k = 2,user_s = 4)
342: DM da
343: PetscInt mx,xs,xe,gxs,gxe
344: PetscInt i,i1,row,col
345: PetscReal k1,k2;
346: PetscScalar val(4)
348: PetscCall(TSGetDM(ts,da,ierr))
349: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
351: i1 = 1
352: k1 = user(user_k+1)
353: k2 = user(user_k+2)
354: do i = xs,xe
355: row = i-gxs
356: col = i-gxs
357: val(1) = shift + k1
358: val(2) = -k2
359: val(3) = -k1
360: val(4) = shift + k2
361: PetscCall(MatSetValuesBlockedLocal(Jpre,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
362: end do
363: PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
364: PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
365: if (J /= Jpre) then
366: PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
367: PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
368: end if
369: end subroutine FormIJacobian
371: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
372: implicit none
373: PetscInt mx,xs,xe,gxs,gxe
374: PetscScalar x(2,xs:xe)
375: PetscReal a(2),k(2),s(2)
376: PetscErrorCode ierr
378: PetscInt i
379: PetscReal one,hx,r,ik
380: one = 1.0
381: hx = one / mx
382: do i=xs,xe
383: r = i*hx
384: if (k(2) .ne. 0.0) then
385: ik = one/k(2)
386: else
387: ik = one
388: end if
389: x(1,i) = one + s(2)*r
390: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
391: end do
392: end subroutine FormInitialSolutionLocal
394: subroutine FormInitialSolution(ts,X,user,ierr)
395: use petscts
396: use petscdmda
397: implicit none
399: TS ts
400: Vec X
401: PetscReal user(6)
402: PetscErrorCode ierr
403: integer user_a,user_k,user_s
404: parameter (user_a = 1,user_k = 3,user_s = 5)
406: DM da
407: PetscInt mx,xs,xe,gxs,gxe
408: PetscScalar,pointer :: xx(:)
410: PetscCall(TSGetDM(ts,da,ierr))
411: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
413: ! Get access to vector data
414: PetscCall(VecGetArrayF90(X,xx,ierr))
416: PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx,user(user_a),user(user_k),user(user_s),ierr))
418: PetscCall(VecRestoreArrayF90(X,xx,ierr))
419: end subroutine FormInitialSolution
421: ! ---------------------------------------------------------------------
422: !
423: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
424: !
425: subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
426: use ex22f_mfmodule
427: implicit none
428: TS ts
429: PetscReal t,shift
430: Vec X,Xdot
431: Mat J,Jpre
432: PetscReal user(6)
433: PetscErrorCode ierr
435: PETSC_SHIFT=shift
436: MFuser=user
438: end subroutine FormIJacobianMF
440: ! -------------------------------------------------------------------
441: !
442: ! MyMult - user provided matrix multiply
443: !
444: ! Input Parameters:
445: !. X - input vector
446: !
447: ! Output Parameter:
448: !. F - function vector
449: !
450: subroutine MyMult(A,X,F,ierr)
451: use ex22f_mfmodule
452: implicit none
454: Mat A
455: Vec X,F
457: PetscErrorCode ierr
458: PetscScalar shift
460: ! Mat J,Jpre
462: PetscReal user(6)
464: integer user_a,user_k,user_s
465: parameter (user_a = 0,user_k = 2,user_s = 4)
467: DM da
468: PetscInt mx,xs,xe,gxs,gxe
469: PetscInt i,i1,row,col
470: PetscReal k1,k2;
471: PetscScalar val(4)
473: shift=PETSC_SHIFT
474: user=MFuser
476: PetscCall(TSGetDM(tscontext,da,ierr))
477: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
479: i1 = 1
480: k1 = user(user_k+1)
481: k2 = user(user_k+2)
483: do i = xs,xe
484: row = i-gxs
485: col = i-gxs
486: val(1) = shift + k1
487: val(2) = -k2
488: val(3) = -k1
489: val(4) = shift + k2
490: PetscCall(MatSetValuesBlockedLocal(Jmat,i1,[row],i1,[col],val,INSERT_VALUES,ierr))
491: end do
493: ! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
494: ! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
495: ! if (J /= Jpre) then
496: PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr))
497: PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr))
498: ! end if
500: PetscCall(MatMult(Jmat,X,F,ierr))
502: end subroutine MyMult
504: !
505: subroutine SaveSolutionToDisk(da,X,gdof,xs,xe)
506: use petscdmda
507: implicit none
509: Vec X
510: DM da
511: PetscInt xs,xe,two
512: PetscInt gdof,i
513: PetscErrorCode ierr
514: PetscScalar data2(2,xs:xe),data(gdof)
515: PetscScalar,pointer :: xx(:)
517: PetscCall(VecGetArrayReadF90(X,xx,ierr))
519: two = 2
520: data2=reshape(xx(gdof:gdof),(/two,xe-xs+1/))
521: data=reshape(data2,(/gdof/))
522: open(1020,file='solution_out_ex22f_mf.txt')
523: do i=1,gdof
524: write(1020,'(e24.16,1x)') data(i)
525: end do
526: close(1020)
528: PetscCall(VecRestoreArrayReadF90(X,xx,ierr))
529: end subroutine SaveSolutionToDisk
531: !/*TEST
532: !
533: ! test:
534: ! args: -da_grid_x 200 -ts_arkimex_type 4
535: ! requires: !single
536: ! output_file: output/ex22f_mf_1.out
537: !
538: !TEST*/