Changeset 5827


Ignore:
Timestamp:
Oct 9, 2008, 12:54:08 PM (17 years ago)
Author:
sudi
Message:

Adding new versions of shallow_water_domain

Location:
anuga_work/development/anuga_1d
Files:
3 added
8 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/analytic_dam_sudi.py

    r5535 r5827  
    2626            else:
    2727                zmax=z
    28 
     28            #print 'func=',func
    2929        if( abs(z) > 99.0):
    3030            print 'no convergence'

  • anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

    r5725 r5827  
    2626//Innermost flux function (using w=z+h)
    2727int _flux_function_wellbalanced(double *q_left, double *q_right,
    28         double z_left, double z_right,
    2928                double normals, double g, double epsilon,
    3029                double *edgeflux, double *max_speed) {
    3130               
    3231                int i;
     32        double z_left, z_right;
    3333                double ql[2], qr[2], flux_left[2], flux_right[2];
    3434                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
    … …  
    3636                double s_max, s_min, denom;
    3737               
     38                w_left = q_left[0];
     39                uh_left = q_left[1];
     40                uh_left = uh_left*normals;
     41        z_left = q_left[2];
     42       
     43                w_right = q_right[0];
     44                uh_right = q_right[1];
     45                uh_right = uh_right*normals;
     46        z_right = q_right[2];
     47         
    3848                if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
    39                 ql[0] = q_left[0];
    40                 ql[1] = q_left[1];
    41                 ql[1] = ql[1]*normals;
    42        
    43                 qr[0] = q_right[0];
    44                 qr[1] = q_right[1];
    45                 qr[1] = qr[1]*normals;
    46          
     49       
    4750                //z = (z_left+z_right)/2.0;
    4851                z_star = max(z_left, z_right);                                                                                                          //equation(7) 
    4952                                   
    5053                // Compute speeds in x-direction
    51                 w_left = ql[0];
    5254                h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1].
    5355                h_left_star = max(0, w_left-z_star);                                                                                                            //(8)
    54                 uh_left = ql[1];
     56   
    5557                if (h_left_star < epsilon) {
    5658                        u_left = 0.0; h_left_star = 0.0;
    … …  
    6062                }
    6163       
    62                 w_right = qr[0];
    6364                h_right = w_right-z_right;                                                                                                                                                                                                     
    6465                h_right_star = max(0, w_right-z_star);                                                                                                          //(9)
    65                 uh_right = qr[1];
     66
    6667                if (h_right_star < epsilon) {
    6768                        u_right = 0.0; h_right_star = 0.0;
    … …  
    134135                double* max_speed_array) {
    135136               
    136                 double flux[2], ql[2], qr[2], edgeflux[2];
     137                double flux[2], ql[3], qr[3], edgeflux[2];
    137138                double zl, zr, max_speed, normal;
    138139                int k, i, ki, n, m, nm=0;
    … …  
    147148                                ql[0] = stage_edge_values[ki];
    148149                                ql[1] = xmom_edge_values[ki];
    149                                 zl = bed_edge_values[ki];
    150                                
     150                                ql[2] = bed_edge_values[ki];
     151                                //ql[3] = velocity_edge_values[ki];
     152                //ql[4] = height_edge_values[ki];
     153               
    151154                                n = neighbours[ki];
    152155                                if (n<0) {
    … …  
    154157                                        qr[0] = stage_boundary_values[m];
    155158                                        qr[1] = xmom_boundary_values[m];
    156                                         zr = zl;
     159                                        qr[2] = ql[2];
     160                   
    157161                                } else {
    158162                                        m = neighbour_vertices[ki];
    … …  
    160164                                        qr[0] = stage_edge_values[nm];
    161165                                        qr[1] = xmom_edge_values[nm];
    162                                         zr = bed_edge_values[nm];                               
     166                                        qr[2] = bed_edge_values[nm];                           
    163167                                }
    164168                               
    165169                                normal = normals[ki];
    166                                 _flux_function_wellbalanced(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
     170                                _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
    167171                                flux[0] -= edgeflux[0];
    168172                                flux[1] -= edgeflux[1];
    … …  
    186190                        max_speed_array[k]=max_speed;
    187191                }
    188                 printf("timestep = %f \n",timestep);
     192                //printf("timestep = %f \n",timestep);
    189193                return timestep;       
    190194        }

  • anuga_work/development/anuga_1d/config.py

    r5741 r5827  
    33
    44epsilon = 1.0e-12
    5 h0 = 1.0e-6
     5h0 = 1.0e-12
    66
    77default_boundary_tag = 'exterior'
    … …  
    3939max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order
    4040
    41 manning = 0.3  #Manning's friction coefficient
     41manning = 0.0  #Manning's friction coefficient
    4242g = 9.8       #Gravity
    4343#g(phi) = 9780313 * (1 + 0.0053024 sin(phi)**2 - 0.000 0059 sin(2*phi)**2) micro m/s**2, where phi is the latitude
    … …  
    109109
    110110#Specific to shallow water W.E.
    111 minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0
    112 maximum_allowed_speed = 100.0
     111minimum_allowed_height = 1.0e-6 #Water depth below which it is considered to be 0
     112maximum_allowed_speed = 1000.0

  • anuga_work/development/anuga_1d/dam_h_elevation.py

    r5743 r5827  
    3131
    3232L=2000.0
    33 N=200
     33N=50
    3434
    3535cell_len=L/N
    … …  
    123123import time
    124124t0=time.time()
    125 yieldstep=45.0 #30.0
    126 finaltime=45.0 #20.0
     125yieldstep=10.0 #30.0
     126finaltime=10.0 #20.0
    127127print "integral", domain.quantities['stage'].get_integral()
    128128for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    … …  
    152152
    153153        plot(X,ElevationQ,X,HeightQ)
    154         #plot1.set_ylim([9.99,10.01])
     154        plot1.set_ylim([-1.0,11.0])
    155155        xlabel('Position')
    156156        ylabel('Stage')

  • anuga_work/development/anuga_1d/dry_dam_sudi.py

    r5587 r5827  
    11import os
    22from math import sqrt, pi
    3 from shallow_water_domain import *
     3from shallow_water_domain_suggestion1 import *
    44from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    55from config import g, epsilon
    … …  
    88   
    99    #t  = 0.0     # time (s)
    10     g  = 9.81    # gravity (m/s^2)
     10    # gravity (m/s^2)
    1111    h1 = 10.0    # depth upstream (m)
    1212    h0 = 0.0     # depth downstream (m)
    … …  
    4343            h[i] = h0
    4444
    45     return h , u*h
     45    return h , u*h, u
    4646
    4747#def newLinePlot(title='Simple Plot'):
    … …  
    6565print "TEST 1D-SOLUTION III -- DRY BED"
    6666
    67 L = 2000.0     # Length of channel (m)
    68 N = 800        # Number of computational cells
    69 cell_len = L/N # Origin = 0.0
    70 
    71 points = zeros(N+1,Float)
    72 for i in range(N+1):
    73     points[i] = i*cell_len
    74 
    75 domain = Domain(points)
    76 
    7767def stage(x):
    7868    h1 = 10.0
    … …  
    9282
    9383finaltime = 20.0
    94 yieldstep = 1.0
     84yieldstep = 20.0
    9585L = 2000.0     # Length of channel (m)
    96 number_of_cells = [200]#,200,500,1000,2000,5000,10000,20000]
     86number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
    9787h_error = zeros(len(number_of_cells),Float)
    9888uh_error = zeros(len(number_of_cells),Float)
    … …  
    110100    domain.set_quantity('stage', stage)
    111101    domain.set_boundary({'exterior': Reflective_boundary(domain)})
    112     domain.default_order = 2
    113     domain.default_time_order = 2
     102    domain.order = 2
     103    domain.set_timestepping_method('rk2')
    114104    domain.cfl = 1.0
    115     domain.limiter = "minmod"
     105    domain.limiter = "vanleer"
    116106
    117107    t0 = time.time()
    … …  
    124114    XmomC = domain.quantities['xmomentum'].centroid_values
    125115    C = domain.centroids
    126     h, uh = analytical_sol(C,domain.time)
     116    h, uh, u = analytical_sol(C,domain.time)
    127117    h_error[k] = 1.0/(N)*sum(abs(h-StageC))
    128118    uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
    … …  
    134124    StageQ = domain.quantities['stage'].vertex_values
    135125    XmomQ = domain.quantities['xmomentum'].vertex_values
    136     h, uh = analytical_sol(X.flat,domain.time)
     126    velQ = domain.quantities['velocity'].vertex_values
     127   
     128    h, uh, u = analytical_sol(X.flat,domain.time)
    137129    x = X.flat
    138130   
    … …  
    151143           'upper right', shadow=True)
    152144    plot2 = subplot(212)
    153     plot(x,uh,x,XmomQ.flat)
     145    plot(x,u,x,velQ.flat)
    154146    plot2.set_ylim([-35,35])
    155147   
    156148    xlabel('Position')
    157     ylabel('Xmomentum')
     149    ylabel('Velocity')
    158150   
    159151    file = "dry_bed_"

  • anuga_work/development/anuga_1d/parabolic_cannal.py

    r5535 r5827  
    11import os
    22from math import sqrt, pi, sin, cos
    3 from shallow_water_1d import *
     3from shallow_water_domain import *
    44from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    55from config import g, epsilon
    6 
     6"""
    77def newLinePlot(title='Simple Plot'):
    88    import Gnuplot
    … …  
    2020    plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1")
    2121    g.plot(plot1,plot2,plot3)
    22 
     22"""
    2323
    2424def analytic_cannal(C,t):
    … …  
    5454
    5555L_x = 2500.0     # Length of channel (m)
    56 N = 8         # Number of compuational cells
     56N = 100         # Number of compuational cells
    5757cell_len = 4*L_x/N   # Origin = 0.0
    5858
    … …  
    6363domain = Domain(points)
    6464
    65 domain.order = 2 #make this unnecessary
    66 domain.default_order = 2
    67 domain.default_time_order = 2
     65domain.order = 2
     66domain.set_timestepping_method('rk2')
    6867domain.cfl = 1.0
    6968domain.beta = 1.0
    … …  
    7271#domain.limiter = "vanalbada"
    7372#domain.limiter = "superbee"
    74 domain.limiter = "steve_minmod"
     73#domain.limiter = "steve_minmod"
     74domain.limiter = "pyvolution"
    7575
    7676def stage(x):
    … …  
    113113domain.set_quantity('elevation',elevation)
    114114#domain.set_quantity('height',height)
    115 
    116115domain.set_boundary({'exterior': Reflective_boundary(domain)})
    117116 
    118117import time
    119118t0 = time.time()
    120 finaltime = 1122.0*0.75
    121 yieldstep = finaltime
    122 yieldstep = 10.0
    123 plot1 = newLinePlot("Stage")
    124 plot2 = newLinePlot("Xmomentum")
     119#finaltime = 1122.0*0.75
     120yieldstep = finaltime = 10.0
     121#yieldstep = 10.0
     122#plot1 = newLinePlot("Stage")
     123#plot2 = newLinePlot("Xmomentum")
    125124C = domain.centroids
    126125X = domain.vertices
    127126StageQ = domain.quantities['stage'].vertex_values
    128127XmomQ = domain.quantities['xmomentum'].vertex_values
     128#Bed = domain.quantities['elevation'].vertex_values
    129129import time
    130130t0 = time.time()
    … …  
    133133    domain.write_time()
    134134    u,h,z,z_b = analytic_cannal(X.flat,t)
    135     linePlot(plot1,X,z,X,z_b,X,StageQ)
    136     linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
     135    #linePlot(plot1,X,z,X,z_b,X,StageQ)
     136    #linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
    137137    HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
    138138    u,hc,z,z_b = analytic_cannal(C,t)
    … …  
    142142    error = 1.0/(N)*sum(abs(hc-HeightQ))
    143143    print 'Error measured at centroids', error
    144     #raw_input()         
    145 
    146 #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    147 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
    148 #import Gnuplot 
    149 X = domain.vertices
    150 u,h,z,z_b = analytic_cannal(X.flat,domain.time)
    151 StageQ = domain.quantities['stage'].vertex_values
    152 XmomQ = domain.quantities['xmomentum'].vertex_values
    153 hold(False)
    154 plot1 = subplot(211)
    155 plot(X,h,X,StageQ)
    156 #plot1.set_ylim([4,11])
    157 #title('Free Surface Elevation of a Dry Dam-Break')
    158 xlabel('Position')
    159 ylabel('Stage')
    160 legend(('Analytical Solution', 'Numerical Solution'),
    161        'upper right', shadow=True)
    162 plot2 = subplot(212)
    163 plot(X,u*h,X,XmomQ)
    164 #plot2.set_ylim([-1,25])
    165 #title('Xmomentum Profile of a Dry Dam-Break')
    166 xlabel('Position')
    167 ylabel('Xmomentum')
    168 show()
     144    # #raw_input()         
     145    # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
     146    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
     147    # #import Gnuplot 
     148    X = domain.vertices
     149    u,h,z,z_b = analytic_cannal(X.flat,domain.time)
     150    StageQ = domain.quantities['stage'].vertex_values
     151    XmomQ = domain.quantities['xmomentum'].vertex_values
     152    hold(False)
     153    plot1 = subplot(211)
     154    plot(X,h,X,StageQ,X,z_b)
     155    #plot1.set_ylim([4,11])
     156    #title('?????????')
     157    xlabel('Position')
     158    ylabel('Stage')
     159    legend(('Analytical Solution', 'Numerical Solution', 'Bed'),
     160           'upper right', shadow=True)
     161    plot2 = subplot(212)
     162    plot(X,u*h,X,XmomQ)
     163    #plot2.set_ylim([-1,25])
     164    #title('?????????????????????????')
     165    xlabel('Position')
     166    ylabel('Xmomentum')
     167    show()

  • anuga_work/development/anuga_1d/quantity.py

    r5743 r5827  
    764764            elif n1 < 0:
    765765                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
    766             elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
    767                 phi = 0.0
     766            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
     767            #    phi = 0.0
    768768            else:
    769769                if limiter == "minmod":

  • anuga_work/development/anuga_1d/steady_flow.py

    r5535 r5827  
    11import os
    2 import Gnuplot
     2#import Gnuplot
    33from math import sqrt, pi
    4 from shallow_water_1d import *
     4from shallow_water_domain import *
    55from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    66from config import g, epsilon
    … …  
    2121h_1 = 0.5
    2222u_1 = 0.6
    23 g = 9.81
     23
    2424
    2525#h_c = (q**2/g)**(1/3)
    … …  
    103103
    104104L = 50.0
    105 N = 50 #800
     105N = 400
    106106cell_len = L/N
    107107points = zeros(N+1,Float)
    … …  
    134134wc, uc, hc = analytical_sol(C)
    135135
    136 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
    137 hold(False)
    138 rc('text', usetex=True)
     136#from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
     137#hold(False)
     138#rc('text', usetex=True)
    139139h_error = zeros(1,Float)
    140140uh_error = zeros(1,Float)
    141141k = 0
    142 yieldstep = 25.0
    143 finaltime = 25.0
     142yieldstep = 1.0
     143finaltime = 1.0
    144144domain.limiter = "vanleer"
    145 #domain.limiter = "steve_minmod"
    146 domain.default_order = 2
    147 domain.default_time_order = 2#check order is not working
     145domain.order = 2
     146domain.set_timestepping_method('rk2')
     147domain.cfl = 1.0
     148print 'THE DOMAIN LIMITER is', domain.limiter
     149import time
     150t0=time.time()
    148151for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    149152    domain.write_time()
    150 print 'Test 1'   
    151 plot1 = subplot(211)
    152 print 'Test 2'
    153 plot(X,w,X,StageQ,X,ElevationQ)
    154 print 'Test 3'
    155 #xlabel('Position')
    156 #ylabel('Stage (m)')
    157 #legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
    158 #       'upper right', shadow=True)
    159 plot1.set_ylim([-0.1,1.0])
    160 print 'Test 4'
    161 plot2 = subplot(212)
    162 print 'Test 5'
    163 plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
    164 print 'Test 6'
    165 #xlabel('x (m)')
    166 #ylabel(r'X-momentum ($m^2/s$)')
    167 h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
    168 uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
    169 print "h_error %.10f" %(h_error[k])
    170 print "uh_error %.10f"% (uh_error[k])
    171 print "h_max %.10f"%max(wc-StageC)
    172 print "uh max  %.10f"%max(uc*hc-XmomC)
    173 filename = "steady_flow_"
    174 filename += domain.limiter
    175 filename += str(400)
    176 filename += ".ps"
    177 #savefig(filename)
    178 show()
    179 
     153    print "integral", domain.quantities['stage'].get_integral()
     154    if t>0.0:
     155        print 'Test 1'
     156        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
     157        #ion()
     158        hold(False)
     159        #rc('text', usetex=True)
     160        clf()
     161       
     162        plot1 = subplot(211)
     163        print 'Test 2'
     164        plot(X,w,X,StageQ,X,ElevationQ)
     165        print 'Test 3'
     166        xlabel('Position')
     167        ylabel('Stage (m)')
     168        legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
     169               'upper right', shadow=True)
     170        plot1.set_ylim([-0.1,1.0])
     171        print 'Test 4'
     172        plot2 = subplot(212)
     173        print 'Test 5'
     174        plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
     175        print 'Test 6'
     176        xlabel('x (m)')
     177        ylabel(r'X-momentum (m^2/s)')
     178        plot2.set_ylim([0.297,0.301])
     179        h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
     180        uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
     181        print "h_error %.10f" %(h_error[k])
     182        print "uh_error %.10f"% (uh_error[k])
     183        print "h_max %.10f"%max(wc-StageC)
     184        print "uh max  %.10f"%max(uc*hc-XmomC)
     185        #filename = "steady_flow_"
     186        #filename += domain.limiter
     187        #filename += str(400)
     188        #filename += ".ps"
     189        #savefig(filename)
     190        show()
     191print 'That took %.2f seconds'%(time.time()-t0)
     192raw_input("Press return key!")
Note: See TracChangeset for help on using the changeset viewer.