Actual source code: ex22f.F
petsc-3.3-p7 2013-05-11
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: implicit none
19: #include <finclude/petscsys.h>
20: #include <finclude/petscvec.h>
21: #include <finclude/petscmat.h>
22: #include <finclude/petscsnes.h>
23: #include <finclude/petscts.h>
24: #include <finclude/petscdmda.h>
25: !
26: ! Create an application context to contain data needed by the
27: ! application-provided call-back routines, FormJacobian() and
28: ! FormFunction(). We use a double precision array with six
29: ! entries, two for each problem parameter a, k, s.
30: !
31: PetscReal user(6)
32: integer user_a,user_k,user_s
33: parameter (user_a = 0,user_k = 2,user_s = 4)
35: external FormRHSFunction,FormIFunction,FormIJacobian
36: external FormInitialSolution
38: TS ts
39: Vec X
40: Mat J
41: PetscInt steps,maxsteps,mx
42: PetscErrorCode ierr
43: DM da
44: PetscReal ftime,dt
45: PetscReal zero,one,pone
46: PetscInt im11,i2
47: PetscBool flg
49: im11 = -11
50: i2 = 2
51: zero = 0.0
52: one = 1.0
53: pone = one / 10
55: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: ! Create distributed array (DMDA) to manage parallel grid and vectors
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,im11,i2,i2, &
61: & PETSC_NULL_INTEGER,da,ierr)
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Extract global vectors from DMDA;
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: call 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: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-a0', &
72: & user(user_a+1),flg,ierr)
73: user(user_a+2) = 0.0
74: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-a1', &
75: & user(user_a+2),flg,ierr)
76: user(user_k+1) = 1000000.0
77: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-k0', &
78: & user(user_k+1),flg,ierr)
79: user(user_k+2) = 2*user(user_k+1)
80: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-k1', &
81: & user(user_k+2),flg,ierr)
82: user(user_s+1) = 0.0
83: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-s0', &
84: & user(user_s+1),flg,ierr)
85: user(user_s+2) = 1.0
86: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-s1', &
87: & user(user_s+2),flg,ierr)
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: ! Create timestepping solver context
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
93: call TSSetDM(ts,da,ierr)
94: call TSSetType(ts,TSARKIMEX,ierr)
95: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction, &
96: & user,ierr)
97: call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
98: call DMCreateMatrix(da,MATAIJ,J,ierr)
99: call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)
101: ftime = 1.0
102: maxsteps = 10000
103: call TSSetDuration(ts,maxsteps,ftime,ierr)
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: ! Set initial conditions
107: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: call FormInitialSolution(ts,X,user,ierr)
109: call TSSetSolution(ts,X,ierr)
110: call VecGetSize(X,mx,ierr)
111: ! Advective CFL, I don't know why it needs so much safety factor.
112: dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
113: call TSSetInitialTimeStep(ts,zero,dt,ierr)
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Set runtime options
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: call TSSetFromOptions(ts,ierr)
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Solve nonlinear system
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: call TSSolve(ts,X,ftime,ierr)
124: call TSGetTimeStepNumber(ts,steps,ierr)
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! Free work space.
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: call MatDestroy(J,ierr)
130: call VecDestroy(X,ierr)
131: call TSDestroy(ts,ierr)
132: call DMDestroy(da,ierr)
133: call PetscFinalize(ierr)
134: end program
136: ! Small helper to extract the layout, result uses 1-based indexing.
137: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
138: implicit none
139: #include <finclude/petscsys.h>
140: #include <finclude/petscdmda.h>
141: DM da
142: PetscInt mx,xs,xe,gxs,gxe
143: PetscErrorCode ierr
144: PetscInt xm,gxm
145: call DMDAGetInfo(da,PETSC_NULL_INTEGER, &
146: & mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
147: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
148: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
149: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
150: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
151: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
152: & xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
153: call DMDAGetGhostCorners(da, &
154: & gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
155: & gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
156: xs = xs + 1
157: gxs = gxs + 1
158: xe = xs + xm - 1
159: gxe = gxs + gxm - 1
160: end subroutine
162: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f, &
163: & a,k,s,ierr)
164: implicit none
165: PetscInt mx,xs,xe,gxs,gxe
166: PetscScalar x(2,xs:xe)
167: PetscScalar xdot(2,xs:xe)
168: PetscScalar f(2,xs:xe)
169: PetscReal a(2),k(2),s(2)
170: PetscErrorCode ierr
171: PetscInt i
172: do 10 i = xs,xe
173: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
174: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
175: 10 continue
176: end subroutine
178: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
179: implicit none
180: #include <finclude/petscsys.h>
181: #include <finclude/petscvec.h>
182: #include <finclude/petscmat.h>
183: #include <finclude/petscsnes.h>
184: #include <finclude/petscts.h>
185: #include <finclude/petscdmda.h>
186: TS ts
187: PetscReal t
188: Vec X,Xdot,F
189: PetscReal user(6)
190: PetscErrorCode ierr
191: integer user_a,user_k,user_s
192: parameter (user_a = 1,user_k = 3,user_s = 5)
194: DM da
195: PetscInt mx,xs,xe,gxs,gxe
196: PetscOffset ixx,ixxdot,iff
197: PetscScalar xx(0:1),xxdot(0:1),ff(0:1)
199: call TSGetDM(ts,da,ierr)
200: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
202: ! Get access to vector data
203: call VecGetArray(X,xx,ixx,ierr)
204: call VecGetArray(Xdot,xxdot,ixxdot,ierr)
205: call VecGetArray(F,ff,iff,ierr)
207: call FormIFunctionLocal(mx,xs,xe,gxs,gxe, &
208: & xx(ixx),xxdot(ixxdot),ff(iff), &
209: & user(user_a),user(user_k),user(user_s),ierr)
211: call VecRestoreArray(X,xx,ixx,ierr)
212: call VecRestoreArray(Xdot,xxdot,ixxdot,ierr)
213: call VecRestoreArray(F,ff,iff,ierr)
214: end subroutine
216: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f, &
217: & a,k,s,ierr)
218: implicit none
219: PetscInt mx,xs,xe,gxs,gxe
220: PetscReal t
221: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
222: PetscReal a(2),k(2),s(2)
223: PetscErrorCode ierr
224: PetscInt i,j
225: PetscReal hx,u0t(2)
226: PetscReal one,two,three,four,six,twelve
227: PetscReal half,third,twothird,sixth
228: PetscReal twelfth
230: one = 1.0
231: two = 2.0
232: three = 3.0
233: four = 4.0
234: six = 6.0
235: twelve = 12.0
236: hx = one / mx
237: u0t(1) = one - sin(twelve*t)**four
238: u0t(2) = 0.0
239: half = one/two
240: third = one / three
241: twothird = two / three
242: sixth = one / six
243: twelfth = one / twelve
244: do 20 i = xs,xe
245: do 10 j = 1,2
246: if (i .eq. 1) then
247: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) &
248: & + sixth*x(j,i+2))
249: else if (i .eq. 2) then
250: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) &
251: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
252: else if (i .eq. mx-1) then
253: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) &
254: & - half*x(j,i) -third*x(j,i+1))
255: else if (i .eq. mx) then
256: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
257: else
258: f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) &
259: & + twothird*x(j,i-1) &
260: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
261: end if
262: 10 continue
263: 20 continue
264: end subroutine
266: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
267: implicit none
268: #include <finclude/petscsys.h>
269: #include <finclude/petscvec.h>
270: #include <finclude/petscmat.h>
271: #include <finclude/petscsnes.h>
272: #include <finclude/petscts.h>
273: #include <finclude/petscdmda.h>
274: TS ts
275: PetscReal t
276: Vec X,F
277: PetscReal user(6)
278: PetscErrorCode ierr
279: integer user_a,user_k,user_s
280: parameter (user_a = 1,user_k = 3,user_s = 5)
281: DM da
282: Vec Xloc
283: PetscInt mx,xs,xe,gxs,gxe
284: PetscOffset ixx,iff
285: PetscScalar xx(0:1),ff(0:1)
287: call TSGetDM(ts,da,ierr)
288: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
290: ! Scatter ghost points to local vector,using the 2-step process
291: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
292: ! By placing code between these two statements, computations can be
293: ! done while messages are in transition.
294: call DMGetLocalVector(da,Xloc,ierr)
295: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
296: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)
298: ! Get access to vector data
299: call VecGetArray(Xloc,xx,ixx,ierr)
300: call VecGetArray(F,ff,iff,ierr)
302: call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe, &
303: & t,xx(ixx),ff(iff), &
304: & user(user_a),user(user_k),user(user_s),ierr)
306: call VecRestoreArray(Xloc,xx,ixx,ierr)
307: call VecRestoreArray(F,ff,iff,ierr)
308: call DMRestoreLocalVector(da,Xloc,ierr)
309: end subroutine
311: ! ---------------------------------------------------------------------
312: !
313: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
314: !
315: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,mstr,user,ierr)
316: implicit none
317: #include <finclude/petscsys.h>
318: #include <finclude/petscvec.h>
319: #include <finclude/petscmat.h>
320: #include <finclude/petscsnes.h>
321: #include <finclude/petscts.h>
322: #include <finclude/petscdmda.h>
323: TS ts
324: PetscReal t,shift
325: Vec X,Xdot
326: Mat J,Jpre
327: MatStructure mstr
328: PetscReal user(6)
329: PetscErrorCode ierr
330: integer user_a,user_k,user_s
331: parameter (user_a = 0,user_k = 2,user_s = 4)
333: DM da
334: PetscInt mx,xs,xe,gxs,gxe
335: PetscInt i,i1,row,col
336: PetscReal k1,k2;
337: PetscScalar val(4)
339: call TSGetDM(ts,da,ierr)
340: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
342: i1 = 1
343: k1 = user(user_k+1)
344: k2 = user(user_k+2)
345: do 10 i = xs,xe
346: row = i-gxs
347: col = i-gxs
348: val(1) = shift + k1
349: val(2) = -k2
350: val(3) = -k1
351: val(4) = shift + k2
352: call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val, &
353: & INSERT_VALUES,ierr)
354: 10 continue
355: call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
356: call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
357: if (J /= Jpre) then
358: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
359: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
360: end if
361: mstr = SAME_NONZERO_PATTERN
362: end subroutine
364: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
365: implicit none
366: PetscInt mx,xs,xe,gxs,gxe
367: PetscScalar x(2,xs:xe)
368: PetscReal a(2),k(2),s(2)
369: PetscErrorCode ierr
371: PetscInt i
372: PetscReal one,hx,r,ik
373: one = 1.0
374: hx = one / mx
375: do 10 i=xs,xe
376: r = i*hx
377: if (k(2) .ne. 0.0) then
378: ik = one/k(2)
379: else
380: ik = one
381: end if
382: x(1,i) = one + s(2)*r
383: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
384: 10 continue
385: end subroutine
387: subroutine FormInitialSolution(ts,X,user,ierr)
388: implicit none
389: #include <finclude/petscsys.h>
390: #include <finclude/petscvec.h>
391: #include <finclude/petscmat.h>
392: #include <finclude/petscsnes.h>
393: #include <finclude/petscts.h>
394: #include <finclude/petscdmda.h>
395: TS ts
396: Vec X
397: PetscReal user(6)
398: PetscErrorCode ierr
399: integer user_a,user_k,user_s
400: parameter (user_a = 1,user_k = 3,user_s = 5)
402: DM da
403: PetscInt mx,xs,xe,gxs,gxe
404: PetscOffset ixx
405: PetscScalar xx(0:1)
407: call TSGetDM(ts,da,ierr)
408: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
410: ! Get access to vector data
411: call VecGetArray(X,xx,ixx,ierr)
413: call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe, &
414: & xx(ixx),user(user_a),user(user_k),user(user_s),ierr)
416: call VecRestoreArray(X,xx,ixx,ierr)
417: end subroutine