Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

#190 Performance changes to linear_wake analysis #348

Merged
merged 4 commits into from
Feb 18, 2021
Merged

#190 Performance changes to linear_wake analysis #348

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
99bcf4d
Some performance changes to linear_wake analysis
Feb 12, 2021
7e48717
fixed minor issues
Feb 12, 2021
10a6106
another bugfix
Feb 12, 2021
851ac0c
Update analysis.py
MathisMewes Feb 12, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 7 additions & 8 deletions examples/linear_wake/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

import matplotlib.pyplot as plt
import scipy.constants as scc
from scipy.stats import norm
import matplotlib
import sys
import numpy as np
Expand Down Expand Up @@ -66,24 +67,22 @@
if (args.gaussian_beam):
sigma_z = 1.41 / kp
peak_density = 0.01*ne
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)
Comment on lines -69 to +70
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I actually just found an even more efficient way, using the same function as the old one.
But the minor differences are still present.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

image
The Old code makes sure to keep the borders at 0. In the symmetric case, shown here, it even leaves two steps

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, if the differences are small then we're good.

else:
nb_array[nz-index_beam_head-beam_length_i:nz-index_beam_head] = 0.01 * ne

# calculating the second derivative of the beam density array
nb_dzdz = np.zeros(nz)
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
Comment on lines -77 to +76
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same question here (although I haven't tested)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


# calculating the theoretical plasma density (see Timon Mehrling's thesis page 41)
n_th = np.zeros(nz)
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
Comment on lines +80 to +85
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same question here (although I haven't tested)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one is exactly equal aside from numerical rounding error

rho_th = n_th * q_e

if args.do_plot:
Expand Down