barbagroup / cfdpython Goto Github PK
View Code? Open in Web Editor NEWA sequence of Jupyter notebooks featuring the "12 Steps to Navier-Stokes" http://lorenabarba.com/
License: Other
A sequence of Jupyter notebooks featuring the "12 Steps to Navier-Stokes" http://lorenabarba.com/
License: Other
In code section "In [4]" BC determination:
p[:, -1] = p[:, -2] ##dp/dy = 0 at x = 2
In fact it is ##dp/dx = 0 at x = 2
I believe the initialization in the Step 7 diffuse() function is incomplete: Only the nodes with an initial value of 2 are reset; all other nodes are left as they were after the last call to diffuse(). For example, if you run diffuse(10) twice in succession and look at statistics such as the maximum and minimum value of u, you will see that the results differ in the two runs. This occurs because the initial values outside of [0.5, 1.0], [0.5, 1.0] have not been reset to 1.0.
This can be remedied by adding "u = numpy.ones((ny, nx)) as the first statement in diffuse().
Thanks for putting these steps to CFD together; they have been very helpful to me. I hope this is a little helpful for you.
Bob
In line #435-#437 in _05_Step_4.ipynb_, when i=0
, it will use u[-1]
as its left-hand-side point. However, with periodic boundary condition, because u[0] = u[-1]
, hence u[-2]
instead of u[-1]
is the correct left-hand-side point of u[0]
.
can I use it for micropolar fluid?
Shouldn't the line:
utoE = lambda u: (u/2)**2
be
utoE = lambda u: (u**2)/2
?
Looking at one of the early notebooks, you have a call to action: What do you observe about the evolution of the hat function under the non-linear convection equation? What happens when you change the numerical parameters and run again?
I wondered if you had considered making an interactive out of this using ipywidgets
? (Or maybe you did and discounted it because it can detract from purposeful exploration and engagement with the code in favour of letting students just play with the interactive shiny bits?!)
For openjournals/jose-reviews#21.
If you plot the pressure profile after convergence of the Navier-Stokes equations coupled with the pressure-poisson equation there is no pressure gradient, is this an artifact of the periodic boundary conditions?
The pressure field does not change at all during iteration...
Now that display()
is a built-in, we don't need calls like from IPython.display import display
See: http://blog.jupyter.org/2017/05/31/release-of-ipython-5-4-6-1-and-rlipython-2/
There have been a few issues around the periodic boundary conditions in Step 4.
At the very least, they need further clarification, especially if we treat u[0] == u[-1]
(which we do, for the moment)
I'd like to move the lesson on defining functions up the stack so student see it sooner.
So we'll introduce 1d linear and non-linear convection, then introduce functions and show off a JSAnimation example. Then when we get to 1d diffusion and burgers, we can embed animations to better show the behavior of the discretizations.
I noticed what I think is a mistake in the derivation for central difference scheme for the 2nd order derivative in Step 3:
It is stated that the 4th order terms are negligible and can be discarded, however this term remains in the next equation. I believe these terms remain until we divide all terms in the equation by dx ^ 2 , thus leaving the remaining 4th order terms to 2nd order terms, yielding the scheme 2nd order accurate.
Hope this makes sense :)
When plotting the numerical solution at step 4, initial conditions of the equation was taken from analytical solution (by using un = u.copy() at the beginning of the FD loop, array u is the analytical solution of function at t=0 in our code). I think we need to use initial conditions below in the picture.
In the fourth code cell in 12_Step_9.ipynb, dy is set to 2 / (ny - 1). To be consistent with the y array, this should be 1 / (ny - 1).
Correcting this has several effects. First, convergence slows. With the current version, laplace2d() stops iterating at 1,133 steps; after applying the correction, it takes 2,043 steps. Second, the revision affects how close laplace2d() comes to the analytic solution. I use sum(p) as a rough comparison among different results from laplace2d(). The current version (2/(ny - 1)) produces a sum of 231.8 when iteration stops at 1,133 steps. The corrected initialization (1/(ny - 1)) produces 220.2 at 2,043 steps. Compare these with the analytic solution sum of 240.25.
I'm translating to Julia as I work through the Steps. My Julia versions of Step 9 come up with the same results as the revised Python version. In addition, it shows that if you iterate the Julia version of laplace2d() for 5,000 steps, you come up pretty close the analytic solution (sum of 239.5 vs. the 240.25 noted above).
Hope this helps. Thanks again for all your work on this. In addition, thank you for keeping it current.
On to Step 10!
Thank you for this great course! I have been using this course to teach CFD-concepts in Belgium.
I'm stuck with the following question. In video lecture 11, a poisson equation is obtained to correct for the divergence in the velocity.
However in notebook 14 the discretised pressure-Poisson equation looks completely different than the one obtained
in the video lecture: https://youtu.be/ZjfxA3qq2Lg?t=1647
My question is: what happened in between the curved brackets in the video lecture and equation 13 in the notebook.
In the 11th step, the solution of the pressure equation is -ve , which is physically invalid. And in the Pressure Poisson equation's discretization, there is no time variable so, how the term containing dt arrives
gca()
with keyword arguments was deprecated in Matplotlib 3.4
usingsubplot()
or add_subplot()
may resolve the question
CC-BY only applies to text and image content. For code, we need to use an open-source code license (MIT license is our preference).
Thank you for this wonderful stuff, professor.
Are the discretization grids used in the cavity flow simulation staggered or collocated ones? From the indexing of the prime variables u, v, p it looks like the collocated grid. But the use of fractional step integrating and pressure correction means staggered gridding...? As a newbie, I am confused.
For openjournals/jose-reviews#21.
The current instructions simply tell users to run the notebook server, but they should probably also include cloning the repository, then moving into the directory, running the server, then opening the first notebook.
The Python module versions in requirements.txt are a bit outdated.
I installed all the dependencies using my system's package manager, and while those packages aren't all at the latest version, some of them are quite far ahead of those in requirements.txt.
These are the versions I have tested:
numpy: 1.21.5
matplotlib: 3.5.2
scipy: 1.8.1
sympy: 1.10.1
There weren't any issues, except for a deprecation in matplotlib 3.6 (changelog is here):
MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
ax = fig.gca(projection='3d')
The fix for this actually quite easy: Replace
ax = fig.gca(projection='3d')
with
ax = fig.add_subplot(projection='3d')
Also, when I checked the matplotlib 2.2.x documentation, it didn't even mention the projection
argument on the gca()
function.
Please update the dependency list to more recent versions.
You could also use compatible versions instead of exact versions. For example:
# requirements.txt
ipywidgets ~=8.0
jupyter ~=1.0.0
numpy ~=1.23
matplotlib ~=3.6.0
scipy ~=1.9
sympy ~=1.11
According to the link below, the L1 norm of a matrix is calculated by summing up absolute values in each column, and then taking the maximum
https://en.wikipedia.org/wiki/Matrix_norm
The notebook for Step 9 calculates the norm by summing up the absolute values of all the elements in the matrix. Also the link provided to read up further on the L1 norm leads to the wiki page for L1 norm of a vector (which I understand is different from the L1 norm of a matrix)
Running 01_Step_1.ipynb I get a division by zero error. The error goes away by setting
dx = 2.0/(nx-1) instead of dx = 2/(nx-1). I am using anaconda Python 2.7.11 :: Anaconda 2.5.0 (x86_64).
Thanks!
It's time.
This CFDPython lessons is the best leasons I 've ever token about CFD. I follow every step until step11. But I got stuck in equation 4 of step11 :
∂2p/∂x2+∂2p/∂y2=−ρ(∂u/∂x∂u/∂x+2∂u/∂y∂v/∂x+∂v/∂y*∂v/∂y)
where did the Poisson equation come from. And How to derive this equation.
I have been confused for month about this equation.
Could be wring, but I don't think the 'dt' term belongs there. The poisson pressure equation doesn't have a 'dt' term.
In linear convection's equation:
`u = numpy.ones(nx)
u[int(.5/dx):int(1 / dx + 1)] = 2`
Line 2 will give a zero error.
Hello, please someone answere me, wherefrom the 1/delta(t)*div(U) term is appearing in Poisson eq? In previous steps was used that div(U) = 0 .... just don't get it
cell[3] :
u[.5/dx : 1/dx+1]=2 #setting u = 2 between 0.5 and 1 as per our I.C.s
Floating point used as indices implies a cast to integer which is probably too much for beginners.
It's minor, but the formula for getting dt is already assuming c=1, which could cause someone to trip over this …
def linearconv(nx):
dx = 2./(nx-1)
nt = 20 #nt is the number of timesteps we want to calculate
c = 1
sigma = .1
dt = sigma*dx
(submitted by email by Richard Galvez, Physics PhD Candidate, Syracuse University)
There is a strange (1 / Delta t) term in the discretized Poisson equation that is present in the written equation, the Python code (in the build_up_b function), and the corresponding Youtube video. This seems to make the solution unstable for small time steps. I think this may be due to a misapplication of discretization to a time derivative of div v. For an incompressible fluid, div v = 0 so this term should vanish. In any case, just tacking on a (1 / Delta t) is not the correct way to discretize a time derivative.
Also need to add semi-colons to suppress output where necessary
Dear BarbaGroud,
Thank you for sharing your work, it's awesome. But I'm having a hard time understanding your 2D implementation. In every of the notebooks I can see that the u is initialised as:
u = np.ones((ny,nx))
This means that the rows of the u array are related to the y direction and the columns are related to the x direction.
Furthermore when (for example) the 2D linear burger equation is solved you are using the following expression:
c_dt/dx_(un[i,j] - un[i-1,j])
To evaluate the derivative in the x-direction. As far as I'm concerned the rows of un, do have the information of the y axis, and thus you should be using dy instead of dx. This problem is not raised up because you are using the same dimensions and subdivisions for the x and y axis, but whenever this is changed there will be some errors.
The problem could be solved changing the u initialisation to:
u = np.ones((nx,ny))
Although the plots should be using the transposed of u in order to get the right axis.
Hope I'm wrong with my comments, but I've been thinking about this for a long time and I just cannot see if I'm making a mistake.
Thank you very much.
And change the appropriate line in the intro notebooks about calling ipython notebook --inline
Hey, when handling the periodic boundary conditions in Step 4, after the second for loop of code snippet #10, you have
u[-1] = un[-1] - un[-1] * dt / dx * (un[-1] - un[-2]) + nu * dt / dx**2 * (un[0] - 2 * un[-1] + un[-2])
where in the last set of parentheses you obtained un[i + 1]
= un[-1 + 1]
= un[0]
.
Given the boundary u(0) = u(2*pi), shouldn't it be un[i + 1]
= un[0 + 1]
= un[1]
instead of un[0]
?
Thanks in advance!
In . There is a little bug in sentence:
u[1:, 1:] = un[1:, 1:] - ((c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -
(c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
The bold left bracket "(" is in a wrong place. It should follow "=" . Correct code should be:
u[1:, 1:] = ( un[1:, 1:] -(c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -
(c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
This little bug make cell5 and cell6 show a different result which suppose to be same.
Splitting hairs, yes, but I believe that's the correct way to combine names of multiple people.
Especially re: pylab inline flags being deprecated
The statement of need in paper.md
does a good job describing the effectiveness of the lessons, but doesn't clearly explain who would benefit from these lessons and why. For example, is this geared towards graduate students in mechanical engineering, physics, math, etc.?
A statement of need: Do the authors clearly state the need for this module and who the target audience is?
The Jupyter links of 12 steps are invalid.
Need to clean up all of these notebooks and get rid of spurious random text. Should also add a quick explanation about the text suppression the first time we use it
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.