Changeset 6928

Timestamp:

Apr 28, 2009, 5:41:15 PM (16 years ago)

Author:

ole

Message:

Added AWI boundary as requested by Nils Goseberg (goseberg@…)
It does not yet have a unit test, but I have verified that all existing tests in ANUGA pass.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (2 diffs)
• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (1 diff)
• shallow_water/__init__.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -19 +19 @@
      import File_boundary
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import AWI_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
      import Dirichlet_boundary
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
@@ -24 +26 @@
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
      import Transmissive_boundary
+
 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
 from anuga.abstract_2d_finite_volumes.region\

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -401 +401 @@
             msg += 'vol_id=%s, edge_id=%s' %(str(vol_id), str(edge_id))
             raise Exception, msg
+
+class AWI_boundary(Boundary):
+    """The AWI_boundary reads values for the conserved
+    quantities (only STAGE) from an sww NetCDF file, and returns interpolated values
+    at the midpoints of each associated boundary segment.
+    Time dependency is interpolated linearly.
+
+    Assumes that file contains a time series and possibly
+    also spatial info. See docstring for File_function in util.py
+    for details about admissible file formats
+
+    AWI boundary must read and interpolate from *smoothed* version
+    as stored in sww and cannot work with the discontinuos triangles.
+
+    Example:
+    Ba = AWI_boundary('source_file.sww', domain)
+
+
+    Note that the resulting solution history is not exactly the same as if
+    the models were coupled as there is no feedback into the source model.
+
+    This was added by Nils Goseberg et al in April 2009
+       
+    """
+
+    def __init__(self, filename, domain, time_thinning=1, 
+                 use_cache=False, verbose=False):
+        import time
+        from Numeric import array, zeros, Float
+        from anuga.config import time_format
+        from anuga.abstract_2d_finite_volumes.util import file_function
+
+        Boundary.__init__(self)
+
+        # Get x,y vertex coordinates for all triangles
+        V = domain.vertex_coordinates
+
+        # Compute midpoint coordinates for all boundary elements
+        # Only a subset may be invoked when boundary is evaluated but
+        # we don't know which ones at this stage since this object can
+        # be attached to
+        # any tagged boundary later on.
+
+        if verbose: print 'Find midpoint coordinates of entire boundary'
+        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
+        boundary_keys = domain.boundary.keys()
+
+
+        xllcorner = domain.geo_reference.get_xllcorner()
+        yllcorner = domain.geo_reference.get_yllcorner()        
+        
+
+        # Make ordering unique #FIXME: should this happen in domain.py?
+        boundary_keys.sort()
+
+        # Record ordering #FIXME: should this also happen in domain.py?
+        self.boundary_indices = {}
+        for i, (vol_id, edge_id) in enumerate(boundary_keys):
+
+            base_index = 3*vol_id
+            x0, y0 = V[base_index, :]
+            x1, y1 = V[base_index+1, :]
+            x2, y2 = V[base_index+2, :]
+            
+            # Compute midpoints
+            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
+            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
+            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
+
+            # Convert to absolute UTM coordinates
+            m[0] += xllcorner
+            m[1] += yllcorner
+            
+            # Register point and index
+            self.midpoint_coordinates[i,:] = m
+
+            # Register index of this boundary edge for use with evaluate
+            self.boundary_indices[(vol_id, edge_id)] = i
+
+
+        if verbose: print 'Initialise file_function'
+        self.F = file_function(filename, domain,
+                               interpolation_points=self.midpoint_coordinates,
+                               time_thinning=time_thinning,
+                               use_cache=use_cache, 
+                               verbose=verbose)
+        self.domain = domain
+
+        # Test
+
+        # Here we'll flag indices outside the mesh as a warning
+        # as suggested by Joaquim Luis in sourceforge posting
+        # November 2007
+        # We won't make it an error as it is conceivable that
+        # only part of mesh boundary is actually used with a given
+        # file boundary sww file. 
+        if hasattr(self.F, 'indices_outside_mesh') and\
+               len(self.F.indices_outside_mesh) > 0:
+            msg = 'WARNING: File_boundary has points outside the mesh '
+            msg += 'given in %s. ' %filename
+            msg += 'See warning message issued by Interpolation_function '
+            msg += 'for details (should appear above somewhere if '
+            msg += 'verbose is True).\n'
+            msg += 'This is perfectly OK as long as the points that are '
+            msg += 'outside aren\'t used on the actual boundary segment.'
+            if verbose is True:            
+                print msg
+            #raise Exception(msg)
+
+        
+        q = self.F(0, point_id=0)
+
+        d = len(domain.conserved_quantities)
+        msg = 'Values specified in file %s must be ' %filename
+        msg += ' a list or an array of length %d' %d
+        assert len(q) == d, msg
+
+
+    def __repr__(self):
+        return 'File boundary'
+
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        """Return linearly interpolated values based on domain.time
+        at midpoint of segment defined by vol_id and edge_id.
+        """
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        t = self.domain.time
+
+        if vol_id is not None and edge_id is not None:
+            i = self.boundary_indices[ vol_id, edge_id ]
+            res = self.F(t, point_id = i)
+
+            if res == NAN:
+                x,y=self.midpoint_coordinates[i,:]
+                msg = 'NAN value found in file_boundary at '
+                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
+
+                if hasattr(self.F, 'indices_outside_mesh') and\
+                       len(self.F.indices_outside_mesh) > 0:
+                    # Check if NAN point is due it being outside
+                    # boundary defined in sww file.
+
+                    if i in self.F.indices_outside_mesh:
+                        msg += 'This point refers to one outside the '
+                        msg += 'mesh defined by the file %s.\n'\
+                               %self.F.filename
+                        msg += 'Make sure that the file covers '
+                        msg += 'the boundary segment it is assigned to '
+                        msg += 'in set_boundary.'
+                    else:
+                        msg += 'This point is inside the mesh defined '
+                        msg += 'the file %s.\n' %self.F.filename
+                        msg += 'Check this file for NANs.'
+                raise Exception, msg
+            
+            q[0] = res[0] # Take stage, leave momentum alone
+            return q
+            #return res 
+        else:
+            # raise 'Boundary call without point_id not implemented'
+            # FIXME: What should the semantics be?
+            return self.F(t)
+
+

anuga_core/source/anuga/shallow_water/__init__.py

@@ -13 +13 @@
      Dirichlet_discharge_boundary,\
      Field_boundary,\
-     Transmissive_stage_zero_momentum_boundary
+     Transmissive_stage_zero_momentum_boundary,\
+     AWI_boundary
 
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -94 +94 @@
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
      import Transmissive_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import AWI_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import AWI_boundary
+
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric