[7559] | 1 | |
---|
| 2 | from anuga.shallow_water.shallow_water_domain import * |
---|
| 3 | from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain |
---|
| 4 | |
---|
[7573] | 5 | from swb_boundary_conditions import Transmissive_boundary |
---|
[7559] | 6 | |
---|
| 7 | ############################################################################## |
---|
| 8 | # Shallow Water Balanced Domain |
---|
| 9 | # |
---|
| 10 | # Uses extra evolved quantities height, elevation, xvelocity, yvelocity |
---|
| 11 | ############################################################################## |
---|
| 12 | |
---|
| 13 | ## |
---|
| 14 | # @brief Class for a shallow water balanced domain. |
---|
| 15 | class Domain(Sww_domain): |
---|
| 16 | |
---|
| 17 | ## |
---|
[7573] | 18 | # @brief Instantiate a shallow water balanced domain. |
---|
[7559] | 19 | # @param coordinates |
---|
| 20 | # @param vertices |
---|
| 21 | # @param boundary |
---|
| 22 | # @param tagged_elements |
---|
| 23 | # @param geo_reference |
---|
| 24 | # @param use_inscribed_circle |
---|
| 25 | # @param mesh_filename |
---|
| 26 | # @param use_cache |
---|
| 27 | # @param verbose |
---|
| 28 | # @param full_send_dict |
---|
| 29 | # @param ghost_recv_dict |
---|
| 30 | # @param processor |
---|
| 31 | # @param numproc |
---|
| 32 | # @param number_of_full_nodes |
---|
| 33 | # @param number_of_full_triangles |
---|
| 34 | def __init__(self, |
---|
| 35 | coordinates=None, |
---|
| 36 | vertices=None, |
---|
| 37 | boundary=None, |
---|
| 38 | tagged_elements=None, |
---|
| 39 | geo_reference=None, |
---|
| 40 | use_inscribed_circle=False, |
---|
| 41 | mesh_filename=None, |
---|
| 42 | use_cache=False, |
---|
| 43 | verbose=False, |
---|
| 44 | full_send_dict=None, |
---|
| 45 | ghost_recv_dict=None, |
---|
| 46 | processor=0, |
---|
| 47 | numproc=1, |
---|
| 48 | number_of_full_nodes=None, |
---|
| 49 | number_of_full_triangles=None): |
---|
| 50 | |
---|
[7573] | 51 | conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
---|
| 52 | |
---|
[7559] | 53 | evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum', \ |
---|
| 54 | 'height', 'elevation', 'xvelocity', 'yvelocity'] |
---|
| 55 | |
---|
| 56 | other_quantities = [ 'friction' ] |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | Sww_domain.__init__(self, |
---|
| 60 | coordinates = coordinates, |
---|
| 61 | vertices = vertices, |
---|
| 62 | boundary = boundary, |
---|
| 63 | tagged_elements = tagged_elements, |
---|
| 64 | geo_reference = geo_reference, |
---|
| 65 | use_inscribed_circle = use_inscribed_circle, |
---|
| 66 | mesh_filename = mesh_filename, |
---|
| 67 | use_cache = use_cache, |
---|
| 68 | verbose = verbose, |
---|
[7573] | 69 | conserved_quantities = conserved_quantities, |
---|
[7559] | 70 | evolved_quantities = evolved_quantities, |
---|
| 71 | other_quantities = other_quantities, |
---|
| 72 | full_send_dict = full_send_dict, |
---|
| 73 | ghost_recv_dict = ghost_recv_dict, |
---|
| 74 | processor = processor, |
---|
| 75 | numproc = numproc, |
---|
| 76 | number_of_full_nodes = number_of_full_nodes, |
---|
| 77 | number_of_full_triangles = number_of_full_triangles) |
---|
| 78 | |
---|
| 79 | #--------------------- |
---|
| 80 | # set some defaults |
---|
| 81 | #--------------------- |
---|
[7573] | 82 | self.set_timestepping_method(1) |
---|
| 83 | self.set_default_order(2) |
---|
[7559] | 84 | self.set_new_mannings_function(True) |
---|
[7573] | 85 | self.set_centroid_transmissive_bc(True) |
---|
[7559] | 86 | |
---|
[7573] | 87 | ## |
---|
| 88 | # @brief Run integrity checks on shallow water balanced domain. |
---|
| 89 | def check_integrity(self): |
---|
| 90 | Sww_domain.check_integrity(self) |
---|
[7559] | 91 | |
---|
[7573] | 92 | #Check that the evolved quantities are correct (order) |
---|
| 93 | msg = 'First evolved quantity must be "stage"' |
---|
| 94 | assert self.evolved_quantities[0] == 'stage', msg |
---|
| 95 | msg = 'Second evolved quantity must be "xmomentum"' |
---|
| 96 | assert self.evolved_quantities[1] == 'xmomentum', msg |
---|
| 97 | msg = 'Third evolved quantity must be "ymomentum"' |
---|
| 98 | assert self.evolved_quantities[2] == 'ymomentum', msg |
---|
| 99 | msg = 'Fourth evolved quantity must be "height"' |
---|
| 100 | assert self.evolved_quantities[3] == 'height', msg |
---|
| 101 | msg = 'Fifth evolved quantity must be "elevation"' |
---|
| 102 | assert self.evolved_quantities[4] == 'elevation', msg |
---|
| 103 | msg = 'Sixth evolved quantity must be "xvelocity"' |
---|
| 104 | assert self.evolved_quantities[5] == 'xvelocity', msg |
---|
| 105 | msg = 'Seventh evolved quantity must be "yvelocity"' |
---|
| 106 | assert self.evolved_quantities[6] == 'yvelocity', msg |
---|
[7559] | 107 | |
---|
[7573] | 108 | msg = 'First other quantity must be "friction"' |
---|
| 109 | assert self.other_quantities[0] == 'friction', msg |
---|
| 110 | |
---|
| 111 | ## |
---|
| 112 | # @brief |
---|
| 113 | def compute_fluxes(self): |
---|
| 114 | #Call correct module function (either from this module or C-extension) |
---|
| 115 | |
---|
[7575] | 116 | from swb_domain_ext import compute_fluxes_c |
---|
| 117 | |
---|
| 118 | #Shortcuts |
---|
| 119 | W = self.quantities['stage'] |
---|
| 120 | UH = self.quantities['xmomentum'] |
---|
| 121 | VH = self.quantities['ymomentum'] |
---|
| 122 | H = self.quantities['height'] |
---|
| 123 | Z = self.quantities['elevation'] |
---|
| 124 | U = self.quantities['xvelocity'] |
---|
| 125 | V = self.quantities['yvelocity'] |
---|
| 126 | |
---|
| 127 | timestep = self.get_evolve_max_timestep() |
---|
| 128 | |
---|
| 129 | self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V) |
---|
| 130 | |
---|
[7573] | 131 | ## |
---|
| 132 | # @brief |
---|
| 133 | def distribute_to_vertices_and_edges(self): |
---|
| 134 | """Distribution from centroids to edges specific to the SWW eqn. |
---|
| 135 | |
---|
| 136 | It will ensure that h (w-z) is always non-negative even in the |
---|
| 137 | presence of steep bed-slopes by taking a weighted average between shallow |
---|
| 138 | and deep cases. |
---|
| 139 | |
---|
| 140 | In addition, all conserved quantities get distributed as per either a |
---|
| 141 | constant (order==1) or a piecewise linear function (order==2). |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | Precondition: |
---|
| 145 | All conserved quantities defined at centroids and bed elevation defined at |
---|
| 146 | edges. |
---|
| 147 | |
---|
| 148 | Postcondition |
---|
| 149 | Evolved quantities defined at vertices and edges |
---|
| 150 | """ |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | #Shortcuts |
---|
[7575] | 154 | W = self.quantities['stage'] |
---|
| 155 | UH = self.quantities['xmomentum'] |
---|
| 156 | VH = self.quantities['ymomentum'] |
---|
| 157 | H = self.quantities['height'] |
---|
| 158 | Z = self.quantities['elevation'] |
---|
| 159 | U = self.quantities['xvelocity'] |
---|
| 160 | V = self.quantities['yvelocity'] |
---|
[7573] | 161 | |
---|
| 162 | #Arrays |
---|
[7575] | 163 | w_C = W.centroid_values |
---|
| 164 | uh_C = UH.centroid_values |
---|
| 165 | vh_C = VH.centroid_values |
---|
| 166 | z_C = Z.centroid_values |
---|
| 167 | h_C = H.centroid_values |
---|
| 168 | u_C = U.centroid_values |
---|
| 169 | v_C = V.centroid_values |
---|
[7573] | 170 | |
---|
| 171 | w_C[:] = num.maximum(w_C, z_C) |
---|
| 172 | |
---|
| 173 | h_C[:] = w_C - z_C |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | assert num.min(h_C) >= 0 |
---|
| 177 | |
---|
| 178 | num.putmask(uh_C, h_C < 1.0e-15, 0.0) |
---|
| 179 | num.putmask(vh_C, h_C < 1.0e-15, 0.0) |
---|
| 180 | num.putmask(h_C, h_C < 1.0e-15, 1.0e-15) |
---|
| 181 | |
---|
| 182 | u_C[:] = uh_C/h_C |
---|
| 183 | v_C[:] = vh_C/h_C |
---|
| 184 | |
---|
[7575] | 185 | for name in [ 'height', 'elevation', 'xvelocity', 'yvelocity' ]: |
---|
[7573] | 186 | Q = self.quantities[name] |
---|
| 187 | if self._order_ == 1: |
---|
| 188 | Q.extrapolate_first_order() |
---|
| 189 | elif self._order_ == 2: |
---|
| 190 | Q.extrapolate_second_order_and_limit_by_edge() |
---|
| 191 | else: |
---|
| 192 | raise 'Unknown order' |
---|
| 193 | |
---|
| 194 | |
---|
[7575] | 195 | w_E = W.edge_values |
---|
| 196 | uh_E = UH.edge_values |
---|
| 197 | vh_E = VH.edge_values |
---|
| 198 | h_E = H.edge_values |
---|
| 199 | z_E = Z.edge_values |
---|
| 200 | u_E = U.edge_values |
---|
| 201 | v_E = V.edge_values |
---|
[7573] | 202 | |
---|
| 203 | |
---|
| 204 | w_E[:] = z_E + h_E |
---|
| 205 | |
---|
| 206 | #num.putmask(u_E, h_temp <= 0.0, 0.0) |
---|
| 207 | #num.putmask(v_E, h_temp <= 0.0, 0.0) |
---|
| 208 | #num.putmask(w_E, h_temp <= 0.0, z_E+h_E) |
---|
| 209 | #num.putmask(h_E, h_E <= 0.0, 0.0) |
---|
| 210 | |
---|
| 211 | uh_E[:] = u_E * h_E |
---|
| 212 | vh_E[:] = v_E * h_E |
---|
| 213 | |
---|
| 214 | """ |
---|
| 215 | print '==========================================================' |
---|
| 216 | print 'Time ', self.get_time() |
---|
| 217 | print h_E |
---|
| 218 | print uh_E |
---|
| 219 | print vh_E |
---|
| 220 | """ |
---|
| 221 | |
---|
| 222 | # Compute vertex values by interpolation |
---|
| 223 | for name in self.evolved_quantities: |
---|
| 224 | Q = self.quantities[name] |
---|
| 225 | Q.interpolate_from_edges_to_vertices() |
---|
| 226 | |
---|
| 227 | |
---|
[7575] | 228 | w_V = W.vertex_values |
---|
| 229 | uh_V = UH.vertex_values |
---|
| 230 | vh_V = VH.vertex_values |
---|
| 231 | z_V = Z.vertex_values |
---|
| 232 | h_V = H.vertex_values |
---|
| 233 | u_V = U.vertex_values |
---|
| 234 | v_V = V.vertex_values |
---|
[7573] | 235 | |
---|
| 236 | |
---|
| 237 | #w_V[:] = z_V + h_V |
---|
| 238 | |
---|
| 239 | #num.putmask(u_V, h_V <= 0.0, 0.0) |
---|
| 240 | #num.putmask(v_V, h_V <= 0.0, 0.0) |
---|
| 241 | #num.putmask(w_V, h_V <= 0.0, z_V) |
---|
| 242 | #num.putmask(h_V, h_V <= 0.0, 0.0) |
---|
| 243 | |
---|
| 244 | uh_V[:] = u_V * h_V |
---|
| 245 | vh_V[:] = v_V * h_V |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | ## |
---|
| 249 | # @brief Code to let us use old shallow water domain BCs |
---|
[7559] | 250 | def conserved_values_to_evolved_values(self, q_cons, q_evol): |
---|
| 251 | """Mapping between conserved quantities and the evolved quantities. |
---|
| 252 | Used where we have a boundary condition which works with conserved |
---|
| 253 | quantities and we now want to use them for the new well balanced |
---|
| 254 | code using the evolved quantities |
---|
| 255 | |
---|
| 256 | Typically the old boundary condition will set the values in q_cons, |
---|
| 257 | |
---|
| 258 | q_evol on input will have the values of the evolved quantities at the |
---|
| 259 | edge point (useful as it provides values for evlevation). |
---|
| 260 | """ |
---|
| 261 | |
---|
| 262 | wc = q_cons[0] |
---|
| 263 | uhc = q_cons[1] |
---|
| 264 | vhc = q_cons[2] |
---|
| 265 | |
---|
| 266 | we = q_evol[0] |
---|
| 267 | uhe = q_evol[1] |
---|
| 268 | vhe = q_evol[2] |
---|
| 269 | |
---|
| 270 | he = q_evol[3] |
---|
| 271 | be = q_evol[4] |
---|
| 272 | ue = q_evol[5] |
---|
| 273 | ve = q_evol[6] |
---|
| 274 | |
---|
| 275 | |
---|
| 276 | hc = wc - be |
---|
| 277 | |
---|
[7573] | 278 | if hc <= 0.0: |
---|
[7559] | 279 | hc = 0.0 |
---|
| 280 | uc = 0.0 |
---|
| 281 | vc = 0.0 |
---|
| 282 | else: |
---|
| 283 | uc = uhc/hc |
---|
| 284 | vc = vhc/hc |
---|
| 285 | |
---|
| 286 | q_evol[0] = wc |
---|
| 287 | q_evol[1] = uhc |
---|
| 288 | q_evol[2] = vhc |
---|
| 289 | |
---|
| 290 | q_evol[3] = hc |
---|
| 291 | q_evol[4] = be |
---|
| 292 | q_evol[5] = uc |
---|
| 293 | q_evol[6] = vc |
---|
| 294 | |
---|
| 295 | |
---|
| 296 | return q_evol |
---|
| 297 | |
---|
| 298 | |
---|