1 | """boundary.py - Classes for implementing boundary conditions |
---|
2 | """ |
---|
3 | |
---|
4 | import numpy |
---|
5 | |
---|
6 | class 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 | |
---|
23 | class 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 | |
---|
51 | class 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 | |
---|
74 | class 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 | |
---|
116 | class 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 | |
---|
160 | class 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 |
---|
167 | admissible file formats |
---|
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) |
---|
269 | class 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 | |
---|
320 | #FIXME: Add periodic boundaries e.g.: |
---|
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 | |
---|
334 | class 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, |
---|
340 | domain to get access to model time |
---|
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 | |
---|
363 | class 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 | |
---|
387 | line = fid.readline() |
---|
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 | |
---|
428 | self.read_time_boundary() #Now read all times in. |
---|
429 | |
---|
430 | |
---|
431 | def read_time_boundary(self): |
---|
432 | from config import time_format |
---|
433 | import time |
---|
434 | |
---|
435 | fid = open(self.filename) |
---|
436 | lines = fid.readlines() |
---|
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 | |
---|