U1_new = U1.copy()
        U2_new = U2.copy()
        for j in range(1, N+1): # finite volume discretization
            U1_new[j] = U1[j] - (dt/dx)*(F1[j] - F1[j-1])       # Continuity Equation update water depth h
            U2_new[j] = U2[j] - (dt/dx)*(F2L[j] - F2R[j-1])     # Momentum Equation h*v
        U1, U2 = U1_new, U2_new
        h = U1
        v = U2 / U1
        eta = h[1:N+1] + b_phys
        h_store.append(eta.copy())
        v_store.append(v[1:N+1].copy())
        t_store.append(t)
        plt.cla()
        plt.subplot(2,1,1)
        plt.plot(x, eta, label="η (free surface)")
        plt.plot(x, b_phys, '--', label="bottom")
        plt.plot(x, eta_initial, ':', label="initial η")
        plt.legend()
        plt.subplot(2,1,2)
        plt.plot(x, v[1:N+1])
        plt.pause(pausetime)
    pbar.close()
    h_array = np.array(h_store)
    v_array = np.array(v_store)
    np.savez("dam_break_data_bottomtopo.npz",
             h_array=h_array,
             v_array=v_array,
             t_store=np.array(t_store),
             x=x)
    print("Saved: dam_break_data_bottomtopo.npz")
    plt.close()
SWE_DamBreak_JA()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

data = np.load("dam_break_data_bottomtopo.npz") # Load data
h_array = data["h_array"]   
v_array = data["v_array"]
t_store = data["t_store"]
x = data["x"]
nt, nx = h_array.shape
print("Loaded:", h_array.shape, v_array.shape)
def bottom(x):
    return np.zeros_like(x) # Reconstruct bottom
b = bottom(x)
h0 = np.min(h_array) # Base water level
fig, axes = plt.subplots(3, 1, figsize=(15, 8), sharex=True)
fig.suptitle("1D Shallow Water - Dam Break (Filled Water + Velocity)", fontsize=12)
ax0, ax1, ax2 = axes
h_min, h_max = np.min(h_array), np.max(h_array) # Axis limit
v_min, v_max = np.min(v_array), np.max(v_array)
b_min = np.min(b)
for ax in (ax0, ax1):
    ax.set_xlim(x.min(), x.max())
    ax.set_ylim(b_min - 1, h_max + 2)
    ax.set_ylabel("Height (h)")
    ax.grid()
ax2.set_xlim(x.min(), x.max())
ax2.set_ylim(v_min * 1.05, v_max * 1.05)
ax2.set_ylabel("Velocity (v)")
ax2.set_xlabel("Distance (m)")
ax2.grid()
ax0.set_title("Initial Water Surface", color="midnightblue")
ax1.set_title("Wave Propagation", color="midnightblue")
ax2.set_title("Velocity", color="midnightblue")