- Timestamp:
- Oct 9, 2008, 12:54:08 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
r5725 r5827 26 26 //Innermost flux function (using w=z+h) 27 27 int _flux_function_wellbalanced(double *q_left, double *q_right, 28 double z_left, double z_right,29 28 double normals, double g, double epsilon, 30 29 double *edgeflux, double *max_speed) { 31 30 32 31 int i; 32 double z_left, z_right; 33 33 double ql[2], qr[2], flux_left[2], flux_right[2]; 34 34 double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left; … … 36 36 double s_max, s_min, denom; 37 37 38 w_left = q_left[0]; 39 uh_left = q_left[1]; 40 uh_left = uh_left*normals; 41 z_left = q_left[2]; 42 43 w_right = q_right[0]; 44 uh_right = q_right[1]; 45 uh_right = uh_right*normals; 46 z_right = q_right[2]; 47 38 48 if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right); 39 ql[0] = q_left[0]; 40 ql[1] = q_left[1]; 41 ql[1] = ql[1]*normals; 42 43 qr[0] = q_right[0]; 44 qr[1] = q_right[1]; 45 qr[1] = qr[1]*normals; 46 49 47 50 //z = (z_left+z_right)/2.0; 48 51 z_star = max(z_left, z_right); //equation(7) 49 52 50 53 // Compute speeds in x-direction 51 w_left = ql[0];52 54 h_left = w_left-z_left; //This is used for computing the edgeflux[1]. 53 55 h_left_star = max(0, w_left-z_star); //(8) 54 uh_left = ql[1]; 56 55 57 if (h_left_star < epsilon) { 56 58 u_left = 0.0; h_left_star = 0.0; … … 60 62 } 61 63 62 w_right = qr[0];63 64 h_right = w_right-z_right; 64 65 h_right_star = max(0, w_right-z_star); //(9) 65 uh_right = qr[1]; 66 66 67 if (h_right_star < epsilon) { 67 68 u_right = 0.0; h_right_star = 0.0; … … 134 135 double* max_speed_array) { 135 136 136 double flux[2], ql[ 2], qr[2], edgeflux[2];137 double flux[2], ql[3], qr[3], edgeflux[2]; 137 138 double zl, zr, max_speed, normal; 138 139 int k, i, ki, n, m, nm=0; … … 147 148 ql[0] = stage_edge_values[ki]; 148 149 ql[1] = xmom_edge_values[ki]; 149 zl = bed_edge_values[ki]; 150 150 ql[2] = bed_edge_values[ki]; 151 //ql[3] = velocity_edge_values[ki]; 152 //ql[4] = height_edge_values[ki]; 153 151 154 n = neighbours[ki]; 152 155 if (n<0) { … … 154 157 qr[0] = stage_boundary_values[m]; 155 158 qr[1] = xmom_boundary_values[m]; 156 zr = zl; 159 qr[2] = ql[2]; 160 157 161 } else { 158 162 m = neighbour_vertices[ki]; … … 160 164 qr[0] = stage_edge_values[nm]; 161 165 qr[1] = xmom_edge_values[nm]; 162 zr= bed_edge_values[nm];166 qr[2] = bed_edge_values[nm]; 163 167 } 164 168 165 169 normal = normals[ki]; 166 _flux_function_wellbalanced(ql, qr, zl, zr,normal, g, epsilon, edgeflux, &max_speed);170 _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed); 167 171 flux[0] -= edgeflux[0]; 168 172 flux[1] -= edgeflux[1]; … … 186 190 max_speed_array[k]=max_speed; 187 191 } 188 printf("timestep = %f \n",timestep);192 //printf("timestep = %f \n",timestep); 189 193 return timestep; 190 194 }
Note: See TracChangeset
for help on using the changeset viewer.