20

Starting with the advection equation is conservative form, $$ \frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t) $$ The Crank-Nicolson method consists of a time averaged centered difference. $$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \...


16

The fundamental quantity in transport is the flux, $\mathbf v u$ for advection. The divergence theorem states that $$\int_\Omega \nabla\cdot (\mathbf v u) = \int_{\partial \Omega} (\mathbf v u) \cdot \mathbf n . $$ An equation is conservative when it is preserves this equality. Dropping to 1D with $\Omega = (a,b)$ and using the equation $u_t + (\mathbf v ...


8

I see several issues: The DFT computed with fft puts the zero mode at the beginning of the array, and if you want to compute the derivative, it is necessary to apply fftshift/ifftshift to the array N to make sure the derivative is correct. It is easy to see for yourself what the correct expression is by working it out with pen and paper, and see also the ...


7

One way or the other, you need to stabilize every discretization for this equation. Traditionally, this was done using methods such as artificial viscosity or its slightly smarter sibling SUPG. But there are many other alternatives to make things work if you'd like to stick with the FEM framework -- e.g., discontinuous Galerkin methods. My summary of this ...


6

Let me give a part of your answer, I would need some more indications from your side to answer you fully. So please read, and write some comments so that I can complete my answer. About notations in numerical methods There are a few mistakes in the way you write your equations. These are only details, but for someone who is used to it, it can a bit ...


6

Numerical diffusion arises from a first-order finite difference approximation to the spatial derivative $\partial u/\partial x$. To see how this is the case, examine the Taylor series expansion for $u_{i+1}$: $$ u(x_{i+1}) = u(x_{i}) + \left.\frac{\partial u}{\partial x}\right|_{x_i} (\Delta x) + \frac{1}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_{...


6

The stationary equation you show transports information from the right to the left via the advection term; it also diffuses slightly. If you switch off the diffusion term altogether, then you only have transport from the right to the left, and you need to also drop the boundary condition at the left: because information is from the right to the left, nothing ...


6

Okay. Let's begin with the first situation. Your equation is: $$ \frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}=0\tag{1}$$ $\textbf{Previous comments}$ There are plenty of webs that can explain you what happens with the error when you does not choose properly the values of your mesh to approximate the derivatives in an equation (Von Neumann ...


5

You solve the 1D-advection equation with c a constant velocity : $$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0~~~~~~~~(1) $$ When you discretize this equation (with an explicit scheme in time and an upwind scheme in space for instance), you get : $$ \frac{u_i^{n+1} - u_i^n}{\Delta t}+\frac{u_i^n-u_{i-1}^{n}}{\Delta x} = 0 ~~~~~~~~(2) $$ ...


5

Your image of the numerical domain of dependence is correct. But try to also draw the analytical domain of dependence, maybe this could help you to better understand what is going on. Note that the analytic solution of $u_t+au_x=0$ is $u(t,x)=u_0(x-at)$. So the slope of the actual dependence is $a$. The CFL condition just says that the numerical dependence ...


5

Eigenvalues with zero eigenvalue correspond to purely oscillatory modes. You can see it by diagonalising the system. Your matrix $A$ can be written as \begin{equation} A = P \Sigma P^{-1} \end{equation} where $\Sigma$ is a diagonal matrix with entries $\lambda_n$ corresponding to the eigenvalues of $Q$. You can now transform your ODE to \begin{equation} dQ/...


5

Look at the cell Péclet number $$\mathrm{Pe}_h = h v / K$$ where $h$ is mesh size, $v$ is the magnitude of velocity, and $K$ is diffusivity. It is analogous to cell Reynolds number for the momentum equation and is small when "thermal diffusivity is large compared to advection". It is common common in macro-scale fluid dynamics that thermal diffusivity $K$...


5

Bluntly speaking, SUPG and alike and RANS are different approaches to different problems that, however, have the same name - instability - and the same phenomenology - the failure of numerical routines. RANS is used to cope with turbulence as an instability of the equation. If a flow is or becomes turbulent the describing equations are instable, e.g. ...


5

Your equation can be written in the following fashion (any spatial derivative approximation is valid), once space is discretised: $$\frac{1}{c}\frac{du_i}{dt}=-\left(\frac{\partial u}{\partial x}\right)_i(t) + v_i(t) \tag{*}$$ Keep in mind that $v_i(t) = v(x_i,t)$. Now the system of equations depends only on time $t$ you can apply Crack Nicholson method to ...


5

There is a difference between the requirements for a hyperbolic pde like $$ u_t + a u_x = 0 $$ and for a purely parabolic pde like $$ u_t = u_{xx} $$ Suppose the solutions are smooth and you approximate them by some finite difference method. Then in case of hyperbolic problem, the maximum error in the numerical solution depends on the time interval of ...


4

In 1-D, the first discretization you've presented is correct. The matrix equation does not look right, though. For starters, it's not clear what your boundary conditions would be or how you would incorporate them. In $N$ dimensions, the advection equation looks like \begin{align} \frac{\partial{u}}{\partial{t}} + \sum_{i=1}^{N}\frac{\partial}{\partial{x^{i}...


4

This is actually a standard modeling problem if you consider the medium that flows through the network to be incompressible (e.g., liquids, or gases at low velocity). Then, you formulate everything in terms of fluxes (liters or kg per second) rather than in discrete parcels. The key realization is that the flux that goes into one end of the pipe equals the ...


4

You stumbled across something very fundamental that is important to understand. For any PDE we can define for each point its "domain of dependence" This is the region of space/time that is able to have an effect on the solution at that point. For your advection problem, it is just the characteristic line through that point. If the point of interest is $x_0,...


4

The difficulty is relative to something, in this case it is relative to diffusion dominated problems. Diffusion dominated aren't "easy" either, they have their own set of problems. I'll start with some favorable qualities of diffusion problems, and then mention why they are not present for advection problems: Discretizations are often symmetric and ...


3

In the time-dependent nonlinear case, if you drop the diffusive term then you have a nonlinear hyperbolic problem. Solutions will naturally generate singularities (discontinuities) in finite time. To extend the solution beyond that time, one must consider weak solutions, and uniqueness is lost. To specify a unique, physically relevant solution one ...


3

In the 2-D case, this corresponds to an elliptic fixed point (an orbital, I believe). You might look into Lyapunov stability, hopefully someone will be able to recommend a good resource on that. I'll try to have another look later to see if I can find resources on this. UPDATE For related material, you might look into limit cycles, or the Poincaré-...


3

The confusion was the misleading variables $F_{j-1/2}$ and $F_{j+1/2}$ with f_left and f_right, which are completely different. f_left and f_right are the interpolated fluxes at a one single face. They must be then upwinded using the advection speed to compute the Flux at a specific cell face. Which means if $C>0$ we take f_left, otherwise we take ...


3

The key steps are to consider the advection equation $u_t + au_x = 0$ where $a=\omega/k$ is the advection speed. Exact solutions to this equation is of the form $u(x,t) = f(x-at)$, where $f(y)$ is an arbitrary function. For example, discretize using a standard Galerkin method we derive the weak form $\int_\Omega v u_t dx + \int_\Omega a u_x = 0$ Assuming ...


3

Your description "the function $h(y,x)$ gives the proportion of the mass $u(y,t)$ at position $y$ that moves to position $x$ in space in a small unit of time" seems to indicate a slight misunderstanding of derivatives. Since the term on the left of your equation ($\partial_t u$) is a time derivative, the "small unit of time" must be infinitesimally small. ...


3

There are multiple questions here, but let's start with the basics. You have written two hyperbolic PDEs; (1) is the continuity equation, which is conservative and (2) is the color equation, which is not conservative. What are the characteristic speeds for these equations? For (1), you have stated correctly that the characteristic speed is $u(x,t)$. ...


3

and it works just fine for c>0 with the exact same boundary conditions, though. That's your issue. Look up the method of characteristics: https://web.stanford.edu/class/math220a/handouts/firstorder.pdf The characteristics for a 1D advection problem flow in one direction. So there are two things involved. One is that the BCs on one side of the interval ...


3

The Numpy function $roll$ performs periodic shift of an array. Using it, the explicit time step for your PDE in a periodic domain can be simply implemented like this: u = u - (1/2)*(c*dt/dx)*(np.roll(u,-1) - np.roll(u,1))


3

The choice of $k$ is restricted also by the discretization of the source term. To see it, rewrite your scheme to \begin{equation} u_m^{n+1} = \left(1 - \frac{k(1-x_m)}{h} - k(1-x_m)\right) u_m^n + \frac{k (1-x_m)}{h} u_{m+1}^n \,. \end{equation} You need $$ 1 - \frac{k(1-x_m)}{h} - k(1-x_m) \ge 0 $$ for all $x_m$. Taking $x_m=0$ (the worst case scenario) you ...


2

User 03161 asserts that the Crank Nicolson method is not appropriate for advection problems, but boyfarrell provides a working code with results visualized in a movie. In fact they are both correct, but neither gives the full perspective. At the beginning of boyfarrell's answer the correct C-N formula for linear advection is written out. In that formula are ...


2

You have discretized an advection equation using a forward difference in time and centered differences in space. You have correctly deduced that this is an unstable discretization; in fact it is unstable even for constant-coefficient advection in one dimension. There are many stable discretizations you could use; the most common (and simplest) is to switch ...


Only top voted, non community-wiki answers of a minimum length are eligible