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: ! This example need rewritten using newly developed TS interface!
12: program main
13: implicit none
14: #include <finclude/petscsys.h>
15: #include <finclude/petscvec.h>
16: #include <finclude/petscmat.h>
17: #include <finclude/petscdmda.h>
18: #include <finclude/petscpc.h>
19: #include <finclude/petscksp.h>
20: #include <finclude/petscsnes.h>
21: #include <finclude/petscts.h>
22: #include <finclude/petscdraw.h>
23: #include <finclude/petscviewer.h>
25: #include "ex1f.h"
27: integer linear_no_matrix,linear_no_time,linear
28: integer nonlinear_no_jacobian,nonlinear
29: parameter (linear_no_matrix = 0,linear_no_time = 1,linear = 2)
30: parameter (nonlinear_no_jacobian = 3,nonlinear = 4)
32: PetscErrorCode ierr
33: PetscInt time_steps,steps
34: PetscMPIInt size
35: integer problem
36: Vec local,global
37: double precision dt,ftime,zero,tmax
38: TS ts
39: Mat A
40: MatStructure A_structure
41: TSProblemType tsproblem
42: PetscDraw draw
43: PetscViewer viewer
44: character*(60) type,tsinfo
45: character*(5) etype
46: PetscInt i1,i60
47: PetscBool flg
48: PetscSizeT five
50: external Monitor,RHSFunctionHeat,RHSMatrixFree,Initial
51: external RHSMatrixHeat,RHSJacobianHeat
53: i1 = 1
54: i60 = 60
55: zero = 0.0
56: time_steps = 100
57: problem = linear_no_matrix
58: A = 0
59: tsproblem = TS_LINEAR
60:
61: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
62: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
64: M = 60
65: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-M',M,flg,ierr)
66: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-time',time_steps, &
67: & flg,ierr)
69: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-nox',flg,ierr)
70: if (flg) then
71: nox = 1
72: else
73: nox = 0
74: endif
75: nrm_2 = 0.0
76: nrm_max = 0.0
78: ! Set up the ghost point communication pattern
80: call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,M,i1,i1, &
81: & PETSC_NULL_INTEGER,da,ierr)
82: call DMCreateGlobalVector(da,global,ierr)
83: call VecGetLocalSize(global,m,ierr)
84: call DMCreateLocalVector(da,local,ierr)
86: ! Set up display to show wave graph
88: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER, &
89: & PETSC_NULL_CHARACTER,80,380,400,160,viewer1,ierr)
90: call PetscViewerDrawGetDraw(viewer1,0,draw,ierr)
91: call PetscDrawSetDoubleBuffer(draw,ierr)
92: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER, &
93: & PETSC_NULL_CHARACTER,80,0,400,160,viewer2,ierr)
94: call PetscViewerDrawGetDraw(viewer2,0,draw,ierr)
95: call PetscDrawSetDoubleBuffer(draw,ierr)
97: ! make work array for evaluating right hand side function
99: call VecDuplicate(local,localwork,ierr)
101: ! make work array for storing exact solution
103: call VecDuplicate(global,csolution,ierr)
105: h = 1.0/(M-1.0)
107: ! set initial conditions
108:
109: call Initial(global,PETSC_NULL_OBJECT,ierr)
110:
111: !
112: ! This example is written to allow one to easily test parts
113: ! of TS, we do not expect users to generally need to use more
114: ! then a single TSProblemType
115: !
116: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-linear_no_matrix',&
117: & flg,ierr)
118: if (flg) then
119: tsproblem = TS_LINEAR
120: problem = linear_no_matrix
121: endif
122: call PetscOptionsHasName(PETSC_NULL_CHARACTER, &
123: & '-linear_constant_matrix',flg,ierr)
124: if (flg) then
125: tsproblem = TS_LINEAR
126: problem = linear_no_time
127: endif
128: call PetscOptionsHasName(PETSC_NULL_CHARACTER, &
129: & '-nonlinear_no_jacobian',flg,ierr)
130: if (flg) then
131: tsproblem = TS_NONLINEAR
132: problem = nonlinear_no_jacobian
133: endif
134: call PetscOptionsHasName(PETSC_NULL_CHARACTER, &
135: & '-nonlinear_jacobian',flg,ierr)
136: if (flg) then
137: tsproblem = TS_NONLINEAR
138: problem = nonlinear
139: endif
140:
141: ! make timestep context
143: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
144: call TSSetProblemType(ts,tsproblem,ierr)
145: call TSMonitorSet(ts,Monitor,PETSC_NULL_OBJECT, &
146: & PETSC_NULL_FUNCTION, ierr)
148: dt = h*h/2.01
150: if (problem .eq. linear_no_matrix) then
151: !
152: ! The user provides the RHS as a Shell matrix.
153: !
154: call MatCreateShell(PETSC_COMM_WORLD,m,M,M,M, &
155: & PETSC_NULL_OBJECT,A,ierr)
156: call MatShellSetOperation(A,MATOP_MULT,RHSMatrixFree,ierr)
157: call TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant, &
158: & PETSC_NULL_OBJECT,ierr)
159: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT, &
160: & TSComputeRHSFunctionLinear,PETSC_NULL_OBJECT,ierr)
161: else if (problem .eq. linear_no_time) then
162: !
163: ! The user provides the RHS as a matrix
164: !
165: call MatCreate(PETSC_COMM_WORLD,A,ierr)
166: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M,ierr)
167: call MatSetFromOptions(A,ierr)
168: call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT, &
169: & ierr)
170: call TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant, &
171: & PETSC_NULL_OBJECT,ierr)
172: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT, &
173: & TSComputeRHSFunctionLinear,PETSC_NULL_OBJECT,ierr)
174: else if (problem .eq. nonlinear_no_jacobian) then
175: !
176: ! The user provides the RHS and a Shell Jacobian
177: !
178: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,RHSFunctionHeat, &
179: & PETSC_NULL_OBJECT,ierr)
180: call MatCreateShell(PETSC_COMM_WORLD,m,M,M,M, &
181: & PETSC_NULL_OBJECT,A,ierr)
182: call MatShellSetOperation(A,MATOP_MULT,RHSMatrixFree,ierr)
183: call TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant, &
184: & PETSC_NULL_OBJECT,ierr)
185: else if (problem .eq. nonlinear) then
186: !
187: ! The user provides the RHS and Jacobian
188: !
189: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,RHSFunctionHeat, &
190: & PETSC_NULL_OBJECT,ierr)
191: call MatCreate(PETSC_COMM_WORLD,A,ierr)
192: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,M,ierr)
193: call MatSetFromOptions(A,ierr)
194: call RHSMatrixHeat(ts,zero,A,A,A_structure,PETSC_NULL_OBJECT, &
195: & ierr)
196: call TSSetRHSJacobian(ts,A,A,RHSJacobianHeat, &
197: & PETSC_NULL_OBJECT,ierr)
198: endif
200: call TSSetFromOptions(ts,ierr)
202: call TSSetInitialTimeStep(ts,zero,dt,ierr)
203: tmax = 100.0
204: call TSSetDuration(ts,time_steps,tmax,ierr)
205: call TSSetSolution(ts,global,ierr)
207: call TSSetUp(ts,ierr)
208: call TSStep(ts,steps,ftime,ierr)
209: call PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,i60,viewer, &
210: & ierr)
211: call TSView(ts,viewer,ierr)
212: call TSGetType(ts,type,ierr)
214: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-test',flg,ierr)
215: if (flg) then
216: !
217: ! copy type to string of length 5 to ensure equality test
218: ! is done correctly
219: !
220: five = 5
221: call PetscStrncpy(etype,type,five,ierr)
222: if (etype .eq. TSEULER) then
223: if (abs(nrm_2/steps - 0.00257244) .gt. 1.e-4) then
224: print*,'Error in Euler method: 2-norm ',nrm_2/steps, &
225: & ' expecting: ',0.00257244
226: endif
227: else
228: if (abs(nrm_2/steps - 0.00506174) .gt. 1.e-4) then
229: print*,'Error in ',tsinfo,': 2-norm ',nrm_2/steps, &
230: & ' expecting: ',0.00506174
231: endif
232: endif
233: else
234: print*,size,' Procs Avg. error 2 norm ',nrm_2/steps, &
235: & nrm_max/steps,tsinfo
236: endif
238: call PetscViewerDestroy(viewer,ierr)
239: call TSDestroy(ts,ierr)
240: call PetscViewerDestroy(viewer1,ierr)
241: call PetscViewerDestroy(viewer2,ierr)
242: call VecDestroy(localwork,ierr)
243: call VecDestroy(csolution,ierr)
244: call VecDestroy(local,ierr)
245: call VecDestroy(global,ierr)
246: call DMDestroy(da,ierr)
247: if (A .ne. 0) then
248: call MatDestroy(A,ierr)
249: endif
251: call PetscFinalize(ierr)
252: end
254: ! -------------------------------------------------------------------
255:
256: subroutine Initial(global,ctx,ierr)
257: implicit none
258: #include <finclude/petscsys.h>
259: #include <finclude/petscvec.h>
260: #include <finclude/petscmat.h>
261: #include <finclude/petscdmda.h>
262: #include <finclude/petscpc.h>
263: #include <finclude/petscksp.h>
264: #include <finclude/petscsnes.h>
265: #include <finclude/petscts.h>
266: #include <finclude/petscviewer.h>
268: #include "ex1f.h"
270: Vec global
271: PetscObject ctx
273: PetscScalar localptr(1)
274: PetscInt i,mybase,myend
275: PetscErrorCode ierr
276: PetscOffset idx
278: ! determine starting point of each processor
280: call VecGetOwnershipRange(global,mybase,myend,ierr)
282: ! Initialize the array
284: call VecGetArray(global,localptr,idx,ierr)
285: do 10, i=mybase,myend-1
286: localptr(i-mybase+1+idx) = sin(PETSC_PI*i*6.*h) + &
287: & 3.*sin(PETSC_PI*i*2.*h)
288: 10 continue
289: call VecRestoreArray(global,localptr,idx,ierr)
290: return
291: end
293: ! ------------------------------------------------------------------------------
294: !
295: ! Exact solution
296: !
297: subroutine Solution(t,sol,ctx)
298: implicit none
299: #include <finclude/petscsys.h>
300: #include <finclude/petscvec.h>
301: #include <finclude/petscmat.h>
302: #include <finclude/petscdmda.h>
303: #include <finclude/petscpc.h>
304: #include <finclude/petscksp.h>
305: #include <finclude/petscsnes.h>
306: #include <finclude/petscts.h>
307: #include <finclude/petscviewer.h>
309: #include "ex1f.h"
311: double precision t
312: Vec sol
313: PetscObject ctx
315: PetscScalar localptr(1),ex1
316: PetscScalar ex2,sc1,sc2
317: PetscInt i,mybase,myend
318: PetscErrorCode ierr
319: PetscOffset idx
321: ! determine starting point of each processor
323: call VecGetOwnershipRange(csolution,mybase,myend,ierr)
325: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t)
326: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t)
327: sc1 = PETSC_PI*6.*h
328: sc2 = PETSC_PI*2.*h
329: call VecGetArray(csolution,localptr,idx,ierr)
330: do 10, i=mybase,myend-1
331: localptr(i-mybase+1+idx) = sin(i*sc1)*ex1 + 3.*sin(i*sc2)*ex2
332: 10 continue
333: call VecRestoreArray(csolution,localptr,idx,ierr)
334: return
335: end
338: ! -----------------------------------------------------------------------------------
340: subroutine Monitor(ts,step,time,global,ctx,ierr)
341: implicit none
342: #include <finclude/petscsys.h>
343: #include <finclude/petscvec.h>
344: #include <finclude/petscmat.h>
345: #include <finclude/petscdmda.h>
346: #include <finclude/petscpc.h>
347: #include <finclude/petscksp.h>
348: #include <finclude/petscsnes.h>
349: #include <finclude/petscts.h>
350: #include <finclude/petscdraw.h>
351: #include <finclude/petscviewer.h>
353: #include "ex1f.h"
354: TS ts
355: PetscInt step
356: PetscObject ctx
357: PetscErrorCode ierr
358: double precision time,lnrm_2,lnrm_max
359: Vec global
360: PetscScalar mone
362: call VecView(global,viewer1,ierr)
364: call Solution(time,csolution,ctx)
365: mone = -1.0
366: call VecAXPY(csolution,mone,global,ierr)
367: call VecNorm(csolution,NORM_2,lnrm_2,ierr)
368: lnrm_2 = sqrt(h)*lnrm_2
369: call VecNorm(csolution,NORM_MAX,lnrm_max,ierr)
371: if (nox .eq. 0) then
372: print*,'timestep ',step,' time ',time,' norm of error ', &
373: & lnrm_2,lnrm_max
374: endif
376: nrm_2 = nrm_2 + lnrm_2
377: nrm_max = nrm_max + lnrm_max
378: call VecView(csolution,viewer2,ierr)
380: return
381: end
383: ! -----------------------------------------------------------------------
385: subroutine RHSMatrixFree(mat,x,y,ierr)
386: implicit none
387: #include <finclude/petscsys.h>
388: #include <finclude/petscvec.h>
389: #include <finclude/petscmat.h>
390: #include <finclude/petscdmda.h>
391: #include <finclude/petscpc.h>
392: #include <finclude/petscksp.h>
393: #include <finclude/petscsnes.h>
394: #include <finclude/petscts.h>
395: #include <finclude/petscviewer.h>
396: Mat mat
397: Vec x,y
398: PetscErrorCode ierr
399: double precision zero
400: TS ts0
402: zero = 0.0
404: ts0 = PETSC_NULL_OBJECT
406: call RHSFunctionHeat(ts0,zero,x,y,PETSC_NULL_OBJECT,ierr)
407: return
408: end
410: ! -------------------------------------------------------------------------
412: subroutine RHSFunctionHeat(ts,t,globalin,globalout,ctx,ierr)
413: implicit none
414: #include <finclude/petscsys.h>
415: #include <finclude/petscvec.h>
416: #include <finclude/petscmat.h>
417: #include <finclude/petscdmda.h>
418: #include <finclude/petscpc.h>
419: #include <finclude/petscksp.h>
420: #include <finclude/petscsnes.h>
421: #include <finclude/petscts.h>
422: #include <finclude/petscviewer.h>
424: #include "ex1f.h"
425: TS ts
426: double precision t
427: Vec globalin,globalout
428: PetscObject ctx
429: Vec local
430: PetscErrorCode ierr
431: PetscInt i,localsize
432: PetscOffset ldx,cdx
433: PetscScalar copyptr(1),localptr(1)
434: PetscScalar sc
436: ! Extract local array
438: call DMCreateLocalVector(da,local,ierr)
439: call DMGlobalToLocalBegin(da,globalin,INSERT_VALUES,local,ierr)
440: call DMGlobalToLocalEnd(da,globalin,INSERT_VALUES,local,ierr)
441: call VecGetArray(local,localptr,ldx,ierr)
443: ! Extract work vector
445: call VecGetArray(localwork,copyptr,cdx,ierr)
447: ! Update Locally - Make array of new values
448: ! Note: For the first and last entry I copy the value
449: ! if this is an interior node it is irrelevant
451: sc = 1.0/(h*h)
452: call VecGetLocalSize(local,localsize,ierr)
453: copyptr(1+cdx) = localptr(1+ldx)
454: do 10, i=1,localsize-2
455: copyptr(i+1+cdx) = sc * (localptr(i+2+ldx) + localptr(i+ldx) - &
456: & 2.0*localptr(i+1+ldx))
457: 10 continue
458: copyptr(localsize-1+1+cdx) = localptr(localsize-1+1+ldx)
459: call VecRestoreArray(localwork,copyptr,cdx,ierr)
460: call VecRestoreArray(local,localptr,ldx,ierr)
461: call VecDestroy(local,ierr)
463: ! Local to Global
465: call DMLocalToGlobalBegin(da,localwork,INSERT_VALUES,globalout, &
466: & ierr)
467: call DMLocalToGlobalEnd(da,localwork,INSERT_VALUES,globalout, &
468: & ierr)
469: return
470: end
472: ! ---------------------------------------------------------------------
474: subroutine RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
475: implicit none
476: #include <finclude/petscsys.h>
477: #include <finclude/petscvec.h>
478: #include <finclude/petscmat.h>
479: #include <finclude/petscdmda.h>
480: #include <finclude/petscpc.h>
481: #include <finclude/petscksp.h>
482: #include <finclude/petscsnes.h>
483: #include <finclude/petscts.h>
484: #include <finclude/petscviewer.h>
486: #include "ex1f.h"
487: Mat AA,BB
488: double precision t
489: TS ts
490: MatStructure str
491: PetscObject ctx
492: PetscErrorCode ierr
494: Mat A
495: PetscInt i,mstart(1),mend(1),idx(3)
496: PetscMPIInt rank,size
497: PetscScalar v(3),stwo,sone
498: PetscInt i1,i3
500: i1 = 1
501: i3 = 3
502: A = AA
503: stwo = -2./(h*h)
504: sone = -.5*stwo
505: str = SAME_NONZERO_PATTERN
507: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
508: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
510: call MatGetOwnershipRange(A,mstart,mend,ierr)
511: if (mstart(1) .eq. 0) then
512: v(1) = 1.0
513: call MatSetValues(A,i1,mstart(1),i1,mstart(1),v,INSERT_VALUES, &
514: & ierr)
515: mstart(1) = mstart(1) + 1
516: endif
517: if (mend(1) .eq. M) then
518: mend(1) = mend(1) - 1
519: v(1) = 1.0
520: call MatSetValues(A,i1,mend,i1,mend,v,INSERT_VALUES,ierr)
521: endif
523: !
524: ! Construct matrice one row at a time
525: !
526: v(1) = sone
527: v(2) = stwo
528: v(3) = sone
529: do 10, i=mstart(1),mend(1)-1
530: idx(1) = i-1
531: idx(2) = i
532: idx(3) = i+1
533: call MatSetValues(A,i1,i,i3,idx,v,INSERT_VALUES,ierr)
534: 10 continue
536: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
537: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
538: return
539: end
541: ! --------------------------------------------------------------------------------------
543: subroutine RHSJacobianHeat(ts,t,x,AA,BB,str,ctx,ierr)
544: implicit none
545: #include <finclude/petscsys.h>
546: #include <finclude/petscvec.h>
547: #include <finclude/petscmat.h>
548: #include <finclude/petscdmda.h>
549: #include <finclude/petscpc.h>
550: #include <finclude/petscksp.h>
551: #include <finclude/petscsnes.h>
552: #include <finclude/petscts.h>
553: #include <finclude/petscviewer.h>
554: TS ts
555: double precision t
556: Vec x
557: Mat AA,BB
558: MatStructure str
559: PetscObject ctx
560: PetscErrorCode ierr
562: call RHSMatrixHeat(ts,t,AA,BB,str,ctx,ierr)
563: return
564: end