/[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 251 by gross, Tue Nov 29 05:54:06 2005 UTC revision 281 by elspeth, Wed Nov 30 09:10:59 2005 UTC
# Line 1  Line 1 
1  # $Id$  #!/usr/bin/env python
2    
3    # $Id$
4    
5  import os  import os, sys
6    import vtk
7  from esys.escript import *  from esys.escript import *
8  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
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-5
16    
# Line 639  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 815  class TsunamiInDeepWater(Model): Line 818  class TsunamiInDeepWater(Model):
818                                   wave_velocity=0.,                                   wave_velocity=0.,
819                                   initial_time_step=None,                                   initial_time_step=None,
820                                   bathymetry=1.,                                   bathymetry=1.,
821                                   safety_factor=0.1)                                   safety_factor=1.)
822    
823         def doInitialization(self):         def doInitialization(self):
824             """             """
# Line 828  class TsunamiInDeepWater(Model): Line 831  class TsunamiInDeepWater(Model):
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             if self.initial_time_step==None: self.initial_time_step=self.__dt             if self.initial_time_step==None: self.initial_time_step=self.__dt
834             self.trace("maximum wave velocity %s m/sec^2"%c_max)             self.trace("maximum wave velocity %s m/sec"%c_max)
835             self.trace("Time step size is %s sec"%self.__dt)             self.trace("Time step size is %s sec"%self.__dt)
836    
837    
# Line 851  class TsunamiInDeepWater(Model): Line 854  class TsunamiInDeepWater(Model):
854             @type dt: C{float}             @type dt: C{float}
855             """             """
856             self.__pde.setValue(X=-self.__c2*grad(self.wave_height))             self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
857    
858             new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()             new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
859    
860             self.wave_velocity=(new_height-self.wave_height)/dt             self.wave_velocity=(new_height-self.wave_height)/dt
861             self.wave_height=new_height             self.wave_height=new_height
862             self.initial_time_step=dt             self.initial_time_step=dt
# Line 877  class SurfMovie(Model): Line 882  class SurfMovie(Model):
882                                   north=5.,                                   north=5.,
883                                   east=3.,                                   east=3.,
884                                   west=15.,                                   west=15.,
885                                     max_height=1.0,
886                                   filename="movie.mpg")                                   filename="movie.mpg")
887    
888               self.paramsFileString = "REFERENCE_FRAME DECODED\n"
889               self.paramsFileString += "FRAME_RATE 24\n"
890               self.paramsFileString += "OUTPUT %s\n" % self.filename
891               self.paramsFileString += "PATTERN IBBPBBPBB\n"
892               self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
893               self.paramsFileString += "GOP_SIZE 20\n"
894               self.paramsFileString += "BSEARCH_ALG CROSS2\n"
895               self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
896               self.paramsFileString += "IQSCALE 10\n"
897               self.paramsFileString += "PQSCALE 11\n"
898               self.paramsFileString += "BQSCALE 16\n"
899               self.paramsFileString += "RANGE 8\n"
900               self.paramsFileString += "SLICES_PER_FRAME 1\n"
901               self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
902               self.paramsFileString += "INPUT_DIR \n"
903               self.paramsFileString += "INPUT_CONVERT *\n"
904               self.paramsFileString += "INPUT\n"
905    
906               self.firstFrame = True
907    
908               self.imageFiles = []
909    
910         def doInitialization(self):         def doInitialization(self):
911            """            """
912            Initializes the time integartion scheme            Initializes the time integration scheme
913            """            """
914            self.__fn=os.tempnam()+".xml"            self.__fn=os.tempnam()+".xml"
915            self.__frame_name=os.tempnam()            self.__frame_name=os.tempnam()
# Line 890  class SurfMovie(Model): Line 918  class SurfMovie(Model):
918            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
919            # wndow(south,west,north,east)            # wndow(south,west,north,east)
920    
921              # the bathymmetry colourmap
922              data = []
923              data.append([-8000, 0,   0,   0])
924              data.append([-7000, 0,   5,   25])
925              data.append([-6000, 0,   10,  50])
926              data.append([-5000, 0,   80,  125])
927              data.append([-4000, 0,   150, 200])
928              data.append([-3000, 86,  197, 184])
929              data.append([-2000, 172, 245, 168])
930              data.append([-1000, 211, 250, 211])
931              data.append([0,     16,  48,  16])
932              data.append([300,   70,  150, 50])
933              data.append([500,   146, 126, 60])
934              data.append([1000,  198, 178, 80])
935              data.append([1250,  250, 230, 100])
936              data.append([1500,  250, 234, 126])
937              data.append([1750,  252, 238, 152])
938              data.append([2000,  252, 243, 177])
939              data.append([2250,  253, 249, 216])
940              data.append([2500,  255, 255, 255])
941    
942              # the amount to scale the data by
943              scale = 255.0
944              numColours = len(data)
945    
946              # convert the colourmap into something vtk is more happy with
947              height = numarray.zeros(numColours, numarray.Float)
948              red = numarray.zeros(numColours, numarray.Float)
949              green = numarray.zeros(numColours, numarray.Float)
950              blue = numarray.zeros(numColours, numarray.Float)
951              for i in range(numColours):
952                  (h, r, g, b) = data[i]
953                  height[i] = float(h)
954                  red[i] = float(r)/scale
955                  green[i] = float(g)/scale
956                  blue[i] = float(b)/scale
957    
958              print "Loading bathymmetry data..."
959              # grab the z data
960              bathZData = self.bathymetry.getData()
961    
962              # get the origin
963              bathOrigin = self.bathymetry.getOrigin()
964              # get the delta_x and delta_y
965              bathSpacing = self.bathymetry.getSpacing()
966    
967              # now construct the x and y data from the spacing and the origin
968              numXPoints = bathZData.shape[1]
969              numYPoints = bathZData.shape[0]
970              numPoints = numXPoints*numYPoints
971    
972              bathXData = numarray.zeros(numXPoints, numarray.Float)
973              bathYData = numarray.zeros(numYPoints, numarray.Float)
974              for i in range(numXPoints):
975                  bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
976    
977              for i in range(numYPoints):
978                  bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
979    
980              # calculate the appropriate window size
981              dataWidth = max(bathXData) - min(bathXData)
982              dataHeight = max(bathYData) - min(bathYData)
983              winHeight = 600
984              winWidth = int(dataWidth*float(winHeight)/dataHeight)
985    
986              print "Done loading bathymmetry data"
987    
988              ### now do the vtk stuff
989    
990              # make sure rendering is offscreen
991              factGraphics = vtk.vtkGraphicsFactory()
992              factGraphics.SetUseMesaClasses(1)
993    
994              factImage = vtk.vtkImagingFactory()
995              factImage.SetUseMesaClasses(1)
996    
997              # make the bathymmetry colourmap
998              transFunc = vtk.vtkColorTransferFunction()
999              transFunc.SetColorSpaceToRGB()
1000              for i in range(numColours):
1001                  transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1002                  h = height[i]
1003                  while i < numColours-1 and h < height[i+1]:
1004                    h += 1
1005                    transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1006    
1007              # set up the lookup table for the wave data
1008              refLut = vtk.vtkLookupTable()
1009              self.lutTrans = vtk.vtkLookupTable()
1010              refLut.Build()
1011              alpha = 0.7   # alpha channel value
1012              for i in range(256):
1013                  (r,g,b,a) = refLut.GetTableValue(255-i)
1014                  if g == 1.0 and b < 0.5 and r < 0.5:
1015                      self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1016                  else:
1017                      self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1018    
1019              print "Generating the bathymmetry vtk object..."
1020    
1021              # set up the renderer and the render window
1022              self.ren = vtk.vtkRenderer()
1023              self.renWin = vtk.vtkRenderWindow()
1024              self.renWin.AddRenderer(self.ren)
1025              self.renWin.SetSize(winWidth, winHeight)
1026              self.renWin.OffScreenRenderingOn()
1027    
1028              # make an unstructured grid
1029              bathGrid = vtk.vtkUnstructuredGrid()
1030    
1031              # make the points for the vtk data
1032              bathPoints = vtk.vtkPoints()
1033              bathPoints.SetNumberOfPoints(numPoints)
1034    
1035              # make the vtk bathymmetry data set
1036              bathData = vtk.vtkFloatArray()
1037              bathData.SetNumberOfComponents(1)
1038              bathData.SetNumberOfTuples(numPoints)
1039    
1040              # add the points and data
1041              count = 0
1042              for i in range(numXPoints):
1043                  for j in range(numYPoints):
1044                      bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1045                      bathData.InsertTuple1(count, bathZData[j,i])
1046                      count += 1
1047            
1048              # set the data to the grid
1049              bathGrid.SetPoints(bathPoints)
1050              bathGrid.GetPointData().SetScalars(bathData)
1051    
1052              # do a delaunay on the grid
1053              bathDelaunay = vtk.vtkDelaunay2D()
1054              bathDelaunay.SetInput(bathGrid)
1055              bathDelaunay.SetTolerance(0.001)
1056              bathDelaunay.SetAlpha(2.5)
1057    
1058              # set up the mapper
1059              bathMapper = vtk.vtkDataSetMapper()
1060              bathMapper.SetInput(bathDelaunay.GetOutput())
1061              bathMapper.SetLookupTable(transFunc)
1062    
1063              # set up the actor
1064              bathActor = vtk.vtkActor()
1065              bathActor.SetMapper(bathMapper)
1066    
1067              self.ren.AddActor(bathActor)
1068    
1069              ### now add the coastline
1070              print "Loading the coastline data..."
1071    
1072              # make the coastline points
1073              coastPoints = vtk.vtkPoints()
1074    
1075              # make the coastline grid
1076              coastGrid = vtk.vtkUnstructuredGrid()
1077    
1078              # now point the points and lines into the grid
1079              totalCoastPoints = 0
1080              for polyline in self.coastline.polylines:
1081                  numPoints = len(polyline)
1082                  coastLine = vtk.vtkPolyLine()
1083                  coastLine.GetPointIds().SetNumberOfIds(numPoints)
1084                  j = 0
1085                  for point in polyline:
1086                      coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1087                      coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1088                      j += 1
1089                  coastGrid.InsertNextCell(coastLine.GetCellType(),
1090                          coastLine.GetPointIds())
1091                  totalCoastPoints += numPoints
1092    
1093              coastGrid.SetPoints(coastPoints)
1094    
1095              # make the coast's mapper
1096              coastMapper = vtk.vtkDataSetMapper()
1097              coastMapper.SetInput(coastGrid)
1098    
1099              # make its actor
1100              coastActor = vtk.vtkActor()
1101              coastActor.SetMapper(coastMapper)
1102              coastActor.GetProperty().SetColor(0,0,0)
1103    
1104              # add the actor to the renderer
1105              self.ren.AddActor(coastActor)
1106    
1107              # set up the actor for the wave
1108              self.waveActor = vtk.vtkActor()
1109    
1110              # add the actor to the renderer
1111              self.ren.AddActor(self.waveActor)
1112    
1113              print "Done loading the coastline data"
1114              
1115         def doStepPostprocessing(self, dt):         def doStepPostprocessing(self, dt):
1116          """          """
1117          Does any necessary postprocessing after each step          Does any necessary postprocessing after each step
# Line 901  class SurfMovie(Model): Line 1123  class SurfMovie(Model):
1123               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1124               # vtkobj=...               # vtkobj=...
1125               # save(self.__frame_name)               # save(self.__frame_name)
1126    
1127                 # make a reader for the data
1128                 waveReader = vtk.vtkXMLUnstructuredGridReader()
1129                 waveReader.SetFileName(self.__fn)
1130                 waveReader.Update()
1131    
1132                 # make the grid
1133                 waveGrid = waveReader.GetOutput()
1134                 waveGrid.Update()
1135    
1136                 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1137                 print "Wave height range %f - %f" % (zMin, zMax)
1138    
1139                 # make a mapper for the grid
1140                 waveMapper = vtk.vtkDataSetMapper()
1141                 waveMapper.SetInput(waveGrid)
1142                 waveMapper.SetLookupTable(self.lutTrans)
1143                 waveMapper.SetScalarRange(-self.max_height, self.max_height)
1144    
1145                 self.waveActor.SetMapper(waveMapper)
1146    
1147                 # do stuff here that can only be done on the first frame
1148                 if self.firstFrame:
1149                     # zoom in a bit
1150                     self.ren.GetActiveCamera().Zoom(2.0)
1151                     self.firstFrame = False
1152    
1153                 # render the window
1154                 self.renWin.Render()
1155    
1156                 # convert the render window to an image
1157                 win2imgFilter = vtk.vtkWindowToImageFilter()
1158                 win2imgFilter.SetInput(self.renWin)
1159    
1160                 # the image name
1161                 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1162      
1163                 # save the image to file
1164                 outWriter = vtk.vtkPNMWriter()
1165                 outWriter.SetInput(win2imgFilter.GetOutput())
1166                 outWriter.SetFileName(imgFname)
1167                 outWriter.Write()
1168                 print "Wrote %s" % imgFname
1169    
1170                 # helpful for debugging:
1171                 #os.system("display %s" % imgFname)
1172    
1173                 self.paramsFileString += "%s\n" % imgFname
1174                 self.imageFiles.append(imgFname)
1175    
1176               self.__next_t+=self.dt               self.__next_t+=self.dt
1177    
1178         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
# Line 920  class SurfMovie(Model): Line 1192  class SurfMovie(Model):
1192            image files.            image files.
1193            """            """
1194            # make the movie into self.filename            # make the movie into self.filename
1195            pass            self.paramsFileString += "END_INPUT\n"
1196              self.paramsFileString += "PIXEL HALF\n"
1197              self.paramsFileString += "ASPECT_RATIO 1\n"
1198    
1199              # write the params string to file
1200              fp = open("%s.params" % self.filename, "w")
1201              fp.write(self.paramsFileString + '\n')
1202              fp.close()
1203    
1204              print "Performing conversion to mpeg"
1205              convertString = "ppmtompeg %s.params" % self.filename
1206              result = os.system(convertString)
1207              if result != 0:
1208                  print "An error occurred in mpeg conversion"
1209    
1210              # now clean up the image files
1211              print "Removing temporary image files"
1212              os.unlink("%s.params" % self.filename)
1213              for fname in self.imageFiles:
1214                  os.unlink(fname)
1215    
1216    def main():
 if __name__=="__main__":  
1217     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
1218     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
1219    
# Line 954  if __name__=="__main__": Line 1244  if __name__=="__main__":
1244     src.end_lat=-12.     src.end_lat=-12.
1245     src.start_long=110.     src.start_long=110.
1246     src.end_long=120.     src.end_long=120.
1247     src.width=0.1     src.width=1.
1248     src.decay_zone=0.01     src.decay_zone=0.01
1249     src.amplitude=1.     src.amplitude=1.
1250    
# Line 974  if __name__=="__main__": Line 1264  if __name__=="__main__":
1264     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
1265     sm.dt=5000.     sm.dt=5000.
1266     sm.filename="movie.mpg"     sm.filename="movie.mpg"
1267       sm.max_height=Link(src,"amplitude")
1268        
1269     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1270     s.writeXML()     s.writeXML()
1271     s.run()     s.run()
1272    
1273    if __name__=="__main__":
1274        from esys.modellib import tsunami
1275        tsunami.main()
1276    
1277    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.251  
changed lines
  Added in v.281

  ViewVC Help
Powered by ViewVC 1.1.26