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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 322 by elspeth, Wed Nov 30 09:10:59 2005 UTC revision 323 by gross, Tue Dec 6 06:18:00 2005 UTC
# Line 12  import math Line 12  import math
12  #import urllib  #import urllib
13  import urllib2  import urllib2
14    
15  EPS=1.e-5  EPS=1.e-8
16    
17  #====================================================================================================================================================  #====================================================================================================================================================
18  class GridData:  class GridData:
# Line 354  class Coastline: Line 354  class Coastline:
354               if pl.getDiameter()>self.region.resolution:               if pl.getDiameter()>self.region.resolution:
355                  short_pl=[pl[0]]                  short_pl=[pl[0]]
356                  for i in range(1,len(pl)):                  for i in range(1,len(pl)):
357                        if pl[i]-short_pl[-1]>self.region.resolution/5:                        if pl[i]-short_pl[-1]>0*EPS+self.region.resolution/10:
358                           short_pl.append(pl[i])                           short_pl.append(pl[i])
359                        elif i==len(pl)-1:                        elif i==len(pl)-1:
360                           short_pl[-1]=pl[i]                           short_pl[-1]=pl[i]
# Line 422  class Coastline: Line 422  class Coastline:
422                  vertices_on_face.append(q)                  vertices_on_face.append(q)
423           vertices_on_face.sort(self.region.comparePointsOnFace)           vertices_on_face.sort(self.region.comparePointsOnFace)
424           index=0           index=0
425           walking_on_water=west_south_is_water           walking_on_water=True
426           l=len(vertices_on_face)           l=len(vertices_on_face)
427           while index<l:           while index<l:
428               p1=vertices_on_face[(index+1)%l]               p1=vertices_on_face[(index+1)%l]
# Line 610  class EarthTriangulation: Line 610  class EarthTriangulation:
610             finley_file.close()             finley_file.close()
611             # get the mesh             # get the mesh
612             out=ReadMesh("%s.msh"%self.fn)             out=ReadMesh("%s.msh"%self.fn)
613             # os.remove("%s.msh"%self.fn)             os.remove("%s.msh"%self.fn)
614             return out             return out
615    
616    
# Line 692  class Bathymetry(Model): Line 692  class Bathymetry(Model):
692             dx=x_grd[1]-ox             dx=x_grd[1]-ox
693             nx=1             nx=1
694             nx=1             nx=1
695             while abs(x_grd[nx]-ox)>1.e-10*diam:             while abs(x_grd[nx]-ox)>EPS*diam:
696                 nx+=1                 nx+=1
697             dy=y_grd[nx]-oy             dy=y_grd[nx]-oy
698             ny=len(x_grd)/nx             ny=len(x_grd)/nx
# Line 789  class TsunamiSource(Model): Line 789  class TsunamiSource(Model):
789             x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat)             x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat)
790             # x_lat-direction             # x_lat-direction
791             s=abs(x_lat_hat)-self.width             s=abs(x_lat_hat)-self.width
792             m=s.whereNegative()             m=whereNegative(s)
793             f1=(1.-m)*exp(-(s*beta)**2)+m             f1=(1.-m)*exp(-(s*beta)**2)+m
794             # x_long-direction             # x_long-direction
795             s=abs(x_long_hat)-dist             s=abs(x_long_hat)-dist
796             m=s.whereNegative()             m=whereNegative(s)
797             f0=(1.-m)*exp(-(s*beta)**2)+m             f0=(1.-m)*exp(-(s*beta)**2)+m
798             self.wave_height=f1*f0*self.amplitude             self.wave_height=f1*f0*self.amplitude
799    
# Line 827  class TsunamiInDeepWater(Model): Line 827  class TsunamiInDeepWater(Model):
827             self.__pde=LinearPDE(self.domain)             self.__pde=LinearPDE(self.domain)
828             self.__pde.setSolverMethod(self.__pde.LUMPING)             self.__pde.setSolverMethod(self.__pde.LUMPING)
829             self.__pde.setValue(D=1.)             self.__pde.setValue(D=1.)
830             self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/RegionOnEarthSurface.RADIUS**2             self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2
831             c_max=math.sqrt(Lsup(self.__c2))             c_max=math.sqrt(Lsup(self.__c2))
832             self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))             self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
833               print c_max,inf(self.domain.getSize())
834             if self.initial_time_step==None: self.initial_time_step=self.__dt             if self.initial_time_step==None: self.initial_time_step=self.__dt
835             self.trace("maximum wave velocity %s m/sec"%c_max)             self.trace("maximum wave velocity %s m/sec"%c_max)
836             self.trace("Time step size is %s sec"%self.__dt)             self.trace("Time step size is %s sec"%self.__dt)
# Line 1119  class SurfMovie(Model): Line 1120  class SurfMovie(Model):
1120          @param dt:          @param dt:
1121          """          """
1122          if self.t>=self.__next_t:          if self.t>=self.__next_t:
1123               print self.t,"write ",Lsup(self.wave_height)               print self.t,"write ",Lsup(self.wave_height)," to ",self.__fn
1124               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1125               # vtkobj=...               # vtkobj=...
1126               # save(self.__frame_name)               # save(self.__frame_name)
# Line 1140  class SurfMovie(Model): Line 1141  class SurfMovie(Model):
1141               waveMapper = vtk.vtkDataSetMapper()               waveMapper = vtk.vtkDataSetMapper()
1142               waveMapper.SetInput(waveGrid)               waveMapper.SetInput(waveGrid)
1143               waveMapper.SetLookupTable(self.lutTrans)               waveMapper.SetLookupTable(self.lutTrans)
1144               waveMapper.SetScalarRange(-self.max_height, self.max_height)               waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1145    
1146               self.waveActor.SetMapper(waveMapper)               self.waveActor.SetMapper(waveMapper)
1147    
# Line 1218  def main(): Line 1219  def main():
1219     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
1220    
1221     oc=OceanRegionCollector()     oc=OceanRegionCollector()
1222     oc.resolution=2     oc.north= 26.7
1223     oc.south=-45.5     oc.west= 102.9
1224     oc.north=-5.4     oc.range360= True
1225     oc.east=161.1     oc.east= 232.6
1226     oc.west=108.2     oc.resolution= 1.
1227     oc.range360=True     oc.south= -71.3
1228    
1229    
1230     b=Bathymetry()     b=Bathymetry()
# Line 1240  def main(): Line 1241  def main():
1241    
1242     src=TsunamiSource()     src=TsunamiSource()
1243     src.domain=Link(oreg,"domain")     src.domain=Link(oreg,"domain")
1244     src.start_lat=-10.     src.decay_zone= 0.01
1245     src.end_lat=-12.     src.end_long= 185.
1246     src.start_long=110.     src.end_lat= -37.
1247     src.end_long=120.     src.width=0.5
1248     src.width=1.     src.start_long= 174.
1249     src.decay_zone=0.01     src.start_lat= -15.
1250     src.amplitude=1.     src.amplitude= 3
1251    
1252     ts=TsunamiInDeepWater()     ts=TsunamiInDeepWater()
1253     ts.domain=Link(oreg,"domain")     ts.domain=Link(oreg,"domain")
# Line 1255  def main(): Line 1256  def main():
1256     ts.bathymetry=Link(oreg,"bathymetry")     ts.bathymetry=Link(oreg,"bathymetry")
1257    
1258     sq=Sequencer()     sq=Sequencer()
1259     sq.t_end=100000.     sq.t_end=15000.
1260    
1261     sm=SurfMovie()     sm=SurfMovie()
1262     sm.bathymetry=Link(b,"bathymetry")     sm.bathymetry=Link(b,"bathymetry")
1263     sm.wave_height=Link(ts,"wave_height")     sm.wave_height=Link(ts,"wave_height")
1264     sm.coastline=Link(oreg,"coastline")     sm.coastline=Link(oreg,"coastline")
1265     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
    sm.dt=5000.  
1266     sm.filename="movie.mpg"     sm.filename="movie.mpg"
1267       sm.north= 8.7
1268       sm.west= 138.9
1269       sm.dt= 50.
1270       sm.east= 196.6
1271       sm.south= -53.3
1272     sm.max_height=Link(src,"amplitude")     sm.max_height=Link(src,"amplitude")
1273        
1274     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
Legend:
Removed from v.322  
changed lines
  Added in v.323
  ViewVC Help
Powered by ViewVC 1.1.26