Changeset 5587


Ignore:
Timestamp:
Jul 30, 2008, 5:03:47 PM (11 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.