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: !
16: #include <petsc/finclude/petscts.h>
18: module ex22f_mfmodule
19: use petscts
20: type AppCtx
21: PetscReal a(2), k(2), s(2)
22: end type AppCtx
24: PetscScalar::PETSC_SHIFT
25: TS::tscontext
26: Mat::Jmat
27: type(AppCtx) MFctx
28: end module ex22f_mfmodule
30: program main
31: use ex22f_mfmodule
32: use petscdm
33: implicit none
35: !
36: ! Create an application context to contain data needed by the
37: ! application-provided call-back routines, FormJacobian() and
38: ! FormFunction(). We use a double precision array with six
39: ! entries, two for each problem parameter a, k, s.
40: !
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
56: type(AppCtx) ctx
58: im11 = 11
59: i2 = 2
60: one = 1.0
61: pone = one/10
63: PetscCallA(PetscInitialize(ierr))
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: ! Create distributed array (DMDA) to manage parallel grid and vectors
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, im11, i2, i2, PETSC_NULL_INTEGER_ARRAY, da, ierr))
69: PetscCallA(DMSetFromOptions(da, ierr))
70: PetscCallA(DMSetUp(da, ierr))
72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: ! Extract global vectors from DMDA
74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: PetscCallA(DMCreateGlobalVector(da, X, ierr))
77: ! Initialize user application context
78: ! Use zero-based indexing for command line parameters to match ex22.c
79: ctx%a(1) = 1.0
80: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr))
81: ctx%a(2) = 0.0
82: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr))
83: ctx%k(1) = 1000000.0
84: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr))
85: ctx%k(2) = 2*ctx%k(1)
86: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr))
87: ctx%s(1) = 0.0
88: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr))
89: ctx%s(2) = 1.0
90: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr))
92: OptionSaveToDisk = .false.
93: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-sdisk', OptionSaveToDisk, flg, ierr))
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Create timestepping solver context
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
98: tscontext = ts
99: PetscCallA(TSSetDM(ts, da, ierr))
100: PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
101: PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr))
103: ! - - - - - - - - -- - - - -
104: ! Matrix free setup
105: PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
106: dof = i2*(xe - xs + 1)
107: gdof = i2*(gxe - gxs + 1)
108: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, dof, dof, gdof, gdof, shell_shift, A, ierr))
110: PetscCallA(MatShellSetOperation(A, MATOP_MULT, MyMult, ierr))
111: ! - - - - - - - - - - - -
113: PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr))
114: PetscCallA(DMSetMatType(da, MATAIJ, ierr))
115: PetscCallA(DMCreateMatrix(da, J, ierr))
117: Jmat = J
119: PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr))
120: PetscCallA(TSSetIJacobian(ts, A, A, FormIJacobianMF, ctx, ierr))
122: ftime = 1.0
123: PetscCallA(TSSetMaxTime(ts, ftime, ierr))
124: PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! Set initial conditions
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: PetscCallA(FormInitialSolution(ts, X, ctx, ierr))
130: PetscCallA(TSSetSolution(ts, X, ierr))
131: PetscCallA(VecGetSize(X, mx, ierr))
132: ! Advective CFL, I don't know why it needs so much safety factor.
133: dt = pone*max(ctx%a(1), ctx%a(2))/mx
134: PetscCallA(TSSetTimeStep(ts, dt, ierr))
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Set runtime options
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: PetscCallA(TSSetFromOptions(ts, ierr))
141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: ! Solve nonlinear system
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: PetscCallA(TSSolve(ts, X, ierr))
146: if (OptionSaveToDisk) then
147: PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
148: dof = i2*(xe - xs + 1)
149: gdof = i2*(gxe - gxs + 1)
150: call SaveSolutionToDisk(da, X, gdof, xs, xe)
151: end if
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Free work space.
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: PetscCallA(MatDestroy(A, ierr))
157: PetscCallA(MatDestroy(J, ierr))
158: PetscCallA(VecDestroy(X, ierr))
159: PetscCallA(TSDestroy(ts, ierr))
160: PetscCallA(DMDestroy(da, ierr))
161: PetscCallA(PetscFinalize(ierr))
162: contains
164: ! Small helper to extract the layout, result uses 1-based indexing.
165: subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
166: use petscdm
167: implicit none
169: DM da
170: PetscInt mx, xs, xe, gxs, gxe
171: PetscErrorCode ierr
172: PetscInt xm, gxm
173: 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))
174: PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
175: PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
176: xs = xs + 1
177: gxs = gxs + 1
178: xe = xs + xm - 1
179: gxe = gxs + gxm - 1
180: end subroutine GetLayout
182: subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
183: implicit none
184: PetscInt mx, xs, xe, gxs, gxe
185: PetscScalar x(2, xs:xe)
186: PetscScalar xdot(2, xs:xe)
187: PetscScalar f(2, xs:xe)
188: PetscReal a(2), k(2), s(2)
189: PetscErrorCode ierr
190: PetscInt i
191: do i = xs, xe
192: f(1, i) = xdot(1, i) + k(1)*x(1, i) - k(2)*x(2, i) - s(1)
193: f(2, i) = xdot(2, i) - k(1)*x(1, i) + k(2)*x(2, i) - s(2)
194: end do
195: end subroutine FormIFunctionLocal
197: subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr)
198: use petscdm
199: use ex22f_mfmodule
200: implicit none
202: TS ts
203: PetscReal t
204: Vec X, Xdot, F
205: PetscErrorCode ierr
206: type(AppCtx) ctx
208: DM da
209: PetscInt mx, xs, xe, gxs, gxe
210: PetscScalar, pointer :: xx(:), xxdot(:), ff(:)
212: PetscCall(TSGetDM(ts, da, ierr))
213: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
215: ! Get access to vector data
216: PetscCall(VecGetArrayRead(X, xx, ierr))
217: PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
218: PetscCall(VecGetArray(F, ff, ierr))
220: PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr))
222: PetscCall(VecRestoreArrayRead(X, xx, ierr))
223: PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
224: PetscCall(VecRestoreArray(F, ff, ierr))
225: end subroutine FormIFunction
227: subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
228: implicit none
229: PetscInt mx, xs, xe, gxs, gxe
230: PetscReal t
231: PetscScalar x(2, gxs:gxe), f(2, xs:xe)
232: PetscReal a(2), k(2), s(2)
233: PetscErrorCode ierr
234: PetscInt i, j
235: PetscReal hx, u0t(2)
236: PetscReal one, two, three, four, six, twelve
237: PetscReal half, third, twothird, sixth
238: PetscReal twelfth
240: one = 1.0
241: two = 2.0
242: three = 3.0
243: four = 4.0
244: six = 6.0
245: twelve = 12.0
246: hx = one/mx
247: u0t(1) = one - sin(twelve*t)**four
248: u0t(2) = 0.0
249: half = one/two
250: third = one/three
251: twothird = two/three
252: sixth = one/six
253: twelfth = one/twelve
254: do i = xs, xe
255: do j = 1, 2
256: if (i == 1) then
257: f(j, i) = a(j)/hx*(third*u0t(j) + half*x(j, i) - x(j, i + 1) + sixth*x(j, i + 2))
258: else if (i == 2) then
259: f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
260: else if (i == mx - 1) then
261: f(j, i) = a(j)/hx*(-sixth*x(j, i - 2) + x(j, i - 1) - half*x(j, i) - third*x(j, i + 1))
262: else if (i == mx) then
263: f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
264: else
265: 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))
266: end if
267: end do
268: end do
270: #ifdef EXPLICIT_INTEGRATOR22
271: do i = xs, xe
272: f(1, i) = f(1, i) - (k(1)*x(1, i) - k(2)*x(2, i) - s(1))
273: f(2, i) = f(2, i) - (-k(1)*x(1, i) + k(2)*x(2, i) - s(2))
274: end do
275: #endif
277: end subroutine FormRHSFunctionLocal
279: subroutine FormRHSFunction(ts, t, X, F, ctx, ierr)
280: use ex22f_mfmodule
281: implicit none
283: TS ts
284: PetscReal t
285: Vec X, F
286: type(AppCtx) ctx
287: PetscErrorCode ierr
288: DM da
289: Vec Xloc
290: PetscInt mx, xs, xe, gxs, gxe
291: PetscScalar, pointer :: xx(:), ff(:)
293: PetscCall(TSGetDM(ts, da, ierr))
294: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
296: ! Scatter ghost points to local vector,using the 2-step process
297: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
298: ! By placing code between these two statements, computations can be
299: ! done while messages are in transition.
300: PetscCall(DMGetLocalVector(da, Xloc, ierr))
301: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
302: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))
304: ! Get access to vector data
305: PetscCall(VecGetArrayRead(Xloc, xx, ierr))
306: PetscCall(VecGetArray(F, ff, ierr))
308: PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr))
310: PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
311: PetscCall(VecRestoreArray(F, ff, ierr))
312: PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
313: end subroutine FormRHSFunction
315: ! ---------------------------------------------------------------------
316: !
317: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
318: !
319: subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
320: use ex22f_mfmodule
321: use petscdm
322: implicit none
324: TS ts
325: PetscReal t, shift
326: Vec X, Xdot
327: Mat J, Jpre
328: type(AppCtx) ctx
329: PetscErrorCode ierr
331: DM da
332: PetscInt mx, xs, xe, gxs, gxe
333: PetscInt i, i1, row, col
334: PetscReal k1, k2
335: PetscScalar val(4)
337: PetscCall(TSGetDM(ts, da, ierr))
338: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
340: i1 = 1
341: k1 = ctx%k(1)
342: k2 = ctx%k(2)
343: do i = xs, xe
344: row = i - gxs
345: col = i - gxs
346: val(1) = shift + k1
347: val(2) = -k2
348: val(3) = -k1
349: val(4) = shift + k2
350: PetscCall(MatSetValuesBlockedLocal(Jpre, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
351: end do
352: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
353: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
354: if (J /= Jpre) then
355: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
356: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
357: end if
358: end subroutine FormIJacobian
360: subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
361: implicit none
362: PetscInt mx, xs, xe, gxs, gxe
363: PetscScalar x(2, xs:xe)
364: PetscReal a(2), k(2), s(2)
365: PetscErrorCode ierr
367: PetscInt i
368: PetscReal one, hx, r, ik
369: one = 1.0
370: hx = one/mx
371: do i = xs, xe
372: r = i*hx
373: if (k(2) /= 0.0) then
374: ik = one/k(2)
375: else
376: ik = one
377: end if
378: x(1, i) = one + s(2)*r
379: x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
380: end do
381: end subroutine FormInitialSolutionLocal
383: subroutine FormInitialSolution(ts, X, ctx, ierr)
384: use ex22f_mfmodule
385: use petscdm
386: implicit none
388: TS ts
389: Vec X
390: type(AppCtx) ctx
391: PetscErrorCode ierr
393: DM da
394: PetscInt mx, xs, xe, gxs, gxe
395: PetscScalar, pointer :: xx(:)
397: PetscCall(TSGetDM(ts, da, ierr))
398: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
400: ! Get access to vector data
401: PetscCall(VecGetArray(X, xx, ierr))
403: PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr))
405: PetscCall(VecRestoreArray(X, xx, ierr))
406: end subroutine FormInitialSolution
408: ! ---------------------------------------------------------------------
409: !
410: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
411: !
412: subroutine FormIJacobianMF(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
413: use ex22f_mfmodule
414: implicit none
415: TS ts
416: PetscReal t, shift
417: Vec X, Xdot
418: Mat J, Jpre
419: type(AppCtx) ctx
420: PetscErrorCode ierr
422: PETSC_SHIFT = shift
423: MFctx = ctx
425: end subroutine FormIJacobianMF
427: ! -------------------------------------------------------------------
428: !
429: ! MyMult - user provided matrix multiply
430: !
431: ! Input Parameters:
432: !. X - input vector
433: !
434: ! Output Parameter:
435: !. F - function vector
436: !
437: subroutine MyMult(A, X, F, ierr)
438: use ex22f_mfmodule
439: implicit none
441: Mat A
442: Vec X, F
444: PetscErrorCode ierr
445: PetscScalar shift
446: type(AppCtx) ctx
448: DM da
449: PetscInt mx, xs, xe, gxs, gxe
450: PetscInt i, i1, row, col
451: PetscReal k1, k2
452: PetscScalar val(4)
454: shift = PETSC_SHIFT
455: ctx = MFctx
457: PetscCall(TSGetDM(tscontext, da, ierr))
458: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
460: i1 = 1
461: k1 = ctx%k(1)
462: k2 = ctx%k(2)
464: do i = xs, xe
465: row = i - gxs
466: col = i - gxs
467: val(1) = shift + k1
468: val(2) = -k2
469: val(3) = -k1
470: val(4) = shift + k2
471: PetscCall(MatSetValuesBlockedLocal(Jmat, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
472: end do
474: ! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
475: ! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
476: ! if (J /= Jpre) then
477: PetscCall(MatAssemblyBegin(Jmat, MAT_FINAL_ASSEMBLY, ierr))
478: PetscCall(MatAssemblyEnd(Jmat, MAT_FINAL_ASSEMBLY, ierr))
479: ! end if
481: PetscCall(MatMult(Jmat, X, F, ierr))
483: end subroutine MyMult
485: !
486: subroutine SaveSolutionToDisk(da, X, gdof, xs, xe)
487: use petscdm
488: implicit none
490: Vec X
491: DM da
492: PetscInt xs, xe, two
493: PetscInt gdof, i
494: PetscErrorCode ierr
495: PetscScalar data2(2, xs:xe), data(gdof)
496: PetscScalar, pointer :: xx(:)
498: PetscCall(VecGetArrayRead(X, xx, ierr))
500: two = 2
501: data2 = reshape(xx(gdof:gdof), (/two, xe - xs + 1/))
502: data = reshape(data2, (/gdof/))
503: open (1020, file='solution_out_ex22f_mf.txt')
504: do i = 1, gdof
505: write (1020, '(e24.16,1x)') data(i)
506: end do
507: close (1020)
509: PetscCall(VecRestoreArrayRead(X, xx, ierr))
510: end subroutine SaveSolutionToDisk
511: end program main
513: !/*TEST
514: !
515: ! test:
516: ! args: -da_grid_x 200 -ts_arkimex_type 4
517: ! requires: !single
518: ! output_file: output/empty.out
519: !
520: !TEST*/