/[escript]/branches/3.4.1/modellib/py_src/tsunami.py.old
ViewVC logotype

Contents of /branches/3.4.1/modellib/py_src/tsunami.py.old

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4587 - (show annotations)
Wed Dec 11 06:17:09 2013 UTC (5 years, 3 months ago) by jfenwick
File size: 50622 byte(s)
Preparation begins

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

  ViewVC Help
Powered by ViewVC 1.1.26