Ignore:
Timestamp:
Jun 14, 2012, 12:31:32 AM (13 years ago)
Author:
davies
Message:

bal_dev: Updates

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_boundary_conditions.py

    r8396 r8446  
    66#####
    77
    8 class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
    9     """ Boundary which operates directly on the fluxes
    10         CAN BE UNSTABLE -- CAREFUL
    11         Sets boundary values of h, uh, vh as in transmissive boundary, but
    12         also ensures:
    13         flux[0] = 0.0,
     8#class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
     9#    """ Boundary which operates directly on the fluxes
     10#        CAN BE UNSTABLE -- CAREFUL
     11#        Sets boundary values of h, uh, vh as in transmissive boundary, but
     12#        also ensures:
     13#        flux[0] = 0.0,
     14#    """
     15#    def __init__(self, domain = None):
     16#        Boundary.__init__(self)
     17#
     18#        if domain is None:
     19#            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
     20#            raise Exception, msg
     21#
     22#        self.domain = domain
     23#
     24#    def __repr__(self):
     25#        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
     26#
     27#    def evaluate(self, vol_id, edge_id):
     28#        """Transmissive boundaries return the edge values
     29#        of the volume they serve.
     30#        """
     31#
     32#        if self.domain.get_centroid_transmissive_bc() :
     33#            q = self.domain.get_evolved_quantities(vol_id)
     34#        else:
     35#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     36#        return q
     37#
     38######
     39#
     40#class  zero_mass_flux_zero_momentum_boundary(Boundary):
     41#    """ Boundary which operates directly on the fluxes
     42#        Sets boundary values of h=neighbour_h, uh=0, vh=0.
     43#        and ensure that:
     44#        flux[0] = 0.0
     45#    """
     46#    def __init__(self, domain = None):
     47#        Boundary.__init__(self)
     48#
     49#        if domain is None:
     50#            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
     51#            raise Exception, msg
     52#
     53#        self.domain = domain
     54#
     55#    def __repr__(self):
     56#        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
     57#
     58#    def evaluate(self, vol_id, edge_id):
     59#        """ Apply chosen values to q[1], q[2]
     60#        """
     61#
     62#        if self.domain.get_centroid_transmissive_bc() :
     63#            q = self.domain.get_evolved_quantities(vol_id)
     64#        else:
     65#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     66#
     67#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     68#
     69#        r = rotate(q, normal, direction = 1)
     70#        #r[1] = -r[1]
     71#        tdamp=0.0
     72#        ndamp=0.0
     73#        r[1]=r[1]*ndamp
     74#        r[2]=r[2]*tdamp
     75#        q = rotate(r, normal, direction = -1)
     76#
     77#        # q[1]=0.
     78#        # q[2]=0.
     79#
     80#        return q
     81#
     82#
     83#class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
     84#    """ Boundary which operates directly on the fluxes
     85#        Sets boundary values of h=neighbour_h, uh=0, vh=vh
     86#        and ensure that:
     87#        flux[0] = 0.0
     88#    """
     89#    def __init__(self, domain = None):
     90#        Boundary.__init__(self)
     91#
     92#        if domain is None:
     93#            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
     94#            raise Exception, msg
     95#
     96#        self.domain = domain
     97#
     98#    def __repr__(self):
     99#        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
     100#
     101#    def evaluate(self, vol_id, edge_id):
     102#        """ Apply chosen values to q[1], q[2]
     103#        """
     104#
     105#        if self.domain.get_centroid_transmissive_bc() :
     106#            q = self.domain.get_evolved_quantities(vol_id)
     107#        else:
     108#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     109#
     110#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     111#
     112#        r = rotate(q, normal, direction = 1)
     113#        #r[1] = -r[1]
     114#        tdamp=1.0
     115#        ndamp=0.0
     116#        r[1]=r[1]*ndamp
     117#        r[2]=r[2]*tdamp
     118#        q = rotate(r, normal, direction = -1)
     119#
     120#        # q[1]=0.
     121#        # q[2]=0.
     122#
     123#        return q
     124#
     125#
     126#class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
     127#    """ Boundary which operates directly on the fluxes
     128#        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
     129#        and ensure that:
     130#        flux[0] = 0.0
     131#    """
     132#    def __init__(self, domain = None):
     133#        Boundary.__init__(self)
     134#
     135#        if domain is None:
     136#            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
     137#            raise Exception, msg
     138#
     139#        self.domain = domain
     140#
     141#    def __repr__(self):
     142#        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
     143#
     144#    def evaluate(self, vol_id, edge_id):
     145#        """ Apply chosen values to q[1], q[2]
     146#        """
     147#
     148#        if self.domain.get_centroid_transmissive_bc() :
     149#            q = self.domain.get_evolved_quantities(vol_id)
     150#        else:
     151#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     152#
     153#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     154#
     155#        r = rotate(q, normal, direction = 1)
     156#        #r[1] = -r[1]
     157#        tdamp=0.0
     158#        ndamp=1.0
     159#        r[1]=r[1]*ndamp
     160#        r[2]=r[2]*tdamp
     161#        q = rotate(r, normal, direction = -1)
     162#
     163#        # q[1]=0.
     164#        # q[2]=0.
     165#
     166#        return q
     167#
     168#
     169class Characteristic_boundary(Boundary):
    14170    """
    15     def __init__(self, domain = None):
     171
     172    Example:
     173
     174    def waveform(t):
     175        return sea_level + normalized_amplitude/cosh(t-25)**2
     176
     177    Bf = Nudge_boundary(domain, waveform)
     178
     179    Underlying domain must be specified when boundary is instantiated
     180    """
     181
     182    def __init__(self, domain=None, function=None):
     183        """ Instantiate a
     184            Nudge_boundary.
     185            domain is the domain containing the boundary
     186            function is the function to apply
     187        """
     188
    16189        Boundary.__init__(self)
    17190
    18191        if domain is None:
    19             msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
     192            msg = 'Domain must be specified for this type boundary'
    20193            raise Exception, msg
    21194
     195        if function is None:
     196            msg = 'Function must be specified for this type boundary'
     197            raise Exception, msg
     198
    22199        self.domain = domain
     200        self.function = function
     201
    23202
    24203    def __repr__(self):
    25         return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
     204        """ Return a representation of this instance. """
     205        msg = 'Nudge_boundary'
     206        msg += '(%s)' % self.domain
     207        return msg
     208
    26209
    27210    def evaluate(self, vol_id, edge_id):
    28         """Transmissive boundaries return the edge values
    29         of the volume they serve.
    30211        """
    31 
    32         if self.domain.get_centroid_transmissive_bc() :
    33             q = self.domain.get_evolved_quantities(vol_id)
     212        """
     213
     214        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
     215        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
     216        depth=max(q[0]-bed,0.0)
     217        dt=self.domain.timestep
     218
     219        normal = self.domain.get_normal(vol_id, edge_id)
     220
     221
     222        t = self.domain.get_time()
     223
     224        if hasattr(self.function, 'time'):
     225            # Roll boundary over if time exceeds
     226            while t > self.function.time[-1]:
     227                msg = 'WARNING: domain time %.2f has exceeded' % t
     228                msg += 'time provided in '
     229                msg += 'Flather_boundary object.\n'
     230                msg += 'I will continue, reusing the object from t==0'
     231                log.critical(msg)
     232                t -= self.function.time[-1]
     233
     234        value = self.function(t)
     235        try:
     236            x = float(value)
     237        except:
     238            x = float(value[0])
     239
     240        ##
     241        ndotq = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
     242
     243        if(depth==0.):
     244            q[0] = x
     245            q[1] = 0.
     246            q[2] = 0.
     247
    34248        else:
    35             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     249
     250            # Asssume sub-critical flow. Set the values of the characteristics as
     251            # appropriate, depending on whether we have inflow or outflow
     252            tmp = (9.8/depth)**0.5
     253            if(ndotq>0.):
     254                # Only need to impose one characteristic
     255                # Assume zero external velocity (dangerous?)
     256                # FIXME: Allow external velocities to be specified in the function
     257                # This would be needed to use the boundary for real inflows (e.g. run_wave.py)
     258                # At the moment, it seems to work okay as a 'passive' outflow boundary (e.g. run_wave.py)
     259                # Here it attenuates the wave at the outflow a bit, but doesn't do too much damage compared
     260                # with the other boundary conditions
     261                w1 = 0. - tmp*x
     262                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
     263                w3 = ndotq/depth + tmp*q[0]
     264               
     265            else:
     266                # Need to set 2 characteristics
     267                # Assume zero external velocity (dangerous?)
     268                # FIXME: Allow this to be specified with a function
     269                w1 = 0. - tmp*x
     270                w2 = 0.
     271                #w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
     272                w3 = ndotq/depth + tmp*q[0]
     273
     274
     275            q[0] = (w3-w1)/(2*tmp)
     276            qperp= (w3+w1)/2. *depth
     277            qpar=  w2*depth
     278
     279            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
     280
     281            q[1] = qperp*normal[0] + qpar*normal[1]
     282            q[2] = qperp*normal[1] -qpar*normal[0]
     283
    36284        return q
    37 
    38 #####
    39 
    40 class  zero_mass_flux_zero_momentum_boundary(Boundary):
    41     """ Boundary which operates directly on the fluxes
    42         Sets boundary values of h=neighbour_h, uh=0, vh=0.
    43         and ensure that:
    44         flux[0] = 0.0
    45     """
    46     def __init__(self, domain = None):
    47         Boundary.__init__(self)
    48 
    49         if domain is None:
    50             msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
    51             raise Exception, msg
    52 
    53         self.domain = domain
    54 
    55     def __repr__(self):
    56         return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
    57 
    58     def evaluate(self, vol_id, edge_id):
    59         """ Apply chosen values to q[1], q[2]
    60         """
    61 
    62         if self.domain.get_centroid_transmissive_bc() :
    63             q = self.domain.get_evolved_quantities(vol_id)
    64         else:
    65             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    66 
    67         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    68 
    69         r = rotate(q, normal, direction = 1)
    70         #r[1] = -r[1]
    71         tdamp=0.0
    72         ndamp=0.0
    73         r[1]=r[1]*ndamp
    74         r[2]=r[2]*tdamp
    75         q = rotate(r, normal, direction = -1)
    76 
    77         # q[1]=0.
    78         # q[2]=0.
    79 
    80         return q
    81 
    82 
    83 class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
    84     """ Boundary which operates directly on the fluxes
    85         Sets boundary values of h=neighbour_h, uh=0, vh=vh
    86         and ensure that:
    87         flux[0] = 0.0
    88     """
    89     def __init__(self, domain = None):
    90         Boundary.__init__(self)
    91 
    92         if domain is None:
    93             msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
    94             raise Exception, msg
    95 
    96         self.domain = domain
    97 
    98     def __repr__(self):
    99         return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
    100 
    101     def evaluate(self, vol_id, edge_id):
    102         """ Apply chosen values to q[1], q[2]
    103         """
    104 
    105         if self.domain.get_centroid_transmissive_bc() :
    106             q = self.domain.get_evolved_quantities(vol_id)
    107         else:
    108             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    109 
    110         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    111 
    112         r = rotate(q, normal, direction = 1)
    113         #r[1] = -r[1]
    114         tdamp=1.0
    115         ndamp=0.0
    116         r[1]=r[1]*ndamp
    117         r[2]=r[2]*tdamp
    118         q = rotate(r, normal, direction = -1)
    119 
    120         # q[1]=0.
    121         # q[2]=0.
    122 
    123         return q
    124 
    125 
    126 class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
    127     """ Boundary which operates directly on the fluxes
    128         Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
    129         and ensure that:
    130         flux[0] = 0.0
    131     """
    132     def __init__(self, domain = None):
    133         Boundary.__init__(self)
    134 
    135         if domain is None:
    136             msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
    137             raise Exception, msg
    138 
    139         self.domain = domain
    140 
    141     def __repr__(self):
    142         return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
    143 
    144     def evaluate(self, vol_id, edge_id):
    145         """ Apply chosen values to q[1], q[2]
    146         """
    147 
    148         if self.domain.get_centroid_transmissive_bc() :
    149             q = self.domain.get_evolved_quantities(vol_id)
    150         else:
    151             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    152 
    153         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    154 
    155         r = rotate(q, normal, direction = 1)
    156         #r[1] = -r[1]
    157         tdamp=0.0
    158         ndamp=1.0
    159         r[1]=r[1]*ndamp
    160         r[2]=r[2]*tdamp
    161         q = rotate(r, normal, direction = -1)
    162 
    163         # q[1]=0.
    164         # q[2]=0.
    165 
    166         return q
    167 
    168 
    169285###
    170286def rotate(q,normal,direction=1):
Note: See TracChangeset for help on using the changeset viewer.