Actual source code: ex22f.F

petsc-3.3-p7 2013-05-11
  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:       Vec            X
 40:       Mat            J
 41:       PetscInt       steps,maxsteps,mx
 42:       PetscErrorCode ierr
 43:       DM             da
 44:       PetscReal      ftime,dt
 45:       PetscReal      zero,one,pone
 46:       PetscInt       im11,i2
 47:       PetscBool      flg

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

 55:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

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

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

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

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

101:       ftime = 1.0
102:       maxsteps = 10000
103:       call TSSetDuration(ts,maxsteps,ftime,ierr)

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

115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: !   Set runtime options
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:       call TSSetFromOptions(ts,ierr)

120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: !  Solve nonlinear system
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:       call TSSolve(ts,X,ftime,ierr)
124:       call TSGetTimeStepNumber(ts,steps,ierr)

126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: !  Free work space.
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:       call MatDestroy(J,ierr)
130:       call VecDestroy(X,ierr)
131:       call TSDestroy(ts,ierr)
132:       call DMDestroy(da,ierr)
133:       call PetscFinalize(ierr)
134:       end program

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

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

178:       subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr)
179:       implicit none
180: #include <finclude/petscsys.h>
181: #include <finclude/petscvec.h>
182: #include <finclude/petscmat.h>
183: #include <finclude/petscsnes.h>
184: #include <finclude/petscts.h>
185: #include <finclude/petscdmda.h>
186:       TS ts
187:       PetscReal t
188:       Vec X,Xdot,F
189:       PetscReal user(6)
190:       PetscErrorCode ierr
191:       integer user_a,user_k,user_s
192:       parameter (user_a = 1,user_k = 3,user_s = 5)

194:       DM             da
195:       PetscInt       mx,xs,xe,gxs,gxe
196:       PetscOffset    ixx,ixxdot,iff
197:       PetscScalar    xx(0:1),xxdot(0:1),ff(0:1)

199:       call TSGetDM(ts,da,ierr)
200:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

202: ! Get access to vector data
203:       call VecGetArray(X,xx,ixx,ierr)
204:       call VecGetArray(Xdot,xxdot,ixxdot,ierr)
205:       call VecGetArray(F,ff,iff,ierr)

207:       call FormIFunctionLocal(mx,xs,xe,gxs,gxe,                         &
208:      &     xx(ixx),xxdot(ixxdot),ff(iff),                               &
209:      &     user(user_a),user(user_k),user(user_s),ierr)

211:       call VecRestoreArray(X,xx,ixx,ierr)
212:       call VecRestoreArray(Xdot,xxdot,ixxdot,ierr)
213:       call VecRestoreArray(F,ff,iff,ierr)
214:       end subroutine

216:       subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,           &
217:      &     a,k,s,ierr)
218:       implicit none
219:       PetscInt mx,xs,xe,gxs,gxe
220:       PetscReal t
221:       PetscScalar x(2,gxs:gxe),f(2,xs:xe)
222:       PetscReal a(2),k(2),s(2)
223:       PetscErrorCode ierr
224:       PetscInt i,j
225:       PetscReal hx,u0t(2)
226:       PetscReal one,two,three,four,six,twelve
227:       PetscReal half,third,twothird,sixth
228:       PetscReal twelfth

230:       one = 1.0
231:       two = 2.0
232:       three = 3.0
233:       four = 4.0
234:       six = 6.0
235:       twelve = 12.0
236:       hx = one / mx
237:       u0t(1) = one - sin(twelve*t)**four
238:       u0t(2) = 0.0
239:       half = one/two
240:       third = one / three
241:       twothird = two / three
242:       sixth = one / six
243:       twelfth = one / twelve
244:       do 20 i = xs,xe
245:          do 10 j = 1,2
246:             if (i .eq. 1) then
247:                f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1)  &
248:      &              + sixth*x(j,i+2))
249:             else if (i .eq. 2) then
250:                f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1)    &
251:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
252:             else if (i .eq. mx-1) then
253:                f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1)             &
254:      &         - half*x(j,i) -third*x(j,i+1))
255:             else if (i .eq. mx) then
256:                f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1))
257:             else
258:                f(j,i) = a(j)/hx*(-twelfth*x(j,i-2)                      &
259:      &              + twothird*x(j,i-1)                                 &
260:      &              - twothird*x(j,i+1) + twelfth*x(j,i+2))
261:             end if
262:  10      continue
263:  20   continue
264:       end subroutine

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

287:       call TSGetDM(ts,da,ierr)
288:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

298: ! Get access to vector data
299:       call VecGetArray(Xloc,xx,ixx,ierr)
300:       call VecGetArray(F,ff,iff,ierr)

302:       call FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,                       &
303:      &     t,xx(ixx),ff(iff),                                           &
304:      &     user(user_a),user(user_k),user(user_s),ierr)

306:       call VecRestoreArray(Xloc,xx,ixx,ierr)
307:       call VecRestoreArray(F,ff,iff,ierr)
308:       call DMRestoreLocalVector(da,Xloc,ierr)
309:       end subroutine

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

333:       DM             da
334:       PetscInt       mx,xs,xe,gxs,gxe
335:       PetscInt       i,i1,row,col
336:       PetscReal      k1,k2;
337:       PetscScalar    val(4)

339:       call TSGetDM(ts,da,ierr)
340:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

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

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

402:       DM             da
403:       PetscInt       mx,xs,xe,gxs,gxe
404:       PetscOffset    ixx
405:       PetscScalar    xx(0:1)

407:       call TSGetDM(ts,da,ierr)
408:       call GetLayout(da,mx,xs,xe,gxs,gxe,ierr)

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

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

416:       call VecRestoreArray(X,xx,ixx,ierr)
417:       end subroutine