| 1 | |
|---|
| 2 | from anuga.shallow_water.shallow_water_domain import * |
|---|
| 3 | from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain |
|---|
| 4 | import numpy |
|---|
| 5 | |
|---|
| 6 | ############################################################################## |
|---|
| 7 | # Shallow Water Balanced Domain -- alternative implementation |
|---|
| 8 | # |
|---|
| 9 | # FIXME: Following the methods in CITE MODSIM PAPER |
|---|
| 10 | # |
|---|
| 11 | ############################################################################## |
|---|
| 12 | |
|---|
| 13 | class Domain(Sww_domain): |
|---|
| 14 | |
|---|
| 15 | def __init__(self, |
|---|
| 16 | coordinates=None, |
|---|
| 17 | vertices=None, |
|---|
| 18 | boundary=None, |
|---|
| 19 | tagged_elements=None, |
|---|
| 20 | geo_reference=None, |
|---|
| 21 | use_inscribed_circle=False, |
|---|
| 22 | mesh_filename=None, |
|---|
| 23 | use_cache=False, |
|---|
| 24 | verbose=False, |
|---|
| 25 | full_send_dict=None, |
|---|
| 26 | ghost_recv_dict=None, |
|---|
| 27 | starttime=0.0, |
|---|
| 28 | processor=0, |
|---|
| 29 | numproc=1, |
|---|
| 30 | number_of_full_nodes=None, |
|---|
| 31 | number_of_full_triangles=None): |
|---|
| 32 | |
|---|
| 33 | conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
|---|
| 34 | |
|---|
| 35 | evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
|---|
| 36 | other_quantities = [ 'elevation', 'friction', 'height', |
|---|
| 37 | 'xvelocity', 'yvelocity', 'x', 'y' ] |
|---|
| 38 | |
|---|
| 39 | |
|---|
| 40 | Sww_domain.__init__(self, |
|---|
| 41 | coordinates = coordinates, |
|---|
| 42 | vertices = vertices, |
|---|
| 43 | boundary = boundary, |
|---|
| 44 | tagged_elements = tagged_elements, |
|---|
| 45 | geo_reference = geo_reference, |
|---|
| 46 | use_inscribed_circle = use_inscribed_circle, |
|---|
| 47 | mesh_filename = mesh_filename, |
|---|
| 48 | use_cache = use_cache, |
|---|
| 49 | verbose = verbose, |
|---|
| 50 | conserved_quantities = conserved_quantities, |
|---|
| 51 | evolved_quantities = evolved_quantities, |
|---|
| 52 | other_quantities = other_quantities, |
|---|
| 53 | full_send_dict = full_send_dict, |
|---|
| 54 | ghost_recv_dict = ghost_recv_dict, |
|---|
| 55 | starttime = starttime, |
|---|
| 56 | processor = processor, |
|---|
| 57 | numproc = numproc, |
|---|
| 58 | number_of_full_nodes = number_of_full_nodes, |
|---|
| 59 | number_of_full_triangles = number_of_full_triangles) |
|---|
| 60 | |
|---|
| 61 | #self.forcing_terms.append(manning_friction_implicit) |
|---|
| 62 | #------------------------------------------------ |
|---|
| 63 | # set some defaults |
|---|
| 64 | # Most of these override the options in config.py |
|---|
| 65 | #------------------------------------------------ |
|---|
| 66 | self.set_CFL(1.0) |
|---|
| 67 | self.set_use_kinematic_viscosity(False) |
|---|
| 68 | self.timestepping_method='rk2'#'rk2'#'euler'#'rk2' |
|---|
| 69 | # The following allows storage of the negative depths associated with this method |
|---|
| 70 | self.minimum_storable_height=-99999999999.0 |
|---|
| 71 | |
|---|
| 72 | self.use_edge_limiter=True |
|---|
| 73 | self.default_order=2 |
|---|
| 74 | self.extrapolate_velocity_second_order=True |
|---|
| 75 | |
|---|
| 76 | # Note that the extrapolation method used in quantity_ext.c (e.g. |
|---|
| 77 | # extrapolate_second_order_and_limit_by_edge) uses a constant value for |
|---|
| 78 | # all the betas. |
|---|
| 79 | self.beta_w=0.0 |
|---|
| 80 | self.beta_w_dry=0.0 |
|---|
| 81 | self.beta_uh=0.0 |
|---|
| 82 | self.beta_uh_dry=0.0 |
|---|
| 83 | self.beta_vh=0.0 |
|---|
| 84 | self.beta_vh_dry=0.0 |
|---|
| 85 | |
|---|
| 86 | #self.optimise_dry_cells=True |
|---|
| 87 | # "self.optimise_dry_cells=False" presently ensures that the stage is |
|---|
| 88 | # always >= minimum_bed_edge value. Actually, we should only need to |
|---|
| 89 | # apply 'False' on the very first time step (to deal with stage values |
|---|
| 90 | # that were initialised below the bed by the user). After this, the |
|---|
| 91 | # algorithm should take care of itself, and 'True' should be okay. |
|---|
| 92 | self.optimise_dry_cells=False |
|---|
| 93 | |
|---|
| 94 | # Because gravity is treated within the flux function, |
|---|
| 95 | # we remove it from the forcing terms. |
|---|
| 96 | #self.forcing_terms.remove(gravity) |
|---|
| 97 | |
|---|
| 98 | # We need the edge_coordinates for the extrapolation |
|---|
| 99 | self.edge_coordinates=self.get_edge_midpoint_coordinates() |
|---|
| 100 | |
|---|
| 101 | # We demand that vertex values are stored uniquely |
|---|
| 102 | self.set_store_vertices_smoothly(False) |
|---|
| 103 | |
|---|
| 104 | self.maximum_allowed_speed=0.0 |
|---|
| 105 | #self.forcing_terms.append(manning_friction_explicit) |
|---|
| 106 | #self.forcing_terms.remove(manning_friction_implicit) |
|---|
| 107 | |
|---|
| 108 | # Keep track of the fluxes through the boundaries |
|---|
| 109 | self.boundary_flux_integral=numpy.ndarray(1) |
|---|
| 110 | self.boundary_flux_integral[0]=0. |
|---|
| 111 | self.boundary_flux_sum=numpy.ndarray(1) |
|---|
| 112 | self.boundary_flux_sum[0]=0. |
|---|
| 113 | |
|---|
| 114 | print '##########################################################################' |
|---|
| 115 | print '#' |
|---|
| 116 | print '# Using shallow_water_balanced2 solver in anuga_work/development/gareth/experimental/balanced_dev/' |
|---|
| 117 | print "#" |
|---|
| 118 | print "# This solver is experimental. Here are some tips on using it" |
|---|
| 119 | print "#" |
|---|
| 120 | print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values" |
|---|
| 121 | print "# , as the latter can be completely misleading near strong gradients in the flow. " |
|---|
| 122 | print "# There is a util.py script in anuga_work/development/gareth/tests which might help you extract" |
|---|
| 123 | print "# quantities at centroid values from sww files." |
|---|
| 124 | print "# Note that to accuractely compute centroid values from sww files, the files need to store " |
|---|
| 125 | print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer" |
|---|
| 126 | print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect" |
|---|
| 127 | print "# ANUGA viewer as well -- I expect this would be lots of work)" |
|---|
| 128 | print "#" |
|---|
| 129 | print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set" |
|---|
| 130 | print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). " |
|---|
| 131 | print "#" |
|---|
| 132 | print "# 3) This solver is not expected to perform well in problems with very" |
|---|
| 133 | print "# shallow water flowing down steep slopes (such that the stage_centroid_value " |
|---|
| 134 | print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests" |
|---|
| 135 | print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) " |
|---|
| 136 | print "#" |
|---|
| 137 | print "# 4) This solver allows the stage_centroid_value to drop to the minimum bed_edge_value" |
|---|
| 138 | print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the" |
|---|
| 139 | print "# bed_centroid_value. This means that triangles store slightly more water than they are" |
|---|
| 140 | print "# typically interpreted to store, which might have significance in some applications." |
|---|
| 141 | print "#" |
|---|
| 142 | print "# You will probably be able to tell this is causing you problems by convergence testing" |
|---|
| 143 | print "#" |
|---|
| 144 | print '# 5) Note that many options in config.py have been overridden by the solver -- I have ' |
|---|
| 145 | print '# deliberately attempted to get the solver to perform well with consistent values of ' |
|---|
| 146 | print '# these parameters -- so I would advise against changing them unless you at least check that ' |
|---|
| 147 | print '# it really does improve things' |
|---|
| 148 | print '#' |
|---|
| 149 | print '##########################################################################' |
|---|
| 150 | |
|---|
| 151 | #------------------------------- |
|---|
| 152 | # Specialisation of Evolve |
|---|
| 153 | #------------------------------- |
|---|
| 154 | # This alters the 'evolve' method from shallow_water_domain.py so that just |
|---|
| 155 | # before the file-output step, the vertex values are replaced with the |
|---|
| 156 | # centroid values. |
|---|
| 157 | # This is useful so that we can unambigously know what the centroid values |
|---|
| 158 | # of various parameters are |
|---|
| 159 | |
|---|
| 160 | def evolve(self, |
|---|
| 161 | yieldstep=None, |
|---|
| 162 | finaltime=None, |
|---|
| 163 | duration=None, |
|---|
| 164 | skip_initial_step=False): |
|---|
| 165 | """Specialisation of basic evolve method from parent class. |
|---|
| 166 | |
|---|
| 167 | Evolve the model by 1 step. |
|---|
| 168 | """ |
|---|
| 169 | |
|---|
| 170 | # Call check integrity here rather than from user scripts |
|---|
| 171 | # self.check_integrity() |
|---|
| 172 | |
|---|
| 173 | #print 'Into Evolve' |
|---|
| 174 | |
|---|
| 175 | msg = 'Attribute self.beta_w must be in the interval [0, 2]' |
|---|
| 176 | assert 0 <= self.beta_w <= 2.0, msg |
|---|
| 177 | |
|---|
| 178 | # Initial update of vertex and edge values before any STORAGE |
|---|
| 179 | # and or visualisation. |
|---|
| 180 | # This is done again in the initialisation of the Generic_Domain |
|---|
| 181 | # evolve loop but we do it here to ensure the values are ok for storage. |
|---|
| 182 | self.distribute_to_vertices_and_edges() |
|---|
| 183 | |
|---|
| 184 | if self.store is True and self.get_time() == self.get_starttime(): |
|---|
| 185 | self.initialise_storage() |
|---|
| 186 | #print 'Into Generic_Domain Evolve' |
|---|
| 187 | # Call basic machinery from parent class |
|---|
| 188 | for t in Generic_Domain.evolve(self, yieldstep=yieldstep, |
|---|
| 189 | finaltime=finaltime, duration=duration, |
|---|
| 190 | skip_initial_step=skip_initial_step): |
|---|
| 191 | #print 'Out of Generic_Domain Evolve' |
|---|
| 192 | # Store model data, e.g. for subsequent visualisation |
|---|
| 193 | if self.store is True: |
|---|
| 194 | # FIXME: Extrapolation is done before writing, because I had |
|---|
| 195 | # trouble correctly computing the centroid values from the |
|---|
| 196 | # vertex values (an 'average' did not seem to work correctly in |
|---|
| 197 | # very shallow cells -- there was a discrepency between |
|---|
| 198 | # domain.quantity['blah'].centroid_values and the values |
|---|
| 199 | # computed from the sww using the vertex averge). There should |
|---|
| 200 | # be a much much more disk-efficient way to store the centroid |
|---|
| 201 | # values than this |
|---|
| 202 | self.extrapolate_second_order_edge_sw() |
|---|
| 203 | |
|---|
| 204 | # Store the timestep |
|---|
| 205 | self.store_timestep() |
|---|
| 206 | |
|---|
| 207 | # Pass control on to outer loop for more specific actions |
|---|
| 208 | yield(t) |
|---|
| 209 | |
|---|
| 210 | #----------------- |
|---|
| 211 | # Flux computation |
|---|
| 212 | #----------------- |
|---|
| 213 | |
|---|
| 214 | ## @brief Compute fluxes and timestep suitable for all volumes in domain. |
|---|
| 215 | # @param domain The domain to calculate fluxes for. |
|---|
| 216 | def compute_fluxes(domain): |
|---|
| 217 | """Compute fluxes and timestep suitable for all volumes in domain. |
|---|
| 218 | |
|---|
| 219 | Compute total flux for each conserved quantity using "flux_function" |
|---|
| 220 | |
|---|
| 221 | Fluxes across each edge are scaled by edgelengths and summed up |
|---|
| 222 | Resulting flux is then scaled by area and stored in |
|---|
| 223 | explicit_update for each of the three conserved quantities |
|---|
| 224 | stage, xmomentum and ymomentum |
|---|
| 225 | |
|---|
| 226 | The maximal allowable speed computed by the flux_function for each volume |
|---|
| 227 | is converted to a timestep that must not be exceeded. The minimum of |
|---|
| 228 | those is computed as the next overall timestep. |
|---|
| 229 | |
|---|
| 230 | Post conditions: |
|---|
| 231 | domain.explicit_update is reset to computed flux values |
|---|
| 232 | domain.timestep is set to the largest step satisfying all volumes. |
|---|
| 233 | |
|---|
| 234 | This wrapper calls the underlying C version of compute fluxes |
|---|
| 235 | """ |
|---|
| 236 | |
|---|
| 237 | import sys |
|---|
| 238 | from swb2_domain_ext import compute_fluxes_ext_central \ |
|---|
| 239 | as compute_fluxes_ext |
|---|
| 240 | |
|---|
| 241 | #print "." |
|---|
| 242 | |
|---|
| 243 | # Shortcuts |
|---|
| 244 | #if(domain.timestepping_method!='rk2'): |
|---|
| 245 | # err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\ |
|---|
| 246 | # 'You need to use rk2 timestepping with this solver, ' +\ |
|---|
| 247 | # ' because there is presently a hack in compute fluxes central. \n'+\ |
|---|
| 248 | # ' The HACK: The timestep will only ' +\ |
|---|
| 249 | # 'be recomputed on every 2nd call to compute_fluxes_central, to support'+\ |
|---|
| 250 | # ' correct treatment of wetting and drying' |
|---|
| 251 | |
|---|
| 252 | # raise Exception, err_mess |
|---|
| 253 | |
|---|
| 254 | if(domain.timestepping_method=='rk2'): |
|---|
| 255 | timestep_order=2 |
|---|
| 256 | elif(domain.timestepping_method=='euler'): |
|---|
| 257 | timestep_order=1 |
|---|
| 258 | #elif(domain.timestepping_method=='rk3'): |
|---|
| 259 | # timestep_order=3 |
|---|
| 260 | else: |
|---|
| 261 | err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver' |
|---|
| 262 | raise Exception, err_mess |
|---|
| 263 | |
|---|
| 264 | #print 'timestep_order=', timestep_order |
|---|
| 265 | |
|---|
| 266 | Stage = domain.quantities['stage'] |
|---|
| 267 | Xmom = domain.quantities['xmomentum'] |
|---|
| 268 | Ymom = domain.quantities['ymomentum'] |
|---|
| 269 | Bed = domain.quantities['elevation'] |
|---|
| 270 | |
|---|
| 271 | timestep = float(sys.maxint) |
|---|
| 272 | |
|---|
| 273 | flux_timestep = compute_fluxes_ext(timestep, |
|---|
| 274 | domain.epsilon, |
|---|
| 275 | domain.minimum_allowed_height, |
|---|
| 276 | domain.g, |
|---|
| 277 | domain.boundary_flux_sum, |
|---|
| 278 | domain.neighbours, |
|---|
| 279 | domain.neighbour_edges, |
|---|
| 280 | domain.normals, |
|---|
| 281 | domain.edgelengths, |
|---|
| 282 | domain.radii, |
|---|
| 283 | domain.areas, |
|---|
| 284 | domain.tri_full_flag, |
|---|
| 285 | Stage.edge_values, |
|---|
| 286 | Xmom.edge_values, |
|---|
| 287 | Ymom.edge_values, |
|---|
| 288 | Bed.edge_values, |
|---|
| 289 | Stage.boundary_values, |
|---|
| 290 | Xmom.boundary_values, |
|---|
| 291 | Ymom.boundary_values, |
|---|
| 292 | domain.boundary_flux_type, |
|---|
| 293 | Stage.explicit_update, |
|---|
| 294 | Xmom.explicit_update, |
|---|
| 295 | Ymom.explicit_update, |
|---|
| 296 | domain.already_computed_flux, |
|---|
| 297 | domain.max_speed, |
|---|
| 298 | int(domain.optimise_dry_cells), |
|---|
| 299 | timestep_order, |
|---|
| 300 | Stage.centroid_values, |
|---|
| 301 | Xmom.centroid_values, |
|---|
| 302 | Ymom.centroid_values, |
|---|
| 303 | Bed.centroid_values, |
|---|
| 304 | Bed.vertex_values) |
|---|
| 305 | |
|---|
| 306 | #import pdb |
|---|
| 307 | #pdb.set_trace() |
|---|
| 308 | #print 'flux timestep: ', flux_timestep, domain.timestep |
|---|
| 309 | if(domain.timestepping_method=='rk2'): |
|---|
| 310 | if(flux_timestep == float(sys.maxint)): |
|---|
| 311 | domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\ |
|---|
| 312 | domain.boundary_flux_sum[0]*domain.timestep*0.5 |
|---|
| 313 | #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum |
|---|
| 314 | domain.boundary_flux_sum[0]=0. |
|---|
| 315 | |
|---|
| 316 | domain.flux_timestep = flux_timestep |
|---|
| 317 | |
|---|
| 318 | |
|---|
| 319 | |
|---|
| 320 | def protect_against_infinitesimal_and_negative_heights(domain): |
|---|
| 321 | """protect against infinitesimal heights and associated high velocities""" |
|---|
| 322 | |
|---|
| 323 | from swb2_domain_ext import protect |
|---|
| 324 | #print'using swb2_protect_against ..' |
|---|
| 325 | |
|---|
| 326 | # shortcuts |
|---|
| 327 | wc = domain.quantities['stage'].centroid_values |
|---|
| 328 | wv = domain.quantities['stage'].vertex_values |
|---|
| 329 | zc = domain.quantities['elevation'].centroid_values |
|---|
| 330 | zv = domain.quantities['elevation'].vertex_values |
|---|
| 331 | xmomc = domain.quantities['xmomentum'].centroid_values |
|---|
| 332 | ymomc = domain.quantities['ymomentum'].centroid_values |
|---|
| 333 | areas = domain.areas |
|---|
| 334 | xc = domain.centroid_coordinates[:,0] |
|---|
| 335 | yc = domain.centroid_coordinates[:,1] |
|---|
| 336 | |
|---|
| 337 | protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, |
|---|
| 338 | domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas, xc, yc) |
|---|
| 339 | |
|---|
| 340 | def conserved_values_to_evolved_values(self, q_cons, q_evol): |
|---|
| 341 | """Mapping between conserved quantities and the evolved quantities. |
|---|
| 342 | Used where we have a boundary condition which works with conserved |
|---|
| 343 | quantities and we now want to use them for the new well balanced |
|---|
| 344 | code using the evolved quantities |
|---|
| 345 | |
|---|
| 346 | Typically the old boundary condition will set the values in q_cons, |
|---|
| 347 | |
|---|
| 348 | q_evol on input will have the values of the evolved quantities at the |
|---|
| 349 | edge point (useful as it provides values for evlevation). |
|---|
| 350 | |
|---|
| 351 | """ |
|---|
| 352 | |
|---|
| 353 | wc = q_cons[0] |
|---|
| 354 | uhc = q_cons[1] |
|---|
| 355 | vhc = q_cons[2] |
|---|
| 356 | |
|---|
| 357 | |
|---|
| 358 | q_evol[0] = wc |
|---|
| 359 | q_evol[1] = uhc |
|---|
| 360 | q_evol[2] = vhc |
|---|
| 361 | |
|---|
| 362 | return q_evol |
|---|
| 363 | |
|---|
| 364 | def distribute_to_vertices_and_edges(self): |
|---|
| 365 | """ Call correct module function """ |
|---|
| 366 | |
|---|
| 367 | if self.use_edge_limiter: |
|---|
| 368 | #distribute_using_edge_limiter(self) |
|---|
| 369 | self.distribute_using_edge_limiter() |
|---|
| 370 | else: |
|---|
| 371 | #distribute_using_vertex_limiter(self) |
|---|
| 372 | self.distribute_using_vertex_limiter() |
|---|
| 373 | |
|---|
| 374 | |
|---|
| 375 | def distribute_using_vertex_limiter(domain): |
|---|
| 376 | """Distribution from centroids to vertices specific to the SWW equation. |
|---|
| 377 | |
|---|
| 378 | Precondition: |
|---|
| 379 | All quantities defined at centroids and bed elevation defined at |
|---|
| 380 | vertices. |
|---|
| 381 | |
|---|
| 382 | Postcondition |
|---|
| 383 | Conserved quantities defined at vertices |
|---|
| 384 | """ |
|---|
| 385 | # Remove very thin layers of water |
|---|
| 386 | domain.protect_against_infinitesimal_and_negative_heights() |
|---|
| 387 | |
|---|
| 388 | # Extrapolate all conserved quantities |
|---|
| 389 | if domain.optimised_gradient_limiter: |
|---|
| 390 | # MH090605 if second order, |
|---|
| 391 | # perform the extrapolation and limiting on |
|---|
| 392 | # all of the conserved quantities |
|---|
| 393 | |
|---|
| 394 | if (domain._order_ == 1): |
|---|
| 395 | for name in domain.conserved_quantities: |
|---|
| 396 | Q = domain.quantities[name] |
|---|
| 397 | Q.extrapolate_first_order() |
|---|
| 398 | elif domain._order_ == 2: |
|---|
| 399 | domain.extrapolate_second_order_sw() |
|---|
| 400 | else: |
|---|
| 401 | raise Exception('Unknown order') |
|---|
| 402 | else: |
|---|
| 403 | # Old code: |
|---|
| 404 | for name in domain.conserved_quantities: |
|---|
| 405 | Q = domain.quantities[name] |
|---|
| 406 | |
|---|
| 407 | if domain._order_ == 1: |
|---|
| 408 | Q.extrapolate_first_order() |
|---|
| 409 | elif domain._order_ == 2: |
|---|
| 410 | Q.extrapolate_second_order_and_limit_by_vertex() |
|---|
| 411 | else: |
|---|
| 412 | raise Exception('Unknown order') |
|---|
| 413 | |
|---|
| 414 | # Take bed elevation into account when water heights are small |
|---|
| 415 | #balance_deep_and_shallow(domain) |
|---|
| 416 | |
|---|
| 417 | # Compute edge values by interpolation |
|---|
| 418 | for name in domain.conserved_quantities: |
|---|
| 419 | Q = domain.quantities[name] |
|---|
| 420 | Q.interpolate_from_vertices_to_edges() |
|---|
| 421 | |
|---|
| 422 | |
|---|
| 423 | def distribute_using_edge_limiter(domain): |
|---|
| 424 | """Distribution from centroids to edges specific to the SWW eqn. |
|---|
| 425 | |
|---|
| 426 | Precondition: |
|---|
| 427 | All quantities defined at centroids and bed elevation defined at |
|---|
| 428 | vertices. |
|---|
| 429 | |
|---|
| 430 | Postcondition |
|---|
| 431 | Conserved quantities defined at vertices |
|---|
| 432 | """ |
|---|
| 433 | #from swb2_domain import protect_against_infinitesimal_and_negative_heights |
|---|
| 434 | |
|---|
| 435 | # Remove very thin layers of water |
|---|
| 436 | #print 'Before protect' |
|---|
| 437 | domain.protect_against_infinitesimal_and_negative_heights() |
|---|
| 438 | #print 'After protect' |
|---|
| 439 | |
|---|
| 440 | #for name in domain.conserved_quantities: |
|---|
| 441 | # Q = domain.quantities[name] |
|---|
| 442 | # if domain._order_ == 1: |
|---|
| 443 | # Q.extrapolate_first_order() |
|---|
| 444 | # elif domain._order_ == 2: |
|---|
| 445 | # #Q.extrapolate_second_order_and_limit_by_edge() |
|---|
| 446 | # # FIXME: This use of 'break' is hacky |
|---|
| 447 | # domain.extrapolate_second_order_edge_sw() |
|---|
| 448 | # break |
|---|
| 449 | # else: |
|---|
| 450 | # raise Exception('Unknown order') |
|---|
| 451 | |
|---|
| 452 | #print 'Before extrapolate' |
|---|
| 453 | domain.extrapolate_second_order_edge_sw() |
|---|
| 454 | #print 'After extrapolate' |
|---|
| 455 | |
|---|
| 456 | #balance_deep_and_shallow(domain) |
|---|
| 457 | |
|---|
| 458 | # Compute vertex values by interpolation |
|---|
| 459 | #for name in domain.conserved_quantities: |
|---|
| 460 | # Q = domain.quantities[name] |
|---|
| 461 | # Q.interpolate_from_edges_to_vertices() |
|---|
| 462 | # Q.interpolate_from_vertices_to_edges() |
|---|
| 463 | |
|---|
| 464 | |
|---|
| 465 | def update_other_quantities(self): |
|---|
| 466 | """ There may be a need to calculates some of the other quantities |
|---|
| 467 | based on the new values of conserved quantities |
|---|
| 468 | |
|---|
| 469 | But that is not needed in this version of the solver. |
|---|
| 470 | """ |
|---|
| 471 | pass |
|---|
| 472 | # The centroid values of height and x and y velocity |
|---|
| 473 | # might not have been setup |
|---|
| 474 | |
|---|
| 475 | #self.update_centroids_of_velocities_and_height() |
|---|
| 476 | # |
|---|
| 477 | #for name in ['height', 'xvelocity', 'yvelocity']: |
|---|
| 478 | # Q = self.quantities[name] |
|---|
| 479 | # if self._order_ == 1: |
|---|
| 480 | # Q.extrapolate_first_order() |
|---|
| 481 | # elif self._order_ == 2: |
|---|
| 482 | # Q.extrapolate_second_order_and_limit_by_edge() |
|---|
| 483 | # else: |
|---|
| 484 | # raise Exception('Unknown order') |
|---|
| 485 | |
|---|
| 486 | |
|---|
| 487 | def update_centroids_of_velocities_and_height(self): |
|---|
| 488 | """Calculate the centroid values of velocities and height based |
|---|
| 489 | on the values of the quantities stage and x and y momentum |
|---|
| 490 | |
|---|
| 491 | Assumes that stage and momentum are up to date |
|---|
| 492 | |
|---|
| 493 | Useful for kinematic viscosity calculations |
|---|
| 494 | """ |
|---|
| 495 | pass |
|---|
| 496 | ## For shallow water we need to update height xvelocity and yvelocity |
|---|
| 497 | |
|---|
| 498 | ##Shortcuts |
|---|
| 499 | #W = self.quantities['stage'] |
|---|
| 500 | #UH = self.quantities['xmomentum'] |
|---|
| 501 | #VH = self.quantities['ymomentum'] |
|---|
| 502 | #H = self.quantities['height'] |
|---|
| 503 | #Z = self.quantities['elevation'] |
|---|
| 504 | #U = self.quantities['xvelocity'] |
|---|
| 505 | #V = self.quantities['yvelocity'] |
|---|
| 506 | |
|---|
| 507 | ##print num.min(W.centroid_values) |
|---|
| 508 | |
|---|
| 509 | ## Make sure boundary values of conserved quantites |
|---|
| 510 | ## are consistent with value of functions at centroids |
|---|
| 511 | ##self.distribute_to_vertices_and_edges() |
|---|
| 512 | #Z.set_boundary_values_from_edges() |
|---|
| 513 | |
|---|
| 514 | ##W.set_boundary_values_from_edges() |
|---|
| 515 | ##UH.set_boundary_values_from_edges() |
|---|
| 516 | ##VH.set_boundary_values_from_edges() |
|---|
| 517 | |
|---|
| 518 | ## Update height values |
|---|
| 519 | #H.set_values(W.centroid_values-Z.centroid_values, location='centroids') |
|---|
| 520 | #H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0, |
|---|
| 521 | # W.boundary_values-Z.boundary_values, 0.0)) |
|---|
| 522 | |
|---|
| 523 | ##assert num.min(H.centroid_values) >= 0 |
|---|
| 524 | ##assert num.min(H.boundary_values) >= 0 |
|---|
| 525 | |
|---|
| 526 | ##Aliases |
|---|
| 527 | #uh_C = UH.centroid_values |
|---|
| 528 | #vh_C = VH.centroid_values |
|---|
| 529 | #h_C = H.centroid_values |
|---|
| 530 | |
|---|
| 531 | #uh_B = UH.boundary_values |
|---|
| 532 | #vh_B = VH.boundary_values |
|---|
| 533 | #h_B = H.boundary_values |
|---|
| 534 | |
|---|
| 535 | #H0 = 1.0e-8 |
|---|
| 536 | # |
|---|
| 537 | #U.set_values(uh_C/(h_C + H0/h_C), location='centroids') |
|---|
| 538 | #V.set_values(vh_C/(h_C + H0/h_C), location='centroids') |
|---|
| 539 | |
|---|
| 540 | #U.set_boundary_values(uh_B/(h_B + H0/h_B)) |
|---|
| 541 | #V.set_boundary_values(vh_B/(h_B + H0/h_B)) |
|---|
| 542 | |
|---|
| 543 | |
|---|
| 544 | |
|---|
| 545 | def extrapolate_second_order_edge_sw(self): |
|---|
| 546 | """Wrapper calling C version of extrapolate_second_order_sw, using |
|---|
| 547 | edges instead of vertices in the extrapolation step. The routine |
|---|
| 548 | then interpolates from edges to vertices. |
|---|
| 549 | The momentum extrapolation / interpolation is based on either |
|---|
| 550 | momentum or velocity, depending on the choice of |
|---|
| 551 | extrapolate_velocity_second_order. |
|---|
| 552 | |
|---|
| 553 | self the domain to operate on |
|---|
| 554 | |
|---|
| 555 | """ |
|---|
| 556 | |
|---|
| 557 | from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2 |
|---|
| 558 | |
|---|
| 559 | # Shortcuts |
|---|
| 560 | Stage = self.quantities['stage'] |
|---|
| 561 | Xmom = self.quantities['xmomentum'] |
|---|
| 562 | Ymom = self.quantities['ymomentum'] |
|---|
| 563 | Elevation = self.quantities['elevation'] |
|---|
| 564 | |
|---|
| 565 | extrapol2(self, |
|---|
| 566 | self.surrogate_neighbours, |
|---|
| 567 | self.neighbour_edges, |
|---|
| 568 | self.number_of_boundaries, |
|---|
| 569 | self.centroid_coordinates, |
|---|
| 570 | Stage.centroid_values, |
|---|
| 571 | Xmom.centroid_values, |
|---|
| 572 | Ymom.centroid_values, |
|---|
| 573 | Elevation.centroid_values, |
|---|
| 574 | self.edge_coordinates, |
|---|
| 575 | Stage.edge_values, |
|---|
| 576 | Xmom.edge_values, |
|---|
| 577 | Ymom.edge_values, |
|---|
| 578 | Elevation.edge_values, |
|---|
| 579 | Stage.vertex_values, |
|---|
| 580 | Xmom.vertex_values, |
|---|
| 581 | Ymom.vertex_values, |
|---|
| 582 | Elevation.vertex_values, |
|---|
| 583 | int(self.optimise_dry_cells), |
|---|
| 584 | int(self.extrapolate_velocity_second_order)) |
|---|
| 585 | |
|---|
| 586 | def set_boundary(self, boundary_map): |
|---|
| 587 | ## Specialisation of set_boundary, which also updates the 'boundary_flux_type' |
|---|
| 588 | Sww_domain.set_boundary(self,boundary_map) |
|---|
| 589 | |
|---|
| 590 | # Add a flag which can be used to distinguish flux boundaries within |
|---|
| 591 | # compute_fluxes_central |
|---|
| 592 | # Initialise to zero (which means 'not a flux_boundary') |
|---|
| 593 | self.boundary_flux_type = self.boundary_edges*0 |
|---|
| 594 | |
|---|
| 595 | # HACK to set the values of domain.boundary_flux |
|---|
| 596 | for i in range(len(self.boundary_objects)): |
|---|
| 597 | # Record the first 10 characters of the name of the boundary. |
|---|
| 598 | # FIXME: There must be a better way than this! |
|---|
| 599 | bndry_name = self.boundary_objects[i][1].__repr__()[0:10] |
|---|
| 600 | |
|---|
| 601 | if(bndry_name=='zero_mass_'): |
|---|
| 602 | # Create flag 'boundary_flux_type' that can be read in |
|---|
| 603 | # compute_fluxes_central |
|---|
| 604 | self.boundary_flux_type[i]=1 |
|---|
| 605 | |
|---|
| 606 | |
|---|
| 607 | #def manning_friction_implicit(domain): |
|---|
| 608 | # """Apply (Manning) friction to water momentum |
|---|
| 609 | # Wrapper for c version |
|---|
| 610 | # """ |
|---|
| 611 | # |
|---|
| 612 | # from shallow_water_ext import manning_friction_flat |
|---|
| 613 | # from shallow_water_ext import manning_friction_sloped |
|---|
| 614 | # |
|---|
| 615 | # xmom = domain.quantities['xmomentum'] |
|---|
| 616 | # ymom = domain.quantities['ymomentum'] |
|---|
| 617 | # |
|---|
| 618 | # x = domain.get_vertex_coordinates() |
|---|
| 619 | # |
|---|
| 620 | # w = domain.quantities['stage'].centroid_values |
|---|
| 621 | # z = domain.quantities['elevation'].vertex_values |
|---|
| 622 | # |
|---|
| 623 | # uh = xmom.centroid_values |
|---|
| 624 | # vh = ymom.centroid_values |
|---|
| 625 | # eta = domain.quantities['friction'].centroid_values |
|---|
| 626 | # |
|---|
| 627 | # xmom_update = xmom.semi_implicit_update |
|---|
| 628 | # ymom_update = ymom.semi_implicit_update |
|---|
| 629 | # |
|---|
| 630 | # eps = domain.epsilon |
|---|
| 631 | # g = domain.g |
|---|
| 632 | # |
|---|
| 633 | # if domain.use_sloped_mannings: |
|---|
| 634 | # manning_friction_sloped(g, eps, x, w, uh, vh, z, eta, xmom_update, \ |
|---|
| 635 | # ymom_update) |
|---|
| 636 | # else: |
|---|
| 637 | # manning_friction_flat(g, eps, w, uh, vh, z, eta, xmom_update, \ |
|---|
| 638 | # #ymom_update) |
|---|
