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: 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: PetscCallA(PetscInitialize(ierr))
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: ! Create distributed array (DMDA) to manage parallel grid and vectors
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr))
60: PetscCallA(DMSetFromOptions(da,ierr))
61: PetscCallA(DMSetUp(da,ierr))
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Extract global vectors from DMDA;
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: PetscCallA(DMCreateGlobalVector(da,X,ierr))
68: ! Initialize user application context
69: ! Use zero-based indexing for command line parameters to match ex22.c
70: user(user_a+1) = 1.0
71: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr))
72: user(user_a+2) = 0.0
73: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr))
74: user(user_k+1) = 1000000.0
75: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0',user(user_k+1),flg,ierr))
76: user(user_k+2) = 2*user(user_k+1)
77: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr))
78: user(user_s+1) = 0.0
79: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr))
80: user(user_s+2) = 1.0
81: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr))
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Create timestepping solver context
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
87: PetscCallA(TSSetDM(ts,da,ierr))
88: PetscCallA(TSSetType(ts,TSARKIMEX,ierr))
89: PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
90: PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr))
91: PetscCallA(DMSetMatType(da,MATAIJ,ierr))
92: PetscCallA(DMCreateMatrix(da,J,ierr))
93: PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr))
95: PetscCallA(TSGetSNES(ts,snes,ierr))
96: PetscCallA(SNESGetLineSearch(snes,linesearch,ierr))
97: PetscCallA(SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr))
99: ftime = 1.0
100: PetscCallA(TSSetMaxTime(ts,ftime,ierr))
101: PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! Set initial conditions
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: PetscCallA(FormInitialSolution(ts,X,user,ierr))
107: PetscCallA(TSSetSolution(ts,X,ierr))
108: PetscCallA(VecGetSize(X,mx,ierr))
109: ! Advective CFL, I don't know why it needs so much safety factor.
110: dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
111: PetscCallA(TSSetTimeStep(ts,dt,ierr))
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! Set runtime options
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: PetscCallA(TSSetFromOptions(ts,ierr))
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: ! Solve nonlinear system
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: PetscCallA(TSSolve(ts,X,ierr))
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! Free work space.
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: PetscCallA(MatDestroy(J,ierr))
127: PetscCallA(VecDestroy(X,ierr))
128: PetscCallA(TSDestroy(ts,ierr))
129: PetscCallA(DMDestroy(da,ierr))
130: PetscCallA(PetscFinalize(ierr))
131: end program
133: ! Small helper to extract the layout, result uses 1-based indexing.
134: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
135: use petscdmda
136: implicit none
138: DM da
139: PetscInt mx,xs,xe,gxs,gxe
140: PetscErrorCode ierr
141: PetscInt xm,gxm
142: 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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
143: PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
144: PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
145: xs = xs + 1
146: gxs = gxs + 1
147: xe = xs + xm - 1
148: gxe = gxs + gxm - 1
149: end subroutine
151: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
152: implicit none
153: PetscInt mx,xs,xe,gxs,gxe
154: PetscScalar x(2,xs:xe)
155: PetscScalar xdot(2,xs:xe)
156: PetscScalar f(2,xs:xe)
157: PetscReal a(2),k(2),s(2)
158: PetscErrorCode ierr
159: PetscInt i
160: do 10 i = xs,xe
161: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
162: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
163: 10 continue
164: end subroutine
166: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
167: use petscts
168: implicit none
170: TS ts
171: PetscReal t
172: Vec X,Xdot,F
173: PetscReal user(6)
174: PetscErrorCode ierr
175: integer user_a,user_k,user_s
176: parameter (user_a = 1,user_k = 3,user_s = 5)
178: DM da
179: PetscInt mx,xs,xe,gxs,gxe
180: PetscOffset ixx,ixxdot,iff
181: PetscScalar xx(0:1),xxdot(0:1),ff(0:1)
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,ixx,ierr))
188: PetscCall(VecGetArrayRead(Xdot,xxdot,ixxdot,ierr))
189: PetscCall(VecGetArray(F,ff,iff,ierr))
191: PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr))
193: PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr))
194: PetscCall(VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr))
195: PetscCall(VecRestoreArray(F,ff,iff,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 20 i = xs,xe
227: do 10 j = 1,2
228: if (i .eq. 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 .eq. 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 .eq. 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 .eq. 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: 10 continue
245: 20 continue
246: end subroutine
248: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
249: use petscts
250: implicit none
252: TS ts
253: PetscReal t
254: Vec X,F
255: PetscReal user(6)
256: PetscErrorCode ierr
257: integer user_a,user_k,user_s
258: parameter (user_a = 1,user_k = 3,user_s = 5)
259: DM da
260: Vec Xloc
261: PetscInt mx,xs,xe,gxs,gxe
262: PetscOffset ixx,iff
263: PetscScalar xx(0:1),ff(0:1)
265: PetscCall(TSGetDM(ts,da,ierr))
266: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
268: ! Scatter ghost points to local vector,using the 2-step process
269: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
270: ! By placing code between these two statements, computations can be
271: ! done while messages are in transition.
272: PetscCall(DMGetLocalVector(da,Xloc,ierr))
273: PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
274: PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))
276: ! Get access to vector data
277: PetscCall(VecGetArrayRead(Xloc,xx,ixx,ierr))
278: PetscCall(VecGetArray(F,ff,iff,ierr))
280: PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx(ixx),ff(iff),user(user_a),user(user_k),user(user_s),ierr))
282: PetscCall(VecRestoreArrayRead(Xloc,xx,ixx,ierr))
283: PetscCall(VecRestoreArray(F,ff,iff,ierr))
284: PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
285: end subroutine
287: ! ---------------------------------------------------------------------
288: !
289: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
290: !
291: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
292: use petscts
293: implicit none
295: TS ts
296: PetscReal t,shift
297: Vec X,Xdot
298: Mat J,Jpre
299: PetscReal user(6)
300: PetscErrorCode ierr
301: integer user_a,user_k,user_s
302: parameter (user_a = 0,user_k = 2,user_s = 4)
304: DM da
305: PetscInt mx,xs,xe,gxs,gxe
306: PetscInt i,i1,row,col
307: PetscReal k1,k2;
308: PetscScalar val(4)
310: PetscCall(TSGetDM(ts,da,ierr))
311: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
313: i1 = 1
314: k1 = user(user_k+1)
315: k2 = user(user_k+2)
316: do 10 i = xs,xe
317: row = i-gxs
318: col = i-gxs
319: val(1) = shift + k1
320: val(2) = -k2
321: val(3) = -k1
322: val(4) = shift + k2
323: PetscCall(MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,INSERT_VALUES,ierr))
324: 10 continue
325: PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
326: PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
327: if (J /= Jpre) then
328: PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
329: PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
330: end if
331: end subroutine
333: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
334: implicit none
335: PetscInt mx,xs,xe,gxs,gxe
336: PetscScalar x(2,xs:xe)
337: PetscReal a(2),k(2),s(2)
338: PetscErrorCode ierr
340: PetscInt i
341: PetscReal one,hx,r,ik
342: one = 1.0
343: hx = one / mx
344: do 10 i=xs,xe
345: r = i*hx
346: if (k(2) .ne. 0.0) then
347: ik = one/k(2)
348: else
349: ik = one
350: end if
351: x(1,i) = one + s(2)*r
352: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
353: 10 continue
354: end subroutine
356: subroutine FormInitialSolution(ts,X,user,ierr)
357: use petscts
358: implicit none
360: TS ts
361: Vec X
362: PetscReal user(6)
363: PetscErrorCode ierr
364: integer user_a,user_k,user_s
365: parameter (user_a = 1,user_k = 3,user_s = 5)
367: DM da
368: PetscInt mx,xs,xe,gxs,gxe
369: PetscOffset ixx
370: PetscScalar xx(0:1)
372: PetscCall(TSGetDM(ts,da,ierr))
373: PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))
375: ! Get access to vector data
376: PetscCall(VecGetArray(X,xx,ixx,ierr))
378: PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx(ixx),user(user_a),user(user_k),user(user_s),ierr))
380: PetscCall(VecRestoreArray(X,xx,ixx,ierr))
381: end subroutine
383: !/*TEST
384: !
385: ! test:
386: ! args: -da_grid_x 200 -ts_arkimex_type 4
387: ! requires: !single
388: !
389: !TEST*/