Changeset 4303

Timestamp:

Mar 14, 2007, 11:22:54 AM (18 years ago)

Author:

duncan

Message:

finishing a momentum check and adding another way of reading txt boundary files.

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (8 diffs)
• test_data_manager.py (modified) (3 diffs)
• test_shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -3030 +3030 @@
 
 
-def timefile2netcdf(filename, quantity_names = None):
+def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
     """Template for converting typical text files with time series to
     NetCDF tms file.
@@ -3043 +3043 @@
       31/08/04 00:00:00, 1.328223 0 0
       31/08/04 00:15:00, 1.292912 0 0
+
+    or time (seconds), value0 value1 value2 ...
+    
+      0.0, 1.328223 0 0
+      0.1, 1.292912 0 0
 
     will provide a time dependent function f(t) with three attributes
@@ -3064 +3069 @@
     assert len(fields) == 2, msg
 
-    try:
-        starttime = calendar.timegm(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 DataTimeError, msg
+    if not time_as_seconds:
+        try:
+            starttime = calendar.timegm(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 DataTimeError, msg
+    else:
+        try:
+            starttime = float(fields[0])
+        except Error:
+            msg = "Bad time format"
+            raise DataTimeError, msg
 
 
@@ -3102 +3114 @@
     for i, line in enumerate(lines):
         fields = line.split(',')
-        realtime = calendar.timegm(time.strptime(fields[0], time_format))
+        if not time_as_seconds:
+            realtime = calendar.timegm(time.strptime(fields[0], time_format))
+        else:
+             realtime = float(fields[0])
         T[i] = realtime - starttime
 
@@ -3126 +3141 @@
     #FIXME: Use Georef
     fid.starttime = starttime
-
 
     # dimension definitions
@@ -4866 +4880 @@
     ymomentum[slice_index] = va*h
 
-def URS2txt(basename_in):
+def urs2txt(basename_in, location_index=None):
     """
     Not finished or tested
@@ -4875 +4889 @@
     quantities = ['HA','UA','VA']
 
+    d = ","
+    
     # instanciate urs_points of the three mux files.
     mux = {}
@@ -4886 +4902 @@
     
     # Convert to utm
-    lat = a_mux.lonlatdep[:,1]
-    long = a_mux.lonlatdep[:,0]
+    latitudes = a_mux.lonlatdep[:,1]
+    longitudes = a_mux.lonlatdep[:,0]
     points_utm, zone = convert_from_latlon_to_utm( \
-        latitudes=lat, longitudes=long)
+        latitudes=latitudes, longitudes=longitudes)
     #print "points_utm", points_utm
     #print "zone", zone
-    elevations = a_mux.lonlatdep[:,2] * -1 #
+    depths = a_mux.lonlatdep[:,2]  #
     
     fid = open(basename_in + '.txt', 'w')
 
-    fid.write("zone" + str(zone) + "\n")
-
-    for elevation, point_utm in map(None, elevations, points_utm):
-        pass
+    fid.write("zone: " + str(zone) + "\n")
+
+    if location_index is not None:
+        #Title
+        li = location_index
+        fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
+                  'Northing' + "\n")
+        fid.write(str(li) +d+ str(latitudes[li])+d+ \
+              str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \
+              str(points_utm[li][01]) + "\n")
+
+    # the non-time dependent stuff
+    #Title
+    fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
+                  'Northing' +d+ 'depth m' + "\n")
+    i = 0
+    for depth, point_utm, lat, long in map(None, depths,
+                                               points_utm, latitudes,
+                                               longitudes):
+        
+        fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \
+                  str(point_utm[01]) +d+ str(depth) + "\n")
+        i +=1
+    #Time dependent
+    if location_index is not None:
+        time_step = a_mux.time_step
+        i = 0
+        #Title
+        fid.write('time' +d+ 'HA depth m'+d+ \
+                      'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \
+                      + "\n")
+        for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
+            fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \
+                      str(UA[location_index]) +d+ str(VA[location_index]) \
+                      + "\n")
+            
+            i +=1
     
 class Urs_points:

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -4939 +4939 @@
         ha=2
         ua=5
-        va=5
+        va=10
         base_name, files = self.write_mux(lat_long_points,
                                           time_step_count, time_step,
@@ -4969 +4969 @@
         elevation = fid.variables['elevation'][:]
         #assert allclose(stage[0,0], e + tide)  #Meters
-
+        #print "xmomentum", xmomentum
+        #print "ymomentum", ymomentum
         #Check the momentums - ua
-        #momentum = velocity*(stage-elevation)
-        #momentum = velocity*(stage+elevation)
-        # -(-elevation) since elevation is inverted in mux files
-        # = n*(e+tide+n) based on how I'm writing these files
-        #answer = n*(e+tide+n)
+        #momentum = velocity*water height
+        #water height = mux_depth + mux_height +tide
+        #water height = mux_depth + mux_height +tide
+        #momentum = velocity*(mux_depth + mux_height +tide)
+        #
+        
+        answer = 115
         actual = xmomentum[0,0]
-        #assert allclose(answer, actual)  #Meters
+        assert allclose(answer, actual)  #Meters^2/ sec
+        answer = 230
+        actual = ymomentum[0,0]
+        assert allclose(answer, actual)  #Meters^2/ sec
 
         # check the stage values, first time step.
@@ -5733 +5739 @@
         # extend this so it interpolates onto the boundary.
         # have it fail if there is NaN
+
+    def not_really_test_urs2txt(self):
+        # not really a test, since it doesn't check the output data
+        
+        #This will write 3 non-gridded mux files, for testing.
+        #If no quantities are passed in,
+        #na and va quantities will be the Easting values.
+        #Depth and ua will be the Northing value.
+        # this was manually checked to be correct
+        
+        tide = 1
+        time_step_count = 3
+        time_step = 2
+
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+
+        # This is gridded
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        base_name, files = self.write_mux(lat_long_points,
+                                          time_step_count, time_step)
+        urs2txt(base_name, 0)
+        print "base_name", base_name
+        
+        self.delete_mux(files)
+        #os.remove(sww_file)
+        # delete the txt file if this becomes automatic
+        
+    def daves_urs2txt(self):
+        # not really a test, since it doesn't check the output data
+        
+        #This will write 3 non-gridded mux files, for testing.
+        #If no quantities are passed in,
+        #na and va quantities will be the Easting values.
+        #Depth and ua will be the Northing value.
+        # this was manually checked to be correct
+        
+        tide = 1
+        time_step_count = 3
+        time_step = 2
+
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+
+        # This is gridded
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        base_name, files = self.write_mux(lat_long_points,
+                                          time_step_count, time_step)
+        urs2txt(base_name, 0)
+        print "base_name", base_name
+        
+        self.delete_mux(files)
+        #os.remove(sww_file)
+        # delete the txt file if this becomes automatic
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Data_Manager,'davids_test_points_urs_ungridded2sww')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2txt')
     #suite = unittest.makeSuite(Test_Data_Manager,'cache_test_URS_points_needed_and_urs_ungridded2sww')
-    #suite = unittest.makeSuite(Test_Data_Manager,'visual_test_URS_points_needed_and_urs_ungridded2sww')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_momentum')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     runner = unittest.TextTestRunner()

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1397 +1397 @@
         from data_manager import timefile2netcdf        
         timefile2netcdf(filename)
+        os.remove(filename + '.txt')
+
+        
+        #Setup wind stress
+        F = file_function(filename + '.tms', quantities = ['Attribute0',
+                                                           'Attribute1'])
+        os.remove(filename + '.tms')
+        
+
+        #print 'F(5)', F(5)
+
+        #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
+        
+        #print dir(F)
+        #print F.T
+        #print F.precomputed_values
+        #
+        #F = file_function(filename + '.txt')        
+        #
+        #print dir(F)
+        #print F.T        
+        #print F.Q
+        
+        W = Wind_stress(F)
+
+        domain.forcing_terms = []
+        domain.forcing_terms.append(W)
+
+        domain.compute_forcing_terms()
+
+        #Compute reference solution
+        const = eta_w*rho_a/rho_w
+
+        N = len(domain) # number_of_triangles
+
+        t = domain.time
+
+        s = speed(t,[1],[0])[0]
+        phi = angle(t,[1],[0])[0]
+
+        #Convert to radians
+        phi = phi*pi/180
+
+
+        #Compute velocity vector (u, v)
+        u = s*cos(phi)
+        v = s*sin(phi)
+
+        #Compute wind stress
+        S = const * sqrt(u**2 + v**2)
+
+        for k in range(N):
+            assert allclose(domain.quantities['stage'].explicit_update[k], 0)
+            assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
+            assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
+
+
+    def test_windfield_from_file_seconds(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin, sqrt
+        from anuga.config import time_format
+        from anuga.abstract_2d_finite_volumes.util import file_function
+        import time
+
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+
+        domain.time = 7 #Take a time that is represented in file (not zero)
+
+        #Write wind stress file (ensure that domain.time is covered)
+        #Take x=1 and y=0
+        filename = 'test_windstress_from_file'
+        start = time.mktime(time.strptime('2000', '%Y'))
+        fid = open(filename + '.txt', 'w')
+        dt = 0.5 #1  #One second interval
+        t = 0.0
+        while t <= 10.0:
+            fid.write('%s, %f %f\n' %(str(t),
+                                      speed(t,[1],[0])[0],
+                                      angle(t,[1],[0])[0]))
+            t += dt
+
+        fid.close()
+
+
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        from data_manager import timefile2netcdf        
+        timefile2netcdf(filename, time_as_seconds=True)
         os.remove(filename + '.txt')