Changeset 5563

Timestamp:

Jul 23, 2008, 4:38:10 PM (17 years ago)

Author:

steve

Message:

Working version of comp_flux_ext

Location:

anuga_work/development/anuga_1d

Files:

• comp_flux_ext.c (modified) (5 diffs)
• domain.py (modified) (1 diff)
• shallow_water_domain.py (modified) (3 diffs)
• test_shallow_water.py (modified) (3 diffs)
• util_ext.h (added)

anuga_work/development/anuga_1d/comp_flux_ext.c

@@ -5 +5 @@
 const double pi = 3.14159265358979;
 
-double max(double a, double b) {
-        double z;
-        z=(a>b)?a:b;
-        return z;}
-        
-double min(double a, double b) {
-        double z;
-        z=(a<b)?a:b;
-        return z;}
+// Shared code snippets
+#include "util_ext.h"
+
+
+/* double max(double a, double b) { */
+/*      double z; */
+/*      z=(a>b)?a:b; */
+/*      return z;} */
+        
+/* double min(double a, double b) { */
+/*      double z; */
+/*      z=(a<b)?a:b; */
+/*      return z;} */
         
 
@@ -154 +158 @@
                                 
                                 normal = normals[ki];
-                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, flux, &max_speed);
+                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
                                 flux[0] -= edgeflux[0];
                                 flux[1] -= edgeflux[1];
@@ -189 +193 @@
 //=========================================================================
 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
-  PyArrayObject *neighbours, 
+  
+    PyArrayObject 
+        *neighbours, 
         *neighbour_vertices,
-    *normals, 
+        *normals, 
         *areas,
-    *stage_edge_values,
-    *xmom_edge_values,
-    *bed_edge_values,
-    *stage_boundary_values,
-    *xmom_boundary_values,
+        *stage_edge_values,
+        *xmom_edge_values,
+        *bed_edge_values,
+        *stage_boundary_values,
+        *xmom_boundary_values,
         *stage_explicit_update,
         *xmom_explicit_update,
@@ -231 +237 @@
   // the explicit update arrays 
   timestep = _compute_fluxes_ext(timestep,
-                                     epsilon,
-                                     g,
-                                     (long*) neighbours -> data,
-                                     (long*) neighbour_vertices -> data,
-                                     (double*) normals -> data,
-                                     (double*) areas -> data,
-                                     (double*) stage_edge_values -> data,
-                                     (double*) xmom_edge_values -> data,
-                                     (double*) bed_edge_values -> data,
-                                     (double*) stage_boundary_values -> data,
-                                     (double*) xmom_boundary_values -> data,
-                                         (double*) stage_explicit_update -> data,
-                                         (double*) xmom_explicit_update -> data,
-                                         number_of_elements,
-                                         (double*) max_speed_array -> data);
+                                 epsilon,
+                                 g,
+                                 (long*) neighbours -> data,
+                                 (long*) neighbour_vertices -> data,
+                                 (double*) normals -> data,
+                                 (double*) areas -> data,
+                                 (double*) stage_edge_values -> data,
+                                 (double*) xmom_edge_values -> data,
+                                 (double*) bed_edge_values -> data,
+                                 (double*) stage_boundary_values -> data,
+                                 (double*) xmom_boundary_values -> data,
+                                 (double*) stage_explicit_update -> data,
+                                 (double*) xmom_explicit_update -> data,
+                                 number_of_elements,
+                                 (double*) max_speed_array -> data);
+
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
 }
 
+
+PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
+  
+   PyObject 
+        *domain,
+        *stage, 
+        *xmom, 
+        *bed;
+
+    PyArrayObject 
+        *neighbours, 
+        *neighbour_vertices,
+        *normals, 
+        *areas,
+        *stage_vertex_values,
+        *xmom_vertex_values,
+        *bed_vertex_values,
+        *stage_boundary_values,
+        *xmom_boundary_values,
+        *stage_explicit_update,
+        *xmom_explicit_update,
+        *max_speed_array;
+    
+  double timestep, epsilon, g;
+  int number_of_elements;
+
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOO",
+                        &timestep,
+                        &domain,
+                        &stage,
+                        &xmom,
+                        &bed)) {
+      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
+      return NULL;
+  }
+
+
+    epsilon           = get_python_double(domain,"epsilon");
+    g                 = get_python_double(domain,"g");
+ 
+    
+    neighbours        = get_consecutive_array(domain, "neighbours");
+    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
+    normals           = get_consecutive_array(domain, "normals");
+    areas             = get_consecutive_array(domain, "areas");    
+    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
+    
+    stage_vertex_values = get_consecutive_array(stage, "vertex_values");    
+    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");    
+    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");    
+
+    stage_boundary_values = get_consecutive_array(stage, "boundary_values");    
+    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");    
+
+
+    stage_explicit_update = get_consecutive_array(stage, "explicit_update");    
+    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");    
+
+
+
+    number_of_elements = stage_vertex_values -> dimensions[0];
+
+
+ 
+    // Call underlying flux computation routine and update 
+    // the explicit update arrays 
+    timestep = _compute_fluxes_ext(timestep,
+                                 epsilon,
+                                 g,
+                                 (long*) neighbours -> data,
+                                 (long*) neighbour_vertices -> data,
+                                 (double*) normals -> data,
+                                 (double*) areas -> data,
+                                 (double*) stage_vertex_values -> data,
+                                 (double*) xmom_vertex_values -> data,
+                                 (double*) bed_vertex_values -> data,
+                                 (double*) stage_boundary_values -> data,
+                                 (double*) xmom_boundary_values -> data,
+                                 (double*) stage_explicit_update -> data,
+                                 (double*) xmom_explicit_update -> data,
+                                 number_of_elements,
+                                 (double*) max_speed_array -> data);
+
+
+  Py_DECREF(neighbours);
+  Py_DECREF(neighbour_vertices);
+  Py_DECREF(normals);
+  Py_DECREF(areas);
+  Py_DECREF(stage_vertex_values);
+  Py_DECREF(xmom_vertex_values);
+  Py_DECREF(bed_vertex_values);
+  Py_DECREF(stage_boundary_values);
+  Py_DECREF(xmom_boundary_values);
+  Py_DECREF(stage_explicit_update);
+  Py_DECREF(xmom_explicit_update);
+  Py_DECREF(max_speed_array);
+
+
+
+
+  // Return updated flux timestep
+  return Py_BuildValue("d", timestep);
+}
+
  
 
@@ -259 +372 @@
 static struct PyMethodDef MethodTable[] = {
   {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
+  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
   {NULL, NULL}
 };

anuga_work/development/anuga_1d/domain.py

@@ -47 +47 @@
         self.centroids = zeros(N, Float)
         self.areas     = zeros(N, Float)
+
+        self.max_speed_array = zeros(N, Float)
 
         self.normals = zeros((N, 2), Float)

anuga_work/development/anuga_1d/shallow_water_domain.py

@@ -90 +90 @@
 
         self.quantities_to_be_stored = ['stage','xmomentum']
+
+
 
 
@@ -155 +157 @@
         #Call correct module function
         #(either from this module or C-extension)
-        compute_fluxes(self)
+        compute_fluxes_C_short(self)
 
     def compute_timestep(self):
@@ -572 +574 @@
     #print domain.quantities['stage'].centroid_values
 
+
+# Compute flux definition
+def compute_fluxes_C_long(domain):
+        from Numeric import zeros, Float
+        import sys
+
+        
+        timestep = float(sys.maxint)
+        #print 'timestep=',timestep
+        #print 'The type of timestep is',type(timestep)
+        
+        epsilon = domain.epsilon
+        #print 'epsilon=',epsilon
+        #print 'The type of epsilon is',type(epsilon)
+        
+        g = domain.g
+        #print 'g=',g
+        #print 'The type of g is',type(g)
+        
+        neighbours = domain.neighbours
+        #print 'neighbours=',neighbours
+        #print 'The type of neighbours is',type(neighbours)
+        
+        neighbour_vertices = domain.neighbour_vertices
+        #print 'neighbour_vertices=',neighbour_vertices
+        #print 'The type of neighbour_vertices is',type(neighbour_vertices)
+        
+        normals = domain.normals
+        #print 'normals=',normals
+        #print 'The type of normals is',type(normals)
+        
+        areas = domain.areas
+        #print 'areas=',areas
+        #print 'The type of areas is',type(areas)
+        
+        stage_edge_values = domain.quantities['stage'].vertex_values
+        #print 'stage_edge_values=',stage_edge_values
+        #print 'The type of stage_edge_values is',type(stage_edge_values)
+        
+        xmom_edge_values = domain.quantities['xmomentum'].vertex_values
+        #print 'xmom_edge_values=',xmom_edge_values
+        #print 'The type of xmom_edge_values is',type(xmom_edge_values)
+        
+        bed_edge_values = domain.quantities['elevation'].vertex_values
+        #print 'bed_edge_values=',bed_edge_values
+        #print 'The type of bed_edge_values is',type(bed_edge_values)
+        
+        stage_boundary_values = domain.quantities['stage'].boundary_values
+        #print 'stage_boundary_values=',stage_boundary_values
+        #print 'The type of stage_boundary_values is',type(stage_boundary_values)
+        
+        xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
+        #print 'xmom_boundary_values=',xmom_boundary_values
+        #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
+        
+        stage_explicit_update = domain.quantities['stage'].explicit_update
+        #print 'stage_explicit_update=',stage_explicit_update
+        #print 'The type of stage_explicit_update is',type(stage_explicit_update)
+        
+        xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
+        #print 'xmom_explicit_update=',xmom_explicit_update
+        #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
+        
+        number_of_elements = len(stage_edge_values)
+        #print 'number_of_elements=',number_of_elements
+        #print 'The type of number_of_elements is',type(number_of_elements)
+        
+        max_speed_array = domain.max_speed_array
+        #print 'max_speed_array=',max_speed_array
+        #print 'The type of max_speed_array is',type(max_speed_array)
+        
+        
+        from comp_flux_ext import compute_fluxes_ext
+        
+        domain.timestep = compute_fluxes_ext(timestep,
+                                  epsilon,
+                                  g,
+                                  neighbours,
+                                  neighbour_vertices,
+                                  normals,
+                                  areas, 
+                                  stage_edge_values,
+                                  xmom_edge_values,
+                                  bed_edge_values,
+                                  stage_boundary_values,
+                                  xmom_boundary_values,
+                                  stage_explicit_update,
+                                  xmom_explicit_update,
+                                  number_of_elements,
+                                  max_speed_array)
+
+
+# Compute flux definition
+def compute_fluxes_C_short(domain):
+        from Numeric import zeros, Float
+        import sys
+
+        
+        timestep = float(sys.maxint)
+
+        stage = domain.quantities['stage']
+        xmom = domain.quantities['xmomentum']
+        bed = domain.quantities['elevation']
+
+
+        from comp_flux_ext import compute_fluxes_ext_short
+
+        domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
+
+
+
 ####################################
+
+
+
+
+
 
 # Module functions for gradient limiting (distribute_to_vertices_and_edges)

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -41 +41 @@
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
         assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([2.0, 8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([-2.0, -8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+    def test_domain_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        #print edgeflux
+
+        from shallow_water_domain import flux_function as domain_flux_function
+
+        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
+
+        #print domainedgeflux
+
+        assert allclose(domainedgeflux, edgeflux)
+        assert allclose(domainmaxspeed, maxspeed)
+
 
 
@@ -128 +173 @@
                 #m = domain.neighbour_edges[k,i]
                 m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
                 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                 qr[0] = stage[n, m]
@@ -169 +215 @@
 
 
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
 
 #-------------------------------------------------------------