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

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

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

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

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

 76: !   Set up the ghost point communication pattern

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

 84: !   Set up display to show wave graph

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

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

 97:       call VecDuplicate(local,localwork,ierr)

 99: !   make work array for storing exact solution

101:       call VecDuplicate(global,csolution,ierr)

103:       h = 1.0/(M-1.0)

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

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

152:       dt = h*h/2.01

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

214:       call TSSetFromOptions(ts,ierr)

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

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

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

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

264:       call PetscFinalize(ierr)
265:       end

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

282: #include "ex1f.h"

284:       Vec         global
285:       PetscObject    ctx

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

292: !   determine starting point of each processor

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

296: !   Initialize the array

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

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

324: #include "ex1f.h"

326:       double precision  t
327:       Vec               sol
328:       PetscObject       ctx

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

336: !   determine starting point of each processor

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

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


353: ! -----------------------------------------------------------------------------------

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

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

378:       call VecView(global,viewer1,ierr)

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

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

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

396:       return
397:       end

399: !  -----------------------------------------------------------------------

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

419:       zero = 0.0

421:       ts0 = PETSC_NULL_OBJECT

423:       call RHSFunctionHeat(ts0,zero,x,y,PETSC_NULL_OBJECT,ierr)
424:       return
425:       end

427: ! -------------------------------------------------------------------------

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

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

454: !  Extract local array

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

461: !  Extract work vector

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

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

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

481: !   Local to Global

483:       call DALocalToGlobal(da,localwork,INSERT_VALUES,globalout,ierr)
484:       return
485:       end

487: !  ---------------------------------------------------------------------

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

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

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

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

523:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
524:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

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

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

552:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
553:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
554:       return
555:       end

557: ! --------------------------------------------------------------------------------------

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

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