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