Changeset 2716 for development/pyvolution-1d/domain.py

Timestamp:

Apr 13, 2006, 4:34:22 PM (18 years ago)

Author:

jakeman

Message:

Adding new files

File:

• development/pyvolution-1d/domain.py (modified) (21 diffs)

development/pyvolution-1d/domain.py

@@ -8 +8 @@
 """
 from generic_boundary_conditions import *
+from coordinate_transforms.geo_reference import Geo_reference
 
 class Domain:
 
     def __init__(self, coordinates, boundary = None,
-                 conserved_quantities = None,
-                 other_quantities = None,
-                 tagged_elements = None):
+                 conserved_quantities = None, other_quantities = None,
+                 tagged_elements = None, geo_reference = None):
         """
         Build 1D elements from x coordinates
@@ -25 +25 @@
         self.coordinates = array(coordinates)
 
+        if geo_reference is None:
+            self.geo_reference = Geo_reference() #Use defaults 
+        else:
+            self.geo_reference = geo_reference
+
         #Register number of Elements
         self.number_of_elements = N = len(self.coordinates)-1
@@ -30 +35 @@
         #Allocate space for neighbour and boundary structures
         self.neighbours = zeros((N, 2), Int)
-        self.neighbour_edges = zeros((N, 2), Int)
+        #self.neighbour_edges = zeros((N, 2), Int)
+        self.neighbour_vertices = zeros((N, 2), Int)
         self.number_of_boundaries = zeros(N, Int)
         self.surrogate_neighbours = zeros((N, 2), Int)
@@ -60 +66 @@
             self.neighbours[i, :] = [-1, -1]
             #Initialise edge ids of neighbours
+            #Initialise vertex ids of neighbours
             #In case of boundaries this slot is not used
-            self.neighbour_edges[i, :] = [-1, -1]
+            #self.neighbour_edges[i, :] = [-1, -1]
+            self.neighbour_vertices[i, :] = [-1, -1]
+
+        self.build_vertexlist()
 
         #Build neighbour structure
@@ -70 +80 @@
 
         #Build boundary dictionary mapping (id, edge) to symbolic tags
+        #Build boundary dictionary mapping (id, vertex) to symbolic tags
         self.build_boundary_dictionary(boundary)
 
         #Build tagged element  dictionary mapping (tag) to array of elements
         self.build_tagged_elements_dictionary(tagged_elements)
-
-
         
         from quantity import Quantity, Conserved_quantity
@@ -105 +114 @@
         self.forcing_terms = []
 
+        #Defaults
+        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
+        self.beta_w = beta_w
+        self.beta_h = beta_h
+        self.epsilon = epsilon
+
+        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
+        #Or maybe get rid of order altogether and use beta_w and beta_h
+        self.default_order = 1
+        self.order = self.default_order
+        self.smallsteps = 0
+        self.max_smallsteps = max_smallsteps
+        self.number_of_steps = 0
+        self.number_of_first_order_steps = 0
+        self.CFL = CFL
+
+        #Model time
+        self.time = 0.0
+        self.finaltime = None
+        self.min_timestep = self.max_timestep = 0.0
+        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        #Checkpointing and storage
+        from config import default_datadir
+        self.datadir = default_datadir
+        self.filename = 'domain'
+        self.checkpoint = False
+
+    def __len__(self):
+        return self.number_of_elements
+
+    def build_vertexlist(self):
+        """Build vertexlist index by vertex ids and for each entry (point id)
+        build a list of (triangles, vertex_id) pairs that use the point
+        as vertex.
+
+        Preconditions:
+          self.coordinates and self.triangles are defined
+
+        Postcondition:
+          self.vertexlist is built
+        """
+        from Numeric import array
+
+        vertexlist = [None]*len(self.coordinates)
+        for i in range(self.number_of_elements):
+
+            #a = self.triangles[i, 0]
+            #b = self.triangles[i, 1]
+            #c = self.triangles[i, 2]
+            a = i
+            b = i + 1
+
+            #Register the vertices v as lists of
+            #(triangle_id, vertex_id) tuples associated with them
+            #This is used for smoothing
+            #for vertex_id, v in enumerate([a,b,c]):
+            for vertex_id, v in enumerate([a,b]):
+                if vertexlist[v] is None:
+                    vertexlist[v] = []
+
+                vertexlist[v].append( (i, vertex_id) )
+
+        self.vertexlist = vertexlist
+
+        
     def build_neighbour_structure(self):
         """Update all registered triangles to point to their neighbours.
@@ -112 +186 @@
         Postconditions:
           neighbours and neighbour_edges is populated
+          neighbours and neighbour_vertices is populated
           number_of_boundaries integer array is defined.
         """
@@ -121 +196 @@
         N = self.number_of_elements
         neighbourdict = {}
-        l_edge = 0
-        r_edge = 1
+        #l_edge = 0
+        #r_edge = 1
+        l_vertex = 0
+        r_vertex = 1
         for i in range(N):
 
@@ -149 +226 @@
             #neighbourdict[a,b] = (i, 1) #(id, edge)
             #neighbourdict[b,a] = (i, 0) #(id, edge)
-            neighbourdict[a,l_edge] = (i, 0) #(id, edge)
-            neighbourdict[b,r_edge] = (i, 1) #(id, edge)
+            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
+            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
+            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
+            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
 
 
@@ -166 +245 @@
             #self.number_of_boundaries[i] = 3
             self.number_of_boundaries[i] = 2
-            if neighbourdict.has_key((b,l_edge)):
-                self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
-                self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
+            #if neighbourdict.has_key((b,l_edge)):
+            if neighbourdict.has_key((b,l_vertex)):
+                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
+                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
+                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
+                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
                 self.number_of_boundaries[i] -= 1
 
-            if neighbourdict.has_key((a,r_edge)):
-                self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
-                self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
+            #if neighbourdict.has_key((a,r_edge)):
+            if neighbourdict.has_key((a,r_vertex)):
+                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
+                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
+                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
+                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
                 self.number_of_boundaries[i] -= 1
+                
             #if neighbourdict.has_key((b,a)):
             #    self.neighbours[i, 1] = neighbourdict[b,a][0]
@@ -231 +317 @@
             for vol_id in range(self.number_of_elements):
                 #for edge_id in range(0, 3):
-                for edge_id in range(0, 2):
-                    if self.neighbours[vol_id, edge_id] < 0:
-                        boundary[(vol_id, edge_id)] = default_boundary_tag
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:   
+                        #boundary[(vol_id, edge_id)] = default_boundary_tag
+                        boundary[(vol_id, vertex_id)] = default_boundary_tag
         else:
             #Check that all keys in given boundary exist
-            for vol_id, edge_id in boundary.keys():
-                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+            #for vol_id, edge_id in boundary.keys():
+            for vol_id, vertex_id in boundary.keys():    
+                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
+                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
                 a, b = self.neighbours.shape
-                assert vol_id < a and edge_id < b, msg
+                #assert vol_id < a and edge_id < b, msg
+                assert vol_id < a and vertex_id < b, msg
 
                 #FIXME: This assert violates internal boundaries (delete it)
@@ -248 +340 @@
             for vol_id in range(self.number_of_elements):
                 #for edge_id in range(0, 3):
-                for edge_id in range(0, 2):    
-                    if self.neighbours[vol_id, edge_id] < 0:
-                        if not boundary.has_key( (vol_id, edge_id) ):
+                #for edge_id in range(0, 2):
+                for vertex_id in range(0, 2):    
+                    #if self.neighbours[vol_id, edge_id] < 0:
+                    if self.neighbours[vol_id, vertex_id] < 0:    
+                        #if not boundary.has_key( (vol_id, edge_id) ):
+                        if not boundary.has_key( (vol_id, vertex_id) ):    
                             msg = 'WARNING: Given boundary does not contain '
-                            msg += 'tags for edge (%d, %d). '\
-                                   %(vol_id, edge_id)
+                            #msg += 'tags for edge (%d, %d). '\
+                            #       %(vol_id, edge_id)
+                            msg += 'tags for vertex (%d, %d). '\
+                                   %(vol_id, vertex_id)
                             msg += 'Assigning default tag (%s).'\
                                    %default_boundary_tag
@@ -264 +361 @@
                             #enable default boundary-tags where
                             #tags a not specified
-                            boundary[ (vol_id, edge_id) ] =\
+                            #boundary[ (vol_id, edge_id) ] =\
+                            boundary[ (vol_id, vertex_id) ] =\
                                       default_boundary_tag
 
@@ -304 +402 @@
         return tags.keys()
         
-
-    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
+    def get_vertex_coordinates(self, obj = False):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+
+        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
+        FIXME, we might make that the default.
+        FIXME Maybe use keyword: continuous = False for this condition?
+
+        
+        """
+
+        if obj is True:
+            from Numeric import concatenate, reshape
+            #V = self.vertex_coordinates
+            V = self.vertices
+            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
+
+            N = V.shape[0]
+            #return reshape(V, (3*N, 2))
+            return reshape(V, (N, 2))
+        else:    
+            #return self.vertex_coordinates
+            return self.vertices
+    
+    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
         """Get conserved quantities at volume vol_id
 
@@ -319 +441 @@
         from Numeric import zeros, Float
 
-        if not (vertex is None or edge is None):
-            msg = 'Values for both vertex and edge was specified.'
-            msg += 'Only one (or none) is allowed.'
-            raise msg
+        #if not (vertex is None):# or edge is None):
+        #    msg = 'Values for both vertex and edge was specified.'
+        #    msg += 'Only one (or none) is allowed.'
+       #     raise msg
 
         q = zeros( len(self.conserved_quantities), Float)
@@ -330 +452 @@
             if vertex is not None:
                 q[i] = Q.vertex_values[vol_id, vertex]
-            elif edge is not None:
-                q[i] = Q.edge_values[vol_id, edge]
+            #elif edge is not None:
+            #    q[i] = Q.edge_values[vol_id, edge]
             else:
                 q[i] = Q.centroid_values[vol_id]
@@ -368 +490 @@
 
         return self.areas[elem_id]
+
+    def get_quantity(self, name, location='vertices', indices = None):
+        """Get values for named quantity
+
+        name: Name of quantity
+
+        In case of location == 'centroids' the dimension values must
+        be a list of a Numerical array of length N, N being the number
+        of elements. Otherwise it must be of dimension Nx3.
+
+        Indices is the set of element ids that the operation applies to.
+
+        The values will be stored in elements following their
+        internal ordering.
+        """
+
+        return self.quantities[name].get_values( location, indices = indices)
 
     def set_quantity(self, name, *args, **kwargs):
@@ -462 +601 @@
         #Loop through edges that lie on the boundary and associate them with
         #callable boundary objects depending on their tags
-        for k, (vol_id, edge_id) in enumerate(x):
-            tag = self.boundary[ (vol_id, edge_id) ]
+        #for k, (vol_id, edge_id) in enumerate(x):
+        for k, (vol_id, vertex_id) in enumerate(x):
+            #tag = self.boundary[ (vol_id, edge_id) ]
+            tag = self.boundary[ (vol_id, vertex_id) ]
 
             if boundary_map.has_key(tag):
@@ -469 +610 @@
 
                 if B is not None:
-                    self.boundary_objects.append( ((vol_id, edge_id), B) )
-                    self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
+                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
+                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
+                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
+                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
                 else:
                     pass
@@ -494 +637 @@
         ##assert hasattr(self, 'boundary_objects')
 
+    def get_name(self):
+        return self.filename
+
+    def set_name(self, name):
+        self.filename = name
+
+    def get_datadir(self):
+        return self.datadir
+
+    def set_datadir(self, name):
+        self.datadir = name
+
+#Main components of evolve
+
+    def evolve(self, yieldstep = None, finaltime = None,
+               skip_initial_step = False):
+        """Evolve model from time=0.0 to finaltime yielding results
+        every yieldstep.
+
+        Internally, smaller timesteps may be taken.
+
+        Evolve is implemented as a generator and is to be called as such, e.g.
+
+        for t in domain.evolve(timestep, yieldstep, finaltime):
+            <Do something with domain and t>
+
+        """
+
+        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 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
+
+        ##self.set_defaults()
+
+        if yieldstep is None:
+            yieldstep = max_timestep
+        else:
+            yieldstep = float(yieldstep)
+
+        self.order = self.default_order
+
+
+        self.yieldtime = 0.0 #Time between 'yields'
+
+        #Initialise interval of timestep sizes (for reporting only)
+        self.min_timestep = max_timestep
+        self.max_timestep = min_timestep
+        self.finaltime = finaltime
+        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()
+
+        #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:
+
+            #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
+            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 abs(self.time - finaltime) < epsilon:
+
+                #FIXME: There is a rare situation where the
+                #final time step is stored twice. Can we make a test?
+
+                # Yield final time and stop
+                yield(self.time)
+                break
+
+
+            if abs(self.yieldtime - yieldstep) < epsilon:
+                # 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
+
+    def distribute_to_vertices_and_edges(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices and edge-midpoints for each volume
+
+        Default implementation is straight first order,
+        i.e. constant values throughout each element and
+        no reference to non-conserved quantities.
+        """
+
+        for name in self.conserved_quantities:
+            Q = self.quantities[name]
+            if self.order == 1:
+                Q.extrapolate_first_order()
+            elif self.order == 2:
+                Q.extrapolate_second_order()
+                Q.limit()
+            else:
+                raise 'Unknown order'
+            Q.interpolate_from_vertices_to_edges()
+
+
     def update_boundary(self):
         """Go through list of boundary objects and update boundary values
@@ -501 +791 @@
         #FIXME: Update only those that change (if that can be worked out)
         #FIXME: Boundary objects should not include ghost nodes.
-        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
-            q = B.evaluate(vol_id, edge_id)
+        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
+        #    q = B.evaluate(vol_id, edge_id)
+        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
+            q = B.evaluate(vol_id, vertex_id)
 
             for j, name in enumerate(self.conserved_quantities):
                 Q = self.quantities[name]
                 Q.boundary_values[i] = q[j]
+
+    def compute_forcing_terms(self):
+        """If there are any forcing functions driving the system
+        they should be defined in Domain subclass and appended to
+        the list self.forcing_terms
+        """
+
+        for f in self.forcing_terms:
+            f(self)