source: inundation/pyvolution/generic_boundary_conditions.py @ 1740

Last change on this file since 1740 was 1704, checked in by ole, 19 years ago

Comments

File size: 14.1 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
170    admissible file formats
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        #Record ordering #FIXME: should this happen in domain.py
200        self.boundary_indices = {}
201
202        for i, (vol_id, edge_id) in enumerate(boundary_keys):
203
204            x0 = V[vol_id, 0]; y0 = V[vol_id, 1]
205            x1 = V[vol_id, 2]; y1 = V[vol_id, 3]
206            x2 = V[vol_id, 4]; y2 = V[vol_id, 5]
207
208            #Midpoints
209            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
210            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
211            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
212            self.midpoint_coordinates[i,:] = m
213
214            self.boundary_indices[(vol_id, edge_id)] = i
215
216
217        if verbose: print 'Initialise file_function'
218        self.F = file_function(filename, domain,
219                               interpolation_points=self.midpoint_coordinates,
220                               verbose = verbose)
221        self.domain = domain
222
223        #Test
224        q = self.F(0, point_id=0)
225
226        d = len(domain.conserved_quantities)
227        msg = 'Values specified in file %s must be ' %filename
228        msg += ' a list or an array of length %d' %d
229        assert len(q) == d, msg
230
231
232    def __repr__(self):
233        return 'File boundary'
234
235
236    def evaluate(self, vol_id=None, edge_id=None):
237        """Return linearly interpolated values based on domain.time
238
239        vol_id and edge_id are ignored
240        """
241
242        t = self.domain.time
243
244        if vol_id is not None and edge_id is not None:
245            i = self.boundary_indices[ vol_id, edge_id ]
246            return self.F(t, point_id = i)
247        else:
248            #raise 'Boundary call without point_id not implemented'
249            #FIXME: What should the semantics be?
250            return self.F(t)
251
252
253
254
255
256
257
258
259#THIS FAR (10/8/4)
260class Connective_boundary(Boundary):
261    """Connective boundary returns values for the
262    conserved quantities from a volume as defined by a connection table
263    mapping between tuples of (volume id, face id) for volumes that
264    have their boundaries connected.
265
266    FIXME: Perhaps include possibility of mapping between
267    different domains as well
268
269    FIXME: In case of shallow water we may want to have a
270    special version that casts this in terms of height rather than stage
271    """
272
273
274    def __init__(self, table):
275        from domain import Volume
276
277        Boundary.__init__(self)
278
279        self.connection_table = table
280        self.Volume = Volume
281
282    def __repr__(self):
283        return 'Connective boundary'
284
285    #FIXME: IF we ever need to get field_values from connected volume,
286    #that method could be overridden here (using same idea as in
287    #get_conserved_quantities
288    #def get_field_values()
289
290    def get_conserved_quantities(self, volume, face=0):
291
292        id = volume.id
293        if self.connection_table.has_key((id, face)):
294            other_id, other_face = self.connection_table[(id, face)]
295
296            other_volume = self.Volume.instances[other_id]
297            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
298            exec(cmd)
299            return q
300        else:
301            msg = 'Volume, face tuple ($d, %d) has not been mapped'\
302                  %(id, face)
303            raise msg
304
305
306
307
308
309#FIXME: Add a boundary with a general function of x,y, and t
310
311#FIXME: Add periodic boundaries e.g.:
312# Attempt at periodic conditions from advection_spik. Remember this
313#
314#first = 2*(N-1)*N
315#for i in range(1,2*N+1,2):
316#    k = first + i-1#
317#
318#    print i,k
319#
320#    domain[i].faces[2].neighbour = domain[k].faces[1]
321#    domain[k].faces[1].neighbour = domain[i].faces[2]
322
323
324
325class General_boundary(Boundary):
326    """General boundary which can compute conserved quantities based on
327    their previous value, conserved quantities of its neighbour and model time.
328
329    Must specify initial conserved quantities,
330    neighbour,
331    domain to get access to model time
332    a function f(q_old, neighbours_q, t) which must return
333    new conserved quantities q as a function time
334
335    FIXME: COMPLETE UNTESTED - JUST AN IDEA
336    """
337
338    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
339        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)
340
341        self.f = f
342        self.domain = domain
343
344
345    def get_conserved_quantities(self, volume=None, face=0):
346        return self.f(self.conserved_quantities,
347                      neighbour.conserved_quantities,
348                      self.domain.time)
349
350
351
352
353
354class File_boundary_old(Boundary):
355    """Boundary values obtained from file and interpolated.
356    conserved quantities as a function of time.
357
358    Assumes that file contains a time series.
359
360    No spatial info assumed.
361    """
362
363
364    def __init__(self, domain = None, filename = None):
365        import time
366        from Numeric import array
367        from config import time_format
368
369        Boundary.__init__(self)
370
371
372        try:
373            fid = open(filename)
374        except Exception, e:
375            msg = 'Boundary file %s could not be opened: %s\n' %(filename, e)
376            raise msg
377
378
379        line = fid.readline()
380        fid.close()
381        fields = line.split(',')
382        msg = 'File %s must have the format date, values'
383        assert len(fields) == 2, msg
384
385        try:
386            starttime = time.mktime(time.strptime(fields[0], time_format))
387        except ValueError:
388            msg = 'First field in file %s must be' %filename
389            msg += ' date-time with format %s.\n' %time_format
390            msg += 'I got %s instead.' %fields[0]
391            raise msg
392
393        #Split values
394        values = []
395        for value in fields[1].split():
396            values.append(float(value))
397
398        q = array(values)
399
400        msg = 'ERROR: File boundary function must return a 1d list or array '
401        assert len(q.shape) == 1, msg
402
403        d = len(domain.conserved_quantities)
404        msg = 'Values specified in file must be a list or an array of length %d' %d
405        assert len(q) == d, msg
406
407        self.filename = filename
408        self.domain = domain
409        self.starttime = starttime
410
411        if domain.starttime is None:
412            domain.starttime = starttime
413        else:
414            msg = 'Start time as specified in domain (%s) is earlier '
415            'than the starttime of file %s: %s'\
416                  %(domain.starttime, self.filename, self.starttime)
417            if self.starttime > domain.starttime:
418                raise msg
419
420        self.read_time_boundary() #Now read all times in.
421
422
423    def read_time_boundary(self):
424        from Numeric import zeros, Float, alltrue
425        from config import time_format
426        import time
427
428        fid = open(self.filename)
429        lines = fid.readlines()
430        fid.close()
431
432        N = len(lines)
433        d = len(self.domain.conserved_quantities)
434
435        T = zeros(N, Float)
436        Q = zeros((N, d), Float)
437
438
439        for i, line in enumerate(lines):
440            fields = line.split(',')
441            real_time = time.mktime(time.strptime(fields[0], time_format))
442
443            T[i] = real_time - self.starttime
444
445
446            for j, value in enumerate(fields[1].split()):
447                Q[i, j] = float(value)
448
449        msg = 'Time boundary must list time as a monotonuosly '
450        msg += 'increasing sequence'
451
452        assert alltrue( T[1:] - T[:-1] > 0 ), msg
453
454        self.T = T     #Time
455        self.Q = Q     #Boundary values
456        self.index = 0 #Initial index
457
458
459    def __repr__(self):
460        return 'File boundary'
461
462    def evaluate(self, vol_id=None, edge_id=None):
463        """Return linearly interpolated values based on domain.time
464        """
465
466        t = self.domain.time
467
468        msg = 'Time given in File boundary does not match model time'
469        if t < self.T[0]: raise msg
470        if t > self.T[-1]: raise msg
471
472        oldindex = self.index
473
474        #Find slot
475        while t > self.T[self.index]: self.index += 1
476        while t < self.T[self.index]: self.index -= 1
477
478        #if oldindex != self.index:
479        #    print 'Changing from %d to %d' %(oldindex, self.index)
480
481
482        #t is now between index and index+1
483        ratio = (t - self.T[self.index])/\
484                (self.T[self.index+1] - self.T[self.index])
485
486        return self.Q[self.index,:] +\
487               ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
488
489
490
491
Note: See TracBrowser for help on using the repository browser.