Changeset 8488

Timestamp:

Aug 1, 2012, 3:45:49 PM (13 years ago)

Author:

steve

Message:

Added in some more unit test forbasic rate operator

Location:

trunk/anuga_core/source/anuga

Files:

• file/test_sww.py (modified) (1 diff)
• operators/rain_operators.py (modified) (1 diff)
• operators/rate_operators.py (modified) (5 diffs)
• operators/test_rate_operators.py (modified) (9 diffs)
• structures/inlet_operator.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/file/test_sww.py

@@ -60 +60 @@
         Bt = Transmissive_boundary(domain)
         Bd = Dirichlet_boundary([0.2,0.,0.])
-        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
+        Bw = Time_boundary(domain=domain,function=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
 
         #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})

trunk/anuga_core/source/anuga/operators/rain_operators.py

@@ -12 +12 @@
 from anuga import Quantity
 import numpy as num
+from warnings import warn
 import anuga.utilities.log as log
 

trunk/anuga_core/source/anuga/operators/rate_operators.py

@@ -14 +14 @@
 from anuga import indent
 import numpy as num
+from warnings import warn
 import anuga.utilities.log as log
 
@@ -51 +52 @@
         # Local variables
         #------------------------------------------
-        self.rate = rate
         self.indices = indices
         self.factor = factor
 
-        #------------------------------------------
-        # Check and store default_rate
-        #------------------------------------------
-        msg = ('Keyword argument default_rate must be either None '
-               'or a function of time.\nI got %s.' % str(default_rate))
-        assert (default_rate is None or
-                isinstance(default_rate, (int, float)) or
-                callable(default_rate)), msg
-
-
-        #------------------------------------------
-        # Allow application longer than data
-        #------------------------------------------
-        if default_rate is not None:
-            # If it is a constant, make it a function
-            if not callable(default_rate):
-                tmp = default_rate
-                default_rate = lambda t: tmp
-
-            # Check that default_rate is a function of one argument
-            try:
-                default_rate(0.0)
-            except:
-                raise Exception(msg)
-
-        self.default_rate = default_rate
+        self.set_rate(rate)
+        self.set_default_rate(default_rate)
+
         self.default_rate_invoked = False    # Flag
-
-
-
 
     def __call__(self):
@@ -109 +83 @@
                          % (self.quantity_name, domain.get_time(), rate))
 
-        if self.indices is None:
-            self.stage_c[:] = self.stage_c[:]  \
-                   + factor*rate*timestep
-        else:
-            self.stage_c[indices] = self.stage_c[indices] \
-                   + factor*rate*timestep
-
+        if rate >= 0.0:
+            if self.indices is None:
+                self.stage_c[:] = self.stage_c[:]  \
+                       + factor*rate*timestep
+            else:
+                self.stage_c[indices] = self.stage_c[indices] \
+                       + factor*rate*timestep
+        else: # Be more careful if rate < 0
+            if self.indices is None:
+                self.stage_c[:] = num.maximun(self.stage_c  \
+                       + factor*rate*timestep, self.elev_c )
+            else:
+                self.stage_c[indices] = num.maximum(self.stage_c[indices] \
+                       + factor*rate*timestep, self.elev_c[indices])
 
     def get_rate(self, t=None):
@@ -142 +123 @@
                     if self.default_rate_invoked is False:
                         # Issue warning the first time
-                        msg = ('%s\n'
+                        msg = ('\n%s'
                            'Instead I will use the default rate: %s\n'
                            'Note: Further warnings will be supressed'
-                           % (str(e), str(self.default_rate)))
+                           % (str(e), self.default_rate(t)))
                         warn(msg)
 
@@ -162 +143 @@
 
         return rate
+
+
+    def set_rate(self, rate):
+        """Ability to change rate while running
+        """
+
+        #------------------------------------------
+        # Check and store rate
+        #------------------------------------------
+        msg = ('Keyword argument rate must be a'
+               'scalar, or a function of time.')
+        assert (isinstance(rate, (int, float)) or
+                callable(rate)), msg
+        self.rate = rate
+
+
+    def set_default_rate(self, default_rate):
+        """
+        Check and store default_rate
+        """
+        msg = ('Default_rate must be either None '
+               'a scalar, or a function of time.\nI got %s.' % str(default_rate))
+        assert (default_rate is None or
+                isinstance(default_rate, (int, float)) or
+                callable(default_rate)), msg
+
+
+        #------------------------------------------
+        # Allow longer than data
+        #------------------------------------------
+        if default_rate is not None:
+            # If it is a constant, make it a function
+            if not callable(default_rate):
+                tmp = default_rate
+                default_rate = lambda t: tmp
+
+            # Check that default_rate is a function of one argument
+            try:
+                default_rate(0.0)
+            except:
+                raise Exception(msg)
+
+        self.default_rate = default_rate
+
 
     def parallel_safe(self):

trunk/anuga_core/source/anuga/operators/test_rate_operators.py

@@ -4 +4 @@
 import unittest, os
 import anuga
-from anuga.shallow_water.shallow_water_domain import Domain
-from boundaries import Reflective_boundary
-from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga import Domain
+from anuga import Reflective_boundary
+from anuga import rectangular_cross_domain
+from anuga import file_function
+
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.file_conversion.file_conversion import timefile2netcdf
-from forcing import *
-from mesh_factory import rectangular
-#from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-from anuga.abstract_2d_finite_volumes.util import file_function
-from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.config import time_format
+
+from rate_operators import *
 
 import numpy as num
 import warnings
-
-
-def scalar_func_list(t, x, y):
-    """Function that returns a scalar.
-
-    Used to test error message when numeric array is expected
-    """
-
-    return [17.7]
-
-
-def speed(t, x, y):
-    """
-    Variable windfield implemented using functions    
-    Large speeds halfway between center and edges
-
-    Low speeds at center and edges
-    """
-
-    from math import exp, cos, pi
-
-    x = num.array(x)
-    y = num.array(y)
-
-    N = len(x)
-    s = 0*x  #New array
-
-    for k in range(N):
-        r = num.sqrt(x[k]**2 + y[k]**2)
-        factor = exp(-(r-0.15)**2)
-        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
-
-    return s
-
-
-def angle(t, x, y):
-    """Rotating field
-    """
-    from math import atan, pi
-
-    x = num.array(x)
-    y = num.array(y)
-
-    N = len(x)
-    a = 0 * x    # New array
-
-    for k in range(N):
-        r = num.sqrt(x[k]**2 + y[k]**2)
-
-        angle = atan(y[k]/x[k])
-
-        if x[k] < 0:
-            angle += pi
-
-        # Take normal direction
-        angle -= pi/2
-
-        # Ensure positive radians
-        if angle < 0:
-            angle += 2*pi
-
-        a[k] = angle/pi*180
-
-    return a
-
-def time_varying_speed(t, x, y):
-    """
-    Variable speed windfield
-    """
-
-    from math import exp, cos, pi
-
-    x = num.array(x,num.float)
-    y = num.array(y,num.float)
-
-    N = len(x)
-    s = 0*x  #New array
-
-    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
-    S=100.
-    for k in range(N):
-        s[k]=S*(1.+t/100.)
-    return s
-
-
-def time_varying_angle(t, x, y):
-    """Rotating field
-    """
-    from math import atan, pi
-
-    x = num.array(x,num.float)
-    y = num.array(y,num.float)
-
-    N = len(x)
-    a = 0 * x    # New array
-
-    phi=135.
-    for k in range(N):
-        a[k]=phi*(1.+t/100.)
-
-    return a
-
-
-def time_varying_pressure(t, x, y):
-    """Rotating field
-    """
-    from math import atan, pi
-
-    x = num.array(x,num.float)
-    y = num.array(y,num.float)
-
-    N = len(x)
-    p = 0 * x    # New array
-
-    p0=1000.
-    for k in range(N):
-        p[k]=p0*(1.-t/100.)
-
-    return p
-
-def spatial_linear_varying_speed(t, x, y):
-    """
-    Variable speed windfield
-    """
-
-    from math import exp, cos, pi
-
-    x = num.array(x)
-    y = num.array(y)
-
-    N = len(x)
-    s = 0*x  #New array
-
-    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
-    s0=250.
-    ymin=num.min(y)
-    xmin=num.min(x)
-    a=0.000025; b=0.0000125;
-    for k in range(N):
-        s[k]=s0*(1+t/100.)+a*x[k]+b*y[k]
-    return s
-
-
-def spatial_linear_varying_angle(t, x, y):
-    """Rotating field
-    """
-    from math import atan, pi
-
-    x = num.array(x)
-    y = num.array(y)
-
-    N = len(x)
-    a = 0 * x    # New array
-
-    phi=135.
-    b1=0.000025; b2=0.00001125;
-    for k in range(N):
-        a[k]=phi*(1+t/100.)+b1*x[k]+b2*y[k]
-    return a
-
-def spatial_linear_varying_pressure(t, x, y):
-    p0=1000;
-    a=0.000025; b=0.0000125;
-
-    x = num.array(x)
-    y = num.array(y)
-
-    N = len(x)
-    p = 0 * x    # New array
-
-    for k in range(N):
-        p[k]=p0*(1.-t/100.)+a*x[k]+b*y[k]
-    return p
-
-
-def grid_1d(x0,dx,nx):
-    x = num.empty(nx,dtype=num.float)
-    for i in range(nx):
-        x[i]=x0+float(i)*dx
-    return x
-    
-
-def ndgrid(x,y):
-    nx = len(x)
-    ny = len(y)
-    X = num.empty(nx*ny,dtype=num.float)
-    Y = num.empty(nx*ny,dtype=num.float)
-    k=0
-    for i in range(nx):
-        for j in range(ny):
-            X[k]=x[i]
-            Y[k]=y[j]
-            k+=1
-    return X,Y
-
-class Test_Forcing(unittest.TestCase):
+import time
+
+
+
+class Test_rate_operators(unittest.TestCase):
     def setUp(self):
         pass
@@ -217 +27 @@
     def tearDown(self):
         pass
-        
-    def write_wind_pressure_field_sts(self,
-                                      field_sts_filename,
-                                      nrows=10,
-                                      ncols=10,
-                                      cellsize=25,
-                                      origin=(0.0,0.0),
-                                      refzone=50,
-                                      timestep=1,
-                                      number_of_timesteps=10,
-                                      angle=135.0,
-                                      speed=100.0,
-                                      pressure=1000.0):
-
-        xllcorner=origin[0]
-        yllcorner=origin[1]
-        starttime = 0; endtime = number_of_timesteps*timestep;
-        no_data = -9999
-
-        time = num.arange(starttime, endtime, timestep, dtype='i')
-
-        x = grid_1d(xllcorner,cellsize,ncols)
-        y = grid_1d(yllcorner,cellsize,nrows)
-        [X,Y] = ndgrid(x,y)
-        number_of_points = nrows*ncols
-
-        wind_speed = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
-        wind_angle = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
-        barometric_pressure = num.empty((number_of_timesteps,nrows*ncols),
-                                        dtype=num.float)
-
-        if ( callable(speed) and callable(angle) and callable(pressure) ):
-            x = num.ones(3, num.float)
-            y = num.ones(3, num.float)
-            try:
-                s = speed(1.0, x=x, y=y)
-                a = angle(1.0, x=x, y=y)
-                p = pressure(1.0, x=x, y=y)
-                use_function=True
-            except Exception, e:
-                msg = 'Function could not be executed.\n'
-                raise Exception, msg
-        else:
-            try :
-                speed=float(speed)
-                angle=float(angle)
-                pressure=float(pressure)
-                use_function=False
-            except:
-                msg = ('Force fields must be a scalar value coercible to float.')
-                raise Exception, msg
-
-        for i,t in enumerate(time):
-            if ( use_function ):
-                wind_speed[i,:] = speed(t,X,Y)
-                wind_angle[i,:] = angle(t,X,Y)
-                barometric_pressure[i,:] = pressure(t,X,Y)
-            else:
-                wind_speed[i,:] = speed
-                wind_angle[i,:] = angle
-                barometric_pressure[i,:] = pressure
-
-        # "Creating the field STS NetCDF file"
-
-        fid = NetCDFFile(field_sts_filename+'.sts', 'w')
-        fid.institution = 'Geoscience Australia'
-        fid.description = "description"
-        fid.starttime = 0.0
-        fid.ncols = ncols
-        fid.nrows = nrows
-        fid.cellsize = cellsize
-        fid.no_data = no_data
-        fid.createDimension('number_of_points', number_of_points)
-        fid.createDimension('number_of_timesteps', number_of_timesteps)
-        fid.createDimension('numbers_in_range', 2)
-
-        fid.createVariable('x', 'd', ('number_of_points',))
-        fid.createVariable('y', 'd', ('number_of_points',))
-        fid.createVariable('time', 'i', ('number_of_timesteps',))
-        fid.createVariable('wind_speed', 'd', ('number_of_timesteps', 
-                                               'number_of_points'))
-        fid.createVariable('wind_speed_range', 'd', ('numbers_in_range', ))
-        fid.createVariable('wind_angle', 'd', ('number_of_timesteps', 
-                                               'number_of_points'))
-        fid.createVariable('wind_angle_range', 'd', ('numbers_in_range',))
-        fid.createVariable('barometric_pressure', 'd', ('number_of_timesteps', 
-                                             'number_of_points'))
-        fid.createVariable('barometric_pressure_range', 'd', ('numbers_in_range',))
-
-
-        fid.variables['wind_speed_range'][:] = num.array([1e+036, -1e+036])
-        fid.variables['wind_angle_range'][:] = num.array([1e+036, -1e+036])
-        fid.variables['barometric_pressure_range'][:] = num.array([1e+036, -1e+036])
-        fid.variables['time'][:] = time
-
-        ws = fid.variables['wind_speed']
-        wa = fid.variables['wind_angle']
-        pr = fid.variables['barometric_pressure']
-
-        for i in xrange(number_of_timesteps):
-            ws[i] = wind_speed[i,:]
-            wa[i] = wind_angle[i,:]
-            pr[i] = barometric_pressure[i,:]
-
-        origin = anuga.coordinate_transforms.geo_reference.Geo_reference(refzone,
-                                                                         xllcorner,
-                                                                         yllcorner)
-        geo_ref = anuga.coordinate_transforms.geo_reference.write_NetCDF_georeference(origin, fid)
-
-        fid.variables['x'][:]=X-geo_ref.get_xllcorner()
-        fid.variables['y'][:]=Y-geo_ref.get_yllcorner()
-
-
-        fid.close()
-
-    def test_constant_wind_stress(self):
+
+
+    def test_rate_operator_simple(self):
         from anuga.config import rho_a, rho_w, eta_w
         from math import pi, cos, sin
@@ -357 +54 @@
         domain.set_boundary({'exterior': Br})
 
-        #Setup only one forcing term, constant wind stress
-        s = 100
-        phi = 135
-        domain.forcing_terms = []
-        domain.forcing_terms.append(Wind_stress(s, phi))
-
-        domain.compute_forcing_terms()
-
-        const = eta_w*rho_a / rho_w
-
-        #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 * num.sqrt(u**2 + v**2)
-
-        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
-
-    def test_variable_wind_stress(self):
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+
+        rate = 1.0
+        factor = 10.0
+        default_rate= 0.0
+
+        operator = Rate_operator(domain, rate=rate, factor=factor, \
+                      indices=indices, default_rate = default_rate)
+        
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        stage_ex = [ 21.,  21.,   1.,  21.]
+
+
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+ 
+    def test_rate_operator_negative_rate(self):
         from anuga.config import rho_a, rho_w, eta_w
         from math import pi, cos, sin
@@ -406 +110 @@
         domain.set_boundary({'exterior': Br})
 
-        domain.time = 5.54    # Take a random time (not zero)
-
-        #Setup only one forcing term, constant wind stress
-        s = 100
-        phi = 135
-        domain.forcing_terms = []
-        domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
-
-        domain.compute_forcing_terms()
-
-        #Compute reference solution
-        const = eta_w*rho_a / rho_w
-
-        N = len(domain)    # number_of_triangles
-
-        xc = domain.get_centroid_coordinates()
-        t = domain.time
-
-        x = xc[:,0]
-        y = xc[:,1]
-        s_vec = speed(t,x,y)
-        phi_vec = angle(t,x,y)
-
-        for k in range(N):
-            # Convert to radians
-            phi = phi_vec[k]*pi / 180
-            s = s_vec[k]
-
-            # Compute velocity vector (u, v)
-            u = s*cos(phi)
-            v = s*sin(phi)
-
-            # Compute wind stress
-            S = const * num.sqrt(u**2 + v**2)
-
-            assert num.allclose(domain.quantities['stage'].explicit_update[k],
-                                0)
-            assert num.allclose(domain.quantities['xmomentum'].\
-                                     explicit_update[k],
-                                S*u)
-            assert num.allclose(domain.quantities['ymomentum'].\
-                                     explicit_update[k],
-                                S*v)
-
-    def test_windfield_from_file(self):
-        import time
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+
+
+
+        #Catchment_Rain_Polygon = read_polygon(join('CatchmentBdy.csv'))
+        #rainfall = file_function(join('1y120m.tms'), quantities=['rainfall'])
+        rate = -1.0
+        factor = 10.0
+        default_rate= 0.0
+
+        operator = Rate_operator(domain, rate=rate, factor=factor, \
+                      indices=indices, default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        stage_ex = [ 0.,  0.,   1.,  0.]
+
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+
+    def test_rate_operator_rate_from_file(self):
         from anuga.config import rho_a, rho_w, eta_w
         from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.abstract_2d_finite_volumes.util import file_function
 
         a = [0.0, 0.0]
@@ -468 +161 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
+
+        #---------------------------------
+        #Typical ASCII file
+        #---------------------------------
+        finaltime = 1200
+        filename = 'test_file_function'
+        fid = open(filename + '.txt', 'w')
+        start = time.mktime(time.strptime('2000', '%Y'))
+        dt = 60  #One minute intervals
+        t = 0.0
+        while t <= finaltime:
+            t_string = time.strftime(time_format, time.gmtime(t+start))
+            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
+            t += dt
+
+        fid.close()
+
+        #Convert ASCII file to NetCDF (Which is what we really like!)
+        timefile2netcdf(filename+'.txt')
+
+
+        #Create file function from time series
+        F = file_function(filename + '.tms',
+                          quantities = ['Attribute0',
+                                        'Attribute1',
+                                        'Attribute2'])
+
+        #Now try interpolation
+        for i in range(20):
+            t = i*10
+            q = F(t)
+
+            #Exact linear intpolation
+            assert num.allclose(q[0], 2*t)
+            if i%6 == 0:
+                assert num.allclose(q[1], t**2)
+                assert num.allclose(q[2], sin(t*pi/600))
+
+        #Check non-exact
+
+        t = 90 #Halfway between 60 and 120
+        q = F(t)
+        assert num.allclose( (120**2 + 60**2)/2, q[1] )
+        assert num.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
+
+
+        t = 100 #Two thirds of the way between between 60 and 120
+        q = F(t)
+        assert num.allclose( 2*120**2/3 + 60**2/3, q[1] )
+        assert num.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
+
+        #os.remove(filename + '.txt')
+        #os.remove(filename + '.tms')
+
+
         domain = Domain(points, vertices)
 
-        # Flat surface with 1m of water
+        #Flat surface with 1m of water
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
@@ -478 +226 @@
         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 = 1    # One second interval
-        t = 0.0
-        while t <= 10.0:
-            t_string = time.strftime(time_format, time.gmtime(t+start))
-
-            fid.write('%s, %f %f\n' %
-                      (t_string, speed(t,[1],[0])[0], angle(t,[1],[0])[0]))
-            t += dt
-
-        fid.close()
-
-        timefile2netcdf(filename + '.txt')
-        os.remove(filename + '.txt')
-
-        # Setup wind stress
-        F = file_function(filename + '.tms',
-                          quantities=['Attribute0', 'Attribute1'])
-        os.remove(filename + '.tms')
-
-        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 * num.sqrt(u**2 + v**2)
-
-        for k in range(N):
-            assert num.allclose(domain.quantities['stage'].explicit_update[k],
-                                0)
-            assert num.allclose(domain.quantities['xmomentum'].\
-                                    explicit_update[k],
-                                S*u)
-            assert num.allclose(domain.quantities['ymomentum'].\
-                                    explicit_update[k],
-                                S*v)
-
-    def test_windfield_from_file_seconds(self):
-        import time
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+
+
+        rate = file_function('test_file_function.tms', quantities=['Attribute1'])
+
+        factor = 1000.0
+        default_rate= 17.7
+
+        operator = Rate_operator(domain, rate=rate, factor=factor, \
+                      indices=indices, default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.set_starttime(360.0)
+        domain.timestep = 1.0
+
+        operator()
+
+
+        d = domain.get_time()**2 * factor + 1.0
+        stage_ex0 = [ d,  d,   1.,  d]
+
+#        print d, domain.get_time(), F(360.0)
+
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex0)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+
+        domain.set_starttime(-10.0)
+        domain.timestep = 1.0
+
+        try:
+            operator()
+        except:
+            pass
+        else:
+            raise Exception('Should have raised an exception, time too early')
+
+
+        domain.set_starttime(1300.0)
+        domain.timestep = 1.0
+
+        operator()
+
+        d = default_rate*factor + d
+        stage_ex1 = [ d,  d,   1.,  d]
+
+#        print domain.quantities['elevation'].centroid_values
+#        print domain.quantities['stage'].centroid_values
+#        print domain.quantities['xmomentum'].centroid_values
+#        print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex1)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+
+    def test_rate_operator_functions_rate_default_rate(self):
         from anuga.config import rho_a, rho_w, eta_w
         from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.abstract_2d_finite_volumes.util import file_function
 
         a = [0.0, 0.0]
@@ -561 +312 @@
         domain = Domain(points, vertices)
 
-        # Flat surface with 1m of water
+        #Flat surface with 1m of water
         domain.set_quantity('elevation', 0)
         domain.set_quantity('stage', 1.0)
@@ -569 +320 @@
         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    # Half 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()
-
-        timefile2netcdf(filename + '.txt', time_as_seconds=True)
-        os.remove(filename + '.txt')
-
-        # Setup wind stress
-        F = file_function(filename + '.tms',
-                          quantities=['Attribute0', 'Attribute1'])
-        os.remove(filename + '.tms')
-
-        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 * num.sqrt(u**2 + v**2)
-
-        for k in range(N):
-            assert num.allclose(domain.quantities['stage'].explicit_update[k],
-                                0)
-            assert num.allclose(domain.quantities['xmomentum'].\
-                                    explicit_update[k],
-                                S*u)
-            assert num.allclose(domain.quantities['ymomentum'].\
-                                    explicit_update[k],
-                                S*v)
-
-    def test_wind_stress_error_condition(self):
-        """Test that windstress reacts properly when forcing functions
-        are wrong - e.g. returns a scalar
-        """
-
-        from math import pi, cos, sin
-        from anuga.config import rho_a, rho_w, eta_w
-
-        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 = 5.54    # Take a random time (not zero)
-
-        # Setup only one forcing term, bad func
-        domain.forcing_terms = []
-
-        try:
-            domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
-                                                    phi=angle))
-        except AssertionError:
-            pass
-        else:
-            msg = 'Should have raised exception'
-            raise Exception, msg
-
-        try:
-            domain.forcing_terms.append(Wind_stress(s=speed, phi=scalar_func))
-        except Exception:
-            pass
-        else:
-            msg = 'Should have raised exception'
-            raise Exception, msg
-
-        try:
-            domain.forcing_terms.append(Wind_stress(s=speed, phi='xx'))
-        except:
-            pass
-        else:
-            msg = 'Should have raised exception'
-            raise Exception, msg
-
-    def test_rainfall(self):
-        from math import pi, cos, sin
-
-        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})
-
-        # Setup only one forcing term, constant rainfall
-        domain.forcing_terms = []
-        domain.forcing_terms.append(Rainfall(domain, rate=2.0))
-
-        domain.compute_forcing_terms()
-        assert num.allclose(domain.quantities['stage'].explicit_update,
-                            2.0/1000)
-
-    def test_rainfall_restricted_by_polygon(self):
-        from math import pi, cos, sin
-
-        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})
-
-        # Setup only one forcing term, constant rainfall
-        # restricted to a polygon enclosing triangle #1 (bce)
-        domain.forcing_terms = []
-        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
-
-        assert num.allclose(R.exchange_area, 2)
-
-        domain.forcing_terms.append(R)
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            2.0/1000)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-    def test_time_dependent_rainfall_restricted_by_polygon(self):
-        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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # restricted to a polygon enclosing triangle #1 (bce)
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=lambda t: 3*t + 7,
-                     polygon = [[1,1], [2,1], [2,2], [1,2]])
-
-        assert num.allclose(R.exchange_area, 2)
+        verbose = False
         
-        domain.forcing_terms.append(R)
-
-        domain.time = 10.
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.time + 7)/1000)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-    def test_time_dependent_rainfall_using_starttime(self):
-        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
-
-        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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # restricted to a polygon enclosing triangle #1 (bce)
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=lambda t: 3*t + 7,
-                     polygon=rainfall_poly)                     
-
-        assert num.allclose(R.exchange_area, 2)
-        
-        domain.forcing_terms.append(R)
-
-        # This will test that time is set to starttime in set_starttime
-        domain.set_starttime(5.0)
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.get_time() + 7)/1000)
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.get_starttime() + 7)/1000)
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-    def test_time_dependent_rainfall_using_georef(self):
-        """test_time_dependent_rainfall_using_georef
-
-        This will also test the General forcing term using georef
-        """
-
-        # Mesh in zone 56 (absolute coords)
-        x0 = 314036.58727982
-        y0 = 6224951.2960092
-
-        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
-        rainfall_poly += [x0, y0]
-
-        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,
-                        geo_reference=Geo_reference(56, x0, y0))
-
-        # 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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # restricted to a polygon enclosing triangle #1 (bce)
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=lambda t: 3*t + 7,
-                     polygon=rainfall_poly)
-
-        assert num.allclose(R.exchange_area, 2)
-        
-        domain.forcing_terms.append(R)
-
-        # This will test that time is set to starttime in set_starttime
-        domain.set_starttime(5.0)
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.get_time() + 7)/1000)
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.get_starttime() + 7)/1000)
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-    def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
-        """
-        Test that default rainfall can be used when given rate runs out of data.
-        """
-
-        import warnings
-        warnings.simplefilter('ignore', UserWarning)
-
-
-        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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # that expires at t==20
-        from anuga.fit_interpolate.interpolate import Modeltime_too_late
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+        factor = 10.0
+
 
         def main_rate(t):
@@ -949 +338 @@
                 raise Modeltime_too_late, msg
             else:
-                return 3*t + 7
-
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=main_rate,
-                     polygon = [[1,1], [2,1], [2,2], [1,2]],
-                     default_rate=5.0)
-
-        assert num.allclose(R.exchange_area, 2)
-        
-        domain.forcing_terms.append(R)
-
-        domain.time = 10.
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            (3*domain.time+7)/1000)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-        domain.time = 100.
-        domain.quantities['stage'].explicit_update[:] = 0.0     # Reset
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1],
-                            5.0/1000) # Default value
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-    def test_rainfall_forcing_with_evolve(self):
-        """test_rainfall_forcing_with_evolve
-
-        Test how forcing terms are called within evolve
-        """
-
-        # FIXME(Ole): This test is just to experiment
-        import warnings
-        warnings.simplefilter('ignore', UserWarning)
-
-
-        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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # that expires at t==20
-        from anuga.fit_interpolate.interpolate import Modeltime_too_late
-
-        def main_rate(t):
-            if t > 20:
-                msg = 'Model time exceeded.'
-                raise Modeltime_too_late, msg
-            else:
-                return 3*t + 7
-
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=main_rate,
-                     polygon=[[1,1], [2,1], [2,2], [1,2]],
-                     default_rate=5.0)
-
-        assert num.allclose(R.exchange_area, 2)
-        
-        domain.forcing_terms.append(R)
-
-        for t in domain.evolve(yieldstep=1, finaltime=25):
-            pass
-            #FIXME(Ole):  A test here is hard because explicit_update also
-            # receives updates from the flux calculation.
-
-
-    def test_rainfall_forcing_with_evolve_1(self):
-        """test_rainfall_forcing_with_evolve
-
-        Test how forcing terms are called within evolve.
-        This test checks that proper exception is thrown when no default_rate is set
-        """
-
-        import warnings
-        warnings.simplefilter('ignore', UserWarning)
-
-
-        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})
-
-        # Setup only one forcing term, time dependent rainfall
-        # that expires at t==20
-        from anuga.fit_interpolate.interpolate import Modeltime_too_late
-
-        def main_rate(t):
-            if t > 20:
-                msg = 'Model time exceeded.'
-                raise Modeltime_too_late, msg
-            else:
-                return 3*t + 7
-
-        domain.forcing_terms = []
-        R = Rainfall(domain,
-                     rate=main_rate,
-                     polygon=[[1,1], [2,1], [2,2], [1,2]])
-
-
-        assert num.allclose(R.exchange_area, 2)
-        
-        domain.forcing_terms.append(R)
-        #for t in domain.evolve(yieldstep=1, finaltime=25):
-        #    pass
-                
-        try:
-            for t in domain.evolve(yieldstep=1, finaltime=25):
-                pass
-        except Modeltime_too_late, e:
-            # Test that error message is as expected
-            assert 'can specify keyword argument default_rate in the forcing function' in str(e)
-        else:
-            raise Exception, 'Should have raised exception'
-
-    def test_constant_wind_stress_from_file(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        nrows=5; ncols = 6;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 6
-        timestep=12*60
-        eps=2e-16
-
-        points, vertices, boundary =rectangular(nrows-2,ncols-2,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        s = 100
-        phi = 135
-        pressure=1000
-        domain.forcing_terms = []
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=10,
-                                      speed=s,
-                                      angle=phi,
-                                      pressure=pressure)
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['wind_speed', 'wind_angle'],
-                          interpolation_points = midpoints)
-
-        W = Wind_stress(F,use_coordinates=False)
-        domain.forcing_terms.append(W)
-        domain.compute_forcing_terms()
-
-        const = eta_w*rho_a / rho_w
-
-        # 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 * num.sqrt(u**2 + v**2)
-
-        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
-
-    def test_variable_windfield_from_file(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        #nrows=25; ncols = 25;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 10
-        timestep=1
-        eps=2.e-16
-        spatial_thinning=1
-
-        points, vertices, boundary =rectangular(nrows-2,ncols-2,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        time=num.arange(0,10,1,num.float)
-        eval_time=time[7];
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-        vertexpoints = domain.get_nodes()
-
-        """
-        x=grid_1d(xllcorner,cellsize,ncols)
-        y=grid_1d(yllcorner,cellsize,nrows)
-        X,Y=num.meshgrid(x,y)
-        interpolation_points=num.empty((X.shape[0]*X.shape[1],2),num.float)
-        k=0
-        for i in range(X.shape[0]):
-            for j in range(X.shape[1]):
-                interpolation_points[k,0]=X[i,j]
-                interpolation_points[k,1]=Y[i,j]
-                k+=1
-
-        z=spatial_linear_varying_speed(eval_time,interpolation_points[:,0],
-                                       interpolation_points[:,1])
-
-        k=0
-        Z=num.empty((X.shape[0],X.shape[1]),num.float)
-        for i in range(X.shape[0]):
-            for j in range(X.shape[1]):
-                Z[i,j]=z[k]
-                k+=1
-
-        Q=num.empty((time.shape[0],points.shape[0]),num.float)
-        for i, t in enumerate(time):
-            Q[i,:]=spatial_linear_varying_speed(t,points[:,0],points[:,1])
-
-        from interpolate import Interpolation_function
-        I  = Interpolation_function(time,Q,
-                                    vertex_coordinates = points,
-                                    triangles = domain.triangles,
-                                    #interpolation_points = midpoints,
-                                    interpolation_points=interpolation_points,
-                                    verbose=False)
-
-        V=num.empty((X.shape[0],X.shape[1]),num.float)
-        for k in range(len(interpolation_points)):
-            assert num.allclose(I(eval_time,k),z[k])
-            V[k/X.shape[1],k%X.shape[1]]=I(eval_time,k)
-
-
-           import mpl_toolkits.mplot3d.axes3d as p3
-           fig=P.figure()
-           ax = p3.Axes3D(fig)
-           ax.plot_surface(X,Y,V)
-           ax.plot_surface(X,Y,Z)
-           P.show()
-
-
-        """
-
-        # Flat surface with 1m of water
-        domain.set_quantity('elevation', 0)
-        domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)
-
-        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
-
-        # Write wind stress file (ensure that domain.time is covered)
-
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=10,
-                                      speed=spatial_linear_varying_speed,
-                                      angle=spatial_linear_varying_angle,
-                                      pressure=spatial_linear_varying_pressure)
-
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
-                     verbose=False)
-
-        # Setup wind stress
-        FW = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['wind_speed', 'wind_angle'],
-                          interpolation_points = midpoints)
-
-        W = Wind_stress(FW,use_coordinates=False)
-
-        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
-
-        xc = domain.get_centroid_coordinates()
-        t = domain.time
-
-        x = xc[:,0]
-        y = xc[:,1]
-        s_vec = spatial_linear_varying_speed(t,x,y)
-        phi_vec = spatial_linear_varying_angle(t,x,y)
-
-        for k in range(N):
-            # Convert to radians
-            phi = phi_vec[k]*pi / 180
-            s = s_vec[k]
-
-            # Compute velocity vector (u, v)
-            u = s*cos(phi)
-            v = s*sin(phi)
-
-            # Compute wind stress
-            S = const * num.sqrt(u**2 + v**2)
-
-            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
-
-            assert num.allclose(domain.quantities['xmomentum'].\
-                                    explicit_update[k],S*u,eps)
-            assert num.allclose(domain.quantities['ymomentum'].\
-                                     explicit_update[k],S*v,eps)
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_variable_pressurefield_from_file(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        #nrows=25; ncols = 25;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 10
-        timestep=1
-        eps=2.e-16
-        spatial_thinning=1
-
-        points, vertices, boundary =rectangular(nrows-2,ncols-2,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        time=num.arange(0,10,1,num.float)
-        eval_time=time[7];
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-        vertexpoints = domain.get_nodes()
-
-        # Flat surface with 1m of water
-        domain.set_quantity('elevation', 0)
-        domain.set_quantity('stage', 1.0)
-        domain.set_quantity('friction', 0)
-
-        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
-
-        # Write wind stress file (ensure that domain.time is covered)
-
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=10,
-                                      speed=spatial_linear_varying_speed,
-                                      angle=spatial_linear_varying_angle,
-                                      pressure=spatial_linear_varying_pressure)
-
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
-                     verbose=False)
-
-        # Setup barometric pressure
-        FP = file_function(field_sts_filename+'.sww', domain,
-                           quantities=['barometric_pressure'],
-                           interpolation_points = vertexpoints)
-
-        P = Barometric_pressure(FP,use_coordinates=False)
-
-
-        domain.forcing_terms = []
-        domain.forcing_terms.append(P)
-
-        domain.compute_forcing_terms()
-
-        N = len(domain)    # number_of_triangles
-
-        xc = domain.get_centroid_coordinates()
-        t = domain.time
-
-        x = xc[:,0]
-        y = xc[:,1]
-        p_vec = spatial_linear_varying_pressure(t,x,y)
-
-        h=1 #depth
-        px=0.000025  #pressure gradient in x-direction
-        py=0.0000125 #pressure gradient in y-direction
-        for k in range(N):
-            # Convert to radians
-            p = p_vec[k]
-
-            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
-
-            assert num.allclose(domain.quantities['xmomentum'].\
-                                    explicit_update[k],h*px/rho_w)
-
-            assert num.allclose(domain.quantities['ymomentum'].\
-                                     explicit_update[k],h*py/rho_w)
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_constant_wind_stress_from_file_evolve(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        nrows=5; ncols = 6;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 27
-        timestep=1
-        eps=2e-16
-
-        points, vertices, boundary =rectangular(nrows-2,ncols-2,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        s = 100
-        phi = 135
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=number_of_timesteps,
-                                      speed=s,
-                                      angle=phi)
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['wind_speed', 'wind_angle'],
-                          interpolation_points = midpoints)
-
-        W = Wind_stress(F,use_coordinates=False)
-        domain.forcing_terms.append(W)
-
-        valuesUsingFunction=num.empty((3,number_of_timesteps+1,midpoints.shape[0]),
-                                      num.float)
-        i=0
-        for t in domain.evolve(yieldstep=1, finaltime=number_of_timesteps*timestep):
-            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
-            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
-            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
-            i+=1
-
-
-        domain_II = Domain(points, vertices, boundary)
-
-        # Flat surface with 1m of water
-        domain_II.set_quantity('elevation', 0)
-        domain_II.set_quantity('stage', 1.0)
-        domain_II.set_quantity('friction', 0)
-
-        Br = Reflective_boundary(domain_II)
-        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        s = 100
-        phi = 135
-        domain_II.forcing_terms = []
-        domain_II.forcing_terms.append(Wind_stress(s, phi))
-
-        i=0;
-        for t in domain_II.evolve(yieldstep=1, 
-                                  finaltime=number_of_timesteps*timestep):
-            assert num.allclose(valuesUsingFunction[0,i],domain_II.quantities['stage'].explicit_update), max(valuesUsingFunction[0,i]-domain_II.quantities['stage'].explicit_update)
-            assert  num.allclose(valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update)
-            assert num.allclose(valuesUsingFunction[2,i],domain_II.quantities['ymomentum'].explicit_update)
-            i+=1
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_temporally_varying_wind_stress_from_file_evolve(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        #nrows=20; ncols = 20;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 28
-        timestep=1.
-        eps=2e-16
-
-        #points, vertices, boundary =rectangular(10,10,
-        points, vertices, boundary =rectangular(5,5,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=number_of_timesteps,
-                                      speed=time_varying_speed,
-                                      angle=time_varying_angle,
-                                      pressure=time_varying_pressure)
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['wind_speed', 'wind_angle'],
-                          interpolation_points = midpoints)
-
-        #W = Wind_stress(F,use_coordinates=False)
-        W = Wind_stress_fast(F,filename=field_sts_filename+'.sww', domain=domain)
-        domain.forcing_terms.append(W)
-
-        valuesUsingFunction=num.empty((3,2*number_of_timesteps,midpoints.shape[0]),
-                                      num.float)
-        i=0
-        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
-            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
-            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
-            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
-            i+=1
-
-
-        domain_II = Domain(points, vertices, boundary)
-
-        # Flat surface with 1m of water
-        domain_II.set_quantity('elevation', 0)
-        domain_II.set_quantity('stage', 1.0)
-        domain_II.set_quantity('friction', 0)
-
-        Br = Reflective_boundary(domain_II)
-        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        domain_II.forcing_terms.append(Wind_stress(s=time_varying_speed, 
-                                                   phi=time_varying_angle))
-
-        i=0;
-        for t in domain_II.evolve(yieldstep=timestep/2., 
-                                  finaltime=(number_of_timesteps-1)*timestep):
-            assert num.allclose(valuesUsingFunction[0,i],
-                                domain_II.quantities['stage'].explicit_update,
-                                eps)
-            #print i,valuesUsingFunction[1,i]
-            assert  num.allclose(valuesUsingFunction[1,i],
-                                 domain_II.quantities['xmomentum'].explicit_update,
-                                 eps),(valuesUsingFunction[1,i]-
-                                 domain_II.quantities['xmomentum'].explicit_update)
-            assert num.allclose(valuesUsingFunction[2,i],
-                                domain_II.quantities['ymomentum'].explicit_update,
-                                eps)
-            #if i==1: assert-1==1
-            i+=1
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_spatially_varying_wind_stress_from_file_evolve(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        nrows=20; ncols = 20;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 28
-        timestep=1.
-        eps=2e-16
-
-        #points, vertices, boundary =rectangular(10,10,
-        points, vertices, boundary =rectangular(5,5,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        midpoints = domain.get_centroid_coordinates()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=number_of_timesteps,
-                                      speed=spatial_linear_varying_speed,
-                                      angle=spatial_linear_varying_angle,
-                                      pressure=spatial_linear_varying_pressure)
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['wind_speed', 'wind_angle'],
-                          interpolation_points = midpoints)
-
-        W = Wind_stress(F,use_coordinates=False)
-        domain.forcing_terms.append(W)
-
-        valuesUsingFunction=num.empty((3,number_of_timesteps,midpoints.shape[0]),
-                                      num.float)
-        i=0
-        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
-            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
-            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
-            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
-            i+=1
-
-
-        domain_II = Domain(points, vertices, boundary)
-
-        # Flat surface with 1m of water
-        domain_II.set_quantity('elevation', 0)
-        domain_II.set_quantity('stage', 1.0)
-        domain_II.set_quantity('friction', 0)
-
-        Br = Reflective_boundary(domain_II)
-        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        domain_II.forcing_terms.append(Wind_stress(s=spatial_linear_varying_speed, 
-                                                   phi=spatial_linear_varying_angle))
-
-        i=0;
-        for t in domain_II.evolve(yieldstep=timestep, 
-                                  finaltime=(number_of_timesteps-1)*timestep):
-            #print valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update
-            assert num.allclose(valuesUsingFunction[0,i],
-                                domain_II.quantities['stage'].explicit_update,
-                                eps)
-            assert  num.allclose(valuesUsingFunction[1,i],
-                                 domain_II.quantities['xmomentum'].explicit_update,
-                                 eps)
-            assert num.allclose(valuesUsingFunction[2,i],
-                                domain_II.quantities['ymomentum'].explicit_update,
-                                eps)
-            i+=1
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_temporally_varying_pressure_stress_from_file_evolve(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        #nrows=20; ncols = 20;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 28
-        timestep=10.
-        eps=2e-16
-
-        #print "Building mesh"
-        #points, vertices, boundary =rectangular(10,10,
-        points, vertices, boundary =rectangular(5,5,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        vertexpoints = domain.get_nodes()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        field_sts_filename = 'wind_field'
-        #print 'Writing pressure field sts file'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=number_of_timesteps,
-                                      speed=time_varying_speed,
-                                      angle=time_varying_angle,
-                                      pressure=time_varying_pressure)
-
-        #print "converting sts to sww"
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        #print 'initialising file_function'
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['barometric_pressure'],
-                          interpolation_points = vertexpoints)
-
-        #P = Barometric_pressure(F,use_coordinates=False)
-        #print 'initialising pressure forcing term'
-        P = Barometric_pressure_fast(p=F,filename=field_sts_filename+'.sww',domain=domain)
-        domain.forcing_terms.append(P)
-
-        valuesUsingFunction=num.empty((3,2*number_of_timesteps,len(domain)),
-                                      num.float)
-        i=0
-        import time as timer
-        t0=timer.time()
-        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
-            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
-            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
-            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
-            i+=1
-            #domain.write_time()
-        t1=timer.time()
-        #print "That took %fs seconds" %(t1-t0)
-
-
-        domain_II = Domain(points, vertices, boundary)
-
-        # Flat surface with 1m of water
-        domain_II.set_quantity('elevation', 0)
-        domain_II.set_quantity('stage', 1.0)
-        domain_II.set_quantity('friction', 0)
-
-        Br = Reflective_boundary(domain_II)
-        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        domain_II.forcing_terms.append(Barometric_pressure(p=time_varying_pressure))
-
-        i=0;
-        for t in domain_II.evolve(yieldstep=timestep/2., 
-                                  finaltime=(number_of_timesteps-1)*timestep):
-            assert num.allclose(valuesUsingFunction[0,i],
-                                domain_II.quantities['stage'].explicit_update,
-                                eps)
-            assert  num.allclose(valuesUsingFunction[1,i],
-                                 domain_II.quantities['xmomentum'].explicit_update,
-                                 eps)
-            assert num.allclose(valuesUsingFunction[2,i],
-                                domain_II.quantities['ymomentum'].explicit_update,
-                                eps)
-            i+=1
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_spatially_varying_pressure_stress_from_file_evolve(self):
-        from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin
-        from anuga.config import time_format
-        from anuga.file_conversion.sts2sww_mesh import sts2sww_mesh
-
-        cellsize = 25
-        #nrows=20; ncols = 20;
-        nrows=10; ncols = 10;
-        refzone=50
-        xllcorner=366000;yllcorner=6369500;
-        number_of_timesteps = 28
-        timestep=1.
-        eps=2e-16
-
-        #points, vertices, boundary =rectangular(10,10,
-        points, vertices, boundary =rectangular(5,5,
-                                                len1=cellsize*(ncols-1),
-                                                len2=cellsize*(nrows-1),
-                                                origin=(xllcorner,yllcorner))
-
-        domain = Domain(points, vertices, boundary)
-        vertexpoints = domain.get_nodes()
-
-        # 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({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        # Setup only one forcing term, constant wind stress
-        field_sts_filename = 'wind_field'
-        self.write_wind_pressure_field_sts(field_sts_filename,
-                                      nrows=nrows,
-                                      ncols=ncols,
-                                      cellsize=cellsize,
-                                      origin=(xllcorner,yllcorner),
-                                      refzone=50,
-                                      timestep=timestep,
-                                      number_of_timesteps=number_of_timesteps,
-                                      speed=spatial_linear_varying_speed,
-                                      angle=spatial_linear_varying_angle,
-                                      pressure=spatial_linear_varying_pressure)
-
-        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
-                     verbose=False)
-
-        # Setup wind stress
-        F = file_function(field_sts_filename+'.sww', domain,
-                          quantities=['barometric_pressure'],
-                          interpolation_points = vertexpoints)
-
-        P = Barometric_pressure(F,use_coordinates=False)
-        domain.forcing_terms.append(P)
-
-        valuesUsingFunction=num.empty((3,number_of_timesteps,len(domain)),
-                                      num.float)
-        i=0
-        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
-            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
-            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
-            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
-            i+=1
-
-
-        domain_II = Domain(points, vertices, boundary)
-
-        # Flat surface with 1m of water
-        domain_II.set_quantity('elevation', 0)
-        domain_II.set_quantity('stage', 1.0)
-        domain_II.set_quantity('friction', 0)
-
-        Br = Reflective_boundary(domain_II)
-        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
-
-        domain_II.forcing_terms.append(Barometric_pressure(p=spatial_linear_varying_pressure))
-
-        i=0;
-        for t in domain_II.evolve(yieldstep=timestep, 
-                                  finaltime=(number_of_timesteps-1)*timestep):
-
-            assert num.allclose(valuesUsingFunction[0,i],
-                                domain_II.quantities['stage'].explicit_update,
-                                eps)
-            assert  num.allclose(valuesUsingFunction[1,i],
-                                 domain_II.quantities['xmomentum'].explicit_update,
-                                 eps)
-            assert num.allclose(valuesUsingFunction[2,i],
-                                domain_II.quantities['ymomentum'].explicit_update,
-                                eps)
-            i+=1
-
-        os.remove(field_sts_filename+'.sts')
-        os.remove(field_sts_filename+'.sww')
-
-    def test_gravity(self):
-        #Assuming no friction
-
-        from anuga.config import g
-
-        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)
-        B = Reflective_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        #Set up for a gradient of (3,0) at mid triangle (bce)
-        def slope(x, y):
-            return 3*x
-
-        h = 0.1
-        def stage(x, y):
-            return slope(x, y) + h
-
-        domain.set_quantity('elevation', slope)
-        domain.set_quantity('stage', stage)
-
-        for name in domain.conserved_quantities:
-            assert num.allclose(domain.quantities[name].explicit_update, 0)
-            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
-
-        domain.update_boundary()
-        domain.compute_fluxes()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
-                            -g*h*3)
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
-
-
-
-    def test_manning_friction_old(self):
-        from anuga.config import g
-
-        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)
-
-        domain.set_sloped_mannings_function(False)
-        
-        B = Reflective_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        #Set up for a gradient of (3,0) at mid triangle (bce)
-        def slope(x, y):
-            return 3*x
-
-        h = 0.1
-        def stage(x, y):
-            return slope(x, y) + h
-
-        eta = 0.07
-        domain.set_quantity('elevation', slope)
-        domain.set_quantity('stage', stage)
-        domain.set_quantity('friction', eta)
-
-        for name in domain.conserved_quantities:
-            assert num.allclose(domain.quantities[name].explicit_update, 0)
-            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
-                            0)
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
-
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
-                            0)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
-                            0)
-
-        #Create some momentum for friction to work with
-        domain.set_quantity('xmomentum', 1)
-        S = -g*eta**2 / h**(7.0/3)
-
-        domain.compute_forcing_terms()
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
-                            S)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
-                            0)
-
-        #A more complex example
-        domain.quantities['stage'].semi_implicit_update[:] = 0.0
-        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
-        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
-
-        domain.set_quantity('xmomentum', 3)
-        domain.set_quantity('ymomentum', 4)
-        # sqrt(3^2 +4^2) = 5
-
-        S = -g*eta**2 / h**(7.0/3)  * 5
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,3*S)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,4*S)
-
-
-    def test_manning_friction_new(self):
-        from anuga.config import g
-        import math
-
-        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)
-        B = Reflective_boundary(domain)
-        domain.set_boundary( {'exterior': B})
-
-        # Use the new function which takes into account the extra
-        # wetted area due to slope of bed
-        domain.set_sloped_mannings_function(True)
-
-        #Set up for a gradient of (3,0) at mid triangle (bce)
-        def slope(x, y):
-            return 3*x
-
-        h = 0.1
-        def stage(x, y):
-            return slope(x, y) + h
-
-        eta = 0.07
-        domain.set_quantity('elevation', slope)
-        domain.set_quantity('stage', stage)
-        domain.set_quantity('friction', eta)
-
-        for name in domain.conserved_quantities:
-            assert num.allclose(domain.quantities[name].explicit_update, 0)
-            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
-
-        domain.compute_forcing_terms()
-
-        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
-                            0)
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
-
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
-                            0)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
-                            0)
-
-        #Create some momentum for friction to work with
-        domain.set_quantity('xmomentum', 1)
-        S = -g*eta**2 / h**(7.0/3) * math.sqrt(10)
-
-        domain.compute_forcing_terms()
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
-                            S)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
-                            0)
-
-        #A more complex example
-        domain.quantities['stage'].semi_implicit_update[:] = 0.0
-        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
-        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
-
-        domain.set_quantity('xmomentum', 3)
-        domain.set_quantity('ymomentum', 4)
-
-        S = -g*eta**2 *5 / h**(7.0/3) * math.sqrt(10.0)
-
-        domain.compute_forcing_terms()
-
-        #print 'S', S
-        #print domain.quantities['xmomentum'].semi_implicit_update
-        #print domain.quantities['ymomentum'].semi_implicit_update
-
-        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,3*S)
-        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,4*S)
-
-
-
-
-
-    def test_inflow_using_circle(self):
-        from math import pi, cos, sin
-
-        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})
-
-        # Setup only one forcing term, constant inflow of 2 m^3/s
-        # on a circle affecting triangles #0 and #1 (bac and bce)
-        domain.forcing_terms = []
-
-        I = Inflow(domain, rate=2.0, center=(1,1), radius=1)
-        domain.forcing_terms.append(I)
-        domain.compute_forcing_terms()
-
-
-        A = I.exchange_area
-        assert num.allclose(A, 4) # Two triangles
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-
-    def test_inflow_using_circle_function(self):
-        from math import pi, cos, sin
-
-        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})
-
-        # Setup only one forcing term, time dependent inflow of 2 m^3/s
-        # on a circle affecting triangles #0 and #1 (bac and bce)
-        domain.forcing_terms = []
-        I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
-        domain.forcing_terms.append(I)
-
-        domain.compute_forcing_terms()
-
-        A = I.exchange_area
-        assert num.allclose(A, 4) # Two triangles
-
-        assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/A)
-        assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/A)
-        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
-
-
-
-
-    def test_inflow_catch_too_few_triangles(self):
-        """
-        Test that exception is thrown if no triangles are covered
-        by the inflow area
-        """
-
-        from math import pi, cos, sin
-
-        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})
-
-        # Setup only one forcing term, constant inflow of 2 m^3/s
-        # on a circle affecting triangles #0 and #1 (bac and bce)
-        try:
-            Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01)
-        except:
-            pass
-        else:
-            msg = 'Should have raised exception'
-            raise Exception, msg
-
+                return 3.0 * t + 7.0
+
+        default_rate = lambda t: 3*t + 7
+
+
+        operator = Rate_operator(domain, rate=main_rate, factor=factor, \
+                      indices=indices, default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        t = operator.get_time()
+        d = operator.get_timestep()*main_rate(t)*factor + 1
+        stage_ex = [ d,  d,   1.,  d]
+
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+
+        domain.set_starttime(30.0)
+        domain.timestep = 1.0
+        operator()
+
+        t = operator.get_time()
+        d = operator.get_timestep()*default_rate(t)*factor + d
+        stage_ex = [ d,  d,   1.,  d]
+
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
 
             
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Forcing, 'test')
+    suite = unittest.makeSuite(Test_rate_operators, 'test')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

trunk/anuga_core/source/anuga/structures/inlet_operator.py

@@ -46 +46 @@
         t = self.domain.get_time()
 
+
+
+        # Need to run global command on all processors
+        current_volume = self.inlet.get_total_water_volume()
+
         Q1 = self.update_Q(t)
         Q2 = self.update_Q(t + timestep)
@@ -52 +57 @@
         volume = Q*timestep
         
-        assert volume >= 0.0, 'Q < 0: Water to be removed from an inlet!'
+        assert current_volume + volume >= 0.0, 'Requesting too much water to be removed from an inlet!'
 
         # store last discharge