Changeset 4712

Timestamp:

Sep 7, 2007, 10:43:46 PM (17 years ago)

Author:

steve

Message:

Working on 2nd order time stepping

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (9 diffs)
• abstract_2d_finite_volumes/quantity.py (modified) (6 diffs)
• abstract_2d_finite_volumes/quantity_ext.c (modified) (4 diffs)
• abstract_2d_finite_volumes/test_domain.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_quantity.py (modified) (1 diff)
• abstract_2d_finite_volumes/util.py (modified) (2 diffs)
• config.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (1 diff)
• shallow_water/test_system.py (modified) (6 diffs)
• utilities/compile.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -31 +31 @@
 
 import types
+
+
 
 class Domain(Mesh):
@@ -176 +178 @@
         from anuga.config import max_smallsteps, beta_w, beta_h, epsilon
         from anuga.config import CFL
+        from anuga.config import timestepping_method
         from anuga.config import protect_against_isolated_degenerate_timesteps
         self.beta_w = beta_w
@@ -182 +185 @@
         self.protect_against_isolated_degenerate_timesteps = protect_against_isolated_degenerate_timesteps
         
+        
 
         # FIXME: Maybe have separate orders for h-limiter and w-limiter?
@@ -194 +198 @@
         self.number_of_first_order_steps = 0
         self.CFL = CFL
-
+        self.set_timestepping_method(timestepping_method)
+        
         self.boundary_map = None  # Will be populated by set_boundary        
         
@@ -954 +959 @@
         return msg
 
+    def get_timestepping_method(self):
+        return self.timestepping_method
+
+    def set_timestepping_method(self,timestepping_method):
+        
+        if timestepping_method in ['euler', 'rk2']:
+            self.timestepping_method = timestepping_method
+            return
+
+        msg = '%s is an incorrect timestepping type'% timestepping_method
+        raise Exception, msg
 
     def get_name(self):
@@ -973 +989 @@
     def set_datadir(self, name):
         self.datadir = name
+
+    def get_starttime(self):
+        return self.starttime
+
+    def set_starttime(self, time):
+        self.starttime = float(time)        
 
 
@@ -1085 +1107 @@
 
         while True:
-            # 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()
+
+            # 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)               
 
             # Update extrema if necessary (for reporting)
             self.update_extrema()            
 
-            # Update boundary values
-            self.update_boundary()
-
-            # Update time
-            self.time += self.timestep
+
             self.yieldtime += self.timestep
             self.number_of_steps += 1
@@ -1145 +1156 @@
 
 
+    def evolve_one_euler_step(self, yieldstep, finaltime):
+        """One Euler Time 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
+
+        
+
+
+    def evolve_one_rk2_step(self,yieldstep, finaltime):
+        """One 2nd order RK timestep"""
+
+        #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
+        #------------------------------------
+
+        self.saxpy_conserved_quantities(0.5, 0.5)
+
+        #------------------------------------
+        #Clean up rk step
+        #------------------------------------
+            
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+
     def evolve_to_end(self, finaltime = 1.0):
         """Iterate evolve all the way to the end
@@ -1153 +1253 @@
 
 
+    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)        
+    
 
     def update_boundary(self):

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1305 +1305 @@
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
+        self.centroid_backup_values = zeros(N, Float)       
+
+        
 
 
@@ -1318 +1321 @@
         return compute_gradients(self)
 
-
     def limit(self):
         #Call correct module function
@@ -1324 +1326 @@
         limit(self)
 
-
     def extrapolate_second_order(self):
         #Call correct module function
@@ -1330 +1331 @@
         extrapolate_second_order(self)
 
+    def backup_centroid_values(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        backup_centroid_values(self)
+
+    def saxpy_centroid_values(self,a,b):
+        #Call correct module function
+        #(either from this module or C-extension)
+        saxpy_centroid_values(self,a,b)
+    
 
 def update(quantity, timestep):
@@ -1417 +1428 @@
 
 
-def extrapolate_second_order(quantity):
-    """Extrapolate conserved quantities from centroid to
-    vertices for each volume using
-    second order scheme.
-    """
-
-    a, b = quantity.compute_gradients()
-
-    X = quantity.domain.get_vertex_coordinates()
+def backup_centroid_values(quantity):
+    """Copy centroid values to backup array"""
+
     qc = quantity.centroid_values
-    qv = quantity.vertex_values
+    qb = quantity.centroid_backup_values
 
     #Check each triangle
     for k in range(len(quantity.domain)):
-        #Centroid coordinates
-        x, y = quantity.domain.centroid_coordinates[k]
-
-        #vertex coordinates
-        x0, y0, x1, y1, x2, y2 = X[k,:]
-
-        #Extrapolate
-        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
-        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
+        qb[k] = qc[k]
+
+
+def saxpy_centroid_values(quantity,a,b):
+    """saxpy operation between centroid value and backup"""
+
+    qc = quantity.centroid_values
+    qb = quantity.centroid_backup_values
+
+    #Check each triangle
+    for k in range(len(quantity.domain)):
+        qc[k] = a*qc[k]+b*qb[k]       
 
 
@@ -1582 +1589 @@
     #Replace python version with c implementations
 
-    from quantity_ext import average_vertex_values
+    from quantity_ext import average_vertex_values, backup_centroid_values, \
+         saxpy_centroid_values
 
     from quantity_ext import compute_gradients, limit,\

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -100 +100 @@
 
 
+
 int _extrapolate(int N,
                  double* centroids,
@@ -162 +163 @@
         return 0;
 }
+
+int _backup_centroid_values(int N,
+                            double* centroid_values,
+                            double* centroid_backup_values) {
+    //Backup centroid values
+
+
+    int k;
+
+    for (k=0; k<N; k++) {
+        centroid_backup_values[k] = centroid_values[k];
+    }
+
+
+    return 0;
+}
+
+
+int _saxpy_centroid_values(int N,
+                           double a,
+                           double b,
+                           double* centroid_values,
+                           double* centroid_backup_values) {
+    //saxby centroid values
+
+
+    int k;
+
+
+    for (k=0; k<N; k++) {
+        centroid_values[k] = a*centroid_values[k] + b*centroid_backup_values[k];
+    }
+
+
+    return 0;
+}
+
 
 int _update(int N,
@@ -305 +343 @@
 
 
+PyObject *backup_centroid_values(PyObject *self, PyObject *args) {
+
+        PyObject *quantity;
+        PyArrayObject *centroid_values, *centroid_backup_values;
+
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "O", &quantity)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: backup_centroid_values could not parse input");
+          return NULL;
+        }
+
+        centroid_values        = get_consecutive_array(quantity, "centroid_values");
+        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
+
+        N = centroid_values -> dimensions[0];
+
+        err = _backup_centroid_values(N,
+                      (double*) centroid_values -> data,
+                      (double*) centroid_backup_values -> data);
+
+
+        //Release and return
+        Py_DECREF(centroid_values);
+        Py_DECREF(centroid_backup_values);
+
+        return Py_BuildValue("");
+}
+
+PyObject *saxpy_centroid_values(PyObject *self, PyObject *args) {
+
+        PyObject *quantity;
+        PyArrayObject *centroid_values, *centroid_backup_values;
+
+        double a,b;
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "Odd", &quantity, &a, &b)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: saxpy_centroid_values could not parse input");
+          return NULL;
+        }
+
+        centroid_values        = get_consecutive_array(quantity, "centroid_values");
+        centroid_backup_values = get_consecutive_array(quantity, "centroid_backup_values");
+
+        N = centroid_values -> dimensions[0];
+
+        err = _saxpy_centroid_values(N,a,b,
+                      (double*) centroid_values -> data,
+                      (double*) centroid_backup_values -> data);
+
+
+        //Release and return
+        Py_DECREF(centroid_values);
+        Py_DECREF(centroid_backup_values);
+
+        return Py_BuildValue("");
+}
+
+
 PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
 
@@ -645 +749 @@
         {"limit", limit, METH_VARARGS, "Print out"},
         {"update", update, METH_VARARGS, "Print out"},
+        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
+        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
         {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
         {"extrapolate_second_order", extrapolate_second_order,

anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

@@ -319 +319 @@
 
 
-
-
-
+    def test_setting_timestepping_method(self):
+        """test_set_quanitities_to_be_monitored
+        """
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        domain = Domain(points, vertices, boundary=None,
+                        conserved_quantities =\
+                        ['stage', 'xmomentum', 'ymomentum'],
+                        other_quantities = ['elevation', 'friction', 'depth'])
+
+
+        domain.timestepping_method = None
+
+
+        # Check that invalid requests are dealt with
+        try:
+            domain.set_timestepping_method('eee')        
+        except:
+            pass
+        else:
+            msg = 'Should have caught illegal method'
+            raise Exception, msg
+
+
+        #Should have no trouble with euler or rk2
+        domain.set_timestepping_method('euler')
+        domain.set_timestepping_method('rk2')
+
+        #test get timestepping method
+        assert domain.get_timestepping_method() == 'rk2'
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -1368 +1368 @@
 
 
+    def test_backup_saxpy_centroid_values(self):
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Set up for a gradient of (3,1), f(x) = 3x+y
+        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
+        d_values = array([1.0, 2.0, 3.0, 4.0])
+        quantity.set_values(c_values, location = 'centroids')
+
+        #Backup
+        quantity.backup_centroid_values()
+
+        #print quantity.vertex_values
+        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
+
+
+        quantity.set_values(d_values, location = 'centroids')
+
+        quantity.saxpy_centroid_values(2.0, 3.0)
+
+        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
+
+
+
     def test_first_order_extrapolator(self):
         quantity = Conserved_quantity(self.mesh4)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -85 +85 @@
             quantities = domain.conserved_quantities
             
-        domain_starttime = domain.starttime
+        domain_starttime = domain.get_starttime()
     else:
         domain_starttime = None
@@ -140 +140 @@
             
             if verbose: print msg
-            domain.starttime = starttime #Modifying model time
+            domain.set_starttime(starttime) #Modifying model time
             if verbose: print 'Domain starttime is now set to %f'\
                %domain.starttime

anuga_core/source/anuga/config.py

@@ -104 +104 @@
            #make changes
 
+# Choose type of timestepping,
+timestepping_method = 'euler' # 1st order euler
+#timestepping_method = 'rk2'   # 2nd Order TVD scheme
+
 # Option to search for signatures where isolated triangles are
 # responsible for a small global timestep.

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -446 +446 @@
 
 
-        #Initial update of vertex and edge values before any storage
+        #Initial update of vertex and edge values before any STORAGE
         #and or visualisation
+        #This is done again in the initialisation of the Generic_Domain
+        #evolve loop but we do it here to ensure the values are ok for storage
         self.distribute_to_vertices_and_edges()
 

anuga_core/source/anuga/shallow_water/test_system.py

@@ -47 +47 @@
         domain.set_name(boundary_name)                 
         domain.set_datadir(dir)          
-        domain.starttime = boundary_starttime
+        domain.set_starttime(boundary_starttime)
         
         # Setup initial conditions
@@ -71 +71 @@
         
         """
+     
         boundary_starttime = 500
         boundary_filename = self.create_sww_boundary(boundary_starttime)
@@ -95 +96 @@
         # Setup boundary conditions
         domain.set_boundary({'exterior': Bf})
-        for t in domain.evolve(yieldstep = 5, finaltime = 10.0):
+
+        
+        for t in domain.evolve(yieldstep = 5.0, finaltime = 10.0):
             pass
-            #domain.write_time()
+            #print domain.write_time()
             #print "domain.time", domain.time
 
@@ -117 +120 @@
         Test that starttime can be set in the middle of a boundary condition
         """
+        
         boundary_starttime = 500
         boundary_filename = self.create_sww_boundary(boundary_starttime)
@@ -134 +138 @@
         domain.set_datadir(dir)
         new_starttime = 510.
-        domain.starttime = new_starttime
+        domain.set_starttime(new_starttime)
 
         # Setup initial conditions
@@ -146 +150 @@
         for t in domain.evolve(yieldstep = 5, finaltime = 9.0):
             pass
-            #domain.write_time()
+            #print domain.write_time()
             #print "domain.time", domain.time
 

anuga_core/source/anuga/utilities/compile.py

@@ -201 +201 @@
     # This will make compile work locally
     utilities_include_dir = '.'
+    utilities_include_dir = '../utilities'