Changeset 8182


Ignore:
Timestamp:
May 27, 2011, 12:14:31 PM (13 years ago)
Author:
paul
Message:

Changed pipe flux function and forcing terms to use dual model formulae

File:
1 edited

Legend:

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

    r8177 r8182  
    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);
     16void compute_flux(double [], double [], double *, double, double, double);
    1717int _flux_function_pipe_naive(double *, double *, double *, double *, double, double, double, double *, double *);
    1818
     
    4141
    4242double 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;
    4444}
    4545
     
    5757// *********************************************************************
    5858// Flux function
    59 void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double denom, double epsilon) {
     59void compute_flux(double flux_m[], double flux_p[], double *edgeflux, double smax, double smin, double epsilon) {
    6060  int i;
    61   if (denom < epsilon) {
     61  if (smax-smin < epsilon) {
    6262    for (i=0; i<2; i++) edgeflux[i] = 0.0;
    6363  } else {
     
    7777// Pressurised forcing terms
    7878double 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);
    8080}
    8181
     
    161161
    162162  // 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);
    164164  fluxtemp0 = edgeflux[0];
    165165  fluxtemp1 = edgeflux[1];
     
    173173
    174174  // Flux computation for right hand side
    175   compute_flux(flux_m, flux_p, edgeflux, s_maxl-s_minl, epsilon);
     175  compute_flux(flux_m, flux_p, edgeflux, s_maxr, s_minr, epsilon);
    176176  edgeflux[0] = edgeflux[0]-fluxtemp0;
    177177  edgeflux[1] = edgeflux[1]-fluxtemp1;
Note: See TracChangeset for help on using the changeset viewer.