Changeset 1463

Timestamp:

May 25, 2005, 5:24:14 PM (19 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• parallel/advection.py (modified) (2 diffs)
• parallel/run_parallel_advection.py (modified) (2 diffs)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/parallel/advection.py

@@ -90 +90 @@
         return flux, max_speed
 
-
     def compute_fluxes(self):
+
+        self.compute_fluxes_python()
+
+    def compute_fluxes_python(self):
         """Compute all fluxes and the timestep suitable for all volumes
         in domain.
@@ -178 +181 @@
         self.timestep = timestep
 
-
+    def compute_fluxes_weave(self):
+        """Compute all fluxes and the timestep suitable for all volumes
+        in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        domain.explicit_update
+
+        The maximal allowable speed computed by the flux_function
+        for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+        domain.explicit_update is reset to computed flux values
+        domain.timestep is set to the largest step satisfying all volumes.
+        """
+
+        import sys
+        from Numeric import zeros, Float
+        from config import max_timestep
+
+        import weave
+        from weave import converters
+
+        N = self.number_of_elements
+
+        local_timestep = float(max_timestep) #FIXME: Get rid of this
+
+        #Shortcuts
+        Stage = self.quantities['stage']
+
+        #Arrays
+        neighbours      = self.neighbours
+        neighbour_edges = self.neighbour_edges
+        normals         = self.normals
+        areas           = self.areas
+        radii           = self.radii
+        edgelengths     = self.edgelengths
+
+        stage_edge      = Stage.edge_values
+        stage_bdry      = Stage.boundary_values
+        stage_update    = Stage.explicit_update
+
+        huge_timestep = float(sys.maxint)
+
+        v1 = self.velocity[0]
+        v2 = self.velocity[1]
+
+        code = """
+        //Loop
+
+        double qr,ql;
+        int m,n;
+        double normal[2];
+        double normal_velocity;
+        double flux, edgeflux;
+        double max_speed;
+        double optimal_timestep;
+        for (int k=0; k<N; k++){
+
+            optimal_timestep = huge_timestep;
+            flux = 0.0;  //Reset work array
+            for (int i=0; i<3; i++){
+                //Quantities inside volume facing neighbour i
+                ql = stage_edge(k, i);
+
+                //Quantities at neighbour on nearest face
+                n = neighbours(k,i);
+                if (n < 0) {
+                    m = -n-1; //Convert neg flag to index
+                    qr = stage_bdry(m);
+                } else {
+                    m = neighbour_edges(k,i);
+                    qr = stage_edge(n, m);
+                }
+
+
+                //Outward pointing normal vector
+                for (int j=0; j<2; j++){
+                    normal[j] = normals(k, 2*i+j);
+                }
+
+
+                //Flux computation using provided function
+                normal_velocity = v1*normal[0] + v2*normal[1];
+
+                if (normal_velocity < 0) {
+                    edgeflux = qr * normal_velocity;
+                } else {
+                    edgeflux = ql * normal_velocity;
+                }
+
+                max_speed = fabs(normal_velocity);
+                flux = flux - edgeflux * edgelengths(k,i);
+
+                //Update optimal_timestep
+                if (max_speed != 0.0) {
+                    optimal_timestep = (optimal_timestep>radii(k)/max_speed) ? radii(k)/max_speed : optimal_timestep;
+                }
+
+            }
+
+            //Normalise by area and store for when all conserved
+            //quantities get updated
+            stage_update(k) = flux/areas(k);
+
+            local_timestep = (local_timestep>optimal_timestep) ? optimal_timestep : local_timestep;
+            std::cout << local_timestep << std::endl;
+        }
+        """
+
+        weave.inline(code, ['stage_edge','stage_bdry','stage_update',
+                             'neighbours','neighbour_edges','normals',
+                             'areas','radii','edgelengths','huge_timestep',
+                             'local_timestep','v1','v2','N'],
+                             type_converters = converters.blitz, compiler='gcc');
+
+        self.timestep = local_timestep
+        print local_timestep
 
     def evolve(self, yieldstep = None, finaltime = None):

inundation/ga/storm_surge/parallel/run_parallel_advection.py

@@ -19 +19 @@
 
 #N = 50
-N = 50
-M = 50
+N = 10
+M = 10
 
 points, vertices, boundary, full_send_dict, ghost_recv_dict = parallel_rectangular(N, M, len1=1.0)
@@ -60 +60 @@
 
 #Check that the boundary value gets propagated to all elements
-for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
+for t in domain.evolve(yieldstep = 0.05, finaltime = 3.0):
     if myid == 0:
         domain.write_time()