Край из-за накопления

Question

Привет всем, я рискну провести симуляцию с python. Мне нужно было бы решить эту проблему: я моделирую движение по траектории канала, и, очевидно, материал накапливается, но у меня есть проблема с условиями на краях.

Может ли кто-нибудь помочь мне в случае, если я захочу сделать канал не таким акцентированным, более "тупым"?

Я оставляю вам код здесь ... любые предложения приветствуются. спасибо

nx, ny = 100,100          
X = np.zeros((nx, ny))    


pmt = 7                  
ppe=0.155                 


X[1:nx-1,1:ny-1] = np.random.randint(8,30, size=(nx-2, ny-2))
X[5:nx-1,1:ny-1] = 6.5*np.ones((nx-6, ny-2))    

plt.figure(figsize=(5,5))
plt.imshow(X, cmap='winter')
plt.colorbar()
plt.show()



#---------------------------------------------------------------------


alfa = 35   
pro = 0.3  
tau_s = 3   
s = 1.5    

#tau = tau_s/s


H_dm = 1/(pro*np.cos(math.radians(alfa))*np.sin(math.radians(alfa))) 
H_cm = H_dm*10          
H_m = (H_dm)*0.1       


#---------------------------------------------------------------------
def snow(X,ix,iy,H_m):                
    snow = H_m                        
    X1 = X.copy()
    edge = True                              
    if ix<2 or ix>nx-2 or iy<2 or iy>ny-2:   
        edge = False                          
        pass  
    if edge:
        X1[ix,iy] += snow                    
    return X1

def avalanche(X,pmt):          
    p_scale_factor = 4         
    X_av = X.copy()
    X_av_app = X.copy()
    for i in range(nx-1):
        if i<nx/5:       
            for j in range(ny-1):
                ix = i+1
                iy = j+1
                p_drop = np.tanh((X[ix,iy]-pmt)/p_scale_factor)  
                rand_drop = np.random.random()
                if rand_drop > 1-p_drop:
                    if j==0:    
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy] = X_av_app[ix+1,iy]+3*av_mass/4
                            X_av[ix+1,iy+1] = X_av_app[ix+1,iy+1]+av_mass/4
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass 

                    elif 0<j<ny-3:
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy-1] = X_av_app[ix+1,iy-1]+av_mass/4
                            X_av[ix+1,iy] = X_av_app[ix+1,iy]+av_mass/2
                            X_av[ix+1,iy+1] = X_av_app[ix+1,iy+1]+av_mass/4
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass
                    elif j==ny-3:
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy] = X_av_app[ix+1,iy]+3*av_mass/4
                            X_av[ix+1,iy-1] = X_av_app[ix+1,iy-1]+av_mass/4
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass 
                    X_av_app = X_av
        else:
            ilim1 = i-nx/5        #condizione che stabilisce da dove inizia il canalone
            ilim2 = nx-(i-nx/5)    #condizione che stabilisce da dove inizia il canalone
            for j in range(ny-1):
                ix = i+1
                iy = j+1
                p_drop = np.tanh((X[ix,iy]-pmt)/p_scale_factor)
                rand_drop = np.random.random()
                if rand_drop > 1-p_drop:
                    if j<ilim1:
                        pass
                    if j==ilim1:
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy+1] = X_av_app[ix+1,iy+1]+av_mass
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass 

                    elif ilim1<j<ilim2:
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy-1] = X_av_app[ix+1,iy-1]+av_mass/4
                            X_av[ix+1,iy] = X_av_app[ix+1,iy]+av_mass/2
                            X_av[ix+1,iy+1] = X_av_app[ix+1,iy+1]+av_mass/4
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass
                    elif j==ilim2:
                        if X[ix,iy] > pmt:
                            av_mass = (X[ix,iy]-pmt)*ppe
                            X_av[ix+1,iy-1] = X_av_app[ix+1,iy-1]+av_mass
                            X_av[ix,iy] = X_av_app[ix,iy]-av_mass 
                    X_av_app = X_av            

    return X_av



X_app = X.copy()
X_ev = np.zeros_like(X)
for i in range(nx-2):
    for j in range(ny-2):
        X_ev = snow(X_app,i+1,j+1,H_m)  
        X_app = X_ev