[7559] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | import unittest, os |
---|
| 4 | import os.path |
---|
| 5 | from math import pi, sqrt |
---|
| 6 | import tempfile |
---|
| 7 | |
---|
| 8 | from anuga.config import g, epsilon |
---|
| 9 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
| 10 | from anuga.utilities.numerical_tools import mean |
---|
| 11 | from anuga.utilities.polygon import is_inside_polygon |
---|
| 12 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
| 13 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
| 14 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
| 15 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
| 16 | |
---|
| 17 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
| 18 | from swb_domain import * |
---|
| 19 | |
---|
| 20 | import numpy as num |
---|
| 21 | |
---|
| 22 | # Get gateway to C implementation of flux function for direct testing |
---|
| 23 | from shallow_water_ext import flux_function_central as flux_function |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | class Test_swb_clean(unittest.TestCase): |
---|
| 29 | def setUp(self): |
---|
| 30 | pass |
---|
| 31 | |
---|
| 32 | def tearDown(self): |
---|
| 33 | pass |
---|
| 34 | |
---|
| 35 | |
---|
| 36 | def test_balance_deep_and_shallow(self): |
---|
| 37 | """Test that balanced limiters preserve conserved quantites. |
---|
| 38 | This test is using old depth based balanced limiters |
---|
| 39 | """ |
---|
| 40 | |
---|
| 41 | import copy |
---|
| 42 | |
---|
| 43 | a = [0.0, 0.0] |
---|
| 44 | b = [0.0, 2.0] |
---|
| 45 | c = [2.0, 0.0] |
---|
| 46 | d = [0.0, 4.0] |
---|
| 47 | e = [2.0, 2.0] |
---|
| 48 | f = [4.0, 0.0] |
---|
| 49 | |
---|
| 50 | points = [a, b, c, d, e, f] |
---|
| 51 | # bac, bce, ecf, dbe |
---|
| 52 | elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 53 | |
---|
| 54 | domain = Domain(points, elements) |
---|
| 55 | domain.check_integrity() |
---|
| 56 | |
---|
| 57 | # Create a deliberate overshoot |
---|
| 58 | domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) |
---|
| 59 | domain.set_quantity('elevation', 0) # Flat bed |
---|
| 60 | stage = domain.quantities['stage'] |
---|
| 61 | |
---|
| 62 | ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy |
---|
| 63 | |
---|
| 64 | # Limit |
---|
| 65 | domain.tight_slope_limiters = 0 |
---|
| 66 | domain.distribute_to_vertices_and_edges() |
---|
| 67 | |
---|
| 68 | # Assert that quantities are conserved |
---|
| 69 | for k in range(len(domain)): |
---|
| 70 | assert num.allclose(ref_centroid_values[k], |
---|
| 71 | num.sum(stage.vertex_values[k,:])/3) |
---|
| 72 | |
---|
| 73 | # Now try with a non-flat bed - closely hugging initial stage in places |
---|
| 74 | # This will create alphas in the range [0, 0.478260, 1] |
---|
| 75 | domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) |
---|
| 76 | domain.set_quantity('elevation', [[0,0,0], |
---|
| 77 | [1.8,1.9,5.9], |
---|
| 78 | [4.6,0,0], |
---|
| 79 | [0,2,4]]) |
---|
| 80 | stage = domain.quantities['stage'] |
---|
| 81 | |
---|
| 82 | ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy |
---|
| 83 | ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy |
---|
| 84 | |
---|
| 85 | # Limit |
---|
| 86 | domain.tight_slope_limiters = 0 |
---|
| 87 | domain.distribute_to_vertices_and_edges() |
---|
| 88 | |
---|
| 89 | # Assert that all vertex quantities have changed |
---|
| 90 | for k in range(len(domain)): |
---|
| 91 | assert not num.allclose(ref_vertex_values[k,:], |
---|
| 92 | stage.vertex_values[k,:]) |
---|
| 93 | # and assert that quantities are still conserved |
---|
| 94 | for k in range(len(domain)): |
---|
| 95 | assert num.allclose(ref_centroid_values[k], |
---|
| 96 | num.sum(stage.vertex_values[k,:])/3) |
---|
| 97 | |
---|
| 98 | # Check actual results |
---|
| 99 | assert num.allclose(stage.vertex_values, |
---|
[7573] | 100 | [[ 2., 2., 2. ], |
---|
| 101 | [ 1.93333333, 2.03333333, 6.03333333], |
---|
| 102 | [ 8.4, 3.8, 3.8 ], |
---|
| 103 | [ 3.33333333, 5.33333333, 7.33333333]]) or \ |
---|
| 104 | num.allclose(stage.vertex_values, |
---|
| 105 | [[ 1.06666667, 3.86666667, 1.06666667], |
---|
| 106 | [ 1.93333333, 2.03333333, 6.03333333], |
---|
| 107 | [ 5.46666667, 3.06666667, 7.46666667], |
---|
| 108 | [ 6.53333333, 4.69333333, 4.77333333]]) |
---|
[7559] | 109 | |
---|
[7573] | 110 | |
---|
| 111 | |
---|
[7559] | 112 | def test_balance_deep_and_shallow_tight_SL(self): |
---|
| 113 | """Test that balanced limiters preserve conserved quantites. |
---|
| 114 | This test is using Tight Slope Limiters |
---|
| 115 | """ |
---|
| 116 | |
---|
| 117 | import copy |
---|
| 118 | |
---|
| 119 | a = [0.0, 0.0] |
---|
| 120 | b = [0.0, 2.0] |
---|
| 121 | c = [2.0, 0.0] |
---|
| 122 | d = [0.0, 4.0] |
---|
| 123 | e = [2.0, 2.0] |
---|
| 124 | f = [4.0, 0.0] |
---|
| 125 | |
---|
| 126 | points = [a, b, c, d, e, f] |
---|
| 127 | # bac, bce, ecf, dbe |
---|
| 128 | elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 129 | |
---|
| 130 | domain = Domain(points, elements) |
---|
| 131 | domain.check_integrity() |
---|
| 132 | |
---|
| 133 | # Create a deliberate overshoot |
---|
| 134 | domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) |
---|
| 135 | domain.set_quantity('elevation', 0) # Flat bed |
---|
| 136 | stage = domain.quantities['stage'] |
---|
| 137 | |
---|
| 138 | ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy |
---|
| 139 | |
---|
| 140 | # Limit |
---|
| 141 | domain.tight_slope_limiters = 1 |
---|
| 142 | domain.distribute_to_vertices_and_edges() |
---|
| 143 | |
---|
| 144 | # Assert that quantities are conserved |
---|
| 145 | for k in range(len(domain)): |
---|
| 146 | assert num.allclose (ref_centroid_values[k], |
---|
| 147 | num.sum(stage.vertex_values[k,:])/3) |
---|
| 148 | |
---|
| 149 | # Now try with a non-flat bed - closely hugging initial stage in places |
---|
| 150 | # This will create alphas in the range [0, 0.478260, 1] |
---|
| 151 | domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) |
---|
| 152 | domain.set_quantity('elevation', [[0,0,0], |
---|
| 153 | [1.8,1.9,5.9], |
---|
| 154 | [4.6,0,0], |
---|
| 155 | [0,2,4]]) |
---|
| 156 | stage = domain.quantities['stage'] |
---|
| 157 | |
---|
| 158 | ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy |
---|
| 159 | ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy |
---|
| 160 | |
---|
| 161 | # Limit |
---|
| 162 | domain.tight_slope_limiters = 1 |
---|
| 163 | domain.distribute_to_vertices_and_edges() |
---|
| 164 | |
---|
| 165 | # Assert that all vertex quantities have changed |
---|
| 166 | for k in range(len(domain)): |
---|
| 167 | assert not num.allclose(ref_vertex_values[k,:], |
---|
| 168 | stage.vertex_values[k,:]) |
---|
| 169 | # and assert that quantities are still conserved |
---|
| 170 | for k in range(len(domain)): |
---|
| 171 | assert num.allclose(ref_centroid_values[k], |
---|
| 172 | num.sum(stage.vertex_values[k,:])/3) |
---|
| 173 | |
---|
| 174 | def test_balance_deep_and_shallow_Froude(self): |
---|
| 175 | """Test that balanced limiters preserve conserved quantites - |
---|
| 176 | and also that excessive Froude numbers are dealt with. |
---|
| 177 | This test is using tight slope limiters. |
---|
| 178 | """ |
---|
| 179 | |
---|
| 180 | import copy |
---|
| 181 | |
---|
| 182 | a = [0.0, 0.0] |
---|
| 183 | b = [0.0, 2.0] |
---|
| 184 | c = [2.0, 0.0] |
---|
| 185 | d = [0.0, 4.0] |
---|
| 186 | e = [2.0, 2.0] |
---|
| 187 | f = [4.0, 0.0] |
---|
| 188 | |
---|
| 189 | points = [a, b, c, d, e, f] |
---|
| 190 | # bac, bce, ecf, dbe |
---|
| 191 | elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 192 | |
---|
| 193 | domain = Domain(points, elements) |
---|
| 194 | domain.check_integrity() |
---|
| 195 | domain.tight_slope_limiters = True |
---|
| 196 | domain.use_centroid_velocities = True |
---|
| 197 | |
---|
| 198 | # Create non-flat bed - closely hugging initial stage in places |
---|
| 199 | # This will create alphas in the range [0, 0.478260, 1] |
---|
| 200 | domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) |
---|
| 201 | domain.set_quantity('elevation', [[0,0,0], |
---|
| 202 | [1.8,1.999,5.999], |
---|
| 203 | [4.6,0,0], |
---|
| 204 | [0,2,4]]) |
---|
| 205 | |
---|
| 206 | # Create small momenta, that nonetheless will generate large speeds |
---|
| 207 | # due to shallow depth at isolated vertices |
---|
| 208 | domain.set_quantity('xmomentum', -0.0058) |
---|
| 209 | domain.set_quantity('ymomentum', 0.0890) |
---|
| 210 | |
---|
| 211 | stage = domain.quantities['stage'] |
---|
| 212 | elevation = domain.quantities['elevation'] |
---|
| 213 | xmomentum = domain.quantities['xmomentum'] |
---|
| 214 | ymomentum = domain.quantities['ymomentum'] |
---|
| 215 | |
---|
| 216 | # Setup triangle #1 to mimick real Froude explosion observed |
---|
| 217 | # in the Onslow example 13 Nov 2007. |
---|
| 218 | stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953] |
---|
| 219 | elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647] |
---|
| 220 | xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066] |
---|
| 221 | ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890] |
---|
| 222 | |
---|
| 223 | xmomentum.interpolate() |
---|
| 224 | ymomentum.interpolate() |
---|
| 225 | stage.interpolate() |
---|
| 226 | elevation.interpolate() |
---|
| 227 | |
---|
| 228 | # Verify interpolation |
---|
| 229 | assert num.allclose(stage.centroid_values[1], 1.5233) |
---|
| 230 | assert num.allclose(elevation.centroid_values[1], 1.2452667) |
---|
| 231 | assert num.allclose(xmomentum.centroid_values[1], -0.0058) |
---|
| 232 | assert num.allclose(ymomentum.centroid_values[1], 0.089) |
---|
| 233 | |
---|
| 234 | # Derived quantities |
---|
| 235 | depth = stage-elevation |
---|
| 236 | u = xmomentum/depth |
---|
| 237 | v = ymomentum/depth |
---|
| 238 | |
---|
| 239 | denom = (depth*g)**0.5 |
---|
| 240 | Fx = u/denom |
---|
| 241 | Fy = v/denom |
---|
| 242 | |
---|
| 243 | # Verify against Onslow example (14 Nov 2007) |
---|
| 244 | assert num.allclose(depth.centroid_values[1], 0.278033) |
---|
| 245 | assert num.allclose(u.centroid_values[1], -0.0208608) |
---|
| 246 | assert num.allclose(v.centroid_values[1], 0.3201055) |
---|
| 247 | |
---|
| 248 | assert num.allclose(denom.centroid_values[1], |
---|
| 249 | num.sqrt(depth.centroid_values[1]*g)) |
---|
| 250 | |
---|
| 251 | assert num.allclose(u.centroid_values[1]/denom.centroid_values[1], |
---|
| 252 | -0.012637746977) |
---|
| 253 | assert num.allclose(Fx.centroid_values[1], |
---|
| 254 | u.centroid_values[1]/denom.centroid_values[1]) |
---|
| 255 | |
---|
| 256 | # Check that Froude numbers are small at centroids. |
---|
| 257 | assert num.allclose(Fx.centroid_values[1], -0.012637746977) |
---|
| 258 | assert num.allclose(Fy.centroid_values[1], 0.193924048435) |
---|
| 259 | |
---|
| 260 | # But Froude numbers are huge at some vertices and edges |
---|
| 261 | assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01, |
---|
| 262 | -1.27775313e+01, |
---|
| 263 | -2.78511420e-03]) |
---|
| 264 | |
---|
| 265 | assert num.allclose(Fx.edge_values[1,:], [-6.89150773e-03, |
---|
| 266 | -7.38672488e-03, |
---|
| 267 | -2.35626238e+01]) |
---|
| 268 | |
---|
| 269 | assert num.allclose(Fy.vertex_values[1,:], [8.99035764e+02, |
---|
| 270 | 2.27440057e+02, |
---|
| 271 | 3.75568430e-02]) |
---|
| 272 | |
---|
| 273 | assert num.allclose(Fy.edge_values[1,:], [1.05748998e-01, |
---|
| 274 | 1.06035244e-01, |
---|
| 275 | 3.88346947e+02]) |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | # The task is now to arrange the limiters such that Froude numbers |
---|
| 279 | # remain under control whil at the same time obeying the conservation |
---|
| 280 | # laws. |
---|
| 281 | ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy |
---|
| 282 | ref_vertex_values = copy.copy(stage.vertex_values[:]) # Copy |
---|
| 283 | |
---|
| 284 | # Limit (and invoke balance_deep_and_shallow) |
---|
| 285 | domain.tight_slope_limiters = 1 |
---|
| 286 | domain.distribute_to_vertices_and_edges() |
---|
| 287 | |
---|
| 288 | # Redo derived quantities |
---|
| 289 | depth = stage - elevation |
---|
| 290 | u = xmomentum/depth |
---|
| 291 | v = ymomentum/depth |
---|
| 292 | |
---|
| 293 | # Assert that all vertex velocities stay within one |
---|
| 294 | # order of magnitude of centroid velocities. |
---|
| 295 | assert num.alltrue(num.absolute(u.vertex_values[1,:]) <= |
---|
| 296 | num.absolute(u.centroid_values[1])*10) |
---|
| 297 | assert num.alltrue(num.absolute(v.vertex_values[1,:]) <= |
---|
| 298 | num.absolute(v.centroid_values[1])*10) |
---|
| 299 | |
---|
| 300 | denom = (depth*g)**0.5 |
---|
| 301 | Fx = u/denom |
---|
| 302 | Fy = v/denom |
---|
| 303 | |
---|
| 304 | # Assert that Froude numbers are less than max value (TBA) |
---|
| 305 | # at vertices, edges and centroids. |
---|
| 306 | from anuga.config import maximum_froude_number |
---|
| 307 | |
---|
| 308 | assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) < |
---|
| 309 | maximum_froude_number) |
---|
| 310 | assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) < |
---|
| 311 | maximum_froude_number) |
---|
| 312 | |
---|
| 313 | # Assert that all vertex quantities have changed |
---|
| 314 | for k in range(len(domain)): |
---|
| 315 | assert not num.allclose(ref_vertex_values[k,:], |
---|
| 316 | stage.vertex_values[k,:]) |
---|
| 317 | |
---|
| 318 | # Assert that quantities are still conserved |
---|
| 319 | for k in range(len(domain)): |
---|
| 320 | assert num.allclose(ref_centroid_values[k], |
---|
| 321 | num.sum(stage.vertex_values[k,:])/3) |
---|
| 322 | |
---|
| 323 | return |
---|
| 324 | |
---|
| 325 | qwidth = 12 |
---|
| 326 | for k in [1]: # range(len(domain)): |
---|
| 327 | print 'Triangle %d (C, V, E)' % k |
---|
| 328 | |
---|
| 329 | print ('stage'.ljust(qwidth), stage.centroid_values[k], |
---|
| 330 | stage.vertex_values[k,:], stage.edge_values[k,:]) |
---|
| 331 | print ('elevation'.ljust(qwidth), elevation.centroid_values[k], |
---|
| 332 | elevation.vertex_values[k,:], elevation.edge_values[k,:]) |
---|
| 333 | print ('depth'.ljust(qwidth), depth.centroid_values[k], |
---|
| 334 | depth.vertex_values[k,:], depth.edge_values[k,:]) |
---|
| 335 | print ('xmomentum'.ljust(qwidth), xmomentum.centroid_values[k], |
---|
| 336 | xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]) |
---|
| 337 | print ('ymomentum'.ljust(qwidth), ymomentum.centroid_values[k], |
---|
| 338 | ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]) |
---|
| 339 | print ('u'.ljust(qwidth),u.centroid_values[k], |
---|
| 340 | u.vertex_values[k,:], u.edge_values[k,:]) |
---|
| 341 | print ('v'.ljust(qwidth), v.centroid_values[k], |
---|
| 342 | v.vertex_values[k,:], v.edge_values[k,:]) |
---|
| 343 | print ('Fx'.ljust(qwidth), Fx.centroid_values[k], |
---|
| 344 | Fx.vertex_values[k,:], Fx.edge_values[k,:]) |
---|
| 345 | print ('Fy'.ljust(qwidth), Fy.centroid_values[k], |
---|
| 346 | Fy.vertex_values[k,:], Fy.edge_values[k,:]) |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | ################################################################################# |
---|
| 350 | |
---|
| 351 | if __name__ == "__main__": |
---|
| 352 | suite = unittest.makeSuite(Test_swb_clean, 'test') |
---|
| 353 | runner = unittest.TextTestRunner(verbosity=1) |
---|
| 354 | runner.run(suite) |
---|