[7839] | 1 | """Class Domain - |
---|
| 2 | 1D interval domains for finite-volume computations of |
---|
| 3 | the shallow water wave equation. |
---|
| 4 | |
---|
| 5 | This module contains a specialisation of class Domain from module domain.py |
---|
| 6 | consisting of methods specific to the Shallow Water Wave Equation |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | U_t + E_x = S |
---|
| 10 | |
---|
| 11 | where |
---|
| 12 | |
---|
| 13 | U = [w, uh] |
---|
| 14 | E = [uh, u^2h + gh^2/2] |
---|
| 15 | S represents source terms forcing the system |
---|
| 16 | (e.g. gravity, friction, wind stress, ...) |
---|
| 17 | |
---|
| 18 | and _t, _x, _y denote the derivative with respect to t, x and y respectiely. |
---|
| 19 | |
---|
| 20 | The quantities are |
---|
| 21 | |
---|
| 22 | symbol variable name explanation |
---|
| 23 | x x horizontal distance from origin [m] |
---|
| 24 | z elevation elevation of bed on which flow is modelled [m] |
---|
| 25 | h height water height above z [m] |
---|
| 26 | w stage absolute water level, w = z+h [m] |
---|
| 27 | u speed in the x direction [m/s] |
---|
| 28 | uh xmomentum momentum in the x direction [m^2/s] |
---|
| 29 | |
---|
| 30 | eta mannings friction coefficient [to appear] |
---|
| 31 | nu wind stress coefficient [to appear] |
---|
| 32 | |
---|
| 33 | The conserved quantities are w, uh |
---|
| 34 | |
---|
| 35 | For details see e.g. |
---|
| 36 | Christopher Zoppou and Stephen Roberts, |
---|
| 37 | Catastrophic Collapse of Water Supply Reservoirs in Urban Areas, |
---|
| 38 | Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999 |
---|
| 39 | |
---|
| 40 | |
---|
| 41 | John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
| 42 | Geoscience Australia, 2006 |
---|
| 43 | """ |
---|
| 44 | |
---|
| 45 | |
---|
[7840] | 46 | from anuga_1d.generic.generic_domain import * |
---|
[7839] | 47 | from sww_boundary_conditions import * |
---|
| 48 | from sww_forcing_terms import * |
---|
| 49 | |
---|
| 50 | #Shallow water domain |
---|
| 51 | class Domain(Generic_domain): |
---|
| 52 | |
---|
| 53 | def __init__(self, coordinates, boundary = None, tagged_elements = None): |
---|
| 54 | |
---|
| 55 | conserved_quantities = ['stage', 'xmomentum'] |
---|
| 56 | evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
| 57 | other_quantities = ['friction'] |
---|
| 58 | Generic_domain.__init__(self, |
---|
| 59 | coordinates = coordinates, |
---|
| 60 | boundary = boundary, |
---|
| 61 | conserved_quantities = conserved_quantities, |
---|
| 62 | evolved_quantities = evolved_quantities, |
---|
| 63 | other_quantities = other_quantities, |
---|
| 64 | tagged_elements = tagged_elements) |
---|
| 65 | |
---|
[7840] | 66 | from anuga_1d.config import minimum_allowed_height, g, h0 |
---|
[7839] | 67 | self.minimum_allowed_height = minimum_allowed_height |
---|
| 68 | self.g = g |
---|
| 69 | self.h0 = h0 |
---|
| 70 | |
---|
| 71 | #forcing terms not included in 1d domain ?WHy? |
---|
| 72 | self.forcing_terms.append(gravity) |
---|
| 73 | #self.forcing_terms.append(manning_friction) |
---|
| 74 | #print "\nI have Removed forcing terms line 64 1dsw" |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | #Stored output |
---|
| 78 | self.store = True |
---|
| 79 | self.format = 'sww' |
---|
| 80 | self.smooth = True |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | #Reduction operation for get_vertex_values |
---|
[7840] | 84 | from anuga_1d.generic.util import mean |
---|
[7839] | 85 | self.reduction = mean |
---|
| 86 | #self.reduction = min #Looks better near steep slopes |
---|
| 87 | |
---|
| 88 | self.set_quantities_to_be_stored(['stage','xmomentum']) |
---|
| 89 | |
---|
| 90 | self.__doc__ = 'sww_vel_domain' |
---|
| 91 | |
---|
| 92 | self.check_integrity() |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | def check_integrity(self): |
---|
| 97 | |
---|
| 98 | #Check that we are solving the shallow water wave equation |
---|
| 99 | |
---|
| 100 | msg = 'First conserved quantity must be "stage"' |
---|
| 101 | assert self.conserved_quantities[0] == 'stage', msg |
---|
| 102 | msg = 'Second conserved quantity must be "xmomentum"' |
---|
| 103 | assert self.conserved_quantities[1] == 'xmomentum', msg |
---|
| 104 | |
---|
| 105 | msg = 'First evolved quantity must be "stage"' |
---|
| 106 | assert self.evolved_quantities[0] == 'stage', msg |
---|
| 107 | msg = 'Second evolved quantity must be "xmomentum"' |
---|
| 108 | assert self.evolved_quantities[1] == 'xmomentum', msg |
---|
| 109 | msg = 'Third evolved quantity must be "elevation"' |
---|
| 110 | assert self.evolved_quantities[2] == 'elevation', msg |
---|
| 111 | msg = 'Fourth evolved quantity must be "height"' |
---|
| 112 | assert self.evolved_quantities[3] == 'height', msg |
---|
| 113 | msg = 'Fifth evolved quantity must be "velocity"' |
---|
| 114 | assert self.evolved_quantities[4] == 'velocity', msg |
---|
| 115 | |
---|
| 116 | Generic_domain.check_integrity(self) |
---|
| 117 | |
---|
| 118 | def compute_fluxes(self): |
---|
| 119 | #Call correct module function |
---|
| 120 | #(either from this module or C-extension) |
---|
| 121 | compute_fluxes_vel(self) |
---|
| 122 | |
---|
| 123 | def distribute_to_vertices_and_edges(self): |
---|
| 124 | #Call correct module function |
---|
| 125 | #(either from this module or C-extension) |
---|
| 126 | distribute_to_vertices_and_edges_limit_w_u(self) |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | #=============== End of Shallow Water Domain =============================== |
---|
| 130 | #----------------------------------- |
---|
| 131 | # Compute fluxes interface |
---|
| 132 | #----------------------------------- |
---|
| 133 | def compute_fluxes(domain): |
---|
| 134 | """ |
---|
| 135 | Python version of compute fluxes (local_compute_fluxes) |
---|
| 136 | is available in test_shallow_water_vel_domain.py |
---|
| 137 | """ |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | from numpy import zeros |
---|
| 141 | import sys |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | timestep = float(sys.maxint) |
---|
| 145 | |
---|
| 146 | stage = domain.quantities['stage'] |
---|
| 147 | xmom = domain.quantities['xmomentum'] |
---|
| 148 | bed = domain.quantities['elevation'] |
---|
| 149 | |
---|
| 150 | |
---|
[7840] | 151 | from anuga_1d.sww_flow.sww_vel_comp_flux_ext import compute_fluxes_ext |
---|
[7839] | 152 | |
---|
| 153 | domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed) |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | #----------------------------------- |
---|
| 157 | # Compute flux definition with vel |
---|
| 158 | #----------------------------------- |
---|
| 159 | def compute_fluxes_vel(domain): |
---|
| 160 | from Numeric import zeros, Float |
---|
| 161 | import sys |
---|
| 162 | |
---|
| 163 | |
---|
| 164 | timestep = float(sys.maxint) |
---|
| 165 | |
---|
| 166 | stage = domain.quantities['stage'] |
---|
| 167 | xmom = domain.quantities['xmomentum'] |
---|
| 168 | bed = domain.quantities['elevation'] |
---|
| 169 | height = domain.quantities['height'] |
---|
| 170 | velocity = domain.quantities['velocity'] |
---|
| 171 | |
---|
| 172 | |
---|
[7840] | 173 | from anuga_1d.sww.sww_vel_comp_flux_ext import compute_fluxes_vel_ext |
---|
[7839] | 174 | |
---|
| 175 | domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity) |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | #-------------------------------------------------------------------------- |
---|
| 180 | def distribute_to_vertices_and_edges_limit_w_u(domain): |
---|
| 181 | """Distribution from centroids to vertices specific to the |
---|
| 182 | shallow water wave |
---|
| 183 | equation. |
---|
| 184 | |
---|
| 185 | It will ensure that h (w-z) is always non-negative even in the |
---|
| 186 | presence of steep bed-slopes by taking a weighted average between shallow |
---|
| 187 | and deep cases. |
---|
| 188 | |
---|
| 189 | In addition, all conserved quantities get distributed as per either a |
---|
| 190 | constant (order==1) or a piecewise linear function (order==2). |
---|
| 191 | |
---|
| 192 | FIXME: more explanation about removal of artificial variability etc |
---|
| 193 | |
---|
| 194 | Precondition: |
---|
| 195 | All quantities defined at centroids and bed elevation defined at |
---|
| 196 | vertices. |
---|
| 197 | |
---|
| 198 | Postcondition |
---|
| 199 | Conserved quantities defined at vertices |
---|
| 200 | |
---|
| 201 | """ |
---|
| 202 | |
---|
| 203 | |
---|
| 204 | #Remove very thin layers of water |
---|
| 205 | #protect_against_infinitesimal_and_negative_heights(domain) |
---|
| 206 | |
---|
| 207 | import sys |
---|
| 208 | from numpy import zeros |
---|
[7840] | 209 | from anuga_1d.config import epsilon, h0 |
---|
[7839] | 210 | |
---|
| 211 | N = domain.number_of_elements |
---|
| 212 | |
---|
| 213 | #Shortcuts |
---|
| 214 | Stage = domain.quantities['stage'] |
---|
| 215 | Xmom = domain.quantities['xmomentum'] |
---|
| 216 | Bed = domain.quantities['elevation'] |
---|
| 217 | Height = domain.quantities['height'] |
---|
| 218 | Velocity = domain.quantities['velocity'] |
---|
| 219 | |
---|
| 220 | #Arrays |
---|
| 221 | w_C = Stage.centroid_values |
---|
| 222 | uh_C = Xmom.centroid_values |
---|
| 223 | z_C = Bed.centroid_values |
---|
| 224 | h_C = Height.centroid_values |
---|
| 225 | u_C = Velocity.centroid_values |
---|
| 226 | |
---|
| 227 | #print id(h_C) |
---|
| 228 | ## for i in range(N): |
---|
| 229 | ## h_C[i] = w_C[i] - z_C[i] |
---|
| 230 | ## if h_C[i] <= 1.0e-12: |
---|
| 231 | ## #print 'h_C[%d]= %15.5e\n' % (i,h_C[i]) |
---|
| 232 | ## h_C[i] = 0.0 |
---|
| 233 | ## w_C[i] = z_C[i] |
---|
| 234 | ## #uh_C[i] = 0.0 |
---|
| 235 | |
---|
| 236 | ## # u_C[i] = 0.0 |
---|
| 237 | ## # else: |
---|
| 238 | ## # u_C[i] = uh_C[i]/h_C[i] |
---|
| 239 | |
---|
| 240 | h0 = 1.0e-12 |
---|
| 241 | for i in range(N): |
---|
| 242 | h_C[i] = w_C[i] - z_C[i] |
---|
| 243 | if h_C[i] < 1.0e-12: |
---|
| 244 | u_C[i] = 0.0 #Could have been negative |
---|
| 245 | h_C[i] = 0.0 |
---|
| 246 | w_C[i] = z_C[i] |
---|
| 247 | else: |
---|
| 248 | #u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i]) |
---|
| 249 | u_C[i] = uh_C[i]/h_C[i] |
---|
| 250 | |
---|
| 251 | for name in [ 'velocity', 'stage' ]: |
---|
| 252 | Q = domain.quantities[name] |
---|
| 253 | if domain.order == 1: |
---|
| 254 | Q.extrapolate_first_order() |
---|
| 255 | elif domain.order == 2: |
---|
| 256 | Q.extrapolate_second_order() |
---|
| 257 | else: |
---|
| 258 | raise 'Unknown order' |
---|
| 259 | |
---|
| 260 | w_V = domain.quantities['stage'].vertex_values |
---|
| 261 | z_V = domain.quantities['elevation'].vertex_values |
---|
| 262 | h_V = domain.quantities['height'].vertex_values |
---|
| 263 | u_V = domain.quantities['velocity'].vertex_values |
---|
| 264 | uh_V = domain.quantities['xmomentum'].vertex_values |
---|
| 265 | |
---|
| 266 | #print w_V |
---|
| 267 | #print z_V |
---|
| 268 | |
---|
| 269 | h_V[:,:] = w_V - z_V |
---|
| 270 | for i in range(len(h_C)): |
---|
| 271 | for j in range(2): |
---|
| 272 | if h_V[i,j] < 0.0 : |
---|
| 273 | #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j]) |
---|
| 274 | dh = h_V[i,j] |
---|
| 275 | h_V[i,j] = 0.0 |
---|
| 276 | w_V[i,j] = z_V[i,j] |
---|
| 277 | h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh |
---|
| 278 | w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh |
---|
| 279 | |
---|
| 280 | uh_V[:,:] = u_V * h_V |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | return |
---|
| 284 | |
---|
| 285 | #--------------------------------------------------------------------------- |
---|
| 286 | def distribute_to_vertices_and_edges_limit_w_uh(domain): |
---|
| 287 | """Distribution from centroids to vertices specific to the |
---|
| 288 | shallow water wave equation. |
---|
| 289 | |
---|
| 290 | In addition, all conserved quantities get distributed as per either a |
---|
| 291 | constant (order==1) or a piecewise linear function (order==2). |
---|
| 292 | |
---|
| 293 | Precondition: |
---|
| 294 | All quantities defined at centroids and bed elevation defined at |
---|
| 295 | vertices. |
---|
| 296 | |
---|
| 297 | Postcondition |
---|
| 298 | Conserved quantities defined at vertices |
---|
| 299 | |
---|
| 300 | """ |
---|
| 301 | |
---|
| 302 | import sys |
---|
| 303 | from numpy import zeros |
---|
[7840] | 304 | from anuga_1d.config import epsilon, h0 |
---|
[7839] | 305 | |
---|
| 306 | N = domain.number_of_elements |
---|
| 307 | |
---|
| 308 | #Shortcuts |
---|
| 309 | Stage = domain.quantities['stage'] |
---|
| 310 | Xmom = domain.quantities['xmomentum'] |
---|
| 311 | Bed = domain.quantities['elevation'] |
---|
| 312 | Height = domain.quantities['height'] |
---|
| 313 | Velocity = domain.quantities['velocity'] |
---|
| 314 | |
---|
| 315 | #Arrays |
---|
| 316 | w_C = Stage.centroid_values |
---|
| 317 | uh_C = Xmom.centroid_values |
---|
| 318 | z_C = Bed.centroid_values |
---|
| 319 | h_C = Height.centroid_values |
---|
| 320 | u_C = Velocity.centroid_values |
---|
| 321 | |
---|
| 322 | |
---|
| 323 | for i in range(N): |
---|
| 324 | h_C[i] = w_C[i] - z_C[i] |
---|
| 325 | if h_C[i] <= 1.0e-6: |
---|
| 326 | #print 'h_C[%d]= %15.5e\n' % (i,h_C[i]) |
---|
| 327 | h_C[i] = 0.0 |
---|
| 328 | w_C[i] = z_C[i] |
---|
| 329 | uh_C[i] = 0.0 |
---|
| 330 | |
---|
| 331 | for name in [ 'stage', 'xmomentum']: |
---|
| 332 | Q = domain.quantities[name] |
---|
| 333 | if domain.order == 1: |
---|
| 334 | Q.extrapolate_first_order() |
---|
| 335 | elif domain.order == 2: |
---|
| 336 | Q.extrapolate_second_order() |
---|
| 337 | else: |
---|
| 338 | raise 'Unknown order' |
---|
| 339 | |
---|
| 340 | w_V = domain.quantities['stage'].vertex_values |
---|
| 341 | z_V = domain.quantities['elevation'].vertex_values |
---|
| 342 | h_V = domain.quantities['height'].vertex_values |
---|
| 343 | u_V = domain.quantities['velocity'].vertex_values |
---|
| 344 | uh_V = domain.quantities['xmomentum'].vertex_values |
---|
| 345 | |
---|
| 346 | h_V[:,:] = w_V - z_V |
---|
| 347 | |
---|
| 348 | for i in range(len(h_C)): |
---|
| 349 | for j in range(2): |
---|
| 350 | if h_V[i,j] < 0.0 : |
---|
| 351 | #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j]) |
---|
| 352 | dh = h_V[i,j] |
---|
| 353 | h_V[i,j] = 0.0 |
---|
| 354 | w_V[i,j] = z_V[i,j] |
---|
| 355 | h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh |
---|
| 356 | w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh |
---|
| 357 | u_V[i,j] = 0.0 |
---|
| 358 | if h_V[i,j] < h0: |
---|
| 359 | u_V[i,j] |
---|
| 360 | else: |
---|
| 361 | u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j]) |
---|
| 362 | |
---|
| 363 | |
---|