Changeset 6454

Timestamp:

Mar 4, 2009, 3:37:20 PM (15 years ago)

Author:

steve

Message:

Added Padarn's C code

Location:

anuga_work/development/anuga_1d

Files:

• channel_domain_ext.c (modified) (9 diffs)
• channel_domain_ext_backup.c (added)
• comp_flux_ext_steve.c (added)

anuga_work/development/anuga_1d/channel_domain_ext.c

@@ -46 +46 @@
 
 
-
+//WELL BALANCED VERSION
 //Innermost flux function (using w=z+h)
-int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
+int _flux_function_channel21(double *q_leftm,double *q_leftp, double *q_rightm,
                            double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed) {
 
@@ -60 +60 @@
     double zphalf,zmhalf,hleftstar,hrightstar;
     double fluxtemp1,fluxtemp0,speedtemp;
-    double batemp;
+    double batemp,bphalf,bmhalf;
 
     zmhalf = max(q_leftm[2],q_leftp[2]);
@@ -100 +100 @@
     hrightstar = q_rightm[3];   
 
+     bphalf = 0.5*(b_rightm+b_rightp);
+     bmhalf = 0.5*(b_leftm+b_leftp);
+     //bphalf = min(b_rightm,b_rightp);
+     //bmhalf = min(b_leftm,b_leftp);
+
+
     soundspeed_leftp = sqrt(g*h_leftp);
     soundspeed_leftm = sqrt(g*h_leftm);
@@ -119 +125 @@
     // Flux formulas for left hand side
     flux_left[0] = d_leftm;
-    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
+    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*bmhalf;
 
     flux_right[0] = d_leftp;
-    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
+    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*bmhalf;
 
     
@@ -132 +138 @@
         edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
 
-        batemp = (h_leftp)*b_leftp-(h_leftm)*b_leftm;
-                edgeflux[0] += s_maxl*s_minl*batemp;
+        batemp = (q_leftp[3]+q_leftp[2])*bmhalf-(q_leftm[3]+q_leftm[2])*bmhalf;
+        edgeflux[0] += s_maxl*s_minl*batemp;
         edgeflux[0] /= denom;
         edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
@@ -147 +153 @@
     // Flux formulas for right hand side
     flux_left[0] = d_rightm;
-    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm;
+    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*bphalf;
 
     flux_right[0] = d_rightp;
-    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp;
+    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*bphalf;
     
    
@@ -161 +167 @@
         edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
 
-        batemp = (h_rightp)*b_rightp-(h_rightm)*b_rightm;
-
-                edgeflux[0] += s_maxr*s_minr*batemp;
+        batemp = (q_rightp[3]+q_rightp[2])*bphalf-(q_rightm[3]+q_rightm[2])*bphalf;
+
+        edgeflux[0] += s_maxr*s_minr*batemp;
         edgeflux[0] /= denom;
         edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
@@ -177 +183 @@
     
    
-     edgeflux[1]-=0.5*g*h_rightm*h_rightm*b_rightm-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*b_leftp;
-
-   
-     //edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
+    edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
+
+    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
+
+    edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
+
+   
+    //edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
       
      //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
@@ -196 +206 @@
     return 0;
 }
-
-
-
+// GOOD BUT NOT WELL BALANCED VERSION
+int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
+                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed){
+
+    int i;
+    double flux_left[2], flux_right[2];
+    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
+     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
+ double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
+ double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
+    double s_maxl, s_minl,s_maxr,s_minr, denom;
+    double zphalf,zmhalf,hleftstar,hrightstar;
+    double fluxtemp1,fluxtemp0,speedtemp;
+    double batemp,bphalf,bmhalf;
+
+    zmhalf = max(q_leftm[2],q_leftp[2]);
+    zphalf = max(q_rightm[2],q_rightp[2]);   
+
+    a_leftm  = q_leftm[0];
+    d_leftm = q_leftm[1];
+    z_leftm  = q_leftm[2];
+    h_leftm  = q_leftm[3];
+    u_leftm  = q_leftm[4];
+    b_leftm  = q_leftm[5];
+    w_leftm  = h_leftm+z_leftm;
+
+    a_leftp  = q_leftp[0];
+    d_leftp = q_leftp[1];
+    z_leftp  = q_leftp[2];
+    h_leftp  = q_leftp[3];
+    u_leftp  = q_leftp[4];
+    b_leftp  = q_leftp[5];
+    w_leftp  = h_leftp+z_leftp;
+
+    a_rightm  = q_rightm[0];
+    d_rightm = q_rightm[1];
+    z_rightm  = q_rightm[2];
+    h_rightm  = q_rightm[3];
+    u_rightm  = q_rightm[4];
+    b_rightm  = q_rightm[5];
+    w_rightm  = h_rightm+z_rightm;
+
+    a_rightp  = q_rightp[0];
+    d_rightp = q_rightp[1];
+    z_rightp  = q_rightp[2];
+    h_rightp  = q_rightp[3];
+    u_rightp  = q_rightp[4];
+    b_rightp  = q_rightp[5];
+    w_rightp  = h_rightp+z_rightp;
+
+    hleftstar = q_leftp[3];
+    hrightstar = q_rightm[3];   
+
+    bphalf = min(b_rightm,b_rightp);
+    bmhalf = min(b_leftm,b_leftp);
+
+    soundspeed_leftp = sqrt(g*h_leftp);
+    soundspeed_leftm = sqrt(g*h_leftm);
+    soundspeed_rightp = sqrt(g*h_rightp);
+    soundspeed_rightm = sqrt(g*h_rightm);
+
+    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
+    if (s_maxl < 0.0) s_maxl = 0.0;
+
+    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
+    if (s_minl > 0.0) s_minl = 0.0;
+
+    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
+    if (s_maxr < 0.0) s_maxr = 0.0;
+
+    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
+    if (s_minr > 0.0) s_minr = 0.0;
+
+    // Flux formulas for left hand side
+    flux_left[0] = d_leftm;
+    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
+
+    flux_right[0] = d_leftp;
+    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
+
+    
+    // Flux computation for left hand side
+    denom = s_maxl-s_minl;
+    if (denom < epsilon) {
+        for (i=0; i<2; i++) edgeflux[i] = 0.0;
+    } else {
+        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
+
+        batemp = (q_leftp[3]+q_leftp[2])*b_leftp-(q_leftm[3]+q_leftm[2])*b_leftm;
+        edgeflux[0] += s_maxl*s_minl*batemp;
+        edgeflux[0] /= denom;
+        edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
+        edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
+        edgeflux[1] /= denom;
+    
+
+    }
+        fluxtemp0 = edgeflux[0];
+        fluxtemp1 = edgeflux[1];
+    
+
+    // Flux formulas for right hand side
+    flux_left[0] = d_rightm;
+    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm;
+
+    flux_right[0] = d_rightp;
+    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp;
+    
+   
+   
+    // Flux computation for right hand side
+    denom = s_maxr-s_minr;
+    if (denom < epsilon) {
+        for (i=0; i<2; i++) edgeflux[i] = 0.0;
+    } else {
+        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
+
+        batemp = (q_rightp[3]+q_rightp[2])*b_rightp-(q_rightm[3]+q_rightm[2])*b_rightm;
+
+        edgeflux[0] += s_maxr*s_minr*batemp;
+        edgeflux[0] /= denom;
+        edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
+        edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
+        edgeflux[1] /= denom;
+    
+
+    }
+    
+   
+    edgeflux[0]=edgeflux[0]-fluxtemp0;
+    edgeflux[1]=edgeflux[1]-fluxtemp1;
+    
+    edgeflux[1]-=-0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightm-z_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightm-b_leftp)*g; 
+
+    //edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
+
+    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
+
+    //edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
+
+   
+    //edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
+      
+     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
+                // Maximal wavespeed
+        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
+            *max_speed = 0.0;
+        }else{
+        speedtemp = max(fabs(s_maxl),fabs(s_minl));
+        speedtemp = max(speedtemp,fabs(s_maxr));
+        speedtemp = max(speedtemp,fabs(s_minr));
+        *max_speed = speedtemp;
+        }
+
+//printf("%f\n",h_right);
+    return 0;
+}
+// NAIEVE VERSION
+int _flux_function_channel2(double *q_leftm,double *q_leftp, double *q_rightm,
+                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed){
+
+    int i;
+    double flux_left[2], flux_right[2];
+    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
+     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
+ double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
+ double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
+    double s_maxl, s_minl,s_maxr,s_minr, denom;
+    double zphalf,zmhalf,hleftstar,hrightstar;
+    double fluxtemp1,fluxtemp0,speedtemp;
+    double batemp,bphalf,bmhalf;
+
+    zmhalf = max(q_leftm[2],q_leftp[2]);
+    zphalf = max(q_rightm[2],q_rightp[2]);   
+
+    a_leftm  = q_leftm[0];
+    d_leftm = q_leftm[1];
+    z_leftm  = q_leftm[2];
+    h_leftm  = q_leftm[3];
+    u_leftm  = q_leftm[4];
+    b_leftm  = q_leftm[5];
+    w_leftm  = h_leftm+z_leftm;
+
+    a_leftp  = q_leftp[0];
+    d_leftp = q_leftp[1];
+    z_leftp  = q_leftp[2];
+    h_leftp  = q_leftp[3];
+    u_leftp  = q_leftp[4];
+    b_leftp  = q_leftp[5];
+    w_leftp  = h_leftp+z_leftp;
+
+    a_rightm  = q_rightm[0];
+    d_rightm = q_rightm[1];
+    z_rightm  = q_rightm[2];
+    h_rightm  = q_rightm[3];
+    u_rightm  = q_rightm[4];
+    b_rightm  = q_rightm[5];
+    w_rightm  = h_rightm+z_rightm;
+
+    a_rightp  = q_rightp[0];
+    d_rightp = q_rightp[1];
+    z_rightp  = q_rightp[2];
+    h_rightp  = q_rightp[3];
+    u_rightp  = q_rightp[4];
+    b_rightp  = q_rightp[5];
+    w_rightp  = h_rightp+z_rightp;
+
+    hleftstar = q_leftp[3];
+    hrightstar = q_rightm[3];   
+
+    bphalf = min(b_rightm,b_rightp);
+    bmhalf = min(b_leftm,b_leftp);
+
+    soundspeed_leftp = sqrt(g*h_leftp);
+    soundspeed_leftm = sqrt(g*h_leftm);
+    soundspeed_rightp = sqrt(g*h_rightp);
+    soundspeed_rightm = sqrt(g*h_rightm);
+
+    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
+    if (s_maxl < 0.0) s_maxl = 0.0;
+
+    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
+    if (s_minl > 0.0) s_minl = 0.0;
+
+    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
+    if (s_maxr < 0.0) s_maxr = 0.0;
+
+    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
+    if (s_minr > 0.0) s_minr = 0.0;
+
+    // Flux formulas for left hand side
+    flux_left[0] = d_leftm;
+    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
+
+    flux_right[0] = d_leftp;
+    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
+
+    
+    // Flux computation for left hand side
+    denom = s_maxl-s_minl;
+    if (denom < epsilon) {
+        for (i=0; i<2; i++) edgeflux[i] = 0.0;
+    } else {
+        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
+
+        batemp = (q_leftp[3]+q_leftp[2])*b_leftp-(q_leftm[3]+q_leftm[2])*b_leftm;
+       
+        edgeflux[0] = 0.5*(flux_left[0]+flux_right[0]);
+        edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]);
+
+    
+
+    }
+        fluxtemp0 = edgeflux[0];
+        fluxtemp1 = edgeflux[1];
+    
+
+    // Flux formulas for right hand side
+    flux_left[0] = d_rightm;
+    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm;
+
+    flux_right[0] = d_rightp;
+    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp;
+    
+   
+   
+    // Flux computation for right hand side
+    denom = s_maxr-s_minr;
+    if (denom < epsilon) {
+        for (i=0; i<2; i++) edgeflux[i] = 0.0;
+    } else {
+        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
+
+        batemp = (q_rightp[3]+q_rightp[2])*b_rightp-(q_rightm[3]+q_rightm[2])*b_rightm;
+
+     
+        edgeflux[0] = 0.5*(flux_right[0]+flux_left[0]);
+
+
+        edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]);
+    
+
+    }
+    
+   
+    edgeflux[0]=edgeflux[0]-fluxtemp0;
+    edgeflux[1]=edgeflux[1]-fluxtemp1;
+    
+    edgeflux[1]-=-0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightm-z_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightm-b_leftp)*g; 
+
+    //edgeflux[1]-=0.5*g*h_rightm*h_rightm*bphalf-0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp-0.5*g*h_leftp*h_leftp*bmhalf;
+
+    // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightm-hleftstar*hleftstar*g*0.5*b_leftp);
+
+    //edgeflux[1]-=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp-2*z_rightm)+hleftstar*(hrightstar+z_leftp-z_rightm))-b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstar-z_leftp+z_rightm)+hleftstar*(hrightstar-2*z_leftp+2*z_rightm)));
+
+   
+    //edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
+      
+     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
+                // Maximal wavespeed
+        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
+            *max_speed = 0.0;
+        }else{
+        speedtemp = max(fabs(s_maxl),fabs(s_minl));
+        speedtemp = max(speedtemp,fabs(s_maxr));
+        speedtemp = max(speedtemp,fabs(s_minr));
+        *max_speed = speedtemp;
+        }
+
+//printf("%f\n",h_right);
+    return 0;
+}