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!
16
17      program main
18#include <petsc/finclude/petscts.h>
19#include <petsc/finclude/petscdmda.h>
20      use petscts
21      implicit none
22
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)
32
33      external FormRHSFunction,FormIFunction,FormIJacobian
34      external FormInitialSolution
35
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
48
49      im11 = 11
50      i2   = 2
51      one = 1.0
52      pone = one / 10
53
54      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
55      if (ierr .ne. 0) then
56         print*,'Unable to initialize PETSc'
57         stop
58      endif
59
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)
67
68! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69!    Extract global vectors from DMDA;
70! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71      call DMCreateGlobalVector(da,X,ierr)
72
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)
93
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)
106
107      call TSGetSNES(ts,snes,ierr)
108      call SNESGetLineSearch(snes,linesearch,ierr)
109      call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)
110
111      ftime = 1.0
112      call TSSetMaxTime(ts,ftime,ierr)
113      call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)
114
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)
124
125! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126!   Set runtime options
127! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128      call TSSetFromOptions(ts,ierr)
129
130! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131!  Solve nonlinear system
132! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133      call TSSolve(ts,X,ierr)
134
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
144
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
149
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
170
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
186
187      subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
188      use petscts
189      implicit none
190
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)
198
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)
203
204      call TSGetDM(ts,da,ierr)
205      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
206
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)
211
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)
215
216      call VecRestoreArrayRead(X,xx,ixx,ierr)
217      call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)
218      call VecRestoreArray(F,ff,iff,ierr)
219      end subroutine
220
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
234
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
271
272      subroutine FormRHSFunction(ts,t,X,F,user,ierr)
273      use petscts
274      implicit none
275
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)
288
289      call TSGetDM(ts,da,ierr)
290      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
291
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)
299
300! Get access to vector data
301      call VecGetArrayRead(Xloc,xx,ixx,ierr)
302      call VecGetArray(F,ff,iff,ierr)
303
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)
307
308      call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
309      call VecRestoreArray(F,ff,iff,ierr)
310      call DMRestoreLocalVector(da,Xloc,ierr)
311      end subroutine
312
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
320
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)
329
330      DM             da
331      PetscInt       mx,xs,xe,gxs,gxe
332      PetscInt       i,i1,row,col
333      PetscReal      k1,k2;
334      PetscScalar    val(4)
335
336      call TSGetDM(ts,da,ierr)
337      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
338
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
359
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
366
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
382
383      subroutine FormInitialSolution(ts,X,user,ierr)
384      use petscts
385      implicit none
386
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)
393
394      DM             da
395      PetscInt       mx,xs,xe,gxs,gxe
396      PetscOffset    ixx
397      PetscScalar    xx(0:1)
398
399      call TSGetDM(ts,da,ierr)
400      call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
401
402! Get access to vector data
403      call VecGetArray(X,xx,ixx,ierr)
404
405      call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,                       &
406     &     xx(ixx),user(user_a),user(user_k),user(user_s),ierr)
407
408      call VecRestoreArray(X,xx,ixx,ierr)
409      end subroutine
410
411!/*TEST
412!
413!    test:
414!      args: -da_grid_x 200 -ts_arkimex_type 4
415!      requires: !single
416!
417!TEST*/
418