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

  ViewVC Help
Powered by ViewVC 1.1.26