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: type AppCtx
23: PetscReal a(2), k(2), s(2)
24: end type AppCtx
26: end module ex22f_modctx
27: program main
28: use ex22f_modctx
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: !
38: TS ts
39: SNES snes
40: SNESLineSearch linesearch
41: Vec X
42: Mat J
43: PetscInt mx
44: PetscErrorCode ierr
45: DM da
46: PetscReal ftime, dt
47: PetscReal one, pone
48: PetscInt im11, i2
49: PetscBool flg
50: type(AppCtx) ctx
52: im11 = 11
53: i2 = 2
54: one = 1.0
55: pone = one/10
57: PetscCallA(PetscInitialize(ierr))
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: ! Create distributed array (DMDA) to manage parallel grid and vectors
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, im11, i2, i2, PETSC_NULL_INTEGER, da, ierr))
63: PetscCallA(DMSetFromOptions(da, ierr))
64: PetscCallA(DMSetUp(da, ierr))
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: ! Extract global vectors from DMDA
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: PetscCallA(DMCreateGlobalVector(da, X, ierr))
71: ! Initialize user application context
72: ! Use zero-based indexing for command line parameters to match ex22.c
73: ctx%a(1) = 1.0
74: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr))
75: ctx%a(2) = 0.0
76: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr))
77: ctx%k(1) = 1000000.0
78: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr))
79: ctx%k(2) = 2*ctx%k(1)
80: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr))
81: ctx%s(1) = 0.0
82: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr))
83: ctx%s(2) = 1.0
84: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr))
86: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: ! Create timestepping solver context
88: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
90: PetscCallA(TSSetDM(ts, da, ierr))
91: PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
92: PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr))
93: PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr))
94: PetscCallA(DMSetMatType(da, MATAIJ, ierr))
95: PetscCallA(DMCreateMatrix(da, J, ierr))
96: PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr))
98: PetscCallA(TSGetSNES(ts, snes, ierr))
99: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
100: PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))
102: ftime = 1.0
103: PetscCallA(TSSetMaxTime(ts, ftime, ierr))
104: PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Set initial conditions
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: PetscCallA(FormInitialSolution(ts, X, ctx, ierr))
110: PetscCallA(TSSetSolution(ts, X, ierr))
111: PetscCallA(VecGetSize(X, mx, ierr))
112: ! Advective CFL, I don't know why it needs so much safety factor.
113: dt = pone*max(ctx%a(1), ctx%a(2))/mx
114: PetscCallA(TSSetTimeStep(ts, dt, ierr))
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: ! Set runtime options
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: PetscCallA(TSSetFromOptions(ts, ierr))
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! Solve nonlinear system
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: PetscCallA(TSSolve(ts, X, ierr))
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! Free work space.
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: PetscCallA(MatDestroy(J, ierr))
130: PetscCallA(VecDestroy(X, ierr))
131: PetscCallA(TSDestroy(ts, ierr))
132: PetscCallA(DMDestroy(da, ierr))
133: PetscCallA(PetscFinalize(ierr))
134: contains
136: ! Small helper to extract the layout, result uses 1-based indexing.
137: subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
138: use petscdm
139: implicit none
141: DM da
142: PetscInt mx, xs, xe, gxs, gxe
143: PetscErrorCode ierr
144: PetscInt xm, gxm
145: 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))
146: PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
147: PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
148: xs = xs + 1
149: gxs = gxs + 1
150: xe = xs + xm - 1
151: gxe = gxs + gxm - 1
152: end subroutine
154: subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
155: implicit none
156: PetscInt mx, xs, xe, gxs, gxe
157: PetscScalar x(2, xs:xe)
158: PetscScalar xdot(2, xs:xe)
159: PetscScalar f(2, xs:xe)
160: PetscReal a(2), k(2), s(2)
161: PetscErrorCode ierr
162: PetscInt i
163: do i = xs, xe
164: f(1, i) = xdot(1, i) + k(1)*x(1, i) - k(2)*x(2, i) - s(1)
165: f(2, i) = xdot(2, i) - k(1)*x(1, i) + k(2)*x(2, i) - s(2)
166: end do
167: end subroutine
169: subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr)
170: use ex22f_modctx
171: implicit none
173: TS ts
174: type(AppCtx) ctx
175: PetscReal t
176: Vec X, Xdot, F
177: PetscErrorCode ierr
179: DM da
180: PetscInt mx, xs, xe, gxs, gxe
181: PetscScalar, pointer :: xx(:), xxdot(:), ff(:)
183: PetscCall(TSGetDM(ts, da, ierr))
184: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
186: ! Get access to vector data
187: PetscCall(VecGetArrayRead(X, xx, ierr))
188: PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
189: PetscCall(VecGetArray(F, ff, ierr))
191: PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr))
193: PetscCall(VecRestoreArrayRead(X, xx, ierr))
194: PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
195: PetscCall(VecRestoreArray(F, ff, ierr))
196: end subroutine
198: subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
199: implicit none
200: PetscInt mx, xs, xe, gxs, gxe
201: PetscReal t
202: PetscScalar x(2, gxs:gxe), f(2, xs:xe)
203: PetscReal a(2), k(2), s(2)
204: PetscErrorCode ierr
205: PetscInt i, j
206: PetscReal hx, u0t(2)
207: PetscReal one, two, three, four, six, twelve
208: PetscReal half, third, twothird, sixth
209: PetscReal twelfth
211: one = 1.0
212: two = 2.0
213: three = 3.0
214: four = 4.0
215: six = 6.0
216: twelve = 12.0
217: hx = one/mx
218: ! The Fortran standard only allows positive base for power functions; Nag compiler fails on this
219: u0t(1) = one - abs(sin(twelve*t))**four
220: u0t(2) = 0.0
221: half = one/two
222: third = one/three
223: twothird = two/three
224: sixth = one/six
225: twelfth = one/twelve
226: do i = xs, xe
227: do j = 1, 2
228: if (i == 1) then
229: f(j, i) = a(j)/hx*(third*u0t(j) + half*x(j, i) - x(j, i + 1) &
230: & + sixth*x(j, i + 2))
231: else if (i == 2) then
232: f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) &
233: & - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
234: else if (i == mx - 1) then
235: f(j, i) = a(j)/hx*(-sixth*x(j, i - 2) + x(j, i - 1) &
236: & - half*x(j, i) - third*x(j, i + 1))
237: else if (i == mx) then
238: f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
239: else
240: f(j, i) = a(j)/hx*(-twelfth*x(j, i - 2) &
241: & + twothird*x(j, i - 1) &
242: & - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
243: end if
244: end do
245: end do
246: end subroutine
248: subroutine FormRHSFunction(ts, t, X, F, ctx, ierr)
249: use ex22f_modctx
250: implicit none
252: type(AppCtx) ctx
253: TS ts
254: PetscReal t
255: Vec X, F
256: PetscErrorCode ierr
257: DM da
258: Vec Xloc
259: PetscInt mx, xs, xe, gxs, gxe
260: PetscScalar, pointer :: xx(:), ff(:)
262: PetscCall(TSGetDM(ts, da, ierr))
263: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
265: ! Scatter ghost points to local vector,using the 2-step process
266: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
267: ! By placing code between these two statements, computations can be
268: ! done while messages are in transition.
269: PetscCall(DMGetLocalVector(da, Xloc, ierr))
270: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
271: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))
273: ! Get access to vector data
274: PetscCall(VecGetArrayRead(Xloc, xx, ierr))
275: PetscCall(VecGetArray(F, ff, ierr))
277: PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr))
279: PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
280: PetscCall(VecRestoreArray(F, ff, ierr))
281: PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
282: end subroutine
284: ! ---------------------------------------------------------------------
285: !
286: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
287: !
288: subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr)
289: use ex22f_modctx
290: implicit none
292: type(AppCtx) ctx
293: TS ts
294: PetscReal t, shift
295: Vec X, Xdot
296: Mat J, Jpre
297: PetscErrorCode ierr
299: DM da
300: PetscInt mx, xs, xe, gxs, gxe
301: PetscInt i, i1, row, col
302: PetscReal k1, k2
303: PetscScalar val(4)
305: PetscCall(TSGetDM(ts, da, ierr))
306: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
308: i1 = 1
309: k1 = ctx%k(1)
310: k2 = ctx%k(2)
311: do i = xs, xe
312: row = i - gxs
313: col = i - gxs
314: val(1) = shift + k1
315: val(2) = -k2
316: val(3) = -k1
317: val(4) = shift + k2
318: PetscCall(MatSetValuesBlockedLocal(Jpre, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
319: end do
320: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
321: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
322: if (J /= Jpre) then
323: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
324: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
325: end if
326: end subroutine
328: subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
329: implicit none
330: PetscInt mx, xs, xe, gxs, gxe
331: PetscScalar x(2, xs:xe)
332: PetscReal a(2), k(2), s(2)
333: PetscErrorCode ierr
335: PetscInt i
336: PetscReal one, hx, r, ik
337: one = 1.0
338: hx = one/mx
339: do i = xs, xe
340: r = i*hx
341: if (k(2) /= 0.0) then
342: ik = one/k(2)
343: else
344: ik = one
345: end if
346: x(1, i) = one + s(2)*r
347: x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
348: end do
349: end subroutine
351: subroutine FormInitialSolution(ts, X, ctx, ierr)
352: use ex22f_modctx
353: implicit none
355: type(AppCtx) ctx
356: TS ts
357: Vec X
358: PetscErrorCode ierr
360: DM da
361: PetscInt mx, xs, xe, gxs, gxe
362: PetscScalar, pointer :: xx(:)
364: PetscCall(TSGetDM(ts, da, ierr))
365: PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
367: ! Get access to vector data
368: PetscCall(VecGetArray(X, xx, ierr))
370: PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr))
372: PetscCall(VecRestoreArray(X, xx, ierr))
373: end subroutine
374: end program
375: !/*TEST
376: !
377: ! test:
378: ! args: -da_grid_x 200 -ts_arkimex_type 4
379: ! requires: !single
380: ! output_file: output/empty.out
381: !
382: !TEST*/