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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 251 - (show annotations)
Tue Nov 29 05:54:06 2005 UTC (13 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 36758 byte(s)
problems with coordinate axis fixed
1 # $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 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 f00=1.-f01
57 f11=(x_lat-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,long=0,lat=0):
80 self.long=long
81 self.lat=lat
82 def __str__(self):
83 return "(%s,%s)"%(self.long,self.lat)
84
85 def __sub__(self,other):
86 return self.dist(other)
87
88 def split(self,p,t):
89 return PointOnEarthSurface(long=self.long+t*(p.long-self.long),lat=self.lat+t*(p.lat-self.lat))
90
91 def midPoint(self,other):
92 return PointOnEarthSurface(long=(self.long+other.long)/2,lat=(self.lat+other.lat)/2)
93
94 def dist(self,other):
95 return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
96
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 def __init__(self,west_south=PointOnEarthSurface(),east_north=PointOnEarthSurface(),resolution=1.):
108 if resolution<=0:
109 raise ValueError,"resolution must be positive"
110 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 raise ValueError,"east south corner must be further down than west north corner"
114 if east_north.lat-west_south.lat < resolution/2:
115 raise ValueError,"latitude length of region must be 2*bigger then resolution"
116 if east_north.long-west_south.long < resolution/2:
117 raise ValueError,"longitude length of region must be 2*bigger then resolution"
118
119 self.west_south=west_south
120 self.east_north=east_north
121 self.resolution=resolution
122
123 def __str__(self):
124 return "RegionOnEarthSurface between %s and %s"%(str(self.west_south),str(self.east_north))
125
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 return self.west_south.long==p.long
131 if face==self.EAST:
132 return p.long==self.east_north.long
133 if face==self.SOUTH:
134 return self.west_south.lat==p.lat
135 if face==self.NORTH:
136 return p.lat==self.east_north.lat
137
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 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
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 if p0.lat<p1.lat:
163 return -1
164 elif p0.lat>p1.lat:
165 return 1
166 else:
167 return 0
168 elif f0 == self.EAST:
169 if p0.lat<p1.lat:
170 return 1
171 elif p0.lat>p1.lat:
172 return -1
173 else:
174 return 0
175 elif f0 == self.NORTH:
176 if p0.long<p1.long:
177 return -1
178 elif p0.long>p1.long:
179 return 1
180 else:
181 return 0
182 else:
183 if p0.long<p1.long:
184 return 1
185 elif p0.long>p1.long:
186 return -1
187 else:
188 return 0
189
190
191 def isInRegion(self,p):
192 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
197 def cutSegment(self,p0,p1):
198 t=None
199 p=None
200 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.west_south.long)
201 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 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.east_north.long)
208 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 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.west_south.lat)
215 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 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.east_north.lat)
222 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 a=p0.long
235 b=p1.long
236 else:
237 a=p0.lat
238 b=p1.lat
239 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 integ+=(p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
323 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 def makeTriangulation(self,west_south_is_water=True,east_south_is_water=True,west_north_is_water=True,east_north_is_water=True):
345 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 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 mid=v_tmp[0].midPoint(v_tmp[1])
376 n_long=-d_lat/l
377 n_lat=d_long/l
378 s=self.region.resolution
379 for seg in s_tmp:
380 p0=vertices[seg[0]]
381 p1=vertices[seg[1]]
382 a_long=p1.long-p0.long
383 a_lat=p1.lat-p0.lat
384 d=a_lat*n_long-a_long*n_lat
385 if abs(d)>0.:
386 t=((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
387 if 0<=t and t<=1:
388 s_tmp=((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
389 if s_tmp>EPS: s=min(s,s_tmp)
390 h=PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
391 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 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
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 walking_on_water=west_south_is_water
423 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 for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].long,vertices[i].lat))
555 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 for i in range(len(holes)): f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
559 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 reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
738 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 segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
747 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 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 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 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 self.trace("source length = %s"%(dist*2.))
785 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 m=s.whereNegative()
790 f1=(1.-m)*exp(-(s*beta)**2)+m
791 # x_long-direction
792 s=abs(x_long_hat)-dist
793 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 oc.resolution=2
932 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 src.end_lat=-12.
955 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 sq.t_end=100000.
969
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 sm.dt=5000.
976 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