/[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 280 by cochrane, Wed Nov 30 06:23:59 2005 UTC revision 918 by gross, Wed Jan 3 06:30:00 2007 UTC
# Line 2  Line 2 
2    
3  # $Id$  # $Id$
4    
5    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
6                        http://www.access.edu.au
7                    Primary Business: Queensland, Australia"""
8    __license__="""Licensed under the Open Software License version 3.0
9                 http://www.opensource.org/licenses/osl-3.0.php"""
10    
11  import os, sys  import os, sys
12  import vtk  import vtk
13  from esys.escript import *  from esys.escript import *
# Line 9  from esys.escript.linearPDEs import Line Line 15  from esys.escript.linearPDEs import Line
15  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
16  import numarray  import numarray
17  import math  import math
18  import urllib  #import urllib
19    import urllib2
20    
21  EPS=1.e-5  EPS=1.e-8
22    
23  #====================================================================================================================================================  #====================================================================================================================================================
24  class GridData:  class GridData:
# Line 353  class Coastline: Line 360  class Coastline:
360               if pl.getDiameter()>self.region.resolution:               if pl.getDiameter()>self.region.resolution:
361                  short_pl=[pl[0]]                  short_pl=[pl[0]]
362                  for i in range(1,len(pl)):                  for i in range(1,len(pl)):
363                        if pl[i]-short_pl[-1]>self.region.resolution/5:                        if pl[i]-short_pl[-1]>0*EPS+self.region.resolution/10:
364                           short_pl.append(pl[i])                           short_pl.append(pl[i])
365                        elif i==len(pl)-1:                        elif i==len(pl)-1:
366                           short_pl[-1]=pl[i]                           short_pl[-1]=pl[i]
# Line 421  class Coastline: Line 428  class Coastline:
428                  vertices_on_face.append(q)                  vertices_on_face.append(q)
429           vertices_on_face.sort(self.region.comparePointsOnFace)           vertices_on_face.sort(self.region.comparePointsOnFace)
430           index=0           index=0
431           walking_on_water=west_south_is_water           walking_on_water=True
432           l=len(vertices_on_face)           l=len(vertices_on_face)
433           while index<l:           while index<l:
434               p1=vertices_on_face[(index+1)%l]               p1=vertices_on_face[(index+1)%l]
# Line 609  class EarthTriangulation: Line 616  class EarthTriangulation:
616             finley_file.close()             finley_file.close()
617             # get the mesh             # get the mesh
618             out=ReadMesh("%s.msh"%self.fn)             out=ReadMesh("%s.msh"%self.fn)
619             # os.remove("%s.msh"%self.fn)             os.remove("%s.msh"%self.fn)
620             return out             return out
621    
622    
# Line 621  class OceanRegionCollector(Model): Line 628  class OceanRegionCollector(Model):
628    
629    
630        """        """
631        def __init__(self,debug=False):        def __init__(self,**kwargs):
632             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
633             self.declareParameter(coastline_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/datafile.txt?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%&range360=%%range360%%&river=false&city=false&citytyp=wcity.dat&volcano=false&hotspot=false&shoreline=true&state=false&WSMdata=false&plate=false",             self.declareParameter(coastline_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/datafile.txt?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%&range360=%%range360%%&river=false&city=false&citytyp=wcity.dat&volcano=false&hotspot=false&shoreline=true&state=false&WSMdata=false&plate=false",
634                                   bathymetry_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/grdfile.xyz?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%",                                   bathymetry_source="http://jamboree.esscc.uq.edu.au/cgi-bin/doreen/grdfile.xyz?west=%%west%%&east=%%east%%&south=%%south%%&north=%%north%%&resolution=%%resolution%%",
635                                   resolution=1.,                                   resolution=1.,
# Line 641  class OceanRegionCollector(Model): Line 648  class OceanRegionCollector(Model):
648             """             """
649             c=self.__mergeParameters(self.coastline_source)             c=self.__mergeParameters(self.coastline_source)
650             b=self.__mergeParameters(self.bathymetry_source)             b=self.__mergeParameters(self.bathymetry_source)
651             self.coastline_stream=urllib.urlopen(c)             self.coastline_stream=urllib2.urlopen(c)
652             self.bathymetry_stream=urllib.urlopen(b)             self.bathymetry_stream=urllib2.urlopen(b)
653    
654        def __mergeParameters(self,txt):        def __mergeParameters(self,txt):
655             return txt.replace("%%west%%",str(self.west))\             return txt.replace("%%west%%",str(self.west))\
# Line 656  class Bathymetry(Model): Line 663  class Bathymetry(Model):
663         """         """
664         generates the bathymetry data within a region on the earth         generates the bathymetry data within a region on the earth
665         """         """
666         def __init__(self,debug=False):         def __init__(self,**kwargs):
667             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
668             self.declareParameter(source="none",             self.declareParameter(source="none",
669                                   bathymetry=1.)                                   bathymetry=1.)
670    
# Line 691  class Bathymetry(Model): Line 698  class Bathymetry(Model):
698             dx=x_grd[1]-ox             dx=x_grd[1]-ox
699             nx=1             nx=1
700             nx=1             nx=1
701             while abs(x_grd[nx]-ox)>1.e-10*diam:             while abs(x_grd[nx]-ox)>EPS*diam:
702                 nx+=1                 nx+=1
703             dy=y_grd[nx]-oy             dy=y_grd[nx]-oy
704             ny=len(x_grd)/nx             ny=len(x_grd)/nx
# Line 705  class OceanRegion(Model): Line 712  class OceanRegion(Model):
712         generates the ocean region with a coast line and a bathymetry         generates the ocean region with a coast line and a bathymetry
713    
714         """         """
715         def __init__(self,debug=False):         def __init__(self,**kwargs):
716             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
717             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
718                                   resolution=1.,                                   resolution=1.,
719                                   south=0.,                                   south=0.,
# Line 744  class OceanRegion(Model): Line 751  class OceanRegion(Model):
751                     name=line[2:]                     name=line[2:]
752                     segs=[]                     segs=[]
753                if not (line=="" or line[0]=="#" or line[0]==">") :                if not (line=="" or line[0]=="#" or line[0]==">") :
754                      print line
755                    x=line.split()                    x=line.split()
756                    segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))                    segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
757                line=f.readline()                line=f.readline()
# Line 758  class TsunamiSource(Model): Line 766  class TsunamiSource(Model):
766         defines a wave in Gaussean form between start and end.         defines a wave in Gaussean form between start and end.
767         """         """
768         GAMMA=0.05         GAMMA=0.05
769         def __init__(self,debug=False):         def __init__(self,**kwargs):
770             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
771             self.declareParameter(domain=None,             self.declareParameter(domain=None,
772                                   start_lat=-10.,                                   start_lat=-10.,
773                                   start_long=110.,                                   start_long=110.,
# Line 788  class TsunamiSource(Model): Line 796  class TsunamiSource(Model):
796             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)
797             # x_lat-direction             # x_lat-direction
798             s=abs(x_lat_hat)-self.width             s=abs(x_lat_hat)-self.width
799             m=s.whereNegative()             m=whereNegative(s)
800             f1=(1.-m)*exp(-(s*beta)**2)+m             f1=(1.-m)*exp(-(s*beta)**2)+m
801             # x_long-direction             # x_long-direction
802             s=abs(x_long_hat)-dist             s=abs(x_long_hat)-dist
803             m=s.whereNegative()             m=whereNegative(s)
804             f0=(1.-m)*exp(-(s*beta)**2)+m             f0=(1.-m)*exp(-(s*beta)**2)+m
805             self.wave_height=f1*f0*self.amplitude             self.wave_height=f1*f0*self.amplitude
806    
# Line 810  class TsunamiInDeepWater(Model): Line 818  class TsunamiInDeepWater(Model):
818         The simulation uses the Verlet scheme.         The simulation uses the Verlet scheme.
819    
820         """         """
821         def __init__(self,debug=False):         def __init__(self,**kwargs):
822             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
823             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
824                                   wave_height=1.,                                   wave_height=1.,
825                                   wave_velocity=0.,                                   wave_velocity=0.,
# Line 826  class TsunamiInDeepWater(Model): Line 834  class TsunamiInDeepWater(Model):
834             self.__pde=LinearPDE(self.domain)             self.__pde=LinearPDE(self.domain)
835             self.__pde.setSolverMethod(self.__pde.LUMPING)             self.__pde.setSolverMethod(self.__pde.LUMPING)
836             self.__pde.setValue(D=1.)             self.__pde.setValue(D=1.)
837             self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/RegionOnEarthSurface.RADIUS**2             self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2
838             c_max=math.sqrt(Lsup(self.__c2))             c_max=math.sqrt(Lsup(self.__c2))
839             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))
840               print c_max,inf(self.domain.getSize())
841             if self.initial_time_step==None: self.initial_time_step=self.__dt             if self.initial_time_step==None: self.initial_time_step=self.__dt
842             self.trace("maximum wave velocity %s m/sec"%c_max)             self.trace("maximum wave velocity %s m/sec"%c_max)
843             self.trace("Time step size is %s sec"%self.__dt)             self.trace("Time step size is %s sec"%self.__dt)
# Line 870  class SurfMovie(Model): Line 879  class SurfMovie(Model):
879         @ivar wave_height: vector data set         @ivar wave_height: vector data set
880         @ivar filename: name of the movie file         @ivar filename: name of the movie file
881         """         """
882         def __init__(self,debug=False):         def __init__(self,**kwargs):
883             Model.__init__(self,debug=debug)             Model.__init__(self,**kwargs)
884             self.declareParameter(bathymetry=1.,             self.declareParameter(bathymetry=1.,
885                                   wave_height=1.,                                   wave_height=1.,
886                                   coastline=None,                                   coastline=None,
# Line 884  class SurfMovie(Model): Line 893  class SurfMovie(Model):
893                                   max_height=1.0,                                   max_height=1.0,
894                                   filename="movie.mpg")                                   filename="movie.mpg")
895    
            self.paramsFileString = "REFERENCE_FRAME DECODED\n"  
            self.paramsFileString += "FRAME_RATE 24\n"  
            self.paramsFileString += "OUTPUT %s\n" % self.filename  
            self.paramsFileString += "PATTERN IBBPBBPBB\n"  
            self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"  
            self.paramsFileString += "GOP_SIZE 20\n"  
            self.paramsFileString += "BSEARCH_ALG CROSS2\n"  
            self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"  
            self.paramsFileString += "IQSCALE 10\n"  
            self.paramsFileString += "PQSCALE 11\n"  
            self.paramsFileString += "BQSCALE 16\n"  
            self.paramsFileString += "RANGE 8\n"  
            self.paramsFileString += "SLICES_PER_FRAME 1\n"  
            self.paramsFileString += "BASE_FILE_FORMAT PNM\n"  
            self.paramsFileString += "INPUT_DIR \n"  
            self.paramsFileString += "INPUT_CONVERT *\n"  
            self.paramsFileString += "INPUT\n"  
   
            self.firstFrame = True  
   
            self.imageFiles = []  
   
896         def doInitialization(self):         def doInitialization(self):
897            """            """
898            Initializes the time integration scheme            Initializes the time integration scheme
# Line 917  class SurfMovie(Model): Line 904  class SurfMovie(Model):
904            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
905            # wndow(south,west,north,east)            # wndow(south,west,north,east)
906    
907              # set up the movie parameters
908              self.paramsFileString = "REFERENCE_FRAME DECODED\n"
909              self.paramsFileString += "FRAME_RATE 24\n"
910              self.paramsFileString += "OUTPUT %s\n" % self.filename
911              self.paramsFileString += "PATTERN IBBPBBPBB\n"
912              self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
913              self.paramsFileString += "GOP_SIZE 20\n"
914              self.paramsFileString += "BSEARCH_ALG CROSS2\n"
915              self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
916              self.paramsFileString += "IQSCALE 10\n"
917              self.paramsFileString += "PQSCALE 11\n"
918              self.paramsFileString += "BQSCALE 16\n"
919              self.paramsFileString += "RANGE 8\n"
920              self.paramsFileString += "SLICES_PER_FRAME 1\n"
921              self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
922              self.paramsFileString += "INPUT_DIR \n"
923              self.paramsFileString += "INPUT_CONVERT *\n"
924              self.paramsFileString += "INPUT\n"
925    
926              self.firstFrame = True
927    
928              self.imageFiles = []
929    
930            # the bathymmetry colourmap            # the bathymmetry colourmap
931            data = []            data = []
932            data.append([-8000, 0,   0,   0])            data.append([-8000, 0,   0,   0])
# Line 1118  class SurfMovie(Model): Line 1128  class SurfMovie(Model):
1128          @param dt:          @param dt:
1129          """          """
1130          if self.t>=self.__next_t:          if self.t>=self.__next_t:
1131               print self.t,"write ",Lsup(self.wave_height)               print self.t,"write ",Lsup(self.wave_height)," to ",self.__fn
1132               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1133               # vtkobj=...               # vtkobj=...
1134               # save(self.__frame_name)               # save(self.__frame_name)
1135    
1136             # check that the file actually exists
1137             if not os.path.exists(self.__fn):
1138             raise IOError, "File not found: %s" % self.__fn
1139    
1140               # make a reader for the data               # make a reader for the data
1141               waveReader = vtk.vtkXMLUnstructuredGridReader()               waveReader = vtk.vtkXMLUnstructuredGridReader()
1142               waveReader.SetFileName(self.__fn)               waveReader.SetFileName(self.__fn)
# Line 1139  class SurfMovie(Model): Line 1153  class SurfMovie(Model):
1153               waveMapper = vtk.vtkDataSetMapper()               waveMapper = vtk.vtkDataSetMapper()
1154               waveMapper.SetInput(waveGrid)               waveMapper.SetInput(waveGrid)
1155               waveMapper.SetLookupTable(self.lutTrans)               waveMapper.SetLookupTable(self.lutTrans)
1156               waveMapper.SetScalarRange(-self.max_height, self.max_height)               waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1157    
1158               self.waveActor.SetMapper(waveMapper)               self.waveActor.SetMapper(waveMapper)
1159    
# Line 1207  class SurfMovie(Model): Line 1221  class SurfMovie(Model):
1221                print "An error occurred in mpeg conversion"                print "An error occurred in mpeg conversion"
1222    
1223            # now clean up the image files            # now clean up the image files
1224            print "Removing temporary image files"            if result == 0:
1225            os.unlink("%s.params" % self.filename)                print "Removing temporary image files"
1226            for fname in self.imageFiles:                os.unlink("%s.params" % self.filename)
1227                os.unlink(fname)                for fname in self.imageFiles:
1228                      os.unlink(fname)
1229    
1230  if __name__=="__main__":  def main():
1231     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
1232     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
1233    
1234     oc=OceanRegionCollector()     oc=OceanRegionCollector()
1235     oc.resolution=2     oc.north= 26.7
1236     oc.south=-45.5     oc.west= 102.9
1237     oc.north=-5.4     oc.range360= True
1238     oc.east=161.1     oc.east= 232.6
1239     oc.west=108.2     oc.resolution= 1.
1240     oc.range360=True     oc.south= -71.3
1241    
1242    
1243     b=Bathymetry()     b=Bathymetry()
# Line 1239  if __name__=="__main__": Line 1254  if __name__=="__main__":
1254    
1255     src=TsunamiSource()     src=TsunamiSource()
1256     src.domain=Link(oreg,"domain")     src.domain=Link(oreg,"domain")
1257     src.start_lat=-10.     src.decay_zone= 0.01
1258     src.end_lat=-12.     src.end_long= 185.
1259     src.start_long=110.     src.end_lat= -37.
1260     src.end_long=120.     src.width=0.5
1261     src.width=1.     src.start_long= 174.
1262     src.decay_zone=0.01     src.start_lat= -15.
1263     src.amplitude=1.     src.amplitude= 3
1264    
1265     ts=TsunamiInDeepWater()     ts=TsunamiInDeepWater()
1266     ts.domain=Link(oreg,"domain")     ts.domain=Link(oreg,"domain")
# Line 1254  if __name__=="__main__": Line 1269  if __name__=="__main__":
1269     ts.bathymetry=Link(oreg,"bathymetry")     ts.bathymetry=Link(oreg,"bathymetry")
1270    
1271     sq=Sequencer()     sq=Sequencer()
1272     sq.t_end=100000.     sq.t_end=15000.
1273    
1274     sm=SurfMovie()     sm=SurfMovie()
1275     sm.bathymetry=Link(b,"bathymetry")     sm.bathymetry=Link(b,"bathymetry")
1276     sm.wave_height=Link(ts,"wave_height")     sm.wave_height=Link(ts,"wave_height")
1277     sm.coastline=Link(oreg,"coastline")     sm.coastline=Link(oreg,"coastline")
1278     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
1279     sm.dt=5000.     sm.filename="mymovie.mpg"
1280     sm.filename="movie.mpg"     sm.north= 8.7
1281       sm.west= 138.9
1282       sm.dt= 50.
1283       sm.east= 196.6
1284       sm.south= -53.3
1285     sm.max_height=Link(src,"amplitude")     sm.max_height=Link(src,"amplitude")
1286        
1287     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1288     # s.writeXML()     s.writeXML()
1289     s.run()     s.run()
1290    
1291    if __name__=="__main__":
1292        from esys.modellib import tsunami
1293        tsunami.main()
1294    
1295  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

Legend:
Removed from v.280  
changed lines
  Added in v.918

  ViewVC Help
Powered by ViewVC 1.1.26