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

Last change on this file since 5521 was 5373, checked in by jakeman, 16 years ago

forgot file in last commit

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