Godunov

flood.py

 from sympy import *
+from matplotlib import pyplot as plt
+from scipy import integrate as sp_int
+import numpy as np
 
 # def wrt_w():
 #     A_sym = Symbol("A")
@@ -26,41 +29,105 @@ from sympy import *
 #     print(simplify(l))
 #     print(simplify(l.subs(w_sym, w)))
 
-def main():
+def plot_expr(f, x_sym, lim, N):
+    F = lambdify(x_sym, f)
+    x = []
+    y = []
+    for i in range(N):
+        xi = i * (lim[1]-lim[0]) / (N-1) + lim[0]
+        x.append(xi)
+        y.append(F(xi))
+    plt.plot(x, y)
+    plt.show()
+
+def godunov_step(Q, delta, c, c_inflow):
+    N = len(c)
+    dc_dt = [-(Q(c[0]) - Q(c_inflow))/delta]
+    for i in range(1, N):
+        dc_dt.append(-(Q(c[i]) - Q(c[i-1]))/delta)
+    return dc_dt
+
+
+    
+def derive_wetted_perim(w, y, A_sym, debug=False):
     # Define symbols and functions
-    A_sym = Symbol("A", positive=True, real=True)
     l = Function("l", positive=True, real=True)(A_sym)
 
     h_sym = h = Symbol("h", positive=True, real=True)
-    y = Symbol("y", positive=True, real=True)
-    
-    # Width of channel wrt vertical coordinate (y)
-    a = Symbol("a", positive=True, real=True)
-    w = sqrt(y/a)
 
     # Calculate length by integrating segment lengths between 0 and h.
     # Add the width at the base to the total length for the case of a flat bottom.
     dx_dy = diff(w/2, y) # x = +- w/2
     dl_dy = 2 * sqrt(dx_dy**2 + 1) # 2* for both left and right side of channel
     l = integrate(dl_dy, (y, 0, h)) + w.subs(y, 0)
-    print("l(h) =", l)
+    if debug:
+        print("l(h) =", l)
 
     # Calculate channel area by integrating width wrt vertical coordinate.
     A = simplify(integrate(w, (y, 0, h)))
-    print("A(h) =", A)
+    if debug:
+        print("A(h) =", A)
 
     # Get h in terms of A (there may be multiple solutions)
     h_of_A_arr = solve(A_sym - A, h)
 
     # For each solution...
     for h_of_A in h_of_A_arr:
-        print("================")
-        print("h(A) =")
-        pprint(h_of_A)
+        if debug:
+            print("================")
+            print("h(A) =")
+            pprint(h_of_A)
 
         # Substitute the expression for h(A) into l(h) to get l(A)
         l_of_A = simplify(l.subs(h, h_of_A))
-        print("l(A) =")
-        pprint(l_of_A)
+        if debug:
+            print("l(A) =")
+            pprint(l_of_A)
+
+        return l_of_A
+
+def godunov_solve(Q, c_init, c_inflow, delta, t):
+    gstep = lambda c, t: godunov_step(Q, delta, c, c_inflow)
+    return sp_int.odeint(gstep, c_init, t)
+    
+def main():
+    g, alpha, f = symbols("g,alpha,f")
+    
+    # Width of channel wrt vertical coordinate (y)
+    y = Symbol("y", positive=True, real=True)
+    a = Symbol("a", positive=True, real=True)
+    w = a #sqrt(y/a)
+
+    A_sym = Symbol("A", positive=True, real=True)
+
+    l_of_A = derive_wetted_perim(w, y, A_sym)
+    print(l_of_A)
+
+    u = simplify(sqrt(A_sym*g*sin(alpha) / (f * l_of_A)))
+    print("u(A) = ")
+    pprint(u)
     
+    Q = (A_sym*u).subs([("g", 9.81), ("alpha", 0.001), ("f", 0.05), (a, 100)])
+    #plot_expr(diff(Q, A_sym), A_sym, [0,30], 1000)
+    #return
+    Q_lamb = lambdify(A_sym, Q)
+
+    delta_x = 0.1
+    L = 100
+
+    delta_t = 0.1
+    T = 200
+    
+    x = np.array([i*delta_x for i in range(int(L/delta_x))])
+    A_init = 1 + 1000/(100+(x-10)**2)
+    
+    A_inflow = 1
+    t = [i*delta_t for i in range(int(T/delta_t))]
+    
+    c_sln = godunov_solve(Q_lamb, A_init, A_inflow, delta_x, t)
+    
+    plt.imshow(c_sln, origin="lower", aspect="auto")
+    plt.show()
+
+
 main()