Changeset 8116

Timestamp:

Jan 31, 2011, 4:33:40 PM (12 years ago)

Author:

jakeman

Message:

jakeman: updated forcing.py to include wind and pressure stress terms. Updated test_forcing.py with associated unit tests

Location:

trunk/anuga_core/source/anuga/shallow_water

Files:

• forcing.py (modified) (10 diffs)
• test_forcing.py (modified) (5 diffs)

trunk/anuga_core/source/anuga/shallow_water/forcing.py

@@ -21 +21 @@
 from warnings import warn
 import numpy as num
-
-
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from copy import copy
 
 def check_forcefield(f):
@@ -62 +63 @@
         except:
             msg = '%s must return vector' % func_msg
-            self.fail(msg)
+            raise Exception, msg
         msg = '%s must return vector of length %d' % (func_msg, N)
         assert result_len == N, msg
@@ -123 +124 @@
         from anuga.config import rho_a, rho_w, eta_w
 
+        self.use_coordinates=True
         if len(args) == 2:
             s = args[0]
@@ -129 +131 @@
             # Assume vector function returning (s, phi)(t,x,y)
             vector_function = args[0]
-            s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
-            phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+            if ( len(kwargs)==1 ):
+                self.use_coordinates=kwargs['use_coordinates']
+            else:
+                self.use_coordinates=True
+            if ( self.use_coordinates ):
+                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+            else:
+                s = lambda t,i: vector_function(t,point_id=i)[0]
+                phi = lambda t,i: vector_function(t,point_id=i)[1]
         else:
            # Assume info is in 2 keyword arguments
@@ -139 +149 @@
                raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
 
-        self.speed = check_forcefield(s)
-        self.phi = check_forcefield(phi)
+        if ( self.use_coordinates ):
+            self.speed = check_forcefield(s)
+            self.phi = check_forcefield(phi)
+        else:
+            self.speed = s
+            self.phi = phi
 
         self.const = eta_w*rho_a/rho_w
-
     ##
     # @brief 'execute' this class instance.
@@ -150 +163 @@
         """Evaluate windfield based on values found in domain"""
 
-        from math import pi, cos, sin, sqrt
-
         xmom_update = domain.quantities['xmomentum'].explicit_update
         ymom_update = domain.quantities['ymomentum'].explicit_update
@@ -160 +171 @@
         if callable(self.speed):
             xc = domain.get_centroid_coordinates()
-            s_vec = self.speed(t, xc[:,0], xc[:,1])
+            if ( self.use_coordinates ):
+                s_vec = self.speed(t, xc[:,0], xc[:,1])
+            else:
+                s_vec=num.empty(N,float)
+                for i in range(N):
+                    s_vec[i]=self.speed(t,i)
         else:
             # Assume s is a scalar
@@ -171 +187 @@
         if callable(self.phi):
             xc = domain.get_centroid_coordinates()
-            phi_vec = self.phi(t, xc[:,0], xc[:,1])
+            if ( self.use_coordinates ):
+                phi_vec = self.phi(t, xc[:,0], xc[:,1])
+            else:
+                phi_vec=num.empty(len(xc),float)
+                for i in range(len(xc)):
+                    phi_vec[i]=self.phi(t,i)
         else:
             # Assume phi is a scalar
@@ -195 +216 @@
                             s_vec, phi_vec, const):
     """Python version of assigning wind field to update vectors.
-    A C version also exists (for speed)
     """
 
@@ -858 +878 @@
         return average_energy
 
+class Barometric_pressure:
+    """ Apply barometric pressure stress to water momentum in terms of
+        barometric pressure p [hPa]. If the pressure data is stored in a file
+        file_function is used to create a callable function. The data file 
+        contains pressure values at a set of possibly arbitrarily located nodes
+        at a set o possibly irregular but increasing times. file_function 
+        interpolates from the file data onto the vertices of the domain.mesh
+        for each time. The file_function is called at every timestep during 
+        the evolve function call.
+    """
+    def __init__(self, *args, **kwargs):
+        """Initialise barometric pressure field from barometric pressure [hPa]
+        Input p can be either scalars or Python functions, e.g.
+
+        P = barometric_pressure(1000)
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def pressure(t,x,y):
+            ...
+            return p
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        P = Barometric_pressure(pressure)
+
+        agruments can also be the ANGUA file_function, e.g.
+        F = file_function(sww_filename,domain,quantities,interpolation_points)
+        The interpolation_points must be the mesh vertices returned by 
+        domain.get_nodes(). Quantities = ['barometric_pressure']
+
+        The file_function is passed using
+
+        P = Barometric_pressure(F, use_coordinates=True/False)
+
+        The instantiated object P can be appended to the list of
+        forcing_terms as in
+
+        domain.forcing_terms.append(P)
+        """
+
+        from anuga.config import rho_a, rho_w, eta_w
+
+        self.use_coordinates=True
+        if len(args) == 1:
+            if ( not callable(args[0]) ):
+                pressure=args[0]
+            else:
+                # Assume vector function returning (pressure)(t,x,y)
+                vector_function = args[0]
+                if ( len(kwargs)==1 ):
+                    self.use_coordinates=kwargs['use_coordinates']
+                else:
+                    self.use_coordinates=True
+
+                if ( self.use_coordinates ):
+                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+                else:
+                    p = lambda t,i: vector_function(t,point_id=i)[0]
+        else:
+           # Assume info is in 1 or 2 keyword arguments
+           if ( len(kwargs) == 1 ):
+               p = kwargs['p']
+           elif ( len(kwargs)==2 ):
+               p = kwargs['p']
+               self.use_coordinates = kwargs['use_coordinates']
+           else:
+               raise Exception, 'Assumes one keyword argument: p=... or two keyword arguments p=...,use_coordinates=...'
+
+        if ( self.use_coordinates ):
+            self.pressure = check_forcefield(p)
+        else:
+            self.pressure = p
+
+    ##
+    # @brief 'execute' this class instance.
+    # @param domain 
+    def __call__(self, domain):
+        """Evaluate pressure field based on values found in domain"""
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        N = domain.get_number_of_nodes()
+        t = domain.time
+
+        if callable(self.pressure):
+            xv = domain.get_nodes()
+            if ( self.use_coordinates ):
+                p_vec = self.pressure(t, xv[:,0], xv[:,1])
+            else:
+                p_vec=num.empty(N,num.float)
+                for i in range(N):
+                    p_vec[i]=self.pressure(t,i)
+        else:
+            # Assume s is a scalar
+            try:
+                p_vec = self.pressure * num.ones(N, num.float)
+            except:
+                msg = 'Pressure must be either callable or a scalar: %s' %self.s
+                raise msg
+
+        stage = domain.quantities['stage']
+        elevation = domain.quantities['elevation']
+
+        #FIXME SR Should avoid allocating memory!
+        height = stage.centroid_values - elevation.centroid_values
+
+        point = domain.get_vertex_coordinates()
+
+        assign_pressure_field_values(height, p_vec, point, domain.triangles,
+                                     xmom_update, ymom_update)
+
+
+##
+# @brief Assign pressure field values
+# @param xmom_update 
+# @param ymom_update
+# @param s_vec 
+# @param phi_vec 
+# @param const 
+def assign_pressure_field_values(height, pressure, x, triangles, 
+                                 xmom_update, ymom_update):
+    """Python version of assigning pressure field to update vectors.
+    """
+
+    from utilities.numerical_tools import gradient
+    from anuga.config import rho_a, rho_w, eta_w
+
+    N = len(height)
+    for k in range(N):
+
+        # Compute pressure slope
+
+        p0 = pressure[triangles[k][0]]
+        p1 = pressure[triangles[k][1]]
+        p2 = pressure[triangles[k][2]]
+
+        k3=3*k
+        x0 = x[k3 + 0][0]
+        y0 = x[k3 + 0][1]
+        x1 = x[k3 + 1][0]
+        y1 = x[k3 + 1][1]
+        x2 = x[k3 + 2][0]
+        y2 = x[k3 + 2][1]
+
+        px,py=gradient(x0, y0, x1, y1, x2, y2, p0, p1, p2)
+
+        xmom_update[k] += height[k]*px/rho_w
+        ymom_update[k] += height[k]*py/rho_w
+
+
+class Barometric_pressure_fast:
+    """ Apply barometric pressure stress to water momentum in terms of
+        barometric pressure p [hPa]. If the pressure data is stored in a file
+        file_function is used to create a callable function. The data file 
+        contains pressure values at a set of possibly arbitrarily located nodes
+        at a set o possibly irregular but increasing times. file_function 
+        interpolates from the file data onto the vertices of the domain.mesh
+        for each time. Two arrays are then stored p0=p(t0,:) and p1=p(t1,:) 
+        where t0<=domain.time<=t1. These arrays are recalculated when necessary
+        i.e t>t1. A linear temporal interpolation is used to approximate 
+        pressure at time t.
+    """
+    def __init__(self, *args, **kwargs):
+        """Initialise barometric pressure field from barometric pressure [hPa]
+        Input p can be either scalars or Python functions, e.g.
+
+        P = barometric_pressure(1000)
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def pressure(t,x,y):
+            ...
+            return p
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        P = Barometric_pressure(pressure)
+
+        Agruments can also be the ANGUA file_function, e.g.
+        F = file_function(sww_filename,domain,quantities,interpolation_points)
+        The interpolation_points must be the mesh vertices returned by 
+        domain.get_nodes(). Quantities = ['barometric_pressure']
+
+        The file_function is passed using
+
+        P = Barometric_pressure(F, filename=swwname, domain=domain)
+
+        The instantiated object P can be appended to the list of
+        forcing_terms as in
+
+        domain.forcing_terms.append(P)
+        """
+
+        from anuga.config import rho_a, rho_w, eta_w
+
+        self.use_coordinates=True
+        if len(args) == 1:
+            if ( not callable(args[0]) ):
+                pressure=args[0]
+            else:
+                # Assume vector function returning (pressure)(t,x,y)
+                vector_function = args[0]
+                if ( len(kwargs)==0 ):
+                    self.usre_coordinates=True
+                elif (len(kwargs)==2):
+                    filename=kwargs['filename']
+                    domain=kwargs['domain']
+                    self.use_coordinates=False
+                else:
+                    raise Exception, 'Assumes zero or two keyword arguments filename=...,domain=...'
+
+                if ( self.use_coordinates ):
+                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+                else:
+                    p = lambda t,i: vector_function(t,point_id=i)[0]
+        else:
+           # Assume info is in 1 or 2 keyword arguments
+           if ( len(kwargs) == 1 ):
+               p = kwargs['p']
+               self.use_coordinates=True
+           elif ( len(kwargs)==3 ):
+               p = kwargs['p']
+               filename = kwargs['filename']
+               domain = kwargs['domain']
+               self.use_coordinates = False
+           else:
+               raise Exception, 'Assumes one keyword argument: p=f(t,x,y,) or three keyword arguments p=f(t,i),filename=...,domain=...'
+
+        if ( self.use_coordinates ):
+            self.pressure = check_forcefield(p)
+        else:
+            self.pressure = p
+
+        if ( callable(self.pressure) and not self.use_coordinates):
+
+            # Open NetCDF file
+            fid = NetCDFFile(filename, netcdf_mode_r)
+            self.file_time = fid.variables['time'][:]
+            fid.close()
+
+            msg = 'pressure_file.starttime > domain.starttime'
+            if (self.file_time[0]>domain.starttime):
+                raise Exception, msg
+
+            msg = 'pressure_file[-1] < domain.starttime'
+            if (self.file_time[-1]<domain.starttime):
+                raise Exception, msg
+
+            msg = 'No pressure values exist for times greater than domain.starttime'
+            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
+                raise Exception, msg
+
+            # FIXME(JJ): How do we check that evolve 
+            # finaltime  < pressure_file.finaltime      
+            
+
+            self.index=0;
+            for i in range(len(self.file_time)):
+                if (self.file_time[i]<domain.starttime):
+                    self.index=i
+                else:
+                    break
+
+            N = domain.get_number_of_nodes()
+            self.prev_pressure_vertex_values=num.empty(N,num.float)
+            self.next_pressure_vertex_values=num.empty(N,num.float)
+            for i in range(N):
+                self.prev_pressure_vertex_values[i]=self.pressure(self.file_time[self.index],i)
+                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i)
+
+        self.p_vec=num.empty(N,num.float)
+
+
+    ##
+    # @brief 'execute' this class instance.
+    # @param domain 
+    def __call__(self, domain):
+        """Evaluate pressure field based on values found in domain"""
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        t = domain.time
+
+        if callable(self.pressure):
+            if ( self.use_coordinates ):
+                xv = domain.get_nodes()
+                self.p_vec = self.pressure(t, xv[:,0], xv[:,1])
+            else:
+                self.update_stored_pressure_values(domain)
+
+                # Linear temporal interpolation of pressure values
+                ratio = (t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index])
+                self.p_vec = self.prev_pressure_vertex_values + ratio*(self.next_pressure_vertex_values - self.prev_pressure_vertex_values)
+
+        else:
+            # Assume s is a scalar
+            try:
+                self.p_vec[:] = self.pressure
+            except:
+                msg = 'Pressure must be either callable function or a scalar: %s' %self.s
+                raise msg
+
+        stage = domain.quantities['stage']
+        elevation = domain.quantities['elevation']
+
+        height = stage.centroid_values - elevation.centroid_values
+
+        point = domain.get_vertex_coordinates()
+
+        assign_pressure_field_values(height, self.p_vec, point, 
+                                     domain.triangles,
+                                     xmom_update, ymom_update)
+
+    def update_stored_pressure_values(self,domain):
+        while (self.file_time[self.index+1]<domain.time):
+            self.index+=1
+            self.prev_pressure_vertex_values=copy(self.next_pressure_vertex_values)
+            for i in range(self.prev_pressure_vertex_values.shape[0]):
+                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i) 
+
+
+class Wind_stress_fast:
+    """ Apply wind stress to water momentum in terms of
+        wind speed [m/s] and wind direction [degrees]. 
+        If the wind data is stored in a file
+        file_function is used to create a callable function. The data file 
+        contains wind speed and direction values at a set of possibly 
+        arbitrarily located nodes
+        at a set of possibly irregular but increasing times. file_function 
+        interpolates from the file data onto the centroids of the domain.mesh
+        for each time. Two arrays for each wind quantity are then stored \
+        q0=q(t0,:) and q1=q(t1,:) 
+        where t0<=domain.time<=t1. These arrays are recalculated when necessary
+        i.e t>t1. A linear temporal interpolation is used to approximate 
+        pressure at time t.
+    """
+    def __init__(self, *args, **kwargs):
+        """Initialise windfield from wind speed s [m/s]
+        and wind direction phi [degrees]
+
+        Inputs v and phi can be either scalars or Python functions, e.g.
+
+        W = Wind_stress(10, 178)
+
+        #FIXME - 'normal' degrees are assumed for now, i.e. the
+        vector (1,0) has zero degrees.
+        We may need to convert from 'compass' degrees later on and also
+        map from True north to grid north.
+
+        Arguments can also be Python functions of t,x,y as in
+
+        def speed(t,x,y):
+            ...
+            return s
+
+        def angle(t,x,y):
+            ...
+            return phi
+
+        where x and y are vectors.
+
+        and then pass the functions in
+
+        W = Wind_stress(speed, angle)
+
+        The instantiated object W can be appended to the list of
+        forcing_terms as in
+
+        Alternatively, one vector valued function for (speed, angle)
+        can be applied, providing both quantities simultaneously.
+        As in
+        W = Wind_stress(F), where returns (speed, angle) for each t.
+
+        domain.forcing_terms.append(W)
+        """
+
+        from anuga.config import rho_a, rho_w, eta_w
+
+        self.use_coordinates=True
+        if len(args) == 2:
+            s = args[0]
+            phi = args[1]
+        elif len(args) == 1:
+            # Assume vector function returning (s, phi)(t,x,y)
+            vector_function = args[0]
+            if ( len(kwargs)==2 ):
+                filename=kwargs['filename']
+                domain=kwargs['domain']
+                self.use_coordinates=False
+            else:
+                self.use_coordinates=True
+            if ( self.use_coordinates ):
+                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
+                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
+            else:
+                s = lambda t,i: vector_function(t,point_id=i)[0]
+                phi = lambda t,i: vector_function(t,point_id=i)[1]
+        else:
+           # Assume info is in 2 keyword arguments
+           if len(kwargs) == 2:
+               s = kwargs['s']
+               phi = kwargs['phi']
+           else:
+               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
+
+        if ( self.use_coordinates ):
+            self.speed = check_forcefield(s)
+            self.phi = check_forcefield(phi)
+        else:
+            self.speed = s
+            self.phi = phi
+
+        N = len(domain)
+        if ( not self.use_coordinates):
+
+            # Open NetCDF file
+            fid = NetCDFFile(filename, netcdf_mode_r)
+            self.file_time = fid.variables['time'][:]
+            fid.close()
+
+            msg = 'wind_file.starttime > domain.starttime'
+            if (self.file_time[0]>domain.starttime):
+                raise Exception, msg
+
+            msg = 'wind_file[-1] < domain.starttime'
+            if (self.file_time[-1]<domain.starttime):
+                raise Exception, msg
+
+            msg = 'No wind values exist for times greater than domain.starttime'
+            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
+                raise Exception, msg
+
+            # FIXME(JJ): How do we check that evolve 
+            # finaltime  < wind_file.finaltime      
+            
+
+            self.index=0;
+            for i in range(len(self.file_time)):
+                if (self.file_time[i]<domain.starttime):
+                    self.index=i
+                else:
+                    break
+
+            self.prev_windspeed_centroid_values=num.empty(N,num.float)
+            self.next_windspeed_centroid_values=num.empty(N,num.float)
+            self.prev_windangle_centroid_values=num.empty(N,num.float)
+            self.next_windangle_centroid_values=num.empty(N,num.float)
+            for i in range(N):
+                self.prev_windspeed_centroid_values[i]=self.speed(self.file_time[self.index],i)
+                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i)
+                self.prev_windangle_centroid_values[i]=self.phi(self.file_time[self.index],i)
+                self.next_windangle_centroid_values[i]=self.phi(self.file_time[self.index+1],i)
+
+        self.s_vec=num.empty(N,num.float)
+        self.phi_vec=num.empty(N,num.float)
+
+        self.const = eta_w*rho_a/rho_w
+    ##
+    # @brief 'execute' this class instance.
+    # @param domain 
+    def __call__(self, domain):
+        """Evaluate windfield based on values found in domain"""
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        N = len(domain)    # number_of_triangles
+        t = domain.time
+
+        if callable(self.speed):
+            if ( self.use_coordinates ):
+                xc = domain.get_centroid_coordinates()
+                self.s_vec = self.speed(t, xc[:,0], xc[:,1])
+            else:
+                self.update_stored_wind_values(domain)
+
+                # Linear temporal interpolation of wind values
+                if t==self.file_time[self.index]:
+                    ratio = 0.
+                else:
+                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
+                self.s_vec = self.prev_windspeed_centroid_values + ratio*(self.next_windspeed_centroid_values - self.prev_windspeed_centroid_values)
+        else:
+            # Assume s is a scalar
+            try:
+                self.s_vec[:] = self.speed
+            except:
+                msg = 'Speed must be either callable or a scalar: %s' %self.s
+                raise msg
+
+        if callable(self.phi):
+            if ( self.use_coordinates ):
+                xc = domain.get_centroid_coordinates()
+                self.phi_vec = self.phi(t, xc[:,0], xc[:,1])
+            else:
+                self.update_stored_wind_values(domain)
+
+                # Linear temporal interpolation of wind values
+                if t==self.file_time[self.index]:
+                    ratio = 0.
+                else:
+                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
+                self.phi_vec = self.prev_windangle_centroid_values + ratio*(self.next_windangle_centroid_values - self.prev_windangle_centroid_values)
+        else:
+            # Assume phi is a scalar
+
+            try:
+                self.phi_vec[:] = self.phi
+            except:
+                msg = 'Angle must be either callable or a scalar: %s' %self.phi
+                raise msg
+
+        assign_windfield_values(xmom_update, ymom_update,
+                                self.s_vec, self.phi_vec, self.const)
+
+    def update_stored_wind_values(self,domain):
+        while (self.file_time[self.index+1]<domain.time):
+            self.index+=1
+            self.prev_windspeed_centroid_values=copy(self.next_windspeed_centroid_values)
+            self.prev_windangle_centroid_values=copy(self.next_windangle_centroid_values)
+            for i in range(self.next_windspeed_centroid_values.shape[0]):
+                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i) 
+                self.next_windangle_centroid_values[i]=self.phi(self.file_time[self.index+1],i) 

trunk/anuga_core/source/anuga/shallow_water/test_forcing.py

@@ -3 +3 @@
 
 import unittest, os
+import anuga
 from anuga.shallow_water.shallow_water_domain import Domain
 from boundaries import Reflective_boundary
@@ -8 +9 @@
 from anuga.file_conversion.file_conversion import timefile2netcdf
 from forcing import *
+from mesh_factory import rectangular
+from file_conversion.sts2sww_mesh import sts2sww_mesh
+from anuga.abstract_2d_finite_volumes.util import file_function
 
 import numpy as num
@@ -76 +80 @@
     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):
@@ -85 +217 @@
         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):
         from anuga.config import rho_a, rho_w, eta_w
@@ -855 +1101 @@
             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
+
+        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
+
+        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
+
+        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
+
+        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
+
+        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
+
+        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
+
+        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
+
+        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')
+
+
             
 if __name__ == "__main__":