Skip to content

Commit

Permalink
soft update for the ion ion beam functional test (#726)
Browse files Browse the repository at this point in the history 
 few comments and new var to make the test more explicit

---------

Co-authored-by: Loïc Darrieumerlou <[email protected]>
  • Loading branch information
@rochSmets
rochSmets and Loïc Darrieumerlou authored Feb 13, 2023
1 parent 774fd61 commit b1a9ec3
Showing 1 changed file with 86 additions and 39 deletions.
125 changes: 86 additions & 39 deletions tests/functional/ionIonBeam/ion_ion_beam1d.py
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,8 @@
#!/usr/bin/env python3

import pyphare.pharein as ph #lgtm [py/import-and-import-from]
from pyphare.pharein import Simulation
from pyphare.pharein import MaxwellianFluidModel
from pyphare.pharein import ElectromagDiagnostics, FluidDiagnostics
from pyphare.pharein import ElectromagDiagnostics, FluidDiagnostics, ParticleDiagnostics
from pyphare.pharein import ElectronModel
from pyphare.simulator.simulator import Simulator
from pyphare.pharein import global_vars as gv
Expand All @@ -18,6 +17,10 @@
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from pyphare.pharesee import plotting
from pathlib import Path
mpl.use('Agg')


Expand All @@ -30,8 +33,8 @@ def config():
Simulation(
smallest_patch_size=20,
largest_patch_size=60,
final_time=50,
time_step=0.0005,
final_time=100,
time_step=0.001,
boundary_types="periodic",
cells=165,
dl=0.2,
Expand Down Expand Up @@ -69,7 +72,7 @@ def vth(x):

vMain = {
"vbulkx": v0, "vbulky": v0, "vbulkz": v0,
"vthx": vth, "vthy": vth, "vthz": vth
"vthx": vth, "vthy": vth, "vthz": vth,
}


Expand All @@ -89,14 +92,16 @@ def vth(x):

sim = ph.global_vars.sim

timestamps = np.arange(0, sim.final_time, 1.)
timestamps = np.arange(0, sim.final_time, 0.1)

for quantity in ["B"]:
for quantity in ["B","E"]:
ElectromagDiagnostics(
quantity=quantity,
write_timestamps=timestamps,
compute_timestamps=timestamps,
)
ParticleDiagnostics(quantity = "domain", population_name = "main" , write_timestamps=timestamps, compute_timestamps=timestamps)
ParticleDiagnostics(quantity = "domain", population_name = "beam" , write_timestamps=timestamps, compute_timestamps=timestamps)



Expand All @@ -106,14 +111,14 @@ def yaebx(x, a, b):
return a*np.exp(np.multiply(b, x))


def growth_b_right_hand(run_path):
def growth_b_right_hand(run_path, time_offset):

file = os.path.join(run_path, "EM_B.h5")
times = get_times_from_h5(file)

dt = times[1]-times[0]
r = Run(run_path)
first_mode = np.array([])

from scipy.optimize import curve_fit

for time in times:
B_hier = r.GetB(time, merged=True, interp="linear")
Expand All @@ -125,16 +130,21 @@ def growth_b_right_hand(run_path):
x = xyz_finest[0][:-1]

by = by_interpolator(x)
bz = by_interpolator(x)
bz = bz_interpolator(x)

# get the mode 1, as it is the most unstable in a box of length 33
mode1 = np.absolute(np.fft.fft(by-1j*bz)[1])
first_mode = np.append(first_mode, mode1)

popt, pcov = curve_fit(yaebx, times, first_mode, p0=[0.08, 0.09])
ioffset = int(time_offset/dt)
imax = find_peaks(first_mode,width = ioffset)[0][0]

# the curve_fit is performed from time index 0 to imax-ioffset as this offset prevent to use
# the final part of the curve which is no more exponential as this is the end of the linear mode
popt, pcov = curve_fit(yaebx, times[:imax-ioffset], first_mode[:imax-ioffset], p0=[0.08, 0.09])

# now the signal is stripped from its exponential part
damped_mode=first_mode*yaebx(times, 1/popt[0], -popt[1])
damped_mode=first_mode[:imax-ioffset]*yaebx(times[:imax-ioffset], 1/popt[0], -popt[1])

# find the omega for which "damped_mode" is the largest :
# this term is twice the one it should be because "mode1" resulting from
Expand All @@ -143,52 +153,89 @@ def growth_b_right_hand(run_path):
# the factor "+1" is because we remove the DC component, so the value
# given by argmax has also to miss this value
omegas = np.fabs(np.fft.fft(damped_mode).real)
omega = 0.5*(omegas[1:omegas.size//2].argmax()+1)*2*np.pi/times[-1]
omega = 0.5*(omegas[1:omegas.size//2].argmax()+1)*2*np.pi/times[imax-1-ioffset]

print(omegas[0:6])
print(omegas[1:omegas.size//2].argmax())
print(2*np.pi/times[-1])
print(0.5*(omegas[1:omegas.size//2].argmax()+1)*2*np.pi/times[-1])

return times, first_mode, popt[0], popt[1], damped_mode, omega


def main():

from pybindlibs.cpp import mpi_rank
time_offset = 10.
# this is an offset so the exponential fit associated to the linear phase is not performed
# until the time at which the B value gets the greater... but is rather before with an offset
# given by time_offset : the exponential fit is performed on [0, times[imax] - time_offset]

config()
simulator = Simulator(gv.sim)
simulator.initialize()

simulator.run()

if mpi_rank() == 0:

times, first_mode, ampl, gamma, damped_mode, omega \
= growth_b_right_hand(os.path.join(os.curdir, "ion_ion_beam1d"))

fig, (ax1, ax2) = plt.subplots(2, 1)

ax1.set_title("Right Hand Resonant mode (Beam instability)")
ax1.stem(times, first_mode, linefmt='-k', basefmt=' ', use_line_collection=True)
ax1.plot(times, yaebx(times, ampl, gamma), color='r', linestyle='-', marker='')
ax1.text(0.04, 0.80, "From Gary et al., 1985 (ApJ : 10.1086/162797)", transform=ax1.transAxes)
ax1.set_ylabel("Most unstable mode")
ax1.set_title("Right Hand Resonant mode (Beam instability)")
ax1.text(0.30, 0.50, "gamma = {:5.3f}... expected 0.09".format(gamma), transform=ax1.transAxes)

ax2.plot(times, damped_mode, color='g', linestyle='', marker='o')
ax2.set_xlabel("Time")
ax2.set_ylabel("Real mode")
ax2.text(0.48, 0.30, "~ 3 periods until t=50", transform=ax2.transAxes)
ax2.text(0.40, 0.20, "omega (real) = {:5.3f}... expected 0.19".format(omega), transform=ax2.transAxes)

times, first_mode, ampl, gamma, damped_mode, omega \
= growth_b_right_hand(os.path.join(os.curdir, "ion_ion_beam1d"), time_offset)

dt = times[1]-times[0]
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(6.5,12), constrained_layout=True)
imax = find_peaks(first_mode, width = int(10/dt))[0][0] # the smallest width of the peak is 10

ax1.set_title("Time evolution of the first right-hand circular mode amplitude")
ax1.plot(times, first_mode, color = 'k',label ="|$\hat{B}_y (m=1,t)-i\hat{B}_z (m=1,t)$|")
ax1.plot(times[:imax], yaebx(times[:imax], ampl, gamma), color='r', linestyle='-', label ="$B_0. \exp(\gamma t), \ with\ \gamma =${:5.5f} (expected 0.09)".format(gamma))
ax1.axvline(0, 0, yaebx(times[imax],ampl,gamma), color = 'red',linestyle = '--')
ax1.axvline(times[imax]-time_offset, 0, yaebx(times[imax],ampl,gamma), color = 'red',linestyle = '--')
ax1.legend()
ax1.set_xlabel("t - Time")

r = Run(os.path.join(os.curdir, "ion_ion_beam1d"))
dt = times[1]-times[0]
ax_t2, ax_t3, ax_t4 = ax2.twinx(), ax3.twinx(), ax4.twinx()
vmin, vmax = -1.3, 5.6
for t in [0,times[imax],times[-1]] :
if t == 0 :
ax, ax_t = ax2, ax_t2
elif t == times[-1]:
ax,ax_t = ax4, ax_t4
else :
ax,ax_t = ax3, ax_t3

ions = r.GetParticles(t, ["main","beam"])


ions.dist_plot(axis=("x", "Vx"),
ax = ax,
norm = 0.4,
finest=True,
gaussian_filter_sigma = (1,1) ,
vmin=vmin,vmax=vmax,
dv=0.05,
title="t = {:.1f}".format(t),
xlabel = "",
ylabel = ""
)
ax.set_xlim((0,33))
ax.axvline(10, vmin, vmax, color = 'k')
ax.axvline(22, vmin, vmax, color = 'k')
ax.axvline(14, vmin, vmax, color = 'red')
ax.axvline(18, vmin, vmax, color = 'red')


E_hier = r.GetE(time = t, merged = True, interp = "linear")
ey_interpolator, xyz_finest = E_hier["Ey"]
ax_t.plot(xyz_finest[0], ey_interpolator(xyz_finest[0]), linewidth = 2, color = 'dimgray')
ax_t.set_ylim((-0.4,0.4))

ax_t3.set_ylabel("Ey(x)")
ax3.set_ylabel("Vx - Velocity")
ax4.set_xlabel("X - Position")
fig.savefig("ion_ion_beam1d.png")

# compare with the values given gary et al. 1985
assert np.fabs(gamma-0.09) < 2e-2
#assert np.fabs(omega-0.19) < 8e-2


if __name__=="__main__":
main()

0 comments on commit b1a9ec3

Please sign in to comment.