Changeset 566

Timestamp:

Nov 16, 2004, 4:52:21 PM (19 years ago)

Author:

ole

Message:

First stab at wind stress from file

Location:

inundation/ga/storm_surge/pyvolution

Files:

• config.py (modified) (1 diff)
• domain.py (modified) (1 diff)
• generic_boundary_conditions.py (modified) (5 diffs)
• pressure_force.py (modified) (2 diffs)
• shallow_water.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/config.py

@@ -40 +40 @@
 
 
+eta_w = 3.0e-3 #Wind stress coefficient
+rho_a = 1.2e-3 #Atmospheric density
+rho_w = 1023   #Fluid density [kg/m^3] (rho_w = 1023 for salt water)
+
+        
+
+
+
 
 beta = 0.9  #]0;1] (e.g. 0.9). 1 allows the steepest gradients,

inundation/ga/storm_surge/pyvolution/domain.py

@@ -64 +64 @@
         self.time = 0.0
         self.finaltime = None
-        self.min_timestep = self.max_timestep = 0.0        
+        self.min_timestep = self.max_timestep = 0.0
+        self.starttime = None
         
         #Checkpointing          

inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

@@ -143 +143 @@
         fid.close()
         fields = line.split(',')
-        msg = 'Boundary file %s must have the format date, values'
+        msg = 'File %s must have the format date, values'
         assert len(fields) == 2, msg
 
         try:
-            start_time = time.mktime(time.strptime(fields[0], time_format))
+            starttime = time.mktime(time.strptime(fields[0], time_format))
         except ValueError:    
-            msg = 'First field in boundary file %s must be' %filename
+            msg = 'First field in file %s must be' %filename
             msg += ' date-time with format %s.\n' %time_format
             msg += 'I got %s instead.' %fields[0] 
@@ -170 +170 @@
         self.filename = filename
         self.domain = domain
-        self.start_time = start_time
-
+        self.starttime = starttime
+
+        if domain.starttime is None:
+            domain.starttime = starttime
+        else:
+            msg = 'Start time as specified in domain (%s) is earlier '
+            'than the starttime of file %s: %s'\
+                  %(domain.starttime, self.filename, self.starttime) 
+            if self.starttime > domain.starttime:
+                raise msg
+        
         self.read_time_boundary() #Now read all times in.
         
@@ -195 +204 @@
             real_time = time.mktime(time.strptime(fields[0], time_format))
 
-            T[i] = real_time - self.start_time
+            T[i] = real_time - self.starttime
 
             
@@ -201 +210 @@
                 Q[i, j] = float(value)
 
-        msg = 'Time boundary must list time as a monotonuosly increasing sequence'          
+        msg = 'Time boundary must list time as a monotonuosly '
+        msg += 'increasing sequence'
+        
         assert alltrue( T[1:] - T[:-1] > 0 ), msg
 
@@ -233 +244 @@
         
         #t is now between index and index+1
-        ratio = (t - self.T[self.index])/(self.T[self.index+1] - self.T[self.index])
-
-        return self.Q[self.index,:] + ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+        ratio = (t - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+        return self.Q[self.index,:] +\
+               ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
 
 

inundation/ga/storm_surge/pyvolution/pressure_force.py

@@ -48 +48 @@
 
 def cyclone(domain):
-    rho = 1023  #Fluid density [kg/m^3] #FIXME: Move to config object
+    from config import rho_w 
     from util import gradient
 
@@ -72 +72 @@
 
         #Update momentum
-        xmom[k] += -px*avg_h/rho
-        ymom[k] += -py*avg_h/rho        
+        xmom[k] += -px*avg_h/rho_w
+        ymom[k] += -py*avg_h/rho_w        
 
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -937 +937 @@
         """
 
+        from config import rho_a, rho_w, eta_w
+        
         self.s = s
         self.phi = phi
 
-        #FIXME: Maybe move to config or read from domain
-        Cw = 3.0e-3 #Wind stress coefficient
-        rho_a = 1.2e-3 #Atmospheric density
-        rho = 1023 #Density of water
-        
-        self.const = Cw*rho_a/rho
+        self.const = eta_w*rho_a/rho_w
         
 
@@ -1007 +1004 @@
         
         
+
+class Wind_stress_from_file:
+    """Apply wind stress read from a file to water momentum
+
+    At this stage the wind field is assumed to depend on time only, i.e
+    no spatial dependency.
+    
+    FIXME: This class may be incorporated in the generic Wind_stress class
+    """
+
+    def __init__(self, filename):
+        """Initialise windfield from a file with the following format
+
+        time [DD/MM/YY hh:mm:ss], speed [m/s] direction [degrees]
+        
+        e.g.
+        
+        03/09/04 19:15:00, 9.53 39
+        
+        domain.forcing_terms.append(W)
+        
+        """
+
+        import time
+        from Numeric import array
+        from config import time_format
+
+
+        from config import rho_a, rho_w, eta_w
+        self.const = eta_w*rho_a/rho_w
+        
+
+        try:
+            fid = open(filename)
+        except Exception, e:
+            msg = 'Wind stress file %s could not be opened: %s\n'\
+                  %(filename, e)
+            raise msg
+
+
+        line = fid.readline()
+        fid.close()
+        fields = line.split(',')
+        msg = 'File %s must have the format date, values'
+        assert len(fields) == 2, msg
+
+        try:
+            starttime = time.mktime(time.strptime(fields[0], time_format))
+        except ValueError:    
+            msg = 'First field in file %s must be' %filename
+            msg += ' date-time with format %s.\n' %time_format
+            msg += 'I got %s instead.' %fields[0] 
+            raise msg
+
+        #Split values
+        values = []
+        for value in fields[1].split():
+            values.append(float(value))
+
+        q = array(values)
+        
+        msg = 'ERROR: File function must return a 1d list or array '
+        assert len(q.shape) == 1, msg
+            
+
+        msg = 'Values specified in file must convert to an array of length 2'
+        assert len(q) == 2, msg
+
+        self.filename = filename
+        self.domain = domain
+        self.starttime = starttime
+
+
+        if domain.starttime is None:
+            domain.starttime = starttime
+        else:
+            msg = 'Start time as specified in domain (%s) is earlier '
+            'than the starttime of file %s: %s'\
+                  %(domain.starttime, self.filename, self.starttime) 
+            if self.starttime > domain.starttime:
+                raise msg
+
+
+        self.read_times() #Now read all times in.
+
+        
+    def read_times(self):
+        from Numeric import zeros, Float, alltrue
+        from config import time_format
+        import time
+        
+        fid = open(self.filename)        
+        lines = fid.readlines()
+        fid.close()
+        
+        N = len(lines)
+
+        T = zeros(N, Float)       #Time
+        Q = zeros((N, 2), Float)  #Wind field (s, phi)
+        
+        for i, line in enumerate(lines):
+            fields = line.split(',')
+            real_time = time.mktime(time.strptime(fields[0], time_format))
+
+            T[i] = real_time - self.start_time
+
+            for j, value in enumerate(fields[1].split()):
+                Q[i, j] = float(value)
+
+        msg = 'File %s must list time as a monotonuosly ' %self.filename
+        msg += 'increasing sequence'
+        assert alltrue( T[1:] - T[:-1] > 0 ), msg
+
+        self.T = T     #Time
+        self.Q = Q     #Windfied
+        self.index = 0 #Initial index
+
+        
+    def __repr__(self):
+        return 'Wind field from file'
+
+        
+
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain
+        """
+
+        from math import pi, cos, sin, sqrt
+        
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update    
+        
+        N = domain.number_of_elements
+        t = self.domain.time
+        
+        msg = 'Time given in File %s does not match model time'\
+              %self.filename
+        
+        if t < self.T[0]: raise msg
+        if t > self.T[-1]: raise msg        
+
+        oldindex = self.index
+
+        #Find slot
+        while t > self.T[self.index]: self.index += 1
+        while t < self.T[self.index]: self.index -= 1
+
+        #t is now between index and index+1
+        ratio = (t - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+        #Compute interpolated values for s and phi
+        q = self.Q[self.index,:] +\
+            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+
+        s = q[0]
+        
+        #Convert to radians
+        phi = q[1]*pi/180
+        
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress for this time step
+        S = self.const * sqrt(u**2 + v**2)
+        xmom_update += S*u
+        ymom_update += S*v        
+        
+
+