/[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 247 - (show annotations)
Tue Nov 29 04:57:38 2005 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 36678 byte(s)
initial checkin 

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 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