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: use petscdm
21: implicit none
22: type AppCtx
23: PetscReal, dimension(2) :: a, k, s
24: end type AppCtx
26: PetscReal :: PETSC_SHIFT
27: TS :: tscontext
28: Mat :: Jmat
29: type(AppCtx) MFctx
31: contains
33: ! Small helper to extract the layout, result uses 1-based indexing.
34: subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
35: DM da
36: PetscInt mx, xs, xe, gxs, gxe
37: PetscErrorCode, intent(out) :: ierr
38: PetscInt xm, gxm
40: 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))
41: PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
42: PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
43: xs = xs + 1
44: gxs = gxs + 1
45: xe = xs + xm - 1
46: gxe = gxs + gxm - 1
47: end subroutine GetLayout
49: subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
50: PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
51: PetscScalar, intent(in) :: x(2, xs:xe)
52: PetscScalar, intent(in) :: xdot(2, xs:xe)
53: PetscReal, intent(in) :: a(2), k(2), s(2)
54: PetscScalar, intent(out) :: f(2, xs:xe)
55: PetscErrorCode, intent(out) :: ierr
57: f(1, :) = xdot(1, :) + k(1)*x(1, :) - k(2)*x(2, :) - s(1)
58: f(2, :) = xdot(2, :) - k(1)*x(1, :) + k(2)*x(2, :) - s(2)
59: ierr = 0
60: end subroutine FormIFunctionLocal
62: subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr)
63: TS ts
64: PetscReal t
65: Vec X, Xdot, F
66: PetscErrorCode, intent(out) :: ierr
67: type(AppCtx) ctx
69: DM da
70: PetscInt mx, xs, xe, gxs, gxe
71: PetscScalar, pointer :: xx(:), xxdot(:), ff(:)
73: PetscCall(TSGetDM(ts, da, ierr))
74: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
76: ! Get access to vector data
77: PetscCall(VecGetArrayRead(X, xx, ierr))
78: PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
79: PetscCall(VecGetArray(F, ff, ierr))
81: PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr))
83: PetscCall(VecRestoreArrayRead(X, xx, ierr))
84: PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
85: PetscCall(VecRestoreArray(F, ff, ierr))
86: end subroutine FormIFunction
88: subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
89: PetscInt mx, xs, xe, gxs, gxe
90: PetscReal, intent(in) :: t
91: PetscScalar x(2, gxs:gxe), f(2, xs:xe)
92: PetscReal a(2), k(2), s(2)
93: PetscErrorCode, intent(out) :: ierr
94: PetscInt i, j
95: PetscReal hx, u0t(2)
96: PetscReal, parameter :: twothird = 2._PETSC_REAL_KIND/3._PETSC_REAL_KIND
98: hx = 1._PETSC_REAL_KIND/mx
99: u0t = [1._PETSC_REAL_KIND - sin(12._PETSC_REAL_KIND*t)**4, 0._PETSC_REAL_KIND]
100: do i = xs, xe
101: do j = 1, 2
102: if (i == 1) then
103: f(j, i) = a(j)/hx*(u0t(j)/3._PETSC_REAL_KIND + 0.5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1) + x(j, i + 2)/6._PETSC_REAL_KIND)
104: else if (i == 2) then
105: f(j, i) = a(j)/hx*(-u0t(j)/12._PETSC_REAL_KIND + twothird*x(j, i - 1) - twothird*x(j, i + 1) + x(j, i + 2)/12._PETSC_REAL_KIND)
106: else if (i == mx - 1) then
107: f(j, i) = a(j)/hx*(-x(j, i - 2)/6._PETSC_REAL_KIND + x(j, i - 1) - 0.5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1)/3._PETSC_REAL_KIND)
108: else if (i == mx) then
109: f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
110: else
111: f(j, i) = a(j)/hx*(-x(j, i - 2)/12._PETSC_REAL_KIND + twothird*x(j, i - 1) - twothird*x(j, i + 1) + x(j, i + 2)/12._PETSC_REAL_KIND)
112: end if
113: end do
114: end do
116: #ifdef EXPLICIT_INTEGRATOR22
117: f(1, :) = f(1, :) - (k(1)*x(1, :) - k(2)*x(2, :) - s(1))
118: f(2, :) = f(2, :) - (-k(1)*x(1, :) + k(2)*x(2, :) - s(2))
119: #endif
120: ierr = 0
121: end subroutine FormRHSFunctionLocal
123: subroutine FormRHSFunction(ts, t, X, F, ctx, ierr)
124: TS ts
125: PetscReal t
126: Vec X, F
127: type(AppCtx) ctx
128: PetscErrorCode, intent(out) :: ierr
129: DM da
130: Vec Xloc
131: PetscInt mx, xs, xe, gxs, gxe
132: PetscScalar, pointer :: xx(:), ff(:)
134: PetscCall(TSGetDM(ts, da, ierr))
135: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
137: ! Scatter ghost points to local vector,using the 2-step process
138: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
139: ! By placing code between these two statements, computations can be
140: ! done while messages are in transition.
141: PetscCall(DMGetLocalVector(da, Xloc, ierr))
142: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
143: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))
145: ! Get access to vector data
146: PetscCall(VecGetArrayRead(Xloc, xx, ierr))
147: PetscCall(VecGetArray(F, ff, ierr))
149: PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr))
151: PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
152: PetscCall(VecRestoreArray(F, ff, ierr))
153: PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
154: end subroutine FormRHSFunction
156: ! ---------------------------------------------------------------------
157: !
158: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
159: !
160: subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
161: TS ts
162: PetscReal, intent(in) :: t, shift
163: Vec X, Xdot
164: Mat J, Jpre
165: type(AppCtx), intent(in) :: ctx
166: PetscErrorCode, intent(out) :: ierr
168: DM da
169: PetscInt mx, xs, xe, gxs, gxe
170: PetscInt i, row, col
171: PetscScalar val(4)
173: PetscCall(TSGetDM(ts, da, ierr))
174: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
176: val = [shift + ctx%k(1), -ctx%k(2), -ctx%k(1), shift + ctx%k(2)]
178: do i = xs, xe
179: row = i - gxs
180: col = i - gxs
181: PetscCall(MatSetValuesBlockedLocal(Jpre, 1_PETSC_INT_KIND, [row], 1_PETSC_INT_KIND, [col], val, INSERT_VALUES, ierr))
182: end do
183: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
184: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
185: if (J /= Jpre) then
186: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
187: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
188: end if
189: end subroutine FormIJacobian
191: subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
192: PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
193: PetscReal, intent(in) :: a(2), k(2), s(2)
194: PetscScalar, intent(out) :: x(2, xs:xe)
195: PetscErrorCode, intent(out) :: ierr
196: PetscInt i
197: PetscReal hx, r, ik
199: hx = 1._PETSC_REAL_KIND/mx
200: do i = xs, xe
201: r = i*hx
202: if (k(2) /= 0.0) then
203: ik = 1._PETSC_REAL_KIND/k(2)
204: else
205: ik = 1._PETSC_REAL_KIND
206: end if
207: x(1, i) = 1._PETSC_REAL_KIND + s(2)*r
208: x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
209: end do
210: ierr = 0
211: end subroutine FormInitialSolutionLocal
213: subroutine FormInitialSolution(ts, X, ctx, ierr)
214: TS ts
215: Vec X
216: type(AppCtx) ctx
217: PetscErrorCode, intent(out) :: ierr
219: DM da
220: PetscInt mx, xs, xe, gxs, gxe
221: PetscScalar, pointer :: xx(:)
223: PetscCall(TSGetDM(ts, da, ierr))
224: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
226: ! Get access to vector data
227: PetscCall(VecGetArray(X, xx, ierr))
229: PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr))
231: PetscCall(VecRestoreArray(X, xx, ierr))
232: end subroutine FormInitialSolution
234: ! ---------------------------------------------------------------------
235: !
236: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
237: !
238: subroutine FormIJacobianMF(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
239: TS ts
240: PetscReal, intent(in) :: t, shift
241: Vec X, Xdot
242: Mat J, Jpre
243: type(AppCtx) ctx
244: PetscErrorCode, intent(out) :: ierr
246: PETSC_SHIFT = shift
247: MFctx = ctx
248: ierr = 0
249: end subroutine FormIJacobianMF
251: ! -------------------------------------------------------------------
252: !
253: ! MyMult - user provided matrix multiply
254: !
255: ! Input Parameters:
256: !. X - input vector
257: !
258: ! Output Parameter:
259: !. F - function vector
260: !
261: subroutine MyMult(A, X, F, ierr)
262: Mat A
263: Vec X, F
264: PetscErrorCode, intent(out) :: ierr
266: DM da
267: PetscInt mx, xs, xe, gxs, gxe
268: PetscInt i, row, col
269: PetscScalar val(4)
271: PetscCall(TSGetDM(tscontext, da, ierr))
272: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
274: val = [PETSC_SHIFT + MFctx%k(1), -MFctx%k(2), -MFctx%k(1), PETSC_SHIFT + MFctx%k(2)]
276: do i = xs, xe
277: row = i - gxs
278: col = i - gxs
279: PetscCall(MatSetValuesBlockedLocal(Jmat, 1_PETSC_INT_KIND, [row], 1_PETSC_INT_KIND, [col], val, INSERT_VALUES, ierr))
280: end do
282: ! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
283: ! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
284: ! if (J /= Jpre) then
285: PetscCall(MatAssemblyBegin(Jmat, MAT_FINAL_ASSEMBLY, ierr))
286: PetscCall(MatAssemblyEnd(Jmat, MAT_FINAL_ASSEMBLY, ierr))
287: ! end if
289: PetscCall(MatMult(Jmat, X, F, ierr))
291: end subroutine MyMult
293: subroutine SaveSolutionToDisk(da, X, gdof, xs, xe)
294: Vec X
295: DM da
296: PetscInt xs, xe
297: PetscInt gdof, i
298: PetscErrorCode ierr
299: PetscScalar data2(2, xs:xe), data(gdof)
300: PetscScalar, pointer :: xx(:)
301: integer :: fhandle
303: PetscCall(VecGetArrayRead(X, xx, ierr))
305: data2 = reshape(xx(gdof:gdof), [2_PETSC_INT_KIND, xe - xs + 1])
306: data = reshape(data2, [gdof])
307: open (newunit=fhandle, file='solution_out_ex22f_mf.txt')
308: do i = 1, gdof
309: write (fhandle, '(e24.16,1x)') data(i)
310: end do
311: close (fhandle)
313: PetscCall(VecRestoreArrayRead(X, xx, ierr))
314: end subroutine SaveSolutionToDisk
315: end module ex22f_mfmodule
317: program main
318: use ex22f_mfmodule
319: use petscdm
320: implicit none
322: !
323: ! Create an application context to contain data needed by the
324: ! application-provided call-back routines, FormJacobian() and
325: ! FormFunction(). We use a double precision array with six
326: ! entries, two for each problem parameter a, k, s.
327: !
328: TS ts
329: Vec X
330: Mat J
331: PetscInt mx
332: PetscBool OptionSaveToDisk
333: PetscErrorCode ierr
334: DM da
335: PetscReal dt
336: PetscReal, parameter :: ftime = 1.0
337: PetscBool flg
339: PetscInt xs, xe, gxs, gxe, dof, gdof
340: PetscScalar shell_shift
341: Mat A
342: type(AppCtx) ctx
344: PetscCallA(PetscInitialize(ierr))
346: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
347: ! Create distributed array (DMDA) to manage parallel grid and vectors
348: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 1_PETSC_INT_KIND, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, da, ierr))
350: PetscCallA(DMSetFromOptions(da, ierr))
351: PetscCallA(DMSetUp(da, ierr))
353: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
354: ! Extract global vectors from DMDA
355: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: PetscCallA(DMCreateGlobalVector(da, X, ierr))
358: ! Initialize user application context
359: ! Use zero-based indexing for command line parameters to match ex22.c
360: ctx%a(1) = 1.0
361: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr))
362: ctx%a(2) = 0.0
363: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr))
364: ctx%k(1) = 1000000.0
365: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr))
366: ctx%k(2) = 2*ctx%k(1)
367: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr))
368: ctx%s(1) = 0.0
369: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr))
370: ctx%s(2) = 1.0
371: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr))
373: OptionSaveToDisk = .false.
374: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-sdisk', OptionSaveToDisk, flg, ierr))
375: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
376: ! Create timestepping solver context
377: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
378: PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
379: tscontext = ts
380: PetscCallA(TSSetDM(ts, da, ierr))
381: PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
382: PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr))
384: ! - - - - - - - - -- - - - -
385: ! Matrix free setup
386: PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
387: dof = 2_PETSC_INT_KIND*(xe - xs + 1)
388: gdof = 2_PETSC_INT_KIND*(gxe - gxs + 1)
389: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, dof, dof, gdof, gdof, shell_shift, A, ierr))
391: PetscCallA(MatShellSetOperation(A, MATOP_MULT, MyMult, ierr))
392: ! - - - - - - - - - - - -
394: PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr))
395: PetscCallA(DMSetMatType(da, MATAIJ, ierr))
396: PetscCallA(DMCreateMatrix(da, J, ierr))
398: Jmat = J
400: PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr))
401: PetscCallA(TSSetIJacobian(ts, A, A, FormIJacobianMF, ctx, ierr))
403: PetscCallA(TSSetMaxTime(ts, ftime, ierr))
404: PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))
406: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407: ! Set initial conditions
408: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409: PetscCallA(FormInitialSolution(ts, X, ctx, ierr))
410: PetscCallA(TSSetSolution(ts, X, ierr))
411: PetscCallA(VecGetSize(X, mx, ierr))
412: ! Advective CFL, I don't know why it needs so much safety factor.
413: dt = max(ctx%a(1), ctx%a(2))/mx/10._PETSC_REAL_KIND
414: PetscCallA(TSSetTimeStep(ts, dt, ierr))
416: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
417: ! Set runtime options
418: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
419: PetscCallA(TSSetFromOptions(ts, ierr))
421: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
422: ! Solve nonlinear system
423: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
424: PetscCallA(TSSolve(ts, X, ierr))
426: if (OptionSaveToDisk) then
427: PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
428: dof = 2_PETSC_INT_KIND*(xe - xs + 1)
429: gdof = 2_PETSC_INT_KIND*(gxe - gxs + 1)
430: call SaveSolutionToDisk(da, X, gdof, xs, xe)
431: end if
433: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434: ! Free work space.
435: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436: PetscCallA(MatDestroy(A, ierr))
437: PetscCallA(MatDestroy(J, ierr))
438: PetscCallA(VecDestroy(X, ierr))
439: PetscCallA(TSDestroy(ts, ierr))
440: PetscCallA(DMDestroy(da, ierr))
441: PetscCallA(PetscFinalize(ierr))
442: end program main
444: !/*TEST
445: !
446: ! test:
447: ! args: -da_grid_x 200 -ts_arkimex_type 4
448: ! requires: !single
449: ! output_file: output/empty.out
450: !
451: !TEST*/