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: program main
18: #include <petsc/finclude/petscts.h>
19: #include <petsc/finclude/petscdmda.h>
20: use petscts
21: implicit none
23: !
24: ! Create an application context to contain data needed by the
25: ! application-provided call-back routines, FormJacobian() and
26: ! FormFunction(). We use a double precision array with six
27: ! entries, two for each problem parameter a, k, s.
28: !
29: PetscReal user(6)
30: integer user_a,user_k,user_s
31: parameter (user_a = 0,user_k = 2,user_s = 4)
33: TS ts
34: SNES snes
35: SNESLineSearch linesearch
36: Vec X
37: Mat J
38: PetscInt mx
39: PetscErrorCode ierr
40: DM da
41: PetscReal ftime,dt
42: PetscReal one,pone
43: PetscInt im11,i2
44: PetscBool flg
46: im11 = 11
47: i2 = 2
48: one = 1.0
49: pone = one / 10
51: PetscCallA(PetscInitialize(ierr))
53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: ! Create distributed array (DMDA) to manage parallel grid and vectors
55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr))
57: PetscCallA(DMSetFromOptions(da,ierr))
58: PetscCallA(DMSetUp(da,ierr))
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Extract global vectors from DMDA
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: PetscCallA(DMCreateGlobalVector(da,X,ierr))
65: ! Initialize user application context
66: ! Use zero-based indexing for command line parameters to match ex22.c
67: user(user_a+1) = 1.0
68: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr))
69: user(user_a+2) = 0.0
70: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr))
71: user(user_k+1) = 1000000.0
72: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0',user(user_k+1),flg,ierr))
73: user(user_k+2) = 2*user(user_k+1)
74: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr))
75: user(user_s+1) = 0.0
76: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr))
77: user(user_s+2) = 1.0
78: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr))
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: ! Create timestepping solver context
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
84: PetscCallA(TSSetDM(ts,da,ierr))
85: PetscCallA(TSSetType(ts,TSARKIMEX,ierr))
86: PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
87: PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr))
88: PetscCallA(DMSetMatType(da,MATAIJ,ierr))
89: PetscCallA(DMCreateMatrix(da,J,ierr))
90: PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr))
92: PetscCallA(TSGetSNES(ts,snes,ierr))
93: PetscCallA(SNESGetLineSearch(snes,linesearch,ierr))
94: PetscCallA(SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr))
96: ftime = 1.0
97: PetscCallA(TSSetMaxTime(ts,ftime,ierr))
98: PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Set initial conditions
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: PetscCallA(FormInitialSolution(ts,X,user,ierr))
104: PetscCallA(TSSetSolution(ts,X,ierr))
105: PetscCallA(VecGetSize(X,mx,ierr))
106: ! Advective CFL, I don't know why it needs so much safety factor.
107: dt = pone * max(user(user_a+1),user(user_a+2)) / mx
108: PetscCallA(TSSetTimeStep(ts,dt,ierr))
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Set runtime options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: PetscCallA(TSSetFromOptions(ts,ierr))
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Solve nonlinear system
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: PetscCallA(TSSolve(ts,X,ierr))
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Free work space.
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: PetscCallA(MatDestroy(J,ierr))
124: PetscCallA(VecDestroy(X,ierr))
125: PetscCallA(TSDestroy(ts,ierr))
126: PetscCallA(DMDestroy(da,ierr))
127: PetscCallA(PetscFinalize(ierr))
128: contains
130: ! Small helper to extract the layout, result uses 1-based indexing.
131: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
132: use petscdm
133: implicit none
135: DM da
136: PetscInt mx,xs,xe,gxs,gxe
137: PetscErrorCode ierr
138: PetscInt xm,gxm
139: 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))
140: PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
141: PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
142: xs = xs + 1
143: gxs = gxs + 1
144: xe = xs + xm - 1
145: gxe = gxs + gxm - 1
146: end subroutine
148: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
149: implicit none
150: PetscInt mx,xs,xe,gxs,gxe
151: PetscScalar x(2,xs:xe)
152: PetscScalar xdot(2,xs:xe)
153: PetscScalar f(2,xs:xe)
154: PetscReal a(2),k(2),s(2)
155: PetscErrorCode ierr
156: PetscInt i
157: do 10 i = xs,xe
158: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
159: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
160: 10 continue
161: end subroutine
163: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
164: use petscts
165: implicit none
167: TS ts
168: PetscReal t
169: Vec X,Xdot,F
170: PetscReal user(6)
171: PetscErrorCode ierr
172: integer user_a,user_k,user_s
173: parameter (user_a = 1,user_k = 3,user_s = 5)
175: DM da
176: PetscInt mx,xs,xe,gxs,gxe
177: PetscScalar,pointer :: xx(:),xxdot(:),ff(:)
179: PetscCall(TSGetDM(ts,da,ierr))
180: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
182: ! Get access to vector data
183: PetscCall(VecGetArrayRead(X,xx,ierr))
184: PetscCall(VecGetArrayRead(Xdot,xxdot,ierr))
185: PetscCall(VecGetArray(F,ff,ierr))
187: PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx,xxdot,ff,user(user_a),user(user_k),user(user_s),ierr))
189: PetscCall(VecRestoreArrayRead(X,xx,ierr))
190: PetscCall(VecRestoreArrayRead(Xdot,xxdot,ierr))
191: PetscCall(VecRestoreArray(F,ff,ierr))
192: end subroutine
194: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
195: implicit none
196: PetscInt mx,xs,xe,gxs,gxe
197: PetscReal t
198: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
199: PetscReal a(2),k(2),s(2)
200: PetscErrorCode ierr
201: PetscInt i,j
202: PetscReal hx,u0t(2)
203: PetscReal one,two,three,four,six,twelve
204: PetscReal half,third,twothird,sixth
205: PetscReal twelfth
207: one = 1.0
208: two = 2.0
209: three = 3.0
210: four = 4.0
211: six = 6.0
212: twelve = 12.0
213: hx = one / mx
214: ! The Fortran standard only allows positive base for power functions; Nag compiler fails on this
215: u0t(1) = one - abs(sin(twelve*t))**four
216: u0t(2) = 0.0
217: half = one/two
218: third = one / three
219: twothird = two / three
220: sixth = one / six
221: twelfth = one / twelve
222: do 20 i = xs,xe
223: do 10 j = 1,2
224: if (i .eq. 1) then
225: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) &
226: & + sixth*x(j,i+2))
227: else if (i .eq. 2) then
228: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) &
229: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
230: else if (i .eq. mx-1) then
231: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) &
232: & - half*x(j,i) -third*x(j,i+1))
233: else if (i .eq. mx) then
234: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
235: else
236: f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) &
237: & + twothird*x(j,i-1) &
238: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
239: end if
240: 10 continue
241: 20 continue
242: end subroutine
244: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
245: use petscts
246: implicit none
248: TS ts
249: PetscReal t
250: Vec X,F
251: PetscReal user(6)
252: PetscErrorCode ierr
253: integer user_a,user_k,user_s
254: parameter (user_a = 1,user_k = 3,user_s = 5)
255: DM da
256: Vec Xloc
257: PetscInt mx,xs,xe,gxs,gxe
258: PetscScalar,pointer :: xx(:),ff(:)
260: PetscCall(TSGetDM(ts,da,ierr))
261: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
263: ! Scatter ghost points to local vector,using the 2-step process
264: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
265: ! By placing code between these two statements, computations can be
266: ! done while messages are in transition.
267: PetscCall(DMGetLocalVector(da,Xloc,ierr))
268: PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
269: PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))
271: ! Get access to vector data
272: PetscCall(VecGetArrayRead(Xloc,xx,ierr))
273: PetscCall(VecGetArray(F,ff,ierr))
275: PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx,ff,user(user_a),user(user_k),user(user_s),ierr))
277: PetscCall(VecRestoreArrayRead(Xloc,xx,ierr))
278: PetscCall(VecRestoreArray(F,ff,ierr))
279: PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
280: end subroutine
282: ! ---------------------------------------------------------------------
283: !
284: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
285: !
286: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
287: use petscts
288: implicit none
290: TS ts
291: PetscReal t,shift
292: Vec X,Xdot
293: Mat J,Jpre
294: PetscReal user(6)
295: PetscErrorCode ierr
296: integer user_a,user_k,user_s
297: parameter (user_a = 0,user_k = 2,user_s = 4)
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 = user(user_k+1)
310: k2 = user(user_k+2)
311: do 10 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: 10 continue
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 10 i=xs,xe
340: r = i*hx
341: if (k(2) .ne. 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: 10 continue
349: end subroutine
351: subroutine FormInitialSolution(ts,X,user,ierr)
352: use petscts
353: implicit none
355: TS ts
356: Vec X
357: PetscReal user(6)
358: PetscErrorCode ierr
359: integer user_a,user_k,user_s
360: parameter (user_a = 1,user_k = 3,user_s = 5)
362: DM da
363: PetscInt mx,xs,xe,gxs,gxe
364: PetscScalar,pointer :: xx(:)
366: PetscCall(TSGetDM(ts,da,ierr))
367: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
369: ! Get access to vector data
370: PetscCall(VecGetArray(X,xx,ierr))
372: PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx,user(user_a),user(user_k),user(user_s),ierr))
374: PetscCall(VecRestoreArray(X,xx,ierr))
375: end subroutine
376: end program
377: !/*TEST
378: !
379: ! test:
380: ! args: -da_grid_x 200 -ts_arkimex_type 4
381: ! requires: !single
382: ! output_file: output/empty.out
383: !
384: !TEST*/