[3293] | 1 | #!/usr/bin/env python |
---|
| 2 | |
---|
| 3 | import unittest |
---|
| 4 | from math import sqrt, pi |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | from quantity import * |
---|
| 8 | #from config import epsilon |
---|
| 9 | from Numeric import allclose, array, zeros, Float |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | class TestCase(unittest.TestCase): |
---|
| 13 | def setUp(self): |
---|
| 14 | from domain import Domain |
---|
| 15 | |
---|
| 16 | a = 0.0 |
---|
| 17 | b = 1.0 |
---|
| 18 | c = 2.0 |
---|
| 19 | d = 2.5 |
---|
| 20 | e = 3.1 |
---|
| 21 | f = 4.0 |
---|
| 22 | |
---|
| 23 | self.points = [a, b, c, d, e, f] |
---|
| 24 | self.centroids = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2] |
---|
| 25 | self.vertex_values = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]] |
---|
| 26 | self.centroid_values = [[1.5, 2.5, 3.5, 4.75, 5.55]] |
---|
| 27 | |
---|
| 28 | self.domain1 = Domain(self.points[0:2]) |
---|
| 29 | self.domain5 = Domain(self.points) |
---|
| 30 | |
---|
| 31 | def tearDown(self): |
---|
| 32 | pass |
---|
| 33 | #print " Tearing down" |
---|
| 34 | |
---|
| 35 | |
---|
| 36 | def test_creation(self): |
---|
| 37 | |
---|
| 38 | quantity = Quantity(self.domain5, self.vertex_values) |
---|
| 39 | assert allclose(quantity.centroid_values, self.centroid_values) |
---|
| 40 | |
---|
| 41 | try: |
---|
| 42 | quantity = Quantity() |
---|
| 43 | except: |
---|
| 44 | pass |
---|
| 45 | else: |
---|
| 46 | raise 'Should have raised empty quantity exception' |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | try: |
---|
| 50 | quantity = Quantity([1,2]) |
---|
| 51 | except AssertionError: |
---|
| 52 | pass |
---|
| 53 | except: |
---|
| 54 | raise 'Should have raised "missing domain object" error' |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | def test_creation_zeros(self): |
---|
| 58 | |
---|
| 59 | quantity = Quantity(self.domain1) |
---|
| 60 | assert allclose(quantity.vertex_values, [[0.,0.]]) |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | quantity = Quantity(self.domain5) |
---|
| 64 | assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.], |
---|
| 65 | [0.,0.], [0.,0.], |
---|
| 66 | [0.,0.]]) |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | def test_interpolation(self): |
---|
| 70 | quantity = Quantity(self.domain1, [[1,2]]) |
---|
| 71 | assert allclose(quantity.centroid_values, 1.5) #Centroid |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | def test_interpolation2(self): |
---|
| 75 | quantity = Quantity(self.domain5, |
---|
| 76 | [[1,2], [5,5], [0,9], [-6, 3], [3,4]]) |
---|
| 77 | assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid |
---|
| 78 | |
---|
| 79 | def test_set_values(self): |
---|
| 80 | quantity = Quantity(self.domain5) |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]], |
---|
| 84 | location = 'vertices') |
---|
| 85 | assert allclose(quantity.vertex_values, |
---|
| 86 | [[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) |
---|
| 87 | assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid |
---|
| 88 | |
---|
| 89 | #Test default |
---|
| 90 | quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) |
---|
| 91 | assert allclose(quantity.vertex_values, |
---|
| 92 | [[1,2], [5,5], [0,0], [-6, 3], [-2,4]]) |
---|
| 93 | assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid |
---|
| 94 | |
---|
| 95 | #Test centroids |
---|
| 96 | quantity.set_values([1,2,3,4,5], location = 'centroids') |
---|
| 97 | assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid |
---|
| 98 | |
---|
| 99 | #Test exceptions |
---|
| 100 | try: |
---|
| 101 | quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]], |
---|
| 102 | location = 'bas kamel tuba') |
---|
| 103 | except: |
---|
| 104 | pass |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | try: |
---|
| 108 | quantity.set_values([[1,2], [0,0]]) |
---|
| 109 | except AssertionError: |
---|
| 110 | pass |
---|
| 111 | except: |
---|
| 112 | raise 'should have raised AssertionError' |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | def test_set_values_const(self): |
---|
| 116 | quantity = Quantity(self.domain5) |
---|
| 117 | |
---|
| 118 | quantity.set_values(1.0, location = 'vertices') |
---|
| 119 | assert allclose(quantity.vertex_values, |
---|
| 120 | [[1,1], [1,1], [1,1], [1,1], [1,1]]) |
---|
| 121 | assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid |
---|
| 122 | |
---|
| 123 | quantity.set_values(2.0, location = 'centroids') |
---|
| 124 | assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2]) |
---|
| 125 | |
---|
| 126 | def test_set_values_func(self): |
---|
| 127 | quantity = Quantity(self.domain5) |
---|
| 128 | |
---|
| 129 | def f(x): |
---|
| 130 | return x**2 |
---|
| 131 | |
---|
| 132 | quantity.set_values(f, location = 'vertices') |
---|
| 133 | assert allclose(quantity.vertex_values, |
---|
| 134 | [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]]) |
---|
| 135 | assert allclose(quantity.centroid_values, |
---|
| 136 | [0.5, 2.5, 5.125, 7.93, 12.805 ]) |
---|
| 137 | |
---|
| 138 | quantity.set_values(f, location = 'centroids') |
---|
| 139 | assert allclose(quantity.centroid_values, |
---|
| 140 | [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2]) |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | def test_gradient(self): |
---|
| 144 | |
---|
| 145 | quantity = Conserved_quantity(self.domain5) |
---|
| 146 | |
---|
| 147 | #Set up for a gradient of (3,0) at mid triangle |
---|
| 148 | quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0], |
---|
| 149 | location = 'centroids') |
---|
| 150 | |
---|
| 151 | #Gradients |
---|
| 152 | quantity.compute_gradients() |
---|
| 153 | N = quantity.gradients.shape[0] |
---|
| 154 | |
---|
| 155 | G = quantity.gradients |
---|
| 156 | G0 = zeros(N, Float) |
---|
| 157 | Q = quantity.centroid_values |
---|
| 158 | X = quantity.domain.centroids |
---|
| 159 | |
---|
| 160 | def grad0(x0,x1,x2,q0,q1,q2): |
---|
| 161 | return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0) |
---|
| 162 | |
---|
| 163 | def grad1(x0,x1,q0,q1): |
---|
| 164 | return (q1 - q0)/(x1 - x0) |
---|
| 165 | |
---|
| 166 | #Check Gradients |
---|
| 167 | G0[0] = grad1(X[0],X[1],Q[0],Q[1]) |
---|
| 168 | for i in range(1,4): |
---|
| 169 | G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1]) |
---|
| 170 | G0[4] = grad1(X[3],X[4],Q[3],Q[4]) |
---|
| 171 | |
---|
| 172 | assert allclose(G,G0) |
---|
| 173 | |
---|
| 174 | quantity.extrapolate_first_order() |
---|
| 175 | |
---|
| 176 | assert allclose(quantity.vertex_values, [[2., 2.], |
---|
| 177 | [4., 4.], |
---|
| 178 | [8., 8.], |
---|
| 179 | [2., 2.], |
---|
| 180 | [5., 5.]]) |
---|
| 181 | |
---|
| 182 | quantity.extrapolate_second_order() |
---|
| 183 | |
---|
| 184 | V = quantity.domain.vertices |
---|
| 185 | for k in range(quantity.domain.number_of_elements): |
---|
| 186 | G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5 |
---|
| 187 | |
---|
| 188 | assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]], |
---|
| 189 | [4.-G0[1], 4.+G0[1]], |
---|
| 190 | [8.-G0[2], 8.+G0[2]], |
---|
| 191 | [2.-G0[3], 2.+G0[3]], |
---|
| 192 | [5.-G0[4], 5.+G0[4]]]) |
---|
| 193 | |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | def test_second_order_extrapolation2(self): |
---|
| 197 | quantity = Conserved_quantity(self.domain5) |
---|
| 198 | |
---|
| 199 | #Set up for a gradient of (2), f(x) = 2x |
---|
| 200 | quantity.set_values([1., 3., 4.5, 5.6, 7.1], |
---|
| 201 | location = 'centroids') |
---|
| 202 | |
---|
| 203 | #Gradients |
---|
| 204 | quantity.compute_gradients() |
---|
| 205 | G = quantity.gradients |
---|
| 206 | |
---|
| 207 | allclose(G, 2.0) |
---|
| 208 | |
---|
| 209 | quantity.extrapolate_second_order() |
---|
| 210 | |
---|
| 211 | V = quantity.domain.vertices |
---|
| 212 | Q = quantity.vertex_values |
---|
| 213 | |
---|
| 214 | assert allclose(Q[:,0], 2.0*V[:,0]) |
---|
| 215 | assert allclose(Q[:,1], 2.0*V[:,1]) |
---|
| 216 | |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | |
---|
| 221 | |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | ## def test_first_order_extrapolator(self): |
---|
| 225 | ## quantity = Conserved_quantity(self.mesh4) |
---|
| 226 | |
---|
| 227 | ## #Test centroids |
---|
| 228 | ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') |
---|
| 229 | ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid |
---|
| 230 | |
---|
| 231 | ## #Extrapolate |
---|
| 232 | ## quantity.extrapolate_first_order() |
---|
| 233 | |
---|
| 234 | ## #Check vertices but not edge values |
---|
| 235 | ## assert allclose(quantity.vertex_values, |
---|
| 236 | ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | def test_second_order_limit_extrapolator(self): |
---|
| 240 | quantity = Conserved_quantity(self.domain5) |
---|
| 241 | |
---|
| 242 | #Set up for a gradient of (3,0) at mid triangle |
---|
| 243 | quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0], |
---|
| 244 | location = 'centroids') |
---|
| 245 | |
---|
| 246 | quantity.extrapolate_second_order() |
---|
| 247 | quantity.limit() |
---|
| 248 | |
---|
| 249 | #Assert that central triangle is limited by neighbours |
---|
| 250 | assert quantity.vertex_values[0,0] >= 2.0 |
---|
| 251 | assert quantity.vertex_values[0,0] <= 4.0 |
---|
| 252 | assert quantity.vertex_values[0,1] >= 2.0 |
---|
| 253 | assert quantity.vertex_values[0,1] <= 4.0 |
---|
| 254 | |
---|
| 255 | assert quantity.vertex_values[1,0] >= 2.0 |
---|
| 256 | assert quantity.vertex_values[1,0] <= 8.0 |
---|
| 257 | assert quantity.vertex_values[1,1] >= 2.0 |
---|
| 258 | assert quantity.vertex_values[1,1] <= 8.0 |
---|
| 259 | |
---|
| 260 | assert quantity.vertex_values[2,0] >= 2.0 |
---|
| 261 | assert quantity.vertex_values[2,0] <= 8.0 |
---|
| 262 | assert quantity.vertex_values[2,1] >= 2.0 |
---|
| 263 | assert quantity.vertex_values[2,1] <= 8.0 |
---|
| 264 | |
---|
| 265 | assert quantity.vertex_values[3,0] >= 0.0 |
---|
| 266 | assert quantity.vertex_values[3,0] <= 8.0 |
---|
| 267 | assert quantity.vertex_values[3,1] >= 0.0 |
---|
| 268 | assert quantity.vertex_values[3,1] <= 8.0 |
---|
| 269 | |
---|
| 270 | assert quantity.vertex_values[4,0] >= 0.0 |
---|
| 271 | assert quantity.vertex_values[4,0] <= 2.0 |
---|
| 272 | assert quantity.vertex_values[4,1] >= 0.0 |
---|
| 273 | assert quantity.vertex_values[4,1] <= 2.0 |
---|
| 274 | |
---|
| 275 | #Assert that quantities are conserved |
---|
| 276 | from Numeric import sum |
---|
| 277 | for k in range(quantity.centroid_values.shape[0]): |
---|
| 278 | assert allclose (quantity.centroid_values[k], |
---|
| 279 | sum(quantity.vertex_values[k,:])/2) |
---|
| 280 | |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | def test_limiter(self): |
---|
| 285 | quantity = Conserved_quantity(self.domain5) |
---|
| 286 | |
---|
| 287 | #Create a deliberate overshoot (e.g. from gradient computation) |
---|
| 288 | quantity.set_values([[3,4], [5,5], [0,0], [-6, 3], [-2,4]], |
---|
| 289 | location = 'vertices') |
---|
| 290 | |
---|
| 291 | #Limit |
---|
| 292 | quantity.limit() |
---|
| 293 | |
---|
| 294 | #Assert that central triangle is limited by neighbours |
---|
| 295 | assert quantity.vertex_values[0,1] >= quantity.centroid_values[0] |
---|
| 296 | assert quantity.vertex_values[1,0] <= quantity.centroid_values[1] |
---|
| 297 | |
---|
| 298 | assert quantity.vertex_values[1,1] >= quantity.centroid_values[1] |
---|
| 299 | |
---|
| 300 | |
---|
| 301 | #Assert that quantities are conserved |
---|
| 302 | from Numeric import sum |
---|
| 303 | for k in range(quantity.centroid_values.shape[0]): |
---|
| 304 | assert allclose (quantity.centroid_values[k], |
---|
| 305 | sum(quantity.vertex_values[k,:])/2) |
---|
| 306 | |
---|
| 307 | |
---|
| 308 | |
---|
| 309 | ## def test_distribute_first_order(self): |
---|
| 310 | ## quantity = Conserved_quantity(self.mesh4) |
---|
| 311 | |
---|
| 312 | ## #Test centroids |
---|
| 313 | ## quantity.set_values([1.,2.,3.,4.], location = 'centroids') |
---|
| 314 | ## assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid |
---|
| 315 | |
---|
| 316 | |
---|
| 317 | ## #Extrapolate |
---|
| 318 | ## quantity.extrapolate_first_order() |
---|
| 319 | |
---|
| 320 | ## #Interpolate |
---|
| 321 | ## quantity.interpolate_from_vertices_to_edges() |
---|
| 322 | |
---|
| 323 | ## assert allclose(quantity.vertex_values, |
---|
| 324 | ## [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) |
---|
| 325 | ## assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], |
---|
| 326 | ## [3,3,3], [4, 4, 4]]) |
---|
| 327 | |
---|
| 328 | |
---|
| 329 | |
---|
| 330 | def test_update_explicit(self): |
---|
| 331 | quantity = Conserved_quantity(self.domain5) |
---|
| 332 | |
---|
| 333 | #Test centroids |
---|
| 334 | quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids') |
---|
| 335 | assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid |
---|
| 336 | |
---|
| 337 | #Set explicit_update |
---|
| 338 | quantity.explicit_update = array( [1.,1.,1.,1.,1.] ) |
---|
| 339 | |
---|
| 340 | #Update with given timestep |
---|
| 341 | quantity.update(0.1) |
---|
| 342 | |
---|
| 343 | x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] ) |
---|
| 344 | assert allclose( quantity.centroid_values, x) |
---|
| 345 | |
---|
| 346 | def test_update_semi_implicit(self): |
---|
| 347 | quantity = Conserved_quantity(self.domain5) |
---|
| 348 | |
---|
| 349 | #Test centroids |
---|
| 350 | quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids') |
---|
| 351 | assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid |
---|
| 352 | |
---|
| 353 | #Set semi implicit update |
---|
| 354 | quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] ) |
---|
| 355 | |
---|
| 356 | #Update with given timestep |
---|
| 357 | quantity.update(0.1) |
---|
| 358 | |
---|
| 359 | x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] ) |
---|
| 360 | assert allclose( quantity.centroid_values, x) |
---|
| 361 | |
---|
| 362 | def test_both_updates(self): |
---|
| 363 | quantity = Conserved_quantity(self.domain5) |
---|
| 364 | |
---|
| 365 | #Test centroids |
---|
| 366 | quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids') |
---|
| 367 | assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid |
---|
| 368 | |
---|
| 369 | #Set explicit_update |
---|
| 370 | quantity.explicit_update = array( [4.,3.,2.,1.,0.0] ) |
---|
| 371 | |
---|
| 372 | #Set semi implicit update |
---|
| 373 | quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] ) |
---|
| 374 | |
---|
| 375 | #Update with given timestep |
---|
| 376 | quantity.update(0.1) |
---|
| 377 | |
---|
| 378 | x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] ) |
---|
| 379 | x /= array( [.9,.9,.9,.9,.9] ) |
---|
| 380 | assert allclose( quantity.centroid_values, x) |
---|
| 381 | |
---|
| 382 | |
---|
| 383 | |
---|
| 384 | #------------------------------------------------------------- |
---|
| 385 | if __name__ == "__main__": |
---|
| 386 | suite = unittest.makeSuite(TestCase,'test') |
---|
| 387 | runner = unittest.TextTestRunner(verbosity=2) |
---|
| 388 | runner.run(suite) |
---|
| 389 | |
---|
| 390 | |
---|