Changeset 267

Timestamp:

Sep 2, 2004, 3:00:07 AM (21 years ago)

Author:

ole

Message:

C implementation of balanced_deep_and_shallow.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• doit.py (modified) (1 diff)
• shallow_water.py (modified) (4 diffs)
• shallow_water_ext.c (modified) (5 diffs)

inundation/ga/storm_surge/pyvolution/doit.py

@@ -1 +1 @@
 import os
+
+#Clean up
+
+for file in os.listdir('.'):
+    if file[-1] == '~':
+        os.remove(file)
+
 
 os.system('python compile.py')

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -511 +511 @@
         hmin = min( hv[k,:] )
 
+        print dz, hmin
+        
         #Create alpha in [0,1], where alpha==0 means using shallow 
         #first order scheme and alpha==1 means using the stage w as
@@ -552 +554 @@
             # xmomc and ymomc (shallow) and momentum
             # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];            
-            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]            
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
+                                    
             
 
@@ -580 +583 @@
 
     from shallow_water_ext import balance_deep_and_shallow
-    balance_deep_and_shallow(domain.number_of_elements,
-                             wc, zc, wv, zv, hv,
+    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
                              xmomc, ymomc, xmomv, ymomv)
 
@@ -985 +987 @@
     gravity = gravity_c
     manning_friction = manning_friction_c
-    #balance_deep_and_shallow = balance_deep_and_shallow_c 
+    balance_deep_and_shallow = balance_deep_and_shallow_c 
     
     #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -176 +176 @@
             
 
-/*
-int _balance_deep_and_shallow()     
-
-    #Computed linear combination between constant levels and and
-    #levels parallel to the bed elevation.     
-    for k in range(domain.number_of_elements):
-        #Compute maximal variation in bed elevation
-        #  This quantitiy is
-        #    dz = max_i abs(z_i - z_c)
-        #  and it is independent of dimension
-        #  In the 1d case zc = (z0+z1)/2
-        #  In the 2d case zc = (z0+z1+z2)/3
-
-        dz = max(abs(zv[k,0]-zc[k]),
-                 abs(zv[k,1]-zc[k]),
-                 abs(zv[k,2]-zc[k]))
-
-        
-        hmin = min( hv[k,:] )
-
-        #Create alpha in [0,1], where alpha==0 means using shallow 
-        #first order scheme and alpha==1 means using the stage w as
-        #computed by the gradient limiter (1st or 2nd order)
-        #
-        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
-        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
+
+int _balance_deep_and_shallow(int N,
+                              double* wc,
+                              double* zc, 
+                              double* hc,                             
+                              double* wv, 
+                              double* zv, 
+                              double* hv,
+                              double* xmomc, 
+                              double* ymomc, 
+                              double* xmomv, 
+                              double* ymomv) { 
+  
+  int k, k3, i;
+  double dz, hmin, alpha;
+  
+  //Compute linear combination between constant levels and and
+  //levels parallel to the bed elevation.     
+  
+  for (k=0; k<N; k++) {
+    // Compute maximal variation in bed elevation
+    // This quantitiy is
+    //     dz = max_i abs(z_i - z_c)
+    // and it is independent of dimension
+    // In the 1d case zc = (z0+z1)/2
+    // In the 2d case zc = (z0+z1+z2)/3
+
+    k3 = 3*k;
+    
+    //FIXME: Try with this one precomputed
+    dz = 0.0;
+    hmin = hv[k3];
+    for (i=0; i<3; i++) {
+      dz = max(dz, fabs(zv[k3+i]-zc[k]));
+      hmin = min(hmin, hv[k3+i]);
+    }
+
+    
+    //Create alpha in [0,1], where alpha==0 means using shallow 
+    //first order scheme and alpha==1 means using the stage w as
+    //computed by the gradient limiter (1st or 2nd order)
+    //
+    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
+    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
+    
+    if (dz > 0.0) 
+      alpha = max( min( 2*hmin/dz, 1.0), 0.0 );
+    else
+      alpha = 1.0;  //Flat bed
+
+
+    //Weighted balance between stage parallel to bed elevation 
+    //(wvi = zvi + hc) and stage as computed by 1st or 2nd 
+    //order gradient limiter
+    //(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
+    //
+    //It follows that the updated wvi is
+    //  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
+    //  zvi + hc + alpha*(hvi - hc)
+    //
+    //Note that hvi = zc+hc-zvi in the first order case (constant).
+
+    if (alpha < 1) {         
+      for (i=0; i<3; i++) {
+        wv[k3+i] = zv[k3+i] + hc[k] + alpha*(hv[k3+i]-hc[k]);
+            
       
-        if dz > 0.0:
-            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
-        else:
-            #Flat bed 
-            alpha = 1.0
-
-
-        #Weighted balance between stage parallel to bed elevation 
-        #(wvi = zvi + hc) and stage as computed by 1st or 2nd 
-        #order gradient limiter
-        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
-        #
-        #It follows that the updated wvi is
-        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
-        #  zvi + hc + alpha*(hvi - hc)
-        #
-        #Note that hvi = zc+hc-zvi in the first order case (constant).
-
-        if alpha < 1:         
-            for i in range(3):
-                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-            
-
-            # Update momentum as a linear combination of 
-            # xmomc and ymomc (shallow) and momentum
-            # from extrapolator xmomv and ymomv (deep).
-            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];            
-            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
-*/            
-            
-        
+        //Update momentum as a linear combination of 
+        //xmomc and ymomc (shallow) and momentum
+        //from extrapolator xmomv and ymomv (deep).
+        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];            
+        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+      } 
+    }
+  }         
+  return 0;
+}
+
 ///////////////////////////////////////////////////////////////////
 // Gateways to Python
@@ -542 +562 @@
   //
   //    balance_deep_and_shallow(domain.number_of_elements,
-  //                             wc, zc, wv, zv, hv,
+  //                             wc, zc, hc, wv, zv, hv,
   //                             xmomc, ymomc, xmomv, ymomv)
   
@@ -549 +569 @@
     *wc,            //Level at centroids
     *zc,            //Elevation at centroids    
+    *hc,            //Height at centroids        
     *wv,            //Level at vertices
     *zv,            //Elevation at vertices
@@ -560 +581 @@
   
   // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv,
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",  
+                        &wc, &zc, &hc, 
+                        &wv, &zv, &hv,
                         &xmomc, &ymomc, &xmomv, &ymomv))
     return NULL;
 
-  
-  /*
+  N = wc -> dimensions[0];
+    
   _balance_deep_and_shallow(N,
                             (double*) wc -> data,
                             (double*) zc -> data, 
+                            (double*) hc -> data,                           
                             (double*) wv -> data, 
                             (double*) zv -> data, 
@@ -576 +600 @@
                             (double*) xmomv -> data, 
                             (double*) ymomv -> data);  
-  */
+  
  
   return Py_BuildValue("");