Changeset 5827 for anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

Timestamp:

Oct 9, 2008, 12:54:08 PM (15 years ago)

Author:

sudi

Message:

Adding new versions of shallow_water_domain

File:

• anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c (modified) (8 diffs)

anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

@@ -26 +26 @@
 //Innermost flux function (using w=z+h)
 int _flux_function_wellbalanced(double *q_left, double *q_right,
-        double z_left, double z_right,
                 double normals, double g, double epsilon,
                 double *edgeflux, double *max_speed) {
                 
                 int i;
+        double z_left, z_right;
                 double ql[2], qr[2], flux_left[2], flux_right[2];
                 double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
@@ -36 +36 @@
                 double s_max, s_min, denom;
                 
+                w_left = q_left[0];
+                uh_left = q_left[1];
+                uh_left = uh_left*normals;
+        z_left = q_left[2];
+        
+                w_right = q_right[0];
+                uh_right = q_right[1];
+                uh_right = uh_right*normals;
+        z_right = q_right[2];
+          
                 if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
-                ql[0] = q_left[0];
-                ql[1] = q_left[1];
-                ql[1] = ql[1]*normals;
-        
-                qr[0] = q_right[0];
-                qr[1] = q_right[1];
-                qr[1] = qr[1]*normals;
-          
+        
                 //z = (z_left+z_right)/2.0;
                 z_star = max(z_left, z_right);                                                                                                          //equation(7)  
                                    
                 // Compute speeds in x-direction
-                w_left = ql[0];
                 h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1].
                 h_left_star = max(0, w_left-z_star);                                                                                                            //(8)
-                uh_left = ql[1];
+    
                 if (h_left_star < epsilon) {
                         u_left = 0.0; h_left_star = 0.0;
@@ -60 +62 @@
                 }
         
-                w_right = qr[0];
                 h_right = w_right-z_right;                                                                                                                                                                                                     
                 h_right_star = max(0, w_right-z_star);                                                                                                          //(9)
-                uh_right = qr[1];
+
                 if (h_right_star < epsilon) {
                         u_right = 0.0; h_right_star = 0.0; 
@@ -134 +135 @@
                 double* max_speed_array) {
                 
-                double flux[2], ql[2], qr[2], edgeflux[2];
+                double flux[2], ql[3], qr[3], edgeflux[2];
                 double zl, zr, max_speed, normal;
                 int k, i, ki, n, m, nm=0;
@@ -147 +148 @@
                                 ql[0] = stage_edge_values[ki];
                                 ql[1] = xmom_edge_values[ki];
-                                zl = bed_edge_values[ki];
-                                
+                                ql[2] = bed_edge_values[ki];
+                                //ql[3] = velocity_edge_values[ki];
+                //ql[4] = height_edge_values[ki];
+                
                                 n = neighbours[ki];
                                 if (n<0) {
@@ -154 +157 @@
                                         qr[0] = stage_boundary_values[m];
                                         qr[1] = xmom_boundary_values[m];
-                                        zr = zl;
+                                        qr[2] = ql[2];
+                    
                                 } else {
                                         m = neighbour_vertices[ki];
@@ -160 +164 @@
                                         qr[0] = stage_edge_values[nm];
                                         qr[1] = xmom_edge_values[nm];
-                                        zr = bed_edge_values[nm];                               
+                                        qr[2] = bed_edge_values[nm];                            
                                 }
                                 
                                 normal = normals[ki];
-                                _flux_function_wellbalanced(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
+                                _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
                                 flux[0] -= edgeflux[0];
                                 flux[1] -= edgeflux[1];
@@ -186 +190 @@
                         max_speed_array[k]=max_speed;
                 }
-                printf("timestep = %f \n",timestep);
+                //printf("timestep = %f \n",timestep);
                 return timestep;        
         }