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$

marked as duplicate by David Ketcheson, Christian Clason, Anton Menshov, Mauro Vanzetto, Community Oct 16 '18 at 21:45

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 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$ – cpraveen 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$ – Herman Jaramillo 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$ – cpraveen 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$ – Herman Jaramillo Oct 15 '18 at 19:46