source: trunk/anuga_work/development/sudi/sw_1d/periodic_waves/johns/generic_boundary_conditions.py @ 7922

Last change on this file since 7922 was 7922, checked in by mungkasi, 14 years ago

Modifying codes so that arrays are compatible with numpy instead of Numeric.

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