[Python]セミラグランジュ法のメソッド

移流方程式解くプログラムがネット上に少なかったのと、メソッド化されているものは日英で探した時に無かった(もしくはバグってるやつ)ので作りました。

#https://github.com/christian512/SemiLagPy/blob/master/advection1D.ipynb
#これを参考にした。

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()