Changeset 8400
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r8390 r8400 631 631 // Innermost flux function (using stage w=z+h) 632 632 633 int _flux_function_central_wb_3(struct edge *E_left, struct edge *E_right, 634 double n1, double n2, 633 int _flux_function_central_wb_3( 634 struct edge *E_left, 635 struct edge *E_right, 636 double n1, 637 double n2, 635 638 double epsilon, 636 639 double h0, … … 657 660 double h1_left, h2_left, h1_right, h2_right; 658 661 double uh1_left, uh2_left, uh1_right, uh2_right; 662 double vh1_left, vh2_left, vh1_right, vh2_right; 659 663 double u1_left, u2_left, u1_right, u2_right; 660 664 double p_left, p_right; … … 674 678 q_left[1] = E_left->uh; 675 679 q_left[2] = E_left->vh; 676 677 680 678 681 … … 736 739 h1_left = E_left->h1; 737 740 uh1_left = E_left->uh1; 741 vh1_left = E_left->vh1; 738 742 u1_left = _compute_speed(&uh1_left, &h1_left, epsilon, h0, limiting_threshold); 739 743 … … 741 745 h2_left = E_left->h2; 742 746 uh2_left = E_left->uh2; 747 vh2_left = E_left->vh2; 743 748 u2_left = _compute_speed(&uh2_left, &h2_left, epsilon, h0, limiting_threshold); 744 749 … … 746 751 h1_right = E_right->h2; 747 752 uh1_right = E_right->uh2; 753 vh1_right = E_right->vh2; 748 754 u1_right = _compute_speed(&uh1_right, &h1_right, epsilon, h0, limiting_threshold); 749 755 … … 751 757 h2_right = E_right->h1; 752 758 uh2_right = E_right->uh1; 759 vh2_right = E_right->vh1; 753 760 u2_right = _compute_speed(&uh2_right, &h2_right, epsilon, h0, limiting_threshold); 754 761 } else { … … 757 764 h1_left = h_left; 758 765 uh1_left = uh_left; 766 vh1_left = vh_left; 759 767 u1_left = u_left; 760 768 … … 762 770 h2_left = h_left; 763 771 uh2_left = uh_left; 772 vh2_left = vh_left; 764 773 u2_left = u_left; 765 774 766 775 // First vertex right edge data (needs interchange) 767 h1_right = h_left; 768 uh1_right = uh_left; 769 u1_right = u_left; 776 h1_right = h_right; 777 uh1_right = uh_right; 778 vh1_right = vh_right; 779 u1_right = u_right; 770 780 771 781 // Second vertex right edge data (needs interchange) 772 h2_right = h_left; 773 uh2_right = uh_left; 774 u2_right = u_left; 775 } 776 777 778 779 780 782 h2_right = h_right; 783 uh2_right = uh_right; 784 vh2_right = vh_right; 785 u2_right = u_right; 786 } 781 787 782 788 //printf("h1_left, h2_left %g %g \n", h1_left, h2_left); 783 789 //printf("h1_right, h2_right %g %g \n", h1_right, h2_right); 784 790 785 786 791 p_left = 0.5*g*h_left*h_left; 787 792 p_right = 0.5*g*h_right*h_right; … … 795 800 796 801 // Flux formulas 797 flux_left[0] = u_left*h_left; 798 flux_left[1] = u_left * uh_left + p_left; 799 flux_left[2] = u_left*vh_left; 800 802 //flux_left[0] = u_left*h_left; 803 //flux_left[1] = u_left * uh_left + p_left; 804 //flux_left[2] = u_left*vh_left; 805 806 807 flux_left[0] = (u1_left*h1_left + 4.0*u_left*h_left + u2_left*h2_left)/6.0; 808 flux_left[1] = (u1_left*uh1_left + 4.0*u_left*uh_left + u2_left*uh2_left)/6.0 + p_left; 809 flux_left[2] = (u1_left*vh1_left + 4.0*u_left*vh_left + u2_left*vh2_left)/6.0; 801 810 802 811 803 812 //printf("flux_left %g %g %g \n",flux_left[0],flux_left[1],flux_left[2]); 804 813 805 flux_right[0] = u_right*h_right; 806 flux_right[1] = u_right*uh_right + p_right; 807 flux_right[2] = u_right*vh_right; 814 //flux_right[0] = u_right*h_right; 815 //flux_right[1] = u_right*uh_right + p_right; 816 //flux_right[2] = u_right*vh_right; 817 818 819 flux_right[0] = (u1_right*h1_right + 4.0*u_right*h_right + u2_right*h2_right)/6.0; 820 flux_right[1] = (u1_right*uh1_right + 4.0*u_right*uh_right + u2_right*uh2_right)/6.0 + p_right; 821 flux_right[2] = (u1_right*vh1_right + 4.0*u_right*vh_right + u2_right*vh2_right)/6.0; 808 822 809 823 //printf("flux_right %g %g %g \n",flux_right[0],flux_right[1],flux_right[2]); … … 3378 3392 // This assumes compute_fluxes called before forcing terms 3379 3393 memset((char*) D->stage_explicit_update, 0, D->number_of_elements * sizeof (double)); 3380 memset((char*) D->xmom_explicit_update, 0, D->number_of_elements * sizeof (double));3381 memset((char*) D->ymom_explicit_update, 0, D->number_of_elements * sizeof (double));3394 memset((char*) D->xmom_explicit_update, 0, D->number_of_elements * sizeof (double)); 3395 memset((char*) D->ymom_explicit_update, 0, D->number_of_elements * sizeof (double)); 3382 3396 3383 3397 // For all triangles … … 3393 3407 } 3394 3408 3395 // Get the "left" values at the edge vertices and midpo nt3409 // Get the "left" values at the edge vertices and midpoint 3396 3410 // from the triangle k, edge i 3397 3411 get_edge_data(&E_left, D, k, i); … … 3425 3439 E_right.z1 = E_left.z; 3426 3440 3427 3428 3441 E_right.w2 = E_right.w; 3429 3442 E_right.h2 = E_right.h; … … 3433 3446 E_right.z2 = E_left.z; 3434 3447 3435 3436 3437 3438 3439 3440 3441 3442 3448 } else { 3443 3449 // Neighbour is a real triangle 3444 3450 m = D->neighbour_edges[k3i]; 3451 n3m = 3*n+m; 3445 3452 3446 3453 get_edge_data(&E_right, D, n, m); … … 3500 3507 // Update triangle k with flux from edge i 3501 3508 D->stage_explicit_update[k] -= edgeflux[0]; 3502 D->xmom_explicit_update[k] -= edgeflux[1];3503 D->ymom_explicit_update[k] -= edgeflux[2];3509 D->xmom_explicit_update[k] -= edgeflux[1]; 3510 D->ymom_explicit_update[k] -= edgeflux[2]; 3504 3511 3505 3512 D->already_computed_flux[k3i] = call; // #k Done … … 3510 3517 if (n >= 0) { 3511 3518 3512 n3m = 3 * n + (m + 1) % 3;3513 3519 D->stage_explicit_update[n] += edgeflux[0]; 3514 D->xmom_explicit_update[n] += edgeflux[1];3515 D->ymom_explicit_update[n] += edgeflux[2];3516 3517 n3m = 3*n+m; 3520 D->xmom_explicit_update[n] += edgeflux[1]; 3521 D->ymom_explicit_update[n] += edgeflux[2]; 3522 3523 3518 3524 D->already_computed_flux[n3m] = call; // #n Done 3519 3525 } … … 3976 3982 stage, xmomentum and ymomentum 3977 3983 3978 The3979 3980 3984 The maximal allowable speed computed by the flux_function for each volume 3981 3985 is converted to a timestep that must not be exceeded. The minimum of -
trunk/anuga_work/shallow_water_balanced_steve/run_sw_merimbula.py
r8395 r8400 96 96 domain.set_timestepping_method('rk2') 97 97 #domain.set_CFL(0.7) 98 #domain.set_beta(1.5)98 domain.set_beta(1.5) 99 99 domain.set_extrapolate_velocity(True) 100 100 domain.set_compute_fluxes_method('wb_3') 101 domain.set_extrapolate_velocity(True) 101 102 domain.set_name('merimbula') 102 103
Note: See TracChangeset
for help on using the changeset viewer.