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!

can I use it for micropolar fluid?

Also need to add semi-colons to suppress output where necessary

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.

It's time.

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

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!

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

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.

Splitting hairs, yes, but I believe that's the correct way to combine names of multiple people.

Looking at 10_Step_8.
In the third set of equations when solving for the next time step and

The last terms for both equations need to change one of the indices to . Like this:

Thank you!

And change the appropriate line in the intro notebooks about calling ipython notebook --inline

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 :)

For openjournals/jose-reviews#21.

Could be wring, but I don't think the 'dt' term belongs there. The poisson pressure equation doesn't have a 'dt' term.

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.

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.

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)

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).

For openjournals/jose-reviews#21.

• Third paragraph in section "Statement of Need": usting >> using

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)

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

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)

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

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?!)

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.

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.

Shouldn't the line:

utoE = lambda u: (u/2)**2

be

utoE = lambda u: (u**2)/2

?

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?

openjournals/jose-reviews#21

Especially re: pylab inline flags being deprecated

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/

The Jupyter links of 12 steps are invalid.

For openjournals/jose-reviews#21.

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.

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

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...

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.

Correctly displaying at Step12 -> copy from there where applicable.

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.

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!

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] .

Hello! I've been trying to watch the videos in the notebooks, but I could not reproduce them, see image below:

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

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.

gca() with keyword arguments was deprecated in Matplotlib 3.4
usingsubplot() or add_subplot()may resolve the question