Changeset 8455

Timestamp:

Jul 5, 2012, 5:17:34 PM (13 years ago)

Author:

steve

Message:

Added in a Set_stage_operator

Location:

trunk/anuga_core/source

Files:

• anuga/__init__.py (modified) (1 diff)
• anuga/abstract_2d_finite_volumes/generic_domain.py (modified) (4 diffs)
• anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (5 diffs)
• anuga/abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
• anuga/file/test_sww.py (modified) (5 diffs)
• anuga/operators/run_set_stage.py (modified) (1 diff)
• anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• anuga/shallow_water/shallow_water_ext.c (modified) (4 diffs)
• anuga_parallel/parallel_api.py (modified) (3 diffs)
• anuga_parallel/parallel_generic_communications.py (modified) (1 diff)
• anuga_parallel/parallel_shallow_water.py (modified) (2 diffs)
• anuga_parallel/run_parallel_sw_merimbula.py (modified) (4 diffs)

trunk/anuga_core/source/anuga/__init__.py

@@ -150 +150 @@
 
 
-
+indent = '   '
 
 #-----------------------------

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -113 +113 @@
 
         self.number_of_boundaries = self.mesh.number_of_boundaries
+        self.boundary_length = self.mesh.boundary_length
+        self.tag_boundary_cells = self.mesh.tag_boundary_cells
         #self.number_of_full_nodes = self.mesh.number_of_full_nodes
         #self.number_of_full_triangles = self.mesh.number_of_full_triangles
@@ -1316 +1318 @@
         self._order_ = self.default_order
 
-        assert finaltime > self.get_starttime(), 'finaltime is less than starttime!'
+        assert finaltime >= self.get_starttime(), 'finaltime is less than starttime!'
         
         if finaltime is not None and duration is not None:
@@ -1691 +1693 @@
         return q_evol
     
-    def update_boundary(self):
+    def update_boundary_old(self):
         """Go through list of boundary objects and update boundary values
         for all conserved quantities on boundary.
@@ -1728 +1730 @@
                     Q = self.quantities[name]
                     Q.boundary_values[i] = q_evol[j]
+
+
+    def update_boundary_old_2(self):
+        """Go through list of boundary objects and update boundary values
+        for all conserved quantities on boundary.
+        It is assumed that the ordering of conserved quantities is
+        consistent between the domain and the boundary object, i.e.
+        the jth element of vector q must correspond to the jth conserved
+        quantity in domain.
+        """
+
+
+        for i in range(self.boundary_length):
+            vol_id  = self.boundary_cells[i]
+            edge_id = self.boundary_edges[i]
+            blah, B = self.boundary_objects[i]
+
+            if B is None:
+                log.critical('WARNING: Ignored boundary segment (None)')
+            else:
+                q_bdry = B.evaluate(vol_id, edge_id)
+
+                if len(q_bdry) == len(self.evolved_quantities):
+                    # conserved and evolved quantities are the same
+                    q_evol = q_bdry
+                elif len(q_bdry) == len(self.conserved_quantities):
+                    # boundary just returns conserved quantities
+                    # Need to calculate all the evolved quantities
+                    # Use default conversion
+
+                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
+
+                    q_evol = self.conserved_values_to_evolved_values \
+                                                            (q_bdry, q_evol)
+                else:
+                    msg = 'Boundary must return array of either conserved'
+                    msg += ' or evolved quantities'
+                    raise Exception(msg)
+
+                for j, name in enumerate(self.evolved_quantities):
+                    Q = self.quantities[name]
+                    Q.boundary_values[i] = q_evol[j]
+
+
+
+#        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+#            if B is None:
+#                log.critical('WARNING: Ignored boundary segment (None)')
+#            else:
+#                q_bdry = B.evaluate(vol_id, edge_id)
+#
+#                if len(q_bdry) == len(self.evolved_quantities):
+#                    # conserved and evolved quantities are the same
+#                    q_evol = q_bdry
+#                elif len(q_bdry) == len(self.conserved_quantities):
+#                    # boundary just returns conserved quantities
+#                    # Need to calculate all the evolved quantities
+#                    # Use default conversion
+#
+#                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#                    q_evol = self.conserved_values_to_evolved_values \
+#                                                            (q_bdry, q_evol)
+#                else:
+#                    msg = 'Boundary must return array of either conserved'
+#                    msg += ' or evolved quantities'
+#                    raise Exception(msg)
+#
+#                for j, name in enumerate(self.evolved_quantities):
+#                    Q = self.quantities[name]
+#                    Q.boundary_values[i] = q_evol[j]
+
+
+    def update_boundary(self):
+        """Go through list of boundary objects and update boundary values
+        for all conserved quantities on boundary.
+        It is assumed that the ordering of conserved quantities is
+        consistent between the domain and the boundary object, i.e.
+        the jth element of vector q must correspond to the jth conserved
+        quantity in domain.
+        """
+
+
+        for tag in self.tag_boundary_cells:
+
+            #print tag
+            
+            B = self.boundary_map[tag]
+            
+            #if B is None:
+            #        log.critical('WARNING: Ignored boundary segment (None)')
+
+            for i in self.tag_boundary_cells[tag]:
+                vol_id  = self.boundary_cells[i]
+                edge_id = self.boundary_edges[i]
+
+                if B is None:
+                    pass
+                else:
+                    q_bdry = B.evaluate(vol_id, edge_id)
+
+                    if len(q_bdry) == len(self.evolved_quantities):
+                        # conserved and evolved quantities are the same
+                        q_evol = q_bdry
+                    elif len(q_bdry) == len(self.conserved_quantities):
+                        # boundary just returns conserved quantities
+                        # Need to calculate all the evolved quantities
+                        # Use default conversion
+
+                        q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
+
+                        q_evol = self.conserved_values_to_evolved_values \
+                                                                (q_bdry, q_evol)
+                    else:
+                        msg = 'Boundary must return array of either conserved'
+                        msg += ' or evolved quantities'
+                        raise Exception(msg)
+
+                    for j, name in enumerate(self.evolved_quantities):
+                        Q = self.quantities[name]
+                        Q.boundary_values[i] = q_evol[j]
+
 
     def compute_fluxes(self):

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -191 +191 @@
         #Update boundary indices FIXME: OBSOLETE
         #self.build_boundary_structure()
+
+
 
         #FIXME check integrity?
@@ -376 +378 @@
 
         self.boundary = boundary
+        self.boundary_length = len(self.boundary)
 
 
@@ -402 +405 @@
         return self.tagged_elements
 
-    def build_boundary_structure(self):
-        """Traverse boundary and
-        enumerate neighbour indices from -1 and
-        counting down.
-
-        Precondition:
-            self.boundary is defined.
-        Post condition:
-            neighbour array has unique negative indices for boundary
-            boundary_segments array imposes an ordering on segments
-            (not otherwise available from the dictionary)
-
-        Note: If a segment is listed in the boundary dictionary
-        it *will* become a boundary - even if there is a neighbouring triangle.
-        This would be the case for internal boundaries
-        """
-
-        #FIXME: Now Obsolete - maybe use some comments from here in
-        #domain.set_boundary
-
-        if self.boundary is None:
-            msg = 'Boundary dictionary must be defined before '
-            msg += 'building boundary structure'
-            raise Exception(msg)
-
-
-        self.boundary_segments = self.boundary.keys()
-        self.boundary_segments.sort()
-
-        index = -1
-        for id, edge in self.boundary_segments:
-
-            #FIXME: One would detect internal boundaries as follows
-            #if self.neighbours[id, edge] > -1:
-            #    log.critical('Internal boundary')
-
-            self.neighbours[id, edge] = index
-
-            self.boundary_enumeration[id,edge] = index
-
-            index -= 1
-
+#    def build_boundary_structure(self):
+#        """Traverse boundary and
+#        enumerate neighbour indices from -1 and
+#        counting down.
+#
+#        Precondition:
+#            self.boundary is defined.
+#        Post condition:
+#            neighbour array has unique negative indices for boundary
+#            boundary_segments array imposes an ordering on segments
+#            (not otherwise available from the dictionary)
+#
+#        Note: If a segment is listed in the boundary dictionary
+#        it *will* become a boundary - even if there is a neighbouring triangle.
+#        This would be the case for internal boundaries
+#        """
+#
+#        #FIXME: Now Obsolete - maybe use some comments from here in
+#        #domain.set_boundary
+#
+#        if self.boundary is None:
+#            msg = 'Boundary dictionary must be defined before '
+#            msg += 'building boundary structure'
+#            raise Exception(msg)
+#
+#
+#        self.boundary_segments = self.boundary.keys()
+#        self.boundary_segments.sort()
+#
+#        index = -1
+#        for id, edge in self.boundary_segments:
+#
+#            #FIXME: One would detect internal boundaries as follows
+#            #if self.neighbours[id, edge] > -1:
+#            #    log.critical('Internal boundary')
+#
+#            self.neighbours[id, edge] = index
+#
+#            self.boundary_enumeration[id,edge] = index
+#
+#            index -= 1
+#
 
 
@@ -467 +470 @@
         self.boundary_enumeration = {}
 
+
+
         X = self.boundary.keys()
         X.sort()
 
+        #print 'X', X
         index = -1
         for id, edge in X:
@@ -484 +490 @@
 
         for id, edge in X:
-
             j = self.boundary_enumeration[id,edge]
             self.boundary_cells[j] = id
             self.boundary_edges[j] = edge
 
+        # For each tag create list of boundary edges
+        self.tag_boundary_cells = {}
+
+        tags = self.get_boundary_tags()
+
+        #print tags
+
+        for tag in tags:
+            self.tag_boundary_cells[tag] = []
+
+
+        for j in range(self.boundary_length):
+            id  = self.boundary_cells[j]
+            edge = self.boundary_edges[j]
+            tag = self.boundary[id, edge]
+            #print tag, id, edge
+            self.tag_boundary_cells[tag].append(j)
+
+
+        #print self.tag_boundary_cells
 
     def get_boundary_tags(self):

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -72 +72 @@
 
         # Allocate space for boundary values
-        self.boundary_length = L = len(domain.boundary)
+        #self.boundary_length = domain.boundary_length
+        self.boundary_length = L = self.domain.boundary_length
         self.boundary_values = num.zeros(L, num.float)
 

trunk/anuga_core/source/anuga/file/test_sww.py

@@ -41 +41 @@
         yiel=0.01
         points, vertices, boundary = rectangular(10,10)
+
+        #print "=============== boundary rect ======================="
+        #print boundary
 
         #Create shallow water domain
@@ -76 +79 @@
             pass
 
+        #print boundary
 
 
@@ -81 +85 @@
         domain2 = load_sww_as_domain(filename, None, fail_if_NaN=False,
                                         verbose=self.verbose)
+
+        # Unfortunately we loss the boundaries top, bottom, left and right,
+        # they are now all lumped into "exterior"
+
+        #print "=============== boundary domain2 ======================="
+        #print domain2.boundary
+        
+
+        #print domain2.get_boundary_tags()
+        
         #points, vertices, boundary = rectangular(15,15)
         #domain2.boundary = boundary
@@ -110 +124 @@
         domain.set_quantity('friction', 0.1)
         domain.store = False
-        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+        domain.set_boundary({'exterior': Bd, 'left' : Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
 
 
@@ -131 +145 @@
         #print 'domain2.boundary'
         #print domain2.boundary
-        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+        domain2.set_boundary({'exterior': Bd, 'left' : Bd,  'right': Bd, 'top': Bd, 'bottom': Bd})
         #domain2.set_boundary({'exterior': Bd})
 

trunk/anuga_core/source/anuga/operators/run_set_stage.py

@@ -60 +60 @@
 from anuga.operators.set_value_operators import Circular_set_stage_operator
 
-cop1 = Circular_set_stage_operator(domain, stage=h0+20.0, center=(0.0, 0.0), radius=100.0 )
-cop2 = Circular_set_stage_operator(domain, stage=h0+20.0, center=(2000.0, 1000.0), radius=100.0 )
+import math
+
+stage1 = lambda t: h0 + 20.0 * math.sin(t/3.0)
+cop1 = Circular_set_stage_operator(domain, stage=stage1, center=(0.0, 0.0), radius=100.0 )
+
+stage2 = lambda t: h0 + 30.0 * math.sin(t/6.0)
+cop2 = Circular_set_stage_operator(domain, stage=stage2, center=(2000.0, 1000.0), radius=100.0 )
 
 #print cop1.statistics()

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -304 +304 @@
         self.store = flag
 
+    def get_store(self):
+        """Get whether data saved to sww file.
+        """
+
+        return self.store
+
         
     def set_sloped_mannings_function(self, flag=True):
@@ -342 +348 @@
         """Get method for computing fluxes.
 
-        Currently
-           original
-           wb_1
+        See set_compute_fluxes_method for possible choices.
         """
 
@@ -1010 +1014 @@
 
 
+    
     def evolve(self,
                yieldstep=None,

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -451 +451 @@
 
     /*Compute fluxes between volumes for the shallow water wave equation
-      cast in terms of the 'stage', 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
+     * cast in terms of the 'stage', 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.
+     *
+     * The pressure term is calculated using simpson's rule so that it
+     * will well balance with the standard ghz_x gravity term
      */
 
@@ -464 +467 @@
     int i;
     //double hl, hr;
-    double uh_left, vh_left, u_left;
-    double uh_right, vh_right, u_right;
+    double uh_left, vh_left, u_left, v_left;
+    double uh_right, vh_right, u_right, v_right;
     double s_min, s_max, soundspeed_left, soundspeed_right;
     double denom, inverse_denominator;
@@ -531 +534 @@
     // Momentum in y-direction
     vh_left = q_left_rotated[2];
+    v_left = _compute_speed(&vh_left, &h_left,
+            epsilon, h0, limiting_threshold);
+
     vh_right = q_right_rotated[2];
+    v_right = _compute_speed(&vh_right, &h_right,
+            epsilon, h0, limiting_threshold);
 
     // Limit y-momentum if necessary
@@ -615 +623 @@
             edgeflux[i] *= inverse_denominator;
         }
+/*
+        if (edgeflux[0] > 0.0) {
+            edgeflux[2] = edgeflux[0]*v_left;
+        } else {
+            edgeflux[2] = edgeflux[0]*v_right;
+        }
+*/
 
         // Maximal wavespeed

trunk/anuga_core/source/anuga_parallel/parallel_api.py

@@ -62 +62 @@
         domain_name = domain.get_name()
         domain_dir = domain.get_datadir()
+        domain_store = domain.get_store()
         georef = domain.geo_reference
         number_of_global_triangles = domain.number_of_triangles
@@ -69 +70 @@
 
         for p in range(1, numprocs):
-            send((domain_name, domain_dir, georef, \
+            send((domain_name, domain_dir, domain_store, georef, \
                   number_of_global_triangles, number_of_global_nodes), p)
     else:
         if verbose: print 'P%d: Receiving domain attributes' %(myid)
 
-        domain_name, domain_dir, georef, \
+        domain_name, domain_dir, domain_store, georef, \
                   number_of_global_triangles, number_of_global_nodes = receive(0)
 
@@ -177 +178 @@
     #------------------------------------------------------------------------
     domain.set_name(domain_name)
-    domain.set_datadir(domain_dir)     
+    domain.set_datadir(domain_dir)
+    domain.set_store(domain_store)
     domain.geo_reference = georef   
 

trunk/anuga_core/source/anuga_parallel/parallel_generic_communications.py

@@ -33 +33 @@
     """
 
-    pypar.barrier()
+
 
     import time

trunk/anuga_core/source/anuga_parallel/parallel_shallow_water.py

@@ -18 +18 @@
 
 import numpy as num
+from os.path import join
 
 
@@ -131 +132 @@
     def sww_merge(self, verbose=False, delete_old=False):
 
-        if self.processor == 0 and self.numproc > 1:
+        if self.processor == 0 and self.numproc > 1 and self.store :
             import anuga.utilities.sww_merge as merge
+
+            global_name = join(self.get_datadir(),self.get_global_name())
             
-            merge.sww_merge_parallel(self.get_global_name(),self.numproc,verbose,delete_old)
+            merge.sww_merge_parallel(global_name,self.numproc,verbose,delete_old)
 
 

trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

@@ -85 +85 @@
     domain = create_domain_from_file(mesh_filename)
     domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
+    domain.set_datadir('Data')
+    domain.set_name('merimbula_new')
+    domain.set_store(False)
     #domain.set_quantity('elevation', Set_Elevation(500.0))
 
@@ -106 +109 @@
 #--------------------------------------------------------------------------
 
+    
 #domain.smooth = False
 domain.set_default_order(2)
@@ -111 +115 @@
 #domain.set_CFL(0.7)
 #domain.set_beta(1.5)
-domain.set_name('merimbula')
+
 
 #for p in range(numprocs):
@@ -144 +148 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
+for t in domain.evolve(yieldstep = yieldtime, finaltime = finaltime):
     if myid == 0:
         domain.write_time()