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:   TS ts
 34:   SNES snes
 35:   SNESLineSearch linesearch
 36:   Vec X
 37:   Mat J
 38:   PetscInt mx
 39:   PetscErrorCode ierr
 40:   DM da
 41:   PetscReal ftime, dt
 42:   PetscReal one, pone
 43:   PetscInt im11, i2
 44:   PetscBool flg

 46:   im11 = 11
 47:   i2 = 2
 48:   one = 1.0
 49:   pone = one/10

 51:   PetscCallA(PetscInitialize(ierr))

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

 60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61: !    Extract global vectors from DMDA
 62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:   PetscCallA(DMCreateGlobalVector(da, X, ierr))

 65: ! Initialize user application context
 66: ! Use zero-based indexing for command line parameters to match ex22.c
 67:   user(user_a + 1) = 1.0
 68:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', user(user_a + 1), flg, ierr))
 69:   user(user_a + 2) = 0.0
 70:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', user(user_a + 2), flg, ierr))
 71:   user(user_k + 1) = 1000000.0
 72:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', user(user_k + 1), flg, ierr))
 73:   user(user_k + 2) = 2*user(user_k + 1)
 74:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', user(user_k + 2), flg, ierr))
 75:   user(user_s + 1) = 0.0
 76:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', user(user_s + 1), flg, ierr))
 77:   user(user_s + 2) = 1.0
 78:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', user(user_s + 2), flg, ierr))

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

 92:   PetscCallA(TSGetSNES(ts, snes, ierr))
 93:   PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
 94:   PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr))

 96:   ftime = 1.0
 97:   PetscCallA(TSSetMaxTime(ts, ftime, ierr))
 98:   PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))

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

110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: !   Set runtime options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:   PetscCallA(TSSetFromOptions(ts, ierr))

115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: !  Solve nonlinear system
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:   PetscCallA(TSSolve(ts, X, ierr))

120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: !  Free work space.
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:   PetscCallA(MatDestroy(J, ierr))
124:   PetscCallA(VecDestroy(X, ierr))
125:   PetscCallA(TSDestroy(ts, ierr))
126:   PetscCallA(DMDestroy(da, ierr))
127:   PetscCallA(PetscFinalize(ierr))
128: contains

130: ! Small helper to extract the layout, result uses 1-based indexing.
131:   subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
132:     use petscdm
133:     implicit none

135:     DM da
136:     PetscInt mx, xs, xe, gxs, gxe
137:     PetscErrorCode ierr
138:     PetscInt xm, gxm
139:     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_ENUM, PETSC_NULL_ENUM, PETSC_NULL_ENUM, PETSC_NULL_ENUM, ierr))
140:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
141:     PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
142:     xs = xs + 1
143:     gxs = gxs + 1
144:     xe = xs + xm - 1
145:     gxe = gxs + gxm - 1
146:   end subroutine

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

163:       subroutine FormIFunction(ts, t, X, Xdot, F, user, ierr)
164:         use petscts
165:         implicit none

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

175:         DM da
176:         PetscInt mx, xs, xe, gxs, gxe
177:         PetscScalar, pointer :: xx(:), xxdot(:), ff(:)

179:         PetscCall(TSGetDM(ts, da, ierr))
180:         PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

182: ! Get access to vector data
183:         PetscCall(VecGetArrayRead(X, xx, ierr))
184:         PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
185:         PetscCall(VecGetArray(F, ff, ierr))

187:         PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, user(user_a), user(user_k), user(user_s), ierr))

189:         PetscCall(VecRestoreArrayRead(X, xx, ierr))
190:         PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
191:         PetscCall(VecRestoreArray(F, ff, ierr))
192:       end subroutine

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

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

244:             subroutine FormRHSFunction(ts, t, X, F, user, ierr)
245:               use petscts
246:               implicit none

248:               TS ts
249:               PetscReal t
250:               Vec X, F
251:               PetscReal user(6)
252:               PetscErrorCode ierr
253:               integer user_a, user_k, user_s
254:               parameter(user_a=1, user_k=3, user_s=5)
255:               DM da
256:               Vec Xloc
257:               PetscInt mx, xs, xe, gxs, gxe
258:               PetscScalar, pointer :: xx(:), ff(:)

260:               PetscCall(TSGetDM(ts, da, ierr))
261:               PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

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

271: ! Get access to vector data
272:               PetscCall(VecGetArrayRead(Xloc, xx, ierr))
273:               PetscCall(VecGetArray(F, ff, ierr))

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

277:               PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
278:               PetscCall(VecRestoreArray(F, ff, ierr))
279:               PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
280:             end subroutine

282: ! ---------------------------------------------------------------------
283: !
284: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
285: !
286:             subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, user, ierr)
287:               use petscts
288:               implicit none

290:               TS ts
291:               PetscReal t, shift
292:               Vec X, Xdot
293:               Mat J, Jpre
294:               PetscReal user(6)
295:               PetscErrorCode ierr
296:               integer user_a, user_k, user_s
297:               parameter(user_a=0, user_k=2, user_s=4)

299:               DM da
300:               PetscInt mx, xs, xe, gxs, gxe
301:               PetscInt i, i1, row, col
302:               PetscReal k1, k2
303:               PetscScalar val(4)

305:               PetscCall(TSGetDM(ts, da, ierr))
306:               PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

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

328:                 subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
329:                   implicit none
330:                   PetscInt mx, xs, xe, gxs, gxe
331:                   PetscScalar x(2, xs:xe)
332:                   PetscReal a(2), k(2), s(2)
333:                   PetscErrorCode ierr

335:                   PetscInt i
336:                   PetscReal one, hx, r, ik
337:                   one = 1.0
338:                   hx = one/mx
339:                   do 10 i = xs, xe
340:                     r = i*hx
341:                     if (k(2) /= 0.0) then
342:                       ik = one/k(2)
343:                     else
344:                       ik = one
345:                     end if
346:                     x(1, i) = one + s(2)*r
347:                     x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
348: 10                  continue
349:                     end subroutine

351:                     subroutine FormInitialSolution(ts, X, user, ierr)
352:                       use petscts
353:                       implicit none

355:                       TS ts
356:                       Vec X
357:                       PetscReal user(6)
358:                       PetscErrorCode ierr
359:                       integer user_a, user_k, user_s
360:                       parameter(user_a=1, user_k=3, user_s=5)

362:                       DM da
363:                       PetscInt mx, xs, xe, gxs, gxe
364:                       PetscScalar, pointer :: xx(:)

366:                       PetscCall(TSGetDM(ts, da, ierr))
367:                       PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

369: ! Get access to vector data
370:                       PetscCall(VecGetArray(X, xx, ierr))

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

374:                       PetscCall(VecRestoreArray(X, xx, ierr))
375:                     end subroutine
376:                     end program
377: !/*TEST
378: !
379: !    test:
380: !      args: -da_grid_x 200 -ts_arkimex_type 4
381: !      requires: !single
382: !      output_file: output/empty.out
383: !
384: !TEST*/