Actual source code: ex22f.F
petsc-3.4.5 2014-06-29
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: SNES snes
40: SNESLineSearch linesearch
41: Vec X
42: Mat J
43: PetscInt steps,maxsteps,mx
44: PetscErrorCode ierr
45: DM da
46: PetscReal ftime,dt
47: PetscReal zero,one,pone
48: PetscInt im11,i2
49: PetscBool flg
51: im11 = -11
52: i2 = 2
53: zero = 0.0
54: one = 1.0
55: pone = one / 10
57: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: ! Create distributed array (DMDA) to manage parallel grid and vectors
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,im11,i2,i2, &
63: & PETSC_NULL_INTEGER,da,ierr)
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: ! Extract global vectors from DMDA;
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: call DMCreateGlobalVector(da,X,ierr)
70: ! Initialize user application context
71: ! Use zero-based indexing for command line parameters to match ex22.c
72: user(user_a+1) = 1.0
73: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-a0', &
74: & user(user_a+1),flg,ierr)
75: user(user_a+2) = 0.0
76: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-a1', &
77: & user(user_a+2),flg,ierr)
78: user(user_k+1) = 1000000.0
79: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-k0', &
80: & user(user_k+1),flg,ierr)
81: user(user_k+2) = 2*user(user_k+1)
82: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-k1', &
83: & user(user_k+2),flg,ierr)
84: user(user_s+1) = 0.0
85: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-s0', &
86: & user(user_s+1),flg,ierr)
87: user(user_s+2) = 1.0
88: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-s1', &
89: & user(user_s+2),flg,ierr)
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Create timestepping solver context
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
95: call TSSetDM(ts,da,ierr)
96: call TSSetType(ts,TSARKIMEX,ierr)
97: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction, &
98: & user,ierr)
99: call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
100: call DMCreateMatrix(da,MATAIJ,J,ierr)
101: call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)
103: call TSGetSNES(ts,snes,ierr)
104: call SNESGetLineSearch(snes,linesearch,ierr)
105: call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)
107: ftime = 1.0
108: maxsteps = 10000
109: call TSSetDuration(ts,maxsteps,ftime,ierr)
111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: ! Set initial conditions
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: call FormInitialSolution(ts,X,user,ierr)
115: call TSSetSolution(ts,X,ierr)
116: call VecGetSize(X,mx,ierr)
117: ! Advective CFL, I don't know why it needs so much safety factor.
118: dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
119: call TSSetInitialTimeStep(ts,zero,dt,ierr)
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! Set runtime options
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: call TSSetFromOptions(ts,ierr)
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! Solve nonlinear system
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: call TSSolve(ts,X,ierr)
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: ! Free work space.
133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: call MatDestroy(J,ierr)
135: call VecDestroy(X,ierr)
136: call TSDestroy(ts,ierr)
137: call DMDestroy(da,ierr)
138: call PetscFinalize(ierr)
139: end program
141: ! Small helper to extract the layout, result uses 1-based indexing.
142: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
143: implicit none
144: #include <finclude/petscsys.h>
145: #include <finclude/petscdmda.h>
146: DM da
147: PetscInt mx,xs,xe,gxs,gxe
148: PetscErrorCode ierr
149: PetscInt xm,gxm
150: call DMDAGetInfo(da,PETSC_NULL_INTEGER, &
151: & mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
152: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
153: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
154: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
155: & PETSC_NULL_INTEGER,ierr)
156: call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
157: & xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
158: call DMDAGetGhostCorners(da, &
159: & gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
160: & gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
161: xs = xs + 1
162: gxs = gxs + 1
163: xe = xs + xm - 1
164: gxe = gxs + gxm - 1
165: end subroutine
167: subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f, &
168: & a,k,s,ierr)
169: implicit none
170: PetscInt mx,xs,xe,gxs,gxe
171: PetscScalar x(2,xs:xe)
172: PetscScalar xdot(2,xs:xe)
173: PetscScalar f(2,xs:xe)
174: PetscReal a(2),k(2),s(2)
175: PetscErrorCode ierr
176: PetscInt i
177: do 10 i = xs,xe
178: f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
179: f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
180: 10 continue
181: end subroutine
183: subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
184: implicit none
185: #include <finclude/petscsys.h>
186: #include <finclude/petscvec.h>
187: #include <finclude/petscmat.h>
188: #include <finclude/petscsnes.h>
189: #include <finclude/petscts.h>
190: #include <finclude/petscdmda.h>
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 VecGetArray(X,xx,ixx,ierr)
209: call VecGetArray(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 VecRestoreArray(X,xx,ixx,ierr)
217: call VecRestoreArray(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: u0t(1) = one - sin(twelve*t)**four
243: u0t(2) = 0.0
244: half = one/two
245: third = one / three
246: twothird = two / three
247: sixth = one / six
248: twelfth = one / twelve
249: do 20 i = xs,xe
250: do 10 j = 1,2
251: if (i .eq. 1) then
252: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) &
253: & + sixth*x(j,i+2))
254: else if (i .eq. 2) then
255: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) &
256: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
257: else if (i .eq. mx-1) then
258: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) &
259: & - half*x(j,i) -third*x(j,i+1))
260: else if (i .eq. mx) then
261: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
262: else
263: f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) &
264: & + twothird*x(j,i-1) &
265: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
266: end if
267: 10 continue
268: 20 continue
269: end subroutine
271: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
272: implicit none
273: #include <finclude/petscsys.h>
274: #include <finclude/petscvec.h>
275: #include <finclude/petscmat.h>
276: #include <finclude/petscsnes.h>
277: #include <finclude/petscts.h>
278: #include <finclude/petscdmda.h>
279: TS ts
280: PetscReal t
281: Vec X,F
282: PetscReal user(6)
283: PetscErrorCode ierr
284: integer user_a,user_k,user_s
285: parameter (user_a = 1,user_k = 3,user_s = 5)
286: DM da
287: Vec Xloc
288: PetscInt mx,xs,xe,gxs,gxe
289: PetscOffset ixx,iff
290: PetscScalar xx(0:1),ff(0:1)
292: call TSGetDM(ts,da,ierr)
293: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
295: ! Scatter ghost points to local vector,using the 2-step process
296: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
297: ! By placing code between these two statements, computations can be
298: ! done while messages are in transition.
299: call DMGetLocalVector(da,Xloc,ierr)
300: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
301: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)
303: ! Get access to vector data
304: call VecGetArray(Xloc,xx,ixx,ierr)
305: call VecGetArray(F,ff,iff,ierr)
307: call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe, &
308: & t,xx(ixx),ff(iff), &
309: & user(user_a),user(user_k),user(user_s),ierr)
311: call VecRestoreArray(Xloc,xx,ixx,ierr)
312: call VecRestoreArray(F,ff,iff,ierr)
313: call DMRestoreLocalVector(da,Xloc,ierr)
314: end subroutine
316: ! ---------------------------------------------------------------------
317: !
318: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
319: !
320: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,mstr,user,ierr)
321: implicit none
322: #include <finclude/petscsys.h>
323: #include <finclude/petscvec.h>
324: #include <finclude/petscmat.h>
325: #include <finclude/petscsnes.h>
326: #include <finclude/petscts.h>
327: #include <finclude/petscdmda.h>
328: TS ts
329: PetscReal t,shift
330: Vec X,Xdot
331: Mat J,Jpre
332: MatStructure mstr
333: PetscReal user(6)
334: PetscErrorCode ierr
335: integer user_a,user_k,user_s
336: parameter (user_a = 0,user_k = 2,user_s = 4)
338: DM da
339: PetscInt mx,xs,xe,gxs,gxe
340: PetscInt i,i1,row,col
341: PetscReal k1,k2;
342: PetscScalar val(4)
344: call TSGetDM(ts,da,ierr)
345: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
347: i1 = 1
348: k1 = user(user_k+1)
349: k2 = user(user_k+2)
350: do 10 i = xs,xe
351: row = i-gxs
352: col = i-gxs
353: val(1) = shift + k1
354: val(2) = -k2
355: val(3) = -k1
356: val(4) = shift + k2
357: call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val, &
358: & INSERT_VALUES,ierr)
359: 10 continue
360: call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
361: call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
362: if (J /= Jpre) then
363: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
364: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
365: end if
366: mstr = SAME_NONZERO_PATTERN
367: end subroutine
369: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
370: implicit none
371: PetscInt mx,xs,xe,gxs,gxe
372: PetscScalar x(2,xs:xe)
373: PetscReal a(2),k(2),s(2)
374: PetscErrorCode ierr
376: PetscInt i
377: PetscReal one,hx,r,ik
378: one = 1.0
379: hx = one / mx
380: do 10 i=xs,xe
381: r = i*hx
382: if (k(2) .ne. 0.0) then
383: ik = one/k(2)
384: else
385: ik = one
386: end if
387: x(1,i) = one + s(2)*r
388: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
389: 10 continue
390: end subroutine
392: subroutine FormInitialSolution(ts,X,user,ierr)
393: implicit none
394: #include <finclude/petscsys.h>
395: #include <finclude/petscvec.h>
396: #include <finclude/petscmat.h>
397: #include <finclude/petscsnes.h>
398: #include <finclude/petscts.h>
399: #include <finclude/petscdmda.h>
400: TS ts
401: Vec X
402: PetscReal user(6)
403: PetscErrorCode ierr
404: integer user_a,user_k,user_s
405: parameter (user_a = 1,user_k = 3,user_s = 5)
407: DM da
408: PetscInt mx,xs,xe,gxs,gxe
409: PetscOffset ixx
410: PetscScalar xx(0:1)
412: call TSGetDM(ts,da,ierr)
413: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
415: ! Get access to vector data
416: call VecGetArray(X,xx,ixx,ierr)
418: call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe, &
419: & xx(ixx),user(user_a),user(user_k),user(user_s),ierr)
421: call VecRestoreArray(X,xx,ixx,ierr)
422: end subroutine