/[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 408 by cochrane, Fri Dec 23 00:56:56 2005 UTC
# Line 9  from esys.escript.linearPDEs import Line Line 9  from esys.escript.linearPDEs import Line
9  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
10  import numarray  import numarray
11  import math  import math
12  import urllib  #import urllib
13    import urllib2
14    
15  EPS=1.e-5  EPS=1.e-8
16    
17  #====================================================================================================================================================  #====================================================================================================================================================
18  class GridData:  class GridData:
# Line 353  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 421  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 609  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 641  class OceanRegionCollector(Model): Line 642  class OceanRegionCollector(Model):
642             """             """
643             c=self.__mergeParameters(self.coastline_source)             c=self.__mergeParameters(self.coastline_source)
644             b=self.__mergeParameters(self.bathymetry_source)             b=self.__mergeParameters(self.bathymetry_source)
645             self.coastline_stream=urllib.urlopen(c)             self.coastline_stream=urllib2.urlopen(c)
646             self.bathymetry_stream=urllib.urlopen(b)             self.bathymetry_stream=urllib2.urlopen(b)
647    
648        def __mergeParameters(self,txt):        def __mergeParameters(self,txt):
649             return txt.replace("%%west%%",str(self.west))\             return txt.replace("%%west%%",str(self.west))\
# Line 691  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 788  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 826  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 881  class SurfMovie(Model): Line 883  class SurfMovie(Model):
883                                   north=5.,                                   north=5.,
884                                   east=3.,                                   east=3.,
885                                   west=15.,                                   west=15.,
886                                     max_height=1.0,
887                                   filename="movie.mpg")                                   filename="movie.mpg")
888    
            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 = []  
   
889         def doInitialization(self):         def doInitialization(self):
890            """            """
891            Initializes the time integration scheme            Initializes the time integration scheme
# Line 916  class SurfMovie(Model): Line 897  class SurfMovie(Model):
897            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
898            # wndow(south,west,north,east)            # wndow(south,west,north,east)
899    
900              # set up the movie parameters
901              self.paramsFileString = "REFERENCE_FRAME DECODED\n"
902              self.paramsFileString += "FRAME_RATE 24\n"
903              self.paramsFileString += "OUTPUT %s\n" % self.filename
904              self.paramsFileString += "PATTERN IBBPBBPBB\n"
905              self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
906              self.paramsFileString += "GOP_SIZE 20\n"
907              self.paramsFileString += "BSEARCH_ALG CROSS2\n"
908              self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
909              self.paramsFileString += "IQSCALE 10\n"
910              self.paramsFileString += "PQSCALE 11\n"
911              self.paramsFileString += "BQSCALE 16\n"
912              self.paramsFileString += "RANGE 8\n"
913              self.paramsFileString += "SLICES_PER_FRAME 1\n"
914              self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
915              self.paramsFileString += "INPUT_DIR \n"
916              self.paramsFileString += "INPUT_CONVERT *\n"
917              self.paramsFileString += "INPUT\n"
918    
919              self.firstFrame = True
920    
921              self.imageFiles = []
922    
923            # the bathymmetry colourmap            # the bathymmetry colourmap
924        data = []            data = []
925        data.append([-8000, 0,   0,   0])            data.append([-8000, 0,   0,   0])
926        data.append([-7000, 0,   5,   25])            data.append([-7000, 0,   5,   25])
927        data.append([-6000, 0,   10,  50])            data.append([-6000, 0,   10,  50])
928        data.append([-5000, 0,   80,  125])            data.append([-5000, 0,   80,  125])
929        data.append([-4000, 0,   150, 200])            data.append([-4000, 0,   150, 200])
930        data.append([-3000, 86,  197, 184])            data.append([-3000, 86,  197, 184])
931        data.append([-2000, 172, 245, 168])            data.append([-2000, 172, 245, 168])
932        data.append([-1000, 211, 250, 211])            data.append([-1000, 211, 250, 211])
933        data.append([0,     16,  48,  16])            data.append([0,     16,  48,  16])
934        data.append([300,   70,  150, 50])            data.append([300,   70,  150, 50])
935        data.append([500,   146, 126, 60])            data.append([500,   146, 126, 60])
936        data.append([1000,  198, 178, 80])            data.append([1000,  198, 178, 80])
937        data.append([1250,  250, 230, 100])            data.append([1250,  250, 230, 100])
938        data.append([1500,  250, 234, 126])            data.append([1500,  250, 234, 126])
939        data.append([1750,  252, 238, 152])            data.append([1750,  252, 238, 152])
940        data.append([2000,  252, 243, 177])            data.append([2000,  252, 243, 177])
941        data.append([2250,  253, 249, 216])            data.append([2250,  253, 249, 216])
942        data.append([2500,  255, 255, 255])            data.append([2500,  255, 255, 255])
943    
944        # the amount to scale the data by            # the amount to scale the data by
945            scale = 255.0            scale = 255.0
946            numColours = len(data)            numColours = len(data)
947    
948        # convert the colourmap into something vtk is more happy with            # convert the colourmap into something vtk is more happy with
949            height = numarray.zeros(numColours, numarray.Float)            height = numarray.zeros(numColours, numarray.Float)
950            red = numarray.zeros(numColours, numarray.Float)            red = numarray.zeros(numColours, numarray.Float)
951            green = numarray.zeros(numColours, numarray.Float)            green = numarray.zeros(numColours, numarray.Float)
# Line 1002  class SurfMovie(Model): Line 1006  class SurfMovie(Model):
1006                  h += 1                  h += 1
1007                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1008    
1009        # set up the lookup table for the wave data            # set up the lookup table for the wave data
1010        refLut = vtk.vtkLookupTable()            refLut = vtk.vtkLookupTable()
1011        self.lutTrans = vtk.vtkLookupTable()            self.lutTrans = vtk.vtkLookupTable()
1012        refLut.Build()            refLut.Build()
1013        alpha = 0.4   # alpha channel value            alpha = 0.7   # alpha channel value
1014        for i in range(256):            for i in range(256):
1015            (r,g,b,a) = refLut.GetTableValue(255-i)                (r,g,b,a) = refLut.GetTableValue(255-i)
1016            if g == 1.0 and b < 0.5 and r < 0.5:                if g == 1.0 and b < 0.5 and r < 0.5:
1017            self.lutTrans.SetTableValue(i, r, g, b, 0.0)                    self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1018            else:                else:
1019            self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)                    self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1020    
1021            print "Generating the bathymmetry vtk object..."            print "Generating the bathymmetry vtk object..."
1022    
# Line 1064  class SurfMovie(Model): Line 1068  class SurfMovie(Model):
1068    
1069            self.ren.AddActor(bathActor)            self.ren.AddActor(bathActor)
1070    
1071        ### now add the coastline            ### now add the coastline
1072            print "Loading the coastline data..."            print "Loading the coastline data..."
1073    
1074        # make the coastline points            # make the coastline points
1075        coastPoints = vtk.vtkPoints()            coastPoints = vtk.vtkPoints()
1076    
1077        # make the coastline grid            # make the coastline grid
1078        coastGrid = vtk.vtkUnstructuredGrid()            coastGrid = vtk.vtkUnstructuredGrid()
1079    
1080        # now point the points and lines into the grid            # now point the points and lines into the grid
1081        totalCoastPoints = 0            totalCoastPoints = 0
1082        for polyline in self.coastline.polylines:            for polyline in self.coastline.polylines:
1083            numPoints = len(polyline)                numPoints = len(polyline)
1084            coastLine = vtk.vtkPolyLine()                coastLine = vtk.vtkPolyLine()
1085            coastLine.GetPointIds().SetNumberOfIds(numPoints)                coastLine.GetPointIds().SetNumberOfIds(numPoints)
1086            j = 0                j = 0
1087            for point in polyline:                for point in polyline:
1088            coastLine.GetPointIds().SetId(j, j+totalCoastPoints)                    coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1089            coastPoints.InsertNextPoint(point.long, point.lat, 0.0)                    coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1090            j += 1                    j += 1
1091            coastGrid.InsertNextCell(coastLine.GetCellType(),                coastGrid.InsertNextCell(coastLine.GetCellType(),
1092                coastLine.GetPointIds())                        coastLine.GetPointIds())
1093            totalCoastPoints += numPoints                totalCoastPoints += numPoints
1094    
1095        coastGrid.SetPoints(coastPoints)            coastGrid.SetPoints(coastPoints)
1096    
1097        # make the coast's mapper            # make the coast's mapper
1098        coastMapper = vtk.vtkDataSetMapper()            coastMapper = vtk.vtkDataSetMapper()
1099        coastMapper.SetInput(coastGrid)            coastMapper.SetInput(coastGrid)
1100    
1101        # make its actor            # make its actor
1102        coastActor = vtk.vtkActor()            coastActor = vtk.vtkActor()
1103        coastActor.SetMapper(coastMapper)            coastActor.SetMapper(coastMapper)
1104        coastActor.GetProperty().SetColor(0,0,0)            coastActor.GetProperty().SetColor(0,0,0)
1105    
1106        # add the actor to the renderer            # add the actor to the renderer
1107        self.ren.AddActor(coastActor)            self.ren.AddActor(coastActor)
1108    
1109            # set up the actor for the wave            # set up the actor for the wave
1110            self.waveActor = vtk.vtkActor()            self.waveActor = vtk.vtkActor()
# Line 1117  class SurfMovie(Model): Line 1121  class SurfMovie(Model):
1121          @param dt:          @param dt:
1122          """          """
1123          if self.t>=self.__next_t:          if self.t>=self.__next_t:
1124               print self.t,"write ",Lsup(self.wave_height)               print self.t,"write ",Lsup(self.wave_height)," to ",self.__fn
1125               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1126               # vtkobj=...               # vtkobj=...
1127               # save(self.__frame_name)               # save(self.__frame_name)
1128    
1129           # make a reader for the data           # check that the file actually exists
1130           waveReader = vtk.vtkXMLUnstructuredGridReader()           if not os.path.exists(self.__fn):
1131           waveReader.SetFileName(self.__fn)           raise IOError, "File not found: %s" % self.__fn
1132           waveReader.Update()  
1133                 # make a reader for the data
1134                 waveReader = vtk.vtkXMLUnstructuredGridReader()
1135                 waveReader.SetFileName(self.__fn)
1136                 waveReader.Update()
1137    
1138               # make the grid               # make the grid
1139               waveGrid = waveReader.GetOutput()               waveGrid = waveReader.GetOutput()
1140           waveGrid.Update()               waveGrid.Update()
1141    
1142               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1143               print "Wave height range %f - %f" % (zMin, zMax)               print "Wave height range %f - %f" % (zMin, zMax)
# Line 1137  class SurfMovie(Model): Line 1145  class SurfMovie(Model):
1145               # make a mapper for the grid               # make a mapper for the grid
1146               waveMapper = vtk.vtkDataSetMapper()               waveMapper = vtk.vtkDataSetMapper()
1147               waveMapper.SetInput(waveGrid)               waveMapper.SetInput(waveGrid)
1148           waveMapper.SetLookupTable(self.lutTrans)               waveMapper.SetLookupTable(self.lutTrans)
1149               waveMapper.SetScalarRange(zMin, zMax)               waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1150    
1151               self.waveActor.SetMapper(waveMapper)               self.waveActor.SetMapper(waveMapper)
1152    
# Line 1165  class SurfMovie(Model): Line 1173  class SurfMovie(Model):
1173               outWriter.Write()               outWriter.Write()
1174               print "Wrote %s" % imgFname               print "Wrote %s" % imgFname
1175    
1176           # helpful for debugging:               # helpful for debugging:
1177           #os.system("display %s" % imgFname)               #os.system("display %s" % imgFname)
1178    
1179               self.paramsFileString += "%s\n" % imgFname               self.paramsFileString += "%s\n" % imgFname
1180           self.imageFiles.append(imgFname)               self.imageFiles.append(imgFname)
1181    
1182               self.__next_t+=self.dt               self.__next_t+=self.dt
1183    
# Line 1205  class SurfMovie(Model): Line 1213  class SurfMovie(Model):
1213            if result != 0:            if result != 0:
1214                print "An error occurred in mpeg conversion"                print "An error occurred in mpeg conversion"
1215    
1216        # now clean up the image files            # now clean up the image files
1217        print "Removing temporary image files"            if result == 0:
1218        os.unlink("%s.params" % self.filename)                print "Removing temporary image files"
1219        for fname in self.imageFiles:                os.unlink("%s.params" % self.filename)
1220            os.unlink(fname)                for fname in self.imageFiles:
1221                      os.unlink(fname)
1222    
1223  if __name__=="__main__":  def main():
1224     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
1225     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
1226    
1227     oc=OceanRegionCollector()     oc=OceanRegionCollector()
1228     oc.resolution=2     oc.north= 26.7
1229     oc.south=-45.5     oc.west= 102.9
1230     oc.north=-5.4     oc.range360= True
1231     oc.east=161.1     oc.east= 232.6
1232     oc.west=108.2     oc.resolution= 1.
1233     oc.range360=True     oc.south= -71.3
1234    
1235    
1236     b=Bathymetry()     b=Bathymetry()
# Line 1238  if __name__=="__main__": Line 1247  if __name__=="__main__":
1247    
1248     src=TsunamiSource()     src=TsunamiSource()
1249     src.domain=Link(oreg,"domain")     src.domain=Link(oreg,"domain")
1250     src.start_lat=-10.     src.decay_zone= 0.01
1251     src.end_lat=-12.     src.end_long= 185.
1252     src.start_long=110.     src.end_lat= -37.
1253     src.end_long=120.     src.width=0.5
1254     src.width=1.     src.start_long= 174.
1255     src.decay_zone=0.01     src.start_lat= -15.
1256     src.amplitude=1.     src.amplitude= 3
1257    
1258     ts=TsunamiInDeepWater()     ts=TsunamiInDeepWater()
1259     ts.domain=Link(oreg,"domain")     ts.domain=Link(oreg,"domain")
# Line 1253  if __name__=="__main__": Line 1262  if __name__=="__main__":
1262     ts.bathymetry=Link(oreg,"bathymetry")     ts.bathymetry=Link(oreg,"bathymetry")
1263    
1264     sq=Sequencer()     sq=Sequencer()
1265     sq.t_end=100000.     sq.t_end=15000.
1266    
1267     sm=SurfMovie()     sm=SurfMovie()
1268     sm.bathymetry=Link(b,"bathymetry")     sm.bathymetry=Link(b,"bathymetry")
1269     sm.wave_height=Link(ts,"wave_height")     sm.wave_height=Link(ts,"wave_height")
1270     sm.coastline=Link(oreg,"coastline")     sm.coastline=Link(oreg,"coastline")
1271     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
1272     sm.dt=5000.     sm.filename="mymovie.mpg"
1273     sm.filename="movie.mpg"     sm.north= 8.7
1274       sm.west= 138.9
1275       sm.dt= 50.
1276       sm.east= 196.6
1277       sm.south= -53.3
1278       sm.max_height=Link(src,"amplitude")
1279        
1280     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1281     # s.writeXML()     s.writeXML()
1282     s.run()     s.run()
1283    
1284    if __name__=="__main__":
1285        from esys.modellib import tsunami
1286        tsunami.main()
1287    
1288  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

Legend:
Removed from v.276  
changed lines
  Added in v.408

  ViewVC Help
Powered by ViewVC 1.1.26