Changeset 2716


Ignore:
Timestamp:
Apr 13, 2006, 4:34:22 PM (18 years ago)
Author:
jakeman
Message:

Adding new files

Location:
development/pyvolution-1d
Files:
8 added
3 edited

Legend:

Unmodified
Added
Removed
  • development/pyvolution-1d/domain.py

    r2705 r2716  
    88"""
    99from generic_boundary_conditions import *
     10from coordinate_transforms.geo_reference import Geo_reference
    1011
    1112class Domain:
    1213
    1314    def __init__(self, coordinates, boundary = None,
    14                  conserved_quantities = None,
    15                  other_quantities = None,
    16                  tagged_elements = None):
     15                 conserved_quantities = None, other_quantities = None,
     16                 tagged_elements = None, geo_reference = None):
    1717        """
    1818        Build 1D elements from x coordinates
     
    2525        self.coordinates = array(coordinates)
    2626
     27        if geo_reference is None:
     28            self.geo_reference = Geo_reference() #Use defaults
     29        else:
     30            self.geo_reference = geo_reference
     31
    2732        #Register number of Elements
    2833        self.number_of_elements = N = len(self.coordinates)-1
     
    3035        #Allocate space for neighbour and boundary structures
    3136        self.neighbours = zeros((N, 2), Int)
    32         self.neighbour_edges = zeros((N, 2), Int)
     37        #self.neighbour_edges = zeros((N, 2), Int)
     38        self.neighbour_vertices = zeros((N, 2), Int)
    3339        self.number_of_boundaries = zeros(N, Int)
    3440        self.surrogate_neighbours = zeros((N, 2), Int)
     
    6066            self.neighbours[i, :] = [-1, -1]
    6167            #Initialise edge ids of neighbours
     68            #Initialise vertex ids of neighbours
    6269            #In case of boundaries this slot is not used
    63             self.neighbour_edges[i, :] = [-1, -1]
     70            #self.neighbour_edges[i, :] = [-1, -1]
     71            self.neighbour_vertices[i, :] = [-1, -1]
     72
     73        self.build_vertexlist()
    6474
    6575        #Build neighbour structure
     
    7080
    7181        #Build boundary dictionary mapping (id, edge) to symbolic tags
     82        #Build boundary dictionary mapping (id, vertex) to symbolic tags
    7283        self.build_boundary_dictionary(boundary)
    7384
    7485        #Build tagged element  dictionary mapping (tag) to array of elements
    7586        self.build_tagged_elements_dictionary(tagged_elements)
    76 
    77 
    7887       
    7988        from quantity import Quantity, Conserved_quantity
     
    105114        self.forcing_terms = []
    106115
     116        #Defaults
     117        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
     118        self.beta_w = beta_w
     119        self.beta_h = beta_h
     120        self.epsilon = epsilon
     121
     122        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
     123        #Or maybe get rid of order altogether and use beta_w and beta_h
     124        self.default_order = 1
     125        self.order = self.default_order
     126        self.smallsteps = 0
     127        self.max_smallsteps = max_smallsteps
     128        self.number_of_steps = 0
     129        self.number_of_first_order_steps = 0
     130        self.CFL = CFL
     131
     132        #Model time
     133        self.time = 0.0
     134        self.finaltime = None
     135        self.min_timestep = self.max_timestep = 0.0
     136        self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
     137        #Checkpointing and storage
     138        from config import default_datadir
     139        self.datadir = default_datadir
     140        self.filename = 'domain'
     141        self.checkpoint = False
     142
     143    def __len__(self):
     144        return self.number_of_elements
     145
     146    def build_vertexlist(self):
     147        """Build vertexlist index by vertex ids and for each entry (point id)
     148        build a list of (triangles, vertex_id) pairs that use the point
     149        as vertex.
     150
     151        Preconditions:
     152          self.coordinates and self.triangles are defined
     153
     154        Postcondition:
     155          self.vertexlist is built
     156        """
     157        from Numeric import array
     158
     159        vertexlist = [None]*len(self.coordinates)
     160        for i in range(self.number_of_elements):
     161
     162            #a = self.triangles[i, 0]
     163            #b = self.triangles[i, 1]
     164            #c = self.triangles[i, 2]
     165            a = i
     166            b = i + 1
     167
     168            #Register the vertices v as lists of
     169            #(triangle_id, vertex_id) tuples associated with them
     170            #This is used for smoothing
     171            #for vertex_id, v in enumerate([a,b,c]):
     172            for vertex_id, v in enumerate([a,b]):
     173                if vertexlist[v] is None:
     174                    vertexlist[v] = []
     175
     176                vertexlist[v].append( (i, vertex_id) )
     177
     178        self.vertexlist = vertexlist
     179
     180       
    107181    def build_neighbour_structure(self):
    108182        """Update all registered triangles to point to their neighbours.
     
    112186        Postconditions:
    113187          neighbours and neighbour_edges is populated
     188          neighbours and neighbour_vertices is populated
    114189          number_of_boundaries integer array is defined.
    115190        """
     
    121196        N = self.number_of_elements
    122197        neighbourdict = {}
    123         l_edge = 0
    124         r_edge = 1
     198        #l_edge = 0
     199        #r_edge = 1
     200        l_vertex = 0
     201        r_vertex = 1
    125202        for i in range(N):
    126203
     
    149226            #neighbourdict[a,b] = (i, 1) #(id, edge)
    150227            #neighbourdict[b,a] = (i, 0) #(id, edge)
    151             neighbourdict[a,l_edge] = (i, 0) #(id, edge)
    152             neighbourdict[b,r_edge] = (i, 1) #(id, edge)
     228            #neighbourdict[a,l_edge] = (i, 0) #(id, edge)
     229            #neighbourdict[b,r_edge] = (i, 1) #(id, edge)
     230            neighbourdict[a,l_vertex] = (i, 0) #(id, vertex)
     231            neighbourdict[b,r_vertex] = (i, 1) #(id, vertex)
    153232
    154233
     
    166245            #self.number_of_boundaries[i] = 3
    167246            self.number_of_boundaries[i] = 2
    168             if neighbourdict.has_key((b,l_edge)):
    169                 self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
    170                 self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
     247            #if neighbourdict.has_key((b,l_edge)):
     248            if neighbourdict.has_key((b,l_vertex)):
     249                #self.neighbours[i, 1] = neighbourdict[b,l_edge][0]
     250                #self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1]
     251                self.neighbours[i, 1] = neighbourdict[b,l_vertex][0]
     252                self.neighbour_vertices[i, 1] = neighbourdict[b,l_vertex][1]
    171253                self.number_of_boundaries[i] -= 1
    172254
    173             if neighbourdict.has_key((a,r_edge)):
    174                 self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
    175                 self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
     255            #if neighbourdict.has_key((a,r_edge)):
     256            if neighbourdict.has_key((a,r_vertex)):
     257                #self.neighbours[i, 0] = neighbourdict[a,r_edge][0]
     258                #self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1]
     259                self.neighbours[i, 0] = neighbourdict[a,r_vertex][0]
     260                self.neighbour_vertices[i, 0] = neighbourdict[a,r_vertex][1]
    176261                self.number_of_boundaries[i] -= 1
     262               
    177263            #if neighbourdict.has_key((b,a)):
    178264            #    self.neighbours[i, 1] = neighbourdict[b,a][0]
     
    231317            for vol_id in range(self.number_of_elements):
    232318                #for edge_id in range(0, 3):
    233                 for edge_id in range(0, 2):
    234                     if self.neighbours[vol_id, edge_id] < 0:
    235                         boundary[(vol_id, edge_id)] = default_boundary_tag
     319                #for edge_id in range(0, 2):
     320                for vertex_id in range(0, 2):   
     321                    #if self.neighbours[vol_id, edge_id] < 0:
     322                    if self.neighbours[vol_id, vertex_id] < 0:   
     323                        #boundary[(vol_id, edge_id)] = default_boundary_tag
     324                        boundary[(vol_id, vertex_id)] = default_boundary_tag
    236325        else:
    237326            #Check that all keys in given boundary exist
    238             for vol_id, edge_id in boundary.keys():
    239                 msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
     327            #for vol_id, edge_id in boundary.keys():
     328            for vol_id, vertex_id in boundary.keys():   
     329                #msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
     330                msg = 'Segment (%d, %d) does not exist' %(vol_id, vertex_id)
    240331                a, b = self.neighbours.shape
    241                 assert vol_id < a and edge_id < b, msg
     332                #assert vol_id < a and edge_id < b, msg
     333                assert vol_id < a and vertex_id < b, msg
    242334
    243335                #FIXME: This assert violates internal boundaries (delete it)
     
    248340            for vol_id in range(self.number_of_elements):
    249341                #for edge_id in range(0, 3):
    250                 for edge_id in range(0, 2):   
    251                     if self.neighbours[vol_id, edge_id] < 0:
    252                         if not boundary.has_key( (vol_id, edge_id) ):
     342                #for edge_id in range(0, 2):
     343                for vertex_id in range(0, 2):   
     344                    #if self.neighbours[vol_id, edge_id] < 0:
     345                    if self.neighbours[vol_id, vertex_id] < 0:   
     346                        #if not boundary.has_key( (vol_id, edge_id) ):
     347                        if not boundary.has_key( (vol_id, vertex_id) ):   
    253348                            msg = 'WARNING: Given boundary does not contain '
    254                             msg += 'tags for edge (%d, %d). '\
    255                                    %(vol_id, edge_id)
     349                            #msg += 'tags for edge (%d, %d). '\
     350                            #       %(vol_id, edge_id)
     351                            msg += 'tags for vertex (%d, %d). '\
     352                                   %(vol_id, vertex_id)
    256353                            msg += 'Assigning default tag (%s).'\
    257354                                   %default_boundary_tag
     
    264361                            #enable default boundary-tags where
    265362                            #tags a not specified
    266                             boundary[ (vol_id, edge_id) ] =\
     363                            #boundary[ (vol_id, edge_id) ] =\
     364                            boundary[ (vol_id, vertex_id) ] =\
    267365                                      default_boundary_tag
    268366
     
    304402        return tags.keys()
    305403       
    306 
    307     def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
     404    def get_vertex_coordinates(self, obj = False):
     405        """Return all vertex coordinates.
     406        Return all vertex coordinates for all triangles as an Nx6 array
     407        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
     408
     409        if obj is True, the x/y pairs are returned in a 3*N x 2 array.
     410        FIXME, we might make that the default.
     411        FIXME Maybe use keyword: continuous = False for this condition?
     412
     413       
     414        """
     415
     416        if obj is True:
     417            from Numeric import concatenate, reshape
     418            #V = self.vertex_coordinates
     419            V = self.vertices
     420            #return concatenate( (V[:,0:2], V[:,2:4], V[:,4:6]), axis=0)
     421
     422            N = V.shape[0]
     423            #return reshape(V, (3*N, 2))
     424            return reshape(V, (N, 2))
     425        else:   
     426            #return self.vertex_coordinates
     427            return self.vertices
     428   
     429    def get_conserved_quantities(self, vol_id, vertex=None):#, edge=None):
    308430        """Get conserved quantities at volume vol_id
    309431
     
    319441        from Numeric import zeros, Float
    320442
    321         if not (vertex is None or edge is None):
    322             msg = 'Values for both vertex and edge was specified.'
    323             msg += 'Only one (or none) is allowed.'
    324             raise msg
     443        #if not (vertex is None):# or edge is None):
     444        #    msg = 'Values for both vertex and edge was specified.'
     445        #    msg += 'Only one (or none) is allowed.'
     446       #     raise msg
    325447
    326448        q = zeros( len(self.conserved_quantities), Float)
     
    330452            if vertex is not None:
    331453                q[i] = Q.vertex_values[vol_id, vertex]
    332             elif edge is not None:
    333                 q[i] = Q.edge_values[vol_id, edge]
     454            #elif edge is not None:
     455            #    q[i] = Q.edge_values[vol_id, edge]
    334456            else:
    335457                q[i] = Q.centroid_values[vol_id]
     
    368490
    369491        return self.areas[elem_id]
     492
     493    def get_quantity(self, name, location='vertices', indices = None):
     494        """Get values for named quantity
     495
     496        name: Name of quantity
     497
     498        In case of location == 'centroids' the dimension values must
     499        be a list of a Numerical array of length N, N being the number
     500        of elements. Otherwise it must be of dimension Nx3.
     501
     502        Indices is the set of element ids that the operation applies to.
     503
     504        The values will be stored in elements following their
     505        internal ordering.
     506        """
     507
     508        return self.quantities[name].get_values( location, indices = indices)
    370509
    371510    def set_quantity(self, name, *args, **kwargs):
     
    462601        #Loop through edges that lie on the boundary and associate them with
    463602        #callable boundary objects depending on their tags
    464         for k, (vol_id, edge_id) in enumerate(x):
    465             tag = self.boundary[ (vol_id, edge_id) ]
     603        #for k, (vol_id, edge_id) in enumerate(x):
     604        for k, (vol_id, vertex_id) in enumerate(x):
     605            #tag = self.boundary[ (vol_id, edge_id) ]
     606            tag = self.boundary[ (vol_id, vertex_id) ]
    466607
    467608            if boundary_map.has_key(tag):
     
    469610
    470611                if B is not None:
    471                     self.boundary_objects.append( ((vol_id, edge_id), B) )
    472                     self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
     612                    #self.boundary_objects.append( ((vol_id, edge_id), B) )
     613                    #self.neighbours[vol_id, edge_id] = -len(self.boundary_objects)
     614                    self.boundary_objects.append( ((vol_id, vertex_id), B) )
     615                    self.neighbours[vol_id, vertex_id] = -len(self.boundary_objects)
    473616                else:
    474617                    pass
     
    494637        ##assert hasattr(self, 'boundary_objects')
    495638
     639    def get_name(self):
     640        return self.filename
     641
     642    def set_name(self, name):
     643        self.filename = name
     644
     645    def get_datadir(self):
     646        return self.datadir
     647
     648    def set_datadir(self, name):
     649        self.datadir = name
     650
     651#Main components of evolve
     652
     653    def evolve(self, yieldstep = None, finaltime = None,
     654               skip_initial_step = False):
     655        """Evolve model from time=0.0 to finaltime yielding results
     656        every yieldstep.
     657
     658        Internally, smaller timesteps may be taken.
     659
     660        Evolve is implemented as a generator and is to be called as such, e.g.
     661
     662        for t in domain.evolve(timestep, yieldstep, finaltime):
     663            <Do something with domain and t>
     664
     665        """
     666
     667        from config import min_timestep, max_timestep, epsilon
     668
     669        #FIXME: Maybe lump into a larger check prior to evolving
     670        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
     671        msg += 'e.g. using the method set_boundary.\n'
     672        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
     673        assert hasattr(self, 'boundary_objects'), msg
     674
     675        ##self.set_defaults()
     676
     677        if yieldstep is None:
     678            yieldstep = max_timestep
     679        else:
     680            yieldstep = float(yieldstep)
     681
     682        self.order = self.default_order
     683
     684
     685        self.yieldtime = 0.0 #Time between 'yields'
     686
     687        #Initialise interval of timestep sizes (for reporting only)
     688        self.min_timestep = max_timestep
     689        self.max_timestep = min_timestep
     690        self.finaltime = finaltime
     691        self.number_of_steps = 0
     692        self.number_of_first_order_steps = 0
     693
     694        #update ghosts
     695        #self.update_ghosts()
     696
     697        #Initial update of vertex and edge values
     698        self.distribute_to_vertices_and_edges()
     699
     700        #Initial update boundary values
     701        self.update_boundary()
     702
     703        #Or maybe restore from latest checkpoint
     704        if self.checkpoint is True:
     705            self.goto_latest_checkpoint()
     706
     707        if skip_initial_step is False:
     708            yield(self.time)  #Yield initial values
     709
     710        while True:
     711
     712            #Compute fluxes across each element edge
     713            self.compute_fluxes()
     714
     715            #Update timestep to fit yieldstep and finaltime
     716            self.update_timestep(yieldstep, finaltime)
     717
     718            #Update conserved quantities
     719            self.update_conserved_quantities()
     720
     721            #update ghosts
     722            #self.update_ghosts()
     723
     724            #Update vertex and edge values
     725            self.distribute_to_vertices_and_edges()
     726
     727            #Update boundary values
     728            self.update_boundary()
     729
     730            #Update time
     731            self.time += self.timestep
     732            self.yieldtime += self.timestep
     733            self.number_of_steps += 1
     734            if self.order == 1:
     735                self.number_of_first_order_steps += 1
     736
     737            #Yield results
     738            if finaltime is not None and abs(self.time - finaltime) < epsilon:
     739
     740                #FIXME: There is a rare situation where the
     741                #final time step is stored twice. Can we make a test?
     742
     743                # Yield final time and stop
     744                yield(self.time)
     745                break
     746
     747
     748            if abs(self.yieldtime - yieldstep) < epsilon:
     749                # Yield (intermediate) time and allow inspection of domain
     750
     751                if self.checkpoint is True:
     752                    self.store_checkpoint()
     753                    self.delete_old_checkpoints()
     754
     755                #Pass control on to outer loop for more specific actions
     756                yield(self.time)
     757
     758                # Reinitialise
     759                self.yieldtime = 0.0
     760                self.min_timestep = max_timestep
     761                self.max_timestep = min_timestep
     762                self.number_of_steps = 0
     763                self.number_of_first_order_steps = 0
     764
     765    def distribute_to_vertices_and_edges(self):
     766        """Extrapolate conserved quantities from centroid to
     767        vertices and edge-midpoints for each volume
     768
     769        Default implementation is straight first order,
     770        i.e. constant values throughout each element and
     771        no reference to non-conserved quantities.
     772        """
     773
     774        for name in self.conserved_quantities:
     775            Q = self.quantities[name]
     776            if self.order == 1:
     777                Q.extrapolate_first_order()
     778            elif self.order == 2:
     779                Q.extrapolate_second_order()
     780                Q.limit()
     781            else:
     782                raise 'Unknown order'
     783            Q.interpolate_from_vertices_to_edges()
     784
     785
    496786    def update_boundary(self):
    497787        """Go through list of boundary objects and update boundary values
     
    501791        #FIXME: Update only those that change (if that can be worked out)
    502792        #FIXME: Boundary objects should not include ghost nodes.
    503         for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
    504             q = B.evaluate(vol_id, edge_id)
     793        #for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
     794        #    q = B.evaluate(vol_id, edge_id)
     795        for i, ((vol_id, vertex_id), B) in enumerate(self.boundary_objects):
     796            q = B.evaluate(vol_id, vertex_id)
    505797
    506798            for j, name in enumerate(self.conserved_quantities):
    507799                Q = self.quantities[name]
    508800                Q.boundary_values[i] = q[j]
     801
     802    def compute_forcing_terms(self):
     803        """If there are any forcing functions driving the system
     804        they should be defined in Domain subclass and appended to
     805        the list self.forcing_terms
     806        """
     807
     808        for f in self.forcing_terms:
     809            f(self)
    509810
    510811
  • development/pyvolution-1d/shallow_water_1d.py

    r2705 r2716  
    4949class Domain(Generic_Domain):
    5050
    51     def __init__(self, coordinates, boundary = None, tagged_elements = None):
     51    def __init__(self, coordinates, boundary = None, tagged_elements = None,
     52                 geo_reference = None):
    5253
    5354        conserved_quantities = ['stage', 'xmomentum']
     
    5556        Generic_Domain.__init__(self, coordinates, boundary,
    5657                                conserved_quantities, other_quantities,
    57                                 tagged_elements)
     58                                tagged_elements, geo_reference)
    5859       
    5960        from config import minimum_allowed_height, g
     
    281282# Flux computation
    282283#def flux_function(normal, ql, qr, zl, zr):
    283 def flux_function(ql, qr, zl, zr):
     284def flux_function(normal, ql, qr, zl, zr):
    284285    """Compute fluxes between volumes for the shallow water wave equation
    285286    cast in terms of w = h+z using the 'central scheme' as described in
     
    305306    #q_right = rotate(qr, normal, direction = 1)
    306307    q_left = ql
     308    q_left[1] = q_left[1]*normal
    307309    q_right = qr
     310    q_right[1] = q_right[1]*normal
    308311
    309312    z = (zl+zr)/2 #Take average of field values
     
    338341
    339342    #Maximal wave speed
    340     #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    341     s_max = max(u_left+soundspeed_left, u_right+soundspeed_right,0)
     343    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    342344   
    343345    #Minimal wave speed
    344     #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    345     s_min = min(u_left-soundspeed_left, u_right-soundspeed_right,0)
     346    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    346347   
    347348    #Flux computation
     
    363364    denom = s_max-s_min
    364365    if denom == 0.0:
    365         edgeflux = array([0.0, 0.0, 0.0])
     366        #edgeflux = array([0.0, 0.0, 0.0])
    366367        edgeflux = array([0.0, 0.0])
    367368        max_speed = 0.0
     
    369370        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
    370371        edgeflux += s_max*s_min*(q_right-q_left)/denom
     372       
     373        edgeflux[1] = -edgeflux[1]*normal
    371374
    372375        #edgeflux = rotate(edgeflux, normal, direction=-1)
     
    407410
    408411    #Arrays
    409     stage = Stage.edge_values
    410     xmom =  Xmom.edge_values
     412    #stage = Stage.edge_values
     413    #xmom =  Xmom.edge_values
    411414 #   ymom =  Ymom.edge_values
    412     bed =   Bed.edge_values
     415    #bed =   Bed.edge_values
     416   
     417    stage = Stage.vertex_values
     418    xmom =  Xmom.vertex_values
     419    bed =   Bed.vertex_values
     420   
     421    print 'stage edge values'
     422    print stage
     423    print 'xmom edge values'
     424    print xmom
     425    print 'bed values'
     426    print bed
    413427
    414428    stage_bdry = Stage.boundary_values
    415429    xmom_bdry =  Xmom.boundary_values
     430    print 'stage_bdry'
     431    print stage_bdry
     432    print 'xmom_bdry'
     433    print xmom_bdry
    416434#    ymom_bdry =  Ymom.boundary_values
    417435
     
    439457                zr = zl #Extend bed elevation to boundary
    440458            else:
    441                 m = domain.neighbour_edges[k,i]
     459                #m = domain.neighbour_edges[k,i]
     460                m = domain.neighbour_vertices[k,i]
    442461                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
    443462                qr = [stage[n, m], xmom[n, m]]
     
    446465
    447466            #Outward pointing normal vector
    448             normal = domain.normals[k, 2*i:2*i+2]
     467            #normal = domain.normals[k, 2*i:2*i+2]
    449468            #CHECK THIS LINE [k, 2*i:2*i+1]
    450469
    451470            #Flux computation using provided function
    452471            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    453             edgeflux, max_speed = flux_function(ql, qr, zl, zr)
    454             flux -= edgeflux * domain.edgelengths[k,i]
     472            print 'ql'
     473            print ql
     474            print 'qr'
     475            print qr
     476           
     477            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     478            print 'edgeflux'
     479            print edgeflux
     480            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
     481            # flux = edgefluxleft - edgefluxright
     482            flux -= edgeflux #* domain.edgelengths[k,i]
    455483
    456484            #Update optimal_timestep
    457485            try:
    458                 timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
     486                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
     487                timestep = min(timestep, 0.5/max_speed)
     488                print 'check time step in compute fluxes is ok'
    459489            except ZeroDivisionError:
    460490                pass
     
    514544    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
    515545                                     domain.neighbours,
    516                                      domain.neighbour_edges,
     546                                     #domain.neighbour_edges,
     547                                     domain.neighbour_vertices,
    517548                                     domain.normals,
    518                                      domain.edgelengths,
    519                                      domain.radii,
     549                                     #domain.edgelengths,
     550                                     #domain.radii,
    520551                                     domain.areas,
    521                                      Stage.edge_values,
    522                                      Xmom.edge_values,
     552                                     #Stage.edge_values,
     553                                     #Xmom.edge_values,
    523554                                     #Ymom.edge_values,
    524                                      Bed.edge_values,
     555                                     #Bed.edge_values,
     556                                     Stage.vertex_values,
     557                                     Xmom.vertex_values,
     558                                     Bed.vertex_values,
    525559                                     Stage.boundary_values,
    526560                                     Xmom.boundary_values,
     
    568602        #MH090605 if second order,
    569603        #perform the extrapolation and limiting on
    570         #all of the conserved quantitie
     604        #all of the conserved quantities
    571605
    572606        if (domain.order == 1):
     
    595629
    596630    #Compute edge values by interpolation
    597     for name in domain.conserved_quantities:
    598         Q = domain.quantities[name]
    599         Q.interpolate_from_vertices_to_edges()
     631    #for name in domain.conserved_quantities:
     632    #    Q = domain.quantities[name]
     633    #    Q.interpolate_from_vertices_to_edges()
    600634
    601635
     
    898932
    899933        #Handy shorthands
    900         self.stage   = domain.quantities['stage'].edge_values
    901         self.xmom    = domain.quantities['xmomentum'].edge_values
     934        #self.stage   = domain.quantities['stage'].edge_values
     935        #self.xmom    = domain.quantities['xmomentum'].edge_values
    902936        #self.ymom    = domain.quantities['ymomentum'].edge_values
    903937        #self.normals = domain.normals
     938        self.stage   = domain.quantities['stage'].vertex_values
     939        self.xmom    = domain.quantities['xmomentum'].vertex_values
    904940
    905941        from Numeric import zeros, Float
     
    931967        r[1] = -q[1]
    932968        q = r
     969        #For start interval there is no outward momentum so do not need to
     970        #reverse direction in this case
    933971
    934972        return q
     
    950988    Stage = domain.quantities['stage']
    951989    Elevation = domain.quantities['elevation']
    952     h = Stage.edge_values - Elevation.edge_values
     990    #h = Stage.edge_values - Elevation.edge_values
     991    h = Stage.vertex_values - Elevation.vertex_values
    953992    v = Elevation.vertex_values
    954993
     
    10051044
    10061045    uh = domain.quantities['xmomentum'].centroid_values
    1007     vh = domain.quantities['ymomentum'].centroid_values
     1046    #vh = domain.quantities['ymomentum'].centroid_values
    10081047    eta = domain.quantities['friction'].centroid_values
    10091048
    10101049    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
    1011     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
     1050    #ymom_update = domain.quantities['ymomentum'].semi_implicit_update
    10121051
    10131052    N = domain.number_of_elements
     
    10181057        if eta[k] >= eps:
    10191058            if h[k] >= eps:
    1020                 S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
     1059                #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
     1060                S = -g * eta[k]**2 * uh[k]
    10211061                S /= h[k]**(7.0/3)
    10221062
    10231063                #Update momentum
    10241064                xmom_update[k] += S*uh[k]
    1025                 ymom_update[k] += S*vh[k]
     1065                #ymom_update[k] += S*vh[k]
    10261066
    10271067
  • development/pyvolution-1d/test_quantity.py

    r2229 r2716  
    112112            raise 'should have raised AssertionError'
    113113
     114
    114115    def test_set_values_const(self):
    115116        quantity = Quantity(self.domain5)
     
    120121        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
    121122
    122 
    123123        quantity.set_values(2.0, location = 'centroids')
    124124        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
    125 
    126125
    127126    def test_set_values_func(self):
Note: See TracChangeset for help on using the changeset viewer.