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/petscksp.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/petscksp.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/petscksp.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/petscksp.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/petscksp.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/petscksp.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/petscksp.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/petscksp.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