Changeset 9671


Ignore:
Timestamp:
Feb 12, 2015, 4:54:53 PM (9 years ago)
Author:
davies
Message:

Implementing implicit discharge evaluation for internal boundary operator

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/anuga/parallel/parallel_internal_boundary_operator.py

    r9669 r9671  
    22import math
    33import numpy
     4from numpy.linalg import solve
    45
    56#from anuga.structures.boyd_box_operator import boyd_box_function
     
    109110        self.smoothing_timescale = smoothing_timescale
    110111
     112    #def __call__(self):
     113    #    """
     114    #        Update for n sub-timesteps
     115    #    """
     116
     117    #    number_of_substeps = 20
     118    #    original_timestep = self.domain.get_timestep()
     119    #    self.domain.timestep = original_timestep/(1.0*number_of_substeps)
     120    #   
     121    #    for i in range(number_of_substeps):
     122    #        anuga.parallel.parallel_structure_operator.Parallel_Structure_operator.__call__(self)
     123
     124    #    self.domain.timestep = original_timestep
     125
    111126    def parallel_safe(self):
    112127
     
    115130
    116131
    117     def discharge_routine(self):
     132    def discharge_routine_old(self):
    118133
    119134        import pypar
     
    252267       
    253268       
     269    def discharge_routine(self):
     270        """
     271            Uses semi-implicit discharge estimation:
     272              Discharge = 0.5*(Q(H0, T0) + Q(H0 + delta_H, T0+delta_T))
     273            where H0 = headwater stage, T0 = tailwater stage, delta_H = change in
     274            headwater stage over a timestep, delta_T = change in tailwater stage over a
     275            timestep.
     276
     277            We can estimate delta_H, delta_T by solving the following system (based on mass conservation):
     278              A0*delta_H = -dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
     279              A1*delta_T =  dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
     280            where A0, A1 are the inlet areas, and dt is the timestep
     281
     282            We linearise the system with the approximation:
     283              Q(H0 + delta_H, T0 + delta_T) ~= Q(H0, T0) + delQ/delH * delta_H + delQ/delT*delta_T
     284            where delQ/delH and delQ/delT are evaluated with numerical finite differences
     285
     286
     287        """
     288
     289        import pypar
     290
     291        local_debug = False
     292       
     293        # If the structure has been closed, then no water gets through
     294        if self.height <= 0.0:
     295            if self.myid == self.master_proc:
     296                Q = 0.0
     297                barrel_velocity = 0.0
     298                outlet_culvert_depth = 0.0
     299                self.case = "Structure is blocked"
     300                self.inflow = self.inlets[0]
     301                self.outflow = self.inlets[1]
     302                return Q, barrel_velocity, outlet_culvert_depth
     303            else:
     304                return None, None, None
     305
     306        #Send attributes of both enquiry points to the master proc
     307        if self.myid == self.master_proc:
     308
     309            if self.myid == self.enquiry_proc[0]:
     310                enq_total_energy0 = self.inlets[0].get_enquiry_total_energy()
     311                enq_stage0 = self.inlets[0].get_enquiry_stage()
     312            else:
     313                enq_total_energy0 = pypar.receive(self.enquiry_proc[0])
     314                enq_stage0 = pypar.receive(self.enquiry_proc[0])
     315
     316
     317            if self.myid == self.enquiry_proc[1]:
     318                enq_total_energy1 = self.inlets[1].get_enquiry_total_energy()
     319                enq_stage1 = self.inlets[1].get_enquiry_stage()
     320            else:
     321                enq_total_energy1 = pypar.receive(self.enquiry_proc[1])
     322                enq_stage1 = pypar.receive(self.enquiry_proc[1])
     323
     324        else:
     325            if self.myid == self.enquiry_proc[0]:
     326                pypar.send(self.inlets[0].get_enquiry_total_energy(), self.master_proc)
     327                pypar.send(self.inlets[0].get_enquiry_stage(), self.master_proc)
     328
     329            if self.myid == self.enquiry_proc[1]:
     330                pypar.send(self.inlets[1].get_enquiry_total_energy(), self.master_proc)
     331                pypar.send(self.inlets[1].get_enquiry_stage(), self.master_proc)
     332
     333        # Send inlet areas to the master proc. FIXME: Inlet areas don't change
     334        # -- perhaps we could just do this once?
     335
     336        # area0
     337        if self.myid in self.inlet_procs[0]:
     338            area0 = self.inlets[0].get_global_area()
     339
     340        if self.myid == self.master_proc:
     341            if self.myid != self.inlet_master_proc[0]:
     342                area0 = pypar.receive(self.inlet_master_proc[0])
     343        elif self.myid == self.inlet_master_proc[0]:
     344            pypar.send(area0, self.master_proc)
     345       
     346        # area1
     347        if self.myid in self.inlet_procs[1]:
     348            area1 = self.inlets[1].get_global_area()
     349
     350        if self.myid == self.master_proc:
     351            if self.myid != self.inlet_master_proc[1]:
     352                area1 = pypar.receive(self.inlet_master_proc[1])
     353        elif self.myid == self.inlet_master_proc[1]:
     354            pypar.send(area1, self.master_proc)
     355
     356        # Compute discharge
     357        if self.myid == self.master_proc:
     358
     359            # Energy or stage as head
     360            if self.use_velocity_head:
     361                E0 = enq_total_energy0
     362                E1 = enq_total_energy1
     363            else:
     364                E0 = enq_stage0
     365                E1 = enq_stage1
     366
     367            # Variables for anuga's structure operator
     368            self.delta_total_energy = E0 - E1
     369            self.driving_energy = max(E0, E1)
     370
     371
     372            Q0 = self.internal_boundary_function(E0, E1)
     373            dt = self.domain.get_timestep()
     374            if dt > 0.:
     375                # Numerical derivatives of Discharge function
     376                dE0 = 0.01
     377                dE1 = 0.01
     378                dQ_dE0 = (self.internal_boundary_function(E0+dE0, E1) - self.internal_boundary_function(E0-dE0, E1))/(2*dE0)
     379                dQ_dE1 = (self.internal_boundary_function(E0, E1+dE1) - self.internal_boundary_function(E0, E1-dE1))/(2*dE1)
     380
     381                # Assemble and solve matrix
     382                hdt = 0.5*dt
     383
     384                M11 = area0 + dQ_dE0*hdt
     385                M12 = dQ_dE1*hdt
     386                M21 = -dQ_dE0*hdt
     387                M22 = area1 - dQ_dE1*hdt
     388
     389                lhs = numpy.array([ [M11, M12], [M21, M22]])
     390                rhs = numpy.array([ -Q0*dt, Q0*dt])
     391                # sol contains delta_E0, delta_E1
     392                sol = solve(lhs, rhs)
     393               
     394                Q = 0.5*(Q0 + ( Q0 + sol[0]*dQ_dE0 + sol[1]*dQ_dE1))
     395            else:
     396                Q = Q0
     397
     398            # Smooth discharge
     399            if dt > 0.:
     400                ts = dt/max(dt, self.smoothing_timescale, 1.0e-30)
     401            else:
     402                # No smoothing
     403                ts = 1.0
     404
     405            self.smooth_Q = self.smooth_Q + ts*(Q - self.smooth_Q)
     406
     407        else:
     408            self.delta_total_energy=numpy.nan
     409            self.driving_energy=numpy.nan
     410
     411
     412        self.inflow_index = 0
     413        self.outflow_index = 1
     414        # master proc orders reversal if applicable
     415        if self.myid == self.master_proc:
     416
     417            # Reverse the inflow and outflow direction?
     418            if Q < 0.:
     419                self.inflow_index = 1
     420                self.outflow_index = 0
     421
     422                for i in self.procs:
     423                    if i == self.master_proc: continue
     424                    pypar.send(True, i)
     425            else:
     426                for i in self.procs:
     427                    if i == self.master_proc: continue
     428                    pypar.send(False, i)
     429
     430        else:
     431            reverse = pypar.receive(self.master_proc)
     432
     433            if reverse:
     434                self.inflow_index = 1
     435                self.outflow_index = 0
     436
     437        # Master proc computes return values
     438        if self.myid == self.master_proc:
     439            # Zero Q if sign's of smooth_Q and Q differ
     440            if numpy.sign(self.smooth_Q) != numpy.sign(Q):
     441                Q = 0.
     442            else:
     443                # Make Q positive (for anuga's structure operator)
     444                Q = min( abs(self.smooth_Q), abs(Q) )
     445            # Variables required by anuga's structure operator which are
     446            # not used
     447            barrel_velocity = numpy.nan
     448            outlet_culvert_depth = numpy.nan
     449            return Q, barrel_velocity, outlet_culvert_depth
     450        else:
     451            return None, None, None
  • trunk/anuga_core/anuga/parallel/parallel_structure_operator.py

    r9658 r9671  
    193193                               enquiry_proc = self.enquiry_proc[0],
    194194                               verbose = self.verbose))
    195        
     195
    196196            if force_constant_inlet_elevations:
    197197                # Try to enforce a constant inlet elevation
    198198                inlet_global_elevation = self.inlets[-1].get_global_average_elevation()
    199199                self.inlets[-1].set_elevations(inlet_global_elevation)
    200             
     200       
    201201        else:
    202202            self.inlets.append(None)
     
    230230                               enquiry_proc = self.enquiry_proc[1],
    231231                               verbose = self.verbose))
     232
     233            if force_constant_inlet_elevations:
     234                # Try to enforce a constant inlet elevation
     235                inlet_global_elevation = self.inlets[-1].get_global_average_elevation()
     236                self.inlets[-1].set_elevations(inlet_global_elevation)
     237           
    232238
    233239        else:
  • trunk/anuga_core/anuga/structures/internal_boundary_operator.py

    r9669 r9671  
    22import math
    33import numpy
     4from numpy.linalg import solve
    45
    56
     
    108109        # Finally, set the smoothing timescale we actually want
    109110        self.smoothing_timescale = smoothing_timescale
     111
     112    #def __call__(self):
     113    #    """
     114    #        Update for n sub-timesteps
     115    #    """
     116
     117    #    number_of_substeps = 1
     118    #    original_timestep = self.domain.get_timestep()
     119    #    self.domain.timestep = original_timestep/(1.0*number_of_substeps)
     120    #   
     121    #    for i in range(number_of_substeps):
     122    #        anuga.Structure_operator.__call__(self)
     123
     124    #    self.domain.timestep = original_timestep
    110125   
    111     ###########################################################################
    112     def discharge_routine(self):
     126    def discharge_routine_old(self):
    113127        """Procedure to determine the inflow and outflow inlets.
    114128        Then use self.internal_boundary_function to do the actual calculation
     
    191205
    192206        return Q, barrel_velocity, outlet_culvert_depth
     207
     208
     209    def discharge_routine(self):
     210        """Uses semi-implicit discharge estimation:
     211          Discharge = 0.5*(Q(H0, T0) + Q(H0 + delta_H, T0+delta_T))
     212        where H0 = headwater stage, T0 = tailwater stage, delta_H = change in
     213        headwater stage over a timestep, delta_T = change in tailwater stage over a
     214        timestep, and Q is the discharge function
     215
     216        We can estimate delta_H, delta_T by solving the following system (based on mass conservation):
     217          A0*delta_H = -dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
     218          A1*delta_T =  dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
     219        where A0, A1 are the inlet areas, and dt is the timestep
     220
     221        We linearise the system with the approximation:
     222          Q(H0 + delta_H, T0 + delta_T) ~= Q(H0, T0) + delQ/delH * delta_H + delQ/delT*delta_T
     223        where delQ/delH and delQ/delT are evaluated with numerical finite differences
     224
     225
     226        """
     227
     228        # Compute energy head or stage at inlets 0 and 1
     229        if self.use_velocity_head:
     230            self.inlet0_energy = self.inlets[0].get_enquiry_total_energy()
     231            self.inlet1_energy = self.inlets[1].get_enquiry_total_energy()
     232        else:
     233            self.inlet0_energy = self.inlets[0].get_enquiry_stage()
     234            self.inlet1_energy = self.inlets[1].get_enquiry_stage()
     235       
     236        # Store these variables for anuga's structure output
     237        self.driving_energy = max(self.inlet0_energy, self.inlet1_energy)
     238        self.delta_total_energy = self.inlet0_energy - self.inlet1_energy
     239
     240        Q0 = self.internal_boundary_function(self.inlet0_energy, self.inlet1_energy)
     241        dt = self.domain.get_timestep()
     242
     243
     244        if dt > 0.:
     245            E0 = self.inlet0_energy
     246            E1 = self.inlet1_energy
     247
     248            # Numerical derivatives
     249            dE0 = 0.01
     250            dE1 = 0.01
     251            dQ_dE0 = (self.internal_boundary_function(E0+dE0, E1) - self.internal_boundary_function(E0-dE0, E1))/(2*dE0)
     252            dQ_dE1 = (self.internal_boundary_function(E0, E1+dE1) - self.internal_boundary_function(E0, E1-dE1))/(2*dE1)
     253
     254            # Assemble and solve matrix
     255            A0 = self.inlets[0].area
     256            A1 = self.inlets[1].area
     257            hdt = 0.5*dt
     258
     259            M11 = A0 + dQ_dE0*hdt
     260            M12 = dQ_dE1*hdt
     261            M21 = -dQ_dE0*hdt
     262            M22 = A1 - dQ_dE1*hdt
     263
     264            lhs = numpy.array([ [M11, M12], [M21, M22]])
     265            rhs = numpy.array([ -Q0*dt, Q0*dt])
     266            # sol contains delta_E0, delta_E1
     267            sol = solve(lhs, rhs)
     268           
     269            Q = 0.5*(Q0 + ( Q0 + sol[0]*dQ_dE0 + sol[1]*dQ_dE1))
     270        else:
     271            Q = Q0
     272
     273
     274        # Use time-smoothed discharge if smoothing_timescale > 0.
     275        if dt > 0.0:
     276            ts = dt/max(dt, self.smoothing_timescale, 1.0e-30)
     277        else:
     278            # No smoothing
     279            ts = 1.0
     280
     281        self.smooth_Q = self.smooth_Q + ts*(Q - self.smooth_Q)
     282
     283        # Define 'inflow' and 'outflow' for anuga's structure operator
     284        if Q >= 0.:
     285            self.inflow  = self.inlets[0]
     286            self.outflow = self.inlets[1]
     287        else:
     288            self.inflow  = self.inlets[1]
     289            self.outflow = self.inlets[0]
     290
     291        # Zero discharge if the sign's of Q and smooth_Q are not the same
     292        if numpy.sign(self.smooth_Q) != numpy.sign(Q):
     293            Q = 0.
     294        else:
     295            # Make Q positive (for anuga's structure operator)
     296            Q = min( abs(self.smooth_Q), abs(Q) )
     297
     298        barrel_velocity = numpy.nan
     299        outlet_culvert_depth = numpy.nan
     300
     301        return Q, barrel_velocity, outlet_culvert_depth
  • trunk/anuga_core/anuga/structures/structure_operator.py

    r9658 r9671  
    130130        else:
    131131            raise Exception, 'Define either exchange_lines or end_points'
    132 
    133        
     132       
     133
    134134        self.inlets = []
    135135        line0 = self.exchange_lines[0] #self.inlet_lines[0]
     
    162162            inlet_global_elevation = self.inlets[-1].get_average_elevation()
    163163            self.inlets[-1].set_elevations(inlet_global_elevation)
    164        
    165164
    166165        tris_0 = self.inlets[0].triangle_indices
     
    200199        self.set_logging(logging)
    201200
     201        if force_constant_inlet_elevations:
     202            # Try to enforce a constant inlet elevation
     203            inlet_global_elevation = self.inlets[-1].get_average_elevation()
     204            self.inlets[-1].set_elevations(inlet_global_elevation)
     205       
    202206
    203207
  • trunk/anuga_work/development/gareth/tests/ras_bridge/channel_floodplain1.py

    r9634 r9671  
    188188    exchange_lines=[bridge_in, bridge_out],
    189189    enquiry_gap=0.01,
    190     use_velocity_head=False,
    191     smoothing_timescale=30.0,
     190    zero_outflow_momentum=False,
     191    smoothing_timescale=0.0,
    192192    logging=True)
    193193   
Note: See TracChangeset for help on using the changeset viewer.