Changeset 5741


Ignore:
Timestamp:
Sep 5, 2008, 9:16:45 PM (14 years ago)
Author:
steve
Message:

Updated timestepping

Location:
anuga_work/development
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/config.py

    r5740 r5741  
    1111
    1212min_timestep = 1.0e-6 #Should be computed based on geometry
    13 max_timestep = 1000
     13max_timestep = 1.0e3
     14max_smallsteps = 50  # Max number of degenerate steps allowed b4 trying first order
    1415#This is how:
    1516#Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)
  • anuga_work/development/anuga_1d/dam_h_elevation.py

    r5738 r5741  
    11import os
    22from math import sqrt, sin, cos, pi, exp
    3 from shallow_water_domain import *
     3from shallow_water_domain_new import *
    44from Numeric import zeros, Float
    55#from analytic_dam_sudi import AnalyticDam
     
    4242
    4343domain.default_order = 2
    44 domain.default_time_order = 1
     44domain.set_timestepping_method('euler')
     45#domain.default_time_order = 2
    4546domain.cfl = 1.0
    4647domain.limiter = "vanleer"
  • anuga_work/development/anuga_1d/domain.py

    r5740 r5741  
    724724    #--------------------------   
    725725
    726     def evolve(self, yieldstep = None, finaltime = None, duration = None,
     726    def evolve(self, yieldstep = None,
     727               finaltime = None,
     728               duration = None,
    727729               skip_initial_step = False):
    728730        """Evolve model through time starting from self.starttime.
     
    759761
    760762        from config import min_timestep, max_timestep, epsilon
    761 
    762         print 'yieldstep', yieldstep
    763         print 'duration', duration
    764763
    765764        # FIXME: Maybe lump into a larger check prior to evolving
     
    13881387    def update_timestep(self, yieldstep, finaltime):
    13891388
    1390         from config import min_timestep
     1389        from config import min_timestep, max_timestep
    13911390
    13921391        # self.timestep is calculated from speed of characteristics
    13931392        # Apply CFL condition here
    1394         timestep = self.CFL*self.timestep
     1393        timestep = min(self.CFL*self.flux_timestep, max_timestep)
    13951394
    13961395        #Record maximal and minimal values of timestep for reporting
     
    14591458   
    14601459    #def update_conserved_quantities(self):
    1461     def update_conserved_quantities(self,timestep):
     1460    def update_conserved_quantities(self):
    14621461        """Update vectors of conserved quantities using previously
    14631462        computed fluxes specified forcing functions.
     
    14691468        d = len(self.conserved_quantities)
    14701469
    1471         #timestep = self.timestep
     1470        timestep = self.timestep
    14721471
    14731472        #Compute forcing terms
    1474         #self.compute_forcing_terms()
     1473        self.compute_forcing_terms()
    14751474
    14761475        #Update conserved_quantities
     
    14791478            Q.update(timestep)
    14801479
    1481             #Clean up
    1482             #Note that Q.explicit_update is reset by compute_fluxes
    1483 
    1484             #MH090605 commented out the following since semi_implicit_update is now re-initialized
    1485             #at the end of the _update function in quantity_ext.c (This is called by the
    1486             #preceeding Q.update(timestep) statement above).
    1487             #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs
    1488             #to 8.35 secs
    1489 
    1490             #Q.semi_implicit_update[:] = 0.0
     1480
    14911481
    14921482if __name__ == "__main__":
  • anuga_work/development/anuga_1d/quantity.py

    r5738 r5741  
    861861        #backup_centroid_values(self)
    862862        N = self.domain.number_of_elements
    863         for i in range(self.N):
    864             quantity.centroid_backup_values[i] = quantity.centroid_values[i]
     863        for i in range(N):
     864            self.centroid_backup_values[i] = self.centroid_values[i]
    865865
    866866    def saxpy_centroid_values(self,a,b):
     
    869869        #saxpy_centroid_values(self,a,b)
    870870        N = self.domain.number_of_elements
    871         for i in range(self.N):
    872             quantity.centroid_backup_values[i] = a*quantity.centroid_values[i] + b*quantity.centroid_backup_values[i]       
     871        for i in range(N):
     872            self.centroid_backup_values[i] = a*self.centroid_values[i] + b*self.centroid_backup_values[i]       
    873873
    874874class Conserved_quantity(Quantity):
  • anuga_work/development/anuga_1d/shallow_water_domain.py

    r5740 r5741  
    575575        #print "flux cell",k,flux[0]
    576576
    577     domain.timestep = timestep
     577    domain.flux_timestep = timestep
    578578    #print domain.quantities['stage'].centroid_values
    579579
     
    652652        from comp_flux_ext import compute_fluxes_ext
    653653       
    654         domain.timestep = compute_fluxes_ext(timestep,
     654        domain.flux_timestep = compute_fluxes_ext(timestep,
    655655                                  epsilon,
    656656                                  g,
     
    685685        from comp_flux_ext import compute_fluxes_ext_short
    686686
    687         domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
     687        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
    688688
    689689
  • anuga_work/development/anuga_1d/shallow_water_domain_new.py

    r5737 r5741  
    165165        distribute_to_vertices_and_edges(self)
    166166       
    167     def evolve(self, yieldstep = None, finaltime = None,
     167    def evolve(self, yieldstep = None,
     168               finaltime = None,
     169               duration = None,
    168170               skip_initial_step = False):
    169171        """Specialisation of basic evolve method from parent class
     
    172174        self.distribute_to_vertices_and_edges()
    173175       
    174         for t in Generic_Domain.evolve(self, yieldstep, finaltime,
     176        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration,
    175177                                       skip_initial_step):
    176178            #Pass control on to outer loop for more specific actions
     
    439441       
    440442
    441     domain.timestep = timestep
     443    domain.flux_timestep = timestep
    442444    #print domain.quantities['stage'].centroid_values
    443445
     
    468470        from comp_flux_ext import compute_fluxes_ext
    469471       
    470         domain.timestep = compute_fluxes_ext(timestep,
     472        domain.flux_timestep = compute_fluxes_ext(timestep,
    471473                                  epsilon,
    472474                                  g,
     
    501503        from comp_flux_ext import compute_fluxes_ext_short
    502504
    503         domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
     505        domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
    504506
    505507
     
    526528        max_speed_array = domain.max_speed_array
    527529               
    528         from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced                  #from comp_flux_ext import compute_fluxes_ext
     530        from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
    529531       
    530         domain.timestep = compute_fluxes_ext_wellbalanced(timestep,
     532        domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
    531533                                  epsilon,
    532534                                  g,
  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5738 r5741  
    55
    66
    7 from shallow_water_domain import *
    8 from shallow_water_domain import flux_function as domain_flux_function
     7from shallow_water_domain_new import *
     8from shallow_water_domain_new import flux_function as domain_flux_function
    99
    1010from Numeric import allclose, array, ones, Float, maximum, zeros
     
    4141        #print domain.quantities['stage'].explicit_update
    4242        #print domain.quantities['xmomentum'].explicit_update
    43         #print stage_ud
     43        print stage_ud
    4444
    4545        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
     
    133133
    134134        domain.default_order = 1
    135         domain.default_time_order = 1
     135        domain.set_timestepping_method('euler')
     136       
    136137        yieldstep=1.0
    137138        finaltime=1.0
     
    167168
    168169        domain.default_order = 2
    169         domain.default_time_order = 1
     170        domain.set_timestepping_method('rk2')
    170171        yieldstep=1.0
    171172        finaltime=1.0
     
    191192
    192193        domain.default_order = 2
    193         domain.default_time_order = 2
     194        domain.set_timestepping_method('rk3')
     195
    194196        yieldstep=1.0
    195197        finaltime=1.0
  • anuga_work/development/demos/netherlands.py

    r5153 r5741  
    2727N = 130 # size = 33800
    2828N = 600 # Size = 720000
    29 N = 100
     29N = 50
    3030
    3131points, elements, boundary = rectangular_cross(N, N)
     
    3434domain.check_integrity()
    3535
     36import sys
     37base = os.path.basename(sys.argv[0])
     38domain.simulation_name = 'netherlands'
    3639domain.set_name(os.path.splitext(__file__)[0])
    37 domain.set_timestepping_method('rk3')
     40domain.set_timestepping_method('rk2')
    3841domain.set_default_order(2)
    3942domain.set_store_vertices_uniquely(True) # Store as internally represented
     
    130133domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
    131134
     135#-------------------------------------------------------------------------------
     136# Copy scripts to time stamped output directory and capture screen
     137# output to file
     138#-------------------------------------------------------------------------------
     139time = strftime('%Y%m%d_%H%M%S',localtime())
     140
     141output_dir = 'dam_break_'+time
     142output_file = 'dam_break'
     143
     144copy_code_files(output_dir,__file__)
     145#start_screen_catcher(output_dir+'_')
     146
    132147
    133148#------------------------------------------------------------------------------
     
    139154    # Initialise real-time visualiser
    140155
    141     pass
    142     #vis = RealtimeVisualiser(domain)
    143     #vis.render_quantity_height("elevation", dynamic=False)
    144     #vis.render_quantity_height("stage", dynamic=True)
    145     #vis.colour_height_quantity('stage', (0.0, 0.0, 0.8))
    146     #vis.start()
     156    #pass
     157    visualise = True
     158    if visualise:
     159        from anuga.visualiser import RealtimeVisualiser
     160        vis = RealtimeVisualiser(domain)
     161        vis.render_quantity_height("elevation", offset=0.01, dynamic=False)
     162        vis.render_quantity_height("stage", dynamic=True)
     163        vis.colour_height_quantity('stage', (0.2, 0.2, 0.8))
     164        vis.start()
     165        import time
     166        time.sleep(2.0)
     167   
     168   
    147169
    148170
     
    151173t0 = time.time()
    152174
    153 for t in domain.evolve(yieldstep = 0.005, finaltime = None):
     175for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
    154176    print domain.timestepping_statistics()
    155     print domain.quantities['stage'].get_values(location='centroids',
    156                                                 indices=[0])
    157                                                
    158     #vis.update()                                               
    159     #time.sleep(0.1)
    160     #raw_input('pause>')
    161     #V.update_quantity('stage')
    162     #rpdb.set_active()
    163     #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
    164 #vis.evolveFinished()
     177    #print domain.quantities['stage'].get_values(location='centroids',
     178                                                #indices=[0])
     179    if visualise:
     180        vis.update()                                           
     181
     182       
     183
     184if visualise: vis.evolveFinished()
    165185
    166186print 'That took %.2f seconds' %(time.time()-t0)
  • anuga_work/development/shallow_water_1d_old/quantity.py

    r4959 r5741  
    5656        #Intialise centroid values
    5757        self.interpolate()
     58
     59
     60        from Numeric import zeros, Float
     61       
     62        #Allocate space for boundary values
     63        #L = len(domain.boundary)
     64        self.boundary_values = zeros(2, Float) #assumes no parrellism
     65
     66        #Allocate space for updates of conserved quantities by
     67        #flux calculations and forcing functions
     68       
     69        N = domain.number_of_elements
     70        self.explicit_update = zeros(N, Float )
     71        self.semi_implicit_update = zeros(N, Float )
     72
     73        self.gradients = zeros(N, Float)
     74        self.qmax = zeros(self.centroid_values.shape, Float)
     75        self.qmin = zeros(self.centroid_values.shape, Float)
     76
     77        self.beta = domain.beta       
    5878
    5979    #Methods for operator overloading
     
    381401        return integral
    382402
    383 class Conserved_quantity(Quantity):
    384     """Class conserved quantity adds to Quantity:
    385 
    386     storage and method for updating, and
    387     methods for extrapolation from centropid to vertices inluding
    388     gradients and limiters
    389     """
    390 
    391     def __init__(self, domain, vertex_values=None):
    392         Quantity.__init__(self, domain, vertex_values)
    393 
    394         from Numeric import zeros, Float
    395        
    396         #Allocate space for boundary values
    397         #L = len(domain.boundary)
    398         self.boundary_values = zeros(2, Float) #assumes no parrellism
    399 
    400         #Allocate space for updates of conserved quantities by
    401         #flux calculations and forcing functions
    402        
    403         N = domain.number_of_elements
    404         self.explicit_update = zeros(N, Float )
    405         self.semi_implicit_update = zeros(N, Float )
    406 
    407         self.gradients = zeros(N, Float)
    408         self.qmax = zeros(self.centroid_values.shape, Float)
    409         self.qmin = zeros(self.centroid_values.shape, Float)
    410 
    411         self.beta = domain.beta
    412403
    413404    def update(self, timestep):
     
    863854                else:
    864855                    qv[k,i] = qc[k]
     856   
     857
     858class Conserved_quantity(Quantity):
     859    """Class conserved quantity adds to Quantity:
     860
     861    storage and method for updating, and
     862    methods for extrapolation from centropid to vertices inluding
     863    gradients and limiters
     864    """
     865
     866    def __init__(self, domain, vertex_values=None):
     867        Quantity.__init__(self, domain, vertex_values)
     868
     869        print "Use Quantity instead of Conserved_quantity"
     870
    865871                   
    866872
Note: See TracChangeset for help on using the changeset viewer.