| 1 | import os.path |
|---|
| 2 | import sys |
|---|
| 3 | |
|---|
| 4 | from anuga.utilities.system_tools import get_pathname_from_package |
|---|
| 5 | from anuga.geometry.polygon_function import Polygon_function |
|---|
| 6 | |
|---|
| 7 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
|---|
| 8 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
|---|
| 9 | from anuga.abstract_2d_finite_volumes.util import file_function |
|---|
| 10 | |
|---|
| 11 | import anuga |
|---|
| 12 | |
|---|
| 13 | #from anuga.structures.boyd_box_operator import Boyd_box_operator |
|---|
| 14 | #from anuga.structures.inlet_operator import Inlet_operator |
|---|
| 15 | |
|---|
| 16 | #from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model |
|---|
| 17 | |
|---|
| 18 | from math import pi, pow, sqrt |
|---|
| 19 | |
|---|
| 20 | import numpy as num |
|---|
| 21 | from parallel_inlet_operator import Parallel_Inlet_operator |
|---|
| 22 | from anuga_parallel import distribute, myid, numprocs, finalize |
|---|
| 23 | from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect |
|---|
| 24 | |
|---|
| 25 | from parallel_operator_factory import Inlet_operator, Boyd_box_operator |
|---|
| 26 | import pypar |
|---|
| 27 | import random |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | """test_that_culvert_runs_rating |
|---|
| 31 | |
|---|
| 32 | This test exercises the culvert and checks values outside rating curve |
|---|
| 33 | are dealt with |
|---|
| 34 | """ |
|---|
| 35 | verbose = True |
|---|
| 36 | path = get_pathname_from_package('anuga.culvert_flows') |
|---|
| 37 | |
|---|
| 38 | length = 40. |
|---|
| 39 | width = 15. |
|---|
| 40 | |
|---|
| 41 | dx = dy = 0.5 # Resolution: Length of subdivisions on both axes |
|---|
| 42 | |
|---|
| 43 | #---------------------------------------------------------------------- |
|---|
| 44 | # Setup initial conditions |
|---|
| 45 | #---------------------------------------------------------------------- |
|---|
| 46 | |
|---|
| 47 | def topography(x, y): |
|---|
| 48 | """Set up a weir |
|---|
| 49 | |
|---|
| 50 | A culvert will connect either side |
|---|
| 51 | """ |
|---|
| 52 | # General Slope of Topography |
|---|
| 53 | z=-x/1000 |
|---|
| 54 | |
|---|
| 55 | N = len(x) |
|---|
| 56 | for i in range(N): |
|---|
| 57 | |
|---|
| 58 | # Sloping Embankment Across Channel |
|---|
| 59 | if 5.0 < x[i] < 10.1: |
|---|
| 60 | # Cut Out Segment for Culvert face |
|---|
| 61 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: |
|---|
| 62 | z[i]=z[i] |
|---|
| 63 | else: |
|---|
| 64 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
|---|
| 65 | if 10.0 < x[i] < 12.1: |
|---|
| 66 | z[i] += 2.5 # Flat Crest of Embankment |
|---|
| 67 | if 12.0 < x[i] < 14.5: |
|---|
| 68 | # Cut Out Segment for Culvert face |
|---|
| 69 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: |
|---|
| 70 | z[i]=z[i] |
|---|
| 71 | else: |
|---|
| 72 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | return z |
|---|
| 76 | |
|---|
| 77 | filename=os.path.join(path, 'example_rating_curve.csv') |
|---|
| 78 | |
|---|
| 79 | line0 = [[10.0, 10.0], [30.0, 10.0]] |
|---|
| 80 | #line0 = [[29.0, 10.0], [30.0, 10.0]] |
|---|
| 81 | line1 = [[0.0, 10.0], [0.0, 15.0]] |
|---|
| 82 | Q0 = file_function('test_hydrograph.tms', quantities=['hydrograph']) |
|---|
| 83 | Q1 = 5.0 |
|---|
| 84 | |
|---|
| 85 | samples = 50 |
|---|
| 86 | |
|---|
| 87 | def run_test(parallel = False, control_data = None, test_points = None, verbose = False): |
|---|
| 88 | success = True |
|---|
| 89 | |
|---|
| 90 | ##----------------------------------------------------------------------- |
|---|
| 91 | ## Setup domain |
|---|
| 92 | ##----------------------------------------------------------------------- |
|---|
| 93 | |
|---|
| 94 | points, vertices, boundary = rectangular_cross(int(length/dx), |
|---|
| 95 | int(width/dy), |
|---|
| 96 | len1=length, |
|---|
| 97 | len2=width) |
|---|
| 98 | |
|---|
| 99 | domain = anuga.Domain(points, vertices, boundary) |
|---|
| 100 | domain.set_name('Test_Parallel_Frac_Op') # Output name |
|---|
| 101 | domain.set_default_order(2) |
|---|
| 102 | |
|---|
| 103 | ##----------------------------------------------------------------------- |
|---|
| 104 | ## Distribute domain |
|---|
| 105 | ##----------------------------------------------------------------------- |
|---|
| 106 | |
|---|
| 107 | if parallel: |
|---|
| 108 | domain = distribute(domain) |
|---|
| 109 | domain.dump_triangulation("frac_op_domain.png") |
|---|
| 110 | |
|---|
| 111 | |
|---|
| 112 | ##----------------------------------------------------------------------- |
|---|
| 113 | ## Setup boundary conditions |
|---|
| 114 | ##----------------------------------------------------------------------- |
|---|
| 115 | |
|---|
| 116 | domain.set_quantity('elevation', topography) |
|---|
| 117 | domain.set_quantity('friction', 0.01) # Constant friction |
|---|
| 118 | domain.set_quantity('stage', |
|---|
| 119 | expression='elevation') # Dry initial condition |
|---|
| 120 | |
|---|
| 121 | |
|---|
| 122 | Bi = anuga.Dirichlet_boundary([5.0, 0.0, 0.0]) |
|---|
| 123 | Br = anuga.Reflective_boundary(domain) # Solid reflective wall |
|---|
| 124 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
|---|
| 125 | |
|---|
| 126 | |
|---|
| 127 | ##----------------------------------------------------------------------- |
|---|
| 128 | ## Determine triangle index coinciding with test points |
|---|
| 129 | ##----------------------------------------------------------------------- |
|---|
| 130 | |
|---|
| 131 | assert(test_points is not None) |
|---|
| 132 | assert(len(test_points) == samples) |
|---|
| 133 | |
|---|
| 134 | tri_ids = [] |
|---|
| 135 | |
|---|
| 136 | for point in test_points: |
|---|
| 137 | try: |
|---|
| 138 | k = domain.get_triangle_containing_point(point) |
|---|
| 139 | if domain.tri_full_flag[k] == 1: |
|---|
| 140 | tri_ids.append(k) |
|---|
| 141 | else: |
|---|
| 142 | tri_ids.append(-1) |
|---|
| 143 | except: |
|---|
| 144 | tri_ids.append(-2) |
|---|
| 145 | |
|---|
| 146 | if verbose: print 'P%d has points = %s' %(myid, tri_ids) |
|---|
| 147 | |
|---|
| 148 | if not parallel: control_data = [] |
|---|
| 149 | |
|---|
| 150 | ################ Define Fractional Operators ########################## |
|---|
| 151 | |
|---|
| 152 | inlet0 = None |
|---|
| 153 | inlet1 = None |
|---|
| 154 | boyd_box0 = None |
|---|
| 155 | |
|---|
| 156 | inlet0 = Inlet_operator(domain, line0, Q0, debug = False) |
|---|
| 157 | inlet1 = Inlet_operator(domain, line1, Q1, debug = False) |
|---|
| 158 | |
|---|
| 159 | # Enquiry point [ 19. 2.5] is contained in two domains in 4 proc case |
|---|
| 160 | |
|---|
| 161 | boyd_box0 = Boyd_box_operator(domain, |
|---|
| 162 | end_points=[[9.0, 2.5],[19.0, 2.5]], |
|---|
| 163 | losses=1.5, |
|---|
| 164 | width=5.0, |
|---|
| 165 | apron=5.0, |
|---|
| 166 | use_momentum_jet=True, |
|---|
| 167 | use_velocity_head=False, |
|---|
| 168 | manning=0.013, |
|---|
| 169 | verbose=False, debug = False) |
|---|
| 170 | |
|---|
| 171 | if inlet0 is not None: inlet0.print_statistics() |
|---|
| 172 | if inlet1 is not None: inlet1.print_statistics() |
|---|
| 173 | if boyd_box0 is not None: boyd_box0.print_statistics() |
|---|
| 174 | |
|---|
| 175 | # if parallel: |
|---|
| 176 | # factory = Parallel_operator_factory(domain, debug = True) |
|---|
| 177 | # |
|---|
| 178 | # inlet0 = factory.inlet_operator_factory(line0, Q0) |
|---|
| 179 | # inlet1 = factory.inlet_operator_factory(line1, Q1) |
|---|
| 180 | # |
|---|
| 181 | # boyd_box0 = factory.boyd_box_operator_factory(end_points=[[9.0, 2.5],[19.0, 2.5]], |
|---|
| 182 | # losses=1.5, |
|---|
| 183 | # width=1.5, |
|---|
| 184 | # apron=5.0, |
|---|
| 185 | # use_momentum_jet=True, |
|---|
| 186 | # use_velocity_head=False, |
|---|
| 187 | # manning=0.013, |
|---|
| 188 | # verbose=False) |
|---|
| 189 | # |
|---|
| 190 | # else: |
|---|
| 191 | # inlet0 = Inlet_operator(domain, line0, Q0) |
|---|
| 192 | # inlet1 = Inlet_operator(domain, line1, Q1) |
|---|
| 193 | # |
|---|
| 194 | # # Enquiry point [ 19. 2.5] is contained in two domains in 4 proc case |
|---|
| 195 | # boyd_box0 = Boyd_box_operator(domain, |
|---|
| 196 | # end_points=[[9.0, 2.5],[19.0, 2.5]], |
|---|
| 197 | # losses=1.5, |
|---|
| 198 | # width=1.5, |
|---|
| 199 | # apron=5.0, |
|---|
| 200 | # use_momentum_jet=True, |
|---|
| 201 | # use_velocity_head=False, |
|---|
| 202 | # manning=0.013, |
|---|
| 203 | # verbose=False) |
|---|
| 204 | |
|---|
| 205 | ####################################################################### |
|---|
| 206 | |
|---|
| 207 | ##----------------------------------------------------------------------- |
|---|
| 208 | ## Evolve system through time |
|---|
| 209 | ##----------------------------------------------------------------------- |
|---|
| 210 | |
|---|
| 211 | for t in domain.evolve(yieldstep = 0.1, finaltime = 38): |
|---|
| 212 | domain.write_time() |
|---|
| 213 | |
|---|
| 214 | #print domain.volumetric_balance_statistics() |
|---|
| 215 | |
|---|
| 216 | stage = domain.get_quantity('stage') |
|---|
| 217 | |
|---|
| 218 | |
|---|
| 219 | if boyd_box0 is not None: boyd_box0.print_timestepping_statistics() |
|---|
| 220 | |
|---|
| 221 | #for i in range(samples): |
|---|
| 222 | # if tri_ids[i] >= 0: |
|---|
| 223 | # if verbose: print 'P%d tri %d, value = %s' %(myid, i, stage.centroid_values[tri_ids[i]]) |
|---|
| 224 | |
|---|
| 225 | sys.stdout.flush() |
|---|
| 226 | |
|---|
| 227 | pass |
|---|
| 228 | |
|---|
| 229 | success = True |
|---|
| 230 | |
|---|
| 231 | ##----------------------------------------------------------------------- |
|---|
| 232 | ## Assign/Test Control data |
|---|
| 233 | ##----------------------------------------------------------------------- |
|---|
| 234 | |
|---|
| 235 | if not parallel: |
|---|
| 236 | stage = domain.get_quantity('stage') |
|---|
| 237 | |
|---|
| 238 | for i in range(samples): |
|---|
| 239 | assert(tri_ids[i] >= 0) |
|---|
| 240 | control_data.append(stage.centroid_values[tri_ids[i]]) |
|---|
| 241 | |
|---|
| 242 | if inlet0 is not None: |
|---|
| 243 | control_data.append(inlet0.inlet.get_average_stage()) |
|---|
| 244 | control_data.append(inlet0.inlet.get_average_xmom()) |
|---|
| 245 | control_data.append(inlet0.inlet.get_average_ymom()) |
|---|
| 246 | control_data.append(inlet0.inlet.get_total_water_volume()) |
|---|
| 247 | control_data.append(inlet0.inlet.get_average_depth()) |
|---|
| 248 | |
|---|
| 249 | if verbose: print 'P%d control_data = %s' %(myid, control_data) |
|---|
| 250 | else: |
|---|
| 251 | stage = domain.get_quantity('stage') |
|---|
| 252 | |
|---|
| 253 | for i in range(samples): |
|---|
| 254 | if tri_ids[i] >= 0: |
|---|
| 255 | local_success = num.allclose(control_data[i], stage.centroid_values[tri_ids[i]]) |
|---|
| 256 | success = success and local_success |
|---|
| 257 | if verbose: |
|---|
| 258 | print 'P%d tri %d, control = %s, actual = %s, Success = %s' %(myid, i, control_data[i], stage.centroid_values[tri_ids[i]], local_success) |
|---|
| 259 | |
|---|
| 260 | |
|---|
| 261 | if inlet0 is not None: |
|---|
| 262 | inlet_master_proc = inlet0.inlet.get_master_proc() |
|---|
| 263 | average_stage = inlet0.inlet.get_global_average_stage() |
|---|
| 264 | average_xmom = inlet0.inlet.get_global_average_xmom() |
|---|
| 265 | average_ymom = inlet0.inlet.get_global_average_ymom() |
|---|
| 266 | average_volume = inlet0.inlet.get_global_total_water_volume() |
|---|
| 267 | average_depth = inlet0.inlet.get_global_average_depth() |
|---|
| 268 | |
|---|
| 269 | if myid == inlet_master_proc: |
|---|
| 270 | if verbose: |
|---|
| 271 | print 'P%d average stage, control = %s, actual = %s' %(myid, control_data[samples], average_stage) |
|---|
| 272 | |
|---|
| 273 | print 'P%d average xmom, control = %s, actual = %s' %(myid, control_data[samples+1], average_xmom) |
|---|
| 274 | |
|---|
| 275 | print 'P%d average ymom, control = %s, actual = %s' %(myid, control_data[samples+2], average_ymom) |
|---|
| 276 | |
|---|
| 277 | print 'P%d average volume, control = %s, actual = %s' %(myid, control_data[samples+3], average_volume) |
|---|
| 278 | |
|---|
| 279 | print 'P%d average depth, control = %s, actual = %s' %(myid, control_data[samples+4], average_depth) |
|---|
| 280 | |
|---|
| 281 | |
|---|
| 282 | #assert(success) |
|---|
| 283 | |
|---|
| 284 | return control_data |
|---|
| 285 | |
|---|
| 286 | |
|---|
| 287 | if __name__=="__main__": |
|---|
| 288 | |
|---|
| 289 | test_points = [] |
|---|
| 290 | |
|---|
| 291 | if myid == 0: |
|---|
| 292 | |
|---|
| 293 | for i in range(samples): |
|---|
| 294 | x = random.randrange(0,1000)/1000.0 * length |
|---|
| 295 | y = random.randrange(0,1000)/1000.0 * width |
|---|
| 296 | point = [x, y] |
|---|
| 297 | test_points.append(point) |
|---|
| 298 | |
|---|
| 299 | for i in range(1,numprocs): |
|---|
| 300 | pypar.send(test_points, i) |
|---|
| 301 | else: |
|---|
| 302 | test_points = pypar.receive(0) |
|---|
| 303 | |
|---|
| 304 | print "Test Points::" |
|---|
| 305 | print test_points |
|---|
| 306 | |
|---|
| 307 | if myid == 0: |
|---|
| 308 | control_data = run_test(parallel=False, test_points = test_points, verbose = True) |
|---|
| 309 | |
|---|
| 310 | for proc in range(1,numprocs): |
|---|
| 311 | pypar.send(control_data, proc) |
|---|
| 312 | else: |
|---|
| 313 | control_data = pypar.receive(0) |
|---|
| 314 | |
|---|
| 315 | pypar.barrier() |
|---|
| 316 | run_test(parallel=True, control_data = control_data, test_points = test_points, verbose = True) |
|---|
| 317 | |
|---|
| 318 | |
|---|
| 319 | finalize() |
|---|
