Changeset 4823
 Timestamp:
 Nov 16, 2007, 5:17:49 PM (16 years ago)
 Location:
 anuga_core/source/anuga/shallow_water
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4815 r4823 464 464 */ 465 465 466 467 468 469 double velocity_balance(double uh_i, 470 double uh, 471 double h_i, 472 double h, 473 double alpha, 474 double epsilon) { 475 // Find alpha such that speed at the vertex is within one 476 // order of magnitude of the centroid speed 477 478 // FIXME(Ole): Work in progress 479 480 double a, b, estimate; 481 double m=10; // One order of magnitude  allow for velocity deviations at vertices 482 483 484 printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n", 485 alpha, uh_i, uh, h_i, h); 486 487 488 489 490 // Shorthands and determine inequality 491 if (fabs(uh) < epsilon) { 492 a = 1.0e10; // Limit 493 } else { 494 a = fabs(uh_i  uh)/fabs(uh); 495 } 496 497 if (h < epsilon) { 498 b = 1.0e10; // Limit 499 } else { 500 b = m*fabs(h_i  h)/h; 501 } 502 503 printf("a %f, b %f\n", a, b); 504 505 if (a > b) { 506 estimate = (m1)/(ab); 507 508 printf("Alpha %f, estimate %f\n", 509 alpha, estimate); 510 511 if (alpha < estimate) { 512 printf("Adjusting alpha from %f to %f\n", 513 alpha, estimate); 514 alpha = estimate; 515 } 516 } else { 517 518 if (h < h_i) { 519 estimate = (m1)/(ab); 520 521 printf("Alpha %f, estimate %f\n", 522 alpha, estimate); 523 524 if (alpha < estimate) { 525 printf("Adjusting alpha from %f to %f\n", 526 alpha, estimate); 527 alpha = estimate; 528 } 529 } 530 // Always fulfilled as alpha and m1 are non negative 531 } 532 533 534 return alpha; 535 } 536 537 466 538 int _balance_deep_and_shallow(int N, 467 539 double beta_h, … … 480 552 481 553 int k, k3, i, excessive_froude_number=0; 554 482 555 double dz, hmin, alpha, h_diff, hc_k; 483 556 double epsilon = 1.0e6; // FIXME: Temporary measure … … 539 612 540 613 // Make alpha as large as possible but still ensure that 541 // final depth is positive 614 // final depth is positive and that velocities at vertices 615 // are controlled 542 616 543 617 if (hmin < H0) { 544 618 alpha = 1.0; 545 619 for (i=0; i<3; i++) { 546 547 // FIXME (Ole): Simplify when (if) hvbar gets retired620 621 // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired 548 622 if (beta_h > epsilon) { 549 h_diff = hvbar[k3+i]  hv[i]; 550 } else { 551 h_diff = hc_k  hv[i]; 623 hc_k = hvbar[k3+i]; // Depth to be used at vertices 552 624 } 553 625 626 h_diff = hc_k  hv[i]; 554 627 if (h_diff <= 0) { 555 628 // Deep water triangle is further away from bed than … … 560 633 // hlimited stage. 561 634 562 // FIXME (Ole): Simplify when (if) hvbar gets retired 563 if (beta_h > epsilon) { 564 alpha = min(alpha, (hvbar[k3+i]  H0)/h_diff); 565 } else { 566 alpha = min(alpha, (hc_k  H0)/h_diff); 567 } 635 alpha = min(alpha, (hc_k  H0)/h_diff); 568 636 } 569 637 } … … 574 642 575 643 } else { 576 // Use wlimited stage exclusively 644 // Use wlimited stage exclusively in deeper water. 577 645 alpha = 1.0; 578 646 } 647 648 649 /* 650 // Experimental code for controlling velocities at vertices. 651 // Adjust alpha (down towards first order) such that 652 // velocities at vertices remain within one order of magnitude 653 // of those at the centroid. 654 655 for (i=0; i<3; i++) { 656 657 // FIXME (Ole): Simplify (remove) when (if) hvbar gets retired 658 if (beta_h > epsilon) { 659 hc_k = hvbar[k3+i]; // Depth to be used at vertices 660 } 661 662 663 alpha = velocity_balance(xmomv[k3+i], xmomc[k], 664 hv[i], hc_k, alpha, epsilon); 665 666 alpha = velocity_balance(ymomv[k3+i], ymomc[k], 667 hv[i], hc_k, alpha, epsilon); 668 } 669 */ 579 670 } 580 581 582 671 583 672 //printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n", 584 673 // k, hmin, dz, alpha, alpha_balance); … … 624 713 625 714 715 //if (0) { // FIXME(Ole): Disabled while testing balancing of velocities above 626 716 if (tight_slope_limiters == 1) { 627 717 
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r4815 r4823 2861 2861 u = xmomentum/depth 2862 2862 v = ymomentum/depth 2863 2864 # Assert that all vertex velocities stay within one 2865 # order of magnitude of centroid velocities. 2866 #print u.vertex_values[1,:] 2867 #print u.centroid_values[1] 2868 2869 assert alltrue(absolute(u.vertex_values[1,:]) <= absolute(u.centroid_values[1])*10) 2870 assert alltrue(absolute(v.vertex_values[1,:]) <= absolute(v.centroid_values[1])*10) 2863 2871 2864 2872 denom = (depth*g)**0.5
Note: See TracChangeset
for help on using the changeset viewer.