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

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

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.

File size: 21.1 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
10import Numeric as num
11
12
13class Boundary:
14    """Base class for boundary conditions.
15       Specific boundary conditions must provide values for
16       the conserved_quantities
17
18       A boundary object has one neighbour; the one it
19       serves.
20    """
21
22    def __init__(self):
23        pass
24
25    def evaluate(self, vol_id=None, edge_id=None):
26        msg = 'Generic class Boundary must be subclassed'
27        raise Exception, msg
28
29
30class Transmissive_boundary(Boundary):
31    """Transmissive boundary returns same conserved quantities as
32    those present in its neighbour volume.
33
34    Underlying domain must be specified when boundary is instantiated
35    """
36
37    def __init__(self, domain = None):
38        Boundary.__init__(self)
39
40        if domain is None:
41            msg = 'Domain must be specified for transmissive boundary'
42            raise Exception, msg
43
44        self.domain = domain
45
46    def __repr__(self):
47        return 'Transmissive_boundary(%s)' %self.domain
48
49    def evaluate(self, vol_id, edge_id):
50        """Transmissive boundaries return the edge values
51        of the volume they serve.
52        """
53
54        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
55        return q
56
57
58class Dirichlet_boundary(Boundary):
59    """Dirichlet boundary returns constant values for the
60    conserved quantities
61    """
62
63
64    def __init__(self, conserved_quantities=None):
65        Boundary.__init__(self)
66
67        if conserved_quantities is None:
68            msg = 'Must specify one value for each conserved quantity'
69            raise Exception, msg
70
71        self.conserved_quantities=num.array(conserved_quantities, num.Float)
72
73    def __repr__(self):
74        return 'Dirichlet boundary (%s)' %self.conserved_quantities
75
76    def evaluate(self, vol_id=None, edge_id=None):
77        return self.conserved_quantities
78
79
80
81class Time_boundary(Boundary):
82    """Time dependent boundary returns values for the
83    conserved quantities as a function of time.
84    Must specify domain to get access to model time and a function of t
85    which must return conserved quantities as a function time.
86   
87    Example:
88      B = Time_boundary(domain,
89                        function=lambda t: [(60<t<3660)*2, 0, 0])
90     
91      This will produce a boundary condition with is a 2m high square wave
92      starting 60 seconds into the simulation and lasting one hour.
93      Momentum applied will be 0 at all times.
94                       
95    """
96
97    # FIXME (Ole): We should rename f to function to be consistent with
98    # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
99    def __init__(self, domain=None,
100                 f=None, # Should be removed and replaced by function below
101                 function=None,
102                 default_boundary=None,
103                 verbose=False):
104        Boundary.__init__(self)
105
106        self.default_boundary = default_boundary
107        self.default_boundary_invoked = False    # Flag
108        self.domain = domain
109        self.verbose = verbose
110
111       
112        # FIXME: Temporary code to deal with both f and function
113        if function is not None and f is not None:
114            raise Exception, 'Specify either function or f to Time_boundary'
115           
116        if function is None and f is None:
117            raise Exception, 'You must specify a function to Time_boundary'
118           
119        if f is None:
120            f = function
121        #####
122       
123        try:
124            q = f(0.0)
125        except Exception, e:
126            msg = 'Function for time boundary could not be executed:\n%s' %e
127            raise msg
128
129
130        try:
131            q = num.array(q, num.Float)
132        except:
133            msg = 'Return value from time boundary function could '
134            msg += 'not be converted into a Numeric array of floats.\n'
135            msg += 'Specified function should return either list or array.\n'
136            msg += 'I got %s' %str(q)
137            raise msg
138
139        msg = 'ERROR: Time boundary function must return a 1d list or array '
140        assert len(q.shape) == 1, msg
141
142        d = len(domain.conserved_quantities)
143        msg = 'Return value for f must be a list or an array of length %d' %d
144        assert len(q) == d, msg
145
146        self.f = f
147        self.domain = domain
148
149    def __repr__(self):
150        return 'Time boundary'
151
152    def evaluate(self, vol_id=None, edge_id=None):
153        # FIXME (Ole): I think this should be get_time(), see ticket:306
154        try:
155            res = self.f(self.domain.time)
156        except Modeltime_too_early, e:
157            raise Modeltime_too_early, e
158        except Modeltime_too_late, e:
159            if self.default_boundary is None:
160                raise Exception, e # Reraise exception
161            else:
162                # Pass control to default boundary
163                res = self.default_boundary.evaluate(vol_id, edge_id)
164               
165                # Ensure that result cannot be manipulated
166                # This is a real danger in case the
167                # default_boundary is a Dirichlet type
168                # for instance.
169                res = res.copy() 
170               
171                if self.default_boundary_invoked is False:
172                    if self.verbose:               
173                        # Issue warning the first time
174                        msg = '%s' %str(e)
175                        msg += 'Instead I will use the default boundary: %s\n'\
176                            %str(self.default_boundary) 
177                        msg += 'Note: Further warnings will be supressed'
178                        print msg
179               
180                    # FIXME (Ole): Replace this crude flag with
181                    # Python's ability to print warnings only once.
182                    # See http://docs.python.org/lib/warning-filter.html
183                    self.default_boundary_invoked = True
184
185        return res
186
187
188
189class File_boundary(Boundary):
190    """The File_boundary reads values for the conserved
191    quantities from an sww NetCDF file, and returns interpolated values
192    at the midpoints of each associated boundary segment.
193    Time dependency is interpolated linearly.
194
195    Assumes that file contains a time series and possibly
196    also spatial info. See docstring for File_function in util.py
197    for details about admissible file formats
198
199    File boundary must read and interpolate from *smoothed* version
200    as stored in sww and cannot work with the discontinuos triangles.
201
202    Example:
203    Bf = File_boundary('source_file.sww', domain)
204
205
206    Note that the resulting solution history is not exactly the same as if
207    the models were coupled as there is no feedback into the source model.
208   
209    Optional keyword argument default_boundary must be either None or
210    an instance of class descending from class Boundary.
211    This will be used in case model time exceeds that available in the
212    underlying data.
213       
214    """
215
216    # FIXME (Ole): I kind of like the name Spatio_Temporal_boundary for this
217    # rather than File_boundary
218
219    def __init__(self, filename, domain, 
220                 time_thinning=1, 
221                 time_limit=None,
222                 boundary_polygon=None,   
223                 default_boundary=None,
224                 use_cache=False, 
225                 verbose=False): 
226
227        import time
228        from anuga.config import time_format
229        from anuga.abstract_2d_finite_volumes.util import file_function
230
231        Boundary.__init__(self)
232
233        # Get x,y vertex coordinates for all triangles
234        V = domain.vertex_coordinates
235
236        # Compute midpoint coordinates for all boundary elements
237        # Only a subset may be invoked when boundary is evaluated but
238        # we don't know which ones at this stage since this object can
239        # be attached to
240        # any tagged boundary later on.
241
242        if verbose: print 'Find midpoint coordinates of entire boundary'
243        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.Float)
244        boundary_keys = domain.boundary.keys()
245
246        xllcorner = domain.geo_reference.get_xllcorner()
247        yllcorner = domain.geo_reference.get_yllcorner()       
248       
249
250        # Make ordering unique #FIXME: should this happen in domain.py?
251        boundary_keys.sort()
252
253        # Record ordering #FIXME: should this also happen in domain.py?
254        self.boundary_indices = {}
255        for i, (vol_id, edge_id) in enumerate(boundary_keys):
256
257            base_index = 3*vol_id
258            x0, y0 = V[base_index, :]
259            x1, y1 = V[base_index+1, :]
260            x2, y2 = V[base_index+2, :]
261           
262            # Compute midpoints
263            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
264            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
265            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
266
267            # Convert to absolute UTM coordinates
268            m[0] += xllcorner
269            m[1] += yllcorner
270           
271            # Register point and index
272            self.midpoint_coordinates[i,:] = m
273
274            # Register index of this boundary edge for use with evaluate
275            self.boundary_indices[(vol_id, edge_id)] = i
276
277        if verbose: print 'Initialise file_function'
278        self.F = file_function(filename,
279                               domain,
280                               quantities=domain.conserved_quantities,
281                               interpolation_points=self.midpoint_coordinates,
282                               time_thinning=time_thinning,
283                               time_limit=time_limit,
284                               use_cache=use_cache, 
285                               verbose=verbose,
286                               boundary_polygon=boundary_polygon)
287                             
288        # Check and store default_boundary
289        msg = 'Keyword argument default_boundary must be either None '
290        msg += 'or a boundary object.\n I got %s' %(str(default_boundary))
291        assert default_boundary is None or\
292            isinstance(default_boundary, Boundary), msg
293        self.default_boundary = default_boundary
294        self.default_boundary_invoked = False    # Flag
295
296        # Store pointer to domain and verbosity
297        self.domain = domain
298        self.verbose = verbose
299
300
301        # Here we'll flag indices outside the mesh as a warning
302        # as suggested by Joaquim Luis in sourceforge posting
303        # November 2007
304        # We won't make it an error as it is conceivable that
305        # only part of mesh boundary is actually used with a given
306        # file boundary sww file.
307        if hasattr(self.F, 'indices_outside_mesh') and\
308               len(self.F.indices_outside_mesh) > 0:
309            msg = 'WARNING: File_boundary has points outside the mesh '
310            msg += 'given in %s. ' %filename
311            msg += 'See warning message issued by Interpolation_function '
312            msg += 'for details (should appear above somewhere if '
313            msg += 'verbose is True).\n'
314            msg += 'This is perfectly OK as long as the points that are '
315            msg += 'outside aren\'t used on the actual boundary segment.'
316            if verbose is True:           
317                print msg
318            #raise Exception(msg)
319
320        # Test that file function can be called
321        q = self.F(0, point_id=0)
322        d = len(domain.conserved_quantities)
323        msg = 'Values specified in file %s must be ' %filename
324        msg += ' a list or an array of length %d' %d
325        assert len(q) == d, msg
326
327
328    def __repr__(self):
329        return 'File boundary'
330
331
332    def evaluate(self, vol_id=None, edge_id=None):
333        """Return linearly interpolated values based on domain.time
334        at midpoint of segment defined by vol_id and edge_id.
335        """
336
337        # FIXME (Ole): I think this should be get_time(), see ticket:306
338        t = self.domain.time
339       
340        if vol_id is not None and edge_id is not None:
341            i = self.boundary_indices[ vol_id, edge_id ]
342           
343            try:
344                res = self.F(t, point_id=i)
345            except Modeltime_too_early, e:
346                raise Modeltime_too_early, e
347            except Modeltime_too_late, e:
348                if self.default_boundary is None:
349                    raise Exception, e # Reraise exception
350                else:
351                    # Pass control to default boundary
352                    res = self.default_boundary.evaluate(vol_id, edge_id)
353                   
354                    # Ensure that result cannot be manipulated
355                    # This is a real danger in case the
356                    # default_boundary is a Dirichlet type
357                    # for instance.
358                    res = res.copy() 
359                   
360                    if self.default_boundary_invoked is False:
361                        # Issue warning the first time
362                        if self.verbose:
363                            msg = '%s' %str(e)
364                            msg += 'Instead I will use the default boundary: %s\n'\
365                                %str(self.default_boundary) 
366                            msg += 'Note: Further warnings will be supressed'
367                            print msg
368                   
369                        # FIXME (Ole): Replace this crude flag with
370                        # Python's ability to print warnings only once.
371                        # See http://docs.python.org/lib/warning-filter.html
372                        self.default_boundary_invoked = True
373                   
374
375            if res == NAN:
376                x,y=self.midpoint_coordinates[i,:]
377                msg = 'NAN value found in file_boundary at '
378                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
379
380                if hasattr(self.F, 'indices_outside_mesh') and\
381                       len(self.F.indices_outside_mesh) > 0:
382                    # Check if NAN point is due it being outside
383                    # boundary defined in sww file.
384
385                    if i in self.F.indices_outside_mesh:
386                        msg += 'This point refers to one outside the '
387                        msg += 'mesh defined by the file %s.\n'\
388                               %self.F.filename
389                        msg += 'Make sure that the file covers '
390                        msg += 'the boundary segment it is assigned to '
391                        msg += 'in set_boundary.'
392                    else:
393                        msg += 'This point is inside the mesh defined '
394                        msg += 'the file %s.\n' %self.F.filename
395                        msg += 'Check this file for NANs.'
396                raise Exception, msg
397           
398            return res
399        else:
400            msg = 'Boundary call without point_id not implemented.\n'
401            msg += 'vol_id=%s, edge_id=%s' %(str(vol_id), str(edge_id))
402            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
Note: See TracBrowser for help on using the repository browser.