barbagroup / CFDPython

A sequence of Jupyter notebooks featuring the "12 Steps to Navier-Stokes" http://lorenabarba.com/
Other
3.43k stars 1.19k forks source link

Maybe there is a bug in code of Array Operations with NumPy #60

Closed fidle closed 5 years ago

fidle commented 5 years ago

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.

mesnardo commented 5 years ago

Thank you for opening this issue @fidle. Correct me if I am wrong, but I think you are referring to In [5] and In [6] in Lesson 06 (06_Array_Operations_with_NumPy.ipynb). If so, the parenthesis in In [6] are indeed misplaced.

Here is a small example that reflects this issue:

> python
Python 3.6.8 |Anaconda custom (64-bit)| (default, Dec 29 2018, 19:04:46)
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
>>> numpy.__version__
'1.15.4'
>>> c = 1.0
>>> dx = 2 / 80
>>> dy = 2 / 80
>>> sigma = 0.2
>>> dt = sigma * dx
>>> un = numpy.random.rand(80, 80)
>>> u1 = un[1:, 1:] - ((c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) - (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
>>> u2 = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) - (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
>>> numpy.allclose(u1, u2)
False
>>>

Since you found the bug @fidle, would you like to provide a patch in a pull-request?

fidle commented 5 years ago

sorry to reply you so late. @mesnardo you are right , the parenthesis in In[6] is misplace. I am sorry I did make it clear in the issue. I hope I could make a pull-request. but I don't know how(I am a fresh man in github). I rewrite this part to make the code output a result, so people can see that both method can get the same result. here is the code `import numpy from matplotlib import pyplot

nx = 81 ny = 81 dx = 2/(nx-1) dy = 2/(ny-1) c = 1 nt = 100 sigma =0.2 dt = sigma*dx

x = numpy.linspace(0,2,nx) y = numpy.linspace(0,2,ny) %%timeit u = numpy.ones((ny, nx)) u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps un = u.copy() row, col = u.shape for j in range(1, row): for i in range(1, col): u[j, i] = (un[j, i] - (c dt / dx (un[j, i] - un[j, i - 1])) - (c dt / dy (un[j, i] - un[j - 1, i]))) u[0, :] = 1 u[-1, :] = 1 u[:, 0] = 1 u[:, -1] = 1 pyplot.pcolormesh(u)`

`%%timeit u = numpy.ones((ny, nx)) u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps un = u.copy() 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:])) u[0, :] = 1 u[-1, :] = 1 u[:, 0] = 1 u[:, -1] = 1 pyplot.pcolormesh(u)`

mesnardo commented 5 years ago

Thank you @fidle. I opened a pull-request to correct the parenthesis. Here are the steps I did to submit a pull-request:

(@fidle, I gave you credit in the commit message 3c692c4.)