From 0efc51180311a841bd607f02b637eaf088fe7ee7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Thu, 18 Nov 2021 10:09:27 +0100 Subject: [PATCH] plasma is working but overflows everytime --- src/scgenerator/data/materials.toml | 10 ++++++++++ src/scgenerator/evaluator.py | 3 ++- src/scgenerator/operators.py | 8 +++++++- src/scgenerator/physics/plasma.py | 10 +++------- src/scgenerator/physics/simulate.py | 19 +++++++++++++------ src/scgenerator/scripts/__init__.py | 2 +- src/scgenerator/solver.py | 6 ++++-- src/scgenerator/spectra.py | 14 ++++++++++++++ src/scgenerator/utils.py | 8 ++++---- 9 files changed, 58 insertions(+), 22 deletions(-) diff --git a/src/scgenerator/data/materials.toml b/src/scgenerator/data/materials.toml index 1eab050..f79930d 100644 --- a/src/scgenerator/data/materials.toml +++ b/src/scgenerator/data/materials.toml @@ -78,6 +78,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field [helium] a = 0.00346 b = 2.38e-5 +ionization_energy = 3.9393356074281e-18 +atomic_number = 2 [helium.sellmeier] B = [2.16463842e-5, 2.10561127e-7, 4.7509272e-5] @@ -131,6 +133,8 @@ source = "Shelton, D. P., & Rice, J. E. (1994). Measurements and calculations of [neon] a = 0.02135 b = 1.709e-5 +ionization_energy = 3.45501365359425e-18 +atomic_number = 10 [neon.sellmeier] B = [1281450000.0, 22048600000.0] @@ -148,6 +152,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field [argon] a = 0.1355 b = 3.201e-5 +ionization_energy = 2.5249661793774e-18 +atomic_number = 18 [argon.sellmeier] B = [2501410000.0, 500283000.0, 52234300000.0] @@ -203,6 +209,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field [krypton] a = 0.2349 b = 3.978e-5 +ionization_energy = 2.2429831039374e-18 +atomic_number = 36 [krypton.sellmeier] B = [2536370000.0, 2736490000.0, 62080200000.0] @@ -220,6 +228,8 @@ source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field [xenon] a = 0.425 b = 5.105e-5 +ionization_energy = 1.94342415157935 +atomic_number = 54 [xenon.sellmeier] B = [3228690000.0, 3553930000.0, 60676400000.0] diff --git a/src/scgenerator/evaluator.py b/src/scgenerator/evaluator.py index 7eb6c38..6c50675 100644 --- a/src/scgenerator/evaluator.py +++ b/src/scgenerator/evaluator.py @@ -375,10 +375,10 @@ default_rules: list[Rule] = [ Rule("n_op", operators.MarcatiliRefractiveIndex), 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), - Rule("conserved_quantity", operators.conserved_quantity), ] envelope_rules = default_rules + [ @@ -419,6 +419,7 @@ envelope_rules = default_rules + [ Rule("dispersion_op", operators.ConstantPolyDispersion), Rule("dispersion_op", operators.DirectDispersion), Rule("linear_operator", operators.EnvelopeLinearOperator), + Rule("conserved_quantity", operators.conserved_quantity), ] full_field_rules = default_rules + [ diff --git a/src/scgenerator/operators.py b/src/scgenerator/operators.py index 111afac..aa6d444 100644 --- a/src/scgenerator/operators.py +++ b/src/scgenerator/operators.py @@ -801,6 +801,7 @@ class ConstantGamma(AbstractGamma): class Plasma(Operator): mat_plasma: plasma.Plasma gas_op: AbstractGas + ionization_fraction = 0.0 def __init__(self, dt: float, gas_op: AbstractGas): self.gas_op = gas_op @@ -812,7 +813,12 @@ class Plasma(Operator): def __call__(self, state: CurrentState) -> np.ndarray: N0 = self.gas_op.number_density(state) - return self.mat_plasma(state.field, N0) + plasma_info = self.mat_plasma(state.field, N0) + self.ionization_fraction = plasma_info.electron_density[-1] / N0 + return plasma_info.polarization + + def values(self) -> dict[str, float]: + return dict(ionization_fraction=self.ionization_fraction) class NoPlasma(NoOpTime, Plasma): diff --git a/src/scgenerator/physics/plasma.py b/src/scgenerator/physics/plasma.py index 155b0d8..8cf35de 100644 --- a/src/scgenerator/physics/plasma.py +++ b/src/scgenerator/physics/plasma.py @@ -9,10 +9,7 @@ from ..math import expm1_int, cumulative_simpson @dataclass class PlasmaInfo: electron_density: np.ndarray - dn_dt: np.ndarray polarization: np.ndarray - loss: np.ndarray - debug: np.ndarray class IonizationRate: @@ -66,12 +63,11 @@ class Plasma: rate = self.rate(field_abs) electron_density = free_electron_density(rate, self.dt, N0) dn_dt = (N0 - electron_density) * rate - out = self.dt * cumulative_simpson( - dn_dt * self.Ip / (field + delta) + polarization = self.dt * cumulative_simpson( + dn_dt * self.Ip / np.where(field_abs < delta, 1e-14, field) + e ** 2 / me * self.dt * cumulative_simpson(electron_density * field) ) - loss = cumulative_simpson(dn_dt * self.Ip / (field + delta)) * self.dt - return PlasmaInfo(electron_density, dn_dt, out, loss, loss) + return PlasmaInfo(electron_density, polarization) 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 52d2f37..adde3d7 100644 --- a/src/scgenerator/physics/simulate.py +++ b/src/scgenerator/physics/simulate.py @@ -164,12 +164,19 @@ class RK4IP: state = self.init_state.copy() yield len(self.stored_spectra) - 1, state - integrator_args = [ - self.params.compute(a) for a in solver.Integrator.factory_args() if a != "init_state" - ] - integrator = solver.Integrator.create( - self.params.integration_scheme, state, *integrator_args - ) + if self.params.adapt_step_size: + integrator_args = [ + self.params.compute(a) + for a in solver.Integrator.factory_args() + if a != "init_state" + ] + integrator = solver.Integrator.create( + self.params.integration_scheme, state, *integrator_args + ) + else: + integrator = solver.ConstantStepIntegrator( + state, self.params.linear_operator, self.params.nonlinear_operator + ) for state in integrator: diff --git a/src/scgenerator/scripts/__init__.py b/src/scgenerator/scripts/__init__.py index e11938b..57f1bc5 100644 --- a/src/scgenerator/scripts/__init__.py +++ b/src/scgenerator/scripts/__init__.py @@ -261,7 +261,7 @@ def finish_plot(fig: plt.Figure, legend_axes: plt.Axes, all_labels: list[str], p def plot_helper(config_path: Path) -> Iterable[tuple[dict, list[str], Parameters]]: cc = cycler(color=[f"C{i}" for i in range(10)]) * cycler(ls=["-", "--"]) - for style, (descriptor, params) in zip(cc, Configuration(config_path)): + for style, (descriptor, params), _ in zip(cc, Configuration(config_path), range(20)): yield style, descriptor.branch.formatted_descriptor(), params diff --git a/src/scgenerator/solver.py b/src/scgenerator/solver.py index 3c16482..c5d007d 100644 --- a/src/scgenerator/solver.py +++ b/src/scgenerator/solver.py @@ -130,7 +130,7 @@ class ConstantStepIntegrator(Integrator, scheme="constant"): linear_operator: LinearOperator, nonlinear_operator: NonLinearOperator, ): - super().__init__(init_state, linear_operator, nonlinear_operator, 0.0) + super().__init__(init_state, linear_operator, nonlinear_operator) def __iter__(self) -> Iterator[CurrentState]: while True: @@ -156,6 +156,7 @@ class AdaptiveIntegrator(Integrator): target_error: float current_error: float = 0.0 steps_rejected: float = 0 + min_step_size = 0.0 def __init__( self, @@ -190,7 +191,8 @@ class AdaptiveIntegrator(Integrator): else: next_h_factor = 2.0 h_next_step = h * next_h_factor - accepted = self.current_error <= 2 * self.target_error + accepted = self.current_error <= 2 * self.target_error or h_next_step < self.min_step_size + h_next_step = max(h_next_step, self.min_step_size) if not accepted: self.steps_rejected += 1 self.logger.info( diff --git a/src/scgenerator/spectra.py b/src/scgenerator/spectra.py index 31f39c5..7ae51b4 100644 --- a/src/scgenerator/spectra.py +++ b/src/scgenerator/spectra.py @@ -17,6 +17,7 @@ from .plotting import ( mean_values_plot, propagation_plot, single_position_plot, + transform_1D_values, transform_2D_propagation, ) from .utils import load_spectrum, simulations_list, load_toml @@ -246,6 +247,19 @@ class SimulationSeries: vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind) return single_position_plot(vals, plot_range, self.params, ax, **kwargs) + def plot_values_1D( + self, + left: float, + right: float, + unit: Union[Callable[[float], float], str], + z_pos: int, + sim_ind: int = 0, + **kwargs, + ) -> tuple[np.ndarray, np.ndarray]: + plot_range = PlotRange(left, right, unit) + vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind) + return transform_1D_values(vals, plot_range, self.params, **kwargs) + def plot_mean( self, left: float, diff --git a/src/scgenerator/utils.py b/src/scgenerator/utils.py index 221edbd..580cd73 100644 --- a/src/scgenerator/utils.py +++ b/src/scgenerator/utils.py @@ -11,7 +11,7 @@ import os import re from collections import defaultdict from dataclasses import dataclass -from functools import cache +from functools import cache, lru_cache from pathlib import Path from string import printable as str_printable from typing import Any, Callable, MutableMapping, Sequence, TypeVar, Set, Union @@ -21,7 +21,7 @@ import pkg_resources as pkg import tomli import tomli_w -from .const import PARAM_FN, PARAM_SEPARATOR, SPEC1_FN, Z_FN, ROOT_PARAMETERS +from .const import PARAM_FN, PARAM_SEPARATOR, SPEC1_FN, SPECN_FN1, Z_FN, ROOT_PARAMETERS from .logger import get_logger from .errors import DuplicateParameterError @@ -134,7 +134,7 @@ def load_previous_spectrum(prev_data_dir: str) -> np.ndarray: return load_spectrum(prev_data_dir / SPEC1_FN.format(num)) -@cache +@lru_cache(20000) def load_spectrum(file: os.PathLike) -> np.ndarray: return np.load(file) @@ -535,7 +535,7 @@ def simulations_list(path: os.PathLike) -> list[Path]: """ paths: list[Path] = [] for pwd, _, files in os.walk(path): - if PARAM_FN in files: + if PARAM_FN in files and SPEC1_FN.format(0) in files: paths.append(Path(pwd)) paths.sort(key=branch_id) return [p for p in paths if p.parent.name == paths[-1].parent.name]