# Changeset 3424

Ignore:
Timestamp:
Jul 26, 2006, 12:12:43 PM (18 years ago)
Message:

minmod limiter for shallow water

Location:
development/pyvolution-1d
Files:
4 edited

Unmodified
Removed
• ## development/pyvolution-1d/analytic_dam.py

 r3401 def __init__(self, h0 = 5.0, h1 = 10.0, L = 2000.0): from math import sqrt self.h0 = h0 # depth upstream (m) self.h1 = h1 # depth downstream (m) self.L  = L  # length of domain def __call__(self, C,t): from Numeric import zeros,Float from math import sqrt, pi #t  = 0.0     # time (s) h0 = self.h0 h1 = self.h1 L = self.L n = len(C)    # number of cells g  = 9.81    # gravity (m/s^2) c0 = sqrt(g*h0) #left celerity c1 = sqrt(g*h1) #right celerity u = zeros(n,Float) h = zeros(n,Float) uh = zeros(n,Float) x = C-L/2 #x = zeros(n,Float) #for i in range(n): #    x[i] = C[i]-1000.0 # Upstream and downstream boundary conditions are set to the intial water # depth for all time. # Calculate Shock Speed #h2 = 7.2692044 #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2)) #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0))) zmin=-100.0 zmax=101.0 if( abs(z) > 99.0): print 'no convergence' self.u2 = u2 self.c0 = c0 self.c1 = c1 self.c2 = c2 self.g = g self.z = z def __call__(self, C,t): from Numeric import zeros,Float from math import sqrt #t  = 0.0     # time (s) h0 = self.h0 h1 = self.h1 L = self.L n = len(C)    # number of cells u2 = self.u2 c0 = self.c0 c1 = self.c1 c2 = self.c2 g = self.g z = self.z u = zeros(n,Float) h = zeros(n,Float) uh = zeros(n,Float) x = C-L/2.0 #x = zeros(n,Float) #for i in range(n): #    x[i] = C[i]-1000.0 # Upstream and downstream boundary conditions are set to the intial water # depth for all time. # Calculate Shock Speed #h2 = 7.2692044 #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2)) #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0))) h2=h0/(1.0-u2/z)
• ## development/pyvolution-1d/dam.py

 r3370 from Numeric import allclose, array, zeros, ones, Float, take, sqrt from config import g, epsilon from analytic_dam import AnalyticDam def analytical_sol(C,t): h0 = 0.1 h1 = 10.0 analytical_sol = AnalyticDam(h0, h1) ## def analytical_sol(C,t): #t  = 0.0     # time (s) g  = 9.81    # gravity (m/s^2) h1 = 10.0    # depth upstream (m) h0 = 5.0     # depth downstream (m) L = 2000.0   # length of stream/domain (m) n = len(C)    # number of cells ##     #t  = 0.0     # time (s) ##     g  = 9.81    # gravity (m/s^2) ##     h1 = 10.0    # depth upstream (m) ##     h0 = 5.0     # depth downstream (m) ##     L = 2000.0   # length of stream/domain (m) ##     n = len(C)    # number of cells u = zeros(n,Float) h = zeros(n,Float) x = C-L/2 #x = zeros(n,Float) #for i in range(n): #    x[i] = C[i]-1000.0 ##     u = zeros(n,Float) ##     h = zeros(n,Float) ##     x = C-L/2 ##     #x = zeros(n,Float) ##     #for i in range(n): ##     #    x[i] = C[i]-1000.0 # Upstream and downstream boundary conditions are set to the intial water # depth for all time. ##     # Upstream and downstream boundary conditions are set to the intial water ##     # depth for all time. # Calculate Shock Speed h2 = 7.2692044 S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2)) u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0))) ##     # Calculate Shock Speed ##     h2 = 7.2692044 ##     S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2)) ##     u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0))) #t=50 #x = (-L/2:L/2) for i in range(n): # Calculate Analytical Solution at time t > 0 u3 = 2/3*(sqrt(g*h1)+x[i]/t) h3 = 4/(9*g)*(sqrt(g*h1)-x[i]/(2*t))*(sqrt(g*h1)-x[i]/(2*t)) ##     #t=50 ##     #x = (-L/2:L/2) ##     for i in range(n): ##         # Calculate Analytical Solution at time t > 0 ##         u3 = 2/3*(sqrt(g*h1)+x[i]/t) ##         h3 = 4/(9*g)*(sqrt(g*h1)-x[i]/(2*t))*(sqrt(g*h1)-x[i]/(2*t)) if ( x[i] <= -t*sqrt(g*h1) ): u[i] = 0.0 h[i] = h1 elif ( x[i] <= t*(u2-sqrt(g*h2)) ): u[i] = u3 h[i] = h3 elif ( x[i] < t*S2 ): u[i] = u2 h[i] = h2 else: u[i] = 0.0 h[i] = h0 ##         if ( x[i] <= -t*sqrt(g*h1) ): ##             u[i] = 0.0 ##             h[i] = h1 ##         elif ( x[i] <= t*(u2-sqrt(g*h2)) ): ##             u[i] = u3 ##             h[i] = h3 ##         elif ( x[i] < t*S2 ): ##             u[i] = u2 ##             h[i] = h2 ##         else: ##             u[i] = 0.0 ##             h[i] = h0 return h , u*h ##     return h , u*h def newLinePlot(title='Simple Plot'): import Gnuplot g = Gnuplot.Gnuplot(persist=1) g = Gnuplot.Gnuplot() g.title(title) g('set data style linespoints') L = 2000.0     # Length of channel (m) N = 100        # Number of compuational cells N = 400        # Number of compuational cells cell_len = L/N # Origin = 0.0 for i in range(len(x)): if x[i]<=1000.0: y[i] = 10.0 y[i] = h1 else: y[i] = 5.0 y[i] = h0 return y domain.set_quantity('stage', stage) domain.order = 2 domain.default_order = 2 domain.cfl = 0.8 #domain.beta = 0.0 print "domain.order", domain.order print "domain.order", domain.default_order if (debug == True): plot1 = newLinePlot("Stage") plot2 = newLinePlot("Momentum") plot3 = newLinePlot("Velocity") import time yieldstep = 1.0 finaltime = 50.0 print domain.neighbour_vertices for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): domain.write_time() linePlot(plot1,X,StageQ,X,y) MomentumQ = domain.quantities['xmomentum'].vertex_values Velocity = MomentumQ/StageQ linePlot(plot2,X,MomentumQ,X,my) linePlot(plot3,X,Velocity,X,my/y) #raw_input('press_return') #pass print C raw_input('press return') #f = file('test_solution_I.out', 'w') #for i in range(N):