Changeset 7868


Ignore:
Timestamp:
Jun 22, 2010, 5:30:32 PM (14 years ago)
Author:
steve
Message:

Added in some more c based limiters

Location:
anuga_work/development/2010-projects/anuga_1d
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py

    r7860 r7868  
    1717                 coordinates,
    1818                 boundary = None,
    19                  conserved_quantities = None,
    20                  evolved_quantities = None,
    21                  other_quantities = None,
     19                 conserved_quantities = [],
     20                 evolved_quantities = [],
     21                 other_quantities = [],
    2222                 tagged_elements = None):
    2323        """
     
    3535        self.number_of_elements = N = len(self.coordinates)-1
    3636
    37         self.beta = 1.0
    38         self.set_limiter("minmod_kurganov")
    3937        self.set_CFL(CFL)
    4038        self.set_timestepping_method(timestepping_method)
     
    869867   
    870868
     869    def get_beta(self):
     870
     871        warn('limiter parameter beta associated with quantity not domain')
     872
     873    def set_beta(self,beta):
     874        """Set the same limiter beta parameter to all evolving quantities
     875        """
     876
     877        for name in self.evolved_quantities:
     878            Q = self.quantities[name]
     879            Q.set_beta(beta)
     880
     881
    871882    def get_limiter(self):
    872         return self.limiter
     883
     884        warn('limiter associated with quantity not domain')
    873885
    874886    def set_limiter(self,limiter):
    875887
    876         possible_limiters = \
    877         ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
    878 
    879         if limiter in possible_limiters:
    880             self.limiter = limiter
    881             return
    882 
    883         msg = '%s is an incorrect limiter type.\n'% limiter
    884         msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
    885         raise Exception, msg
     888        for name in self.evolved_quantities:
     889            Q = self.quantities[name]
     890            Q.set_limiter(limiter)
     891       
    886892
    887893       
  • anuga_work/development/2010-projects/anuga_1d/base/generic_mesh.py

    r7852 r7868  
    1818
    1919
    20     points  = x_0 + (x_1 - x_0)*numpy.arange(n+1,dtype=numpy.float)
     20    points  = x_0 + (x_1 - x_0)*numpy.arange(n+1,dtype=numpy.float)/n
    2121    boundary = {(0, 0): 'left', (n-1, 1): 'right'}
    2222
  • anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py

    r7855 r7868  
    3737def vanleer(a,b):
    3838
    39     from numpy import abs, where
     39    from numpy import fabs, where
    4040
    41     return where((abs(a) + abs(b) >= 1.0e-12), (a*abs(b)+abs(a)*b)/(abs(a)+abs(b)), 0.0)
     41    return where((fabs(a) + fabs(b) >= 1.0e-12), (a*fabs(b)+fabs(a)*b)/(fabs(a)+fabs(b)), 0.0)
    4242
    4343
  • anuga_work/development/2010-projects/anuga_1d/base/quantity.py

    r7855 r7868  
    7070        #flux calculations and forcing functions
    7171       
    72         N = domain.number_of_elements
     72        self.N = domain.number_of_elements
     73        N = self.N
    7374        self.explicit_update = numpy.zeros(N, numpy.float )
    7475        self.semi_implicit_update = numpy.zeros(N, numpy.float )
     
    7879        self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
    7980
    80         self.beta = domain.beta
     81        #These are taken from domain but can be set for each quantity
     82        # if so desired
     83        self.beta = 0.0
     84        self.limiter = 'vanleer'
    8185
    8286
     
    604608        return integral
    605609
     610    def get_beta(self,beta):
     611        """Get limiting parameter
     612        """
     613       
     614        return self.beta
     615
     616    def set_beta(self,beta):
     617        """Set limiting parameter
     618        """
     619
     620        #Probably should test that it is not too large
     621        self.beta = beta
     622
     623
     624    def get_limiter(self):
     625        return self.limiter
     626
     627    def set_limiter(self,limiter):
     628
     629        possible_limiters = \
     630        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
     631
     632        if limiter in possible_limiters:
     633            self.limiter = limiter
     634            return
     635
     636        msg = '%s is an incorrect limiter type.\n'% limiter
     637        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
     638        raise Exception, msg
    606639
    607640    def update(self, timestep):
     
    795828
    796829        #Gradient
    797         G[1:-1] = xmic( self.domain.beta, d1, d2 )     
     830        G[1:-1] = xmic( self.beta, d1, d2 )     
    798831
    799832       
     
    818851        """
    819852
    820         if self.domain.limiter == "pyvolution":
     853        if self.limiter == "pyvolution":
    821854            self.limit_pyvolution()
    822855
    823         elif self.domain.limiter == "minmod_steve":
     856        elif self.limiter == "minmod_steve":
    824857            self.limit_minmod()
    825858
     
    9851018        from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
    9861019
    987         limiter = self.domain.limiter
     1020        limiter = self.get_limiter()
    9881021        #print limiter
    9891022       
    9901023        #print 'limit_range'
    991         N = self.domain.number_of_elements
     1024        N = self.N
    9921025        qc = self.centroid_values
    9931026        qv = self.vertex_values
     
    10131046        phi[N-1] = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2])
    10141047
     1048
    10151049        if limiter == "minmod":
    10161050            phi[1:-1] = minmod(beta_p[1:-1],beta_m[1:-1])
     
    10231057
    10241058        elif limiter == "minmod_kurganov":
    1025             theta = 2.0
     1059            theta = self.beta
    10261060            phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
    10271061
     
    10341068            msg = 'Unknown limiter'
    10351069            raise Exception, msg
    1036    
     1070
     1071
     1072       
    10371073        qv[:,0] = qc + phi*dx[:,0]
    10381074        qv[:,1] = qc + phi*dx[:,1]
    10391075
    10401076
    1041         #Phi limiter
    1042 #        for k in range(N):
    1043 #            n0 = self.domain.neighbours[k,0]
    1044 #            n1 = self.domain.neighbours[k,1]
    1045 #            if n0 < 0:
    1046 #                phi = (qc[k+1] - qc[k])/(xc[k+1] - xc[k])
    1047 #            elif n1 < 0:
    1048 #                phi = (qc[k] - qc[k-1])/(xc[k] - xc[k-1])
    1049 #            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
    1050 #            #    phi = 0.0
    1051 #            else:
    1052 #                if limiter == "minmod":
    1053 #                    phi = minmod(beta_p[k],beta_m[k])
    1054 #
    1055 #                elif limiter == "minmod_kurganov":#Change this
    1056 #                    # Also known as monotonized central difference limiter
    1057 #                    # if theta = 2.0
    1058 #                    theta = 2.0
    1059 #                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
    1060 #
    1061 #                elif limiter == "superbee":
    1062 #                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
    1063 #                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
    1064 #                    phi = maxmod(slope1,slope2)
    1065 #
    1066 #                elif limiter == "vanleer":
    1067 #                    phi = vanleer(beta_p[k],beta_m[k])
    1068 #
    1069 #                elif limiter == "vanalbada":
    1070 #                    phi = vanalbada(beta_m[k],beta_p[k])
    1071 #
    1072 #            for i in range(2):
    1073 #                qv[k,i] = qc[k] + phi*dx[k,i]
     1077
    10741078
    10751079    def limit_steve_slope(self):   
  • anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

    r7860 r7868  
    376376
    377377        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
    378        
     378
     379        quantity.beta = 2.0
     380        quantity.limiter = "minmod_kurganov"
    379381        quantity.extrapolate_second_order()
    380382
     
    405407
    406408        #Work out the others
    407 
     409        quantity.beta = 2.0
     410        quantity.limiter = "minmod_kurganov"
    408411        quantity.extrapolate_second_order()
    409412       
     
    496499
    497500
     501
     502
     503    def test_limit_minmod(self):
     504        quantity = Quantity(self.domain4)
     505
     506        #Create a deliberate overshoot (e.g. from gradient computation)
     507        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     508
     509        #Limit
     510        quantity.limiter = 'minmod'
     511        quantity.limit_range()
     512
     513
     514
     515
     516        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
     517                                                 [  2.4,   9.6],
     518                                                 [ 10.6,  11.4],
     519                                                 [ 11.4,  14.6]] )
     520
     521        from anuga_1d.base.quantity_ext import limit_minmod_ext
     522
     523        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     524
     525        limit_minmod_ext(quantity)
     526
     527
     528
     529
     530        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
     531                                                 [  2.4,   9.6],
     532                                                 [ 10.6,  11.4],
     533                                                 [ 11.4,  14.6]] )
     534
     535
     536
     537    def test_limit_minmod_kurganov(self):
     538        quantity = Quantity(self.domain4)
     539
     540        #Create a deliberate overshoot (e.g. from gradient computation)
     541        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     542
     543        #Limit
     544        quantity.limiter = 'minmod_kurganov'
     545        quantity.beta = 0.5
     546        quantity.limit_range()
     547
     548
     549        #print quantity.vertex_values
     550
     551        assert allclose(quantity.vertex_values, [[ -2.4,   2.4],
     552                                                 [  4.2,   7.8],
     553                                                 [ 10.8,  11.2],
     554                                                 [ 11.4,  14.6]])
     555
     556        from anuga_1d.base.quantity_ext import limit_minmod_kurganov_ext
     557
     558        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     559
     560        limit_minmod_kurganov_ext(quantity)
     561
     562
     563        #print quantity.vertex_values
     564
     565
     566    def test_limit_vanleer(self):
     567        quantity = Quantity(self.domain4)
     568
     569        #Create a deliberate overshoot (e.g. from gradient computation)
     570        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     571
     572        #Limit
     573        quantity.set_limiter('vanleer')
     574        quantity.limit_range()
     575
     576
     577        #print quantity.vertex_values
     578
     579        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
     580                                                 [  2.32653061,   9.67346939],
     581                                                 [ 10.39393939,  11.60606061],
     582                                                 [ 11.4,         14.6       ]] )
     583
     584        from anuga_1d.base.quantity_ext import limit_vanleer_ext
     585
     586        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
     587
     588        limit_vanleer_ext(quantity)
     589
     590
     591        #print quantity.vertex_values
     592
     593        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
     594                                                 [  2.32653061,   9.67346939],
     595                                                 [ 10.39393939,  11.60606061],
     596                                                 [ 11.4,         14.6       ]] )
     597#
    498598
    499599    def test_find_qmax_qmin(self):
     
    542642
    543643        #Extrapolate
     644        quantity.beta = 2.0
     645        quantity.limiter = "minmod_kurganov"
    544646        quantity.extrapolate_second_order()
    545647
  • anuga_work/development/2010-projects/anuga_1d/config.py

    r7860 r7868  
    22"""
    33
    4 epsilon = 1.0e-6
    5 h0 = 1.0e-6
     4epsilon = 1.0e-12
     5h0 = 1.0e-12
    66
    77default_boundary_tag = 'exterior'
  • anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

    r7860 r7868  
    1212
    1313h1 = 10.0
    14 h0 = 0.0
     14h0 = 0.1
    1515
    1616def analytical_sol(C,t):
     
    8181k = 0
    8282
    83 N = 800
     83N = 3200
    8484print "Evaluating domain with %d cells" %N
    85 domain = Domain(*uniform_mesh(N))
     85domain = Domain(*uniform_mesh(N,x_1=L))
    8686   
    8787domain.set_quantity('stage', stage)
     
    9090
    9191domain.set_boundary({'left': Br, 'right' : Br})
    92 domain.order = 2
    93 domain.set_timestepping_method('euler')
     92domain.set_spatial_order(2)
     93domain.set_timestepping_method('rk2')
    9494domain.set_CFL(1.0)
    95 domain.set_limiter("minmod")
     95domain.set_limiter("minmod_kurganov")
    9696#domain.h0=0.0001
    9797
     
    136136
    137137pylab.plot(X.flat,VelocityV.flat)
    138 plot2.set_ylim([-20,10])
     138plot2.set_ylim([-15,15])
    139139
    140140pylab.xlabel('Position')
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

    r7852 r7868  
    9090        self.quantities_to_be_stored = ['stage','xmomentum']
    9191
    92         self.__doc__ = 'shallow_water_domain'
     92        self.__doc__ = 'sww_domain'
    9393
    9494
     
    178178            yield(t)
    179179
    180     def initialise_storage(self):
    181         """Create and initialise self.writer object for storing data.
    182         Also, save x and bed elevation
    183         """
    184 
    185         import data_manager
    186 
    187         #Initialise writer
    188         self.writer = data_manager.get_dataobject(self, mode = 'w')
    189 
    190         #Store vertices and connectivity
    191         self.writer.store_connectivity()
    192 
    193 
    194     def store_timestep(self, name):
    195         """Store named quantity and time.
    196 
    197         Precondition:
    198            self.write has been initialised
    199         """
    200         self.writer.store_timestep(name)
     180
     181
    201182
    202183
     
    253234    u_C   = Velocity.centroid_values
    254235       
    255     #print id(h_C)
    256     #FIXME replace with numpy.where
    257     for i in range(N):
    258         h_C[i] = w_C[i] - z_C[i]                                               
    259         if h_C[i] <= 0.0:
    260             #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
    261             h_C[i] = 1.0e-15
    262             w_C[i] = z_C[i]
    263             uh_C[i] = 0.0
    264                
    265    
    266 ##    for i in range(len(h_C)):
    267 ##      if h_C[i] < epsilon:
    268 ##          u_C[i] = 0.0  #Could have been negative
    269 ##          h_C[i] = 0.0
    270 ##      else:
    271 
    272     u_C[:,]  = uh_C/(h_C + h0/h_C)
     236    #Calculate height (and fix negatives)better be non-negative!
     237    h_C[:] = w_C - z_C
     238    u_C[:]  = uh_C/(h_C + h0/h_C)
    273239       
    274240    for name in [ 'velocity', 'stage' ]:
     
    284250            raise 'Unknown order'
    285251
     252
    286253    stage_V = domain.quantities['stage'].vertex_values                 
    287254    bed_V   = domain.quantities['elevation'].vertex_values     
     
    291258
    292259    h_V[:,:]    = stage_V - bed_V
    293     for i in range(len(h_C)):
    294         for j in range(2):
    295             if h_V[i,j] < 0.0 :
    296                 #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
    297                 dh = h_V[i,j]
    298                 h_V[i,j] = 0.0
    299                 stage_V[i,j] = bed_V[i,j]
    300                 h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
    301                 stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
    302                
     260
     261#    # protect from edge values going negative
     262#    h_V[:,1] = numpy.where(h_V[:,0] < 0.0 , h_V[:,1]-h_V[:,0], h_V[:,1])
     263#    h_V[:,0] = numpy.where(h_V[:,0] < 0.0 , 0.0, h_V[:,0])
     264#
     265#    h_V[:,0] = numpy.where(h_V[:,1] < 0.0 , h_V[:,0]-h_V[:,1], h_V[:,0])
     266#    h_V[:,1] = numpy.where(h_V[:,1] < 0.0 , 0.0, h_V[:,1])
     267#
     268#
     269#    stage_V[:,:] = bed_V + h_V
    303270    xmom_V[:,:] = u_V * h_V
    304271   
  • anuga_work/development/2010-projects/anuga_1d/utilities/util_ext.h

    r7839 r7868  
    3131
    3232
     33double sign(double x) {
     34  //Return sign of a double
     35
     36  if (x>0.0) return 1.0;
     37  else if (x<0.0) return -1.0;
     38  else return 0.0;
     39}
     40
     41
    3342int _gradient(double x0, double y0,
    3443              double x1, double y1,
Note: See TracChangeset for help on using the changeset viewer.