Changeset 8400


Ignore:
Timestamp:
Apr 12, 2012, 9:27:36 PM (12 years ago)
Author:
steve
Message:

new well balanced scheme

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8390 r8400  
    631631// Innermost flux function (using stage w=z+h)
    632632
    633 int _flux_function_central_wb_3(struct edge *E_left, struct edge *E_right,
    634         double n1, double n2,
     633int _flux_function_central_wb_3(
     634        struct edge *E_left,
     635        struct edge *E_right,
     636        double n1,
     637        double n2,
    635638        double epsilon,
    636639        double h0,
     
    657660    double h1_left, h2_left, h1_right, h2_right;
    658661    double uh1_left, uh2_left, uh1_right, uh2_right;
     662    double vh1_left, vh2_left, vh1_right, vh2_right;
    659663    double u1_left, u2_left, u1_right, u2_right;
    660664    double p_left, p_right;
     
    674678    q_left[1] = E_left->uh;
    675679    q_left[2] = E_left->vh;
    676 
    677680
    678681
     
    736739        h1_left = E_left->h1;
    737740        uh1_left = E_left->uh1;
     741        vh1_left = E_left->vh1;
    738742        u1_left = _compute_speed(&uh1_left, &h1_left, epsilon, h0, limiting_threshold);
    739743
     
    741745        h2_left = E_left->h2;
    742746        uh2_left = E_left->uh2;
     747        vh2_left = E_left->vh2;
    743748        u2_left = _compute_speed(&uh2_left, &h2_left, epsilon, h0, limiting_threshold);
    744749
     
    746751        h1_right = E_right->h2;
    747752        uh1_right = E_right->uh2;
     753        vh1_right = E_right->vh2;
    748754        u1_right = _compute_speed(&uh1_right, &h1_right, epsilon, h0, limiting_threshold);
    749755
     
    751757        h2_right = E_right->h1;
    752758        uh2_right = E_right->uh1;
     759        vh2_right = E_right->vh1;
    753760        u2_right = _compute_speed(&uh2_right, &h2_right, epsilon, h0, limiting_threshold);
    754761    } else {
     
    757764        h1_left  = h_left;
    758765        uh1_left = uh_left;
     766        vh1_left = vh_left;
    759767        u1_left  = u_left;
    760768
     
    762770        h2_left  = h_left;
    763771        uh2_left = uh_left;
     772        vh2_left = vh_left;
    764773        u2_left  = u_left;
    765774
    766775        // 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;
    770780
    771781        // 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    }
    781787
    782788    //printf("h1_left, h2_left %g %g \n", h1_left, h2_left);
    783789    //printf("h1_right, h2_right %g %g \n", h1_right, h2_right);
    784790
    785 
    786791    p_left  = 0.5*g*h_left*h_left;
    787792    p_right = 0.5*g*h_right*h_right;
     
    795800
    796801    // 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;
    801810
    802811
    803812    //printf("flux_left %g %g %g \n",flux_left[0],flux_left[1],flux_left[2]);
    804813
    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;
    808822
    809823    //printf("flux_right %g %g %g \n",flux_right[0],flux_right[1],flux_right[2]);
     
    33783392    // This assumes compute_fluxes called before forcing terms
    33793393    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));
    33823396
    33833397    // For all triangles
     
    33933407            }
    33943408
    3395             // Get the "left" values at the edge vertices and midpont
     3409            // Get the "left" values at the edge vertices and midpoint
    33963410            // from the triangle k, edge i
    33973411            get_edge_data(&E_left, D, k, i);
     
    34253439                E_right.z1  = E_left.z;
    34263440               
    3427 
    34283441                E_right.w2  = E_right.w;
    34293442                E_right.h2  = E_right.h;
     
    34333446                E_right.z2  = E_left.z;
    34343447
    3435                
    3436 
    3437 
    3438 
    3439                
    3440 
    3441 
    34423448            } else {
    34433449                // Neighbour is a real triangle
    34443450                m = D->neighbour_edges[k3i];
     3451                n3m  = 3*n+m;
    34453452
    34463453                get_edge_data(&E_right, D, n, m);
     
    35003507            // Update triangle k with flux from edge i
    35013508            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];
    35043511
    35053512            D->already_computed_flux[k3i] = call; // #k Done
     
    35103517            if (n >= 0) {
    35113518
    3512                 n3m  = 3 * n + (m + 1) % 3;
    35133519                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
    35183524                D->already_computed_flux[n3m] = call; // #n Done
    35193525            }
     
    39763982      stage, xmomentum and ymomentum
    39773983
    3978       The
    3979 
    39803984      The maximal allowable speed computed by the flux_function for each volume
    39813985      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  
    9696domain.set_timestepping_method('rk2')
    9797#domain.set_CFL(0.7)
    98 #domain.set_beta(1.5)
     98domain.set_beta(1.5)
    9999domain.set_extrapolate_velocity(True)
    100100domain.set_compute_fluxes_method('wb_3')
     101domain.set_extrapolate_velocity(True)
    101102domain.set_name('merimbula')
    102103
Note: See TracChangeset for help on using the changeset viewer.