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/petscsys.h
13: #include finclude/petscvec.h
14: #include finclude/petscmat.h
15: #include finclude/petscda.h
16: #include finclude/petscpc.h
17: #include finclude/petscksp.h
18: #include finclude/petscsnes.h
19: #include finclude/petscts.h
20: #include finclude/petscdraw.h
21: #include finclude/petscviewer.h
23: #include "ex1f.h"
25: integer linear_no_matrix,linear_no_time,linear
26: integer nonlinear_no_jacobian,nonlinear
27: parameter (linear_no_matrix = 0,linear_no_time = 1,linear = 2)
28: parameter (nonlinear_no_jacobian = 3,nonlinear = 4)
30: PetscErrorCode ierr
31: PetscInt time_steps,steps
32: PetscMPIInt size
33: integer problem
34: Vec local,global
35: double precision dt,ftime,zero,tmax
36: TS ts
37: Mat A
38: MatStructure A_structure
39: TSProblemType tsproblem
40: PetscDraw draw
41: PetscViewer viewer
42: character*(60) type,tsinfo
43: character*(5) etype
44: PetscInt i1,i60
45: PetscTruth flg
46: PetscSizeT five
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) 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) 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) 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) 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) 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) 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) then
230: !
231: ! copy type to string of length 5 to ensure equality test
232: ! is done correctly
233: !
234: five = 5
235: call PetscStrncpy(etype,type,five,ierr)
236: if (etype .eq. TSEULER) then
237: if (abs(nrm_2/steps - 0.00257244) .gt. 1.e-4) then
238: print*,'Error in Euler method: 2-norm ',nrm_2/steps, &
239: & ' expecting: ',0.00257244
240: endif
241: else
242: if (abs(nrm_2/steps - 0.00506174) .gt. 1.e-4) then
243: print*,'Error in ',tsinfo,': 2-norm ',nrm_2/steps, &
244: & ' expecting: ',0.00506174
245: endif
246: endif
247: else
248: print*,size,' Procs Avg. error 2 norm ',nrm_2/steps, &
249: & nrm_max/steps,tsinfo
250: endif
252: call PetscViewerDestroy(viewer,ierr)
253: call TSDestroy(ts,ierr)
254: call PetscViewerDestroy(viewer1,ierr)
255: call PetscViewerDestroy(viewer2,ierr)
256: call VecDestroy(localwork,ierr)
257: call VecDestroy(csolution,ierr)
258: call VecDestroy(local,ierr)
259: call VecDestroy(global,ierr)
260: call DADestroy(da,ierr)
261: if (A .ne. 0) then
262: call MatDestroy(A,ierr)
263: endif
265: call PetscFinalize(ierr)
266: end
268: ! -------------------------------------------------------------------
269:
270: subroutine Initial(global,ctx,ierr)
271: implicit none
272: #include finclude/petscsys.h
273: #include finclude/petscvec.h
274: #include finclude/petscmat.h
275: #include finclude/petscda.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/petscsys.h
314: #include finclude/petscvec.h
315: #include finclude/petscmat.h
316: #include finclude/petscda.h
317: #include finclude/petscpc.h
318: #include finclude/petscksp.h
319: #include finclude/petscsnes.h
320: #include finclude/petscts.h
321: #include finclude/petscviewer.h
323: #include "ex1f.h"
325: double precision t
326: Vec sol
327: PetscObject ctx
329: PetscScalar localptr(1),ex1
330: PetscScalar ex2,sc1,sc2
331: PetscInt i,mybase,myend
332: PetscErrorCode ierr
333: PetscOffset idx
335: ! determine starting point of each processor
337: call VecGetOwnershipRange(csolution,mybase,myend,ierr)
339: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t)
340: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t)
341: sc1 = PETSC_PI*6.*h
342: sc2 = PETSC_PI*2.*h
343: call VecGetArray(csolution,localptr,idx,ierr)
344: do 10, i=mybase,myend-1
345: localptr(i-mybase+1+idx) = sin(i*sc1)*ex1 + 3.*sin(i*sc2)*ex2
346: 10 continue
347: call VecRestoreArray(csolution,localptr,idx,ierr)
348: return
349: end
352: ! -----------------------------------------------------------------------------------
354: subroutine Monitor(ts,step,time,global,ctx,ierr)
355: implicit none
356: #include finclude/petscsys.h
357: #include finclude/petscvec.h
358: #include finclude/petscmat.h
359: #include finclude/petscda.h
360: #include finclude/petscpc.h
361: #include finclude/petscksp.h
362: #include finclude/petscsnes.h
363: #include finclude/petscts.h
364: #include finclude/petscdraw.h
365: #include finclude/petscviewer.h
367: #include "ex1f.h"
368: TS ts
369: PetscInt step
370: PetscObject ctx
371: PetscErrorCode ierr
372: double precision time,lnrm_2,lnrm_max
373: Vec global
374: PetscScalar mone
376: call VecView(global,viewer1,ierr)
378: call Solution(time,csolution,ctx)
379: mone = -1.0
380: call VecAXPY(csolution,mone,global,ierr)
381: call VecNorm(csolution,NORM_2,lnrm_2,ierr)
382: lnrm_2 = sqrt(h)*lnrm_2
383: call VecNorm(csolution,NORM_MAX,lnrm_max,ierr)
385: if (nox .eq. 0) then
386: print*,'timestep ',step,' time ',time,' norm of error ', &
387: & lnrm_2,lnrm_max
388: endif
390: nrm_2 = nrm_2 + lnrm_2
391: nrm_max = nrm_max + lnrm_max
392: call VecView(csolution,viewer2,ierr)
394: return
395: end
397: ! -----------------------------------------------------------------------
399: subroutine RHSMatrixFree(mat,x,y,ierr)
400: implicit none
401: #include finclude/petscsys.h
402: #include finclude/petscvec.h
403: #include finclude/petscmat.h
404: #include finclude/petscda.h
405: #include finclude/petscpc.h
406: #include finclude/petscksp.h
407: #include finclude/petscsnes.h
408: #include finclude/petscts.h
409: #include finclude/petscviewer.h
410: Mat mat
411: Vec x,y
412: PetscErrorCode ierr
413: double precision zero
414: TS ts0
416: zero = 0.0
418: ts0 = PETSC_NULL_OBJECT
420: call RHSFunctionHeat(ts0,zero,x,y,PETSC_NULL_OBJECT,ierr)
421: return
422: end
424: ! -------------------------------------------------------------------------
426: subroutine RHSFunctionHeat(ts,t,globalin,globalout,ctx,ierr)
427: implicit none
428: #include finclude/petscsys.h
429: #include finclude/petscvec.h
430: #include finclude/petscmat.h
431: #include finclude/petscda.h
432: #include finclude/petscpc.h
433: #include finclude/petscksp.h
434: #include finclude/petscsnes.h
435: #include finclude/petscts.h
436: #include finclude/petscviewer.h
438: #include "ex1f.h"
439: TS ts
440: double precision t
441: Vec globalin,globalout
442: PetscObject ctx
443: Vec local
444: PetscErrorCode ierr
445: PetscInt i,localsize
446: PetscOffset ldx,cdx
447: PetscScalar copyptr(1),localptr(1)
448: PetscScalar sc
450: ! Extract local array
452: call DACreateLocalVector(da,local,ierr)
453: call DAGlobalToLocalBegin(da,globalin,INSERT_VALUES,local,ierr)
454: call DAGlobalToLocalEnd(da,globalin,INSERT_VALUES,local,ierr)
455: call VecGetArray(local,localptr,ldx,ierr)
457: ! Extract work vector
459: call VecGetArray(localwork,copyptr,cdx,ierr)
461: ! Update Locally - Make array of new values
462: ! Note: For the first and last entry I copy the value
463: ! if this is an interior node it is irrelevant
465: sc = 1.0/(h*h)
466: call VecGetLocalSize(local,localsize,ierr)
467: copyptr(1+cdx) = localptr(1+ldx)
468: do 10, i=1,localsize-2
469: copyptr(i+1+cdx) = sc * (localptr(i+2+ldx) + localptr(i+ldx) - &
470: & 2.0*localptr(i+1+ldx))
471: 10 continue
472: copyptr(localsize-1+1+cdx) = localptr(localsize-1+1+ldx)
473: call VecRestoreArray(localwork,copyptr,cdx,ierr)
474: call VecRestoreArray(local,localptr,ldx,ierr)
475: call VecDestroy(local,ierr)
477: ! Local to Global
479: call DALocalToGlobal(da,localwork,INSERT_VALUES,globalout,ierr)
480: return
481: end
483: ! ---------------------------------------------------------------------
485: subroutine RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
486: implicit none
487: #include finclude/petscsys.h
488: #include finclude/petscvec.h
489: #include finclude/petscmat.h
490: #include finclude/petscda.h
491: #include finclude/petscpc.h
492: #include finclude/petscksp.h
493: #include finclude/petscsnes.h
494: #include finclude/petscts.h
495: #include finclude/petscviewer.h
497: #include "ex1f.h"
498: Mat AA,BB
499: double precision t
500: TS ts
501: MatStructure str
502: PetscObject ctx
503: PetscErrorCode ierr
505: Mat A
506: PetscInt i,mstart(1),mend(1),idx(3)
507: PetscMPIInt rank,size
508: PetscScalar v(3),stwo,sone
509: PetscInt i1,i3
511: i1 = 1
512: i3 = 3
513: A = AA
514: stwo = -2./(h*h)
515: sone = -.5*stwo
516: str = SAME_NONZERO_PATTERN
518: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
519: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
521: call MatGetOwnershipRange(A,mstart,mend,ierr)
522: if (mstart(1) .eq. 0) then
523: v(1) = 1.0
524: call MatSetValues(A,i1,mstart(1),i1,mstart(1),v,INSERT_VALUES, &
525: & ierr)
526: mstart(1) = mstart(1) + 1
527: endif
528: if (mend(1) .eq. M) then
529: mend(1) = mend(1) - 1
530: v(1) = 1.0
531: call MatSetValues(A,i1,mend,i1,mend,v,INSERT_VALUES,ierr)
532: endif
534: !
535: ! Construct matrice one row at a time
536: !
537: v(1) = sone
538: v(2) = stwo
539: v(3) = sone
540: do 10, i=mstart(1),mend(1)-1
541: idx(1) = i-1
542: idx(2) = i
543: idx(3) = i+1
544: call MatSetValues(A,i1,i,i3,idx,v,INSERT_VALUES,ierr)
545: 10 continue
547: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
548: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
549: return
550: end
552: ! --------------------------------------------------------------------------------------
554: subroutine RHSJacobianHeat(ts,t,x,AA,BB,str,ctx,ierr)
555: implicit none
556: #include finclude/petscsys.h
557: #include finclude/petscvec.h
558: #include finclude/petscmat.h
559: #include finclude/petscda.h
560: #include finclude/petscpc.h
561: #include finclude/petscksp.h
562: #include finclude/petscsnes.h
563: #include finclude/petscts.h
564: #include finclude/petscviewer.h
565: TS ts
566: double precision t
567: Vec x
568: Mat AA,BB
569: MatStructure str
570: PetscObject ctx
571: PetscErrorCode ierr
573: call RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
574: return
575: end