dedebenui/scgenerator
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Some more cleanups, especially tests

Browse Source
This commit is contained in:
Benoît Sierro
2023-07-24 14:28:41 +02:00
parent 41072ea1f2
commit be0a9b8c20
71 changed files with 142 additions and 1529 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
examples/Optica_PM2000D/Optica_PM2000D.toml Executable file
View File

@@ -0,0 +1,18 @@
name = "PM2000D"
mean_power = 0.23
field_file = "Pos30000New.npz"
repetition_rate = 40e6
wavelength = 1546e-9
dt = 1e-15
t_num = 8192
tolerated_error = 1e-6
quantum_noise = true
raman_type = "agrawal"
z_num = 128
length = 0.3
dispersion_file = "PM2000D_2 extrapolated 4 0.npz"
interpolation_degree = 12
A_eff_file = "PM2000D_A_eff_marcuse.npz"
n2 = 4.5e-20

BIN
examples/Optica_PM2000D/PM2000D_2 extrapolated 4 0.npz Normal file
View File

Binary file not shown.

BIN
examples/Optica_PM2000D/PM2000D_A_eff_marcuse.npz Normal file
View File

Binary file not shown.

BIN
examples/Optica_PM2000D/Pos30000New.npz Normal file
View File

Binary file not shown.

70
examples/Optica_PM2000D/test_Optica_PM2000D.py Normal file
View File

@@ -0,0 +1,70 @@
"""
May 2023
Testing the new solver / operators structure
using parameters from the 2022 Optica paper
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import scgenerator as sc
import scgenerator.math as math
import scgenerator.physics.units as units
import scgenerator.plotting as plot
import scgenerator.solver as sol
def main():
params = sc.Parameters(**sc.open_single_config("./tests/Optica_PM2000D/Optica_PM2000D.toml"))
# print(params.nonlinear_operator)
# print(params.compute("dispersion_op"))
# print(params.linear_operator)
# print(params.spec_0)
# print(params.compute("gamma_op"))
#
# plt.plot(params.w, params.linear_operator(0).imag)
# plt.show()
breakpoint()
res = sol.integrate(
params.spec_0, params.length, params.linear_operator, params.nonlinear_operator
)
new_z = np.linspace(0, params.length, 256)
specs2 = math.abs2(res.spectra)
specs2 = units.to_WL(specs2, params.l)
x = params.l
# x = units.THz.inv(w)
# new_x = np.linspace(100, 2200, 1024)
new_x = np.linspace(800e-9, 2000e-9, 1024)
solution = interp1d(res.z, specs2, axis=0)(new_z)
solution = interp1d(x, solution)(new_x)
solution = units.to_log2D(solution)
plt.imshow(
solution,
origin="lower",
aspect="auto",
extent=plot.get_extent(1e9 * new_x, new_z * 1e2),
vmin=-30,
)
plt.show()
fields = np.fft.irfft(res.spectra)
solution = math.abs2(fields)
solution = interp1d(res.z, solution, axis=0)(new_z)
solution.T[:] /= solution.max(axis=1)
plt.imshow(
solution,
origin="lower",
aspect="auto",
extent=plot.get_extent(params.t * 1e15, new_z * 1e2),
)
plt.show()
if __name__ == "__main__":
main()

15
examples/Travers/Travers.toml Executable file
View File

@@ -0,0 +1,15 @@
name = "Travers"
repetition_rate = 1e3
wavelength = 800e-9
tolerated_error = 1e-6
z_num = 128
t_num = 4096
dt = 50e-18
core_radius = 125e-6
gas_name = "helium"
pressure = 400e2
energy = 0.4e-3
width = 10e-15
length = 3
full_field = true
shape = "sech"

69
examples/Travers/Travers2019.py Normal file
View File

@@ -0,0 +1,69 @@
"""
Testing the the solver / operator mechanism with
parameters from the 2019 Travers paper
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import scgenerator as sc
import scgenerator.math as math
import scgenerator.physics.units as units
import scgenerator.plotting as plot
import scgenerator.solver as sol
def main():
params = sc.Parameters(**sc.open_single_config("./tests/Travers/Travers.toml"))
# print(params.nonlinear_operator)
# print(params.compute("dispersion_op"))
# print(params.linear_operator)
# print(params.spec_0)
# print(params.compute("gamma_op"))
#
# plt.plot(params.w, params.linear_operator(0).imag)
# plt.show()
breakpoint()
res = sol.integrate(
params.spec_0, params.length, params.linear_operator, params.nonlinear_operator
)
new_z = np.linspace(0, params.length, 256)
specs2 = math.abs2(res.spectra)
specs2 = units.to_WL(specs2, params.l)
x = params.l
# x = units.THz.inv(w)
# new_x = np.linspace(100, 2200, 1024)
new_x = np.linspace(100e-9, 1200e-9, 1024)
solution = interp1d(res.z, specs2, axis=0)(new_z)
solution = interp1d(x, solution)(new_x)
solution = units.to_log2D(solution)
plt.imshow(
solution,
origin="lower",
aspect="auto",
extent=plot.get_extent(1e9 * new_x, new_z * 1e2),
vmin=-30,
)
plt.show()
fields = np.fft.irfft(res.spectra)
solution = math.abs2(fields)
solution = interp1d(res.z, solution, axis=0)(new_z)
solution.T[:] /= solution.max(axis=1)
plt.imshow(
solution,
origin="lower",
aspect="auto",
extent=plot.get_extent(params.t * 1e15, new_z * 1e2),
)
plt.show()
if __name__ == "__main__":
main()

0
examples/test_rate.py → examples/ion_rate_plot.py
View File

151
examples/test_integrator.py
View File

@@ -1,151 +0,0 @@
from collections import defaultdict
import matplotlib.pyplot as plt
import numpy as np
import pytest
from scipy.interpolate import interp1d
import scgenerator as sc
import scgenerator.operators as op
import scgenerator.solver as so
def test_rk43_absorbtion_only():
n = 129
end = 1.0
h = 2**-3
w_c = np.linspace(-5, 5, n)
ind = np.argsort(w_c)
spec0 = np.exp(-(w_c**2))
init_state = op.SimulationState(spec0, end, h)
lin = op.envelope_linear_operator(
op.no_op_freq(n),
op.constant_array_operator(np.ones(n) * np.log(2)),
)
non_lin = op.no_op_freq(n)
judge = so.adaptive_judge(1e-6, 4)
stepper = so.ERK43(lin, non_lin)
for state in so.integrate(stepper, init_state, h, judge, max_step_size=0.125):
if state.z >= end:
break
assert np.max(state.spectrum2) == pytest.approx(0.5)
def test_rk43_soliton(plot=False):
n = 1024
b2 = sc.fiber.D_to_beta2(sc.units.D_ps_nm_km(24), 835e-9)
gamma = 0.08
t0_fwhm = 50e-15
p0 = 1.26e3
t0 = sc.pulse.width_to_t0(t0_fwhm, "sech")
print(np.sqrt(t0**2 / np.abs(b2) * gamma * p0))
disp_len = t0**2 / np.abs(b2)
end = disp_len * 0.5 * np.pi
targets = np.linspace(0, end, 32)
print(end)
h = 2**-6
t = np.linspace(-200e-15, 200e-15, n)
w_c = np.pi * 2 * np.fft.fftfreq(n, t[1] - t[0])
ind = np.argsort(w_c)
field0 = sc.pulse.sech_pulse(t, t0, p0)
init_state = op.SimulationState(np.fft.fft(field0), end, h)
no_op = op.no_op_freq(n)
lin = op.envelope_linear_operator(
op.constant_polynomial_dispersion([b2], w_c),
no_op,
# op.constant_array_operator(np.ones(n) * np.log(2)),
)
non_lin = op.envelope_nonlinear_operator(
op.constant_array_operator(np.ones(n) * gamma),
no_op,
op.envelope_spm(0),
no_op,
)
# new_state = init_state.copy()
# plt.plot(t, init_state.field2)
# new_state.set_spectrum(non_lin(init_state))
# plt.plot(t, new_state.field2)
# new_state.set_spectrum(lin(init_state))
# plt.plot(t, new_state.field2)
# print(new_state.spectrum2.max())
# plt.show()
# return
judge = so.adaptive_judge(1e-6, 4)
stepper = so.ERKIP43Stepper(lin, non_lin)
# stepper.set_state(init_state)
# state, error = stepper.take_step(init_state, 1e-3, True)
# print(error)
# plt.plot(t, stepper.fine.field2)
# plt.plot(t, stepper.coarse.field2)
# plt.show()
# return
target = 0
stats = defaultdict(list)
saved = []
zs = []
for state in so.integrate(stepper, init_state, h, judge, max_step_size=0.125):
# print(f"z = {state.z*100:.2f}")
saved.append(state.spectrum2[ind])
zs.append(state.z)
for k, v in state.stats.items():
stats[k].append(v)
if state.z > end:
break
print(len(saved))
if plot:
interp = interp1d(zs, saved, axis=0)
plt.imshow(sc.units.to_log(interp(targets)), origin="lower", aspect="auto", vmin=-40)
plt.show()
plt.plot(stats["z"][1:], np.diff(stats["z"]))
plt.show()
def test_simple_euler():
n = 129
end = 1.0
h = 2**-3
w_c = np.linspace(-5, 5, n)
ind = np.argsort(w_c)
spec0 = np.exp(-(w_c**2))
init_state = op.SimulationState(spec0, end, h)
lin = op.envelope_linear_operator(
op.no_op_freq(n),
op.constant_array_operator(np.ones(n) * np.log(2)),
)
euler = so.ConstantEuler(lin)
target = 0
end = 1.0
h = 2**-6
for state in so.integrate(euler, init_state, h):
if state.z >= target:
target += 0.125
plt.plot(w_c[ind], state.spectrum2[ind], label=f"z={state.z:.3f}")
if target > end:
print(np.max(state.spectrum2))
break
plt.title(f"{h = }")
plt.legend()
plt.show()
def benchmark():
for _ in range(50):
test_rk43_soliton()
if __name__ == "__main__":
test_rk43_soliton()
benchmark()