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

Last change on this file since 5991 was 5991, checked in by ole, 15 years ago

Work on startime and comments re ticket:306

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