This video is inspired by Lorena Barba's amazing course on CFD in Python that you can access either here: lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/ or here: github.com/barbagroup/CFDPython BTW: Did you know that her code was featured in Nvidia's recent keynote event: ruclips.net/video/jhDiaUL_RaM/видео.html where they showed this code scaled really well with their drop-in replacement for NumPy accelerating CFD to multiple GPU nodes at the same time. Just mind-blowing!
You're very welcome 😊 that's great to hear. Definitely, doing CFD computations in python has a different application/interest than doing them in established and sophisticated tools like openfoam by being more experimental or simpler, but I find them extremely insightful. Especially for someone trying to learn the concepts.
You're very welcome 🤗 Thanks for the suggestion. It's more of a nichy topics. For now, I want to keep the videos rather general to address a larger audience
Could you make a video on this example but in 3d and with no slip boundary conditions on only five of the walls, (the wall with slip condition is the one where the 'stream' is that causes the swirls.) I really love this video and the pressure displays.
Hi, thanks for commenting. :) That would indeed be a great extension to 3D. Though, I am unsure whether it would be worthwhile to implement similarly to this 2D video due to a considerable increase in computational requirements. Maybe that would be a great task for the Julia Programming language. I will note down your idea, but I can't promise to have a video in the near future :) Of course, if you come up with a 3D extension, feel free to open up a pull request on the GitHub repo: github.com/Ceyron/machine-learning-and-simulation
@@hasnaouiacademy7899 Will definitely come in the far future ;) As of now, I think Python or Julia is better because the overhead and boilerplate is lower
@@MachineLearningSimulation it's just to get some knowledge about how to programme with CFD. Do you know how to work with machine learning to solve PDE equations. I am very interested in this field.
@@hasnaouiacademy7899 There are a lot of flavors to using ML for PDE solving. You can have full surrogates like PINNs or just data-driven models like turbulence closures etc. Maybe a good first pointer is this review article: www.annualreviews.org/doi/abs/10.1146/annurev-fluid-010719-060214
You're very welcome 🤗 FVM is definitely on my agenda, at the moment I am focusing on methods that are more relevant to my PhD: finite differences and spectral methods. But I will definitely come back to other approaches in the future. In 1d and on structured grids FVM is also very similar to FDM. The projection method (in the sense of making a flow field incompressible) is similar for all discretization methods in that you end up solving a pressure poisson problem.
Hi, thx a lot for this tutorial, very helpful ^^. I am new to CFD and I came across the existence of different type of discretization schemes (upwind, linear upwind, and central). I guess u are making use of the central one in this video. Also, why did u choose this one specifically? I guess in such a simple case it would not make any difference, but I suppose there are scenarios in which one gives better results than the others. Are there some rules depending on the domain structure? Thx again!
Hi, thanks for the kind feedback 😊 I'm happy the video was of help to you. Indeed, you correctly observed that the convection/advection part of the Navier Stokes equations was discretized using a central scheme, in this case the second order accurate central scheme for the first derivative. For many scenarios, this is a bad choice. Here it is only viable for making the tutorial simpler. If you look at the prototypical hyperbolic PDE: the 1D advection equation and you discretizatized it with an ftcs (forward-in-time central-in-space) scheme (like we did in this video here), you would even find that it is unconditionally unstable. In other words: you can make your time step as small as possible and the simulation will still be unstable (it develops unphysical oscillations). The usual remedy is the use of a first order upwind scheme. If your advection speed is postive, you discretize the first derivative with a backward stencil. Vice versa, if the advection speed was negative, you would use a forward stencil. The reason, we are able to make usage of a central scheme in this video was because we have enough diffusion. In other words, the kinematic viscosity was chosen high enough. There is also a dimensionless number for this: the Peclet number (closely related to the Reynolds number). We know that for an advection-diffusion equation a central scheme for the advection part is only stable if we can guarantee that the Peclet number in each part of the domain is smaller than 2. This is the scenario for which we are simulating here. In general, for any realistic CFD, you will mostly rely on upwinding biased scheme. This is a general term and just means that your stencil points for the derivative approximation are not centered around your desired degree of freedom, but favor one side over the over, which side is typically dependent on the flow direction (similar to the first-order upwind). In the CFD community (especially for openfoam) I believe that the word "upwind" just refers to the first order upwind, as described before whereas "Liner upwind" refers to a second order upwinding. (Despite the latter being more desirable it has also some disadvantages). Hope that helped 😊 Let me know if something is still unclear.
which editor are you using? i am using spyder and i absolutely hate how slow it performs on my computer also i have 3d ascii cfd data generated from a fortran code. is there a way to plot it on python? (x, y, z, pressure)
Hi, I am using VS Code together with the Python extension. I like this combination a lot and it performs reasonably fast. Regarding your second question: There are probably multiple approaches. If your data comes in a standard format like CSV, there are a lot of reader libraries that can read the files and, for instance, create numpy arrays out of it. Maybe check out this tutorial: numpy.org/devdocs/user/how-to-io.html It can also be helpful to look into Pandas.
Thank you very much for the video. I am wondering if you may teach me how to adjust the boundary conditions such that the flow goes around a cylinder. Thank you.
You're very welcome 🤗 thanks for the comment. I believe you are interested in a setting different from this example in that you also have an inflow and an outflow. In the example of the lid driven cavity, there domain is actually closed. Take a look at the other videos in the series: 🌊 Simulations simply implemented in Python or Julia: ruclips.net/p/PLISXH-iEM4JmgBfU_QU262MQTYa7DoJK0 For example the one with the transient pipe flow that shows how to prescribe inflow and outflow. I also have a video on the van karman vortex street (flow behind a cylinder), but it uses an lbm approach. For a FDM approach, the easiest you could do is to create a mask for the nodes that belong to the cylinder then force zero Dirichlet BCs on them for the velocity and and zero neumann on the pressure.
Sure, you can manipulate the kinematic viscosity. That should result in a change in the Reynold's number ;) But be aware of the possible stability issues I noted in the pinned comment.
Certainly, some form of residual check would be necessary to judge about the quality of the solver. Though, I think it's a bit beyond the purpose of the video as this just aims to be the simplest (arguably also the crudest) way to solve the NS equations 😅 There will most likely be a huge residual in the continuity as the pressure poisson problem is just solved with a Jacobi smoother for a fixed number of iterations, no convergence check is made. Also, for the momentum only one Picard iteration is performed, also without convergence check. Would be interesting to see how that affects the long-term solution since it converges to a steady-state 😁 Long story short: if you have a nice way of adding a residual print-out feel free to open a PR over on GitHub :)
Hello, after updating the velocities (applying the pressure correction), don't we need to check if du/dx + dv/dy = 0, and if it is not, keep repeating the process for that time step till it is?
Hi, great question! :) It depends on the context of what you are simulating for. If you do a true engineering CFD simulation, you have to ensure both global and local continuity (at least to within a sufficient tolerance). As you mention, one way would be to repeat solving the pressure correction problem with subsequent projection until this tolerance is reached. One should then also put emphasis on having the correct BCs for the poisson problem (here we incorrectly used homogeneous Dirichlet BC at the top for convenience) as well as employing a better linear solver: the Jacobi solver with a fixed amount of iterations is arguably a bag choice. There are many approaches as to how to select outer, inner and mid-iterations giving, e.g., rise to algorithms commonly used in engineering like SIMPLE or PISO. I can totally recommend the book by Stefan Turek called "efficient solvers for incomprehensible flow problems" as it unifies many of the approaches from the perspective of solving a coupled saddle point system. Initially, I said it depends because even a bad incompressibility projection (i.e. one that does not reduce the divergence enough) is still creating swirly motion. We want a vortex to get the lid driven cavity. You can also try this yourself, by commenting out the projection step, the flow will not get swirly. Hope that cleared it up a bit :). As a final remark: the intent of the video is educational, it's not a rock solid CFD code. I don't think you can write such a code in this limited amount of time 😁 Let me know if sth is unclear :)
@@MachineLearningSimulation Thanks for the detailed reply! Makes sense. I was just asking because I was running similar problems with a different method and a slightly different approach, however I was getting a program that is very slow. I was trying to compare what the difference was between the two programs, and figured that it was the incompressibility check, so I just wanted to know if there was a way to get around it.
@@atharvnaik that's an interesting observation. 😊 My experience is that the pressure correction problems are often the toughest parts in each CFD simulation (not only from the analytical perspective with deriving the correct BCs, but especially for computing them numerically). It commonly becomes the part of the code the most computational time and effort is spent. In other words: with taking a more lax approach as in the video we can simplify many things and "speed them up" 😁
In my understanding, the Courant number (or CFL number) is a dimensionless quantity defining how far information propagates within the next time step. It is therefore defined as u_max * dt / dx. It is the stability condition of the Upwind Finite Differencing scheme when analyzed with von Neumann stability analysis. However, in the video we used a stability constraint for the FTCS scheme of the diffusion equation because we do not have a balanced discretization in space and time. Hence, I would say the introduced "stability safety factor" serves a similar purpose as the CFL number, but is not the number. What do you think?
@@MachineLearningSimulation without a doubt, it's not the Courant number... I referred to it because that was the interpretation I gave for the solution you presented, regarding the issue of stability. So I figured it was something similar. I understand perfectly.
Important note, that I forgot in the video: Towards the end, we implemented a stability criterion for the diffusive part of the equation. Be aware that this does not mean you can select the kinematic viscosity to be arbitrarily small. If you do so, the (local) Peclet Number might be larger than 2 (en.wikipedia.org/wiki/P%C3%A9clet_number) leading to the unconditional instability of the Forward-Time-Central-Space discretization of the nonlinear convection. This is because the FTCS scheme only works for advection-diffusion problems (like the Navier-Stokes equations are) if diffusive transport is dominant. If you want to avoid this problem, you need to discretize the nonlinear convection term in means of an Upwind Difference. If you choose to use the First Order Upwind, then the time step has to be chosen in accordance to the CFL condition (en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) which can be a tougher restriction than the constraint we have from the FTCS scheme. In order to enforce the CFL condition, you might also have to adapt the time step throughout the simulation, as it then depends on the maximum velocity in the domain.
Additionally, be aware that at even lower kinematic viscosities you might see spurious oscillations in the pressure field. This is the so-called checkerboard pattern that is caused by using a central difference for the pressure gradient on a co-located grid. A co-located grid means that all variables (u, v and p) are stored at the same points in the domain. There are two remedies: 1) Use a staggered grid. This will have pressure, u velocity and v velocity being saved at different points in the mesh. 2) Use a corrective interpolation, e.g. the one proposed by Rhie & Chow.
How come I get this when running the code? RuntimeWarning: invalid value encountered in add u_prev * d_u_prev__d_x + v_prev * d_u_prev__d_y RuntimeWarning: overflow encountered in multiply
That would probably refer to some of your values turning NaN (not a number) during the simulation. Often, this happens when the computation becomes instable. Maybe you have changed some of the parameters (like number of discretization points, kinematic viscosity or the time step size). Try to lower the time step size ;)
Is there any source that I can receive detailed explanation regarding the central difference scheme that is provided in 8.00 (row 157)? Actually I could not umderstand that part clearly.
Sure, I think you can find more details on that in introductory books go CFD like the one by Ferziger and Peric or the one by Versteeg. :) It is a second order approximation to the first derivative by (u_{i+1,j} - u_{i-1,j})/(2 * dx). Just Note that due to the array layout the "i-indexing" is in the column dimension of the matrix, and not in the row dimension as one would expect by lexicographical ordering. Hope that helps :). Let me know if that is still unclear.
Very interesting and very well explained! That's the kind of video that helps students like me. Simple question: does the fluid inside the cavity can "flow out"? What if the "lid" is another fluid, is it possible to modify your code to compute the mixing that will occur?
Thanks for the kind feedback 😊 There is a minor mistake in the video, in that I said that the top boundary is something like the interface to a free stream that is passing by. Actually, the lid-driven cavity scenario is closed off; three walls are stationary, and one (the top) is moving to the right. But back to your question: if the top was open and a free stream with a different fluid passed by: generally speaking, yes it is possible to simulate the mixture of two fluids. However, this is way beyond the complexity of this video. There are approaches like the volume of fluid method but I have to refer you to a text book on intermediate topics in CFD for that (maybe it is covered in fertiger & Peric). If one had the simpler case of a passive species being transported (like smoke) in the liquid, it is possible to do so by a simple additional advection equation (that is easier to solve than the NS equations). There is also the option that the smoke asserts a force back into the liquid based on its buoyancy. I also have video on that for which I used the phiflow package: ruclips.net/video/KMfcF9XvVio/видео.htmlsi=Y065eoSm6WiBG-jr Hope that answered your question. 😊
Hello, I was working on a problem of comparable difficulty, however I have used loops for everything instead of directly using the arrays. The code takes ages to run if I try making the grid finer ( splitting it into even 10 parts for a length of 1 takes more than 40 mins). Any tips on how to speed it up?
Hi, there are a lot of factors. Generally speaking, expressing operations in (double-) for loops in Python is a common slowdown factor. It could yield a significant speedup to translate your operations into efficient NumPy array calls as that avoids the notoriously slow Python looping. However, that depends a bit on the size of problems. My rule of thumb is that NumPy array operations tend to be relevant for everything with more than 200 or 300 degrees of freedom, which in 2D would already be achieved with a ~ 15x15 discretization. Usually, the most time-critical aspect in incompressible Navier-Stokes simulations is the elliptic pressure correction problem. In this video, I use a fixed number of Jacobi iterations. Arguably, that's a terrible way to solve a linear system, as I am guaranteed to run into troubles at finer discretizations (due to the growing condition number of the problem). Maybe, if you are using an off-the-shelve linear solver (like Conjugate gradient), it could be helpful to reduce its tolerances such that it requires fewer costly iterations to converge. For larger problems (let's say anything larger than 300x300 you might also want to use a preconditioner - one of the most efficient preconditioner for elliptic problems are multigrid methods). Maybe also try translating your code to Julia :) Good luck with the fantastic world of PDE solvers. :)
@@atharvnaik Thanks for the kind feedback :). I am happy that the videos are of help. It also took me quite a while to get my head around it. I try to make the videos in a style I would have loved to had myself, i.e., to use a top-down approach by starting code-first and then investigate all the nitty details.
Hey, thanks for the kind words, :) Really appreciated. Regarding your question: I think it would not be too difficult to extend to 3D. The Laplace operator would get a third term, and also in the Central Differences discretization of the convection we would have a third one. As long as you keep the degrees of freedom reasonably low, it should also run in an acceptable time. Of course, make sure the simulation is stable. It might be a cool challenge. If you want to share something: feel free to open a Pull Request over on the channel's GitHub repo: github.com/Ceyron/machine-learning-and-simulation Regarding non-isothermal: I would say that depends on the concrete modeling. In the easiest case, you just have to solve another advection-diffusion equation in each time iteration. That's a cool idea for a future video. I will note that down. ;)
Hi, thanks for the question. Unfortunately, I do not have much experience yet on working with Navier-Stokes in cylindrical coordinates. Though, more generally speaking, projection methods (at least P1 projection methods like Chorin's splitting) work such that you ignore the pressure gradient term in the momentum equation, then solve an elliptic PDE for a (pseudo-) pressure correction and then correct the velocities. I assume this might be similarly applicable to cylindrical coordinates. I think this Wikipedia Article might be a good first start: en.wikipedia.org/wiki/Projection_method_(fluid_dynamics) What is the concrete problem you are solving?
@@MachineLearningSimulation Thank you, will try again keeping this in mind. Thanks a lot for the wonderful video, found it helpful to get started with code. I am trying to solve startup laminar flow of liquid in a pipe using projection method till the fluid reaches steady state flow. Initial velocity of the flow is given.
Where does the external flow come from? Is the lid being moved and therefore causing flow? Or are you programming some arbitrary flow to force vorticity?
Hi, that was an error in my understanding of the lid-driven cavity at that time. Correct would be that all four boundaries are walls, and the top wall is moving to the right. Hence, there should be no "external flow". Though, I think both are mathematically similar since we have a prescribed Dirichlet boundary on the velocity. Probably, the pressure handling changes. Do you know more?
Why not apply the boundary conditions to the previous velocity? Now, the very first timestep iteration, the calculation of the tentative velocity does nothing but create another zero matrix. Seems inefficient.
Well, yes, there is always room for improvements in computational efficiency, but that's not the main goal of the video ;). I believe that what you propose is just the tip of the iceberg. The involved Jacobi solver for the pressure Poisson problem is the biggest weak-point of the implementation. There are many other linear solves that are mathematically more efficient (by a large amount). You find many more such things in the implementation :D.
@@MachineLearningSimulation So it is the Jacobi solver? Good to know, because I have been trying to figure out how this algorithm would converge, but now I know the name.
I have looked into it and cannot figure out how you solved the pressure Poisson equation. From what I understand, you are solving the laplace(p_prev) - rho/delta t * gradient(u) = 0. in every interior position in p_next, you write laplace(p_prev) - rho/delta t * gradient(u). But that is formula, not an equation. At what point is laplace(p_prev) - rho/delta t * gradient(u) = 0 actually solved? What is the x value you are looking for? I guess the real question is: This video gives the equation you solve when using the Jacobi iteration method ruclips.net/video/VH0TZlkZPRo/видео.html&ab_channel=EMPossible I can't figure out how this codes translate into solving this equation. Which part is the vector b, which is the matrix A and what is x in the code. Could you please elaborate?
@@MachineLearningSimulation Hi! To my knowledge, we need an additional padding of grid on every sides to implement the boundary conditions. But the code doesn't have any. Please let me know about this.
I think you are referring to a cell-centered discretization which is common when using finite volumes. Here, however, we have a vertex-centered discretization, mostly (but necessarily) due to usage of Finite differences. Hence, we have degrees of freedom that lie directly on the boundary which allows us to prescribe both Dirichlet BCs for the Velocity and Neumann BCs for the pressure. The latter are unfortunately only of first order consistency. Is that what you were referring to? ☺️
@@MachineLearningSimulation Yes, this is what I am referring to. I thought that for pressure projection, we need to discretize using approach similar to finite volume, with pressure being cell centered (same as your pipe flow video), otherwise it might create pressure oscillations.
Sir, can you please explain why there is a kind of phase separation in the streamline plot???? Will it persist in this particular problem of lid driven cavity throughout the run time ???
@@MachineLearningSimulation Sir , actually I just want to know the explanation of what that different colour layers does mean which are appearing at time of simulation ??? Physical explanation
Do you refer to the central difference approximation of the Laplace operator or of the first derivative? :) In case of the Laplace operator, "-4" in front of the central point should be correct since we are taking a second order approximation of the second derivative in both x and y-axis. Hence, the stencil is 0 1 0 1 -4 1 0 1 0 For the central difference approximation of the first derivative (which we used to evaluate the convection/advection), the second order stencil does not include the point itself. Hence, it is -1 0 1 (wrt x) or 1 0 -1 (wrt y) What do you think? Btw: the stability check is not for the CFL condition, but for the stability of the FTCS scheme of the heat equation (en.wikipedia.org/wiki/FTCS_scheme#Stability ). This is because we can only simulate diffusion dominant flows with explicit central difference approximations to the advection. For higher Reynolds numbers, one would observe pressure oscillations. For these high diffusion cases (as simulated here), the FTCS stability criterion is often stricter. Hope I got you right :). There is of course a chance that I made a mistake. Please let me know, then I can fix the file on GitHub and leave a remark under the video.
I think that there are definitely resemblances to what the original PISO was doing, but I would not coin the approach here PISO. More generally speaking, the PISO algorithm (and also the SIMPLE algorithm) are just some special treatments of the coupled NS problem, see e.g. "Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach" by Stefan Turek. The approach done here is inspired by Lorena Barba's "12 steps to Navier Stokes course", but uses the simpler incompressibility projection by Helmholtz decomposition.
Hi, it seems to me that your simulation turned unstable. As such, some (and eventually all) your floating point values in the arrays turned NaN (not a number) which what you're error messages are telling you. This is caused by a too high time step. If you lower it again, it should be fine. The simulation here is only conditionally stable since we are using both an explicit convection and an explicit diffusion treatment. Hence, there will be some condition for the time step length in relation to the spatial discretization length (and some other quantities) that has to be fulfilled.
![State of Play | September 24, 2024 [ENGLISH]](http://i.ytimg.com/vi/fVs0suPHAX8/mqdefault.jpg)
This video is inspired by Lorena Barba's amazing course on CFD in Python that you can access either here: lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/ or here: github.com/barbagroup/CFDPython
BTW: Did you know that her code was featured in Nvidia's recent keynote event: ruclips.net/video/jhDiaUL_RaM/видео.html where they showed this code scaled really well with their drop-in replacement for NumPy accelerating CFD to multiple GPU nodes at the same time. Just mind-blowing!
Do you know which minute of Nvidia's recent keynote event? I couldn't find it in the "Accelerating Science by a Million-X" or any other.
@@josealvarezcabrera561 It's briefly mentioned at 25:46 in ruclips.net/video/jhDiaUL_RaM/видео.html
Thank you for these videos! I am trying to improve at python (particularly for CFD) after mostly using openfoam, and these have been perfect.
You're very welcome 😊 that's great to hear.
Definitely, doing CFD computations in python has a different application/interest than doing them in established and sophisticated tools like openfoam by being more experimental or simpler, but I find them extremely insightful. Especially for someone trying to learn the concepts.
Brilliant tutorial, man. Such a clear explanation
Thanks for the feedback :). Much appreciated!
Thank you for this wonderful video.
Please can you also do another one for a vertical flow for two phase liquid and gas
You're very welcome 🤗
Thanks for the suggestion. It's more of a nichy topics. For now, I want to keep the videos rather general to address a larger audience
I highly appreciate your efforts on these excellent videos! Please keep it up!
Appreciate the generous donations 😊 Thanks a lot.
More videos on PDEs (especially their interplay with) neural networks will come in the future. :)
Could you make a video on this example but in 3d and with no slip boundary conditions on only five of the walls, (the wall with slip condition is the one where the 'stream' is that causes the swirls.) I really love this video and the pressure displays.
Hi,
thanks for commenting. :)
That would indeed be a great extension to 3D. Though, I am unsure whether it would be worthwhile to implement similarly to this 2D video due to a considerable increase in computational requirements. Maybe that would be a great task for the Julia Programming language.
I will note down your idea, but I can't promise to have a video in the near future :)
Of course, if you come up with a 3D extension, feel free to open up a pull request on the GitHub repo: github.com/Ceyron/machine-learning-and-simulation
Nice Explanation! Which text editor are you using btw?
You're welcome 😊
The editor is VS code.
@@MachineLearningSimulation what font do you use?
@@josetjaw8161 Should be the deault theme for VS code under Linux.
Great work, I like it, especially when you use tqdm
Thanks. ❤️
Yes, I think tqdm is a neat tool to show the transient simulation progress.
@@MachineLearningSimulation please if you can make a simulation with c++
thank's
@@hasnaouiacademy7899 Will definitely come in the far future ;)
As of now, I think Python or Julia is better because the overhead and boilerplate is lower
@@MachineLearningSimulation it's just to get some knowledge about how to programme with CFD.
Do you know how to work with machine learning to solve PDE equations.
I am very interested in this field.
@@hasnaouiacademy7899 There are a lot of flavors to using ML for PDE solving. You can have full surrogates like PINNs or just data-driven models like turbulence closures etc. Maybe a good first pointer is this review article: www.annualreviews.org/doi/abs/10.1146/annurev-fluid-010719-060214
Thank you .
Can u make a tutorial for FVM , projection method ? That would be so helpful..
You're very welcome 🤗
FVM is definitely on my agenda, at the moment I am focusing on methods that are more relevant to my PhD: finite differences and spectral methods. But I will definitely come back to other approaches in the future. In 1d and on structured grids FVM is also very similar to FDM. The projection method (in the sense of making a flow field incompressible) is similar for all discretization methods in that you end up solving a pressure poisson problem.
Thank You for the explanation.
You're very welcome. :)
Hi, thx a lot for this tutorial, very helpful ^^.
I am new to CFD and I came across the existence of different type of discretization schemes (upwind, linear upwind, and central). I guess u are making use of the central one in this video. Also, why did u choose this one specifically?
I guess in such a simple case it would not make any difference, but I suppose there are scenarios in which one gives better results than the others. Are there some rules depending on the domain structure? Thx again!
Hi,
thanks for the kind feedback 😊 I'm happy the video was of help to you.
Indeed, you correctly observed that the convection/advection part of the Navier Stokes equations was discretized using a central scheme, in this case the second order accurate central scheme for the first derivative. For many scenarios, this is a bad choice. Here it is only viable for making the tutorial simpler.
If you look at the prototypical hyperbolic PDE: the 1D advection equation and you discretizatized it with an ftcs (forward-in-time central-in-space) scheme (like we did in this video here), you would even find that it is unconditionally unstable. In other words: you can make your time step as small as possible and the simulation will still be unstable (it develops unphysical oscillations).
The usual remedy is the use of a first order upwind scheme. If your advection speed is postive, you discretize the first derivative with a backward stencil. Vice versa, if the advection speed was negative, you would use a forward stencil.
The reason, we are able to make usage of a central scheme in this video was because we have enough diffusion. In other words, the kinematic viscosity was chosen high enough. There is also a dimensionless number for this: the Peclet number (closely related to the Reynolds number). We know that for an advection-diffusion equation a central scheme for the advection part is only stable if we can guarantee that the Peclet number in each part of the domain is smaller than 2. This is the scenario for which we are simulating here.
In general, for any realistic CFD, you will mostly rely on upwinding biased scheme. This is a general term and just means that your stencil points for the derivative approximation are not centered around your desired degree of freedom, but favor one side over the over, which side is typically dependent on the flow direction (similar to the first-order upwind).
In the CFD community (especially for openfoam) I believe that the word "upwind" just refers to the first order upwind, as described before whereas "Liner upwind" refers to a second order upwinding.
(Despite the latter being more desirable it has also some disadvantages).
Hope that helped 😊 Let me know if something is still unclear.
which editor are you using? i am using spyder and i absolutely hate how slow it performs on my computer
also i have 3d ascii cfd data generated from a fortran code. is there a way to plot it on python? (x, y, z, pressure)
Hi,
I am using VS Code together with the Python extension. I like this combination a lot and it performs reasonably fast.
Regarding your second question: There are probably multiple approaches. If your data comes in a standard format like CSV, there are a lot of reader libraries that can read the files and, for instance, create numpy arrays out of it. Maybe check out this tutorial: numpy.org/devdocs/user/how-to-io.html
It can also be helpful to look into Pandas.
Thank you very much for the video. I am wondering if you may teach me how to adjust the boundary conditions such that the flow goes around a cylinder. Thank you.
You're very welcome 🤗 thanks for the comment.
I believe you are interested in a setting different from this example in that you also have an inflow and an outflow. In the example of the lid driven cavity, there domain is actually closed. Take a look at the other videos in the series: 🌊 Simulations simply implemented in Python or Julia: ruclips.net/p/PLISXH-iEM4JmgBfU_QU262MQTYa7DoJK0
For example the one with the transient pipe flow that shows how to prescribe inflow and outflow. I also have a video on the van karman vortex street (flow behind a cylinder), but it uses an lbm approach. For a FDM approach, the easiest you could do is to create a mask for the nodes that belong to the cylinder then force zero Dirichlet BCs on them for the velocity and and zero neumann on the pressure.
@@MachineLearningSimulation Thank you so so much!!!!
You're welcome:)
you are excellent man
Thanks 🙏
sir can we manipulate this code for different Reynolds number?
Sure, you can manipulate the kinematic viscosity. That should result in a change in the Reynold's number ;)
But be aware of the possible stability issues I noted in the pinned comment.
What about adding residual error?
Certainly, some form of residual check would be necessary to judge about the quality of the solver. Though, I think it's a bit beyond the purpose of the video as this just aims to be the simplest (arguably also the crudest) way to solve the NS equations 😅 There will most likely be a huge residual in the continuity as the pressure poisson problem is just solved with a Jacobi smoother for a fixed number of iterations, no convergence check is made. Also, for the momentum only one Picard iteration is performed, also without convergence check. Would be interesting to see how that affects the long-term solution since it converges to a steady-state 😁
Long story short: if you have a nice way of adding a residual print-out feel free to open a PR over on GitHub :)
Hello, after updating the velocities (applying the pressure correction), don't we need to check if du/dx + dv/dy = 0, and if it is not, keep repeating the process for that time step till it is?
Hi, great question! :)
It depends on the context of what you are simulating for. If you do a true engineering CFD simulation, you have to ensure both global and local continuity (at least to within a sufficient tolerance). As you mention, one way would be to repeat solving the pressure correction problem with subsequent projection until this tolerance is reached. One should then also put emphasis on having the correct BCs for the poisson problem (here we incorrectly used homogeneous Dirichlet BC at the top for convenience) as well as employing a better linear solver: the Jacobi solver with a fixed amount of iterations is arguably a bag choice.
There are many approaches as to how to select outer, inner and mid-iterations giving, e.g., rise to algorithms commonly used in engineering like SIMPLE or PISO. I can totally recommend the book by Stefan Turek called "efficient solvers for incomprehensible flow problems" as it unifies many of the approaches from the perspective of solving a coupled saddle point system.
Initially, I said it depends because even a bad incompressibility projection (i.e. one that does not reduce the divergence enough) is still creating swirly motion. We want a vortex to get the lid driven cavity. You can also try this yourself, by commenting out the projection step, the flow will not get swirly.
Hope that cleared it up a bit :). As a final remark: the intent of the video is educational, it's not a rock solid CFD code. I don't think you can write such a code in this limited amount of time 😁
Let me know if sth is unclear :)
@@MachineLearningSimulation Thanks for the detailed reply! Makes sense. I was just asking because I was running similar problems with a different method and a slightly different approach, however I was getting a program that is very slow. I was trying to compare what the difference was between the two programs, and figured that it was the incompressibility check, so I just wanted to know if there was a way to get around it.
@@atharvnaik that's an interesting observation. 😊
My experience is that the pressure correction problems are often the toughest parts in each CFD simulation (not only from the analytical perspective with deriving the correct BCs, but especially for computing them numerically). It commonly becomes the part of the code the most computational time and effort is spent.
In other words: with taking a more lax approach as in the video we can simplify many things and "speed them up" 😁
Can you provide the link for the next video you were talking about in the start of the video?
Hey,
it's unfortunately not released yet (not even produced). You can expect it later this year.
26:28... like Courant Number?
In my understanding, the Courant number (or CFL number) is a dimensionless quantity defining how far information propagates within the next time step. It is therefore defined as u_max * dt / dx. It is the stability condition of the Upwind Finite Differencing scheme when analyzed with von Neumann stability analysis.
However, in the video we used a stability constraint for the FTCS scheme of the diffusion equation because we do not have a balanced discretization in space and time. Hence, I would say the introduced "stability safety factor" serves a similar purpose as the CFL number, but is not the number.
What do you think?
@@MachineLearningSimulation without a doubt, it's not the Courant number... I referred to it because that was the interpretation I gave for the solution you presented, regarding the issue of stability. So I figured it was something similar. I understand perfectly.
Mister, these are incredible videos. Keep up the good work. Best Regards from Germany.
Thank so much ❤️
I didn't understood shouldn't u_tent [0, :] be top boundary? Row indexing starting from upper row? 13:21
Important note, that I forgot in the video: Towards the end, we implemented a stability criterion for the diffusive part of the equation. Be aware that this does not mean you can select the kinematic viscosity to be arbitrarily small. If you do so, the (local) Peclet Number might be larger than 2 (en.wikipedia.org/wiki/P%C3%A9clet_number) leading to the unconditional instability of the Forward-Time-Central-Space discretization of the nonlinear convection.
This is because the FTCS scheme only works for advection-diffusion problems (like the Navier-Stokes equations are) if diffusive transport is dominant.
If you want to avoid this problem, you need to discretize the nonlinear convection term in means of an Upwind Difference. If you choose to use the First Order Upwind, then the time step has to be chosen in accordance to the CFL condition (en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) which can be a tougher restriction than the constraint we have from the FTCS scheme. In order to enforce the CFL condition, you might also have to adapt the time step throughout the simulation, as it then depends on the maximum velocity in the domain.
Additionally, be aware that at even lower kinematic viscosities you might see spurious oscillations in the pressure field. This is the so-called checkerboard pattern that is caused by using a central difference for the pressure gradient on a co-located grid. A co-located grid means that all variables (u, v and p) are stored at the same points in the domain.
There are two remedies:
1) Use a staggered grid. This will have pressure, u velocity and v velocity being saved at different points in the mesh.
2) Use a corrective interpolation, e.g. the one proposed by Rhie & Chow.
How come I get this when running the code?
RuntimeWarning: invalid value encountered in add
u_prev * d_u_prev__d_x + v_prev * d_u_prev__d_y
RuntimeWarning: overflow encountered in multiply
That would probably refer to some of your values turning NaN (not a number) during the simulation. Often, this happens when the computation becomes instable. Maybe you have changed some of the parameters (like number of discretization points, kinematic viscosity or the time step size). Try to lower the time step size ;)
@@MachineLearningSimulation Thanks, it runs.
@@antares6998 You're welcome.
Is there any source that I can receive detailed explanation regarding the central difference scheme that is provided in 8.00 (row 157)? Actually I could not umderstand that part clearly.
Sure, I think you can find more details on that in introductory books go CFD like the one by Ferziger and Peric or the one by Versteeg. :)
It is a second order approximation to the first derivative by (u_{i+1,j} - u_{i-1,j})/(2 * dx). Just Note that due to the array layout the "i-indexing" is in the column dimension of the matrix, and not in the row dimension as one would expect by lexicographical ordering.
Hope that helps :). Let me know if that is still unclear.
Very interesting and very well explained! That's the kind of video that helps students like me.
Simple question: does the fluid inside the cavity can "flow out"? What if the "lid" is another fluid, is it possible to modify your code to compute the mixing that will occur?
Thanks for the kind feedback 😊
There is a minor mistake in the video, in that I said that the top boundary is something like the interface to a free stream that is passing by. Actually, the lid-driven cavity scenario is closed off; three walls are stationary, and one (the top) is moving to the right.
But back to your question: if the top was open and a free stream with a different fluid passed by: generally speaking, yes it is possible to simulate the mixture of two fluids. However, this is way beyond the complexity of this video. There are approaches like the volume of fluid method but I have to refer you to a text book on intermediate topics in CFD for that (maybe it is covered in fertiger & Peric).
If one had the simpler case of a passive species being transported (like smoke) in the liquid, it is possible to do so by a simple additional advection equation (that is easier to solve than the NS equations). There is also the option that the smoke asserts a force back into the liquid based on its buoyancy. I also have video on that for which I used the phiflow package: ruclips.net/video/KMfcF9XvVio/видео.htmlsi=Y065eoSm6WiBG-jr
Hope that answered your question. 😊
Nice. Adding this to my projects
Awesome video! Thank you!
Great, thanks for enjoying these simulation scripts videos :)
Hello, I was working on a problem of comparable difficulty, however I have used loops for everything instead of directly using the arrays. The code takes ages to run if I try making the grid finer ( splitting it into even 10 parts for a length of 1 takes more than 40 mins). Any tips on how to speed it up?
Hi,
there are a lot of factors. Generally speaking, expressing operations in (double-) for loops in Python is a common slowdown factor. It could yield a significant speedup to translate your operations into efficient NumPy array calls as that avoids the notoriously slow Python looping. However, that depends a bit on the size of problems. My rule of thumb is that NumPy array operations tend to be relevant for everything with more than 200 or 300 degrees of freedom, which in 2D would already be achieved with a ~ 15x15 discretization.
Usually, the most time-critical aspect in incompressible Navier-Stokes simulations is the elliptic pressure correction problem. In this video, I use a fixed number of Jacobi iterations. Arguably, that's a terrible way to solve a linear system, as I am guaranteed to run into troubles at finer discretizations (due to the growing condition number of the problem). Maybe, if you are using an off-the-shelve linear solver (like Conjugate gradient), it could be helpful to reduce its tolerances such that it requires fewer costly iterations to converge. For larger problems (let's say anything larger than 300x300 you might also want to use a preconditioner - one of the most efficient preconditioner for elliptic problems are multigrid methods).
Maybe also try translating your code to Julia :)
Good luck with the fantastic world of PDE solvers. :)
Thanks man! Your videos are really amazing for someone trying to figure out this stuff on their own
@@atharvnaik Thanks for the kind feedback :). I am happy that the videos are of help. It also took me quite a while to get my head around it. I try to make the videos in a style I would have loved to had myself, i.e., to use a top-down approach by starting code-first and then investigate all the nitty details.
How difficult would this be to expand this solution strategy to a 3D system? Or a nonisothermal system?
Also this was an awesome video. Thanks for the great explanations :).
Hey, thanks for the kind words, :) Really appreciated.
Regarding your question: I think it would not be too difficult to extend to 3D. The Laplace operator would get a third term, and also in the Central Differences discretization of the convection we would have a third one. As long as you keep the degrees of freedom reasonably low, it should also run in an acceptable time. Of course, make sure the simulation is stable. It might be a cool challenge. If you want to share something: feel free to open a Pull Request over on the channel's GitHub repo: github.com/Ceyron/machine-learning-and-simulation
Regarding non-isothermal: I would say that depends on the concrete modeling. In the easiest case, you just have to solve another advection-diffusion equation in each time iteration. That's a cool idea for a future video. I will note that down. ;)
Could you please help on how to use projection method to solve Navier-Stokes equation in cylindrical coordinates?
Hi,
thanks for the question. Unfortunately, I do not have much experience yet on working with Navier-Stokes in cylindrical coordinates. Though, more generally speaking, projection methods (at least P1 projection methods like Chorin's splitting) work such that you ignore the pressure gradient term in the momentum equation, then solve an elliptic PDE for a (pseudo-) pressure correction and then correct the velocities. I assume this might be similarly applicable to cylindrical coordinates. I think this Wikipedia Article might be a good first start: en.wikipedia.org/wiki/Projection_method_(fluid_dynamics)
What is the concrete problem you are solving?
@@MachineLearningSimulation Thank you, will try again keeping this in mind. Thanks a lot for the wonderful video, found it helpful to get started with code.
I am trying to solve startup laminar flow of liquid in a pipe using projection method till the fluid reaches steady state flow. Initial velocity of the flow is given.
You're welcome 😊
Maybe you also find this video helpful: ruclips.net/video/rV8tD2nQfkk/видео.html
Where does the external flow come from? Is the lid being moved and therefore causing flow? Or are you programming some arbitrary flow to force vorticity?
Hi,
that was an error in my understanding of the lid-driven cavity at that time. Correct would be that all four boundaries are walls, and the top wall is moving to the right. Hence, there should be no "external flow".
Though, I think both are mathematically similar since we have a prescribed Dirichlet boundary on the velocity. Probably, the pressure handling changes. Do you know more?
that is GREAT, Well Done
Thanks a lot ❤️
Why not apply the boundary conditions to the previous velocity? Now, the very first timestep iteration, the calculation of the tentative velocity does nothing but create another zero matrix. Seems inefficient.
Well, yes, there is always room for improvements in computational efficiency, but that's not the main goal of the video ;).
I believe that what you propose is just the tip of the iceberg. The involved Jacobi solver for the pressure Poisson problem is the biggest weak-point of the implementation. There are many other linear solves that are mathematically more efficient (by a large amount).
You find many more such things in the implementation :D.
@@MachineLearningSimulation So it is the Jacobi solver? Good to know, because I have been trying to figure out how this algorithm would converge, but now I know the name.
I have looked into it and cannot figure out how you solved the pressure Poisson equation. From what I understand, you are solving the laplace(p_prev) - rho/delta t * gradient(u) = 0. in every interior position in p_next, you write laplace(p_prev) - rho/delta t * gradient(u). But that is formula, not an equation. At what point is laplace(p_prev) - rho/delta t * gradient(u) = 0 actually solved? What is the x value you are looking for?
I guess the real question is:
This video gives the equation you solve when using the Jacobi iteration method
ruclips.net/video/VH0TZlkZPRo/видео.html&ab_channel=EMPossible
I can't figure out how this codes translate into solving this equation. Which part is the vector b, which is the matrix A and what is x in the code. Could you please elaborate?
Hello!
Did you used staggered grid or collocated grid for this simulation?
Hi,
that's a colocated grid here.
@@MachineLearningSimulation Hi!
To my knowledge, we need an additional padding of grid on every sides to implement the boundary conditions. But the code doesn't have any. Please let me know about this.
I think you are referring to a cell-centered discretization which is common when using finite volumes. Here, however, we have a vertex-centered discretization, mostly (but necessarily) due to usage of Finite differences. Hence, we have degrees of freedom that lie directly on the boundary which allows us to prescribe both Dirichlet BCs for the Velocity and Neumann BCs for the pressure. The latter are unfortunately only of first order consistency.
Is that what you were referring to? ☺️
@@MachineLearningSimulation Yes, this is what I am referring to. I thought that for pressure projection, we need to discretize using approach similar to finite volume, with pressure being cell centered (same as your pipe flow video), otherwise it might create pressure oscillations.
@@MachineLearningSimulation If we have a cell based discretization, will the results alter?
Sir, can you please explain why there is a kind of phase separation in the streamline plot???? Will it persist in this particular problem of lid driven cavity throughout the run time ???
Hi,
thanks for the comment :).
Can you elaborate what you mean by phase separation?
@@MachineLearningSimulation Sir , actually I just want to know the explanation of what that different colour layers does mean which are appearing at time of simulation ??? Physical explanation
The contour plot is the pseudo-pressure. The streamlines are created by the velocity field.
Is that what you refer to?
@@MachineLearningSimulation Can you please explain the changing of pressure physically
The central difference is wrong, that's why it's not stable even when the CFL condition is met. It's -2 * central point, not -4.
Do you refer to the central difference approximation of the Laplace operator or of the first derivative? :) In case of the Laplace operator, "-4" in front of the central point should be correct since we are taking a second order approximation of the second derivative in both x and y-axis. Hence, the stencil is
0 1 0
1 -4 1
0 1 0
For the central difference approximation of the first derivative (which we used to evaluate the convection/advection), the second order stencil does not include the point itself. Hence, it is
-1 0 1 (wrt x)
or
1
0
-1 (wrt y)
What do you think?
Btw: the stability check is not for the CFL condition, but for the stability of the FTCS scheme of the heat equation (en.wikipedia.org/wiki/FTCS_scheme#Stability ). This is because we can only simulate diffusion dominant flows with explicit central difference approximations to the advection. For higher Reynolds numbers, one would observe pressure oscillations. For these high diffusion cases (as simulated here), the FTCS stability criterion is often stricter.
Hope I got you right :). There is of course a chance that I made a mistake. Please let me know, then I can fix the file on GitHub and leave a remark under the video.
Isn't it PISO algorithm?
I think that there are definitely resemblances to what the original PISO was doing, but I would not coin the approach here PISO. More generally speaking, the PISO algorithm (and also the SIMPLE algorithm) are just some special treatments of the coupled NS problem, see e.g. "Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach" by Stefan Turek.
The approach done here is inspired by Lorena Barba's "12 steps to Navier Stokes course", but uses the simpler incompressibility projection by Helmholtz decomposition.
@@MachineLearningSimulation thank you for the helpful explanation
I am using VS Code to run the program, I have varied the TIME_STEP_LENGTH, but I keep getting the following errors:
" 0%| | 0/500 [00:00
Hi,
it seems to me that your simulation turned unstable. As such, some (and eventually all) your floating point values in the arrays turned NaN (not a number) which what you're error messages are telling you.
This is caused by a too high time step. If you lower it again, it should be fine.
The simulation here is only conditionally stable since we are using both an explicit convection and an explicit diffusion treatment. Hence, there will be some condition for the time step length in relation to the spatial discretization length (and some other quantities) that has to be fulfilled.