Actual source code: ex22f.F

  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:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 55:       if (ierr .ne. 0) then
 56:          print*,'Unable to initialize PETSc'
 57:          stop
 58:       endif

 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)

 68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69: !    Extract global vectors from DMDA;
 70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:       call DMCreateGlobalVector(da,X,ierr)

 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)

 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)

107:       call TSGetSNES(ts,snes,ierr)
108:       call SNESGetLineSearch(snes,linesearch,ierr)
109:       call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)

111:       ftime = 1.0
112:       call TSSetMaxTime(ts,ftime,ierr)
113:       call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)

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)

125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: !   Set runtime options
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:       call TSSetFromOptions(ts,ierr)

130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: !  Solve nonlinear system
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:       call TSSolve(ts,X,ierr)

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

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

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:       use petscts
189:       implicit none

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 VecGetArrayRead(X,xx,ixx,ierr)
209:       call VecGetArrayRead(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 VecRestoreArrayRead(X,xx,ixx,ierr)
217:       call VecRestoreArrayRead(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: !     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

272:       subroutine FormRHSFunction(ts,t,X,F,user,ierr)
273:       use petscts
274:       implicit none

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)

289:       call TSGetDM(ts,da,ierr)
290:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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)

300: ! Get access to vector data
301:       call VecGetArrayRead(Xloc,xx,ixx,ierr)
302:       call VecGetArray(F,ff,iff,ierr)

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)

308:       call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
309:       call VecRestoreArray(F,ff,iff,ierr)
310:       call DMRestoreLocalVector(da,Xloc,ierr)
311:       end subroutine

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

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)

330:       DM             da
331:       PetscInt       mx,xs,xe,gxs,gxe
332:       PetscInt       i,i1,row,col
333:       PetscReal      k1,k2;
334:       PetscScalar    val(4)

336:       call TSGetDM(ts,da,ierr)
337:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

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

383:       subroutine FormInitialSolution(ts,X,user,ierr)
384:       use petscts
385:       implicit none

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)

394:       DM             da
395:       PetscInt       mx,xs,xe,gxs,gxe
396:       PetscOffset    ixx
397:       PetscScalar    xx(0:1)

399:       call TSGetDM(ts,da,ierr)
400:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

402: ! Get access to vector data
403:       call VecGetArray(X,xx,ixx,ierr)

405:       call FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,                       &
406:      &     xx(ixx),user(user_a),user(user_k),user(user_s),ierr)

408:       call VecRestoreArray(X,xx,ixx,ierr)
409:       end subroutine

411: !/*TEST
412: !
413: !    test:
414: !      args: -da_grid_x 200 -ts_arkimex_type 4
415: !      requires: !single
416: !
417: !TEST*/