/[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 2067 - (show annotations)
Thu Nov 20 06:22:56 2008 UTC (10 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 50580 byte(s)
Moved 'make dir' operation to the correct spot and extended the number
of timesteps a bit.


1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 import os, sys
23 import vtk
24 from esys.escript import *
25 from esys.escript.linearPDEs import LinearPDE
26 from esys.escript.modelframe import Model
27 import numarray
28 import math
29 from tempfile import mkstemp
30
31 WORKDIR="./work/"
32 EPS=1.e-8
33
34 #=============================================================================
35
36 class GridData:
37 """
38 This object is used to store data on a grid.
39 It will be replaced by Bruce at a later stage.
40
41 data[i,j] are x=j*s[0]+o[0] and y=i*s[1]+o[1]
42
43 for 0<=j<n[0] and 0<=i<n[1]
44 """
45 def __init__(self, s, o, n, data):
46 self.s=tuple(s)
47 self.o=tuple(o)
48 self.n=tuple(n)
49 self.data=data
50
51 def getOrigin(self):
52 return self.o
53
54 def getSpacing(self):
55 return self.s
56
57 def getData(self):
58 return self.data
59
60 def interpolate(self,x):
61 if hasattr(x, "getNumberOfDataPoints"):
62 x_shape0 = x.getNumberOfDataPoints()
63 return_data_object = True
64 else:
65 x_array = numarray.array(x)
66 x_shape0 = x_array.shape[0]
67 return_data_object = False
68
69 data=numarray.zeros(x_shape0, numarray.Float)
70 ox,oy = self.getOrigin()
71 dx,dy = self.getSpacing()
72 data_array = self.getData()
73 i_dx = 1
74 i_dy = 1
75 for i in range(x_shape0):
76 if return_data_object:
77 x_array = x.getValueOfDataPoint(i)
78 x_long = x_array[0]-ox
79 x_lat = x_array[1]-oy
80 else:
81 x_long = x_array[i,0]-ox
82 x_lat = x_array[i,1]-oy
83 j0 = min(max(int(x_long/dx),0), data_array.shape[1]-1-i_dy)
84 j1 = min(max(int(x_lat/dy),0), data_array.shape[0]-1-i_dx)
85 f01 = (x_long-j0*dx)/dx
86 f00 = 1.-f01
87 f11 = (x_lat-j1*dy)/dy
88 f10 = 1.-f11
89 H00 = data_array[j1,j0]
90 H01 = data_array[j1,j0+i_dx]
91 H11 = data_array[j1+i_dy,j0+i_dx]
92 H10 = data_array[j1+i_dy,j0]
93 data[i] = (H00*f00+H01*f01)*f10 + (H10*f00+H11*f01)*f11
94
95 if return_data_object:
96 dataObj = Scalar(0, x.getFunctionSpace(), True)
97 for i in range(x_shape0):
98 dataObj.setValueOfDataPoint(i, data[i])
99 return dataObj
100 else:
101 return data
102
103
104 class PointOnEarthSurface:
105 """
106 Coordinates of a point on the surface of the Earth
107 """
108 def __init__(self, long=0, lat=0):
109 self.long=long
110 self.lat=lat
111
112 def __str__(self):
113 return "(%s,%s)"%(self.long,self.lat)
114
115 def __sub__(self, other):
116 return self.dist(other)
117
118 def split(self,p,t):
119 return PointOnEarthSurface(long=self.long+t*(p.long-self.long),
120 lat=self.lat+t*(p.lat-self.lat) )
121
122 def midPoint(self,other):
123 return PointOnEarthSurface(long=(self.long+other.long)/2,
124 lat=(self.lat+other.lat)/2 )
125
126 def dist(self,other):
127 return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
128
129
130 class RegionOnEarthSurface:
131 """
132 Defines a region by a south-east and north-west corner
133 """
134 RADIUS=6378.8e3
135 GRAVITY=9.81
136 WEST=0
137 NORTH=1
138 EAST=2
139 SOUTH=3
140 def __init__(self, west_south=PointOnEarthSurface(), east_north=PointOnEarthSurface(), resolution=1.):
141 if resolution <= 0:
142 raise ValueError, "resolution must be positive"
143 if west_south.long >= east_north.long:
144 raise ValueError, "south-west corner must be west of north-east corner"
145 if west_south.lat >= east_north.lat:
146 raise ValueError, "south-east corner must be south of north-west corner"
147 if east_north.lat-west_south.lat < resolution/2:
148 raise ValueError, "latitude length of region must be at least 2*larger than the resolution"
149 if east_north.long-west_south.long < resolution/2:
150 raise ValueError, "longitude length of region must be at least 2*larger than the resolution"
151
152 self.west_south = west_south
153 self.east_north = east_north
154 self.resolution = resolution
155
156 def __str__(self):
157 return "RegionOnEarthSurface between %s and %s" % (str(self.west_south), str(self.east_north))
158
159 def isOnFace(self, p):
160 return self.isOnThisFace(p,self.SOUTH) or \
161 self.isOnThisFace(p,self.NORTH) or \
162 self.isOnThisFace(p,self.WEST) or \
163 self.isOnThisFace(p,self.EAST)
164
165 def isOnThisFace(self, p, face):
166 if face==self.WEST:
167 return self.west_south.long==p.long
168 if face==self.EAST:
169 return p.long==self.east_north.long
170 if face==self.SOUTH:
171 return self.west_south.lat==p.lat
172 if face==self.NORTH:
173 return p.lat==self.east_north.lat
174
175 def isBoxVertex(self, p):
176 return ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.SOUTH) ) or \
177 ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.NORTH) ) or \
178 ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.NORTH) ) or \
179 ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.SOUTH) )
180
181 def getFace(self, p):
182 # order is critical
183 if self.west_south.long == p.long: return self.WEST
184 if p.lat == self.east_north.lat: return self.NORTH
185 if p.long == self.east_north.long: return self.EAST
186 if self.west_south.lat == p.lat: return self.SOUTH
187
188 def comparePointsOnFace(self, p0, p1):
189 f0 = self.getFace(p0)
190 f1 = self.getFace(p1)
191
192 if f0 < f1:
193 return -1
194 elif f0 > f1:
195 return 1
196 else:
197 if f0 == self.WEST:
198 if p0.lat < p1.lat:
199 return -1
200 elif p0.lat > p1.lat:
201 return 1
202 else:
203 return 0
204 elif f0 == self.EAST:
205 if p0.lat < p1.lat:
206 return 1
207 elif p0.lat > p1.lat:
208 return -1
209 else:
210 return 0
211 elif f0 == self.NORTH:
212 if p0.long < p1.long:
213 return -1
214 elif p0.long > p1.long:
215 return 1
216 else:
217 return 0
218 else:
219 if p0.long < p1.long:
220 return 1
221 elif p0.long > p1.long:
222 return -1
223 else:
224 return 0
225
226 def isInRegion(self, p):
227 return self.west_south.long <= p.long \
228 and p.long <= self.east_north.long \
229 and self.west_south.lat <= p.lat \
230 and p.lat <= self.east_north.lat
231
232 def cutSegment(self, p0, p1):
233 t = None
234 p = None
235 tmp = self.interceptOfSegment(p0, p1, d=0, v=self.west_south.long)
236 if not tmp == None:
237 p_tmp = p0.split(p1, tmp)
238 if self.isInRegion(p_tmp) and t < tmp:
239 t = tmp
240 p = p_tmp
241
242 tmp = self.interceptOfSegment(p0, p1, d=0, v=self.east_north.long)
243 if not tmp == None:
244 p_tmp = p0.split(p1, tmp)
245 if self.isInRegion(p_tmp) and t < tmp:
246 t = tmp
247 p = p_tmp
248
249 tmp = self.interceptOfSegment(p0, p1, d=1, v=self.west_south.lat)
250 if not tmp == None:
251 p_tmp = p0.split(p1, tmp)
252 if self.isInRegion(p_tmp) and t < tmp:
253 t = tmp
254 p = p_tmp
255
256 tmp = self.interceptOfSegment(p0, p1, d=1, v=self.east_north.lat)
257 if not tmp == None:
258 p_tmp = p0.split(p1, tmp)
259 if self.isInRegion(p_tmp) and t < tmp:
260 t = tmp
261 p = p_tmp
262 return p
263
264 def interceptOfSegment(self, p0, p1, d=0, v=0.):
265 """
266 Finds t in [0,1] such that p0+t*(p0-p1) cuts x[d]=v. If it does not
267 exist None is returned.
268 """
269 if d == 0:
270 a = p0.long
271 b = p1.long
272 else:
273 a = p0.lat
274 b = p1.lat
275 if b == a:
276 if a == v:
277 t = 0
278 else:
279 t = None
280 else:
281 t = (v-a)/(b-a)
282 if not (0<=t and t<=1):
283 t = None
284 return t
285
286
287 class Polyline:
288 """
289 Defines a set of segments through a list of coordinates
290 """
291 def __init__(self, list_of_coordinates=[], name="none"):
292 c=[]
293 if len(list_of_coordinates)>0:
294 for i in range(len(list_of_coordinates)-1):
295 if list_of_coordinates[i]-list_of_coordinates[i+1] > 0:
296 c.append(list_of_coordinates[i])
297 c.append(list_of_coordinates[-1])
298 self.list_of_coordinates = c
299 self.name = name
300
301 def getDiameter(self):
302 out = 0.
303 for p in self.list_of_coordinates:
304 for q in self.list_of_coordinates:
305 out = max(out, p-q)
306 return out
307
308 def isLoop(self):
309 if len(self) > 0:
310 return not self[0]-self[-1]>EPS
311 else:
312 return True
313
314 def insert(self, index, coordinate):
315 """
316 Inserts point before index
317 """
318 if self.list_of_coordinates[index]-coordinate < EPS:
319 return index
320 elif self.list_of_coordinates[index-1]-coordinate < EPS:
321 return index-1
322 else:
323 self.list_of_coordinates.insert(index, coordinate)
324 return index
325
326 def split(self, index):
327 """
328 Splits the polyline at point index
329 """
330 return Polyline(self.list_of_coordinates[:index], self.name), \
331 Polyline(self.list_of_coordinates[index:], self.name)
332
333 def __getitem__(self, s):
334 return self.list_of_coordinates.__getitem__(s)
335
336 def __iter__(self):
337 return iter(self.list_of_coordinates)
338
339 def __str__(self):
340 if self.isLoop():
341 out = "loop["
342 else:
343 out = "["
344 for i in self:
345 if out[-1] == "[":
346 out += "%s" % str(i)
347 else:
348 out += ",%s" % str(i)
349 return out+"]"
350
351 def __len__(self):
352 return len(self.list_of_coordinates)
353
354 def orientation(self):
355 """
356 Returns the orientation of the polyline
357 """
358 if not self.isLoop():
359 raise TypeError, "polyline is not a loop"
360
361 integ = 0.
362 for i in range(len(self.list_of_coordinates)-1):
363 p0 = self.list_of_coordinates[i]
364 p1 = self.list_of_coordinates[i+1]
365 integ += (p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
366 if integ >= 0.:
367 return 1.
368 else:
369 return -1.
370
371 def givePositiveOrientation(self):
372 if self.orientation() < 0:
373 self.list_of_coordinates.reverse()
374
375
376 class Coastline:
377 """
378 Defines a coast line by a Polyline within a RegionOnEarthSurface
379 """
380 def __init__(self, region, name="none"):
381 self.region = region
382 self.name = name
383 self.polylines = []
384
385 def __str__(self):
386 out = self.name + " in " + str(self.region)
387 for pl in self.polylines:
388 out += "\n" + str(pl)
389 return out
390
391 def makeTriangulation(self, west_south_is_water=True, east_south_is_water=True, west_north_is_water=True, east_north_is_water=True):
392 self.clean()
393 vertices=[]
394 segments=[]
395 holes=[]
396 vertices_on_face=[]
397 for pl in self.polylines:
398 if pl.getDiameter() > self.region.resolution:
399 short_pl = [pl[0]]
400 for i in range(1, len(pl)):
401 if pl[i]-short_pl[-1] > 0*EPS+self.region.resolution/10:
402 short_pl.append(pl[i])
403 elif i == len(pl)-1:
404 short_pl[-1] = pl[i]
405 if pl.isLoop():
406 if len(short_pl) > 3:
407 offset = len(vertices)
408 v_tmp = [short_pl[0]]
409 s_tmp = []
410 for i in range(1, len(short_pl)):
411 if i == len(short_pl)-1:
412 s_tmp.append((len(v_tmp)-1+offset, offset))
413 else:
414 s_tmp.append((len(v_tmp)-1+offset, len(v_tmp)+offset))
415 v_tmp.append(short_pl[i])
416 vertices += v_tmp
417 segments += s_tmp
418 # find a point in the loop:
419 d_long = v_tmp[1].long-v_tmp[0].long
420 d_lat = v_tmp[1].lat-v_tmp[0].lat
421 l = math.sqrt(d_long**2+d_lat**2)
422 mid = v_tmp[0].midPoint(v_tmp[1])
423 n_long = -d_lat/l
424 n_lat = d_long/l
425 s = self.region.resolution
426 for seg in s_tmp:
427 p0 = vertices[seg[0]]
428 p1 = vertices[seg[1]]
429 a_long = p1.long-p0.long
430 a_lat = p1.lat-p0.lat
431 d = a_lat*n_long-a_long*n_lat
432 if abs(d) > 0.:
433 t = ((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
434 if 0<=t and t<=1:
435 s_tmp = ((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
436 if s_tmp > EPS:
437 s = min(s,s_tmp)
438 h = PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
439 holes.append(h)
440 else:
441 if len(short_pl) > 1:
442 if self.region.isOnFace(short_pl[0]):
443 vertices_on_face.append(short_pl[0])
444 if self.region.isOnFace(short_pl[-1]):
445 vertices_on_face.append(short_pl[-1])
446 vertices.append(short_pl[0])
447 for i in range(1, len(short_pl)):
448 segments.append((len(vertices)-1, len(vertices)))
449 vertices.append(short_pl[i])
450 # put into the bounding box:
451 new_vertices = []
452 if west_south_is_water:
453 new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long, lat=self.region.west_south.lat))
454 if east_south_is_water:
455 new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long, lat=self.region.west_south.lat))
456 if west_north_is_water:
457 new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long, lat=self.region.east_north.lat))
458 if east_north_is_water:
459 new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long, lat=self.region.east_north.lat))
460
461 # add new vertices if they don't exist yet
462 for q in new_vertices:
463 for p2 in vertices_on_face:
464 if p2-q < EPS:
465 q = None
466 raise ValueError, "coast line crosses boundary box vertex. This case is currently not supported."
467 if not q == None:
468 vertices.append(q)
469 vertices_on_face.append(q)
470 vertices_on_face.sort(self.region.comparePointsOnFace)
471 index = 0
472 walking_on_water = True
473 l = len(vertices_on_face)
474 while index < l:
475 p1 = vertices_on_face[(index+1)%l]
476 p0 = vertices_on_face[index]
477 if walking_on_water:
478 segments.append((vertices.index(p0), vertices.index(p1)))
479 walking_on_water = False
480 else:
481 if self.region.isBoxVertex(p0):
482 segments.append((vertices.index(p0), vertices.index(p1)))
483 else:
484 walking_on_water = True
485 index += 1
486 return EarthTriangulation(vertices, segments, holes, self.region.resolution)
487
488 def clean(self):
489 """
490 Cleans up the coastline by joining polylines to loops or connecting
491 faces of the region
492 """
493 # find poylines that are linked
494 while True:
495 k0 = None
496 for pl in self.polylines:
497 if not pl.isLoop():
498 for k in [0,-1]:
499 for ql in self.polylines:
500 if not (ql==pl or ql.isLoop()):
501 for k2 in [0,-1]:
502 if ql[k2]-pl[k] < EPS:
503 pl0 = pl
504 pl1 = ql
505 k0 = k
506 k1 = k2
507 break
508 if not k0==None: break # ql
509 if not k0==None: break # k
510 if not k0==None: break # pl
511
512 if k0 == None:
513 break
514 else:
515 self.polylines.remove(pl0)
516 self.polylines.remove(pl1)
517 pl0c = pl0.list_of_coordinates
518 pl1c = pl1.list_of_coordinates
519 if k0 == 0: pl0c.reverse()
520 if k1 == -1: pl1c.reverse()
521 pl = Polyline(pl0c+pl1c[1:], pl0.name+" + "+pl1.name)
522 self.append(pl)
523
524 # find a polyline that is not a loop and has an end or start point
525 # not on the face of the region:
526 while True:
527 pl = None
528 k = None
529 for pl2 in self.polylines:
530 if not pl2.isLoop():
531 pl = pl2
532 if not self.region.isOnFace(pl[0]): k=0
533 if not self.region.isOnFace(pl[-1]): k=-1
534 if not k == None: break
535 if k == None: break
536 self.polylines.remove(pl)
537 d_min = 50000.
538 k_min = None
539 for pl2 in self.polylines:
540 if not pl2.isLoop():
541 for k2 in [0,-1]:
542 if not self.region.isOnFace(pl2[k2]):
543 d2 = pl2[k2]-pl[k]
544 if d2 < d_min:
545 d_min = d2
546 pl_min = pl2
547 k_min = k2
548 if k_min == None:
549 raise ValueError, "cannot link coastline %s to any other coastline." % pl.name
550 plc = pl.list_of_coordinates
551 plc_min = pl_min.list_of_coordinates
552 if k == 0: plc.reverse()
553 if k_min == -1: plc_min.reverse()
554 if d_min < EPS:
555 new_pl = Polyline(plc+plc_min[1:], pl.name+" + "+pl_min.name)
556 else:
557 new_pl = Polyline(plc+plc_min, pl.name+" + "+pl_min.name)
558 self.polylines.remove(pl_min)
559 self.append(new_pl)
560 # give positive orientation to loops:
561 for pl in self.polylines:
562 if pl.isLoop(): pl.givePositiveOrientation()
563
564 def append(self, polyline=Polyline()):
565 """
566 Appends a polyline
567 """
568 if len(polyline) > 1:
569 pl = []
570 outside_region = None
571 for i in range(len(polyline)):
572 if not self.region.isInRegion(polyline[i]):
573 outside_region = i
574 break
575 #pl.append(self.region.nudgeToFace(polyline[i]))
576 pl.append(polyline[i])
577 if not outside_region == None:
578 if outside_region == 0:
579 for i in range(outside_region+1, len(polyline)):
580 if self.region.isInRegion(polyline[i]):
581 polyline.insert(i,self.region.cutSegment( \
582 polyline[i-1], polyline[i]))
583 pl1 = polyline.split(i)[1]
584 self.append(pl1)
585 break
586 else:
587 # split polyline in two parts: first is fully within the
588 # region, the other starts with point outside the region
589 c = self.region.cutSegment(polyline[outside_region-1], polyline[outside_region])
590 i = polyline.insert(outside_region, c)
591 pl0,pl1 = polyline.split(i+1)
592 self.append(pl0)
593 self.append(pl1)
594 else:
595 if len(pl) > 1:
596 pply = Polyline(pl,polyline.name)
597 self.polylines.append(pply)
598
599
600 class EarthTriangulation:
601 """
602 Generates earth mesh and triangulates it
603 """
604 GENERATOR = "triangle -pqa%g %s"
605
606 def __init__(self,vertices=[], segments=[], holes=[], resolution=1.):
607 (fd, self.fn) = mkstemp(suffix=".poly")
608 self.fn = self.fn.replace(".poly", "")
609 f=os.fdopen(fd, "w")
610 # write triangle input file
611 f.writelines("%d %d %d %d\n" % (len(vertices),2,0,0))
612 for i in range(len(vertices)):
613 f.writelines("%d %e %e\n" % (i,vertices[i].long,vertices[i].lat))
614 f.writelines("%d %d\n"%(len(segments),0))
615 for i in range(len(segments)):
616 f.writelines("%d %d %d\n" % (i,segments[i][0],segments[i][1]))
617 f.writelines("%d\n"%(len(holes)))
618 for i in range(len(holes)):
619 f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
620 f.close()
621
622 # create mesh using generator tool
623 os.system(self.GENERATOR % (resolution**2, self.fn))
624
625 # read the resulting mesh file
626 self.node_coordinates=[]
627 self.node_tags=[]
628 self.node_ids=[]
629 self.triangles_nodes=[]
630 self.triangles_id=[]
631 node_file=open("%s.1.node" % self.fn, "r")
632 nodes=node_file.readline().strip().split()
633 nn=int(nodes[0])
634 for i in range(nn):
635 nodes=node_file.readline().strip().split()
636 self.node_coordinates.append((float(nodes[1]),float(nodes[2])))
637 self.node_tags.append(int(nodes[3]))
638 self.node_ids.append(int(nodes[0]))
639 node_file.close()
640 ele_file=open("%s.1.ele" % self.fn, "r")
641 elem=ele_file.readline().strip().split()
642 ne=int(elem[0])
643 for i in range(ne):
644 elem=ele_file.readline().strip().split()
645 self.triangles_id.append(int(elem[0]))
646 self.triangles_nodes.append((int(elem[1]),int(elem[2]),int(elem[3])))
647 ele_file.close()
648
649 # clean up
650 os.remove("%s.poly" % self.fn)
651 os.remove("%s.1.poly" % self.fn)
652 os.remove("%s.1.node" % self.fn)
653 os.remove("%s.1.ele" % self.fn)
654
655 def getFinleyDomain(self):
656 from esys.finley import ReadMesh
657 finley_file = open("%s.msh" % self.fn, "w")
658 finley_file.writelines("%s\n2D-Nodes %d\n" % (self.fn,
659 len(self.node_coordinates)))
660 for i in range(len(self.node_coordinates)):
661 finley_file.writelines("%s %s %s %e %e\n" % (self.node_ids[i],
662 self.node_ids[i],self.node_tags[i],
663 self.node_coordinates[i][0],
664 self.node_coordinates[i][1]))
665
666 finley_file.writelines("Tri3 %d\n" % len(self.triangles_nodes))
667 for i in range(len(self.triangles_nodes)):
668 finley_file.writelines("%s 0 %s %s %s\n" % (
669 self.triangles_id[i],
670 self.triangles_nodes[i][0],
671 self.triangles_nodes[i][1],
672 self.triangles_nodes[i][2]))
673 finley_file.writelines("Line2 0\n")
674 finley_file.writelines("Line2_Contact 0\n")
675 finley_file.writelines("Point1 0\n")
676 finley_file.close()
677 # read the mesh with finley
678 out=ReadMesh("%s.msh" % self.fn)
679 os.remove("%s.msh" % self.fn)
680 return out
681
682
683 #=================================
684 # Model interfaces #
685 #=================================
686 class OceanRegionCollector(Model):
687 """
688 Opens input streams for coastline and bathymetry data (generated by GMT)
689 """
690
691 def __init__(self,**kwargs):
692 Model.__init__(self,**kwargs)
693 self.declareParameter(coastline_source = "coastline_default.dat",
694 bathymetry_source = "bathymetry_default.dat",
695 resolution = 1.,
696 south = 0.,
697 north = 10.,
698 east = 0.,
699 west = 20.,
700 range360 = True,
701 coastline_stream = None,
702 bathymetry_stream = None)
703
704 def doInitialization(self):
705 """
706 Initializes the ocean region by reading coordinate files
707 """
708 # We presume that %%'s in the source name requires calling a command
709 # otherwise it is a filename containing the required data
710 if self.coastline_source.find("%%") != -1:
711 c=self.__mergeParameters(self.coastline_source)
712 self.coastline_stream=os.popen(c).read()
713 else:
714 self.coastline_stream=self.coastline_source
715
716 if self.bathymetry_source.find("%%") != -1:
717 b=self.__mergeParameters(self.bathymetry_source)
718 self.bathymetry_stream=os.popen(b).read()
719 else:
720 self.bathymetry_stream=self.bathymetry_source
721
722 def __mergeParameters(self,txt):
723 return txt.replace("%%west%%", str(self.west)) \
724 .replace("%%east%%", str(self.east)) \
725 .replace("%%south%%", str(self.south)) \
726 .replace("%%north%%", str(self.north)) \
727 .replace("%%resolution%%", str(self.resolution)) \
728 .replace("%%range360%%", str(self.range360).lower())
729
730
731 #=============================================================================
732 class Bathymetry(Model):
733 """
734 Generates the bathymetry data within a region on the earth
735 """
736 def __init__(self,**kwargs):
737 Model.__init__(self,**kwargs)
738 self.declareParameter(source = "none",
739 bathymetry = 1.)
740
741 def doInitialization(self):
742 """
743 Initializes the bathymetry grid
744 """
745
746 print "Intializing bathymetry..."
747 if hasattr(self.source,"readline"):
748 f=self.source
749 else:
750 f=open(self.source,"r")
751 x_grd_list=[]
752 y_grd_list=[]
753 data_grd_list=[]
754 line=f.readline().strip()
755 while line != "":
756 v=line.split()
757 x_grd_list.append(float(v[0]))
758 y_grd_list.append(float(v[1]))
759 data_grd_list.append(float(v[2]))
760 line=f.readline().strip()
761 self.trace("%s lines have been read from %s." % (len(data_grd_list), self.source))
762 data_grd=numarray.array(data_grd_list)
763 x_grd=numarray.array(x_grd_list)
764 y_grd=numarray.array(y_grd_list)
765 if len(x_grd)<2:
766 raise ValueError,"%s: data base is too small"%str(self)
767 ox=x_grd[0]
768 oy=y_grd[0]
769 diam=max(abs(x_grd[len(x_grd)-1]-ox),abs(y_grd[len(y_grd)-1]-oy))
770 dx=x_grd[1]-ox
771 nx=1
772 nx=1
773 while abs(x_grd[nx]-ox)>EPS*diam:
774 nx+=1
775 dy=y_grd[nx]-oy
776 ny=len(x_grd)/nx
777 data_grd.resize((ny,nx))
778 self.bathymetry=GridData(s=[dx,dy],o=[ox,oy],n=[nx,ny],data=data_grd)
779 self.trace("%sx%s grid with %sx%s spacing." % (nx,ny,dx,dy))
780
781
782 #=============================================================================
783 class OceanRegion(Model):
784 """
785 Generates the ocean region with a coast line and a bathymetry
786 """
787 def __init__(self,**kwargs):
788 Model.__init__(self,**kwargs)
789 self.declareParameter(domain = None,
790 resolution = 1.,
791 south = 0.,
792 north = 10.,
793 east = 0.,
794 west = 20.,
795 bathymetry = None,
796 bathymetry_data = None,
797 coastline = None,
798 source = "none")
799
800 def doInitialization(self):
801 """
802 Initializes the ocean region
803 """
804
805 # create output directory if necessary
806 if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)
807
808 print "Initializing ocean region..."
809 if hasattr(self.source,"readline"):
810 f=self.source
811 data_name=f.geturl()
812 else:
813 f=open(self.source,"r")
814 data_name=self.source
815
816 segs=[]
817 name=""
818 line=f.readline()
819 while not line=="":
820 line=line.strip()
821 if line[:7]=="#region":
822 data=line[line.index("[")+1:line.index("]")].split(",")
823 reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
824 self.coastline=Coastline(region=reg,name=data_name)
825 elif line.find("Shore Bin")>-1:
826 self.coastline.append(Polyline(segs,name))
827 name=line[2:]
828 segs=[]
829 if not (line=="" or line[0]=="#" or line[0]==">") :
830 x=line.split()
831 segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
832 line=f.readline()
833 self.coastline.append(Polyline(segs,name))
834 d=self.bathymetry_data.interpolate([[self.east,self.south],
835 [self.west,self.south],
836 [self.east,self.north],
837 [self.west,self.north]])
838 self.domain = self.coastline.makeTriangulation(
839 east_south_is_water=d[0]<=0,
840 west_south_is_water=d[1]<=0,
841 east_north_is_water=d[2]<=0,
842 west_north_is_water=d[3]<=0
843 ).getFinleyDomain()
844 self.domain.dump(os.path.join(WORKDIR, "coastline.nc"))
845 self.bathymetry = maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
846
847
848 #=============================================================================
849 class TsunamiSource(Model):
850 """
851 Defines a wave in Gaussian form between start and end.
852 """
853 GAMMA=0.05
854 def __init__(self,**kwargs):
855 Model.__init__(self,**kwargs)
856 self.declareParameter(domain = None,
857 start_lat = -10.,
858 start_long = 110.,
859 end_lat = -12.,
860 end_long = 120.,
861 width = 5.,
862 decay_zone = 0.1,
863 amplitude = 1.,
864 wave_height = 1.)
865
866 def doInitialization(self):
867 """
868 Sets the initial wave form
869 """
870 beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
871 x=self.domain.getX()
872 x_long = x[0]
873 x_lat = x[1]
874 mid_long = (self.start_long+self.end_long)/2
875 mid_lat = (self.start_lat+self.end_lat)/2
876 dist = math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
877 a = (self.start_long-mid_long)/dist
878 b = (self.start_lat-mid_lat)/dist
879 self.trace("source length = %s" % (dist*2.))
880 x_long_hat = a*(x_long-mid_long)+b*(x_lat-mid_lat)
881 x_lat_hat = -b*(x_long-mid_long)+a*(x_lat-mid_lat)
882 # x_lat-direction
883 s = abs(x_lat_hat)-self.width
884 m = whereNegative(s)
885 f1 = (1.-m)*exp(-(s*beta)**2)+m
886
887 # x_long-direction
888 s = abs(x_long_hat)-dist
889 m = whereNegative(s)
890 f0 = (1.-m)*exp(-(s*beta)**2)+m
891 self.wave_height = f1*f0*self.amplitude
892
893
894 #=============================================================================
895 class TsunamiInDeepWater(Model):
896 """
897 Runs the deep water tsunami model based on a simplified form of the
898 shallow water equation.
899
900 M{d^2 h/dt^2 =div(c grad(h)) }
901
902 where h is the wave height above sea level, and c=sqrt(g*H),
903 with H - depth of the water level, g - gravity constant
904
905 The simulation uses the Verlet scheme.
906 """
907 def __init__(self,**kwargs):
908 Model.__init__(self,**kwargs)
909 self.declareParameter(domain = None,
910 wave_height = 1.,
911 wave_velocity = 0.,
912 initial_time_step = None,
913 bathymetry = 1.,
914 safety_factor = 1.)
915
916 def doInitialization(self):
917 """
918 Initializes the time integration scheme
919 """
920
921 self.__pde = LinearPDE(self.domain)
922 self.__pde.setSolverMethod(self.__pde.LUMPING)
923 self.__pde.setValue(D=1.)
924 self.__c2 = RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2
925 c_max = math.sqrt(Lsup(self.__c2))
926 self.__dt = self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
927 if self.initial_time_step==None:
928 self.initial_time_step=self.__dt
929 self.trace("Maximum wave velocity %s m/sec" % c_max)
930 self.trace("Time step size is %s sec" % self.__dt)
931
932 def getSafeTimeStepSize(self, dt):
933 """
934 Returns new step size
935
936 @param dt: last time step size used
937 @type dt: C{float}
938 @return: time step size that can safely be used
939 @rtype: C{float}
940 """
941 return self.__dt
942
943 def doStepPostprocessing(self,dt):
944 """
945 Performs the time step using the Verlet scheme
946
947 @param dt: time step size to be used
948 @type dt: C{float}
949 """
950 self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
951
952 new_height = self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
953
954 self.wave_velocity = (new_height-self.wave_height)/dt
955 self.wave_height = new_height
956 self.initial_time_step = dt
957 self.trace("Wave height range is %e %e" % (inf(self.wave_height), sup(self.wave_height)))
958
959
960 #=============================================================================
961 class SurfMovie(Model):
962 """
963 movie from a wave propagation on the sea
964
965 @ivar time: current time
966 @ivar bathymetry: scalar data set
967 @ivar wave_height: vector data set
968 @ivar filename: name of the movie file
969 """
970 def __init__(self,**kwargs):
971 Model.__init__(self,**kwargs)
972 self.declareParameter(bathymetry = 1.,
973 wave_height = 1.,
974 coastline = None,
975 t = 0.,
976 dt = 1.,
977 south = 2.,
978 north = 5.,
979 east = 3.,
980 west = 15.,
981 max_height = 1.,
982 filename = "tsunami.mpg")
983
984 def doInitialization(self):
985 """
986 Initializes the time integration scheme
987 """
988 self.__frame_name="tsunami"
989 self.__next_t=self.dt
990 self.__fileIndex=0
991
992 self.firstFrame = True
993 self.imageFiles = []
994
995 factGraphics = vtk.vtkGraphicsFactory()
996 factGraphics.SetUseMesaClasses(1)
997 factImage = vtk.vtkImagingFactory()
998 factImage.SetUseMesaClasses(1)
999
1000 print "Preparing bathymetry data..."
1001 # get depth (z), origin, dx and dy of the bathymetry
1002 bathZData = self.bathymetry.getData()
1003 bathOrigin = self.bathymetry.getOrigin()
1004 bathSpacing = self.bathymetry.getSpacing()
1005
1006 # now construct the x and y data from the spacing and the origin
1007 numXPoints = bathZData.shape[1]
1008 numYPoints = bathZData.shape[0]
1009 numPoints = numXPoints*numYPoints
1010
1011 bathXData = numarray.zeros(numXPoints, numarray.Float)
1012 bathYData = numarray.zeros(numYPoints, numarray.Float)
1013 for i in range(numXPoints):
1014 bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
1015
1016 for i in range(numYPoints):
1017 bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
1018
1019 # prepare the bathymetry data set
1020 bathPoints = vtk.vtkPoints()
1021 bathPoints.SetNumberOfPoints(numPoints)
1022 bathData = vtk.vtkFloatArray()
1023 bathData.SetNumberOfComponents(1)
1024 bathData.SetNumberOfTuples(numPoints)
1025
1026 # add the points and data values
1027 count = 0
1028 for i in range(numXPoints):
1029 for j in range(numYPoints):
1030 bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1031 bathData.InsertTuple1(count, bathZData[j,i])
1032 count += 1
1033
1034 # create an unstructured grid for the bathymetry
1035 bathGrid = vtk.vtkUnstructuredGrid()
1036 bathGrid.SetPoints(bathPoints)
1037 bathGrid.GetPointData().SetScalars(bathData)
1038
1039 # do a delaunay on the grid
1040 bathDelaunay = vtk.vtkDelaunay2D()
1041 bathDelaunay.SetInput(bathGrid)
1042 bathDelaunay.SetTolerance(0.001)
1043 bathDelaunay.SetAlpha(2.5)
1044
1045 # save static bathymetry file
1046 writer = vtk.vtkPolyDataWriter()
1047 writer.SetInput(bathDelaunay.GetOutput())
1048 writer.SetScalarsName("bathymetry")
1049 writer.SetFileName(os.path.join(WORKDIR, "bathymetry.vtk"))
1050 writer.Write()
1051
1052 # create the bathymetry colourmap
1053 data = []
1054 data.append([-8000, 0, 0, 0])
1055 data.append([-7000, 0, 5, 25])
1056 data.append([-6000, 0, 10, 50])
1057 data.append([-5000, 0, 80, 125])
1058 data.append([-4000, 0, 150, 200])
1059 data.append([-3000, 86, 197, 184])
1060 data.append([-2000, 172, 245, 168])
1061 data.append([-1000, 211, 250, 211])
1062 data.append([0, 16, 48, 16])
1063 data.append([300, 70, 150, 50])
1064 data.append([500, 146, 126, 60])
1065 data.append([1000, 198, 178, 80])
1066 data.append([1250, 250, 230, 100])
1067 data.append([1500, 250, 234, 126])
1068 data.append([1750, 252, 238, 152])
1069 data.append([2000, 252, 243, 177])
1070 data.append([2250, 253, 249, 216])
1071 data.append([2500, 255, 255, 255])
1072
1073 # the amount to scale the data by
1074 scale = 255.0
1075 numColours = len(data)
1076
1077 # convert the colourmap into something vtk is more happy with
1078 height = numarray.zeros(numColours, numarray.Float)
1079 red = numarray.zeros(numColours, numarray.Float)
1080 green = numarray.zeros(numColours, numarray.Float)
1081 blue = numarray.zeros(numColours, numarray.Float)
1082 for i in range(numColours):
1083 (h, r, g, b) = data[i]
1084 height[i] = float(h)
1085 red[i] = float(r)/scale
1086 green[i] = float(g)/scale
1087 blue[i] = float(b)/scale
1088
1089 transFunc = vtk.vtkColorTransferFunction()
1090 transFunc.SetColorSpaceToRGB()
1091 for i in range(numColours):
1092 transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1093 h = height[i]
1094 while i < numColours-1 and h < height[i+1]:
1095 h += 1
1096 transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1097
1098 # set up the mapper
1099 bathMapper = vtk.vtkDataSetMapper()
1100 bathMapper.SetInput(bathDelaunay.GetOutput())
1101 bathMapper.SetLookupTable(transFunc)
1102
1103 # set up the actor
1104 self.bathActor = vtk.vtkActor()
1105 self.bathActor.SetMapper(bathMapper)
1106
1107 print "Done preparing bathymetry data"
1108
1109 ### prepare and add the coastline ###
1110
1111 print "Preparing coastline data..."
1112 # create the coastline grid
1113 coastGrid = vtk.vtkUnstructuredGrid()
1114 coastPoints = vtk.vtkPoints()
1115 totalCoastPoints = 0
1116 for polyline in self.coastline.polylines:
1117 numPoints = len(polyline)
1118 coastLine = vtk.vtkPolyLine()
1119 coastLine.GetPointIds().SetNumberOfIds(numPoints)
1120 j = 0
1121 for point in polyline:
1122 coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1123 coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1124 j += 1
1125 coastGrid.InsertNextCell(coastLine.GetCellType(),
1126 coastLine.GetPointIds())
1127 totalCoastPoints += numPoints
1128
1129 coastGrid.SetPoints(coastPoints)
1130 coastMapper = vtk.vtkDataSetMapper()
1131 coastMapper.SetInput(coastGrid)
1132
1133 # create an actor for the coastline
1134 self.coastActor = vtk.vtkActor()
1135 self.coastActor.SetMapper(coastMapper)
1136 self.coastActor.GetProperty().SetColor(0,0,0)
1137
1138 print "Done preparing coastline data"
1139
1140 # set up the lookup table for the wave data
1141 refLut = vtk.vtkLookupTable()
1142 self.lutTrans = vtk.vtkLookupTable()
1143 refLut.Build()
1144 alpha = 0.7 # alpha channel value
1145 for i in range(256):
1146 (r,g,b,a) = refLut.GetTableValue(255-i)
1147 if g == 1.0 and b < 0.5 and r < 0.5:
1148 self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1149 else:
1150 self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1151
1152 # create an actor for the wave
1153 self.waveActor = vtk.vtkActor()
1154
1155 # calculate window size
1156 dataWidth = max(bathXData) - min(bathXData)
1157 dataHeight = max(bathYData) - min(bathYData)
1158 self._winHeight = 600
1159 self._winWidth = int(dataWidth*float(self._winHeight)/dataHeight)
1160
1161 def saveImage(self, prefix):
1162 """
1163 Saves current state as PNM image by saving a VTK file with the
1164 wave height first and reading it back afterwards.
1165 """
1166
1167 # do stuff here that should only be done once
1168 if self.firstFrame:
1169 ### set up the renderer and the render window ###
1170 self.ren = vtk.vtkRenderer()
1171 self.renWin = vtk.vtkRenderWindow()
1172 self.renWin.AddRenderer(self.ren)
1173 self.renWin.SetSize(self._winWidth, self._winHeight)
1174 # add actors to the renderer
1175 self.ren.AddActor(self.bathActor)
1176 self.ren.AddActor(self.coastActor)
1177 self.ren.AddActor(self.waveActor)
1178 self.ren.GetActiveCamera().Zoom(2.0)
1179
1180 # use interactor...
1181 self.iren = vtk.vtkRenderWindowInteractor()
1182 self.iren.SetRenderWindow(self.renWin)
1183 self.iren.Initialize()
1184 self.iren.Start()
1185
1186 # ...or make sure rendering is offscreen
1187 #self.renWin.OffScreenRenderingOn()
1188
1189 self.firstFrame = False
1190
1191 vtuFile = prefix+".vtu"
1192 imgFile = prefix+".pnm"
1193
1194 # save the VTK file first
1195 print "Writing", vtuFile
1196 saveVTK(vtuFile, h=self.wave_height)
1197
1198 # make a reader for the data
1199 waveReader = vtk.vtkXMLUnstructuredGridReader()
1200 waveReader.SetFileName(vtuFile)
1201 waveReader.Update()
1202
1203 waveGrid = waveReader.GetOutput()
1204 waveGrid.Update()
1205
1206 waveMapper = vtk.vtkDataSetMapper()
1207 waveMapper.SetInput(waveGrid)
1208 waveMapper.SetLookupTable(self.lutTrans)
1209 waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1210
1211 self.waveActor.SetMapper(waveMapper)
1212
1213 # now render the window and save the image file
1214 print "Rendering to", imgFile
1215 self.renWin.Render()
1216
1217 # convert the render window to an image
1218 win2imgFilter = vtk.vtkWindowToImageFilter()
1219 win2imgFilter.SetInput(self.renWin)
1220
1221 # save the image to file
1222 outWriter = vtk.vtkPNMWriter()
1223 outWriter.SetInput(win2imgFilter.GetOutput())
1224 outWriter.SetFileName(imgFile)
1225 outWriter.Write()
1226 self.imageFiles.append(imgFile)
1227
1228 def doStepPostprocessing(self, dt):
1229 """
1230 Does any necessary postprocessing after each step
1231
1232 @param dt:
1233 """
1234 if self.t >= self.__next_t:
1235 prefix = os.path.join(WORKDIR, \
1236 "%s%04d" % (self.__frame_name, self.__fileIndex))
1237
1238 zMin = inf(self.wave_height)
1239 zMax = sup(self.wave_height)
1240 print "T = %7.1f: Wave height range=%f...%f" % (self.t, zMin, zMax)
1241
1242 #self.saveImage(prefix)
1243 self.wave_height.dump(os.path.join(WORKDIR, \
1244 "waveheight%04d.nc" % self.__fileIndex))
1245
1246 self.__next_t += self.dt
1247 self.__fileIndex += 1
1248
1249 def getSafeTimeStepSize(self,dt):
1250 """
1251 returns new step size
1252
1253 @param dt: last time step size used
1254 @type dt: C{float}
1255 @return: time step size that can savely be used
1256 @rtype: C{float}
1257 """
1258 return self.__next_t-self.t
1259
1260 def doFinalization(self):
1261 """
1262 Finalises the visualisation. For instance, makes a movie of the
1263 image files.
1264 """
1265
1266 if len(self.imageFiles) == 0:
1267 return
1268
1269 # set up the movie parameters
1270 paramsFileString = "REFERENCE_FRAME DECODED\n"
1271 paramsFileString += "FRAME_RATE 24\n"
1272 paramsFileString += "OUTPUT %s\n" % self.filename
1273 paramsFileString += "PATTERN IBBPBBPBB\n"
1274 paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
1275 paramsFileString += "GOP_SIZE 20\n"
1276 paramsFileString += "BSEARCH_ALG CROSS2\n"
1277 paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
1278 paramsFileString += "IQSCALE 10\n"
1279 paramsFileString += "PQSCALE 11\n"
1280 paramsFileString += "BQSCALE 16\n"
1281 paramsFileString += "RANGE 8\n"
1282 paramsFileString += "SLICES_PER_FRAME 1\n"
1283 paramsFileString += "BASE_FILE_FORMAT PNM\n"
1284 paramsFileString += "INPUT_DIR .\n"
1285 paramsFileString += "INPUT_CONVERT *\n"
1286 paramsFileString += "INPUT\n"
1287 for fname in self.imageFiles:
1288 paramsFileString += "%s\n" % fname
1289 paramsFileString += "END_INPUT\n"
1290 paramsFileString += "PIXEL HALF\n"
1291 paramsFileString += "ASPECT_RATIO 1\n"
1292
1293 # write the parameter file
1294 fp = open("%s.params" % self.filename, "w")
1295 fp.write(paramsFileString)
1296 fp.close()
1297
1298 print "Performing conversion to mpeg"
1299 convertString = "ppmtompeg %s.params" % self.filename
1300 # write the movie into self.filename
1301 result = os.system(convertString)
1302 if result != 0:
1303 print "An error occurred in mpeg conversion"
1304
1305 # now clean up the image files
1306 if result == 0:
1307 print "Removing temporary image files"
1308 os.unlink("%s.params" % self.filename)
1309 for fname in self.imageFiles:
1310 os.unlink(fname)
1311
1312
1313 def main():
1314 from esys.escript.modelframe import Link,Simulation
1315 from esys.modellib.input import Sequencer
1316
1317 oc = OceanRegionCollector()
1318 oc.range360 = True
1319 oc.west = 102.9
1320 oc.east = 232.6
1321 oc.south = -71.3
1322 oc.north = 26.7
1323 oc.resolution = 0.25
1324
1325 # Uncomment and use the following if GMT is available and the scripts work
1326 #oc.coastline_source="tsunami_create_coast.py 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 earthquake=false"
1327 #oc.bathymetry_source = "tsunami_create_bath.py west=%%west%% east=%%east%% south=%%south%% north=%%north%% resolution=%%resolution%%"
1328
1329 # ...otherwise use pre-created data files (range parameters are ignored)
1330 oc.coastline_source = "coastline_default.dat"
1331 oc.bathymetry_source = "bathymetry_default.dat"
1332
1333 b = Bathymetry()
1334 b.source = Link(oc,"bathymetry_stream")
1335
1336 oreg = OceanRegion()
1337 oreg.source = Link(oc,"coastline_stream")
1338 oreg.resolution = Link(oc, "resolution")
1339 oreg.south = Link(oc, "south")
1340 oreg.north = Link(oc, "north")
1341 oreg.east = Link(oc, "east")
1342 oreg.west = Link(oc, "west")
1343 oreg.bathymetry_data = Link(b, "bathymetry")
1344
1345 src = TsunamiSource()
1346 src.domain = Link(oreg, "domain")
1347 src.decay_zone = 0.01
1348 src.end_long = 185.
1349 src.end_lat = -37.
1350 src.width = 0.5
1351 src.start_long = 174.
1352 src.start_lat = -15.
1353 src.amplitude = 3
1354
1355 ts = TsunamiInDeepWater()
1356 ts.domain = Link(oreg, "domain")
1357 ts.wave_height = Link(src, "wave_height")
1358 ts.wave_velocity = 0.
1359 ts.bathymetry = Link(oreg, "bathymetry")
1360
1361 sq = Sequencer()
1362 sq.t_end = 30000.
1363
1364 sm = SurfMovie()
1365 sm.bathymetry = Link(b, "bathymetry")
1366 sm.wave_height = Link(ts, "wave_height")
1367 sm.coastline = Link(oreg, "coastline")
1368 sm.t = Link(sq,"t")
1369 sm.filename = "tsunami.mpg"
1370 sm.north = 8.7
1371 sm.west = 138.9
1372 sm.dt = 50.
1373 sm.east = 196.6
1374 sm.south = -53.3
1375 sm.max_height = Link(src, "amplitude")
1376
1377 s = Simulation([sq,oc,b,oreg,src,ts,sm])
1378 #s.writeXML()
1379 s.run()
1380
1381 if __name__=="__main__":
1382 #from esys.modellib import tsunami
1383 main()
1384
1385 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26