Changeset 515

Timestamp:

Nov 9, 2004, 10:41:43 PM (20 years ago)

Author:

ole

Message:

Identified potential optimisation in Manning friction using Taylor expansion

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity_ext.c (modified) (1 diff)
• run_profile.py (modified) (1 diff)
• shallow_water.py (modified) (3 diffs)
• shallow_water_ext.c (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -169 +169 @@
   double denominator, x;
 
+  
   //Divide semi_implicit update by conserved quantity  
   for (k=0; k<N; k++) {  

inundation/ga/storm_surge/pyvolution/run_profile.py

@@ -20 +20 @@
 #
 
-N = 32
+N = 16
 
 print 'Creating domain'

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -819 +819 @@
     
     for k in range(N):
-        if h[k] >= eps:
-            S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
-            S /= h[k]**(7.0/3)
-
-            #Update momentum
-            xmom_update[k] += S*uh[k]
-            ymom_update[k] += S*vh[k]
+        if eta[k] >= eps:
+            if h[k] >= eps:
+                S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
+                S /= h[k]**(7.0/3)
+
+                #Update momentum
+                xmom_update[k] += S*uh[k]
+                ymom_update[k] += S*vh[k]
             
         
@@ -838 +839 @@
     w = domain.quantities['level'].centroid_values
     z = domain.quantities['elevation'].centroid_values
-    h = w-z
+    
     uh = xmom.centroid_values
     vh = ymom.centroid_values
@@ -851 +852 @@
 
     from shallow_water_ext import manning_friction
-    manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
+    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
         
 

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -156 +156 @@
         
 void _manning_friction(double g, double eps, int N, 
-                       double* h, double* uh, double* vh, 
+                       double* w, double* z, 
+                       double* uh, double* vh, 
                        double* eta, double* xmom, double* ymom) {     
 
   int k;
-  double S;
+  double S, h;
   
   for (k=0; k<N; k++) {
-    if (h[k] >= eps) {
-      S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-      S /= pow(h[k], 7.0/3);      
-
-      //Update momentum
-      xmom[k] += S*uh[k];
-      ymom[k] += S*vh[k];
+    if (eta[k] > eps) {
+      h = w[k]-z[k];
+      if (h >= eps) {
+        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 7.0/3);      //Expensive
+        //S /= h*h*(1 + h/3.0 - h*h/9.0);  //FIXME: Use a Taylor expansion
+        //FIXME: This will save a lot of time
+
+        //Update momentum
+        xmom[k] += S*uh[k];
+        ymom[k] += S*vh[k];
+      }
     }
   }
-  
 }
             
@@ -344 +349 @@
   
   
-  PyArrayObject *h, *uh, *vh, *eta, *xmom, *ymom;
+  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
   int N;
   double g, eps;
     
-  if (!PyArg_ParseTuple(args, "ddOOOOOO",
-                        &g, &eps, &h, &uh, &vh, &eta, 
+  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
+                        &g, &eps, &w, &z, &uh, &vh, &eta, 
                         &xmom, &ymom))  
     return NULL;  
 
-  N = h -> dimensions[0];    
+  N = w -> dimensions[0];    
   _manning_friction(g, eps, N,
-                    (double*) h -> data,
+                    (double*) w -> data,
+                    (double*) z -> data,                    
                     (double*) uh -> data, 
                     (double*) vh -> data,