[Python] Méthode Semi-Lagrange

Je l'ai fait parce qu'il y avait peu de programmes sur le net pour résoudre l'équation de translocation, et il n'y en avait pas de méthodique quand je l'ai recherché en japonais et en anglais (ou un buggy).

#https://github.com/christian512/SemiLagPy/blob/master/advection1D.ipynb
#J'en ai parlé.

import numpy as np
import matplotlib.pyplot as plt

class SemiLagrangian1D:
    """
    SemiLagrangian scheme for periodic 1D problems.
    """
    def __init__(self,x,y,u,dt):
        """
        Initizalizes the SemiLagrangian object given initial configuration
        x : x positions (must be equidistant)
        y : variable at all x positions at time t
        u : speed of advection
        dt: time step
        """
        # Check dimensions
        assert x.shape[0] == y.shape[0], "x and y have different length"
        # Get the x stepsize
        dx = x[1] - x[0]
        # check if x is equidistant
        assert x[-1] == x[0] + dx*(x.shape[0]-1), "x array not equidistant"
        # Set class variables
        self.x = x
        self.y = y
        self.u = u
        self.dt = dt
        self.dx = dx

    def cubic_evolve(self,nt=1):
        """
        Propagates the variable using a cubic interpolation scheme.
        nt: number of time stepsize
        """
        # loop through time steps
        for l in range(nt):
            # temporary array
            y_temp = np.zeros(self.y.shape[0])
            # loop through array
            for i in range(self.y.shape[0]):
                # idx left to departure point
                x_dep = self.x[i]-self.u*self.dt
                j = int(np.floor(x_dep/self.dx))
                # alpha
                a = (self.x[i]-self.u*self.dt - j*self.dx)/self.dx
                # calculate next time step
                f = lambda x: x % self.y.shape[0] if x >= self.y.shape[0] else x
                y_temp[i] = - a * (1-a)*(2-a)/6 * self.y[f(j-1)]
                y_temp[i] += (1-a**2)*(2-a)/2 * self.y[f(j)]
                y_temp[i] += a*(1+a)*(2-a)/2 * self.y[f(j+1)]
                y_temp[i] -= a*(1-a**2)/6 * self.y[f(j+2)]
            self.y = np.copy(y_temp)
        return self.y

    def linear_interporation(self,nt=1):
        for l in range(nt):
            # temporary array
            y_tmp = np.zeros(self.y.shape[0])
            for i in range(1,NX-1):
                x_tmp = self.x[0]+i*self.dx
                x_tmp -= self.u*self.dt

                ib = int((x_tmp-self.x[0])/self.dx)
                if ib < 0: ib = 0
                if ib >= len(self.x): ib = NX-2
                distance = (x_tmp - ib*self.dx)/DX

                y_tmp[i] = (1-distance)*self.y[ib]+distance*self.y[ib+1]

        self.y = np.copy(y_tmp)
        return self.y





L = 5 # Domain length
NX = 1000 # Number of grid points
X = np.linspace(0,L,num=NX) # Array of grid points
DX = X[1] - X[0] # spatial grid point distance
DT = 0.01 # time step size
NT = 1000 # number of time steps
CA = (1 + 0.02*np.sqrt(130)) # advection speed
print('Courant= ' + str(CA*DT/DX))

# Sine function
F_ini2 = np.zeros(X.shape[0])
F_ini2[0:int(X.shape[0]/5)] = np.sin(np.pi*X[0:int(X.shape[0]/5)])

# Plotting
# plt.plot(X,F_ini2,label='Sine function')
# plt.legend()
# plt.show()

# Create array for storing evolution of F(x,t)
F_arr = np.zeros([NT,X.shape[0]])
F_arr[0,:] = np.copy(F_ini2)
# Create SemiLagrangian object
sl = SemiLagrangian1D(X,F_arr[0,:],CA,DT)

# Calculate the variable for all time steps using semi-lagrangian
for i in range(1,NT):
    F_arr[i,:] = sl.cubic_evolve(nt=1)
    #F_arr[i,:] = sl.linear_interporation(nt=i)
#ans = sl.linear_interporation(nt=300)

# Plot initial vs. evolution
# print(F_arr[300,:])
print("f(0.5,t_0) = {}".format(F_arr[0,int(X.shape[0]/10)]))
print("f(x1,t1) = {}".format(F_arr[300,int(X.shape[0]/10) + int(-CA*X.shape[0]/5)]))
plt.plot(X,F_arr[0,:],label='Initial')
plt.plot(X,F_arr[300,:],label='Evolved')
plt.legend()
plt.show()