Changeset 8116


Ignore:
Timestamp:
Jan 31, 2011, 4:33:40 PM (13 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:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/shallow_water/forcing.py

    r8073 r8116  
    2121from warnings import warn
    2222import numpy as num
    23 
    24 
     23from Scientific.IO.NetCDF import NetCDFFile
     24from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
     25from copy import copy
    2526
    2627def check_forcefield(f):
     
    6263        except:
    6364            msg = '%s must return vector' % func_msg
    64             self.fail(msg)
     65            raise Exception, msg
    6566        msg = '%s must return vector of length %d' % (func_msg, N)
    6667        assert result_len == N, msg
     
    123124        from anuga.config import rho_a, rho_w, eta_w
    124125
     126        self.use_coordinates=True
    125127        if len(args) == 2:
    126128            s = args[0]
     
    129131            # Assume vector function returning (s, phi)(t,x,y)
    130132            vector_function = args[0]
    131             s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
    132             phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
     133            if ( len(kwargs)==1 ):
     134                self.use_coordinates=kwargs['use_coordinates']
     135            else:
     136                self.use_coordinates=True
     137            if ( self.use_coordinates ):
     138                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
     139                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
     140            else:
     141                s = lambda t,i: vector_function(t,point_id=i)[0]
     142                phi = lambda t,i: vector_function(t,point_id=i)[1]
    133143        else:
    134144           # Assume info is in 2 keyword arguments
     
    139149               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
    140150
    141         self.speed = check_forcefield(s)
    142         self.phi = check_forcefield(phi)
     151        if ( self.use_coordinates ):
     152            self.speed = check_forcefield(s)
     153            self.phi = check_forcefield(phi)
     154        else:
     155            self.speed = s
     156            self.phi = phi
    143157
    144158        self.const = eta_w*rho_a/rho_w
    145 
    146159    ##
    147160    # @brief 'execute' this class instance.
     
    150163        """Evaluate windfield based on values found in domain"""
    151164
    152         from math import pi, cos, sin, sqrt
    153 
    154165        xmom_update = domain.quantities['xmomentum'].explicit_update
    155166        ymom_update = domain.quantities['ymomentum'].explicit_update
     
    160171        if callable(self.speed):
    161172            xc = domain.get_centroid_coordinates()
    162             s_vec = self.speed(t, xc[:,0], xc[:,1])
     173            if ( self.use_coordinates ):
     174                s_vec = self.speed(t, xc[:,0], xc[:,1])
     175            else:
     176                s_vec=num.empty(N,float)
     177                for i in range(N):
     178                    s_vec[i]=self.speed(t,i)
    163179        else:
    164180            # Assume s is a scalar
     
    171187        if callable(self.phi):
    172188            xc = domain.get_centroid_coordinates()
    173             phi_vec = self.phi(t, xc[:,0], xc[:,1])
     189            if ( self.use_coordinates ):
     190                phi_vec = self.phi(t, xc[:,0], xc[:,1])
     191            else:
     192                phi_vec=num.empty(len(xc),float)
     193                for i in range(len(xc)):
     194                    phi_vec[i]=self.phi(t,i)
    174195        else:
    175196            # Assume phi is a scalar
     
    195216                            s_vec, phi_vec, const):
    196217    """Python version of assigning wind field to update vectors.
    197     A C version also exists (for speed)
    198218    """
    199219
     
    858878        return average_energy
    859879
     880class Barometric_pressure:
     881    """ Apply barometric pressure stress to water momentum in terms of
     882        barometric pressure p [hPa]. If the pressure data is stored in a file
     883        file_function is used to create a callable function. The data file
     884        contains pressure values at a set of possibly arbitrarily located nodes
     885        at a set o possibly irregular but increasing times. file_function
     886        interpolates from the file data onto the vertices of the domain.mesh
     887        for each time. The file_function is called at every timestep during
     888        the evolve function call.
     889    """
     890    def __init__(self, *args, **kwargs):
     891        """Initialise barometric pressure field from barometric pressure [hPa]
     892        Input p can be either scalars or Python functions, e.g.
     893
     894        P = barometric_pressure(1000)
     895
     896        Arguments can also be Python functions of t,x,y as in
     897
     898        def pressure(t,x,y):
     899            ...
     900            return p
     901
     902        where x and y are vectors.
     903
     904        and then pass the functions in
     905
     906        P = Barometric_pressure(pressure)
     907
     908        agruments can also be the ANGUA file_function, e.g.
     909        F = file_function(sww_filename,domain,quantities,interpolation_points)
     910        The interpolation_points must be the mesh vertices returned by
     911        domain.get_nodes(). Quantities = ['barometric_pressure']
     912
     913        The file_function is passed using
     914
     915        P = Barometric_pressure(F, use_coordinates=True/False)
     916
     917        The instantiated object P can be appended to the list of
     918        forcing_terms as in
     919
     920        domain.forcing_terms.append(P)
     921        """
     922
     923        from anuga.config import rho_a, rho_w, eta_w
     924
     925        self.use_coordinates=True
     926        if len(args) == 1:
     927            if ( not callable(args[0]) ):
     928                pressure=args[0]
     929            else:
     930                # Assume vector function returning (pressure)(t,x,y)
     931                vector_function = args[0]
     932                if ( len(kwargs)==1 ):
     933                    self.use_coordinates=kwargs['use_coordinates']
     934                else:
     935                    self.use_coordinates=True
     936
     937                if ( self.use_coordinates ):
     938                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
     939                else:
     940                    p = lambda t,i: vector_function(t,point_id=i)[0]
     941        else:
     942           # Assume info is in 1 or 2 keyword arguments
     943           if ( len(kwargs) == 1 ):
     944               p = kwargs['p']
     945           elif ( len(kwargs)==2 ):
     946               p = kwargs['p']
     947               self.use_coordinates = kwargs['use_coordinates']
     948           else:
     949               raise Exception, 'Assumes one keyword argument: p=... or two keyword arguments p=...,use_coordinates=...'
     950
     951        if ( self.use_coordinates ):
     952            self.pressure = check_forcefield(p)
     953        else:
     954            self.pressure = p
     955
     956    ##
     957    # @brief 'execute' this class instance.
     958    # @param domain
     959    def __call__(self, domain):
     960        """Evaluate pressure field based on values found in domain"""
     961
     962        xmom_update = domain.quantities['xmomentum'].explicit_update
     963        ymom_update = domain.quantities['ymomentum'].explicit_update
     964
     965        N = domain.get_number_of_nodes()
     966        t = domain.time
     967
     968        if callable(self.pressure):
     969            xv = domain.get_nodes()
     970            if ( self.use_coordinates ):
     971                p_vec = self.pressure(t, xv[:,0], xv[:,1])
     972            else:
     973                p_vec=num.empty(N,num.float)
     974                for i in range(N):
     975                    p_vec[i]=self.pressure(t,i)
     976        else:
     977            # Assume s is a scalar
     978            try:
     979                p_vec = self.pressure * num.ones(N, num.float)
     980            except:
     981                msg = 'Pressure must be either callable or a scalar: %s' %self.s
     982                raise msg
     983
     984        stage = domain.quantities['stage']
     985        elevation = domain.quantities['elevation']
     986
     987        #FIXME SR Should avoid allocating memory!
     988        height = stage.centroid_values - elevation.centroid_values
     989
     990        point = domain.get_vertex_coordinates()
     991
     992        assign_pressure_field_values(height, p_vec, point, domain.triangles,
     993                                     xmom_update, ymom_update)
     994
     995
     996##
     997# @brief Assign pressure field values
     998# @param xmom_update
     999# @param ymom_update
     1000# @param s_vec
     1001# @param phi_vec
     1002# @param const
     1003def assign_pressure_field_values(height, pressure, x, triangles,
     1004                                 xmom_update, ymom_update):
     1005    """Python version of assigning pressure field to update vectors.
     1006    """
     1007
     1008    from utilities.numerical_tools import gradient
     1009    from anuga.config import rho_a, rho_w, eta_w
     1010
     1011    N = len(height)
     1012    for k in range(N):
     1013
     1014        # Compute pressure slope
     1015
     1016        p0 = pressure[triangles[k][0]]
     1017        p1 = pressure[triangles[k][1]]
     1018        p2 = pressure[triangles[k][2]]
     1019
     1020        k3=3*k
     1021        x0 = x[k3 + 0][0]
     1022        y0 = x[k3 + 0][1]
     1023        x1 = x[k3 + 1][0]
     1024        y1 = x[k3 + 1][1]
     1025        x2 = x[k3 + 2][0]
     1026        y2 = x[k3 + 2][1]
     1027
     1028        px,py=gradient(x0, y0, x1, y1, x2, y2, p0, p1, p2)
     1029
     1030        xmom_update[k] += height[k]*px/rho_w
     1031        ymom_update[k] += height[k]*py/rho_w
     1032
     1033
     1034class Barometric_pressure_fast:
     1035    """ Apply barometric pressure stress to water momentum in terms of
     1036        barometric pressure p [hPa]. If the pressure data is stored in a file
     1037        file_function is used to create a callable function. The data file
     1038        contains pressure values at a set of possibly arbitrarily located nodes
     1039        at a set o possibly irregular but increasing times. file_function
     1040        interpolates from the file data onto the vertices of the domain.mesh
     1041        for each time. Two arrays are then stored p0=p(t0,:) and p1=p(t1,:)
     1042        where t0<=domain.time<=t1. These arrays are recalculated when necessary
     1043        i.e t>t1. A linear temporal interpolation is used to approximate
     1044        pressure at time t.
     1045    """
     1046    def __init__(self, *args, **kwargs):
     1047        """Initialise barometric pressure field from barometric pressure [hPa]
     1048        Input p can be either scalars or Python functions, e.g.
     1049
     1050        P = barometric_pressure(1000)
     1051
     1052        Arguments can also be Python functions of t,x,y as in
     1053
     1054        def pressure(t,x,y):
     1055            ...
     1056            return p
     1057
     1058        where x and y are vectors.
     1059
     1060        and then pass the functions in
     1061
     1062        P = Barometric_pressure(pressure)
     1063
     1064        Agruments can also be the ANGUA file_function, e.g.
     1065        F = file_function(sww_filename,domain,quantities,interpolation_points)
     1066        The interpolation_points must be the mesh vertices returned by
     1067        domain.get_nodes(). Quantities = ['barometric_pressure']
     1068
     1069        The file_function is passed using
     1070
     1071        P = Barometric_pressure(F, filename=swwname, domain=domain)
     1072
     1073        The instantiated object P can be appended to the list of
     1074        forcing_terms as in
     1075
     1076        domain.forcing_terms.append(P)
     1077        """
     1078
     1079        from anuga.config import rho_a, rho_w, eta_w
     1080
     1081        self.use_coordinates=True
     1082        if len(args) == 1:
     1083            if ( not callable(args[0]) ):
     1084                pressure=args[0]
     1085            else:
     1086                # Assume vector function returning (pressure)(t,x,y)
     1087                vector_function = args[0]
     1088                if ( len(kwargs)==0 ):
     1089                    self.usre_coordinates=True
     1090                elif (len(kwargs)==2):
     1091                    filename=kwargs['filename']
     1092                    domain=kwargs['domain']
     1093                    self.use_coordinates=False
     1094                else:
     1095                    raise Exception, 'Assumes zero or two keyword arguments filename=...,domain=...'
     1096
     1097                if ( self.use_coordinates ):
     1098                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
     1099                else:
     1100                    p = lambda t,i: vector_function(t,point_id=i)[0]
     1101        else:
     1102           # Assume info is in 1 or 2 keyword arguments
     1103           if ( len(kwargs) == 1 ):
     1104               p = kwargs['p']
     1105               self.use_coordinates=True
     1106           elif ( len(kwargs)==3 ):
     1107               p = kwargs['p']
     1108               filename = kwargs['filename']
     1109               domain = kwargs['domain']
     1110               self.use_coordinates = False
     1111           else:
     1112               raise Exception, 'Assumes one keyword argument: p=f(t,x,y,) or three keyword arguments p=f(t,i),filename=...,domain=...'
     1113
     1114        if ( self.use_coordinates ):
     1115            self.pressure = check_forcefield(p)
     1116        else:
     1117            self.pressure = p
     1118
     1119        if ( callable(self.pressure) and not self.use_coordinates):
     1120
     1121            # Open NetCDF file
     1122            fid = NetCDFFile(filename, netcdf_mode_r)
     1123            self.file_time = fid.variables['time'][:]
     1124            fid.close()
     1125
     1126            msg = 'pressure_file.starttime > domain.starttime'
     1127            if (self.file_time[0]>domain.starttime):
     1128                raise Exception, msg
     1129
     1130            msg = 'pressure_file[-1] < domain.starttime'
     1131            if (self.file_time[-1]<domain.starttime):
     1132                raise Exception, msg
     1133
     1134            msg = 'No pressure values exist for times greater than domain.starttime'
     1135            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
     1136                raise Exception, msg
     1137
     1138            # FIXME(JJ): How do we check that evolve
     1139            # finaltime  < pressure_file.finaltime     
     1140           
     1141
     1142            self.index=0;
     1143            for i in range(len(self.file_time)):
     1144                if (self.file_time[i]<domain.starttime):
     1145                    self.index=i
     1146                else:
     1147                    break
     1148
     1149            N = domain.get_number_of_nodes()
     1150            self.prev_pressure_vertex_values=num.empty(N,num.float)
     1151            self.next_pressure_vertex_values=num.empty(N,num.float)
     1152            for i in range(N):
     1153                self.prev_pressure_vertex_values[i]=self.pressure(self.file_time[self.index],i)
     1154                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i)
     1155
     1156        self.p_vec=num.empty(N,num.float)
     1157
     1158
     1159    ##
     1160    # @brief 'execute' this class instance.
     1161    # @param domain
     1162    def __call__(self, domain):
     1163        """Evaluate pressure field based on values found in domain"""
     1164
     1165        xmom_update = domain.quantities['xmomentum'].explicit_update
     1166        ymom_update = domain.quantities['ymomentum'].explicit_update
     1167
     1168        t = domain.time
     1169
     1170        if callable(self.pressure):
     1171            if ( self.use_coordinates ):
     1172                xv = domain.get_nodes()
     1173                self.p_vec = self.pressure(t, xv[:,0], xv[:,1])
     1174            else:
     1175                self.update_stored_pressure_values(domain)
     1176
     1177                # Linear temporal interpolation of pressure values
     1178                ratio = (t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index])
     1179                self.p_vec = self.prev_pressure_vertex_values + ratio*(self.next_pressure_vertex_values - self.prev_pressure_vertex_values)
     1180
     1181        else:
     1182            # Assume s is a scalar
     1183            try:
     1184                self.p_vec[:] = self.pressure
     1185            except:
     1186                msg = 'Pressure must be either callable function or a scalar: %s' %self.s
     1187                raise msg
     1188
     1189        stage = domain.quantities['stage']
     1190        elevation = domain.quantities['elevation']
     1191
     1192        height = stage.centroid_values - elevation.centroid_values
     1193
     1194        point = domain.get_vertex_coordinates()
     1195
     1196        assign_pressure_field_values(height, self.p_vec, point,
     1197                                     domain.triangles,
     1198                                     xmom_update, ymom_update)
     1199
     1200    def update_stored_pressure_values(self,domain):
     1201        while (self.file_time[self.index+1]<domain.time):
     1202            self.index+=1
     1203            self.prev_pressure_vertex_values=copy(self.next_pressure_vertex_values)
     1204            for i in range(self.prev_pressure_vertex_values.shape[0]):
     1205                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i)
     1206
     1207
     1208class Wind_stress_fast:
     1209    """ Apply wind stress to water momentum in terms of
     1210        wind speed [m/s] and wind direction [degrees].
     1211        If the wind data is stored in a file
     1212        file_function is used to create a callable function. The data file
     1213        contains wind speed and direction values at a set of possibly
     1214        arbitrarily located nodes
     1215        at a set of possibly irregular but increasing times. file_function
     1216        interpolates from the file data onto the centroids of the domain.mesh
     1217        for each time. Two arrays for each wind quantity are then stored \
     1218        q0=q(t0,:) and q1=q(t1,:)
     1219        where t0<=domain.time<=t1. These arrays are recalculated when necessary
     1220        i.e t>t1. A linear temporal interpolation is used to approximate
     1221        pressure at time t.
     1222    """
     1223    def __init__(self, *args, **kwargs):
     1224        """Initialise windfield from wind speed s [m/s]
     1225        and wind direction phi [degrees]
     1226
     1227        Inputs v and phi can be either scalars or Python functions, e.g.
     1228
     1229        W = Wind_stress(10, 178)
     1230
     1231        #FIXME - 'normal' degrees are assumed for now, i.e. the
     1232        vector (1,0) has zero degrees.
     1233        We may need to convert from 'compass' degrees later on and also
     1234        map from True north to grid north.
     1235
     1236        Arguments can also be Python functions of t,x,y as in
     1237
     1238        def speed(t,x,y):
     1239            ...
     1240            return s
     1241
     1242        def angle(t,x,y):
     1243            ...
     1244            return phi
     1245
     1246        where x and y are vectors.
     1247
     1248        and then pass the functions in
     1249
     1250        W = Wind_stress(speed, angle)
     1251
     1252        The instantiated object W can be appended to the list of
     1253        forcing_terms as in
     1254
     1255        Alternatively, one vector valued function for (speed, angle)
     1256        can be applied, providing both quantities simultaneously.
     1257        As in
     1258        W = Wind_stress(F), where returns (speed, angle) for each t.
     1259
     1260        domain.forcing_terms.append(W)
     1261        """
     1262
     1263        from anuga.config import rho_a, rho_w, eta_w
     1264
     1265        self.use_coordinates=True
     1266        if len(args) == 2:
     1267            s = args[0]
     1268            phi = args[1]
     1269        elif len(args) == 1:
     1270            # Assume vector function returning (s, phi)(t,x,y)
     1271            vector_function = args[0]
     1272            if ( len(kwargs)==2 ):
     1273                filename=kwargs['filename']
     1274                domain=kwargs['domain']
     1275                self.use_coordinates=False
     1276            else:
     1277                self.use_coordinates=True
     1278            if ( self.use_coordinates ):
     1279                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
     1280                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
     1281            else:
     1282                s = lambda t,i: vector_function(t,point_id=i)[0]
     1283                phi = lambda t,i: vector_function(t,point_id=i)[1]
     1284        else:
     1285           # Assume info is in 2 keyword arguments
     1286           if len(kwargs) == 2:
     1287               s = kwargs['s']
     1288               phi = kwargs['phi']
     1289           else:
     1290               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
     1291
     1292        if ( self.use_coordinates ):
     1293            self.speed = check_forcefield(s)
     1294            self.phi = check_forcefield(phi)
     1295        else:
     1296            self.speed = s
     1297            self.phi = phi
     1298
     1299        N = len(domain)
     1300        if ( not self.use_coordinates):
     1301
     1302            # Open NetCDF file
     1303            fid = NetCDFFile(filename, netcdf_mode_r)
     1304            self.file_time = fid.variables['time'][:]
     1305            fid.close()
     1306
     1307            msg = 'wind_file.starttime > domain.starttime'
     1308            if (self.file_time[0]>domain.starttime):
     1309                raise Exception, msg
     1310
     1311            msg = 'wind_file[-1] < domain.starttime'
     1312            if (self.file_time[-1]<domain.starttime):
     1313                raise Exception, msg
     1314
     1315            msg = 'No wind values exist for times greater than domain.starttime'
     1316            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
     1317                raise Exception, msg
     1318
     1319            # FIXME(JJ): How do we check that evolve
     1320            # finaltime  < wind_file.finaltime     
     1321           
     1322
     1323            self.index=0;
     1324            for i in range(len(self.file_time)):
     1325                if (self.file_time[i]<domain.starttime):
     1326                    self.index=i
     1327                else:
     1328                    break
     1329
     1330            self.prev_windspeed_centroid_values=num.empty(N,num.float)
     1331            self.next_windspeed_centroid_values=num.empty(N,num.float)
     1332            self.prev_windangle_centroid_values=num.empty(N,num.float)
     1333            self.next_windangle_centroid_values=num.empty(N,num.float)
     1334            for i in range(N):
     1335                self.prev_windspeed_centroid_values[i]=self.speed(self.file_time[self.index],i)
     1336                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i)
     1337                self.prev_windangle_centroid_values[i]=self.phi(self.file_time[self.index],i)
     1338                self.next_windangle_centroid_values[i]=self.phi(self.file_time[self.index+1],i)
     1339
     1340        self.s_vec=num.empty(N,num.float)
     1341        self.phi_vec=num.empty(N,num.float)
     1342
     1343        self.const = eta_w*rho_a/rho_w
     1344    ##
     1345    # @brief 'execute' this class instance.
     1346    # @param domain
     1347    def __call__(self, domain):
     1348        """Evaluate windfield based on values found in domain"""
     1349
     1350        xmom_update = domain.quantities['xmomentum'].explicit_update
     1351        ymom_update = domain.quantities['ymomentum'].explicit_update
     1352
     1353        N = len(domain)    # number_of_triangles
     1354        t = domain.time
     1355
     1356        if callable(self.speed):
     1357            if ( self.use_coordinates ):
     1358                xc = domain.get_centroid_coordinates()
     1359                self.s_vec = self.speed(t, xc[:,0], xc[:,1])
     1360            else:
     1361                self.update_stored_wind_values(domain)
     1362
     1363                # Linear temporal interpolation of wind values
     1364                if t==self.file_time[self.index]:
     1365                    ratio = 0.
     1366                else:
     1367                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
     1368                self.s_vec = self.prev_windspeed_centroid_values + ratio*(self.next_windspeed_centroid_values - self.prev_windspeed_centroid_values)
     1369        else:
     1370            # Assume s is a scalar
     1371            try:
     1372                self.s_vec[:] = self.speed
     1373            except:
     1374                msg = 'Speed must be either callable or a scalar: %s' %self.s
     1375                raise msg
     1376
     1377        if callable(self.phi):
     1378            if ( self.use_coordinates ):
     1379                xc = domain.get_centroid_coordinates()
     1380                self.phi_vec = self.phi(t, xc[:,0], xc[:,1])
     1381            else:
     1382                self.update_stored_wind_values(domain)
     1383
     1384                # Linear temporal interpolation of wind values
     1385                if t==self.file_time[self.index]:
     1386                    ratio = 0.
     1387                else:
     1388                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
     1389                self.phi_vec = self.prev_windangle_centroid_values + ratio*(self.next_windangle_centroid_values - self.prev_windangle_centroid_values)
     1390        else:
     1391            # Assume phi is a scalar
     1392
     1393            try:
     1394                self.phi_vec[:] = self.phi
     1395            except:
     1396                msg = 'Angle must be either callable or a scalar: %s' %self.phi
     1397                raise msg
     1398
     1399        assign_windfield_values(xmom_update, ymom_update,
     1400                                self.s_vec, self.phi_vec, self.const)
     1401
     1402    def update_stored_wind_values(self,domain):
     1403        while (self.file_time[self.index+1]<domain.time):
     1404            self.index+=1
     1405            self.prev_windspeed_centroid_values=copy(self.next_windspeed_centroid_values)
     1406            self.prev_windangle_centroid_values=copy(self.next_windangle_centroid_values)
     1407            for i in range(self.next_windspeed_centroid_values.shape[0]):
     1408                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i)
     1409                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

    r8068 r8116  
    33
    44import unittest, os
     5import anuga
    56from anuga.shallow_water.shallow_water_domain import Domain
    67from boundaries import Reflective_boundary
     
    89from anuga.file_conversion.file_conversion import timefile2netcdf
    910from forcing import *
     11from mesh_factory import rectangular
     12from file_conversion.sts2sww_mesh import sts2sww_mesh
     13from anuga.abstract_2d_finite_volumes.util import file_function
    1014
    1115import numpy as num
     
    7680    return a
    7781
     82def time_varying_speed(t, x, y):
     83    """
     84    Variable speed windfield
     85    """
     86
     87    from math import exp, cos, pi
     88
     89    x = num.array(x,num.float)
     90    y = num.array(y,num.float)
     91
     92    N = len(x)
     93    s = 0*x  #New array
     94
     95    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
     96    S=100.
     97    for k in range(N):
     98        s[k]=S*(1.+t/100.)
     99    return s
     100
     101
     102def time_varying_angle(t, x, y):
     103    """Rotating field
     104    """
     105    from math import atan, pi
     106
     107    x = num.array(x,num.float)
     108    y = num.array(y,num.float)
     109
     110    N = len(x)
     111    a = 0 * x    # New array
     112
     113    phi=135.
     114    for k in range(N):
     115        a[k]=phi*(1.+t/100.)
     116
     117    return a
     118
     119
     120def time_varying_pressure(t, x, y):
     121    """Rotating field
     122    """
     123    from math import atan, pi
     124
     125    x = num.array(x,num.float)
     126    y = num.array(y,num.float)
     127
     128    N = len(x)
     129    p = 0 * x    # New array
     130
     131    p0=1000.
     132    for k in range(N):
     133        p[k]=p0*(1.-t/100.)
     134
     135    return p
     136
     137def spatial_linear_varying_speed(t, x, y):
     138    """
     139    Variable speed windfield
     140    """
     141
     142    from math import exp, cos, pi
     143
     144    x = num.array(x)
     145    y = num.array(y)
     146
     147    N = len(x)
     148    s = 0*x  #New array
     149
     150    #dx=x[-1]-x[0]; dy = y[-1]-y[0]
     151    s0=250.
     152    ymin=num.min(y)
     153    xmin=num.min(x)
     154    a=0.000025; b=0.0000125;
     155    for k in range(N):
     156        s[k]=s0*(1+t/100.)+a*x[k]+b*y[k]
     157    return s
     158
     159
     160def spatial_linear_varying_angle(t, x, y):
     161    """Rotating field
     162    """
     163    from math import atan, pi
     164
     165    x = num.array(x)
     166    y = num.array(y)
     167
     168    N = len(x)
     169    a = 0 * x    # New array
     170
     171    phi=135.
     172    b1=0.000025; b2=0.00001125;
     173    for k in range(N):
     174        a[k]=phi*(1+t/100.)+b1*x[k]+b2*y[k]
     175    return a
     176
     177def spatial_linear_varying_pressure(t, x, y):
     178    p0=1000;
     179    a=0.000025; b=0.0000125;
     180
     181    x = num.array(x)
     182    y = num.array(y)
     183
     184    N = len(x)
     185    p = 0 * x    # New array
     186
     187    for k in range(N):
     188        p[k]=p0*(1.-t/100.)+a*x[k]+b*y[k]
     189    return p
     190
     191
     192def grid_1d(x0,dx,nx):
     193    x = num.empty(nx,dtype=num.float)
     194    for i in range(nx):
     195        x[i]=x0+float(i)*dx
     196    return x
    78197   
     198
     199def ndgrid(x,y):
     200    nx = len(x)
     201    ny = len(y)
     202    X = num.empty(nx*ny,dtype=num.float)
     203    Y = num.empty(nx*ny,dtype=num.float)
     204    k=0
     205    for i in range(nx):
     206        for j in range(ny):
     207            X[k]=x[i]
     208            Y[k]=y[j]
     209            k+=1
     210    return X,Y
    79211
    80212class Test_Forcing(unittest.TestCase):
     
    85217        pass
    86218       
     219    def write_wind_pressure_field_sts(self,
     220                                      field_sts_filename,
     221                                      nrows=10,
     222                                      ncols=10,
     223                                      cellsize=25,
     224                                      origin=(0.0,0.0),
     225                                      refzone=50,
     226                                      timestep=1,
     227                                      number_of_timesteps=10,
     228                                      angle=135.0,
     229                                      speed=100.0,
     230                                      pressure=1000.0):
     231
     232        xllcorner=origin[0]
     233        yllcorner=origin[1]
     234        starttime = 0; endtime = number_of_timesteps*timestep;
     235        no_data = -9999
     236
     237        time = num.arange(starttime, endtime, timestep, dtype='i')
     238
     239        x = grid_1d(xllcorner,cellsize,ncols)
     240        y = grid_1d(yllcorner,cellsize,nrows)
     241        [X,Y] = ndgrid(x,y)
     242        number_of_points = nrows*ncols
     243
     244        wind_speed = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
     245        wind_angle = num.empty((number_of_timesteps,nrows*ncols),dtype=num.float)
     246        barometric_pressure = num.empty((number_of_timesteps,nrows*ncols),
     247                                        dtype=num.float)
     248
     249        if ( callable(speed) and callable(angle) and callable(pressure) ):
     250            x = num.ones(3, num.float)
     251            y = num.ones(3, num.float)
     252            try:
     253                s = speed(1.0, x=x, y=y)
     254                a = angle(1.0, x=x, y=y)
     255                p = pressure(1.0, x=x, y=y)
     256                use_function=True
     257            except Exception, e:
     258                msg = 'Function could not be executed.\n'
     259                raise Exception, msg
     260        else:
     261            try :
     262                speed=float(speed)
     263                angle=float(angle)
     264                pressure=float(pressure)
     265                use_function=False
     266            except:
     267                msg = ('Force fields must be a scalar value coercible to float.')
     268                raise Exception, msg
     269
     270        for i,t in enumerate(time):
     271            if ( use_function ):
     272                wind_speed[i,:] = speed(t,X,Y)
     273                wind_angle[i,:] = angle(t,X,Y)
     274                barometric_pressure[i,:] = pressure(t,X,Y)
     275            else:
     276                wind_speed[i,:] = speed
     277                wind_angle[i,:] = angle
     278                barometric_pressure[i,:] = pressure
     279
     280        # "Creating the field STS NetCDF file"
     281
     282        fid = NetCDFFile(field_sts_filename+'.sts', 'w')
     283        fid.institution = 'Geoscience Australia'
     284        fid.description = "description"
     285        fid.starttime = 0.0
     286        fid.ncols = ncols
     287        fid.nrows = nrows
     288        fid.cellsize = cellsize
     289        fid.no_data = no_data
     290        fid.createDimension('number_of_points', number_of_points)
     291        fid.createDimension('number_of_timesteps', number_of_timesteps)
     292        fid.createDimension('numbers_in_range', 2)
     293
     294        fid.createVariable('x', 'd', ('number_of_points',))
     295        fid.createVariable('y', 'd', ('number_of_points',))
     296        fid.createVariable('time', 'i', ('number_of_timesteps',))
     297        fid.createVariable('wind_speed', 'd', ('number_of_timesteps',
     298                                               'number_of_points'))
     299        fid.createVariable('wind_speed_range', 'd', ('numbers_in_range', ))
     300        fid.createVariable('wind_angle', 'd', ('number_of_timesteps',
     301                                               'number_of_points'))
     302        fid.createVariable('wind_angle_range', 'd', ('numbers_in_range',))
     303        fid.createVariable('barometric_pressure', 'd', ('number_of_timesteps',
     304                                             'number_of_points'))
     305        fid.createVariable('barometric_pressure_range', 'd', ('numbers_in_range',))
     306
     307
     308        fid.variables['wind_speed_range'][:] = num.array([1e+036, -1e+036])
     309        fid.variables['wind_angle_range'][:] = num.array([1e+036, -1e+036])
     310        fid.variables['barometric_pressure_range'][:] = num.array([1e+036, -1e+036])
     311        fid.variables['time'][:] = time
     312
     313        ws = fid.variables['wind_speed']
     314        wa = fid.variables['wind_angle']
     315        pr = fid.variables['barometric_pressure']
     316
     317        for i in xrange(number_of_timesteps):
     318            ws[i] = wind_speed[i,:]
     319            wa[i] = wind_angle[i,:]
     320            pr[i] = barometric_pressure[i,:]
     321
     322        origin = anuga.coordinate_transforms.geo_reference.Geo_reference(refzone,
     323                                                                         xllcorner,
     324                                                                         yllcorner)
     325        geo_ref = anuga.coordinate_transforms.geo_reference.write_NetCDF_georeference(origin, fid)
     326
     327        fid.variables['x'][:]=X-geo_ref.get_xllcorner()
     328        fid.variables['y'][:]=Y-geo_ref.get_yllcorner()
     329
     330
     331        fid.close()
     332
    87333    def test_constant_wind_stress(self):
    88334        from anuga.config import rho_a, rho_w, eta_w
     
    8551101            raise Exception, 'Should have raised exception'
    8561102
     1103    def test_constant_wind_stress_from_file(self):
     1104        from anuga.config import rho_a, rho_w, eta_w
     1105        from math import pi, cos, sin
     1106
     1107        cellsize = 25
     1108        nrows=5; ncols = 6;
     1109        refzone=50
     1110        xllcorner=366000;yllcorner=6369500;
     1111        number_of_timesteps = 6
     1112        timestep=12*60
     1113        eps=2e-16
     1114
     1115        points, vertices, boundary =rectangular(nrows-2,ncols-2,
     1116                                                len1=cellsize*(ncols-1),
     1117                                                len2=cellsize*(nrows-1),
     1118                                                origin=(xllcorner,yllcorner))
     1119
     1120        domain = Domain(points, vertices, boundary)
     1121        midpoints = domain.get_centroid_coordinates()
     1122
     1123        # Flat surface with 1m of water
     1124        domain.set_quantity('elevation', 0)
     1125        domain.set_quantity('stage', 1.0)
     1126        domain.set_quantity('friction', 0)
     1127
     1128        Br = Reflective_boundary(domain)
     1129        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1130
     1131        # Setup only one forcing term, constant wind stress
     1132        s = 100
     1133        phi = 135
     1134        pressure=1000
     1135        domain.forcing_terms = []
     1136        field_sts_filename = 'wind_field'
     1137        self.write_wind_pressure_field_sts(field_sts_filename,
     1138                                      nrows=nrows,
     1139                                      ncols=ncols,
     1140                                      cellsize=cellsize,
     1141                                      origin=(xllcorner,yllcorner),
     1142                                      refzone=50,
     1143                                      timestep=timestep,
     1144                                      number_of_timesteps=10,
     1145                                      speed=s,
     1146                                      angle=phi,
     1147                                      pressure=pressure)
     1148
     1149        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1150                     verbose=False)
     1151
     1152        # Setup wind stress
     1153        F = file_function(field_sts_filename+'.sww', domain,
     1154                          quantities=['wind_speed', 'wind_angle'],
     1155                          interpolation_points = midpoints)
     1156
     1157        W = Wind_stress(F,use_coordinates=False)
     1158        domain.forcing_terms.append(W)
     1159        domain.compute_forcing_terms()
     1160
     1161        const = eta_w*rho_a / rho_w
     1162
     1163        # Convert to radians
     1164        phi = phi*pi / 180
     1165
     1166        # Compute velocity vector (u, v)
     1167        u = s*cos(phi)
     1168        v = s*sin(phi)
     1169
     1170        # Compute wind stress
     1171        S = const * num.sqrt(u**2 + v**2)
     1172
     1173        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
     1174        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
     1175        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
     1176
     1177    def test_variable_windfield_from_file(self):
     1178        from anuga.config import rho_a, rho_w, eta_w
     1179        from math import pi, cos, sin
     1180        from anuga.config import time_format
     1181
     1182        cellsize = 25
     1183        #nrows=25; ncols = 25;
     1184        nrows=10; ncols = 10;
     1185        refzone=50
     1186        xllcorner=366000;yllcorner=6369500;
     1187        number_of_timesteps = 10
     1188        timestep=1
     1189        eps=2.e-16
     1190        spatial_thinning=1
     1191
     1192        points, vertices, boundary =rectangular(nrows-2,ncols-2,
     1193                                                len1=cellsize*(ncols-1),
     1194                                                len2=cellsize*(nrows-1),
     1195                                                origin=(xllcorner,yllcorner))
     1196
     1197        time=num.arange(0,10,1,num.float)
     1198        eval_time=time[7];
     1199
     1200        domain = Domain(points, vertices, boundary)
     1201        midpoints = domain.get_centroid_coordinates()
     1202        vertexpoints = domain.get_nodes()
     1203
     1204        """
     1205        x=grid_1d(xllcorner,cellsize,ncols)
     1206        y=grid_1d(yllcorner,cellsize,nrows)
     1207        X,Y=num.meshgrid(x,y)
     1208        interpolation_points=num.empty((X.shape[0]*X.shape[1],2),num.float)
     1209        k=0
     1210        for i in range(X.shape[0]):
     1211            for j in range(X.shape[1]):
     1212                interpolation_points[k,0]=X[i,j]
     1213                interpolation_points[k,1]=Y[i,j]
     1214                k+=1
     1215
     1216        z=spatial_linear_varying_speed(eval_time,interpolation_points[:,0],
     1217                                       interpolation_points[:,1])
     1218
     1219        k=0
     1220        Z=num.empty((X.shape[0],X.shape[1]),num.float)
     1221        for i in range(X.shape[0]):
     1222            for j in range(X.shape[1]):
     1223                Z[i,j]=z[k]
     1224                k+=1
     1225
     1226        Q=num.empty((time.shape[0],points.shape[0]),num.float)
     1227        for i, t in enumerate(time):
     1228            Q[i,:]=spatial_linear_varying_speed(t,points[:,0],points[:,1])
     1229
     1230        from interpolate import Interpolation_function
     1231        I  = Interpolation_function(time,Q,
     1232                                    vertex_coordinates = points,
     1233                                    triangles = domain.triangles,
     1234                                    #interpolation_points = midpoints,
     1235                                    interpolation_points=interpolation_points,
     1236                                    verbose=False)
     1237
     1238        V=num.empty((X.shape[0],X.shape[1]),num.float)
     1239        for k in range(len(interpolation_points)):
     1240            assert num.allclose(I(eval_time,k),z[k])
     1241            V[k/X.shape[1],k%X.shape[1]]=I(eval_time,k)
     1242
     1243
     1244           import mpl_toolkits.mplot3d.axes3d as p3
     1245           fig=P.figure()
     1246           ax = p3.Axes3D(fig)
     1247           ax.plot_surface(X,Y,V)
     1248           ax.plot_surface(X,Y,Z)
     1249           P.show()
     1250
     1251
     1252        """
     1253
     1254        # Flat surface with 1m of water
     1255        domain.set_quantity('elevation', 0)
     1256        domain.set_quantity('stage', 1.0)
     1257        domain.set_quantity('friction', 0)
     1258
     1259        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
     1260
     1261        # Write wind stress file (ensure that domain.time is covered)
     1262
     1263        field_sts_filename = 'wind_field'
     1264        self.write_wind_pressure_field_sts(field_sts_filename,
     1265                                      nrows=nrows,
     1266                                      ncols=ncols,
     1267                                      cellsize=cellsize,
     1268                                      origin=(xllcorner,yllcorner),
     1269                                      refzone=50,
     1270                                      timestep=timestep,
     1271                                      number_of_timesteps=10,
     1272                                      speed=spatial_linear_varying_speed,
     1273                                      angle=spatial_linear_varying_angle,
     1274                                      pressure=spatial_linear_varying_pressure)
     1275
     1276
     1277        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
     1278                     verbose=False)
     1279
     1280        # Setup wind stress
     1281        FW = file_function(field_sts_filename+'.sww', domain,
     1282                          quantities=['wind_speed', 'wind_angle'],
     1283                          interpolation_points = midpoints)
     1284
     1285        W = Wind_stress(FW,use_coordinates=False)
     1286
     1287        domain.forcing_terms = []
     1288        domain.forcing_terms.append(W)
     1289
     1290        domain.compute_forcing_terms()
     1291
     1292        # Compute reference solution
     1293        const = eta_w*rho_a / rho_w
     1294
     1295        N = len(domain)    # number_of_triangles
     1296
     1297        xc = domain.get_centroid_coordinates()
     1298        t = domain.time
     1299
     1300        x = xc[:,0]
     1301        y = xc[:,1]
     1302        s_vec = spatial_linear_varying_speed(t,x,y)
     1303        phi_vec = spatial_linear_varying_angle(t,x,y)
     1304
     1305        for k in range(N):
     1306            # Convert to radians
     1307            phi = phi_vec[k]*pi / 180
     1308            s = s_vec[k]
     1309
     1310            # Compute velocity vector (u, v)
     1311            u = s*cos(phi)
     1312            v = s*sin(phi)
     1313
     1314            # Compute wind stress
     1315            S = const * num.sqrt(u**2 + v**2)
     1316
     1317            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
     1318
     1319            assert num.allclose(domain.quantities['xmomentum'].\
     1320                                    explicit_update[k],S*u,eps)
     1321            assert num.allclose(domain.quantities['ymomentum'].\
     1322                                     explicit_update[k],S*v,eps)
     1323
     1324        os.remove(field_sts_filename+'.sts')
     1325        os.remove(field_sts_filename+'.sww')
     1326
     1327    def test_variable_pressurefield_from_file(self):
     1328        from anuga.config import rho_a, rho_w, eta_w
     1329        from math import pi, cos, sin
     1330        from anuga.config import time_format
     1331
     1332        cellsize = 25
     1333        #nrows=25; ncols = 25;
     1334        nrows=10; ncols = 10;
     1335        refzone=50
     1336        xllcorner=366000;yllcorner=6369500;
     1337        number_of_timesteps = 10
     1338        timestep=1
     1339        eps=2.e-16
     1340        spatial_thinning=1
     1341
     1342        points, vertices, boundary =rectangular(nrows-2,ncols-2,
     1343                                                len1=cellsize*(ncols-1),
     1344                                                len2=cellsize*(nrows-1),
     1345                                                origin=(xllcorner,yllcorner))
     1346
     1347        time=num.arange(0,10,1,num.float)
     1348        eval_time=time[7];
     1349
     1350        domain = Domain(points, vertices, boundary)
     1351        midpoints = domain.get_centroid_coordinates()
     1352        vertexpoints = domain.get_nodes()
     1353
     1354        # Flat surface with 1m of water
     1355        domain.set_quantity('elevation', 0)
     1356        domain.set_quantity('stage', 1.0)
     1357        domain.set_quantity('friction', 0)
     1358
     1359        domain.time = 7*timestep    # Take a time that is represented in file (not zero)
     1360
     1361        # Write wind stress file (ensure that domain.time is covered)
     1362
     1363        field_sts_filename = 'wind_field'
     1364        self.write_wind_pressure_field_sts(field_sts_filename,
     1365                                      nrows=nrows,
     1366                                      ncols=ncols,
     1367                                      cellsize=cellsize,
     1368                                      origin=(xllcorner,yllcorner),
     1369                                      refzone=50,
     1370                                      timestep=timestep,
     1371                                      number_of_timesteps=10,
     1372                                      speed=spatial_linear_varying_speed,
     1373                                      angle=spatial_linear_varying_angle,
     1374                                      pressure=spatial_linear_varying_pressure)
     1375
     1376
     1377        sts2sww_mesh(field_sts_filename,spatial_thinning=spatial_thinning,
     1378                     verbose=False)
     1379
     1380        # Setup barometric pressure
     1381        FP = file_function(field_sts_filename+'.sww', domain,
     1382                           quantities=['barometric_pressure'],
     1383                           interpolation_points = vertexpoints)
     1384
     1385        P = Barometric_pressure(FP,use_coordinates=False)
     1386
     1387
     1388        domain.forcing_terms = []
     1389        domain.forcing_terms.append(P)
     1390
     1391        domain.compute_forcing_terms()
     1392
     1393        N = len(domain)    # number_of_triangles
     1394
     1395        xc = domain.get_centroid_coordinates()
     1396        t = domain.time
     1397
     1398        x = xc[:,0]
     1399        y = xc[:,1]
     1400        p_vec = spatial_linear_varying_pressure(t,x,y)
     1401
     1402        h=1 #depth
     1403        px=0.000025  #pressure gradient in x-direction
     1404        py=0.0000125 #pressure gradient in y-direction
     1405        for k in range(N):
     1406            # Convert to radians
     1407            p = p_vec[k]
     1408
     1409            assert num.allclose(domain.quantities['stage'].explicit_update[k],0)
     1410
     1411            assert num.allclose(domain.quantities['xmomentum'].\
     1412                                    explicit_update[k],h*px/rho_w)
     1413
     1414            assert num.allclose(domain.quantities['ymomentum'].\
     1415                                     explicit_update[k],h*py/rho_w)
     1416
     1417        os.remove(field_sts_filename+'.sts')
     1418        os.remove(field_sts_filename+'.sww')
     1419
     1420    def test_constant_wind_stress_from_file_evolve(self):
     1421        from anuga.config import rho_a, rho_w, eta_w
     1422        from math import pi, cos, sin
     1423        from anuga.config import time_format
     1424
     1425        cellsize = 25
     1426        nrows=5; ncols = 6;
     1427        refzone=50
     1428        xllcorner=366000;yllcorner=6369500;
     1429        number_of_timesteps = 27
     1430        timestep=1
     1431        eps=2e-16
     1432
     1433        points, vertices, boundary =rectangular(nrows-2,ncols-2,
     1434                                                len1=cellsize*(ncols-1),
     1435                                                len2=cellsize*(nrows-1),
     1436                                                origin=(xllcorner,yllcorner))
     1437
     1438        domain = Domain(points, vertices, boundary)
     1439        midpoints = domain.get_centroid_coordinates()
     1440
     1441        # Flat surface with 1m of water
     1442        domain.set_quantity('elevation', 0)
     1443        domain.set_quantity('stage', 1.0)
     1444        domain.set_quantity('friction', 0)
     1445
     1446        Br = Reflective_boundary(domain)
     1447        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1448
     1449        # Setup only one forcing term, constant wind stress
     1450        s = 100
     1451        phi = 135
     1452        field_sts_filename = 'wind_field'
     1453        self.write_wind_pressure_field_sts(field_sts_filename,
     1454                                      nrows=nrows,
     1455                                      ncols=ncols,
     1456                                      cellsize=cellsize,
     1457                                      origin=(xllcorner,yllcorner),
     1458                                      refzone=50,
     1459                                      timestep=timestep,
     1460                                      number_of_timesteps=number_of_timesteps,
     1461                                      speed=s,
     1462                                      angle=phi)
     1463
     1464        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1465                     verbose=False)
     1466
     1467        # Setup wind stress
     1468        F = file_function(field_sts_filename+'.sww', domain,
     1469                          quantities=['wind_speed', 'wind_angle'],
     1470                          interpolation_points = midpoints)
     1471
     1472        W = Wind_stress(F,use_coordinates=False)
     1473        domain.forcing_terms.append(W)
     1474
     1475        valuesUsingFunction=num.empty((3,number_of_timesteps+1,midpoints.shape[0]),
     1476                                      num.float)
     1477        i=0
     1478        for t in domain.evolve(yieldstep=1, finaltime=number_of_timesteps*timestep):
     1479            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
     1480            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
     1481            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
     1482            i+=1
     1483
     1484
     1485        domain_II = Domain(points, vertices, boundary)
     1486
     1487        # Flat surface with 1m of water
     1488        domain_II.set_quantity('elevation', 0)
     1489        domain_II.set_quantity('stage', 1.0)
     1490        domain_II.set_quantity('friction', 0)
     1491
     1492        Br = Reflective_boundary(domain_II)
     1493        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1494
     1495        s = 100
     1496        phi = 135
     1497        domain_II.forcing_terms = []
     1498        domain_II.forcing_terms.append(Wind_stress(s, phi))
     1499
     1500        i=0;
     1501        for t in domain_II.evolve(yieldstep=1,
     1502                                  finaltime=number_of_timesteps*timestep):
     1503            assert num.allclose(valuesUsingFunction[0,i],domain_II.quantities['stage'].explicit_update), max(valuesUsingFunction[0,i]-domain_II.quantities['stage'].explicit_update)
     1504            assert  num.allclose(valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update)
     1505            assert num.allclose(valuesUsingFunction[2,i],domain_II.quantities['ymomentum'].explicit_update)
     1506            i+=1
     1507
     1508        os.remove(field_sts_filename+'.sts')
     1509        os.remove(field_sts_filename+'.sww')
     1510
     1511    def test_temporally_varying_wind_stress_from_file_evolve(self):
     1512        from anuga.config import rho_a, rho_w, eta_w
     1513        from math import pi, cos, sin
     1514        from anuga.config import time_format
     1515
     1516        cellsize = 25
     1517        #nrows=20; ncols = 20;
     1518        nrows=10; ncols = 10;
     1519        refzone=50
     1520        xllcorner=366000;yllcorner=6369500;
     1521        number_of_timesteps = 28
     1522        timestep=1.
     1523        eps=2e-16
     1524
     1525        #points, vertices, boundary =rectangular(10,10,
     1526        points, vertices, boundary =rectangular(5,5,
     1527                                                len1=cellsize*(ncols-1),
     1528                                                len2=cellsize*(nrows-1),
     1529                                                origin=(xllcorner,yllcorner))
     1530
     1531        domain = Domain(points, vertices, boundary)
     1532        midpoints = domain.get_centroid_coordinates()
     1533
     1534        # Flat surface with 1m of water
     1535        domain.set_quantity('elevation', 0)
     1536        domain.set_quantity('stage', 1.0)
     1537        domain.set_quantity('friction', 0)
     1538
     1539        Br = Reflective_boundary(domain)
     1540        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1541
     1542        # Setup only one forcing term, constant wind stress
     1543        field_sts_filename = 'wind_field'
     1544        self.write_wind_pressure_field_sts(field_sts_filename,
     1545                                      nrows=nrows,
     1546                                      ncols=ncols,
     1547                                      cellsize=cellsize,
     1548                                      origin=(xllcorner,yllcorner),
     1549                                      refzone=50,
     1550                                      timestep=timestep,
     1551                                      number_of_timesteps=number_of_timesteps,
     1552                                      speed=time_varying_speed,
     1553                                      angle=time_varying_angle,
     1554                                      pressure=time_varying_pressure)
     1555
     1556        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1557                     verbose=False)
     1558
     1559        # Setup wind stress
     1560        F = file_function(field_sts_filename+'.sww', domain,
     1561                          quantities=['wind_speed', 'wind_angle'],
     1562                          interpolation_points = midpoints)
     1563
     1564        #W = Wind_stress(F,use_coordinates=False)
     1565        W = Wind_stress_fast(F,filename=field_sts_filename+'.sww', domain=domain)
     1566        domain.forcing_terms.append(W)
     1567
     1568        valuesUsingFunction=num.empty((3,2*number_of_timesteps,midpoints.shape[0]),
     1569                                      num.float)
     1570        i=0
     1571        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
     1572            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
     1573            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
     1574            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
     1575            i+=1
     1576
     1577
     1578        domain_II = Domain(points, vertices, boundary)
     1579
     1580        # Flat surface with 1m of water
     1581        domain_II.set_quantity('elevation', 0)
     1582        domain_II.set_quantity('stage', 1.0)
     1583        domain_II.set_quantity('friction', 0)
     1584
     1585        Br = Reflective_boundary(domain_II)
     1586        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1587
     1588        domain_II.forcing_terms.append(Wind_stress(s=time_varying_speed,
     1589                                                   phi=time_varying_angle))
     1590
     1591        i=0;
     1592        for t in domain_II.evolve(yieldstep=timestep/2.,
     1593                                  finaltime=(number_of_timesteps-1)*timestep):
     1594            assert num.allclose(valuesUsingFunction[0,i],
     1595                                domain_II.quantities['stage'].explicit_update,
     1596                                eps)
     1597            #print i,valuesUsingFunction[1,i]
     1598            assert  num.allclose(valuesUsingFunction[1,i],
     1599                                 domain_II.quantities['xmomentum'].explicit_update,
     1600                                 eps),(valuesUsingFunction[1,i]-
     1601                                 domain_II.quantities['xmomentum'].explicit_update)
     1602            assert num.allclose(valuesUsingFunction[2,i],
     1603                                domain_II.quantities['ymomentum'].explicit_update,
     1604                                eps)
     1605            #if i==1: assert-1==1
     1606            i+=1
     1607
     1608        os.remove(field_sts_filename+'.sts')
     1609        os.remove(field_sts_filename+'.sww')
     1610
     1611    def test_spatially_varying_wind_stress_from_file_evolve(self):
     1612        from anuga.config import rho_a, rho_w, eta_w
     1613        from math import pi, cos, sin
     1614        from anuga.config import time_format
     1615
     1616        cellsize = 25
     1617        nrows=20; ncols = 20;
     1618        nrows=10; ncols = 10;
     1619        refzone=50
     1620        xllcorner=366000;yllcorner=6369500;
     1621        number_of_timesteps = 28
     1622        timestep=1.
     1623        eps=2e-16
     1624
     1625        #points, vertices, boundary =rectangular(10,10,
     1626        points, vertices, boundary =rectangular(5,5,
     1627                                                len1=cellsize*(ncols-1),
     1628                                                len2=cellsize*(nrows-1),
     1629                                                origin=(xllcorner,yllcorner))
     1630
     1631        domain = Domain(points, vertices, boundary)
     1632        midpoints = domain.get_centroid_coordinates()
     1633
     1634        # Flat surface with 1m of water
     1635        domain.set_quantity('elevation', 0)
     1636        domain.set_quantity('stage', 1.0)
     1637        domain.set_quantity('friction', 0)
     1638
     1639        Br = Reflective_boundary(domain)
     1640        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1641
     1642        # Setup only one forcing term, constant wind stress
     1643        field_sts_filename = 'wind_field'
     1644        self.write_wind_pressure_field_sts(field_sts_filename,
     1645                                      nrows=nrows,
     1646                                      ncols=ncols,
     1647                                      cellsize=cellsize,
     1648                                      origin=(xllcorner,yllcorner),
     1649                                      refzone=50,
     1650                                      timestep=timestep,
     1651                                      number_of_timesteps=number_of_timesteps,
     1652                                      speed=spatial_linear_varying_speed,
     1653                                      angle=spatial_linear_varying_angle,
     1654                                      pressure=spatial_linear_varying_pressure)
     1655
     1656        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1657                     verbose=False)
     1658
     1659        # Setup wind stress
     1660        F = file_function(field_sts_filename+'.sww', domain,
     1661                          quantities=['wind_speed', 'wind_angle'],
     1662                          interpolation_points = midpoints)
     1663
     1664        W = Wind_stress(F,use_coordinates=False)
     1665        domain.forcing_terms.append(W)
     1666
     1667        valuesUsingFunction=num.empty((3,number_of_timesteps,midpoints.shape[0]),
     1668                                      num.float)
     1669        i=0
     1670        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
     1671            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
     1672            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
     1673            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
     1674            i+=1
     1675
     1676
     1677        domain_II = Domain(points, vertices, boundary)
     1678
     1679        # Flat surface with 1m of water
     1680        domain_II.set_quantity('elevation', 0)
     1681        domain_II.set_quantity('stage', 1.0)
     1682        domain_II.set_quantity('friction', 0)
     1683
     1684        Br = Reflective_boundary(domain_II)
     1685        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1686
     1687        domain_II.forcing_terms.append(Wind_stress(s=spatial_linear_varying_speed,
     1688                                                   phi=spatial_linear_varying_angle))
     1689
     1690        i=0;
     1691        for t in domain_II.evolve(yieldstep=timestep,
     1692                                  finaltime=(number_of_timesteps-1)*timestep):
     1693            #print valuesUsingFunction[1,i],domain_II.quantities['xmomentum'].explicit_update
     1694            assert num.allclose(valuesUsingFunction[0,i],
     1695                                domain_II.quantities['stage'].explicit_update,
     1696                                eps)
     1697            assert  num.allclose(valuesUsingFunction[1,i],
     1698                                 domain_II.quantities['xmomentum'].explicit_update,
     1699                                 eps)
     1700            assert num.allclose(valuesUsingFunction[2,i],
     1701                                domain_II.quantities['ymomentum'].explicit_update,
     1702                                eps)
     1703            i+=1
     1704
     1705        os.remove(field_sts_filename+'.sts')
     1706        os.remove(field_sts_filename+'.sww')
     1707
     1708    def test_temporally_varying_pressure_stress_from_file_evolve(self):
     1709        from anuga.config import rho_a, rho_w, eta_w
     1710        from math import pi, cos, sin
     1711        from anuga.config import time_format
     1712
     1713        cellsize = 25
     1714        #nrows=20; ncols = 20;
     1715        nrows=10; ncols = 10;
     1716        refzone=50
     1717        xllcorner=366000;yllcorner=6369500;
     1718        number_of_timesteps = 28
     1719        timestep=10.
     1720        eps=2e-16
     1721
     1722        #print "Building mesh"
     1723        #points, vertices, boundary =rectangular(10,10,
     1724        points, vertices, boundary =rectangular(5,5,
     1725                                                len1=cellsize*(ncols-1),
     1726                                                len2=cellsize*(nrows-1),
     1727                                                origin=(xllcorner,yllcorner))
     1728
     1729        domain = Domain(points, vertices, boundary)
     1730        vertexpoints = domain.get_nodes()
     1731
     1732        # Flat surface with 1m of water
     1733        domain.set_quantity('elevation', 0)
     1734        domain.set_quantity('stage', 1.0)
     1735        domain.set_quantity('friction', 0)
     1736
     1737        Br = Reflective_boundary(domain)
     1738        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1739
     1740        # Setup only one forcing term, constant wind stress
     1741        field_sts_filename = 'wind_field'
     1742        #print 'Writing pressure field sts file'
     1743        self.write_wind_pressure_field_sts(field_sts_filename,
     1744                                      nrows=nrows,
     1745                                      ncols=ncols,
     1746                                      cellsize=cellsize,
     1747                                      origin=(xllcorner,yllcorner),
     1748                                      refzone=50,
     1749                                      timestep=timestep,
     1750                                      number_of_timesteps=number_of_timesteps,
     1751                                      speed=time_varying_speed,
     1752                                      angle=time_varying_angle,
     1753                                      pressure=time_varying_pressure)
     1754
     1755        #print "converting sts to sww"
     1756        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1757                     verbose=False)
     1758
     1759        #print 'initialising file_function'
     1760        # Setup wind stress
     1761        F = file_function(field_sts_filename+'.sww', domain,
     1762                          quantities=['barometric_pressure'],
     1763                          interpolation_points = vertexpoints)
     1764
     1765        #P = Barometric_pressure(F,use_coordinates=False)
     1766        #print 'initialising pressure forcing term'
     1767        P = Barometric_pressure_fast(p=F,filename=field_sts_filename+'.sww',domain=domain)
     1768        domain.forcing_terms.append(P)
     1769
     1770        valuesUsingFunction=num.empty((3,2*number_of_timesteps,len(domain)),
     1771                                      num.float)
     1772        i=0
     1773        import time as timer
     1774        t0=timer.time()
     1775        for t in domain.evolve(yieldstep=timestep/2., finaltime=(number_of_timesteps-1)*timestep):
     1776            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
     1777            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
     1778            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
     1779            i+=1
     1780            #domain.write_time()
     1781        t1=timer.time()
     1782        #print "That took %fs seconds" %(t1-t0)
     1783
     1784
     1785        domain_II = Domain(points, vertices, boundary)
     1786
     1787        # Flat surface with 1m of water
     1788        domain_II.set_quantity('elevation', 0)
     1789        domain_II.set_quantity('stage', 1.0)
     1790        domain_II.set_quantity('friction', 0)
     1791
     1792        Br = Reflective_boundary(domain_II)
     1793        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1794
     1795        domain_II.forcing_terms.append(Barometric_pressure(p=time_varying_pressure))
     1796
     1797        i=0;
     1798        for t in domain_II.evolve(yieldstep=timestep/2.,
     1799                                  finaltime=(number_of_timesteps-1)*timestep):
     1800            assert num.allclose(valuesUsingFunction[0,i],
     1801                                domain_II.quantities['stage'].explicit_update,
     1802                                eps)
     1803            assert  num.allclose(valuesUsingFunction[1,i],
     1804                                 domain_II.quantities['xmomentum'].explicit_update,
     1805                                 eps)
     1806            assert num.allclose(valuesUsingFunction[2,i],
     1807                                domain_II.quantities['ymomentum'].explicit_update,
     1808                                eps)
     1809            i+=1
     1810
     1811        os.remove(field_sts_filename+'.sts')
     1812        os.remove(field_sts_filename+'.sww')
     1813
     1814    def test_spatially_varying_pressure_stress_from_file_evolve(self):
     1815        from anuga.config import rho_a, rho_w, eta_w
     1816        from math import pi, cos, sin
     1817        from anuga.config import time_format
     1818
     1819        cellsize = 25
     1820        #nrows=20; ncols = 20;
     1821        nrows=10; ncols = 10;
     1822        refzone=50
     1823        xllcorner=366000;yllcorner=6369500;
     1824        number_of_timesteps = 28
     1825        timestep=1.
     1826        eps=2e-16
     1827
     1828        #points, vertices, boundary =rectangular(10,10,
     1829        points, vertices, boundary =rectangular(5,5,
     1830                                                len1=cellsize*(ncols-1),
     1831                                                len2=cellsize*(nrows-1),
     1832                                                origin=(xllcorner,yllcorner))
     1833
     1834        domain = Domain(points, vertices, boundary)
     1835        vertexpoints = domain.get_nodes()
     1836
     1837        # Flat surface with 1m of water
     1838        domain.set_quantity('elevation', 0)
     1839        domain.set_quantity('stage', 1.0)
     1840        domain.set_quantity('friction', 0)
     1841
     1842        Br = Reflective_boundary(domain)
     1843        domain.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1844
     1845        # Setup only one forcing term, constant wind stress
     1846        field_sts_filename = 'wind_field'
     1847        self.write_wind_pressure_field_sts(field_sts_filename,
     1848                                      nrows=nrows,
     1849                                      ncols=ncols,
     1850                                      cellsize=cellsize,
     1851                                      origin=(xllcorner,yllcorner),
     1852                                      refzone=50,
     1853                                      timestep=timestep,
     1854                                      number_of_timesteps=number_of_timesteps,
     1855                                      speed=spatial_linear_varying_speed,
     1856                                      angle=spatial_linear_varying_angle,
     1857                                      pressure=spatial_linear_varying_pressure)
     1858
     1859        sts2sww_mesh(field_sts_filename,spatial_thinning=1,
     1860                     verbose=False)
     1861
     1862        # Setup wind stress
     1863        F = file_function(field_sts_filename+'.sww', domain,
     1864                          quantities=['barometric_pressure'],
     1865                          interpolation_points = vertexpoints)
     1866
     1867        P = Barometric_pressure(F,use_coordinates=False)
     1868        domain.forcing_terms.append(P)
     1869
     1870        valuesUsingFunction=num.empty((3,number_of_timesteps,len(domain)),
     1871                                      num.float)
     1872        i=0
     1873        for t in domain.evolve(yieldstep=timestep, finaltime=(number_of_timesteps-1)*timestep):
     1874            valuesUsingFunction[0,i]=domain.quantities['stage'].explicit_update
     1875            valuesUsingFunction[1,i]=domain.quantities['xmomentum'].explicit_update
     1876            valuesUsingFunction[2,i]=domain.quantities['ymomentum'].explicit_update
     1877            i+=1
     1878
     1879
     1880        domain_II = Domain(points, vertices, boundary)
     1881
     1882        # Flat surface with 1m of water
     1883        domain_II.set_quantity('elevation', 0)
     1884        domain_II.set_quantity('stage', 1.0)
     1885        domain_II.set_quantity('friction', 0)
     1886
     1887        Br = Reflective_boundary(domain_II)
     1888        domain_II.set_boundary({'top': Br, 'bottom' :Br, 'left': Br, 'right': Br})
     1889
     1890        domain_II.forcing_terms.append(Barometric_pressure(p=spatial_linear_varying_pressure))
     1891
     1892        i=0;
     1893        for t in domain_II.evolve(yieldstep=timestep,
     1894                                  finaltime=(number_of_timesteps-1)*timestep):
     1895
     1896            assert num.allclose(valuesUsingFunction[0,i],
     1897                                domain_II.quantities['stage'].explicit_update,
     1898                                eps)
     1899            assert  num.allclose(valuesUsingFunction[1,i],
     1900                                 domain_II.quantities['xmomentum'].explicit_update,
     1901                                 eps)
     1902            assert num.allclose(valuesUsingFunction[2,i],
     1903                                domain_II.quantities['ymomentum'].explicit_update,
     1904                                eps)
     1905            i+=1
     1906
     1907        os.remove(field_sts_filename+'.sts')
     1908        os.remove(field_sts_filename+'.sww')
     1909
     1910
    8571911           
    8581912if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.