#190 Performance changes to linear_wake analysis #348
Conversation
There was a problem hiding this comment.
Thanks for this cleaning PR @MathisMewes! The code looks nicer without these loops indeed. You will find comments & questions below.
| for i in range( int(nz/2) -1): | ||
| nb_array[int(nz/2)-i ] = peak_density * np.exp(-0.5*((i*dzeta)/sigma_z)**2 ) | ||
| nb_array[int(nz/2)+i ] = peak_density * np.exp(-0.5*((i*dzeta)/sigma_z)**2 ) | ||
| nb_array = peak_density*np.sqrt(2*np.pi)*norm.pdf(np.linspace(-nz/2,nz/2,nz)*dzeta/sigma_z) |
There was a problem hiding this comment.
Did you print both arrays and check that they are equal? I tried this with random values of nz etc. and they seem to differ.
There was a problem hiding this comment.
There are 2 minor differences, but they should be functionally equivalent:
The old code shifts the center of the gauss function away from 0, if nz is even.
The new code keeps the gaussfunction centered around 0.
This shift could be added to the new method aswell, but it seems to be a redundant feature of the discretization.
I also noticed a difference for very small nz (e.g. 5). The numerical methods seem to diverge a bit there.
There was a problem hiding this comment.
I actually just found an even more efficient way, using the same function as the old one.
But the minor differences are still present.
There was a problem hiding this comment.
I actually found a small mistake in my approach and fixed it. Now there are just two optional differences to the old one remaining, which i could also adopt in the new code if you prefer.
The Old code makes sure to keep the borders at 0. In the symmetric case, shown here, it even leaves two steps
As seen here the old code shifts the center of the distribution to the nearest discretized step, while the new one keeps it at the real center, which is between two steps.
There was a problem hiding this comment.
OK, if the differences are small then we're good.
| for i in range(nz-1): | ||
| nb_dzdz[i] = (nb_array[i-1] -2*nb_array[i] + nb_array[i+1] )/dzeta**2 | ||
| nb_dzdz[1:nz-1] = (nb_array[0:nz-2] - 2*nb_array[1:nz-1] + nb_array[2:nz])/dzeta**2 |
There was a problem hiding this comment.
Same question here (although I haven't tested)
There was a problem hiding this comment.
This takes the second derivative, so the first and last element become useless.
The old method actually mistakenly calculated it for the first element aswell, while leaving the last at 0.
It used the value from the last element in the process (nb_array[-1]), which can lead to some unexpected results.
The new Method simply keeps both first and last element at 0. Could be included again without trouble, if you want.
The rest is exactly the same.
| tmp = np.zeros([nz,nz],dtype=float) | ||
| for i in np.arange(nz-1,-1,-1): | ||
| tmp = 0. | ||
| for j in range(nz-i): | ||
| tmp += 1./kp*math.sin(kp*dzeta*(i-(nz-1-j)))*nb_dzdz[nz-1-j] | ||
| n_th[i] = tmp*dzeta + nb_array[i] | ||
| tmp[i,j]= i-(nz-1-j) | ||
| tmp = (dzeta/kp*np.sin(kp*dzeta*tmp) * np.full([nz,nz],1) * nb_dzdz[np.linspace(nz-1,0,nz,dtype=int)]) | ||
| n_th = np.sum(tmp,axis = 1) + nb_array |
There was a problem hiding this comment.
Same question here (although I haven't tested)
There was a problem hiding this comment.
This one is exactly equal aside from numerical rounding error
| for i in range( int(nz/2) -1): | ||
| nb_array[int(nz/2)-i ] = peak_density * np.exp(-0.5*((i*dzeta)/sigma_z)**2 ) | ||
| nb_array[int(nz/2)+i ] = peak_density * np.exp(-0.5*((i*dzeta)/sigma_z)**2 ) | ||
| nb_array = peak_density*np.sqrt(2*np.pi)*norm.pdf(np.linspace(-nz/2,nz/2,nz)*dzeta/sigma_z) |
There was a problem hiding this comment.
OK, if the differences are small then we're good.
Working on #190 :
analysis of linear_wake example:
Implemented more effective code using numpy commands instead of for loops.
Compared results between new and old scripts applied to the same simulation data. Outcome seems to be consistent.
Small enough (< few 100s of lines), otherwise it should probably be split into smaller PRs
Tested (describe the tests in the PR description)
Runs on GPU (basic: the code compiles and run well with the new module)
Contains an automated test (checksum and/or comparison with theory)
Documented: all elements (classes and their members, functions, namespaces, etc.) are documented
Constified (All that can be
constisconst)Code is clean (no unwanted comments, )
Style and code conventions are respected at the bottom of https://github.com/Hi-PACE/hipace
Proper label and GitHub project, if applicable