source: trunk/anuga_core/source/anuga/shallow_water/forcing.py @ 8780

Last change on this file since 8780 was 8780, checked in by steve, 11 years ago

Some changes to allow netcdf4 use
File size: 47.4 KB
Line 
1"""
2Environmental forcing functions, such as wind and rainfall.
3
4Constraints: See GPL license in the user guide
5Version: 1.0 ($Revision: 7731 $)
6ModifiedBy:
7    $Author: hudson $
8    $Date: 2010-05-18 14:54:05 +1000 (Tue, 18 May 2010) $
9"""
10
11
12from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
13from anuga.utilities.numerical_tools import ensure_numeric
14from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
15                                              Modeltime_too_late
16from anuga.geometry.polygon import is_inside_polygon, inside_polygon, \
17                                    polygon_area
18from anuga.geospatial_data.geospatial_data import ensure_geospatial
19
20from warnings import warn
21import numpy as num
22from anuga.file.netcdf import NetCDFFile
23from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
24from copy import copy
25
26def check_forcefield(f):
27    """Check that force object is as expected.
28   
29    Check that f is either:
30    1: a callable object f(t,x,y), where x and y are vectors
31       and that it returns an array or a list of same length
32       as x and y
33    2: a scalar
34    """
35
36    if callable(f):
37        N = 3
38        x = num.ones(3, num.float)
39        y = num.ones(3, num.float)
40        try:
41            q = f(1.0, x=x, y=y)
42        except Exception, e:
43            msg = 'Function %s could not be executed:\n%s' %(f, e)
44            # FIXME: Reconsider this semantics
45            raise Exception(msg)
46
47        try:
48            q = num.array(q, num.float)
49        except:
50            msg = ('Return value from vector function %s could not '
51                   'be converted into a numeric array of floats.\nSpecified '
52                   'function should return either list or array.' % f)
53            raise Exception(msg)
54
55        # Is this really what we want?
56        # info is "(func name, filename, defining line)"
57        func_info = (f.func_name, f.func_code.co_filename,
58                     f.func_code.co_firstlineno)
59        func_msg = 'Function %s (defined in %s, line %d)' % func_info
60        try:
61            result_len = len(q)
62        except:
63            msg = '%s must return vector' % func_msg
64            raise Exception(msg)
65        msg = '%s must return vector of length %d' % (func_msg, N)
66        assert result_len == N, msg
67    else:
68        try:
69            f = float(f)
70        except:
71            msg = ('Force field %s must be a scalar value coercible to float.'
72                   % str(f))
73            raise Exception(msg)
74
75    return f
76
77
78
79class Wind_stress:
80    """Apply wind stress to water momentum in terms of
81    wind speed [m/s] and wind direction [degrees]
82    """
83    def __init__(self, *args, **kwargs):
84        """Initialise windfield from wind speed s [m/s]
85        and wind direction phi [degrees]
86
87        Inputs v and phi can be either scalars or Python functions, e.g.
88
89        W = Wind_stress(10, 178)
90
91        #FIXME - 'normal' degrees are assumed for now, i.e. the
92        vector (1,0) has zero degrees.
93        We may need to convert from 'compass' degrees later on and also
94        map from True north to grid north.
95
96        Arguments can also be Python functions of t,x,y as in
97
98        def speed(t,x,y):
99            ...
100            return s
101
102        def angle(t,x,y):
103            ...
104            return phi
105
106        where x and y are vectors.
107
108        and then pass the functions in
109
110        W = Wind_stress(speed, angle)
111
112        The instantiated object W can be appended to the list of
113        forcing_terms as in
114
115        Alternatively, one vector valued function for (speed, angle)
116        can be applied, providing both quantities simultaneously.
117        As in
118        W = Wind_stress(F), where returns (speed, angle) for each t.
119
120        domain.forcing_terms.append(W)
121        """
122
123        from anuga.config import rho_a, rho_w, eta_w
124
125        self.use_coordinates=True
126        if len(args) == 2:
127            s = args[0]
128            phi = args[1]
129        elif len(args) == 1:
130            # Assume vector function returning (s, phi)(t,x,y)
131            vector_function = args[0]
132            if ( len(kwargs)==1 ):
133                self.use_coordinates=kwargs['use_coordinates']
134            else:
135                self.use_coordinates=True
136            if ( self.use_coordinates ):
137                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
138                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
139            else:
140                s = lambda t,i: vector_function(t,point_id=i)[0]
141                phi = lambda t,i: vector_function(t,point_id=i)[1]
142        else:
143           # Assume info is in 2 keyword arguments
144           if len(kwargs) == 2:
145               s = kwargs['s']
146               phi = kwargs['phi']
147           else:
148               raise Exception('Assumes two keyword arguments: s=..., phi=....')
149
150        if ( self.use_coordinates ):
151            self.speed = check_forcefield(s)
152            self.phi = check_forcefield(phi)
153        else:
154            self.speed = s
155            self.phi = phi
156
157        self.const = eta_w*rho_a/rho_w
158
159    def __call__(self, domain):
160        """Evaluate windfield based on values found in domain"""
161
162        xmom_update = domain.quantities['xmomentum'].explicit_update
163        ymom_update = domain.quantities['ymomentum'].explicit_update
164
165        N = len(domain)    # number_of_triangles
166        t = domain.time
167
168        if callable(self.speed):
169            xc = domain.get_centroid_coordinates()
170            if ( self.use_coordinates ):
171                s_vec = self.speed(t, xc[:,0], xc[:,1])
172            else:
173                s_vec=num.empty(N,float)
174                for i in range(N):
175                    s_vec[i]=self.speed(t,i)
176        else:
177            # Assume s is a scalar
178            try:
179                s_vec = self.speed * num.ones(N, num.float)
180            except:
181                msg = 'Speed must be either callable or a scalar: %s' %self.s
182                raise Exception(msg)
183
184        if callable(self.phi):
185            xc = domain.get_centroid_coordinates()
186            if ( self.use_coordinates ):
187                phi_vec = self.phi(t, xc[:,0], xc[:,1])
188            else:
189                phi_vec=num.empty(len(xc),float)
190                for i in range(len(xc)):
191                    phi_vec[i]=self.phi(t,i)
192        else:
193            # Assume phi is a scalar
194
195            try:
196                phi_vec = self.phi * num.ones(N, num.float)
197            except:
198                msg = 'Angle must be either callable or a scalar: %s' %self.phi
199                raise Exception(msg)
200
201        assign_windfield_values(xmom_update, ymom_update,
202                                s_vec, phi_vec, self.const)
203
204
205def assign_windfield_values(xmom_update, ymom_update,
206                            s_vec, phi_vec, const):
207    """Python version of assigning wind field to update vectors.
208    """
209
210    from math import pi, cos, sin, sqrt
211
212    N = len(s_vec)
213    for k in range(N):
214        s = s_vec[k]
215        phi = phi_vec[k]
216
217        # Convert to radians
218        phi = phi*pi/180
219
220        # Compute velocity vector (u, v)
221        u = s*cos(phi)
222        v = s*sin(phi)
223
224        # Compute wind stress
225        S = const * sqrt(u**2 + v**2)
226        xmom_update[k] += S*u
227        ymom_update[k] += S*v
228
229
230class General_forcing:
231    """General explicit forcing term for update of quantity
232
233    This is used by Inflow and Rainfall for instance
234
235
236    General_forcing(quantity_name, rate, center, radius, polygon)
237
238    domain:     ANUGA computational domain
239    quantity_name: Name of quantity to update.
240                   It must be a known conserved quantity.
241
242    rate [?/s]: Total rate of change over the specified area.
243                This parameter can be either a constant or a
244                function of time. Positive values indicate increases,
245                negative values indicate decreases.
246                Rate can be None at initialisation but must be specified
247                before forcing term is applied (i.e. simulation has started).
248
249    center [m]: Coordinates at center of flow point
250    radius [m]: Size of circular area
251    polygon:    Arbitrary polygon
252    default_rate: Rate to be used if rate fails (e.g. if model time exceeds its data)
253                  Admissible types: None, constant number or function of t
254
255
256    Either center, radius or polygon can be specified but not both.
257    If neither are specified the entire domain gets updated.
258    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
259
260    Inflow or Rainfall for examples of use
261    """
262
263
264    # FIXME (AnyOne) : Add various methods to allow spatial variations
265
266    def __init__(self,
267                 domain,
268                 quantity_name,
269                 rate=0.0,
270                 center=None,
271                 radius=None,
272                 polygon=None,
273                 default_rate=None,
274                 verbose=False):
275
276        from math import pi, cos, sin
277
278        if domain.numproc > 1:
279            msg = 'Not implemented to run in parallel'
280            assert self.parallel_safe(), msg
281
282        if center is None:
283            msg = 'I got radius but no center.'
284            assert radius is None, msg
285
286        if radius is None:
287            msg += 'I got center but no radius.'
288            assert center is None, msg
289
290        self.domain = domain
291        self.quantity_name = quantity_name
292        self.rate = rate
293        self.center = ensure_numeric(center)
294        self.radius = radius
295        self.polygon = polygon
296        self.verbose = verbose
297        self.value = 0.0    # Can be used to remember value at
298                            # previous timestep in order to obtain rate
299
300        # Get boundary (in absolute coordinates)
301        bounding_polygon = domain.get_boundary_polygon()
302
303        # Update area if applicable
304        if center is not None and radius is not None:
305            assert len(center) == 2
306            msg = 'Polygon cannot be specified when center and radius are'
307            assert polygon is None, msg
308
309            # Check that circle center lies within the mesh.
310            msg = 'Center %s specified for forcing term did not' % str(center)
311            msg += 'fall within the domain boundary.'
312            assert is_inside_polygon(center, bounding_polygon), msg
313
314            # Check that circle periphery lies within the mesh.
315            N = 100
316            periphery_points = []
317            for i in range(N):
318                theta = 2*pi*i/100
319
320                x = center[0] + radius*cos(theta)
321                y = center[1] + radius*sin(theta)
322
323                periphery_points.append([x,y])
324
325            for point in periphery_points:
326                msg = 'Point %s on periphery for forcing term' % str(point)
327                msg += ' did not fall within the domain boundary.'
328                assert is_inside_polygon(point, bounding_polygon), msg
329
330        if polygon is not None:
331            # Check that polygon lies within the mesh.
332            for point in self.polygon:
333                msg = 'Point %s in polygon for forcing term' % str(point)
334                msg += ' did not fall within the domain boundary.'
335                assert is_inside_polygon(point, bounding_polygon), msg
336
337        # Pointer to update vector
338        self.update = domain.quantities[self.quantity_name].explicit_update
339
340        # Determine indices in flow area
341        N = len(domain)
342        points = domain.get_centroid_coordinates(absolute=True)
343
344        # Calculate indices in exchange area for this forcing term
345        self.exchange_indices = None
346        if self.center is not None and self.radius is not None:
347            # Inlet is circular
348            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
349
350            self.exchange_indices = []
351            for k in range(N):
352                x, y = points[k,:]    # Centroid
353
354                c = self.center
355                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
356                    self.exchange_indices.append(k)
357
358        if self.polygon is not None:
359            # Inlet is polygon
360            self.exchange_indices = inside_polygon(points, self.polygon)
361
362        if self.exchange_indices is None:
363            self.exchange_area = polygon_area(bounding_polygon)
364        else:   
365            if len(self.exchange_indices) == 0:
366                msg = 'No triangles have been identified in '
367                msg += 'specified region: %s' % inlet_region
368                raise Exception(msg)
369
370            # Compute exchange area as the sum of areas of triangles identified
371            # by circle or polygon
372            self.exchange_area = 0.0
373            for i in self.exchange_indices:
374                self.exchange_area += domain.areas[i]
375           
376
377        msg = 'Exchange area in forcing term'
378        msg += ' has area = %f' %self.exchange_area
379        assert self.exchange_area > 0.0           
380           
381               
382
383           
384        # Check and store default_rate
385        msg = ('Keyword argument default_rate must be either None '
386               'or a function of time.\nI got %s.' % str(default_rate))
387        assert (default_rate is None or
388                isinstance(default_rate, (int, float)) or
389                callable(default_rate)), msg
390
391        if default_rate is not None:
392            # If it is a constant, make it a function
393            if not callable(default_rate):
394                tmp = default_rate
395                default_rate = lambda t: tmp
396
397            # Check that default_rate is a function of one argument
398            try:
399                default_rate(0.0)
400            except:
401                raise Exception(msg)
402
403        self.default_rate = default_rate
404        self.default_rate_invoked = False    # Flag
405
406    def __call__(self, domain):
407        """Apply inflow function at time specified in domain, update stage"""
408
409        # Call virtual method allowing local modifications
410        t = domain.get_time()
411        try:
412            rate = self.update_rate(t)
413        except Modeltime_too_early, e:
414            raise Modeltime_too_early(e)
415        except Modeltime_too_late, e:
416            if self.default_rate is None:
417                msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
418                msg += 'You can specify keyword argument default_rate in the '
419                msg += 'forcing function to tell it what to do in the absence of time data.'
420                raise Modeltime_too_late(msg)
421            else:
422                # Pass control to default rate function
423                rate = self.default_rate(t)
424
425                if self.default_rate_invoked is False:
426                    # Issue warning the first time
427                    msg = ('%s\n'
428                           'Instead I will use the default rate: %s\n'
429                           'Note: Further warnings will be supressed'
430                           % (str(e), str(self.default_rate)))
431                    warn(msg)
432
433                    # FIXME (Ole): Replace this crude flag with
434                    # Python's ability to print warnings only once.
435                    # See http://docs.python.org/lib/warning-filter.html
436                    self.default_rate_invoked = True
437
438        if rate is None:
439            msg = ('Attribute rate must be specified in General_forcing '
440                   'or its descendants before attempting to call it')
441            raise Exception(msg)
442
443        # Now rate is a number
444        if self.verbose is True:
445            log.critical('Rate of %s at time = %.2f = %f'
446                         % (self.quantity_name, domain.get_time(), rate))
447
448        if self.exchange_indices is None:
449            self.update[:] += rate
450        else:
451            # Brute force assignment of restricted rate
452            for k in self.exchange_indices:
453                self.update[k] += rate
454
455    def update_rate(self, t):
456        """Virtual method allowing local modifications by writing an
457        overriding version in descendant
458        """
459
460        if callable(self.rate):
461            rate = self.rate(t)
462        else:
463            rate = self.rate
464
465        return rate
466
467    def get_quantity_values(self, quantity_name=None):
468        """Return values for specified quantity restricted to opening
469
470        Optionally a quantity name can be specified if values from another
471        quantity is sought
472        """
473
474        if quantity_name is None:
475            quantity_name = self.quantity_name
476
477        q = self.domain.quantities[quantity_name]
478        return q.get_values(location='centroids',
479                            indices=self.exchange_indices)
480
481    def set_quantity_values(self, val, quantity_name=None):
482        """Set values for specified quantity restricted to opening
483
484        Optionally a quantity name can be specified if values from another
485        quantity is sought
486        """
487
488        if quantity_name is None:
489            quantity_name = self.quantity_name
490
491        q = self.domain.quantities[self.quantity_name]
492        q.set_values(val,
493                     location='centroids',
494                     indices=self.exchange_indices)
495
496
497    def parallel_safe(self):
498        """
499        These forcing terms only work on individual processors, the polygon
500        had better not stride over multiple sub meshes
501        """
502        return True
503
504class Rainfall(General_forcing):
505    """Class Rainfall - general 'rain over entire domain' forcing term.
506
507    Used for implementing Rainfall over the entire domain.
508
509        Current Limited to only One Gauge..
510
511        Need to add Spatial Varying Capability
512        (This module came from copying and amending the Inflow Code)
513
514    Rainfall(rain)
515
516    domain
517    rain [mm/s]:  Total rain rate over the specified domain.
518                  NOTE: Raingauge Data needs to reflect the time step.
519                  IE: if Gauge is mm read at a time step, then the input
520                  here is as mm/(timeStep) so 10mm in 5minutes becomes
521                  10/(5x60) = 0.0333mm/s.
522
523                  This parameter can be either a constant or a
524                  function of time. Positive values indicate inflow,
525                  negative values indicate outflow.
526                  (and be used for Infiltration - Write Seperate Module)
527                  The specified flow will be divided by the area of
528                  the inflow region and then applied to update the
529                  stage quantity.
530
531    polygon: Specifies a polygon to restrict the rainfall.
532
533    Examples
534    How to put them in a run File...
535
536    #------------------------------------------------------------------------
537    # Setup specialised forcing terms
538    #------------------------------------------------------------------------
539    # This is the new element implemented by Ole and Rudy to allow direct
540    # input of Rainfall in mm/s
541
542    catchmentrainfall = Rainfall(rate=file_function('Q100_2hr_Rain.tms'))
543                        # Note need path to File in String.
544                        # Else assumed in same directory
545
546    domain.forcing_terms.append(catchmentrainfall)
547    """
548
549    def __init__(self,
550                 domain,
551                 rate=0.0,
552                 center=None,
553                 radius=None,
554                 polygon=None,
555                 default_rate=None,
556                 verbose=False):
557
558        # Converting mm/s to m/s to apply in ANUGA)
559        if callable(rate):
560            rain = lambda t: rate(t)/1000.0
561        else:
562            rain = rate/1000.0
563
564        if default_rate is not None:
565            if callable(default_rate):
566                default_rain = lambda t: default_rate(t)/1000.0
567            else:
568                default_rain = default_rate/1000.0
569        else:
570            default_rain = None
571
572           
573           
574        General_forcing.__init__(self,
575                                 domain,
576                                 'stage',
577                                 rate=rain,
578                                 center=center,
579                                 radius=radius,
580                                 polygon=polygon,
581                                 default_rate=default_rain,
582                                 verbose=verbose)
583
584
585class Inflow(General_forcing):
586    """Class Inflow - general 'rain and drain' forcing term.
587
588    Useful for implementing flows in and out of the domain.
589
590    Inflow(flow, center, radius, polygon)
591
592    domain
593    rate [m^3/s]: Total flow rate over the specified area.
594                  This parameter can be either a constant or a
595                  function of time. Positive values indicate inflow,
596                  negative values indicate outflow.
597                  The specified flow will be divided by the area of
598                  the inflow region and then applied to update stage.
599    center [m]: Coordinates at center of flow point
600    radius [m]: Size of circular area
601    polygon:    Arbitrary polygon.
602
603    Either center, radius or polygon must be specified
604
605    Examples
606
607    # Constant drain at 0.003 m^3/s.
608    # The outflow area is 0.07**2*pi=0.0154 m^2
609    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
610    #
611    Inflow((0.7, 0.4), 0.07, -0.003)
612
613
614    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
615    # The inflow area is 0.03**2*pi = 0.00283 m^2
616    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
617    # over the specified area
618    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
619
620
621    #------------------------------------------------------------------------
622    # Setup specialised forcing terms
623    #------------------------------------------------------------------------
624    # This is the new element implemented by Ole to allow direct input
625    # of Inflow in m^3/s
626
627    hydrograph = Inflow(center=(320, 300), radius=10,
628                        rate=file_function('Q/QPMF_Rot_Sub13.tms'))
629
630    domain.forcing_terms.append(hydrograph)
631    """
632
633    def __init__(self,
634                 domain,
635                 rate=0.0,
636                 center=None,
637                 radius=None,
638                 polygon=None,
639                 default_rate=None,
640                 verbose=False):
641        """Create an instance of the class
642
643        domain        Domain of interest
644        rate          Total rain rate over the specified domain (mm/s)
645        center
646        radius
647        polygon       Polygon to restrict rainfall
648        default_rate
649        verbose       True if this instance is to be verbose
650        """
651
652        # Create object first to make area is available
653        General_forcing.__init__(self,
654                                 domain,
655                                 'stage',
656                                 rate=rate,
657                                 center=center,
658                                 radius=radius,
659                                 polygon=polygon,
660                                 default_rate=default_rate,
661                                 verbose=verbose)
662
663    def update_rate(self, t):
664        """Virtual method allowing local modifications by writing an
665        overriding version in descendant
666
667        t  New rate object
668
669        This one converts m^3/s to m/s which can be added directly
670        to 'stage' in ANUGA
671        """
672
673        if callable(self.rate):
674            _rate = self.rate(t)/self.exchange_area
675        else:
676            _rate = self.rate/self.exchange_area
677
678        return _rate
679
680
681class Cross_section:
682    """Class Cross_section - a class to setup a cross section from
683    which you can then calculate flow and energy through cross section
684
685    Cross_section(domain, polyline)
686
687    domain:
688    polyline: Representation of desired cross section - it may contain
689              multiple sections allowing for complex shapes. Assume
690              absolute UTM coordinates.
691              Format [[x0, y0], [x1, y1], ...]
692    verbose:
693    """
694
695    def __init__(self,
696                 domain,
697                 polyline=None,
698                 verbose=False):
699        """Create an instance of Cross_section.
700
701        domain    domain of interest
702        polyline  polyline defining cross section
703        verbose   True if this instance is to be verbose
704        """
705       
706        self.domain = domain
707        self.polyline = polyline
708        self.verbose = verbose
709       
710        # Find all intersections and associated triangles.
711        self.segments = self.domain.get_intersecting_segments(self.polyline,
712                                                              use_cache=True,
713                                                              verbose=self.verbose)
714       
715        # Get midpoints
716        self.midpoints = segment_midpoints(self.segments)
717
718        # Make midpoints Geospatial instances
719        self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
720
721    def set_verbose(self,verbose=True):
722        """Set verbose mode true or flase"""
723
724        self.verbose=verbose
725
726    def get_flow_through_cross_section(self):
727        """ Output: Total flow [m^3/s] across cross section.
728        """
729
730        # Get interpolated values
731        xmomentum = self.domain.get_quantity('xmomentum')
732        ymomentum = self.domain.get_quantity('ymomentum')
733
734        uh = xmomentum.get_values(interpolation_points=self.midpoints,
735                                  use_cache=True)
736        vh = ymomentum.get_values(interpolation_points=self.midpoints,
737                                  use_cache=True)
738
739        # Compute and sum flows across each segment
740        total_flow = 0
741        for i in range(len(uh)):
742            # Inner product of momentum vector with segment normal [m^2/s]
743            normal = self.segments[i].normal
744            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
745
746            # Flow across this segment [m^3/s]
747            segment_flow = normal_momentum*self.segments[i].length
748
749            # Accumulate
750            total_flow += segment_flow
751
752        return total_flow
753 
754
755    def get_energy_through_cross_section(self, kind='total'):
756        """Obtain average energy head [m] across specified cross section.
757
758        Output:
759            E: Average energy [m] across given segments for all stored times.
760
761        The average velocity is computed for each triangle intersected by
762        the polyline and averaged weighted by segment lengths.
763
764        The typical usage of this function would be to get average energy of
765        flow in a channel, and the polyline would then be a cross section
766        perpendicular to the flow.
767
768        #FIXME (Ole) - need name for this energy reflecting that its dimension
769        is [m].
770        """
771
772        from anuga.config import g, epsilon, velocity_protection as h0
773       
774        # Get interpolated values
775        stage = self.domain.get_quantity('stage')
776        elevation = self.domain.get_quantity('elevation')
777        xmomentum = self.domain.get_quantity('xmomentum')
778        ymomentum = self.domain.get_quantity('ymomentum')
779
780        w = stage.get_values(interpolation_points=self.midpoints, use_cache=True)
781        z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True)
782        uh = xmomentum.get_values(interpolation_points=self.midpoints,
783                                  use_cache=True)
784        vh = ymomentum.get_values(interpolation_points=self.midpoints,
785                                  use_cache=True)
786        h = w-z                # Depth
787
788        # Compute total length of polyline for use with weighted averages
789        total_line_length = 0.0
790        for segment in self.segments:
791            total_line_length += segment.length
792
793        # Compute and sum flows across each segment
794        average_energy = 0.0
795        for i in range(len(w)):
796            # Average velocity across this segment
797            if h[i] > epsilon:
798                # Use protection against degenerate velocities
799                u = uh[i]/(h[i] + h0/h[i])
800                v = vh[i]/(h[i] + h0/h[i])
801            else:
802                u = v = 0.0
803
804            speed_squared = u*u + v*v
805            kinetic_energy = 0.5*speed_squared/g
806
807            if kind == 'specific':
808                segment_energy = h[i] + kinetic_energy
809            elif kind == 'total':
810                segment_energy = w[i] + kinetic_energy
811            else:
812                msg = 'Energy kind must be either "specific" or "total".'
813                msg += ' I got %s' %kind
814
815            # Add to weighted average
816            weigth = self.segments[i].length/total_line_length
817            average_energy += segment_energy*weigth
818
819        return average_energy
820
821class Barometric_pressure:
822    """ Apply barometric pressure stress to water momentum in terms of
823        barometric pressure p [hPa]. If the pressure data is stored in a file
824        file_function is used to create a callable function. The data file
825        contains pressure values at a set of possibly arbitrarily located nodes
826        at a set o possibly irregular but increasing times. file_function
827        interpolates from the file data onto the vertices of the domain.mesh
828        for each time. The file_function is called at every timestep during
829        the evolve function call.
830    """
831    def __init__(self, *args, **kwargs):
832        """Initialise barometric pressure field from barometric pressure [hPa]
833        Input p can be either scalars or Python functions, e.g.
834
835        P = barometric_pressure(1000)
836
837        Arguments can also be Python functions of t,x,y as in
838
839        def pressure(t,x,y):
840            ...
841            return p
842
843        where x and y are vectors.
844
845        and then pass the functions in
846
847        P = Barometric_pressure(pressure)
848
849        agruments can also be the ANGUA file_function, e.g.
850        F = file_function(sww_filename,domain,quantities,interpolation_points)
851        The interpolation_points must be the mesh vertices returned by
852        domain.get_nodes(). Quantities = ['barometric_pressure']
853
854        The file_function is passed using
855
856        P = Barometric_pressure(F, use_coordinates=True/False)
857
858        The instantiated object P can be appended to the list of
859        forcing_terms as in
860
861        domain.forcing_terms.append(P)
862        """
863
864        from anuga.config import rho_a, rho_w, eta_w
865
866        self.use_coordinates=True
867        if len(args) == 1:
868            if ( not callable(args[0]) ):
869                pressure=args[0]
870            else:
871                # Assume vector function returning (pressure)(t,x,y)
872                vector_function = args[0]
873                if ( len(kwargs)==1 ):
874                    self.use_coordinates=kwargs['use_coordinates']
875                else:
876                    self.use_coordinates=True
877
878                if ( self.use_coordinates ):
879                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
880                else:
881                    p = lambda t,i: vector_function(t,point_id=i)[0]
882        else:
883           # Assume info is in 1 or 2 keyword arguments
884           if ( len(kwargs) == 1 ):
885               p = kwargs['p']
886           elif ( len(kwargs)==2 ):
887               p = kwargs['p']
888               self.use_coordinates = kwargs['use_coordinates']
889           else:
890               raise Exception('Assumes one keyword argument: p=... or two '
891                               'keyword arguments p=...,use_coordinates=...')
892
893        if ( self.use_coordinates ):
894            self.pressure = check_forcefield(p)
895        else:
896            self.pressure = p
897
898    def __call__(self, domain):
899        """Evaluate pressure field based on values found in domain"""
900
901        xmom_update = domain.quantities['xmomentum'].explicit_update
902        ymom_update = domain.quantities['ymomentum'].explicit_update
903
904        N = domain.get_number_of_nodes()
905        t = domain.time
906
907        if callable(self.pressure):
908            xv = domain.get_nodes()
909            if ( self.use_coordinates ):
910                p_vec = self.pressure(t, xv[:,0], xv[:,1])
911            else:
912                p_vec=num.empty(N,num.float)
913                for i in range(N):
914                    p_vec[i]=self.pressure(t,i)
915        else:
916            # Assume s is a scalar
917            try:
918                p_vec = self.pressure * num.ones(N, num.float)
919            except:
920                msg = 'Pressure must be either callable or a scalar: %s' %self.s
921                raise Exception(msg)
922
923        stage = domain.quantities['stage']
924        elevation = domain.quantities['elevation']
925
926        #FIXME SR Should avoid allocating memory!
927        height = stage.centroid_values - elevation.centroid_values
928
929        point = domain.get_vertex_coordinates()
930
931        assign_pressure_field_values(height, p_vec, point, domain.triangles,
932                                     xmom_update, ymom_update)
933
934
935def assign_pressure_field_values(height, pressure, x, triangles, 
936                                 xmom_update, ymom_update):
937    """Python version of assigning pressure field to update vectors.
938    """
939
940    from utilities.numerical_tools import gradient
941    from anuga.config import rho_a, rho_w, eta_w
942
943    N = len(height)
944    for k in range(N):
945
946        # Compute pressure slope
947
948        p0 = pressure[triangles[k][0]]
949        p1 = pressure[triangles[k][1]]
950        p2 = pressure[triangles[k][2]]
951
952        k3=3*k
953        x0 = x[k3 + 0][0]
954        y0 = x[k3 + 0][1]
955        x1 = x[k3 + 1][0]
956        y1 = x[k3 + 1][1]
957        x2 = x[k3 + 2][0]
958        y2 = x[k3 + 2][1]
959
960        px,py=gradient(x0, y0, x1, y1, x2, y2, p0, p1, p2)
961
962        xmom_update[k] += height[k]*px/rho_w
963        ymom_update[k] += height[k]*py/rho_w
964
965
966class Barometric_pressure_fast:
967    """ Apply barometric pressure stress to water momentum in terms of
968        barometric pressure p [hPa]. If the pressure data is stored in a file
969        file_function is used to create a callable function. The data file
970        contains pressure values at a set of possibly arbitrarily located nodes
971        at a set o possibly irregular but increasing times. file_function
972        interpolates from the file data onto the vertices of the domain.mesh
973        for each time. Two arrays are then stored p0=p(t0,:) and p1=p(t1,:)
974        where t0<=domain.time<=t1. These arrays are recalculated when necessary
975        i.e t>t1. A linear temporal interpolation is used to approximate
976        pressure at time t.
977    """
978    def __init__(self, *args, **kwargs):
979        """Initialise barometric pressure field from barometric pressure [hPa]
980        Input p can be either scalars or Python functions, e.g.
981
982        P = barometric_pressure(1000)
983
984        Arguments can also be Python functions of t,x,y as in
985
986        def pressure(t,x,y):
987            ...
988            return p
989
990        where x and y are vectors.
991
992        and then pass the functions in
993
994        P = Barometric_pressure(pressure)
995
996        Agruments can also be the ANGUA file_function, e.g.
997        F = file_function(sww_filename,domain,quantities,interpolation_points)
998        The interpolation_points must be the mesh vertices returned by
999        domain.get_nodes(). Quantities = ['barometric_pressure']
1000
1001        The file_function is passed using
1002
1003        P = Barometric_pressure(F, filename=swwname, domain=domain)
1004
1005        The instantiated object P can be appended to the list of
1006        forcing_terms as in
1007
1008        domain.forcing_terms.append(P)
1009        """
1010
1011        from anuga.config import rho_a, rho_w, eta_w
1012
1013        self.use_coordinates=True
1014        if len(args) == 1:
1015            if ( not callable(args[0]) ):
1016                pressure=args[0]
1017            else:
1018                # Assume vector function returning (pressure)(t,x,y)
1019                vector_function = args[0]
1020                if ( len(kwargs)==0 ):
1021                    self.usre_coordinates=True
1022                elif (len(kwargs)==2):
1023                    filename=kwargs['filename']
1024                    domain=kwargs['domain']
1025                    self.use_coordinates=False
1026                else:
1027                    raise Exception('Assumes zero or two keyword arguments '
1028                                    'filename=...,domain=...')
1029
1030                if ( self.use_coordinates ):
1031                    p = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1032                else:
1033                    p = lambda t,i: vector_function(t,point_id=i)[0]
1034        else:
1035           # Assume info is in 1 or 2 keyword arguments
1036           if ( len(kwargs) == 1 ):
1037               p = kwargs['p']
1038               self.use_coordinates=True
1039           elif ( len(kwargs)==3 ):
1040               p = kwargs['p']
1041               filename = kwargs['filename']
1042               domain = kwargs['domain']
1043               self.use_coordinates = False
1044           else:
1045               raise Exception('Assumes one keyword argument: p=f(t,x,y,) or '
1046                               'three keyword arguments '
1047                               'p=f(t,i),filename=...,domain=...')
1048
1049        if ( self.use_coordinates ):
1050            self.pressure = check_forcefield(p)
1051        else:
1052            self.pressure = p
1053
1054        if ( callable(self.pressure) and not self.use_coordinates):
1055
1056            # Open NetCDF file
1057            fid = NetCDFFile(filename, netcdf_mode_r)
1058            self.file_time = fid.variables['time'][:]
1059            fid.close()
1060
1061            msg = 'pressure_file.starttime > domain.starttime'
1062            if (self.file_time[0]>domain.starttime):
1063                raise Exception(msg)
1064
1065            msg = 'pressure_file[-1] < domain.starttime'
1066            if (self.file_time[-1]<domain.starttime):
1067                raise Exception(msg)
1068
1069            msg = 'No pressure values exist for times greater than domain.starttime'
1070            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
1071                raise Exception, msg
1072
1073            # FIXME(JJ): How do we check that evolve
1074            # finaltime  < pressure_file.finaltime     
1075           
1076
1077            self.index=0;
1078            for i in range(len(self.file_time)):
1079                if (self.file_time[i]<domain.starttime):
1080                    self.index=i
1081                else:
1082                    break
1083
1084            N = domain.get_number_of_nodes()
1085            self.prev_pressure_vertex_values=num.empty(N,num.float)
1086            self.next_pressure_vertex_values=num.empty(N,num.float)
1087            for i in range(N):
1088                self.prev_pressure_vertex_values[i]=self.pressure(self.file_time[self.index],i)
1089                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i)
1090
1091        self.p_vec=num.empty(N,num.float)
1092
1093
1094    def __call__(self, domain):
1095        """Evaluate pressure field based on values found in domain"""
1096
1097        xmom_update = domain.quantities['xmomentum'].explicit_update
1098        ymom_update = domain.quantities['ymomentum'].explicit_update
1099
1100        t = domain.time
1101
1102        if callable(self.pressure):
1103            if ( self.use_coordinates ):
1104                xv = domain.get_nodes()
1105                self.p_vec = self.pressure(t, xv[:,0], xv[:,1])
1106            else:
1107                self.update_stored_pressure_values(domain)
1108
1109                # Linear temporal interpolation of pressure values
1110                ratio = (t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index])
1111                self.p_vec = self.prev_pressure_vertex_values + ratio*(self.next_pressure_vertex_values - self.prev_pressure_vertex_values)
1112
1113        else:
1114            # Assume s is a scalar
1115            try:
1116                self.p_vec[:] = self.pressure
1117            except:
1118                msg = 'Pressure must be either callable function or a scalar: %s' %self.s
1119                raise Exception(msg)
1120
1121        stage = domain.quantities['stage']
1122        elevation = domain.quantities['elevation']
1123
1124        height = stage.centroid_values - elevation.centroid_values
1125
1126        point = domain.get_vertex_coordinates()
1127
1128        assign_pressure_field_values(height, self.p_vec, point, 
1129                                     domain.triangles,
1130                                     xmom_update, ymom_update)
1131
1132    def update_stored_pressure_values(self,domain):
1133        while (self.file_time[self.index+1]<domain.time):
1134            self.index+=1
1135            self.prev_pressure_vertex_values=copy(self.next_pressure_vertex_values)
1136            for i in range(self.prev_pressure_vertex_values.shape[0]):
1137                self.next_pressure_vertex_values[i]=self.pressure(self.file_time[self.index+1],i) 
1138
1139
1140class Wind_stress_fast:
1141    """ Apply wind stress to water momentum in terms of
1142        wind speed [m/s] and wind direction [degrees].
1143        If the wind data is stored in a file
1144        file_function is used to create a callable function. The data file
1145        contains wind speed and direction values at a set of possibly
1146        arbitrarily located nodes
1147        at a set of possibly irregular but increasing times. file_function
1148        interpolates from the file data onto the centroids of the domain.mesh
1149        for each time. Two arrays for each wind quantity are then stored \
1150        q0=q(t0,:) and q1=q(t1,:)
1151        where t0<=domain.time<=t1. These arrays are recalculated when necessary
1152        i.e t>t1. A linear temporal interpolation is used to approximate
1153        pressure at time t.
1154    """
1155    def __init__(self, *args, **kwargs):
1156        """Initialise windfield from wind speed s [m/s]
1157        and wind direction phi [degrees]
1158
1159        Inputs v and phi can be either scalars or Python functions, e.g.
1160
1161        W = Wind_stress(10, 178)
1162
1163        #FIXME - 'normal' degrees are assumed for now, i.e. the
1164        vector (1,0) has zero degrees.
1165        We may need to convert from 'compass' degrees later on and also
1166        map from True north to grid north.
1167
1168        Arguments can also be Python functions of t,x,y as in
1169
1170        def speed(t,x,y):
1171            ...
1172            return s
1173
1174        def angle(t,x,y):
1175            ...
1176            return phi
1177
1178        where x and y are vectors.
1179
1180        and then pass the functions in
1181
1182        W = Wind_stress(speed, angle)
1183
1184        The instantiated object W can be appended to the list of
1185        forcing_terms as in
1186
1187        Alternatively, one vector valued function for (speed, angle)
1188        can be applied, providing both quantities simultaneously.
1189        As in
1190        W = Wind_stress(F), where returns (speed, angle) for each t.
1191
1192        domain.forcing_terms.append(W)
1193        """
1194
1195        from anuga.config import rho_a, rho_w, eta_w
1196
1197        self.use_coordinates=True
1198        if len(args) == 2:
1199            s = args[0]
1200            phi = args[1]
1201        elif len(args) == 1:
1202            # Assume vector function returning (s, phi)(t,x,y)
1203            vector_function = args[0]
1204            if ( len(kwargs)==2 ):
1205                filename=kwargs['filename']
1206                domain=kwargs['domain']
1207                self.use_coordinates=False
1208            else:
1209                self.use_coordinates=True
1210            if ( self.use_coordinates ):
1211                s = lambda t,x,y: vector_function(t,x=x,y=y)[0]
1212                phi = lambda t,x,y: vector_function(t,x=x,y=y)[1]
1213            else:
1214                s = lambda t,i: vector_function(t,point_id=i)[0]
1215                phi = lambda t,i: vector_function(t,point_id=i)[1]
1216        else:
1217           # Assume info is in 2 keyword arguments
1218           if len(kwargs) == 2:
1219               s = kwargs['s']
1220               phi = kwargs['phi']
1221           else:
1222               raise Exception('Assumes two keyword arguments: s=...,phi=....')
1223
1224        if ( self.use_coordinates ):
1225            self.speed = check_forcefield(s)
1226            self.phi = check_forcefield(phi)
1227        else:
1228            self.speed = s
1229            self.phi = phi
1230
1231        N = len(domain)
1232        if ( not self.use_coordinates):
1233
1234            # Open NetCDF file
1235            fid = NetCDFFile(filename, netcdf_mode_r)
1236            self.file_time = fid.variables['time'][:]
1237            fid.close()
1238
1239            msg = 'wind_file.starttime > domain.starttime'
1240            if (self.file_time[0]>domain.starttime):
1241                raise Exception(msg)
1242
1243            msg = 'wind_file[-1] < domain.starttime'
1244            if (self.file_time[-1]<domain.starttime):
1245                raise Exception(msg)
1246
1247            msg = 'No wind values exist for times greater than domain.starttime'
1248            if (self.file_time[-2]<domain.starttime and self.file_time[-1]>domain.starttime):
1249                raise Exception(msg)
1250
1251            # FIXME(JJ): How do we check that evolve
1252            # finaltime  < wind_file.finaltime     
1253           
1254
1255            self.index=0;
1256            for i in range(len(self.file_time)):
1257                if (self.file_time[i]<domain.starttime):
1258                    self.index=i
1259                else:
1260                    break
1261
1262            self.prev_windspeed_centroid_values=num.empty(N,num.float)
1263            self.next_windspeed_centroid_values=num.empty(N,num.float)
1264            self.prev_windangle_centroid_values=num.empty(N,num.float)
1265            self.next_windangle_centroid_values=num.empty(N,num.float)
1266            for i in range(N):
1267                self.prev_windspeed_centroid_values[i]=self.speed(self.file_time[self.index],i)
1268                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i)
1269                self.prev_windangle_centroid_values[i]=self.phi(self.file_time[self.index],i)
1270                self.next_windangle_centroid_values[i]=self.phi(self.file_time[self.index+1],i)
1271
1272        self.s_vec=num.empty(N,num.float)
1273        self.phi_vec=num.empty(N,num.float)
1274
1275        self.const = eta_w*rho_a/rho_w
1276
1277    def __call__(self, domain):
1278        """Evaluate windfield based on values found in domain"""
1279
1280        xmom_update = domain.quantities['xmomentum'].explicit_update
1281        ymom_update = domain.quantities['ymomentum'].explicit_update
1282
1283        N = len(domain)    # number_of_triangles
1284        t = domain.time
1285
1286        if callable(self.speed):
1287            if ( self.use_coordinates ):
1288                xc = domain.get_centroid_coordinates()
1289                self.s_vec = self.speed(t, xc[:,0], xc[:,1])
1290            else:
1291                self.update_stored_wind_values(domain)
1292
1293                # Linear temporal interpolation of wind values
1294                if t==self.file_time[self.index]:
1295                    ratio = 0.
1296                else:
1297                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
1298                self.s_vec = self.prev_windspeed_centroid_values + ratio*(self.next_windspeed_centroid_values - self.prev_windspeed_centroid_values)
1299        else:
1300            # Assume s is a scalar
1301            try:
1302                self.s_vec[:] = self.speed
1303            except:
1304                msg = 'Speed must be either callable or a scalar: %s' %self.s
1305                raise Exception(msg)
1306
1307        if callable(self.phi):
1308            if ( self.use_coordinates ):
1309                xc = domain.get_centroid_coordinates()
1310                self.phi_vec = self.phi(t, xc[:,0], xc[:,1])
1311            else:
1312                self.update_stored_wind_values(domain)
1313
1314                # Linear temporal interpolation of wind values
1315                if t==self.file_time[self.index]:
1316                    ratio = 0.
1317                else:
1318                    ratio = ((t - self.file_time[self.index]) / (self.file_time[self.index+1]-self.file_time[self.index]))
1319                self.phi_vec = self.prev_windangle_centroid_values + ratio*(self.next_windangle_centroid_values - self.prev_windangle_centroid_values)
1320        else:
1321            # Assume phi is a scalar
1322
1323            try:
1324                self.phi_vec[:] = self.phi
1325            except:
1326                msg = 'Angle must be either callable or a scalar: %s' %self.phi
1327                raise Exception(msg)
1328
1329        assign_windfield_values(xmom_update, ymom_update,
1330                                self.s_vec, self.phi_vec, self.const)
1331
1332    def update_stored_wind_values(self,domain):
1333        while (self.file_time[self.index+1]<domain.time):
1334            self.index+=1
1335            self.prev_windspeed_centroid_values=copy(self.next_windspeed_centroid_values)
1336            self.prev_windangle_centroid_values=copy(self.next_windangle_centroid_values)
1337            for i in range(self.next_windspeed_centroid_values.shape[0]):
1338                self.next_windspeed_centroid_values[i]=self.speed(self.file_time[self.index+1],i) 
1339                self.next_windangle_centroid_values[i]=self.phi(self.file_time[self.index+1],i) 
Note: See TracBrowser for help on using the repository browser.