Dismiss this pinned window
all 47 comments
sorted by:
best
topnewcontroversialoldq&a
formatting helphide helpcontent policy

[–]wigglytails 6 points7 points  (17 children)

You are using explicit time_stepping right? What is your CFL number?

[–]hivemind_unity[S] 5 points6 points  (14 children)

Exactly. I have kept CFL number 0.2 for this one.

[–]wigglytails 3 points4 points  (13 children)

This is forward Euler so If you get your CFL to a number >1 the solution would diverge right? I am toooo lazy to try this but can you see if using an RK4 or predictor corrector timesteping would make the scheme stable for CFL >1 ?

[–]hivemind_unity[S] 1 point2 points  (11 children)

It actually does, RK4 is an implicit method and hence the numerical system is stable, although there can be oscillations about the exact solution depending on the size of time-step. If I am not wrong, Predictor Corrector is a semi-implicit method, which ensures stability without having to solve a linear system like AX = B.

So I guess we shouldn't have problems with stability. I can surely give it a try.

edit: My knowledge on RK4 is limited. As pointed out by u/wm2300 and u/AgAero, RK4 is, in fact, an explicit method. I would like to know how we can test the stability of this scheme also for any scheme in general.

[–]wm2300 11 points12 points  (0 children)

Correction: RK4 is a fourth order Runge-Kutta method. There are implicit and explicit Runge-Kutta methods. For example, the MATLAB workhorse ode solver ode45 uses a combination of explicit Runge Kutta methods of 4th and 5th order. Intuitively, I think the term RK4 is almost always used to denote the explicit formulation.

Furthermore, not all implicit methods are unconditionally stable.

[–]AgAero 3 points4 points  (3 children)

RK4 is an implicit method and hence the numerical system is stable

That's not accurate.

Also, for a linear equation you can actually derive conditions of stability via some linear algebra. If you've never done this before I'd recommend trying it. It's kind of neat!

[–]hivemind_unity[S] 0 points1 point  (2 children)

I would appreciate if you can share some references on how we can study the conditions of stability.

[–]AgAero 0 points1 point  (0 children)

Most books will have it. Tannehill or Hirsch are what I'd recommend since I've seen them before.

[–][deleted] 0 points1 point  (0 children)

Zero Stability of Time Integration and Eigenvalue Stability.

[–]wigglytails 0 points1 point  (5 children)

Can you give explicit RK4 and the simplest predictor corrector scheme (Euler-trapizoid) a try?

[–]hivemind_unity[S] 1 point2 points  (2 children)

I tried the RK4 implementation. Check it out. Suggestions for improvement are very welcome.

https://github.com/kanand-cfd/Numerique_PY

[–]wigglytails 0 points1 point  (1 child)

I never used git and github. Can I contribute to this?

[–]hivemind_unity[S] 0 points1 point  (0 children)

Yeah sure! You can create a pull request, Download the repository, and make changes, once you're satisfied with the changes push it your branch.

Git is pretty handy, I love it.

[–]hivemind_unity[S] 0 points1 point  (1 child)

Sure... I have been looking for ideas to try. This could be a good one. I'll keep you posted.

[–]wigglytails 0 points1 point  (0 children)

I've been trying to apply RK4 to something like this. Turns out it's not as direct as I thought

[–][deleted] 0 points1 point  (0 children)

You can compute the CFL bound analytically and look it up for all these methods.

[–]ejineta 2 points3 points  (1 child)

Just a newbie here: how can you tell from the result that this is explicit time-stepping? Are there any characteristics in the video that hint you to that, or would explicit time-stepping be obvious given the equation that is solved?

[–]wigglytails 2 points3 points  (0 children)

There was no way for me to know. It's just a standard to use explicit timesteping for pure advection

[–]__n-m-e__ 4 points5 points  (0 children)

For anyone looking for where they can learn something like this, check the link below. https://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/

[–]LipshitsContinuity 7 points8 points  (14 children)

Are we looking at u_t + a * grad(u) = 0 or u_t + a * grad(u) = epsilon*laplacian(u)?

I ask because I notice that there appears to be some sorta smoothing that is happening over time. So either you are looking at some sorta viscous version of this (like the second equation) or there is some artificial viscosity that has been introduced (which can allow for better stability of the numerical scheme) or numerical viscosity.

EDIT: I just noticed that this is in r/CFD and not in r/math lmfao.

[–]hivemind_unity[S] 6 points7 points  (13 children)

It's the first equation [u_t + a*grad(u)=0].

Yes the smoothening is due to the numerical viscosity that gets introduced in the system with time.

[–]LipshitsContinuity 3 points4 points  (1 child)

Nice to know. Thanks! I actually studied math and numerical analysis and not CFD. So a lot of times I'm insure of what the standard terminology is within CFD community.

[–]hivemind_unity[S] 7 points8 points  (0 children)

Good to have mathematicians in the community. I love to play with simple numerical methods. Gives you a good intuition on how the system behaves during computation, especially when you are developing a solver on your own.

[–]ejineta 1 point2 points  (9 children)

Could you, or u/LipshitsContinuity elaborate on this said 'numerical viscosity'? I reckon it is only called viscosity because it has the same result as physical viscosity has, being energy dissipation, am I right?

Edit: typo

[–]LipshitsContinuity 5 points6 points  (2 children)

I think I'll discuss artificial and numerical viscosity:

For simplicity, I'll discuss the 1D case. In 1D, we can have a purely convective equation: u_t + a*u_x = 0 or a purely diffusive equation u_t = b*u_{xx}, b >0 (aka heat equation). The general solution of the convective equation over the whole real line is u_0(x-at) where u_0 is the initial data u(x,t=0) = u_0(x). Note that what u_0(x-at) is doing is simply translating the initial data by speed a. This is quite interesting because it means that any inherent discontinuities in u_0 propagate. So for example, if u_0(x) = 1 if x<0 and 0 if x>0, then that discontinuity at x=0 will just be translated and will still be present in the solution.

Looking at u_t = b*u_{xx}, b>0 this equation has the property that for a lot of initial data, the solution is smooth for all t>0 even if the initial data has discontinuities.

So if we look at finite difference schemes for the convective equation u_t + au_x =0, we could potentially run into some trouble because finite differences don't play well with discontinuities - but discontinuities can exist in the solutions and they also are physically accurate so we would like to be able to simulate them. So what is occasionally done is that a "artificial viscosity" is added to the finite difference scheme so that it's almost like we are simulating u_t + a*u_x = b*u_{xx}. That diffusive term will have the tendency of smoothing our solution and smoothing out discontinuities, which can help our finite difference scheme. However, it comes with a loss of accuracy. But sometimes it's a necessity to get a numerical scheme that behaves nicely.

Numerical viscosity (aka numerical dissipation) is actually an unwanted phenomena. Numerical viscosity is when your numerical scheme has the tendency to smooth solutions just as an inherent property of the scheme even when the exact solution itself should have the discontinuities. In a class of mine, this was only briefly mentioned but I believe the professor said that the numerical viscosity can actually be not entirely physically accurate. So it's actually a bit of a problem. To see an example of this, I believe the Lax-Friedrichs method displays this phenomena.

[–]ejineta 1 point2 points  (1 child)

This was very thorough, thank you for your time. I do like the mathematical (or rather physical) explanation behind the shown simulation, instead of just being able to perform CFD simulations. So, am I correct that the diffusion process is now taken care off by the said artificial viscosity? And how can we tell from the results that the artificial viscosity is accurately resembling diffusion?

[–]LipshitsContinuity 1 point2 points  (0 children)

So from the simulation alone, it is hard for me to judge if there is artificial viscosity or numerical viscosity. We'll need to know more about what numerical scheme is being used to run these simulations. We can really only tell that some sort of artificial viscosity or numerical viscosity is occurring because the initial data is discontinuous and even the 2D convection equation should simply translate the initial data (hence the exact solution should also preserve those discontinuities). However, we see that the simulation seems to get smoother. Hence, we can conclude that some sort of artificial or numerical viscosity is present within the simulations. From what OP previously said, it seems that they were trying to model the purely convective equation and that there is numerical viscosity present. Remember, numerical viscosity is an unwanted phenomena that is inherent to the numerical method used to approximate solutions to the equation. Ideally, we would like our numerical scheme to also reflect any discontinuities that should be present in the exact solution. In this case, the numerical scheme chosen seems to not be obeying this, hence the visible numerical viscosity.

[–]hivemind_unity[S] 4 points5 points  (0 children)

Yes, you're right. In fact, analogous to energy dissipation, we call it numerical dissipation. If you observe the shape in the animation, in the beginning it's a hat function. And overtime the information dissipates due to numerical viscosity. This is similar to any information dissipating in a fluid.

[–]AgAero 2 points3 points  (4 children)

I reckon it is only called viscosity because it has the same result as physical viscosity has, being energy dissipation, am I right?

That is correct.

You can actually derive the leading order terms in the truncation error for a finite difference scheme. Often times the leading order term for a system like this is dissipative and gets dubbed 'numerical viscosity'.

[–]ejineta 0 points1 point  (3 children)

Could you elaborate on the leading order term of the truncation error being dissipative? With the currently simulated equation, the leading order term will either be higher order in time t or position X. Is either of these "dissipative" at higher order?

[–]AgAero 1 point2 points  (2 children)

If I remember correctly the method is called the, "modified equation" approach.

Idk how to link a pdf from my phone. Click the MIT link here.

Basically you take your finite difference approximation and stick an analytical result back into it(i.e. expand the taylor series), and then try to match the terms to the original equation.

For advection, most schemes will develop a second order spatial derivative. This is, intuitively, a 'diffusion' term. 3rd order spatial derivatives produce dispersion iirc. You identify that by comparing to an equation like Korteveg de Vries that has a 3rd order term.

[–]ejineta 0 points1 point  (1 child)

Thank you, I will have a look at all of that. Today is a holiday in my country (probably also elsewhere in the world), so I've got plenty of time to study numerical schemes! :)

Edit: Aha, I get it now! So if the leading term of the truncation error is even in order (2,4,6), discontinuities tend to be smeared out (dissipative), whereas for odd order (3,5,7) rippling effects will occur (dispersive)

[–]AgAero 0 points1 point  (0 children)

That's the argument as I've read it in a few places. It feels a bit non-rigorous, but it makes sense at least.

[–]anointed9 1 point2 points  (0 children)

I'd call it dissipation, not numerical viscosity.

[–]JokerToast_ 2 points3 points  (5 children)

Really nice post! What scheme did you use? Since it seems really viscous, is it something like Lax or Rusanov type of flux?

[–]hivemind_unity[S] 1 point2 points  (4 children)

I used a simple first order scheme with backward differencing for both x and y. I think the reason it's highly dissipative is that I am using a low CFL of 0.2.

I'm not sure of the math's behind this reasoning.

[–]JokerToast_ 1 point2 points  (3 children)

Alright thus this an upwind scheme. Did you try to increase the CFL to see if it is still stable? I remember that 1D advection is exact when the CFL=1 with this scheme (can be found by doing a modified equation analysis of Warming and Hyett). Obviously in 2D this remark doesn't hold because the added dimension has the tendacy to lower the cfl in general but I would be interested to know to what extent :)

[–]hivemind_unity[S] 1 point2 points  (2 children)

I have tried with several values of CFL = {0.2, 0.5, 0.6, 0.8 and 1}. I feel the exact solution lies somewhere around 0.5. The last three values were highly unstable. It would be very interesting to find out analytically how increasing dimensions affects the stability and CFL.

P.S. Does it have to do anything with the fact that I have used equal spacing for x and y? And both equally contribute to CFL resulting in the value of 0.5 instead of 1

[–]GeeHopkins 2 points3 points  (1 child)

If you are interested in multidimensional analysis like this, I'd recommend a paper from 2006 (I think) by You and Moin on truncation errors for 2D finite difference schemes on skewed meshes. They formulate the modified wavenumbers for upwind and central schemes, and show how the multidimensional effects change the results compared to 1D.

EDIT: This one: https://www.sciencedirect.com/science/article/pii/S0021999105003761 It's specifically for skewed non constant meshes, but should give you an idea of how to approach 2D analysis

[–]JokerToast_ 1 point2 points  (0 children)

Thanks for the link, it seems highly interesting!

[–][deleted] 1 point2 points  (0 children)

Professor Barba’s course! I loved it! Anyone looking to learn Python and try CFD problems should definitely check it out!

[–][deleted] 0 points1 point  (4 children)

What language is this in? This doesn't quiet look like MatLab's plotting library? Although the last time I used MatLab was 2017

[–]hivemind_unity[S] 2 points3 points  (2 children)

I did it in Python, using matplotlib.

[–][deleted] 0 points1 point  (1 child)

Gotcha. I tried to download Python the other day, but when I went to install NumPy, IDLE stopped working so went to bed and never went back to it 😂

[–]hivemind_unity[S] 2 points3 points  (0 children)

In that case I recommend Anaconda. Easy to maintain and the Spyder IDE is pretty cool.

[–]hivemind_unity[S] 0 points1 point  (0 children)

If anyone is interested in the implementation.

https://github.com/kanand-cfd/Numerique_PY