# Changeset 5827

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

Adding new versions of shallow_water_domain

Location:
anuga_work/development/anuga_1d
Files:
8 edited

Unmodified
Removed
• ## anuga_work/development/anuga_1d/analytic_dam_sudi.py

 r5535 else: zmax=z #print 'func=',func if( abs(z) > 99.0): print 'no convergence'
• ## anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

 r5725 //Innermost flux function (using w=z+h) int _flux_function_wellbalanced(double *q_left, double *q_right, double z_left, double z_right, double normals, double g, double epsilon, double *edgeflux, double *max_speed) { int i; double z_left, z_right; double ql[2], qr[2], flux_left[2], flux_right[2]; double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left; double s_max, s_min, denom; w_left = q_left[0]; uh_left = q_left[1]; uh_left = uh_left*normals; z_left = q_left[2]; w_right = q_right[0]; uh_right = q_right[1]; uh_right = uh_right*normals; z_right = q_right[2]; if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right); ql[0] = q_left[0]; ql[1] = q_left[1]; ql[1] = ql[1]*normals; qr[0] = q_right[0]; qr[1] = q_right[1]; qr[1] = qr[1]*normals; //z = (z_left+z_right)/2.0; z_star = max(z_left, z_right);                                                                                                          //equation(7) // Compute speeds in x-direction w_left = ql[0]; h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1]. h_left_star = max(0, w_left-z_star);                                                                                                            //(8) uh_left = ql[1]; if (h_left_star < epsilon) { u_left = 0.0; h_left_star = 0.0; } w_right = qr[0]; h_right = w_right-z_right; h_right_star = max(0, w_right-z_star);                                                                                                          //(9) uh_right = qr[1]; if (h_right_star < epsilon) { u_right = 0.0; h_right_star = 0.0; double* max_speed_array) { double flux[2], ql[2], qr[2], edgeflux[2]; double flux[2], ql[3], qr[3], edgeflux[2]; double zl, zr, max_speed, normal; int k, i, ki, n, m, nm=0; ql[0] = stage_edge_values[ki]; ql[1] = xmom_edge_values[ki]; zl = bed_edge_values[ki]; ql[2] = bed_edge_values[ki]; //ql[3] = velocity_edge_values[ki]; //ql[4] = height_edge_values[ki]; n = neighbours[ki]; if (n<0) { qr[0] = stage_boundary_values[m]; qr[1] = xmom_boundary_values[m]; zr = zl; qr[2] = ql[2]; } else { m = neighbour_vertices[ki]; qr[0] = stage_edge_values[nm]; qr[1] = xmom_edge_values[nm]; zr = bed_edge_values[nm]; qr[2] = bed_edge_values[nm]; } normal = normals[ki]; _flux_function_wellbalanced(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed); _flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed); flux[0] -= edgeflux[0]; flux[1] -= edgeflux[1]; max_speed_array[k]=max_speed; } printf("timestep = %f \n",timestep); //printf("timestep = %f \n",timestep); return timestep; }
• ## anuga_work/development/anuga_1d/config.py

 r5741 epsilon = 1.0e-12 h0 = 1.0e-6 h0 = 1.0e-12 default_boundary_tag = 'exterior' max_smallsteps = 50  #Max number of degenerate steps allowed b4 trying first order manning = 0.3  #Manning's friction coefficient manning = 0.0  #Manning's friction coefficient g = 9.8       #Gravity #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 #Specific to shallow water W.E. minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0 maximum_allowed_speed = 100.0 minimum_allowed_height = 1.0e-6 #Water depth below which it is considered to be 0 maximum_allowed_speed = 1000.0
• ## anuga_work/development/anuga_1d/dam_h_elevation.py

 r5743 L=2000.0 N=200 N=50 cell_len=L/N import time t0=time.time() yieldstep=45.0 #30.0 finaltime=45.0 #20.0 yieldstep=10.0 #30.0 finaltime=10.0 #20.0 print "integral", domain.quantities['stage'].get_integral() for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): plot(X,ElevationQ,X,HeightQ) #plot1.set_ylim([9.99,10.01]) plot1.set_ylim([-1.0,11.0]) xlabel('Position') ylabel('Stage')
• ## anuga_work/development/anuga_1d/dry_dam_sudi.py

 r5587 import os from math import sqrt, pi from shallow_water_domain import * from shallow_water_domain_suggestion1 import * from Numeric import allclose, array, zeros, ones, Float, take, sqrt from config import g, epsilon #t  = 0.0     # time (s) g  = 9.81    # gravity (m/s^2) # gravity (m/s^2) h1 = 10.0    # depth upstream (m) h0 = 0.0     # depth downstream (m) h[i] = h0 return h , u*h return h , u*h, u #def newLinePlot(title='Simple Plot'): print "TEST 1D-SOLUTION III -- DRY BED" L = 2000.0     # Length of channel (m) N = 800        # Number of computational cells cell_len = L/N # Origin = 0.0 points = zeros(N+1,Float) for i in range(N+1): points[i] = i*cell_len domain = Domain(points) def stage(x): h1 = 10.0 finaltime = 20.0 yieldstep = 1.0 yieldstep = 20.0 L = 2000.0     # Length of channel (m) number_of_cells = [200]#,200,500,1000,2000,5000,10000,20000] number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000] h_error = zeros(len(number_of_cells),Float) uh_error = zeros(len(number_of_cells),Float) domain.set_quantity('stage', stage) domain.set_boundary({'exterior': Reflective_boundary(domain)}) domain.default_order = 2 domain.default_time_order = 2 domain.order = 2 domain.set_timestepping_method('rk2') domain.cfl = 1.0 domain.limiter = "minmod" domain.limiter = "vanleer" t0 = time.time() XmomC = domain.quantities['xmomentum'].centroid_values C = domain.centroids h, uh = analytical_sol(C,domain.time) h, uh, u = analytical_sol(C,domain.time) h_error[k] = 1.0/(N)*sum(abs(h-StageC)) uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC)) StageQ = domain.quantities['stage'].vertex_values XmomQ = domain.quantities['xmomentum'].vertex_values h, uh = analytical_sol(X.flat,domain.time) velQ = domain.quantities['velocity'].vertex_values h, uh, u = analytical_sol(X.flat,domain.time) x = X.flat 'upper right', shadow=True) plot2 = subplot(212) plot(x,uh,x,XmomQ.flat) plot(x,u,x,velQ.flat) plot2.set_ylim([-35,35]) xlabel('Position') ylabel('Xmomentum') ylabel('Velocity') file = "dry_bed_"
• ## anuga_work/development/anuga_1d/parabolic_cannal.py

 r5535 import os from math import sqrt, pi, sin, cos from shallow_water_1d import * from shallow_water_domain import * from Numeric import allclose, array, zeros, ones, Float, take, sqrt from config import g, epsilon """ def newLinePlot(title='Simple Plot'): import Gnuplot plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1") g.plot(plot1,plot2,plot3) """ def analytic_cannal(C,t): L_x = 2500.0     # Length of channel (m) N = 8         # Number of compuational cells N = 100         # Number of compuational cells cell_len = 4*L_x/N   # Origin = 0.0 domain = Domain(points) domain.order = 2 #make this unnecessary domain.default_order = 2 domain.default_time_order = 2 domain.order = 2 domain.set_timestepping_method('rk2') domain.cfl = 1.0 domain.beta = 1.0 #domain.limiter = "vanalbada" #domain.limiter = "superbee" domain.limiter = "steve_minmod" #domain.limiter = "steve_minmod" domain.limiter = "pyvolution" def stage(x): domain.set_quantity('elevation',elevation) #domain.set_quantity('height',height) domain.set_boundary({'exterior': Reflective_boundary(domain)}) import time t0 = time.time() finaltime = 1122.0*0.75 yieldstep = finaltime yieldstep = 10.0 plot1 = newLinePlot("Stage") plot2 = newLinePlot("Xmomentum") #finaltime = 1122.0*0.75 yieldstep = finaltime = 10.0 #yieldstep = 10.0 #plot1 = newLinePlot("Stage") #plot2 = newLinePlot("Xmomentum") C = domain.centroids X = domain.vertices StageQ = domain.quantities['stage'].vertex_values XmomQ = domain.quantities['xmomentum'].vertex_values #Bed = domain.quantities['elevation'].vertex_values import time t0 = time.time() domain.write_time() u,h,z,z_b = analytic_cannal(X.flat,t) linePlot(plot1,X,z,X,z_b,X,StageQ) linePlot(plot2,X,u*h,X,u*h,X,XmomQ) #linePlot(plot1,X,z,X,z_b,X,StageQ) #linePlot(plot2,X,u*h,X,u*h,X,XmomQ) HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values u,hc,z,z_b = analytic_cannal(C,t) error = 1.0/(N)*sum(abs(hc-HeightQ)) print 'Error measured at centroids', error #raw_input() #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot #import Gnuplot X = domain.vertices u,h,z,z_b = analytic_cannal(X.flat,domain.time) StageQ = domain.quantities['stage'].vertex_values XmomQ = domain.quantities['xmomentum'].vertex_values hold(False) plot1 = subplot(211) plot(X,h,X,StageQ) #plot1.set_ylim([4,11]) #title('Free Surface Elevation of a Dry Dam-Break') xlabel('Position') ylabel('Stage') legend(('Analytical Solution', 'Numerical Solution'), 'upper right', shadow=True) plot2 = subplot(212) plot(X,u*h,X,XmomQ) #plot2.set_ylim([-1,25]) #title('Xmomentum Profile of a Dry Dam-Break') xlabel('Position') ylabel('Xmomentum') show() # #raw_input() # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot # #import Gnuplot X = domain.vertices u,h,z,z_b = analytic_cannal(X.flat,domain.time) StageQ = domain.quantities['stage'].vertex_values XmomQ = domain.quantities['xmomentum'].vertex_values hold(False) plot1 = subplot(211) plot(X,h,X,StageQ,X,z_b) #plot1.set_ylim([4,11]) #title('?????????') xlabel('Position') ylabel('Stage') legend(('Analytical Solution', 'Numerical Solution', 'Bed'), 'upper right', shadow=True) plot2 = subplot(212) plot(X,u*h,X,XmomQ) #plot2.set_ylim([-1,25]) #title('?????????????????????????') xlabel('Position') ylabel('Xmomentum') show()
• ## anuga_work/development/anuga_1d/quantity.py

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