Changeset 9681 for trunk/anuga_core


Ignore:
Timestamp:
Feb 23, 2015, 9:21:46 AM (10 years ago)
Author:
davies
Message:

Major changes + bugfix for internal_boundary_operator and related

Location:
trunk/anuga_core/anuga
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/anuga/parallel/parallel_boyd_box_operator.py

    r9657 r9681  
    6363                                          use_momentum_jet=use_momentum_jet,
    6464                                          zero_outflow_momentum=(not use_momentum_jet),
    65                                           use_old_momentum_method=True,
    66                                           force_constant_inlet_elevations=False,
     65                                          use_old_momentum_method=True,
     66                                          always_use_Q_wetdry_adjustment=True,
     67                                          force_constant_inlet_elevations=False,
    6768                                          description=description,
    6869                                          label=label,

  • trunk/anuga_core/anuga/parallel/parallel_boyd_pipe_operator.py

    r9657 r9681  
    6161                                          zero_outflow_momentum=(not use_momentum_jet),
    6262                                          use_old_momentum_method=True,
     63                                          always_use_Q_wetdry_adjustment=True,
    6364                                          force_constant_inlet_elevations=False,
    6465                                          description=description,

  • trunk/anuga_core/anuga/parallel/parallel_internal_boundary_operator.py

    r9672 r9681  
    33import numpy
    44from numpy.linalg import solve
     5import scipy
     6import scipy.optimize as sco
    57
    68#from anuga.structures.boyd_box_operator import boyd_box_function
    … …  
    6870                                          zero_outflow_momentum=zero_outflow_momentum,
    6971                                          use_old_momentum_method=False,
     72                                          always_use_Q_wetdry_adjustment=False,
    7073                                          force_constant_inlet_elevations=force_constant_inlet_elevations,
    7174                                          description=description,
    … …  
    288291        """
    289292            Uses semi-implicit discharge estimation:
    290               Discharge = 0.5*(Q(H0, T0) + Q(H0 + delta_H, T0+delta_T))
     293              Discharge = (1-theta)*Q(H0, T0) + theta*Q(H0 + delta_H, T0+delta_T))
    291294            where H0 = headwater stage, T0 = tailwater stage, delta_H = change in
    292295            headwater stage over a timestep, delta_T = change in tailwater stage over a
    293             timestep.
    294 
    295             We can estimate delta_H, delta_T by solving the following system (based on mass conservation):
    296               A0*delta_H = -dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
    297               A1*delta_T =  dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
    298             where A0, A1 are the inlet areas, and dt is the timestep
    299 
    300             We linearise the system with the approximation:
    301               Q(H0 + delta_H, T0 + delta_T) ~= Q(H0, T0) + delQ/delH * delta_H + delQ/delT*delta_T
    302             where delQ/delH and delQ/delT are evaluated with numerical finite differences
    303 
     296            timestep, and Q is the discharge function, and theta is a constant in
     297            [0,1] determining the degree of implicitness (currently hardcoded).
     298
     299            Note this is effectively assuming:
     300            1) Q is a function of stage, not energy (so we can relate mass change directly to delta_H, delta_T). We
     301               could generalise it to the energy case ok.
     302            2) The stage is computed on the exchange line (or the change in
     303                stage at the enquiry point is effectively the same as that on the exchange
     304                line)
    304305
    305306        """
    … …  
    390391            Q0 = self.internal_boundary_function(E0, E1)
    391392            dt = self.domain.get_timestep()
     393           
    392394            if dt > 0.:
    393                 # Numerical derivatives of Discharge function
    394                 dE0 = 0.01
    395                 dE1 = 0.01
    396                 dQ_dE0 = (self.internal_boundary_function(E0+dE0, E1) - self.internal_boundary_function(E0-dE0, E1))/(2*dE0)
    397                 dQ_dE1 = (self.internal_boundary_function(E0, E1+dE1) - self.internal_boundary_function(E0, E1-dE1))/(2*dE1)
    398 
    399                 # Assemble and solve matrix
    400                 hdt = 0.5*dt
    401 
    402                 M11 = area0 + dQ_dE0*hdt
    403                 M12 = dQ_dE1*hdt
    404                 M21 = -dQ_dE0*hdt
    405                 M22 = area1 - dQ_dE1*hdt
    406 
    407                 lhs = numpy.array([ [M11, M12], [M21, M22]])
    408                 rhs = numpy.array([ -Q0*dt, Q0*dt])
    409                 # sol contains delta_E0, delta_E1
    410                 sol = solve(lhs, rhs)
    411                
    412                 #Q = 0.5*(Q0 + ( Q0 + sol[0]*dQ_dE0 + sol[1]*dQ_dE1))
    413                 Q1 = self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
    414                 Q = 0.5*(Q0 + Q1)
     395                # Key constants for iterative solution
     396                theta = 1.0
     397                sol = numpy.array([0., 0.]) # estimate of (delta_H, delta_T)
     398                areas = numpy.array([area0, area1])
     399
     400                # Use scipy root finding
     401                def F_to_solve(sol):
     402                    Q1 =  self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
     403                    discharge = (1-theta)*Q0 + theta*Q1
     404                    output = sol*areas - discharge*dt*numpy.array([-1., 1.])
     405                    return(output)
     406
     407                final_sol = sco.root(F_to_solve, sol, method='lm').x
     408                Q1 =  self.internal_boundary_function(E0 + final_sol[0], E1 + final_sol[1])
     409                Q = (1.0-theta)*Q0 + theta*Q1
     410
    415411            else:
    416412                Q = Q0

  • trunk/anuga_core/anuga/parallel/parallel_structure_operator.py

    r9671 r9681  
    5252                 zero_outflow_momentum,
    5353                 use_old_momentum_method,
     54                 always_use_Q_wetdry_adjustment,
    5455                 force_constant_inlet_elevations,
    5556                 description,
    … …  
    125126            raise Exception(msg)
    126127        self.use_old_momentum_method = use_old_momentum_method
     128        self.always_use_Q_wetdry_adjustment = always_use_Q_wetdry_adjustment
    127129
    128130        if description == None:
    … …  
    150152        self.accumulated_flow = 0.0
    151153        self.discharge = 0.0
     154        self.discharge_function_value = 0.0
    152155        self.velocity = 0.0
    153156        self.outlet_depth = 0.0
    … …  
    280283        # Master proc of structure only
    281284        if self.myid == self.master_proc:
    282             #if old_inflow_depth > 0.0 :
    283             #    Q_star = Q/old_inflow_depth
    284             #else:
    285             #    Q_star = 0.0
    286285            if old_inflow_depth > 0.0 :
    287286                dt_Q_on_d = timestep*Q/old_inflow_depth
    … …  
    289288                dt_Q_on_d = 0.0
    290289
     290            # Check whether we should use the wet-dry Q adjustment (where Q is
     291            # multiplied by new_inflow_depth/old_inflow_depth)
     292            always_use_Q_wetdry_adjustment = self.always_use_Q_wetdry_adjustment
     293            # Always use it if we are near wet-dry
     294            use_Q_wetdry_adjustment = ((always_use_Q_wetdry_adjustment) |\
     295                (old_inflow_depth*inflow_area <= Q*timestep))
     296
    291297            factor = 1.0/(1.0 + dt_Q_on_d/inflow_area)
    292             new_inflow_depth = old_inflow_depth*factor
     298       
     299            if use_Q_wetdry_adjustment:
     300                new_inflow_depth = old_inflow_depth*factor
     301                if old_inflow_depth > 0.:
     302                    timestep_star = timestep*new_inflow_depth/old_inflow_depth
     303                else:
     304                    timestep_star = 0.
     305            else:
     306                new_inflow_depth = old_inflow_depth - timestep*Q/inflow_area
     307                timestep_star = timestep
    293308
    294309            #new_inflow_xmom = old_inflow_xmom*factor
    … …  
    301316
    302317            else:
    303                 # For the momentum balance, note that Q also advects the momentum,
    304                 # which has an average value of new_inflow_mom (or old_inflow_mom). For
    305                 # consistency we keep using the (new_inflow_depth/old_inflow_depth)
    306                 # factor for discharge:
     318                # For the momentum balance, note that Q also transports the velocity,
     319                # which has an average value of new_inflow_mom/depth (or old_inflow_mom/depth).
    307320                #
    308321                #     new_inflow_xmom*inflow_area =
    309322                #     old_inflow_xmom*inflow_area -
    310                 #     timestep*Q*(new_inflow_depth/old_inflow_depth)*new_inflow_xmom
     323                #     timestep*Q*(new_inflow_xmom/old_inflow_depth)
    311324                # and:
    312325                #     new_inflow_ymom*inflow_area =
    313326                #     old_inflow_ymom*inflow_area -
    314                 #     timestep*Q*(new_inflow_depth/old_inflow_depth)*new_inflow_ymom
     327                #     timestep*Q*(new_inflow_ymom/old_inflow_depth)
    315328                #
    316                 # The choice of new_inflow_mom in the final term at the end might be
    317                 # replaced with old_inflow_mom
     329                # The choice of new_inflow_mom in the final term might be
     330                # replaced with old_inflow_mom.
    318331                #
    319                 factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/inflow_area)
     332                # The units balance: (m^2/s)*(m^2) = (m^2/s)*(m^2) - s*(m^3/s)*(m^2/s)*(m^(-1))
     333                #
     334                if old_inflow_depth > 0.:
     335                    if use_Q_wetdry_adjustment:
     336                        factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/(old_inflow_depth*inflow_area))
     337                    else:
     338                        factor2 = 1.0/(1.0 + timestep*Q/(old_inflow_depth*inflow_area))
     339                else:
     340                    factor2 = 0.
     341
    320342                new_inflow_xmom = old_inflow_xmom*factor2
    321343                new_inflow_ymom = old_inflow_ymom*factor2
    … …  
    370392
    371393            # set outflow
    372             if old_inflow_depth > 0.0 :
    373                 timestep_star = timestep*new_inflow_depth/old_inflow_depth
    374             else:
    375                 timestep_star = 0.0
    376 
    377394            outflow_extra_depth = Q*timestep_star/outflow_area
    378395            outflow_direction = - outflow_outward_culvert_vector
    … …  
    383400            # Update Stats
    384401            self.discharge  = Q*timestep_star/timestep #outflow_extra_depth*self.outflow.get_area()/timestep
     402            self.discharge_function_value = Q
    385403            self.velocity = barrel_speed #self.discharge/outlet_depth/self.width
    386404
    … …  
    612630            message += 'Type: %s\n' % self.structure_type
    613631            message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
     632            message += 'Discharge function value [m^3/s]: %.2f\n' % self.discharge_function_value
    614633            message += 'Velocity  [m/s]: %.2f\n' % self.velocity
    615634            message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
    … …  
    631650            self.log_filename = self.domain.get_datadir() + '/' + self.label + '.log'
    632651            log_to_file(self.log_filename, stats, mode='w')
    633             log_to_file(self.log_filename, 'time,discharge,velocity,driving_energy,delta_total_energy')
     652            log_to_file(self.log_filename, 'time,discharge,discharge_function_value,velocity,driving_energy,delta_total_energy')
    634653
    635654            #log_to_file(self.log_filename, self.culvert_type)
    … …  
    655674        message  = '%.5f, ' % self.domain.get_time()
    656675        message += '%.5f, ' % self.discharge
     676        message += '%.5f, ' % self.discharge_function_value
    657677        message += '%.5f, ' % self.velocity
    658678        message += '%.5f, ' % self.driving_energy

  • trunk/anuga_core/anuga/parallel/parallel_weir_orifice_trapezoid_operator.py

    r9657 r9681  
    6262                                          zero_outflow_momentum=(not use_momentum_jet),
    6363                                          use_old_momentum_method=True,
     64                                          always_use_Q_wetdry_adjustment=True,
    6465                                          force_constant_inlet_elevations=False,
    6566                                          description=description,

  • trunk/anuga_core/anuga/structures/boyd_box_operator_Amended3.py

    r8125 r9681  
    3232                 verbose=False):
    3333
    34                      
     34        raise Exception('(Feb 2015) This code seems broken -- the structure operator call is incorrect -- retire?')
     35
    3536        anuga.Structure_operator.__init__(self,
    3637                                          domain,

  • trunk/anuga_core/anuga/structures/internal_boundary_operator.py

    r9672 r9681  
    33import numpy
    44from numpy.linalg import solve
     5import scipy
     6import scipy.optimize as sco
    57
    68
    … …  
    7577                                          zero_outflow_momentum=zero_outflow_momentum,
    7678                                          use_old_momentum_method=False,
     79                                          always_use_Q_wetdry_adjustment=False,
    7780                                          force_constant_inlet_elevations=force_constant_inlet_elevations,
    7881                                          description=description,
    … …  
    227230    def discharge_routine_implicit(self):
    228231        """Uses semi-implicit discharge estimation:
    229           Discharge = 0.5*(Q(H0, T0) + Q(H0 + delta_H, T0+delta_T))
     232          Discharge = (1-theta)*Q(H0, T0) + theta*Q(H0 + delta_H, T0+delta_T))
    230233        where H0 = headwater stage, T0 = tailwater stage, delta_H = change in
    231234        headwater stage over a timestep, delta_T = change in tailwater stage over a
    232         timestep, and Q is the discharge function
    233 
    234         We can estimate delta_H, delta_T by solving the following system (based on mass conservation):
    235           A0*delta_H = -dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
    236           A1*delta_T =  dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
    237         where A0, A1 are the inlet areas, and dt is the timestep
    238 
    239         We linearise the system with the approximation:
    240           Q(H0 + delta_H, T0 + delta_T) ~= Q(H0, T0) + delQ/delH * delta_H + delQ/delT*delta_T
    241         where delQ/delH and delQ/delT are evaluated with numerical finite differences
    242 
     235        timestep, and Q is the discharge function, and theta is a constant in
     236        [0,1] determining the degree of implicitness (currently hard-coded).
     237
     238        Note this is effectively assuming:
     239        1) Q is a function of stage, not energy (so we can relate mass change directly to delta_H, delta_T). We
     240           could generalise it to the energy case ok.
     241        2) The stage is computed on the exchange line (or the change in
     242            stage at the enquiry point is effectively the same as that on the exchange
     243            line)
    243244
    244245        """
    … …  
    263264            E0 = self.inlet0_energy
    264265            E1 = self.inlet1_energy
    265 
    266             # Numerical derivatives
    267             dE0 = 0.01
    268             dE1 = 0.01
    269             dQ_dE0 = (self.internal_boundary_function(E0+dE0, E1) - self.internal_boundary_function(E0-dE0, E1))/(2*dE0)
    270             dQ_dE1 = (self.internal_boundary_function(E0, E1+dE1) - self.internal_boundary_function(E0, E1-dE1))/(2*dE1)
    271 
    272             # Assemble and solve matrix
    273266            A0 = self.inlets[0].area
    274267            A1 = self.inlets[1].area
    275             hdt = 0.5*dt
    276 
    277             M11 = A0 + dQ_dE0*hdt
    278             M12 = dQ_dE1*hdt
    279             M21 = -dQ_dE0*hdt
    280             M22 = A1 - dQ_dE1*hdt
    281 
    282             lhs = numpy.array([ [M11, M12], [M21, M22]])
    283             rhs = numpy.array([ -Q0*dt, Q0*dt])
    284             # sol contains delta_E0, delta_E1
    285             sol = solve(lhs, rhs)
    286            
    287             #Q = 0.5*(Q0 + ( Q0 + sol[0]*dQ_dE0 + sol[1]*dQ_dE1))
    288             Q1 =  self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
    289             Q = 0.5*(Q0 + Q1)
     268            theta = 1.0
     269
     270            sol = numpy.array([0., 0.]) # estimate of (delta_H, delta_T)
     271            areas = numpy.array([A0, A1])
     272
     273            # Use scipy root finding
     274            def F_to_solve(sol):
     275                Q1 =  self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
     276                discharge = (1.0-theta)*Q0 + theta*Q1
     277                output = sol*areas - discharge*dt*numpy.array([-1., 1.])
     278                return(output)
     279
     280            final_sol = sco.root(F_to_solve, sol, method='lm').x
     281            Q1 =  self.internal_boundary_function(E0 + final_sol[0], E1 + final_sol[1])
     282            Q = (1.0-theta)*Q0 + theta*Q1
     283
    290284        else:
    291285            Q = Q0

  • trunk/anuga_core/anuga/structures/structure_operator.py

    r9671 r9681  
    3939                 zero_outflow_momentum=True,
    4040                 use_old_momentum_method=True,
     41                 always_use_Q_wetdry_adjustment=True,
    4142                 force_constant_inlet_elevations=False,
    4243                 description=None,
    … …  
    9495            raise Exception(msg)
    9596        self.use_old_momentum_method = use_old_momentum_method
     97        self.always_use_Q_wetdry_adjustment = always_use_Q_wetdry_adjustment
    9698
    9799
    … …  
    119121        self.accumulated_flow = 0.0
    120122        self.discharge = 0.0
     123        self.discharge_function_value = 0.0
    121124        self.velocity = 0.0
    122125        self.outlet_depth = 0.0
    … …  
    218221        old_inflow_ymom = self.inflow.get_average_ymom()
    219222
    220         #semi_implicit = True
    221         #if semi_implicit:   
    222 
    223         # FIXME: This update replaces Q with Q*new_inflow_depth/old_inflow_depth
    224         # This has good wetting and drying properties but means that
    225         # discharge != Q.
    226         # We should be able to check if old_inflow_depth*old_inflow_area > Q*dt,
    227         # and in that case we don't need this implicit trick, and could
    228         # have the discharge = Q (whereas in the nearly-dry case we want
    229         # the trick).
    230 
    231223        # Implement the update of flow over a timestep by
    232224        # using a semi-implict update. This ensures that
    … …  
    237229            dt_Q_on_d = 0.0
    238230
    239         # The depth update is:
     231        always_use_Q_wetdry_adjustment = self.always_use_Q_wetdry_adjustment
     232        use_Q_wetdry_adjustment = ((always_use_Q_wetdry_adjustment) |\
     233            (old_inflow_depth*self.inflow.get_area() <= Q*timestep))
     234        # If we use the 'Q_adjustment', then the discharge is rescaled so that
     235        # the depth update is:
    240236        #    new_inflow_depth*inflow_area =
    241237        #    old_inflow_depth*inflow_area -
    242238        #    timestep*Q*(new_inflow_depth/old_inflow_depth)
    243         # The last term in () is a wet-dry improvement trick
    244         # Note inflow_area is the area of all triangles in the inflow
    245         # region -- so this is a volumetric balance equation
     239        # The last term in () is a wet-dry improvement trick (which rescales Q)
     240        #
     241        # Before Feb 2015 this rescaling was always done (even if the flow was
     242        # not near wet-dry). Now if use_old_Q_adjustment=False, the rescaling
     243        # is only done if required to avoid drying
     244        #
    246245        #
    247246        factor = 1.0/(1.0 + dt_Q_on_d/self.inflow.get_area())
    248         new_inflow_depth = old_inflow_depth*factor
     247
     248
     249        if use_Q_wetdry_adjustment:
     250            # FIXME:
     251            new_inflow_depth = old_inflow_depth*factor
     252
     253            if(old_inflow_depth>0.):
     254                timestep_star = timestep*new_inflow_depth/old_inflow_depth
     255            else:
     256                timestep_star = 0.
     257
     258        else:
     259            new_inflow_depth = old_inflow_depth - timestep*Q/self.inflow.get_area()
     260            timestep_star = timestep
    249261
    250262        if(self.use_old_momentum_method):
    … …  
    256268        else:
    257269            # For the momentum balance, note that Q also advects the momentum,
    258             # The volumetric momentum flux should be Q*momentum, where
     270            # The volumetric momentum flux should be Q*momentum/depth, where
    259271            # momentum has an average value of new_inflow_mom (or
    260             # old_inflow_mom). For consistency we keep using the
    261             # (new_inflow_depth/old_inflow_depth) factor for discharge:
     272            # old_inflow_mom).  We use old_inflow_depth for depth
    262273            #
    263274            #     new_inflow_xmom*inflow_area =
    264275            #     old_inflow_xmom*inflow_area -
    265             #     [timestep*Q*(new_inflow_depth/old_inflow_depth)]*new_inflow_xmom
     276            #     [timestep*Q]*new_inflow_xmom/old_inflow_depth
    266277            # and:
    267278            #     new_inflow_ymom*inflow_area =
    268279            #     old_inflow_ymom*inflow_area -
    269             #     [timestep*Q*(new_inflow_depth/old_inflow_depth)]*new_inflow_ymom
     280            #     [timestep*Q]*new_inflow_ymom/old_inflow_depth
    270281            #
    271             # The choice of new_inflow_mom in the final term at the end might be
    272             # replaced with old_inflow_mom
     282            # The units balance: m^2/s*m^2 = m^2/s*m^2 - s*m^3/s*m^2/s *m^(-1)
    273283            #
    274             factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/self.inflow.get_area())
     284            if old_inflow_depth > 0.:
     285                if use_Q_wetdry_adjustment:
     286                    # Replace dt*Q with dt*Q*new_inflow_depth/old_inflow_depth = dt_Q_on_d*new_inflow_depth
     287                    factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/(old_inflow_depth*self.inflow.get_area()))
     288                else:
     289                    factor2 = 1.0/(1.0 + timestep*Q/(old_inflow_depth*self.inflow.get_area()))
     290            else:
     291                factor2 = 0.
     292
    275293            new_inflow_xmom = old_inflow_xmom*factor2
    276294            new_inflow_ymom = old_inflow_ymom*factor2
    … …  
    289307
    290308        # set outflow
    291         if old_inflow_depth > 0.0 :
    292             timestep_star = timestep*new_inflow_depth/old_inflow_depth
    293         else:
    294             timestep_star = 0.0
    295 
    296309        outflow_extra_depth = Q*timestep_star/self.outflow.get_area()
    297310        outflow_direction = - self.outflow.outward_culvert_vector
    … …  
    306319        self.accumulated_flow += gain
    307320        self.discharge  = Q*timestep_star/timestep
     321        self.discharge_function_value = Q
    308322        self.velocity =   barrel_speed
    309323        self.outlet_depth = outlet_depth
    … …  
    541555       
    542556        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
     557        message += 'Discharge_function_value [m^3/s]: %.2f\n' % self.discharge_function_value
    543558        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
    544559        message += 'Outlet Depth  [m]: %.2f\n' % self.outlet_depth
    … …  
    562577            self.log_filename = self.domain.get_datadir() + '/' + self.label + '.log'
    563578            log_to_file(self.log_filename, self.statistics(), mode='w')
    564             log_to_file(self.log_filename, 'time, discharge, velocity, accumulated_flow, driving_energy, delta_total_energy')
     579            log_to_file(self.log_filename, 'time, discharge, discharge_function_value, velocity, accumulated_flow, driving_energy, delta_total_energy')
    565580
    566581            #log_to_file(self.log_filename, self.culvert_type)
    … …  
    571586        message  = '%.5f, ' % self.domain.get_time()
    572587        message += '%.5f, ' % self.discharge
     588        message += '%.5f, ' % self.discharge_function_value
    573589        message += '%.5f, ' % self.velocity
    574590        message += '%.5f, ' % self.accumulated_flow
Note: See TracChangeset for help on using the changeset viewer.