2
$\begingroup$

I am implementing the advection equation $u_x+(1/c)u_t=0$ following a Crank-Nicholson finite difference scheme. The equation for this is \begin{eqnarray*} -\frac{\gamma}{4} w_{n-3 j+1} + w_{n-2 j+1} + \frac{\gamma}{4} w_{n-1 j+1} = \frac{\gamma}{4} w_{n-3 j} + w_{n-2 j} - \frac{\gamma}{4} w_{n-1 j}. \end{eqnarray*}

In matrix form it can be written as:

\begin{eqnarray*} A \bf{w}_{j+1} = \it B {\bf{w}_j } + \bf{b} \end{eqnarray*} where

\begin{eqnarray*} A = \begin{pmatrix} 1 & \frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\ -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} \\ 0 & 0 & \cdots & \cdots & -\frac{\gamma}{4} & 1 \end{pmatrix} \; , \; \bf{w}_{j+1} = \begin{pmatrix} w_{1 j+1} \\ w_{2 j+1} \\ \vdots \\ \vdots \\ \vdots \\ w_{n-2 j+1} \\ \end{pmatrix} \end{eqnarray*}

\begin{eqnarray*} B = \begin{pmatrix} 1 & -\frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\ \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} \\ 0 & 0 & \cdots & \cdots & \frac{\gamma}{4} & 1 \end{pmatrix} \; , \; \bf{w}_{j} = \begin{pmatrix} w_{1 j} \\ w_{2 j} \\ \vdots \\ \vdots \\ \vdots \\ w_{n-2 j} \\ \end{pmatrix} \; , \; \bf{b} = \begin{pmatrix} \frac{\gamma}{4} w_{0 j+1} + \frac{\gamma}{4} w_{0j} \\ 0 \\ \vdots \\ 0 \\ -\frac{\gamma}{4} w_{n-1 j+1} - \frac{\gamma}{4} w_{n-1 j} \end{pmatrix} \end{eqnarray*}

As an initial condition I have a Gaussian pulse $u(x,0)=\mathrm{e}^{-400(x-1/2)^2}$. The spacial domain is $x \in [0,1]$. Constant velocity $c=1$, $dt=0.001$, $dx=0.002$ and $\gamma=c dt/dx$.

My implementation is in Python and I can provide it upon request.

I attach a few screen captures of the "wavefield".

Time 0.3 seconds

time 0.4 seconds

time 0.6 seconds

time 0.7 seconds

Can someone please help me understand what is going on after the pulse hits the left boundary? I am using Dirichlet boundary conditions with 0 on the boundary.

Thanks.

*****Update*****

I found someone else asking the same question and some interesting answers which I am still trying to digest. Still I would like to have more discussion on this issue.

Here is the link on this

$\endgroup$
6
  • 2
    $\begingroup$ If your speed $c > 0$ then you can use Dirichlet bc only at left boundary. What are you doing in your implementation ? Also it is standard to use a superscript to denote time level, and I think thats a good convention. $\endgroup$
    – cfdlab
    Oct 15 '18 at 13:50
  • $\begingroup$ @PraveenChandrashekar Thanks for your input. Yes, I think you are right about the boundary condition. Mathematically sounds right, but numerically, how do I deal with data beyond the right boundary on my centered (CN) scheme? You are right also about the convention of the superscript for time. Thanks. $\endgroup$ Oct 15 '18 at 13:57
  • $\begingroup$ There is no unique way to deal with the right boundary, which is an artificial boundary. You can use backward difference for the last point instead of central difference. Or put a ghost point beyond right boundary, extrapolate solution from the last two points to ghost point, and then you have enough stencil to apply your scheme. For $c>0$ the initial profile you put will just move out of the right boundary and after some time you are just left with a constant solution, unless you put a time dependent bc at the left. $\endgroup$
    – cfdlab
    Oct 15 '18 at 15:06
  • $\begingroup$ You need to use a ghost cell or point which it a value outside your domain to help you setup the boundary conditions. Normally when doing so the actual value of the ghost point cancels out, leaving the problem well posed within the domain. $\endgroup$
    – boyfarrell
    Oct 15 '18 at 18:35
  • $\begingroup$ @DavidKetcheson: You are right about that. Even there is another question posted on the StackExchange computational Science that is equivalent. scicomp.stackexchange.com/questions/21734/… While the questions are equivalent we would like to keep all answers and comments. What is the best way to achieve this? Thanks. $\endgroup$ Oct 15 '18 at 19:46