Changeset 6042


Ignore:
Timestamp:
Dec 2, 2008, 9:46:45 AM (16 years ago)
Author:
steve
Message:

Found bug associated with changing max_speed_array to an integer

Location:
anuga_work/development/anuga_1d
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/channel_domain.py

    r6038 r6042  
    9292
    9393
    94     def set_quantities_to_be_stored(self, q):
    95         """Specify which quantities will be stored in the sww file.
    96 
    97         q must be either:
    98           - the name of a quantity
    99           - a list of quantity names
    100           - None
    101 
    102         In the two first cases, the named quantities will be stored at each
    103         yieldstep
    104         (This is in addition to the quantities elevation and friction) 
    105         If q is None, storage will be switched off altogether.
    106         """
    107 
    108 
    109         if q is None:
    110             self.quantities_to_be_stored = []           
    111             self.store = False
    112             return
    113 
    114         if isinstance(q, basestring):
    115             q = [q] # Turn argument into a list
    116 
    117         #Check correcness   
    118         for quantity_name in q:
    119             msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
    120             assert quantity_name in self.conserved_quantities, msg
    121        
    122         self.quantities_to_be_stored = q
    123        
    124 
    12594    def check_integrity(self):
    12695
    … …  
    156125       
    157126
    158 #=============== End of Shallow Water Domain ===============================
     127#=============== End of Channel Domain ===============================
    159128#-----------------------------------
    160129# Compute fluxes interface

  • anuga_work/development/anuga_1d/comp_flux_vel_ext.c

    r5999 r6042  
    145145    double max_speed, normal;
    146146    int k, i, ki, n, m, nm=0;
    147    
     147
     148    //printf("Inside _comp\n");   
    148149   
    149150    for (k=0; k<number_of_elements; k++) {
    150151        flux[0] = 0.0;
    151152        flux[1] = 0.0;
    152        
     153        //printf("k = %d\n",k);
     154
    153155        for (i=0; i<2; i++) {
    154156            ki = k*2+i;
    … …  
    526528                        &velocity)) {
    527529      PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
     530      printf("comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
    528531      return NULL;
    529532  }
    … …  
    560563    number_of_elements = stage_vertex_values -> dimensions[0];
    561564 
     565
     566    //printf("comp_flux_vel 1 N = %d %15.5e %15.5e %15.5e %15.5e %15.5e \n",number_of_elements,cfl,timestep,epsilon,g,h0);
     567
     568
     569
    562570    // Call underlying flux computation routine and update
    563571    // the explicit update arrays
    … …  
    585593                                       number_of_elements,
    586594                                       (double*) max_speed_array -> data);
    587    
     595
     596    //printf("comp_flux_vel 2 \n");   
    588597   
    589598    Py_DECREF(neighbours);
    … …  
    605614    Py_DECREF(max_speed_array);
    606615   
    607    
     616    //printf("comp_flux_vel 3 \n");
     617
    608618    // Return updated flux timestep
    609619    return Py_BuildValue("d", timestep);

  • anuga_work/development/anuga_1d/domain.py

    r5844 r6042  
    431431        self.tagged_elements = tagged_elements
    432432
     433
     434    def set_quantities_to_be_stored(self, q):
     435        """Specify which quantities will be stored in the sww file.
     436
     437        q must be either:
     438          - the name of a quantity
     439          - a list of quantity names
     440          - None
     441
     442        In the two first cases, the named quantities will be stored at each
     443        yieldstep
     444        (This is in addition to the quantities elevation and friction) 
     445        If q is None, storage will be switched off altogether.
     446        """
     447
     448
     449        if q is None:
     450            self.quantities_to_be_stored = []           
     451            self.store = False
     452            return
     453
     454        if isinstance(q, basestring):
     455            q = [q] # Turn argument into a list
     456
     457        #Check correcness   
     458        for quantity_name in q:
     459            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
     460            assert quantity_name in self.conserved_quantities, msg
     461       
     462        self.quantities_to_be_stored = q
     463       
     464
     465
     466
     467
    433468    def get_boundary_tags(self):
    434469        """Return list of available boundary tags
    … …  
    909944            elif self.get_timestepping_method() == 'rk3':
    910945                self.evolve_one_rk3_step(yieldstep,finaltime)               
     946
    911947           
    912948            # Update extrema if necessary (for reporting)
    … …  
    914950
    915951
     952           
    916953            self.yieldtime += self.timestep
    917954            self.number_of_steps += 1
    918955            if self._order_ == 1:
    919956                self.number_of_first_order_steps += 1
     957
    920958
    921959            # Yield results
    … …  
    936974                break
    937975
    938 
    939976            if self.yieldtime >= yieldstep:
    940977                # Yield (intermediate) time and allow inspection of domain
    … …  
    945982
    946983                # Pass control on to outer loop for more specific actions
     984
    947985                yield(self.time)
    948986
    … …  
    953991                self.number_of_steps = 0
    954992                self.number_of_first_order_steps = 0
    955                 self.max_speed_array = 0.0
     993                #self.max_speed_array = 0.0
    956994
    957995
    … …  
    961999        Q^{n+1} = E(h) Q^n
    9621000        """
     1001
    9631002
    9641003        # Compute fluxes across each element edge

  • anuga_work/development/anuga_1d/dry_dam_vel.py

    r5845 r6042  
    6565import time
    6666
    67 finaltime = 10.0
    68 yieldstep = finaltime
     67finaltime = 2.0
     68yieldstep = 0.1
    6969L = 2000.0     # Length of channel (m)
    70 number_of_cells = [800]#,200,500,1000,2000,5000,10000,20000]
    71 h_error = zeros(len(number_of_cells),Float)
    72 uh_error = zeros(len(number_of_cells),Float)
    73 u_error = zeros(len(number_of_cells),Float)
     70
    7471k = 0
    75 for i in range(len(number_of_cells)):
    76     N = int(number_of_cells[i])
    77     print "Evaluating domain with %d cells" %N
    78     cell_len = L/N # Origin = 0.0
    79     points = zeros(N+1,Float)
    80     for j in range(N+1):
    81         points[j] = j*cell_len
     72
     73N = 800
     74print "Evaluating domain with %d cells" %N
     75cell_len = L/N # Origin = 0.0
     76points = zeros(N+1,Float)
     77
     78for j in range(N+1):
     79    points[j] = j*cell_len
    8280       
    83     domain = Domain(points)
     81domain = Domain(points)
    8482   
    85     domain.set_quantity('stage', stage)
    86     domain.set_boundary({'exterior': Reflective_boundary(domain)})
    87     domain.order = 2
    88     domain.set_timestepping_method('rk2')
    89     domain.set_CFL(1.0)
    90     domain.set_limiter("vanleer")
    91     #domain.h0=0.0001
     83domain.set_quantity('stage', stage)
     84domain.set_boundary({'exterior': Reflective_boundary(domain)})
     85domain.order = 2
     86domain.set_timestepping_method('euler')
     87domain.set_CFL(1.0)
     88domain.set_limiter("vanleer")
     89#domain.h0=0.0001
    9290
    93     t0 = time.time()
     91t0 = time.time()
    9492
    95     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    96         domain.write_time()
     93for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     94    domain.write_time()
    9795
    98     N = float(N)
    99     StageC = domain.quantities['stage'].centroid_values
    100     XmomC = domain.quantities['xmomentum'].centroid_values
    101     u_C = domain.quantities['velocity'].centroid_values
    102     C = domain.centroids
    103     h, uh, u = analytical_sol(C,domain.time)
    104     h_error[k] = 1.0/(N)*sum(abs(h-StageC))
    105     uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
    106     u_error[k] = 1.0/(N)*sum(abs(u-u_C))
    107     print "h_error %.10f" %(h_error[k])
    108     print "uh_error %.10f"% (uh_error[k])
    109     print "u_error %.10f"% (u_error[k])
    110     k = k+1
    111     print 'That took %.2f seconds' %(time.time()-t0)
    112     X = domain.vertices
    113     StageQ = domain.quantities['stage'].vertex_values
    114     XmomQ = domain.quantities['xmomentum'].vertex_values
    115     velQ = domain.quantities['velocity'].vertex_values
    116    
    117     h, uh, u = analytical_sol(X.flat,domain.time)
    118     x = X.flat
    119    
    120     from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    121     print 'test 2'
    122     hold(False)
    123     print 'test 3'
    124     plot1 = subplot(211)
    125     print 'test 4'
    126     plot(x,h,x,StageQ.flat)
    127     print 'test 5'
    128     plot1.set_ylim([-1,11])
    129     xlabel('Position')
    130     ylabel('Stage')
    131     legend(('Analytical Solution', 'Numerical Solution'),
    132            'upper right', shadow=True)
    133     plot2 = subplot(212)
    134     plot(x,u,x,velQ.flat)
    135     plot2.set_ylim([-35,35])
    136    
    137     xlabel('Position')
    138     ylabel('Velocity')
    139    
    140     file = "dry_bed_"
    141     file += str(number_of_cells[i])
    142     file += ".eps"
    143     #savefig(file)
    144     show()
    145    
    146 print "Error in height", h_error
    147 print "Error in xmom", uh_error
     96print 'end'
     97
     98

  • anuga_work/development/anuga_1d/shallow_water_vel_domain.py

    r5845 r6042  
    9292
    9393
    94     def set_quantities_to_be_stored(self, q):
    95         """Specify which quantities will be stored in the sww file.
    96 
    97         q must be either:
    98           - the name of a quantity
    99           - a list of quantity names
    100           - None
    101 
    102         In the two first cases, the named quantities will be stored at each
    103         yieldstep
    104         (This is in addition to the quantities elevation and friction) 
    105         If q is None, storage will be switched off altogether.
    106         """
    107 
    108 
    109         if q is None:
    110             self.quantities_to_be_stored = []           
    111             self.store = False
    112             return
    113 
    114         if isinstance(q, basestring):
    115             q = [q] # Turn argument into a list
    116 
    117         #Check correcness   
    118         for quantity_name in q:
    119             msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
    120             assert quantity_name in self.conserved_quantities, msg
    121        
    122         self.quantities_to_be_stored = q
    123        
    12494
    12595    def check_integrity(self):
    … …  
    205175
    206176
     177
    207178#--------------------------------------------------------------------------
    208179def distribute_to_vertices_and_edges_limit_w_u(domain):
    … …  
    255226       
    256227    #print id(h_C)
    257     for i in range(N):
    258         h_C[i] = w_C[i] - z_C[i]                                               
    259         if h_C[i] <= 1.0e-12:
    260             #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
    261             h_C[i] = 0.0
    262             w_C[i] = z_C[i]
    263             #uh_C[i] = 0.0
     228##     for i in range(N):
     229##      h_C[i] = w_C[i] - z_C[i]                                               
     230##         if h_C[i] <= 1.0e-12:
     231##          #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
     232##          h_C[i] = 0.0
     233##          w_C[i] = z_C[i]
     234##          #uh_C[i] = 0.0
    264235           
    265  #           u_C[i] = 0.0
    266  #       else:
    267  #           u_C[i] = uh_C[i]/h_C[i]
     236## #           u_C[i] = 0.0
     237## #       else:
     238## #           u_C[i] = uh_C[i]/h_C[i]
    268239               
    269240    h0 = 1.0e-12   
    270     for i in range(len(h_C)):
     241    for i in range(N):
     242        h_C[i] = w_C[i] - z_C[i]                                               
    271243        if h_C[i] < 1.0e-12:
    272244            u_C[i] = 0.0  #Could have been negative
    273245            h_C[i] = 0.0
     246            w_C[i] = z_C[i]
    274247        else:
    275             u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
    276             #u_C[i]  = uh_C[i]/h_C[i]
     248            #u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])
     249            u_C[i]  = uh_C[i]/h_C[i]
    277250       
    278251    for name in [ 'velocity', 'stage' ]:

  • anuga_work/development/anuga_1d/test_shallow_water_vel_domain.py

    r5845 r6042  
    191191
    192192
    193 ##     def test_evolve_first_order(self):
    194 ##         """
    195 ##         Compare still lake solution for various versions of shallow_water_domain
    196 ##         """
    197        
    198 ##         def slope_square(x):
    199 ##             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
    200 
    201 ##         domain = Domain(self.points2)
    202 ##         domain.set_quantity('stage',10.0)
    203 ##         domain.set_quantity('elevation',slope_square)
    204 ##         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    205 
    206 ##         domain.default_order = 1
    207 ##         domain.set_timestepping_method('euler')
    208        
    209 ##         yieldstep=1.0
    210 ##         finaltime=1.0
    211 
    212 ##         for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    213 ##             pass
    214        
    215 ##         ##   print domain.quantities['stage'].vertex_values
    216 ##         ##   print domain.quantities['elevation'].vertex_values
    217 ##         ##        print domain.quantities['xmomentum'].vertex_values
    218 ##         ## 
    219 ##         ##
    220 ##         ##        print domain.quantities['stage'].centroid_values 
    221 ##         ##   print domain.quantities['elevation'].centroid_values   
    222 ##         ##        print domain.quantities['xmomentum'].centroid_values   
    223 
    224 ##         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
    225 ##         assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
     193    def test_evolve_first_order(self):
     194        """
     195        Compare still lake solution for various versions of shallow_water_domain
     196        """
     197       
     198        def slope_square(x):
     199            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
     200
     201        domain = Domain(self.points2)
     202        domain.set_quantity('stage',10.0)
     203        domain.set_quantity('elevation',slope_square)
     204        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
     205
     206        domain.default_order = 1
     207        domain.set_timestepping_method('euler')
     208       
     209        yieldstep=0.25
     210        finaltime=1.0
     211
     212        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
     213            pass
     214       
     215        ##      print domain.quantities['stage'].vertex_values
     216        ##      print domain.quantities['elevation'].vertex_values
     217        ##        print domain.quantities['xmomentum'].vertex_values
     218        ## 
     219        ##
     220        ##        print domain.quantities['stage'].centroid_values 
     221        ##      print domain.quantities['elevation'].centroid_values   
     222        ##        print domain.quantities['xmomentum'].centroid_values   
     223
     224        #assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
     225        #assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
    226226
    227227

  • anuga_work/development/anuga_1d/util_ext.h

    r5742 r6042  
    261261   
    262262  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
     263
     264  //printf("B = %p\n",(void*)B);
    263265  if (!B) {
    264     PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
     266    printf("util_ext.h: get_consecutive_array could not obtain python object");
     267    printf(" %s\n",name);
     268    fflush(stdout);
     269    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain python object");
    265270    return NULL;
    266271  }     
    … …  
    273278 
    274279  if (!A) {
     280    printf("util_ext.h: get_consecutive_array could not obtain array object");
     281    printf(" %s \n",name);
     282    fflush(stdout);
    275283    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array");
    276284    return NULL;
    277285  }
     286
     287
    278288  return A;
    279289}
Note: See TracChangeset for help on using the changeset viewer.