Changeset 3335


Ignore:
Timestamp:
Jul 14, 2006, 2:11:20 PM (18 years ago)
Author:
jakeman
Message:

New file dam.py tests numerical solution against analytical solution of
Stoker 57. Numerical results match Stoker's roughly but oscillations
present. Most likely limiting is not working correctly. Q.limit does
have an effect but unsure whether this effect is the right one.

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

Legend:

Unmodified
Added
Removed

  • development/pyvolution-1d/domain.py

    r3295 r3335  
    557557            Q = self.create_quantity_from_expression(expression)
    558558            kwargs['quantity'] = Q
    559 
    560 
    561         print "self.vertices1"
    562         print self.vertices
    563559       
    564560        #Assign values
    565561        self.quantities[name].set_values(*args, **kwargs)
    566 
    567         print "self.vertices2"
    568         print self.vertices
    569562
    570563    def set_boundary(self, boundary_map):
    … …  
    663656        ##assert hasattr(self, 'boundary_objects')
    664657
     658    def write_time(self):
     659        print self.timestepping_statistics()
     660
     661    def timestepping_statistics(self):
     662        """Return string with time stepping statistics for printing or logging
     663        """
     664
     665        msg = ''
     666        if self.min_timestep == self.max_timestep:
     667            msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
     668                   %(self.time, self.min_timestep, self.number_of_steps,
     669                     self.number_of_first_order_steps)
     670        elif self.min_timestep > self.max_timestep:
     671            msg += 'Time = %.4f, steps=%d (%d)'\
     672                   %(self.time, self.number_of_steps,
     673                     self.number_of_first_order_steps)
     674        else:
     675            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
     676                   %(self.time, self.min_timestep,
     677                     self.max_timestep, self.number_of_steps,
     678                     self.number_of_first_order_steps)
     679
     680        return msg
     681
    665682    def get_name(self):
    666683        return self.filename
    … …  
    676693
    677694#Main components of evolve
    678 
    679695    def evolve(self, yieldstep = None, finaltime = None,
    680696               skip_initial_step = False):
    … …  
    721737        #update ghosts
    722738        #self.update_ghosts()
    723 
     739   
    724740        #Initial update of vertex and edge values
    725741        self.distribute_to_vertices_and_edges()
    726 
     742       
    727743        #Initial update boundary values
    728744        self.update_boundary()
    … …  
    861877                else:
    862878                    #Try to overcome situation by switching to 1 order
     879                    print "changing Order 1"
    863880                    self.order = 1
    864881

  • development/pyvolution-1d/quantity.py

    r3322 r3335  
    134134                  Default is "vertices"
    135135        """
    136 
    137         print "In set_function value"
    138         print self.domain.vertices
    139136       
    140137        if location == 'centroids':
    141 
    142             print "In set at centroid"           
     138           
    143139            P = self.domain.centroids
    144140            self.set_values(f(P), location)
    … …  
    146142            #Vertices
    147143            P = self.domain.get_vertices()
    148 
    149             print "P"
    150             print P
    151             print P[:,0]
    152             print f(P[:,0])
    153 
    154             print P[:,1]
    155             print f(P[:,1])
     144           
    156145            for i in range(2):
    157                 self.vertex_values[:,i] = f(P[:,i])
    158 
    159 
    160            
    161 
    162             print "In set at vertices"
    163             print self.vertex_values
    164 
    165         print "Out set_function value"
    166         print self.domain.vertices       
    167 
     146                self.vertex_values[:,i] = f(P[:,i])       
    168147
    169148    def set_array_values(self, values, location='vertices'):
    … …  
    564543        vertex values are updated
    565544        """
    566 
     545        print "in Q.limit"
    567546        from Numeric import zeros, Float
    568547
    … …  
    706685    print Q2.centroid_values
    707686    print Qc
    708 
     687    raw_input('press_return')
     688   
    709689    g3 = newLinePlot('plot 3')
    710690    linePlot(g3,Xc,Qc)
    … …  
    720700    raw_input('press return')
    721701
    722 
    723 
    724 
    725 
    726 
    727 

  • development/pyvolution-1d/shallow_water_1d.py

    r3293 r3335  
    6363
    6464        #forcing terms not included in 1d domain ?WHy?
    65         #self.forcing_terms.append(gravity)
    66         #self.forcing_terms.append(manning_friction)
     65        self.forcing_terms.append(gravity)
     66        self.forcing_terms.append(manning_friction)
    6767        #print "\nI have Removed forcing terms line 64 1dsw"
    6868
    … …  
    174174        #and or visualisation
    175175        self.distribute_to_vertices_and_edges()
    176 
     176       
    177177        #Initialise real time viz if requested
    178178        #if self.visualise is True and self.time == 0.0:
    … …  
    432432    #Loop
    433433    timestep = float(sys.maxint)
     434    enter = True
    434435    for k in range(N):
    435436
    … …  
    473474            # flux = edgefluxleft - edgefluxright
    474475            flux -= edgeflux #* domain.edgelengths[k,i]
    475 
    476476            #Update optimal_timestep
    477477            try:
    478478                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
    479479                #timestep = 0.01
     480           
    480481                timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
    481                 if timestep < 0.00001:
    482                     #print 'max_speed', max_speed
     482                if (timestep < 1e-6) & (enter == True):
     483                    #print "domain.order", domain.order
     484                    #domain.write_time()
     485                    print "cell number", k
     486                    print "timestep", timestep
     487                    print 'max_speed', max_speed
    483488                    s = domain.quantities['stage']
    484489                    s = s.centroid_values
    485490                    xm = domain.quantities['xmomentum']
    486491                    xm = xm.centroid_values
    487                     #print 'h', s[k]
    488                     #print 'xm', xm[k]
    489                     #print 'u', xm[k]/s[k]
    490                     #break
    491                 #timestep = 0.01
    492                 #print 'areas', domain.areas[k]
    493                 #print "timestep", timestep
     492                    print 'h', s[k]
     493                    print 'xm', xm[k]
     494                    print 'u', xm[k]/s[k]
     495                    enter = False
     496           
    494497            except ZeroDivisionError:
    495498                pass
    … …  
    497500        #Normalise by area and store for when all conserved
    498501        #quantities get updated
    499         #flux /= domain.areas[k]
     502        flux /= domain.areas[k]
    500503        # ADD ABOVE LINE AGAIN
    501504        Stage.explicit_update[k] = flux[0]
    … …  
    599602    """
    600603
    601     from config import optimised_gradient_limiter
     604    #from config import optimised_gradient_limiter
    602605
    603606    #Remove very thin layers of water
Note: See TracChangeset for help on using the changeset viewer.