source: trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py @ 7770

Last change on this file since 7770 was 7770, checked in by hudson, 14 years ago

Broke up data_manager into more modules - some failing tests.

File size: 27.6 KB
Line 
1'''
2    Operations to extract information from an SWW file.
3'''
4
5##
6# @brief Get values for quantities interpolated to polyline midpoints from SWW.
7# @param filename Path to file to read.
8# @param quantity_names Quantity names to get.
9# @param polyline Representation of desired cross-section.
10# @param verbose True if this function is to be verbose.
11# @return (segments, i_func) where segments is a list of Triangle_intersection
12#         instances and i_func is an instance of Interpolation_function.
13# @note For 'polyline' assume absolute UTM coordinates.
14def get_interpolated_quantities_at_polyline_midpoints(filename,
15                                                      quantity_names=None,
16                                                      polyline=None,
17                                                      verbose=False):
18    """Get values for quantities interpolated to polyline midpoints from SWW
19
20    Input:
21        filename - Name of sww file
22        quantity_names - Names of quantities to load
23        polyline: Representation of desired cross section - it may contain
24                  multiple sections allowing for complex shapes. Assume
25                  absolute UTM coordinates.
26                  Format [[x0, y0], [x1, y1], ...]
27
28    Output:
29        segments: list of instances of class Triangle_intersection
30        interpolation_function: Instance of class Interpolation_function
31
32
33      This function is used by get_flow_through_cross_section and
34      get_energy_through_cross_section
35    """
36
37    from anuga.fit_interpolate.interpolate import Interpolation_function
38
39    # Get mesh and quantities from sww file
40    X = get_mesh_and_quantities_from_file(filename,
41                                          quantities=quantity_names,
42                                          verbose=verbose)
43    mesh, quantities, time = X
44
45    # Find all intersections and associated triangles.
46    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
47
48    # Get midpoints
49    interpolation_points = segment_midpoints(segments)
50
51    # Interpolate
52    if verbose:
53        log.critical('Interpolating - total number of interpolation points = %d'
54                     % len(interpolation_points))
55
56    I = Interpolation_function(time,
57                               quantities,
58                               quantity_names=quantity_names,
59                               vertex_coordinates=mesh.nodes,
60                               triangles=mesh.triangles,
61                               interpolation_points=interpolation_points,
62                               verbose=verbose)
63
64    return segments, I
65
66
67##
68# @brief Obtain flow (m^3/s) perpendicular to specified cross section.
69# @param filename Path to file to read.
70# @param polyline Representation of desired cross-section.
71# @param verbose Trie if this function is to be verbose.
72# @return (time, Q) where time and Q are lists of time and flow respectively.
73def get_flow_through_cross_section(filename, polyline, verbose=False):
74    """Obtain flow (m^3/s) perpendicular to specified cross section.
75
76    Inputs:
77        filename: Name of sww file
78        polyline: Representation of desired cross section - it may contain
79                  multiple sections allowing for complex shapes. Assume
80                  absolute UTM coordinates.
81                  Format [[x0, y0], [x1, y1], ...]
82
83    Output:
84        time: All stored times in sww file
85        Q: Hydrograph of total flow across given segments for all stored times.
86
87    The normal flow is computed for each triangle intersected by the polyline
88    and added up.  Multiple segments at different angles are specified the
89    normal flows may partially cancel each other.
90
91    The typical usage of this function would be to get flow through a channel,
92    and the polyline would then be a cross section perpendicular to the flow.
93    """
94
95    quantity_names =['elevation',
96                     'stage',
97                     'xmomentum',
98                     'ymomentum']
99
100    # Get values for quantities at each midpoint of poly line from sww file
101    X = get_interpolated_quantities_at_polyline_midpoints(filename,
102                                                          quantity_names=\
103                                                              quantity_names,
104                                                          polyline=polyline,
105                                                          verbose=verbose)
106    segments, interpolation_function = X
107
108    # Get vectors for time and interpolation_points
109    time = interpolation_function.time
110    interpolation_points = interpolation_function.interpolation_points
111
112    if verbose: log.critical('Computing hydrograph')
113
114    # Compute hydrograph
115    Q = []
116    for t in time:
117        total_flow = 0
118        for i in range(len(interpolation_points)):
119            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
120            normal = segments[i].normal
121
122            # Inner product of momentum vector with segment normal [m^2/s]
123            normal_momentum = uh*normal[0] + vh*normal[1]
124
125            # Flow across this segment [m^3/s]
126            segment_flow = normal_momentum * segments[i].length
127
128            # Accumulate
129            total_flow += segment_flow
130
131        # Store flow at this timestep
132        Q.append(total_flow)
133
134
135    return time, Q
136
137
138##
139# @brief Get average energy across a cross-section.
140# @param filename Path to file of interest.
141# @param polyline Representation of desired cross-section.
142# @param kind Select energy to compute: 'specific' or 'total'.
143# @param verbose True if this function is to be verbose.
144# @return (time, E) where time and E are lists of timestep and energy.
145def get_energy_through_cross_section(filename,
146                                     polyline,
147                                     kind='total',
148                                     verbose=False):
149    """Obtain average energy head [m] across specified cross section.
150
151    Inputs:
152        polyline: Representation of desired cross section - it may contain
153                  multiple sections allowing for complex shapes. Assume
154                  absolute UTM coordinates.
155                  Format [[x0, y0], [x1, y1], ...]
156        kind:     Select which energy to compute.
157                  Options are 'specific' and 'total' (default)
158
159    Output:
160        E: Average energy [m] across given segments for all stored times.
161
162    The average velocity is computed for each triangle intersected by
163    the polyline and averaged weighted by segment lengths.
164
165    The typical usage of this function would be to get average energy of
166    flow in a channel, and the polyline would then be a cross section
167    perpendicular to the flow.
168
169    #FIXME (Ole) - need name for this energy reflecting that its dimension
170    is [m].
171    """
172
173    from anuga.config import g, epsilon, velocity_protection as h0
174
175    quantity_names =['elevation',
176                     'stage',
177                     'xmomentum',
178                     'ymomentum']
179
180    # Get values for quantities at each midpoint of poly line from sww file
181    X = get_interpolated_quantities_at_polyline_midpoints(filename,
182                                                          quantity_names=\
183                                                              quantity_names,
184                                                          polyline=polyline,
185                                                          verbose=verbose)
186    segments, interpolation_function = X
187
188    # Get vectors for time and interpolation_points
189    time = interpolation_function.time
190    interpolation_points = interpolation_function.interpolation_points
191
192    if verbose: log.critical('Computing %s energy' % kind)
193
194    # Compute total length of polyline for use with weighted averages
195    total_line_length = 0.0
196    for segment in segments:
197        total_line_length += segment.length
198
199    # Compute energy
200    E = []
201    for t in time:
202        average_energy = 0.0
203        for i, p in enumerate(interpolation_points):
204            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
205
206            # Depth
207            h = depth = stage-elevation
208
209            # Average velocity across this segment
210            if h > epsilon:
211                # Use protection against degenerate velocities
212                u = uh / (h + h0/h)
213                v = vh / (h + h0/h)
214            else:
215                u = v = 0.0
216
217            speed_squared = u*u + v*v
218            kinetic_energy = 0.5 * speed_squared / g
219
220            if kind == 'specific':
221                segment_energy = depth + kinetic_energy
222            elif kind == 'total':
223                segment_energy = stage + kinetic_energy
224            else:
225                msg = 'Energy kind must be either "specific" or "total". '
226                msg += 'I got %s' % kind
227
228            # Add to weighted average
229            weigth = segments[i].length / total_line_length
230            average_energy += segment_energy * weigth
231
232        # Store energy at this timestep
233        E.append(average_energy)
234
235    return time, E
236
237
238##
239# @brief Return highest elevation where depth > 0.
240# @param filename Path to SWW file of interest.
241# @param polygon If specified resrict to points inside this polygon.
242# @param time_interval If specified resrict to within the time specified.
243# @param verbose True if this function is  to be verbose.
244def get_maximum_inundation_elevation(filename,
245                                     polygon=None,
246                                     time_interval=None,
247                                     verbose=False):
248    """Return highest elevation where depth > 0
249
250    Usage:
251    max_runup = get_maximum_inundation_elevation(filename,
252                                                 polygon=None,
253                                                 time_interval=None,
254                                                 verbose=False)
255
256    filename is a NetCDF sww file containing ANUGA model output.
257    Optional arguments polygon and time_interval restricts the maximum
258    runup calculation
259    to a points that lie within the specified polygon and time interval.
260
261    If no inundation is found within polygon and time_interval the return value
262    is None signifying "No Runup" or "Everything is dry".
263
264    See general function get_maximum_inundation_data for details.
265    """
266
267    runup, _ = get_maximum_inundation_data(filename,
268                                           polygon=polygon,
269                                           time_interval=time_interval,
270                                           verbose=verbose)
271    return runup
272
273
274##
275# @brief Return location of highest elevation where h > 0
276# @param filename Path to SWW file to read.
277# @param polygon If specified resrict to points inside this polygon.
278# @param time_interval If specified resrict to within the time specified.
279# @param verbose True if this function is  to be verbose.
280def get_maximum_inundation_location(filename,
281                                    polygon=None,
282                                    time_interval=None,
283                                    verbose=False):
284    """Return location of highest elevation where h > 0
285
286    Usage:
287    max_runup_location = get_maximum_inundation_location(filename,
288                                                         polygon=None,
289                                                         time_interval=None,
290                                                         verbose=False)
291
292    filename is a NetCDF sww file containing ANUGA model output.
293    Optional arguments polygon and time_interval restricts the maximum
294    runup calculation
295    to a points that lie within the specified polygon and time interval.
296
297    If no inundation is found within polygon and time_interval the return value
298    is None signifying "No Runup" or "Everything is dry".
299
300    See general function get_maximum_inundation_data for details.
301    """
302
303    _, max_loc = get_maximum_inundation_data(filename,
304                                             polygon=polygon,
305                                             time_interval=time_interval,
306                                             verbose=verbose)
307    return max_loc
308
309
310##
311# @brief Compute maximum run up height from SWW file.
312# @param filename Path to SWW file to read.
313# @param polygon If specified resrict to points inside this polygon.
314# @param time_interval If specified resrict to within the time specified.
315# @param use_centroid_values
316# @param verbose True if this function is to be verbose.
317# @return (maximal_runup, maximal_runup_location)
318def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
319                                use_centroid_values=False,
320                                verbose=False):
321    """Compute maximum run up height from sww file.
322
323    Usage:
324    runup, location = get_maximum_inundation_data(filename,
325                                                  polygon=None,
326                                                  time_interval=None,
327                                                  verbose=False)
328
329    Algorithm is as in get_maximum_inundation_elevation from
330    shallow_water_domain except that this function works with the sww file and
331    computes the maximal runup height over multiple timesteps.
332
333    Optional arguments polygon and time_interval restricts the maximum runup
334    calculation to a points that lie within the specified polygon and time
335    interval.
336
337    Polygon is assumed to be in (absolute) UTM coordinates in the same zone
338    as domain.
339
340    If no inundation is found within polygon and time_interval the return value
341    is None signifying "No Runup" or "Everything is dry".
342    """
343
344    # We are using nodal values here as that is what is stored in sww files.
345
346    # Water depth below which it is considered to be 0 in the model
347    # FIXME (Ole): Allow this to be specified as a keyword argument as well
348
349    from anuga.geometry.polygon import inside_polygon
350    from anuga.config import minimum_allowed_height
351    from Scientific.IO.NetCDF import NetCDFFile
352
353    dir, base = os.path.split(filename)
354
355    iterate_over = get_all_swwfiles(dir, base)
356
357    # Read sww file
358    if verbose: log.critical('Reading from %s' % filename)
359    # FIXME: Use general swwstats (when done)
360
361    maximal_runup = None
362    maximal_runup_location = None
363
364    for file, swwfile in enumerate (iterate_over):
365        # Read sww file
366        filename = join(dir, swwfile+'.sww')
367
368        if verbose: log.critical('Reading from %s' % filename)
369        # FIXME: Use general swwstats (when done)
370
371        fid = NetCDFFile(filename)
372
373        # Get geo_reference
374        # sww files don't have to have a geo_ref
375        try:
376            geo_reference = Geo_reference(NetCDFObject=fid)
377        except AttributeError, e:
378            geo_reference = Geo_reference() # Default georef object
379
380        xllcorner = geo_reference.get_xllcorner()
381        yllcorner = geo_reference.get_yllcorner()
382        zone = geo_reference.get_zone()
383
384        # Get extent
385        volumes = fid.variables['volumes'][:]
386        x = fid.variables['x'][:] + xllcorner
387        y = fid.variables['y'][:] + yllcorner
388
389        # Get the relevant quantities (Convert from single precison)
390        elevation = num.array(fid.variables['elevation'][:], num.float)
391        stage = num.array(fid.variables['stage'][:], num.float)
392
393        # Here's where one could convert nodal information to centroid
394        # information but is probably something we need to write in C.
395        # Here's a Python thought which is NOT finished!!!
396        if use_centroid_values is True:
397            x = get_centroid_values(x, volumes)
398            y = get_centroid_values(y, volumes)
399            elevation = get_centroid_values(elevation, volumes)
400
401        # Spatial restriction
402        if polygon is not None:
403            msg = 'polygon must be a sequence of points.'
404            assert len(polygon[0]) == 2, msg
405            # FIXME (Ole): Make a generic polygon input check in polygon.py
406            # and call it here
407            points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
408                                                            y[:,num.newaxis]),
409                                                            axis=1))
410            point_indices = inside_polygon(points, polygon)
411
412            # Restrict quantities to polygon
413            elevation = num.take(elevation, point_indices, axis=0)
414            stage = num.take(stage, point_indices, axis=1)
415
416            # Get info for location of maximal runup
417            points_in_polygon = num.take(points, point_indices, axis=0)
418
419            x = points_in_polygon[:,0]
420            y = points_in_polygon[:,1]
421        else:
422            # Take all points
423            point_indices = num.arange(len(x))
424
425        # Temporal restriction
426        time = fid.variables['time'][:]
427        all_timeindices = num.arange(len(time))
428        if time_interval is not None:
429            msg = 'time_interval must be a sequence of length 2.'
430            assert len(time_interval) == 2, msg
431            msg = 'time_interval %s must not be decreasing.' % time_interval
432            assert time_interval[1] >= time_interval[0], msg
433            msg = 'Specified time interval [%.8f:%.8f] ' % tuple(time_interval)
434            msg += 'must does not match model time interval: [%.8f, %.8f]\n' \
435                   % (time[0], time[-1])
436            if time_interval[1] < time[0]: raise ValueError(msg)
437            if time_interval[0] > time[-1]: raise ValueError(msg)
438
439            # Take time indices corresponding to interval (& is bitwise AND)
440            timesteps = num.compress((time_interval[0] <= time) \
441                                     & (time <= time_interval[1]),
442                                     all_timeindices)
443
444            msg = 'time_interval %s did not include any model timesteps.' \
445                  % time_interval
446            assert not num.alltrue(timesteps == 0), msg
447        else:
448            # Take them all
449            timesteps = all_timeindices
450
451        fid.close()
452
453        # Compute maximal runup for each timestep
454        #maximal_runup = None
455        #maximal_runup_location = None
456        #maximal_runups = [None]
457        #maximal_runup_locations = [None]
458
459        for i in timesteps:
460            if use_centroid_values is True:
461                stage_i = get_centroid_values(stage[i,:], volumes)
462            else:
463                stage_i = stage[i,:]
464
465            depth = stage_i - elevation
466
467            # Get wet nodes i.e. nodes with depth>0 within given region
468            # and timesteps
469            wet_nodes = num.compress(depth > minimum_allowed_height,
470                                     num.arange(len(depth)))
471
472            if num.alltrue(wet_nodes == 0):
473                runup = None
474            else:
475                # Find maximum elevation among wet nodes
476                wet_elevation = num.take(elevation, wet_nodes, axis=0)
477                runup_index = num.argmax(wet_elevation)
478                runup = max(wet_elevation)
479                assert wet_elevation[runup_index] == runup       # Must be True
480
481            if runup > maximal_runup:
482                maximal_runup = runup      # works even if maximal_runup is None
483
484                # Record location
485                wet_x = num.take(x, wet_nodes, axis=0)
486                wet_y = num.take(y, wet_nodes, axis=0)
487                maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]
488
489    return maximal_runup, maximal_runup_location
490
491
492
493    def test_get_maximum_inundation_from_sww(self):
494        """test_get_maximum_inundation_from_sww(self)
495
496        Test of get_maximum_inundation_elevation()
497        and get_maximum_inundation_location() from data_manager.py
498
499        This is based on test_get_maximum_inundation_3(self) but works with the
500        stored results instead of with the internal data structure.
501
502        This test uses the underlying get_maximum_inundation_data for tests
503        """
504
505        from anuga.abstract_2d_finite_volumes.mesh_factory \
506                import rectangular_cross
507        from data_manager import get_maximum_inundation_elevation
508        from data_manager import get_maximum_inundation_location
509        from data_manager import get_maximum_inundation_data
510
511        initial_runup_height = -0.4
512        final_runup_height = -0.3
513
514        #--------------------------------------------------------------
515        # Setup computational domain
516        #--------------------------------------------------------------
517        N = 10
518        points, vertices, boundary = rectangular_cross(N, N)
519        domain = Domain(points, vertices, boundary)
520        domain.set_name('runup_test')
521        domain.set_maximum_allowed_speed(1.0)
522
523        # FIXME: This works better with old limiters so far
524        domain.tight_slope_limiters = 0
525
526        #--------------------------------------------------------------
527        # Setup initial conditions
528        #--------------------------------------------------------------
529        def topography(x, y):
530            return -x/2                             # linear bed slope
531
532        # Use function for elevation
533        domain.set_quantity('elevation', topography)
534        domain.set_quantity('friction', 0.)                # Zero friction
535        # Constant negative initial stage
536        domain.set_quantity('stage', initial_runup_height)
537
538        #--------------------------------------------------------------
539        # Setup boundary conditions
540        #--------------------------------------------------------------
541        Br = Reflective_boundary(domain)                       # Reflective wall
542        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
543
544        # All reflective to begin with (still water)
545        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
546
547        #--------------------------------------------------------------
548        # Test initial inundation height
549        #--------------------------------------------------------------
550        indices = domain.get_wet_elements()
551        z = domain.get_quantity('elevation').\
552                get_values(location='centroids', indices=indices)
553        assert num.alltrue(z < initial_runup_height)
554
555        q_ref = domain.get_maximum_inundation_elevation()
556        # First order accuracy
557        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
558
559        #--------------------------------------------------------------
560        # Let triangles adjust
561        #--------------------------------------------------------------
562        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
563            pass
564
565        #--------------------------------------------------------------
566        # Test inundation height again
567        #--------------------------------------------------------------
568        q_ref = domain.get_maximum_inundation_elevation()
569        q = get_maximum_inundation_elevation('runup_test.sww')
570        msg = 'We got %f, should have been %f' % (q, q_ref)
571        assert num.allclose(q, q_ref, rtol=1.0/N), msg
572
573        q = get_maximum_inundation_elevation('runup_test.sww')
574        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
575        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
576
577        # Test error condition if time interval is out
578        try:
579            q = get_maximum_inundation_elevation('runup_test.sww',
580                                                 time_interval=[2.0, 3.0])
581        except ValueError:
582            pass
583        else:
584            msg = 'should have caught wrong time interval'
585            raise Exception, msg
586
587        # Check correct time interval
588        q, loc = get_maximum_inundation_data('runup_test.sww',
589                                             time_interval=[0.0, 3.0])
590        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
591        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
592        assert num.allclose(-loc[0]/2, q)    # From topography formula
593
594        #--------------------------------------------------------------
595        # Update boundary to allow inflow
596        #--------------------------------------------------------------
597        domain.set_boundary({'right': Bd})
598
599        #--------------------------------------------------------------
600        # Evolve system through time
601        #--------------------------------------------------------------
602        q_max = None
603        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
604                               skip_initial_step = True):
605            q = domain.get_maximum_inundation_elevation()
606            if q > q_max:
607                q_max = q
608
609        #--------------------------------------------------------------
610        # Test inundation height again
611        #--------------------------------------------------------------
612        indices = domain.get_wet_elements()
613        z = domain.get_quantity('elevation').\
614                get_values(location='centroids', indices=indices)
615
616        assert num.alltrue(z < final_runup_height)
617
618        q = domain.get_maximum_inundation_elevation()
619        # First order accuracy
620        assert num.allclose(q, final_runup_height, rtol=1.0/N)
621
622        q, loc = get_maximum_inundation_data('runup_test.sww',
623                                             time_interval=[3.0, 3.0])
624        msg = 'We got %f, should have been %f' % (q, final_runup_height)
625        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
626        assert num.allclose(-loc[0]/2, q)    # From topography formula
627
628        q = get_maximum_inundation_elevation('runup_test.sww')
629        loc = get_maximum_inundation_location('runup_test.sww')
630        msg = 'We got %f, should have been %f' % (q, q_max)
631        assert num.allclose(q, q_max, rtol=1.0/N), msg
632        assert num.allclose(-loc[0]/2, q)    # From topography formula
633
634        q = get_maximum_inundation_elevation('runup_test.sww',
635                                             time_interval=[0, 3])
636        msg = 'We got %f, should have been %f' % (q, q_max)
637        assert num.allclose(q, q_max, rtol=1.0/N), msg
638
639        # Check polygon mode
640        # Runup region
641        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
642        q = get_maximum_inundation_elevation('runup_test.sww',
643                                             polygon = polygon,
644                                             time_interval=[0, 3])
645        msg = 'We got %f, should have been %f' % (q, q_max)
646        assert num.allclose(q, q_max, rtol=1.0/N), msg
647
648        # Offshore region
649        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
650        q, loc = get_maximum_inundation_data('runup_test.sww',
651                                             polygon = polygon,
652                                             time_interval=[0, 3])
653        msg = 'We got %f, should have been %f' % (q, -0.475)
654        assert num.allclose(q, -0.475, rtol=1.0/N), msg
655        assert is_inside_polygon(loc, polygon)
656        assert num.allclose(-loc[0]/2, q)    # From topography formula
657
658        # Dry region
659        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
660        q, loc = get_maximum_inundation_data('runup_test.sww',
661                                             polygon = polygon,
662                                             time_interval=[0, 3])
663        msg = 'We got %s, should have been None' % (q)
664        assert q is None, msg
665        msg = 'We got %s, should have been None' % (loc)
666        assert loc is None, msg
667
668        # Check what happens if no time point is within interval
669        try:
670            q = get_maximum_inundation_elevation('runup_test.sww',
671                                                 time_interval=[2.75, 2.75])
672        except AssertionError:
673            pass
674        else:
675            msg = 'Time interval should have raised an exception'
676            raise Exception, msg
677
678        # Cleanup
679        try:
680            os.remove(domain.get_name() + '.sww')
681        except:
682            pass
683            #FIXME(Ole): Windows won't allow removal of this
684
Note: See TracBrowser for help on using the repository browser.