[6145] | 1 | #!/usr/bin/env python |
---|
[7780] | 2 | """ |
---|
| 3 | Generic boundary conditions for a domain. |
---|
| 4 | |
---|
| 5 | A boundary represents the edge of the model, where inflow, outflow, and |
---|
| 6 | reflection can take place. |
---|
| 7 | |
---|
| 8 | The boundaries in this model can be applied universally across all |
---|
| 9 | domain models, without being tied to a particular implementation. |
---|
| 10 | """ |
---|
[6145] | 11 | |
---|
| 12 | import unittest |
---|
| 13 | from math import sqrt, pi |
---|
| 14 | |
---|
| 15 | from generic_boundary_conditions import * |
---|
[7737] | 16 | from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain |
---|
[6145] | 17 | from anuga.config import epsilon |
---|
| 18 | |
---|
[7276] | 19 | import numpy as num |
---|
[6145] | 20 | |
---|
| 21 | |
---|
| 22 | class Test_Generic_Boundary_Conditions(unittest.TestCase): |
---|
| 23 | def setUp(self): |
---|
| 24 | pass |
---|
| 25 | #print " Setting up" |
---|
| 26 | |
---|
| 27 | def tearDown(self): |
---|
| 28 | pass |
---|
| 29 | #print " Tearing down" |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | def test_generic(self): |
---|
| 33 | b = Boundary() |
---|
| 34 | |
---|
| 35 | try: |
---|
| 36 | b.evaluate() |
---|
| 37 | except: |
---|
| 38 | pass |
---|
| 39 | else: |
---|
| 40 | raise 'Should have raised exception' |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | def test_dirichlet_empty(self): |
---|
| 44 | |
---|
| 45 | try: |
---|
| 46 | Bd = Dirichlet_boundary() |
---|
| 47 | except: |
---|
| 48 | pass |
---|
| 49 | else: |
---|
| 50 | raise 'Should have raised exception' |
---|
| 51 | |
---|
| 52 | def test_dirichlet(self): |
---|
| 53 | x = [3.14,0,0.1] |
---|
| 54 | Bd = Dirichlet_boundary(x) |
---|
| 55 | |
---|
| 56 | q = Bd.evaluate() |
---|
| 57 | assert num.allclose(q, x) |
---|
| 58 | |
---|
| 59 | |
---|
[7505] | 60 | def test_time(self): |
---|
| 61 | from quantity import Quantity |
---|
| 62 | |
---|
| 63 | a = [0.0, 0.0] |
---|
| 64 | b = [0.0, 2.0] |
---|
| 65 | c = [2.0,0.0] |
---|
| 66 | d = [0.0, 4.0] |
---|
| 67 | e = [2.0, 2.0] |
---|
| 68 | f = [4.0,0.0] |
---|
| 69 | |
---|
| 70 | points = [a, b, c, d, e, f] |
---|
| 71 | |
---|
| 72 | #bac, bce, ecf, dbe |
---|
| 73 | elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 74 | |
---|
[7737] | 75 | domain = Generic_Domain(points, elements) |
---|
[7505] | 76 | domain.check_integrity() |
---|
| 77 | |
---|
| 78 | domain.conserved_quantities = ['stage', 'ymomentum'] |
---|
[7519] | 79 | domain.evolved_quantities = ['stage', 'ymomentum'] |
---|
[7505] | 80 | domain.quantities['stage'] =\ |
---|
| 81 | Quantity(domain, [[1,2,3], [5,5,5], |
---|
| 82 | [0,0,9], [-6, 3, 3]]) |
---|
| 83 | |
---|
| 84 | domain.quantities['ymomentum'] =\ |
---|
| 85 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 86 | [0,0,9], [-6, 3, 3]]) |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | domain.check_integrity() |
---|
| 90 | |
---|
| 91 | #Test time bdry, you need to provide a domain and function |
---|
| 92 | try: |
---|
| 93 | T = Time_boundary(domain) |
---|
| 94 | except: |
---|
| 95 | pass |
---|
| 96 | else: |
---|
| 97 | raise 'Should have raised exception' |
---|
| 98 | |
---|
| 99 | #Test time bdry, you need to provide a function |
---|
| 100 | try: |
---|
| 101 | T = Time_boundary() |
---|
| 102 | except: |
---|
| 103 | pass |
---|
| 104 | else: |
---|
| 105 | raise 'Should have raised exception' |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | def function(t): |
---|
| 109 | return [1.0, 0.0] |
---|
| 110 | |
---|
| 111 | T = Time_boundary(domain, function) |
---|
| 112 | |
---|
| 113 | from anuga.config import default_boundary_tag |
---|
| 114 | domain.set_boundary( {default_boundary_tag: T} ) |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | #FIXME: should not necessarily be true always. |
---|
| 118 | #E.g. with None as a boundary object. |
---|
| 119 | assert len(domain.boundary) == len(domain.boundary_objects) |
---|
| 120 | |
---|
| 121 | q = T.evaluate(0, 2) #Vol=0, edge=2 |
---|
| 122 | |
---|
| 123 | assert num.allclose(q, [1.0, 0.0]) |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | def test_time_space_boundary(self): |
---|
| 127 | from quantity import Quantity |
---|
| 128 | |
---|
| 129 | a = [0.0, 0.0] |
---|
| 130 | b = [0.0, 2.0] |
---|
| 131 | c = [2.0,0.0] |
---|
| 132 | d = [0.0, 4.0] |
---|
| 133 | e = [2.0, 2.0] |
---|
| 134 | f = [4.0,0.0] |
---|
| 135 | |
---|
| 136 | points = [a, b, c, d, e, f] |
---|
| 137 | |
---|
| 138 | #bac, bce, ecf, dbe |
---|
| 139 | elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 140 | |
---|
[7737] | 141 | domain = Generic_Domain(points, elements) |
---|
[7505] | 142 | domain.check_integrity() |
---|
| 143 | |
---|
| 144 | domain.conserved_quantities = ['stage', 'ymomentum'] |
---|
[7519] | 145 | domain.evolved_quantities = ['stage', 'ymomentum'] |
---|
[7505] | 146 | domain.quantities['stage'] =\ |
---|
| 147 | Quantity(domain, [[1,2,3], [5,5,5], |
---|
| 148 | [0,0,9], [-6, 3, 3]]) |
---|
| 149 | |
---|
| 150 | domain.quantities['ymomentum'] =\ |
---|
| 151 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 152 | [0,0,9], [-6, 3, 3]]) |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | domain.check_integrity() |
---|
| 156 | |
---|
| 157 | #Test time space bdry, you need to provide a domain and function |
---|
| 158 | try: |
---|
| 159 | T = Time_space_boundary(domain) |
---|
| 160 | except: |
---|
| 161 | pass |
---|
| 162 | else: |
---|
| 163 | raise 'Should have raised exception' |
---|
| 164 | |
---|
| 165 | #Test time bdry, you need to provide a function |
---|
| 166 | try: |
---|
| 167 | T = Time_space_boundary() |
---|
| 168 | except: |
---|
| 169 | pass |
---|
| 170 | else: |
---|
| 171 | raise 'Should have raised exception' |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | def function(t,x,y): |
---|
| 175 | return [x,y] |
---|
| 176 | |
---|
| 177 | T = Time_space_boundary(domain, function) |
---|
| 178 | |
---|
| 179 | from anuga.config import default_boundary_tag |
---|
| 180 | domain.set_boundary( {default_boundary_tag: T} ) |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | #FIXME: should not necessarily be true always. |
---|
| 184 | #E.g. with None as a boundary object. |
---|
| 185 | assert len(domain.boundary) == len(domain.boundary_objects) |
---|
| 186 | |
---|
| 187 | q = T.evaluate(0, 2) #Vol=0, edge=2 |
---|
| 188 | assert num.allclose(q, domain.get_edge_midpoint_coordinate(0,2)) |
---|
| 189 | |
---|
| 190 | |
---|
[7519] | 191 | q = T.evaluate(1, 1) #Vol=1, edge=1 |
---|
[7505] | 192 | assert num.allclose(q, domain.get_edge_midpoint_coordinate(1,1)) |
---|
| 193 | |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | |
---|
| 197 | |
---|
[6145] | 198 | def test_transmissive(self): |
---|
| 199 | from quantity import Quantity |
---|
| 200 | |
---|
| 201 | a = [0.0, 0.0] |
---|
| 202 | b = [0.0, 2.0] |
---|
| 203 | c = [2.0,0.0] |
---|
| 204 | d = [0.0, 4.0] |
---|
| 205 | e = [2.0, 2.0] |
---|
| 206 | f = [4.0,0.0] |
---|
| 207 | |
---|
| 208 | points = [a, b, c, d, e, f] |
---|
| 209 | |
---|
| 210 | #bac, bce, ecf, dbe |
---|
| 211 | elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 212 | |
---|
[7737] | 213 | domain = Generic_Domain(points, elements) |
---|
[6145] | 214 | domain.check_integrity() |
---|
| 215 | |
---|
| 216 | domain.conserved_quantities = ['stage', 'ymomentum'] |
---|
[7519] | 217 | domain.evolved_quantities = ['stage', 'ymomentum'] |
---|
[6145] | 218 | domain.quantities['stage'] =\ |
---|
| 219 | Quantity(domain, [[1,2,3], [5,5,5], |
---|
| 220 | [0,0,9], [-6, 3, 3]]) |
---|
| 221 | |
---|
| 222 | domain.quantities['ymomentum'] =\ |
---|
| 223 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 224 | [0,0,9], [-6, 3, 3]]) |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | domain.check_integrity() |
---|
| 228 | |
---|
| 229 | #Test transmissve bdry |
---|
| 230 | try: |
---|
| 231 | T = Transmissive_boundary() |
---|
| 232 | except: |
---|
| 233 | pass |
---|
| 234 | else: |
---|
| 235 | raise 'Should have raised exception' |
---|
| 236 | |
---|
| 237 | T = Transmissive_boundary(domain) |
---|
| 238 | |
---|
| 239 | from anuga.config import default_boundary_tag |
---|
| 240 | domain.set_boundary( {default_boundary_tag: T} ) |
---|
| 241 | |
---|
| 242 | |
---|
| 243 | #FIXME: should not necessarily be true always. |
---|
| 244 | #E.g. with None as a boundary object. |
---|
| 245 | assert len(domain.boundary) == len(domain.boundary_objects) |
---|
| 246 | |
---|
| 247 | q = T.evaluate(0, 2) #Vol=0, edge=2 |
---|
| 248 | |
---|
| 249 | assert num.allclose(q, [1.5, 2.5]) |
---|
| 250 | |
---|
| 251 | |
---|
[7573] | 252 | # Now set the centroid_transmissive_bc flag to true |
---|
| 253 | domain.set_centroid_transmissive_bc(True) |
---|
| 254 | |
---|
| 255 | q = T.evaluate(0, 2) #Vol=0, edge=2 |
---|
| 256 | |
---|
| 257 | assert num.allclose(q, [2.0 ,3.0]) # centroid value |
---|
| 258 | |
---|
| 259 | |
---|
| 260 | |
---|
| 261 | |
---|
| 262 | |
---|
[6145] | 263 | def NOtest_fileboundary_time_only(self): |
---|
| 264 | """Test that boundary values can be read from file and interpolated |
---|
| 265 | This is using the .tms file format |
---|
| 266 | |
---|
| 267 | See also test_util for comprenhensive testing of the underlying |
---|
| 268 | file_function and also tests in test_datamanager which tests |
---|
| 269 | file_function using the sts format |
---|
| 270 | """ |
---|
| 271 | #FIXME (Ole): This test was disabled 18 August 2008 as no |
---|
| 272 | # need for this was found. Rather I implemented an Exception |
---|
| 273 | # to catch possible errors in the model setup |
---|
| 274 | |
---|
| 275 | from quantity import Quantity |
---|
| 276 | import time, os |
---|
| 277 | from math import sin, pi |
---|
| 278 | from anuga.config import time_format |
---|
| 279 | |
---|
| 280 | a = [0.0, 0.0] |
---|
| 281 | b = [0.0, 2.0] |
---|
| 282 | c = [2.0, 0.0] |
---|
| 283 | d = [0.0, 4.0] |
---|
| 284 | e = [2.0, 2.0] |
---|
| 285 | f = [4.0, 0.0] |
---|
| 286 | |
---|
| 287 | points = [a, b, c, d, e, f] |
---|
| 288 | |
---|
| 289 | #bac, bce, ecf, dbe |
---|
| 290 | elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 291 | |
---|
[7737] | 292 | domain = Generic_Domain(points, elements) |
---|
[6145] | 293 | domain.conserved_quantities = ['stage', 'ymomentum'] |
---|
| 294 | domain.quantities['stage'] =\ |
---|
| 295 | Quantity(domain, [[1,2,3], [5,5,5], |
---|
| 296 | [0,0,9], [-6, 3, 3]]) |
---|
| 297 | |
---|
| 298 | domain.quantities['ymomentum'] =\ |
---|
| 299 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 300 | [0,0,9], [-6, 3, 3]]) |
---|
| 301 | |
---|
| 302 | domain.check_integrity() |
---|
| 303 | |
---|
| 304 | |
---|
| 305 | #Write file |
---|
| 306 | filename = 'boundarytest' + str(time.time()) |
---|
| 307 | fid = open(filename + '.txt', 'w') |
---|
| 308 | start = time.mktime(time.strptime('2000', '%Y')) |
---|
| 309 | dt = 5*60 #Five minute intervals |
---|
| 310 | for i in range(10): |
---|
| 311 | t = start + i*dt |
---|
| 312 | t_string = time.strftime(time_format, time.gmtime(t)) |
---|
| 313 | |
---|
| 314 | fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) |
---|
| 315 | fid.close() |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | #Convert ASCII file to NetCDF (Which is what we really like!) |
---|
| 319 | |
---|
| 320 | from anuga.shallow_water.data_manager import timefile2netcdf |
---|
| 321 | |
---|
| 322 | timefile2netcdf(filename, quantity_names = ['stage', 'ymomentum']) |
---|
| 323 | |
---|
| 324 | |
---|
| 325 | |
---|
| 326 | F = File_boundary(filename + '.tms', domain) |
---|
| 327 | |
---|
| 328 | |
---|
| 329 | os.remove(filename + '.txt') |
---|
| 330 | os.remove(filename + '.tms') |
---|
| 331 | |
---|
| 332 | |
---|
| 333 | |
---|
| 334 | |
---|
| 335 | #Check that midpoint coordinates at boundary are correctly computed |
---|
| 336 | assert num.allclose( F.midpoint_coordinates, |
---|
| 337 | [[1.0, 0.0], [0.0, 1.0], [3.0, 0.0], |
---|
| 338 | [3.0, 1.0], [1.0, 3.0], [0.0, 3.0]]) |
---|
| 339 | |
---|
| 340 | #assert allclose(F.midpoint_coordinates[(3,2)], [0.0, 3.0]) |
---|
| 341 | #assert allclose(F.midpoint_coordinates[(3,1)], [1.0, 3.0]) |
---|
| 342 | #assert allclose(F.midpoint_coordinates[(0,2)], [0.0, 1.0]) |
---|
| 343 | #assert allclose(F.midpoint_coordinates[(0,0)], [1.0, 0.0]) |
---|
| 344 | #assert allclose(F.midpoint_coordinates[(2,0)], [3.0, 0.0]) |
---|
| 345 | #assert allclose(F.midpoint_coordinates[(2,1)], [3.0, 1.0]) |
---|
| 346 | |
---|
| 347 | |
---|
| 348 | #Check time interpolation |
---|
| 349 | from anuga.config import default_boundary_tag |
---|
| 350 | domain.set_boundary( {default_boundary_tag: F} ) |
---|
| 351 | |
---|
| 352 | domain.time = 5*30/2 #A quarter way through first step |
---|
| 353 | q = F.evaluate() |
---|
| 354 | assert num.allclose(q, [1.0/4, sin(2*pi/10)/4]) |
---|
| 355 | |
---|
| 356 | |
---|
| 357 | domain.time = 2.5*5*60 #Half way between steps 2 and 3 |
---|
| 358 | q = F.evaluate() |
---|
| 359 | assert num.allclose(q, [2.5, (sin(2*2*pi/10) + sin(3*2*pi/10))/2]) |
---|
| 360 | |
---|
| 361 | |
---|
| 362 | |
---|
| 363 | def test_fileboundary_exception(self): |
---|
| 364 | """Test that boundary object complains if number of |
---|
| 365 | conserved quantities are wrong |
---|
| 366 | """ |
---|
| 367 | |
---|
| 368 | from quantity import Quantity |
---|
| 369 | import time, os |
---|
| 370 | from math import sin, pi |
---|
| 371 | from anuga.config import time_format |
---|
| 372 | |
---|
| 373 | a = [0.0, 0.0] |
---|
| 374 | b = [0.0, 2.0] |
---|
| 375 | c = [2.0,0.0] |
---|
| 376 | d = [0.0, 4.0] |
---|
| 377 | e = [2.0, 2.0] |
---|
| 378 | f = [4.0,0.0] |
---|
| 379 | |
---|
| 380 | points = [a, b, c, d, e, f] |
---|
| 381 | |
---|
| 382 | #bac, bce, ecf, dbe |
---|
| 383 | elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] |
---|
| 384 | |
---|
[7737] | 385 | domain = Generic_Domain(points, elements) |
---|
[6145] | 386 | domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
[7519] | 387 | domain.evolved_quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
[6145] | 388 | domain.quantities['stage'] =\ |
---|
| 389 | Quantity(domain, [[1,2,3], [5,5,5], |
---|
| 390 | [0,0,9], [-6, 3, 3]]) |
---|
| 391 | |
---|
| 392 | domain.quantities['xmomentum'] =\ |
---|
| 393 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 394 | [0,0,9], [-6, 3, 3]]) |
---|
| 395 | domain.quantities['ymomentum'] =\ |
---|
| 396 | Quantity(domain, [[2,3,4], [5,5,5], |
---|
| 397 | [0,0,9], [-6, 3, 3]]) |
---|
| 398 | |
---|
| 399 | domain.check_integrity() |
---|
| 400 | |
---|
| 401 | #Write file (with only two values) |
---|
| 402 | filename = 'boundarytest' + str(time.time()) |
---|
| 403 | fid = open(filename + '.txt', 'w') |
---|
| 404 | start = time.mktime(time.strptime('2000', '%Y')) |
---|
| 405 | dt = 5*60 #Five minute intervals |
---|
| 406 | for i in range(10): |
---|
| 407 | t = start + i*dt |
---|
| 408 | t_string = time.strftime(time_format, time.gmtime(t)) |
---|
| 409 | |
---|
| 410 | fid.write('%s,%f %f\n' %(t_string, 1.0*i, sin(i*2*pi/10))) |
---|
| 411 | fid.close() |
---|
| 412 | |
---|
| 413 | |
---|
| 414 | #Convert ASCII file to NetCDF (Which is what we really like!) |
---|
[7780] | 415 | from anuga.file_conversion.file_conversion import timefile2netcdf |
---|
[6145] | 416 | |
---|
| 417 | timefile2netcdf(filename, quantity_names = ['stage', 'xmomentum']) |
---|
| 418 | |
---|
| 419 | |
---|
| 420 | try: |
---|
| 421 | F = File_boundary(filename + '.tms', |
---|
| 422 | domain) |
---|
| 423 | except: |
---|
| 424 | pass |
---|
| 425 | else: |
---|
| 426 | raise 'Should have raised an exception' |
---|
| 427 | |
---|
| 428 | os.remove(filename + '.txt') |
---|
| 429 | os.remove(filename + '.tms') |
---|
| 430 | |
---|
| 431 | |
---|
[7276] | 432 | #------------------------------------------------------------- |
---|
[6145] | 433 | |
---|
| 434 | if __name__ == "__main__": |
---|
| 435 | suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test') |
---|
| 436 | runner = unittest.TextTestRunner() |
---|
| 437 | runner.run(suite) |
---|