Changeset 566 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Timestamp:

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

Author:

ole

Message:

First stab at wind stress from file

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (2 diffs)

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        
+        
+
+