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

Last change on this file since 4712 was 4618, checked in by ole, 17 years ago

Added comment

File size: 11.3 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
84    def __init__(self, domain = None, f = None):
85        Boundary.__init__(self)
86
87        try:
88            q = f(0.0)
89        except Exception, e:
90            msg = 'Function for time boundary could not be executed:\n%s' %e
91            raise msg
92
93
94        from Numeric import array, Float
95        try:
96            q = array(q).astype(Float)
97        except:
98            msg = 'Return value from time boundary function could '
99            msg += 'not be converted into a Numeric array of floats.\n'
100            msg += 'Specified function should return either list or array.'
101            raise msg
102
103        msg = 'ERROR: Time boundary function must return a 1d list or array '
104        assert len(q.shape) == 1, msg
105
106        d = len(domain.conserved_quantities)
107        msg = 'Return value for f must be a list or an array of length %d' %d
108        assert len(q) == d, msg
109
110        self.f = f
111        self.domain = domain
112
113    def __repr__(self):
114        return 'Time boundary'
115
116    def evaluate(self, vol_id=None, edge_id=None):
117        return self.f(self.domain.time)
118
119
120class File_boundary_time(Boundary):
121    """Boundary values obtained from file and interpolated.
122    conserved quantities as a function of time.
123
124    Assumes that file contains a time series.
125
126    No spatial info assumed.
127    """
128
129    #FIXME: Is this still necessary
130
131    def __init__(self, filename, domain):
132        import time
133        from Numeric import array
134        from anuga.config import time_format
135        from anuga.abstract_2d_finite_volumes.util import File_function
136
137        Boundary.__init__(self)
138
139        self.F = File_function(filename, domain)
140        self.domain = domain
141
142        #Test
143        q = self.F(0)
144
145        d = len(domain.conserved_quantities)
146        msg = 'Values specified in file must be a list or an array of length %d' %d
147        assert len(q) == d, msg
148
149
150    def __repr__(self):
151        return 'File boundary'
152
153    def evaluate(self, vol_id=None, edge_id=None):
154        """Return linearly interpolated values based on domain.time
155
156        vol_id and edge_id are ignored
157        """
158
159        t = self.domain.time
160        return self.F(t)
161
162
163
164
165class File_boundary(Boundary):
166    """The File_boundary reads values for the conserved
167    quantities from an sww NetCDF file, and returns interpolated values
168    at the midpoints of each associated boundary segment.
169    Time dependency is interpolated linearly.
170
171    Assumes that file contains a time series and possibly
172    also spatial info. See docstring for File_function in util.py
173    for details about admissible file formats
174
175    File boundary must read and interpolate from *smoothed* version
176    as stored in sww and cannot work with the discontinuos triangles.
177
178    Example:
179    Bf = File_boundary('source_file.sww', domain)
180
181
182    Note that the resulting solution history is not exactly the same as if
183    the models were coupled as there is no feedback into the source model.
184       
185    """
186
187    # FIXME (Ole): I kind of like the name Spatio_Temporal_boundary for this
188    # rather than File_boundary
189
190    def __init__(self, filename, domain, time_thinning=1, 
191                 use_cache=False, verbose=False):
192        import time
193        from Numeric import array, zeros, Float
194        from anuga.config import time_format
195        from anuga.abstract_2d_finite_volumes.util import file_function
196
197        Boundary.__init__(self)
198
199        #Get x,y vertex coordinates for all triangles
200        V = domain.vertex_coordinates
201
202        #Compute midpoint coordinates for all boundary elements
203        #Only a subset may be invoked when boundary is evaluated but
204        #we don't know which ones at this stage since this object can
205        #be attached to
206        #any tagged boundary later on.
207
208        if verbose: print 'Find midpoint coordinates of entire boundary'
209        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
210        boundary_keys = domain.boundary.keys()
211
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
245        if verbose: print 'Initialise file_function'
246        self.F = file_function(filename, domain,
247                               interpolation_points=self.midpoint_coordinates,
248                               time_thinning=time_thinning,
249                               use_cache=use_cache, 
250                               verbose=verbose)
251        self.domain = domain
252
253        #Test
254        q = self.F(0, point_id=0)
255
256        d = len(domain.conserved_quantities)
257        msg = 'Values specified in file %s must be ' %filename
258        msg += ' a list or an array of length %d' %d
259        assert len(q) == d, msg
260
261
262    def __repr__(self):
263        return 'File boundary'
264
265
266    def evaluate(self, vol_id=None, edge_id=None):
267        """Return linearly interpolated values based on domain.time
268        at midpoint of segment defined by vol_id and edge_id.
269        """
270
271        t = self.domain.time
272
273        if vol_id is not None and edge_id is not None:
274            i = self.boundary_indices[ vol_id, edge_id ]
275            res = self.F(t, point_id = i)
276
277            if res == NAN:
278                x,y=self.midpoint_coordinates[i,:]
279                msg = 'NAN value found in file_boundary at '
280                msg += 'point id #%d: (%.2f, %.2f)' %(i, x, y)
281                #print msg
282                raise Exception, msg
283           
284            return res
285        else:
286            #raise 'Boundary call without point_id not implemented'
287            #FIXME: What should the semantics be?
288            return self.F(t)
289
290
291
292
293
294
295
296
297#THIS FAR (10/8/4)
298class Connective_boundary(Boundary):
299    """Connective boundary returns values for the
300    conserved quantities from a volume as defined by a connection table
301    mapping between tuples of (volume id, face id) for volumes that
302    have their boundaries connected.
303
304    FIXME: Perhaps include possibility of mapping between
305    different domains as well
306
307    FIXME: In case of shallow water we may want to have a
308    special version that casts this in terms of height rather than stage
309    """
310
311
312    def __init__(self, table):
313        from domain import Volume
314
315        Boundary.__init__(self)
316
317        self.connection_table = table
318        self.Volume = Volume
319
320    def __repr__(self):
321        return 'Connective boundary'
322
323    #FIXME: IF we ever need to get field_values from connected volume,
324    #that method could be overridden here (using same idea as in
325    #get_conserved_quantities
326    #def get_field_values()
327
328    def get_conserved_quantities(self, volume, face=0):
329
330        id = volume.id
331        if self.connection_table.has_key((id, face)):
332            other_id, other_face = self.connection_table[(id, face)]
333
334            other_volume = self.Volume.instances[other_id]
335            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
336            exec(cmd)
337            return q
338        else:
339            msg = 'Volume, face tuple ($d, %d) has not been mapped'\
340                  %(id, face)
341            raise msg
342
343
344
345
346
347#FIXME: Add a boundary with a general function of x,y, and t
348
349#FIXME: Add periodic boundaries e.g.:
350# Attempt at periodic conditions from advection_spik. Remember this
351#
352#first = 2*(N-1)*N
353#for i in range(1,2*N+1,2):
354#    k = first + i-1#
355#
356#    print i,k
357#
358#    domain[i].faces[2].neighbour = domain[k].faces[1]
359#    domain[k].faces[1].neighbour = domain[i].faces[2]
360
361
362
363class General_boundary(Boundary):
364    """General boundary which can compute conserved quantities based on
365    their previous value, conserved quantities of its neighbour and model time.
366
367    Must specify initial conserved quantities,
368    neighbour,
369    domain to get access to model time
370    a function f(q_old, neighbours_q, t) which must return
371    new conserved quantities q as a function time
372
373    FIXME: COMPLETE UNTESTED - JUST AN IDEA
374    """
375
376    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
377        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)
378
379        self.f = f
380        self.domain = domain
381
382
383    def get_conserved_quantities(self, volume=None, face=0):
384        return self.f(self.conserved_quantities,
385                      neighbour.conserved_quantities,
386                      self.domain.time)
387
388
389
390
Note: See TracBrowser for help on using the repository browser.