1 | from anuga_1d.sqpipe.sqpipe_domain import Sqpipe_domain |
---|
2 | from anuga_1d.sqpipe.sqpipe_forcing_terms import * |
---|
3 | from anuga_1d.sqpipe.sqpipe_boundary_conditions import * |
---|
4 | |
---|
5 | |
---|
6 | class Domain(Sqpipe_domain): |
---|
7 | |
---|
8 | def __init__(self, coordinates, boundary = None, tagged_elements = None, bulk_modulus = 1.0): |
---|
9 | conserved_quantities = ['area', 'discharge'] |
---|
10 | forcing_terms = [] |
---|
11 | state = numpy.zeros(len(coordinates)-1, numpy.int) |
---|
12 | Sqpipe_domain.__init__(self, |
---|
13 | coordinates = coordinates, |
---|
14 | conserved_quantities = conserved_quantities, |
---|
15 | forcing_terms = forcing_terms, |
---|
16 | boundary = boundary, |
---|
17 | state = state, |
---|
18 | bulk_modulus = bulk_modulus, |
---|
19 | # update_state = False, |
---|
20 | tagged_elements = tagged_elements) |
---|
21 | |
---|
22 | |
---|
23 | self.__doc__ = 'sqpipe_area_discharge_domain' |
---|
24 | |
---|
25 | def distribute_to_vertices_and_edges(self): |
---|
26 | # Shortcuts |
---|
27 | h0 = self.h0 |
---|
28 | epsilon = self.epsilon |
---|
29 | |
---|
30 | area = self.quantities['area'] |
---|
31 | discharge = self.quantities['discharge'] |
---|
32 | bed = self.quantities['elevation'] |
---|
33 | height = self.quantities['height'] |
---|
34 | velocity = self.quantities['velocity'] |
---|
35 | width = self.quantities['width'] |
---|
36 | top = self.quantities['top'] |
---|
37 | stage = self.quantities['stage'] |
---|
38 | |
---|
39 | #Arrays |
---|
40 | a_C = area.centroid_values |
---|
41 | d_C = discharge.centroid_values |
---|
42 | z_C = bed.centroid_values |
---|
43 | h_C = height.centroid_values |
---|
44 | u_C = velocity.centroid_values |
---|
45 | b_C = width.centroid_values |
---|
46 | t_C = top.centroid_values |
---|
47 | w_C = stage.centroid_values |
---|
48 | |
---|
49 | #Calculate height (and fix negatives)better be non-negative! |
---|
50 | h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon? |
---|
51 | w_C[:] = z_C + h_C |
---|
52 | u_C[:] = d_C/(h_C + h0/(h_C + h0)) |
---|
53 | |
---|
54 | for name in ['velocity', 'stage' ]: |
---|
55 | Q = self.quantities[name] |
---|
56 | if self.order == 1: |
---|
57 | Q.extrapolate_first_order() |
---|
58 | elif self.order == 2: |
---|
59 | Q.extrapolate_second_order() |
---|
60 | else: |
---|
61 | raise 'Unknown order' |
---|
62 | |
---|
63 | # These have been extrapolated |
---|
64 | w_V = self.quantities['stage'].vertex_values |
---|
65 | u_V = self.quantities['velocity'].vertex_values |
---|
66 | |
---|
67 | # These are the given geometry and remain fixed |
---|
68 | z_V = bed.vertex_values |
---|
69 | b_V = width.vertex_values |
---|
70 | t_V = top.vertex_values |
---|
71 | |
---|
72 | # These need to be updated |
---|
73 | a_V = area.vertex_values |
---|
74 | d_V = discharge.vertex_values |
---|
75 | h_V = height.vertex_values |
---|
76 | |
---|
77 | h_V[:,:] = w_V - z_V |
---|
78 | |
---|
79 | # Fix up small heights |
---|
80 | h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0]) |
---|
81 | h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1]) |
---|
82 | |
---|
83 | h_V[:,0] = h_0 |
---|
84 | h_V[:,1] = h_1 |
---|
85 | |
---|
86 | h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0]) |
---|
87 | h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1]) |
---|
88 | |
---|
89 | h_V[:,0] = h_0 |
---|
90 | h_V[:,1] = h_1 |
---|
91 | |
---|
92 | # It may still be possible that h is small |
---|
93 | # If we set h to zero, we should also set u to 0 |
---|
94 | h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V) |
---|
95 | u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V) |
---|
96 | |
---|
97 | # Reconstruct conserved quantities and make everything consistent |
---|
98 | # with new h values |
---|
99 | w_V[:,:] = z_V + h_V |
---|
100 | a_V[:,:] = h_V * b_V |
---|
101 | d_V[:,:] = u_V * a_V |
---|
102 | |
---|
103 | |
---|
104 | return |
---|
105 | |
---|
106 | |
---|