Changeset 6928


Ignore:
Timestamp:
Apr 28, 2009, 5:41:15 PM (15 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:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r6717 r6928  
    1919     import File_boundary
    2020from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     21     import AWI_boundary
     22from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    2123     import Dirichlet_boundary
    2224from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     
    2426from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    2527     import Transmissive_boundary
     28
    2629from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
    2730from anuga.abstract_2d_finite_volumes.region\
  • anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6889 r6928  
    401401            msg += 'vol_id=%s, edge_id=%s' %(str(vol_id), str(edge_id))
    402402            raise Exception, msg
     403
     404class AWI_boundary(Boundary):
     405    """The AWI_boundary reads values for the conserved
     406    quantities (only STAGE) from an sww NetCDF file, and returns interpolated values
     407    at the midpoints of each associated boundary segment.
     408    Time dependency is interpolated linearly.
     409
     410    Assumes that file contains a time series and possibly
     411    also spatial info. See docstring for File_function in util.py
     412    for details about admissible file formats
     413
     414    AWI boundary must read and interpolate from *smoothed* version
     415    as stored in sww and cannot work with the discontinuos triangles.
     416
     417    Example:
     418    Ba = AWI_boundary('source_file.sww', domain)
     419
     420
     421    Note that the resulting solution history is not exactly the same as if
     422    the models were coupled as there is no feedback into the source model.
     423
     424    This was added by Nils Goseberg et al in April 2009
     425       
     426    """
     427
     428    def __init__(self, filename, domain, time_thinning=1,
     429                 use_cache=False, verbose=False):
     430        import time
     431        from Numeric import array, zeros, Float
     432        from anuga.config import time_format
     433        from anuga.abstract_2d_finite_volumes.util import file_function
     434
     435        Boundary.__init__(self)
     436
     437        # Get x,y vertex coordinates for all triangles
     438        V = domain.vertex_coordinates
     439
     440        # Compute midpoint coordinates for all boundary elements
     441        # Only a subset may be invoked when boundary is evaluated but
     442        # we don't know which ones at this stage since this object can
     443        # be attached to
     444        # any tagged boundary later on.
     445
     446        if verbose: print 'Find midpoint coordinates of entire boundary'
     447        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
     448        boundary_keys = domain.boundary.keys()
     449
     450
     451        xllcorner = domain.geo_reference.get_xllcorner()
     452        yllcorner = domain.geo_reference.get_yllcorner()       
     453       
     454
     455        # Make ordering unique #FIXME: should this happen in domain.py?
     456        boundary_keys.sort()
     457
     458        # Record ordering #FIXME: should this also happen in domain.py?
     459        self.boundary_indices = {}
     460        for i, (vol_id, edge_id) in enumerate(boundary_keys):
     461
     462            base_index = 3*vol_id
     463            x0, y0 = V[base_index, :]
     464            x1, y1 = V[base_index+1, :]
     465            x2, y2 = V[base_index+2, :]
     466           
     467            # Compute midpoints
     468            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
     469            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
     470            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
     471
     472            # Convert to absolute UTM coordinates
     473            m[0] += xllcorner
     474            m[1] += yllcorner
     475           
     476            # Register point and index
     477            self.midpoint_coordinates[i,:] = m
     478
     479            # Register index of this boundary edge for use with evaluate
     480            self.boundary_indices[(vol_id, edge_id)] = i
     481
     482
     483        if verbose: print 'Initialise file_function'
     484        self.F = file_function(filename, domain,
     485                               interpolation_points=self.midpoint_coordinates,
     486                               time_thinning=time_thinning,
     487                               use_cache=use_cache,
     488                               verbose=verbose)
     489        self.domain = domain
     490
     491        # Test
     492
     493        # Here we'll flag indices outside the mesh as a warning
     494        # as suggested by Joaquim Luis in sourceforge posting
     495        # November 2007
     496        # We won't make it an error as it is conceivable that
     497        # only part of mesh boundary is actually used with a given
     498        # file boundary sww file.
     499        if hasattr(self.F, 'indices_outside_mesh') and\
     500               len(self.F.indices_outside_mesh) > 0:
     501            msg = 'WARNING: File_boundary has points outside the mesh '
     502            msg += 'given in %s. ' %filename
     503            msg += 'See warning message issued by Interpolation_function '
     504            msg += 'for details (should appear above somewhere if '
     505            msg += 'verbose is True).\n'
     506            msg += 'This is perfectly OK as long as the points that are '
     507            msg += 'outside aren\'t used on the actual boundary segment.'
     508            if verbose is True:           
     509                print msg
     510            #raise Exception(msg)
     511
     512       
     513        q = self.F(0, point_id=0)
     514
     515        d = len(domain.conserved_quantities)
     516        msg = 'Values specified in file %s must be ' %filename
     517        msg += ' a list or an array of length %d' %d
     518        assert len(q) == d, msg
     519
     520
     521    def __repr__(self):
     522        return 'File boundary'
     523
     524
     525    def evaluate(self, vol_id=None, edge_id=None):
     526        """Return linearly interpolated values based on domain.time
     527        at midpoint of segment defined by vol_id and edge_id.
     528        """
     529        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
     530        t = self.domain.time
     531
     532        if vol_id is not None and edge_id is not None:
     533            i = self.boundary_indices[ vol_id, edge_id ]
     534            res = self.F(t, point_id = i)
     535
     536            if res == NAN:
     537                x,y=self.midpoint_coordinates[i,:]
     538                msg = 'NAN value found in file_boundary at '
     539                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
     540
     541                if hasattr(self.F, 'indices_outside_mesh') and\
     542                       len(self.F.indices_outside_mesh) > 0:
     543                    # Check if NAN point is due it being outside
     544                    # boundary defined in sww file.
     545
     546                    if i in self.F.indices_outside_mesh:
     547                        msg += 'This point refers to one outside the '
     548                        msg += 'mesh defined by the file %s.\n'\
     549                               %self.F.filename
     550                        msg += 'Make sure that the file covers '
     551                        msg += 'the boundary segment it is assigned to '
     552                        msg += 'in set_boundary.'
     553                    else:
     554                        msg += 'This point is inside the mesh defined '
     555                        msg += 'the file %s.\n' %self.F.filename
     556                        msg += 'Check this file for NANs.'
     557                raise Exception, msg
     558           
     559            q[0] = res[0] # Take stage, leave momentum alone
     560            return q
     561            #return res
     562        else:
     563            # raise 'Boundary call without point_id not implemented'
     564            # FIXME: What should the semantics be?
     565            return self.F(t)
     566
     567
  • anuga_core/source/anuga/shallow_water/__init__.py

    r6265 r6928  
    1313     Dirichlet_discharge_boundary,\
    1414     Field_boundary,\
    15      Transmissive_stage_zero_momentum_boundary
     15     Transmissive_stage_zero_momentum_boundary,\
     16     AWI_boundary
    1617
    1718
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6654 r6928  
    9494from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    9595     import Transmissive_boundary
     96from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     97     import AWI_boundary
     98from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     99     import AWI_boundary
     100
    96101from anuga.pmesh.mesh_interface import create_mesh_from_regions
    97102from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
Note: See TracChangeset for help on using the changeset viewer.