/[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 323 by gross, Tue Dec 6 06:18:00 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    
889             self.paramsFileString = "REFERENCE_FRAME DECODED\n"             self.paramsFileString = "REFERENCE_FRAME DECODED\n"
# Line 903  class SurfMovie(Model): Line 906  class SurfMovie(Model):
906    
907             self.firstFrame = True             self.firstFrame = True
908    
909         self.imageFiles = []             self.imageFiles = []
910    
911         def doInitialization(self):         def doInitialization(self):
912            """            """
# Line 917  class SurfMovie(Model): Line 920  class SurfMovie(Model):
920            # wndow(south,west,north,east)            # wndow(south,west,north,east)
921    
922            # the bathymmetry colourmap            # the bathymmetry colourmap
923        data = []            data = []
924        data.append([-8000, 0,   0,   0])            data.append([-8000, 0,   0,   0])
925        data.append([-7000, 0,   5,   25])            data.append([-7000, 0,   5,   25])
926        data.append([-6000, 0,   10,  50])            data.append([-6000, 0,   10,  50])
927        data.append([-5000, 0,   80,  125])            data.append([-5000, 0,   80,  125])
928        data.append([-4000, 0,   150, 200])            data.append([-4000, 0,   150, 200])
929        data.append([-3000, 86,  197, 184])            data.append([-3000, 86,  197, 184])
930        data.append([-2000, 172, 245, 168])            data.append([-2000, 172, 245, 168])
931        data.append([-1000, 211, 250, 211])            data.append([-1000, 211, 250, 211])
932        data.append([0,     16,  48,  16])            data.append([0,     16,  48,  16])
933        data.append([300,   70,  150, 50])            data.append([300,   70,  150, 50])
934        data.append([500,   146, 126, 60])            data.append([500,   146, 126, 60])
935        data.append([1000,  198, 178, 80])            data.append([1000,  198, 178, 80])
936        data.append([1250,  250, 230, 100])            data.append([1250,  250, 230, 100])
937        data.append([1500,  250, 234, 126])            data.append([1500,  250, 234, 126])
938        data.append([1750,  252, 238, 152])            data.append([1750,  252, 238, 152])
939        data.append([2000,  252, 243, 177])            data.append([2000,  252, 243, 177])
940        data.append([2250,  253, 249, 216])            data.append([2250,  253, 249, 216])
941        data.append([2500,  255, 255, 255])            data.append([2500,  255, 255, 255])
942    
943        # the amount to scale the data by            # the amount to scale the data by
944            scale = 255.0            scale = 255.0
945            numColours = len(data)            numColours = len(data)
946    
947        # convert the colourmap into something vtk is more happy with            # convert the colourmap into something vtk is more happy with
948            height = numarray.zeros(numColours, numarray.Float)            height = numarray.zeros(numColours, numarray.Float)
949            red = numarray.zeros(numColours, numarray.Float)            red = numarray.zeros(numColours, numarray.Float)
950            green = numarray.zeros(numColours, numarray.Float)            green = numarray.zeros(numColours, numarray.Float)
# Line 1002  class SurfMovie(Model): Line 1005  class SurfMovie(Model):
1005                  h += 1                  h += 1
1006                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])                  transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1007    
1008        # set up the lookup table for the wave data            # set up the lookup table for the wave data
1009        refLut = vtk.vtkLookupTable()            refLut = vtk.vtkLookupTable()
1010        self.lutTrans = vtk.vtkLookupTable()            self.lutTrans = vtk.vtkLookupTable()
1011        refLut.Build()            refLut.Build()
1012        alpha = 0.4   # alpha channel value            alpha = 0.7   # alpha channel value
1013        for i in range(256):            for i in range(256):
1014            (r,g,b,a) = refLut.GetTableValue(255-i)                (r,g,b,a) = refLut.GetTableValue(255-i)
1015            if g == 1.0 and b < 0.5 and r < 0.5:                if g == 1.0 and b < 0.5 and r < 0.5:
1016            self.lutTrans.SetTableValue(i, r, g, b, 0.0)                    self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1017            else:                else:
1018            self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)                    self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1019    
1020            print "Generating the bathymmetry vtk object..."            print "Generating the bathymmetry vtk object..."
1021    
# Line 1064  class SurfMovie(Model): Line 1067  class SurfMovie(Model):
1067    
1068            self.ren.AddActor(bathActor)            self.ren.AddActor(bathActor)
1069    
1070        ### now add the coastline            ### now add the coastline
1071            print "Loading the coastline data..."            print "Loading the coastline data..."
1072    
1073        # make the coastline points            # make the coastline points
1074        coastPoints = vtk.vtkPoints()            coastPoints = vtk.vtkPoints()
1075    
1076        # make the coastline grid            # make the coastline grid
1077        coastGrid = vtk.vtkUnstructuredGrid()            coastGrid = vtk.vtkUnstructuredGrid()
1078    
1079        # now point the points and lines into the grid            # now point the points and lines into the grid
1080        totalCoastPoints = 0            totalCoastPoints = 0
1081        for polyline in self.coastline.polylines:            for polyline in self.coastline.polylines:
1082            numPoints = len(polyline)                numPoints = len(polyline)
1083            coastLine = vtk.vtkPolyLine()                coastLine = vtk.vtkPolyLine()
1084            coastLine.GetPointIds().SetNumberOfIds(numPoints)                coastLine.GetPointIds().SetNumberOfIds(numPoints)
1085            j = 0                j = 0
1086            for point in polyline:                for point in polyline:
1087            coastLine.GetPointIds().SetId(j, j+totalCoastPoints)                    coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1088            coastPoints.InsertNextPoint(point.long, point.lat, 0.0)                    coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1089            j += 1                    j += 1
1090            coastGrid.InsertNextCell(coastLine.GetCellType(),                coastGrid.InsertNextCell(coastLine.GetCellType(),
1091                coastLine.GetPointIds())                        coastLine.GetPointIds())
1092            totalCoastPoints += numPoints                totalCoastPoints += numPoints
1093    
1094        coastGrid.SetPoints(coastPoints)            coastGrid.SetPoints(coastPoints)
1095    
1096        # make the coast's mapper            # make the coast's mapper
1097        coastMapper = vtk.vtkDataSetMapper()            coastMapper = vtk.vtkDataSetMapper()
1098        coastMapper.SetInput(coastGrid)            coastMapper.SetInput(coastGrid)
1099    
1100        # make its actor            # make its actor
1101        coastActor = vtk.vtkActor()            coastActor = vtk.vtkActor()
1102        coastActor.SetMapper(coastMapper)            coastActor.SetMapper(coastMapper)
1103        coastActor.GetProperty().SetColor(0,0,0)            coastActor.GetProperty().SetColor(0,0,0)
1104    
1105        # add the actor to the renderer            # add the actor to the renderer
1106        self.ren.AddActor(coastActor)            self.ren.AddActor(coastActor)
1107    
1108            # set up the actor for the wave            # set up the actor for the wave
1109            self.waveActor = vtk.vtkActor()            self.waveActor = vtk.vtkActor()
# Line 1117  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)
1127    
1128           # make a reader for the data               # make a reader for the data
1129           waveReader = vtk.vtkXMLUnstructuredGridReader()               waveReader = vtk.vtkXMLUnstructuredGridReader()
1130           waveReader.SetFileName(self.__fn)               waveReader.SetFileName(self.__fn)
1131           waveReader.Update()               waveReader.Update()
1132    
1133               # make the grid               # make the grid
1134               waveGrid = waveReader.GetOutput()               waveGrid = waveReader.GetOutput()
1135           waveGrid.Update()               waveGrid.Update()
1136    
1137               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()               (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1138               print "Wave height range %f - %f" % (zMin, zMax)               print "Wave height range %f - %f" % (zMin, zMax)
# Line 1137  class SurfMovie(Model): Line 1140  class SurfMovie(Model):
1140               # make a mapper for the grid               # make a mapper for the grid
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(zMin, zMax)               waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1145    
1146               self.waveActor.SetMapper(waveMapper)               self.waveActor.SetMapper(waveMapper)
1147    
# Line 1165  class SurfMovie(Model): Line 1168  class SurfMovie(Model):
1168               outWriter.Write()               outWriter.Write()
1169               print "Wrote %s" % imgFname               print "Wrote %s" % imgFname
1170    
1171           # helpful for debugging:               # helpful for debugging:
1172           #os.system("display %s" % imgFname)               #os.system("display %s" % imgFname)
1173    
1174               self.paramsFileString += "%s\n" % imgFname               self.paramsFileString += "%s\n" % imgFname
1175           self.imageFiles.append(imgFname)               self.imageFiles.append(imgFname)
1176    
1177               self.__next_t+=self.dt               self.__next_t+=self.dt
1178    
# Line 1205  class SurfMovie(Model): Line 1208  class SurfMovie(Model):
1208            if result != 0:            if result != 0:
1209                print "An error occurred in mpeg conversion"                print "An error occurred in mpeg conversion"
1210    
1211        # now clean up the image files            # now clean up the image files
1212        print "Removing temporary image files"            print "Removing temporary image files"
1213        os.unlink("%s.params" % self.filename)            os.unlink("%s.params" % self.filename)
1214        for fname in self.imageFiles:            for fname in self.imageFiles:
1215            os.unlink(fname)                os.unlink(fname)
1216    
1217  if __name__=="__main__":  def main():
1218     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
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 1238  if __name__=="__main__": Line 1241  if __name__=="__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 1253  if __name__=="__main__": Line 1256  if __name__=="__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")
1273        
1274     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1275     # s.writeXML()     s.writeXML()
1276     s.run()     s.run()
1277    
1278    if __name__=="__main__":
1279        from esys.modellib import tsunami
1280        tsunami.main()
1281    
1282  # vim: expandtab shiftwidth=4:  # vim: expandtab shiftwidth=4:

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

  ViewVC Help
Powered by ViewVC 1.1.26