Changeset 8184


Ignore:
Timestamp:
Jun 3, 2011, 12:51:05 PM (13 years ago)
Author:
paul
Message:

Added Godunov method for flux computation of pipes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c

    r8182 r8184  
    1414double forcingterm_press(double, double, double, double, double, double, double, double, double);
    1515double 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);
     17int _flux_function_pipe(double *, double *, double *, double *, double, double, double, double *, double *);
    1818
    1919// Shared code snippets
     
    2121
    2222
    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
    2528
    2629
    2730// *********************************************************************
    2831// 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 }
    3632
    3733// Formulae for pressurised fluxes
     
    5753// *********************************************************************
    5854// Flux function
    59 void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double smax, double smin, double epsilon) {
     55
     56// Naive version
     57void 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
    6059  int i;
     60
    6161  if (smax-smin < epsilon) {
    6262    for (i=0; i<2; i++) edgeflux[i] = 0.0;
     
    6666  }
    6767}
     68
     69// Godunov method
     70void 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
     89void 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
    68107// *********************************************************************
    69108
     
    88127
    89128// *********************************************************************
    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){
     129int _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) {
    93130
    94131  double flux_m[2], flux_p[2];
     
    161198
    162199  // 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);
    164201  fluxtemp0 = edgeflux[0];
    165202  fluxtemp1 = edgeflux[1];
     
    173210
    174211  // 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);
    176213  edgeflux[0] = edgeflux[0]-fluxtemp0;
    177214  edgeflux[1] = edgeflux[1]-fluxtemp1;
Note: See TracChangeset for help on using the changeset viewer.