Changeset 266

Timestamp:

Sep 1, 2004, 6:31:42 PM (21 years ago)

Author:

ole

Message:

 

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (3 diffs)
• shallow_water_ext.c (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -474 +474 @@
 
 def balance_deep_and_shallow(domain):
-    """Compute linear combination between stage as computed by gradient-limiters and
-    stage computed as constant height above bed.
-    The former takes precedence when heights are large compared to the bed slope while the latter
-    takes precedence when heights are relatively small.
-    Anything in between is computed as a balanced linear combination in order to avoid numerical
-    disturbances which would otherwise appear as a result of hard switching between modes.
+    """Compute linear combination between stage as computed by
+    gradient-limiters and stage computed as constant height above bed.
+    The former takes precedence when heights are large compared to the
+    bed slope while the latter takes precedence when heights are
+    relatively small.  Anything in between is computed as a balanced
+    linear combination in order to avoid numerical disturbances which
+    would otherwise appear as a result of hard switching between
+    modes.
     """
     
@@ -553 +555 @@
             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
             
+
+def balance_deep_and_shallow_c(domain):
+    """Wrapper for C implementation
+    """
+    
+    #Shortcuts
+    wc = domain.quantities['level'].centroid_values
+    zc = domain.quantities['elevation'].centroid_values
+    hc = wc - zc
+    
+    wv = domain.quantities['level'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    hv = wv-zv
+
+    #Momentums at centroids
+    xmomc = domain.quantities['xmomentum'].centroid_values
+    ymomc = domain.quantities['ymomentum'].centroid_values        
+
+    #Momentums at vertices
+    xmomv = domain.quantities['xmomentum'].vertex_values
+    ymomv = domain.quantities['ymomentum'].vertex_values         
+
+    
+
+    from shallow_water_ext import balance_deep_and_shallow
+    balance_deep_and_shallow(domain.number_of_elements,
+                             wc, zc, wv, zv, hv,
+                             xmomc, ymomc, xmomv, ymomv)
+
     
 
@@ -954 +985 @@
     gravity = gravity_c
     manning_friction = manning_friction_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.
+      
+        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,:];                     
+*/            
+            
         
 ///////////////////////////////////////////////////////////////////
@@ -482 +539 @@
 
 
+PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
+  //
+  //    balance_deep_and_shallow(domain.number_of_elements,
+  //                             wc, zc, wv, zv, hv,
+  //                             xmomc, ymomc, xmomv, ymomv)
+  
+
+  PyArrayObject 
+    *wc,            //Level at centroids
+    *zc,            //Elevation at centroids    
+    *wv,            //Level at vertices
+    *zv,            //Elevation at vertices
+    *hv,            //Heights at vertices    
+    *xmomc,         //Momentums at centroids and vertices
+    *ymomc, 
+    *xmomv, 
+    *ymomv;    
+    
+  int N; //, err;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv,
+                        &xmomc, &ymomc, &xmomv, &ymomv))
+    return NULL;
+
+  
+  /*
+  _balance_deep_and_shallow(N,
+                            (double*) wc -> data,
+                            (double*) zc -> data, 
+                            (double*) wv -> data, 
+                            (double*) zv -> data, 
+                            (double*) hv -> data,
+                            (double*) xmomc -> data, 
+                            (double*) ymomc -> data, 
+                            (double*) xmomv -> data, 
+                            (double*) ymomv -> data);  
+  */
+ 
+  return Py_BuildValue("");  
+}
 
 
@@ -499 +597 @@
   {"gravity", gravity, METH_VARARGS, "Print out"},      
   {"manning_friction", manning_friction, METH_VARARGS, "Print out"},        
+  {"balance_deep_and_shallow", balance_deep_and_shallow, 
+   METH_VARARGS, "Print out"},          
   //{"distribute_to_vertices_and_edges", 
   // distribute_to_vertices_and_edges, METH_VARARGS},