Discrete transparent boundary conditions for leap-frog scheme applied to transport equation

$$\partial_t u +c_x\partial_x u + c_y\partial_y u=0.$$
In [1]:
#%matplotlib notebook
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
In [2]:
plt.rcParams['font.size'] = 18.0         # font size
plt.rcParams['mathtext.fontset'] = 'cm'  # computer moder math font
plt.rcParams['text.usetex'] = True       # use tex engine for everything (useful for the axes labels)
plt.rcParams['figure.figsize']=[9,6]

1D transport equation

Jump to 2D transport equation

Data definition

  • Velocity $c$
  • $x\in[x_\ell,x_r]$
  • final time $T$
  • number of grid points $J$
  • initial datum $u_0$
In [3]:
c=1
xl=-3
xr=3
T=10
J=1001
u0=lambda x: np.exp(-10*x**2)

Computation of some parameters

  • mesh size $\delta x$
  • CFL condition $\mu \leq cfl <1$, with $\mu=c \delta t/\delta x$
In [4]:
x=np.linspace(xl,xr,J)
dx=x[1]-x[0]
CFL=5./6.
dt=CFL*dx/np.abs(c)
N=np.floor(T/dt)
if np.mod(N,2): N+=1
N=np.int(N)
mu=c*dt/dx
if mu>=1 : print('Warning: CFL condition is not satisfied')

Computation of two initial states

The state $u^1$ is computed from $u^0$ with second order Lax-Wendroff scheme $$ \forall j=1,\cdots,J, \quad u_j^1=u_j^0-\frac{\mu}{2}(u_{j+1}^0-u_{j-1}^0)+\frac{\mu^2}{2}(u_{j+1}^0-2u_j^0+u_{j-1}^0). $$

In [5]:
un0=u0(x)

unjm1=np.zeros(np.shape(un0))
unjp1=np.zeros(np.shape(un0))
unjm1[1:]=un0[0:J-1]
unjp1[0:J-1]=un0[1:J]
un1=un0-(mu/2)*(unjp1-unjm1)+(mu**2/2)*(unjp1-2*un0+unjm1)

fig, axes = plt.subplots()
axes.plot(x,un0)
axes.set_xlabel('$x$')
axes.set_ylabel('initial datum $u^0$')
plt.show()

Leap-frog scheme with DTBC

The leap-frog scheme is given by $$ \forall j=1,\cdots,J, \quad u_j^{n+2}=u_j^n-\mu(u_{j+1}^{n+1}-u_{j-1}^{n+1}). $$ complemented with the discrete transparent boundary conditions $$ u_{0}^{n+2}=-\sum_{0\leq m\leq (n+1)/2} s_m^0u_1^{n+1-2m}, \qquad u_{J+1}^{n+2}=\sum_{0\leq m\leq (n+1)/2} s_m^0u_J^{n+1-2m}. $$

Computation of convolution kernel coefficients $(s_m^0)_m$

In [6]:
s=np.zeros(N+2)
s[0]=mu
s[1]=mu*(1-mu**2)
for j in range(3,N+3):
    s[j-1]=(2*j-3)*(1-2*mu**2)*s[j-2]/(j)-(j-3)*s[j-3]/(j)
alp=np.zeros(N+2)
alp[::2]=s[:N//2+1]

fig, axes = plt.subplots()
axes.plot(np.log10(np.abs(s)),'.')
axes.set_xlabel('$m$')
axes.set_ylabel('$\log_{10}|s_m^0|$')
plt.show()

Definition of some array to store:

  • Evolution of the solution at boundary nodes: u1, uJ
  • Boundary values: CLl, CLr
  • Evolution of $l^\infty$ norm: Linf, $l^\infty$ error: errLinf
  • Global evolution of the solution: pl
In [7]:
uJ=np.zeros(N+2); uJ[0]=0; uJ[1]=0;
u1=np.zeros(N+2); u1[0]=0; u1[1]=0;
CLr=0;CLl=0;
Linf=np.zeros(N+2); vt=Linf.copy(); errLinf=Linf.copy();
errLinf[0]=0; vt[0]=0;
errLinf[1]=np.max(np.abs(un1-u0(x-c*dt))); vt[1]=dt;
Linf[0]=np.max(np.abs(un0)); Linf[1]=np.max(np.abs(un1));
pl=np.zeros((N+2,J))
pl[0,:]=un0.copy()
pl[1,:]=un1.copy()

Main procedure

In [8]:
tloc=0;n=1;
fig, axes = plt.subplots()
axes.plot(x,un0)
unj=un0.copy(); unp1j=un1.copy();
#-------------------------------------------------------------------------#
# Main procedure
while tloc<T-dt: 
    #unp1jp1=np.concatenate((unp1j[1:J],np.array([CLr])))
    #unp1jm1=np.concatenate((np.array([CLl]),unp1j[:J-1]))
    unp1jp1=np.append(unp1j[1:J],CLr)
    unp1jm1=np.append(CLl,unp1j[:J-1])
    unp2j=unj-mu*(unp1jp1-unp1jm1)
    
    uJ[n+1]=unp2j[J-1]
    u1[n+1]=unp2j[0]
    
    v2=np.flip(alp[:n+1])
    v1=uJ[:n+1]
    CLr=np.dot(v1,v2)
    
    v1=u1[:n+1]
    CLl=-np.dot(v1,v2)
    
    errLinf[n+1]=np.max(np.abs(unp2j-u0(x-c*(n+1)*dt)))
    Linf[n+1]=np.max(np.abs(unp2j))    
    vt[n+1]=(n+1)*dt
    pl[n+1,:]=unp2j.copy()
    
    if np.mod(n,100)==0 :
#        plt.plot(x,unp2j-u0(x-(n+1)*dt))
        axes.plot(x,unp2j)
#        print('Time ',(n+1)*dt,'/',T)
    n=n+1
    tloc=tloc+dt    
    
    unj=unp1j.copy()
    unp1j=unp2j.copy()
    
axes.set_title('Evolution of $(u_j^n)$')
axes.set_xlabel('x')
plt.show()
In [9]:
fig, axes = plt.subplots(1,2)
fig.set_figwidth(16)
im0=axes[0].imshow(np.abs(pl),interpolation="bicubic",\
                   origin="lower",extent=[xl,xr,0,T],\
                   cmap=plt.get_cmap('jet'))
fig.colorbar(im0,ax=axes[0])
#im0.set_clim(-16,0)
axes[0].set_aspect(0.5)
axes[0].set_title('Evolution of $u_j^n$ with DTBC')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$t$')
im1=axes[1].imshow(np.log10(abs(pl)),interpolation="bicubic",\
                   origin="lower",extent=[xl,xr,0,T],\
                   cmap=plt.get_cmap('jet'))
fig.colorbar(im1,ax=axes[1])
im1.set_clim(-16,0)
axes[1].set_aspect(0.5)
axes[1].set_title('Evolution of $\log_{10}(u_j^n)$ with DTBC')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('$t$')
plt.show()
In [10]:
fig, axes = plt.subplots(1,2)
fig.set_figwidth(16)
axes[0].plot(vt,Linf)
axes[0].set_xlabel('$t$'); 
axes[0].set_title('Evolution of $\|u_j^n\|_{l^\infty}$')
axes[1].semilogy(vt,errLinf)
axes[1].set_xlabel('$t$'); 
axes[1].set_title('Evolution of $\|u_j^n-u(t^n,x_j)\|_{l^\infty}$')
plt.show()

Leap-frog scheme with Approximate DTBC (sums of exponentials)

The convolution kernel coefficients $(s_k^0)_k$ are approximated by the sum of exponentials $$ s_k^0 \approx \tilde{s}_k^0=\sum_{m=1}^M b_m q_m^{-k} $$ where $b_m$ and $q_m$ are computed with $P_N/Q_M$ Pade approximant of the power series $\sum_{k\in \mathbb{N}} s_k^0 x^k$.

Computation of parameters $(b,q)$ to approximate kernel coefficients with $(M,N)$ Pade approximants

We use the floating-point arithmetic with arbitrary accuracy Python library mpmath:

  • dps denotes the decimal accuracy

Definition of calc_sum_exp_coef function

In [11]:
from mpmath import mp
#mp.dps = 120
mp.dps = 80

def calc_sum_exp_coef(arg_list):
    # arg_list = [c,xl,xr,T,J,Nt,M,N]
    c=mp.mpf(arg_list[0])
    xl=mp.mpf(arg_list[1])
    xr=mp.mpf(arg_list[2])
    T=mp.mpf(arg_list[3])    
    J=int(arg_list[4])
    Nt=int(arg_list[5])
    M=int(arg_list[6])
    N=int(arg_list[7])
    L=M
    dx=(xr-xl)/(J-1)
    dt=T/Nt
    mu=c*dt/dx

    sL=[2*mu**2]
    sL=sL+[2*mu**2*(1-mu**2)]
    for j in range(3,Nt+1):
        stmp=(2*j-3)*(1-2*mu**2)*sL[j-1-1]/j-(j-3)*sL[j-1-2]/j
        sL=sL+[stmp]

    for j in range(Nt):
        sL[j]=sL[j]/2/mu

    s=mp.matrix(sL)
    
    Mat = mp.matrix(N+M)
    for k in range(M):
        for j in range(N+M-k):
            Mat[j+k,k]=s[j]

    for k in range(M,N+M):
        Mat[k-M,k]=mp.mpf(-1)

    A=Mat[0:M,0:M]
    B=Mat[0:M,M:N+M]
    C=Mat[M:N+M,0:M]

    rhs=mp.matrix(M,1)
    rhs[0]=mp.mpf(1)
    x = mp.lu_solve(A, rhs)
    Am=mp.matrix(M)
    for k in range(M):
        for j in range(M-k):
            Am[j+k,k]=x[j]

    rhs=-s[1:N+M+1]
    y1=rhs[0:M]
    y2=rhs[M:N+M]
    if N>0:
        CAmB=C*Am*B
        Pr=mp.lu_solve(CAmB,C*Am*y1-y2)
        Qr=Am*(y1-B*Pr)
    else:
        Qr=Am*y1

    P=mp.matrix(N+1,1)
    Q=mp.matrix(M+1,1)

    P[N]=s[0]
    Q[M]=1

    if N>0:
        for j in range(N):
            P[j]=Pr[N-1-j]

    for j in range(M):
        Q[j]=Qr[M-1-j]

    Qprime=mp.matrix(M,1)
    for j in range(M):
        Qprime[j]=(M-j)*Q[j]

    Qf=np.zeros(M)
    for j in range(M):
        Qf[j]=float(mp.nstr(Q[j],n=20))
    qf=np.roots(Qf)
    qf_init=mp.matrix(qf)

    #q=mp.polyroots(Q,maxsteps=500,cleanup=True, extraprec=240,roots_init=qf_init)
    q=mp.polyroots(Q,maxsteps=500,cleanup=True, extraprec=150,roots_init=qf_init)
    #q=mp.polyroots(Q,maxsteps=500,cleanup=True, extraprec=100)
    b=mp.matrix(M,1)
    for j in range(M):
        b[j]=(-mp.polyval(P,q[j])/mp.polyval(Qprime,q[j]))*q[j]**0
        
    # Test roots
    module_q=mp.matrix(M,1)
    for j in range(M):
        module_q[j]=mp.sqrt(q[j].real**2+q[j].imag**2)

    mqmin=module_q[0]
    for j in range(1,M):
        if module_q[j]<mqmin:
            mqmin=module_q[j]

    if mqmin<=1:
        print('There exists a root with modulus <=1')

    bf=np.zeros(M,dtype=complex)
    qf=np.zeros(M,dtype=complex)
    for j in range(M):
        bf[j]=float(mp.nstr(b[j].real,n=50))+1j*float(mp.nstr(b[j].imag,n=50))
        qf[j]=float(mp.nstr(q[j].real,n=50))+1j*float(mp.nstr(q[j].imag,n=50))
        
    return [bf,qf]

Computation of coefficients $(b_m,q_m)$ with Pade approximant $P_{N}/Q_{M}$

In [12]:
Ml=50;Nl=4;
[b,q]=calc_sum_exp_coef([c,xl,xr,T,J,N,Ml,Nl])
In [13]:
fig, axes = plt.subplots(1,2)
fig.set_figwidth(16)
th=np.linspace(0,2*np.pi,90)
axes[0].plot(np.cos(th),np.sin(th))
axes[0].plot(np.real(q),np.imag(q),'.')
axes[0].set_aspect(1.0)
axes[0].set_title('Roots of polynomial $Q_M$')
axes[0].set_xlabel('Re($q_m$)')
axes[0].set_ylabel('Im($q_m$)')

# Computation of approximate kernel coefficients
st=np.zeros(s.shape)
Q = np.array(q, dtype=np.complex256)
for n in range(st.size):
    for j in range(Ml):
#        st[n]+=np.real(b[j]/q[j]**(n+1))
        st[n]+=np.real(b[j]/Q[j])
    Q=Q*q
#axes[1].semilogy(np.abs(s),'.')
#axes[1].semilogy(np.abs(st),'o')
axes[1].plot(np.log10(np.abs(s)),'.',label='$\log_{10}(s^0_n)$')
axes[1].plot(np.log10(np.abs(st)),'o',label=r'$\log_{10}(\tilde{s}^0_n)$')
axes[1].set_aspect(80)
axes[1].set_xlabel('$n$')
axes[1].legend()
#axes[1].legend=(['$\log_{10}(s^0_n)$','$\log_{10}(\tilde{s}^0_n)$'])
plt.show()

Definition of various arrays

In [14]:
cof_l_even=np.zeros(Ml); cof_r_even=np.zeros(Ml)
cof_l_odd=np.zeros(Ml);  cof_r_odd=np.zeros(Ml)

CLr=0; CLl=0
pl=np.zeros((N+2,J))
pl[0,:]=un0.copy()
pl[1,:]=un1.copy()

Linf=np.zeros(N+2); vt=Linf.copy(); errLinf=Linf.copy();
errLinf[0]=0; vt[0]=0;
errLinf[1]=np.max(np.abs(un1-u0(x-c*dt))); vt[1]=dt;
Linf[0]=np.max(np.abs(un0)); Linf[1]=np.max(np.abs(un1));

Main procedure

In [15]:
tloc=0;n=1;
fig, axes = plt.subplots()
unj=un0.copy(); unp1j=un1.copy();
#-------------------------------------------------------------------------#
# Programme principal
for n in range(N):
    #unp1jp1=np.concatenate((unp1j[1:J],np.array([CLr])))
    #unp1jm1=np.concatenate((np.array([CLl]),unp1j[:J-1]))
    unp1jp1=np.append(unp1j[1:J],CLr)
    unp1jm1=np.append(CLl,unp1j[:J-1])
    unp2j=unj-mu*(unp1jp1-unp1jm1)
    
    if np.mod(n,2)==1:
        v1=cof_r_odd/q+unp1j[J-1]*b/q
        v2=cof_l_odd/q+unp1j[0]*b/q
        CLr=np.sum(np.real(v1))
        CLl=-np.sum(np.real(v2))
        cof_r_odd=v1
        cof_l_odd=v2
    else :
        v1=cof_r_even/q+unp1j[J-1]*b/q
        v2=cof_l_even/q+unp1j[0]*b/q
        CLr=np.sum(np.real(v1))
        CLl=-np.sum(np.real(v2))
        cof_r_even=v1
        cof_l_even=v2
    
    pl[n+2,:]=unp2j.copy()
    errLinf[n+2]=np.max(np.abs(unp2j-u0(x-c*(n+2)*dt)))
    Linf[n+2]=np.max(np.abs(unp2j))    
    vt[n+2]=(n+2)*dt    
    
    if np.mod(n,100)==0 :
#        plt.plot(x,unp2j-u0(x-(n+1)*dt))
        axes.plot(x,unp2j)
#        print('Temps ',(n+1)*dt,'/',T)

    unj=unp1j.copy()
    unp1j=unp2j.copy()
axes.set_title('Evolution of $(u_j^n)$')
axes.set_xlabel('x')
plt.show()
In [16]:
fig, axes = plt.subplots(1,2)
fig.set_figwidth(16)
im0=axes[0].imshow(np.abs(pl),interpolation="bicubic",\
                   origin="lower",extent=[xl,xr,0,T],cmap=plt.get_cmap('jet'))
fig.colorbar(im0,ax=axes[0])
#im0.set_clim(-16,0)
axes[0].set_aspect(0.5)
axes[0].set_title('Evolution of $u_j^n$ with ADTBC')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$t$')
im1=axes[1].imshow(np.log10(abs(pl)),interpolation="bicubic",\
                   origin="lower",extent=[xl,xr,0,T],cmap=plt.get_cmap('jet'))
fig.colorbar(im1,ax=axes[1])
im1.set_clim(-16,0)
axes[1].set_aspect(0.5)
axes[1].set_title('Evolution of $\log_{10}(u_j^n)$ with ADTBC')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('$t$')
plt.show()
In [17]:
fig, axes = plt.subplots(1,2)
fig.set_figwidth(16)
axes[0].plot(vt,Linf)
axes[0].set_xlabel('$t$'); 
axes[0].set_title('Evolution of $\|u_j^n\|_{l^\infty}$')
axes[1].semilogy(vt,errLinf)
axes[1].set_xlabel('$t$'); 
axes[1].set_title('Evolution of $\|u_j^n-u(t^n,x_j)\|_{l^\infty}$')
plt.show()

2D transport equation

Jump to 1D transport equation

Data definition

  • Velocity $\mathbf{c}=(c_x,c_y)$
  • $x\in[x_\ell,x_r]$, $y\in[y_b,y_t]$
  • final time $T$
  • number of grid points $J\times K$
  • mesh size $\delta x$, $\delta y$
  • CFL condition $|\mu_x|+|\mu_y| \leq cfl <1$, with $\mu_x=c_x \delta t/\delta x$ and $\mu_y=c_y \delta t/\delta y$
  • initial datum $u_0$
In [18]:
# Velocity
cx=1;cy=0.1;
# Domain [xl,xr] x [yb,yt] x [0,T]
xl=-3; xr=3;
yb=-2; yt=2;
T=8;
# Mesh, with CFL
J=301;K=201;
x=np.linspace(xl,xr,J)
y=np.linspace(yb,yt,K)
dx=x[1]-x[0];dy=y[1]-y[0];
CFL=0.5;
dt=CFL*(dx*dy)/(abs(cx)*dy+abs(cy)*dx);
X,Y=np.meshgrid(x,y);
mux=cx*dt/dx;
muy=cy*dt/dy;
# Number of time steps
N=np.floor(T/dt)
if np.mod(N,2): N+=1
N=np.int(N)

# Initial datum
u0=np.exp(-5*X**2-5*Y**2)
un=u0.copy()

Definition of discrete derivative operators, $\forall (j,k)\in \mathbb{Z}^2$ $$ \mathrm{Dxp}(u_{j,k})=(u_{j+1,k}-u_{j,k})/\delta x, \quad \mathrm{Dxm}(u_{j,k})=(u_{j,k}-u_{j-1,k})/\delta x,\quad \mathrm{Dx0}(u_{j,k})=(u_{j+1,k}-u_{j-1,k})/2\delta x $$ $$ \mathrm{Dyp}(u_{j,k})=(u_{j,k+1}-u_{j,k})/\delta y, \quad \mathrm{Dym}(u_{j,k})=(u_{j,k}-u_{j,k-1})/\delta y,\quad \mathrm{Dy0}(u_{j,k})=(u_{j,k+1}-u_{j,k-1})/2\delta y $$

In [19]:
# Periodic finite difference operators
Dxp=lambda Z: (np.roll(Z,-1,axis=1)-Z);
Dxm=lambda Z: (Z-np.roll(Z,+1,axis=1));
Dx0=lambda Z: (np.roll(Z,-1,axis=1)-np.roll(Z,+1,axis=1));

Dyp=lambda Z: (np.roll(Z,-1,axis=0)-Z);
Dym=lambda Z: (Z-np.roll(Z,+1,axis=0));
Dy0=lambda Z: (np.roll(Z,-1,axis=0)-np.roll(Z,+1,axis=0));

# Dirichlet finite difference operators
# Centered derivative with respect to x, Zl at left boundary node
# and Zr at right boundary node
def Dx0D(Z,Zl,Zr):
    Zjp2k=np.roll(Z,-1,axis=1); Zjp2k[:,-1]=Zr;
    Zjm1k=np.roll(Z,+1,axis=1); Zjm1k[:,0]=Zl;
    Res=Zjp2k-Zjm1k;
    return Res

# Centered derivative with respect to y, Zb at bottom boundary node
# and Zt at top boundary node
def Dy0D(Z,Zb,Zt):
    Zjkp2=np.roll(Z,-1,axis=0); Zjkp2[-1,:]=Zt;
    Zjkm1=np.roll(Z,+1,axis=0); Zjkm1[0,:]=Zb;
    Res=Zjkp2-Zjkm1;
    return Res

Computation of two initial states

The state $u^1$ is computed from $u^0$ with second order Lax-Wendroff scheme $\forall j=1,\cdots,J$ and $\forall k=1,\cdots,K$ $$ \begin{array}{cl} u_{j,k}^1=u_{j,k}^0&- \frac{\mu_x}{2}(u_{j+1,k}^0-u_{j-1,k}^0) -\frac{\mu_y}{2}(u_{j,k+1}^0-u_{j,k-1}^0)\\ & +\frac{\mu_x^2}{2}(u_{j+1,k}^0-2u_{j,k}^0+u_{j-1,k}^0) +\frac{\mu_y^2}{2}(u_{j,k+1}^0-2u_{j,k}^0+u_{j,k-1}^0)\\ & +\frac{\mu_x\, \mu_y}{4} \, (u_{j+1,k+1}^0-u_{j+1,k-1}^0-u_{j-1,k+1}^0+u_{j-1,k-1}^0) \end{array} $$

In [20]:
# Step 1: Computation of u^{n+1} with Lax Wendroff 
zerX=np.zeros(K)
zerY=np.zeros(J)
DxO=Dx0(un); DyO=Dy0(un);
D2x=Dxp(Dxm(un)); D2y=Dyp(Dym(un));
Dxy=Dx0(Dy0(un));
# Homogeneous Dirichlet BC wrt X and Y
DxO[:,1-1]=zerX.copy(); DyO[1-1,:]=zerY.copy();
DxO[:,J-1]=zerX.copy(); DyO[K-1,:]=zerY.copy();
D2x[:,1-1]=zerX.copy(); D2y[1-1,:]=zerY.copy();
D2x[:,J-1]=zerX.copy(); D2y[K-1,:]=zerY.copy();
Dxy[:,1-1]=zerX.copy(); Dxy[1-1,:]=zerY.copy();
Dxy[:,J-1]=zerX.copy(); Dxy[K-1,:]=zerY.copy();

unp1=un-(mux/2)*(DxO)+(mux**2/2)*(D2x)-(muy/2)*(DyO)+(muy**2/2)*(D2y)+(mux*muy/4)*Dxy;
u1=unp1.copy()

Leap-frog scheme with DTBC

The leap-frog scheme is given by $$ u_{j,k}^{n+2}-u_{j,k}^n +\mu_x \, \big( u_{j+1,k}^{n+1}-u_{j-1,k}^{n+1} \big) +\mu_y \, \big( u_{j,k+1}^{n+1}-u_{j,k-1}^{n+1} \big) \, = \, 0 \, $$ which holds for $1 \le j \le J$ and $1 \le k \le K$ complemented with the approximate discrete transparent boundary conditions of various order

  • Order 0 $$ \begin{align} u_{0,k}^{n+2} \, &= \, -\sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{1,k}^{n+1-2\, m}, \\ u_{J+1,k}^{n+2} \, & = \, \sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{J,k}^{n+1-2\, m}, \\ u_{j,0}^{n+2} \, & = \, -\sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,1}^{n+1-2\, m}, \\ u_{j,K+1}^{n+2} \, &= \, \sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,K}^{n+1-2\, m} \end{align} $$

  • Order 1 $$ \begin{align} u_{0,k}^{n+2} \, &= \, -\sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{1,k}^{n+1-2\, m} -\sum_{1 \le m \le (n+2)/2} \, s_m^1 \, (u_{1,k+1}^{n+2-2\, m}-u_{1,k-1}^{n+2-2\, m}), \\ u_{J+1,k}^{n+2} \, & = \, \sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{J,k}^{n+1-2\, m} +\sum_{1 \le m \le (n+2)/2} \, s_m^1 \, (u_{J,k+1}^{n+2-2\, m}-u_{J,k-1}^{n+2-2\, m}), \\ u_{j,0}^{n+2} \, & = \, -\sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,1}^{n+1-2\, m} -\sum_{1 \le m \le (n+2)/2} \, t_m^1 \, (u_{j+1,1}^{n+2-2\, m}-u_{j-1,1}^{n+2-2\, m}, \\ u_{j,K+1}^{n+2} \, &= \, \sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,K}^{n+1-2\, m} +\sum_{1 \le m \le (n+2)/2} \, t_m^1 \, (u_{j+1,K}^{n+2-2\, m}-u_{j-1,K}^{n+2-2\, m}) \end{align} $$

  • Order 2 $$ \begin{align} u_{0,k}^{n+2} \, &= \,-\sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{1,k}^{n+1-2\, m} -\sum_{1 \le m \le (n+2)/2} \, s_m^1 \, (u_{1,k+1}^{n+2-2\, m}-u_{1,k-1}^{n+2-2\, m}) \notag \\ & \qquad \qquad \qquad -\sum_{1 \le m \le (n+1)/2} \, s_m^2 \, \big( u_{1,k+1}^{n+1-2\, m}-2\, u_{1,k}^{n+1-2\, m}+u_{1,k-1}^{n+1-2\, m} \big) , \\ u_{J+1,k}^{n+2} \, & = \,\sum_{0 \le m \le (n+1)/2} \, s_m^0 \, u_{J,k}^{n+1-2\, m} +\sum_{1 \le m \le (n+2)/2} \, s_m^1 \, (u_{J,k+1}^{n+2-2\, m}-u_{J,k-1}^{n+2-2\, m}) \notag \\ & \qquad \qquad \qquad +\sum_{1 \le m \le (n+1)/2} \, s_m^2 \, \big( u_{J,k+1}^{n+1-2\, m}-2\, u_{J,k}^{n+1-2\, m}+u_{J,k-1}^{n+1-2\, m} \big) , \\ u_{j,0}^{n+2} \, & = \, -\sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,1}^{n+1-2\, m} -\sum_{1 \le m \le (n+2)/2} \, t_m^1 \, (u_{j+1,1}^{n+2-2\, m}-u_{j-1,1}^{n+2-2\, m}) \notag \\ & \qquad \qquad \qquad -\sum_{1 \le m \le (n+1)/2} \, t_m^2 \, \big( u_{j+1,1}^{n+1-2\, m}-2\, u_{j,1}^{n+1-2\, m}+u_{j-1,1}^{n+1-2\, m} \big), \\ u_{j,K+1}^{n+2} \, &= \, \sum_{0 \le m \le (n+1)/2} \, t_m^0 \, u_{j,K}^{n+1-2\, m} +\sum_{1 \le m \le (n+2)/2} \, t_m^1 \, (u_{j+1,K}^{n+2-2\, m}-u_{j-1,K}^{n+2-2\, m}) \notag \\ & \qquad \qquad \qquad +\sum_{1 \le m \le (n+1)/2} \, t_m^2 \, \big( u_{j+1,K}^{n+1-2\, m}-2\, u_{j,K}^{n+1-2\, m}+u_{j-1,K}^{n+1-2\, m} \big) \end{align} $$

Computation of convolution kernel coefficients $(s_m^0)_m$, $(s_m^1)_m$, $(s_m^2)_m$, $(t_m^0)_m$, $(t_m^1)_m$ and $(t_m^2)_m$

In [21]:
sx=np.zeros(N+2); #s_m^0
sx[0]=mux
for j in range(N+1):
    v=sx[0:j+1]
    sx[j+1]=sx[j]-mux*np.dot(np.flip(v),v)

sy=np.zeros(N+2); #t_m^0
sy[0]=muy
for j in range(N+1):
    v=sy[0:j+1]
    sy[j+1]=sy[j]-muy*np.dot(np.flip(v),v)

tx=np.zeros(N+2); #s_m^1
for j in range(N+1):
    v=tx[0:j+1]
    w=sx[0:j+1]
    tx[j+1]=tx[j]-2*mux*np.dot(np.flip(w),v)-muy*sx[j]

ty=np.zeros(N+2); #t_m^1
for j in range(N+1):
    v=ty[0:j+1]
    w=sy[0:j+1]
    ty[j+1]=ty[j]-2*muy*np.dot(np.flip(w),v)-mux*sy[j]

rx=np.zeros(N+2); #s_m^2
rx[0]=0
rx[1]=4*mux*muy**2
for p in range(N):
    v=rx[0:p+2]
    w=sx[0:p+2]
    ww=tx[0:p+3]
    rx[p+2]=rx[p+1]-2*mux*np.dot(np.flip(w),v)\
            -4*muy*tx[p+2]\
            -4*mux*np.dot(np.flip(ww),ww)

ry=np.zeros(N+2); #t_m^2
ry[0]=0
ry[1]=4*muy*mux**2
for p in range(N):
    v=ry[0:p+2]
    w=sy[0:p+2]
    ww=ty[0:p+3]
    ry[p+2]=ry[p+1]-2*muy*np.dot(np.flip(w),v)\
            -4*mux*ty[p+2]\
            -4*muy*np.dot(np.flip(ww),ww)

fig, axes = plt.subplots(3,2)
fig.set_figwidth(16)
fig.set_figheight(16)
axes[0,0].plot(np.log10(np.abs(sx)),'.')
axes[0,1].plot(np.log10(np.abs(sy)),'.')
axes[1,0].plot(np.log10(np.abs(tx[1:])),'.')
axes[1,1].plot(np.log10(np.abs(ty[1:])),'.')
axes[2,0].plot(np.log10(np.abs(rx[2:])),'.')
axes[2,1].plot(np.log10(np.abs(ry[2:])),'.')
axes[0,0].set_xlabel('$m$'); axes[0,0].set_title('$s_m^0$')
axes[0,1].set_xlabel('$m$'); axes[0,1].set_title('$t_m^0$')
axes[1,0].set_xlabel('$m$'); axes[1,0].set_title('$s_m^1$')
axes[1,1].set_xlabel('$m$'); axes[1,1].set_title('$t_m^1$')
axes[2,0].set_xlabel('$m$'); axes[2,0].set_title('$s_m^2$')
axes[2,1].set_xlabel('$m$'); axes[2,1].set_title('$t_m^2$')
plt.subplots_adjust(top = 0.99, bottom=0.01, hspace=0.3, wspace=0.15)
plt.show()

Definition of some local arrays

In [22]:
# Arrays to save the boundary values along iterations
uJk=np.zeros((N+2,K)); uJk[0,:]=zerX.copy(); uJk[1,:]=zerX.copy();
u1k=np.zeros((N+2,K)); u1k[0,:]=zerX.copy(); u1k[1,:]=zerX.copy();
ujK=np.zeros((N+2,J)); ujK[0,:]=zerY.copy(); ujK[1,:]=zerY.copy();
uj1=np.zeros((N+2,J)); uj1[0,:]=zerY.copy(); uj1[1,:]=zerY.copy();

# Arrays to store the boundary conditions
CLr=zerX.copy();CLl=zerX.copy();
CLb=zerY.copy();CLt=zerY.copy();
# 
#          0 1                         J J+1     <---j
# k = K+1    |                         |
# k = K  ----|----------CLt------------|-----
#            |                         |
#            |                         |
#           CLl                       CLr
#            |                         |
#            |                         |
# k = 1  ----|----------CLb------------|-----
# k = 0      |                         |

# Arrays to store the values of ghost nodes
GhostT=np.zeros((N+2,2))
GhostB=np.zeros((N+2,2))
GhostL=np.zeros((N+2,2))
GhostR=np.zeros((N+2,2))

# Arrays for global evolution plots
pl=np.zeros((N+2,K,J))
pl[0,:,:]=un.copy()
pl[1,:,:]=unp1.copy()

Main procedure

In [23]:
# Order = 0;
# Order = 1;
Order = 2; # (2,1)
# Order = 3; # (2,2) unstable

alpx=np.zeros(N+2); alpx[::2]=sx[:N//2+1]
alpy=np.zeros(N+2); alpy[::2]=sy[:N//2+1]
betx=np.zeros(N+2); betx[::2]=tx[1:N//2+2]
bety=np.zeros(N+2); bety[::2]=ty[1:N//2+2]
if Order == 2:
    gamx=np.zeros(N+2); gamx[::2]=rx[1:N//2+2]
    gamy=np.zeros(N+2);
elif Order == 3:
    gamx=np.zeros(N+2); gamx[::2]=rx[1:N//2+2]
    gamy=np.zeros(N+2); gamy[::2]=ry[1:N//2+2]

tloc=0;cpt=1;
un=u0.copy()
unp1=u1.copy()
# Main procedure
while tloc<T-dt:
    Dx0np1=Dx0D(unp1,CLl,CLr)
    Dy0np1=Dy0D(unp1,CLb,CLt)
    
    unp2=un-mux*Dx0np1-muy*Dy0np1
    
    #/----\____/----\____/----\____/----\____/----\____/----\____/----\____
    # Save of boundary traces
    uJk[cpt+2-1,:]=unp2[:,J-1];
    u1k[cpt+2-1,:]=unp2[:,1-1];
    ujK[cpt+2-1,:]=unp2[K-1,:];
    uj1[cpt+2-1,:]=unp2[1-1,:];
    
    #/----\____/----\____/----\____/----\____/----\____/----\____/----\____
    # Computation of convolutions
    # Right boundary
    v1=alpx[:cpt+1]
    for k in range(K):
        v2=uJk[:cpt+1,k]
        CLr[k]=np.dot(np.flip(v2),v1)
    if Order>=1:
        v1=betx[:cpt]
        for k in range(1,K-1):
            v2=uJk[:cpt,k+1]-uJk[:cpt,k-1]
            CLr[k]+=np.dot(np.flip(v2),v1)
        v2=uJk[:cpt,1]-GhostB[:cpt,1]
        CLr[0]+=np.dot(np.flip(v2),v1)
        v2=GhostT[:cpt,1]-uJk[:cpt,K-2]
        CLr[K-1]+=np.dot(np.flip(v2),v1)
        if Order >=2:
            v1=gamx[:cpt-1]
            for k in range(1,K-1):
                v2=uJk[:cpt-1,k+1]-2*uJk[:cpt-1,k]+uJk[:cpt-1,k-1]
                CLr[k]+=np.dot(np.flip(v2),v1)
            v2=uJk[:cpt-1,1]-2*uJk[:cpt-1,0]+GhostB[:cpt-1,1]
            CLr[0]+=np.dot(np.flip(v2),v1)
            v2=GhostT[:cpt-1,1]-2*uJk[:cpt-1,K-1]+uJk[:cpt-1,K-2]
            CLr[K-1]+=np.dot(np.flip(v2),v1)
    
    # Left boundary   
    v1=-alpx[:cpt+1]
    for k in range(K):
        v2=u1k[:cpt+1,k]
        CLl[k]=np.dot(np.flip(v2),v1)
    if Order>=1:
        v1=-betx[:cpt]
        for k in range(1,K-1):
            v2=u1k[:cpt,k+1]-u1k[:cpt,k-1]
            CLl[k]+=np.dot(np.flip(v2),v1)
        v2=u1k[:cpt,1]-GhostB[:cpt,0]
        CLl[0]+=np.dot(np.flip(v2),v1)
        v2=GhostT[:cpt,0]-u1k[:cpt,K-2]
        CLl[K-1]+=np.dot(np.flip(v2),v1)
        if Order >=2:
            v1=-gamx[:cpt-1]
            for k in range(1,K-1):
                v2=u1k[:cpt-1,k+1]-2*u1k[:cpt-1,k]+u1k[:cpt-1,k-1]
                CLl[k]+=np.dot(np.flip(v2),v1)
            v2=u1k[:cpt-1,1]-2*u1k[:cpt-1,0]+GhostB[:cpt-1,0]
            CLl[0]+=np.dot(np.flip(v2),v1)
            v2=GhostT[:cpt-1,0]-2*u1k[:cpt-1,K-1]+u1k[:cpt-1,K-2]
            CLl[K-1]+=np.dot(np.flip(v2),v1)
    
    # Top boundary   
    v1=alpy[:cpt+1]
    for j in range(J):
        v2=ujK[:cpt+1,j]
        CLt[j]=np.dot(np.flip(v2),v1)   
    if Order>=1:
        v1=bety[:cpt]
        for j in range(1,J-1):
            v2=ujK[:cpt,j+1]-ujK[:cpt,j-1]
            CLt[j]+=np.dot(np.flip(v2),v1)
        v2=ujK[:cpt,1]-GhostL[:cpt,1]
        CLt[0]+=np.dot(np.flip(v2),v1)
        v2=GhostR[:cpt,1]-ujK[:cpt,J-2]
        CLt[J-1]+=np.dot(np.flip(v2),v1)       
        if Order >=2:
            v1=gamy[:cpt-1]
            for j in range(1,J-1):
                v2=ujK[:cpt-1,j+1]-2*ujK[:cpt-1,j]+ujK[:cpt-1,j-1]
                CLt[j]+=np.dot(np.flip(v2),v1)
            v2=ujK[:cpt-1,1]-2*ujK[:cpt-1,0]+GhostL[:cpt-1,1]
            CLt[0]+=np.dot(np.flip(v2),v1)
            v2=GhostR[:cpt-1,1]-2*ujK[:cpt-1,J-1]+ujK[:cpt-1,J-2]
            CLt[J-1]+=np.dot(np.flip(v2),v1)
    
    # Bottom boundary    
    v1=-alpy[:cpt+1]
    for j in range(J):
        v2=uj1[:cpt+1,j]
        CLb[j]=np.dot(np.flip(v2),v1)        
    if Order>=1:
        v1=-bety[:cpt]
        for j in range(1,J-1):
            v2=uj1[:cpt,j+1]-uj1[:cpt,j-1]
            CLb[j]+=np.dot(np.flip(v2),v1)
        v2=uj1[:cpt,1]-GhostL[:cpt,0]
        CLb[0]+=np.dot(np.flip(v2),v1)
        v2=GhostR[:cpt,0]-uj1[:cpt,J-2]
        CLb[J-1]+=np.dot(np.flip(v2),v1)     
        if Order >=2:
            v1=-gamy[:cpt-1]
            for j in range(1,J-1):
                v2=uj1[:cpt-1,j+1]-2*uj1[:cpt-1,j]+uj1[:cpt-1,j-1]
                CLb[j]+=np.dot(np.flip(v2),v1)
            v2=uj1[:cpt-1,1]-2*uj1[:cpt-1,0]+GhostL[:cpt-1,0]
            CLb[0]+=np.dot(np.flip(v2),v1)
            v2=GhostR[:cpt-1,0]-2*uj1[:cpt-1,J-1]+uj1[:cpt-1,J-2]
            CLb[J-1]+=np.dot(np.flip(v2),v1)
    
    cpt=cpt+1;
    tloc=tloc+dt;
    pl[cpt+1-1,:,:]=unp2.copy();

    GhostT[cpt,0]=CLt[0]; GhostT[cpt,1]=CLt[J-1]
    GhostB[cpt,0]=CLb[0]; GhostB[cpt,1]=CLb[J-1]
    GhostL[cpt,0]=CLl[0]; GhostL[cpt,1]=CLl[K-1]
    GhostR[cpt,0]=CLr[0]; GhostR[cpt,1]=CLr[K-1]           
    un=unp1.copy()
    unp1=unp2.copy()
In [24]:
vt=np.arange(0,T+2*dt,dt)
fig, axes = plt.subplots(3,2)
fig.set_figwidth(16); fig.set_figheight(16)
fig.suptitle('Evolution of $u_{j,k}^n$')
j=0;
for k in range(3):
    for l in range(2):
        if (j==0):
            lab=0
        else:
            lab=abs(np.ceil((vt.size-1)*j/(6-1))-1); lab=int(lab)
        KK=pl[lab,:,:]
        im=axes[k,l].imshow(np.abs(KK),interpolation="bicubic",\
                            origin="lower",extent=[xl,xr,yb,yt],cmap=plt.get_cmap('jet'))
        axes[k,l].set_aspect(1)
        axes[k,l].set_title('$t = {:.3}$'.format(vt[lab]))
        im.set_clim(0,1)
        fig.colorbar(im,ax=axes[k,l])
        j+=1
In [25]:
vt=np.arange(0,T+2*dt,dt)
fig, axes = plt.subplots(3,2)
fig.set_figwidth(16); fig.set_figheight(16)
fig.suptitle('Evolution of $|u_{j,k}^n|$')
j=0;
for k in range(3):
    for l in range(2):
        if (j==0):
            lab=0
        else:
            lab=abs(np.ceil((vt.size-1)*j/(6-1))-1); lab=int(lab)
        KK=pl[lab,:,:]
        im=axes[k,l].imshow(np.abs(KK),interpolation="bicubic",\
                            origin="lower",extent=[xl,xr,yb,yt],cmap=plt.get_cmap('jet'))
        axes[k,l].set_aspect(1)
        axes[k,l].set_title('$t = {:.3}$'.format(vt[lab]))
        fig.colorbar(im,ax=axes[k,l],format='%.0e')
        j+=1
In [26]:
vt=np.arange(0,T+2*dt,dt)
fig, axes = plt.subplots(3,2)
fig.set_figwidth(16); fig.set_figheight(16)
fig.suptitle('Evolution of $\log_{10}(|u_{j,k}^n|)$')
j=0;
for k in range(3):
    for l in range(2):
        if (j==0):
            lab=0
        else:
            lab=abs(np.ceil((vt.size-1)*j/(6-1))-1); lab=int(lab)
        KK=np.log10(np.abs(pl[lab,:,:]))
        im=axes[k,l].imshow(KK,interpolation="bicubic",\
                            origin="lower",extent=[xl,xr,yb,yt],cmap=plt.get_cmap('jet'))
        axes[k,l].set_aspect(1)
        axes[k,l].set_title('$t = {:.3}$'.format(vt[lab]))
        fig.colorbar(im,ax=axes[k,l])
        im.set_clim(-16,0)
        j+=1
In [ ]:
In [ ]:
In [ ]: