Changeset 205
- Timestamp:
- Aug 23, 2004, 2:45:42 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 8 edited
-
config.py (modified) (1 diff)
-
domain.py (modified) (6 diffs)
-
quantity.py (modified) (7 diffs)
-
shallow_water.py (modified) (9 diffs)
-
test_quantity.py (modified) (6 diffs)
-
test_shallow_water.py (modified) (3 diffs)
-
test_util.py (modified) (4 diffs)
-
util_ext.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/config.py
r198 r205 57 57 58 58 use_extensions = True #Try to use C-extensions 59 use_extensions = False #Do not use C-extensions59 #use_extensions = False #Do not use C-extensions 60 60 61 61 use_psyco = True #Use psyco optimisations -
inundation/ga/storm_surge/pyvolution/domain.py
r196 r205 36 36 for name in self.conserved_quantities + other_quantities: 37 37 self.quantities[name] = Quantity(self) 38 39 40 #Create an empty list for explicit forcing terms 41 # 42 # Explicit terms must have the form 43 # 44 # G(q, t) 45 # 46 # and explicit scheme is 47 # 48 # q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) 49 # 50 # 51 # FIXME: How to call and how function should look 52 53 self.explicit_forcing_terms = [] 54 55 56 #Create an empty list for semi implicit forcing terms if any 57 # 58 # Semi implicit forcing terms are assumed to have the form 59 # 60 # G(q, t) = H(q, t) q 61 # 62 # and the semi implicit scheme will then be 63 # 64 # q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) 65 66 self.semi_implicit_forcing_terms = [] 67 38 68 39 69 #Defaults … … 62 92 self.store=False 63 93 self.format = 'dat' 64 #self.smooth = False65 94 self.smooth = True 66 95 … … 386 415 they should be defined in Domain subclass 387 416 """ 388 pass 417 418 for f in self.explicit_forcing_terms: 419 420 pass 421 389 422 390 423 … … 395 428 396 429 from Numeric import ones, sum, equal, Float 397 398 430 399 431 N = self.number_of_elements … … 402 434 timestep = self.timestep 403 435 436 #Reset semi_implicit_updates 437 #(explicit updates are already set by compute_fluyzes) 438 for i, name in enumerate(self.conserved_quantities): 439 Q = self.quantities[name] 440 Q.semi_implicit_update[:] = 0. 441 442 #Compute forcing terms 404 443 self.compute_forcing_terms() 405 444 … … 421 460 for name in self.conserved_quantities: 422 461 Q = self.quantities[name] 423 Q.distribute_to_vertices_and_edges() 424 462 if self.order == 1: 463 Q.first_order_extrapolator() 464 elif self.order == 2: 465 Q.second_order_extrapolator() 466 Q.limiter() 467 else: 468 raise 'Unknown order' 469 Q.interpolate_from_vertices_to_edges() 425 470 426 471 -
inundation/ga/storm_surge/pyvolution/quantity.py
r198 r205 106 106 self.centroid_values /= denominator 107 107 108 108 109 def compute_gradients(self): 109 110 """Compute gradients of triangle surfaces defined by centroids of … … 138 139 139 140 k0 = k #Self 140 141 141 #Find index of the one neighbour 142 142 for k1 in self.domain.neighbours[k,:]: … … 144 144 break 145 145 assert k1 != k0 146 assert k1 >= 0 147 146 assert k1 >= 0 147 148 148 #Get data 149 149 q0 = self.centroid_values[k0] … … 159 159 b[k] = (x0*q1 - x1*q0)/det 160 160 161 161 162 else: 162 163 #One or zero missing neighbours … … 180 181 #Gradient 181 182 a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) 182 #array([q0]), array([q1]), array([q2]))183 183 184 184 return a, b 185 185 186 186 187 def second_order_limiter(self): 188 """Extrapolate conserved quantities from centroid to 189 vertices for each volume using second order scheme. 190 187 def limiter(self): 188 """Limit slopes for each volume to eliminate artificial variance 189 introduced by e.g. second order extrapolator 190 191 This is an unsophisticated limiter as it does not take into 192 account dependencies among quantities. 193 191 194 precondition: 192 195 vertex values are estimated from gradient … … 242 245 243 246 244 def first_order_ limiter(self):247 def first_order_extrapolator(self): 245 248 """Extrapolate conserved quantities from centroid to 246 249 vertices for each volume using … … 254 257 qv[:,i] = qc 255 258 256 257 def distribute_to_vertices_and_edges(self, order=1):259 260 def second_order_extrapolator(self): 258 261 """Extrapolate conserved quantities from centroid to 259 vertices and edge-midpoints for each volume using 260 specified order. 261 262 This is an unsophisticated limiter as it does not take into 263 account dependencies among quantities. 264 """ 265 266 if order == 1: 267 self.first_order_limiter() 268 elif order == 2: 269 a, b = self.compute_gradients() 270 271 self.second_order_limiter() 272 else: 273 msg = 'ERROR: unknown order %d' %order 274 raise msg 275 276 self.interpolate_from_vertices_to_edges() 277 262 vertices for each volume using 263 second order scheme. 264 """ 265 266 a, b = self.compute_gradients() 267 268 V = self.domain.get_vertex_coordinates() 269 qc = self.centroid_values 270 qv = self.vertex_values 271 272 #Check each triangle 273 for k in range(self.domain.number_of_elements): 274 #Centroid coordinates 275 x, y = self.domain.centroids[k] 276 277 #vertex coordinates 278 x0, y0, x1, y1, x2, y2 = V[k,:] 279 280 #Extrapolate 281 qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y) 282 qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y) 283 qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y) 284 285 278 286 279 287 def interpolate_from_vertices_to_edges(self): -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r198 r205 71 71 compute_fluxes(self) 72 72 73 def first_order_limiter(self):73 def distribute_to_vertices_and_edges(self): 74 74 #Call correct module function 75 first_order_limiter(self) 75 #(either from this module or C-extension) 76 distribute_to_vertices_and_edges(self) 76 77 77 78 #################################### … … 259 260 Xmom.explicit_update[k] = flux[1] 260 261 Ymom.explicit_update[k] = flux[2] 261 262 262 263 timestep = min(timestep, optimal_timestep) 263 264 264 265 domain.timestep = timestep 266 265 267 266 268 267 269 #################################### 268 # Functions for gradient limiting (distribute_to_vertices_and_edges) 270 # Module functions for gradient limiting (distribute_to_vertices_and_edges) 271 272 def distribute_to_vertices_and_edges(domain): 273 274 protect_against_infinitesimal_heights_centroid(domain) 275 if domain.order == 1: 276 first_order_extrapolator(domain) 277 elif domain.order == 2: 278 second_order_extrapolator(domain) 279 else: 280 raise 'Unknown order' 269 281 270 282 def protect_against_infinitesimal_heights_centroid(domain): … … 297 309 298 310 299 300 301 def first_order_limiter(domain): 302 """First order limiter function, specific to the shallow water wave 303 equation. 311 312 313 314 def first_order_extrapolator(domain): 315 """First order extrapolator function, specific 316 to the shallow water wave equation. 304 317 305 318 It will ensure that h (w-z) is always non-negative even in the … … 315 328 """ 316 329 317 #FIXME: Might go elsewhere318 protect_against_infinitesimal_heights_centroid(domain)319 330 320 331 #Update conserved quantities using straight first order 321 332 for name in domain.conserved_quantities: 322 333 Q = domain.quantities[name] 323 Q.first_order_limiter() 324 325 334 Q.first_order_extrapolator() 335 Q.interpolate_from_vertices_to_edges() 336 337 326 338 #Water levels at centroids 327 339 wc = domain.quantities['level'].centroid_values … … 377 389 alpha = 1.0 378 390 379 #Update stage391 #Update water levels 380 392 for i in range(3): 381 393 #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i]) … … 384 396 385 397 386 def second_order_ limiter(domain):398 def second_order_extrapolator(domain): 387 399 """Second order limiter function, specific to the shallow water wave 388 400 equation. … … 406 418 407 419 """ 408 409 #FIXME: Might go elsewhere 410 protect_against_infinitesimal_heights_centroid(domain) 420 421 #FIXME: first and second order migh merge 422 423 from Numeric import minimum, maximum 411 424 412 425 #Update conserved quantities using straight second order 413 426 for name in domain.conserved_quantities: 414 427 Q = domain.quantities[name] 415 Q.second_order_limiter() 416 417 428 429 Q.second_order_extrapolator() 430 Q.limiter() 431 Q.interpolate_from_vertices_to_edges() 432 433 418 434 #Water levels at centroids 419 435 wc = domain.quantities['level'].centroid_values … … 425 441 hc = wc - zc 426 442 427 #Momentums at centroids428 xmomc = domain.quantities['xmomentum'].centroid_values429 ymomc = domain.quantities['ymomentum'].centroid_values430 431 443 #Water levels at vertices 432 #wv = domain.quantities['level'].vertex_values433 # 444 wv = domain.quantities['level'].vertex_values 445 434 446 #Bed elevations at vertices 435 #zv = domain.quantities['elevation'].vertex_values 436 # 437 #Momentums at vertices 438 #xmomv = domain.quantities['xmomentum'].vertex_values 439 #ymomv = domain.quantities['ymomentum'].vertex_values 440 447 zv = domain.quantities['elevation'].vertex_values 448 449 #Water depths at vertices 450 hv = wv-zv 451 452 453 454 #Computed linear combination between constant levels and and 455 #levels parallel to the bed elevation. 456 for k in range(domain.number_of_elements): 457 458 #Compute maximal variation in bed elevation 459 # This quantitiy is 460 # dz = max_i abs(z_i - z_c) 461 # and it is independent of dimension 462 # In the 1d case zc = (z0+z1)/2 463 # In the 2d case zc = (z0+z1+z2)/3 464 # 465 466 z_range = max(abs(zv[k,0]-zc[k]), 467 abs(zv[k,1]-zc[k]), 468 abs(zv[k,2]-zc[k])) 469 470 471 hmin = minimum( hv[k, :] ) 472 473 #Create alpha in [0,1], where alpha==0 means using shallow 474 #first order scheme and alpha==1 means using the limited 475 #second order scheme for the stage w 476 #If hmin < 0 then alpha=0 reverrting to first order. 477 478 if z_range > 0.0: 479 alpha = max( min( 2*hmin/z_range, 1.0), 0.0) 480 else: 481 alpha = 1.0 482 483 #Update water levels 484 # (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 485 # zvi + hc + alpha*(hvi - hc) 486 for i in range(3): 487 wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) 488 489 490 #Momentums at centroids 491 xmomc = domain.quantities['xmomentum'].centroid_values 492 ymomc = domain.quantities['ymomentum'].centroid_values 493 494 #Momentums at vertices 495 xmomv = domain.quantities['xmomentum'].vertex_values 496 ymomv = domain.quantities['ymomentum'].vertex_values 497 498 # Update momentum as a linear combination of 499 # xmomc and ymomc (shallow) and momentum 500 # from extrapolator xmomv and ymomv (deep). 501 502 503 ##f##or i in range(3): 504 ##### xmomv[k,i] = (1-alpha)*xmomc[k] + alpha*xmomv[k,i]; 505 506 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]; 507 508 509 #Finally, protect against infinitesimal heights and high speeds 510 #Water depths at vertices 511 hv = wv-zv 512 hc = wc-zc 441 513 for k in range(domain.number_of_elements): 442 #Protect against infinitesimal heights and high velocities 443 if hc[k] < domain.minimum_allowed_height: 444 #Control level and height 514 hmax = maximum(hv[k,:]) 515 516 if hmax < minimum_allowed_height: 517 #Reset negative heights to bed elevation 445 518 if hc[k] < 0.0: 446 wc[k] = zc[k]; hc[k] = 0.0 447 448 #Control momentum as well! 449 xmomc[k] = ymomc[k] = 0.0 450 451 452 #Compute second order limiter for each conserved quantity 453 for name in domain.conserved_quantities: 454 domain.quantities[name].second_order_limiter() 455 456 457 # //Adjusted water heights at vertices 458 # hv0 = qv0[0]-zv0; 459 # hv1 = qv1[0]-zv1; 460 # hv2 = qv2[0]-zv2; 461 462 # hmin = min( min(hv0, hv1), hv2 ); 463 464 # // ----------------------------------- 465 # // Second run - Linear combination with 466 # // first order scheme for shallow depths 467 # // ----------------------------------- 468 469 470 # //Compute maximal variation in bed elevation 471 # // 472 # //This quantitiy is 473 # // dz = max_i abs(z_i - z_c) 474 # //and it is independent of dimension 475 # //In the 1d case zc = (z0+z1)/2 476 # //In the 2d case zc = (z0+z1+z2)/3 477 # // 478 # z_range = max( max(fabs(zv0-zc), fabs(zv1-zc)), fabs(zv2-zc) ); 479 480 # //Create alpha in [0,1], where alpha==0 means using shallow 481 # //first order scheme and alpha==1 means using the limited 482 # //second order scheme for the stage w 483 # //If hmin < 0 then alpha=0 and first order method is 484 # if (z_range>0.0) { 485 # alpha = max( min(2*hmin/z_range, 1.0), 0.0); 486 # } else { 487 # alpha = 1.0; 488 # } 489 490 # if (alpha < 1) { 491 # //Update stage 492 # // 493 # //(1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 494 # //zvi + hc + alpha*(hvi - hc) 495 # // 496 # qv0[0] = zv0 + hc + alpha*(hv0-hc); 497 # qv1[0] = zv1 + hc + alpha*(hv1-hc); 498 # qv2[0] = zv2 + hc + alpha*(hv2-hc); 499 500 501 # //Update momentum as a linear combination of 502 # //qc[j] (shallow) and momentum from gradient limiter (deep). 503 # for (j=1; j<number_of_conserved_quantities; j++) { 504 # qv0[j] = (1-alpha)*qc[j] + alpha*qv0[j]; 505 # qv1[j] = (1-alpha)*qc[j] + alpha*qv1[j]; 506 # qv2[j] = (1-alpha)*qc[j] + alpha*qv2[j]; 507 # } 508 # } 509 510 511 # //Finally, protect against infinitesimal heights and high speeds 512 # hv0 = qv0[0]-zv0; hv1 = qv1[0]-zv1; hv2 = qv2[0]-zv2; 513 # hmax = max(hv0, max(hv1, hv2)); 514 # if (hmax < minimum_allowed_height) { 515 516 # //Set negative heights to bed elevation 517 # if (hc < 0.0) {qc[0] = zc; hc = 0.0;} 518 # if (hv0 < 0.0) {qv0[0] = zv0; hv0 = 0.0;} 519 # if (hv1 < 0.0) {qv1[0] = zv1; hv1 = 0.0;} 520 # if (hv2 < 0.0) {qv2[0] = zv2; hv2 = 0.0;} 521 522 523 # //FIXME: Do we need this desperately?? 524 # //FIXME: We definitely do encounter negative heights, 525 # //but what about momentum? 526 527 # //Control momentum 528 # //qc[1] = 0.0; qc[2] = 0.0; 529 # //qv0[1] = 0.0; qv0[2] = 0.0; 530 # //qv1[1] = 0.0; qv1[2] = 0.0; 531 # //qv2[1] = 0.0; qv2[2] = 0.0; 532 # } 533 # } 534 535 536 537 538 def limiter_org(k): 539 """limiter function, specific to the shallow water wave equation. 540 It will ensure that h is always non-negative. 541 This version works on the consecutive data structure 542 """ 543 544 #This version should be used to reproduce the bad limiting when levels 545 #are close to bed elevation 546 547 import sys 548 from config import beta 549 from Numeric import ones, Float 550 551 #conserved quantities at centroid 552 qc = Volume.conserved_quantities_centroid[k, :] 553 N = len(qc) 554 555 zc = Volume.field_values_centroid[k,0] #bed elevation at centroid 556 wc = qc[0] #Stage at centroid 557 hc = wc - zc #Water depth at centroid 558 559 #conserved quant. at vertices 560 qv0 = Volume.conserved_quantities_vertex0[k,:] 561 qv1 = Volume.conserved_quantities_vertex1[k,:] 562 qv2 = Volume.conserved_quantities_vertex2[k,:] 563 564 565 #Bed elevation at vertices 566 zv0 = Volume.field_values_vertex0[k,0] 567 zv1 = Volume.field_values_vertex1[k,0] 568 zv2 = Volume.field_values_vertex2[k,0] 569 570 #deltas 571 dq0 = qv0-qc 572 dq1 = qv1-qc 573 dq2 = qv2-qc 574 575 576 #------------------------------------- 577 hmin = hc 578 qmin = qmax = qc 579 for i in range(3): 580 n = Volume.neighbours[k,i] 581 if n >= 0: 582 qn = Volume.conserved_quantities_centroid[n,:] 583 zn = Volume.field_values_centroid[n,0] 584 hn = qn[0]-zn 585 586 hmin = min(hmin, hn) 587 qmin = minimum(qmin, qn) 588 qmax = maximum(qmax, qn) 589 590 hmin = max(hmin,0.0) 591 592 593 #----------------------------------- 594 # First run - limit on conserved quantities 595 # based on neighbouring heights 596 #----------------------------------- 597 598 dqmax = qmax - qc 599 dqmin = qmin - qc 600 phi = ones(len(qc), Float) 601 for j in range(N): 602 #First find the gradient limiter (phi) 603 604 #Vertex0: 605 r = 1.0 606 if dq0[j] > 0: r = dqmax[j]/dq0[j] 607 if dq0[j] < 0: r = dqmin[j]/dq0[j] 608 phi[j] = min( min(r*beta, 1), phi[j] ) 609 610 #Vertex1: 611 r = 1.0 612 if dq1[j] > 0: r = dqmax[j]/dq1[j] 613 if dq1[j] < 0: r = dqmin[j]/dq1[j] 614 phi[j] = min( min(r*beta, 1), phi[j] ) 615 616 #Vertex2: 617 r = 1.0 618 if dq2[j] > 0: r = dqmax[j]/dq2[j] 619 if dq2[j] < 0: r = dqmin[j]/dq2[j] 620 phi[j] = min( min(r*beta, 1), phi[j] ) 621 622 #The update using phi limiter 623 qv0[j] = qc[j] + phi[j]*dq0[j] #Vertex0 624 qv1[j] = qc[j] + phi[j]*dq1[j] #Vertex1 625 qv2[j] = qc[j] + phi[j]*dq2[j] #Vertex2 626 627 628 629 #----------------------------------- 630 # Second run - lower limit on height 631 # based on neighbouring heights for 632 # surfaces constructed below bed 633 #----------------------------------- 634 635 #Adjusted water heights at vertices 636 hv0 = qv0[0]-zv0 637 hv1 = qv1[0]-zv1 638 hv2 = qv2[0]-zv2 639 640 #Deltas 641 dhmin = hmin - hc 642 dh0 = hv0 - hc 643 dh1 = hv1 - hc 644 dh2 = hv2 - hc 645 646 647 #Work out if update is needed and calculate phi limiter 648 phi = 1.0 649 update_needed = False 650 if hv0 < 0.0: 651 update_needed = True 652 r = 1.0 653 if dh0 < 0: r = dhmin/dh0 654 phi = min( min(beta*r, 1), phi ) 655 656 if hv1 < 0.0: 657 update_needed = True 658 r = 1.0 659 if dh1 < 0: r = dhmin/dh1 660 phi = min( min(beta*r, 1), phi ) 661 662 if hv2 < 0.0: 663 update_needed = True 664 r = 1.0 665 if dh2 < 0: r = dhmin/dh2 666 phi = min( min(beta*r, 1), phi ) 667 668 669 #Update stage at vertices 670 if update_needed: 671 hv0 = hc + phi*dh0 672 if hv0 < 0.0: raise 'hv0 = %12f hc = %12f ' % (hv0 , hc) 673 qv0[0] = zv0 + hv0 674 675 hv1 = hc + phi*dh1 676 if hv1 < 0.0: raise 'hv1 = %12f hc = %12f ' % (hv1 , hc) 677 qv1[0] = zv1 + hv1 678 679 hv2 = hc + phi*dh2 680 if hv2 < 0.0: raise 'hv2 = %12f hc = %12f ' % (hv2 , hc) 681 qv2[0] = zv2 + hv2 682 683 684 #In either case protect against infinitesimal heights 685 dry(k) 686 687 def dry(k): 688 """Protect agains infinitesimal heights (and high speeds) 689 There are no states for any conserved quantity when h is 690 between minimum allowed height and 0. 691 """ 692 693 #volume = Volume.instances[k] 694 695 maxhv = 0.0 696 697 wv0 = Volume.conserved_quantities_vertex0[k,0] 698 zv0 = Volume.field_values_vertex0[k,0] 699 hv0 = wv0-zv0 700 701 wv1 = Volume.conserved_quantities_vertex1[k,0] 702 zv1 = Volume.field_values_vertex1[k,0] 703 hv1 = wv1-zv1 704 705 wv2 = Volume.conserved_quantities_vertex2[k,0] 706 zv2 = Volume.field_values_vertex2[k,0] 707 hv2 = wv2-zv2 708 709 maxhv = max(hv0, hv1, hv2) 710 711 if maxhv < minimum_allowed_height: 712 #Set water level to bed elevation and zero momentums. 713 Volume.conserved_quantities_vertex0[k,:] =\ 714 array([Volume.field_values_vertex0[k,0], 0, 0]) 715 716 Volume.conserved_quantities_vertex1[k,:] =\ 717 array([Volume.field_values_vertex1[k,0], 0, 0]) 718 719 Volume.conserved_quantities_vertex2[k,:] =\ 720 array([Volume.field_values_vertex2[k,0], 0, 0]) 721 722 #Set only momentums to zero 723 #Volume.conserved_quantities_vertex0[k,1:] = array([0., 0.]) 724 #Volume.conserved_quantities_vertex1[k,1:] = array([0., 0.]) 725 #Volume.conserved_quantities_vertex2[k,1:] = array([0., 0.]) 726 727 728 729 730 def distribute_to_vertices_and_edges(domain): 731 """Extrapolate conserved quantities from centroid to 732 vertices and edge-midpoints for each volume 733 """ 734 735 from config import minimum_allowed_height 736 737 from mesh import Point 738 order = domain.order 739 740 limiter = domain.limiter 741 if limiter is None: 742 from shallow_water import limit_on_w_then_h as limiter 743 744 745 for k in range(Volume.number_of_instances): 746 747 #conserved quantities at centroid 748 q = Volume.conserved_quantities_centroid[k,:] 749 750 if order == 1: 751 ######################## 752 #NOTE: This is now specific to the shallow water wave equation 753 754 zc = Volume.field_values_centroid[k,0] #bed elevation at centroid 755 wc = q[0] #Stage at centroid 756 hc = wc - zc #Water depth at centroid 757 758 #Protect against infinitesimal heights and high velocities 759 if hc < minimum_allowed_height: 760 #Set momentums to zero 761 q[0] = zc 762 q[1] = 0.0 763 q[2] = 0.0 764 765 766 #Bed elevation at vertices 767 zv0 = Volume.field_values_vertex0[k,0] 768 zv1 = Volume.field_values_vertex1[k,0] 769 zv2 = Volume.field_values_vertex2[k,0] 770 771 #Update conserved quant. at vertices 772 Volume.conserved_quantities_vertex0[k,1:] = q[1:] #Momentum 773 Volume.conserved_quantities_vertex1[k,1:] = q[1:] #Momentum 774 Volume.conserved_quantities_vertex2[k,1:] = q[1:] #Momentum 775 776 #Stage 777 z_range = max( max(abs(zv0-zv1), abs(zv0-zv2)), abs(zv1-zv2) ) 778 if z_range>0: 779 alpha = min(hc/z_range, 1.0) 780 else: 781 alpha = 1 782 783 zz = hc+alpha*zc 784 Volume.conserved_quantities_vertex0[k,0] = zz + (1-alpha)*zv0 785 Volume.conserved_quantities_vertex1[k,0] = zz + (1-alpha)*zv1 786 Volume.conserved_quantities_vertex2[k,0] = zz + (1-alpha)*zv2 787 788 elif order == 2: 789 a, b = compute_gradient(k) 790 791 #Compute extrapolated values at vertices 792 x,y = Point.coordinates[Volume.points[k,3]] #Centroid 793 794 #Get vertex coordinates 795 i = Volume.points[k,0] 796 x0 = Point.coordinates[i,0] 797 y0 = Point.coordinates[i,1] 798 799 i = Volume.points[k,1] 800 x1 = Point.coordinates[i,0] 801 y1 = Point.coordinates[i,1] 802 803 i = Volume.points[k,2] 804 x2 = Point.coordinates[i,0] 805 y2 = Point.coordinates[i,1] 806 807 808 #Extrapolate 809 Volume.conserved_quantities_vertex0[k, :] = q + a*(x0-x) + b*(y0-y) 810 Volume.conserved_quantities_vertex1[k, :] = q + a*(x1-x) + b*(y1-y) 811 Volume.conserved_quantities_vertex2[k, :] = q + a*(x2-x) + b*(y2-y) 812 813 814 if limiter is not None: 815 limiter(k) 816 else: 817 raise 'Unknown order' 818 819 820 #Compute conserved quantities as interpolated 821 #quantities at each edge midpoint 822 q0 = Volume.conserved_quantities_vertex0[k,:] 823 q1 = Volume.conserved_quantities_vertex1[k,:] 824 q2 = Volume.conserved_quantities_vertex2[k,:] 825 826 Volume.conserved_quantities_face0[k,:] = 0.5*(q1 + q2) 827 Volume.conserved_quantities_face1[k,:] = 0.5*(q0 + q2) 828 Volume.conserved_quantities_face2[k,:] = 0.5*(q0 + q1) 519 wc[k] = zc[k] 520 ###hc[k] = 0.0 521 for i in range(3): 522 if hv[k,i] < 0.0: 523 wv[k,i] = zv[k,i] 524 ##hv0 = 0.0;} 525 526 829 527 830 528 … … 869 567 870 568 569 ######################### 570 #Standard forcing terms: 571 572 def gravity(domain): 573 """Implement forcing function for bed slope working with 574 consecutive data structures of class Volume 575 """ 576 577 from config import g 578 from util import gradient 579 from Numeric import zeros, Float, array, sum 580 581 Level = domain.quantities['level'] 582 Xmom = domain.quantities['xmomentum'] 583 Ymom = domain.quantities['ymomentum'] 584 585 Elevation = domain.quantities['elevation'] 586 h = Level.edge_values - Elevation.edge_values 587 V = domain.get_vertex_coordinates() 588 589 for k in range(domain.number_of_elements): 590 avg_h = sum( h[k,:] )/3 591 592 #Compute bed slope 593 x0, y0, x1, y1, x2, y2 = V[k,:] 594 z0, z1, z2 = Elevation.vertex_values[k,:] 595 zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 596 597 #Update momentum 598 Xmom.explicit_update[k] += -g*zx*avg_h 599 Ymom.explicit_update[k] += -g*zy*avg_h 600 601 602 def manning_friction(domain): 603 """Apply (Manning) friction to water momentum 604 """ 605 606 import Numeric 607 from math import sqrt 608 from config import g, minimum_allowed_height 609 610 Level = domain.quantities['level'] 611 Xmom = domain.quantities['xmomentum'] 612 Ymom = domain.quantities['ymomentum'] 613 614 Friction = domain.quantities['friction'] 615 616 for k in range(domain.number_of_elements): 617 w = Level.centroid_values[k] 618 uh = Xmom.centroid_values[k] 619 vh = Ymom.centroid_values[k] 620 manning = Friction.centroid_values[k] 621 622 if w >= minimum_allowed_height: 623 S = -g * manning**2 * sqrt((uh**2 + vh**2)) 624 S /= w**(7.0/3) 625 626 #Update momentum 627 Xmom.explicit_update[k] += S*uh 628 Ymom.explicit_update[k] += S*vh 629 871 630 872 631 ############################################## -
inundation/ga/storm_surge/pyvolution/test_quantity.py
r198 r205 148 148 quantity = Quantity(self.mesh4) 149 149 150 #Set up for a gradient of ( 1,0)151 quantity.set_values([2.0 /3, 4.0/3, 8.0/3, 2.0/3],150 #Set up for a gradient of (3,0) at mid triangle 151 quantity.set_values([2.0, 4.0, 8.0, 2.0], 152 152 location = 'centroids') 153 153 154 154 #Gradients 155 quantity.compute_gradients() 156 157 #FIXME: 158 159 #HERTIL 160 161 #Gradient of fitted pwl surface 162 #a, b = compute_gradient(v2.id) 163 164 #assert a == 1.0 165 #assert b == 0.0 166 167 168 #And now for the second order stuff 169 # - the full second order extrapolation 170 #domain.order = 2 171 #domain.limiter = dummy_limiter 172 #distribute_to_vertices_and_edges(domain) 173 174 #assert v2.conserved_quantities_vertex0 == 0.0 175 #assert v2.conserved_quantities_vertex1 == 2.0 176 #assert v2.conserved_quantities_vertex2 == 2.0 155 a, b = quantity.compute_gradients() 156 157 158 #gradient bewteen t0 and t1 is undefined as det==0 159 assert a[0] == 0.0 160 assert b[0] == 0.0 161 #The others are OK 162 for i in range(1,4): 163 assert a[i] == 3.0 164 assert b[i] == 0.0 165 166 167 quantity.second_order_extrapolator() 168 169 assert allclose(quantity.vertex_values, [[2., 2., 2.], 170 [0., 6., 6.], 171 [6., 6., 12.], 172 [0., 0., 6.]]) 177 173 178 174 … … 355 351 356 352 357 def test_first_order_ limiter(self):353 def test_first_order_extrapolator(self): 358 354 quantity = Quantity(self.mesh4) 359 355 … … 363 359 364 360 #Extrapolate 365 quantity.first_order_ limiter()361 quantity.first_order_extrapolator() 366 362 367 363 #Check vertices but not edge values … … 370 366 371 367 372 def test_second_order_limiter(self): 373 quantity = Quantity(self.mesh4) 374 375 #Create a deliberate overshoot (e.g. from gradient computation) 376 quantity.set_values([[0,3,3], [6,2,2], [3,8,5], [3,5,8]]) 377 378 #Limit and extrapolate 379 quantity.second_order_limiter() 368 def test_second_order_extrapolator(self): 369 quantity = Quantity(self.mesh4) 370 371 #Set up for a gradient of (3,0) at mid triangle 372 quantity.set_values([2.0, 4.0, 8.0, 2.0], 373 location = 'centroids') 374 375 376 377 quantity.second_order_extrapolator() 378 quantity.limiter() 380 379 381 380 382 381 #Assert that central triangle is limited by neighbours 383 assert quantity.vertex_values[1,0] <= quantity.vertex_values[2,2]384 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]385 386 assert quantity.vertex_values[1,1] <= quantity.vertex_values[ 3,0]382 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 383 assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1] 384 385 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 387 386 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 388 387 389 388 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 390 assert quantity.vertex_values[1,2] >= quantity.vertex_values[ 0,1]389 assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1] 391 390 392 391 … … 398 397 399 398 400 #Check vertices but not edge values 401 #assert allclose(quantity.vertex_values, 402 # [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) 399 400 401 402 def test_limiter(self): 403 quantity = Quantity(self.mesh4) 404 405 #Create a deliberate overshoot (e.g. from gradient computation) 406 quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 407 408 409 #Limit 410 quantity.limiter() 411 412 #Assert that central triangle is limited by neighbours 413 assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0] 414 assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1] 415 416 assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1] 417 assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2] 418 419 assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0] 420 assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1] 421 422 423 424 #Assert that quantities are conserved 425 from Numeric import sum 426 for k in range(quantity.centroid_values.shape[0]): 427 assert allclose (quantity.centroid_values[k], 428 sum(quantity.vertex_values[k,:])/3) 429 403 430 404 431 … … 412 439 413 440 #Extrapolate 414 quantity.distribute_to_vertices_and_edges(order=1) 441 quantity.first_order_extrapolator() 442 443 #Interpolate 444 quantity.interpolate_from_vertices_to_edges() 415 445 416 446 assert allclose(quantity.vertex_values, -
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r198 r205 536 536 537 537 538 def test_first_order_ limiter_const_z(self):538 def test_first_order_extrapolator_const_z(self): 539 539 540 540 a = [0.0, 0.0] … … 564 564 565 565 566 domain.first_order_limiter() 566 domain.order = 1 567 domain.distribute_to_vertices_and_edges() 567 568 568 569 #Check that centroid values were distributed to vertices … … 608 609 #- so that the limiter has something to work with 609 610 assert not alltrue(alltrue(greater_equal(L,E-epsilon))) 610 611 domain.first_order_limiter() 611 612 domain.order = 1 613 domain.distribute_to_vertices_and_edges() 612 614 613 615 #Check that all levels are above elevation (within eps) -
inundation/ga/storm_surge/pyvolution/test_util.py
r195 r205 22 22 x2 = 0.0; y2 = 1.0; z2 = 0.0 23 23 24 zx, zy = gradient(x0, y0, x1, y1, x2, y2, array([z0]), array([z1]), array([z2]))24 zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 25 25 26 26 assert zx == -1.0 … … 50 50 x2 = 2.0/3; y2 = 8.0/3 51 51 52 q0 = array([2.0+2.0/3])53 q1 = array([8.0+2.0/3])54 q2 = array([2.0+8.0/3])52 q0 = 2.0+2.0/3 53 q1 = 8.0+2.0/3 54 q2 = 2.0+8.0/3 55 55 56 56 #Gradient of fitted pwl surface 57 57 a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) 58 58 59 assert abs(a [0]- 3.0) < epsilon60 assert abs(b [0]- 1.0) < epsilon59 assert abs(a - 3.0) < epsilon 60 assert abs(b - 1.0) < epsilon 61 61 62 62 … … 84 84 x2 = 2.0/3; y2 = 8.0/3 85 85 86 q0 = array([2.0+2.0/3])87 q1 = array([8.0+2.0/3])88 q2 = array([2.0+8.0/3])86 q0 = 2.0+2.0/3 87 q1 = 8.0+2.0/3 88 q2 = 2.0+8.0/3 89 89 90 90 #Gradient of fitted pwl surface 91 91 a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2) 92 92 93 assert abs(a[0] - 3.0) < epsilon 94 assert abs(b[0] - 1.0) < epsilon 95 96 def test_gradient_C_extension_scalars(self): 97 from util_ext import gradient as gradient_c 98 99 from RandomArray import uniform, seed 100 seed(17, 53) 101 102 x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6) 103 104 q0 = 10.0 105 q1 = 3.0 106 q2 = 7.0 107 108 109 #Gradient of fitted pwl surface 110 from util import gradient_python 111 a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2) 112 113 #print a_ref, b_ref 114 try: 115 a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2) 116 except TypeError: 117 pass 118 else: 119 raise 'gradient c-extension should throw exception when q''s are scalars' 93 assert abs(a - 3.0) < epsilon 94 assert abs(b - 1.0) < epsilon 120 95 121 96 … … 132 107 q2 = uniform(7.0, 20.0, 4) 133 108 134 135 #Gradient of fitted pwl surface136 from util import gradient_python137 a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)138 109 139 #print a_ref, b_ref 140 a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2) 110 for i in range(4): 111 #Gradient of fitted pwl surface 112 from util import gradient_python 113 a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, 114 q0[i], q1[i], q2[i]) 141 115 142 #print a, a_ref, b, b_ref 143 for i in range(4): 144 assert abs(a[i] - a_ref[i]) < epsilon 145 assert abs(b[i] - b_ref[i]) < epsilon 116 #print a_ref, b_ref 117 a, b = gradient_c(x0, y0, x1, y1, x2, y2, 118 q0[i], q1[i], q2[i]) 119 120 #print a, a_ref, b, b_ref 121 assert abs(a - a_ref) < epsilon 122 assert abs(b - b_ref) < epsilon 146 123 147 124 -
inundation/ga/storm_surge/pyvolution/util_ext.c
r198 r205 18 18 double x1, double y1, 19 19 double x2, double y2, 20 double *q0, double *q1, double *q2, 21 double *a, double *b, 22 int N) { 20 double q0, double q1, double q2, 21 double *a, double *b) { 23 22 24 23 double det; 25 int i;26 24 27 25 det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0); 28 29 for (i=0; i<N; i++) {30 a[i] = (y2-y0)*(q1[i]-q0[i]) - (y1-y0)*(q2[i]-q0[i]);31 a[i] /= det;32 26 33 b[i] = (x1-x0)*(q2[i]-q0[i]) - (x2-x0)*(q1[i]-q0[i]); 34 b[i] /= det; 35 } 27 *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0); 28 *a /= det; 29 30 *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0); 31 *b /= det; 36 32 37 33 return 0; … … 81 77 // 82 78 83 double x0, y0, x1, y1, x2, y2; 84 PyObject *Q0, *Q1, *Q2, *result; 85 PyArrayObject *q0, *q1, *q2, *a, *b; 86 int dimensions[1]; 87 int N; 88 79 double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b; 80 PyObject *result; 89 81 90 82 // Convert Python arguments to C 91 if (!PyArg_ParseTuple(args, "dddddd OOO", &x0, &y0, &x1, &y1, &x2, &y2,92 & Q0, &Q1, &Q2))83 if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2, 84 &q0, &q1, &q2)) 93 85 return NULL; 94 95 //Input check96 if PyArray_Check(Q0) {97 q0 = (PyArrayObject *)98 PyArray_ContiguousFromObject(Q0, PyArray_DOUBLE, 0, 0);99 } else {100 PyErr_SetString(PyExc_TypeError, "q0 must be a numeric array");101 return NULL;102 }103 104 if PyArray_Check(Q1) {105 q1 = (PyArrayObject *)106 PyArray_ContiguousFromObject(Q1, PyArray_DOUBLE, 0, 0);107 } else {108 PyErr_SetString(PyExc_TypeError, "q1 must be a numeric array");109 return NULL;110 }111 112 if PyArray_Check(Q2) {113 q2 = (PyArrayObject *)114 PyArray_ContiguousFromObject(Q2, PyArray_DOUBLE, 0, 0);115 } else {116 PyErr_SetString(PyExc_TypeError, "q2 must be a numeric array");117 return NULL;118 }119 120 //Allocate space for return vectors a and b121 122 N = q0 -> dimensions[0];123 dimensions[0] = N;124 a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);125 b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);126 86 127 87 128 88 // Call underlying routine 129 89 _gradient(x0, y0, x1, y1, x2, y2, 130 (double *) q0 -> data, 131 (double *) q1 -> data, 132 (double *) q2 -> data, 133 (double *) a -> data, 134 (double *) b-> data, N); 90 q0, q1, q2, &a, &b); 135 91 136 Py_DECREF(q0); 137 Py_DECREF(q1); 138 Py_DECREF(q2); 139 140 // Return arrays a and b 141 result = Py_BuildValue("OO", a, b); 142 Py_DECREF(a); 143 Py_DECREF(b); 92 // Return values a and b 93 result = Py_BuildValue("dd", a, b); 144 94 return result; 145 95 }
Note: See TracChangeset
for help on using the changeset viewer.
