Actual source code: ex1f.F

  1: !
  2: !       Formatted test for TS routines.
  3: !
  4: !          Solves U_t = U_xx
  5: !     F(t,u) = (u_i+1 - 2u_i + u_i-1)/h^2
  6: !       using several different schemes.
  7: !
  8: !23456789012345678901234567890123456789012345678901234567890123456789012

 10:       program main
 11:       implicit none
 12:  #include finclude/petsc.h
 13:  #include finclude/petscvec.h
 14:  #include finclude/petscmat.h
 15:  #include finclude/petscda.h
 16:  #include finclude/petscsys.h
 17:  #include finclude/petscpc.h
 18:  #include finclude/petscksp.h
 19:  #include finclude/petscsnes.h
 20:  #include finclude/petscts.h
 21:  #include finclude/petscdraw.h
 22:  #include finclude/petscviewer.h

 24: #include "ex1f.h"

 26:       integer   linear_no_matrix,linear_no_time,linear
 27:       integer   nonlinear_no_jacobian,nonlinear
 28:       parameter (linear_no_matrix = 0,linear_no_time = 1,linear = 2)
 29:       parameter (nonlinear_no_jacobian = 3,nonlinear = 4)

 31:       PetscErrorCode  ierr
 32:       PetscInt time_steps,steps
 33:       PetscMPIInt size
 34:       integer problem
 35:       Vec              local,global
 36:       double precision dt,ftime,zero,tmax
 37:       TS               ts
 38:       Mat              A
 39:       MatStructure     A_structure
 40:       TSProblemType    tsproblem
 41:       PetscDraw        draw
 42:       PetscViewer      viewer
 43:       character*(60)   type,tsinfo
 44:       character*(5)    etype
 45:       PetscInt         i1,i60
 46:       PetscTruth flg
 47:       PetscSizeT five

 49:       external Monitor,RHSFunctionHeat,RHSMatrixFree,Initial
 50:       external RHSMatrixHeat,RHSJacobianHeat

 52:       i1         = 1
 53:       i60        = 60
 54:       zero       = 0.0
 55:       time_steps = 100
 56:       problem    = linear_no_matrix
 57:       A          = 0
 58:       tsproblem  = TS_LINEAR
 59: 
 60:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 61:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

 63:       M = 60
 64:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-M',M,flg,ierr)
 65:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-time',time_steps,   &
 66:      &     flg,ierr)

 68:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-nox',flg,ierr)
 69:       if (flg) then
 70:         nox = 1
 71:       else
 72:         nox = 0
 73:       endif
 74:       nrm_2   = 0.0
 75:       nrm_max = 0.0

 77: !   Set up the ghost point communication pattern

 79:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,M,i1,i1,            &
 80:      &     PETSC_NULL_INTEGER,da,ierr)
 81:       call DACreateGlobalVector(da,global,ierr)
 82:       call VecGetLocalSize(global,m,ierr)
 83:       call DACreateLocalVector(da,local,ierr)

 85: !   Set up display to show wave graph

 87:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,        &
 88:      &     PETSC_NULL_CHARACTER,80,380,400,160,viewer1,ierr)
 89:       call PetscViewerDrawGetDraw(viewer1,0,draw,ierr)
 90:       call PetscDrawSetDoubleBuffer(draw,ierr)
 91:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,        &
 92:      &     PETSC_NULL_CHARACTER,80,0,400,160,viewer2,ierr)
 93:       call PetscViewerDrawGetDraw(viewer2,0,draw,ierr)
 94:       call PetscDrawSetDoubleBuffer(draw,ierr)

 96: !   make work array for evaluating right hand side function

 98:       call VecDuplicate(local,localwork,ierr)

100: !   make work array for storing exact solution

102:       call VecDuplicate(global,csolution,ierr)

104:       h = 1.0/(M-1.0)

106: !   set initial conditions
107: 
108:       call Initial(global,PETSC_NULL_OBJECT,ierr)
109: 
110: !
111: !     This example is written to allow one to easily test parts
112: !    of TS, we do not expect users to generally need to use more
113: !    then a single TSProblemType
114: !
115:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-linear_no_matrix',&
116:      &                    flg,ierr)
117:       if (flg) then
118:         tsproblem = TS_LINEAR
119:         problem   = linear_no_matrix
120:       endif
121:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,                    &
122:      &     '-linear_constant_matrix',flg,ierr)
123:       if (flg) then
124:         tsproblem = TS_LINEAR
125:         problem   = linear_no_time
126:       endif
127:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,                    &
128:      &     '-linear_variable_matrix',flg,ierr)
129:       if (flg) then
130:         tsproblem = TS_LINEAR
131:         problem   = linear
132:       endif
133:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,                    &
134:      &     '-nonlinear_no_jacobian',flg,ierr)
135:       if (flg) then
136:         tsproblem = TS_NONLINEAR
137:         problem   = nonlinear_no_jacobian
138:       endif
139:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,                    &
140:      &     '-nonlinear_jacobian',flg,ierr)
141:       if (flg) then
142:         tsproblem = TS_NONLINEAR
143:         problem   = nonlinear
144:       endif
145: 
146: !   make timestep context

148:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
149:       call TSSetProblemType(ts,tsproblem,ierr)
150:       call TSMonitorSet(ts,Monitor,PETSC_NULL_OBJECT,                   &
151:      &                  PETSC_NULL_FUNCTION, ierr)

153:       dt = h*h/2.01

155:       if (problem .eq. linear_no_matrix) then
156: !
157: !         The user provides the RHS as a Shell matrix.
158: !
159:          call MatCreateShell(PETSC_COMM_WORLD,m,M,M,M,                  &
160:      &        PETSC_NULL_OBJECT,A,ierr)
161:          call MatShellSetOperation(A,MATOP_MULT,RHSMatrixFree,ierr)
162:          call TSSetMatrices(ts,A,PETSC_NULL_FUNCTION,                    &
163:      &           PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,                  &
164:      &           DIFFERENT_NONZERO_PATTERN,PETSC_NULL_OBJECT,ierr)
165:       else if (problem .eq. linear_no_time) then
166: !
167: !         The user provides the RHS as a matrix
168: !
169:          call MatCreate(PETSC_COMM_WORLD,A,ierr)
170:          call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M,ierr)
171:          call MatSetFromOptions(A,ierr)
172:          call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT,  &
173:      &        ierr)
174:          call TSSetMatrices(ts,A,PETSC_NULL_FUNCTION,                    &
175:      &           PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,                  &
176:      &           DIFFERENT_NONZERO_PATTERN,PETSC_NULL_OBJECT,ierr)
177:       else if (problem .eq. linear) then
178: !
179: !         The user provides the RHS as a time dependent matrix
180: !
181:          call MatCreate(PETSC_COMM_WORLD,A,ierr)
182:          call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M,ierr)
183:          call MatSetFromOptions(A,ierr)
184:          call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT,  &
185:      &        ierr)
186:          call TSSetMatrices(ts,A,RHSMatrixHeat,                          &
187:      &           PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,                  &
188:      &           DIFFERENT_NONZERO_PATTERN,PETSC_NULL_OBJECT,ierr)
189:       else if (problem .eq. nonlinear_no_jacobian) then
190: !
191: !         The user provides the RHS and a Shell Jacobian
192: !
193:          call TSSetRHSFunction(ts,RHSFunctionHeat,PETSC_NULL_OBJECT,    &
194:      &        ierr)
195:          call MatCreateShell(PETSC_COMM_WORLD,m,M,M,M,                  &
196:      &        PETSC_NULL_OBJECT,A,ierr)
197:          call MatShellSetOperation(A,MATOP_MULT,RHSMatrixFree,ierr)
198:          call TSSetRHSJacobian(ts,A,A,PETSC_NULL_FUNCTION,               &
199:      &        PETSC_NULL_OBJECT,ierr)
200:       else if (problem .eq. nonlinear) then
201: !
202: !         The user provides the RHS and Jacobian
203: !
204:          call TSSetRHSFunction(ts,RHSFunctionHeat,PETSC_NULL_OBJECT,     &
205:      &                         ierr)
206:          call MatCreate(PETSC_COMM_WORLD,A,ierr)
207:          call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M,ierr)
208:          call MatSetFromOptions(A,ierr)
209:          call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT,   &
210:      &        ierr)
211:          call TSSetRHSJacobian(ts,A,A,RHSJacobianHeat,                   &
212:      &PETSC_NULL_OBJECT,ierr)
213:       endif

215:       call TSSetFromOptions(ts,ierr)

217:       call TSSetInitialTimeStep(ts,zero,dt,ierr)
218:       tmax = 100.0
219:       call TSSetDuration(ts,time_steps,tmax,ierr)
220:       call TSSetSolution(ts,global,ierr)

222:       call TSSetUp(ts,ierr)
223:       call TSStep(ts,steps,ftime,ierr)
224:       call PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,i60,viewer,         &
225:      &                           ierr)
226:       call TSView(ts,viewer,ierr)
227:       call TSGetType(ts,type,ierr)

229:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-test',flg,ierr)
230:       if (flg) then
231: !
232: !         copy type to string of length 5 to ensure equality test
233: !         is done correctly
234: !
235:         five = 5
236:         call PetscStrncpy(etype,type,five,ierr)
237:         if (etype .eq. TS_EULER) then
238:           if (abs(nrm_2/steps - 0.00257244) .gt. 1.e-4) then
239:             print*,'Error in Euler method: 2-norm ',nrm_2/steps,        &
240:      &            ' expecting: ',0.00257244
241:           endif
242:         else
243:           if (abs(nrm_2/steps - 0.00506174) .gt. 1.e-4) then
244:             print*,'Error in ',tsinfo,': 2-norm ',nrm_2/steps,          &
245:      &             ' expecting: ',0.00506174
246:           endif
247:         endif
248:       else
249:         print*,size,' Procs Avg. error 2 norm ',nrm_2/steps,            &
250:      &              nrm_max/steps,tsinfo
251:       endif

253:       call PetscViewerDestroy(viewer,ierr)
254:       call TSDestroy(ts,ierr)
255:       call PetscViewerDestroy(viewer1,ierr)
256:       call PetscViewerDestroy(viewer2,ierr)
257:       call VecDestroy(localwork,ierr)
258:       call VecDestroy(csolution,ierr)
259:       call VecDestroy(local,ierr)
260:       call VecDestroy(global,ierr)
261:       call DADestroy(da,ierr)
262:       if (A .ne. 0) then
263:         call MatDestroy(A,ierr)
264:       endif

266:       call PetscFinalize(ierr)
267:       end

269: !  -------------------------------------------------------------------
270: 
271:       subroutine Initial(global,ctx,ierr)
272:       implicit none
273:  #include finclude/petsc.h
274:  #include finclude/petscvec.h
275:  #include finclude/petscmat.h
276:  #include finclude/petscda.h
277:  #include finclude/petscsys.h
278:  #include finclude/petscpc.h
279:  #include finclude/petscksp.h
280:  #include finclude/petscsnes.h
281:  #include finclude/petscts.h
282:  #include finclude/petscviewer.h

284: #include "ex1f.h"

286:       Vec         global
287:       PetscObject    ctx

289:       PetscScalar localptr(1)
290:       PetscInt     i,mybase,myend
291:       PetscErrorCode ierr
292:       PetscOffset idx

294: !   determine starting point of each processor

296:       call VecGetOwnershipRange(global,mybase,myend,ierr)

298: !   Initialize the array

300:       call VecGetArray(global,localptr,idx,ierr)
301:       do 10, i=mybase,myend-1
302:         localptr(i-mybase+1+idx) = sin(PETSC_PI*i*6.*h) +               &
303:      &                             3.*sin(PETSC_PI*i*2.*h)
304:  10   continue
305:       call VecRestoreArray(global,localptr,idx,ierr)
306:       return
307:       end

309: ! ------------------------------------------------------------------------------
310: !
311: !       Exact solution
312: !
313:       subroutine Solution(t,sol,ctx)
314:       implicit none
315:  #include finclude/petsc.h
316:  #include finclude/petscvec.h
317:  #include finclude/petscmat.h
318:  #include finclude/petscda.h
319:  #include finclude/petscsys.h
320:  #include finclude/petscpc.h
321:  #include finclude/petscksp.h
322:  #include finclude/petscsnes.h
323:  #include finclude/petscts.h
324:  #include finclude/petscviewer.h

326: #include "ex1f.h"

328:       double precision  t
329:       Vec               sol
330:       PetscObject       ctx

332:       PetscScalar localptr(1),ex1
333:       PetscScalar ex2,sc1,sc2
334:       PetscInt    i,mybase,myend
335:       PetscErrorCode ierr
336:       PetscOffset       idx

338: !   determine starting point of each processor

340:       call VecGetOwnershipRange(csolution,mybase,myend,ierr)

342:       ex1 = exp(-36.*PETSC_PI*PETSC_PI*t)
343:       ex2 = exp(-4.*PETSC_PI*PETSC_PI*t)
344:       sc1 = PETSC_PI*6.*h
345:       sc2 = PETSC_PI*2.*h
346:       call VecGetArray(csolution,localptr,idx,ierr)
347:       do 10, i=mybase,myend-1
348:         localptr(i-mybase+1+idx) = sin(i*sc1)*ex1 + 3.*sin(i*sc2)*ex2
349:  10   continue
350:       call VecRestoreArray(csolution,localptr,idx,ierr)
351:       return
352:       end


355: ! -----------------------------------------------------------------------------------

357:       subroutine Monitor(ts,step,time,global,ctx,ierr)
358:       implicit none
359:  #include finclude/petsc.h
360:  #include finclude/petscvec.h
361:  #include finclude/petscmat.h
362:  #include finclude/petscda.h
363:  #include finclude/petscsys.h
364:  #include finclude/petscpc.h
365:  #include finclude/petscksp.h
366:  #include finclude/petscsnes.h
367:  #include finclude/petscts.h
368:  #include finclude/petscdraw.h
369:  #include finclude/petscviewer.h

371: #include "ex1f.h"
372:       TS                ts
373:       PetscInt           step
374:       PetscObject     ctx
375:       PetscErrorCode ierr
376:       double precision  time,lnrm_2,lnrm_max
377:       Vec               global
378:       PetscScalar       mone

380:       call VecView(global,viewer1,ierr)

382:       call Solution(time,csolution,ctx)
383:       mone = -1.0
384:       call VecAXPY(csolution,mone,global,ierr)
385:       call VecNorm(csolution,NORM_2,lnrm_2,ierr)
386:       lnrm_2 = sqrt(h)*lnrm_2
387:       call VecNorm(csolution,NORM_MAX,lnrm_max,ierr)

389:       if (nox .eq. 0) then
390:         print*,'timestep ',step,' time ',time,' norm of error ',        &
391:      &         lnrm_2,lnrm_max
392:       endif

394:       nrm_2   = nrm_2 + lnrm_2
395:       nrm_max = nrm_max +  lnrm_max
396:       call VecView(csolution,viewer2,ierr)

398:       return
399:       end

401: !  -----------------------------------------------------------------------

403:       subroutine RHSMatrixFree(mat,x,y,ierr)
404:       implicit none
405:  #include finclude/petsc.h
406:  #include finclude/petscvec.h
407:  #include finclude/petscmat.h
408:  #include finclude/petscda.h
409:  #include finclude/petscsys.h
410:  #include finclude/petscpc.h
411:  #include finclude/petscksp.h
412:  #include finclude/petscsnes.h
413:  #include finclude/petscts.h
414:  #include finclude/petscviewer.h
415:       Mat              mat
416:       Vec              x,y
417:       PetscErrorCode          ierr
418:       double precision zero
419:       TS               ts0

421:       zero = 0.0

423:       ts0 = PETSC_NULL_OBJECT

425:       call RHSFunctionHeat(ts0,zero,x,y,PETSC_NULL_OBJECT,ierr)
426:       return
427:       end

429: ! -------------------------------------------------------------------------

431:       subroutine RHSFunctionHeat(ts,t,globalin,globalout,ctx,ierr)
432:       implicit none
433:  #include finclude/petsc.h
434:  #include finclude/petscvec.h
435:  #include finclude/petscmat.h
436:  #include finclude/petscda.h
437:  #include finclude/petscsys.h
438:  #include finclude/petscpc.h
439:  #include finclude/petscksp.h
440:  #include finclude/petscsnes.h
441:  #include finclude/petscts.h
442:  #include finclude/petscviewer.h

444: #include "ex1f.h"
445:       TS               ts
446:       double precision t
447:       Vec globalin,globalout
448:       PetscObject ctx
449:       Vec local
450:       PetscErrorCode  ierr
451:       PetscInt i,localsize
452:       PetscOffset ldx,cdx
453:       PetscScalar copyptr(1),localptr(1)
454:       PetscScalar sc

456: !  Extract local array

458:       call DACreateLocalVector(da,local,ierr)
459:       call DAGlobalToLocalBegin(da,globalin,INSERT_VALUES,local,ierr)
460:       call DAGlobalToLocalEnd(da,globalin,INSERT_VALUES,local,ierr)
461:       call VecGetArray(local,localptr,ldx,ierr)

463: !  Extract work vector

465:       call VecGetArray(localwork,copyptr,cdx,ierr)

467: !   Update Locally - Make array of new values
468: !   Note: For the first and last entry I copy the value
469: !   if this is an interior node it is irrelevant

471:       sc = 1.0/(h*h)
472:       call VecGetLocalSize(local,localsize,ierr)
473:       copyptr(1+cdx) = localptr(1+ldx)
474:       do 10, i=1,localsize-2
475:         copyptr(i+1+cdx) = sc * (localptr(i+2+ldx) + localptr(i+ldx) -  &
476:      &                     2.0*localptr(i+1+ldx))
477:  10   continue
478:       copyptr(localsize-1+1+cdx) = localptr(localsize-1+1+ldx)
479:       call VecRestoreArray(localwork,copyptr,cdx,ierr)
480:       call VecRestoreArray(local,localptr,ldx,ierr)
481:       call VecDestroy(local,ierr)

483: !   Local to Global

485:       call DALocalToGlobal(da,localwork,INSERT_VALUES,globalout,ierr)
486:       return
487:       end

489: !  ---------------------------------------------------------------------

491:       subroutine RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
492:       implicit none
493:  #include finclude/petsc.h
494:  #include finclude/petscvec.h
495:  #include finclude/petscmat.h
496:  #include finclude/petscda.h
497:  #include finclude/petscsys.h
498:  #include finclude/petscpc.h
499:  #include finclude/petscksp.h
500:  #include finclude/petscsnes.h
501:  #include finclude/petscts.h
502:  #include finclude/petscviewer.h

504: #include "ex1f.h"
505:       Mat              AA,BB
506:       double precision t
507:       TS               ts
508:       MatStructure     str
509:       PetscObject          ctx
510:       PetscErrorCode ierr

512:       Mat              A
513:       PetscInt    i,mstart(1),mend(1),idx(3)
514:       PetscMPIInt rank,size
515:       PetscScalar v(3),stwo,sone
516:       PetscInt    i1,i3

518:       i1 = 1
519:       i3 = 3
520:       A    = AA
521:       stwo = -2./(h*h)
522:       sone = -.5*stwo
523:       str  = SAME_NONZERO_PATTERN

525:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
526:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

528:       call MatGetOwnershipRange(A,mstart,mend,ierr)
529:       if (mstart(1) .eq. 0) then
530:         v(1) = 1.0
531:         call MatSetValues(A,i1,mstart(1),i1,mstart(1),v,INSERT_VALUES,  &
532:      &       ierr)
533:         mstart(1) = mstart(1) + 1
534:       endif
535:       if (mend(1) .eq. M) then
536:         mend(1) = mend(1) - 1
537:         v(1) = 1.0
538:         call MatSetValues(A,i1,mend,i1,mend,v,INSERT_VALUES,ierr)
539:       endif

541: !
542: !     Construct matrice one row at a time
543: !
544:       v(1) = sone
545:       v(2) = stwo
546:       v(3) = sone
547:       do 10, i=mstart(1),mend(1)-1
548:         idx(1) = i-1
549:         idx(2) = i
550:         idx(3) = i+1
551:         call MatSetValues(A,i1,i,i3,idx,v,INSERT_VALUES,ierr)
552:  10   continue

554:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
555:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
556:       return
557:       end

559: ! --------------------------------------------------------------------------------------

561:       subroutine RHSJacobianHeat(ts,t,x,AA,BB,str,ctx,ierr)
562:       implicit none
563:  #include finclude/petsc.h
564:  #include finclude/petscvec.h
565:  #include finclude/petscmat.h
566:  #include finclude/petscda.h
567:  #include finclude/petscsys.h
568:  #include finclude/petscpc.h
569:  #include finclude/petscksp.h
570:  #include finclude/petscsnes.h
571:  #include finclude/petscts.h
572:  #include finclude/petscviewer.h
573:       TS               ts
574:       double precision t
575:       Vec              x
576:       Mat              AA,BB
577:       MatStructure     str
578:       PetscObject      ctx
579:       PetscErrorCode ierr

581:       call RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
582:       return
583:       end