Actual source code: ex22f.F

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

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

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

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

 58:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

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

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

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

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

105:       call TSGetSNES(ts,snes,ierr)
106:       call SNESGetLineSearch(snes,linesearch,ierr)
107:       call SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)

109:       ftime = 1.0
110:       maxsteps = 10000
111:       call TSSetDuration(ts,maxsteps,ftime,ierr)

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

123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: !   Set runtime options
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:       call TSSetFromOptions(ts,ierr)

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

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

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

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

186:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
187:       implicit none
188: #include <petsc/finclude/petscsys.h>
189: #include <petsc/finclude/petscvec.h>
190: #include <petsc/finclude/petscmat.h>
191: #include <petsc/finclude/petscsnes.h>
192: #include <petsc/finclude/petscts.h>
193: #include <petsc/finclude/petscdm.h>
194: #include <petsc/finclude/petscdmda.h>
195:       TS ts
196:       PetscReal t
197:       Vec X,Xdot,F
198:       PetscReal user(6)
199:       PetscErrorCode ierr
200:       integer user_a,user_k,user_s
201:       parameter (user_a = 1,user_k = 3,user_s = 5)

203:       DM             da
204:       PetscInt       mx,xs,xe,gxs,gxe
205:       PetscOffset    ixx,ixxdot,iff
206:       PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

208:       call TSGetDM(ts,da,ierr)
209:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

211: ! Get access to vector data
212:       call VecGetArrayRead(X,xx,ixx,ierr)
213:       call VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)
214:       call VecGetArray(F,ff,iff,ierr)

216:       call FormIFunctionLocal(mx,xs,xe,gxs,gxe,                         &
217:      &     xx(ixx),xxdot(ixxdot),ff(iff),                               &
218:      &     user(user_a),user(user_k),user(user_s),ierr)

220:       call VecRestoreArrayRead(X,xx,ixx,ierr)
221:       call VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)
222:       call VecRestoreArray(F,ff,iff,ierr)
223:       end subroutine

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

239:       one = 1.0
240:       two = 2.0
241:       three = 3.0
242:       four = 4.0
243:       six = 6.0
244:       twelve = 12.0
245:       hx = one / mx
246:       u0t(1) = one - sin(twelve*t)**four
247:       u0t(2) = 0.0
248:       half = one/two
249:       third = one / three
250:       twothird = two / three
251:       sixth = one / six
252:       twelfth = one / twelve
253:       do 20 i = xs,xe
254:          do 10 j = 1,2
255:             if (i .eq. 1) then
256:                f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
257:      &              + sixth*x(j,i+2))
258:             else if (i .eq. 2) then
259:                f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
260:      &              - 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)             &
263:      &         - half*x(j,i) -third*x(j,i+1))
264:             else if (i .eq. mx) then
265:                f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
266:             else
267:                f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
268:      &              + twothird*x(j,i-1)                                 &
269:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
270:             end if
271:  10      continue
272:  20   continue
273:       end subroutine

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

297:       call TSGetDM(ts,da,ierr)
298:       call 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:       call DMGetLocalVector(da,Xloc,ierr)
305:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)
306:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)

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

312:       call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,                       &
313:      &     t,xx(ixx),ff(iff),                                           &
314:      &     user(user_a),user(user_k),user(user_s),ierr)

316:       call VecRestoreArrayRead(Xloc,xx,ixx,ierr)
317:       call VecRestoreArray(F,ff,iff,ierr)
318:       call DMRestoreLocalVector(da,Xloc,ierr)
319:       end subroutine

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

343:       DM             da
344:       PetscInt       mx,xs,xe,gxs,gxe
345:       PetscInt       i,i1,row,col
346:       PetscReal      k1,k2;
347:       PetscScalar    val(4)

349:       call TSGetDM(ts,da,ierr)
350:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

373:       subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr)
374:       implicit none
375:       PetscInt mx,xs,xe,gxs,gxe
376:       PetscScalar x(2,xs:xe)
377:       PetscReal a(2),k(2),s(2)
378:       PetscErrorCode ierr

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

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

412:       DM             da
413:       PetscInt       mx,xs,xe,gxs,gxe
414:       PetscOffset    ixx
415:       PetscScalar    xx(0:1)

417:       call TSGetDM(ts,da,ierr)
418:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

426:       call VecRestoreArray(X,xx,ixx,ierr)
427:       end subroutine