Actual source code: ex22f.F
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: external FormRHSFunction,FormIFunction,FormIJacobian
34: external FormInitialSolution
36: TS ts
37: SNES snes
38: SNESLineSearch linesearch
39: Vec X
40: Mat J
41: PetscInt mx
42: PetscErrorCode ierr
43: DM da
44: PetscReal ftime,dt
45: PetscReal one,pone
46: PetscInt im11,i2
47: PetscBool flg
49: im11 = 11
50: i2 = 2
51: one = 1.0
52: pone = one / 10
54: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
55: if (ierr .ne. 0) then
56: print*,'Unable to initialize PETSc'
57: stop
58: endif
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Create distributed array (DMDA) to manage parallel grid and vectors
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2, &
64: & PETSC_NULL_INTEGER,da,ierr)
65: call DMSetFromOptions(da,ierr)
66: call DMSetUp(da,ierr)
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: ! Extract global vectors from DMDA;
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: call DMCreateGlobalVector(da,X,ierr)
73: ! Initialize user application context
74: ! Use zero-based indexing for command line parameters to match ex22.c
75: user(user_a+1) = 1.0
76: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
77: & '-a0',user(user_a+1),flg,ierr)
78: user(user_a+2) = 0.0
79: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
80: & '-a1',user(user_a+2),flg,ierr)
81: user(user_k+1) = 1000000.0
82: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
83: & '-k0',user(user_k+1),flg,ierr)
84: user(user_k+2) = 2*user(user_k+1)
85: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
86: & '-k1', user(user_k+2),flg,ierr)
87: user(user_s+1) = 0.0
88: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
89: & '-s0',user(user_s+1),flg,ierr)
90: user(user_s+2) = 1.0
91: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
92: & '-s1',user(user_s+2),flg,ierr)
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Create timestepping solver context
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
98: call TSSetDM(ts,da,ierr)
99: call TSSetType(ts,TSARKIMEX,ierr)
100: call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction, &
101: & user,ierr)
102: call TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr)
103: call DMSetMatType(da,MATAIJ,ierr)
104: call DMCreateMatrix(da,J,ierr)
105: call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)
107: call TSGetSNES(ts,snes,ierr)
108: call SNESGetLineSearch(snes,linesearch,ierr)
109: call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)
111: ftime = 1.0
112: call TSSetMaxTime(ts,ftime,ierr)
113: call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Set initial conditions
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: call FormInitialSolution(ts,X,user,ierr)
119: call TSSetSolution(ts,X,ierr)
120: call VecGetSize(X,mx,ierr)
121: ! Advective CFL, I don't know why it needs so much safety factor.
122: dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
123: call TSSetTimeStep(ts,dt,ierr)
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Set runtime options
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: call TSSetFromOptions(ts,ierr)
130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: ! Solve nonlinear system
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: call TSSolve(ts,X,ierr)
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! Free work space.
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: call MatDestroy(J,ierr)
139: call VecDestroy(X,ierr)
140: call TSDestroy(ts,ierr)
141: call DMDestroy(da,ierr)
142: call PetscFinalize(ierr)
143: end program
145: ! Small helper to extract the layout, result uses 1-based indexing.
146: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
147: use petscdmda
148: implicit none
150: DM da
151: PetscInt mx,xs,xe,gxs,gxe
152: PetscErrorCode ierr
153: PetscInt xm,gxm
154: call DMDAGetInfo(da,PETSC_NULL_INTEGER, &
155: & mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
156: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
157: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
158: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
159: & PETSC_NULL_INTEGER,ierr)
160: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
161: & xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
162: call DMDAGetGhostCorners(da, &
163: & gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
164: & gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
165: xs = xs + 1
166: gxs = gxs + 1
167: xe = xs + xm - 1
168: gxe = gxs + gxm - 1
169: end subroutine
171: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f, &
172: & a,k,s,ierr)
173: implicit none
174: PetscInt mx,xs,xe,gxs,gxe
175: PetscScalar x(2,xs:xe)
176: PetscScalar xdot(2,xs:xe)
177: PetscScalar f(2,xs:xe)
178: PetscReal a(2),k(2),s(2)
179: PetscErrorCode ierr
180: PetscInt i
181: do 10 i = xs,xe
182: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
183: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
184: 10 continue
185: end subroutine
187: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
188: use petscts
189: implicit none
191: TS ts
192: PetscReal t
193: Vec X,Xdot,F
194: PetscReal user(6)
195: PetscErrorCode ierr
196: integer user_a,user_k,user_s
197: parameter (user_a = 1,user_k = 3,user_s = 5)
199: DM da
200: PetscInt mx,xs,xe,gxs,gxe
201: PetscOffset ixx,ixxdot,iff
202: PetscScalar xx(0:1),xxdot(0:1),ff(0:1)
204: call TSGetDM(ts,da,ierr)
205: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
207: ! Get access to vector data
208: call VecGetArrayRead(X,xx,ixx,ierr)
209: call VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)
210: call VecGetArray(F,ff,iff,ierr)
212: call FormIFunctionLocal(mx,xs,xe,gxs,gxe, &
213: & xx(ixx),xxdot(ixxdot),ff(iff), &
214: & user(user_a),user(user_k),user(user_s),ierr)
216: call VecRestoreArrayRead(X,xx,ixx,ierr)
217: call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)
218: call VecRestoreArray(F,ff,iff,ierr)
219: end subroutine
221: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f, &
222: & a,k,s,ierr)
223: implicit none
224: PetscInt mx,xs,xe,gxs,gxe
225: PetscReal t
226: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
227: PetscReal a(2),k(2),s(2)
228: PetscErrorCode ierr
229: PetscInt i,j
230: PetscReal hx,u0t(2)
231: PetscReal one,two,three,four,six,twelve
232: PetscReal half,third,twothird,sixth
233: PetscReal twelfth
235: one = 1.0
236: two = 2.0
237: three = 3.0
238: four = 4.0
239: six = 6.0
240: twelve = 12.0
241: hx = one / mx
242: ! The Fortran standard only allows positive base for power functions; Nag compiler fails on this
243: u0t(1) = one - abs(sin(twelve*t))**four
244: u0t(2) = 0.0
245: half = one/two
246: third = one / three
247: twothird = two / three
248: sixth = one / six
249: twelfth = one / twelve
250: do 20 i = xs,xe
251: do 10 j = 1,2
252: if (i .eq. 1) then
253: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) &
254: & + sixth*x(j,i+2))
255: else if (i .eq. 2) then
256: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) &
257: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
258: else if (i .eq. mx-1) then
259: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) &
260: & - half*x(j,i) -third*x(j,i+1))
261: else if (i .eq. mx) then
262: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
263: else
264: f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) &
265: & + twothird*x(j,i-1) &
266: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
267: end if
268: 10 continue
269: 20 continue
270: end subroutine
272: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
273: use petscts
274: implicit none
276: TS ts
277: PetscReal t
278: Vec X,F
279: PetscReal user(6)
280: PetscErrorCode ierr
281: integer user_a,user_k,user_s
282: parameter (user_a = 1,user_k = 3,user_s = 5)
283: DM da
284: Vec Xloc
285: PetscInt mx,xs,xe,gxs,gxe
286: PetscOffset ixx,iff
287: PetscScalar xx(0:1),ff(0:1)
289: call TSGetDM(ts,da,ierr)
290: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
292: ! Scatter ghost points to local vector,using the 2-step process
293: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
294: ! By placing code between these two statements, computations can be
295: ! done while messages are in transition.
296: call DMGetLocalVector(da,Xloc,ierr)
297: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
298: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)
300: ! Get access to vector data
301: call VecGetArrayRead(Xloc,xx,ixx,ierr)
302: call VecGetArray(F,ff,iff,ierr)
304: call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe, &
305: & t,xx(ixx),ff(iff), &
306: & user(user_a),user(user_k),user(user_s),ierr)
308: call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
309: call VecRestoreArray(F,ff,iff,ierr)
310: call DMRestoreLocalVector(da,Xloc,ierr)
311: end subroutine
313: ! ---------------------------------------------------------------------
314: !
315: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
316: !
317: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
318: use petscts
319: implicit none
321: TS ts
322: PetscReal t,shift
323: Vec X,Xdot
324: Mat J,Jpre
325: PetscReal user(6)
326: PetscErrorCode ierr
327: integer user_a,user_k,user_s
328: parameter (user_a = 0,user_k = 2,user_s = 4)
330: DM da
331: PetscInt mx,xs,xe,gxs,gxe
332: PetscInt i,i1,row,col
333: PetscReal k1,k2;
334: PetscScalar val(4)
336: call TSGetDM(ts,da,ierr)
337: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
339: i1 = 1
340: k1 = user(user_k+1)
341: k2 = user(user_k+2)
342: do 10 i = xs,xe
343: row = i-gxs
344: col = i-gxs
345: val(1) = shift + k1
346: val(2) = -k2
347: val(3) = -k1
348: val(4) = shift + k2
349: call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val, &
350: & INSERT_VALUES,ierr)
351: 10 continue
352: call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
353: call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
354: if (J /= Jpre) then
355: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
356: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
357: end if
358: end subroutine
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 10 i=xs,xe
372: r = i*hx
373: if (k(2) .ne. 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: 10 continue
381: end subroutine
383: subroutine FormInitialSolution(ts,X,user,ierr)
384: use petscts
385: implicit none
387: TS ts
388: Vec X
389: PetscReal user(6)
390: PetscErrorCode ierr
391: integer user_a,user_k,user_s
392: parameter (user_a = 1,user_k = 3,user_s = 5)
394: DM da
395: PetscInt mx,xs,xe,gxs,gxe
396: PetscOffset ixx
397: PetscScalar xx(0:1)
399: call TSGetDM(ts,da,ierr)
400: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
402: ! Get access to vector data
403: call VecGetArray(X,xx,ixx,ierr)
405: call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe, &
406: & xx(ixx),user(user_a),user(user_k),user(user_s),ierr)
408: call VecRestoreArray(X,xx,ixx,ierr)
409: end subroutine
411: !/*TEST
412: !
413: ! test:
414: ! args: -da_grid_x 200 -ts_arkimex_type 4
415: ! requires: !single
416: !
417: !TEST*/