source:anuga_work/development/shallow_water_1d/generic_boundary_conditions.py@4960

Last change on this file since 4960 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 14.5 KB
Line
1"""boundary.py - Classes for implementing boundary conditions
2"""
3
4
5class Boundary:
6    """Base class for boundary conditions.
7       Specific boundary conditions must provide values for
8       the conserved_quantities
9
10       A boundary object has one neighbour; the one it
11       serves.
12    """
13
14    def __init__(self):
15        pass
16
17    def evaluate(self, vol_id=None, edge_id=None):
18        msg = 'Generic class Boundary must be subclassed'
19        raise msg
20
21
22class Transmissive_boundary(Boundary):
23    """Transmissive boundary returns same conserved quantities as
24    those present in its neighbour volume.
25
26    Underlying domain must be specified when boundary is instantiated
27    """
28
29    def __init__(self, domain = None):
30        Boundary.__init__(self)
31
32        if domain is None:
33            msg = 'Domain must be specified for transmissive boundary'
34            raise msg
35
36        self.domain = domain
37
38    def __repr__(self):
39        return 'Transmissive_boundary(%s)' %self.domain
40
41    def evaluate(self, vol_id, edge_id):
42        """Transmissive boundaries return the edge values
43        of the volume they serve.
44        """
45
46        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
47        return q
48
49
50class Dirichlet_boundary(Boundary):
51    """Dirichlet boundary returns constant values for the
52    conserved quantities
53    """
54
55
56    def __init__(self, conserved_quantities=None):
57        Boundary.__init__(self)
58
59        if conserved_quantities is None:
60            msg = 'Must specify one value for each conserved quantity'
61            raise msg
62
63        from Numeric import array, Float
64        self.conserved_quantities=array(conserved_quantities).astype(Float)
65
66    def __repr__(self):
67        return 'Dirichlet boundary (%s)' %self.conserved_quantities
68
69    def evaluate(self, vol_id=None, edge_id=None):
70        return self.conserved_quantities
71
72
73
74class Time_boundary(Boundary):
75    """Time dependent boundary returns values for the
76    conserved quantities as a function of time.
77    Must specify domain to get access to model time and a function f(t)
78    which must return conserved quantities as a function time
79    """
80
81
82    def __init__(self, domain = None, f = None):
83        Boundary.__init__(self)
84
85        try:
86            q = f(0.0)
87        except Exception, e:
88            msg = 'Function for time boundary could not be executed:\n%s' %e
89            raise msg
90
91
92        from Numeric import array, Float
93        try:
94            q = array(q).astype(Float)
95        except:
96            msg = 'Return value from time boundary function could '
97            msg += 'not be converted into a Numeric array of floats.\n'
98            msg += 'Specified function should return either list or array.'
99            raise msg
100
101        msg = 'ERROR: Time boundary function must return a 1d list or array '
102        assert len(q.shape) == 1, msg
103
104        d = len(domain.conserved_quantities)
105        msg = 'Return value for f must be a list or an array of length %d' %d
106        assert len(q) == d, msg
107
108        self.f = f
109        self.domain = domain
110
111    def __repr__(self):
112        return 'Time boundary'
113
114    def evaluate(self, vol_id=None, edge_id=None):
115        return self.f(self.domain.time)
116
117
118class File_boundary_time(Boundary):
119    """Boundary values obtained from file and interpolated.
120    conserved quantities as a function of time.
121
122    Assumes that file contains a time series.
123
124    No spatial info assumed.
125    """
126
127    #FIXME: Is this still necessary
128
129    def __init__(self, filename, domain):
130        import time
131        from Numeric import array
132        from config import time_format
133        from util import File_function
134
135        Boundary.__init__(self)
136
137        self.F = File_function(filename, domain)
138        self.domain = domain
139
140        #Test
141        q = self.F(0)
142
143        d = len(domain.conserved_quantities)
144        msg = 'Values specified in file must be a list or an array of length %d' %d
145        assert len(q) == d, msg
146
147
148    def __repr__(self):
149        return 'File boundary'
150
151    def evaluate(self, vol_id=None, edge_id=None):
152        """Return linearly interpolated values based on domain.time
153
154        vol_id and edge_id are ignored
155        """
156
157        t = self.domain.time
158        return self.F(t)
159
160
161
162
163class File_boundary(Boundary):
164    """Boundary values obtained from file and interpolated.
165    conserved quantities as a function of time.
166
167    Assumes that file contains a time series and possibly
168    also spatial info.
169    See docstring for File_function in util.py for details about
171
172    The full solution history is not exactly the same as if
173    the models were coupled:
174    File boundary must read and interpolate from *smoothed* version
175    as stored in sww and cannot work with the discontinuos triangles.
176
177    """
178
179    def __init__(self, filename, domain, verbose = False):
180        import time
181        from Numeric import array, zeros, Float
182        from config import time_format
183        from util import file_function
184
185        Boundary.__init__(self)
186
187        #Get x,y vertex coordinates for all triangles
188        V = domain.vertex_coordinates
189
190        #Compute midpoint coordinates for all boundary elements
191        #Only a subset may be invoked when boundary is evaluated but
192        #we don't know which ones at this stage since this object can be attached to
193        #any tagged boundary later on.
194
195        if verbose: print 'Find midpoint coordinates of entire boundary'
196        self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
197        boundary_keys = domain.boundary.keys()
198
199
200        xllcorner = domain.geo_reference.get_xllcorner()
201        yllcorner = domain.geo_reference.get_yllcorner()
202
203
204        #Make ordering unique #FIXME: should this happen in domain.py?
205        boundary_keys.sort()
206
207
208        #Record ordering #FIXME: should this also happen in domain.py?
209        self.boundary_indices = {}
210        for i, (vol_id, edge_id) in enumerate(boundary_keys):
211
212            x0 = V[vol_id, 0]; y0 = V[vol_id, 1]
213            x1 = V[vol_id, 2]; y1 = V[vol_id, 3]
214            x2 = V[vol_id, 4]; y2 = V[vol_id, 5]
215
216            #Compute midpoints
217            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
218            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
219            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
220
221            #Convert to absolute UTM coordinates
222            m[0] += xllcorner
223            m[1] += yllcorner
224
225            #Register point and index
226            self.midpoint_coordinates[i,:] = m
227            self.boundary_indices[(vol_id, edge_id)] = i
228
229
230        if verbose: print 'Initialise file_function'
231        self.F = file_function(filename, domain,
232                               interpolation_points=self.midpoint_coordinates,
233                               verbose = verbose)
234        self.domain = domain
235
236        #Test
237        q = self.F(0, point_id=0)
238
239        d = len(domain.conserved_quantities)
240        msg = 'Values specified in file %s must be ' %filename
241        msg += ' a list or an array of length %d' %d
242        assert len(q) == d, msg
243
244
245    def __repr__(self):
246        return 'File boundary'
247
248
249    def evaluate(self, vol_id=None, edge_id=None):
250        """Return linearly interpolated values based on domain.time
251
252        vol_id and edge_id are ignored
253        """
254
255        t = self.domain.time
256
257        if vol_id is not None and edge_id is not None:
258            i = self.boundary_indices[ vol_id, edge_id ]
259            return self.F(t, point_id = i)
260        else:
261            #raise 'Boundary call without point_id not implemented'
262            #FIXME: What should the semantics be?
263            return self.F(t)
264
265
266
267
268
269
270
271
272#THIS FAR (10/8/4)
273class Connective_boundary(Boundary):
274    """Connective boundary returns values for the
275    conserved quantities from a volume as defined by a connection table
276    mapping between tuples of (volume id, face id) for volumes that
277    have their boundaries connected.
278
279    FIXME: Perhaps include possibility of mapping between
280    different domains as well
281
282    FIXME: In case of shallow water we may want to have a
283    special version that casts this in terms of height rather than stage
284    """
285
286
287    def __init__(self, table):
288        from domain import Volume
289
290        Boundary.__init__(self)
291
292        self.connection_table = table
293        self.Volume = Volume
294
295    def __repr__(self):
296        return 'Connective boundary'
297
298    #FIXME: IF we ever need to get field_values from connected volume,
299    #that method could be overridden here (using same idea as in
300    #get_conserved_quantities
301    #def get_field_values()
302
303    def get_conserved_quantities(self, volume, face=0):
304
305        id = volume.id
306        if self.connection_table.has_key((id, face)):
307            other_id, other_face = self.connection_table[(id, face)]
308
309            other_volume = self.Volume.instances[other_id]
310            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
311            exec(cmd)
312            return q
313        else:
314            msg = 'Volume, face tuple (\$d, %d) has not been mapped'\
315                  %(id, face)
316            raise msg
317
318
319
320
321
322#FIXME: Add a boundary with a general function of x,y, and t
323
325# Attempt at periodic conditions from advection_spik. Remember this
326#
327#first = 2*(N-1)*N
328#for i in range(1,2*N+1,2):
329#    k = first + i-1#
330#
331#    print i,k
332#
333#    domain[i].faces[2].neighbour = domain[k].faces[1]
334#    domain[k].faces[1].neighbour = domain[i].faces[2]
335
336
337
338class General_boundary(Boundary):
339    """General boundary which can compute conserved quantities based on
340    their previous value, conserved quantities of its neighbour and model time.
341
342    Must specify initial conserved quantities,
343    neighbour,
345    a function f(q_old, neighbours_q, t) which must return
346    new conserved quantities q as a function time
347
348    FIXME: COMPLETE UNTESTED - JUST AN IDEA
349    """
350
351    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
352        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)
353
354        self.f = f
355        self.domain = domain
356
357
358    def get_conserved_quantities(self, volume=None, face=0):
359        return self.f(self.conserved_quantities,
360                      neighbour.conserved_quantities,
361                      self.domain.time)
362
363
364
365
366
367class File_boundary_old(Boundary):
368    """Boundary values obtained from file and interpolated.
369    conserved quantities as a function of time.
370
371    Assumes that file contains a time series.
372
373    No spatial info assumed.
374    """
375
376
377    def __init__(self, domain = None, filename = None):
378        import time
379        from Numeric import array
380        from config import time_format
381
382        Boundary.__init__(self)
383
384
385        try:
386            fid = open(filename)
387        except Exception, e:
388            msg = 'Boundary file %s could not be opened: %s\n' %(filename, e)
389            raise msg
390
391
393        fid.close()
394        fields = line.split(',')
395        msg = 'File %s must have the format date, values'
396        assert len(fields) == 2, msg
397
398        try:
399            starttime = time.mktime(time.strptime(fields[0], time_format))
400        except ValueError:
401            msg = 'First field in file %s must be' %filename
402            msg += ' date-time with format %s.\n' %time_format
403            msg += 'I got %s instead.' %fields[0]
404            raise msg
405
406        #Split values
407        values = []
408        for value in fields[1].split():
409            values.append(float(value))
410
411        q = array(values)
412
413        msg = 'ERROR: File boundary function must return a 1d list or array '
414        assert len(q.shape) == 1, msg
415
416        d = len(domain.conserved_quantities)
417        msg = 'Values specified in file must be a list or an array of length %d' %d
418        assert len(q) == d, msg
419
420        self.filename = filename
421        self.domain = domain
422        self.starttime = starttime
423
424        if domain.starttime is None:
425            domain.starttime = starttime
426        else:
427            msg = 'Start time as specified in domain (%s) is earlier '
428            'than the starttime of file %s: %s'\
429                  %(domain.starttime, self.filename, self.starttime)
430            if self.starttime > domain.starttime:
431                raise msg
432
434
435
437        from Numeric import zeros, Float, alltrue
438        from config import time_format
439        import time
440
441        fid = open(self.filename)
443        fid.close()
444
445        N = len(lines)
446        d = len(self.domain.conserved_quantities)
447
448        T = zeros(N, Float)
449        Q = zeros((N, d), Float)
450
451
452        for i, line in enumerate(lines):
453            fields = line.split(',')
454            real_time = time.mktime(time.strptime(fields[0], time_format))
455
456            T[i] = real_time - self.starttime
457
458
459            for j, value in enumerate(fields[1].split()):
460                Q[i, j] = float(value)
461
462        msg = 'Time boundary must list time as a monotonuosly '
463        msg += 'increasing sequence'
464
465        assert alltrue( T[1:] - T[:-1] > 0 ), msg
466
467        self.T = T     #Time
468        self.Q = Q     #Boundary values
469        self.index = 0 #Initial index
470
471
472    def __repr__(self):
473        return 'File boundary'
474
475    def evaluate(self, vol_id=None, edge_id=None):
476        """Return linearly interpolated values based on domain.time
477        """
478
479        t = self.domain.time
480
481        msg = 'Time given in File boundary does not match model time'
482        if t < self.T[0]: raise msg
483        if t > self.T[-1]: raise msg
484
485        oldindex = self.index
486
487        #Find slot
488        while t > self.T[self.index]: self.index += 1
489        while t < self.T[self.index]: self.index -= 1
490
491        #if oldindex != self.index:
492        #    print 'Changing from %d to %d' %(oldindex, self.index)
493
494
495        #t is now between index and index+1
496        ratio = (t - self.T[self.index])/\
497                (self.T[self.index+1] - self.T[self.index])
498
499        return self.Q[self.index,:] +\
500               ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
501
502
503
504
Note: See TracBrowser for help on using the repository browser.