Actual source code: ex22f.F
petsc-3.7.5 2017-01-01
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 <petsc/finclude/petscsys.h>
20: #include <petsc/finclude/petscvec.h>
21: #include <petsc/finclude/petscmat.h>
22: #include <petsc/finclude/petscsnes.h>
23: #include <petsc/finclude/petscts.h>
24: #include <petsc/finclude/petscdm.h>
25: #include <petsc/finclude/petscdmda.h>
26: !
27: ! Create an application context to contain data needed by the
28: ! application-provided call-back routines, FormJacobian() and
29: ! FormFunction(). We use a double precision array with six
30: ! entries, two for each problem parameter a, k, s.
31: !
32: PetscReal user(6)
33: integer user_a,user_k,user_s
34: parameter (user_a = 0,user_k = 2,user_s = 4)
36: external FormRHSFunction,FormIFunction,FormIJacobian
37: external FormInitialSolution
39: TS ts
40: SNES snes
41: SNESLineSearch linesearch
42: Vec X
43: Mat J
44: PetscInt maxsteps,mx
45: PetscErrorCode ierr
46: DM da
47: PetscReal ftime,dt
48: PetscReal zero,one,pone
49: PetscInt im11,i2
50: PetscBool flg
52: im11 = -11
53: i2 = 2
54: zero = 0.0
55: one = 1.0
56: pone = one / 10
58: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
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)
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: ! Extract global vectors from DMDA;
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: call DMCreateGlobalVector(da,X,ierr)
71: ! Initialize user application context
72: ! Use zero-based indexing for command line parameters to match ex22.c
73: user(user_a+1) = 1.0
74: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
75: & '-a0',user(user_a+1),flg,ierr)
76: user(user_a+2) = 0.0
77: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
78: & '-a1',user(user_a+2),flg,ierr)
79: user(user_k+1) = 1000000.0
80: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
81: & '-k0',user(user_k+1),flg,ierr)
82: user(user_k+2) = 2*user(user_k+1)
83: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
84: & '-k1', user(user_k+2),flg,ierr)
85: user(user_s+1) = 0.0
86: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
87: & '-s0',user(user_s+1),flg,ierr)
88: user(user_s+2) = 1.0
89: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
90: & '-s1',user(user_s+2),flg,ierr)
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: ! Create timestepping solver context
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
96: call TSSetDM(ts,da,ierr)
97: call TSSetType(ts,TSARKIMEX,ierr)
98: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction, &
99: & user,ierr)
100: call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
101: call DMSetMatType(da,MATAIJ,ierr)
102: call DMCreateMatrix(da,J,ierr)
103: call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)
105: call TSGetSNES(ts,snes,ierr)
106: call SNESGetLineSearch(snes,linesearch,ierr)
107: call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)
109: ftime = 1.0
110: maxsteps = 10000
111: call TSSetDuration(ts,maxsteps,ftime,ierr)
112: call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)
113:
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Set initial conditions
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: call FormInitialSolution(ts,X,user,ierr)
118: call TSSetSolution(ts,X,ierr)
119: call VecGetSize(X,mx,ierr)
120: ! Advective CFL, I don't know why it needs so much safety factor.
121: dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
122: call TSSetInitialTimeStep(ts,zero,dt,ierr)
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Set runtime options
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: call TSSetFromOptions(ts,ierr)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: ! Solve nonlinear system
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: call TSSolve(ts,X,ierr)
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: ! Free work space.
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: call MatDestroy(J,ierr)
138: call VecDestroy(X,ierr)
139: call TSDestroy(ts,ierr)
140: call DMDestroy(da,ierr)
141: call PetscFinalize(ierr)
142: end program
144: ! Small helper to extract the layout, result uses 1-based indexing.
145: subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
146: implicit none
147: #include <petsc/finclude/petscsys.h>
148: #include <petsc/finclude/petscdm.h>
149: #include <petsc/finclude/petscdmda.h>
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: implicit none
189: #include <petsc/finclude/petscsys.h>
190: #include <petsc/finclude/petscvec.h>
191: #include <petsc/finclude/petscmat.h>
192: #include <petsc/finclude/petscsnes.h>
193: #include <petsc/finclude/petscts.h>
194: #include <petsc/finclude/petscdm.h>
195: #include <petsc/finclude/petscdmda.h>
196: TS ts
197: PetscReal t
198: Vec X,Xdot,F
199: PetscReal user(6)
200: PetscErrorCode ierr
201: integer user_a,user_k,user_s
202: parameter (user_a = 1,user_k = 3,user_s = 5)
204: DM da
205: PetscInt mx,xs,xe,gxs,gxe
206: PetscOffset ixx,ixxdot,iff
207: PetscScalar xx(0:1),xxdot(0:1),ff(0:1)
209: call TSGetDM(ts,da,ierr)
210: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
212: ! Get access to vector data
213: call VecGetArrayRead(X,xx,ixx,ierr)
214: call VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)
215: call VecGetArray(F,ff,iff,ierr)
217: call FormIFunctionLocal(mx,xs,xe,gxs,gxe, &
218: & xx(ixx),xxdot(ixxdot),ff(iff), &
219: & user(user_a),user(user_k),user(user_s),ierr)
221: call VecRestoreArrayRead(X,xx,ixx,ierr)
222: call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)
223: call VecRestoreArray(F,ff,iff,ierr)
224: end subroutine
226: subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f, &
227: & a,k,s,ierr)
228: implicit none
229: PetscInt mx,xs,xe,gxs,gxe
230: PetscReal t
231: PetscScalar x(2,gxs:gxe),f(2,xs:xe)
232: PetscReal a(2),k(2),s(2)
233: PetscErrorCode ierr
234: PetscInt i,j
235: PetscReal hx,u0t(2)
236: PetscReal one,two,three,four,six,twelve
237: PetscReal half,third,twothird,sixth
238: PetscReal twelfth
240: one = 1.0
241: two = 2.0
242: three = 3.0
243: four = 4.0
244: six = 6.0
245: twelve = 12.0
246: hx = one / mx
247: u0t(1) = one - sin(twelve*t)**four
248: u0t(2) = 0.0
249: half = one/two
250: third = one / three
251: twothird = two / three
252: sixth = one / six
253: twelfth = one / twelve
254: do 20 i = xs,xe
255: do 10 j = 1,2
256: if (i .eq. 1) then
257: f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) &
258: & + sixth*x(j,i+2))
259: else if (i .eq. 2) then
260: f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) &
261: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
262: else if (i .eq. mx-1) then
263: f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) &
264: & - half*x(j,i) -third*x(j,i+1))
265: else if (i .eq. mx) then
266: f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
267: else
268: f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) &
269: & + twothird*x(j,i-1) &
270: & - twothird*x(j,i+1) + twelfth*x(j,i+2))
271: end if
272: 10 continue
273: 20 continue
274: end subroutine
276: subroutine FormRHSFunction(ts,t,X,F,user,ierr)
277: implicit none
278: #include <petsc/finclude/petscsys.h>
279: #include <petsc/finclude/petscvec.h>
280: #include <petsc/finclude/petscmat.h>
281: #include <petsc/finclude/petscsnes.h>
282: #include <petsc/finclude/petscts.h>
283: #include <petsc/finclude/petscdm.h>
284: #include <petsc/finclude/petscdmda.h>
285: TS ts
286: PetscReal t
287: Vec X,F
288: PetscReal user(6)
289: PetscErrorCode ierr
290: integer user_a,user_k,user_s
291: parameter (user_a = 1,user_k = 3,user_s = 5)
292: DM da
293: Vec Xloc
294: PetscInt mx,xs,xe,gxs,gxe
295: PetscOffset ixx,iff
296: PetscScalar xx(0:1),ff(0:1)
298: call TSGetDM(ts,da,ierr)
299: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
301: ! Scatter ghost points to local vector,using the 2-step process
302: ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
303: ! By placing code between these two statements, computations can be
304: ! done while messages are in transition.
305: call DMGetLocalVector(da,Xloc,ierr)
306: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
307: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)
309: ! Get access to vector data
310: call VecGetArrayRead(Xloc,xx,ixx,ierr)
311: call VecGetArray(F,ff,iff,ierr)
313: call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe, &
314: & t,xx(ixx),ff(iff), &
315: & user(user_a),user(user_k),user(user_s),ierr)
317: call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
318: call VecRestoreArray(F,ff,iff,ierr)
319: call DMRestoreLocalVector(da,Xloc,ierr)
320: end subroutine
322: ! ---------------------------------------------------------------------
323: !
324: ! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
325: !
326: subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
327: implicit none
328: #include <petsc/finclude/petscsys.h>
329: #include <petsc/finclude/petscvec.h>
330: #include <petsc/finclude/petscmat.h>
331: #include <petsc/finclude/petscsnes.h>
332: #include <petsc/finclude/petscts.h>
333: #include <petsc/finclude/petscdm.h>
334: #include <petsc/finclude/petscdmda.h>
335: TS ts
336: PetscReal t,shift
337: Vec X,Xdot
338: Mat J,Jpre
339: PetscReal user(6)
340: PetscErrorCode ierr
341: integer user_a,user_k,user_s
342: parameter (user_a = 0,user_k = 2,user_s = 4)
344: DM da
345: PetscInt mx,xs,xe,gxs,gxe
346: PetscInt i,i1,row,col
347: PetscReal k1,k2;
348: PetscScalar val(4)
350: call TSGetDM(ts,da,ierr)
351: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
353: i1 = 1
354: k1 = user(user_k+1)
355: k2 = user(user_k+2)
356: do 10 i = xs,xe
357: row = i-gxs
358: col = i-gxs
359: val(1) = shift + k1
360: val(2) = -k2
361: val(3) = -k1
362: val(4) = shift + k2
363: call MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val, &
364: & INSERT_VALUES,ierr)
365: 10 continue
366: call MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)
367: call MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)
368: if (J /= Jpre) then
369: call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)
370: call MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)
371: end if
372: end subroutine
374: subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
375: implicit none
376: PetscInt mx,xs,xe,gxs,gxe
377: PetscScalar x(2,xs:xe)
378: PetscReal a(2),k(2),s(2)
379: PetscErrorCode ierr
381: PetscInt i
382: PetscReal one,hx,r,ik
383: one = 1.0
384: hx = one / mx
385: do 10 i=xs,xe
386: r = i*hx
387: if (k(2) .ne. 0.0) then
388: ik = one/k(2)
389: else
390: ik = one
391: end if
392: x(1,i) = one + s(2)*r
393: x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
394: 10 continue
395: end subroutine
397: subroutine FormInitialSolution(ts,X,user,ierr)
398: implicit none
399: #include <petsc/finclude/petscsys.h>
400: #include <petsc/finclude/petscvec.h>
401: #include <petsc/finclude/petscmat.h>
402: #include <petsc/finclude/petscsnes.h>
403: #include <petsc/finclude/petscts.h>
404: #include <petsc/finclude/petscdm.h>
405: #include <petsc/finclude/petscdmda.h>
406: TS ts
407: Vec X
408: PetscReal user(6)
409: PetscErrorCode ierr
410: integer user_a,user_k,user_s
411: parameter (user_a = 1,user_k = 3,user_s = 5)
413: DM da
414: PetscInt mx,xs,xe,gxs,gxe
415: PetscOffset ixx
416: PetscScalar xx(0:1)
418: call TSGetDM(ts,da,ierr)
419: call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
421: ! Get access to vector data
422: call VecGetArray(X,xx,ixx,ierr)
424: call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe, &
425: & xx(ixx),user(user_a),user(user_k),user(user_s),ierr)
427: call VecRestoreArray(X,xx,ixx,ierr)
428: end subroutine