Changeset 4733


Ignore:
Timestamp:
Sep 12, 2007, 4:13:11 PM (15 years ago)
Author:
ole
Message:

Further optimisation when beta_h == 0 (gained about 5s bring the
total time of 30s down to 25s in the okushiri performance test).

Also cleared out some obsolete code

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/config.py

    r4732 r4733  
    8181beta_h      = 0.2
    8282
    83 # beta_h can be safely put to zero esp if we are using tight_slope_limiters = 1. This will
    84 # also speed things up.
     83# beta_h can be safely put to zero esp if we are using
     84# tight_slope_limiters = 1. This will
     85# also speed things up in general
    8586beta_h = 0.0
    8687
     
    131132
    132133use_psyco = True  #Use psyco optimisations
    133 se_psyco = False  #Do not use psyco optimisations
     134#use_psyco = False  #Do not use psyco optimisations
    134135
    135136
  • anuga_core/source/anuga/shallow_water/__init__.py

    r4360 r4733  
    1212    Transmissive_Momentum_Set_Stage_boundary,\
    1313    Dirichlet_Discharge_boundary,\
    14     Constant_stage, Constant_height,\
    1514    Field_boundary
    1615
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4713 r4733  
    206206        self.smooth = not flag
    207207
    208         #Reduction operation for get_vertex_values
     208        # Reduction operation for get_vertex_values
    209209        if reduction is None:
    210210            self.reduction = mean
     
    506506
    507507
    508 #Rotation of momentum vector
     508# Rotation of momentum vector
    509509def rotate(q, normal, direction = 1):
    510510    """Rotate the momentum component q (q[1], q[2])
     
    12211221    """
    12221222
    1223     #Shortcuts
     1223    # FIXME (Ole): I reckon this can be simplified significantly:
     1224    #
     1225    # Always use beta_h == 0, and phase it out.
     1226    # Compute hc and hv in the c-code
     1227    # Omit updating xmomv
     1228    #
     1229    from shallow_water_ext import balance_deep_and_shallow
     1230   
     1231    # Shortcuts
    12241232    wc = domain.quantities['stage'].centroid_values
    12251233    zc = domain.quantities['elevation'].centroid_values
    1226     hc = wc - zc
     1234    #hc = wc - zc
    12271235
    12281236    wv = domain.quantities['stage'].vertex_values
    12291237    zv = domain.quantities['elevation'].vertex_values
    1230     hv = wv - zv
    1231 
    1232     #Momentums at centroids
     1238    #hv = wv - zv
     1239
     1240    # Momentums at centroids
    12331241    xmomc = domain.quantities['xmomentum'].centroid_values
    12341242    ymomc = domain.quantities['ymomentum'].centroid_values
    12351243
    1236     #Momentums at vertices
     1244    # Momentums at vertices
    12371245    xmomv = domain.quantities['xmomentum'].vertex_values
    12381246    ymomv = domain.quantities['ymomentum'].vertex_values
     
    12411249    if domain.beta_h > 0:
    12421250        hvbar = h_limiter(domain)
     1251       
     1252        balance_deep_and_shallow(domain, domain.beta_h,
     1253                                 wc, zc, wv, zv, hvbar,
     1254                                 xmomc, ymomc, xmomv, ymomv)       
    12431255    else:
    12441256        # print 'Using first order h-limiter'
    1245      
     1257        # FIXME: Pass wc in for now - it will be ignored.
    12461258       
    12471259        #This is how one would make a first order h_limited value
    12481260        #as in the old balancer (pre 17 Feb 2005):
    12491261        # If we wish to hard wire this, one should modify the C-code
    1250         from Numeric import zeros, Float
    1251         hvbar = zeros( (len(hc), 3), Float)
    1252         for i in range(3):
    1253             hvbar[:,i] = hc[:]
    1254 
    1255     from shallow_water_ext import balance_deep_and_shallow
    1256     balance_deep_and_shallow(domain, wc, zc, hc, wv, zv, hv, hvbar,
    1257                              xmomc, ymomc, xmomv, ymomv)
     1262        #from Numeric import zeros, Float
     1263        #hvbar = zeros( (len(wc), 3), Float)
     1264        #for i in range(3):
     1265        #    hvbar[:,i] = wc[:] - zc[:]
     1266
     1267        balance_deep_and_shallow(domain, domain.beta_h,
     1268                                 wc, zc, wv, zv, wc,
     1269                                 xmomc, ymomc, xmomv, ymomv)
     1270
    12581271
    12591272
     
    20052018
    20062019
    2007 ##############################
    2008 #OBSOLETE STUFF
    2009 
    2010 def balance_deep_and_shallow_old(domain):
    2011     """Compute linear combination between stage as computed by
    2012     gradient-limiters and stage computed as constant height above bed.
    2013     The former takes precedence when heights are large compared to the
    2014     bed slope while the latter takes precedence when heights are
    2015     relatively small.  Anything in between is computed as a balanced
    2016     linear combination in order to avoid numerical disturbances which
    2017     would otherwise appear as a result of hard switching between
    2018     modes.
    2019     """
    2020 
    2021     #OBSOLETE
    2022 
    2023     #Shortcuts
    2024     wc = domain.quantities['stage'].centroid_values
    2025     zc = domain.quantities['elevation'].centroid_values
    2026     hc = wc - zc
    2027 
    2028     wv = domain.quantities['stage'].vertex_values
    2029     zv = domain.quantities['elevation'].vertex_values
    2030     hv = wv-zv
    2031 
    2032 
    2033     #Computed linear combination between constant stages and and
    2034     #stages parallel to the bed elevation.
    2035     for k in range(len(domain)):
    2036         #Compute maximal variation in bed elevation
    2037         #  This quantitiy is
    2038         #    dz = max_i abs(z_i - z_c)
    2039         #  and it is independent of dimension
    2040         #  In the 1d case zc = (z0+z1)/2
    2041         #  In the 2d case zc = (z0+z1+z2)/3
    2042 
    2043         dz = max(abs(zv[k,0]-zc[k]),
    2044                  abs(zv[k,1]-zc[k]),
    2045                  abs(zv[k,2]-zc[k]))
    2046 
    2047 
    2048         hmin = min( hv[k,:] )
    2049 
    2050         #Create alpha in [0,1], where alpha==0 means using shallow
    2051         #first order scheme and alpha==1 means using the stage w as
    2052         #computed by the gradient limiter (1st or 2nd order)
    2053         #
    2054         #If hmin > dz/2 then alpha = 1 and the bed will have no effect
    2055         #If hmin < 0 then alpha = 0 reverting to constant height above bed.
    2056 
    2057         if dz > 0.0:
    2058             alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
    2059         else:
    2060             #Flat bed
    2061             alpha = 1.0
    2062 
    2063         #Weighted balance between stage parallel to bed elevation
    2064         #(wvi = zvi + hc) and stage as computed by 1st or 2nd
    2065         #order gradient limiter
    2066         #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
    2067         #
    2068         #It follows that the updated wvi is
    2069         #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
    2070         #  zvi + hc + alpha*(hvi - hc)
    2071         #
    2072         #Note that hvi = zc+hc-zvi in the first order case (constant).
    2073 
    2074         if alpha < 1:
    2075             for i in range(3):
    2076                 wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
    2077 
    2078 
    2079             #Momentums at centroids
    2080             xmomc = domain.quantities['xmomentum'].centroid_values
    2081             ymomc = domain.quantities['ymomentum'].centroid_values
    2082 
    2083             #Momentums at vertices
    2084             xmomv = domain.quantities['xmomentum'].vertex_values
    2085             ymomv = domain.quantities['ymomentum'].vertex_values
    2086 
    2087             # Update momentum as a linear combination of
    2088             # xmomc and ymomc (shallow) and momentum
    2089             # from extrapolator xmomv and ymomv (deep).
    2090             xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]
    2091             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
    2092 
    2093 
    2094 
    2095 
    2096 
    2097 ###########################
    2098 ###########################
    2099 #Geometries
    2100 
    2101 
    2102 #FIXME: Rethink this way of creating values.
    2103 
    2104 
    2105 class Weir:
    2106     """Set a bathymetry for weir with a hole and a downstream gutter
    2107     x,y are assumed to be in the unit square
    2108     """
    2109 
    2110     def __init__(self, stage):
    2111         self.inflow_stage = stage
    2112 
    2113     def __call__(self, x, y):
    2114         from Numeric import zeros, Float
    2115         from math import sqrt
    2116 
    2117         N = len(x)
    2118         assert N == len(y)
    2119 
    2120         z = zeros(N, Float)
    2121         for i in range(N):
    2122             z[i] = -x[i]/2  #General slope
    2123 
    2124             #Flattish bit to the left
    2125             if x[i] < 0.3:
    2126                 z[i] = -x[i]/10
    2127 
    2128             #Weir
    2129             if x[i] >= 0.3 and x[i] < 0.4:
    2130                 z[i] = -x[i]+0.9
    2131 
    2132             #Dip
    2133             x0 = 0.6
    2134             #depth = -1.3
    2135             depth = -1.0
    2136             #plateaux = -0.9
    2137             plateaux = -0.6
    2138             if y[i] < 0.7:
    2139                 if x[i] > x0 and x[i] < 0.9:
    2140                     z[i] = depth
    2141 
    2142                 #RHS plateaux
    2143                 if x[i] >= 0.9:
    2144                     z[i] = plateaux
    2145 
    2146 
    2147             elif y[i] >= 0.7 and y[i] < 1.5:
    2148                 #Restrict and deepen
    2149                 if x[i] >= x0 and x[i] < 0.8:
    2150                     z[i] = depth-(y[i]/3-0.3)
    2151                     #z[i] = depth-y[i]/5
    2152                     #z[i] = depth
    2153                 elif x[i] >= 0.8:
    2154                     #RHS plateaux
    2155                     z[i] = plateaux
    2156 
    2157             elif y[i] >= 1.5:
    2158                 if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
    2159                     #Widen up and stay at constant depth
    2160                     z[i] = depth-1.5/5
    2161                 elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
    2162                     #RHS plateaux
    2163                     z[i] = plateaux
    2164 
    2165 
    2166             #Hole in weir (slightly higher than inflow condition)
    2167             if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    2168                 z[i] = -x[i]+self.inflow_stage + 0.02
    2169 
    2170             #Channel behind weir
    2171             x0 = 0.5
    2172             if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
    2173                 z[i] = -x[i]+self.inflow_stage + 0.02
    2174 
    2175             if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
    2176                 #Flatten it out towards the end
    2177                 z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
    2178 
    2179             #Hole to the east
    2180             x0 = 1.1; y0 = 0.35
    2181             #if x[i] < -0.2 and y < 0.5:
    2182             if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    2183                 z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
    2184 
    2185             #Tiny channel draining hole
    2186             if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
    2187                 z[i] = -0.9 #North south
    2188 
    2189             if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
    2190                 z[i] = -1.0 + (x[i]-0.9)/3 #East west
    2191 
    2192 
    2193 
    2194             #Stuff not in use
    2195 
    2196             #Upward slope at inlet to the north west
    2197             #if x[i] < 0.0: # and y[i] > 0.5:
    2198             #    #z[i] = -y[i]+0.5  #-x[i]/2
    2199             #    z[i] = x[i]/4 - y[i]**2 + 0.5
    2200 
    2201             #Hole to the west
    2202             #x0 = -0.4; y0 = 0.35 # center
    2203             #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    2204             #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
    2205 
    2206 
    2207 
    2208 
    2209 
    2210         return z/2
    2211 
    2212 class Weir_simple:
    2213     """Set a bathymetry for weir with a hole and a downstream gutter
    2214     x,y are assumed to be in the unit square
    2215     """
    2216 
    2217     def __init__(self, stage):
    2218         self.inflow_stage = stage
    2219 
    2220     def __call__(self, x, y):
    2221         from Numeric import zeros, Float
    2222 
    2223         N = len(x)
    2224         assert N == len(y)
    2225 
    2226         z = zeros(N, Float)
    2227         for i in range(N):
    2228             z[i] = -x[i]  #General slope
    2229 
    2230             #Flat bit to the left
    2231             if x[i] < 0.3:
    2232                 z[i] = -x[i]/10  #General slope
    2233 
    2234             #Weir
    2235             if x[i] > 0.3 and x[i] < 0.4:
    2236                 z[i] = -x[i]+0.9
    2237 
    2238             #Dip
    2239             if x[i] > 0.6 and x[i] < 0.9:
    2240                 z[i] = -x[i]-0.5  #-y[i]/5
    2241 
    2242             #Hole in weir (slightly higher than inflow condition)
    2243             if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    2244                 z[i] = -x[i]+self.inflow_stage + 0.05
    2245 
    2246 
    2247         return z/2
    2248 
    2249 
    2250 
    2251 class Constant_stage:
    2252     """Set an initial condition with constant stage
    2253     """
    2254     def __init__(self, s):
    2255         self.s = s
    2256 
    2257     def __call__(self, x, y):
    2258         return self.s
    2259 
    2260 class Constant_height:
    2261     """Set an initial condition with constant water height, e.g
    2262     stage s = z+h
    2263     """
    2264 
    2265     def __init__(self, W, h):
    2266         self.W = W
    2267         self.h = h
    2268 
    2269     def __call__(self, x, y):
    2270         if self.W is None:
    2271             from Numeric import ones, Float
    2272             return self.h*ones(len(x), Float)
    2273         else:
    2274             return self.W(x,y) + self.h
    2275 
    2276 
    2277 
    2278 
    2279 def compute_fluxes_python(domain):
    2280     """Compute all fluxes and the timestep suitable for all volumes
    2281     in domain.
    2282 
    2283     Compute total flux for each conserved quantity using "flux_function"
    2284 
    2285     Fluxes across each edge are scaled by edgelengths and summed up
    2286     Resulting flux is then scaled by area and stored in
    2287     explicit_update for each of the three conserved quantities
    2288     stage, xmomentum and ymomentum
    2289 
    2290     The maximal allowable speed computed by the flux_function for each volume
    2291     is converted to a timestep that must not be exceeded. The minimum of
    2292     those is computed as the next overall timestep.
    2293 
    2294     Post conditions:
    2295       domain.explicit_update is reset to computed flux values
    2296       domain.timestep is set to the largest step satisfying all volumes.
    2297     """
    2298 
    2299     import sys
    2300     from Numeric import zeros, Float
    2301 
    2302     N = len(domain) # number_of_triangles
    2303 
    2304     #Shortcuts
    2305     Stage = domain.quantities['stage']
    2306     Xmom = domain.quantities['xmomentum']
    2307     Ymom = domain.quantities['ymomentum']
    2308     Bed = domain.quantities['elevation']
    2309 
    2310     #Arrays
    2311     stage = Stage.edge_values
    2312     xmom =  Xmom.edge_values
    2313     ymom =  Ymom.edge_values
    2314     bed =   Bed.edge_values
    2315 
    2316     stage_bdry = Stage.boundary_values
    2317     xmom_bdry =  Xmom.boundary_values
    2318     ymom_bdry =  Ymom.boundary_values
    2319 
    2320     flux = zeros((N,3), Float) #Work array for summing up fluxes
    2321 
    2322     #Loop
    2323     timestep = float(sys.maxint)
    2324     for k in range(N):
    2325 
    2326         for i in range(3):
    2327             #Quantities inside volume facing neighbour i
    2328             ql = [stage[k, i], xmom[k, i], ymom[k, i]]
    2329             zl = bed[k, i]
    2330 
    2331             #Quantities at neighbour on nearest face
    2332             n = domain.neighbours[k,i]
    2333             if n < 0:
    2334                 m = -n-1 #Convert negative flag to index
    2335                 qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
    2336                 zr = zl #Extend bed elevation to boundary
    2337             else:
    2338                 m = domain.neighbour_edges[k,i]
    2339                 qr = [stage[n, m], xmom[n, m], ymom[n, m]]
    2340                 zr = bed[n, m]
    2341 
    2342 
    2343             #Outward pointing normal vector
    2344             normal = domain.normals[k, 2*i:2*i+2]
    2345 
    2346             #Flux computation using provided function
    2347             edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    2348 
    2349             flux[k,:] = edgeflux
    2350 
    2351     return flux
    2352 
    2353 
    2354 
    2355 
    2356 
    2357 
    2358 
    2359 ##############################################
    2360 #Initialise module
    2361 
     2020#------------------
     2021# Initialise module
     2022#------------------
    23622023
    23632024from anuga.utilities import compile
    23642025if compile.can_use_C_extension('shallow_water_ext.c'):
    2365     #Replace python version with c implementations
     2026    # Replace python version with c implementations
    23662027
    23672028    from shallow_water_ext import rotate, assign_windfield_values
     
    23802041
    23812042
    2382 #Optimisation with psyco
     2043# Optimisation with psyco
    23832044from anuga.config import use_psyco
    23842045if use_psyco:
     
    24012062    pass
    24022063
    2403 # Profiling stuff
    2404 #import profile
    2405 #profiler = profile.Profile()
     2064
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4731 r4733  
    1818#include <stdio.h>
    1919
    20 //Shared code snippets
     20// Shared code snippets
    2121#include "util_ext.h"
    2222
     
    316316
    317317// Computational function for flux computation (using stage w=z+h)
    318 // FIXME (Ole): Is this used anywhere??
    319318int flux_function_kinetic(double *q_left, double *q_right,
    320319                  double z_left, double z_right,
     
    471470
    472471int _balance_deep_and_shallow(int N,
     472                              double beta_h,
    473473                              double* wc,
    474474                              double* zc,
    475                               double* hc,
    476475                              double* wv,
    477476                              double* zv,
    478                               double* hv,
    479                               double* hvbar,
     477                              double* hvbar, // Retire this
    480478                              double* xmomc,
    481479                              double* ymomc,
     
    487485
    488486  int k, k3, i;
    489   double dz, hmin, alpha, h_diff;
     487  double dz, hmin, alpha, h_diff, hc_k;
     488  double epsilon = 1.0e-6; // Temporary measure
     489  double hv[3]; // Depths at vertices
    490490
    491491  // Compute linear combination between w-limited stages and
     
    501501
    502502    k3 = 3*k;
    503 
     503    hc_k = wc[k] - zc[k]; // Centroid value at triangle k
    504504   
    505505    dz = 0.0;
     
    511511    }
    512512
     513    // Calculate depth at vertices
     514    hv[0] = wv[k3] -   zv[k3];
     515    hv[1] = wv[k3+1] - zv[k3+1];
     516    hv[2] = wv[k3+2] - zv[k3+2];       
     517   
    513518    // Calculate minimal depth across all three vertices
    514     hmin = min(hv[k3], min(hv[k3+1], hv[k3+2]));
    515    
     519    hmin = min(hv[0], min(hv[1], hv[2]));
    516520   
    517521
     
    544548        for (i=0; i<3; i++) {
    545549
    546           h_diff = hvbar[k3+i] - hv[k3+i];
     550          // FIXME (Ole): Simplify when (if) hvbar gets retired
     551          if (beta_h > epsilon) {
     552            h_diff = hvbar[k3+i] - hv[i];
     553          } else {
     554            h_diff = hc_k - hv[i];       
     555          }
    547556       
    548557          if (h_diff <= 0) {
     
    554563            // h-limited stage.
    555564           
    556             alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
     565            // FIXME (Ole): Simplify when (if) hvbar gets retired           
     566            if (beta_h > epsilon) {       
     567              alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
     568            } else {
     569              alpha = min(alpha, (hc_k - H0)/h_diff);       
     570            }
    557571          }
    558572        }
     
    594608    if (alpha < 1) {
    595609      for (i=0; i<3; i++) {
    596         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
     610     
     611        // FIXME (Ole): Simplify when (if) hvbar gets retired       
     612        if (beta_h > epsilon) {   
     613          wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[i];
     614        } else {
     615          wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];   
     616        }
    597617
    598618        // Update momentum as a linear combination of
     
    647667        //New code: Adjust momentum to guarantee speeds are physical
    648668        //          ensure h is non negative
    649         //FIXME (Ole): This is only implemented in this C extension and
    650         //             has no Python equivalent
    651669             
    652670        if (hc <= 0.0) {
     
    14511469
    14521470
    1453 
    1454 // FIXME (Ole): This function is obsolete as of 12 Sep 2007
    1455 PyObject *extrapolate_second_order_sw_original(PyObject *self, PyObject *args) {
    1456   /*Compute the vertex values based on a linear reconstruction on each triangle
    1457     These values are calculated as follows:
    1458     1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
    1459     formed by the centroids of its three neighbours.
    1460     2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
    1461     of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
    1462     of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
    1463     original triangle has gradient (a,b).
    1464     3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions
    1465     find_qmin_and_qmax and limit_gradient
    1466 
    1467     Python call:
    1468     extrapolate_second_order_sw(domain.surrogate_neighbours,
    1469                                 domain.number_of_boundaries
    1470                                 domain.centroid_coordinates,
    1471                                 Stage.centroid_values
    1472                                 Xmom.centroid_values
    1473                                 Ymom.centroid_values
    1474                                 domain.vertex_coordinates,
    1475                                 Stage.vertex_values,
    1476                                 Xmom.vertex_values,
    1477                                 Ymom.vertex_values)
    1478 
    1479     Post conditions:
    1480             The vertices of each triangle have values from a limited linear reconstruction
    1481             based on centroid values
    1482 
    1483   */
    1484   PyArrayObject *surrogate_neighbours,
    1485     *number_of_boundaries,
    1486     *centroid_coordinates,
    1487     *stage_centroid_values,
    1488     *xmom_centroid_values,
    1489     *ymom_centroid_values,
    1490         *elevation_centroid_values,
    1491     *vertex_coordinates,
    1492     *stage_vertex_values,
    1493     *xmom_vertex_values,
    1494     *ymom_vertex_values,
    1495         *elevation_vertex_values;
    1496   PyObject *domain, *Tmp;
    1497   //Local variables
    1498   double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
    1499   int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
    1500   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
    1501   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    1502   double dqv[3], qmin, qmax, hmin, hmax;
    1503   double hc, h0, h1, h2;
    1504   double epsilon=1.0e-12;
    1505   int optimise_dry_cells=1; // Optimisation flag   
    1506   double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
    1507   double minimum_allowed_height;
    1508  
    1509   // Provisional jumps from centroids to v'tices and safety factor re limiting
    1510   // by which these jumps are limited
    1511   // Convert Python arguments to C
    1512   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    1513                         &domain,
    1514                         &surrogate_neighbours,
    1515                         &number_of_boundaries,
    1516                         &centroid_coordinates,
    1517                         &stage_centroid_values,
    1518                         &xmom_centroid_values,
    1519                         &ymom_centroid_values,
    1520                         &elevation_centroid_values,
    1521                         &vertex_coordinates,
    1522                         &stage_vertex_values,
    1523                         &xmom_vertex_values,
    1524                         &ymom_vertex_values,
    1525                         &elevation_vertex_values,
    1526                         &optimise_dry_cells)) {                 
    1527                        
    1528     PyErr_SetString(PyExc_RuntimeError, "Input arguments to extrapolate_second_order_sw failed");
    1529     return NULL;
    1530   }
    1531 
    1532  
    1533   // Get the safety factor beta_w, set in the config.py file. This is used in the limiting process
    1534   Tmp = PyObject_GetAttrString(domain, "beta_w");
    1535   if (!Tmp) {
    1536     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
    1537     return NULL;
    1538   } 
    1539   beta_w = PyFloat_AsDouble(Tmp);
    1540   Py_DECREF(Tmp);
    1541  
    1542   Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
    1543   if (!Tmp) {
    1544     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
    1545     return NULL;
    1546   } 
    1547   beta_w_dry = PyFloat_AsDouble(Tmp);
    1548   Py_DECREF(Tmp);
    1549  
    1550   Tmp = PyObject_GetAttrString(domain, "beta_uh");
    1551   if (!Tmp) {
    1552     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
    1553     return NULL;
    1554   } 
    1555   beta_uh = PyFloat_AsDouble(Tmp);
    1556   Py_DECREF(Tmp);
    1557  
    1558   Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
    1559   if (!Tmp) {
    1560     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
    1561     return NULL;
    1562   } 
    1563   beta_uh_dry = PyFloat_AsDouble(Tmp);
    1564   Py_DECREF(Tmp);
    1565 
    1566   Tmp = PyObject_GetAttrString(domain, "beta_vh");
    1567   if (!Tmp) {
    1568     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
    1569     return NULL;
    1570   } 
    1571   beta_vh = PyFloat_AsDouble(Tmp);
    1572   Py_DECREF(Tmp);
    1573  
    1574   Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
    1575   if (!Tmp) {
    1576     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
    1577     return NULL;
    1578   } 
    1579   beta_vh_dry = PyFloat_AsDouble(Tmp);
    1580   Py_DECREF(Tmp);
    1581  
    1582   Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
    1583   if (!Tmp) {
    1584     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
    1585     return NULL;
    1586   } 
    1587   minimum_allowed_height = PyFloat_AsDouble(Tmp);
    1588   Py_DECREF(Tmp); 
    1589 
    1590   Tmp = PyObject_GetAttrString(domain, "epsilon");
    1591   if (!Tmp) {
    1592     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
    1593     return NULL;
    1594   } 
    1595   epsilon = PyFloat_AsDouble(Tmp);
    1596   Py_DECREF(Tmp); 
    1597    
    1598   number_of_elements = stage_centroid_values -> dimensions[0];
    1599   for (k=0; k<number_of_elements; k++) {
    1600     k3=k*3;
    1601     k6=k*6;
    1602 
    1603    
    1604     if (((long *) number_of_boundaries->data)[k]==3){
    1605       // No neighbours, set gradient on the triangle to zero*/
    1606       ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
    1607       ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
    1608       ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
    1609       ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
    1610       ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
    1611       ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
    1612       ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
    1613       ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
    1614       ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
    1615       continue;
    1616     }
    1617     else {
    1618       // Triangle k has one or more neighbours.
    1619       // Get centroid and vertex coordinates of the triangle
    1620      
    1621       // Get the vertex coordinates
    1622       xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
    1623       xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
    1624       xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
    1625      
    1626       // Get the centroid coordinates
    1627       coord_index=2*k;
    1628       x=((double *)centroid_coordinates->data)[coord_index];
    1629       y=((double *)centroid_coordinates->data)[coord_index+1];
    1630      
    1631       // Store x- and y- differentials for the vertices of the FV triangle relative to the centroid
    1632       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    1633       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
    1634     }
    1635 
    1636 
    1637    
    1638    
    1639    
    1640            
    1641     if (((long *)number_of_boundaries->data)[k]<=1){
    1642    
    1643       //==============================================
    1644       // Number of boundaries <= 1
    1645       //==============================================   
    1646    
    1647    
    1648       // If no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
    1649       // If one boundary, auxiliary triangle is formed from this centroid and its two neighbours
    1650       k0=((long *)surrogate_neighbours->data)[k3];
    1651       k1=((long *)surrogate_neighbours->data)[k3+1];
    1652       k2=((long *)surrogate_neighbours->data)[k3+2];
    1653      
    1654       // Get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
    1655       coord_index=2*k0;
    1656       x0=((double *)centroid_coordinates->data)[coord_index];
    1657       y0=((double *)centroid_coordinates->data)[coord_index+1];
    1658       coord_index=2*k1;
    1659       x1=((double *)centroid_coordinates->data)[coord_index];
    1660       y1=((double *)centroid_coordinates->data)[coord_index+1];
    1661       coord_index=2*k2;
    1662       x2=((double *)centroid_coordinates->data)[coord_index];
    1663       y2=((double *)centroid_coordinates->data)[coord_index+1];
    1664      
    1665       // Store x- and y- differentials for the vertices of the auxiliary triangle
    1666       dx1=x1-x0; dx2=x2-x0;
    1667       dy1=y1-y0; dy2=y2-y0;
    1668      
    1669       // Calculate 2*area of the auxiliary triangle
    1670       area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
    1671      
    1672       // If the mesh is 'weird' near the boundary, the triangle might be flat or clockwise:
    1673       if (area2<=0) {
    1674         PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: negative triangle area encountered");
    1675         return NULL;
    1676       } 
    1677      
    1678       // Calculate heights of neighbouring cells
    1679       hc = ((double *)stage_centroid_values->data)[k]  - ((double *)elevation_centroid_values->data)[k];
    1680       h0 = ((double *)stage_centroid_values->data)[k0] - ((double *)elevation_centroid_values->data)[k0];
    1681       h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
    1682       h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
    1683       hmin = min(hc,min(h0,min(h1,h2)));  // FIXME Don't need to include hc
    1684      
    1685      
    1686       if (optimise_dry_cells) {     
    1687         // Check if linear reconstruction is necessary for triangle k
    1688         // This check will exclude dry cells.
    1689 
    1690         hmax = max(h0,max(h1,h2));     
    1691         if (hmax < epsilon) {
    1692           continue;
    1693         }
    1694       }
    1695 
    1696            
    1697       //-----------------------------------
    1698       // stage
    1699       //-----------------------------------     
    1700      
    1701       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1702       dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
    1703      
    1704       // Calculate differentials between the vertices of the auxiliary triangle
    1705       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
    1706       dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
    1707      
    1708       // Calculate the gradient of stage on the auxiliary triangle
    1709       a = dy2*dq1 - dy1*dq2;
    1710       a /= area2;
    1711       b = dx1*dq2 - dx2*dq1;
    1712       b /= area2;
    1713      
    1714       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1715       dqv[0]=a*dxv0+b*dyv0;
    1716       dqv[1]=a*dxv1+b*dyv1;
    1717       dqv[2]=a*dxv2+b*dyv2;
    1718      
    1719       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1720       // and compute jumps from the centroid to the min and max
    1721       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1722      
    1723       // Playing with dry wet interface
    1724       hmin = qmin;
    1725       beta_tmp = beta_w;
    1726       if (hmin<minimum_allowed_height)
    1727         beta_tmp = beta_w_dry;
    1728       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1729       for (i=0;i<3;i++)
    1730         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    1731      
    1732      
    1733       //-----------------------------------
    1734       // xmomentum
    1735       //-----------------------------------           
    1736 
    1737       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1738       dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
    1739      
    1740       // Calculate differentials between the vertices of the auxiliary triangle
    1741       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
    1742       dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
    1743      
    1744       // Calculate the gradient of xmom on the auxiliary triangle
    1745       a = dy2*dq1 - dy1*dq2;
    1746       a /= area2;
    1747       b = dx1*dq2 - dx2*dq1;
    1748       b /= area2;
    1749      
    1750       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1751       dqv[0]=a*dxv0+b*dyv0;
    1752       dqv[1]=a*dxv1+b*dyv1;
    1753       dqv[2]=a*dxv2+b*dyv2;
    1754      
    1755       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1756       // and compute jumps from the centroid to the min and max
    1757       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1758       beta_tmp = beta_uh;
    1759       if (hmin<minimum_allowed_height)
    1760         beta_tmp = beta_uh_dry;
    1761       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1762       for (i=0;i<3;i++)
    1763         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    1764      
    1765      
    1766       //-----------------------------------
    1767       // ymomentum
    1768       //-----------------------------------                 
    1769 
    1770       // Calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    1771       dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
    1772      
    1773       // Calculate differentials between the vertices of the auxiliary triangle
    1774       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
    1775       dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
    1776      
    1777       // Calculate the gradient of xmom on the auxiliary triangle
    1778       a = dy2*dq1 - dy1*dq2;
    1779       a /= area2;
    1780       b = dx1*dq2 - dx2*dq1;
    1781       b /= area2;
    1782      
    1783       // Calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    1784       dqv[0]=a*dxv0+b*dyv0;
    1785       dqv[1]=a*dxv1+b*dyv1;
    1786       dqv[2]=a*dxv2+b*dyv2;
    1787      
    1788       // Now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    1789       // and compute jumps from the centroid to the min and max
    1790       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    1791       beta_tmp = beta_vh;
    1792       if (hmin<minimum_allowed_height)
    1793         beta_tmp = beta_vh_dry;
    1794       limit_gradient(dqv,qmin,qmax,beta_tmp);//the gradient will be limited
    1795       for (i=0;i<3;i++)
    1796         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    1797     } // End number_of_boundaries <=1
    1798     else{
    1799 
    1800       //==============================================
    1801       // Number of boundaries == 2
    1802       //==============================================       
    1803        
    1804       // One internal neighbour and gradient is in direction of the neighbour's centroid
    1805      
    1806       // Find the only internal neighbour
    1807       for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
    1808         if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
    1809           break;
    1810       }
    1811       if ((k2==k3+3)) {//if we didn't find an internal neighbour
    1812         PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: Internal neighbour not found");     
    1813         return NULL;//error
    1814       }
    1815      
    1816       k1=((long *)surrogate_neighbours->data)[k2];
    1817      
    1818       // The coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
    1819       coord_index=2*k1;
    1820       x1=((double *)centroid_coordinates->data)[coord_index];
    1821       y1=((double *)centroid_coordinates->data)[coord_index+1];
    1822      
    1823       // Compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
    1824       dx1=x1-x; dy1=y1-y;
    1825      
    1826       // Set area2 as the square of the distance
    1827       area2=dx1*dx1+dy1*dy1;
    1828      
    1829       // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
    1830       // respectively correspond to the x- and y- gradients of the conserved quantities
    1831       dx2=1.0/area2;
    1832       dy2=dx2*dy1;
    1833       dx2*=dx1;
    1834      
    1835      
    1836      
    1837       //-----------------------------------
    1838       // stage
    1839       //-----------------------------------           
    1840 
    1841       // Compute differentials
    1842       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
    1843      
    1844       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1845       a=dq1*dx2;
    1846       b=dq1*dy2;
    1847      
    1848       // Calculate provisional vertex jumps, to be limited
    1849       dqv[0]=a*dxv0+b*dyv0;
    1850       dqv[1]=a*dxv1+b*dyv1;
    1851       dqv[2]=a*dxv2+b*dyv2;
    1852      
    1853       // Now limit the jumps
    1854       if (dq1>=0.0){
    1855         qmin=0.0;
    1856         qmax=dq1;
    1857       }
    1858       else{
    1859         qmin=dq1;
    1860         qmax=0.0;
    1861       }
    1862      
    1863      
    1864       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1865       for (i=0;i<3;i++)
    1866         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    1867      
    1868       //-----------------------------------
    1869       // xmomentum
    1870       //-----------------------------------                       
    1871      
    1872       // Compute differentials
    1873       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
    1874      
    1875       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1876       a=dq1*dx2;
    1877       b=dq1*dy2;
    1878      
    1879       // Calculate provisional vertex jumps, to be limited
    1880       dqv[0]=a*dxv0+b*dyv0;
    1881       dqv[1]=a*dxv1+b*dyv1;
    1882       dqv[2]=a*dxv2+b*dyv2;
    1883      
    1884       // Now limit the jumps
    1885       if (dq1>=0.0){
    1886         qmin=0.0;
    1887         qmax=dq1;
    1888       }
    1889       else{
    1890         qmin=dq1;
    1891         qmax=0.0;
    1892       }
    1893       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1894       for (i=0;i<3;i++)
    1895         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    1896      
    1897       //-----------------------------------
    1898       // ymomentum
    1899       //-----------------------------------                       
    1900 
    1901       // Compute differentials
    1902       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
    1903      
    1904       // Calculate the gradient between the centroid of the FV triangle and that of its neighbour
    1905       a=dq1*dx2;
    1906       b=dq1*dy2;
    1907      
    1908       // Calculate provisional vertex jumps, to be limited
    1909       dqv[0]=a*dxv0+b*dyv0;
    1910       dqv[1]=a*dxv1+b*dyv1;
    1911       dqv[2]=a*dxv2+b*dyv2;
    1912      
    1913       // Now limit the jumps
    1914       if (dq1>=0.0){
    1915         qmin=0.0;
    1916         qmax=dq1;
    1917       }
    1918       else{
    1919         qmin=dq1;
    1920         qmax=0.0;
    1921       }
    1922       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    1923       for (i=0;i<3;i++)
    1924         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    1925     }//else [number_of_boundaries==2]
    1926   }//for k=0 to number_of_elements-1
    1927  
    1928   return Py_BuildValue("");
    1929 }//extrapolate_second-order_sw
    19301471
    19311472
     
    22951836}
    22961837
    2297 
    2298 
    2299 
    2300 // THIS FUNCTION IS NOW OBSOLETE
    2301 PyObject *compute_fluxes_ext_central_original(PyObject *self, PyObject *args) {
    2302   /*Compute all fluxes and the timestep suitable for all volumes
    2303     in domain.
    2304 
    2305     Compute total flux for each conserved quantity using "flux_function_central"
    2306 
    2307     Fluxes across each edge are scaled by edgelengths and summed up
    2308     Resulting flux is then scaled by area and stored in
    2309     explicit_update for each of the three conserved quantities
    2310     stage, xmomentum and ymomentum
    2311 
    2312     The maximal allowable speed computed by the flux_function for each volume
    2313     is converted to a timestep that must not be exceeded. The minimum of
    2314     those is computed as the next overall timestep.
    2315 
    2316     Python call:
    2317     domain.timestep = compute_fluxes(timestep,
    2318                                      domain.epsilon,
    2319                                      domain.H0,
    2320                                      domain.g,
    2321                                      domain.neighbours,
    2322                                      domain.neighbour_edges,
    2323                                      domain.normals,
    2324                                      domain.edgelengths,
    2325                                      domain.radii,
    2326                                      domain.areas,
    2327                                      tri_full_flag,
    2328                                      Stage.edge_values,
    2329                                      Xmom.edge_values,
    2330                                      Ymom.edge_values,
    2331                                      Bed.edge_values,
    2332                                      Stage.boundary_values,
    2333                                      Xmom.boundary_values,
    2334                                      Ymom.boundary_values,
    2335                                      Stage.explicit_update,
    2336                                      Xmom.explicit_update,
    2337                                      Ymom.explicit_update,
    2338                                      already_computed_flux,
    2339                                      optimise_dry_cells)                                     
    2340 
    2341 
    2342     Post conditions:
    2343       domain.explicit_update is reset to computed flux values
    2344       domain.timestep is set to the largest step satisfying all volumes.
    2345 
    2346 
    2347   */
    2348 
    2349 
    2350   PyArrayObject *neighbours, *neighbour_edges,
    2351     *normals, *edgelengths, *radii, *areas,
    2352     *tri_full_flag,
    2353     *stage_edge_values,
    2354     *xmom_edge_values,
    2355     *ymom_edge_values,
    2356     *bed_edge_values,
    2357     *stage_boundary_values,
    2358     *xmom_boundary_values,
    2359     *ymom_boundary_values,
    2360     *stage_explicit_update,
    2361     *xmom_explicit_update,
    2362     *ymom_explicit_update,
    2363     *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    2364     *max_speed_array; //Keeps track of max speeds for each triangle
    2365 
    2366 
    2367   // Local variables
    2368   double timestep, max_speed, epsilon, g, H0, length, area;
    2369   int optimise_dry_cells=0; // Optimisation flag 
    2370   double normal[2], ql[3], qr[3], zl, zr;
    2371   double edgeflux[3]; // Work array for summing up fluxes
    2372 
    2373   int number_of_elements, k, i, m, n;
    2374 
    2375   int ki, nm=0, ki2; // Index shorthands
    2376   static long call=1; // Static local variable flagging already computed flux
    2377 
    2378 
    2379   // Convert Python arguments to C
    2380   if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    2381                         &timestep,
    2382                         &epsilon,
    2383                         &H0,
    2384                         &g,
    2385                         &neighbours,
    2386                         &neighbour_edges,
    2387                         &normals,
    2388                         &edgelengths, &radii, &areas,
    2389                         &tri_full_flag,
    2390                         &stage_edge_values,
    2391                         &xmom_edge_values,
    2392                         &ymom_edge_values,
    2393                         &bed_edge_values,
    2394                         &stage_boundary_values,
    2395                         &xmom_boundary_values,
    2396                         &ymom_boundary_values,
    2397                         &stage_explicit_update,
    2398                         &xmom_explicit_update,
    2399                         &ymom_explicit_update,
    2400                         &already_computed_flux,
    2401                         &max_speed_array,
    2402                         &optimise_dry_cells)) {
    2403     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    2404     return NULL;
    2405   }
    2406  
    2407  
    2408   number_of_elements = stage_edge_values -> dimensions[0];
    2409 
    2410   call++; // Flag 'id' of flux calculation for this timestep
    2411  
    2412   // Set explicit_update to zero for all conserved_quantities.
    2413   // This assumes compute_fluxes called before forcing terms
    2414   for (k=0; k<number_of_elements; k++) {
    2415     ((double *) stage_explicit_update -> data)[k]=0.0;
    2416     ((double *) xmom_explicit_update -> data)[k]=0.0;
    2417     ((double *) ymom_explicit_update -> data)[k]=0.0; 
    2418   }
    2419  
    2420   // For all triangles
    2421   for (k=0; k<number_of_elements; k++) {
    2422  
    2423     // Loop through neighbours and compute edge flux for each 
    2424     for (i=0; i<3; i++) {
    2425       ki = k*3+i; // Linear index (triangle k, edge i)
    2426      
    2427       if (((long *) already_computed_flux->data)[ki] == call)
    2428         // We've already computed the flux across this edge
    2429         continue;
    2430        
    2431        
    2432       ql[0] = ((double *) stage_edge_values -> data)[ki];
    2433       ql[1] = ((double *) xmom_edge_values -> data)[ki];
    2434       ql[2] = ((double *) ymom_edge_values -> data)[ki];
    2435       zl =    ((double *) bed_edge_values -> data)[ki];
    2436 
    2437       // Quantities at neighbour on nearest face
    2438       n = ((long *) neighbours -> data)[ki];
    2439       if (n < 0) {
    2440         m = -n-1; // Convert negative flag to index
    2441        
    2442         qr[0] = ((double *) stage_boundary_values -> data)[m];
    2443         qr[1] = ((double *) xmom_boundary_values -> data)[m];
    2444         qr[2] = ((double *) ymom_boundary_values -> data)[m];
    2445         zr = zl; //Extend bed elevation to boundary
    2446       } else {
    2447         m = ((long *) neighbour_edges -> data)[ki];
    2448         nm = n*3+m; // Linear index (triangle n, edge m)
    2449        
    2450         qr[0] = ((double *) stage_edge_values -> data)[nm];
    2451         qr[1] = ((double *) xmom_edge_values -> data)[nm];
    2452         qr[2] = ((double *) ymom_edge_values -> data)[nm];
    2453         zr =    ((double *) bed_edge_values -> data)[nm];
    2454       }
    2455      
    2456      
    2457       if (optimise_dry_cells) {     
    2458         // Check if flux calculation is necessary across this edge
    2459         // This check will exclude dry cells.
    2460         // This will also optimise cases where zl != zr as
    2461         // long as both are dry
    2462 
    2463         if ( fabs(ql[0] - zl) < epsilon &&
    2464              fabs(qr[0] - zr) < epsilon ) {
    2465           // Cell boundary is dry
    2466          
    2467           ((long *) already_computed_flux -> data)[ki] = call; // #k Done       
    2468           if (n>=0)
    2469             ((long *) already_computed_flux -> data)[nm] = call; // #n Done
    2470        
    2471           max_speed = 0.0;
    2472           continue;
    2473         }
    2474       }
    2475      
    2476            
    2477       // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    2478       ki2 = 2*ki; //k*6 + i*2
    2479       normal[0] = ((double *) normals -> data)[ki2];
    2480       normal[1] = ((double *) normals -> data)[ki2+1];
    2481      
    2482 
    2483       // Edge flux computation (triangle k, edge i)
    2484       flux_function_central(ql, qr, zl, zr,
    2485                             normal[0], normal[1],
    2486                             epsilon, H0, g,
    2487                             edgeflux, &max_speed);
    2488      
    2489      
    2490       // Multiply edgeflux by edgelength
    2491       length = ((double *) edgelengths -> data)[ki];
    2492       edgeflux[0] *= length;           
    2493       edgeflux[1] *= length;           
    2494       edgeflux[2] *= length;                       
    2495      
    2496      
    2497       // Update triangle k with flux from edge i
    2498       ((double *) stage_explicit_update -> data)[k] -= edgeflux[0];
    2499       ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1];
    2500       ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2];
    2501      
    2502       ((long *) already_computed_flux -> data)[ki] = call; // #k Done
    2503      
    2504      
    2505       // Update neighbour n with same flux but reversed sign
    2506       if (n>=0){
    2507         ((double *) stage_explicit_update -> data)[n] += edgeflux[0];
    2508         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1];
    2509         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2];
    2510        
    2511         ((long *) already_computed_flux -> data)[nm] = call; // #n Done
    2512       }
    2513      
    2514      
    2515       // Update timestep based on edge i and possibly neighbour n
    2516       if ( ((long *) tri_full_flag -> data)[k] == 1) {
    2517         if (max_speed > epsilon) {
    2518           timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    2519           if (n>=0)
    2520             timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    2521         }
    2522       }
    2523      
    2524     } // End edge i
    2525    
    2526    
    2527     // Normalise triangle k by area and store for when all conserved
    2528     // quantities get updated
    2529     area = ((double *) areas -> data)[k];
    2530     ((double *) stage_explicit_update -> data)[k] /= area;
    2531     ((double *) xmom_explicit_update -> data)[k] /= area;
    2532     ((double *) ymom_explicit_update -> data)[k] /= area;
    2533    
    2534    
    2535     // Keep track of maximal speeds
    2536     ((double *) max_speed_array -> data)[k] = max_speed;   
    2537    
    2538   } // End triangle k
    2539  
    2540   return Py_BuildValue("d", timestep);
    2541 }
    25421838
    25431839
     
    27652061PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
    27662062  //
    2767   //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
     2063  //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
    27682064  //                             xmomc, ymomc, xmomv, ymomv)
    27692065
     
    27722068    *wc,            //Stage at centroids
    27732069    *zc,            //Elevation at centroids
    2774     *hc,            //Height at centroids
    27752070    *wv,            //Stage at vertices
    27762071    *zv,            //Elevation at vertices
    2777     *hv,            //Depths at vertices
    27782072    *hvbar,         //h-Limited depths at vertices
    27792073    *xmomc,         //Momentums at centroids and vertices
     
    27852079   
    27862080  double alpha_balance = 2.0;
    2787   double H0;
     2081  double H0, beta_h;
    27882082
    27892083  int N, tight_slope_limiters; //, err;
    27902084
    27912085  // Convert Python arguments to C
    2792   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOO",
     2086  if (!PyArg_ParseTuple(args, "OdOOOOOOOOO",
    27932087                        &domain,
    2794                         &wc, &zc, &hc,
    2795                         &wv, &zv, &hv, &hvbar,
     2088                        &beta_h,
     2089                        &wc, &zc,
     2090                        &wv, &zv, &hvbar,
    27962091                        &xmomc, &ymomc, &xmomv, &ymomv)) {
    2797     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     2092    PyErr_SetString(PyExc_RuntimeError,
     2093                    "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
    27982094    return NULL;
    27992095  } 
     
    28322128 
    28332129 
    2834      
    2835   //alpha_balance = 2.0;
    2836   //H0 = 0.001;
    2837   //tight_slope_limiters = 1;
    2838  
    28392130  N = wc -> dimensions[0];
    2840 
    28412131  _balance_deep_and_shallow(N,
     2132                            beta_h,
    28422133                            (double*) wc -> data,
    28432134                            (double*) zc -> data,
    2844                             (double*) hc -> data,
    28452135                            (double*) wv -> data,
    28462136                            (double*) zv -> data,
    2847                             (double*) hv -> data,
    28482137                            (double*) hvbar -> data,
    28492138                            (double*) xmomc -> data,
     
    30842373  Py_InitModule("shallow_water_ext", MethodTable);
    30852374
    3086   import_array();     //Necessary for handling of NumPY structures
    3087 }
     2375  import_array(); // Necessary for handling of NumPY structures
     2376}
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4731 r4733  
    1111from shallow_water_domain import *
    1212from shallow_water_domain import flux_function_central as flux_function
     13
     14class Weir:
     15    """Set a bathymetry for weir with a hole and a downstream gutter
     16    x,y are assumed to be in the unit square
     17    """
     18
     19    def __init__(self, stage):
     20        self.inflow_stage = stage
     21
     22    def __call__(self, x, y):
     23        from Numeric import zeros, Float
     24        from math import sqrt
     25
     26        N = len(x)
     27        assert N == len(y)
     28
     29        z = zeros(N, Float)
     30        for i in range(N):
     31            z[i] = -x[i]/2  #General slope
     32
     33            #Flattish bit to the left
     34            if x[i] < 0.3:
     35                z[i] = -x[i]/10
     36
     37            #Weir
     38            if x[i] >= 0.3 and x[i] < 0.4:
     39                z[i] = -x[i]+0.9
     40
     41            #Dip
     42            x0 = 0.6
     43            #depth = -1.3
     44            depth = -1.0
     45            #plateaux = -0.9
     46            plateaux = -0.6
     47            if y[i] < 0.7:
     48                if x[i] > x0 and x[i] < 0.9:
     49                    z[i] = depth
     50
     51                #RHS plateaux
     52                if x[i] >= 0.9:
     53                    z[i] = plateaux
     54
     55
     56            elif y[i] >= 0.7 and y[i] < 1.5:
     57                #Restrict and deepen
     58                if x[i] >= x0 and x[i] < 0.8:
     59                    z[i] = depth-(y[i]/3-0.3)
     60                    #z[i] = depth-y[i]/5
     61                    #z[i] = depth
     62                elif x[i] >= 0.8:
     63                    #RHS plateaux
     64                    z[i] = plateaux
     65
     66            elif y[i] >= 1.5:
     67                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
     68                    #Widen up and stay at constant depth
     69                    z[i] = depth-1.5/5
     70                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
     71                    #RHS plateaux
     72                    z[i] = plateaux
     73
     74
     75            #Hole in weir (slightly higher than inflow condition)
     76            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
     77                z[i] = -x[i]+self.inflow_stage + 0.02
     78
     79            #Channel behind weir
     80            x0 = 0.5
     81            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
     82                z[i] = -x[i]+self.inflow_stage + 0.02
     83
     84            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
     85                #Flatten it out towards the end
     86                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
     87
     88            #Hole to the east
     89            x0 = 1.1; y0 = 0.35
     90            #if x[i] < -0.2 and y < 0.5:
     91            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
     92                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
     93
     94            #Tiny channel draining hole
     95            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
     96                z[i] = -0.9 #North south
     97
     98            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
     99                z[i] = -1.0 + (x[i]-0.9)/3 #East west
     100
     101
     102
     103            #Stuff not in use
     104
     105            #Upward slope at inlet to the north west
     106            #if x[i] < 0.0: # and y[i] > 0.5:
     107            #    #z[i] = -y[i]+0.5  #-x[i]/2
     108            #    z[i] = x[i]/4 - y[i]**2 + 0.5
     109
     110            #Hole to the west
     111            #x0 = -0.4; y0 = 0.35 # center
     112            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
     113            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
     114
     115
     116
     117
     118
     119        return z/2
     120
     121class Weir_simple:
     122    """Set a bathymetry for weir with a hole and a downstream gutter
     123    x,y are assumed to be in the unit square
     124    """
     125
     126    def __init__(self, stage):
     127        self.inflow_stage = stage
     128
     129    def __call__(self, x, y):
     130        from Numeric import zeros, Float
     131
     132        N = len(x)
     133        assert N == len(y)
     134
     135        z = zeros(N, Float)
     136        for i in range(N):
     137            z[i] = -x[i]  #General slope
     138
     139            #Flat bit to the left
     140            if x[i] < 0.3:
     141                z[i] = -x[i]/10  #General slope
     142
     143            #Weir
     144            if x[i] > 0.3 and x[i] < 0.4:
     145                z[i] = -x[i]+0.9
     146
     147            #Dip
     148            if x[i] > 0.6 and x[i] < 0.9:
     149                z[i] = -x[i]-0.5  #-y[i]/5
     150
     151            #Hole in weir (slightly higher than inflow condition)
     152            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
     153                z[i] = -x[i]+self.inflow_stage + 0.05
     154
     155
     156        return z/2
     157
     158
     159
    13160
    14161#Variable windfield implemented using functions
     
    32663413
    32673414        #Initial condition
    3268         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3415        #domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3416        domain.set_quantity('stage', expression='elevation+0.05')
    32693417        domain.check_integrity()
    32703418
     
    33103458
    33113459        #Initial condition
    3312         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3460        domain.set_quantity('stage', expression='elevation+0.05')
    33133461        domain.check_integrity()
    33143462
     
    33823530
    33833531        #Initial condition
    3384         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3532        domain.set_quantity('stage', expression='elevation+0.05')
    33853533        domain.check_integrity()
    33863534
     
    34793627
    34803628        #Initial condition
    3481         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3629        domain.set_quantity('stage', expression='elevation+0.05')
    34823630        domain.check_integrity()
    34833631
     
    35813729
    35823730        #Initial condition
    3583         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     3731        domain.set_quantity('stage', expression='elevation+0.05')
    35843732        domain.check_integrity()
    35853733
     
    39424090
    39434091        #Initial condition
    3944         domain.set_quantity('stage', Constant_height(x_slope, 0.05))
     4092        domain.set_quantity('stage', expression='elevation+0.05')
    39454093        domain.check_integrity()
    39464094
     
    39844132        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
    39854133
    3986         domain.set_quantity('stage', Constant_height(Z, 0.))
     4134        domain.set_quantity('stage', expression='elevation')
    39874135
    39884136        for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2):
Note: See TracChangeset for help on using the changeset viewer.