/[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 1989 by caltinay, Fri Nov 7 04:11:33 2008 UTC revision 2625 by jfenwick, Fri Aug 21 06:30:25 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  import os, sys  import os, sys
23  import vtk  import vtk
24  from esys.escript import *  from esys.escript import *
25  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
26  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
27  import numarray  import numpy
28  import math  import math
29  from tempfile import mkstemp  from tempfile import mkstemp
30    
# Line 62  class GridData: Line 62  class GridData:
62              x_shape0 = x.getNumberOfDataPoints()              x_shape0 = x.getNumberOfDataPoints()
63              return_data_object = True              return_data_object = True
64          else:          else:
65              x_array = numarray.array(x)              x_array = numpy.array(x)
66              x_shape0 = x_array.shape[0]              x_shape0 = x_array.shape[0]
67              return_data_object = False              return_data_object = False
68    
69          data=numarray.zeros(x_shape0, numarray.Float)          data=numpy.zeros(x_shape0, numpy.float_)
70          ox,oy = self.getOrigin()          ox,oy = self.getOrigin()
71          dx,dy = self.getSpacing()          dx,dy = self.getSpacing()
72          data_array = self.getData()          data_array = self.getData()
# Line 74  class GridData: Line 74  class GridData:
74          i_dy = 1          i_dy = 1
75          for i in range(x_shape0):          for i in range(x_shape0):
76              if return_data_object:              if return_data_object:
77                  x_array = x.getValueOfDataPoint(i)                  x_array = x.getTupleForDataPoint(i)
78                  x_long = x_array[0]-ox                  x_long = x_array[0]-ox
79                  x_lat = x_array[1]-oy                  x_lat = x_array[1]-oy
80              else:              else:
# Line 759  class Bathymetry(Model): Line 759  class Bathymetry(Model):
759              data_grd_list.append(float(v[2]))              data_grd_list.append(float(v[2]))
760              line=f.readline().strip()              line=f.readline().strip()
761          self.trace("%s lines have been read from %s." % (len(data_grd_list), self.source))          self.trace("%s lines have been read from %s." % (len(data_grd_list), self.source))
762          data_grd=numarray.array(data_grd_list)          data_grd=numpy.array(data_grd_list)
763          x_grd=numarray.array(x_grd_list)          x_grd=numpy.array(x_grd_list)
764          y_grd=numarray.array(y_grd_list)          y_grd=numpy.array(y_grd_list)
765          if len(x_grd)<2:          if len(x_grd)<2:
766              raise ValueError,"%s: data base is too small"%str(self)              raise ValueError,"%s: data base is too small"%str(self)
767          ox=x_grd[0]          ox=x_grd[0]
# Line 802  class OceanRegion(Model): Line 802  class OceanRegion(Model):
802          Initializes the ocean region          Initializes the ocean region
803          """          """
804    
805            # create output directory if necessary
806            if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)
807    
808          print "Initializing ocean region..."          print "Initializing ocean region..."
809          if hasattr(self.source,"readline"):          if hasattr(self.source,"readline"):
810              f=self.source              f=self.source
# Line 894  class TsunamiInDeepWater(Model): Line 897  class TsunamiInDeepWater(Model):
897      Runs the deep water tsunami model based on a simplified form of the      Runs the deep water tsunami model based on a simplified form of the
898      shallow water equation.      shallow water equation.
899    
900      M{d^2 h/dt^2 =div(c grad(h)) }      *d^2 h/dt^2 =div(c grad(h))*
901    
902      where h is the wave height above sea level, and c=sqrt(g*H),      where h is the wave height above sea level, and c=sqrt(g*H),
903      with H - depth of the water level, g - gravity constant      with H - depth of the water level, g - gravity constant
# Line 918  class TsunamiInDeepWater(Model): Line 921  class TsunamiInDeepWater(Model):
921          self.__pde = LinearPDE(self.domain)          self.__pde = LinearPDE(self.domain)
922          self.__pde.setSolverMethod(self.__pde.LUMPING)          self.__pde.setSolverMethod(self.__pde.LUMPING)
923          self.__pde.setValue(D=1.)          self.__pde.setValue(D=1.)
924          self.__c2 = RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2          self.__c2 = RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numpy.pi/360.)**2
925          c_max = math.sqrt(Lsup(self.__c2))          c_max = math.sqrt(Lsup(self.__c2))
926          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))
927          if self.initial_time_step==None:          if self.initial_time_step==None:
# Line 930  class TsunamiInDeepWater(Model): Line 933  class TsunamiInDeepWater(Model):
933          """          """
934          Returns new step size          Returns new step size
935    
936          @param dt: last time step size used          :param dt: last time step size used
937          @type dt: C{float}          :type dt: ``float``
938          @return: time step size that can safely be used          :return: time step size that can safely be used
939          @rtype: C{float}          :rtype: ``float``
940          """          """
941          return self.__dt          return self.__dt
942    
# Line 941  class TsunamiInDeepWater(Model): Line 944  class TsunamiInDeepWater(Model):
944          """          """
945          Performs the time step using the Verlet scheme          Performs the time step using the Verlet scheme
946    
947          @param dt: time step size to be used          :param dt: time step size to be used
948          @type dt: C{float}          :type dt: ``float``
949          """          """
950          self.__pde.setValue(X=-self.__c2*grad(self.wave_height))          self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
951    
# Line 959  class SurfMovie(Model): Line 962  class SurfMovie(Model):
962      """      """
963      movie from a wave propagation on the sea      movie from a wave propagation on the sea
964    
965      @ivar time: current time      :ivar time: current time
966      @ivar bathymetry: scalar data set      :ivar bathymetry: scalar data set
967      @ivar wave_height: vector data set      :ivar wave_height: vector data set
968      @ivar filename: name of the movie file      :ivar filename: name of the movie file
969      """      """
970      def __init__(self,**kwargs):      def __init__(self,**kwargs):
971          Model.__init__(self,**kwargs)          Model.__init__(self,**kwargs)
# Line 986  class SurfMovie(Model): Line 989  class SurfMovie(Model):
989          self.__next_t=self.dt          self.__next_t=self.dt
990          self.__fileIndex=0          self.__fileIndex=0
991    
         # create output directory if necessary  
         if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)  
   
992          self.firstFrame = True          self.firstFrame = True
993          self.imageFiles = []          self.imageFiles = []
994    
# Line 1008  class SurfMovie(Model): Line 1008  class SurfMovie(Model):
1008          numYPoints = bathZData.shape[0]          numYPoints = bathZData.shape[0]
1009          numPoints = numXPoints*numYPoints          numPoints = numXPoints*numYPoints
1010    
1011          bathXData = numarray.zeros(numXPoints, numarray.Float)          bathXData = numpy.zeros(numXPoints, numpy.float_)
1012          bathYData = numarray.zeros(numYPoints, numarray.Float)          bathYData = numpy.zeros(numYPoints, numpy.float_)
1013          for i in range(numXPoints):          for i in range(numXPoints):
1014              bathXData[i] = bathOrigin[0] + i*bathSpacing[0]              bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
1015    
# Line 1075  class SurfMovie(Model): Line 1075  class SurfMovie(Model):
1075          numColours = len(data)          numColours = len(data)
1076    
1077          # convert the colourmap into something vtk is more happy with          # convert the colourmap into something vtk is more happy with
1078          height = numarray.zeros(numColours, numarray.Float)          height = numpy.zeros(numColours, numpy.float_)
1079          red = numarray.zeros(numColours, numarray.Float)          red = numpy.zeros(numColours, numpy.float_)
1080          green = numarray.zeros(numColours, numarray.Float)          green = numpy.zeros(numColours, numpy.float_)
1081          blue = numarray.zeros(numColours, numarray.Float)          blue = numpy.zeros(numColours, numpy.float_)
1082          for i in range(numColours):          for i in range(numColours):
1083              (h, r, g, b) = data[i]              (h, r, g, b) = data[i]
1084              height[i] = float(h)              height[i] = float(h)
# Line 1229  class SurfMovie(Model): Line 1229  class SurfMovie(Model):
1229          """          """
1230          Does any necessary postprocessing after each step          Does any necessary postprocessing after each step
1231    
1232          @param dt:          :param dt:
1233          """          """
1234          if self.t >= self.__next_t:          if self.t >= self.__next_t:
1235              prefix = os.path.join(WORKDIR, \              prefix = os.path.join(WORKDIR, \
# Line 1250  class SurfMovie(Model): Line 1250  class SurfMovie(Model):
1250          """          """
1251          returns new step size          returns new step size
1252    
1253          @param dt: last time step size used          :param dt: last time step size used
1254          @type dt: C{float}          :type dt: ``float``
1255          @return: time step size that can savely be used          :return: time step size that can savely be used
1256          @rtype: C{float}          :rtype: ``float``
1257          """          """
1258          return self.__next_t-self.t          return self.__next_t-self.t
1259    
# Line 1359  def main(): Line 1359  def main():
1359      ts.bathymetry = Link(oreg, "bathymetry")      ts.bathymetry = Link(oreg, "bathymetry")
1360    
1361      sq = Sequencer()      sq = Sequencer()
1362      sq.t_end = 20000.      sq.t_end = 30000.
1363    
1364      sm = SurfMovie()      sm = SurfMovie()
1365      sm.bathymetry = Link(b, "bathymetry")      sm.bathymetry = Link(b, "bathymetry")

Legend:
Removed from v.1989  
changed lines
  Added in v.2625

  ViewVC Help
Powered by ViewVC 1.1.26