Actual source code: ex22f.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: #include <petsc/finclude/petscts.h>
18: #include <petsc/finclude/petscdmda.h>
20: module ex22f_modctx
21: use petscts
22: use petscdm
23: implicit none
24: type AppCtx
25: PetscReal a(2), k(2), s(2)
26: end type AppCtx
27: contains
29: ! Small helper to extract the layout, result uses 1-based indexing.
30: subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
31: DM da
32: PetscInt mx, xs, xe, gxs, gxe
33: PetscErrorCode, intent(out) :: ierr
34: PetscInt xm, gxm
35: 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))
36: PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
37: PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
38: xs = xs + 1
39: gxs = gxs + 1
40: xe = xs + xm - 1
41: gxe = gxs + gxm - 1
42: end subroutine
44: subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
45: PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
46: PetscScalar, dimension(2, xs:xe), intent(in) :: x, xdot
47: PetscScalar, dimension(2, xs:xe), intent(out) :: f
48: PetscReal, dimension(2) :: a, k, s
49: PetscErrorCode, intent(out) :: ierr
51: f(1, :) = xdot(1, :) + k(1)*x(1, :) - k(2)*x(2, :) - s(1)
52: f(2, :) = xdot(2, :) - k(1)*x(1, :) + k(2)*x(2, :) - s(2)
53: ierr = 0
54: end subroutine
56: subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr)
57: TS ts
58: type(AppCtx) ctx
59: PetscReal t
60: Vec X, Xdot, F
61: PetscErrorCode, intent(out) :: ierr
63: DM da
64: PetscInt mx, xs, xe, gxs, gxe
65: PetscScalar, pointer :: xx(:), xxdot(:), ff(:)
67: PetscCall(TSGetDM(ts, da, ierr))
68: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
70: ! Get access to vector data
71: PetscCall(VecGetArrayRead(X, xx, ierr))
72: PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
73: PetscCall(VecGetArray(F, ff, ierr))
75: PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr))
77: PetscCall(VecRestoreArrayRead(X, xx, ierr))
78: PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
79: PetscCall(VecRestoreArray(F, ff, ierr))
80: end subroutine
82: subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
83: PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
84: PetscReal, intent(in) :: t
85: PetscScalar x(2, gxs:gxe), f(2, xs:xe)
86: PetscReal, dimension(2) :: a, k, s, u0t
87: PetscErrorCode ierr
88: PetscInt i, j
89: PetscReal hx
90: PetscReal, parameter :: twelfth = 1._PETSC_REAL_KIND/12._PETSC_REAL_KIND, twothird = 2._PETSC_REAL_KIND/3._PETSC_REAL_KIND
92: hx = 1.0_PETSC_REAL_KIND/mx
93: ! The Fortran standard only allows positive base for power functions; Nag compiler fails on this
94: u0t = [1.0_PETSC_REAL_KIND - abs(sin(12._PETSC_REAL_KIND*t))**4, 0._PETSC_REAL_KIND]
95: do i = xs, xe
96: do j = 1, 2
97: if (i == 1) then
98: f(j, i) = a(j)/hx*(u0t(j)/3._PETSC_REAL_KIND + .5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1) + x(j, i + 2)/6._PETSC_REAL_KIND)
99: else if (i == 2) then
100: f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
101: else if (i == mx - 1) then
102: f(j, i) = a(j)/hx*(-x(j, i - 2)/6._PETSC_REAL_KIND + x(j, i - 1) - .5_PETSC_REAL_KIND*x(j, i) - x(j, i + 1)/3._PETSC_REAL_KIND)
103: else if (i == mx) then
104: f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
105: else
106: 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))
107: end if
108: end do
109: end do
110: end subroutine
112: subroutine FormRHSFunction(ts, t, X, F, ctx, ierr)
113: type(AppCtx) ctx
114: TS ts
115: PetscReal t
116: Vec X, F
117: PetscErrorCode, intent(out) :: ierr
118: DM da
119: Vec Xloc
120: PetscInt mx, xs, xe, gxs, gxe
121: PetscScalar, pointer :: xx(:), ff(:)
123: PetscCall(TSGetDM(ts, da, ierr))
124: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
126: ! Scatter ghost points to local vector,using the 2-step process
127: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
128: ! By placing code between these two statements, computations can be
129: ! done while messages are in transition.
130: PetscCall(DMGetLocalVector(da, Xloc, ierr))
131: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
132: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))
134: ! Get access to vector data
135: PetscCall(VecGetArrayRead(Xloc, xx, ierr))
136: PetscCall(VecGetArray(F, ff, ierr))
138: PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr))
140: PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
141: PetscCall(VecRestoreArray(F, ff, ierr))
142: PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
143: end subroutine
145: ! ---------------------------------------------------------------------
146: !
147: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
148: !
149: subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
150: type(AppCtx) ctx
151: TS ts
152: PetscReal t, shift
153: Vec X, Xdot
154: Mat J, Jpre
155: PetscErrorCode ierr
157: DM da
158: PetscInt mx, xs, xe, gxs, gxe
159: PetscInt i, row, col
160: PetscReal k1, k2
161: PetscScalar val(4)
163: PetscCall(TSGetDM(ts, da, ierr))
164: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
166: k1 = ctx%k(1)
167: k2 = ctx%k(2)
168: do i = xs, xe
169: row = i - gxs
170: col = i - gxs
171: val = [shift + k1, -k2, -k1, shift + k2]
172: PetscCall(MatSetValuesBlockedLocal(Jpre, 1_PETSC_INT_KIND, [row], 1_PETSC_INT_KIND, [col], val, INSERT_VALUES, ierr))
173: end do
174: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
175: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
176: if (J /= Jpre) then
177: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
178: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
179: end if
180: end subroutine
182: subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
183: PetscInt, intent(in) :: mx, xs, xe, gxs, gxe
184: PetscScalar, intent(out) :: x(2, xs:xe)
185: PetscReal, dimension(2) :: a, k, s
186: PetscErrorCode, intent(out) :: ierr
187: PetscInt i
188: PetscReal hx, r, ik
190: hx = 1._PETSC_REAL_KIND/mx
191: do i = xs, xe
192: r = i*hx
193: if (k(2) /= 0.0) then
194: ik = 1._PETSC_REAL_KIND/k(2)
195: else
196: ik = 1._PETSC_REAL_KIND
197: end if
198: x(1, i) = 1._PETSC_REAL_KIND + s(2)*r
199: x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
200: end do
201: ierr = 0
202: end subroutine
204: subroutine FormInitialSolution(ts, X, ctx, ierr)
205: type(AppCtx) ctx
206: TS ts
207: Vec X
208: PetscErrorCode, intent(out) :: ierr
209: DM da
210: PetscInt mx, xs, xe, gxs, gxe
211: PetscScalar, pointer :: xx(:)
213: PetscCall(TSGetDM(ts, da, ierr))
214: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
216: ! Get access to vector data
217: PetscCall(VecGetArray(X, xx, ierr))
219: PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr))
221: PetscCall(VecRestoreArray(X, xx, ierr))
222: end subroutine
224: end module ex22f_modctx
225: program main
226: use ex22f_modctx
227: implicit none
229: !
230: ! Create an application context to contain data needed by the
231: ! application-provided call-back routines, FormJacobian() and
232: ! FormFunction(). We use a double precision array with six
233: ! entries, two for each problem parameter a, k, s.
234: !
236: TS ts
237: SNES snes
238: SNESLineSearch linesearch
239: Vec X
240: Mat J
241: PetscInt mx
242: PetscErrorCode ierr
243: DM da
244: PetscReal, parameter :: ftime = 1.0
245: PetscReal dt
246: PetscBool flg
247: type(AppCtx) ctx
249: PetscCallA(PetscInitialize(ierr))
251: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: ! Create distributed array (DMDA) to manage parallel grid and vectors
253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11_PETSC_INT_KIND, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER, da, ierr))
255: PetscCallA(DMSetFromOptions(da, ierr))
256: PetscCallA(DMSetUp(da, ierr))
258: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259: ! Extract global vectors from DMDA
260: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
261: PetscCallA(DMCreateGlobalVector(da, X, ierr))
263: ! Initialize user application context
264: ! Use zero-based indexing for command line parameters to match ex22.c
265: ctx%a(1) = 1.0
266: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr))
267: ctx%a(2) = 0.0
268: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr))
269: ctx%k(1) = 1000000.0
270: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr))
271: ctx%k(2) = 2*ctx%k(1)
272: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr))
273: ctx%s(1) = 0.0
274: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr))
275: ctx%s(2) = 1.0
276: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr))
278: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
279: ! Create timestepping solver context
280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
282: PetscCallA(TSSetDM(ts, da, ierr))
283: PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
284: PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr))
285: PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr))
286: PetscCallA(DMSetMatType(da, MATAIJ, ierr))
287: PetscCallA(DMCreateMatrix(da, J, ierr))
288: PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr))
290: PetscCallA(TSGetSNES(ts, snes, ierr))
291: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
292: PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))
294: PetscCallA(TSSetMaxTime(ts, ftime, ierr))
295: PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))
297: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
298: ! Set initial conditions
299: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300: PetscCallA(FormInitialSolution(ts, X, ctx, ierr))
301: PetscCallA(TSSetSolution(ts, X, ierr))
302: PetscCallA(VecGetSize(X, mx, ierr))
303: ! Advective CFL, I don't know why it needs so much safety factor.
304: dt = .1_PETSC_REAL_KIND*max(ctx%a(1), ctx%a(2))/mx
305: PetscCallA(TSSetTimeStep(ts, dt, ierr))
307: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308: ! Set runtime options
309: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310: PetscCallA(TSSetFromOptions(ts, ierr))
312: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313: ! Solve nonlinear system
314: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
315: PetscCallA(TSSolve(ts, X, ierr))
317: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318: ! Free work space.
319: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: PetscCallA(MatDestroy(J, ierr))
321: PetscCallA(VecDestroy(X, ierr))
322: PetscCallA(TSDestroy(ts, ierr))
323: PetscCallA(DMDestroy(da, ierr))
324: PetscCallA(PetscFinalize(ierr))
325: end program
326: !/*TEST
327: !
328: ! test:
329: ! args: -da_grid_x 200 -ts_arkimex_type 4
330: ! requires: !single
331: ! output_file: output/empty.out
332: !
333: !TEST*/