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