diff --git a/src/scgenerator/evaluator.py b/src/scgenerator/evaluator.py index cc44e64..e617046 100644 --- a/src/scgenerator/evaluator.py +++ b/src/scgenerator/evaluator.py @@ -386,7 +386,6 @@ default_rules: list[Rule] = [ Rule("n_op", operators.MarcatiliAdjustedRefractiveIndex), Rule("n_op", operators.HasanRefractiveIndex), Rule("gas_op", operators.ConstantGas), - Rule("raman_op", operators.NoRaman, priorities=-1), Rule("loss_op", operators.NoLoss, priorities=-1), Rule("conserved_quantity", operators.NoConservedQuantity, priorities=-1), ] @@ -422,6 +421,7 @@ envelope_rules = default_rules + [ Rule("spm_op", operators.NoEnvelopeSPM, priorities=-1), Rule("spm_op", operators.EnvelopeSPM), Rule("raman_op", operators.EnvelopeRaman), + Rule("raman_op", operators.NoEnvelopeRaman, priorities=-1), Rule("nonlinear_operator", operators.EnvelopeNonLinearOperator), Rule("loss_op", operators.CustomLoss, priorities=3), Rule("loss_op", operators.CapillaryLoss, priorities=2, conditions=dict(loss="capillary")), @@ -454,5 +454,6 @@ full_field_rules = default_rules + [ ), Rule("plasma_op", operators.Plasma, conditions=dict(photoionization=True)), Rule("plasma_op", operators.NoPlasma, priorities=-1), + Rule("raman_op", operators.NoFullFieldRaman, priorities=-1), Rule("nonlinear_operator", operators.FullFieldNonLinearOperator), ] diff --git a/src/scgenerator/operators.py b/src/scgenerator/operators.py index 21fb67f..1b13bfb 100644 --- a/src/scgenerator/operators.py +++ b/src/scgenerator/operators.py @@ -630,7 +630,7 @@ class CustomLoss(ConstantLoss): class AbstractRaman(Operator): - f_r: float = 0.0 + fraction: float = 0.0 @abstractmethod def __call__(self, state: CurrentState) -> np.ndarray: @@ -647,27 +647,34 @@ class AbstractRaman(Operator): """ -class NoRaman(NoOpTime, AbstractRaman): +class NoEnvelopeRaman(NoOpTime, AbstractRaman): + pass + + +class NoFullFieldRaman(NoOpFreq, AbstractRaman): pass class EnvelopeRaman(AbstractRaman): def __init__(self, raman_type: str, t: np.ndarray): self.hr_w = fiber.delayed_raman_w(t, raman_type) - self.f_r = 0.245 if raman_type == "agrawal" else 0.18 + self.fraction = 0.245 if raman_type == "agrawal" else 0.18 def __call__(self, state: CurrentState) -> np.ndarray: - return self.f_r * np.fft.ifft(self.hr_w * np.fft.fft(state.field2)) + return self.fraction * np.fft.ifft(self.hr_w * np.fft.fft(state.field2)) class FullFieldRaman(AbstractRaman): - def __init__(self, raman_type: str, t: np.ndarray, chi3: float): + def __init__(self, raman_type: str, t: np.ndarray, w: np.ndarray, chi3: float): self.hr_w = fiber.delayed_raman_w(t, raman_type) - self.f_r = 0.245 if raman_type == "agrawal" else 0.18 - self.multiplier = units.epsilon0 * chi3 * self.f_r + self.fraction = 0.245 if raman_type == "agrawal" else 0.18 + self.factor_in = units.epsilon0 * chi3 * self.fraction + self.factor_out = 1j * w ** 2 / (2.0 * units.c ** 2 * units.epsilon0) def __call__(self, state: CurrentState) -> np.ndarray: - return self.multiplier * np.fft.ifft(np.fft.fft(state.field2) * self.hr_w) + return self.factor_out * np.fft.rfft( + self.factor_in * state.field * np.fft.irfft(self.hr_w * np.fft.rfft(state.field2)) + ) ################################################## @@ -703,19 +710,20 @@ class NoFullFieldSPM(NoOpTime, AbstractSPM): class EnvelopeSPM(AbstractSPM): def __init__(self, raman_op: AbstractRaman): - self.fraction = 1 - raman_op.f_r + self.fraction = 1 - raman_op.fraction def __call__(self, state: CurrentState) -> np.ndarray: return self.fraction * state.field2 class FullFieldSPM(AbstractSPM): - def __init__(self, raman_op: AbstractRaman, chi3: float): - self.fraction = 1 - raman_op.f_r - self.factor = self.fraction * chi3 * units.epsilon0 + def __init__(self, raman_op: AbstractRaman, w: np.ndarray, chi3: float): + self.fraction = 1 - raman_op.fraction + self.factor_out = 1j * w ** 2 / (2.0 * units.c ** 2 * units.epsilon0) + self.factor_in = self.fraction * chi3 * units.epsilon0 def __call__(self, state: CurrentState) -> np.ndarray: - return self.factor * state.field2 * state.field + return self.factor_out * np.fft.rfft(self.factor_in * state.field2 * state.field) ################################################## @@ -804,19 +812,20 @@ class Plasma(Operator): gas_op: AbstractGas ionization_fraction = 0.0 - def __init__(self, dt: float, gas_op: AbstractGas): + def __init__(self, dt: float, w: np.ndarray, gas_op: AbstractGas): self.gas_op = gas_op self.mat_plasma = plasma.Plasma( dt, self.gas_op.material_dico["ionization_energy"], plasma.IonizationRateADK(self.gas_op.material_dico["ionization_energy"]), ) + self.factor_out = -w / (2.0 * units.c ** 2 * units.epsilon0) def __call__(self, state: CurrentState) -> np.ndarray: N0 = self.gas_op.number_density(state) plasma_info = self.mat_plasma(state.field, N0) self.ionization_fraction = plasma_info.electron_density[-1] / N0 - return plasma_info.polarization + return self.factor_out * np.fft.rfft(plasma_info.dp_dt) def values(self) -> dict[str, float]: return dict(ionization_fraction=self.ionization_fraction) @@ -902,7 +911,7 @@ def conserved_quantity( if not adapt_step_size: return NoConservedQuantity() logger = get_logger(__name__) - raman = not isinstance(raman_op, NoRaman) + raman = not isinstance(raman_op, NoEnvelopeRaman) loss = not isinstance(loss_op, NoLoss) if raman and loss: logger.debug("Conserved quantity : photon number with loss") @@ -1023,4 +1032,4 @@ class FullFieldNonLinearOperator(NonLinearOperator): def __call__(self, state: CurrentState) -> np.ndarray: total_nonlinear = self.spm_op(state) + self.raman_op(state) + self.plasma_op(state) - return self.factor / self.beta_op(state) * np.fft.rfft(total_nonlinear) + return total_nonlinear / self.beta_op(state) diff --git a/src/scgenerator/physics/plasma.py b/src/scgenerator/physics/plasma.py index 4f794de..8276843 100644 --- a/src/scgenerator/physics/plasma.py +++ b/src/scgenerator/physics/plasma.py @@ -12,7 +12,7 @@ T_a = TypeVar("T_a", np.floating, np.ndarray) @dataclass class PlasmaInfo: electron_density: np.ndarray - polarization: np.ndarray + dp_dt: np.ndarray rate: np.ndarray dn_dt: np.ndarray debug: np.ndarray @@ -79,18 +79,15 @@ class Plasma: rate = self.rate(field_abs) electron_density = free_electron_density(rate, self.dt, N0) dn_dt: np.ndarray = (N0 - electron_density) * rate - integrand = np.zeros_like(field) - integrand[nzi] = dn_dt[nzi] * self.ionization_energy / field[nzi] - energy_loss = self.dt * cumulative_simpson(integrand) - added_phase = ( - self.dt ** 2 - * e ** 2 - / me - * cumulative_simpson(cumulative_simpson(electron_density * field)) - ) - polarization = energy_loss + added_phase - return PlasmaInfo(electron_density, polarization, rate, dn_dt, (energy_loss, added_phase)) + loss_term = np.zeros_like(field) + loss_term[nzi] = dn_dt[nzi] * self.ionization_energy / field[nzi] + + phase_term = self.dt * e ** 2 / me * cumulative_simpson(electron_density * field) + + # polarization = -loss_integrated + phase_integrated + dp_dt = loss_term + return PlasmaInfo(electron_density, dp_dt, rate, dn_dt, phase_term) def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray: diff --git a/src/scgenerator/solver.py b/src/scgenerator/solver.py index 60e9586..c6262e2 100644 --- a/src/scgenerator/solver.py +++ b/src/scgenerator/solver.py @@ -19,6 +19,11 @@ from .operators import ( from .utils import get_arg_names import warnings + +class IntegratorError(Exception): + pass + + # warnings.filterwarnings("error")