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

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

Better error message in Time_boundary

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