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