1 | |
---|
2 | from anuga.shallow_water.shallow_water_domain import * |
---|
3 | from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain |
---|
4 | |
---|
5 | |
---|
6 | ############################################################################## |
---|
7 | # Shallow Water Balanced Domain -- alternative implementation |
---|
8 | # |
---|
9 | # FIXME: Following the methods in CITE MODSIM PAPER |
---|
10 | # |
---|
11 | ############################################################################## |
---|
12 | |
---|
13 | class Domain(Sww_domain): |
---|
14 | |
---|
15 | def __init__(self, |
---|
16 | coordinates=None, |
---|
17 | vertices=None, |
---|
18 | boundary=None, |
---|
19 | tagged_elements=None, |
---|
20 | geo_reference=None, |
---|
21 | use_inscribed_circle=False, |
---|
22 | mesh_filename=None, |
---|
23 | use_cache=False, |
---|
24 | verbose=False, |
---|
25 | full_send_dict=None, |
---|
26 | ghost_recv_dict=None, |
---|
27 | starttime=0.0, |
---|
28 | processor=0, |
---|
29 | numproc=1, |
---|
30 | number_of_full_nodes=None, |
---|
31 | number_of_full_triangles=None): |
---|
32 | |
---|
33 | conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
---|
34 | |
---|
35 | evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum'] |
---|
36 | other_quantities = [ 'elevation', 'friction', 'height', |
---|
37 | 'xvelocity', 'yvelocity', 'x', 'y' ] |
---|
38 | |
---|
39 | |
---|
40 | Sww_domain.__init__(self, |
---|
41 | coordinates = coordinates, |
---|
42 | vertices = vertices, |
---|
43 | boundary = boundary, |
---|
44 | tagged_elements = tagged_elements, |
---|
45 | geo_reference = geo_reference, |
---|
46 | use_inscribed_circle = use_inscribed_circle, |
---|
47 | mesh_filename = mesh_filename, |
---|
48 | use_cache = use_cache, |
---|
49 | verbose = verbose, |
---|
50 | conserved_quantities = conserved_quantities, |
---|
51 | evolved_quantities = evolved_quantities, |
---|
52 | other_quantities = other_quantities, |
---|
53 | full_send_dict = full_send_dict, |
---|
54 | ghost_recv_dict = ghost_recv_dict, |
---|
55 | starttime = starttime, |
---|
56 | processor = processor, |
---|
57 | numproc = numproc, |
---|
58 | number_of_full_nodes = number_of_full_nodes, |
---|
59 | number_of_full_triangles = number_of_full_triangles) |
---|
60 | |
---|
61 | #------------------------------------------------ |
---|
62 | # set some defaults |
---|
63 | # Most of these override the options in config.py |
---|
64 | #------------------------------------------------ |
---|
65 | self.set_CFL(1.0) |
---|
66 | self.set_use_kinematic_viscosity(False) |
---|
67 | self.timestepping_method='rk2' |
---|
68 | # The following allows storage of the negative depths associated with this method |
---|
69 | self.minimum_storable_height=-99999999999.0 |
---|
70 | |
---|
71 | self.use_edge_limiter=True |
---|
72 | self.default_order=2 |
---|
73 | self.extrapolate_velocity_second_order=True |
---|
74 | |
---|
75 | # Note that the extrapolation method used in quantity_ext.c (e.g. |
---|
76 | # extrapolate_second_order_and_limit_by_edge) uses a constant value for |
---|
77 | # all the betas. |
---|
78 | self.beta_w=1.0 |
---|
79 | self.beta_w_dry=1.0 |
---|
80 | self.beta_uh=1.0 |
---|
81 | self.beta_uh_dry=1.0 |
---|
82 | self.beta_vh=1.0 |
---|
83 | self.beta_vh_dry=1.0 |
---|
84 | |
---|
85 | #self.optimise_dry_cells=True |
---|
86 | # "self.optimise_dry_cells=False" presently ensures that the stage is |
---|
87 | # always >= minimum_bed_edge value. Actually, we should only need to |
---|
88 | # apply 'False' on the very first time step (to deal with stage values |
---|
89 | # that were initialised below the bed by the user). After this, the |
---|
90 | # algorithm should take care of itself, and 'True' should be okay. |
---|
91 | self.optimise_dry_cells=False |
---|
92 | |
---|
93 | # Because gravity is treated within the flux function, |
---|
94 | # we remove it from the forcing terms. |
---|
95 | #self.forcing_terms.remove(gravity) |
---|
96 | |
---|
97 | # We need the edge_coordinates for the extrapolation |
---|
98 | self.edge_coordinates=self.get_edge_midpoint_coordinates() |
---|
99 | |
---|
100 | # We demand that vertex values are stored uniquely |
---|
101 | self.set_store_vertices_smoothly(False) |
---|
102 | |
---|
103 | self.maximum_allowed_speed=0.0 |
---|
104 | |
---|
105 | |
---|
106 | print '##########################################################################' |
---|
107 | print '#' |
---|
108 | print '# Using shallow_water_balanced2 solver in anuga_work/development/gareth/experimental/balanced_dev/' |
---|
109 | print "#" |
---|
110 | print "# This solver is experimental. Here are some tips on using it" |
---|
111 | print "#" |
---|
112 | print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values" |
---|
113 | print "# , as the latter can be completely misleading near strong gradients in the flow. " |
---|
114 | print "# There is a util.py script in anuga_work/development/gareth/tests which might help you extract" |
---|
115 | print "# quantities at centroid values from sww files." |
---|
116 | print "# Note that to accuractely compute centroid values from sww files, the files need to store " |
---|
117 | print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer" |
---|
118 | print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect" |
---|
119 | print "# ANUGA viewer as well -- I expect this would be lots of work)" |
---|
120 | print "#" |
---|
121 | print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set" |
---|
122 | print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). With this setting, velocity" |
---|
123 | print "# seems to be very well behaved" |
---|
124 | print "#" |
---|
125 | print "# 3) This solver is not expected to perform well in problems with very" |
---|
126 | print "# shallow water flowing down steep slopes (such that the stage_centroid_value " |
---|
127 | print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests" |
---|
128 | print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) " |
---|
129 | print "#" |
---|
130 | print "# 4) This solver allows the stage_centroid_value to drop to the minimum bed_edge_value" |
---|
131 | print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the" |
---|
132 | print "# bed_centroid_value. This means that triangles store slightly more water than they are" |
---|
133 | print "# typically interpreted to store, which might have significance in some applications." |
---|
134 | print "#" |
---|
135 | print "# You will probably be able to tell this is causing you problems by convergence testing" |
---|
136 | print "#" |
---|
137 | print '# 5) Note that many options in config.py have been overridden by the solver -- I have ' |
---|
138 | print '# deliberately attempted to get the solver to perform well with consistent values of ' |
---|
139 | print '# these parameters -- so I would advise against changing them unless you at least check that ' |
---|
140 | print '# it really does improve things' |
---|
141 | print '#' |
---|
142 | print '##########################################################################' |
---|
143 | |
---|
144 | #------------------------------- |
---|
145 | # Specialisation of Evolve |
---|
146 | #------------------------------- |
---|
147 | # This alters the 'evolve' method from shallow_water_domain.py so that just |
---|
148 | # before the file-output step, the vertex values are replaced with the |
---|
149 | # centroid values. |
---|
150 | # This is useful so that we can unambigously know what the centroid values |
---|
151 | # of various parameters are |
---|
152 | |
---|
153 | def evolve(self, |
---|
154 | yieldstep=None, |
---|
155 | finaltime=None, |
---|
156 | duration=None, |
---|
157 | skip_initial_step=False): |
---|
158 | """Specialisation of basic evolve method from parent class. |
---|
159 | |
---|
160 | Evolve the model by 1 step. |
---|
161 | """ |
---|
162 | |
---|
163 | # Call check integrity here rather than from user scripts |
---|
164 | # self.check_integrity() |
---|
165 | |
---|
166 | #print 'Into Evolve' |
---|
167 | |
---|
168 | msg = 'Attribute self.beta_w must be in the interval [0, 2]' |
---|
169 | assert 0 <= self.beta_w <= 2.0, msg |
---|
170 | |
---|
171 | # Initial update of vertex and edge values before any STORAGE |
---|
172 | # and or visualisation. |
---|
173 | # This is done again in the initialisation of the Generic_Domain |
---|
174 | # evolve loop but we do it here to ensure the values are ok for storage. |
---|
175 | self.distribute_to_vertices_and_edges() |
---|
176 | |
---|
177 | if self.store is True and self.get_time() == self.get_starttime(): |
---|
178 | self.initialise_storage() |
---|
179 | #print 'Into Generic_Domain Evolve' |
---|
180 | # Call basic machinery from parent class |
---|
181 | for t in Generic_Domain.evolve(self, yieldstep=yieldstep, |
---|
182 | finaltime=finaltime, duration=duration, |
---|
183 | skip_initial_step=skip_initial_step): |
---|
184 | #print 'Out of Generic_Domain Evolve' |
---|
185 | # Store model data, e.g. for subsequent visualisation |
---|
186 | if self.store is True: |
---|
187 | # FIXME: Extrapolation is done before writing, because I had |
---|
188 | # trouble correctly computing the centroid values from the |
---|
189 | # vertex values (an 'average' did not seem to work correctly in |
---|
190 | # very shallow cells -- there was a discrepency between |
---|
191 | # domain.quantity['blah'].centroid_values and the values |
---|
192 | # computed from the sww using the vertex averge). There should |
---|
193 | # be a much much more disk-efficient way to store the centroid |
---|
194 | # values than this |
---|
195 | self.extrapolate_second_order_edge_sw() |
---|
196 | |
---|
197 | # Store the timestep |
---|
198 | self.store_timestep() |
---|
199 | |
---|
200 | # Pass control on to outer loop for more specific actions |
---|
201 | yield(t) |
---|
202 | |
---|
203 | #----------------- |
---|
204 | # Flux computation |
---|
205 | #----------------- |
---|
206 | |
---|
207 | ## @brief Compute fluxes and timestep suitable for all volumes in domain. |
---|
208 | # @param domain The domain to calculate fluxes for. |
---|
209 | def compute_fluxes(domain): |
---|
210 | """Compute fluxes and timestep suitable for all volumes in domain. |
---|
211 | |
---|
212 | Compute total flux for each conserved quantity using "flux_function" |
---|
213 | |
---|
214 | Fluxes across each edge are scaled by edgelengths and summed up |
---|
215 | Resulting flux is then scaled by area and stored in |
---|
216 | explicit_update for each of the three conserved quantities |
---|
217 | stage, xmomentum and ymomentum |
---|
218 | |
---|
219 | The maximal allowable speed computed by the flux_function for each volume |
---|
220 | is converted to a timestep that must not be exceeded. The minimum of |
---|
221 | those is computed as the next overall timestep. |
---|
222 | |
---|
223 | Post conditions: |
---|
224 | domain.explicit_update is reset to computed flux values |
---|
225 | domain.timestep is set to the largest step satisfying all volumes. |
---|
226 | |
---|
227 | This wrapper calls the underlying C version of compute fluxes |
---|
228 | """ |
---|
229 | |
---|
230 | import sys |
---|
231 | from swb2_domain_ext import compute_fluxes_ext_central \ |
---|
232 | as compute_fluxes_ext |
---|
233 | |
---|
234 | #print "." |
---|
235 | |
---|
236 | # Shortcuts |
---|
237 | Stage = domain.quantities['stage'] |
---|
238 | Xmom = domain.quantities['xmomentum'] |
---|
239 | Ymom = domain.quantities['ymomentum'] |
---|
240 | Bed = domain.quantities['elevation'] |
---|
241 | |
---|
242 | timestep = float(sys.maxint) |
---|
243 | |
---|
244 | flux_timestep = compute_fluxes_ext(timestep, |
---|
245 | domain.epsilon, |
---|
246 | domain.H0, |
---|
247 | domain.g, |
---|
248 | domain.neighbours, |
---|
249 | domain.neighbour_edges, |
---|
250 | domain.normals, |
---|
251 | domain.edgelengths, |
---|
252 | domain.radii, |
---|
253 | domain.areas, |
---|
254 | domain.tri_full_flag, |
---|
255 | Stage.edge_values, |
---|
256 | Xmom.edge_values, |
---|
257 | Ymom.edge_values, |
---|
258 | Bed.edge_values, |
---|
259 | Stage.boundary_values, |
---|
260 | Xmom.boundary_values, |
---|
261 | Ymom.boundary_values, |
---|
262 | domain.boundary_flux_type, |
---|
263 | Stage.explicit_update, |
---|
264 | Xmom.explicit_update, |
---|
265 | Ymom.explicit_update, |
---|
266 | domain.already_computed_flux, |
---|
267 | domain.max_speed, |
---|
268 | int(domain.optimise_dry_cells), |
---|
269 | Stage.centroid_values, |
---|
270 | Bed.centroid_values, |
---|
271 | Bed.vertex_values) |
---|
272 | |
---|
273 | #import pdb |
---|
274 | #pdb.set_trace() |
---|
275 | |
---|
276 | domain.flux_timestep = flux_timestep |
---|
277 | |
---|
278 | |
---|
279 | def protect_against_infinitesimal_and_negative_heights(domain): |
---|
280 | """protect against infinitesimal heights and associated high velocities""" |
---|
281 | |
---|
282 | from swb2_domain_ext import protect |
---|
283 | #print'using swb2_protect_against ..' |
---|
284 | |
---|
285 | # shortcuts |
---|
286 | wc = domain.quantities['stage'].centroid_values |
---|
287 | wv = domain.quantities['stage'].vertex_values |
---|
288 | zc = domain.quantities['elevation'].centroid_values |
---|
289 | zv = domain.quantities['elevation'].vertex_values |
---|
290 | xmomc = domain.quantities['xmomentum'].centroid_values |
---|
291 | ymomc = domain.quantities['ymomentum'].centroid_values |
---|
292 | |
---|
293 | protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, |
---|
294 | domain.epsilon, wc, wv, zc,zv, xmomc, ymomc) |
---|
295 | |
---|
296 | def conserved_values_to_evolved_values(self, q_cons, q_evol): |
---|
297 | """Mapping between conserved quantities and the evolved quantities. |
---|
298 | Used where we have a boundary condition which works with conserved |
---|
299 | quantities and we now want to use them for the new well balanced |
---|
300 | code using the evolved quantities |
---|
301 | |
---|
302 | Typically the old boundary condition will set the values in q_cons, |
---|
303 | |
---|
304 | q_evol on input will have the values of the evolved quantities at the |
---|
305 | edge point (useful as it provides values for evlevation). |
---|
306 | |
---|
307 | """ |
---|
308 | |
---|
309 | wc = q_cons[0] |
---|
310 | uhc = q_cons[1] |
---|
311 | vhc = q_cons[2] |
---|
312 | |
---|
313 | |
---|
314 | q_evol[0] = wc |
---|
315 | q_evol[1] = uhc |
---|
316 | q_evol[2] = vhc |
---|
317 | |
---|
318 | return q_evol |
---|
319 | |
---|
320 | def distribute_to_vertices_and_edges(self): |
---|
321 | """ Call correct module function """ |
---|
322 | |
---|
323 | if self.use_edge_limiter: |
---|
324 | #distribute_using_edge_limiter(self) |
---|
325 | self.distribute_using_edge_limiter() |
---|
326 | else: |
---|
327 | #distribute_using_vertex_limiter(self) |
---|
328 | self.distribute_using_vertex_limiter() |
---|
329 | |
---|
330 | |
---|
331 | def distribute_using_vertex_limiter(domain): |
---|
332 | """Distribution from centroids to vertices specific to the SWW equation. |
---|
333 | |
---|
334 | Precondition: |
---|
335 | All quantities defined at centroids and bed elevation defined at |
---|
336 | vertices. |
---|
337 | |
---|
338 | Postcondition |
---|
339 | Conserved quantities defined at vertices |
---|
340 | """ |
---|
341 | # Remove very thin layers of water |
---|
342 | domain.protect_against_infinitesimal_and_negative_heights() |
---|
343 | |
---|
344 | # Extrapolate all conserved quantities |
---|
345 | if domain.optimised_gradient_limiter: |
---|
346 | # MH090605 if second order, |
---|
347 | # perform the extrapolation and limiting on |
---|
348 | # all of the conserved quantities |
---|
349 | |
---|
350 | if (domain._order_ == 1): |
---|
351 | for name in domain.conserved_quantities: |
---|
352 | Q = domain.quantities[name] |
---|
353 | Q.extrapolate_first_order() |
---|
354 | elif domain._order_ == 2: |
---|
355 | domain.extrapolate_second_order_sw() |
---|
356 | else: |
---|
357 | raise Exception('Unknown order') |
---|
358 | else: |
---|
359 | # Old code: |
---|
360 | for name in domain.conserved_quantities: |
---|
361 | Q = domain.quantities[name] |
---|
362 | |
---|
363 | if domain._order_ == 1: |
---|
364 | Q.extrapolate_first_order() |
---|
365 | elif domain._order_ == 2: |
---|
366 | Q.extrapolate_second_order_and_limit_by_vertex() |
---|
367 | else: |
---|
368 | raise Exception('Unknown order') |
---|
369 | |
---|
370 | # Take bed elevation into account when water heights are small |
---|
371 | #balance_deep_and_shallow(domain) |
---|
372 | |
---|
373 | # Compute edge values by interpolation |
---|
374 | for name in domain.conserved_quantities: |
---|
375 | Q = domain.quantities[name] |
---|
376 | Q.interpolate_from_vertices_to_edges() |
---|
377 | |
---|
378 | |
---|
379 | def distribute_using_edge_limiter(domain): |
---|
380 | """Distribution from centroids to edges specific to the SWW eqn. |
---|
381 | |
---|
382 | Precondition: |
---|
383 | All quantities defined at centroids and bed elevation defined at |
---|
384 | vertices. |
---|
385 | |
---|
386 | Postcondition |
---|
387 | Conserved quantities defined at vertices |
---|
388 | """ |
---|
389 | #from swb2_domain import protect_against_infinitesimal_and_negative_heights |
---|
390 | |
---|
391 | # Remove very thin layers of water |
---|
392 | #print 'Before protect' |
---|
393 | domain.protect_against_infinitesimal_and_negative_heights() |
---|
394 | #print 'After protect' |
---|
395 | |
---|
396 | #for name in domain.conserved_quantities: |
---|
397 | # Q = domain.quantities[name] |
---|
398 | # if domain._order_ == 1: |
---|
399 | # Q.extrapolate_first_order() |
---|
400 | # elif domain._order_ == 2: |
---|
401 | # #Q.extrapolate_second_order_and_limit_by_edge() |
---|
402 | # # FIXME: This use of 'break' is hacky |
---|
403 | # domain.extrapolate_second_order_edge_sw() |
---|
404 | # break |
---|
405 | # else: |
---|
406 | # raise Exception('Unknown order') |
---|
407 | |
---|
408 | #print 'Before extrapolate' |
---|
409 | domain.extrapolate_second_order_edge_sw() |
---|
410 | #print 'After extrapolate' |
---|
411 | |
---|
412 | #balance_deep_and_shallow(domain) |
---|
413 | |
---|
414 | # Compute vertex values by interpolation |
---|
415 | #for name in domain.conserved_quantities: |
---|
416 | # Q = domain.quantities[name] |
---|
417 | # Q.interpolate_from_edges_to_vertices() |
---|
418 | # #Q.interpolate_from_vertices_to_edges() |
---|
419 | |
---|
420 | |
---|
421 | def update_other_quantities(self): |
---|
422 | """ There may be a need to calculates some of the other quantities |
---|
423 | based on the new values of conserved quantities |
---|
424 | |
---|
425 | But that is not needed in this version of the solver. |
---|
426 | """ |
---|
427 | pass |
---|
428 | # The centroid values of height and x and y velocity |
---|
429 | # might not have been setup |
---|
430 | |
---|
431 | #self.update_centroids_of_velocities_and_height() |
---|
432 | # |
---|
433 | #for name in ['height', 'xvelocity', 'yvelocity']: |
---|
434 | # Q = self.quantities[name] |
---|
435 | # if self._order_ == 1: |
---|
436 | # Q.extrapolate_first_order() |
---|
437 | # elif self._order_ == 2: |
---|
438 | # Q.extrapolate_second_order_and_limit_by_edge() |
---|
439 | # else: |
---|
440 | # raise Exception('Unknown order') |
---|
441 | |
---|
442 | |
---|
443 | def update_centroids_of_velocities_and_height(self): |
---|
444 | """Calculate the centroid values of velocities and height based |
---|
445 | on the values of the quantities stage and x and y momentum |
---|
446 | |
---|
447 | Assumes that stage and momentum are up to date |
---|
448 | |
---|
449 | Useful for kinematic viscosity calculations |
---|
450 | """ |
---|
451 | pass |
---|
452 | ## For shallow water we need to update height xvelocity and yvelocity |
---|
453 | |
---|
454 | ##Shortcuts |
---|
455 | #W = self.quantities['stage'] |
---|
456 | #UH = self.quantities['xmomentum'] |
---|
457 | #VH = self.quantities['ymomentum'] |
---|
458 | #H = self.quantities['height'] |
---|
459 | #Z = self.quantities['elevation'] |
---|
460 | #U = self.quantities['xvelocity'] |
---|
461 | #V = self.quantities['yvelocity'] |
---|
462 | |
---|
463 | ##print num.min(W.centroid_values) |
---|
464 | |
---|
465 | ## Make sure boundary values of conserved quantites |
---|
466 | ## are consistent with value of functions at centroids |
---|
467 | ##self.distribute_to_vertices_and_edges() |
---|
468 | #Z.set_boundary_values_from_edges() |
---|
469 | |
---|
470 | ##W.set_boundary_values_from_edges() |
---|
471 | ##UH.set_boundary_values_from_edges() |
---|
472 | ##VH.set_boundary_values_from_edges() |
---|
473 | |
---|
474 | ## Update height values |
---|
475 | #H.set_values(W.centroid_values-Z.centroid_values, location='centroids') |
---|
476 | #H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0, |
---|
477 | # W.boundary_values-Z.boundary_values, 0.0)) |
---|
478 | |
---|
479 | ##assert num.min(H.centroid_values) >= 0 |
---|
480 | ##assert num.min(H.boundary_values) >= 0 |
---|
481 | |
---|
482 | ##Aliases |
---|
483 | #uh_C = UH.centroid_values |
---|
484 | #vh_C = VH.centroid_values |
---|
485 | #h_C = H.centroid_values |
---|
486 | |
---|
487 | #uh_B = UH.boundary_values |
---|
488 | #vh_B = VH.boundary_values |
---|
489 | #h_B = H.boundary_values |
---|
490 | |
---|
491 | #H0 = 1.0e-8 |
---|
492 | # |
---|
493 | #U.set_values(uh_C/(h_C + H0/h_C), location='centroids') |
---|
494 | #V.set_values(vh_C/(h_C + H0/h_C), location='centroids') |
---|
495 | |
---|
496 | #U.set_boundary_values(uh_B/(h_B + H0/h_B)) |
---|
497 | #V.set_boundary_values(vh_B/(h_B + H0/h_B)) |
---|
498 | |
---|
499 | |
---|
500 | |
---|
501 | def extrapolate_second_order_edge_sw(self): |
---|
502 | """Wrapper calling C version of extrapolate_second_order_sw, using |
---|
503 | edges instead of vertices in the extrapolation step. The routine |
---|
504 | then interpolates from edges to vertices. |
---|
505 | The momentum extrapolation / interpolation is based on either |
---|
506 | momentum or velocity, depending on the choice of |
---|
507 | extrapolate_velocity_second_order. |
---|
508 | |
---|
509 | self the domain to operate on |
---|
510 | |
---|
511 | """ |
---|
512 | |
---|
513 | from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2 |
---|
514 | |
---|
515 | # Shortcuts |
---|
516 | Stage = self.quantities['stage'] |
---|
517 | Xmom = self.quantities['xmomentum'] |
---|
518 | Ymom = self.quantities['ymomentum'] |
---|
519 | Elevation = self.quantities['elevation'] |
---|
520 | |
---|
521 | extrapol2(self, |
---|
522 | self.surrogate_neighbours, |
---|
523 | self.number_of_boundaries, |
---|
524 | self.centroid_coordinates, |
---|
525 | Stage.centroid_values, |
---|
526 | Xmom.centroid_values, |
---|
527 | Ymom.centroid_values, |
---|
528 | Elevation.centroid_values, |
---|
529 | self.edge_coordinates, |
---|
530 | Stage.edge_values, |
---|
531 | Xmom.edge_values, |
---|
532 | Ymom.edge_values, |
---|
533 | Elevation.edge_values, |
---|
534 | Stage.vertex_values, |
---|
535 | Xmom.vertex_values, |
---|
536 | Ymom.vertex_values, |
---|
537 | Elevation.vertex_values, |
---|
538 | int(self.optimise_dry_cells), |
---|
539 | int(self.extrapolate_velocity_second_order)) |
---|
540 | |
---|
541 | def set_boundary(self, boundary_map): |
---|
542 | ## Specialisation of set_boundary, which also updates the 'boundary_flux_type' |
---|
543 | Sww_domain.set_boundary(self,boundary_map) |
---|
544 | |
---|
545 | # Add a flag which can be used to distinguish flux boundaries within |
---|
546 | # compute_fluxes_central |
---|
547 | # Initialise to zero (which means 'not a flux_boundary') |
---|
548 | self.boundary_flux_type = self.boundary_edges*0 |
---|
549 | |
---|
550 | # HACK to set the values of domain.boundary_flux |
---|
551 | for i in range(len(self.boundary_objects)): |
---|
552 | # Record the first 10 characters of the name of the boundary. |
---|
553 | # FIXME: There must be a better way than this! |
---|
554 | bndry_name = self.boundary_objects[i][1].__repr__()[0:10] |
---|
555 | |
---|
556 | if(bndry_name=='zero_mass_'): |
---|
557 | # Create flag 'boundary_flux_type' that can be read in |
---|
558 | # compute_fluxes_central |
---|
559 | self.boundary_flux_type[i]=1 |
---|
560 | |
---|
561 | |
---|