The course is offered under CC-BY, which is fantastic. This is not obvious from the repository (unless you dig into the iPython notebooks). A top-level license file makes this obvious both to people browsing and to search engines.

Fix #40

Use Latex split environment. Remove line breaks to display equation correctly.

You can view the modified notebooks with Nbviewer to check the equations:

• Step 7: before and after;
• Step 8: before and after;
• Step 10: before and after;
• Step 11: before and after;
• Step 12: before and after.

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!

The u[0] calculation had terms that were un[-2] in both spatial derivatives. Those should be un[-1].

Additionally it is worth noting that the periodic conditions can all be handled correctly by the u[i] calculation as written by simply changing the range to for i in range(-1, nx-1) as below. This makes use of python negative indexing which wraps to the end of the array, effectively computing u[-1] and u[0] first.

for n in range(nt):
    un = u.copy()
    for i in range(-1, nx-1):
        u[i] = un[i] - un[i] * dt/dx *(un[i] - un[i-1]) + nu*dt/dx**2*\
                (un[i+1]-2*un[i]+un[i-1])

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.

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!

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

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.

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)

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

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.

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

Previously, there has been a typo throughout the code implementation in Step 12 in the periodic boundary conditions. As a quick example, If u_{i,j} was at the right boundary, then u_{i+1,j} was took to be at the left boundary i.e u_{0,j} (instead of u_{1,j}). This PR fixes those issues, although hardly changes the results.

Bumps numpy from 1.14.3 to 1.22.0.

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting @dependabot rebase.

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.