Changeset 5587


Ignore:
Timestamp:
Jul 30, 2008, 5:03:47 PM (16 years ago)
Author:
steve
Message:

Added in minimum height

Location:
anuga_work/development/anuga_1d
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/comp_flux_ext.c

    r5564 r5587  
    2121       
    2222
     23
     24
     25// Function to obtain speed from momentum and depth.
     26// This is used by flux functions
     27// Input parameters uh and h may be modified by this function.
     28double _compute_speed(double *uh,
     29                      double *h,
     30                      double epsilon,
     31                      double h0) {
     32 
     33  double u;
     34 
     35  if (*h < epsilon) {
     36    *h = 0.0;  //Could have been negative
     37     u = 0.0;
     38  } else {
     39    u = *uh/(*h + h0/ *h);   
     40  }
     41 
     42
     43  // Adjust momentum to be consistent with speed
     44  *uh = u * *h;
     45 
     46  return u;
     47}
     48
     49
     50
    2351//Innermost flux function (using w=z+h)
    2452int _flux_function(double *q_left, double *q_right,
    2553        double z_left, double z_right,
    26                 double normals, double g, double epsilon,
     54                   double normals, double g, double epsilon, double h0,
    2755                double *edgeflux, double *max_speed) {
    2856               
    … …  
    3361                double s_max, s_min, denom;
    3462               
    35                
     63                //printf("h0 = %f \n",h0);
    3664                ql[0] = q_left[0];
    3765                ql[1] = q_left[1];
    … …  
    5482                h_left = w_left-z;
    5583                uh_left = ql[1];
    56                 if (h_left < epsilon) {
    57                         u_left = 0.0; h_left = 0.0;
    58                 }
    59                 else {
    60                         u_left = uh_left/h_left;
    61                 }
    62        
     84
     85                u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
     86
    6387                w_right = qr[0];
    6488                h_right = w_right-z;
    6589                uh_right = qr[1];
    66                 if (h_right < epsilon) {
    67                         u_right = 0.0; h_right = 0.0;
    68                 }
    69                 else {
    70                         u_right = uh_right/h_right;
    71                 }
     90
     91                u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
    7292 
    7393                soundspeed_left = sqrt(g*h_left);
    … …  
    113133// Computational function for flux computation
    114134double _compute_fluxes_ext(double timestep,
    115                 double epsilon,
    116                 double g,
    117                 long* neighbours,
    118                 long* neighbour_vertices,
    119                 double* normals,
    120                 double* areas,
    121                 double* stage_edge_values,
    122                 double* xmom_edge_values,
    123                 double* bed_edge_values,
    124                 double* stage_boundary_values,
    125                 double* xmom_boundary_values,
    126                 double* stage_explicit_update,
    127                 double* xmom_explicit_update,
    128                 int number_of_elements,
    129                 double* max_speed_array) {
    130                
    131                 double flux[2], ql[2], qr[2], edgeflux[2];
     135                           double epsilon,
     136                           double g,
     137                           double h0,
     138                           long* neighbours,
     139                           long* neighbour_vertices,
     140                           double* normals,
     141                           double* areas,
     142                           double* stage_edge_values,
     143                           double* xmom_edge_values,
     144                           double* bed_edge_values,
     145                           double* stage_boundary_values,
     146                           double* xmom_boundary_values,
     147                           double* stage_explicit_update,
     148                           double* xmom_explicit_update,
     149                           int number_of_elements,
     150                           double* max_speed_array) {
     151               
     152                double flux[2], ql[2], qr[2], edgeflux[2];
    132153                double zl, zr, max_speed, normal;
    133154                int k, i, ki, n, m, nm=0;
    … …  
    159180                               
    160181                                normal = normals[ki];
    161                                 _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
     182                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
    162183                                flux[0] -= edgeflux[0];
    163184                                flux[1] -= edgeflux[1];
    … …  
    209230        *max_speed_array;
    210231   
    211   double timestep, epsilon, g;
     232  double timestep, epsilon, g, h0;
    212233  int number_of_elements;
    213234   
    214235  // Convert Python arguments to C
    215   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
     236  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO",
    216237                        &timestep,
    217238                        &epsilon,
    218239                        &g,
     240                        &h0,
    219241                        &neighbours,
    220242                        &neighbour_vertices,
    … …  
    240262                                 epsilon,
    241263                                 g,
     264                                 h0,
    242265                                 (long*) neighbours -> data,
    243266                                 (long*) neighbour_vertices -> data,
    … …  
    281304        *max_speed_array;
    282305   
    283   double timestep, epsilon, g;
     306  double timestep, epsilon, g, h0;
    284307  int number_of_elements;
    285308
    … …  
    299322    epsilon           = get_python_double(domain,"epsilon");
    300323    g                 = get_python_double(domain,"g");
     324    h0                = get_python_double(domain,"h0");
    301325 
    302326   
    … …  
    327351    // the explicit update arrays
    328352    timestep = _compute_fluxes_ext(timestep,
    329                                  epsilon,
    330                                  g,
    331                                  (long*) neighbours -> data,
    332                                  (long*) neighbour_vertices -> data,
    333                                  (double*) normals -> data,
    334                                  (double*) areas -> data,
    335                                  (double*) stage_vertex_values -> data,
    336                                  (double*) xmom_vertex_values -> data,
    337                                  (double*) bed_vertex_values -> data,
    338                                  (double*) stage_boundary_values -> data,
    339                                  (double*) xmom_boundary_values -> data,
    340                                  (double*) stage_explicit_update -> data,
    341                                  (double*) xmom_explicit_update -> data,
    342                                  number_of_elements,
    343                                  (double*) max_speed_array -> data);
     353                                   epsilon,
     354                                   g,
     355                                   h0,
     356                                   (long*) neighbours -> data,
     357                                   (long*) neighbour_vertices -> data,
     358                                   (double*) normals -> data,
     359                                   (double*) areas -> data,
     360                                   (double*) stage_vertex_values -> data,
     361                                   (double*) xmom_vertex_values -> data,
     362                                   (double*) bed_vertex_values -> data,
     363                                   (double*) stage_boundary_values -> data,
     364                                   (double*) xmom_boundary_values -> data,
     365                                   (double*) stage_explicit_update -> data,
     366                                   (double*) xmom_explicit_update -> data,
     367                                   number_of_elements,
     368                                   (double*) max_speed_array -> data);
    344369
    345370

  • anuga_work/development/anuga_1d/config.py

    r5535 r5587  
    1 """Module where global model parameters are set
     1"""Module where global model parameters are set for anuga_1d
    22"""
    33
    44epsilon = 1.0e-12
     5h0 = 1.0e-6
    56
    67default_boundary_tag = 'exterior'
    … …  
    7576
    7677
    77 beta_w = 0.9
     78beta_w = 1.5
    7879beta_h = 0.2
    7980CFL = 1.0  #FIXME (ole): Is this in use yet??

  • anuga_work/development/anuga_1d/dam_h_sudi.py

    r5535 r5587  
    11import os
    22from math import sqrt
    3 from shallow_water_h import *
     3#from shallow_water_h import *
     4from shallow_water_domain import *
    45from Numeric import zeros, Float
    56from analytic_dam_sudi import AnalyticDam
    … …  
    3132
    3233L=2000.0
    33 N=100
     34N=400
    3435
    3536cell_len=L/N
    … …  
    4142domain=Domain(points)
    4243
    43 domain.default_order = 1
     44domain.default_order = 2
    4445domain.default_time_order = 1
    45 #domain.cfl = 1.0
    46 #domain.limiter = "minmod"
     46domain.cfl = 1.0
     47domain.limiter = "vanleer"
    4748
    4849
    4950
    50 
    51 
    52    
    53    
    5451def height(x):
    5552    y=zeros(len(x), Float)
    … …  
    6360    return y
    6461
    65 domain.set_quantity('height', height)
     62domain.set_quantity('stage',height) #('height', height)
    6663domain.order=domain.default_order
    6764print "domain order", domain.order
    … …  
    7875t0=time.time()
    7976yieldstep=30.0
    80 finaltime=30.0
    81 print "integral", domain.quantities['height'].get_integral()
     77finaltime=20.0
     78print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
    8279for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    8380    domain.write_time()
    84     print "integral", domain.quantities['height'].get_integral()
     81    print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
    8582    if t>0.0:
    86         HeightQ=domain.quantities['height'].vertex_values
     83        HeightQ=domain.quantities['stage'].vertex_values #['height'].vertex_values
    8784        MomentumQ=domain.quantities['xmomentum'].vertex_values
    8885        h, uh=analytical_sol(X.flat, domain.time)

  • anuga_work/development/anuga_1d/dry_dam_sudi.py

    r5565 r5587  
    2222        # Calculate Analytical Solution at time t > 0
    2323        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t)
    24         h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
     24        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
     25        u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1))
     26        h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1))
    2527
    26         if ( x[i] <= -t*sqrt(g*h1) ):
     28        if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)):
     29            u[i] = 0.0
     30            h[i] = h0
     31        elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)):
     32            u[i] = u3_
     33            h[i] = h3_
     34
     35        elif ( x[i] <= -t*sqrt(g*h1) ):
    2736            u[i] = 0.0
    2837            h[i] = h1
    … …  
    5766
    5867L = 2000.0     # Length of channel (m)
    59 N = 100        # Number of compuational cells
     68N = 800        # Number of computational cells
    6069cell_len = L/N # Origin = 0.0
    6170
    … …  
    8594yieldstep = 1.0
    8695L = 2000.0     # Length of channel (m)
    87 number_of_cells = [100]#,200,500,1000,2000,5000,10000,20000]
     96number_of_cells = [200]#,200,500,1000,2000,5000,10000,20000]
    8897h_error = zeros(len(number_of_cells),Float)
    8998uh_error = zeros(len(number_of_cells),Float)
    … …  
    102111    domain.set_boundary({'exterior': Reflective_boundary(domain)})
    103112    domain.default_order = 2
    104     domain.default_time_order = 1
     113    domain.default_time_order = 2
    105114    domain.cfl = 1.0
    106     domain.limiter = "vanleer"
     115    domain.limiter = "minmod"
    107116
    108117    t0 = time.time()
    109118
    110119    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    111         pass
     120        domain.write_time()
    112121
    113122    N = float(N)

  • anuga_work/development/anuga_1d/run_profile_dry.py

    r5565 r5587  
    114114
    115115S = pstats.Stats(FN)
    116 s = S.sort_stats('cumulative').print_stats(20)
     116s = S.sort_stats('cumulative').print_stats(30)
    117117
    118118print s   

  • anuga_work/development/anuga_1d/shallow_water_domain.py

    r5563 r5587  
    5959                                tagged_elements)
    6060       
    61         from config import minimum_allowed_height, g
     61        from config import minimum_allowed_height, g, h0
    6262        self.minimum_allowed_height = minimum_allowed_height
    6363        self.g = g
     64        self.h0 = h0
    6465
    6566        #forcing terms not included in 1d domain ?WHy?
    … …  
    308309    """
    309310
    310     from config import g, epsilon
     311    from config import g, epsilon, h0
    311312    from math import sqrt
    312313    from Numeric import array
    … …  
    333334        h_left = 0.0
    334335    else:
    335         u_left  = uh_left/h_left
     336        u_left  = uh_left/(h_left +  h0/h_left)
     337
     338
     339    uh_left = u_left*h_left
    336340
    337341
    … …  
    339343    h_right  = w_right-z
    340344    uh_right = q_right[1]
    341 
    342345
    343346    if h_right < epsilon:
    … …  
    345348        h_right = 0.0
    346349    else:
    347         u_right  = uh_right/h_right
     350        u_right  = uh_right/(h_right + h0/h_right)
     351
     352    uh_right = u_right*h_right
    348353
    349354    #vh_left  = q_left[2]

  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5563 r5587  
    5050        zr = 0.0
    5151
     52        #This assumes h0 = 1.0e-3!!
    5253        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
    53 
    54         assert allclose(array([2.0, 8.9],Float), edgeflux)
    55         assert allclose(5.1304951685, maxspeed)
     54        #print maxspeed
     55        #print edgeflux
     56       
     57        assert allclose(array([1.998002, 8.89201198],Float), edgeflux)
     58        assert allclose(5.1284971665, maxspeed)
    5659
    5760        normal = -1.0
    … …  
    6366        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
    6467
    65         assert allclose(array([-2.0, -8.9],Float), edgeflux)
    66         assert allclose(5.1304951685, maxspeed)
     68
     69        #print maxspeed
     70        #print edgeflux       
     71       
     72        assert allclose(array([-1.998002, -8.89201198],Float), edgeflux)
     73        assert allclose(5.1284971665, maxspeed)
    6774
    6875    def test_domain_flux_function(self):
    … …  
    231238    """
    232239
    233     from config import g, epsilon
     240    from config import g, epsilon, h0
    234241    from math import sqrt
    235242    from Numeric import array
    … …  
    256263        h_left = 0.0
    257264    else:
    258         u_left  = uh_left/h_left
    259 
     265        u_left  = uh_left/(h_left +  h0/h_left)
     266
     267
     268    uh_left = u_left*h_left
    260269
    261270    w_right  = q_right[0]  #w=h+z
    … …  
    268277        h_right = 0.0
    269278    else:
    270         u_right  = uh_right/h_right
    271 
     279        u_right  = uh_right/(h_right + h0/h_right)
     280
     281    uh_right = u_right*h_right
     282   
    272283    #vh_left  = q_left[2]
    273284    #vh_right = q_right[2]
Note: See TracChangeset for help on using the changeset viewer.