Changeset 5738

Timestamp:

Sep 5, 2008, 3:06:01 PM (15 years ago)

Author:

steve

Message:

Copied evolve code from anuga_2d

Location:

anuga_work/development/anuga_1d

Files:

• compile_all.py (modified) (1 diff)
• dam_h_elevation.py (modified) (6 diffs)
• domain.py (modified) (5 diffs)
• quantity.py (modified) (3 diffs)
• test_shallow_water.py (modified) (5 diffs)

anuga_work/development/anuga_1d/compile_all.py

@@ -6 +6 @@
 
 anuga_dir = '..' + os.sep + '..' + os.sep + '..'
-#os.chdir(anuga_dir)
-print os.getcwd()  
-
 utilities_dir = anuga_dir + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities'
 
-print anuga_dir
-print utilities_dir
 execfile( utilities_dir + os.sep + 'compile.py')
 
 
+os.chdir(buildroot)    
 
-os.chdir(buildroot)    
-#execfile('test_all.py')
     

anuga_work/development/anuga_1d/dam_h_elevation.py

@@ -42 +42 @@
 
 domain.default_order = 2
-domain.default_time_order = 2
+domain.default_time_order = 1
 domain.cfl = 1.0
 domain.limiter = "vanleer"
@@ -61 +61 @@
     for i in range(len(x)):
         if 1100.0 <= x[i] <=1200.0:
-            y[i]=10.001
+            y[i]=10.0
         else:
             y[i]=10.0
@@ -140 +140 @@
         from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
         #print 'Test1'
+        ion()
         hold(False)
         #clf()
@@ -150 +151 @@
 
         plot(X,ElevationQ,X,HeightQ)
-        plot1.set_ylim([-1,12])
+        plot1.set_ylim([9.99,10.01])
         xlabel('Position')
         ylabel('Stage')
@@ -161 +162 @@
         plot2 = subplot(212)
         plot(X,MomentumQ)
-        plot2.set_ylim([-1,1])
+        #plot2.set_ylim([-1,1])
         
         #legend( ('Numerical Solution', 'for momentum'),
@@ -173 +174 @@
         file += ".eps"
         #savefig(file)
-        show()
+        #show()
         
 print 'That took %.2f seconds'%(time.time()-t0)
+
+raw_input("Press any key")

anuga_work/development/anuga_1d/domain.py

@@ -697 +697 @@
         self.datadir = name
 
-#Main components of evolve
-    def evolve(self, yieldstep = None, finaltime = None,
+
+    #--------------------------
+    # Main components of evolve
+    #--------------------------    
+
+    def evolve(self,
+               yieldstep = None,
+               finaltime = None,
+               duration = None,
+               skip_initial_step = False):
+        """Evolve model through time starting from self.starttime.
+
+
+        yieldstep: Interval between yields where results are stored,
+                   statistics written and domain inspected or
+                   possibly modified. If omitted the internal predefined
+                   max timestep is used.
+                   Internally, smaller timesteps may be taken.
+
+        duration: Duration of simulation
+
+        finaltime: Time where simulation should end. This is currently
+        relative time.  So it's the same as duration.
+
+        If both duration and finaltime are given an exception is thrown.
+
+
+        skip_initial_step: Boolean flag that decides whether the first
+        yield step is skipped or not. This is useful for example to avoid
+        duplicate steps when multiple evolve processes are dove tailed.
+
+
+        Evolve is implemented as a generator and is to be called as such, e.g.
+
+        for t in domain.evolve(yieldstep, finaltime):
+            <Do something with domain and t>
+
+
+        All times are given in seconds
+
+        """
+
+        from config import min_timestep, max_timestep, epsilon
+
+        # FIXME: Maybe lump into a larger check prior to evolving
+        msg = 'Boundary tags must be bound to boundary objects before '
+        msg += 'evolving system, '
+        msg += 'e.g. using the method set_boundary.\n'
+        msg += 'This system has the boundary tags %s '\
+               %self.get_boundary_tags()
+        assert hasattr(self, 'boundary_objects'), msg
+
+
+        if yieldstep is None:
+            yieldstep = max_timestep
+        else:
+            yieldstep = float(yieldstep)
+
+        self._order_ = self.default_order
+
+
+        if finaltime is not None and duration is not None:
+            # print 'F', finaltime, duration
+            msg = 'Only one of finaltime and duration may be specified'
+            raise msg
+        else:
+            if finaltime is not None:
+                self.finaltime = float(finaltime)
+            if duration is not None:
+                self.finaltime = self.starttime + float(duration)
+
+
+
+        N = len(self) # Number of triangles
+        self.yieldtime = 0.0 # Track time between 'yields'
+
+        # Initialise interval of timestep sizes (for reporting only)
+        self.min_timestep = max_timestep
+        self.max_timestep = min_timestep
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Initial update of vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update extrema if necessary (for reporting)
+        self.update_extrema()
+        
+        # Initial update boundary values
+        self.update_boundary()
+
+        # Or maybe restore from latest checkpoint
+        if self.checkpoint is True:
+            self.goto_latest_checkpoint()
+
+        if skip_initial_step is False:
+            yield(self.time)  # Yield initial values
+
+        while True:
+
+            # Evolve One Step, using appropriate timestepping method
+            if self.get_timestepping_method() == 'euler':
+                self.evolve_one_euler_step(yieldstep,finaltime)
+                
+            elif self.get_timestepping_method() == 'rk2':
+                self.evolve_one_rk2_step(yieldstep,finaltime)
+
+            elif self.get_timestepping_method() == 'rk3':
+                self.evolve_one_rk3_step(yieldstep,finaltime)               
+            
+            # Update extrema if necessary (for reporting)
+            self.update_extrema()            
+
+
+            self.yieldtime += self.timestep
+            self.number_of_steps += 1
+            if self._order_ == 1:
+                self.number_of_first_order_steps += 1
+
+            # Yield results
+            if finaltime is not None and self.time >= finaltime-epsilon:
+
+                if self.time > finaltime:
+                    # FIXME (Ole, 30 April 2006): Do we need this check?
+                    # Probably not (Ole, 18 September 2008). Now changed to
+                    # Exception
+                    msg = 'WARNING (domain.py): time overshot finaltime. '
+                    msg += 'Contact Ole.Nielsen@ga.gov.au'
+                    raise Exception, msg
+                    
+
+                # Yield final time and stop
+                self.time = finaltime
+                yield(self.time)
+                break
+
+
+            if self.yieldtime >= yieldstep:
+                # Yield (intermediate) time and allow inspection of domain
+
+                if self.checkpoint is True:
+                    self.store_checkpoint()
+                    self.delete_old_checkpoints()
+
+                # Pass control on to outer loop for more specific actions
+                yield(self.time)
+
+                # Reinitialise
+                self.yieldtime = 0.0
+                self.min_timestep = max_timestep
+                self.max_timestep = min_timestep
+                self.number_of_steps = 0
+                self.number_of_first_order_steps = 0
+                self.max_speed = zeros(N, Float)
+
+
+    def evolve_one_euler_step(self, yieldstep, finaltime):
+        """
+        One Euler Time Step
+        Q^{n+1} = E(h) Q^n
+        """
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        
+
+
+    def evolve_one_rk2_step(self, yieldstep, finaltime):
+        """
+        One 2nd order RK timestep
+        Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n
+        """
+
+        # Save initial initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        #--------------------------------------
+        # First euler step
+        #--------------------------------------
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        #------------------------------------
+        # Second Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        # Combine initial and final values
+        # of conserved quantities and cleanup
+        #------------------------------------
+        
+        # Combine steps
+        self.saxpy_conserved_quantities(0.5, 0.5)
+
+        #-----------------------------------
+        # clean up vertex values
+        #-----------------------------------
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+
+
+    def evolve_one_rk3_step(self, yieldstep, finaltime):
+        """
+        One 3rd order RK timestep
+        Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n  (at time t^n + h/2)
+        Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1})
+        """
+
+        # Save initial initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        initial_time = self.time
+        
+        #--------------------------------------
+        # First euler step
+        #--------------------------------------
+
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Update time
+        self.time += self.timestep
+
+        #------------------------------------
+        # Second Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine steps to obtain intermediate
+        #solution at time t^n + 0.5 h
+        #------------------------------------
+
+        # Combine steps
+        self.saxpy_conserved_quantities(0.25, 0.75)
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Set substep time
+        self.time = initial_time + self.timestep*0.5
+
+        #------------------------------------
+        # Third Euler step
+        #------------------------------------
+            
+        # Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        # Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        # Combine final and initial values
+        # and cleanup
+        #------------------------------------
+        
+        # Combine steps
+        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
+ 
+        # Update ghosts
+        self.update_ghosts()
+
+        # Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        # Update boundary values
+        self.update_boundary()
+
+        # Set new time
+        self.time = initial_time + self.timestep       
+        
+
+    def backup_conserved_quantities(self):
+        N = len(self) # Number_of_triangles
+
+        # Backup conserved_quantities centroid values
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.backup_centroid_values()        
+
+    def saxpy_conserved_quantities(self,a,b):
+        N = len(self) #number_of_triangles
+
+        # Backup conserved_quantities centroid values
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            Q.saxpy_centroid_values(a,b)        
+
+
+    #==============================
+    # John Jakeman's old evolve code
+    #=============================
+
+    def evolve_john(self, yieldstep = None, finaltime = None,
                skip_initial_step = False):
         """Evolve model from time=0.0 to finaltime yielding results
@@ -728 +1103 @@
 
         self.order = self.default_order
+        
         self.time_order = self.default_time_order
 
@@ -741 +1117 @@
         
         #update ghosts
-        #self.update_ghosts()
+        self.update_ghosts()
     
         #Initial update of vertex and edge values
@@ -969 +1345 @@
 
 
+    def update_ghosts(self):
+        pass
+    
     def update_boundary(self):
         """Go through list of boundary objects and update boundary values
@@ -1035 +1414 @@
 
         self.timestep = timestep
-    
+
+    def update_extrema(self):
+        pass
 
     def compute_forcing_terms(self):

anuga_work/development/anuga_1d/quantity.py

@@ -51 +51 @@
         #Allocate space for other quantities
         self.centroid_values = zeros(N, Float)
+        self.centroid_backup_values = zeros(N, Float)
         #self.edge_values = zeros((N, 2), Float)
         #edge values are values of the ends of each interval
@@ -854 +855 @@
                 else:
                     qv[k,i] = qc[k]
-    
+
+    def backup_centroid_values(self):
+        # Call correct module function
+        # (either from this module or C-extension)
+        #backup_centroid_values(self)
+        N = self.domain.number_of_elements
+        for i in range(self.N):
+            quantity.centroid_backup_values[i] = quantity.centroid_values[i]
+
+    def saxpy_centroid_values(self,a,b):
+        # Call correct module function
+        # (either from this module or C-extension)
+        #saxpy_centroid_values(self,a,b)
+        N = self.domain.number_of_elements
+        for i in range(self.N):
+            quantity.centroid_backup_values[i] = a*quantity.centroid_values[i] + b*quantity.centroid_backup_values[i]        
 
 class Conserved_quantity(Quantity):
@@ -869 +885 @@
         print "Use Quantity instead of Conserved_quantity"
 
-                    
+
 
 def newLinePlot(title='Simple Plot'):

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -41 +41 @@
         #print domain.quantities['stage'].explicit_update
         #print domain.quantities['xmomentum'].explicit_update
-        print stage_ud
+        #print stage_ud
 
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
@@ -134 +134 @@
         domain.default_order = 1
         domain.default_time_order = 1
-        yieldstep=10.0
-        finaltime=10.0
+        yieldstep=1.0
+        finaltime=1.0
 
         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
@@ -168 +168 @@
         domain.default_order = 2
         domain.default_time_order = 1
-        yieldstep=10.0
-        finaltime=10.0
+        yieldstep=1.0
+        finaltime=1.0
 
         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
@@ -177 +177 @@
         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
 
-    def test_evolve_euler_second_order_space_time(self):
+    def test_evolve_second_order_space_time(self):
         """
         Compare still lake solution for various versions of shallow_water_domain
@@ -192 +192 @@
         domain.default_order = 2
         domain.default_time_order = 2
-        yieldstep=10.0
-        finaltime=10.0
+        yieldstep=1.0
+        finaltime=1.0
 
         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):