Changeset 245
- Timestamp:
- Aug 30, 2004, 6:44:45 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/quantity.py
r242 r245 333 333 334 334 def limit(self): 335 """Limit slopes for each volume to eliminate artificial variance 336 introduced by e.g. second order extrapolator 337 338 This is an unsophisticated limiter as it does not take into 339 account dependencies among quantities. 340 341 precondition: 342 vertex values are estimated from gradient 343 postcondition: 344 vertex values are updated 345 """ 346 347 from Numeric import zeros, Float 348 349 beta = self.domain.beta 350 351 qc = self.centroid_values 352 qv = self.vertex_values 353 354 #Deltas between vertex and centroid values 355 dq = zeros( qv.shape, Float) 356 for i in range(3): 357 dq[:,i] = qv[:,i] - qc 358 359 #Find min and max of this and neighbour's centroid values 360 qmax = zeros(qc.shape, Float) 361 qmin = zeros(qc.shape, Float) 362 363 for k in range(self.domain.number_of_elements): 364 qmax[k] = qmin[k] = qc[k] 365 for i in range(3): 366 n = self.domain.neighbours[k,i] 367 if n >= 0: 368 qn = qc[n] #Neighbour's centroid value 369 370 qmin[k] = min(qmin[k], qn) 371 qmax[k] = max(qmax[k], qn) 372 373 374 #Diffences between centroids and maxima/minima 375 dqmax = qmax - qc 376 dqmin = qmin - qc 377 378 for k in range(self.domain.number_of_elements): 379 380 #Find the gradient limiter (phi) across vertices 381 phi = 1.0 382 for i in range(3): 383 r = 1.0 384 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] 385 if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] 386 387 phi = min( min(r*beta, 1), phi ) 388 389 #Then update using phi limiter 390 for i in range(3): 391 qv[k,i] = qc[k] + phi*dq[k,i] 392 335 #Call correct module function 336 #(either from this module or C-extension) 337 limit(self) 338 393 339 394 340 def extrapolate_first_order(self): … … 431 377 432 378 379 380 def limit(quantity): 381 """Limit slopes for each volume to eliminate artificial variance 382 introduced by e.g. second order extrapolator 383 384 This is an unsophisticated limiter as it does not take into 385 account dependencies among quantities. 386 387 precondition: 388 vertex values are estimated from gradient 389 postcondition: 390 vertex values are updated 391 """ 392 393 from Numeric import zeros, Float 394 395 N = quantity.domain.number_of_elements 396 397 beta = quantity.domain.beta 398 399 qc = quantity.centroid_values 400 qv = quantity.vertex_values 401 402 #Find min and max of this and neighbour's centroid values 403 qmax = zeros(qc.shape, Float) 404 qmin = zeros(qc.shape, Float) 405 406 for k in range(N): 407 qmax[k] = qmin[k] = qc[k] 408 for i in range(3): 409 n = quantity.domain.neighbours[k,i] 410 if n >= 0: 411 qn = qc[n] #Neighbour's centroid value 412 413 qmin[k] = min(qmin[k], qn) 414 qmax[k] = max(qmax[k], qn) 415 416 417 #Diffences between centroids and maxima/minima 418 dqmax = qmax - qc 419 dqmin = qmin - qc 420 421 #Deltas between vertex and centroid values 422 dq = zeros(qv.shape, Float) 423 for i in range(3): 424 dq[:,i] = qv[:,i] - qc 425 426 #Phi limiter 427 for k in range(N): 428 429 #Find the gradient limiter (phi) across vertices 430 phi = 1.0 431 for i in range(3): 432 r = 1.0 433 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] 434 if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] 435 436 phi = min( min(r*beta, 1), phi ) 437 438 #Then update using phi limiter 439 for i in range(3): 440 qv[k,i] = qc[k] + phi*dq[k,i] 441 442 443 444 import compile 445 if compile.can_use_C_extension('quantity_ext.c'): 446 #Replace python version with c implementations 447 448 pass 449 #from quantity_ext import limit 450 -
inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
r240 r245 17 17 #include "math.h" 18 18 19 //FIXME: Should live in util_ext.c 19 20 double max(double x, double y) { 20 21 //Return maximum of two doubles … … 33 34 34 35 // Computational function for rotation 35 // FIXME: Should probably go into shallow_water_ext.c36 36 int _rotate(double *q, double n1, double n2) { 37 37 /*Rotate the momentum component q (q[1], q[2])
Note: See TracChangeset
for help on using the changeset viewer.