Actual source code: ex1f.F

  1: !
  2: ! "$Id: ex1f.F,v 1.39 2001/08/07 03:04:24 balay Exp $";
  3: !
  4: !       Formatted test for TS routines.
  5: !
  6: !          Solves U_t = U_xx
  7: !     F(t,u) = (u_i+1 - 2u_i + u_i-1)/h^2
  8: !       using several different schemes.
  9: !
 10: !23456789012345678901234567890123456789012345678901234567890123456789012

 12:       program main
 13:       implicit none
 14:  #include finclude/petsc.h
 15:  #include finclude/petscvec.h
 16:  #include finclude/petscmat.h
 17:  #include finclude/petscda.h
 18:  #include finclude/petscsys.h
 19:  #include finclude/petscpc.h
 20:  #include finclude/petscsles.h
 21:  #include finclude/petscsnes.h
 22:  #include finclude/petscts.h
 23:  #include finclude/petscdraw.h
 24:  #include finclude/petscviewer.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:       integer          ierr,time_steps,steps,flg,size,problem
 32:       Vec              local,global
 33:       double precision dt,ftime,zero,tmax
 34:       TS               ts
 35:       Mat              A
 36:       MatStructure     A_structure
 37:       TSProblemType    tsproblem
 38:       PetscDraw        draw
 39:       PetscViewer      viewer
 40:       character*(60)   type,tsinfo
 41:       character*(5)    etype

 43: !
 44: !    Application context data, stored in common block
 45: !
 46:       Vec              localwork,csolution
 47:       DA               da
 48:       PetscViewer           viewer1,viewer2
 49:       integer          M
 50:       double precision h,nrm_2,nrm_max,nox
 51:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
 52:      &               ,viewer2,M

 54:       external Monitor,RHSFunctionHeat,RHSMatrixFree,Initial
 55:       external RHSMatrixHeat,RHSJacobianHeat

 57:       zero       = 0.0
 58:       time_steps = 100
 59:       problem    = linear_no_matrix
 60:       A          = 0
 61:       tsproblem  = TS_LINEAR
 62: 
 63:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 64:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

 66:       M = 60
 67:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-M',M,flg,ierr)
 68:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-time',time_steps,   &
 69:      &     flg,ierr)

 71:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-nox',flg,ierr)
 72:       if (flg .eq. 1) then
 73:         nox = 1
 74:       else
 75:         nox = 0
 76:       endif
 77:       nrm_2   = 0.0
 78:       nrm_max = 0.0

 80: !   Set up the ghost point communication pattern

 82:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,M,1,1,            &
 83:      &     PETSC_NULL_INTEGER,da,ierr)
 84:       call DACreateGlobalVector(da,global,ierr)
 85:       call VecGetLocalSize(global,m,ierr)
 86:       call DACreateLocalVector(da,local,ierr)

 88: !   Set up display to show wave graph

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

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

101:       call VecDuplicate(local,localwork,ierr)

103: !   make work array for storing exact solution

105:       call VecDuplicate(global,csolution,ierr)

107:       h = 1.0/(M-1.0)

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

151:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
152:       call TSSetProblemType(ts,tsproblem,ierr)
153:       call TSSetMonitor(ts,Monitor,PETSC_NULL_OBJECT,                   &
154:      &                  PETSC_NULL_FUNCTION, ierr)

156:       dt = h*h/2.01

158:       if (problem .eq. linear_no_matrix) then
159: !
160: !         The user provides the RHS as a Shell matrix.
161: !
162:          call MatCreateShell(PETSC_COMM_WORLD,m,M,M,M,                  &
163:      &        PETSC_NULL_OBJECT,A,ierr)
164:          call MatShellSetOperation(A,MATOP_MULT,RHSMatrixFree,ierr)
165:          call TSSetRHSMatrix(ts,A,A,PETSC_NULL_FUNCTION,                &
166:      &        PETSC_NULL_OBJECT,ierr)
167:       else if (problem .eq. linear_no_time) then
168: !
169: !         The user provides the RHS as a matrix
170: !
171:          call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,M, &
172:      &                  A,ierr)
173:          call MatSetFromOptions(A,ierr)
174:          call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT,  &
175:      &        ierr)
176:          call TSSetRHSMatrix(ts,A,A,PETSC_NULL_FUNCTION,                &
177:      &        PETSC_NULL_OBJECT,ierr)
178:       else if (problem .eq. linear) then
179: !
180: !         The user provides the RHS as a time dependent matrix
181: !
182:          call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,M,M, &
183:      &                  A,ierr)
184:          call MatSetFromOptions(A,ierr)
185:          call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT,  &
186:      &        ierr)
187:          call TSSetRHSMatrix(ts,A,A,RHSMatrixHeat,PETSC_NULL_OBJECT,    &
188:      &        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,PETSC_DECIDE,PETSC_DECIDE,M,M,A &
207:      &                  ,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,60,viewer,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/petscsles.h
278:  #include finclude/petscsnes.h
279:  #include finclude/petscts.h
280:  #include finclude/petscviewer.h
281:       Vec         global
282:       integer     ctx

284:       PetscScalar localptr(1)
285:       integer     i,mybase,myend,ierr
286:       PetscOffset idx

288: !
289: !    Application context data, stored in common block
290: !
291:       Vec              localwork,csolution
292:       DA               da
293:       PetscViewer           viewer1,viewer2
294:       integer          M
295:       double precision h,nrm_2,nrm_max,nox
296:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
297:      &               ,viewer2,M

299: !   determine starting point of each processor

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

303: !   Initialize the array

305:       call VecGetArray(global,localptr,idx,ierr)
306:       do 10, i=mybase,myend-1
307:         localptr(i-mybase+1+idx) = sin(PETSC_PI*i*6.*h) +               &
308:      &                             3.*sin(PETSC_PI*i*2.*h)
309:  10   continue
310:       call VecRestoreArray(global,localptr,idx,ierr)
311:       return
312:       end

314: ! ------------------------------------------------------------------------------
315: !
316: !       Exact solution
317: !
318:       subroutine Solution(t,sol,ctx)
319:       implicit none
320:  #include finclude/petsc.h
321:  #include finclude/petscvec.h
322:  #include finclude/petscmat.h
323:  #include finclude/petscda.h
324:  #include finclude/petscsys.h
325:  #include finclude/petscpc.h
326:  #include finclude/petscsles.h
327:  #include finclude/petscsnes.h
328:  #include finclude/petscts.h
329:  #include finclude/petscviewer.h
330:       double precision  t
331:       Vec               sol
332:       integer           ctx

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

338: !
339: !    Application context data, stored in common block
340: !
341:       Vec               localwork,csolution
342:       DA                da
343:       PetscViewer            viewer1,viewer2
344:       integer           M
345:       double precision  h,nrm_2,nrm_max,nox
346:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
347:      &               ,viewer2,M

349: !   determine starting point of each processor

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

353:       ex1 = exp(-36.*PETSC_PI*PETSC_PI*t)
354:       ex2 = exp(-4.*PETSC_PI*PETSC_PI*t)
355:       sc1 = PETSC_PI*6.*h
356:       sc2 = PETSC_PI*2.*h
357:       call VecGetArray(csolution,localptr,idx,ierr)
358:       do 10, i=mybase,myend-1
359:         localptr(i-mybase+1+idx) = sin(i*sc1)*ex1 + 3.*sin(i*sc2)*ex2
360:  10   continue
361:       call VecRestoreArray(csolution,localptr,idx,ierr)
362:       return
363:       end


366: ! -----------------------------------------------------------------------------------

368:       subroutine Monitor(ts,step,time,global,ctx,ierr)
369:       implicit none
370:  #include finclude/petsc.h
371:  #include finclude/petscvec.h
372:  #include finclude/petscmat.h
373:  #include finclude/petscda.h
374:  #include finclude/petscsys.h
375:  #include finclude/petscpc.h
376:  #include finclude/petscsles.h
377:  #include finclude/petscsnes.h
378:  #include finclude/petscts.h
379:  #include finclude/petscdraw.h
380:  #include finclude/petscviewer.h
381:       TS                ts
382:       integer           step,ctx,ierr
383:       double precision  time,lnrm_2,lnrm_max
384:       Vec               global
385:       PetscScalar       mone

387: !
388: !    Application context data, stored in common block
389: !
390:       Vec              localwork,csolution
391:       DA               da
392:       PetscViewer           viewer1,viewer2
393:       integer          M
394:       double precision h,nrm_2,nrm_max,nox
395:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
396:      &               ,viewer2,M

398:       call VecView(global,viewer1,ierr)

400:       call Solution(time,csolution,ctx)
401:       mone = -1.0
402:       call VecAXPY(mone,global,csolution,ierr)
403:       call VecNorm(csolution,NORM_2,lnrm_2,ierr)
404:       lnrm_2 = sqrt(h)*lnrm_2
405:       call VecNorm(csolution,NORM_MAX,lnrm_max,ierr)

407:       if (nox .eq. 0) then
408:         print*,'timestep ',step,' time ',time,' norm of error ',        &
409:      &         lnrm_2,lnrm_max
410:       endif

412:       nrm_2   = nrm_2 + lnrm_2
413:       nrm_max = nrm_max +  lnrm_max
414:       call VecView(csolution,viewer2,ierr)

416:       return
417:       end

419: !  -----------------------------------------------------------------------

421:       subroutine RHSMatrixFree(mat,x,y,ierr)
422:       implicit none
423:  #include finclude/petsc.h
424:  #include finclude/petscvec.h
425:  #include finclude/petscmat.h
426:  #include finclude/petscda.h
427:  #include finclude/petscsys.h
428:  #include finclude/petscpc.h
429:  #include finclude/petscsles.h
430:  #include finclude/petscsnes.h
431:  #include finclude/petscts.h
432:  #include finclude/petscviewer.h
433:       Mat              mat
434:       Vec              x,y
435:       integer          ierr
436:       double precision zero

438:       zero = 0.0

440:       call RHSFunctionHeat(0,zero,x,y,0,ierr)
441:       return
442:       end

444: ! -------------------------------------------------------------------------

446:       subroutine RHSFunctionHeat(ts,t,globalin,globalout,ctx,ierr)
447:       implicit none
448:  #include finclude/petsc.h
449:  #include finclude/petscvec.h
450:  #include finclude/petscmat.h
451:  #include finclude/petscda.h
452:  #include finclude/petscsys.h
453:  #include finclude/petscpc.h
454:  #include finclude/petscsles.h
455:  #include finclude/petscsnes.h
456:  #include finclude/petscts.h
457:  #include finclude/petscviewer.h
458:       TS               ts
459:       double precision t
460:       Vec              globalin,globalout
461:       integer          ctx
462:       Vec              local
463:       integer          ierr,i,localsize
464:       PetscOffset      ldx,cdx
465:       PetscScalar      copyptr(1),localptr(1),sc
466: !
467: !    Application context data, stored in common block
468: !
469:       Vec              localwork,csolution
470:       DA               da
471:       PetscViewer           viewer1,viewer2
472:       integer          M
473:       double precision h,nrm_2,nrm_max,nox
474:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
475:      &               ,viewer2,M

477: !  Extract local array

479:       call DACreateLocalVector(da,local,ierr)
480:       call DAGlobalToLocalBegin(da,globalin,INSERT_VALUES,local,ierr)
481:       call DAGlobalToLocalEnd(da,globalin,INSERT_VALUES,local,ierr)
482:       call VecGetArray(local,localptr,ldx,ierr)

484: !  Extract work vector

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

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

492:       sc = 1.0/(h*h)
493:       call VecGetLocalSize(local,localsize,ierr)
494:       copyptr(1+cdx) = localptr(1+ldx)
495:       do 10, i=1,localsize-2
496:         copyptr(i+1+cdx) = sc * (localptr(i+2+ldx) + localptr(i+ldx) -  &
497:      &                     2.0*localptr(i+1+ldx))
498:  10   continue
499:       copyptr(localsize-1+1+cdx) = localptr(localsize-1+1+ldx)
500:       call VecRestoreArray(localwork,copyptr,cdx,ierr)
501:       call VecRestoreArray(local,localptr,ldx,ierr)
502:       call VecDestroy(local,ierr)

504: !   Local to Global

506:       call DALocalToGlobal(da,localwork,INSERT_VALUES,globalout,ierr)
507:       return
508:       end

510: !  ---------------------------------------------------------------------

512:       subroutine RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
513:       implicit none
514:  #include finclude/petsc.h
515:  #include finclude/petscvec.h
516:  #include finclude/petscmat.h
517:  #include finclude/petscda.h
518:  #include finclude/petscsys.h
519:  #include finclude/petscpc.h
520:  #include finclude/petscsles.h
521:  #include finclude/petscsnes.h
522:  #include finclude/petscts.h
523:  #include finclude/petscviewer.h
524:       Mat              AA,BB
525:       double precision t
526:       TS               ts
527:       MatStructure     str
528:       integer          ctx,ierr

530:       Mat              A
531:       integer          i,mstart,mend,rank,size,idx(3)
532:       PetscScalar      v(3),stwo,sone

534: !
535: !    Application context data, stored in common block
536: !
537:       Vec              localwork,csolution
538:       DA               da
539:       PetscViewer           viewer1,viewer2
540:       integer          M
541:       double precision h,nrm_2,nrm_max,nox
542:       common /tsctx/ h,nrm_2,nrm_max,nox,localwork,csolution,da,viewer1 &
543:      &               ,viewer2,M

545:       A    = AA
546:       stwo = -2./(h*h)
547:       sone = -.5*stwo
548:       str  = SAME_NONZERO_PATTERN

550:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
551:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

553:       call MatGetOwnershipRange(A,mstart,mend,ierr)
554:       if (mstart .eq. 0) then
555:         v(1) = 1.0
556:         call MatSetValues(A,1,mstart,1,mstart,v,INSERT_VALUES,ierr)
557:         mstart = mstart + 1
558:       endif
559:       if (mend .eq. M) then
560:         mend = mend - 1
561:         v(1) = 1.0
562:         call MatSetValues(A,1,mend,1,mend,v,INSERT_VALUES,ierr)
563:       endif

565: !
566: !     Construct matrice one row at a time
567: !
568:       v(1) = sone
569:       v(2) = stwo
570:       v(3) = sone
571:       do 10, i=mstart,mend-1
572:         idx(1) = i-1
573:         idx(2) = i
574:         idx(3) = i+1
575:         call MatSetValues(A,1,i,3,idx,v,INSERT_VALUES,ierr)
576:  10   continue

578:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
579:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
580:       return
581:       end

583: ! --------------------------------------------------------------------------------------

585:       subroutine RHSJacobianHeat(ts,t,x,AA,BB,str,ctx,ierr)
586:       implicit none
587:  #include finclude/petsc.h
588:  #include finclude/petscvec.h
589:  #include finclude/petscmat.h
590:  #include finclude/petscda.h
591:  #include finclude/petscsys.h
592:  #include finclude/petscpc.h
593:  #include finclude/petscsles.h
594:  #include finclude/petscsnes.h
595:  #include finclude/petscts.h
596:  #include finclude/petscviewer.h
597:       TS               ts
598:       double precision t
599:       Vec              x
600:       Mat              AA,BB
601:       MatStructure     str
602:       integer          ctx,ierr

604:       call RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
605:       return
606:       end