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 .eq. 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 .eq. 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 .eq. 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 .eq. 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) .ne. 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*/