Actual source code: ex22f.F90

  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:       PetscCallA(PetscInitialize(ierr))

 56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57: !  Create distributed array (DMDA) to manage parallel grid and vectors
 58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 59:       PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr))
 60:       PetscCallA(DMSetFromOptions(da,ierr))
 61:       PetscCallA(DMSetUp(da,ierr))

 63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64: !    Extract global vectors from DMDA;
 65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66:       PetscCallA(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:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr))
 72:       user(user_a+2) = 0.0
 73:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr))
 74:       user(user_k+1) = 1000000.0
 75:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0',user(user_k+1),flg,ierr))
 76:       user(user_k+2) = 2*user(user_k+1)
 77:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr))
 78:       user(user_s+1) = 0.0
 79:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr))
 80:       user(user_s+2) = 1.0
 81:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr))

 83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84: !    Create timestepping solver context
 85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86:       PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
 87:       PetscCallA(TSSetDM(ts,da,ierr))
 88:       PetscCallA(TSSetType(ts,TSARKIMEX,ierr))
 89:       PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr))
 90:       PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr))
 91:       PetscCallA(DMSetMatType(da,MATAIJ,ierr))
 92:       PetscCallA(DMCreateMatrix(da,J,ierr))
 93:       PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr))

 95:       PetscCallA(TSGetSNES(ts,snes,ierr))
 96:       PetscCallA(SNESGetLineSearch(snes,linesearch,ierr))
 97:       PetscCallA(SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr))

 99:       ftime = 1.0
100:       PetscCallA(TSSetMaxTime(ts,ftime,ierr))
101:       PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))

103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: !  Set initial conditions
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:       PetscCallA(FormInitialSolution(ts,X,user,ierr))
107:       PetscCallA(TSSetSolution(ts,X,ierr))
108:       PetscCallA(VecGetSize(X,mx,ierr))
109: !  Advective CFL, I don't know why it needs so much safety factor.
110:       dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
111:       PetscCallA(TSSetTimeStep(ts,dt,ierr))

113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: !   Set runtime options
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:       PetscCallA(TSSetFromOptions(ts,ierr))

118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: !  Solve nonlinear system
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:       PetscCallA(TSSolve(ts,X,ierr))

123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: !  Free work space.
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:       PetscCallA(MatDestroy(J,ierr))
127:       PetscCallA(VecDestroy(X,ierr))
128:       PetscCallA(TSDestroy(ts,ierr))
129:       PetscCallA(DMDestroy(da,ierr))
130:       PetscCallA(PetscFinalize(ierr))
131:       end program

133: ! Small helper to extract the layout, result uses 1-based indexing.
134:       subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
135:       use petscdmda
136:       implicit none

138:       DM da
139:       PetscInt mx,xs,xe,gxs,gxe
140:       PetscErrorCode ierr
141:       PetscInt xm,gxm
142:       PetscCall(DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
143:       PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
144:       PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
145:       xs = xs + 1
146:       gxs = gxs + 1
147:       xe = xs + xm - 1
148:       gxe = gxs + gxm - 1
149:       end subroutine

151:       subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr)
152:       implicit none
153:       PetscInt mx,xs,xe,gxs,gxe
154:       PetscScalar x(2,xs:xe)
155:       PetscScalar xdot(2,xs:xe)
156:       PetscScalar f(2,xs:xe)
157:       PetscReal a(2),k(2),s(2)
158:       PetscErrorCode ierr
159:       PetscInt i
160:       do 10 i = xs,xe
161:          f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
162:          f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
163:  10   continue
164:       end subroutine

166:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
167:       use petscts
168:       implicit none

170:       TS ts
171:       PetscReal t
172:       Vec X,Xdot,F
173:       PetscReal user(6)
174:       PetscErrorCode ierr
175:       integer user_a,user_k,user_s
176:       parameter (user_a = 1,user_k = 3,user_s = 5)

178:       DM             da
179:       PetscInt       mx,xs,xe,gxs,gxe
180:       PetscOffset    ixx,ixxdot,iff
181:       PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

183:       PetscCall(TSGetDM(ts,da,ierr))
184:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

186: ! Get access to vector data
187:       PetscCall(VecGetArrayRead(X,xx,ixx,ierr))
188:       PetscCall(VecGetArrayRead(Xdot,xxdot,ixxdot,ierr))
189:       PetscCall(VecGetArray(F,ff,iff,ierr))

191:       PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr))

193:       PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr))
194:       PetscCall(VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr))
195:       PetscCall(VecRestoreArray(F,ff,iff,ierr))
196:       end subroutine

198:       subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr)
199:       implicit none
200:       PetscInt mx,xs,xe,gxs,gxe
201:       PetscReal t
202:       PetscScalar x(2,gxs:gxe),f(2,xs:xe)
203:       PetscReal a(2),k(2),s(2)
204:       PetscErrorCode ierr
205:       PetscInt i,j
206:       PetscReal hx,u0t(2)
207:       PetscReal one,two,three,four,six,twelve
208:       PetscReal half,third,twothird,sixth
209:       PetscReal twelfth

211:       one = 1.0
212:       two = 2.0
213:       three = 3.0
214:       four = 4.0
215:       six = 6.0
216:       twelve = 12.0
217:       hx = one / mx
218: !     The Fortran standard only allows positive base for power functions; Nag compiler fails on this
219:       u0t(1) = one - abs(sin(twelve*t))**four
220:       u0t(2) = 0.0
221:       half = one/two
222:       third = one / three
223:       twothird = two / three
224:       sixth = one / six
225:       twelfth = one / twelve
226:       do 20 i = xs,xe
227:          do 10 j = 1,2
228:             if (i .eq. 1) then
229:                f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
230:      &              + sixth*x(j,i+2))
231:             else if (i .eq. 2) then
232:                f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
233:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
234:             else if (i .eq. mx-1) then
235:                f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1)             &
236:      &         - half*x(j,i) -third*x(j,i+1))
237:             else if (i .eq. mx) then
238:                f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
239:             else
240:                f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
241:      &              + twothird*x(j,i-1)                                 &
242:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
243:             end if
244:  10      continue
245:  20   continue
246:       end subroutine

248:       subroutine FormRHSFunction(ts,t,X,F,user,ierr)
249:       use petscts
250:       implicit none

252:       TS ts
253:       PetscReal t
254:       Vec X,F
255:       PetscReal user(6)
256:       PetscErrorCode ierr
257:       integer user_a,user_k,user_s
258:       parameter (user_a = 1,user_k = 3,user_s = 5)
259:       DM             da
260:       Vec            Xloc
261:       PetscInt       mx,xs,xe,gxs,gxe
262:       PetscOffset    ixx,iff
263:       PetscScalar    xx(0:1),ff(0:1)

265:       PetscCall(TSGetDM(ts,da,ierr))
266:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

268: !     Scatter ghost points to local vector,using the 2-step process
269: !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
270: !     By placing code between these two statements, computations can be
271: !     done while messages are in transition.
272:       PetscCall(DMGetLocalVector(da,Xloc,ierr))
273:       PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr))
274:       PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr))

276: ! Get access to vector data
277:       PetscCall(VecGetArrayRead(Xloc,xx,ixx,ierr))
278:       PetscCall(VecGetArray(F,ff,iff,ierr))

280:       PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx(ixx),ff(iff),user(user_a),user(user_k),user(user_s),ierr))

282:       PetscCall(VecRestoreArrayRead(Xloc,xx,ixx,ierr))
283:       PetscCall(VecRestoreArray(F,ff,iff,ierr))
284:       PetscCall(DMRestoreLocalVector(da,Xloc,ierr))
285:       end subroutine

287: ! ---------------------------------------------------------------------
288: !
289: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
290: !
291:       subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr)
292:       use petscts
293:       implicit none

295:       TS ts
296:       PetscReal t,shift
297:       Vec X,Xdot
298:       Mat J,Jpre
299:       PetscReal user(6)
300:       PetscErrorCode ierr
301:       integer user_a,user_k,user_s
302:       parameter (user_a = 0,user_k = 2,user_s = 4)

304:       DM             da
305:       PetscInt       mx,xs,xe,gxs,gxe
306:       PetscInt       i,i1,row,col
307:       PetscReal      k1,k2;
308:       PetscScalar    val(4)

310:       PetscCall(TSGetDM(ts,da,ierr))
311:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

313:       i1 = 1
314:       k1 = user(user_k+1)
315:       k2 = user(user_k+2)
316:       do 10 i = xs,xe
317:          row = i-gxs
318:          col = i-gxs
319:          val(1) = shift + k1
320:          val(2) = -k2
321:          val(3) = -k1
322:          val(4) = shift + k2
323:          PetscCall(MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,INSERT_VALUES,ierr))
324:  10   continue
325:       PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
326:       PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
327:       if (J /= Jpre) then
328:          PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr))
329:          PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr))
330:       end if
331:       end subroutine

333:       subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
334:       implicit none
335:       PetscInt mx,xs,xe,gxs,gxe
336:       PetscScalar x(2,xs:xe)
337:       PetscReal a(2),k(2),s(2)
338:       PetscErrorCode ierr

340:       PetscInt i
341:       PetscReal one,hx,r,ik
342:       one = 1.0
343:       hx = one / mx
344:       do 10 i=xs,xe
345:          r = i*hx
346:          if (k(2) .ne. 0.0) then
347:             ik = one/k(2)
348:          else
349:             ik = one
350:          end if
351:          x(1,i) = one + s(2)*r
352:          x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
353:  10   continue
354:       end subroutine

356:       subroutine FormInitialSolution(ts,X,user,ierr)
357:       use petscts
358:       implicit none

360:       TS ts
361:       Vec X
362:       PetscReal user(6)
363:       PetscErrorCode ierr
364:       integer user_a,user_k,user_s
365:       parameter (user_a = 1,user_k = 3,user_s = 5)

367:       DM             da
368:       PetscInt       mx,xs,xe,gxs,gxe
369:       PetscOffset    ixx
370:       PetscScalar    xx(0:1)

372:       PetscCall(TSGetDM(ts,da,ierr))
373:       PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr))

375: ! Get access to vector data
376:       PetscCall(VecGetArray(X,xx,ixx,ierr))

378:       PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx(ixx),user(user_a),user(user_k),user(user_s),ierr))

380:       PetscCall(VecRestoreArray(X,xx,ixx,ierr))
381:       end subroutine

383: !/*TEST
384: !
385: !    test:
386: !      args: -da_grid_x 200 -ts_arkimex_type 4
387: !      requires: !single
388: !
389: !TEST*/