source: anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py @ 5657

Last change on this file since 5657 was 5657, checked in by ole, 16 years ago

Implemented default_boundary option in File_boundary and Field_boundary as
per ticket:293 and added a note in the user manual.

File size: 15.6 KB
Line 
1"""boundary.py - Classes for implementing boundary conditions
2"""
3
4from warnings import warn
5
6from anuga.utilities.numerical_tools import NAN   
7from anuga.fit_interpolate.interpolate import Modeltime_too_late
8from anuga.fit_interpolate.interpolate import Modeltime_too_early
9
10
11class Boundary:
12    """Base class for boundary conditions.
13       Specific boundary conditions must provide values for
14       the conserved_quantities
15
16       A boundary object has one neighbour; the one it
17       serves.
18    """
19
20    def __init__(self):
21        pass
22
23    def evaluate(self, vol_id=None, edge_id=None):
24        msg = 'Generic class Boundary must be subclassed'
25        raise Exception, msg
26
27
28class Transmissive_boundary(Boundary):
29    """Transmissive boundary returns same conserved quantities as
30    those present in its neighbour volume.
31
32    Underlying domain must be specified when boundary is instantiated
33    """
34
35    def __init__(self, domain = None):
36        Boundary.__init__(self)
37
38        if domain is None:
39            msg = 'Domain must be specified for transmissive boundary'
40            raise Exception, msg
41
42        self.domain = domain
43
44    def __repr__(self):
45        return 'Transmissive_boundary(%s)' %self.domain
46
47    def evaluate(self, vol_id, edge_id):
48        """Transmissive boundaries return the edge values
49        of the volume they serve.
50        """
51
52        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
53        return q
54
55
56class Dirichlet_boundary(Boundary):
57    """Dirichlet boundary returns constant values for the
58    conserved quantities
59    """
60
61
62    def __init__(self, conserved_quantities=None):
63        Boundary.__init__(self)
64
65        if conserved_quantities is None:
66            msg = 'Must specify one value for each conserved quantity'
67            raise Exception, msg
68
69        from Numeric import array, Float
70        self.conserved_quantities=array(conserved_quantities).astype(Float)
71
72    def __repr__(self):
73        return 'Dirichlet boundary (%s)' %self.conserved_quantities
74
75    def evaluate(self, vol_id=None, edge_id=None):
76        return self.conserved_quantities
77
78
79
80class Time_boundary(Boundary):
81    """Time dependent boundary returns values for the
82    conserved quantities as a function of time.
83    Must specify domain to get access to model time and a function f(t)
84    which must return conserved quantities as a function time
85    """
86
87    # FIXME (Ole): We should rename f to function to be consistent with
88    # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
89    def __init__(self, domain = None, f = None):
90        Boundary.__init__(self)
91
92        try:
93            q = f(0.0)
94        except Exception, e:
95            msg = 'Function for time boundary could not be executed:\n%s' %e
96            raise msg
97
98
99        from Numeric import array, Float
100        try:
101            q = array(q).astype(Float)
102        except:
103            msg = 'Return value from time boundary function could '
104            msg += 'not be converted into a Numeric array of floats.\n'
105            msg += 'Specified function should return either list or array.'
106            raise msg
107
108        msg = 'ERROR: Time boundary function must return a 1d list or array '
109        assert len(q.shape) == 1, msg
110
111        d = len(domain.conserved_quantities)
112        msg = 'Return value for f must be a list or an array of length %d' %d
113        assert len(q) == d, msg
114
115        self.f = f
116        self.domain = domain
117
118    def __repr__(self):
119        return 'Time boundary'
120
121    def evaluate(self, vol_id=None, edge_id=None):
122        return self.f(self.domain.time)
123
124
125class File_boundary_time(Boundary):
126    """Boundary values obtained from file and interpolated.
127    conserved quantities as a function of time.
128
129    Assumes that file contains a time series.
130
131    No spatial info assumed.
132    """
133
134    #FIXME: Is this still necessary
135
136    def __init__(self, filename, domain):
137        import time
138        from Numeric import array
139        from anuga.config import time_format
140        from anuga.abstract_2d_finite_volumes.util import File_function
141
142        Boundary.__init__(self)
143
144        self.F = File_function(filename, domain)
145        self.domain = domain
146
147        #Test
148        q = self.F(0)
149
150        d = len(domain.conserved_quantities)
151        msg = 'Values specified in file must be a list or an array of length %d' %d
152        assert len(q) == d, msg
153
154
155    def __repr__(self):
156        return 'File boundary'
157
158    def evaluate(self, vol_id=None, edge_id=None):
159        """Return linearly interpolated values based on domain.time
160
161        vol_id and edge_id are ignored
162        """
163
164        t = self.domain.time
165        return self.F(t)
166
167
168
169
170class File_boundary(Boundary):
171    """The File_boundary reads values for the conserved
172    quantities from an sww NetCDF file, and returns interpolated values
173    at the midpoints of each associated boundary segment.
174    Time dependency is interpolated linearly.
175
176    Assumes that file contains a time series and possibly
177    also spatial info. See docstring for File_function in util.py
178    for details about admissible file formats
179
180    File boundary must read and interpolate from *smoothed* version
181    as stored in sww and cannot work with the discontinuos triangles.
182
183    Example:
184    Bf = File_boundary('source_file.sww', domain)
185
186
187    Note that the resulting solution history is not exactly the same as if
188    the models were coupled as there is no feedback into the source model.
189   
190    Optional keyword argument default_boundary must be either None or
191    an instance of class descending from class Boundary.
192    This will be used in case model time exceeds that available in the
193    underlying data.
194       
195    """
196
197    # FIXME (Ole): I kind of like the name Spatio_Temporal_boundary for this
198    # rather than File_boundary
199
200    def __init__(self, filename, domain, 
201                 time_thinning=1, 
202                 boundary_polygon=None,   
203                 default_boundary=None,
204                 use_cache=False, 
205                 verbose=False): 
206
207        import time
208        from Numeric import array, zeros, Float
209        from anuga.config import time_format
210        from anuga.abstract_2d_finite_volumes.util import file_function
211
212        Boundary.__init__(self)
213
214        # Get x,y vertex coordinates for all triangles
215        V = domain.vertex_coordinates
216
217        # Compute midpoint coordinates for all boundary elements
218        # Only a subset may be invoked when boundary is evaluated but
219        # we don't know which ones at this stage since this object can
220        # be attached to
221        # any tagged boundary later on.
222
223        if verbose: print 'Find midpoint coordinates of entire boundary'
224        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
225        boundary_keys = domain.boundary.keys()
226
227        xllcorner = domain.geo_reference.get_xllcorner()
228        yllcorner = domain.geo_reference.get_yllcorner()       
229       
230
231        # Make ordering unique #FIXME: should this happen in domain.py?
232        boundary_keys.sort()
233
234        # Record ordering #FIXME: should this also happen in domain.py?
235        self.boundary_indices = {}
236        for i, (vol_id, edge_id) in enumerate(boundary_keys):
237
238            base_index = 3*vol_id
239            x0, y0 = V[base_index, :]
240            x1, y1 = V[base_index+1, :]
241            x2, y2 = V[base_index+2, :]
242           
243            # Compute midpoints
244            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
245            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
246            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
247
248            # Convert to absolute UTM coordinates
249            m[0] += xllcorner
250            m[1] += yllcorner
251           
252            # Register point and index
253            self.midpoint_coordinates[i,:] = m
254
255            # Register index of this boundary edge for use with evaluate
256            self.boundary_indices[(vol_id, edge_id)] = i
257
258        if verbose: print 'Initialise file_function'
259        self.F = file_function(filename,
260                               domain,
261                               quantities=domain.conserved_quantities,
262                               interpolation_points=self.midpoint_coordinates,
263                               time_thinning=time_thinning,
264                               use_cache=use_cache, 
265                               verbose=verbose,
266                               boundary_polygon=boundary_polygon)
267                             
268        # Check and store default_boundary
269        msg = 'Keyword argument default_boundary must be either None '
270        msg += 'or a boundary object.\n I got %s' %(str(default_boundary))
271        assert default_boundary is None or isinstance(default_boundary, Boundary), msg
272        self.default_boundary = default_boundary
273        self.default_boundary_invoked = False    # Flag
274
275        # Store pointer to domain
276        self.domain = domain
277       
278        self.verbose = verbose
279
280        # Test
281
282        # Here we'll flag indices outside the mesh as a warning
283        # as suggested by Joaquim Luis in sourceforge posting
284        # November 2007
285        # We won't make it an error as it is conceivable that
286        # only part of mesh boundary is actually used with a given
287        # file boundary sww file.
288        if hasattr(self.F, 'indices_outside_mesh') and\
289               len(self.F.indices_outside_mesh) > 0:
290            msg = 'WARNING: File_boundary has points outside the mesh '
291            msg += 'given in %s. ' %filename
292            msg += 'See warning message issued by Interpolation_function '
293            msg += 'for details (should appear above somewhere if '
294            msg += 'verbose is True).\n'
295            msg += 'This is perfectly OK as long as the points that are '
296            msg += 'outside aren\'t used on the actual boundary segment.'
297            if verbose is True:           
298                print msg
299            #raise Exception(msg)
300
301        # Test that file function can be called
302        q = self.F(0, point_id=0)
303        d = len(domain.conserved_quantities)
304        msg = 'Values specified in file %s must be ' %filename
305        msg += ' a list or an array of length %d' %d
306        assert len(q) == d, msg
307
308
309    def __repr__(self):
310        return 'File boundary'
311
312
313    def evaluate(self, vol_id=None, edge_id=None):
314        """Return linearly interpolated values based on domain.time
315        at midpoint of segment defined by vol_id and edge_id.
316        """
317
318        t = self.domain.time
319       
320        if vol_id is not None and edge_id is not None:
321            i = self.boundary_indices[ vol_id, edge_id ]
322           
323           
324            #res = self.F(t, point_id = i)           
325            try:
326                res = self.F(t, point_id = i)
327            except (Modeltime_too_late, Modeltime_too_early), e:
328                if self.default_boundary is None:
329                    raise Exception, e # Reraise exception
330                else:
331                    if self.default_boundary_invoked is False:
332                        # Issue warning the first time
333                        msg = '%s' %str(e)
334                        msg += 'Instead I will using the default boundary: %s\n'\
335                            %str(self.default_boundary) 
336                        msg += 'Note: Further warnings will be supressed'
337                        warn(msg)
338                   
339                        # FIXME (Ole): Replace this crude flag with
340                        # Python's ability to print warnings only once.
341                        # See http://docs.python.org/lib/warning-filter.html
342                        self.default_boundary_invoked = True
343                   
344                    # Pass control to default boundary
345                    res = self.default_boundary.evaluate(vol_id, edge_id)
346                   
347               
348            if res == NAN:
349                x,y=self.midpoint_coordinates[i,:]
350                msg = 'NAN value found in file_boundary at '
351                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
352
353                if hasattr(self.F, 'indices_outside_mesh') and\
354                       len(self.F.indices_outside_mesh) > 0:
355                    # Check if NAN point is due it being outside
356                    # boundary defined in sww file.
357
358                    if i in self.F.indices_outside_mesh:
359                        msg += 'This point refers to one outside the '
360                        msg += 'mesh defined by the file %s.\n'\
361                               %self.F.filename
362                        msg += 'Make sure that the file covers '
363                        msg += 'the boundary segment it is assigned to '
364                        msg += 'in set_boundary.'
365                    else:
366                        msg += 'This point is inside the mesh defined '
367                        msg += 'the file %s.\n' %self.F.filename
368                        msg += 'Check this file for NANs.'
369                raise Exception, msg
370           
371            return res
372        else:
373            # raise 'Boundary call without point_id not implemented'
374            # FIXME: What should the semantics be?
375            return self.F(t)
376
377
378
379
380
381
382
383
384#THIS FAR (10/8/4)
385class Connective_boundary(Boundary):
386    """Connective boundary returns values for the
387    conserved quantities from a volume as defined by a connection table
388    mapping between tuples of (volume id, face id) for volumes that
389    have their boundaries connected.
390
391    FIXME: Perhaps include possibility of mapping between
392    different domains as well
393
394    FIXME: In case of shallow water we may want to have a
395    special version that casts this in terms of height rather than stage
396    """
397
398
399    def __init__(self, table):
400        from domain import Volume
401
402        Boundary.__init__(self)
403
404        self.connection_table = table
405        self.Volume = Volume
406
407    def __repr__(self):
408        return 'Connective boundary'
409
410    #FIXME: IF we ever need to get field_values from connected volume,
411    #that method could be overridden here (using same idea as in
412    #get_conserved_quantities
413    #def get_field_values()
414
415    def get_conserved_quantities(self, volume, face=0):
416
417        id = volume.id
418        if self.connection_table.has_key((id, face)):
419            other_id, other_face = self.connection_table[(id, face)]
420
421            other_volume = self.Volume.instances[other_id]
422            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
423            exec(cmd)
424            return q
425        else:
426            msg = 'Volume, face tuple ($d, %d) has not been mapped'\
427                  %(id, face)
428            raise msg
429
430
431
432
433
434#FIXME: Add a boundary with a general function of x,y, and t
435
436#FIXME: Add periodic boundaries e.g.:
437# Attempt at periodic conditions from advection_spik. Remember this
438#
439#first = 2*(N-1)*N
440#for i in range(1,2*N+1,2):
441#    k = first + i-1#
442#
443#    print i,k
444#
445#    domain[i].faces[2].neighbour = domain[k].faces[1]
446#    domain[k].faces[1].neighbour = domain[i].faces[2]
447
448
449
450class General_boundary(Boundary):
451    """General boundary which can compute conserved quantities based on
452    their previous value, conserved quantities of its neighbour and model time.
453
454    Must specify initial conserved quantities,
455    neighbour,
456    domain to get access to model time
457    a function f(q_old, neighbours_q, t) which must return
458    new conserved quantities q as a function time
459
460    FIXME: COMPLETE UNTESTED - JUST AN IDEA
461    """
462
463    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
464        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)
465
466        self.f = f
467        self.domain = domain
468
469
470    def get_conserved_quantities(self, volume=None, face=0):
471        return self.f(self.conserved_quantities,
472                      neighbour.conserved_quantities,
473                      self.domain.time)
474
475
476
477
Note: See TracBrowser for help on using the repository browser.