Actual source code: ex22f.F

petsc-3.4.4 2014-03-13
  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:       implicit none
 19: #include <finclude/petscsys.h>
 20: #include <finclude/petscvec.h>
 21: #include <finclude/petscmat.h>
 22: #include <finclude/petscsnes.h>
 23: #include <finclude/petscts.h>
 24: #include <finclude/petscdmda.h>
 25: !
 26: !  Create an application context to contain data needed by the
 27: !  application-provided call-back routines, FormJacobian() and
 28: !  FormFunction(). We use a double precision array with six
 29: !  entries, two for each problem parameter a, k, s.
 30: !
 31:       PetscReal user(6)
 32:       integer user_a,user_k,user_s
 33:       parameter (user_a = 0,user_k = 2,user_s = 4)

 35:       external FormRHSFunction,FormIFunction,FormIJacobian
 36:       external FormInitialSolution

 38:       TS             ts
 39:       SNES           snes
 40:       SNESLineSearch linesearch
 41:       Vec            X
 42:       Mat            J
 43:       PetscInt       steps,maxsteps,mx
 44:       PetscErrorCode ierr
 45:       DM             da
 46:       PetscReal      ftime,dt
 47:       PetscReal      zero,one,pone
 48:       PetscInt       im11,i2
 49:       PetscBool      flg

 51:       im11 = -11
 52:       i2   = 2
 53:       zero = 0.0
 54:       one = 1.0
 55:       pone = one / 10

 57:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 60: !  Create distributed array (DMDA) to manage parallel grid and vectors
 61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:       call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,im11,i2,i2, &
 63:      &     PETSC_NULL_INTEGER,da,ierr)

 65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66: !    Extract global vectors from DMDA;
 67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68:       call DMCreateGlobalVector(da,X,ierr)

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

 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92: !    Create timestepping solver context
 93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
 95:       call TSSetDM(ts,da,ierr)
 96:       call TSSetType(ts,TSARKIMEX,ierr)
 97:       call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormRHSFunction,       &
 98:      &     user,ierr)
 99:       call TSSetIFunction(ts,PETSC_NULL_OBJECT,FormIFunction,user,ierr)
100:       call DMCreateMatrix(da,MATAIJ,J,ierr)
101:       call TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)

103:       call TSGetSNES(ts,snes,ierr)
104:       call SNESGetLineSearch(snes,linesearch,ierr)
105:       call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)

107:       ftime = 1.0
108:       maxsteps = 10000
109:       call TSSetDuration(ts,maxsteps,ftime,ierr)

111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: !  Set initial conditions
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:       call FormInitialSolution(ts,X,user,ierr)
115:       call TSSetSolution(ts,X,ierr)
116:       call VecGetSize(X,mx,ierr)
117: !  Advective CFL, I don't know why it needs so much safety factor.
118:       dt = pone * max(user(user_a+1),user(user_a+2)) / mx;
119:       call TSSetInitialTimeStep(ts,zero,dt,ierr)

121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: !   Set runtime options
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:       call TSSetFromOptions(ts,ierr)

126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: !  Solve nonlinear system
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:       call TSSolve(ts,X,ierr)

131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: !  Free work space.
133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:       call MatDestroy(J,ierr)
135:       call VecDestroy(X,ierr)
136:       call TSDestroy(ts,ierr)
137:       call DMDestroy(da,ierr)
138:       call PetscFinalize(ierr)
139:       end program

141: ! Small helper to extract the layout, result uses 1-based indexing.
142:       subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr)
143:       implicit none
144: #include <finclude/petscsys.h>
145: #include <finclude/petscdmda.h>
146:       DM da
147:       PetscInt mx,xs,xe,gxs,gxe
148:       PetscErrorCode ierr
149:       PetscInt xm,gxm
150:             call DMDAGetInfo(da,PETSC_NULL_INTEGER,                     &
151:      &     mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                    &
152:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
153:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                       &
154:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,    &
155:      &     PETSC_NULL_INTEGER,ierr)
156:       call DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,  &
157:      &     xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
158:       call DMDAGetGhostCorners(da,                                      &
159:      &     gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                   &
160:      &     gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
161:       xs = xs + 1
162:       gxs = gxs + 1
163:       xe = xs + xm - 1
164:       gxe = gxs + gxm - 1
165:       end subroutine

167:       subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,          &
168:      &     a,k,s,ierr)
169:       implicit none
170:       PetscInt mx,xs,xe,gxs,gxe
171:       PetscScalar x(2,xs:xe)
172:       PetscScalar xdot(2,xs:xe)
173:       PetscScalar f(2,xs:xe)
174:       PetscReal a(2),k(2),s(2)
175:       PetscErrorCode ierr
176:       PetscInt i
177:       do 10 i = xs,xe
178:          f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1)
179:          f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2)
180:  10   continue
181:       end subroutine

183:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
184:       implicit none
185: #include <finclude/petscsys.h>
186: #include <finclude/petscvec.h>
187: #include <finclude/petscmat.h>
188: #include <finclude/petscsnes.h>
189: #include <finclude/petscts.h>
190: #include <finclude/petscdmda.h>
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 VecGetArray(X,xx,ixx,ierr)
209:       call VecGetArray(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 VecRestoreArray(X,xx,ixx,ierr)
217:       call VecRestoreArray(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:       u0t(1) = one - sin(twelve*t)**four
243:       u0t(2) = 0.0
244:       half = one/two
245:       third = one / three
246:       twothird = two / three
247:       sixth = one / six
248:       twelfth = one / twelve
249:       do 20 i = xs,xe
250:          do 10 j = 1,2
251:             if (i .eq. 1) then
252:                f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
253:      &              + sixth*x(j,i+2))
254:             else if (i .eq. 2) then
255:                f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
256:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
257:             else if (i .eq. mx-1) then
258:                f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1)             &
259:      &         - half*x(j,i) -third*x(j,i+1))
260:             else if (i .eq. mx) then
261:                f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
262:             else
263:                f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
264:      &              + twothird*x(j,i-1)                                 &
265:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
266:             end if
267:  10      continue
268:  20   continue
269:       end subroutine

271:       subroutine FormRHSFunction(ts,t,X,F,user,ierr)
272:       implicit none
273: #include <finclude/petscsys.h>
274: #include <finclude/petscvec.h>
275: #include <finclude/petscmat.h>
276: #include <finclude/petscsnes.h>
277: #include <finclude/petscts.h>
278: #include <finclude/petscdmda.h>
279:       TS ts
280:       PetscReal t
281:       Vec X,F
282:       PetscReal user(6)
283:       PetscErrorCode ierr
284:       integer user_a,user_k,user_s
285:       parameter (user_a = 1,user_k = 3,user_s = 5)
286:       DM             da
287:       Vec            Xloc
288:       PetscInt       mx,xs,xe,gxs,gxe
289:       PetscOffset    ixx,iff
290:       PetscScalar    xx(0:1),ff(0:1)

292:       call TSGetDM(ts,da,ierr)
293:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

295: !     Scatter ghost points to local vector,using the 2-step process
296: !        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
297: !     By placing code between these two statements, computations can be
298: !     done while messages are in transition.
299:       call DMGetLocalVector(da,Xloc,ierr)
300:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
301:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)

303: ! Get access to vector data
304:       call VecGetArray(Xloc,xx,ixx,ierr)
305:       call VecGetArray(F,ff,iff,ierr)

307:       call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,                       &
308:      &     t,xx(ixx),ff(iff),                                           &
309:      &     user(user_a),user(user_k),user(user_s),ierr)

311:       call VecRestoreArray(Xloc,xx,ixx,ierr)
312:       call VecRestoreArray(F,ff,iff,ierr)
313:       call DMRestoreLocalVector(da,Xloc,ierr)
314:       end subroutine

316: ! ---------------------------------------------------------------------
317: !
318: !  IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot
319: !
320:       subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,mstr,user,ierr)
321:       implicit none
322: #include <finclude/petscsys.h>
323: #include <finclude/petscvec.h>
324: #include <finclude/petscmat.h>
325: #include <finclude/petscsnes.h>
326: #include <finclude/petscts.h>
327: #include <finclude/petscdmda.h>
328:       TS ts
329:       PetscReal t,shift
330:       Vec X,Xdot
331:       Mat J,Jpre
332:       MatStructure mstr
333:       PetscReal user(6)
334:       PetscErrorCode ierr
335:       integer user_a,user_k,user_s
336:       parameter (user_a = 0,user_k = 2,user_s = 4)

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

344:       call TSGetDM(ts,da,ierr)
345:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

376:       PetscInt i
377:       PetscReal one,hx,r,ik
378:       one = 1.0
379:       hx = one / mx
380:       do 10 i=xs,xe
381:          r = i*hx
382:          if (k(2) .ne. 0.0) then
383:             ik = one/k(2)
384:          else
385:             ik = one
386:          end if
387:          x(1,i) = one + s(2)*r
388:          x(2,i) = k(1)*ik*x(1,i) + s(2)*ik
389:  10   continue
390:       end subroutine

392:       subroutine FormInitialSolution(ts,X,user,ierr)
393:       implicit none
394: #include <finclude/petscsys.h>
395: #include <finclude/petscvec.h>
396: #include <finclude/petscmat.h>
397: #include <finclude/petscsnes.h>
398: #include <finclude/petscts.h>
399: #include <finclude/petscdmda.h>
400:       TS ts
401:       Vec X
402:       PetscReal user(6)
403:       PetscErrorCode ierr
404:       integer user_a,user_k,user_s
405:       parameter (user_a = 1,user_k = 3,user_s = 5)

407:       DM             da
408:       PetscInt       mx,xs,xe,gxs,gxe
409:       PetscOffset    ixx
410:       PetscScalar    xx(0:1)

412:       call TSGetDM(ts,da,ierr)
413:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

421:       call VecRestoreArray(X,xx,ixx,ierr)
422:       end subroutine