Changeset 8182
- Timestamp:
- May 27, 2011, 12:14:31 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c
r8177 r8182 14 14 double forcingterm_press(double, double, double, double, double, double, double, double, double); 15 15 double forcingterm_fs(double, double, double, double, double, double, double, double, double); 16 void compute_flux(double [], double [], double *, double, double );16 void compute_flux(double [], double [], double *, double, double, double); 17 17 int _flux_function_pipe_naive(double *, double *, double *, double *, double, double, double, double *, double *); 18 18 … … 41 41 42 42 double flux_formula_press1(double d, double u, double h, double b, double t, double g) { 43 return u*d + h*b;43 return u*d + (h-t)*b + 0.5*g*t*t*b; 44 44 } 45 45 … … 57 57 // ********************************************************************* 58 58 // Flux function 59 void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double denom, double epsilon) {59 void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double smax, double smin, double epsilon) { 60 60 int i; 61 if ( denom< epsilon) {61 if (smax-smin < epsilon) { 62 62 for (i=0; i<2; i++) edgeflux[i] = 0.0; 63 63 } else { … … 77 77 // Pressurised forcing terms 78 78 double forcingterm_press(double h_m, double h_p, double b_m, double b_p, double z_m, double z_p, double t_m, double t_p, double g) { 79 return g * ( 1/(0.5 * (t_m + t_p))) * ((h_m - h_p) * 0.5 * (b_m + b_p) + (b_m - b_p) * 0.5 * (h_m + h_p));79 return g * (0.5*(h_m+h_p) - 0.5*(t_m+t_p)) * (b_m-b_p) + g * (0.5*(h_m+h_p) - 0.5*(t_m+t_p)) * 0.5*(b_m+b_p) * (t_m-t_p) * 1/(0.5*(t_m+t_p)) + forcingterm_fs(t_m, t_p, b_m, b_p, z_m, z_p, t_m, t_p, g); 80 80 } 81 81 … … 161 161 162 162 // Flux computation for left hand side 163 compute_flux(flux_m, flux_p, edgeflux, s_maxl -s_minl, epsilon);163 compute_flux(flux_m, flux_p, edgeflux, s_maxl, s_minl, epsilon); 164 164 fluxtemp0 = edgeflux[0]; 165 165 fluxtemp1 = edgeflux[1]; … … 173 173 174 174 // Flux computation for right hand side 175 compute_flux(flux_m, flux_p, edgeflux, s_max l-s_minl, epsilon);175 compute_flux(flux_m, flux_p, edgeflux, s_maxr, s_minr, epsilon); 176 176 edgeflux[0] = edgeflux[0]-fluxtemp0; 177 177 edgeflux[1] = edgeflux[1]-fluxtemp1;
Note: See TracChangeset
for help on using the changeset viewer.