/[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 276 by cochrane, Wed Nov 30 04:31:40 2005 UTC revision 628 by elspeth, Thu Mar 23 02:27:57 2006 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 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 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 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 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 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 881  class SurfMovie(Model): Line 890  class SurfMovie(Model):
890                                   north=5.,                                   north=5.,
891                                   east=3.,                                   east=3.,
892                                   west=15.,                                   west=15.,
893                                     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 916  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])
933        data.append([-7000, 0,   5,   25])            data.append([-7000, 0,   5,   25])
934        data.append([-6000, 0,   10,  50])            data.append([-6000, 0,   10,  50])
935        data.append([-5000, 0,   80,  125])            data.append([-5000, 0,   80,  125])
936        data.append([-4000, 0,   150, 200])            data.append([-4000, 0,   150, 200])
937        data.append([-3000, 86,  197, 184])            data.append([-3000, 86,  197, 184])
938        data.append([-2000, 172, 245, 168])            data.append([-2000, 172, 245, 168])
939        data.append([-1000, 211, 250, 211])            data.append([-1000, 211, 250, 211])
940        data.append([0,     16,  48,  16])            data.append([0,     16,  48,  16])
941        data.append([300,   70,  150, 50])            data.append([300,   70,  150, 50])
942        data.append([500,   146, 126, 60])            data.append([500,   146, 126, 60])
943        data.append([1000,  198, 178, 80])            data.append([1000,  198, 178, 80])
944        data.append([1250,  250, 230, 100])            data.append([1250,  250, 230, 100])
945        data.append([1500,  250, 234, 126])            data.append([1500,  250, 234, 126])
946        data.append([1750,  252, 238, 152])            data.append([1750,  252, 238, 152])
947        data.append([2000,  252, 243, 177])            data.append([2000,  252, 243, 177])
948        data.append([2250,  253, 249, 216])            data.append([2250,  253, 249, 216])
949        data.append([2500,  255, 255, 255])            data.append([2500,  255, 255, 255])
950    
951        # the amount to scale the data by            # the amount to scale the data by
952            scale = 255.0            scale = 255.0
953            numColours = len(data)            numColours = len(data)
954    
955        # convert the colourmap into something vtk is more happy with            # convert the colourmap into something vtk is more happy with
956            height = numarray.zeros(numColours, numarray.Float)            height = numarray.zeros(numColours, numarray.Float)
957            red = numarray.zeros(numColours, numarray.Float)            red = numarray.zeros(numColours, numarray.Float)
958            green = numarray.zeros(numColours, numarray.Float)            green = numarray.zeros(numColours, numarray.Float)
# Line 1002  class SurfMovie(Model): Line 1013  class SurfMovie(Model):
1013                  h += 1                  h += 1
1014                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1015    
1016        # set up the lookup table for the wave data            # set up the lookup table for the wave data
1017        refLut = vtk.vtkLookupTable()            refLut = vtk.vtkLookupTable()
1018        self.lutTrans = vtk.vtkLookupTable()            self.lutTrans = vtk.vtkLookupTable()
1019        refLut.Build()            refLut.Build()
1020        alpha = 0.4   # alpha channel value            alpha = 0.7   # alpha channel value
1021        for i in range(256):            for i in range(256):
1022            (r,g,b,a) = refLut.GetTableValue(255-i)                (r,g,b,a) = refLut.GetTableValue(255-i)
1023            if g == 1.0 and b < 0.5 and r < 0.5:                if g == 1.0 and b < 0.5 and r < 0.5:
1024            self.lutTrans.SetTableValue(i, r, g, b, 0.0)                    self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1025            else:                else:
1026            self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)                    self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1027    
1028            print "Generating the bathymmetry vtk object..."            print "Generating the bathymmetry vtk object..."
1029    
# Line 1064  class SurfMovie(Model): Line 1075  class SurfMovie(Model):
1075    
1076            self.ren.AddActor(bathActor)            self.ren.AddActor(bathActor)
1077    
1078        ### now add the coastline            ### now add the coastline
1079            print "Loading the coastline data..."            print "Loading the coastline data..."
1080    
1081        # make the coastline points            # make the coastline points
1082        coastPoints = vtk.vtkPoints()            coastPoints = vtk.vtkPoints()
1083    
1084        # make the coastline grid            # make the coastline grid
1085        coastGrid = vtk.vtkUnstructuredGrid()            coastGrid = vtk.vtkUnstructuredGrid()
1086    
1087        # now point the points and lines into the grid            # now point the points and lines into the grid
1088        totalCoastPoints = 0            totalCoastPoints = 0
1089        for polyline in self.coastline.polylines:            for polyline in self.coastline.polylines:
1090            numPoints = len(polyline)                numPoints = len(polyline)
1091            coastLine = vtk.vtkPolyLine()                coastLine = vtk.vtkPolyLine()
1092            coastLine.GetPointIds().SetNumberOfIds(numPoints)                coastLine.GetPointIds().SetNumberOfIds(numPoints)
1093            j = 0                j = 0
1094            for point in polyline:                for point in polyline:
1095            coastLine.GetPointIds().SetId(j, j+totalCoastPoints)                    coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1096            coastPoints.InsertNextPoint(point.long, point.lat, 0.0)                    coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1097            j += 1                    j += 1
1098            coastGrid.InsertNextCell(coastLine.GetCellType(),                coastGrid.InsertNextCell(coastLine.GetCellType(),
1099                coastLine.GetPointIds())                        coastLine.GetPointIds())
1100            totalCoastPoints += numPoints                totalCoastPoints += numPoints
1101    
1102        coastGrid.SetPoints(coastPoints)            coastGrid.SetPoints(coastPoints)
1103    
1104        # make the coast's mapper            # make the coast's mapper
1105        coastMapper = vtk.vtkDataSetMapper()            coastMapper = vtk.vtkDataSetMapper()
1106        coastMapper.SetInput(coastGrid)            coastMapper.SetInput(coastGrid)
1107    
1108        # make its actor            # make its actor
1109        coastActor = vtk.vtkActor()            coastActor = vtk.vtkActor()
1110        coastActor.SetMapper(coastMapper)            coastActor.SetMapper(coastMapper)
1111        coastActor.GetProperty().SetColor(0,0,0)            coastActor.GetProperty().SetColor(0,0,0)
1112    
1113        # add the actor to the renderer            # add the actor to the renderer
1114        self.ren.AddActor(coastActor)            self.ren.AddActor(coastActor)
1115    
1116            # set up the actor for the wave            # set up the actor for the wave
1117            self.waveActor = vtk.vtkActor()            self.waveActor = vtk.vtkActor()
# Line 1117  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           # make a reader for the data           # check that the file actually exists
1137           waveReader = vtk.vtkXMLUnstructuredGridReader()           if not os.path.exists(self.__fn):
1138           waveReader.SetFileName(self.__fn)           raise IOError, "File not found: %s" % self.__fn
1139           waveReader.Update()  
1140                 # make a reader for the data
1141                 waveReader = vtk.vtkXMLUnstructuredGridReader()
1142                 waveReader.SetFileName(self.__fn)
1143                 waveReader.Update()
1144    
1145               # make the grid               # make the grid
1146               waveGrid = waveReader.GetOutput()               waveGrid = waveReader.GetOutput()
1147           waveGrid.Update()               waveGrid.Update()
1148    
1149               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1150               print "Wave height range %f - %f" % (zMin, zMax)               print "Wave height range %f - %f" % (zMin, zMax)
# Line 1137  class SurfMovie(Model): Line 1152  class SurfMovie(Model):
1152               # make a mapper for the grid               # make a mapper for the grid
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(zMin, zMax)               waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1157    
1158               self.waveActor.SetMapper(waveMapper)               self.waveActor.SetMapper(waveMapper)
1159    
# Line 1165  class SurfMovie(Model): Line 1180  class SurfMovie(Model):
1180               outWriter.Write()               outWriter.Write()
1181               print "Wrote %s" % imgFname               print "Wrote %s" % imgFname
1182    
1183           # helpful for debugging:               # helpful for debugging:
1184           #os.system("display %s" % imgFname)               #os.system("display %s" % imgFname)
1185    
1186               self.paramsFileString += "%s\n" % imgFname               self.paramsFileString += "%s\n" % imgFname
1187           self.imageFiles.append(imgFname)               self.imageFiles.append(imgFname)
1188    
1189               self.__next_t+=self.dt               self.__next_t+=self.dt
1190    
# Line 1205  class SurfMovie(Model): Line 1220  class SurfMovie(Model):
1220            if result != 0:            if result != 0:
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 1238  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 1253  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")
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.276  
changed lines
  Added in v.628

  ViewVC Help
Powered by ViewVC 1.1.26