<jupyter_start><jupyter_text>---- Two-sided Ritter dam break$x$-direction<jupyter_code># Finite volume with Rusanov flux scheme
# start time
time   = 0
dtplot = 0.20
tplot  = dtplot
# initialisation
x   = np.linspace(xmin,xmax,nx)
y   = np.linspace(ymin,ymax,ny)
X,Y = np.meshgrid(x,y,indexing='ij')
h   = 0*X
#h   = np.exp(-((X-mu)**2+(Y-mu)**2)/sigma**2)
#h   = np.where(X <= xmax/2, 2., 0.01)
h   = np.where(np.sqrt((X-mu)**2)<= radius,2.,0.01)
hu  = 0*X
hv  = 0*X
# add ghost cells
X  = addGhostCells2D(X,xmin,xmax,ymin,ymax)
Y  = addGhostCells2D(Y,xmin,xmax,ymin,ymax)
h  = addGhostCells2D(h,0,0,0,0)
hu = addGhostCells2D(hu,0,0,0,0)
hv = addGhostCells2D(hv,0,0,0,0)
# initial plot
YM = Y[int(ny/2),int(ny/2)]
print(int(ny/2),int(ny/2),YM)
plt.figure(figsize=(10,6))
plt.ylim(0,2.1)
plt.xlabel('x [m]')
plt.ylabel('Water height [m]')
plt.title('Ritter dam break')
#plt.contour(X[1:-1,1:-1],Y[1:-1,1:-1],h[1:-1,1:-1],levels=[0.5,0.8,1.5])
plt.plot(X[1:-1,1:-1][Y[1:-1,1:-1]==YM],
         h[1:-1,1:-1][Y[1:-1,1:-1]==YM],linewidth=4,color='gray',label='h$_0$')
# solution
while (time < T):
    # set boundary conditions
    h  = bcNeumann2D(h)
    hu = bcNeumann2D(hu)
    hv = bcNeumann2D(hv)
    
    h = np.maximum(0.01,h)
    c = np.sqrt(g*h)
    u = hu/h #uu*h/h
    v = hv/h #vv*h/h
    
    maxeigenu = np.max(np.maximum(np.abs(u-c),np.abs(u+c)))
    maxeigenv = np.max(np.maximum(np.abs(v-c),np.abs(v+c)))
    maxeigen = max(maxeigenu,maxeigenv)
    dt = dtCourant(dx,dy,maxeigenu,maxeigenv,time,T,CFL)
    #print('dt: ',dx,dt,time,maxeigen)
    # advance time
    time  = time+dt
    
    Fhup = h[1:,:]*u[1:,:]
    Fhum = h[:-1,:]*u[:-1,:]
    Fhvp = h[:,1:]*v[:,1:]
    Fhvm = h[:,:-1]*v[:,:-1]
    
    Fhu2p = h[1:,:]*u[1:,:]**2 + g*h[1:,:]**2/2
    Fhu2m = h[:-1,:]*u[:-1,:]**2 + g*h[:-1,:]**2/2
    Fhuvp = h[:,1:]*u[:,1:]*v[:,1:]
    Fhuvm = h[:,:-1]*u[:,:-1]*v[:,:-1]
    
    Fhuv2p = h[1:,:]*u[1:,:]*v[1:,:]
    Fhuv2m = h[:-1,:]*u[:-1,:]*v[:-1,:]
    Fhv2p = h[:,1:]*v[:,1:]**2 + g*h[:,1:]**2/2
    Fhv2m = h[:,:-1]*v[:,:-1]**2 + g*h[:,:-1]**2/2
    
    # R = Lax-Friedrichs Flux
    Rhx  = (Fhup+Fhum)/2 - dx/dt/2*(h[1:,:]-h[:-1,:])
    Rhy  = (Fhvp+Fhvm)/2 - dy/dt/2*(h[:,1:]-h[:,:-1])
    
    Rhux  = (Fhu2p+Fhu2m)/2 - dx/dt/2*(hu[1:,:]-hu[:-1,:])
    Rhuy  = (Fhuvp+Fhuvm)/2 - dy/dt/2*(hu[:,1:]-hu[:,:-1])
    
    Rhvx  = (Fhuv2p+Fhuv2m)/2 - dx/dt/2*(hv[1:,:]-hv[:-1,:])
    Rhvy  = (Fhv2p+Fhv2m)/2 - dy/dt/2*(hv[:,1:]-hv[:,:-1])
    
    # update inner points
    h[1:-1,1:-1] = h[1:-1,1:-1] - dt/dx*(Rhx[1:,1:-1]-Rhx[:-1,1:-1]) \
                                - dt/dy*(Rhy[1:-1,1:]-Rhy[1:-1,:-1])
    hu[1:-1,1:-1] = hu[1:-1,1:-1] - dt/dx*(Rhux[1:,1:-1]-Rhux[:-1,1:-1]) \
                                - dt/dy*(Rhuy[1:-1,1:]-Rhuy[1:-1,:-1])
    hv[1:-1,1:-1] = hv[1:-1,1:-1] - dt/dx*(Rhvx[1:,1:-1]-Rhvx[:-1,1:-1]) \
                                - dt/dy*(Rhvy[1:-1,1:]-Rhvy[1:-1,:-1])
    if (time > tplot):
        #plt.contour(X,Y,h,levels=[0.5,0.8,1.5])
        plt.plot(X[1:-1,1:-1][Y[1:-1,1:-1]==YM],
                 h[1:-1,1:-1][Y[1:-1,1:-1]==YM],label='h(t='+str(round(time,2))+' s)')
        plt.legend()
        tplot = tplot + dtplot<jupyter_output>100 100 19.8<jupyter_text>---- Two-sided Ritter dam break$y$-direction<jupyter_code># Finite volume with Rusanov flux scheme
# start time
time   = 0
dtplot = 0.20
tplot  = dtplot
# initialisation
x   = np.linspace(xmin,xmax,nx)
y   = np.linspace(ymin,ymax,ny)
X,Y = np.meshgrid(x,y,indexing='ij')
h   = 0*X
h   = np.where(np.sqrt((Y-mu)**2)<= radius,2.,0.01)
hu  = 0*X
hv  = 0*X
# add ghost cells
X  = addGhostCells2D(X,xmin,xmax,ymin,ymax)
Y  = addGhostCells2D(Y,xmin,xmax,ymin,ymax)
h  = addGhostCells2D(h,0,0,0,0)
hu = addGhostCells2D(hu,0,0,0,0)
hv = addGhostCells2D(hv,0,0,0,0)
# initial plot
XM = X[int(nx/2),int(nx/2)]
print(int(nx/2),int(nx/2),XM)
plt.figure(figsize=(10,6))
plt.ylim(0,2.1)
plt.xlabel('y [m]')
plt.ylabel('Water height [m]')
plt.title('Ritter dam break')
#plt.contour(X[1:-1,1:-1],Y[1:-1,1:-1],h[1:-1,1:-1],levels=[0.5,0.8,1.5])
plt.plot(Y[1:-1,1:-1][X[1:-1,1:-1]==XM],
         h[1:-1,1:-1][X[1:-1,1:-1]==XM],linewidth=4,color='gray',label='h$_0$')
# solution
while (time < T):
    # set boundary conditions
    h  = bcNeumann2D(h)
    hu = bcNeumann2D(hu)
    hv = bcNeumann2D(hv)
    
    h = np.maximum(0.01,h)
    c = np.sqrt(g*h)
    u = hu/h #uu*h/h
    v = hv/h #vv*h/h
    
    maxeigenu = np.max(np.maximum(np.abs(u-c),np.abs(u+c)))
    maxeigenv = np.max(np.maximum(np.abs(v-c),np.abs(v+c)))
    maxeigen  = max(maxeigenu,maxeigenv)
    dt = dtCourant(dx,dy,maxeigenu,maxeigenv,time,T,CFL)
    #print('dt: ',dx,dt,time,maxeigen)
    # advance time
    time  = time+dt
    
    Fhup = h[1:,:]*u[1:,:]
    Fhum = h[:-1,:]*u[:-1,:]
    Fhvp = h[:,1:]*v[:,1:]
    Fhvm = h[:,:-1]*v[:,:-1]
    
    Fhu2p = hu[1:,:]*u[1:,:] + 0.5*g*h[1:,:]**2
    Fhu2m = hu[:-1,:]*u[:-1,:] + 0.5*g*h[:-1,:]**2
    Fhuvp = h[:,1:]*u[:,1:]*v[:,1:]
    Fhuvm = h[:,:-1]*u[:,:-1]*v[:,:-1]
    
    Fhuv2p = h[1:,:]*u[1:,:]*v[1:,:]
    Fhuv2m = h[:-1,:]*u[:-1,:]*v[:-1,:]
    Fhv2p = hv[:,1:]*v[:,1:] + 0.5*g*h[:,1:]**2
    Fhv2m = hv[:,:-1]*v[:,:-1] + 0.5*g*h[:,:-1]**2
    
    # R = Lax-Friedrichs Flux
    Rhx  = (Fhup+Fhum)/2 - dx/dt/2*(h[1:,:]-h[:-1,:])
    Rhy  = (Fhvp+Fhvm)/2 - dy/dt/2*(h[:,1:]-h[:,:-1])
    
    Rhux  = (Fhu2p+Fhu2m)/2 - dx/dt/2*(hu[1:,:]-hu[:-1,:])
    Rhuy  = (Fhuvp+Fhuvm)/2 - dy/dt/2*(hu[:,1:]-hu[:,:-1])
    
    Rhvx  = (Fhuv2p+Fhuv2m)/2 - dx/dt/2*(hv[1:,:]-hv[:-1,:])
    Rhvy  = (Fhv2p+Fhv2m)/2 - dy/dt/2*(hv[:,1:]-hv[:,:-1])
    
    # update inner points
    h[1:-1,1:-1] = h[1:-1,1:-1] - dt/dx*(Rhx[1:,1:-1]-Rhx[:-1,1:-1]) \
                                - dt/dy*(Rhy[1:-1,1:]-Rhy[1:-1,:-1])
    hu[1:-1,1:-1] = hu[1:-1,1:-1] - dt/dx*(Rhux[1:,1:-1]-Rhux[:-1,1:-1]) \
                                - dt/dy*(Rhuy[1:-1,1:]-Rhuy[1:-1,:-1])
    hv[1:-1,1:-1] = hv[1:-1,1:-1] - dt/dx*(Rhvx[1:,1:-1]-Rhvx[:-1,1:-1]) \
                                - dt/dy*(Rhvy[1:-1,1:]-Rhvy[1:-1,:-1])
    if (time > tplot):
        #plt.contour(X,Y,h,levels=[0.5,0.8,1.5])
        plt.plot(Y[1:-1,1:-1][X[1:-1,1:-1]==XM],
                 h[1:-1,1:-1][X[1:-1,1:-1]==XM],label='h(t='+str(round(time,2))+' s)')
        plt.legend()
        tplot = tplot + dtplot<jupyter_output>100 100 19.8<jupyter_text>---- Two-sided Ritter dam break$x$- and $y$-directions<jupyter_code># Finite volume with Rusanov flux scheme
# start time
time   = 0
dtplot = 0.20
tplot  = dtplot
# initialisation
x   = np.linspace(xmin,xmax,nx)
y   = np.linspace(ymin,ymax,ny)
X,Y = np.meshgrid(x,y,indexing='ij')
h   = 0*X
#h   = np.where(np.sqrt((Y-mu)**2)<= radius,2.,0.01)
h   = np.where((np.sqrt((X-mu)**2)<= radius) & (np.sqrt((Y-mu)**2)<= radius),2.,0.01)
hu  = 0*X
hv  = 0*X
# add ghost cells
X  = addGhostCells2D(X,xmin,xmax,ymin,ymax)
Y  = addGhostCells2D(Y,xmin,xmax,ymin,ymax)
h  = addGhostCells2D(h,0,0,0,0)
hu = addGhostCells2D(hu,0,0,0,0)
hv = addGhostCells2D(hv,0,0,0,0)
# initial plot
XM = X[int(nx/2),int(nx/2)]
print(int(nx/2),int(nx/2),XM)
plt.figure(figsize=(10,6))
plt.ylim(0,2.1)
plt.xlabel('y [m]')
plt.ylabel('Water height [m]')
plt.title('Ritter dam break')
#plt.contour(X[1:-1,1:-1],Y[1:-1,1:-1],h[1:-1,1:-1],levels=[0.5,0.8,1.5])
plt.plot(Y[1:-1,1:-1][X[1:-1,1:-1]==XM],
         h[1:-1,1:-1][X[1:-1,1:-1]==XM],linewidth=4,color='gray',label='h$_0$')
# solution
while (time < T):
    # set boundary conditions
    h  = bcNeumann2D(h)
    hu = bcNeumann2D(hu)
    hv = bcNeumann2D(hv)
    
    h = np.maximum(0.1,h)
    c = np.sqrt(g*h)
    u = hu/h #uu*h/h
    v = hv/h #vv*h/h
    hu = np.where(hu < 0.1, 0, hu)
    hv = np.where(hv < 0.1, 0, hv)
    
    maxeigenu = np.max(np.maximum(np.abs(u-c),np.abs(u+c)))
    maxeigenv = np.max(np.maximum(np.abs(v-c),np.abs(v+c)))
    maxeigen  = max(maxeigenu,maxeigenv)
    dt = dtCourant(dx,dy,maxeigenu,maxeigenv,time,T,CFL)
    #print('dt: ',dx,dt,time,maxeigen)
    # advance time
    time  = time+dt
    
    Fhup = hu[1:,:]
    Fhum = hu[:-1,:]
    Fhvp = hv[:,1:]
    Fhvm = hv[:,:-1]
    
    Fhu2p = hu[1:,:]*hu[1:,:]/h[1:,:] + 0.5*g*h[1:,:]**2
    Fhu2m = hu[:-1,:]*hu[:-1,:]/h[:-1,:] + 0.5*g*h[:-1,:]**2
    Fhuvp = hu[:,1:]*hv[:,1:]/h[:,1:]
    Fhuvm = hu[:,:-1]*hv[:,:-1]/h[:,:-1]
    
    Fhuv2p = hu[1:,:]*hv[1:,:]/h[1:,:]
    Fhuv2m = hu[:-1,:]*hv[:-1,:]/h[:-1,:]
    Fhv2p  = hv[:,1:]*hv[:,1:]/h[:,1:] + 0.5*g*h[:,1:]**2
    Fhv2m  = hv[:,:-1]*hv[:,:-1]/h[:,:-1] + 0.5*g*h[:,:-1]**2
    
    # R = Lax-Friedrichs Flux
    Rhx  = (Fhup+Fhum)/2 - dx/dt/2*(h[1:,:]-h[:-1,:])
    Rhy  = (Fhvp+Fhvm)/2 - dy/dt/2*(h[:,1:]-h[:,:-1])
    
    Rhux = (Fhu2p+Fhu2m)/2 - dx/dt/2*(hu[1:,:]-hu[:-1,:])
    Rhuy = (Fhuvp+Fhuvm)/2 - dy/dt/2*(hu[:,1:]-hu[:,:-1])
    
    Rhvx = (Fhuv2p+Fhuv2m)/2 - dx/dt/2*(hv[1:,:]-hv[:-1,:])
    Rhvy = (Fhv2p+Fhv2m)/2 - dy/dt/2*(hv[:,1:]-hv[:,:-1])
    
    # update inner points
    h[1:-1,1:-1]  = h[1:-1,1:-1] - dt/dx*(Rhx[1:,1:-1]-Rhx[:-1,1:-1]) \
                                 - dt/dy*(Rhy[1:-1,1:]-Rhy[1:-1,:-1])
    hu[1:-1,1:-1] = hu[1:-1,1:-1] - dt/dx*(Rhux[1:,1:-1]-Rhux[:-1,1:-1]) \
                                  - dt/dy*(Rhuy[1:-1,1:]-Rhuy[1:-1,:-1])
    hv[1:-1,1:-1] = hv[1:-1,1:-1] - dt/dx*(Rhvx[1:,1:-1]-Rhvx[:-1,1:-1]) \
                                  - dt/dy*(Rhvy[1:-1,1:]-Rhvy[1:-1,:-1])
    if (time > tplot):
        #plt.contour(X,Y,h,levels=[0.5,0.8,1.5])
        plt.plot(Y[1:-1,1:-1][X[1:-1,1:-1]==XM],
                 h[1:-1,1:-1][X[1:-1,1:-1]==XM],label='h(t='+str(round(time,2))+' s)')
        plt.legend()
        tplot = tplot + dtplot<jupyter_output>100 100 19.8
{'fraction_non_alphanumeric': 0.3060524126455907, 'fraction_numerical': 0.08090682196339434, 'mean_word_length': 2.76399217221135, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '331163', 'license': 'MIT', 'n_tokens_mistral': 5489, 'n_tokens_neox': 5257, 'n_words': 562, 'path': '.ipynb_checkpoints/Dynamics_lab09_SWE2DwithLocalLaxFriedrich_extended-checkpoint.ipynb', 'repo_name': 'georgkaufmann/Lecture_Dynamics', 'size': '9616'}