Changeset 8274 for trunk

Timestamp:

Dec 8, 2011, 11:38:18 AM (13 years ago)

Author:

paul

Message:

Fixed loss of mass in sqpipe and rearranged code a little

Location:

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe

Files:

• parabolic_canal.py (modified) (3 diffs)
• sqpipe.h (modified) (14 diffs)
• sqpipe_comp_flux.c (modified) (2 diffs)
• sqpipe_default.h (modified) (10 diffs)
• sqpipe_domain.py (modified) (1 diff)
• sqpipe_python.h (modified) (2 diffs)
• test_sqpipe_domain.py (modified) (1 diff)

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py

@@ -94 +94 @@
     pylab.ion()
 
-    x = domain.get_centroids()
-    z = domain.get_quantity('elevation', 'centroids')
-    w = domain.get_quantity('stage', 'centroids')
-    h = domain.get_quantity('height', 'centroids')
-    v = domain.get_quantity('velocity', 'centroids')
-    t = domain.get_quantity('top', 'centroids')
-    s = domain.state
-    m = domain.get_mass()
-    M = m.sum() * numpy.ones_like(x)
+    x, z, w, h, v, t, s, m, M = get_quantities(domain)
 
-    plot1 = pylab.subplot(511)
-    zplot, = pylab.plot(x, z)
-    wplot, = pylab.plot(x, w)
-    ztplot, = pylab.plot(x, z+t)    
-    plot1.set_ylim([-1,50])
-    pylab.xlabel('Position')
-    pylab.ylabel('Stage')
-
-    plot2 = pylab.subplot(512)
-    hplot, = pylab.plot(x, h)
-    tplot, = pylab.plot(x, t)
-    plot2.set_ylim([-1,20])
-    pylab.xlabel('Position')
-    pylab.ylabel('Height')
-
-    plot3 = pylab.subplot(513)
-    vplot, = pylab.plot(x, v)
-    plot3.set_ylim([-30,30])
-    pylab.xlabel('Position')
-    pylab.ylabel('Velocity')
-
-    plot4 = pylab.subplot(514)
-    splot, = pylab.plot(x, s)
-    plot4.set_ylim([-1,2])
-    pylab.xlabel('Position')
-    pylab.ylabel('State')
-
-    plot5 = pylab.subplot(515)
-    Mplot, = pylab.plot(x, m)
-    mplot, = pylab.plot(x, M)
-    plot5.set_ylim([-1,45000])
-    pylab.xlabel('Position')
-    pylab.ylabel('Mass')
-
+    zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot = make_plots(x, z, w, h, v, t, s, m, M)
 
     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
-        z = domain.get_quantity('elevation', 'centroids')
-        w = domain.get_quantity('stage', 'centroids')
-        h = domain.get_quantity('height', 'centroids')
-        t = domain.get_quantity('top', 'centroids')
-        v = domain.get_quantity('velocity', 'centroids')
-        s = domain.state
-        m = domain.get_mass()
-        M = m.sum() * numpy.ones_like(x)
+        x, z, w, h, v, t, s, m, M = get_quantities(domain)
         
         zplot.set_ydata(z)
@@ -166 +118 @@
     import pylab
 
-    x = domain.get_centroids()
-    z = domain.get_quantity('elevation', 'centroids')
-    w = domain.get_quantity('stage', 'centroids')
-    h = domain.get_quantity('height', 'centroids')
-    t = domain.get_quantity('top', 'centroids')
-    u = domain.get_quantity('velocity', 'centroids')
+    x, z, w, h, v, t, s, m, M = get_quantities(domain)
 
     pylab.ioff()
     pylab.hold(False)
 
-    plot1 = pylab.subplot(311)
-    pylab.plot(x, z, x, w, x, z+t)
-    plot1.set_ylim([-1,40])
-    pylab.xlabel('Position')
-    pylab.ylabel('Stage')
-
-    plot2 = pylab.subplot(312)
-    pylab.plot(x, h, x, t)
-    plot2.set_ylim([-1,10])
-    pylab.xlabel('Position')
-    pylab.ylabel('Height')
-
-    plot3 = pylab.subplot(313)
-    pylab.plot(x, u)
-    plot3.set_ylim([-15,15])
-    pylab.xlabel('Position')
-    pylab.ylabel('Velocity')
-
+    make_plots(x, z, w, h, v, t, s, m, M)
+    
     pylab.show()
 
@@ -205 +136 @@
         f.write("%s %s %s\n" % (x[i], z[i], w[i]))
     f.close()
+
+def get_quantities(domain):
+    x = domain.get_centroids()
+    z = domain.get_quantity('elevation', 'centroids')
+    w = domain.get_quantity('stage', 'centroids')
+    h = domain.get_quantity('height', 'centroids')
+    v = domain.get_quantity('velocity', 'centroids')
+    t = domain.get_quantity('top', 'centroids')
+    s = domain.state
+    m = domain.get_mass()
+    M = m.sum() * numpy.ones_like(x)
+
+    return x, z, w, h, v, t, s, m, M
+
+def make_plots(x, z, w, h, v, t, s, m, M):
+    import pylab
+    
+    plot1 = pylab.subplot(321)
+    zplot, = pylab.plot(x, z)
+    wplot, = pylab.plot(x, w)
+    ztplot, = pylab.plot(x, z+t)    
+    plot1.set_ylim([-1,50])
+    pylab.xlabel('Position')
+    pylab.ylabel('Stage')
+
+    plot2 = pylab.subplot(322)
+    hplot, = pylab.plot(x, h)
+    tplot, = pylab.plot(x, t)
+    plot2.set_ylim([-1,20])
+    pylab.xlabel('Position')
+    pylab.ylabel('Height')
+
+    plot3 = pylab.subplot(323)
+    vplot, = pylab.plot(x, v)
+    plot3.set_ylim([-30,30])
+    pylab.xlabel('Position')
+    pylab.ylabel('Velocity')
+
+    plot4 = pylab.subplot(324)
+    splot, = pylab.plot(x, s)
+    plot4.set_ylim([-1,2])
+    pylab.xlabel('Position')
+    pylab.ylabel('State')
+
+    plot5 = pylab.subplot(325)
+    mplot, = pylab.plot(x, m)
+    plot5.set_ylim([-1,1000])
+    pylab.xlabel('Position')
+    pylab.ylabel('Mass')
+
+    plot6 = pylab.subplot(326) Mplot, = pylab.plot(x, M)
+    plot6.set_ylim([-1,45000]) pylab.xlabel('Position')
+    pylab.ylabel('Total Mass')
+
+    return zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h

@@ -69 +69 @@
   double b;
   double t;
+  long state;
 };
 
@@ -79 +80 @@
 
 // Methods
+double* sqpipe_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
+double sqpipe_quantity_sound_speed (struct quantity *q, struct params *p);
+double sqpipe_quantity_get_conserved (struct quantity *q, int k, double normal);
+
+// Helper methods
 double* sqpipe_quantity_free_surface_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
 double* sqpipe_quantity_pressurised_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
@@ -85 +91 @@
 
 // Implementation
+double* sqpipe_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
+  if (((struct sqpipe_quantity *)q)->state == 0) {
+    return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
+  } else {
+    return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
+  }
+}
+
+double sqpipe_quantity_sound_speed (struct quantity *q, struct params *p) {
+  if (((struct sqpipe_quantity *)q)->state == 0) {
+    return sqpipe_quantity_free_surface_sound_speed(q, p);
+  } else {
+    return sqpipe_quantity_pressurised_sound_speed(q, p);
+  }
+}
+
+double sqpipe_quantity_get_conserved (struct quantity *q, int k, double normal) {
+  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
+  double c;
+  
+  switch (k) {
+  case 0:
+    c = p->a;
+    break;
+  case 1:
+    // This should be normal^2 p->d and normal is +/- 1
+    c = p->d;
+    break;
+  default:
+    c = 0;
+  }
+  return c;
+}
+
+
 // Flux formula taking into account the normal
 // Note u, d should be normal*u, normal*d but since normal is +/- 1
@@ -121 +162 @@
 void sqpipe_quantity_init(struct sqpipe_quantity *q) {
   static struct quantity_vtable vtable = {
-    NULL,
-    NULL,
-    NULL,
+    &sqpipe_quantity_flux_formula,
+    &sqpipe_quantity_sound_speed,
+    &sqpipe_quantity_get_conserved,
     &sqpipe_quantity_destroy
   };
@@ -136 +177 @@
   q->b = 0.0;
   q->t = 0.0;
+  q->state = 0;
 }
 
@@ -237 +279 @@
 struct sqpipe_cell {
   struct cell super;
+  
+  long state;
 };
 
@@ -248 +292 @@
 // Methods
 double sqpipe_cell_extreme_sound_speeds (struct cell *c, struct params *p);
+double* sqpipe_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
+
+// Helper methods
 double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce);
 double* sqpipe_cell_pressurised_forcing_terms (struct cell *c, struct params *p, double *cellforce);
@@ -270 +317 @@
 }
 
+double* sqpipe_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
+  if (((struct sqpipe_cell *)c)->state == 0) {
+    return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);    
+  } else {
+    return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
+  }
+}
+
 double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
   struct sqpipe_quantity *qin = (struct sqpipe_quantity *) c->edges[0]->qout;
@@ -297 +352 @@
     &cell_flux_function_generic, 
     &sqpipe_cell_extreme_sound_speeds,
-    NULL,
+        &sqpipe_cell_forcing_terms,
     &sqpipe_cell_destroy
   };
@@ -309 +364 @@
 
   c->super.num_edges = 2;
+  ((struct sqpipe_cell *)c)->state = 0;
 }
 
@@ -365 +421 @@
 // Methods
 struct quantity* sqpipe_quantities_get_quantity(struct quantities *qs, int i, struct quantity *q);
+void sqpipe_quantities_update (struct quantities *qs, double *flux, int k);
 
 // Implementation
@@ -383 +440 @@
 }
 
+void sqpipe_quantities_update (struct quantities *qs, double *flux, int k) {
+  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
+  ps->a[k] = flux[0];
+  ps->d[k] = flux[1];
+}
+
+
 void sqpipe_quantities_init(struct sqpipe_quantities *q) {
   static struct quantities_vtable vtable = {
     &sqpipe_quantities_get_quantity,
-    NULL,
+    &sqpipe_quantities_update,
     &sqpipe_quantities_destroy
   };
@@ -426 +490 @@
 struct sqpipe_domain {
   struct domain super;
+  long *state;
 };
 
 // Construction
-void sqpipe_domain_init(struct sqpipe_domain *q);
+void sqpipe_domain_init(struct sqpipe_domain *d);
 struct sqpipe_domain* sqpipe_domain_new();
 
 // Destruction
-struct domain* sqpipe_domain_destroy(struct domain *D);
+struct domain* sqpipe_domain_destroy(struct domain *d);
 
 // Methods
-void sqpipe_domain_init(struct sqpipe_domain *D);
+struct cell* sqpipe_domain_get_cell(struct domain *D, int k, struct cell *c);
 
 // Implementation
+struct cell* sqpipe_domain_get_cell(struct domain *D, int k, struct cell *c) {
+  struct sqpipe_cell *sc = (struct sqpipe_cell *) c;
+  struct sqpipe_domain *sD = (struct sqpipe_domain *) D;
+  int i, ki, n;
+  long state_in, state_out, state;
+  
+  // First get a generic cell
+  c = domain_get_cell_generic(D, k, c);
+
+  // Set state
+  // The state should be set to state_in for qin and state_out for qout
+  // Temporarily, I set the state of the edge, aribtrarily choosing
+  // pressurised for transition points
+  // This will be fixed when the flux is computed at the edge rather than
+  // at the quantity
+  sc->state = sD->state[k];
+  for (i = 0; i<c->num_edges; i++) {
+        // The inner state is that of this cell
+        state_in = sc->state;
+        // The outer state is that of the corresponding neighbour cell
+        ki = k*(c->num_edges) + i;
+        n = D->neighbours[ki];
+        state_out = sD->state[n];
+        
+        // If both in and out state are the same, choose that state
+        // Otherwise, choose pressurised
+        if (state_in == state_out) {
+                state = state_in;
+        } else {
+                state = 1;
+        }
+    ((struct sqpipe_quantity *)(c->edges[i]->qin))->state = state;      
+    ((struct sqpipe_quantity *)(c->edges[i]->qout))->state = state;
+  }
+  
+  return c;
+}
+
+
 void sqpipe_domain_init(struct sqpipe_domain *d) {
   static struct domain_vtable vtable = {
     NULL,
     NULL,
-    &domain_compute_fluxes_generic,
-    &domain_get_cell_generic,
+    &domain_compute_fluxes_cells_generic,
+    &sqpipe_domain_get_cell,
     &sqpipe_domain_destroy
   };
@@ -463 +567 @@
   ((struct domain *)d)->max_speed_array = NULL;
 
+  d->state = NULL;
   // Can't set params to NULL here since sqpipe_domain_new() has allocated 
   // memory

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_comp_flux.c

@@ -25 +25 @@
   struct domain *D = (struct domain *) sqpipe_default_domain_new();
   double timestep;
-  PyArrayObject *state;
   PyObject *domain;
 
@@ -34 +33 @@
   }
 
+  // Generic sqpipe domain
   D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
 
   // We need the state 
-  state = get_consecutive_array(domain, "state");
-  ((struct sqpipe_default_domain *)D)->state = (long *) state->data;
+  //state = get_consecutive_array(domain, "state");
+  //((struct sqpipe_default_domain *)D)->state = (long *) state->data;
 
+  // Compute the fluxes and return the new timestep
   timestep = domain_compute_fluxes(D);
   domain_destroy(D);

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_default.h

@@ -7 +7 @@
 struct sqpipe_default_quantity {
   struct sqpipe_quantity super;
-  long state;
 };
 
@@ -17 +16 @@
 
 // Methods
-double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
-double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p);
-double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal);
 
 // Implementation
-double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
-  if (((struct sqpipe_default_quantity *)q)->state == 0) {
-    return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
-  } else {
-    return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
-  }
-}
-
-double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p) {
-  if (((struct sqpipe_default_quantity *)q)->state == 0) {
-    return sqpipe_quantity_free_surface_sound_speed(q, p);
-  } else {
-    return sqpipe_quantity_pressurised_sound_speed(q, p);
-  }
-}
-
-double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal) {
-  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
-  double c;
-  
-  switch (k) {
-  case 0:
-    c = p->a;
-    break;
-  case 1:
-    // This should be normal^2 p->d and normal is +/- 1
-    c = p->d;
-    break;
-  default:
-    c = 0;
-  }
-  return c;
-}
-
 void sqpipe_default_quantity_init (struct sqpipe_default_quantity *q) {
   struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
 
-  sqpipe_quantity_init(p);
-
-  static struct quantity_vtable vtable = {
-    &sqpipe_default_quantity_flux_formula,
-    &sqpipe_default_quantity_sound_speed,
-    &sqpipe_default_quantity_get_conserved,
-    &sqpipe_quantity_destroy
-  };
-  p->super.vtable = &vtable;
-
-  ((struct sqpipe_default_quantity *)p)->state = 0;
+  sqpipe_quantity_init(p);  
 }
 
@@ -98 +50 @@
 
 // Implementation
-
 void sqpipe_default_edge_init (struct sqpipe_default_edge *q) {
   sqpipe_edge_init((struct sqpipe_edge *)q);
@@ -122 +73 @@
 struct sqpipe_default_cell {
   struct sqpipe_cell super;
-
-  long state;
 };
 
@@ -133 +82 @@
 
 // Methods
-double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
 
 // Implementation
-double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
-  if (((struct sqpipe_default_cell *)c)->state == 0) {
-    return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);    
-  } else {
-    return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
-  }
-}
-
-void sqpipe_default_cell_init (struct sqpipe_default_cell *q) {
-  struct sqpipe_cell *p = (struct sqpipe_cell *) q;
+void sqpipe_default_cell_init (struct sqpipe_default_cell *c) {
+  struct sqpipe_cell *p = (struct sqpipe_cell *) c;
 
   sqpipe_cell_init(p);
-
-  static struct cell_vtable vtable = {
-    &cell_flux_function_generic,
-    &sqpipe_cell_extreme_sound_speeds,
-    &sqpipe_default_cell_forcing_terms,
-    &sqpipe_cell_destroy
-  };
-  p->super.vtable = &vtable;
-
-  ((struct sqpipe_default_cell *)p)->state = 0;
 }
 
@@ -193 +123 @@
 
 // Methods
-void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k);
-
 
 // Implementation
-void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k) {
-  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
-  ps->a[k] = flux[0];
-  ps->d[k] = flux[1];
-}
-
 void sqpipe_default_quantities_init (struct sqpipe_default_quantities *q) {
   struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
 
   sqpipe_quantities_init(p);
-
-  static struct quantities_vtable vtable = {
-    &sqpipe_quantities_get_quantity,
-    &sqpipe_default_quantities_update,
-    &sqpipe_quantities_destroy
-  };
-  p->super.vtable = &vtable;
 }
 
@@ -232 +147 @@
 struct sqpipe_default_domain {
   struct sqpipe_domain super;
-  long *state;
 };
 
@@ -244 +158 @@
 struct edge* sqpipe_default_domain_new_edge(struct domain *d);
 struct cell* sqpipe_default_domain_new_cell(struct domain *d);
-struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c);
 
 // Implementation
@@ -258 +171 @@
 }
 
-struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c) {
-  struct sqpipe_default_cell *sc = (struct sqpipe_default_cell *) c;
-  struct sqpipe_default_domain *sD = (struct sqpipe_default_domain *) D;
-  int i;
-  
-  c = domain_get_cell_generic(D, k, c);
-
-  // Set state
-  sc->state = sD->state[k];
-  for (i = 0; i<c->num_edges; i++) {
-    ((struct sqpipe_default_quantity *)(c->edges[i]->qin))->state = sc->state;
-    ((struct sqpipe_default_quantity *)(c->edges[i]->qout))->state = sc->state;
-  }
-  
-  return c;
-}
-
 void sqpipe_default_domain_init (struct sqpipe_default_domain *d) {
   struct sqpipe_domain *p = (struct sqpipe_domain *)d;
@@ -283 +179 @@
     &sqpipe_default_domain_new_edge,
     &sqpipe_default_domain_new_cell,
-    &domain_compute_fluxes_generic,
-    &sqpipe_default_domain_get_cell,
+    &domain_compute_fluxes_cells_generic,
+    &sqpipe_domain_get_cell,
     &sqpipe_domain_destroy
   };
   p->super.vtable = &vtable;
-
-  ((struct sqpipe_default_domain *)p)->state = NULL;
 }
 

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

@@ -270 +270 @@
         """
 
-        # Equivalent area A is defined by A = (\rho/\rho_0) A_{pipe}
-        # A_{pipe} = width * top
-        # A = width * height
-        # \rho_0 = 1
-        # Therefore \rho = height/top for pressurised states
-        # For free surface \rho = 1
-
-        # Coarse approximation of mass M in a cell, assume rho is constant
-        # M = \Delta x * \rho * A
-        #   = \Delta_x * A * height/top (pressurised)
-        #   = \Delta_x * A (un-pressurised)
-
-        area      = self.quantities['area']       
-        height    = self.quantities['height']        
-        top       = self.quantities['top']
-        
-        #Arrays   
-        a   = area.centroid_values        
-        h   = height.centroid_values
-        t   = top.centroid_values
-        
-        M = numpy.where(self.state == 0, self.areas * a, self.areas * a * h/t)
-
-        return M
+        # The equivalent area is the conserved mass quantity
+        # It is equal to \rho/\rho_0 A where A is the cross sectional
+        # area of the fluid where $\rho_0 is some fixed reference density
+        area = self.quantities['area']
+        a = area.centroid_values
+        
+        return a * self.areas
+        
 
 # Auxillary methods

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_python.h

@@ -4 +4 @@
 struct sqpipe_domain* sqpipe_domain_python_get(struct sqpipe_domain *sD, PyObject *domain, double timestep) {
   PyObject *quantities;
+  PyArrayObject *state;
   struct domain *D = (struct domain *)sD;
   struct sqpipe_params *sp = (struct sqpipe_params *) D->params;
-
+  
   // Get a generic domain
   D = get_python_domain(D, domain, timestep);
@@ -15 +16 @@
   D->boundary_values = sqpipe_quantities_python_get(D->boundary_values, quantities, "boundary_values");
   D->explicit_update = sqpipe_quantities_python_get(D->explicit_update, quantities, "explicit_update");
+
+  // Get cell level quantities
+  state = get_consecutive_array(domain, "state");
+  sD->state = (long *) state->data;
 
   // Get the specific parameters for square pipes

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py

@@ -5 +5 @@
 
 domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
-finaltime = 10000.0
+finaltime = 1000.0
 yieldstep = 10.0
+#m0 = (domain.get_mass()).sum()
+#for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+#       m = (domain.get_mass()).sum()
+#       if abs(m - m0) > 0.5:
+#               print m0 - m
+#               
+#       m0 = m
+        
+#outfile="out"
+#anuga_1d.sqpipe.parabolic_canal.write_domain(domain, outfile)
+
 anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
+print "finished"
 anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)