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 | class 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 | |
---|
65 | class 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 | |
---|
93 | class 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 zeros |
---|
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 | |
---|
118 | class 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 | |
---|
165 | class 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 | |
---|
210 | class 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 zeros, array |
---|
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) |
---|
321 | class 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 | |
---|
386 | class 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 | |
---|
415 | class 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 | |
---|