Changeset 5162

Timestamp:

Mar 11, 2008, 8:43:22 PM (17 years ago)

Author:

steve

Message:

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py (modified) (3 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (26 diffs)
• anuga_core/source/anuga/advection/advection_ext.c (modified) (6 diffs)
• anuga_core/source/anuga/config.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (8 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (6 diffs)
• anuga_core/source/anuga/shallow_water/test_data_manager.py (modified) (2 diffs)
• anuga_core/source/anuga/utilities/util_ext.h (modified) (4 diffs)
• anuga_core/source/anuga/utilities/xml_tools.py (modified) (1 diff)
• anuga_validation/convergence_study/animatesww2d_alt.py (modified) (2 diffs)
• anuga_validation/convergence_study/convergence_structured.py (modified) (5 diffs)
• anuga_validation/convergence_study/dam_break.py (added)
• anuga_validation/convergence_study/wave.py (added)
• anuga_work/development/analytical_solutions/oblique_shock.py (modified) (5 diffs)
• anuga_work/development/near_shore_PMD/run_beach.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -10 +10 @@
 from Numeric import allclose, argmax, zeros, Float
 from anuga.config import epsilon
+from anuga.config import beta_euler, beta_rk2
 
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
@@ -301 +302 @@
         self.default_order = n
         self._order_ = self.default_order
+
+        if self.default_order == 1:
+            self.set_timestepping_method('euler')
+            #self.set_all_limiters(beta_euler)
+
+        if self.default_order == 2:
+            self.set_timestepping_method('rk2')
+            #self.set_all_limiters(beta_rk2)
         
 
@@ -1550 +1559 @@
             elif self._order_ == 2:
                 Q.extrapolate_second_order()
-                Q.limit()
+                #Q.limit()
             else:
                 raise 'Unknown order'
-            Q.interpolate_from_vertices_to_edges()
+            #Q.interpolate_from_vertices_to_edges()
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -62 +62 @@
         #Allocate space for Gradient
         self.x_gradient = zeros(N, Float)
-        self.y_gradient = zeros(N, Float)        
+        self.y_gradient = zeros(N, Float)
+
+        #Allocate space for Limiter Phi
+        self.phi = zeros(N, Float)        
 
         #Intialise centroid and edge_values
@@ -74 +77 @@
         #flux calculations and forcing functions
 
-        N = len(domain) # number_of_triangles
+        #Allocate space for update fields
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
-        self.centroid_backup_values = zeros(N, Float)       
+        self.centroid_backup_values = zeros(N, Float)
+
+        self.beta = 1.0
 
 
@@ -1403 +1408 @@
         compute_gradients(self)
         extrapolate_from_gradient(self)
-
+        
     def extrapolate_second_order_and_limit(self):
         #Call correct module function
         #(either from this module or C-extension)
-        extrapolate_second_order_and_limit(self)        
+        extrapolate_second_order_and_limit(self)
+
+    def bound_vertices_below_by_constant(self, bound):
+        #Call correct module function
+        #(either from this module or C-extension)
+        bound_vertices_below_by_constant(self, bound)
+
+    def bound_vertices_below_by_quantity(self, quantity):
+        #Call correct module function
+        #(either from this module or C-extension)
+
+        #check consistency
+        assert self.domain == quantity.domain
+        bound_vertices_below_by_quantity(self, quantity)                        
 
     def backup_centroid_values(self):
@@ -1450 +1468 @@
          limit_gradient_by_neighbour,\
          extrapolate_from_gradient,\
+         extrapolate_second_order_and_limit,\
+         bound_vertices_below_by_constant,\
+         bound_vertices_below_by_quantity,\
          interpolate_from_vertices_to_edges,\
          interpolate_from_edges_to_vertices,\

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -145 +145 @@
                                          double* vertex_values,
                                          double* edge_values,
+                                         double* phi,
                                          double* x_gradient,
                                          double* y_gradient) {
@@ -151 +152 @@
         double x, y, x0, y0, x1, y1, x2, y2;
         long n;
-        double qmin, qmax, qn, qc;
-        double dq, dqa[3], phi, r;
-
+        double qmin, qmax, qc;
+        double qn[3];
+        double dq, dqa[3], r;
 
         for (k=0; k<N; k++){
@@ -161 +162 @@
 
                 // Centroid coordinates
-                x = centroids[k2]; y = centroids[k2+1];
+                x = centroids[k2+0]; 
+                y = centroids[k2+1];
 
                 // vertex coordinates
                 // x0, y0, x1, y1, x2, y2 = X[k,:]
-                x0 = vertex_coordinates[k6 + 0];
+                x0 = vertex_coordinates[k6 + 0];
                 y0 = vertex_coordinates[k6 + 1];
                 x1 = vertex_coordinates[k6 + 2];
@@ -172 +174 @@
                 y2 = vertex_coordinates[k6 + 5];
 
-                //Centroid Value
-                qc = centroid_values[k];
-
                 // Extrapolate to Vertices
                 vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
@@ -184 +183 @@
                 edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
                 edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
-
-
-                phi = 1.0;
+        }
+
+
+
+        for (k=0; k<N; k++){
+                k6 = 6*k;
+                k3 = 3*k;
+                k2 = 2*k;
+
+
+                qc = centroid_values[k];
+                
+                qmin = qc;
+                qmax = qc;
+                
+                for (i=0; i<3; i++) {
+                    n = neighbours[k3+i];
+                    if (n < 0) {
+                        qn[i] = qc;
+                    } else {
+                        qn[i] = centroid_values[n];
+                    }
+
+                    qmin = min(qmin, qn[i]);
+                    qmax = max(qmax, qn[i]);
+                }
+
+                //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
+                //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
+
+/*              for (i=0; i<3; i++) { */
+/*                  n = neighbours[k3+i]; */
+/*                  if (n < 0) { */
+/*                      qn[i] = qc; */
+/*                      qmin[i] = qtmin; */
+/*                      qmax[i] = qtmax; */
+/*                  }  */
+/*              } */
+                        
+                phi[k] = 1.0;
 
                 for (i=0; i<3; i++) {
@@ -192 +228 @@
                     dqa[i] = dq;                      //Save dq for use in updating vertex values
 
-                    n = neighbours[k3+i];
-                    if (n >= 0) {
-                        qn = centroid_values[n];     //Neighbour's centroid value
-
-                        qmin = min(qc, qn);
-                        qmax = max(qc, qn);
-
-
-                        r = 1.0;
+                    r = 1.0;
       
-                        if (dq > 0.0) r = (qmax - qc)/dq;
-                        if (dq < 0.0) r = (qmin - qc)/dq;      
+                    if (dq > 0.0) r = (qmax - qc)/dq;
+                    if (dq < 0.0) r = (qmin - qc)/dq;
                     
-                        phi = min( min(r*beta, 1.0), phi);    
-                    }
+                    phi[k] = min( min(r*beta, 1.0), phi[k]);
+                    
                 }
 
 
+
                 //Update gradient, edge and vertex values using phi limiter
-                x_gradient[k] = x_gradient[k]*phi;
-                y_gradient[k] = y_gradient[k]*phi;
-
-                edge_values[k3+0] = qc + phi*dqa[0];
-                edge_values[k3+1] = qc + phi*dqa[1];
-                edge_values[k3+2] = qc + phi*dqa[2];
+                x_gradient[k] = x_gradient[k]*phi[k];
+                y_gradient[k] = y_gradient[k]*phi[k];
+
+                edge_values[k3+0] = qc + phi[k]*dqa[0];
+                edge_values[k3+1] = qc + phi[k]*dqa[1];
+                edge_values[k3+2] = qc + phi[k]*dqa[2];
 
 
@@ -222 +251 @@
                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+                
 
         }
@@ -290 +320 @@
 }
 
+
+
+
 int _limit_edges_by_all_neighbours(int N, double beta,
-                     double* centroid_values,
-                     double* vertex_values,
-                     double* edge_values,
-                     long*   neighbours) {
+                                   double* centroid_values,
+                                   double* vertex_values,
+                                   double* edge_values,
+                                   long*   neighbours,
+                                   double* x_gradient,
+                                   double* y_gradient) {
 
         int i, k, k2, k3, k6;
@@ -334 +369 @@
                 }
     
-                //Update vertex and edge values using phi limiter
+                //Update gradient, vertex and edge values using phi limiter
+                x_gradient[k] = x_gradient[k]*phi;
+                y_gradient[k] = y_gradient[k]*phi;
+
                 edge_values[k3+0] = qc + phi*dqa[0];
                 edge_values[k3+1] = qc + phi*dqa[1];
@@ -353 +391 @@
                      double* vertex_values,
                      double* edge_values,
+                     long*   neighbours) {
+
+        int i, k, k2, k3, k6;
+        long n;
+        double qmin, qmax, qn, qc;
+        double dq, dqa[3], phi, r;
+
+        for (k=0; k<N; k++){
+                k6 = 6*k;
+                k3 = 3*k;
+                k2 = 2*k;
+
+                qc = centroid_values[k];
+                phi = 1.0;
+
+                for (i=0; i<3; i++) {
+                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+                    dqa[i] = dq;                      //Save dqa for use in updating vertex values
+
+                    n = neighbours[k3+i];
+                    qn = qc;
+                    if (n >= 0)  qn = centroid_values[n]; //Neighbour's centroid value
+
+                    qmin = min(qc, qn);
+                    qmax = max(qc, qn);
+
+                    r = 1.0;
+      
+                    if (dq > 0.0) r = (qmax - qc)/dq;
+                    if (dq < 0.0) r = (qmin - qc)/dq;      
+                    
+                    phi = min( min(r*beta, 1.0), phi);    
+                   
+                }
+
+
+                //Update edge and vertex values using phi limiter
+                edge_values[k3+0] = qc + phi*dqa[0];
+                edge_values[k3+1] = qc + phi*dqa[1];
+                edge_values[k3+2] = qc + phi*dqa[2];
+                
+                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
+                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
+                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+
+        }
+
+        return 0;
+}
+
+
+int _limit_gradient_by_neighbour(int N, double beta,
+                     double* centroid_values,
+                     double* vertex_values,
+                     double* edge_values,
+                     double* x_gradient,
+                     double* y_gradient,
                      long*   neighbours) {
 
@@ -403 +498 @@
 }
 
-
-int _limit_gradient_by_neighbour(int N, double beta,
+int _bound_vertices_below_by_constant(int N, double bound,
                      double* centroid_values,
                      double* vertex_values,
                      double* edge_values,
                      double* x_gradient,
-                     double* y_gradient,
-                     long*   neighbours) {
+                     double* y_gradient) {
 
         int i, k, k2, k3, k6;
-        long n;
-        double qmin, qmax, qn, qc;
+        double qmin, qc;
         double dq, dqa[3], phi, r;
 
@@ -423 +515 @@
 
                 qc = centroid_values[k];
+                qmin = bound;
+
+
                 phi = 1.0;
-
-                for (i=0; i<3; i++) {
-                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+                for (i=0; i<3; i++) {    
+                    r = 1.0;
+      
+                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
                     dqa[i] = dq;                      //Save dq for use in updating vertex values
-
-                    n = neighbours[k3+i];
-                    if (n >= 0) {
-                        qn = centroid_values[n]; //Neighbour's centroid value
-
-                        qmin = min(qc, qn);
-                        qmax = max(qc, qn);
-
-                        r = 1.0;
       
-                        if (dq > 0.0) r = (qmax - qc)/dq;
-                        if (dq < 0.0) r = (qmin - qc)/dq;      
+                    if (dq < 0.0) r = (qmin - qc)/dq;      
+  
                     
-                        phi = min( min(r*beta, 1.0), phi);    
-                    }
+                    phi = min( min(r, 1.0), phi);    
                 }
-
-
-                //Update edge and vertex values using phi limiter
-                edge_values[k3+0] = qc + phi*dqa[0];
-                edge_values[k3+1] = qc + phi*dqa[1];
-                edge_values[k3+2] = qc + phi*dqa[2];
+    
+
+                //Update gradient, vertex and edge values using phi limiter
+                x_gradient[k] = x_gradient[k]*phi;
+                y_gradient[k] = y_gradient[k]*phi;
+
+                vertex_values[k3+0] = qc + phi*dqa[0];
+                vertex_values[k3+1] = qc + phi*dqa[1];
+                vertex_values[k3+2] = qc + phi*dqa[2];
                 
-                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
-                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
-                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
+                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
+                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
 
         }
@@ -460 +549 @@
 }
 
-
+int _bound_vertices_below_by_quantity(int N,
+                                      double* bound_vertex_values,
+                                      double* centroid_values,
+                                      double* vertex_values,
+                                      double* edge_values,
+                                      double* x_gradient,
+                                      double* y_gradient) {
+
+        int i, k, k2, k3, k6;
+        double qmin, qc;
+        double dq, dqa[3], phi, r;
+
+        for (k=0; k<N; k++){
+                k6 = 6*k;
+                k3 = 3*k;
+                k2 = 2*k;
+
+                qc = centroid_values[k];
+
+                phi = 1.0;
+                for (i=0; i<3; i++) {    
+                    r = 1.0;
+      
+                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
+                    dqa[i] = dq;                      //Save dq for use in updating vertex values
+      
+                    qmin = bound_vertex_values[k3+i];
+                    if (dq < 0.0) r = (qmin - qc)/dq;      
+  
+                    
+                    phi = min( min(r, 1.0), phi);    
+                }
+    
+
+                //Update gradient, vertex and edge values using phi limiter
+                x_gradient[k] = x_gradient[k]*phi;
+                y_gradient[k] = y_gradient[k]*phi;
+
+                vertex_values[k3+0] = qc + phi*dqa[0];
+                vertex_values[k3+1] = qc + phi*dqa[1];
+                vertex_values[k3+2] = qc + phi*dqa[2];
+                
+                edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
+                edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
+                edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
+
+        }
+
+        return 0;
+}
 
 int _interpolate_from_vertices_to_edges(int N,
@@ -906 +1044 @@
 
 
-/* PyObject *compute_gradients(PyObject *self, PyObject *args) { */
-/*   //"""Compute gradients of triangle surfaces defined by centroids of */
-/*   //neighbouring volumes. */
-/*   //If one edge is on the boundary, use own centroid as neighbour centroid. */
-/*   //If two or more are on the boundary, fall back to first order scheme. */
-/*   //""" */
-
-
-/*      PyObject *quantity, *domain, *R; */
-/*      PyArrayObject */
-/*              *centroids,            //Coordinates at centroids */
-/*              *centroid_values,      //Values at centroids */
-/*              *number_of_boundaries, //Number of boundaries for each triangle */
-/*              *surrogate_neighbours, //True neighbours or - if one missing - self */
-/*              *a, *b;                //Return values */
-
-/*      int dimensions[1], N, err; */
-
-/*      // Convert Python arguments to C */
-/*      if (!PyArg_ParseTuple(args, "O", &quantity)) { */
-/*        PyErr_SetString(PyExc_RuntimeError,  */
-/*                        "quantity_ext.c: compute_gradients could not parse input");    */
-/*        return NULL; */
-/*      } */
-
-/*      domain = PyObject_GetAttrString(quantity, "domain"); */
-/*      if (!domain) { */
-/*        PyErr_SetString(PyExc_RuntimeError,  */
-/*                        "compute_gradients could not obtain domain object from quantity"); */
-/*        return NULL; */
-/*      } */
-
-/*      // Get pertinent variables */
-
-/*      centroids = get_consecutive_array(domain, "centroid_coordinates"); */
-/*      centroid_values = get_consecutive_array(quantity, "centroid_values"); */
-/*      surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); */
-/*      number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); */
-
-/*      N = centroid_values -> dimensions[0]; */
-
-/*      // Release */
-/*      Py_DECREF(domain); */
-
-/*      // Allocate space for return vectors a and b (don't DECREF) */
-/*      dimensions[0] = N; */
-/*      a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
-/*      b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
-
-
-
-/*      err = _compute_gradients(N, */
-/*                      (double*) centroids -> data, */
-/*                      (double*) centroid_values -> data, */
-/*                      (long*) number_of_boundaries -> data, */
-/*                      (long*) surrogate_neighbours -> data, */
-/*                      (double*) a -> data, */
-/*                      (double*) b -> data); */
-
-/*      if (err != 0) { */
-/*        PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
-/*        return NULL; */
-/*      } */
-
-/*      // Release */
-/*      Py_DECREF(centroids); */
-/*      Py_DECREF(centroid_values); */
-/*      Py_DECREF(number_of_boundaries); */
-/*      Py_DECREF(surrogate_neighbours); */
-
-/*      // Build result, release and return */
-/*      R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); */
-/*      Py_DECREF(a); */
-/*      Py_DECREF(b); */
-/*      return R; */
-/* } */
-
-
-
 PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
 
@@ -1087 +1146 @@
         *quantity_vertex_values,        //Values at vertices
         *quantity_edge_values,          //Values at edges
+        *quantity_phi,                  //limiter phi values
         *quantity_x_gradient,           //x gradient
         *quantity_y_gradient;           //y gradient
@@ -1097 +1157 @@
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOd", 
-                        &domain,
-                        &quantity,
-                        &beta)) {
+  if (!PyArg_ParseTuple(args, "O",&quantity)) {
       PyErr_SetString(PyExc_RuntimeError, 
                       "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");      
+      return NULL;
+  }
+
+  domain = PyObject_GetAttrString(quantity, "domain");
+  if (!domain) {
+      PyErr_SetString(PyExc_RuntimeError, 
+                      "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
       return NULL;
   }
@@ -1117 +1181 @@
   quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
   quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_phi                  = get_consecutive_array(quantity, "phi");
   quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
   quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
 
+  beta = get_python_double(quantity,"beta");
+
   ntri = quantity_centroid_values -> dimensions[0];
-
-
 
   err = _compute_gradients(ntri,
@@ -1139 +1204 @@
 
 
-  err = _extrapolate_and_limit_from_gradient(ntri, beta,
-                                             (double*) domain_centroids -> data,
-                                             (long*)   domain_neighbours -> data,
-                                             (double*) quantity_centroid_values -> data,
-                                             (double*) domain_vertex_coordinates -> data,
-                                             (double*) quantity_vertex_values -> data,
-                                             (double*) quantity_edge_values -> data,
-                                             (double*) quantity_x_gradient -> data,
-                                             (double*) quantity_y_gradient -> data);
+  err = _extrapolate_from_gradient(ntri, 
+                                   (double*) domain_centroids -> data,
+                                   (double*) quantity_centroid_values -> data,
+                                   (double*) domain_vertex_coordinates -> data,
+                                   (double*) quantity_vertex_values -> data,
+                                   (double*) quantity_edge_values -> data,
+                                   (double*) quantity_x_gradient -> data,
+                                   (double*) quantity_y_gradient -> data);
 
   if (err != 0) {
       PyErr_SetString(PyExc_RuntimeError,
-                          "quantity_ext.c: Internal function _extrapolate_and_limit_from_gradient failed");
+                          "quantity_ext.c: Internal function _extrapolate_from_gradient failed");
+      return NULL;
+  }
+
+
+  err = _limit_edges_by_all_neighbours(ntri, beta,
+                                       (double*) quantity_centroid_values -> data,
+                                       (double*) quantity_vertex_values -> data,
+                                       (double*) quantity_edge_values -> data,
+                                       (long*)   domain_neighbours -> data,
+                                       (double*) quantity_x_gradient -> data,
+                                       (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _limit_edges_by_all_neighbours failed");
       return NULL;
   }
@@ -1165 +1244 @@
   Py_DECREF(quantity_vertex_values);
   Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_phi);
   Py_DECREF(quantity_x_gradient);
   Py_DECREF(quantity_y_gradient);
@@ -1445 +1525 @@
   //
 
-        PyObject *quantity, *domain, *Tmp;
+        PyObject *quantity, *domain;
         PyArrayObject
             *vertex_values,   //Conserved quantities at vertices
             *centroid_values, //Conserved quantities at centroids
             *edge_values,     //Conserved quantities at edges
+            *x_gradient,
+            *y_gradient,
             *neighbours;
 
@@ -1463 +1545 @@
         }
 
-        domain = PyObject_GetAttrString(quantity, "domain");
-        if (!domain) {
-          PyErr_SetString(PyExc_RuntimeError, 
-                          "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity");                        
-          
-          return NULL;
-        }
-
-        // Get safety factor beta_w
-        Tmp = PyObject_GetAttrString(domain, "beta_w");
-        if (!Tmp) {
-          PyErr_SetString(PyExc_RuntimeError, 
-                          "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain");                  
-          
-          return NULL;
-        }       
+        domain = get_python_object(quantity, "domain");
 
 
@@ -1486 +1553 @@
         vertex_values    = get_consecutive_array(quantity, "vertex_values");
         edge_values      = get_consecutive_array(quantity, "edge_values");
-        beta_w           = PyFloat_AsDouble(Tmp);
-
+        x_gradient       = get_consecutive_array(quantity, "x_gradient");
+        y_gradient       = get_consecutive_array(quantity, "y_gradient");
+        beta_w           = get_python_double(domain,"beta_w");
+
+        Py_DECREF(domain);
 
         N = centroid_values -> dimensions[0];
@@ -1495 +1565 @@
                                              (double*) vertex_values -> data,
                                              (double*) edge_values -> data,
-                                             (long*)   neighbours -> data);
+                                             (long*)   neighbours -> data,
+                                             (double*) x_gradient -> data,
+                                             (double*) y_gradient -> data);
 
         if (err != 0) {
           PyErr_SetString(PyExc_RuntimeError,
-                          "Internal function _limit_by_vertex failed");
+                          "quantity_ect.c: limit_edges_by_all_neighbours internal function _limit_edges_by_all_neighbours failed");
           return NULL;
         }       
@@ -1509 +1581 @@
         Py_DECREF(vertex_values);
         Py_DECREF(edge_values);
-        Py_DECREF(Tmp);
+
+
+
+        return Py_BuildValue("");
+}
+
+PyObject *bound_vertices_below_by_constant(PyObject *self, PyObject *args) {
+  //Bound a quantity below by a contant (useful for ensuring positivity
+  //precondition:
+  //  vertex values are already calulated, gradient consistent
+  //postcondition:
+  //  gradient, vertex and edge values are updated
+  //
+
+        PyObject *quantity, *domain;
+        PyArrayObject
+            *vertex_values,   //Conserved quantities at vertices
+            *centroid_values, //Conserved quantities at centroids
+            *edge_values,     //Conserved quantities at edges
+            *x_gradient,
+            *y_gradient;
+
+        double bound; //Safety factor
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "Od", &quantity, &bound)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: bound_vertices_below_by_constant could not parse input");
+          return NULL;
+        }
+
+        domain = get_python_object(quantity, "domain");
+
+        // Get pertinent variables
+        centroid_values  = get_consecutive_array(quantity, "centroid_values");
+        vertex_values    = get_consecutive_array(quantity, "vertex_values");
+        edge_values      = get_consecutive_array(quantity, "edge_values");
+        x_gradient       = get_consecutive_array(quantity, "x_gradient");
+        y_gradient       = get_consecutive_array(quantity, "y_gradient");
+
+
+
+        Py_DECREF(domain);
+
+        N = centroid_values -> dimensions[0];
+
+        err = _bound_vertices_below_by_constant(N, bound,
+                                             (double*) centroid_values -> data,
+                                             (double*) vertex_values -> data,
+                                             (double*) edge_values -> data,
+                                             (double*) x_gradient -> data,
+                                             (double*) y_gradient -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ect.c: bound_vertices_below_by_constant internal function _bound_vertices_below_by_constant failed");
+          return NULL;
+        }       
+
+
+        // Release
+        Py_DECREF(centroid_values);
+        Py_DECREF(vertex_values);
+        Py_DECREF(edge_values);
+        Py_DECREF(x_gradient);
+        Py_DECREF(y_gradient);
+
+
+
+        return Py_BuildValue("");
+}
+
+
+PyObject *bound_vertices_below_by_quantity(PyObject *self, PyObject *args) {
+  //Bound a quantity below by a contant (useful for ensuring positivity
+  //precondition:
+  //  vertex values are already calulated, gradient consistent
+  //postcondition:
+  //  gradient, vertex and edge values are updated
+  //
+
+        PyObject *quantity, *bounding_quantity, *domain;
+        PyArrayObject
+            *vertex_values,   //Conserved quantities at vertices
+            *centroid_values, //Conserved quantities at centroids
+            *edge_values,     //Conserved quantities at edges
+            *x_gradient,
+            *y_gradient,
+            *bound_vertex_values;
+
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "OO", &quantity, &bounding_quantity)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: bound_vertices_below_by_quantity could not parse input");
+          return NULL;
+        }
+
+        domain = get_python_object(quantity, "domain");
+
+        // Get pertinent variables
+        centroid_values     = get_consecutive_array(quantity, "centroid_values");
+        vertex_values       = get_consecutive_array(quantity, "vertex_values");
+        edge_values         = get_consecutive_array(quantity, "edge_values");
+        x_gradient          = get_consecutive_array(quantity, "x_gradient");
+        y_gradient          = get_consecutive_array(quantity, "y_gradient");
+        bound_vertex_values = get_consecutive_array(bounding_quantity, "vertex_values");
+
+
+
+        Py_DECREF(domain);
+
+        N = centroid_values -> dimensions[0];
+
+        err = _bound_vertices_below_by_quantity(N, 
+                                                (double*) bound_vertex_values -> data,
+                                                (double*) centroid_values -> data,
+                                                (double*) vertex_values -> data,
+                                                (double*) edge_values -> data,
+                                                (double*) x_gradient -> data,
+                                                (double*) y_gradient -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ect.c: bound_vertices_below_by_quantity internal function _bound_vertices_below_by_quantity failed");
+          return NULL;
+        }       
+
+
+        // Release
+        Py_DECREF(centroid_values);
+        Py_DECREF(vertex_values);
+        Py_DECREF(edge_values);
+        Py_DECREF(x_gradient);
+        Py_DECREF(y_gradient);
+        Py_DECREF(bound_vertex_values);
+
 
 
@@ -1589 +1801 @@
 
         // Release
+        Py_DECREF(domain);
         Py_DECREF(neighbours);
         Py_DECREF(centroid_values);
@@ -1700 +1913 @@
         {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
         {"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
+        {"bound_vertices_below_by_constant", bound_vertices_below_by_constant, METH_VARARGS, "Print out"},
+        {"bound_vertices_below_by_quantity", bound_vertices_below_by_quantity, METH_VARARGS, "Print out"},
         {"update", update, METH_VARARGS, "Print out"},
         {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},

anuga_core/source/anuga/advection/advection_ext.c

@@ -2 +2 @@
 //
 // To compile (Python2.X):
-// use python ../utilities/compile.py
-//
-// See the module advection.py
+// use python ../utilities/compile.py to
+// compile all C files within a directory
 //
 //
@@ -19 +18 @@
 
 
+//-------------------------------------------
+// Low level routines (called from wrappers)
+//------------------------------------------
 
 double _compute_fluxes(
@@ -38 +40 @@
 
  
-        //Loop
+        //Local Variables
 
         double qr,ql;
@@ -51 +53 @@
         int k3;
         
+        //Loop through triangles
+
         timestep = max_timestep;
 
@@ -59 +63 @@
             for (i=0; i<3; i++){
                 k_i = k3+i;
-                //Quantities inside volume facing neighbour i
+                //Quantities inside triangle facing neighbour i
                 ql = quantity_edge[k_i];
 
@@ -144 +148 @@
 
   PyObject *domain, *quantity;
+ 
   PyArrayObject 
     * quantity_update,

anuga_core/source/anuga/config.py

@@ -100 +100 @@
 #timestepping_method = 'rk2'   # 2nd Order TVD scheme
 
+# rk2 is a little more stable than euler, so rk2 timestepping
+# can deal with a larger beta when slope limiting the reconstructed
+# solution. The large beta is needed if solving problems sensitive
+# to numerical diffusion, like a small forced wave in an ocean
+beta_euler = 1.0
+beta_rk2   = 1.6
+
+
+
 # Option to search for signatures where isolated triangles are
 # responsible for a small global timestep.

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -103 +103 @@
 from anuga.config import alpha_balance
 from anuga.config import optimise_dry_cells
+from anuga.config import optimised_gradient_limiter
 
 #---------------------
@@ -176 +177 @@
         self.minimum_storable_height = minimum_storable_height
         self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
-        
+
+        # Limiters
+        self.use_old_limiter = True
+
+        self.optimised_gradient_limiter = optimised_gradient_limiter
                 
 
@@ -187 +192 @@
         self.beta_w      = beta
         self.beta_w_dry  = beta
+        self.quantities['stage'].beta = beta
+        
         self.beta_uh     = beta
         self.beta_uh_dry = beta
+        self.quantities['xmomentum'].beta = beta
+        
         self.beta_vh     = beta
         self.beta_vh_dry = beta
+        self.quantities['ymomentum'].beta = beta
+        
         self.beta_h      = beta
         
@@ -384 +395 @@
         # Call correct module function
         # (either from this module or C-extension)
-        distribute_to_vertices_and_edges(self)
+        if self.use_old_limiter:
+            distribute_to_vertices_and_edges(self)
+        else:
+            for name in self.conserved_quantities:
+                Q = self.quantities[name]
+                if self._order_ == 1:
+                    Q.extrapolate_first_order()
+                elif self._order_ == 2:
+                    Q.extrapolate_second_order_and_limit()
+                    if name == 'stage':
+                        Q.bound_vertices_below_by_quantity(self.quantities['elevation'])
+                else:
+                    raise 'Unknown order'
+ 
 
 
@@ -399 +423 @@
         # self.check_integrity()
 
-        msg = 'Parameter beta_h must be in the interval [0, 1['
-        assert 0 <= self.beta_h <= 1.0, msg
-        msg = 'Parameter beta_w must be in the interval [0, 1['
-        assert 0 <= self.beta_w <= 1.0, msg
+        msg = 'Parameter beta_h must be in the interval [0, 2['
+        assert 0 <= self.beta_h <= 2.0, msg
+        msg = 'Parameter beta_w must be in the interval [0, 2['
+        assert 0 <= self.beta_w <= 2.0, msg
 
 
@@ -737 +761 @@
     """
 
-    from anuga.config import optimised_gradient_limiter
+    
 
     # Remove very thin layers of water
@@ -743 +767 @@
 
     # Extrapolate all conserved quantities
-    if optimised_gradient_limiter:
+    if domain.optimised_gradient_limiter:
         # MH090605 if second order,
         # perform the extrapolation and limiting on
@@ -764 +788 @@
                 Q.extrapolate_first_order()
             elif domain._order_ == 2:
-
-                # Experiment
-                #if name == 'stage':
-                #    #print name, 'second'
-                #    Q.extrapolate_second_order()
-                #    Q.limit()
-                #else:
-                #    #print name, 'first'                
-                #    Q.extrapolate_first_order()
-                #    #Q.extrapolate_second_order()
-                #    #Q.limit()                
-                
                 Q.extrapolate_second_order()
                 Q.limit()

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1088 +1088 @@
   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
   double dqv[3], qmin, qmax, hmin, hmax;
-  double hc, h0, h1, h2, beta_tmp;
+  double hc, h0, h1, h2, beta_tmp, hfactor;
 
         
@@ -1185 +1185 @@
       h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
       h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-      hmin = min(h0,min(h1,h2));
-      
+      hmin = min(min(h0,min(h1,h2)),hc);
+      //hfactor = hc/(hc + 1.0);
+
+      hfactor = 0.0;
+      if (hmin > 0.1 ) {
+          hfactor = (hmin-0.1)/(hmin+0.4);
+      }
       
       if (optimise_dry_cells) {      
@@ -1230 +1235 @@
       
       // Playing with dry wet interface
-      hmin = qmin;
-      beta_tmp = beta_w;
-      if (hmin<minimum_allowed_height)
-        beta_tmp = beta_w_dry;
+      //hmin = qmin;
+      //beta_tmp = beta_w_dry;
+      //if (hmin>minimum_allowed_height)
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
         
+      //printf("min_alled_height = %f\n",minimum_allowed_height);
+      //printf("hmin = %f\n",hmin);
+      //printf("beta_w = %f\n",beta_w);
+      //printf("beta_tmp = %f\n",beta_tmp);
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp); 
@@ -1271 +1280 @@
       // from the centroid to the min and max
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      beta_tmp = beta_uh;
-      if (hmin<minimum_allowed_height)
-        beta_tmp = beta_uh_dry;
-        
+      //beta_tmp = beta_uh;
+      //if (hmin<minimum_allowed_height)
+      //beta_tmp = beta_uh_dry;
+      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp);
@@ -1311 +1321 @@
       // from the centroid to the min and max
       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
-      beta_tmp = beta_vh;
-      
-      if (hmin<minimum_allowed_height)
-        beta_tmp = beta_vh_dry;
-        
+      
+      //beta_tmp = beta_vh;
+      //
+      //if (hmin<minimum_allowed_height)
+      //beta_tmp = beta_vh_dry;
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
+
       // Limit the gradient
       limit_gradient(dqv,qmin,qmax,beta_tmp);
@@ -1563 +1575 @@
   
 
-  beta_w      = get_double(domain,"beta_w");
-  beta_w_dry  = get_double(domain,"beta_w_dry");
-  beta_uh     = get_double(domain,"beta_uh");
-  beta_uh_dry = get_double(domain,"beta_uh_dry");
-  beta_vh     = get_double(domain,"beta_vh");
-  beta_vh_dry = get_double(domain,"beta_vh_dry");  
-
-  minimum_allowed_height = get_double(domain,"minimum_allowed_height");
-  epsilon = get_double(domain,"epsilon");
+  beta_w                 = get_python_double(domain,"beta_w");
+  beta_w_dry             = get_python_double(domain,"beta_w_dry");
+  beta_uh                = get_python_double(domain,"beta_uh");
+  beta_uh_dry            = get_python_double(domain,"beta_uh_dry");
+  beta_vh                = get_python_double(domain,"beta_vh");
+  beta_vh_dry            = get_python_double(domain,"beta_vh_dry");  
+
+  minimum_allowed_height = get_python_double(domain,"minimum_allowed_height");
+  epsilon                = get_python_double(domain,"epsilon");
 
   number_of_elements = stage_centroid_values -> dimensions[0];  

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -248 +248 @@
         
         range = fid.variables['xmomentum_range'][:]
-        assert allclose(range,[0,0.4695096]) or\
-               allclose(range,[0,0.47790655]) # Old slope limiters
+
+
+
+        assert allclose(range,[0,0.4695096]) or \
+               allclose(range,[0,0.47790655]) or\
+               allclose(range,[0,0.46941957]) 
         
         range = fid.variables['ymomentum_range'][:]
         #assert allclose(range,[0,0.02174380])
+        
         assert allclose(range,[0,0.02174439]) or\
-               allclose(range,[0,0.02283983]) # Old slope limiters 
+               allclose(range,[0,0.02283983]) or\
+               allclose(range,[0,0.0217342]) # Old slope limiters 
         
         fid.close()
@@ -324 +330 @@
 
         extrema = fid.variables['xmomentum.extrema'][:]
-        assert allclose(extrema,[-0.06062178, 0.47886313])
+        assert allclose(extrema,[-0.06062178, 0.47873023])
         
         extrema = fid.variables['ymomentum.extrema'][:]
-        assert allclose(extrema,[0.00, 0.06241221])        
+        assert allclose(extrema,[0.00, 0.0625786])        
 
         time_interval = fid.variables['extrema.time_interval'][:]

anuga_core/source/anuga/utilities/util_ext.h

@@ -261 +261 @@
     
   B = (PyArrayObject*) PyObject_GetAttrString(O, name);
-  if (!B) return NULL;     
+  if (!B) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }     
   
   //Convert to consecutive array
@@ -269 +272 @@
   Py_DECREF(B); //FIXME: Is this really needed??
   
-  if (!A) return NULL;
+  if (!A) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
+    return NULL;
+  }
   return A;
 }
 
-double get_double(PyObject *O, char *name) {
+double get_python_double(PyObject *O, char *name) {
   PyObject *TObject;
   double tmp;
@@ -281 +287 @@
   TObject = PyObject_GetAttrString(O, name);
   if (!TObject) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_double could not obtain double from object");
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
     return 0.0;
   }  
@@ -293 +299 @@
 
 
-
+PyObject *get_python_object(PyObject *O, char *name) {
+  PyObject *Oout;
+
+  Oout = PyObject_GetAttrString(O, name);
+  if (!Oout) {
+    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object");
+    return NULL;
+  }
+
+  return Oout;
+
+}
+
+
+

anuga_core/source/anuga/utilities/xml_tools.py

@@ -239 +239 @@
         fid = xml
 
-    try:
-        dom = parse(fid)
-    except Exception, e:
-        # Throw filename into dom exception
-        msg = 'XML file "%s" could not be parsed: ' %fid.name
-        msg += str(e)
-        raise Exception, msg
+    dom = parse(fid)
+
+##     try:
+##         dom = parse(fid)
+##     except Exception, e:
+##         # Throw filename into dom exception
+##         msg = 'XML file "%s" could not be parsed: ' %fid.name
+##         msg += str(e)
+##         raise Exception, msg
 
     return dom2object(dom)

anuga_validation/convergence_study/animatesww2d_alt.py

@@ -1 @@
-def animatesww2d_alt(sww_filename, movie_filename, range):
-    """Plot cross section of model output
-    """
+## def animatesww2d_alt(sww_filename, movie_filename, range):
+##     """Plot cross section of model output
+##     """
 
-    # Read SWW file     
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(sww_filename, 'r')
+##     # Read SWW file  
+##     from Scientific.IO.NetCDF import NetCDFFile
+##     fid = NetCDFFile(sww_filename, 'r')
     
-    x = fid.variables['x']
-    y = fid.variables['y']
-    volumes = fid.variables['volumes']
+##     x = fid.variables['x']
+##     y = fid.variables['y']
+##     volumes = fid.variables['volumes']
     
-    elevation = fid.variables['elevation']
-    time = fid.variables['time']
-    stage = fid.variables['stage']
-    xmomentum = fid.variables['xmomentum']
-    ymomentum = fid.variables['ymomentum']
+##     elevation = fid.variables['elevation']
+##     time = fid.variables['time']
+##     stage = fid.variables['stage']
+##     xmomentum = fid.variables['xmomentum']
+##     ymomentum = fid.variables['ymomentum']
     
 def animatesww2d(swwfile):
@@ -34 +34 @@
     # interpolation points are for y = 0
     # number of increments in x
-    m = 100 
+    m = 1000 
     max_x = 100000.
     x = 0

anuga_validation/convergence_study/convergence_structured.py

@@ -17 +17 @@
 from anuga.shallow_water import Transmissive_boundary
 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from anuga.geospatial_data.geospatial_data import *
+#from anuga.geospatial_data.geospatial_data import *
 from math import cos
+from Numeric import zeros, Float
 
 #------------------------------------------------------------------------------
 # Setup computational domain
 #------------------------------------------------------------------------------
-dx = 1000.
+dx = 200.
 dy = dx
 L = 100000.
@@ -33 +34 @@
 domain = Domain(points, vertices, boundary) 
 
-domain.set_timestepping_method('euler')
+domain.set_timestepping_method('rk2')
 domain.set_default_order(2)
 domain.set_name('myexample9')                
 domain.set_datadir('.')                     # Use current directory for output
 
-domain.beta_w      = 1.0
-domain.beta_w_dry  = 0.2
-domain.beta_uh     = 1.0
-domain.beta_uh_dry = 0.2
-domain.beta_vh     = 1.0
-domain.beta_vh_dry = 0.2
-domain.beta_h      = 1.0
+domain.set_all_limiters(0.9)
+
+print domain.beta_w
+domain.use_old_limiter = False
+domain.CFL = 1.0
 
 #------------------------------------------------------------------------------
@@ -50 +49 @@
 #------------------------------------------------------------------------------
 #domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('elevation',-100)
+domain.set_quantity('elevation',0.0)
 domain.set_quantity('friction', 0.00)
-domain.set_quantity('stage', 0.0)            
+
+h0 = 10.0
+h1 = 1.0
+
+def height(x,y):
+    z = zeros(len(x),Float)
+    for i in range(len(x)):
+        if x[i]<=50000.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+
+
+domain.set_quantity('stage', height)
+#domain.set_quantity('stage', 0.0)            
 
 #-----------------------------------------------------------------------------
@@ -65 +79 @@
 Bw = Time_boundary(domain=domain,     # Time dependent boundary  
 ## Sine wave
-                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+#                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
+## Single wave
+                   f=lambda t: [h0, 0.0, 0.0])
 ## Sawtooth?
 #                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
@@ -92 +108 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 20.0, finaltime = 10*40*60.):
+for t in domain.evolve(yieldstep = 50.0, finaltime = 10*40*60.):
     domain.write_time()
     vis.update()

anuga_work/development/analytical_solutions/oblique_shock.py

@@ -29 +29 @@
 
 
-from shallow_water import Domain, Constant_height
-from shallow_water import Transmissive_boundary, Reflective_boundary,\
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
      Dirichlet_boundary
-
+from anuga.visualiser import RealtimeVisualiser
 from math import sqrt, cos, sin, pi
-from mesh_factory import oblique
+from mesh_factory import oblique_cross
 
 
@@ -44 +44 @@
 leny = 30.
 lenx = 40.
-n = 50
-m = 60
+n = 100
+m = 120
 
-points, elements, boundary = oblique(m, n, lenx, leny)
+points, elements, boundary = oblique_cross(m, n, lenx, leny)
 domain = Domain(points, elements, boundary)
 
@@ -54 +54 @@
 
 # Store output
-domain.store=True
+#domain.store=True
 
 # Output format
-domain.format="sww" #NET.CDF binary format
+#domain.format="sww" #NET.CDF binary format
                     # "dat" for ASCII
 
@@ -87 +87 @@
 #Initial condition
 h = 0.5
-domain.set_quantity('stage', Constant_height(x_slope, h) )
+domain.set_quantity('stage', expression = 'elevation + %d'%h )
 
+#---------------------------------
+# Setup visualization
+#---------------------------------
+vis = RealtimeVisualiser(domain)
+vis.render_quantity_height("elevation", dynamic=False)
+vis.render_quantity_height("stage", zScale=10, dynamic=True)
+vis.colour_height_quantity('stage', (0.75, 0.5, 0.5))
+vis.start()
 
 ######################
@@ -96 +104 @@
 for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
     domain.write_time()
+    vis.update()
 
 print 'That took %.2f seconds' %(time.time()-t0)
+
+
+vis.evolveFinished()
+vis.join()
 
 #FIXME: Compute average water depth on either side of shock and compare

anuga_work/development/near_shore_PMD/run_beach.py

@@ -167 +167 @@
     domain.set_boundary( {'wall': Br, 'wave': Bwp_velocity} )
 
+
     #-------------------------------------------------------------------------
     # Evolve system through time
@@ -174 +175 @@
     for t in domain.evolve(yieldstep, finaltime):    
         domain.write_time()
+
     print 'That took %.2f seconds' %(time.time()-t0)
     print 'finished'