Changeset 5738


Ignore:
Timestamp:
Sep 5, 2008, 3:06:01 PM (14 years ago)
Author:
steve
Message:

Copied evolve code from anuga_2d

Location:
anuga_work/development/anuga_1d
Files:
5 edited

Legend:

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

    r5728 r5738  
    66
    77anuga_dir = '..' + os.sep + '..' + os.sep + '..'
    8 #os.chdir(anuga_dir)
    9 print os.getcwd() 
    10 
    118utilities_dir = anuga_dir + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities'
    129
    13 print anuga_dir
    14 print utilities_dir
    1510execfile( utilities_dir + os.sep + 'compile.py')
    1611
    1712
     13os.chdir(buildroot)   
    1814
    19 os.chdir(buildroot)   
    20 #execfile('test_all.py')
    2115   
  • anuga_work/development/anuga_1d/dam_h_elevation.py

    r5737 r5738  
    4242
    4343domain.default_order = 2
    44 domain.default_time_order = 2
     44domain.default_time_order = 1
    4545domain.cfl = 1.0
    4646domain.limiter = "vanleer"
     
    6161    for i in range(len(x)):
    6262        if 1100.0 <= x[i] <=1200.0:
    63             y[i]=10.001
     63            y[i]=10.0
    6464        else:
    6565            y[i]=10.0
     
    140140        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
    141141        #print 'Test1'
     142        ion()
    142143        hold(False)
    143144        #clf()
     
    150151
    151152        plot(X,ElevationQ,X,HeightQ)
    152         plot1.set_ylim([-1,12])
     153        plot1.set_ylim([9.99,10.01])
    153154        xlabel('Position')
    154155        ylabel('Stage')
     
    161162        plot2 = subplot(212)
    162163        plot(X,MomentumQ)
    163         plot2.set_ylim([-1,1])
     164        #plot2.set_ylim([-1,1])
    164165       
    165166        #legend( ('Numerical Solution', 'for momentum'),
     
    173174        file += ".eps"
    174175        #savefig(file)
    175         show()
     176        #show()
    176177       
    177178print 'That took %.2f seconds'%(time.time()-t0)
     179
     180raw_input("Press any key")
  • anuga_work/development/anuga_1d/domain.py

    r5737 r5738  
    697697        self.datadir = name
    698698
    699 #Main components of evolve
    700     def evolve(self, yieldstep = None, finaltime = None,
     699
     700    #--------------------------
     701    # Main components of evolve
     702    #--------------------------   
     703
     704    def evolve(self,
     705               yieldstep = None,
     706               finaltime = None,
     707               duration = None,
     708               skip_initial_step = False):
     709        """Evolve model through time starting from self.starttime.
     710
     711
     712        yieldstep: Interval between yields where results are stored,
     713                   statistics written and domain inspected or
     714                   possibly modified. If omitted the internal predefined
     715                   max timestep is used.
     716                   Internally, smaller timesteps may be taken.
     717
     718        duration: Duration of simulation
     719
     720        finaltime: Time where simulation should end. This is currently
     721        relative time.  So it's the same as duration.
     722
     723        If both duration and finaltime are given an exception is thrown.
     724
     725
     726        skip_initial_step: Boolean flag that decides whether the first
     727        yield step is skipped or not. This is useful for example to avoid
     728        duplicate steps when multiple evolve processes are dove tailed.
     729
     730
     731        Evolve is implemented as a generator and is to be called as such, e.g.
     732
     733        for t in domain.evolve(yieldstep, finaltime):
     734            <Do something with domain and t>
     735
     736
     737        All times are given in seconds
     738
     739        """
     740
     741        from config import min_timestep, max_timestep, epsilon
     742
     743        # FIXME: Maybe lump into a larger check prior to evolving
     744        msg = 'Boundary tags must be bound to boundary objects before '
     745        msg += 'evolving system, '
     746        msg += 'e.g. using the method set_boundary.\n'
     747        msg += 'This system has the boundary tags %s '\
     748               %self.get_boundary_tags()
     749        assert hasattr(self, 'boundary_objects'), msg
     750
     751
     752        if yieldstep is None:
     753            yieldstep = max_timestep
     754        else:
     755            yieldstep = float(yieldstep)
     756
     757        self._order_ = self.default_order
     758
     759
     760        if finaltime is not None and duration is not None:
     761            # print 'F', finaltime, duration
     762            msg = 'Only one of finaltime and duration may be specified'
     763            raise msg
     764        else:
     765            if finaltime is not None:
     766                self.finaltime = float(finaltime)
     767            if duration is not None:
     768                self.finaltime = self.starttime + float(duration)
     769
     770
     771
     772        N = len(self) # Number of triangles
     773        self.yieldtime = 0.0 # Track time between 'yields'
     774
     775        # Initialise interval of timestep sizes (for reporting only)
     776        self.min_timestep = max_timestep
     777        self.max_timestep = min_timestep
     778        self.number_of_steps = 0
     779        self.number_of_first_order_steps = 0
     780
     781
     782        # Update ghosts
     783        self.update_ghosts()
     784
     785        # Initial update of vertex and edge values
     786        self.distribute_to_vertices_and_edges()
     787
     788        # Update extrema if necessary (for reporting)
     789        self.update_extrema()
     790       
     791        # Initial update boundary values
     792        self.update_boundary()
     793
     794        # Or maybe restore from latest checkpoint
     795        if self.checkpoint is True:
     796            self.goto_latest_checkpoint()
     797
     798        if skip_initial_step is False:
     799            yield(self.time)  # Yield initial values
     800
     801        while True:
     802
     803            # Evolve One Step, using appropriate timestepping method
     804            if self.get_timestepping_method() == 'euler':
     805                self.evolve_one_euler_step(yieldstep,finaltime)
     806               
     807            elif self.get_timestepping_method() == 'rk2':
     808                self.evolve_one_rk2_step(yieldstep,finaltime)
     809
     810            elif self.get_timestepping_method() == 'rk3':
     811                self.evolve_one_rk3_step(yieldstep,finaltime)               
     812           
     813            # Update extrema if necessary (for reporting)
     814            self.update_extrema()           
     815
     816
     817            self.yieldtime += self.timestep
     818            self.number_of_steps += 1
     819            if self._order_ == 1:
     820                self.number_of_first_order_steps += 1
     821
     822            # Yield results
     823            if finaltime is not None and self.time >= finaltime-epsilon:
     824
     825                if self.time > finaltime:
     826                    # FIXME (Ole, 30 April 2006): Do we need this check?
     827                    # Probably not (Ole, 18 September 2008). Now changed to
     828                    # Exception
     829                    msg = 'WARNING (domain.py): time overshot finaltime. '
     830                    msg += 'Contact Ole.Nielsen@ga.gov.au'
     831                    raise Exception, msg
     832                   
     833
     834                # Yield final time and stop
     835                self.time = finaltime
     836                yield(self.time)
     837                break
     838
     839
     840            if self.yieldtime >= yieldstep:
     841                # Yield (intermediate) time and allow inspection of domain
     842
     843                if self.checkpoint is True:
     844                    self.store_checkpoint()
     845                    self.delete_old_checkpoints()
     846
     847                # Pass control on to outer loop for more specific actions
     848                yield(self.time)
     849
     850                # Reinitialise
     851                self.yieldtime = 0.0
     852                self.min_timestep = max_timestep
     853                self.max_timestep = min_timestep
     854                self.number_of_steps = 0
     855                self.number_of_first_order_steps = 0
     856                self.max_speed = zeros(N, Float)
     857
     858
     859    def evolve_one_euler_step(self, yieldstep, finaltime):
     860        """
     861        One Euler Time Step
     862        Q^{n+1} = E(h) Q^n
     863        """
     864
     865        # Compute fluxes across each element edge
     866        self.compute_fluxes()
     867
     868        # Update timestep to fit yieldstep and finaltime
     869        self.update_timestep(yieldstep, finaltime)
     870
     871        # Update conserved quantities
     872        self.update_conserved_quantities()
     873
     874        # Update ghosts
     875        self.update_ghosts()
     876
     877        # Update vertex and edge values
     878        self.distribute_to_vertices_and_edges()
     879
     880        # Update boundary values
     881        self.update_boundary()
     882
     883        # Update time
     884        self.time += self.timestep
     885
     886       
     887
     888
     889    def evolve_one_rk2_step(self, yieldstep, finaltime):
     890        """
     891        One 2nd order RK timestep
     892        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
     893        """
     894
     895        # Save initial initial conserved quantities values
     896        self.backup_conserved_quantities()           
     897
     898        #--------------------------------------
     899        # First euler step
     900        #--------------------------------------
     901
     902        # Compute fluxes across each element edge
     903        self.compute_fluxes()
     904
     905        # Update timestep to fit yieldstep and finaltime
     906        self.update_timestep(yieldstep, finaltime)
     907
     908        # Update conserved quantities
     909        self.update_conserved_quantities()
     910
     911        # Update ghosts
     912        self.update_ghosts()
     913
     914        # Update vertex and edge values
     915        self.distribute_to_vertices_and_edges()
     916
     917        # Update boundary values
     918        self.update_boundary()
     919
     920        # Update time
     921        self.time += self.timestep
     922
     923        #------------------------------------
     924        # Second Euler step
     925        #------------------------------------
     926           
     927        # Compute fluxes across each element edge
     928        self.compute_fluxes()
     929
     930        # Update conserved quantities
     931        self.update_conserved_quantities()
     932
     933        #------------------------------------
     934        # Combine initial and final values
     935        # of conserved quantities and cleanup
     936        #------------------------------------
     937       
     938        # Combine steps
     939        self.saxpy_conserved_quantities(0.5, 0.5)
     940
     941        #-----------------------------------
     942        # clean up vertex values
     943        #-----------------------------------
     944 
     945        # Update ghosts
     946        self.update_ghosts()
     947
     948        # Update vertex and edge values
     949        self.distribute_to_vertices_and_edges()
     950
     951        # Update boundary values
     952        self.update_boundary()
     953
     954
     955
     956    def evolve_one_rk3_step(self, yieldstep, finaltime):
     957        """
     958        One 3rd order RK timestep
     959        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
     960        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
     961        """
     962
     963        # Save initial initial conserved quantities values
     964        self.backup_conserved_quantities()           
     965
     966        initial_time = self.time
     967       
     968        #--------------------------------------
     969        # First euler step
     970        #--------------------------------------
     971
     972        # Compute fluxes across each element edge
     973        self.compute_fluxes()
     974
     975        # Update timestep to fit yieldstep and finaltime
     976        self.update_timestep(yieldstep, finaltime)
     977
     978        # Update conserved quantities
     979        self.update_conserved_quantities()
     980
     981        # Update ghosts
     982        self.update_ghosts()
     983
     984        # Update vertex and edge values
     985        self.distribute_to_vertices_and_edges()
     986
     987        # Update boundary values
     988        self.update_boundary()
     989
     990        # Update time
     991        self.time += self.timestep
     992
     993        #------------------------------------
     994        # Second Euler step
     995        #------------------------------------
     996           
     997        # Compute fluxes across each element edge
     998        self.compute_fluxes()
     999
     1000        # Update conserved quantities
     1001        self.update_conserved_quantities()
     1002
     1003        #------------------------------------
     1004        #Combine steps to obtain intermediate
     1005        #solution at time t^n + 0.5 h
     1006        #------------------------------------
     1007
     1008        # Combine steps
     1009        self.saxpy_conserved_quantities(0.25, 0.75)
     1010 
     1011        # Update ghosts
     1012        self.update_ghosts()
     1013
     1014        # Update vertex and edge values
     1015        self.distribute_to_vertices_and_edges()
     1016
     1017        # Update boundary values
     1018        self.update_boundary()
     1019
     1020        # Set substep time
     1021        self.time = initial_time + self.timestep*0.5
     1022
     1023        #------------------------------------
     1024        # Third Euler step
     1025        #------------------------------------
     1026           
     1027        # Compute fluxes across each element edge
     1028        self.compute_fluxes()
     1029
     1030        # Update conserved quantities
     1031        self.update_conserved_quantities()
     1032
     1033        #------------------------------------
     1034        # Combine final and initial values
     1035        # and cleanup
     1036        #------------------------------------
     1037       
     1038        # Combine steps
     1039        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
     1040 
     1041        # Update ghosts
     1042        self.update_ghosts()
     1043
     1044        # Update vertex and edge values
     1045        self.distribute_to_vertices_and_edges()
     1046
     1047        # Update boundary values
     1048        self.update_boundary()
     1049
     1050        # Set new time
     1051        self.time = initial_time + self.timestep       
     1052       
     1053
     1054    def backup_conserved_quantities(self):
     1055        N = len(self) # Number_of_triangles
     1056
     1057        # Backup conserved_quantities centroid values
     1058        for name in self.conserved_quantities:
     1059            Q = self.quantities[name]
     1060            Q.backup_centroid_values()       
     1061
     1062    def saxpy_conserved_quantities(self,a,b):
     1063        N = len(self) #number_of_triangles
     1064
     1065        # Backup conserved_quantities centroid values
     1066        for name in self.conserved_quantities:
     1067            Q = self.quantities[name]
     1068            Q.saxpy_centroid_values(a,b)       
     1069
     1070
     1071    #==============================
     1072    # John Jakeman's old evolve code
     1073    #=============================
     1074
     1075    def evolve_john(self, yieldstep = None, finaltime = None,
    7011076               skip_initial_step = False):
    7021077        """Evolve model from time=0.0 to finaltime yielding results
     
    7281103
    7291104        self.order = self.default_order
     1105       
    7301106        self.time_order = self.default_time_order
    7311107
     
    7411117       
    7421118        #update ghosts
    743         #self.update_ghosts()
     1119        self.update_ghosts()
    7441120   
    7451121        #Initial update of vertex and edge values
     
    9691345
    9701346
     1347    def update_ghosts(self):
     1348        pass
     1349   
    9711350    def update_boundary(self):
    9721351        """Go through list of boundary objects and update boundary values
     
    10351414
    10361415        self.timestep = timestep
    1037    
     1416
     1417    def update_extrema(self):
     1418        pass
    10381419
    10391420    def compute_forcing_terms(self):
  • anuga_work/development/anuga_1d/quantity.py

    r5536 r5738  
    5151        #Allocate space for other quantities
    5252        self.centroid_values = zeros(N, Float)
     53        self.centroid_backup_values = zeros(N, Float)
    5354        #self.edge_values = zeros((N, 2), Float)
    5455        #edge values are values of the ends of each interval
     
    854855                else:
    855856                    qv[k,i] = qc[k]
    856    
     857
     858    def backup_centroid_values(self):
     859        # Call correct module function
     860        # (either from this module or C-extension)
     861        #backup_centroid_values(self)
     862        N = self.domain.number_of_elements
     863        for i in range(self.N):
     864            quantity.centroid_backup_values[i] = quantity.centroid_values[i]
     865
     866    def saxpy_centroid_values(self,a,b):
     867        # Call correct module function
     868        # (either from this module or C-extension)
     869        #saxpy_centroid_values(self,a,b)
     870        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]       
    857873
    858874class Conserved_quantity(Quantity):
     
    869885        print "Use Quantity instead of Conserved_quantity"
    870886
    871                    
     887
    872888
    873889def newLinePlot(title='Simple Plot'):
  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5731 r5738  
    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 )
     
    134134        domain.default_order = 1
    135135        domain.default_time_order = 1
    136         yieldstep=10.0
    137         finaltime=10.0
     136        yieldstep=1.0
     137        finaltime=1.0
    138138
    139139        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
     
    168168        domain.default_order = 2
    169169        domain.default_time_order = 1
    170         yieldstep=10.0
    171         finaltime=10.0
     170        yieldstep=1.0
     171        finaltime=1.0
    172172
    173173        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
     
    177177        assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
    178178
    179     def test_evolve_euler_second_order_space_time(self):
     179    def test_evolve_second_order_space_time(self):
    180180        """
    181181        Compare still lake solution for various versions of shallow_water_domain
     
    192192        domain.default_order = 2
    193193        domain.default_time_order = 2
    194         yieldstep=10.0
    195         finaltime=10.0
     194        yieldstep=1.0
     195        finaltime=1.0
    196196
    197197        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
Note: See TracChangeset for help on using the changeset viewer.