[229] | 1 | """Class Quantity - Implements values at each triangular element |
---|
| 2 | |
---|
| 3 | To create: |
---|
| 4 | |
---|
| 5 | Quantity(domain, vertex_values) |
---|
| 6 | |
---|
| 7 | domain: Associated domain structure. Required. |
---|
| 8 | |
---|
| 9 | vertex_values: N x 3 array of values at each vertex for each element. |
---|
| 10 | Default None |
---|
| 11 | |
---|
| 12 | If vertex_values are None Create array of zeros compatible with domain. |
---|
| 13 | Otherwise check that it is compatible with dimenions of domain. |
---|
| 14 | Otherwise raise an exception |
---|
| 15 | """ |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | class Quantity: |
---|
| 19 | |
---|
| 20 | def __init__(self, domain, vertex_values=None): |
---|
| 21 | |
---|
[242] | 22 | from mesh import Mesh |
---|
[229] | 23 | from Numeric import array, zeros, Float |
---|
| 24 | |
---|
| 25 | msg = 'First argument in Quantity.__init__ ' |
---|
| 26 | msg += 'must be of class Mesh (or a subclass thereof)' |
---|
| 27 | assert isinstance(domain, Mesh), msg |
---|
| 28 | |
---|
| 29 | if vertex_values is None: |
---|
| 30 | N = domain.number_of_elements |
---|
| 31 | self.vertex_values = zeros((N, 3), Float) |
---|
| 32 | else: |
---|
[265] | 33 | self.vertex_values = array(vertex_values).astype(Float) |
---|
[229] | 34 | |
---|
| 35 | N, V = self.vertex_values.shape |
---|
| 36 | assert V == 3,\ |
---|
| 37 | 'Three vertex values per element must be specified' |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | msg = 'Number of vertex values (%d) must be consistent with'\ |
---|
| 41 | %N |
---|
| 42 | msg += 'number of elements in specified domain (%d).'\ |
---|
| 43 | %domain.number_of_elements |
---|
| 44 | |
---|
| 45 | assert N == domain.number_of_elements, msg |
---|
| 46 | |
---|
| 47 | self.domain = domain |
---|
| 48 | |
---|
| 49 | #Allocate space for other quantities |
---|
| 50 | self.centroid_values = zeros(N, Float) |
---|
| 51 | self.edge_values = zeros((N, 3), Float) |
---|
| 52 | |
---|
| 53 | #Intialise centroid and edge_values |
---|
| 54 | self.interpolate() |
---|
| 55 | |
---|
[275] | 56 | def __len__(self): |
---|
| 57 | return self.centroid_values.shape[0] |
---|
[389] | 58 | |
---|
[229] | 59 | def interpolate(self): |
---|
| 60 | """Compute interpolated values at edges and centroid |
---|
| 61 | Pre-condition: vertex_values have been set |
---|
| 62 | """ |
---|
| 63 | |
---|
| 64 | N = self.vertex_values.shape[0] |
---|
| 65 | for i in range(N): |
---|
| 66 | v0 = self.vertex_values[i, 0] |
---|
| 67 | v1 = self.vertex_values[i, 1] |
---|
| 68 | v2 = self.vertex_values[i, 2] |
---|
| 69 | |
---|
| 70 | self.centroid_values[i] = (v0 + v1 + v2)/3 |
---|
[242] | 71 | |
---|
| 72 | self.interpolate_from_vertices_to_edges() |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | def interpolate_from_vertices_to_edges(self): |
---|
[265] | 76 | #Call correct module function |
---|
| 77 | #(either from this module or C-extension) |
---|
| 78 | interpolate_from_vertices_to_edges(self) |
---|
[229] | 79 | |
---|
[242] | 80 | |
---|
[461] | 81 | def set_values(self, X, location='vertices', indexes = None): |
---|
[242] | 82 | """Set values for quantity |
---|
[229] | 83 | |
---|
[242] | 84 | X: Compatible list, Numeric array (see below), constant or function |
---|
| 85 | location: Where values are to be stored. |
---|
| 86 | Permissible options are: vertices, edges, centroid |
---|
| 87 | Default is "vertices" |
---|
| 88 | |
---|
| 89 | In case of location == 'centroid' the dimension values must |
---|
| 90 | be a list of a Numerical array of length N, N being the number |
---|
[590] | 91 | of elements. Otherwise it must be of dimension Nx3 |
---|
[242] | 92 | |
---|
| 93 | The values will be stored in elements following their |
---|
| 94 | internal ordering. |
---|
| 95 | |
---|
| 96 | If values are described a function, it will be evaluated at specified points |
---|
[590] | 97 | Indexes is the set of element ids that the operation applies to. |
---|
[546] | 98 | |
---|
[242] | 99 | If selected location is vertices, values for centroid and edges |
---|
| 100 | will be assigned interpolated values. |
---|
| 101 | In any other case, only values for the specified locations |
---|
| 102 | will be assigned and the others will be left undefined. |
---|
| 103 | """ |
---|
| 104 | |
---|
| 105 | if location not in ['vertices', 'centroids', 'edges']: |
---|
| 106 | msg = 'Invalid location: %s' %location |
---|
| 107 | raise msg |
---|
| 108 | |
---|
| 109 | if X is None: |
---|
| 110 | msg = 'Given values are None' |
---|
| 111 | raise msg |
---|
| 112 | |
---|
| 113 | import types |
---|
| 114 | |
---|
| 115 | if callable(X): |
---|
| 116 | #Use function specific method |
---|
[590] | 117 | self.set_function_values(X, location, indexes = indexes) |
---|
[242] | 118 | elif type(X) in [types.FloatType, types.IntType, types.LongType]: |
---|
| 119 | if location == 'centroids': |
---|
[517] | 120 | if (indexes == None): |
---|
| 121 | self.centroid_values[:] = X |
---|
| 122 | else: |
---|
| 123 | #Brute force |
---|
[546] | 124 | for i in indexes: |
---|
[517] | 125 | self.centroid_values[i,:] = X |
---|
| 126 | |
---|
[242] | 127 | elif location == 'edges': |
---|
[517] | 128 | if (indexes == None): |
---|
| 129 | self.edge_values[:] = X |
---|
| 130 | else: |
---|
| 131 | #Brute force |
---|
[546] | 132 | for i in indexes: |
---|
[517] | 133 | self.edge_values[i,:] = X |
---|
[242] | 134 | else: |
---|
[517] | 135 | if (indexes == None): |
---|
| 136 | self.vertex_values[:] = X |
---|
| 137 | else: |
---|
| 138 | #Brute force |
---|
[546] | 139 | for i_vertex in indexes: |
---|
| 140 | self.vertex_values[i_vertex,:] = X |
---|
[242] | 141 | |
---|
| 142 | else: |
---|
| 143 | #Use array specific method |
---|
[461] | 144 | self.set_array_values(X, location, indexes = indexes) |
---|
[242] | 145 | |
---|
| 146 | if location == 'vertices': |
---|
| 147 | #Intialise centroid and edge_values |
---|
| 148 | self.interpolate() |
---|
| 149 | |
---|
[531] | 150 | |
---|
| 151 | def get_values(self, location='vertices', indexes = None): |
---|
| 152 | """get values for quantity |
---|
[242] | 153 | |
---|
[531] | 154 | return X, Compatible list, Numeric array (see below) |
---|
| 155 | location: Where values are to be stored. |
---|
| 156 | Permissible options are: vertices, edges, centroid |
---|
| 157 | Default is "vertices" |
---|
[242] | 158 | |
---|
[531] | 159 | In case of location == 'centroid' the dimension values must |
---|
| 160 | be a list of a Numerical array of length N, N being the number |
---|
[590] | 161 | of elements. Otherwise it must be of dimension Nx3 |
---|
[531] | 162 | |
---|
[590] | 163 | The returned values with be a list the length of indexes |
---|
| 164 | (N if indexes = None). Each value will be a list of the three |
---|
| 165 | vertex values for this quantity. |
---|
| 166 | |
---|
| 167 | Indexes is the set of element ids that the operation applies to. |
---|
| 168 | |
---|
[531] | 169 | """ |
---|
| 170 | from Numeric import take |
---|
| 171 | |
---|
| 172 | if location not in ['vertices', 'centroids', 'edges']: |
---|
| 173 | msg = 'Invalid location: %s' %location |
---|
| 174 | raise msg |
---|
| 175 | |
---|
| 176 | if (indexes == None): |
---|
| 177 | indexes = range(len(self)) |
---|
| 178 | |
---|
| 179 | if location == 'centroids': |
---|
| 180 | return take(self.centroid_values,indexes) |
---|
| 181 | elif location == 'edges': |
---|
| 182 | return take(self.edge_values,indexes) |
---|
| 183 | else: |
---|
| 184 | return take(self.vertex_values,indexes) |
---|
| 185 | |
---|
| 186 | |
---|
[590] | 187 | def set_function_values(self, f, location='vertices', indexes = None): |
---|
| 188 | """Set values for quantity using specified function |
---|
| 189 | |
---|
| 190 | f: x, y -> z Function where x, y and z are arrays |
---|
| 191 | location: Where values are to be stored. |
---|
[593] | 192 | Permissible options are: vertices, centroid |
---|
[590] | 193 | Default is "vertices" |
---|
| 194 | """ |
---|
| 195 | from Numeric import take |
---|
| 196 | |
---|
| 197 | if (indexes == None): |
---|
| 198 | indexes = range(len(self)) |
---|
| 199 | is_subset = False |
---|
| 200 | else: |
---|
| 201 | is_subset = True |
---|
| 202 | if location == 'centroids': |
---|
| 203 | P = take(self.domain.centroid_coordinates,indexes) |
---|
| 204 | if is_subset: |
---|
| 205 | self.set_values(f(P[:,0], P[:,1]), location, indexes = indexes) |
---|
| 206 | else: |
---|
| 207 | self.set_values(f(P[:,0], P[:,1]), location) |
---|
| 208 | elif location == 'vertices': |
---|
| 209 | P = self.domain.vertex_coordinates |
---|
| 210 | if is_subset: |
---|
| 211 | #Brute force |
---|
| 212 | for e in indexes: |
---|
| 213 | for i in range(3): |
---|
| 214 | self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1]) |
---|
| 215 | else: |
---|
| 216 | for i in range(3): |
---|
| 217 | self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1]) |
---|
| 218 | else: |
---|
| 219 | raise 'Not implemented: %s' %location |
---|
[591] | 220 | |
---|
[242] | 221 | |
---|
[459] | 222 | def set_array_values(self, values, location='vertices', indexes = None): |
---|
[242] | 223 | """Set values for quantity |
---|
| 224 | |
---|
| 225 | values: Numeric array |
---|
| 226 | location: Where values are to be stored. |
---|
| 227 | Permissible options are: vertices, edges, centroid |
---|
| 228 | Default is "vertices" |
---|
[461] | 229 | |
---|
| 230 | indexes - if this action is carried out on a subset of elements |
---|
| 231 | The element indexes are specified here. |
---|
| 232 | |
---|
[242] | 233 | In case of location == 'centroid' the dimension values must |
---|
| 234 | be a list of a Numerical array of length N, N being the number |
---|
[590] | 235 | of elements. |
---|
[242] | 236 | |
---|
[590] | 237 | Otherwise it must be of dimension Nx3 |
---|
| 238 | |
---|
[242] | 239 | The values will be stored in elements following their |
---|
| 240 | internal ordering. |
---|
| 241 | |
---|
| 242 | If selected location is vertices, values for centroid and edges |
---|
| 243 | will be assigned interpolated values. |
---|
| 244 | In any other case, only values for the specified locations |
---|
| 245 | will be assigned and the others will be left undefined. |
---|
| 246 | """ |
---|
| 247 | |
---|
[459] | 248 | from Numeric import array, Float, Int |
---|
[242] | 249 | |
---|
| 250 | values = array(values).astype(Float) |
---|
| 251 | |
---|
[459] | 252 | if (indexes <> None): |
---|
| 253 | indexes = array(indexes).astype(Int) |
---|
| 254 | msg = 'Number of values must match number of indexes' |
---|
| 255 | assert values.shape[0] == indexes.shape[0], msg |
---|
| 256 | |
---|
[242] | 257 | N = self.centroid_values.shape[0] |
---|
| 258 | |
---|
| 259 | if location == 'centroids': |
---|
| 260 | assert len(values.shape) == 1, 'Values array must be 1d' |
---|
[344] | 261 | |
---|
[459] | 262 | if indexes == None: |
---|
| 263 | msg = 'Number of values must match number of elements' |
---|
| 264 | assert values.shape[0] == N, msg |
---|
[344] | 265 | |
---|
[459] | 266 | self.centroid_values = values |
---|
| 267 | else: |
---|
| 268 | msg = 'Number of values must match number of indexes' |
---|
| 269 | assert values.shape[0] == indexes.shape[0], msg |
---|
| 270 | |
---|
| 271 | #Brute force |
---|
| 272 | for i in range(len(indexes)): |
---|
| 273 | self.centroid_values[indexes[i]] = values[i] |
---|
| 274 | |
---|
[242] | 275 | elif location == 'edges': |
---|
| 276 | assert len(values.shape) == 2, 'Values array must be 2d' |
---|
[344] | 277 | |
---|
| 278 | msg = 'Number of values must match number of elements' |
---|
| 279 | assert values.shape[0] == N, msg |
---|
| 280 | |
---|
[242] | 281 | msg = 'Array must be N x 3' |
---|
| 282 | assert values.shape[1] == 3, msg |
---|
| 283 | |
---|
| 284 | self.edge_values = values |
---|
| 285 | else: |
---|
[344] | 286 | if len(values.shape) == 1: |
---|
[242] | 287 | |
---|
[459] | 288 | if indexes == None: |
---|
| 289 | #Values are being specified once for each unique vertex |
---|
| 290 | msg = 'Number of values must match number of vertices' |
---|
| 291 | assert values.shape[0] == self.domain.coordinates.shape[0], msg |
---|
| 292 | self.set_vertex_values(values) |
---|
| 293 | else: |
---|
| 294 | for element_index, value in map(None, indexes, values): |
---|
[517] | 295 | self.vertex_values[element_index, :] = value |
---|
| 296 | |
---|
[344] | 297 | elif len(values.shape) == 2: |
---|
| 298 | #Vertex values are given as a triplet for each triangle |
---|
| 299 | |
---|
| 300 | msg = 'Array must be N x 3' |
---|
| 301 | assert values.shape[1] == 3, msg |
---|
[459] | 302 | |
---|
| 303 | if indexes == None: |
---|
| 304 | self.vertex_values = values |
---|
| 305 | else: |
---|
| 306 | for element_index, value in map(None, indexes, values): |
---|
| 307 | self.vertex_values[element_index] = value |
---|
[344] | 308 | else: |
---|
| 309 | msg = 'Values array must be 1d or 2d' |
---|
| 310 | raise msg |
---|
[459] | 311 | |
---|
| 312 | |
---|
[389] | 313 | # FIXME have a get_vertex_values as well, so the 'level' quantity can be |
---|
| 314 | # set, based on the elevation |
---|
[324] | 315 | def set_vertex_values(self, A): |
---|
| 316 | """Set vertex values for all triangles based on input array A |
---|
| 317 | which is assumed to have one entry per (unique) vertex, i.e. |
---|
| 318 | one value for each row in array self.domain.coordinates. |
---|
| 319 | """ |
---|
[242] | 320 | |
---|
[324] | 321 | from Numeric import array, Float |
---|
| 322 | |
---|
| 323 | #Assert that A can be converted to a Numeric array of appropriate dim |
---|
| 324 | A = array(A, Float) |
---|
| 325 | |
---|
| 326 | assert len(A.shape) == 1 |
---|
| 327 | assert A.shape[0] == self.domain.coordinates.shape[0] |
---|
| 328 | |
---|
| 329 | N = A.shape[0] |
---|
| 330 | |
---|
| 331 | #Go through list of unique vertices |
---|
| 332 | for k in range(N): |
---|
| 333 | L = self.domain.vertexlist[k] |
---|
| 334 | |
---|
| 335 | if L is None: continue #In case there are unused points |
---|
| 336 | |
---|
| 337 | #Go through all triangle, vertex pairs |
---|
| 338 | #touching vertex k and set corresponding vertex value |
---|
| 339 | for triangle_id, vertex_id in L: |
---|
| 340 | self.vertex_values[triangle_id, vertex_id] = A[k] |
---|
| 341 | |
---|
| 342 | #Intialise centroid and edge_values |
---|
| 343 | self.interpolate() |
---|
[460] | 344 | |
---|
[459] | 345 | |
---|
[274] | 346 | def smooth_vertex_values(self, value_array='field_values', |
---|
| 347 | precision = None): |
---|
| 348 | """ Smooths field_values or conserved_quantities data. |
---|
| 349 | TODO: be able to smooth individual fields |
---|
| 350 | NOTE: This function does not have a test. |
---|
[283] | 351 | FIXME: NOT DONE - do we need it? |
---|
[529] | 352 | FIXME: this function isn't called by anything. |
---|
| 353 | Maybe it should be removed..-DSG |
---|
[274] | 354 | """ |
---|
| 355 | |
---|
| 356 | from Numeric import concatenate, zeros, Float, Int, array, reshape |
---|
| 357 | |
---|
| 358 | |
---|
| 359 | A,V = self.get_vertex_values(xy=False, |
---|
| 360 | value_array=value_array, |
---|
| 361 | smooth = True, |
---|
| 362 | precision = precision) |
---|
| 363 | |
---|
| 364 | #Set some field values |
---|
| 365 | for volume in self: |
---|
| 366 | for i,v in enumerate(volume.vertices): |
---|
| 367 | if value_array == 'field_values': |
---|
| 368 | volume.set_field_values('vertex', i, A[v,:]) |
---|
| 369 | elif value_array == 'conserved_quantities': |
---|
| 370 | volume.set_conserved_quantities('vertex', i, A[v,:]) |
---|
| 371 | |
---|
| 372 | if value_array == 'field_values': |
---|
| 373 | self.precompute() |
---|
| 374 | elif value_array == 'conserved_quantities': |
---|
| 375 | Volume.interpolate_conserved_quantities() |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | #Method for outputting model results |
---|
[288] | 379 | #FIXME: Split up into geometric and numeric stuff. |
---|
| 380 | #FIXME: Geometric (X,Y,V) should live in mesh.py |
---|
[292] | 381 | #FIXME: STill remember to move XY to mesh |
---|
[274] | 382 | def get_vertex_values(self, |
---|
| 383 | xy=True, |
---|
| 384 | smooth = None, |
---|
| 385 | precision = None, |
---|
| 386 | reduction = None): |
---|
| 387 | """Return vertex values like an OBJ format |
---|
| 388 | |
---|
| 389 | The vertex values are returned as one sequence in the 1D float array A. |
---|
| 390 | If requested the coordinates will be returned in 1D arrays X and Y. |
---|
| 391 | |
---|
| 392 | The connectivity is represented as an integer array, V, of dimension |
---|
| 393 | M x 3, where M is the number of volumes. Each row has three indices |
---|
| 394 | into the X, Y, A arrays defining the triangle. |
---|
| 395 | |
---|
| 396 | if smooth is True, vertex values corresponding to one common |
---|
| 397 | coordinate set will be smoothed according to the given |
---|
| 398 | reduction operator. In this case vertex coordinates will be |
---|
| 399 | de-duplicated. |
---|
| 400 | |
---|
| 401 | If no smoothings is required, vertex coordinates and values will |
---|
[292] | 402 | be aggregated as a concatenation of values at |
---|
[274] | 403 | vertices 0, vertices 1 and vertices 2 |
---|
| 404 | |
---|
| 405 | |
---|
| 406 | Calling convention |
---|
| 407 | if xy is True: |
---|
| 408 | X,Y,A,V = get_vertex_values |
---|
| 409 | else: |
---|
| 410 | A,V = get_vertex_values |
---|
| 411 | |
---|
| 412 | """ |
---|
| 413 | |
---|
| 414 | from Numeric import concatenate, zeros, Float, Int, array, reshape |
---|
| 415 | |
---|
| 416 | |
---|
| 417 | if smooth is None: |
---|
| 418 | smooth = self.domain.smooth |
---|
| 419 | |
---|
| 420 | if precision is None: |
---|
| 421 | precision = Float |
---|
| 422 | |
---|
| 423 | if reduction is None: |
---|
| 424 | reduction = self.domain.reduction |
---|
[291] | 425 | |
---|
| 426 | #Create connectivity |
---|
| 427 | |
---|
[274] | 428 | if smooth == True: |
---|
[275] | 429 | |
---|
[529] | 430 | V = self.domain.get_vertices() |
---|
[275] | 431 | N = len(self.domain.vertexlist) |
---|
| 432 | A = zeros(N, precision) |
---|
[291] | 433 | |
---|
[274] | 434 | #Smoothing loop |
---|
[275] | 435 | for k in range(N): |
---|
| 436 | L = self.domain.vertexlist[k] |
---|
| 437 | |
---|
| 438 | #Go through all triangle, vertex pairs |
---|
| 439 | #contributing to vertex k and register vertex value |
---|
[297] | 440 | |
---|
| 441 | if L is None: continue #In case there are unused points |
---|
| 442 | |
---|
| 443 | contributions = [] |
---|
[275] | 444 | for volume_id, vertex_id in L: |
---|
| 445 | v = self.vertex_values[volume_id, vertex_id] |
---|
| 446 | contributions.append(v) |
---|
[274] | 447 | |
---|
[275] | 448 | A[k] = reduction(contributions) |
---|
[274] | 449 | |
---|
| 450 | |
---|
| 451 | if xy is True: |
---|
[281] | 452 | X = self.domain.coordinates[:,0].astype(precision) |
---|
| 453 | Y = self.domain.coordinates[:,1].astype(precision) |
---|
| 454 | |
---|
[274] | 455 | return X, Y, A, V |
---|
| 456 | else: |
---|
| 457 | return A, V |
---|
| 458 | else: |
---|
| 459 | #Don't smooth |
---|
| 460 | |
---|
[529] | 461 | # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]] |
---|
| 462 | # These vert_id's will relate to the verts created bellow |
---|
| 463 | m = len(self.domain) #Number of volumes |
---|
| 464 | M = 3*m #Total number of unique vertices |
---|
| 465 | V = reshape(array(range(M)).astype(Int), (m,3)) |
---|
| 466 | |
---|
[281] | 467 | A = self.vertex_values.flat |
---|
[274] | 468 | |
---|
| 469 | #Do vertex coordinates |
---|
| 470 | if xy is True: |
---|
[275] | 471 | C = self.domain.get_vertex_coordinates() |
---|
| 472 | |
---|
[282] | 473 | X = C[:,0:6:2].copy() |
---|
| 474 | Y = C[:,1:6:2].copy() |
---|
[275] | 475 | |
---|
| 476 | return X.flat, Y.flat, A, V |
---|
[274] | 477 | else: |
---|
| 478 | return A, V |
---|
| 479 | |
---|
| 480 | |
---|
| 481 | |
---|
| 482 | |
---|
| 483 | |
---|
[242] | 484 | class Conserved_quantity(Quantity): |
---|
| 485 | """Class conserved quantity adds to Quantity: |
---|
| 486 | |
---|
| 487 | boundary values, storage and method for updating, and |
---|
| 488 | methods for extrapolation from centropid to vertices inluding |
---|
| 489 | gradients and limiters |
---|
| 490 | """ |
---|
| 491 | |
---|
| 492 | def __init__(self, domain, vertex_values=None): |
---|
| 493 | Quantity.__init__(self, domain, vertex_values) |
---|
| 494 | |
---|
| 495 | from Numeric import zeros, Float |
---|
| 496 | |
---|
| 497 | #Allocate space for boundary values |
---|
| 498 | L = len(domain.boundary) |
---|
| 499 | self.boundary_values = zeros(L, Float) |
---|
| 500 | |
---|
| 501 | #Allocate space for updates of conserved quantities by |
---|
| 502 | #flux calculations and forcing functions |
---|
| 503 | |
---|
| 504 | N = domain.number_of_elements |
---|
| 505 | self.explicit_update = zeros(N, Float ) |
---|
| 506 | self.semi_implicit_update = zeros(N, Float ) |
---|
| 507 | |
---|
| 508 | |
---|
[229] | 509 | def update(self, timestep): |
---|
[272] | 510 | #Call correct module function |
---|
| 511 | #(either from this module or C-extension) |
---|
| 512 | return update(self, timestep) |
---|
[229] | 513 | |
---|
| 514 | |
---|
| 515 | def compute_gradients(self): |
---|
[260] | 516 | #Call correct module function |
---|
| 517 | #(either from this module or C-extension) |
---|
| 518 | return compute_gradients(self) |
---|
| 519 | |
---|
[229] | 520 | |
---|
| 521 | def limit(self): |
---|
[245] | 522 | #Call correct module function |
---|
| 523 | #(either from this module or C-extension) |
---|
| 524 | limit(self) |
---|
[229] | 525 | |
---|
| 526 | |
---|
| 527 | def extrapolate_first_order(self): |
---|
| 528 | """Extrapolate conserved quantities from centroid to |
---|
| 529 | vertices for each volume using |
---|
| 530 | first order scheme. |
---|
| 531 | """ |
---|
| 532 | |
---|
| 533 | qc = self.centroid_values |
---|
| 534 | qv = self.vertex_values |
---|
| 535 | |
---|
| 536 | for i in range(3): |
---|
| 537 | qv[:,i] = qc |
---|
| 538 | |
---|
| 539 | |
---|
| 540 | def extrapolate_second_order(self): |
---|
[255] | 541 | #Call correct module function |
---|
| 542 | #(either from this module or C-extension) |
---|
| 543 | extrapolate_second_order(self) |
---|
| 544 | |
---|
| 545 | |
---|
[272] | 546 | def update(quantity, timestep): |
---|
| 547 | """Update centroid values based on values stored in |
---|
| 548 | explicit_update and semi_implicit_update as well as given timestep |
---|
[476] | 549 | |
---|
| 550 | Function implementing forcing terms must take on argument |
---|
| 551 | which is the domain and they must update either explicit |
---|
| 552 | or implicit updates, e,g,: |
---|
| 553 | |
---|
| 554 | def gravity(domain): |
---|
| 555 | .... |
---|
| 556 | domain.quantities['xmomentum'].explicit_update = ... |
---|
| 557 | domain.quantities['ymomentum'].explicit_update = ... |
---|
| 558 | |
---|
| 559 | |
---|
| 560 | |
---|
| 561 | Explicit terms must have the form |
---|
| 562 | |
---|
| 563 | G(q, t) |
---|
| 564 | |
---|
| 565 | and explicit scheme is |
---|
| 566 | |
---|
| 567 | q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) |
---|
| 568 | |
---|
| 569 | |
---|
| 570 | Semi implicit forcing terms are assumed to have the form |
---|
| 571 | |
---|
| 572 | G(q, t) = H(q, t) q |
---|
| 573 | |
---|
| 574 | and the semi implicit scheme will then be |
---|
| 575 | |
---|
| 576 | q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) |
---|
| 577 | |
---|
| 578 | |
---|
[272] | 579 | """ |
---|
| 580 | |
---|
| 581 | from Numeric import sum, equal, ones, Float |
---|
| 582 | |
---|
| 583 | N = quantity.centroid_values.shape[0] |
---|
[458] | 584 | |
---|
| 585 | |
---|
[476] | 586 | #Divide H by conserved quantity to obtain G (see docstring above) |
---|
[514] | 587 | |
---|
| 588 | |
---|
[458] | 589 | for k in range(N): |
---|
| 590 | x = quantity.centroid_values[k] |
---|
| 591 | if x == 0.0: |
---|
[514] | 592 | #FIXME: Is this right |
---|
[458] | 593 | quantity.semi_implicit_update[k] = 0.0 |
---|
| 594 | else: |
---|
| 595 | quantity.semi_implicit_update[k] /= x |
---|
| 596 | |
---|
[272] | 597 | #Explicit updates |
---|
| 598 | quantity.centroid_values += timestep*quantity.explicit_update |
---|
| 599 | |
---|
| 600 | #Semi implicit updates |
---|
| 601 | denominator = ones(N, Float)-timestep*quantity.semi_implicit_update |
---|
[265] | 602 | |
---|
[272] | 603 | if sum(equal(denominator, 0.0)) > 0.0: |
---|
| 604 | msg = 'Zero division in semi implicit update. Call Stephen :-)' |
---|
| 605 | raise msg |
---|
| 606 | else: |
---|
| 607 | #Update conserved_quantities from semi implicit updates |
---|
| 608 | quantity.centroid_values /= denominator |
---|
| 609 | |
---|
| 610 | |
---|
[265] | 611 | def interpolate_from_vertices_to_edges(quantity): |
---|
| 612 | """Compute edge values from vertex values using linear interpolation |
---|
| 613 | """ |
---|
| 614 | |
---|
| 615 | for k in range(quantity.vertex_values.shape[0]): |
---|
| 616 | q0 = quantity.vertex_values[k, 0] |
---|
| 617 | q1 = quantity.vertex_values[k, 1] |
---|
| 618 | q2 = quantity.vertex_values[k, 2] |
---|
| 619 | |
---|
| 620 | quantity.edge_values[k, 0] = 0.5*(q1+q2) |
---|
| 621 | quantity.edge_values[k, 1] = 0.5*(q0+q2) |
---|
| 622 | quantity.edge_values[k, 2] = 0.5*(q0+q1) |
---|
| 623 | |
---|
| 624 | |
---|
| 625 | |
---|
| 626 | def extrapolate_second_order(quantity): |
---|
[255] | 627 | """Extrapolate conserved quantities from centroid to |
---|
| 628 | vertices for each volume using |
---|
| 629 | second order scheme. |
---|
| 630 | """ |
---|
[229] | 631 | |
---|
[265] | 632 | a, b = quantity.compute_gradients() |
---|
[229] | 633 | |
---|
[265] | 634 | X = quantity.domain.get_vertex_coordinates() |
---|
| 635 | qc = quantity.centroid_values |
---|
| 636 | qv = quantity.vertex_values |
---|
[255] | 637 | |
---|
| 638 | #Check each triangle |
---|
[265] | 639 | for k in range(quantity.domain.number_of_elements): |
---|
[255] | 640 | #Centroid coordinates |
---|
[305] | 641 | x, y = quantity.domain.centroid_coordinates[k] |
---|
[229] | 642 | |
---|
[255] | 643 | #vertex coordinates |
---|
[260] | 644 | x0, y0, x1, y1, x2, y2 = X[k,:] |
---|
[255] | 645 | |
---|
| 646 | #Extrapolate |
---|
| 647 | qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y) |
---|
| 648 | qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y) |
---|
| 649 | qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y) |
---|
[229] | 650 | |
---|
[245] | 651 | |
---|
[260] | 652 | def compute_gradients(quantity): |
---|
| 653 | """Compute gradients of triangle surfaces defined by centroids of |
---|
| 654 | neighbouring volumes. |
---|
| 655 | If one edge is on the boundary, use own centroid as neighbour centroid. |
---|
| 656 | If two or more are on the boundary, fall back to first order scheme. |
---|
| 657 | """ |
---|
| 658 | |
---|
| 659 | from Numeric import zeros, Float |
---|
| 660 | from util import gradient |
---|
| 661 | |
---|
[305] | 662 | centroid_coordinates = quantity.domain.centroid_coordinates |
---|
[260] | 663 | surrogate_neighbours = quantity.domain.surrogate_neighbours |
---|
| 664 | centroid_values = quantity.centroid_values |
---|
| 665 | number_of_boundaries = quantity.domain.number_of_boundaries |
---|
| 666 | |
---|
| 667 | N = centroid_values.shape[0] |
---|
| 668 | |
---|
| 669 | a = zeros(N, Float) |
---|
| 670 | b = zeros(N, Float) |
---|
| 671 | |
---|
| 672 | for k in range(N): |
---|
| 673 | if number_of_boundaries[k] < 2: |
---|
| 674 | #Two or three true neighbours |
---|
| 675 | |
---|
| 676 | #Get indices of neighbours (or self when used as surrogate) |
---|
| 677 | k0, k1, k2 = surrogate_neighbours[k,:] |
---|
| 678 | |
---|
[261] | 679 | #Get data |
---|
[260] | 680 | q0 = centroid_values[k0] |
---|
| 681 | q1 = centroid_values[k1] |
---|
| 682 | q2 = centroid_values[k2] |
---|
| 683 | |
---|
[305] | 684 | x0, y0 = centroid_coordinates[k0] #V0 centroid |
---|
| 685 | x1, y1 = centroid_coordinates[k1] #V1 centroid |
---|
| 686 | x2, y2 = centroid_coordinates[k2] #V2 centroid |
---|
[260] | 687 | |
---|
| 688 | #Gradient |
---|
| 689 | a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) |
---|
| 690 | |
---|
| 691 | elif number_of_boundaries[k] == 2: |
---|
| 692 | #One true neighbour |
---|
| 693 | |
---|
| 694 | #Get index of the one neighbour |
---|
| 695 | for k0 in surrogate_neighbours[k,:]: |
---|
| 696 | if k0 != k: break |
---|
| 697 | assert k0 != k |
---|
| 698 | |
---|
| 699 | k1 = k #self |
---|
| 700 | |
---|
| 701 | #Get data |
---|
| 702 | q0 = centroid_values[k0] |
---|
| 703 | q1 = centroid_values[k1] |
---|
| 704 | |
---|
[305] | 705 | x0, y0 = centroid_coordinates[k0] #V0 centroid |
---|
| 706 | x1, y1 = centroid_coordinates[k1] #V1 centroid |
---|
[260] | 707 | |
---|
| 708 | #Gradient |
---|
| 709 | det = x0*y1 - x1*y0 |
---|
| 710 | if det != 0.0: |
---|
| 711 | a[k] = (y1*q0 - y0*q1)/det |
---|
| 712 | b[k] = (x0*q1 - x1*q0)/det |
---|
| 713 | |
---|
| 714 | else: |
---|
| 715 | #No true neighbours - |
---|
| 716 | #Fall back to first order scheme |
---|
| 717 | pass |
---|
| 718 | |
---|
| 719 | |
---|
| 720 | return a, b |
---|
| 721 | |
---|
| 722 | |
---|
| 723 | |
---|
[245] | 724 | def limit(quantity): |
---|
| 725 | """Limit slopes for each volume to eliminate artificial variance |
---|
| 726 | introduced by e.g. second order extrapolator |
---|
| 727 | |
---|
| 728 | This is an unsophisticated limiter as it does not take into |
---|
| 729 | account dependencies among quantities. |
---|
| 730 | |
---|
| 731 | precondition: |
---|
| 732 | vertex values are estimated from gradient |
---|
| 733 | postcondition: |
---|
| 734 | vertex values are updated |
---|
| 735 | """ |
---|
| 736 | |
---|
| 737 | from Numeric import zeros, Float |
---|
| 738 | |
---|
| 739 | N = quantity.domain.number_of_elements |
---|
| 740 | |
---|
| 741 | beta = quantity.domain.beta |
---|
| 742 | |
---|
| 743 | qc = quantity.centroid_values |
---|
| 744 | qv = quantity.vertex_values |
---|
| 745 | |
---|
| 746 | #Find min and max of this and neighbour's centroid values |
---|
| 747 | qmax = zeros(qc.shape, Float) |
---|
| 748 | qmin = zeros(qc.shape, Float) |
---|
| 749 | |
---|
| 750 | for k in range(N): |
---|
| 751 | qmax[k] = qmin[k] = qc[k] |
---|
| 752 | for i in range(3): |
---|
| 753 | n = quantity.domain.neighbours[k,i] |
---|
| 754 | if n >= 0: |
---|
| 755 | qn = qc[n] #Neighbour's centroid value |
---|
| 756 | |
---|
| 757 | qmin[k] = min(qmin[k], qn) |
---|
| 758 | qmax[k] = max(qmax[k], qn) |
---|
| 759 | |
---|
| 760 | |
---|
| 761 | #Diffences between centroids and maxima/minima |
---|
| 762 | dqmax = qmax - qc |
---|
| 763 | dqmin = qmin - qc |
---|
| 764 | |
---|
| 765 | #Deltas between vertex and centroid values |
---|
| 766 | dq = zeros(qv.shape, Float) |
---|
| 767 | for i in range(3): |
---|
| 768 | dq[:,i] = qv[:,i] - qc |
---|
| 769 | |
---|
| 770 | #Phi limiter |
---|
| 771 | for k in range(N): |
---|
| 772 | |
---|
| 773 | #Find the gradient limiter (phi) across vertices |
---|
| 774 | phi = 1.0 |
---|
| 775 | for i in range(3): |
---|
| 776 | r = 1.0 |
---|
| 777 | if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] |
---|
| 778 | if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] |
---|
| 779 | |
---|
| 780 | phi = min( min(r*beta, 1), phi ) |
---|
| 781 | |
---|
| 782 | #Then update using phi limiter |
---|
| 783 | for i in range(3): |
---|
| 784 | qv[k,i] = qc[k] + phi*dq[k,i] |
---|
| 785 | |
---|
| 786 | |
---|
| 787 | |
---|
| 788 | import compile |
---|
| 789 | if compile.can_use_C_extension('quantity_ext.c'): |
---|
| 790 | #Replace python version with c implementations |
---|
[259] | 791 | |
---|
[262] | 792 | from quantity_ext import limit, compute_gradients,\ |
---|
[272] | 793 | extrapolate_second_order, interpolate_from_vertices_to_edges, update |
---|
[265] | 794 | |
---|