Changeset 5587 for anuga_work/development/anuga_1d/comp_flux_ext.c

Timestamp:

Jul 30, 2008, 5:03:47 PM (16 years ago)

Author:

steve

Message:

Added in minimum height

File:

• anuga_work/development/anuga_1d/comp_flux_ext.c (modified) (10 diffs)

anuga_work/development/anuga_1d/comp_flux_ext.c

@@ -21 +21 @@
         
 
+
+
+// Function to obtain speed from momentum and depth.
+// This is used by flux functions
+// Input parameters uh and h may be modified by this function.
+double _compute_speed(double *uh, 
+                      double *h, 
+                      double epsilon, 
+                      double h0) {
+  
+  double u;
+  
+  if (*h < epsilon) {
+    *h = 0.0;  //Could have been negative
+     u = 0.0;
+  } else {
+    u = *uh/(*h + h0/ *h);    
+  }
+  
+
+  // Adjust momentum to be consistent with speed
+  *uh = u * *h;
+  
+  return u;
+}
+
+
+
 //Innermost flux function (using w=z+h)
 int _flux_function(double *q_left, double *q_right,
         double z_left, double z_right,
-                double normals, double g, double epsilon,
+                   double normals, double g, double epsilon, double h0,
                 double *edgeflux, double *max_speed) {
                 
@@ -33 +61 @@
                 double s_max, s_min, denom;
                 
-                
+                //printf("h0 = %f \n",h0);
                 ql[0] = q_left[0];
                 ql[1] = q_left[1];
@@ -54 +82 @@
                 h_left = w_left-z;
                 uh_left = ql[1];
-                if (h_left < epsilon) {
-                        u_left = 0.0; h_left = 0.0;
-                }
-                else {
-                        u_left = uh_left/h_left;
-                }
-        
+
+                u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
+
                 w_right = qr[0];
                 h_right = w_right-z;
                 uh_right = qr[1];
-                if (h_right < epsilon) {
-                        u_right = 0.0; h_right = 0.0; 
-                }
-                else {
-                        u_right = uh_right/h_right;
-                }
+
+                u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
   
                 soundspeed_left = sqrt(g*h_left);
@@ -113 +133 @@
 // Computational function for flux computation
 double _compute_fluxes_ext(double timestep,
-                double epsilon,
-                double g,
-                long* neighbours,
-                long* neighbour_vertices,
-                double* normals,
-                double* areas,
-                double* stage_edge_values,
-                double* xmom_edge_values,
-                double* bed_edge_values,
-                double* stage_boundary_values,
-                double* xmom_boundary_values,
-                double* stage_explicit_update,
-                double* xmom_explicit_update,
-                int number_of_elements,
-                double* max_speed_array) {
-                
-                double flux[2], ql[2], qr[2], edgeflux[2];
+                           double epsilon,
+                           double g,
+                           double h0,
+                           long* neighbours,
+                           long* neighbour_vertices,
+                           double* normals,
+                           double* areas,
+                           double* stage_edge_values,
+                           double* xmom_edge_values,
+                           double* bed_edge_values,
+                           double* stage_boundary_values,
+                           double* xmom_boundary_values,
+                           double* stage_explicit_update,
+                           double* xmom_explicit_update,
+                           int number_of_elements,
+                           double* max_speed_array) {
+                
+                double flux[2], ql[2], qr[2], edgeflux[2];
                 double zl, zr, max_speed, normal;
                 int k, i, ki, n, m, nm=0;
@@ -159 +180 @@
                                 
                                 normal = normals[ki];
-                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
+                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
                                 flux[0] -= edgeflux[0];
                                 flux[1] -= edgeflux[1];
@@ -209 +230 @@
         *max_speed_array;
     
-  double timestep, epsilon, g;
+  double timestep, epsilon, g, h0;
   int number_of_elements;
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO",
                         &timestep,
                         &epsilon,
                         &g,
+                        &h0,
                         &neighbours,
                         &neighbour_vertices,
@@ -240 +262 @@
                                  epsilon,
                                  g,
+                                 h0,
                                  (long*) neighbours -> data,
                                  (long*) neighbour_vertices -> data,
@@ -281 +304 @@
         *max_speed_array;
     
-  double timestep, epsilon, g;
+  double timestep, epsilon, g, h0;
   int number_of_elements;
 
@@ -299 +322 @@
     epsilon           = get_python_double(domain,"epsilon");
     g                 = get_python_double(domain,"g");
+    h0                = get_python_double(domain,"h0");
  
     
@@ -327 +351 @@
     // the explicit update arrays 
     timestep = _compute_fluxes_ext(timestep,
-                                 epsilon,
-                                 g,
-                                 (long*) neighbours -> data,
-                                 (long*) neighbour_vertices -> data,
-                                 (double*) normals -> data,
-                                 (double*) areas -> data,
-                                 (double*) stage_vertex_values -> data,
-                                 (double*) xmom_vertex_values -> data,
-                                 (double*) bed_vertex_values -> data,
-                                 (double*) stage_boundary_values -> data,
-                                 (double*) xmom_boundary_values -> data,
-                                 (double*) stage_explicit_update -> data,
-                                 (double*) xmom_explicit_update -> data,
-                                 number_of_elements,
-                                 (double*) max_speed_array -> data);
+                                   epsilon,
+                                   g,
+                                   h0,
+                                   (long*) neighbours -> data,
+                                   (long*) neighbour_vertices -> data,
+                                   (double*) normals -> data,
+                                   (double*) areas -> data,
+                                   (double*) stage_vertex_values -> data,
+                                   (double*) xmom_vertex_values -> data,
+                                   (double*) bed_vertex_values -> data,
+                                   (double*) stage_boundary_values -> data,
+                                   (double*) xmom_boundary_values -> data,
+                                   (double*) stage_explicit_update -> data,
+                                   (double*) xmom_explicit_update -> data,
+                                   number_of_elements,
+                                   (double*) max_speed_array -> data);