diff --git a/src/scgenerator/__init__.py b/src/scgenerator/__init__.py index 29b9853..0eeca53 100644 --- a/src/scgenerator/__init__.py +++ b/src/scgenerator/__init__.py @@ -4,7 +4,7 @@ from .evaluator import Evaluator from .legacy import convert_sim_folder from .math import abs2, argclosest, normalized, span from .parameter import Configuration, Parameters -from .physics import fiber, materials, pulse, simulate, units +from .physics import fiber, materials, pulse, simulate, units, plasma from .physics.simulate import RK4IP, parallel_RK4IP, run_simulation from .physics.units import PlotRange from .plotting import ( @@ -18,5 +18,5 @@ from .plotting import ( transform_mean_values, ) from .spectra import SimulationSeries, Spectrum -from .utils import Paths, _open_config, open_single_config +from .utils import Paths, _open_config, open_single_config, simulations_list from .variationer import DescriptorDict, VariationDescriptor, Variationer, VariationSpecsError diff --git a/src/scgenerator/const.py b/src/scgenerator/const.py index 268643f..ea4eee4 100644 --- a/src/scgenerator/const.py +++ b/src/scgenerator/const.py @@ -2,6 +2,12 @@ __version__ = "0.2.4dev" from typing import Any +ONE_2 = 1 / 2 +ONE_3 = 1 / 3 +ONE_FOURTH = 1 / 4 +ONE_FIFTH = 1 / 5 +ONE_6 = 1 / 6 + def pbar_format(worker_id: int) -> dict[str, Any]: if worker_id == 0: diff --git a/src/scgenerator/math.py b/src/scgenerator/math.py index 1213030..756d40e 100644 --- a/src/scgenerator/math.py +++ b/src/scgenerator/math.py @@ -1,7 +1,6 @@ from typing import Union import numpy as np -from scipy.integrate import cumulative_trapezoid from scipy.special import jn_zeros import numba @@ -11,10 +10,9 @@ pi = np.pi c = 299792458.0 -def inverse_integral_exponential(y: np.ndarray, dx: float) -> np.ndarray: - """evaluates exp(-func(y)) from x=-inf to x""" - # return np.exp(-cumulative_trapezoid(y, dx=dx, initial=0)) - return np.exp(-cumulative_simpson(y) * dx) +def expm1_int(y: np.ndarray, dx: float) -> np.ndarray: + """evaluates 1 - exp( -∫func(y(x))dx ) from x=-inf to x""" + return -np.expm1(-cumulative_simpson(y) * dx) def span(*vec): diff --git a/src/scgenerator/operators.py b/src/scgenerator/operators.py index 3889372..eb80db3 100644 --- a/src/scgenerator/operators.py +++ b/src/scgenerator/operators.py @@ -69,6 +69,16 @@ class CurrentState: class Operator(ABC): + def values(self) -> dict[str, float]: + return {} + + def get_values(self) -> dict[str, float]: + out = self.values() + for operator in self.__dict__.values(): + if isinstance(operator, Operator): + out |= operator.get_values() + return out + def __repr__(self) -> str: value_pair_list = list(self.__dict__.items()) if len(value_pair_list) == 0: diff --git a/src/scgenerator/physics/plasma.py b/src/scgenerator/physics/plasma.py index 290684f..cad4f09 100644 --- a/src/scgenerator/physics/plasma.py +++ b/src/scgenerator/physics/plasma.py @@ -3,7 +3,7 @@ import numpy as np import scipy.special from .units import e, hbar, me -from ..math import inverse_integral_exponential, cumulative_simpson +from ..math import expm1_int, cumulative_simpson @dataclass @@ -12,7 +12,7 @@ class PlasmaInfo: dn_dt: np.ndarray polarization: np.ndarray loss: np.ndarray - phase_effect: np.ndarray + debug: np.ndarray class IonizationRate: @@ -64,17 +64,15 @@ class Plasma: field_abs = np.abs(field) delta = 1e-14 * field_abs.max() rate = self.rate(field_abs) - exp_int = inverse_integral_exponential(rate, self.dt) - electron_density = N0 * (1 - exp_int) - dn_dt = N0 * rate * exp_int + exp_int = expm1_int(rate, self.dt) + electron_density = N0 * exp_int + dn_dt = (N0 - electron_density) * rate out = self.dt * cumulative_simpson( dn_dt * self.Ip / (field + delta) + e ** 2 / me * self.dt * cumulative_simpson(electron_density * field) ) loss = cumulative_simpson(dn_dt * self.Ip / (field + delta)) * self.dt - phase_effect = e ** 2 / me * self.dt * cumulative_simpson(electron_density * field) - phase_effect = exp_int - return PlasmaInfo(electron_density, dn_dt, out, loss, phase_effect) + return PlasmaInfo(electron_density, dn_dt, out, loss, loss) def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray: diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py index 90582ee..b8fc1fb 100644 --- a/src/scgenerator/physics/simulate.py +++ b/src/scgenerator/physics/simulate.py @@ -13,6 +13,7 @@ from ..logger import get_logger from ..operators import CurrentState from ..parameter import Configuration, Parameters from ..pbar import PBars, ProgressBarActor, progress_worker +from ..const import ONE_2, ONE_3, ONE_6 try: import ray @@ -248,14 +249,16 @@ class RK4IP: while not keep: h = h_next_step - expD = np.exp(h / 2 * self.params.linear_operator(self.state)) + expD = np.exp(h * ONE_2 * self.params.linear_operator(self.state)) A_I = expD * self.state.spectrum k1 = expD * (h * self.params.nonlinear_operator(self.state)) - k2 = h * self.params.nonlinear_operator(self.state.replace(A_I + k1 / 2)) - k3 = h * self.params.nonlinear_operator(self.state.replace(A_I + k2 / 2)) + k2 = h * self.params.nonlinear_operator(self.state.replace(A_I + k1 * ONE_2)) + k3 = h * self.params.nonlinear_operator(self.state.replace(A_I + k2 * ONE_2)) k4 = h * self.params.nonlinear_operator(self.state.replace(expD * (A_I + k3))) - new_state = self.state.replace(expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6) + new_state = self.state.replace( + expD * (A_I + k1 * ONE_6 + k2 * ONE_3 + k3 * ONE_3) + k4 * ONE_6 + ) self.cons_qty[step] = self.params.conserved_quantity(new_state) if self.params.adapt_step_size: @@ -266,8 +269,8 @@ class RK4IP: progress_str = f"step {step} rejected with h = {h:.4e}, doing over" self.logger.debug(progress_str) keep = False - h_next_step = h / 2 - elif cons_qty_change_ok < curr_p_change <= 2 * cons_qty_change_ok: + h_next_step = h * ONE_2 + elif cons_qty_change_ok < curr_p_change <= 2.0 * cons_qty_change_ok: keep = True h_next_step = h / self.size_fac elif curr_p_change < 0.1 * cons_qty_change_ok: