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

  ViewVC Help
Powered by ViewVC 1.1.26