Changeset 8187 for trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c

Timestamp:

Jun 19, 2011, 3:51:16 PM (14 years ago)

Author:

paul

Message:

Rewrote structure of forcing term computation

File:

• trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c (modified) (4 diffs)

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

@@ -12 +12 @@
 void compute_flux_godunov(double[], double[], double, double, double, double, double[]);
 void compute_flux_wb(double[], double[], double, double, double, double, double[]);
-double forcingterm_press_naive(double[], double[], double);
-double forcingterm_fs_naive(double[], double[], double);
+double forcingterm_formula_press_naive(double[], double[], double);
+double forcingterm_formula_fs_naive(double[], double[], double);
+double compute_forcingterm(double[], double[], double);
 int _flux_function_pipe(double *, double *, double *, double *, double, double, double, double *, double *);
 
@@ -23 +24 @@
 #define flux_formula1 flux_formula_press1
 #define compute_flux compute_flux_godunov
-#define forcingterm forcingterm_press_naive
-
+#define forcingterm_formula_press forcingterm_formula_press_naive
+#define forcingterm_formula_fs forcingterm_formula_fs_naive
+#define compute_forcingterm compute_forcingterm_press
 
 // Some formulae
@@ -97 +99 @@
 // *********************************************************************
 
-
-// *********************************************************************
-// Forcing terms
+// *********************************************************************
+// Forcing term formulae
 
 // Pressurised forcing terms
 // Naive version
-double forcingterm_press_naive(double *q_m, double *q_p, double g) {
+double forcingterm_formula_press_naive(double *q_m, double *q_p, double g) {
+  
+  return 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]));
+}
+
+// Free surface forcing terms
+// Naive version
+double forcingterm_formula_fs_naive(double *q_m, double *q_p, double g) {
+  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]);
+}
+// *********************************************************************
+
+// *********************************************************************
+// Forcing term computation
+
+// Pressurised forcing terms computation
+double compute_forcingterm_press(double *q_m, double *q_p, double g) {
   double ret;
-
-  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]));
-
-  // Add in the free surface pressure term with h=t
-  // Does this overwite q_p,q_m in the caller? If so, we don't want that!
+  double h_m, h_p;
+
+  // Compute pressurised forcing term
+  ret = forcingterm_formula_press(q_m, q_p, g);
+
+  // Add in free surface forcing term with h = t
+  // This is a hack: We temporarily change h (q[3]) to t (q[6])
+  // then call fs forcing term, then change h back.
+  h_p = q_p[3];
+  h_m = q_m[3];
   q_p[3] = q_p[6];
   q_m[3] = q_m[6];
-  ret += forcingterm_fs_naive(q_m, q_p, g);
+  ret += forcingterm_formula_fs(q_m, q_p, g);
+  q_p[3] = h_p;
+  q_m[3] = h_m;
 
   return ret;
 }
 
-// Free surface forcing terms
-// Naive version
-double forcingterm_fs_naive(double *q_m, double *q_p, double g) {
-  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]);
+// Free surface forcing terms computation
+double compute_forcingterm_fs(double *q_m, double *q_p, double g) {
+  return forcingterm_formula_fs(q_m, q_p, g);
 }
 // *********************************************************************
@@ -174 +197 @@
 
   // Add in Forcing terms  
-  edgeflux[1] -= forcingterm(q_rightm, q_leftp, g);
+  edgeflux[1] -= compute_forcingterm(q_rightm, q_leftp, g);
 
   // Maximal wavespeed