source: anuga_work/development/sudi/sw_1d/avalanche/B_momentum/generic_boundary_conditions.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

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