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

Added in minimum height

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.