/[escript]/trunk/modellib/py_src/tsunami.py
ViewVC logotype

Annotation of /trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 247 - (hide annotations)
Tue Nov 29 04:57:38 2005 UTC (13 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 36678 byte(s)
initial checkin 

1 gross 247 # $Id:$
2    
3    
4     import os
5     from esys.escript import *
6     from esys.escript.linearPDEs import LinearPDE
7     from esys.escript.modelframe import Model
8     import numarray
9     import math
10     import urllib
11    
12     EPS=1.e-5
13    
14     #====================================================================================================================================================
15     class GridData:
16     """
17     this object is used to store data on grid.
18     it will be replaced by Bruce at a later stage.
19    
20     data[i,j] are data are x=j*s[0]+o[0] and y=i*s[1]+o[1]
21    
22     for 0<=j<n[0] and 0<=i<n[1]
23     """
24     def __init__(self,s,o,n,data):
25     self.s=tuple(s)
26     self.o=tuple(o)
27     self.n=tuple(n)
28     self.data=data
29    
30     def getOrigin(self):
31     return self.o
32     def getSpacing(self):
33     return self.s
34     def getData(self):
35     return self.data
36    
37     def interpolate(self,x):
38     if hasattr(x,"convertToNumArray"):
39     x_array=x.convertToNumArray()
40     return_data_object=True
41     else:
42     x_array=numarray.array(x)
43     return_data_object=False
44     data=numarray.zeros(x_array.shape[0],numarray.Float)
45     ox,oy=self.getOrigin()
46     dx,dy=self.getSpacing()
47     data_array=self.getData()
48     i_dx=1
49     i_dy=1
50     for i in range(x_array.shape[0]):
51     x0=x_array[i,0]-ox
52     x1=x_array[i,1]-oy
53     j0=min(max(int(x0/dx),0),data_array.shape[1]-1-i_dy)
54     j1=min(max(int(x1/dy),0),data_array.shape[0]-1-i_dx)
55     f01=(x0-j0*dx)/dx
56     f00=1.-f01
57     f11=(x1-j1*dy)/dy
58     f10=1.-f11
59     H00=data_array[j1,j0]
60     H01=data_array[j1,j0+i_dx]
61     H11=data_array[j1+i_dy,j0+i_dx]
62     H10=data_array[j1+i_dy,j0]
63     data[i]=(H00*f00+H01*f01)*f10+(H10*f00+H11*f01)*f11
64     if return_data_object:
65     out=Scalar(0,x.getFunctionSpace())
66     out.fillFromNumArray(data)
67     return out
68     else:
69     return data
70    
71     def getVTK(self):
72     pass
73    
74    
75     class PointOnEarthSurface:
76     """
77     coordinates of a point on the surface of the Earth
78     """
79     def __init__(self,phi=0,theta=0):
80     self.phi=phi
81     self.theta=theta
82     def __str__(self):
83     return "(%s,%s)"%(self.phi,self.theta)
84    
85     def __sub__(self,other):
86     return self.dist(other)
87    
88     def split(self,p,t):
89     phi=self.phi+t*(p.phi-self.phi)
90     theta=self.theta+t*(p.theta-self.theta)
91     return PointOnEarthSurface(phi,theta)
92    
93     def midPoint(self,other):
94     return PointOnEarthSurface((self.phi+other.phi)/2,(self.theta+other.theta)/2)
95    
96     def dist(self,other):
97     return math.sqrt((self.phi-other.phi)**2+(self.theta-other.theta)**2)
98    
99     class RegionOnEarthSurface:
100     """
101     defines an region by a south,east and north,west corner
102     """
103     RADIUS=6378.8e3
104     GRAVITY=9.81
105     WEST=0
106     NORTH=1
107     EAST=2
108     SOUTH=3
109     def __init__(self,south_east=PointOnEarthSurface(0.,0.),north_west=PointOnEarthSurface(10.,10.),resolution=1.):
110     if resolution<=0:
111     raise ValueError,"resolution must be positive"
112     if south_east.phi>=north_west.phi:
113     raise ValueError,"east south corner must be further east than west north corner"
114     if south_east.theta>=north_west.theta:
115     raise ValueError,"east south corner must be further down than west north corner"
116     if north_west.theta-south_east.theta < resolution/2:
117     raise ValueError,"latitude length of region must be 2*bigger then resolution"
118     if north_west.phi-south_east.phi < resolution/2:
119     raise ValueError,"longitude length of region must be 2*bigger then resolution"
120    
121     self.south_east=south_east
122     self.north_west=north_west
123     self.resolution=resolution
124    
125     def __str__(self):
126     return "RegionOnEarthSurface between %s and %s"%(str(self.south_east),str(self.north_west))
127    
128     def isOnFace(self,p):
129     return self.isOnThisFace(p,self.SOUTH) or self.isOnThisFace(p,self.NORTH) or self.isOnThisFace(p,self.WEST) or self.isOnThisFace(p,self.EAST)
130     def isOnThisFace(self,p,face):
131     if face==self.WEST:
132     return self.south_east.phi==p.phi
133     if face==self.EAST:
134     return p.phi==self.north_west.phi
135     if face==self.SOUTH:
136     return self.south_east.theta==p.theta
137     if face==self.NORTH:
138     return p.theta==self.north_west.theta
139    
140     def isBoxVertex(self,p):
141     return ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.SOUTH) ) or \
142     ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.NORTH) ) or \
143     ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.NORTH) ) or \
144     ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.SOUTH) )
145    
146    
147     def getFace(self,p):
148     # order is critical
149     if self.south_east.phi==p.phi: return self.WEST
150     if p.theta==self.north_west.theta: return self.NORTH
151     if p.phi==self.north_west.phi: return self.EAST
152     if self.south_east.theta==p.theta: return self.SOUTH
153    
154     def comparePointsOnFace(self,p0,p1):
155     f0=self.getFace(p0)
156     f1=self.getFace(p1)
157    
158     if f0<f1:
159     return -1
160     elif f0>f1:
161     return 1
162     else:
163     if f0 == self.WEST:
164     if p0.theta<p1.theta:
165     return -1
166     elif p0.theta>p1.theta:
167     return 1
168     else:
169     return 0
170     elif f0 == self.EAST:
171     if p0.theta<p1.theta:
172     return 1
173     elif p0.theta>p1.theta:
174     return -1
175     else:
176     return 0
177     elif f0 == self.NORTH:
178     if p0.phi<p1.phi:
179     return -1
180     elif p0.phi>p1.phi:
181     return 1
182     else:
183     return 0
184     else:
185     if p0.phi<p1.phi:
186     return 1
187     elif p0.phi>p1.phi:
188     return -1
189     else:
190     return 0
191    
192    
193     def isInRegion(self,p):
194     return self.south_east.phi<=p.phi \
195     and p.phi<=self.north_west.phi \
196     and self.south_east.theta<=p.theta \
197     and p.theta<=self.north_west.theta
198    
199     def cutSegment(self,p0,p1):
200     t=None
201     p=None
202     tmp=self.interseptOfSegment(p0,p1,d=0,v=self.south_east.phi)
203     if not tmp==None:
204     p_tmp=p0.split(p1,tmp)
205     if self.isInRegion(p_tmp) and t<tmp:
206     t=tmp
207     p=p_tmp
208    
209     tmp=self.interseptOfSegment(p0,p1,d=0,v=self.north_west.phi)
210     if not tmp==None:
211     p_tmp=p0.split(p1,tmp)
212     if self.isInRegion(p_tmp) and t<tmp:
213     t=tmp
214     p=p_tmp
215    
216     tmp=self.interseptOfSegment(p0,p1,d=1,v=self.south_east.theta)
217     if not tmp==None:
218     p_tmp=p0.split(p1,tmp)
219     if self.isInRegion(p_tmp) and t<tmp:
220     t=tmp
221     p=p_tmp
222    
223     tmp=self.interseptOfSegment(p0,p1,d=1,v=self.north_west.theta)
224     if not tmp==None:
225     p_tmp=p0.split(p1,tmp)
226     if self.isInRegion(p_tmp) and t<tmp:
227     t=tmp
228     p=p_tmp
229     return p
230    
231     def interseptOfSegment(self,p0,p1,d=0,v=0.):
232     """
233     find t isn [0,1] such that p0+t*(p0-p1) cuts x[d]=v. if it does nit exist None is returned.
234     """
235     if d==0:
236     a=p0.phi
237     b=p1.phi
238     else:
239     a=p0.theta
240     b=p1.theta
241     if b==a:
242     if a==v:
243     t=0
244     else:
245     t=None
246     else:
247     t=(v-a)/(b-a)
248     if not (0<=t and t<=1): t=None
249     return t
250    
251     class Polyline:
252     """
253     defines set of segments through a list of coordinates
254     """
255     def __init__(self,list_of_coordinates=[],name="none"):
256     c=[]
257     if len(list_of_coordinates)>0:
258     for i in range(len(list_of_coordinates)-1):
259     if list_of_coordinates[i]-list_of_coordinates[i+1]>0: c.append(list_of_coordinates[i])
260     c.append(list_of_coordinates[-1])
261     self.list_of_coordinates=c
262     self.name=name
263     def getDiameter(self):
264     out=0.
265     for p in self.list_of_coordinates:
266     for q in self.list_of_coordinates:
267     out=max(out,p-q)
268     return out
269     def isLoop(self):
270     if len(self)>0:
271     return not self[0]-self[-1]>EPS
272     else:
273     return True
274    
275     def insert(self,index,coordinate):
276     """
277     insert point before index
278     """
279     if self.list_of_coordinates[index]-coordinate<EPS:
280     return index
281     elif self.list_of_coordinates[index-1]-coordinate<EPS:
282     return index-1
283     else:
284     self.list_of_coordinates.insert(index,coordinate)
285     return index
286    
287     def split(self,index):
288     """
289     splits the polyline at point index
290     """
291     return Polyline(self.list_of_coordinates[:index],self.name),Polyline(self.list_of_coordinates[index:],self.name)
292    
293     def __getitem__(self,s):
294     return self.list_of_coordinates.__getitem__(s)
295     def __iter__(self):
296     return iter(self.list_of_coordinates)
297    
298     def __str__(self):
299     if self.isLoop():
300     out="loop["
301     else:
302     out="["
303     for i in self:
304     if out[-1]=="[":
305     out+="%s"%str(i)
306     else:
307     out+=",%s"%str(i)
308     return out+"]"
309    
310     def __len__(self):
311     return len(self.list_of_coordinates)
312    
313     def orientation(self):
314     """
315     returns the orientation of the polyline
316     """
317     if not self.isLoop():
318     raise TypeError,"polyline is not a loop"
319    
320     integ=0.
321     for i in range(len(self.list_of_coordinates)-1):
322     p0=self.list_of_coordinates[i]
323     p1=self.list_of_coordinates[i+1]
324     integ+=(p1.theta-p0.theta)*(p1.phi+p0.phi)-(p1.phi-p0.phi)*(p1.theta+p0.theta)
325     if integ>=0.:
326     return 1.
327     else:
328     return -1.
329    
330     def givePositiveOrientation(self):
331     if self.orientation()<0: self.list_of_coordinates.reverse()
332    
333     class Coastline:
334     """
335     defines a coast line by a Polyline within a RegionOnEarthSurface
336     """
337     def __init__(self,region,name="none"):
338     self.region=region
339     self.name=name
340     self.polylines=[]
341     def __str__(self):
342     out=self.name+" in "+str(self.region)
343     for pl in self.polylines:
344     out+="\n"+str(pl)
345     return out
346     def makeTriangulation(self,south_east_is_water=True,south_west_is_water=True,north_west_is_water=True,north_east_is_water=True):
347     self.clean()
348     vertices=[]
349     segments=[]
350     holes=[]
351     vertices_on_face=[]
352     for pl in self.polylines:
353     if pl.getDiameter()>self.region.resolution:
354     short_pl=[pl[0]]
355     for i in range(1,len(pl)):
356     if pl[i]-short_pl[-1]>self.region.resolution/5:
357     short_pl.append(pl[i])
358     elif i==len(pl)-1:
359     short_pl[-1]=pl[i]
360     if pl.isLoop():
361     if len(short_pl)>3:
362     offset=len(vertices)
363     v_tmp=[short_pl[0]]
364     s_tmp=[]
365     for i in range(1,len(short_pl)):
366     if i==len(short_pl)-1:
367     s_tmp.append((len(v_tmp)-1+offset,offset))
368     else:
369     s_tmp.append((len(v_tmp)-1+offset,len(v_tmp)+offset))
370     v_tmp.append(short_pl[i])
371     vertices+=v_tmp
372     segments+=s_tmp
373     # find a point in the loop:
374     d_phi=v_tmp[1].phi-v_tmp[0].phi
375     d_theta=v_tmp[1].theta-v_tmp[0].theta
376     l=math.sqrt(d_phi**2+d_theta**2)
377     mid=v_tmp[0].midPoint(v_tmp[1])
378     n_phi=-d_theta/l
379     n_theta=d_phi/l
380     s=self.region.resolution
381     for seg in s_tmp:
382     p0=vertices[seg[0]]
383     p1=vertices[seg[1]]
384     a_phi=p1.phi-p0.phi
385     a_theta=p1.theta-p0.theta
386     d=a_theta*n_phi-a_phi*n_theta
387     if abs(d)>0.:
388     t=((mid.theta-p0.theta)*n_phi-(mid.phi-p0.phi)*n_theta)/d
389     if 0<=t and t<=1:
390     s_tmp=((mid.theta-p0.theta)*a_phi-(mid.phi-p0.phi)*a_theta)/d
391     if s_tmp>EPS: s=min(s,s_tmp)
392     h=PointOnEarthSurface(mid.phi+s/2*n_phi,mid.theta+s/2*n_theta)
393     holes.append(h)
394     else:
395     if len(short_pl)>1:
396     if self.region.isOnFace(short_pl[0]): vertices_on_face.append(short_pl[0])
397     if self.region.isOnFace(short_pl[-1]): vertices_on_face.append(short_pl[-1])
398     vertices.append(short_pl[0])
399     for i in range(1,len(short_pl)):
400     segments.append((len(vertices)-1,len(vertices)))
401     vertices.append(short_pl[i])
402     # put into the bounding box:
403     new_vertices=[]
404     if south_east_is_water:
405     new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.south_east.theta))
406     if north_east_is_water:
407     new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.north_west.theta))
408     if north_west_is_water:
409     new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.north_west.theta))
410     if south_west_is_water:
411     new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.south_east.theta))
412    
413     # add new vertices if they don't exist yet
414     for q in new_vertices:
415     for p2 in vertices_on_face:
416     if p2-q<EPS:
417     q=None
418     raise ValueError,"coast line crosses boundary box vertex. This case is currenrly not supported."
419     if not q==None:
420     vertices.append(q)
421     vertices_on_face.append(q)
422     vertices_on_face.sort(self.region.comparePointsOnFace)
423     index=0
424     walking_on_water=south_east_is_water
425     l=len(vertices_on_face)
426     while index<l:
427     p1=vertices_on_face[(index+1)%l]
428     p0=vertices_on_face[index]
429     if walking_on_water:
430     segments.append((vertices.index(p0),vertices.index(p1)))
431     walking_on_water=False
432     else:
433     if self.region.isBoxVertex(p0):
434     segments.append((vertices.index(p0),vertices.index(p1)))
435     else:
436     walking_on_water=True
437     index+=1
438     return EarthTriangulation(vertices,segments,holes,self.region.resolution)
439    
440     def clean(self):
441     """
442     cleans up the coast line by joining polylines to loops or connecting faces of the region
443     """
444     # find a poylines that are linked
445     while True:
446     k0=None
447     for pl in self.polylines:
448     if not pl.isLoop():
449     for k in [0,-1]:
450     for ql in self.polylines:
451     if not (ql==pl or ql.isLoop()):
452     for k2 in [0,-1]:
453     if ql[k2]-pl[k]<EPS:
454     pl0=pl
455     pl1=ql
456     k0=k
457     k1=k2
458     break
459     if not k0==None: break # ql
460     if not k0==None: break # k
461     if not k0==None: break # pl
462    
463     if k0==None:
464     break
465     else:
466     self.polylines.remove(pl0)
467     self.polylines.remove(pl1)
468     pl0c=pl0.list_of_coordinates
469     pl1c=pl1.list_of_coordinates
470     if k0==0: pl0c.reverse()
471     if k1==-1: pl1c.reverse()
472     pl=Polyline(pl0c+pl1c[1:],pl0.name+" + "+pl1.name)
473     self.append(pl)
474    
475     # find a polyline that is not a loop and has an end or start point not on the face of the region:
476     while True:
477     pl=None
478     k=None
479     for pl2 in self.polylines:
480     if not pl2.isLoop():
481     pl=pl2
482     if not self.region.isOnFace(pl[0]): k=0
483     if not self.region.isOnFace(pl[-1]): k=-1
484     if not k==None: break
485     if k==None: break
486     self.polylines.remove(pl)
487     d_min=50000.
488     k_min=None
489     for pl2 in self.polylines:
490     if not pl2.isLoop():
491     for k2 in [0,-1]:
492     if not self.region.isOnFace(pl2[k2]):
493     d2=pl2[k2]-pl[k]
494     if d2<d_min:
495     d_min=d2
496     pl_min=pl2
497     k_min=k2
498     if k_min==None:
499     raise ValueError,"cannot link coastline %s to any other coastline."%pl.name
500     plc=pl.list_of_coordinates
501     plc_min=pl_min.list_of_coordinates
502     if k==0: plc.reverse()
503     if k_min==-1: plc_min.reverse()
504     if d_min<EPS:
505     new_pl=Polyline(plc+plc_min[1:],pl.name+" + "+pl_min.name)
506     else:
507     new_pl=Polyline(plc+plc_min,pl.name+" + "+pl_min.name)
508     self.polylines.remove(pl_min)
509     self.append(new_pl)
510     # give positive orientation to loops:
511     for pl in self.polylines:
512     if pl.isLoop(): pl.givePositiveOrientation()
513    
514     def append(self,polyline=Polyline()):
515     """append a polyline """
516     if len(polyline)>1:
517     pl=[]
518     outside_region=None
519     for i in range(len(polyline)):
520     if not self.region.isInRegion(polyline[i]):
521     outside_region=i
522     break
523     # pl.append(self.region.nudgeToFace(polyline[i]))
524     pl.append(polyline[i])
525     if not outside_region==None:
526     if outside_region==0:
527     for i in range(outside_region+1,len(polyline)):
528     if self.region.isInRegion(polyline[i]):
529     polyline.insert(i,self.region.cutSegment(polyline[i-1],\
530     polyline[i]))
531     pl1=polyline.split(i)[1]
532     self.append(pl1)
533     break
534     else:
535     # split polyline in two part first one is fully within the region the other starts with
536     # point outside the region
537     c=self.region.cutSegment(polyline[outside_region-1],polyline[outside_region])
538     i=polyline.insert(outside_region,c)
539     pl0,pl1=polyline.split(i+1)
540     self.append(pl0)
541     self.append(pl1)
542     else:
543     if len(pl)>1:
544     pply= Polyline(pl,polyline.name)
545     self.polylines.append(pply)
546    
547     class EarthTriangulation:
548     GENERATOR="triangle -pqa%g %s"
549    
550     def __init__(self,vertices=[],segments=[],holes=[],resolution=1.):
551     self.fn=os.tempnam()
552     # write triangle input file
553     poly_file=self.fn+".poly"
554     f=open(poly_file,"w")
555     f.writelines("%d %d %d %d\n"%(len(vertices),2,0,0))
556     for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].phi,vertices[i].theta))
557     f.writelines("%d %d\n"%(len(segments),0))
558     for i in range(len(segments)): f.writelines("%d %d %d\n"%(i,segments[i][0],segments[i][1]))
559     f.writelines("%d\n"%(len(holes)))
560     for i in range(len(holes)): f.writelines("%d %e %e\n"%(i,holes[i].phi,holes[i].theta))
561     f.close()
562     # start mesh generator:
563     os.system(self.GENERATOR%(resolution**2,poly_file))
564     # read mesh file:
565     self.node_coordinates=[]
566     self.node_tags=[]
567     self.node_ids=[]
568     self.triangles_nodes=[]
569     self.triangles_id=[]
570     node_file=open("%s.1.node"%self.fn,"r")
571     nodes=node_file.readline().strip().split()
572     nn=int(nodes[0])
573     for i in range(nn):
574     nodes=node_file.readline().strip().split()
575     self.node_coordinates.append((float(nodes[1]),float(nodes[2])))
576     self.node_tags.append(int(nodes[3]))
577     self.node_ids.append(int(nodes[0]))
578     node_file.close()
579     ele_file=open("%s.1.ele"%self.fn,"r")
580     elem=ele_file.readline().strip().split()
581     ne=int(elem[0])
582     for i in range(ne):
583     elem=ele_file.readline().strip().split()
584     self.triangles_id.append(int(elem[0]))
585     self.triangles_nodes.append((int(elem[1]),int(elem[2]),int(elem[3])))
586    
587     ele_file.close()
588     # os.remove("%s.1.node"%self.fn)
589     # os.remove("%s.1.ele"%self.fn)
590     # os.remove("%s.ploy"%self.fn)
591    
592     def getFinleyDomain(self):
593     from esys.finley import ReadMesh
594     finley_file=open("%s.msh"%self.fn,"w")
595     finley_file.writelines("%s\n2D-Nodes %d\n"%(self.fn,len(self.node_coordinates)))
596     for i in range(len(self.node_coordinates)):
597     finley_file.writelines("%s %s %s %e %e\n"%(self.node_ids[i],self.node_ids[i],self.node_tags[i],\
598     self.node_coordinates[i][0],self.node_coordinates[i][1]))
599    
600     finley_file.writelines("Tri3 %d\n"%len(self.triangles_nodes))
601     for i in range(len(self.triangles_nodes)):
602     finley_file.writelines("%s 0 %s %s %s\n"%(self.triangles_id[i], \
603     self.triangles_nodes[i][0], \
604     self.triangles_nodes[i][1], \
605     self.triangles_nodes[i][2]))
606     finley_file.writelines("Line2 %d\n"%0)
607     finley_file.writelines("Line2_Contact %d\n"%0)
608     finley_file.writelines("Point1 %d\n"%0)
609     finley_file.close()
610     # get the mesh
611     out=ReadMesh("%s.msh"%self.fn)
612     # os.remove("%s.msh"%self.fn)
613     return out
614    
615    
616     #=================================
617     # Model interfaces:
618     #================================
619     class OceanRegionCollector(Model):
620     """
621    
622    
623     """
624     def __init__(self,debug=False):
625     Model.__init__(self,debug=debug)
626     self.declareParameter(coastline_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/datafile.txt?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%&range360=%%range360%%&river=false&city=false&citytyp=wcity.dat&volcano=false&hotspot=false&shoreline=true&state=false&WSMdata=false&plate=false",
627     bathymetry_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/grdfile.xyz?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%",
628     resolution=1.,
629     south=0.,
630     north=10.,
631     east=0.,
632     west=20.,
633     range360=True,
634     coastline_stream=None,
635     bathymetry_stream=None)
636    
637    
638     def doInitialization(self):
639     """
640     Initializes the ocean region
641     """
642     c=self.__mergeParameters(self.coastline_source)
643     b=self.__mergeParameters(self.bathymetry_source)
644     self.coastline_stream=urllib.urlopen(c)
645     self.bathymetry_stream=urllib.urlopen(b)
646    
647     def __mergeParameters(self,txt):
648     return txt.replace("%%west%%",str(self.west))\
649     .replace("%%east%%",str(self.east))\
650     .replace("%%south%%",str(self.south))\
651     .replace("%%north%%",str(self.north))\
652     .replace("%%resolution%%",str(self.resolution)) \
653     .replace("%%range360%%",str(self.range360).lower())
654    
655     class Bathymetry(Model):
656     """
657     generates the bathymetry data within a region on the earth
658     """
659     def __init__(self,debug=False):
660     Model.__init__(self,debug=debug)
661     self.declareParameter(source="none",
662     bathymetry=1.)
663    
664     def doInitialization(self):
665     """
666     Initializes the
667     """
668     if hasattr(self.source,"readline"):
669     f=self.source
670     else:
671     f=open(filename,"r")
672     x_grd_list=[]
673     y_grd_list=[]
674     data_grd_list=[]
675     line=f.readline().strip()
676     while line!="":
677     v=line.split()
678     x_grd_list.append(float(v[0]))
679     y_grd_list.append(float(v[1]))
680     data_grd_list.append(float(v[2]))
681     line=f.readline().strip()
682     self.trace("%s data have been read from %s."%(len(data_grd_list),self.source))
683     data_grd=numarray.array(data_grd_list)
684     x_grd=numarray.array(x_grd_list)
685     y_grd=numarray.array(y_grd_list)
686     if len(x_grd)<2:
687     raise ValueError,"%s: data base is two small"%str(self)
688     ox=x_grd[0]
689     oy=y_grd[0]
690     diam=max(abs(x_grd[len(x_grd)-1]-ox),abs(y_grd[len(y_grd)-1]-oy))
691     dx=x_grd[1]-ox
692     nx=1
693     nx=1
694     while abs(x_grd[nx]-ox)>1.e-10*diam:
695     nx+=1
696     dy=y_grd[nx]-oy
697     ny=len(x_grd)/nx
698     data_grd.resize((ny,nx))
699     self.bathymetry=GridData(s=[dx,dy],o=[ox,oy],n=[nx,ny],data=data_grd)
700     self.trace("%s x %s grid with %s x %s spacing."%(nx,ny,dx,dy))
701    
702    
703     class OceanRegion(Model):
704     """
705     generates the ocean region with a coast line and a bathymetry
706    
707     """
708     def __init__(self,debug=False):
709     Model.__init__(self,debug=debug)
710     self.declareParameter(domain=None, \
711     resolution=1.,
712     south=0.,
713     north=10.,
714     east=0.,
715     west=20.,
716     bathymetry=None,
717     bathymetry_data=None,
718     coastline=None,
719     source="none")
720    
721     def doInitialization(self):
722     """
723     Initializes the ocean region
724     """
725     if hasattr(self.source,"readline"):
726     f=self.source
727     data_name=f.geturl()
728     else:
729     f=open(self.source,"r")
730     data_name=self.source
731    
732     segs=[]
733     name=""
734     line=f.readline()
735     while not line=="":
736     line=line.strip()
737     if line[:7]=="#region":
738     data=line[line.index("[")+1:line.index("]")].split(",")
739     reg=RegionOnEarthSurface(PointOnEarthSurface(self.south,self.west),PointOnEarthSurface(self.north,self.east),self.resolution)
740     self.trace(str(reg))
741     self.coastline=Coastline(region=reg,name=data_name)
742     elif line.find("Shore Bin")>-1:
743     self.coastline.append(Polyline(segs,name))
744     name=line[2:]
745     segs=[]
746     if not (line=="" or line[0]=="#" or line[0]==">") :
747     x=line.split()
748     segs.append(PointOnEarthSurface(float(x[0]),float(x[1])))
749     line=f.readline()
750     self.coastline.append(Polyline(segs,name))
751     d=self.bathymetry_data.interpolate([[self.east,self.south],[self.west,self.south],[self.east,self.north],[self.west,self.north]])
752     self.domain=self.coastline.makeTriangulation(south_east_is_water=d[0]<0,south_west_is_water=d[1]<0,north_east_is_water=d[2]<0,north_west_is_water=d[3]<0).getFinleyDomain()
753     self.bathymetry=maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
754    
755    
756     class TsunamiSource(Model):
757     """
758     defines a wave in Gaussean form between start and end.
759     """
760     GAMMA=0.05
761     def __init__(self,debug=False):
762     Model.__init__(self,debug=debug)
763     self.declareParameter(domain=None,
764     start_lat=-10.,
765     start_long=110.,
766     end_lat=-12.,
767     end_long=120.,
768     width=5.,
769     decay_zone=0.1,
770     amplitude=1.,
771     wave_height=1.)
772    
773     def doInitialization(self):
774     """
775     set initial wave form
776     """
777     beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
778     x=self.domain.getX()
779     x0=x[0]
780     x1=x[1]
781     mid0=(self.start_long+self.end_long)/2
782     mid1=(self.start_lat+self.end_lat)/2
783     dist=math.sqrt((mid0-self.end_long)**2+(mid1-self.end_lat)**2)
784     a=(self.start_long-mid0)/dist
785     b=(self.end_long-mid1)/dist
786     self.trace("source length = %s"%(dist*2.))
787    
788     x0_hat= a*(x0-mid0)+b*(x1-mid1)
789     x1_hat=-b*(x0-mid0)+a*(x1-mid1)
790     # x1-direction
791     s=abs(x1_hat)-self.width
792     m=s.whereNegative()
793     f1=(1.-m)*exp(-(s*beta)**2)+m
794     # x0-direction
795     s=abs(x0_hat)-dist
796     m=s.whereNegative()
797     f0=(1.-m)*exp(-(s*beta)**2)+m
798     self.wave_height=f1*f0*self.amplitude
799    
800     #====================================================================================================================================================
801     class TsunamiInDeepWater(Model):
802     """
803     Runs the deep water tsunami model based on a simplfied form of the shallow water equation.
804    
805    
806     M{d^2 h/dt^2 =div(c grad(h)) }
807    
808     where h is the wave height above sea level, and with H the depth of the water level and g is the gravity constant
809     c=sqrt(g*H).
810    
811     The simulation uses the Verlet scheme.
812    
813     """
814     def __init__(self,debug=False):
815     Model.__init__(self,debug=debug)
816     self.declareParameter(domain=None, \
817     wave_height=1.,
818     wave_velocity=0.,
819     initial_time_step=None,
820     bathymetry=1.,
821     safety_factor=0.1)
822    
823     def doInitialization(self):
824     """
825     Initializes the time integartion scheme
826     """
827     self.__pde=LinearPDE(self.domain)
828     self.__pde.setSolverMethod(self.__pde.LUMPING)
829     self.__pde.setValue(D=1.)
830     self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/RegionOnEarthSurface.RADIUS**2
831     c_max=math.sqrt(Lsup(self.__c2))
832     self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
833     if self.initial_time_step==None: self.initial_time_step=self.__dt
834     self.trace("maximum wave velocity %s m/sec^2"%c_max)
835     self.trace("Time step size is %s sec"%self.__dt)
836    
837    
838     def getSafeTimeStepSize(self,dt):
839     """
840     returns new step size
841    
842     @param dt: last time step size used
843     @type dt: C{float}
844     @return: time step size that can savely be used
845     @rtype: C{float}
846     """
847     return self.__dt
848    
849     def doStepPostprocessing(self,dt):
850     """
851     perform the time step using the Valet scheme
852    
853     @param dt: time step size to be used
854     @type dt: C{float}
855     """
856     self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
857     new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
858     self.wave_velocity=(new_height-self.wave_height)/dt
859     self.wave_height=new_height
860     self.initial_time_step=dt
861     self.trace("Wave height range is %e %e"%(inf(self.wave_height),sup(self.wave_height)))
862    
863     class SurfMovie(Model):
864     """
865     movie from a wave propagation on the sea
866    
867     @ivar time: current time
868     @ivar bathymetry: scalar data set
869     @ivar wave_height: vector data set
870     @ivar filename: name of the movie file
871     """
872     def __init__(self,debug=False):
873     Model.__init__(self,debug=debug)
874     self.declareParameter(bathymetry=1.,
875     wave_height=1.,
876     coastline=None,
877     t=0.,
878     dt=1.,
879     south=2.,
880     north=5.,
881     east=3.,
882     west=15.,
883     filename="movie.mpg")
884    
885     def doInitialization(self):
886     """
887     Initializes the time integartion scheme
888     """
889     self.__fn=os.tempnam()+".xml"
890     self.__frame_name=os.tempnam()
891     self.__next_t=self.dt
892     # self.coastline.getVTK()
893     # self.bathymetry.getVTK()
894     # wndow(south,west,north,east)
895    
896     def doStepPostprocessing(self, dt):
897     """
898     Does any necessary postprocessing after each step
899    
900     @param dt:
901     """
902     if self.t>=self.__next_t:
903     print self.t,"write ",Lsup(self.wave_height)
904     saveVTK(self.__fn,h=self.wave_height)
905     # vtkobj=...
906     # save(self.__frame_name)
907     self.__next_t+=self.dt
908    
909     def getSafeTimeStepSize(self,dt):
910     """
911     returns new step size
912    
913     @param dt: last time step size used
914     @type dt: C{float}
915     @return: time step size that can savely be used
916     @rtype: C{float}
917     """
918     return self.__next_t-self.t
919    
920     def doFinalization(self):
921     """
922     Finalises the visualisation. For instance, makes a movie of the
923     image files.
924     """
925     # make the movie into self.filename
926     pass
927    
928    
929     if __name__=="__main__":
930     from esys.escript.modelframe import Link,Simulation
931     from esys.modellib.input import Sequencer
932    
933     oc=OceanRegionCollector()
934     oc.resolution=0.5
935     oc.south=-45.5
936     oc.north=-5.4
937     oc.east=161.1
938     oc.west=108.2
939     oc.range360=True
940    
941    
942     b=Bathymetry()
943     b.source=Link(oc,"bathymetry_stream")
944    
945     oreg=OceanRegion()
946     oreg.source=Link(oc,"coastline_stream")
947     oreg.resolution=Link(oc,"resolution")
948     oreg.south=Link(oc,"south")
949     oreg.north=Link(oc,"north")
950     oreg.east=Link(oc,"east")
951     oreg.west=Link(oc,"west")
952     oreg.bathymetry_data=Link(b,"bathymetry")
953    
954     src=TsunamiSource()
955     src.domain=Link(oreg,"domain")
956     src.start_lat=-10.
957     src.start_long=110.
958     src.end_lat=-12.
959     src.end_long=120.
960     src.width=0.1
961     src.decay_zone=0.01
962     src.amplitude=1.
963    
964     ts=TsunamiInDeepWater()
965     ts.domain=Link(oreg,"domain")
966     ts.wave_height=Link(src,"wave_height")
967     ts.wave_velocity=0.
968     ts.bathymetry=Link(oreg,"bathymetry")
969    
970     sq=Sequencer()
971     sq.t_end=1000.
972    
973     sm=SurfMovie()
974     sm.bathymetry=Link(b,"bathymetry")
975     sm.wave_height=Link(ts,"wave_height")
976     sm.coastline=Link(oreg,"coastline")
977     sm.t=Link(sq,"t")
978     sm.dt=100.
979     sm.filename="movie.mpg"
980    
981     s=Simulation([sq,oc,b,oreg,src,ts,sm])
982     s.writeXML()
983     s.run()

  ViewVC Help
Powered by ViewVC 1.1.26