Changeset 8184
- Timestamp:
- Jun 3, 2011, 12:51:05 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c
r8182 r8184 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, double);17 int _flux_function_pipe _naive(double *, double *, double *, double *, double, double, double, double *, double *);16 //void compute_flux(double [], double [], double *, double, double, double); 17 int _flux_function_pipe(double *, double *, double *, double *, double, double, double, double *, double *); 18 18 19 19 // Shared code snippets … … 21 21 22 22 23 // Hard coded flux computation choice 24 #define _flux_function_pipe _flux_function_pipe_naive 23 // Hard coded flux computation choices 24 // Temporary hack 25 #define flux_formula0 flux_formula_press0 26 #define flux_formula1 flux_formula_press1 27 #define compute_flux compute_flux_godunov 25 28 26 29 27 30 // ********************************************************************* 28 31 // Formulae for fluxes 29 double flux_formula0(double d, double u, double h, double b, double t, double g) {30 return flux_formula_press0(d, u, h, b, t, g);31 }32 33 double flux_formula1(double d, double u, double h, double b, double t, double g) {34 return flux_formula_press1(d, u, h, b, t, g);35 }36 32 37 33 // Formulae for pressurised fluxes … … 57 53 // ********************************************************************* 58 54 // Flux function 59 void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double smax, double smin, double epsilon) { 55 56 // Naive version 57 void compute_flux_naive(double flux_m[], double flux_p[], double a_m, double a_p, double d_m, double d_p, double *edgeflux, double smax, double smin, double epsilon) { 58 60 59 int i; 60 61 61 if (smax-smin < epsilon) { 62 62 for (i=0; i<2; i++) edgeflux[i] = 0.0; … … 66 66 } 67 67 } 68 69 // Godunov method 70 void compute_flux_godunov(double flux_m[], double flux_p[], double a_m, double a_p, double d_m, double d_p, double *edgeflux, double smax, double smin, double epsilon) { 71 72 int i; 73 double denom; 74 75 denom = smax - smin; 76 if (denom < epsilon) { 77 for (i=0; i<2; i++) edgeflux[i] = 0.0; 78 } else { 79 edgeflux[0] = smax*flux_m[0] - smin*flux_p[0]; 80 edgeflux[0] += smax*smin * (a_p - a_m); 81 edgeflux[0] /= denom; 82 edgeflux[1] = smax*flux_m[1] - smin*flux_p[1]; 83 edgeflux[1] += smax*smin * (d_p - d_m); 84 edgeflux[1] /= denom; 85 } 86 } 87 88 // Godunov method with well balancing 89 void compute_flux_wb(double flux_m[], double flux_p[], double a_m, double a_p, double d_m, double d_p, double *edgeflux, double smax, double smin, double epsilon) { 90 // TODO 91 int i; 92 double denom; 93 94 denom = smax - smin; 95 if (denom < epsilon) { 96 for (i=0; i<2; i++) edgeflux[i] = 0.0; 97 } else { 98 edgeflux[0] = smax*flux_m[0] - smin*flux_p[0]; 99 edgeflux[0] += smax*smin * (a_p - a_m); 100 edgeflux[0] /= denom; 101 edgeflux[1] = smax*flux_m[1] - smin*flux_p[1]; 102 edgeflux[1] += smax*smin * (d_p - d_m); 103 edgeflux[1] /= denom; 104 } 105 } 106 68 107 // ********************************************************************* 69 108 … … 88 127 89 128 // ********************************************************************* 90 // NAIVE VERSION 91 int _flux_function_pipe_naive(double *q_leftm,double *q_leftp, double *q_rightm, 92 double *q_rightp, double g, double epsilon, double h0, double *edgeflux, double *max_speed){ 129 int _flux_function_pipe(double *q_leftm,double *q_leftp, double *q_rightm, double *q_rightp, double g, double epsilon, double h0, double *edgeflux, double *max_speed) { 93 130 94 131 double flux_m[2], flux_p[2]; … … 161 198 162 199 // Flux computation for left hand side 163 compute_flux(flux_m, flux_p, edgeflux, s_maxl, s_minl, epsilon);200 compute_flux(flux_m, flux_p, a_leftm, a_leftp, d_leftm, d_leftp, edgeflux, s_maxl, s_minl, epsilon); 164 201 fluxtemp0 = edgeflux[0]; 165 202 fluxtemp1 = edgeflux[1]; … … 173 210 174 211 // Flux computation for right hand side 175 compute_flux(flux_m, flux_p, edgeflux, s_maxr, s_minr, epsilon);212 compute_flux(flux_m, flux_p, a_rightm, a_rightp, d_rightm, d_rightp, edgeflux, s_maxr, s_minr, epsilon); 176 213 edgeflux[0] = edgeflux[0]-fluxtemp0; 177 214 edgeflux[1] = edgeflux[1]-fluxtemp1;
Note: See TracChangeset
for help on using the changeset viewer.