Actual source code: ex22f_mf.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: module ex22f_mfmodule
 18: #include <petsc/finclude/petscts.h>
 19:   use petscts
 20:   PetscScalar::PETSC_SHIFT
 21:   TS::tscontext
 22:   Mat::Jmat
 23:   PetscReal::MFuser(6)
 24: end module ex22f_mfmodule

 26: program main
 27:   use ex22f_mfmodule
 28:   use petscdm
 29:   implicit none

 31:   !
 32:   !     Create an application context to contain data needed by the
 33:   !     application-provided call-back routines, FormJacobian() and
 34:   !     FormFunction(). We use a double precision array with six
 35:   !     entries, two for each problem parameter a, k, s.
 36:   !
 37:   PetscReal user(6)
 38:   integer user_a, user_k, user_s
 39:   parameter(user_a=0, user_k=2, user_s=4)

 41:   TS ts
 42:   Vec X
 43:   Mat J
 44:   PetscInt mx
 45:   PetscBool OptionSaveToDisk
 46:   PetscErrorCode ierr
 47:   DM da
 48:   PetscReal ftime, dt
 49:   PetscReal one, pone
 50:   PetscInt im11, i2
 51:   PetscBool flg

 53:   PetscInt xs, xe, gxs, gxe, dof, gdof
 54:   PetscScalar shell_shift
 55:   Mat A

 57:   im11 = 11
 58:   i2 = 2
 59:   one = 1.0
 60:   pone = one/10

 62:   PetscCallA(PetscInitialize(ierr))

 64:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:   !  Create distributed array (DMDA) to manage parallel grid and vectors
 66:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, im11, i2, i2, PETSC_NULL_INTEGER_ARRAY, da, ierr))
 68:   PetscCallA(DMSetFromOptions(da, ierr))
 69:   PetscCallA(DMSetUp(da, ierr))

 71:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 72:   !    Extract global vectors from DMDA
 73:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:   PetscCallA(DMCreateGlobalVector(da, X, ierr))

 76:   ! Initialize user application context
 77:   ! Use zero-based indexing for command line parameters to match ex22.c
 78:   user(user_a + 1) = 1.0
 79:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', user(user_a + 1), flg, ierr))
 80:   user(user_a + 2) = 0.0
 81:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', user(user_a + 2), flg, ierr))
 82:   user(user_k + 1) = 1000000.0
 83:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', user(user_k + 1), flg, ierr))
 84:   user(user_k + 2) = 2*user(user_k + 1)
 85:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', user(user_k + 2), flg, ierr))
 86:   user(user_s + 1) = 0.0
 87:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', user(user_s + 1), flg, ierr))
 88:   user(user_s + 2) = 1.0
 89:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', user(user_s + 2), flg, ierr))

 91:   OptionSaveToDisk = .FALSE.
 92:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-sdisk', OptionSaveToDisk, flg, ierr))
 93:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:   !    Create timestepping solver context
 95:   !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:   PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
 97:   tscontext = ts
 98:   PetscCallA(TSSetDM(ts, da, ierr))
 99:   PetscCallA(TSSetType(ts, TSARKIMEX, ierr))
100:   PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, user, ierr))

102:   ! - - - - - - - - -- - - - -
103:   !   Matrix free setup
104:   PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
105:   dof = i2*(xe - xs + 1)
106:   gdof = i2*(gxe - gxs + 1)
107:   PetscCallA(MatCreateShell(PETSC_COMM_WORLD, dof, dof, gdof, gdof, shell_shift, A, ierr))

109:   PetscCallA(MatShellSetOperation(A, MATOP_MULT, MyMult, ierr))
110:   ! - - - - - - - - - - - -

112:   PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, user, ierr))
113:   PetscCallA(DMSetMatType(da, MATAIJ, ierr))
114:   PetscCallA(DMCreateMatrix(da, J, ierr))

116:   Jmat = J

118:   PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, user, ierr))
119:   PetscCallA(TSSetIJacobian(ts, A, A, FormIJacobianMF, user, ierr))

121:   ftime = 1.0
122:   PetscCallA(TSSetMaxTime(ts, ftime, ierr))
123:   PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))

125:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:   !  Set initial conditions
127:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:   PetscCallA(FormInitialSolution(ts, X, user, ierr))
129:   PetscCallA(TSSetSolution(ts, X, ierr))
130:   PetscCallA(VecGetSize(X, mx, ierr))
131:   !  Advective CFL, I don't know why it needs so much safety factor.
132:   dt = pone*max(user(user_a + 1), user(user_a + 2))/mx
133:   PetscCallA(TSSetTimeStep(ts, dt, ierr))

135:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:   !   Set runtime options
137:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:   PetscCallA(TSSetFromOptions(ts, ierr))

140:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:   !  Solve nonlinear system
142:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:   PetscCallA(TSSolve(ts, X, ierr))

145:   if (OptionSaveToDisk) then
146:     PetscCallA(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))
147:     dof = i2*(xe - xs + 1)
148:     gdof = i2*(gxe - gxs + 1)
149:     call SaveSolutionToDisk(da, X, gdof, xs, xe)
150:   end if

152:   ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:   !  Free work space.
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:   PetscCallA(MatDestroy(A, ierr))
156:   PetscCallA(MatDestroy(J, ierr))
157:   PetscCallA(VecDestroy(X, ierr))
158:   PetscCallA(TSDestroy(ts, ierr))
159:   PetscCallA(DMDestroy(da, ierr))
160:   PetscCallA(PetscFinalize(ierr))
161: contains

163: ! Small helper to extract the layout, result uses 1-based indexing.
164:   subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr)
165:     use petscdm
166:     implicit none

168:     DM da
169:     PetscInt mx, xs, xe, gxs, gxe
170:     PetscErrorCode ierr
171:     PetscInt xm, gxm
172:     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))
173:     PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
174:     PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
175:     xs = xs + 1
176:     gxs = gxs + 1
177:     xe = xs + xm - 1
178:     gxe = gxs + gxm - 1
179:   end subroutine GetLayout

181:   subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr)
182:     implicit none
183:     PetscInt mx, xs, xe, gxs, gxe
184:     PetscScalar x(2, xs:xe)
185:     PetscScalar xdot(2, xs:xe)
186:     PetscScalar f(2, xs:xe)
187:     PetscReal a(2), k(2), s(2)
188:     PetscErrorCode ierr
189:     PetscInt i
190:     do i = xs, xe
191:       f(1, i) = xdot(1, i) + k(1)*x(1, i) - k(2)*x(2, i) - s(1)
192:       f(2, i) = xdot(2, i) - k(1)*x(1, i) + k(2)*x(2, i) - s(2)
193:     end do
194:   end subroutine FormIFunctionLocal

196:   subroutine FormIFunction(ts, t, X, Xdot, F, user, ierr)
197:     use petscdm
198:     use petscts
199:     implicit none

201:     TS ts
202:     PetscReal t
203:     Vec X, Xdot, F
204:     PetscReal user(6)
205:     PetscErrorCode ierr
206:     integer user_a, user_k, user_s
207:     parameter(user_a=1, user_k=3, user_s=5)

209:     DM da
210:     PetscInt mx, xs, xe, gxs, gxe
211:     PetscScalar, pointer :: xx(:), xxdot(:), ff(:)

213:     PetscCall(TSGetDM(ts, da, ierr))
214:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

216:     ! Get access to vector data
217:     PetscCall(VecGetArrayRead(X, xx, ierr))
218:     PetscCall(VecGetArrayRead(Xdot, xxdot, ierr))
219:     PetscCall(VecGetArray(F, ff, ierr))

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

223:     PetscCall(VecRestoreArrayRead(X, xx, ierr))
224:     PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr))
225:     PetscCall(VecRestoreArray(F, ff, ierr))
226:   end subroutine FormIFunction

228:   subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr)
229:     implicit none
230:     PetscInt mx, xs, xe, gxs, gxe
231:     PetscReal t
232:     PetscScalar x(2, gxs:gxe), f(2, xs:xe)
233:     PetscReal a(2), k(2), s(2)
234:     PetscErrorCode ierr
235:     PetscInt i, j
236:     PetscReal hx, u0t(2)
237:     PetscReal one, two, three, four, six, twelve
238:     PetscReal half, third, twothird, sixth
239:     PetscReal twelfth

241:     one = 1.0
242:     two = 2.0
243:     three = 3.0
244:     four = 4.0
245:     six = 6.0
246:     twelve = 12.0
247:     hx = one/mx
248:     u0t(1) = one - sin(twelve*t)**four
249:     u0t(2) = 0.0
250:     half = one/two
251:     third = one/three
252:     twothird = two/three
253:     sixth = one/six
254:     twelfth = one/twelve
255:     do i = xs, xe
256:       do j = 1, 2
257:         if (i == 1) then
258:           f(j, i) = a(j)/hx*(third*u0t(j) + half*x(j, i) - x(j, i + 1) + sixth*x(j, i + 2))
259:         else if (i == 2) then
260:           f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
261:         else if (i == mx - 1) then
262:           f(j, i) = a(j)/hx*(-sixth*x(j, i - 2) + x(j, i - 1) - half*x(j, i) - third*x(j, i + 1))
263:         else if (i == mx) then
264:           f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1))
265:         else
266:           f(j, i) = a(j)/hx*(-twelfth*x(j, i - 2) + twothird*x(j, i - 1) - twothird*x(j, i + 1) + twelfth*x(j, i + 2))
267:         end if
268:       end do
269:     end do

271: #ifdef EXPLICIT_INTEGRATOR22
272:     do i = xs, xe
273:       f(1, i) = f(1, i) - (k(1)*x(1, i) - k(2)*x(2, i) - s(1))
274:       f(2, i) = f(2, i) - (-k(1)*x(1, i) + k(2)*x(2, i) - s(2))
275:     end do
276: #endif

278:   end subroutine FormRHSFunctionLocal

280:   subroutine FormRHSFunction(ts, t, X, F, user, ierr)
281:     use petscts
282:     use petscdm
283:     implicit none

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:     PetscScalar, pointer :: xx(:), ff(:)

297:     PetscCall(TSGetDM(ts, da, ierr))
298:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

300:     !     Scatter ghost points to local vector,using the 2-step process
301:     !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
302:     !     By placing code between these two statements, computations can be
303:     !     done while messages are in transition.
304:     PetscCall(DMGetLocalVector(da, Xloc, ierr))
305:     PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr))
306:     PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr))

308:     ! Get access to vector data
309:     PetscCall(VecGetArrayRead(Xloc, xx, ierr))
310:     PetscCall(VecGetArray(F, ff, ierr))

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

314:     PetscCall(VecRestoreArrayRead(Xloc, xx, ierr))
315:     PetscCall(VecRestoreArray(F, ff, ierr))
316:     PetscCall(DMRestoreLocalVector(da, Xloc, ierr))
317:   end subroutine FormRHSFunction

319: ! ---------------------------------------------------------------------
320: !
321: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
322: !
323:   subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, user, ierr)
324:     use petscts
325:     use petscdm
326:     implicit none

328:     TS ts
329:     PetscReal t, shift
330:     Vec X, Xdot
331:     Mat J, Jpre
332:     PetscReal user(6)
333:     PetscErrorCode ierr
334:     integer user_a, user_k, user_s
335:     parameter(user_a=0, user_k=2, user_s=4)

337:     DM da
338:     PetscInt mx, xs, xe, gxs, gxe
339:     PetscInt i, i1, row, col
340:     PetscReal k1, k2
341:     PetscScalar val(4)

343:     PetscCall(TSGetDM(ts, da, ierr))
344:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

346:     i1 = 1
347:     k1 = user(user_k + 1)
348:     k2 = user(user_k + 2)
349:     do i = xs, xe
350:       row = i - gxs
351:       col = i - gxs
352:       val(1) = shift + k1
353:       val(2) = -k2
354:       val(3) = -k1
355:       val(4) = shift + k2
356:       PetscCall(MatSetValuesBlockedLocal(Jpre, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
357:     end do
358:     PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr))
359:     PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr))
360:     if (J /= Jpre) then
361:       PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
362:       PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
363:     end if
364:   end subroutine FormIJacobian

366:   subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr)
367:     implicit none
368:     PetscInt mx, xs, xe, gxs, gxe
369:     PetscScalar x(2, xs:xe)
370:     PetscReal a(2), k(2), s(2)
371:     PetscErrorCode ierr

373:     PetscInt i
374:     PetscReal one, hx, r, ik
375:     one = 1.0
376:     hx = one/mx
377:     do i = xs, xe
378:       r = i*hx
379:       if (k(2) /= 0.0) then
380:         ik = one/k(2)
381:       else
382:         ik = one
383:       end if
384:       x(1, i) = one + s(2)*r
385:       x(2, i) = k(1)*ik*x(1, i) + s(2)*ik
386:     end do
387:   end subroutine FormInitialSolutionLocal

389:   subroutine FormInitialSolution(ts, X, user, ierr)
390:     use petscts
391:     use petscdm
392:     implicit none

394:     TS ts
395:     Vec X
396:     PetscReal user(6)
397:     PetscErrorCode ierr
398:     integer user_a, user_k, user_s
399:     parameter(user_a=1, user_k=3, user_s=5)

401:     DM da
402:     PetscInt mx, xs, xe, gxs, gxe
403:     PetscScalar, pointer :: xx(:)

405:     PetscCall(TSGetDM(ts, da, ierr))
406:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

408:     ! Get access to vector data
409:     PetscCall(VecGetArray(X, xx, ierr))

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

413:     PetscCall(VecRestoreArray(X, xx, ierr))
414:   end subroutine FormInitialSolution

416: ! ---------------------------------------------------------------------
417: !
418: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
419: !
420:   subroutine FormIJacobianMF(ts, t, X, Xdot, shift, J, Jpre, user, ierr)
421:     use ex22f_mfmodule
422:     implicit none
423:     TS ts
424:     PetscReal t, shift
425:     Vec X, Xdot
426:     Mat J, Jpre
427:     PetscReal user(6)
428:     PetscErrorCode ierr

430:     PETSC_SHIFT = shift
431:     MFuser = user

433:   end subroutine FormIJacobianMF

435: ! -------------------------------------------------------------------
436: !
437: !   MyMult - user provided matrix multiply
438: !
439: !   Input Parameters:
440: !.  X - input vector
441: !
442: !   Output Parameter:
443: !.  F - function vector
444: !
445:   subroutine MyMult(A, X, F, ierr)
446:     use ex22f_mfmodule
447:     implicit none

449:     Mat A
450:     Vec X, F

452:     PetscErrorCode ierr
453:     PetscScalar shift

455: !  Mat J,Jpre

457:     PetscReal user(6)

459:     integer user_a, user_k, user_s
460:     parameter(user_a=0, user_k=2, user_s=4)

462:     DM da
463:     PetscInt mx, xs, xe, gxs, gxe
464:     PetscInt i, i1, row, col
465:     PetscReal k1, k2
466:     PetscScalar val(4)

468:     shift = PETSC_SHIFT
469:     user = MFuser

471:     PetscCall(TSGetDM(tscontext, da, ierr))
472:     PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr))

474:     i1 = 1
475:     k1 = user(user_k + 1)
476:     k2 = user(user_k + 2)

478:     do i = xs, xe
479:       row = i - gxs
480:       col = i - gxs
481:       val(1) = shift + k1
482:       val(2) = -k2
483:       val(3) = -k1
484:       val(4) = shift + k2
485:       PetscCall(MatSetValuesBlockedLocal(Jmat, i1, [row], i1, [col], val, INSERT_VALUES, ierr))
486:     end do

488: !  PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr))
489: !  PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr))
490: !  if (J /= Jpre) then
491:     PetscCall(MatAssemblyBegin(Jmat, MAT_FINAL_ASSEMBLY, ierr))
492:     PetscCall(MatAssemblyEnd(Jmat, MAT_FINAL_ASSEMBLY, ierr))
493: !  end if

495:     PetscCall(MatMult(Jmat, X, F, ierr))

497:   end subroutine MyMult

499: !
500:   subroutine SaveSolutionToDisk(da, X, gdof, xs, xe)
501:     use petscdm
502:     implicit none

504:     Vec X
505:     DM da
506:     PetscInt xs, xe, two
507:     PetscInt gdof, i
508:     PetscErrorCode ierr
509:     PetscScalar data2(2, xs:xe), data(gdof)
510:     PetscScalar, pointer :: xx(:)

512:     PetscCall(VecGetArrayRead(X, xx, ierr))

514:     two = 2
515:     data2 = reshape(xx(gdof:gdof), (/two, xe - xs + 1/))
516:     data = reshape(data2, (/gdof/))
517:     open (1020, file='solution_out_ex22f_mf.txt')
518:     do i = 1, gdof
519:       write (1020, '(e24.16,1x)') data(i)
520:     end do
521:     close (1020)

523:     PetscCall(VecRestoreArrayRead(X, xx, ierr))
524:   end subroutine SaveSolutionToDisk
525: end program main

527: !/*TEST
528: !
529: !    test:
530: !      args: -da_grid_x 200 -ts_arkimex_type 4
531: !      requires: !single
532: !      output_file: output/empty.out
533: !
534: !TEST*/