# source:trunk/anuga_work/development/2010-projects/anuga_1d/base/generic_boundary_conditions.py@7884

Last change on this file since 7884 was 7884, checked in by steve, 13 years ago

Moving 2010 project

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