Ignore:
Timestamp:
Jun 19, 2011, 3:43:16 PM (14 years ago)
Author:
paul
Message:

Rewrote structure of pipe flux computation.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r8184 r8186  
    55const double pi = 3.14159265358979;
    66
    7 double flux_formula0(double, double, double, double, double, double);
    8 double flux_formula1(double, double, double, double, double, double);
    9 double flux_formula_press0(double, double, double, double, double, double);
    10 double flux_formula_press1(double, double, double, double, double, double);
    11 double flux_formula_fs0(double, double, double, double, double, double);
    12 double flux_formula_fs1(double, double, double, double, double, double);
    13 double forcingterm(double, double, double, double, double, double, double, double, double);
    14 double forcingterm_press(double, double, double, double, double, double, double, double, double);
    15 double forcingterm_fs(double, double, double, double, double, double, double, double, double);
    16 //void compute_flux(double [], double [], double *, double, double, double);
     7double flux_formula_press0(double*, double);
     8double flux_formula_press1(double*, double);
     9double flux_formula_fs0(double*, double);
     10double flux_formula_fs1(double*, double);
     11void compute_flux_naive(double[], double[], double, double, double, double, double[]);
     12void compute_flux_godunov(double[], double[], double, double, double, double, double[]);
     13void compute_flux_wb(double[], double[], double, double, double, double, double[]);
     14double forcingterm_press_naive(double[], double[], double);
     15double forcingterm_fs_naive(double[], double[], double);
    1716int _flux_function_pipe(double *, double *, double *, double *, double, double, double, double *, double *);
    1817
    … …  
    2019#include "util_ext.h"
    2120
    22 
    23 // Hard coded flux computation choices
    24 // Temporary hack
     21// Hard coded choices of computational methods
    2522#define flux_formula0 flux_formula_press0
    2623#define flux_formula1 flux_formula_press1
    2724#define compute_flux compute_flux_godunov
     25#define forcingterm forcingterm_press_naive
     26
     27
     28// Some formulae
     29#define SOUNDSPEED(a, b) sqrt((a) * (b))
     30#define POSMAX(a, b) max(max(a, b), 0)
     31#define NEGMIN(a, b) min(min(a, b), 0)
    2832
    2933
    … …  
    3236
    3337// Formulae for pressurised fluxes
    34 double flux_formula_press0(double d, double u, double h, double b, double t, double g) {
    35   return d;
    36 }
    37 
    38 double flux_formula_press1(double d, double u, double h, double b, double t, double g) {
    39   return u*d + (h-t)*b + 0.5*g*t*t*b;
     38double flux_formula_press0(double *q, double g) {
     39  return q[1];
     40}
     41
     42double flux_formula_press1(double *q, double g) {
     43  return q[4]*q[1] + (q[3]-q[6])*q[5] + 0.5*g*q[6]*q[6]*q[5];
    4044}
    4145
    4246// Formulae for free surface
    43 double flux_formula_fs0(double d, double u, double h, double b, double t, double g) {
    44   return d;
    45 }
    46 
    47 double flux_formula_fs1(double d, double u, double h, double b, double t, double g) {
    48   return u*d + 0.5*g*h*h*b;
    49 }
    50 // *********************************************************************
    51 
    52 
    53 // *********************************************************************
    54 // Flux function
     47double flux_formula_fs0(double *q, double g) {
     48  return q[1];
     49}
     50
     51double flux_formula_fs1(double *q, double g) {
     52  return q[4]*q[1] + 0.5*g*q[3]*q[3]*q[5];
     53}
     54// *********************************************************************
     55
     56
     57// *********************************************************************
     58// Flux computation method
    5559
    5660// 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) {
     61void compute_flux_naive(double *q_m, double *q_p, double g, double smax, double smin, double epsilon, double *edgeflux){
    5862
    5963  int i;
    … …  
    6266    for (i=0; i<2; i++) edgeflux[i] = 0.0;
    6367  } else {
    64     edgeflux[0] = 0.5*(flux_m[0]+flux_p[0]);
    65     edgeflux[1] = 0.5*(flux_m[1]+flux_p[1]);
     68    edgeflux[0] = 0.5*(flux_formula0(q_m, g) + flux_formula0(q_p, g));
     69    edgeflux[1] = 0.5*(flux_formula1(q_m, g) + flux_formula1(q_p, g));
    6670  }
    6771}
    6872
    6973// 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) {
     74void compute_flux_godunov(double *q_m, double *q_p, double g, double smax, double smin, double epsilon, double *edgeflux) {
    7175
    7276  int i;
    7377  double denom;
     78  double a_p, a_m;
    7479
    7580  denom = smax - smin;
    … …  
    7782    for (i=0; i<2; i++) edgeflux[i] = 0.0;
    7883  } else {
    79     edgeflux[0] = smax*flux_m[0] - smin*flux_p[0];
     84    edgeflux[0] = smax*flux_formula0(q_m, g) - smin*flux_formula0(q_p, g);
     85    a_p = (q_p[3]+q_p[2])*q_p[5];
     86    a_m = (q_m[3]+q_m[2])*q_m[5];
    8087    edgeflux[0] += smax*smin * (a_p - a_m);
    8188    edgeflux[0] /= denom;
    82     edgeflux[1] = smax*flux_m[1] - smin*flux_p[1];
    83     edgeflux[1] += smax*smin * (d_p - d_m);
     89
     90    edgeflux[1] = smax*flux_formula1(q_m, g) - smin*flux_formula1(q_p, g);
     91    edgeflux[1] += smax*smin * (q_p[4]*a_p - q_m[4]*a_m);
    8492    edgeflux[1] /= denom;
    8593    }
    8694}
    8795
    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 
    107 // *********************************************************************
    108 
    109 
    110 // *********************************************************************
    111 // Compute forcing terms
    112 double forcingterm(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) {
    113   return forcingterm_press(h_m, h_p, b_m, b_p, z_m, z_p, t_m, t_p, g);
    114 }
     96
     97// *********************************************************************
     98
     99
     100// *********************************************************************
     101// Forcing terms
    115102
    116103// Pressurised forcing terms
    117 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) {
    118   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);
     104// Naive version
     105double forcingterm_press_naive(double *q_m, double *q_p, double g) {
     106  double ret;
     107
     108  ret = g * (0.5*(q_m[3]+q_p[3]) - 0.5*(q_m[6]+q_p[6])) * (q_m[5]-q_p[5]) + g * (0.5*(q_m[3]+q_p[3]) - 0.5*(q_m[6]+q_p[6])) * 0.5*(q_m[5]+q_p[5]) * (q_m[6]-q_p[6]) * 1/(0.5*(q_m[6]+q_p[6]));
     109
     110  // Add in the free surface pressure term with h=t
     111  // Does this overwite q_p,q_m in the caller? If so, we don't want that!
     112  q_p[3] = q_p[6];
     113  q_m[3] = q_m[6];
     114  ret += forcingterm_fs_naive(q_m, q_p, g);
     115
     116  return ret;
    119117}
    120118
    121119// Free surface forcing terms
    122 double forcingterm_fs(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) {
    123   return -0.5*0.5 * g * (h_m + h_p) * (b_m + b_p) * (z_m - z_p) + 0.5 * g * (h_m + h_p) * (h_m + h_p) * 0.5*0.5 * (b_m - b_p);
    124 }
    125 // *********************************************************************
    126 
    127 
    128 // *********************************************************************
    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) {
    130 
    131   double flux_m[2], flux_p[2];
    132   double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, t_leftm, soundspeed_leftm;
    133   double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, t_leftp, soundspeed_leftp;
    134   double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, t_rightm, soundspeed_rightm;
    135   double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, t_rightp, soundspeed_rightp;
    136   double s_maxl, s_minl,s_maxr,s_minr;
    137   double fluxtemp1,fluxtemp0,speedtemp;
    138 
    139   a_leftm  = q_leftm[0];
    140   d_leftm = q_leftm[1];
    141   z_leftm  = q_leftm[2];
    142   h_leftm  = q_leftm[3];
    143   u_leftm  = q_leftm[4];
    144   b_leftm  = q_leftm[5];
    145   t_leftm  = q_leftm[6];
    146   w_leftm  = h_leftm+z_leftm;
    147 
    148   a_leftp  = q_leftp[0];
    149   d_leftp = q_leftp[1];
    150   z_leftp  = q_leftp[2];
    151   h_leftp  = q_leftp[3];
    152   u_leftp  = q_leftp[4];
    153   b_leftp  = q_leftp[5];
    154   t_leftp  = q_leftp[6];
    155   w_leftp  = h_leftp+z_leftp;
    156 
    157   a_rightm  = q_rightm[0];
    158   d_rightm = q_rightm[1];
    159   z_rightm  = q_rightm[2];
    160   h_rightm  = q_rightm[3];
    161   u_rightm  = q_rightm[4];
    162   b_rightm  = q_rightm[5];
    163   t_rightm  = q_rightm[6];
    164   w_rightm  = h_rightm+z_rightm;
    165 
    166   a_rightp  = q_rightp[0];
    167   d_rightp = q_rightp[1];
    168   z_rightp  = q_rightp[2];
    169   h_rightp  = q_rightp[3];
    170   u_rightp  = q_rightp[4];
    171   b_rightp  = q_rightp[5];
    172   t_rightp  = q_rightp[6];
    173   w_rightp  = h_rightp+z_rightp;
    174 
    175   soundspeed_leftp = sqrt(g*h_leftp);
    176   soundspeed_leftm = sqrt(g*h_leftm);
    177   soundspeed_rightp = sqrt(g*h_rightp);
    178   soundspeed_rightm = sqrt(g*h_rightm);
    179 
    180   s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
    181   if (s_maxl < 0.0) s_maxl = 0.0;
    182 
    183   s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
    184   if (s_minl > 0.0) s_minl = 0.0;
    185 
    186   s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
    187   if (s_maxr < 0.0) s_maxr = 0.0;
    188 
    189   s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
    190   if (s_minr > 0.0) s_minr = 0.0;
    191 
    192   // Flux formulas for left hand side
    193   flux_m[0] = flux_formula0(d_leftm, u_leftm, h_leftm, b_leftm, t_leftm, g);
    194   flux_m[1] = flux_formula1(d_leftm, u_leftm, h_leftm, b_leftm, t_leftm, g);
    195 
    196   flux_p[0] = flux_formula0(d_leftp, u_leftp, h_leftp, b_leftp, t_leftp, g);
    197   flux_p[1] = flux_formula1(d_leftp, u_leftp, h_leftp, b_leftp, t_leftp, g);
    198 
    199   // Flux computation for left hand side
    200   compute_flux(flux_m, flux_p, a_leftm, a_leftp, d_leftm, d_leftp, edgeflux, s_maxl, s_minl, epsilon);
     120// Naive version
     121double forcingterm_fs_naive(double *q_m, double *q_p, double g) {
     122  return -0.5*0.5 * g * (q_m[3] + q_p[3]) * (q_m[5] + q_p[5]) * (q_m[2] - q_p[2]) + 0.5 * g * (q_m[3] + q_p[3]) * (q_m[3] + q_p[3]) * 0.5*0.5 * (q_m[5] - q_p[5]);
     123}
     124// *********************************************************************
     125
     126
     127// *********************************************************************
     128/* Function to compute the flux to return to python
     129   q_... are arrays of the quantities such as bed, height, velocity etc.
     130   left refers to the left interface of the cell.
     131   right refers to the left interface of the cell.
     132   m refers to the values just to the left of the interface.
     133   p refers to the values just to the right of the interface.
     134
     135   The quantities for each index are
     136
     137   0 Area (a)
     138   1 Discharge (d)
     139   2 Elevation (z)
     140   3 Height (h)
     141   4 Velocity (u)
     142   5 Width (b)
     143   6 Top (t)
     144*/
     145int _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) {
     146
     147  double soundspeed_leftm, soundspeed_leftp, soundspeed_rightm, soundspeed_rightp;
     148  double s_maxl, s_minl, s_maxr, s_minr;
     149  double fluxtemp1, fluxtemp0, speedtemp;
     150
     151  // Compute sound speeds
     152  soundspeed_leftp = SOUNDSPEED(g, q_leftp[3]);
     153  soundspeed_leftm = SOUNDSPEED(g, q_leftm[3]);
     154  soundspeed_rightp = SOUNDSPEED(g, q_rightp[3]);
     155  soundspeed_rightm = SOUNDSPEED(g, q_rightm[3]);
     156
     157  // Compute max and min sound speeds
     158  s_maxl = POSMAX(q_leftm[4] + soundspeed_leftm, q_leftp[4] + soundspeed_leftp);
     159  s_minl = NEGMIN(q_leftm[4] - soundspeed_leftm, q_leftp[4] - soundspeed_leftp);
     160  s_maxr = POSMAX(q_rightm[4] + soundspeed_rightm, q_rightp[4] + soundspeed_rightp);
     161  s_minr = NEGMIN(q_rightm[4] - soundspeed_rightm, q_rightp[4] - soundspeed_rightp);
     162
     163  // Flux computation for left hand side of the current cell
     164  compute_flux(q_leftm, q_leftp, g, s_maxl, s_minl, epsilon, edgeflux);
    201165  fluxtemp0 = edgeflux[0];
    202166  fluxtemp1 = edgeflux[1];
    203167
    204   // Flux formulas for right hand side
    205   flux_m[0] = flux_formula0(d_rightm, u_rightm, h_rightm, b_rightm, t_rightm, g);
    206   flux_m[1] = flux_formula1(d_rightm, u_rightm, h_rightm, b_rightm, t_rightm, g);
    207 
    208   flux_p[0] = flux_formula0(d_rightp, u_rightp, h_rightp, b_rightp, t_rightm, g);
    209   flux_p[1] = flux_formula1(d_rightp, u_rightp, h_rightp, b_rightp, t_rightm, g);
    210 
    211   // Flux computation for right hand side
    212   compute_flux(flux_m, flux_p, a_rightm, a_rightp, d_rightm, d_rightp, edgeflux, s_maxr, s_minr, epsilon);
    213   edgeflux[0] = edgeflux[0]-fluxtemp0;
    214   edgeflux[1] = edgeflux[1]-fluxtemp1;
    215 
    216   // Forcing terms 
    217   edgeflux[1] -= forcingterm(h_rightm, h_leftp, b_rightm, b_leftp, z_rightm, z_leftp, t_rightm, t_leftp, g);
     168  // Flux computation for right hand side of the current cell
     169  compute_flux(q_rightm, q_rightp, g, s_maxr, s_minr, epsilon, edgeflux);
     170
     171  // Total flux is the difference of the left and right fluxes
     172  edgeflux[0] = edgeflux[0] - fluxtemp0;
     173  edgeflux[1] = edgeflux[1] - fluxtemp1;
     174
     175  // Add in Forcing terms 
     176  edgeflux[1] -= forcingterm(q_rightm, q_leftp, g);
    218177
    219178  // Maximal wavespeed
    … …  
    221180    *max_speed = 0.0;
    222181  }else{
    223     speedtemp = max(fabs(s_maxl),fabs(s_minl));
    224     speedtemp = max(speedtemp,fabs(s_maxr));
    225     speedtemp = max(speedtemp,fabs(s_minr));
     182    speedtemp = max(fabs(s_maxl), fabs(s_minl));
     183    speedtemp = max(speedtemp, fabs(s_maxr));
     184    speedtemp = max(speedtemp, fabs(s_minr));
    226185    *max_speed = speedtemp;
    227186  }
Note: See TracChangeset for help on using the changeset viewer.