[7839] | 1 | import quantity |
---|
| 2 | #!/usr/bin/env python |
---|
| 3 | |
---|
| 4 | import unittest |
---|
| 5 | from math import sqrt, pi |
---|
| 6 | |
---|
| 7 | |
---|
[7855] | 8 | from anuga_1d.base.generic_domain import Generic_domain as Domain |
---|
[7839] | 9 | #from shallow_water_domain import flux_function as domain_flux_function |
---|
| 10 | |
---|
[7855] | 11 | from anuga_1d.base.quantity import * |
---|
[7839] | 12 | |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | from numpy import allclose, array, ones, zeros |
---|
| 16 | import numpy |
---|
| 17 | |
---|
| 18 | |
---|
| 19 | class Test_Quantity(unittest.TestCase): |
---|
| 20 | def setUp(self): |
---|
| 21 | self.points3 = [0.0, 1.0, 2.0, 3.0] |
---|
| 22 | self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]] |
---|
| 23 | self.domain3 = Domain(self.points3) |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | self.points4 = [0.0, 1.0, 2.5, 3.0, 5.0] |
---|
| 28 | self.vertex_values4 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]] |
---|
| 29 | self.centroid_values4 = [1.5, 4.5, 0.5, 4.5] |
---|
| 30 | self.boundary4 = {(0, 0): 'left', (3, 1): 'right'} |
---|
| 31 | self.domain4 = Domain(self.points4,self.boundary4) |
---|
| 32 | |
---|
| 33 | self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] |
---|
| 34 | self.domain10 = Domain(self.points10) |
---|
| 35 | |
---|
| 36 | self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] |
---|
| 37 | self.domain6 = Domain(self.points6) |
---|
| 38 | |
---|
| 39 | def tearDown(self): |
---|
| 40 | pass |
---|
| 41 | #print " Tearing down" |
---|
| 42 | |
---|
| 43 | |
---|
| 44 | def test_creat_with_boundary(self): |
---|
| 45 | |
---|
| 46 | assert self.domain4.boundary == {(0, 0): 'left', (3, 1): 'right'} |
---|
[7855] | 47 | |
---|
[7839] | 48 | |
---|
| 49 | def test_creation(self): |
---|
| 50 | |
---|
| 51 | quantity = Quantity(self.domain3) |
---|
| 52 | assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]]) |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | try: |
---|
| 56 | quantity = Quantity() |
---|
| 57 | except: |
---|
| 58 | pass |
---|
| 59 | else: |
---|
| 60 | raise 'Should have raised empty quantity exception' |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | try: |
---|
| 64 | quantity = Quantity([1,2,3]) |
---|
| 65 | except AssertionError: |
---|
| 66 | pass |
---|
[7860] | 67 | else: |
---|
[7839] | 68 | raise 'Should have raised "mising domain object" error' |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | def test_creation_zeros(self): |
---|
| 72 | |
---|
| 73 | quantity = Quantity(self.domain3) |
---|
| 74 | assert allclose(quantity.centroid_values, [[0.,0.,0.]]) |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | quantity = Quantity(self.domain4) |
---|
| 78 | assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.], |
---|
| 79 | [0.,0.], [0.,0.]]) |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | def test_interpolation(self): |
---|
| 83 | quantity = Quantity(self.domain4, self.vertex_values4) |
---|
| 84 | assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | |
---|
| 88 | def test_interpolation2(self): |
---|
| 89 | quantity = Quantity(self.domain4, self.vertex_values4) |
---|
| 90 | assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid |
---|
| 91 | |
---|
| 92 | quantity.extrapolate_second_order() |
---|
| 93 | |
---|
| 94 | #print quantity.vertex_values |
---|
| 95 | assert allclose(quantity.vertex_values,[[ 0.3, 2.7], |
---|
| 96 | [ 4.5, 4.5], |
---|
| 97 | [ 0.5, 0.5], |
---|
| 98 | [ 1.3, 7.7]]) |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | def test_boundary_allocation(self): |
---|
| 104 | quantity = Quantity(self.domain4, |
---|
| 105 | [[1,2], [5,5], [0,9], [-6, 3]]) |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | assert quantity.boundary_values.shape[0] == len(self.domain4.boundary) |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | def test_set_values(self): |
---|
| 112 | quantity = Quantity(self.domain4) |
---|
| 113 | |
---|
| 114 | # Test vertices |
---|
| 115 | quantity.set_values([[1,2], [5,5], [0,9], [-6, 3]], location = 'vertices') |
---|
| 116 | assert allclose(quantity.vertex_values, |
---|
| 117 | [[1,2], [5,5], [0,9], [-6, 3]]) |
---|
| 118 | assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5]) #Centroid |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | # Test unique_vertices |
---|
| 123 | quantity.set_values([1,2,4,-5,6], location='unique_vertices') |
---|
| 124 | assert allclose(quantity.vertex_values, |
---|
| 125 | [[1,2], [2,4], [4,-5], [-5,6]]) |
---|
| 126 | assert allclose(quantity.centroid_values, [1.5, 3., -.5, .5]) #Centroid |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | # Test centroids |
---|
| 130 | quantity.set_values([1,2,3,4], location = 'centroids') |
---|
| 131 | assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid |
---|
| 132 | |
---|
| 133 | # Test exceptions |
---|
| 134 | try: |
---|
| 135 | quantity.set_values([[1,3], [5,5], [0,9], [-6, 3]], |
---|
| 136 | location = 'bas kamel tuba') |
---|
| 137 | except: |
---|
| 138 | pass |
---|
| 139 | |
---|
| 140 | |
---|
| 141 | try: |
---|
| 142 | quantity.set_values([[1,3], [0,9]]) |
---|
| 143 | except AssertionError: |
---|
| 144 | pass |
---|
| 145 | except: |
---|
| 146 | raise 'should have raised Assertionerror' |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | def test_set_values_const(self): |
---|
| 151 | quantity = Quantity(self.domain4) |
---|
| 152 | |
---|
| 153 | quantity.set_values(1.0, location = 'vertices') |
---|
| 154 | |
---|
| 155 | assert allclose(quantity.vertex_values, [[1,1], [1,1], [1,1], [1, 1]]) |
---|
| 156 | assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | quantity.set_values(2.0, location = 'centroids') |
---|
| 160 | |
---|
| 161 | assert allclose(quantity.centroid_values, [2, 2, 2, 2]) |
---|
| 162 | |
---|
| 163 | |
---|
| 164 | def test_set_values_func(self): |
---|
| 165 | quantity = Quantity(self.domain4) |
---|
| 166 | |
---|
| 167 | def f(x): |
---|
| 168 | return x*x |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | |
---|
| 172 | quantity.set_values(f, location = 'vertices') |
---|
| 173 | |
---|
| 174 | |
---|
| 175 | assert allclose(quantity.vertex_values, |
---|
| 176 | [[0,1], [1,6.25], [6.25,9], [9,25]]) |
---|
| 177 | |
---|
| 178 | assert allclose(quantity.centroid_values, |
---|
| 179 | [0.5, 3.625, 7.625, 34*0.5]) |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | quantity.set_values(f, location = 'centroids') |
---|
| 183 | |
---|
| 184 | |
---|
| 185 | assert allclose(quantity.centroid_values, |
---|
| 186 | [0.25, 3.0625, 7.5625, 16.0]) |
---|
| 187 | |
---|
| 188 | |
---|
| 189 | def test_integral(self): |
---|
| 190 | quantity = Quantity(self.domain4) |
---|
| 191 | |
---|
| 192 | #Try constants first |
---|
| 193 | const = 5 |
---|
| 194 | quantity.set_values(const, location = 'vertices') |
---|
| 195 | #print 'Q', quantity.get_integral() |
---|
| 196 | |
---|
| 197 | assert allclose(quantity.get_integral(), self.domain4.get_area() * const) |
---|
| 198 | |
---|
| 199 | #Try with a linear function |
---|
| 200 | def f(x): |
---|
| 201 | return x |
---|
| 202 | |
---|
| 203 | quantity.set_values(f, location = 'vertices') |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | assert allclose (quantity.centroid_values, |
---|
| 207 | [ 0.5, 1.75, 2.75, 4. ]) |
---|
| 208 | |
---|
| 209 | assert allclose (quantity.vertex_values, [[ 0., 1. ], |
---|
| 210 | [ 1., 2.5], |
---|
| 211 | [ 2.5, 3. ], |
---|
| 212 | [ 3., 5. ]]) |
---|
| 213 | |
---|
| 214 | |
---|
| 215 | ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0 |
---|
| 216 | |
---|
| 217 | assert allclose (quantity.get_integral(), ref_integral) |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | |
---|
| 221 | |
---|
| 222 | |
---|
| 223 | def test_set_values_from_quantity(self): |
---|
| 224 | |
---|
| 225 | quantity1 = Quantity(self.domain4) |
---|
| 226 | quantity1.set_values([0,1,2,3,4], location='unique_vertices') |
---|
| 227 | |
---|
| 228 | assert allclose(quantity1.vertex_values, |
---|
| 229 | [[0,1], [1,2], [2,3], [3,4]]) |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | quantity2 = Quantity(self.domain4) |
---|
| 233 | quantity2.set_values(quantity1) |
---|
| 234 | assert allclose(quantity2.vertex_values, |
---|
| 235 | [[0,1], [1,2], [2,3], [3,4]]) |
---|
| 236 | |
---|
| 237 | quantity2.set_values(2*quantity1) |
---|
| 238 | |
---|
| 239 | assert allclose(quantity2.vertex_values, |
---|
| 240 | [[0,2], [2,4], [4,6], [6,8]]) |
---|
| 241 | |
---|
| 242 | quantity2.set_values(2*quantity1 + 3) |
---|
| 243 | assert allclose(quantity2.vertex_values, |
---|
| 244 | [[3,5], [5,7], [7,9], [9,11]]) |
---|
| 245 | |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | def test_overloading(self): |
---|
| 249 | |
---|
| 250 | quantity1 = Quantity(self.domain4) |
---|
| 251 | quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]], |
---|
| 252 | location = 'vertices') |
---|
| 253 | |
---|
| 254 | assert allclose(quantity1.vertex_values, |
---|
| 255 | [[0,1], [1,2], [2,3], [3,4]]) |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | quantity2 = Quantity(self.domain4) |
---|
| 259 | quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]], |
---|
| 260 | location = 'vertices') |
---|
| 261 | |
---|
| 262 | |
---|
| 263 | |
---|
| 264 | quantity3 = Quantity(self.domain4) |
---|
| 265 | quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]], |
---|
| 266 | location = 'vertices') |
---|
| 267 | |
---|
| 268 | |
---|
| 269 | # Negation |
---|
| 270 | Q = -quantity1 |
---|
| 271 | assert allclose(Q.vertex_values, -quantity1.vertex_values) |
---|
| 272 | assert allclose(Q.centroid_values, -quantity1.centroid_values) |
---|
| 273 | |
---|
| 274 | |
---|
| 275 | # Addition |
---|
| 276 | Q = quantity1 + 7 |
---|
| 277 | assert allclose(Q.vertex_values, quantity1.vertex_values + 7) |
---|
| 278 | assert allclose(Q.centroid_values, quantity1.centroid_values + 7) |
---|
| 279 | |
---|
| 280 | Q = 7 + quantity1 |
---|
| 281 | assert allclose(Q.vertex_values, quantity1.vertex_values + 7) |
---|
| 282 | assert allclose(Q.centroid_values, quantity1.centroid_values + 7) |
---|
| 283 | |
---|
| 284 | Q = quantity1 + quantity2 |
---|
| 285 | assert allclose(Q.vertex_values, |
---|
| 286 | quantity1.vertex_values + quantity2.vertex_values) |
---|
| 287 | assert allclose(Q.centroid_values, |
---|
| 288 | quantity1.centroid_values + quantity2.centroid_values) |
---|
| 289 | |
---|
| 290 | Q = quantity1 + quantity2 - 3 |
---|
| 291 | assert allclose(Q.vertex_values, |
---|
| 292 | quantity1.vertex_values + quantity2.vertex_values - 3) |
---|
| 293 | |
---|
| 294 | Q = quantity1 - quantity2 |
---|
| 295 | assert allclose(Q.vertex_values, |
---|
| 296 | quantity1.vertex_values - quantity2.vertex_values) |
---|
| 297 | |
---|
| 298 | #Scaling |
---|
| 299 | Q = quantity1*3 |
---|
| 300 | assert allclose(Q.vertex_values, quantity1.vertex_values*3) |
---|
| 301 | assert allclose(Q.centroid_values, quantity1.centroid_values*3) |
---|
| 302 | |
---|
| 303 | Q = 3*quantity1 |
---|
| 304 | assert allclose(Q.vertex_values, quantity1.vertex_values*3) |
---|
| 305 | |
---|
| 306 | #Multiplication |
---|
| 307 | Q = quantity1 * quantity2 |
---|
| 308 | assert allclose(Q.vertex_values, |
---|
| 309 | quantity1.vertex_values * quantity2.vertex_values) |
---|
| 310 | |
---|
| 311 | #Linear combinations |
---|
| 312 | Q = 4*quantity1 + 2 |
---|
| 313 | assert allclose(Q.vertex_values, |
---|
| 314 | 4*quantity1.vertex_values + 2) |
---|
| 315 | |
---|
| 316 | Q = quantity1*quantity2 + 2 |
---|
| 317 | assert allclose(Q.vertex_values, |
---|
| 318 | quantity1.vertex_values * quantity2.vertex_values + 2) |
---|
| 319 | |
---|
| 320 | Q = quantity1*quantity2 + quantity3 |
---|
| 321 | assert allclose(Q.vertex_values, |
---|
| 322 | quantity1.vertex_values * quantity2.vertex_values + |
---|
| 323 | quantity3.vertex_values) |
---|
| 324 | Q = quantity1*quantity2 + 3*quantity3 |
---|
| 325 | assert allclose(Q.vertex_values, |
---|
| 326 | quantity1.vertex_values * quantity2.vertex_values + |
---|
| 327 | 3*quantity3.vertex_values) |
---|
| 328 | Q = quantity1*quantity2 + 3*quantity3 + 5.0 |
---|
| 329 | assert allclose(Q.vertex_values, |
---|
| 330 | quantity1.vertex_values * quantity2.vertex_values + |
---|
| 331 | 3*quantity3.vertex_values + 5) |
---|
| 332 | |
---|
| 333 | Q = quantity1*quantity2 - quantity3 |
---|
| 334 | assert allclose(Q.vertex_values, |
---|
| 335 | quantity1.vertex_values * quantity2.vertex_values - |
---|
| 336 | quantity3.vertex_values) |
---|
| 337 | Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0 |
---|
| 338 | assert allclose(Q.vertex_values, |
---|
| 339 | 1.5*quantity1.vertex_values * quantity2.vertex_values - |
---|
| 340 | 3*quantity3.vertex_values + 5) |
---|
| 341 | |
---|
| 342 | #Try combining quantities and arrays and scalars |
---|
| 343 | Q = 1.5*quantity1*quantity2.vertex_values -\ |
---|
| 344 | 3*quantity3.vertex_values + 5.0 |
---|
| 345 | assert allclose(Q.vertex_values, |
---|
| 346 | 1.5*quantity1.vertex_values * quantity2.vertex_values - |
---|
| 347 | 3*quantity3.vertex_values + 5) |
---|
| 348 | |
---|
| 349 | |
---|
| 350 | #Powers |
---|
| 351 | Q = quantity1**2 |
---|
| 352 | assert allclose(Q.vertex_values, quantity1.vertex_values**2) |
---|
| 353 | |
---|
| 354 | Q = quantity1**2 +quantity2**2 |
---|
| 355 | assert allclose(Q.vertex_values, |
---|
| 356 | quantity1.vertex_values**2 + \ |
---|
| 357 | quantity2.vertex_values**2) |
---|
| 358 | |
---|
| 359 | Q = (quantity1**2 +quantity2**2)**0.5 |
---|
| 360 | assert allclose(Q.vertex_values, |
---|
| 361 | (quantity1.vertex_values**2 + \ |
---|
| 362 | quantity2.vertex_values**2)**0.5) |
---|
| 363 | |
---|
| 364 | def test_compute_gradient(self): |
---|
| 365 | quantity = Quantity(self.domain6) |
---|
| 366 | |
---|
| 367 | #Set up for a gradient of (2,0) at mid triangle |
---|
| 368 | quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0], |
---|
| 369 | location = 'centroids') |
---|
| 370 | |
---|
| 371 | |
---|
| 372 | #Gradients |
---|
| 373 | quantity.compute_gradients() |
---|
| 374 | |
---|
| 375 | a = quantity.gradients |
---|
| 376 | |
---|
| 377 | assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5]) |
---|
[7868] | 378 | |
---|
| 379 | quantity.beta = 2.0 |
---|
| 380 | quantity.limiter = "minmod_kurganov" |
---|
[7839] | 381 | quantity.extrapolate_second_order() |
---|
| 382 | |
---|
| 383 | |
---|
| 384 | assert allclose(quantity.vertex_values, [[ 1., 3. ], |
---|
| 385 | [ 4., 4. ], |
---|
| 386 | [ 4., 4. ], |
---|
| 387 | [ 4., 6.], |
---|
| 388 | [ 8.25, 11.75], |
---|
| 389 | [ 11., 13. ]]) |
---|
| 390 | |
---|
| 391 | |
---|
| 392 | |
---|
| 393 | def test_second_order_extrapolation2(self): |
---|
| 394 | quantity = Quantity(self.domain4) |
---|
| 395 | |
---|
| 396 | #Set up for a gradient of (3,1), f(x) = 3x+y |
---|
| 397 | quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3], |
---|
| 398 | location = 'centroids') |
---|
| 399 | |
---|
| 400 | #Gradients |
---|
| 401 | quantity.compute_gradients() |
---|
| 402 | |
---|
| 403 | a = quantity.gradients |
---|
| 404 | |
---|
| 405 | |
---|
| 406 | assert allclose(a[1], 2.8) |
---|
| 407 | |
---|
| 408 | #Work out the others |
---|
[7868] | 409 | quantity.beta = 2.0 |
---|
| 410 | quantity.limiter = "minmod_kurganov" |
---|
[7839] | 411 | quantity.extrapolate_second_order() |
---|
| 412 | |
---|
| 413 | |
---|
| 414 | assert allclose(quantity.vertex_values[1,0], 3.33333333) |
---|
| 415 | assert allclose(quantity.vertex_values[1,1], 7.33333333) |
---|
| 416 | |
---|
| 417 | assert allclose(quantity.centroid_values[1], 0.5*(7.33333333+3.33333333) ) |
---|
| 418 | |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | |
---|
| 422 | def test_backup_saxpy_centroid_values(self): |
---|
| 423 | quantity = Quantity(self.domain4) |
---|
| 424 | |
---|
| 425 | #Set up for a gradient of (3,1), f(x) = 3x+y |
---|
| 426 | c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3]) |
---|
| 427 | d_values = array([1.0, 2.0, 3.0, 4.0]) |
---|
| 428 | quantity.set_values(c_values, location = 'centroids') |
---|
| 429 | |
---|
| 430 | #Backup |
---|
| 431 | quantity.backup_centroid_values() |
---|
| 432 | |
---|
| 433 | #print quantity.vertex_values |
---|
| 434 | assert allclose(quantity.centroid_values, quantity.centroid_backup_values) |
---|
| 435 | |
---|
| 436 | |
---|
| 437 | quantity.set_values(d_values, location = 'centroids') |
---|
| 438 | |
---|
| 439 | quantity.saxpy_centroid_values(2.0, 3.0) |
---|
| 440 | |
---|
| 441 | assert allclose(quantity.centroid_values, 2.0*d_values + 3.0*c_values) |
---|
| 442 | |
---|
| 443 | |
---|
| 444 | |
---|
| 445 | def test_first_order_extrapolator(self): |
---|
| 446 | quantity = Quantity(self.domain4) |
---|
| 447 | |
---|
| 448 | centroid_values = array([1.,2.,3.,4.]) |
---|
| 449 | quantity.set_values(centroid_values, location = 'centroids') |
---|
| 450 | assert allclose(quantity.centroid_values, centroid_values) #Centroid |
---|
| 451 | |
---|
| 452 | #Extrapolate |
---|
| 453 | quantity.extrapolate_first_order() |
---|
| 454 | |
---|
| 455 | #Check that gradient is zero |
---|
| 456 | a = quantity.gradients |
---|
| 457 | assert allclose(a, [0,0,0,0]) |
---|
| 458 | |
---|
| 459 | |
---|
| 460 | #Check vertices but not edge values |
---|
| 461 | assert allclose(quantity.vertex_values, |
---|
| 462 | [[1,1], [2,2], [3,3], [4,4]]) |
---|
| 463 | |
---|
| 464 | |
---|
| 465 | def test_second_order_extrapolator(self): |
---|
| 466 | quantity = Quantity(self.domain4) |
---|
| 467 | |
---|
| 468 | #Set up for a gradient of (3,0) at mid triangle |
---|
| 469 | quantity.set_values([2.0, 4.0, 8.0, 2.0], |
---|
| 470 | location = 'centroids') |
---|
| 471 | |
---|
| 472 | |
---|
| 473 | |
---|
| 474 | quantity.extrapolate_second_order() |
---|
| 475 | |
---|
| 476 | |
---|
| 477 | #Assert that quantities are conserved |
---|
| 478 | for k in range(quantity.centroid_values.shape[0]): |
---|
| 479 | assert allclose (quantity.centroid_values[k], |
---|
| 480 | numpy.sum(quantity.vertex_values[k,:])/2.0) |
---|
| 481 | |
---|
| 482 | |
---|
| 483 | def test_limit(self): |
---|
| 484 | quantity = Quantity(self.domain4) |
---|
| 485 | |
---|
| 486 | #Create a deliberate overshoot (e.g. from gradient computation) |
---|
| 487 | quantity.set_values([[0,0], [2,20], [-20,3], [8,3]]) |
---|
| 488 | |
---|
| 489 | #Limit |
---|
| 490 | quantity.limit_minmod() |
---|
| 491 | |
---|
| 492 | |
---|
| 493 | #cells 1 and 2 are local max and min |
---|
| 494 | assert quantity.vertex_values[1][0] == quantity.centroid_values[1] |
---|
| 495 | assert quantity.vertex_values[1][1] == quantity.centroid_values[1] |
---|
| 496 | |
---|
| 497 | assert quantity.vertex_values[2][0] == quantity.centroid_values[2] |
---|
| 498 | assert quantity.vertex_values[2][1] == quantity.centroid_values[2] |
---|
| 499 | |
---|
| 500 | |
---|
| 501 | |
---|
[7868] | 502 | |
---|
| 503 | def test_limit_minmod(self): |
---|
| 504 | quantity = Quantity(self.domain4) |
---|
| 505 | |
---|
| 506 | #Create a deliberate overshoot (e.g. from gradient computation) |
---|
| 507 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 508 | |
---|
| 509 | #Limit |
---|
| 510 | quantity.limiter = 'minmod' |
---|
| 511 | quantity.limit_range() |
---|
| 512 | |
---|
| 513 | |
---|
| 514 | |
---|
| 515 | |
---|
| 516 | assert allclose(quantity.vertex_values, [[ -2.4, 2.4], |
---|
| 517 | [ 2.4, 9.6], |
---|
| 518 | [ 10.6, 11.4], |
---|
| 519 | [ 11.4, 14.6]] ) |
---|
| 520 | |
---|
| 521 | from anuga_1d.base.quantity_ext import limit_minmod_ext |
---|
| 522 | |
---|
| 523 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 524 | |
---|
| 525 | limit_minmod_ext(quantity) |
---|
| 526 | |
---|
| 527 | |
---|
| 528 | |
---|
| 529 | |
---|
| 530 | assert allclose(quantity.vertex_values, [[ -2.4, 2.4], |
---|
| 531 | [ 2.4, 9.6], |
---|
| 532 | [ 10.6, 11.4], |
---|
| 533 | [ 11.4, 14.6]] ) |
---|
| 534 | |
---|
| 535 | |
---|
| 536 | |
---|
| 537 | def test_limit_minmod_kurganov(self): |
---|
| 538 | quantity = Quantity(self.domain4) |
---|
| 539 | |
---|
| 540 | #Create a deliberate overshoot (e.g. from gradient computation) |
---|
| 541 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 542 | |
---|
| 543 | #Limit |
---|
| 544 | quantity.limiter = 'minmod_kurganov' |
---|
| 545 | quantity.beta = 0.5 |
---|
| 546 | quantity.limit_range() |
---|
| 547 | |
---|
| 548 | |
---|
| 549 | #print quantity.vertex_values |
---|
| 550 | |
---|
| 551 | assert allclose(quantity.vertex_values, [[ -2.4, 2.4], |
---|
| 552 | [ 4.2, 7.8], |
---|
| 553 | [ 10.8, 11.2], |
---|
| 554 | [ 11.4, 14.6]]) |
---|
| 555 | |
---|
| 556 | from anuga_1d.base.quantity_ext import limit_minmod_kurganov_ext |
---|
| 557 | |
---|
| 558 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 559 | |
---|
| 560 | limit_minmod_kurganov_ext(quantity) |
---|
| 561 | |
---|
| 562 | |
---|
| 563 | #print quantity.vertex_values |
---|
| 564 | |
---|
| 565 | |
---|
| 566 | def test_limit_vanleer(self): |
---|
| 567 | quantity = Quantity(self.domain4) |
---|
| 568 | |
---|
| 569 | #Create a deliberate overshoot (e.g. from gradient computation) |
---|
| 570 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 571 | |
---|
| 572 | #Limit |
---|
| 573 | quantity.set_limiter('vanleer') |
---|
| 574 | quantity.limit_range() |
---|
| 575 | |
---|
| 576 | |
---|
| 577 | #print quantity.vertex_values |
---|
| 578 | |
---|
| 579 | assert allclose(quantity.vertex_values, [[ -2.4, 2.4 ], |
---|
| 580 | [ 2.32653061, 9.67346939], |
---|
| 581 | [ 10.39393939, 11.60606061], |
---|
| 582 | [ 11.4, 14.6 ]] ) |
---|
| 583 | |
---|
| 584 | from anuga_1d.base.quantity_ext import limit_vanleer_ext |
---|
| 585 | |
---|
| 586 | quantity.set_values([[0,0], [2,10], [9,13], [12,14]]) |
---|
| 587 | |
---|
| 588 | limit_vanleer_ext(quantity) |
---|
| 589 | |
---|
| 590 | |
---|
| 591 | #print quantity.vertex_values |
---|
| 592 | |
---|
| 593 | assert allclose(quantity.vertex_values, [[ -2.4, 2.4 ], |
---|
| 594 | [ 2.32653061, 9.67346939], |
---|
| 595 | [ 10.39393939, 11.60606061], |
---|
| 596 | [ 11.4, 14.6 ]] ) |
---|
| 597 | # |
---|
| 598 | |
---|
[7855] | 599 | def test_find_qmax_qmin(self): |
---|
| 600 | quantity = Quantity(self.domain4) |
---|
| 601 | |
---|
| 602 | |
---|
| 603 | quantity.set_values([1.0, 4.0, 8.0, 2.0], |
---|
| 604 | location = 'centroids') |
---|
| 605 | |
---|
| 606 | |
---|
| 607 | |
---|
| 608 | quantity.find_qmax_qmin() |
---|
| 609 | |
---|
| 610 | |
---|
| 611 | assert allclose(quantity.qmax, [4.0, 8.0, 8.0, 8.0]) |
---|
| 612 | assert allclose(quantity.qmin, [1.0, 1.0, 2.0, 2.0]) |
---|
| 613 | |
---|
| 614 | |
---|
[7839] | 615 | def test_distribute_first_order(self): |
---|
| 616 | quantity = Quantity(self.domain4) |
---|
| 617 | |
---|
| 618 | #Test centroids |
---|
| 619 | centroid_values = array([1.,2.,3.,4.]) |
---|
| 620 | quantity.set_values(centroid_values, location = 'centroids') |
---|
| 621 | assert allclose(quantity.centroid_values, centroid_values) #Centroid |
---|
| 622 | |
---|
| 623 | |
---|
| 624 | #Extrapolate from centroid to vertices and edges |
---|
| 625 | quantity.extrapolate_first_order() |
---|
| 626 | |
---|
| 627 | assert allclose(quantity.vertex_values,[[ 1., 1.], |
---|
| 628 | [ 2., 2.], |
---|
| 629 | [ 3., 3.], |
---|
| 630 | [ 4., 4.]]) |
---|
| 631 | |
---|
| 632 | |
---|
| 633 | |
---|
| 634 | def test_distribute_second_order(self): |
---|
| 635 | quantity = Quantity(self.domain4) |
---|
| 636 | |
---|
| 637 | #Test centroids |
---|
| 638 | centroid_values = array([2.,4.,8.,2.]) |
---|
| 639 | quantity.set_values(centroid_values, location = 'centroids') |
---|
| 640 | assert allclose(quantity.centroid_values, centroid_values) #Centroid |
---|
| 641 | |
---|
| 642 | |
---|
| 643 | #Extrapolate |
---|
[7868] | 644 | quantity.beta = 2.0 |
---|
| 645 | quantity.limiter = "minmod_kurganov" |
---|
[7839] | 646 | quantity.extrapolate_second_order() |
---|
| 647 | |
---|
| 648 | assert allclose(quantity.vertex_values, [[ 1.2, 2.8], |
---|
| 649 | [ 2., 6. ], |
---|
| 650 | [ 8., 8. ], |
---|
| 651 | [ 6.8, -2.8]]) |
---|
| 652 | |
---|
| 653 | |
---|
| 654 | def test_update_explicit(self): |
---|
| 655 | quantity = Quantity(self.domain4) |
---|
| 656 | |
---|
| 657 | #Test centroids |
---|
| 658 | quantity.set_values([1.,2.,3.,4.], location = 'centroids') |
---|
| 659 | assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid |
---|
| 660 | |
---|
| 661 | #Set explicit_update |
---|
| 662 | quantity.explicit_update = array( [1.,1.,1.,1.] ) |
---|
| 663 | |
---|
| 664 | #Update with given timestep |
---|
| 665 | quantity.update(0.1) |
---|
| 666 | |
---|
| 667 | x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] ) |
---|
| 668 | assert allclose( quantity.centroid_values, x) |
---|
| 669 | |
---|
| 670 | def test_update_semi_implicit(self): |
---|
| 671 | quantity = Quantity(self.domain4) |
---|
| 672 | |
---|
| 673 | #Test centroids |
---|
| 674 | quantity.set_values([1.,2.,3.,4.], location = 'centroids') |
---|
| 675 | assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid |
---|
| 676 | |
---|
| 677 | #Set semi implicit update |
---|
| 678 | quantity.semi_implicit_update = array([1.,1.,1.,1.]) |
---|
| 679 | |
---|
| 680 | #Update with given timestep |
---|
| 681 | timestep = 0.1 |
---|
| 682 | quantity.update(timestep) |
---|
| 683 | |
---|
| 684 | sem = array([1.,1.,1.,1.]) |
---|
| 685 | denom = ones(4, numpy.float)-timestep*sem |
---|
| 686 | |
---|
| 687 | x = array([1, 2, 3, 4])/denom |
---|
| 688 | assert allclose( quantity.centroid_values, x) |
---|
| 689 | |
---|
| 690 | |
---|
| 691 | def test_both_updates(self): |
---|
| 692 | quantity = Quantity(self.domain4) |
---|
| 693 | |
---|
| 694 | #Test centroids |
---|
| 695 | centroid_values = array( [1, 2, 3, 4] ) |
---|
| 696 | quantity.set_values(centroid_values, location = 'centroids') |
---|
| 697 | assert allclose(quantity.centroid_values, centroid_values) #Centroid |
---|
| 698 | |
---|
| 699 | #Set explicit_update |
---|
| 700 | explicit_update = array( [4.,3.,2.,1.] ) |
---|
| 701 | quantity.explicit_update[:,] = explicit_update |
---|
| 702 | |
---|
| 703 | #Set semi implicit update |
---|
| 704 | semi_implicit_update = array( [1.,1.,1.,1.] ) |
---|
| 705 | quantity.semi_implicit_update[:,] = semi_implicit_update |
---|
| 706 | |
---|
| 707 | #Update with given timestep |
---|
| 708 | timestep = 0.1 |
---|
| 709 | quantity.update(0.1) |
---|
| 710 | |
---|
| 711 | denom = 1.0-timestep*semi_implicit_update |
---|
| 712 | x = (centroid_values + timestep*explicit_update)/denom |
---|
| 713 | |
---|
| 714 | assert allclose( quantity.centroid_values, x) |
---|
| 715 | |
---|
| 716 | |
---|
| 717 | |
---|
| 718 | #------------------------------------------------------------- |
---|
| 719 | if __name__ == "__main__": |
---|
| 720 | suite = unittest.makeSuite(Test_Quantity, 'test') |
---|
| 721 | runner = unittest.TextTestRunner() |
---|
| 722 | runner.run(suite) |
---|