From 3cfd29821fe845f4706a49a9cb5c7b1e3e0ba666 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 10 Jun 2026 19:26:06 +0000 Subject: [PATCH] Centralize algorithm parameters in grouped Settings with config-file loading Extract magic numbers (grid resolutions, solver tolerances, search heuristics, spline parameters, plotting and IO writer defaults) from equilibrium.py, equi_tools.py, surfaces.py, field_line_tracers.py, plotting.py, geqdsk.py and omas.py into pleque.config.Settings, organized into per-topic sections (grid, splines, critical_points, lcfs, field_line_tracing, flux_surfaces, plotting, io). Function arguments still override settings per call; settings are now read at use time instead of being captured at import time. Settings load from the first existing of: $PLEQUE_CONFIG_FILE, ./pleque.toml, the per-user config dir (%APPDATA%\pleque on Windows, $XDG_CONFIG_HOME/pleque elsewhere), or a [tool.pleque] table in the nearest pyproject.toml, with PLEQUE_-prefixed environment variables taking precedence over files. The winning file is recorded on the settings object (config_source, config_files_consulted) and reload_settings() re-reads the configuration at runtime. The old flat fields moved into sections: npsi_grid -> flux_surfaces.n_psi, psin0 -> flux_surfaces.psi_n_min, nr_grid/nz_grid -> lcfs.search_grid_nr/nz (the None sentinel is gone; defaults are 700/1200 directly). Also fix lcfs_field_line ignoring its vect_no and xp_shift arguments (they were hardcoded in the track_plasma_boundary call). Add tests for defaults, precedence, source tracking, env vars and reload; document the feature in docs/source/configuration.rst (kept in sync by a test) and README. pydantic-settings gains the toml extra for Python 3.10 support. https://claude.ai/code/session_01HRFoBXP3uhqR5cZuURzdhS --- README.md | 18 ++ docs/source/configuration.rst | 269 ++++++++++++++++++++ docs/source/index.rst | 1 + pleque/config/__init__.py | 21 ++ pleque/config/settings.py | 391 +++++++++++++++++++++++++++-- pleque/core/equilibrium.py | 127 ++++++---- pleque/io/geqdsk.py | 20 +- pleque/io/omas.py | 6 +- pleque/utils/equi_tools.py | 76 ++++-- pleque/utils/field_line_tracers.py | 15 +- pleque/utils/plotting.py | 38 +-- pleque/utils/surfaces.py | 26 +- poetry.lock | 21 +- pyproject.toml | 2 +- tests/test_settings.py | 196 +++++++++++++++ 15 files changed, 1095 insertions(+), 132 deletions(-) create mode 100644 docs/source/configuration.rst create mode 100644 tests/test_settings.py diff --git a/README.md b/README.md index 71ac532..a3e6142 100644 --- a/README.md +++ b/README.md @@ -102,6 +102,24 @@ the same spatial shape as the requested coordinates for scalar quantities. * Passing mesh-shaped `R` and `Z` arrays with `grid=False` is treated as elementwise evaluation and preserves the mesh shape. +## Configuration + +All tunable algorithm parameters (grid resolutions, solver tolerances, search heuristics, ...) +have built-in defaults that can be overridden by a `pleque.toml` file — in the working +directory, in the user config directory (`~/.config/pleque/` or `%APPDATA%\pleque\`), or as a +`[tool.pleque]` table in `pyproject.toml` — or by `PLEQUE_`-prefixed environment variables: + +```toml +[flux_surfaces] +n_psi = 300 + +[lcfs] +search_grid_nr = 1000 +``` + +See the [configuration documentation](https://pleque.readthedocs.io/en/latest/configuration.html) +for the file lookup order and the full settings reference. + ## Authors * **Lukáš Kripner** - [kripnerl](https://github.com/kripnerl) diff --git a/docs/source/configuration.rst b/docs/source/configuration.rst new file mode 100644 index 0000000..e647032 --- /dev/null +++ b/docs/source/configuration.rst @@ -0,0 +1,269 @@ +Configuration +============= + +PLEQUE's algorithms are controlled by a set of tunable parameters — grid +resolutions, solver tolerances, search heuristics, plotting defaults. All of +them have sensible built-in defaults, and all of them can be overridden at +three levels: + +1. **per call** — most functions accept the parameter as a keyword argument, +2. **per process** — environment variables with the ``PLEQUE_`` prefix, +3. **persistently** — a ``pleque.toml`` configuration file (or a + ``[tool.pleque]`` table in ``pyproject.toml``). + +Configuration file +------------------ + +On first use PLEQUE looks for a configuration file in the following locations +and loads the **first one that exists** (configuration files are *not* merged — +the remaining locations are ignored): + +.. list-table:: + :header-rows: 1 + :widths: 10 50 40 + + * - Priority + - Location + - Notes + * - 1 + - File pointed to by the ``PLEQUE_CONFIG_FILE`` environment variable + - Raises ``FileNotFoundError`` if set but missing. + * - 2 + - ``./pleque.toml`` + - Current working directory. + * - 3 + - ``pleque.toml`` in the user configuration directory + - Linux/macOS: ``$XDG_CONFIG_HOME/pleque/pleque.toml`` + (``~/.config/pleque/pleque.toml`` by default). + Windows: ``%APPDATA%\pleque\pleque.toml``. + * - 4 + - ``[tool.pleque]`` table in the nearest ``pyproject.toml`` + - Found by walking up from the current working directory; only counts + if the ``[tool.pleque]`` table is actually present. + * - 5 + - Built-in defaults + - The values listed in the reference below. + +Format example +-------------- + +``pleque.toml`` contains one TOML table per settings section; any subset of +keys may be given, everything else keeps its default: + +.. code-block:: toml + + default_cocos = 3 + + [flux_surfaces] + n_psi = 300 + psi_n_min = 0.005 + + [lcfs] + search_grid_nr = 1000 + search_grid_nz = 1500 + + [field_line_tracing] + rtol = 1e-9 + +The equivalent form inside a project's ``pyproject.toml``: + +.. code-block:: toml + + [tool.pleque] + default_cocos = 3 + + [tool.pleque.flux_surfaces] + n_psi = 300 + +.. note:: + Unknown keys in a configuration file are silently ignored (this keeps old + PLEQUE versions compatible with configs written for newer ones) — double + check the spelling of the keys you set. + +Environment variables +--------------------- + +Every setting can also be set through an environment variable with the +``PLEQUE_`` prefix; nested sections are addressed with a double underscore. +Environment variables take precedence over values from configuration files: + +.. code-block:: bash + + export PLEQUE_DEFAULT_COCOS=3 + export PLEQUE_FLUX_SURFACES__N_PSI=300 + export PLEQUE_LCFS__REFINEMENT_TOLERANCE=1e-8 + +Programmatic access +------------------- + +.. code-block:: python + + from pleque.config import get_settings, reload_settings + + settings = get_settings() # cached, process-wide instance + settings.flux_surfaces.n_psi # -> 200 + + settings.config_source # Path of the loaded file, or None + settings.config_files_consulted # all locations checked, in order + + settings = reload_settings() # re-read files/env (clears the cache) + +``get_settings()`` is cached (``functools.lru_cache``), so PLEQUE reads the +configuration once per process; call ``reload_settings()`` after changing a +configuration file or environment variable at runtime. A fresh, uncached +instance can be built with ``pleque.config.load_settings()``, and ad-hoc +overrides can be passed directly: ``Settings(flux_surfaces={"n_psi": 50})``. + +Settings reference +------------------ + +Top level +^^^^^^^^^ + +=================== ======= ========= ========================================================== +Name Type Default Description +=================== ======= ========= ========================================================== +``default_cocos`` int 3 COCOS convention assumed when the input does not specify one. +=================== ======= ========= ========================================================== + +``[grid]`` — default computational grids +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +==================================== ======= ========= ========================================================== +Name Type Default Description +==================================== ======= ========= ========================================================== +``default_nr`` int 1000 Number of R points of the default ``Equilibrium.grid()`` grid. +``default_nz`` int 2000 Number of Z points of the default ``Equilibrium.grid()`` grid. +``synthetic_wall_margin`` float 0.01 Fractional inset of the rectangular first wall synthesized when no wall is given. +``synthetic_wall_points_per_side`` int 20 Number of points per side of the synthesized rectangular first wall. +==================================== ======= ========= ========================================================== + +``[splines]`` — spline construction +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +==================== ======= ========= ========================================================== +Name Type Default Description +==================== ======= ========= ========================================================== +``psi_order`` int 3 Order of the 2D psi(R, Z) spline. +``psi_smooth`` float 0.0 Smoothing factor ``s`` of the 2D psi(R, Z) spline. +``profile_order`` int 3 Order ``k`` of 1D profile splines (F, pressure, q, ...). +``profile_smooth`` float 0.0 Smoothing factor ``s`` of 1D profile splines. +==================== ======= ========= ========================================================== + +``[critical_points]`` — magnetic axis and X-point search +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +====================================== ======= ========= ========================================================== +Name Type Default Description +====================================== ======= ========= ========================================================== +``find_extremes_order`` int 20 Number of neighbouring grid points used by ``argrelmin`` when locating extremes of \|grad psi\|. +``find_extremes_max_iter`` int 10 Maximum number of retries (with reduced order) when no O-point candidate is found. +``gradient_threshold`` float 1.0 Upper bound on \|grad psi\|^2 for a grid point to count as a critical-point candidate. +``vicinity_radius`` float 0.1 Half-size [m] of the (R, Z) box used when refining a critical point by local minimization. +``minimizer_xtol`` float 1e-7 ``xtol`` passed to the Powell/TNC minimizers refining extremes. +``relocation_threshold`` float 1e-2 If unbounded minimization moves a candidate point by more than this, the bounded minimizer is used instead. +``axis_vertical_weight`` float 5.0 Weight of the radial distance when ranking O-point candidates. +``axis_out_of_wall_penalty`` float 1e-3 Score multiplier penalizing O-point candidates outside the first wall. +``x_point_sort_epsilon`` float 1e-3 Regularization added to X-point ranking terms to avoid zero scores. +``monotonicity_test_points`` int 10 Number of test points for the psi-monotonicity check between magnetic axis and X-point. +``limiter_monotonicity_test_points`` int 50 Number of test points for the psi-monotonicity check between magnetic axis and limiter point. +====================================== ======= ========= ========================================================== + +``[lcfs]`` — last closed flux surface search +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +============================ ======= ========= ========================================================== +Name Type Default Description +============================ ======= ========= ========================================================== +``search_grid_nr`` int 700 Number of R points of the grid used to locate the LCFS contour. +``search_grid_nz`` int 1200 Number of Z points of the grid used to locate the LCFS contour. +``refinement_tolerance`` float 1e-10 Relative psi error under which the iterative LCFS refinement stops. +``inner_psi_n_offset`` float 1e-5 Offset used as psi_n = 1 - offset when contouring the LCFS as an inner flux surface. +``initial_psi_offset`` float 1e-4 Fraction of (psi_lcfs - psi_axis) used to shift psi inwards for the initial LCFS candidate. +``surface_step_damping`` float 0.99 Damping factor of the gradient step in the iterative flux-surface refinement. +``x_point_shift`` float 1e-6 Distance [m] by which the separatrix tracing start point is shifted from the X-point. +============================ ======= ========= ========================================================== + +``[field_line_tracing]`` — field-line tracing (ODE solver) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +============================== ======= =========== ========================================================== +Name Type Default Description +============================== ======= =========== ========================================================== +``atol`` float 1e-6 Absolute tolerance of the field-line ODE solver. +``rtol`` float 1e-8 Relative tolerance of the field-line ODE solver. +``x_point_atol_scale`` float 1e-3 For X-point plasmas the absolute tolerance is reduced to ``min(atol, distance_to_x_point * x_point_atol_scale)``. +``max_step`` float 1e-2 Maximum integration step in toroidal angle phi [rad]. +``max_toroidal_turns`` int 50 Maximum number of toroidal turns integrated when tracing. +``poloidal_stop_resolution`` float pi/1024 Poloidal-angle resolution [rad] of the stopper used by ``Equilibrium.trace_field_line``. +``default_stop_resolution`` float pi/360 Default poloidal-angle resolution [rad] of ``poloidal_angle_stopper_factory``. +``target_stopper_atol`` float 1e-6 Distance tolerance [m] of the target-point stopper used when tracing flux surfaces. +============================== ======= =========== ========================================================== + +``[flux_surfaces]`` — flux-surface generation and 1D profiles +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +========================= ======= ========= ========================================================== +Name Type Default Description +========================= ======= ========= ========================================================== +``n_psi`` int 200 Default number of psi_n levels used to evaluate flux surfaces. +``psi_n_min`` float 0.01 Minimum psi_n used for flux-surface averaging and volume integration. +``contour_step`` float 1e-3 Spatial resolution [m] of flux-surface contours found on the (R, Z) grid. +``contour_grid_nr`` int 100 Default number of R points of the grid used for flux-surface contouring. +``contour_grid_nz`` int 100 Default number of Z points of the grid used for flux-surface contouring. +``trace_step`` float 1e-3 Maximum step [m] along the flux surface used by ``Equilibrium.trace_flux_surface``. +``trace_max_turns`` int 4 Maximum number of toroidal turns integrated when tracing a flux surface. +``q_psi_n_min`` float 0.01 Lowest psi_n of the grid on which the q profile is evaluated. +``q_psi_n_step`` float 0.005 psi_n step of the grid on which the q profile is evaluated. +``q_search_psi_n_max`` float 0.95 Default upper psi_n bound when searching a flux surface with a given q. +``q_search_psi_n_cap`` float 0.99 Hard upper psi_n bound of the root finder used when searching a flux surface with a given q. +``midplane_map_points`` int 100 Number of points of the midplane r_mid -> psi mapping spline. +========================= ======= ========= ========================================================== + +``[plotting]`` — plotting helpers +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +============================ ======= ========= ========================================================== +Name Type Default Description +============================ ======= ========= ========================================================== +``grid_nr`` int 400 Number of R points of grids used by plotting routines. +``grid_nz`` int 600 Number of Z points of grids used by plotting routines. +``debug_contour_levels`` int 60 Number of psi contour levels in the debug overview plot. +``psi_contour_levels`` int 20 Default number of contour levels of ``plot_psi_contours``. +``rational_q_psi_n_max`` float 0.95 Upper psi_n bound when looking up rational q surfaces for plotting. +``sol_spacing`` float 2e-3 Default radial spacing [m] of near-SOL contours. +``axis_margin`` float 1/12 Margin around the first wall, as a fraction of its size, used for axis limits. +============================ ======= ========= ========================================================== + +``[io]`` — file writers +^^^^^^^^^^^^^^^^^^^^^^^ + +==================== ======= ========= ========================================================== +Name Type Default Description +==================== ======= ========= ========================================================== +``geqdsk_nx`` int 64 Default number of R points of a written G-EQDSK file. +``geqdsk_ny`` int 128 Default number of Z points of a written G-EQDSK file. +``geqdsk_nbdry`` int 200 Default number of boundary points of a written G-EQDSK file. +``omas_n_psi`` int 200 Default number of 1D profile points written to an OMAS data structure. +``omas_grid_step`` float 1e-3 Default (R, Z) grid step [m] of 2D profiles written to an OMAS data structure. +==================== ======= ========= ========================================================== + +Migration from older versions +----------------------------- + +Up to PLEQUE 0.0.8 the (undocumented) ``Settings`` class had four flat fields. +They moved into sections: + +==================== ================================ +Old name New name +==================== ================================ +``npsi_grid`` ``flux_surfaces.n_psi`` +``psin0`` ``flux_surfaces.psi_n_min`` +``nr_grid`` ``lcfs.search_grid_nr`` +``nz_grid`` ``lcfs.search_grid_nz`` +==================== ================================ + +``nr_grid``/``nz_grid`` previously defaulted to ``None`` ("use the fallback +700/1200"); the new fields default to 700/1200 directly. Unprefixed +environment variables (e.g. ``NPSI_GRID``) are no longer read — use the +``PLEQUE_`` prefix. diff --git a/docs/source/index.rst b/docs/source/index.rst index 68689f4..72b52ab 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -15,6 +15,7 @@ PLEQUE is a code which allows easy and quick access to tokamak equilibria obtain coordinates flux_expansion naming_convention + configuration API Reference diff --git a/pleque/config/__init__.py b/pleque/config/__init__.py index e69de29..a31ccdc 100644 --- a/pleque/config/__init__.py +++ b/pleque/config/__init__.py @@ -0,0 +1,21 @@ +from pleque.config.settings import ( + CONFIG_FILE_ENV_VAR, + CONFIG_FILE_NAME, + Settings, + get_settings, + load_settings, + reload_settings, + resolve_config_file, + user_config_path, +) + +__all__ = [ + "CONFIG_FILE_ENV_VAR", + "CONFIG_FILE_NAME", + "Settings", + "get_settings", + "load_settings", + "reload_settings", + "resolve_config_file", + "user_config_path", +] diff --git a/pleque/config/settings.py b/pleque/config/settings.py index fc04492..7604349 100644 --- a/pleque/config/settings.py +++ b/pleque/config/settings.py @@ -1,28 +1,387 @@ +""" +Central configuration for PLEQUE. + +All tunable algorithm parameters (grid resolutions, solver tolerances, search +heuristics, plotting defaults, ...) live in the :class:`Settings` model below, +grouped into per-topic sections. Values are resolved with the following +precedence (highest first): + +1. keyword arguments passed directly to ``Settings(...)``, +2. environment variables (``PLEQUE_`` prefix, ``__`` as section delimiter, + e.g. ``PLEQUE_FLUX_SURFACES__N_PSI=300``), +3. a single configuration file -- the first one found of: + + a. the file pointed to by the ``PLEQUE_CONFIG_FILE`` environment variable, + b. ``pleque.toml`` in the current working directory, + c. ``pleque.toml`` in the user configuration directory + (``%APPDATA%\\pleque\\`` on Windows, ``$XDG_CONFIG_HOME/pleque/`` or + ``~/.config/pleque/`` elsewhere), + d. a ``[tool.pleque]`` table in the nearest ``pyproject.toml`` found by + walking up from the current working directory, + +4. the hard-coded defaults defined in this module. + +Configuration files are *not* merged: the first existing file wins and the +remaining locations are ignored. The winning file (if any) is recorded on the +settings object as :attr:`Settings.config_source`. + +See ``docs/source/configuration.rst`` for the full reference. +""" + +import math +import os +import sys from functools import lru_cache +from pathlib import Path + +from pydantic import BaseModel, Field, PrivateAttr +from pydantic_settings import ( + BaseSettings, + PydanticBaseSettingsSource, + SettingsConfigDict, + TomlConfigSettingsSource, +) + +if sys.version_info >= (3, 11): + import tomllib +else: + import tomli as tomllib + +CONFIG_FILE_NAME = "pleque.toml" +CONFIG_FILE_ENV_VAR = "PLEQUE_CONFIG_FILE" + + +class GridSettings(BaseModel): + """Default computational (R, Z) grids.""" + + default_nr: int = Field(1000, description="Number of R points of the default `Equilibrium.grid()` grid.") + default_nz: int = Field(2000, description="Number of Z points of the default `Equilibrium.grid()` grid.") + synthetic_wall_margin: float = Field( + 0.01, + description="Fractional inset of the rectangular first wall synthesized when no wall is given, " + "relative to the computational grid size.", + ) + synthetic_wall_points_per_side: int = Field( + 20, description="Number of points per side of the synthesized rectangular first wall." + ) + + +class SplineSettings(BaseModel): + """Spline construction parameters.""" + + psi_order: int = Field(3, description="Order of the 2D psi(R, Z) `RectBivariateSpline`.") + psi_smooth: float = Field(0.0, description="Smoothing factor `s` of the 2D psi(R, Z) spline.") + profile_order: int = Field( + 3, description="Order `k` of 1D profile splines (F, pressure, q, ... as functions of psi_n)." + ) + profile_smooth: float = Field(0.0, description="Smoothing factor `s` of 1D profile splines.") + + +class CriticalPointSettings(BaseModel): + """Search for critical points of psi (magnetic axis, X-points).""" + + find_extremes_order: int = Field( + 20, description="Number of neighbouring grid points used by `argrelmin` when locating extremes of |grad psi|." + ) + find_extremes_max_iter: int = Field( + 10, description="Maximum number of retries (with reduced order) when no O-point candidate is found." + ) + gradient_threshold: float = Field( + 1.0, description="Upper bound on |grad psi|^2 for a grid point to count as a critical-point candidate." + ) + vicinity_radius: float = Field( + 0.1, description="Half-size [m] of the (R, Z) box used when refining a critical point by local minimization." + ) + minimizer_xtol: float = Field(1e-7, description="`xtol` passed to the Powell/TNC minimizers refining extremes.") + relocation_threshold: float = Field( + 1e-2, + description="If unbounded minimization moves a candidate point by more than this (in psi-grid units), " + "the bounded minimizer is used instead.", + ) + axis_vertical_weight: float = Field( + 5.0, + description="Weight of the radial distance (relative to the vertical one) when ranking O-point candidates " + "by distance from the expected magnetic-axis position.", + ) + axis_out_of_wall_penalty: float = Field( + 1e-3, description="Score multiplier penalizing O-point candidates located outside the first wall." + ) + x_point_sort_epsilon: float = Field( + 1e-3, description="Regularization added to X-point ranking terms to avoid zero scores." + ) + monotonicity_test_points: int = Field( + 10, description="Number of test points for the psi-monotonicity check between magnetic axis and X-point." + ) + limiter_monotonicity_test_points: int = Field( + 50, description="Number of test points for the psi-monotonicity check between magnetic axis and limiter point." + ) + + +class LcfsSettings(BaseModel): + """Last closed flux surface search and refinement.""" + + search_grid_nr: int = Field(700, description="Number of R points of the grid used to locate the LCFS contour.") + search_grid_nz: int = Field(1200, description="Number of Z points of the grid used to locate the LCFS contour.") + refinement_tolerance: float = Field( + 1e-10, description="Relative psi error under which the iterative LCFS refinement stops." + ) + inner_psi_n_offset: float = Field( + 1e-5, description="Offset used as psi_n = 1 - offset when contouring the LCFS as an inner flux surface." + ) + initial_psi_offset: float = Field( + 1e-4, + description="Fraction of (psi_lcfs - psi_axis) used to shift psi inwards when searching " + "for the initial closed LCFS candidate.", + ) + surface_step_damping: float = Field( + 0.99, description="Damping factor of the gradient step in the iterative flux-surface refinement." + ) + x_point_shift: float = Field( + 1e-6, description="Distance [m] by which the separatrix tracing start point is shifted from the X-point." + ) + + +class FieldLineTracingSettings(BaseModel): + """Field-line tracing (`scipy.integrate.solve_ivp`) parameters.""" + + atol: float = Field(1e-6, description="Absolute tolerance of the field-line ODE solver.") + rtol: float = Field(1e-8, description="Relative tolerance of the field-line ODE solver.") + x_point_atol_scale: float = Field( + 1e-3, + description="For X-point plasmas the absolute tolerance is reduced to " + "min(atol, distance_to_x_point * x_point_atol_scale).", + ) + max_step: float = Field(1e-2, description="Maximum integration step in toroidal angle phi [rad].") + max_toroidal_turns: int = Field(50, description="Maximum number of toroidal turns integrated when tracing.") + poloidal_stop_resolution: float = Field( + math.pi / 1024, + description="Poloidal-angle resolution [rad] of the stopper used by `Equilibrium.trace_field_line`.", + ) + default_stop_resolution: float = Field( + math.pi / 360, description="Default poloidal-angle resolution [rad] of `poloidal_angle_stopper_factory`." + ) + target_stopper_atol: float = Field( + 1e-6, description="Distance tolerance [m] of the target-point stopper used when tracing flux surfaces." + ) + -from pydantic import Field -from pydantic_settings import BaseSettings +class FluxSurfaceSettings(BaseModel): + """Flux-surface generation, tracing and derived 1D profiles.""" + + n_psi: int = Field(200, description="Default number of psi_n levels used to evaluate flux surfaces.") + psi_n_min: float = Field(0.01, description="Minimum psi_n used for flux-surface averaging and volume integration.") + contour_step: float = Field( + 1e-3, description="Spatial resolution [m] of flux-surface contours found on the (R, Z) grid." + ) + contour_grid_nr: int = Field( + 100, description="Default number of R points of the grid used for flux-surface contouring." + ) + contour_grid_nz: int = Field( + 100, description="Default number of Z points of the grid used for flux-surface contouring." + ) + trace_step: float = Field( + 1e-3, description="Maximum step [m] along the flux surface used by `Equilibrium.trace_flux_surface`." + ) + trace_max_turns: int = Field( + 4, description="Maximum number of toroidal turns integrated when tracing a flux surface." + ) + q_psi_n_min: float = Field(0.01, description="Lowest psi_n of the grid on which the q profile is evaluated.") + q_psi_n_step: float = Field(0.005, description="psi_n step of the grid on which the q profile is evaluated.") + q_search_psi_n_max: float = Field( + 0.95, description="Default upper psi_n bound when searching a flux surface with a given q." + ) + q_search_psi_n_cap: float = Field( + 0.99, description="Hard upper psi_n bound of the root finder used when searching a flux surface with a given q." + ) + midplane_map_points: int = Field(100, description="Number of points of the midplane r_mid -> psi mapping spline.") + + +class PlottingSettings(BaseModel): + """Default parameters of plotting helpers.""" + + grid_nr: int = Field(400, description="Number of R points of grids used by plotting routines.") + grid_nz: int = Field(600, description="Number of Z points of grids used by plotting routines.") + debug_contour_levels: int = Field(60, description="Number of psi contour levels in the debug overview plot.") + psi_contour_levels: int = Field(20, description="Default number of contour levels of `plot_psi_contours`.") + rational_q_psi_n_max: float = Field( + 0.95, description="Upper psi_n bound when looking up rational q surfaces for plotting." + ) + sol_spacing: float = Field(2e-3, description="Default radial spacing [m] of near-SOL contours.") + axis_margin: float = Field( + 1 / 12, description="Margin around the first wall, as a fraction of its size, used for axis limits." + ) + + +class IOSettings(BaseModel): + """Default parameters of file writers.""" + + geqdsk_nx: int = Field(64, description="Default number of R points of a written G-EQDSK file.") + geqdsk_ny: int = Field(128, description="Default number of Z points of a written G-EQDSK file.") + geqdsk_nbdry: int = Field(200, description="Default number of boundary points of a written G-EQDSK file.") + omas_n_psi: int = Field(200, description="Default number of 1D profile points written to an OMAS data structure.") + omas_grid_step: float = Field( + 1e-3, description="Default (R, Z) grid step [m] of 2D profiles written to an OMAS data structure." + ) + + +def _windows_config_dir() -> Path: + """Per-user config dir on Windows: ``%APPDATA%\\pleque`` (``~\\pleque`` if ``APPDATA`` is unset).""" + appdata = os.environ.get("APPDATA") + base = Path(appdata) if appdata else Path.home() + return base / "pleque" + + +def _posix_config_dir() -> Path: + """Per-user config dir elsewhere: ``$XDG_CONFIG_HOME/pleque``, defaulting to ``~/.config/pleque``.""" + xdg = os.environ.get("XDG_CONFIG_HOME") + base = Path(xdg) if xdg else Path.home() / ".config" + return base / "pleque" + + +def user_config_path() -> Path: + """Return the per-user config file path (not necessarily existing).""" + config_dir = _windows_config_dir() if os.name == "nt" else _posix_config_dir() + return config_dir / CONFIG_FILE_NAME + + +def _load_pyproject_table(path: Path) -> dict: + """Return the [tool.pleque] table of a pyproject.toml file ({} if absent or unreadable).""" + try: + with path.open("rb") as f: + data = tomllib.load(f) + except (OSError, tomllib.TOMLDecodeError): + return {} + return data.get("tool", {}).get("pleque", {}) + + +def _find_project_pyproject() -> Path | None: + """Walk up from cwd to the first pyproject.toml containing a [tool.pleque] table.""" + cwd = Path.cwd() + for directory in (cwd, *cwd.parents): + candidate = directory / "pyproject.toml" + if candidate.is_file(): + return candidate if _load_pyproject_table(candidate) else None + return None + + +class _TomlTableSettingsSource(PydanticBaseSettingsSource): + """Settings source backed by an in-memory mapping (a TOML table read beforehand).""" + + def __init__(self, settings_cls: type[BaseSettings], data: dict): + super().__init__(settings_cls) + self._data = data + + def get_field_value(self, field, field_name): + return self._data.get(field_name), field_name, False + + def __call__(self) -> dict: + return dict(self._data) + + +def resolve_config_file() -> tuple[Path | None, tuple[Path, ...]]: + """Resolve the configuration file to load. + + Returns ``(winner, consulted)`` where ``winner`` is the first existing + config file in precedence order (or ``None`` if none exists) and + ``consulted`` lists all locations that were checked. + + Raises :class:`FileNotFoundError` if ``PLEQUE_CONFIG_FILE`` is set but + does not point to an existing file. + """ + explicit = os.environ.get(CONFIG_FILE_ENV_VAR) + if explicit: + path = Path(explicit) + if not path.is_file(): + raise FileNotFoundError(f"{CONFIG_FILE_ENV_VAR} points to a non-existent file: {path}") + return path, (path,) + + consulted = [Path.cwd() / CONFIG_FILE_NAME, user_config_path()] + winner = next((p for p in consulted if p.is_file()), None) + if winner is None: + pyproject = _find_project_pyproject() + if pyproject is not None: + consulted.append(pyproject) + winner = pyproject + return winner, tuple(consulted) class Settings(BaseSettings): """ - Configuration settings for PLEQUE application. + PLEQUE configuration, grouped into per-topic sections. + All defaults can be overridden by a ``pleque.toml`` file, a + ``[tool.pleque]`` table in ``pyproject.toml``, environment variables + (``PLEQUE_`` prefix, ``__`` section delimiter) or constructor keyword + arguments — see the module docstring for the precedence rules. """ - npsi_grid: int = Field(200, description="Default number of flux labels used for evaluation of flux surfaces.") - nr_grid: int | None = Field(None, - description= - """Default number of radial points of rectangular grid used for various calculations. - If None, the input grid is used.""" - ) - nz_grid: int | None = Field(None, - description= - """Default number of vertical points of rectangular grid used for various calculations. - If None, the input grid is used.""" + + model_config = SettingsConfigDict( + env_prefix="PLEQUE_", + env_nested_delimiter="__", ) - psin0: float = Field(0.01, description="Minimum flux surface psi_n used for flux surface averaging and volume integration.") + + default_cocos: int = Field(3, description="COCOS convention assumed when the input does not specify one.") + grid: GridSettings = Field(default_factory=GridSettings) + splines: SplineSettings = Field(default_factory=SplineSettings) + critical_points: CriticalPointSettings = Field(default_factory=CriticalPointSettings) + lcfs: LcfsSettings = Field(default_factory=LcfsSettings) + field_line_tracing: FieldLineTracingSettings = Field(default_factory=FieldLineTracingSettings) + flux_surfaces: FluxSurfaceSettings = Field(default_factory=FluxSurfaceSettings) + plotting: PlottingSettings = Field(default_factory=PlottingSettings) + io: IOSettings = Field(default_factory=IOSettings) + + _config_source: Path | None = PrivateAttr(default=None) + _config_files_consulted: tuple[Path, ...] = PrivateAttr(default=()) + + @property + def config_source(self) -> Path | None: + """The configuration file values were loaded from, or None for built-in defaults.""" + return self._config_source + + @property + def config_files_consulted(self) -> tuple[Path, ...]: + """All configuration file locations checked during loading, in precedence order.""" + return self._config_files_consulted + + @classmethod + def settings_customise_sources( + cls, + settings_cls: type[BaseSettings], + init_settings: PydanticBaseSettingsSource, + env_settings: PydanticBaseSettingsSource, + dotenv_settings: PydanticBaseSettingsSource, + file_secret_settings: PydanticBaseSettingsSource, + ) -> tuple[PydanticBaseSettingsSource, ...]: + winner, _ = resolve_config_file() + sources: list[PydanticBaseSettingsSource] = [init_settings, env_settings] + if winner is not None: + if winner.name == "pyproject.toml": + sources.append(_TomlTableSettingsSource(settings_cls, _load_pyproject_table(winner))) + else: + sources.append(TomlConfigSettingsSource(settings_cls, winner)) + sources.append(file_secret_settings) + return tuple(sources) + + +def load_settings() -> Settings: + """Build a fresh :class:`Settings` instance (bypassing the cache), recording its source.""" + # The file is resolved a second time inside settings_customise_sources; both + # resolutions are a few stat() calls and see the same file except for a + # negligible race window. + winner, consulted = resolve_config_file() + settings = Settings() + settings._config_source = winner + settings._config_files_consulted = consulted + return settings @lru_cache -def get_settings(): - return Settings() +def get_settings() -> Settings: + """Return the process-wide cached :class:`Settings` instance.""" + return load_settings() + + +def reload_settings() -> Settings: + """Clear the settings cache and reload from configuration files and environment.""" + get_settings.cache_clear() + return get_settings() diff --git a/pleque/core/equilibrium.py b/pleque/core/equilibrium.py index af074eb..dc17f2b 100644 --- a/pleque/core/equilibrium.py +++ b/pleque/core/equilibrium.py @@ -17,7 +17,6 @@ from pleque.utils.surfaces import track_plasma_boundary from pleque.utils.tools import arglis, xp_sections -settings = get_settings() class Equilibrium: """ @@ -44,9 +43,9 @@ def __init__(self, x_points=None, strike_points=None, init_method="hints", - spline_order=3, - spline_smooth=0, - find_extremes_order=20, + spline_order=None, + spline_smooth=None, + find_extremes_order=None, cocos=None, verbose=False, ): @@ -68,10 +67,10 @@ def __init__(self, If "hints" module use given optional arguments as a help with initialization. If "fast-forward" module use given optional arguments as final and doesn't try to correct. *Note:* Only "hints" method is currently tested. - :param spline_order: - :param spline_smooth: + :param spline_order: Order of the 2D psi spline. Defaults to the value from PLEQUE settings. + :param spline_smooth: Smoothing factor of the 2D psi spline. Defaults to the value from PLEQUE settings. :param find_extremes_order: Define number of points on internal grid used for identifying - magnetic axis and x-points. + magnetic axis and x-points. Defaults to the value from PLEQUE settings. :param cocos: At the moment module assume cocos to be 3 (no other option). The implemetnation is not fully working. Be aware of signs in the module! :param verbose: @@ -83,13 +82,21 @@ def __init__(self, print('Equilibrium module initialization') print('---------------------------------') + cfg = get_settings() + if spline_order is None: + spline_order = cfg.splines.psi_order + if spline_smooth is None: + spline_smooth = cfg.splines.psi_smooth + if find_extremes_order is None: + find_extremes_order = cfg.critical_points.find_extremes_order + if cocos is None: if "cocos" in basedata.attrs: cocos = basedata.attrs["cocos"] elif "cocos" in basedata: cocos = basedata["cocos"] else: - cocos = 3 + cocos = cfg.default_cocos self._basedata = basedata self._verbose = verbose @@ -173,10 +180,11 @@ def _load_spatial_data(self, basedata: xarray.Dataset, first_wall) -> tuple: # todo: remove this if possible # lets reduce the wall a bit to be have some plasma behind the wall - rwall_min += dr / 100 - rwall_max -= dr / 100 - zwall_min += dz / 100 - zwall_max -= dz / 100 + wall_cfg = get_settings().grid + rwall_min += dr * wall_cfg.synthetic_wall_margin + rwall_max -= dr * wall_cfg.synthetic_wall_margin + zwall_min += dz * wall_cfg.synthetic_wall_margin + zwall_max -= dz * wall_cfg.synthetic_wall_margin corners = np.array( [[rwall_min, zwall_max], [rwall_max, zwall_max], [rwall_max, zwall_min], @@ -184,8 +192,8 @@ def _load_spatial_data(self, basedata: xarray.Dataset, first_wall) -> tuple: newwall_r = [] newwall_z = [] for i in range(-1, 3): - rs = np.linspace(corners[i, 0], corners[i + 1, 0], 20) - zs = np.linspace(corners[i, 1], corners[i + 1, 1], 20) + rs = np.linspace(corners[i, 0], corners[i + 1, 0], wall_cfg.synthetic_wall_points_per_side) + zs = np.linspace(corners[i, 1], corners[i + 1, 1], wall_cfg.synthetic_wall_points_per_side) newwall_r += list(rs) newwall_z += list(zs) self._first_wall = np.stack((newwall_r, newwall_z)).T @@ -304,8 +312,9 @@ def _find_critical_points(self, r, z, find_extremes_order: int): def _setup_plasma_boundary(self): """Find LCFS contour and strike points using the recognized plasma type.""" - nr = settings.nr_grid if settings.nr_grid is not None else 700 - nz = settings.nz_grid if settings.nz_grid is not None else 1200 + lcfs_cfg = get_settings().lcfs + nr = lcfs_cfg.search_grid_nr + nz = lcfs_cfg.search_grid_nz rs = np.linspace(self.R_min, self.R_max, nr) zs = np.linspace(self.Z_min, self.Z_max, nz) @@ -327,7 +336,7 @@ def _setup_plasma_boundary(self): close_lcfs = eq_tools.find_close_lcfs(self._psi_lcfs, rs, zs, self._spl_psi, self._mg_axis, self._psi_axis) - while surf.fluxsurf_error(self._spl_psi, close_lcfs, self._psi_lcfs) > 1e-10: + while surf.fluxsurf_error(self._spl_psi, close_lcfs, self._psi_lcfs) > lcfs_cfg.refinement_tolerance: close_lcfs = eq_tools.find_surface_step(self._spl_psi, self._psi_lcfs, close_lcfs) if self._verbose: @@ -361,22 +370,23 @@ def _setup_1d_profiles(self, psi_n, F, FFprime, pressure, pprime): F = np.zeros_like(psi_n) self._vacuum = True - self._fpol_spl = UnivariateSpline(psi_n, F, k=3, s=0) + spl_cfg = get_settings().splines + self._fpol_spl = UnivariateSpline(psi_n, F, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) if FFprime is None: self._df_dpsin_spl = self._fpol_spl.derivative() Fprime = self._df_dpsin_spl(psi_n) / (self._psi_lcfs - self._psi_axis) FFprime = F * Fprime - self._pressure_spl = UnivariateSpline(psi_n, pressure, k=3, s=0) + self._pressure_spl = UnivariateSpline(psi_n, pressure, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) if pprime is None: self._dp_dpsin_spl = self._pressure_spl.derivative() pprime = self._dp_dpsin_spl(psi_n) / (self._psi_lcfs - self._psi_axis) - self._pprime_spl = UnivariateSpline(psi_n, pprime, k=3, s=0) - self._Fprime_spl = UnivariateSpline(psi_n, Fprime, k=3, s=0) - self._FFprime_spl = UnivariateSpline(psi_n, FFprime, k=3, s=0) + self._pprime_spl = UnivariateSpline(psi_n, pprime, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) + self._Fprime_spl = UnivariateSpline(psi_n, Fprime, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) + self._FFprime_spl = UnivariateSpline(psi_n, FFprime, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) self.fluxfuncs.add_flux_func('F', F, psi_n=psi_n) self.fluxfuncs.add_flux_func('FFprime', FFprime, psi_n=psi_n) @@ -605,9 +615,12 @@ def Bvec_norm(self, *coordinates, swap_order=False, R=None, Z=None, coord_type=N # XXXXXX TODO TODO TODO @deprecated('The structure and behaviour of this function will change soon!\n' 'to keep the same behaviour use `_flux_surface` instead.') - def flux_surface(self, *coordinates, resolution=(1e-3, 1e-3), dim="step", + def flux_surface(self, *coordinates, resolution=None, dim="step", closed=True, inlcfs=True, R=None, Z=None, psi_n=None, coord_type=None, **coords): + if resolution is None: + contour_step = get_settings().flux_surfaces.contour_step + resolution = (contour_step, contour_step) return self._flux_surface(*coordinates, resolution=resolution, dim=dim, closed=closed, inlcfs=inlcfs, R=R, Z=Z, psi_n=psi_n, coord_type=coord_type, **coords) @@ -650,7 +663,7 @@ def _flux_surface(self, *coordinates, resolution=None, dim="step", # todo: to get lcfs, here is small trick. This should be handled better # otherwise it may return crossed loop if np.isclose(coordinates.psi_n[0], 1) and inlcfs: - psi_n = 1 - 1e-5 + psi_n = 1 - get_settings().lcfs.inner_psi_n_offset else: psi_n = coordinates.psi_n[0] @@ -1161,7 +1174,8 @@ def grid(self, resolution=None, dim="step"): :param resolution: Iterable of size 2 or a number. If a number is passed, R and Z dimensions will have the same size or step (depending on dim parameter). Different R and Z resolutions or dimension sizes can be required by passing an iterable of size 2. - If None, default grid of size (1000, 2000) is returned. + If None, the default grid given by the PLEQUE settings (`grid.default_nr`, + `grid.default_nz`) is returned. :param dim: iterable of size 2 or string ('step', 'size'). Default is "step", determines the meaning of the resolution. If "step" used, values in resolution are interpreted as step length in psi poloidal map. If "size" is used, @@ -1173,8 +1187,9 @@ def grid(self, resolution=None, dim="step"): if resolution is None: if not hasattr(self, '_default_grid'): # TODO THIS is slow now. Decrease resolution and then use find_fluxsurface_step (!!!) - R = np.linspace(self._basedata.R.min().item(), self._basedata.R.max().item(), 1000) - Z = np.linspace(self._basedata.Z.min().item(), self._basedata.Z.max().item(), 2000) + grid_cfg = get_settings().grid + R = np.linspace(self._basedata.R.min().item(), self._basedata.R.max().item(), grid_cfg.default_nr) + Z = np.linspace(self._basedata.Z.min().item(), self._basedata.Z.max().item(), grid_cfg.default_nz) self._default_grid = self.coordinates(R=R, Z=Z, grid=True) return self._default_grid else: @@ -1841,11 +1856,12 @@ def trace_field_line(self, *coordinates, R: np.array = None, Z: np.array = None, else: phi0 = coords.phi[i] - atol = 1e-6 + flt_cfg = get_settings().field_line_tracing + atol = flt_cfg.atol if self.is_xpoint_plasma: xp = self._x_point xp_dist = np.sqrt(np.sum((xp - y0) ** 2)) - atol = np.minimum(xp_dist * 1e-3, atol) + atol = np.minimum(xp_dist * flt_cfg.x_point_atol_scale, atol) if self._verbose: print(f'>>> tracing from: {y0[0]:3f},{y0[1]:3f},{phi0:3f}') @@ -1874,7 +1890,7 @@ def trace_field_line(self, *coordinates, R: np.array = None, Z: np.array = None, stopper = flt.poloidal_angle_stopper_factory(y0, self.magnetic_axis.as_array()[0], dphidtheta * direction, - stop_res=np.pi / 1024) + stop_res=flt_cfg.poloidal_stop_resolution) else: if self._verbose: print('>>> z-lim stopper is used') @@ -1890,18 +1906,18 @@ def trace_field_line(self, *coordinates, R: np.array = None, Z: np.array = None, dphidtheta = np.sign(self.F0) * self._cocosdic['sigma_pol'] * self._cocosdic['sigma_cyl'] stopper = flt.poloidal_angle_stopper_factory(y0, self.magnetic_axis.as_array()[0], dphidtheta * direction, - stop_res=np.pi / 1024) + stop_res=flt_cfg.poloidal_stop_resolution) # todo: define somehow sufficient tolerances sol = solve_ivp(dphifunc, - (phi0, direction * sigma_B0 * (2 * np.pi * 50 + phi0)), + (phi0, direction * sigma_B0 * (2 * np.pi * flt_cfg.max_toroidal_turns + phi0)), y0, # method='RK45', method='LSODA', events=stopper, - max_step=1e-2, # we want high phi resolution + max_step=flt_cfg.max_step, # we want high phi resolution atol=atol, - rtol=1e-8, + rtol=flt_cfg.rtol, ) if self._verbose: @@ -1956,7 +1972,7 @@ def trace_field_line(self, *coordinates, R: np.array = None, Z: np.array = None, return res - def lcfs_field_line(self, vect_no=0, xp_shift=1e-6, phi0: float = 0.0): + def lcfs_field_line(self, vect_no=0, xp_shift=None, phi0: float = 0.0): """ Computes (some) field line laying on last closed flux surface. @@ -1965,7 +1981,7 @@ def lcfs_field_line(self, vect_no=0, xp_shift=1e-6, phi0: float = 0.0): Index of eigenvector determing the direction of integration. Default is 0. xp_shift: float A small positional adjustment for the x-point in the plasma boundary tracking. - Default is 1e-6. + Defaults to the value from PLEQUE settings (`lcfs.x_point_shift`). phi0: float Toroidal angle on which is field line initiated. @@ -1974,10 +1990,12 @@ def lcfs_field_line(self, vect_no=0, xp_shift=1e-6, phi0: float = 0.0): A representation of the LCFS field line as a result of the plasma boundary tracking method. """ - lcfs = track_plasma_boundary(self, self._x_point, vect_no=0, xp_shift=1e-6, phi_0=phi0) + if xp_shift is None: + xp_shift = get_settings().lcfs.x_point_shift + lcfs = track_plasma_boundary(self, self._x_point, vect_no=vect_no, xp_shift=xp_shift, phi_0=phi0) return lcfs - def trace_flux_surface(self, *coordinates, s_resolution=1e-3, R=None, + def trace_flux_surface(self, *coordinates, s_resolution=None, R=None, Z=None, psi_n=None, coord_type=None, **coords): """ Find a closed flux surface inside LCFS with requested values of psi or psi-normalized. @@ -1992,13 +2010,18 @@ def trace_flux_surface(self, *coordinates, s_resolution=1e-3, R=None, :param coordinates: specifies flux surface to search for (by spatial point or values of psi or psi normalised). If coordinates is spatial point (dim=2) then the trace starts at the midplane. Coordinates.grid must be False. - :param s_resolution: max_step in the distance along the flux surface contour + :param s_resolution: max_step in the distance along the flux surface contour. + Defaults to the value from PLEQUE settings (`flux_surfaces.trace_step`). :return: FluxSurface """ from scipy.integrate import solve_ivp import pleque.utils.field_line_tracers as flt + fs_cfg = get_settings().flux_surfaces + if s_resolution is None: + s_resolution = fs_cfg.trace_step + coords = self.coordinates(*coordinates, R=R, Z=Z, coord_type=coord_type, **coords) if coords.dim == 1: coords = self.coordinates(R=self.magnetic_axis.R + coords.r_mid, Z=0) @@ -2010,7 +2033,7 @@ def trace_flux_surface(self, *coordinates, s_resolution=1e-3, R=None, atol=s_resolution ** 2) sol = solve_ivp(ds_func, - (0, coords.r_mid * 2 * np.pi * 4), + (0, coords.r_mid * 2 * np.pi * fs_cfg.trace_max_turns), y0, method='LSODA', events=stopper, @@ -2035,13 +2058,15 @@ def surfacefuncs(self): self._surfacefunc = SurfaceFunctions(self) # filters out methods from self return self._surfacefunc - def to_geqdsk(self, file, nx=64, ny=128, q_positive=True, use_basedata=False, cocos_out=3): + def to_geqdsk(self, file, nx=None, ny=None, q_positive=True, use_basedata=False, cocos_out=3): """ Write a GEQDSK/g-file equilibrium file. :param file: str, file name - :param nx: int, number radial points and profiles points - :param ny: int, number of vertical points + :param nx: int, number radial points and profiles points. + Defaults to the value from PLEQUE settings (`io.geqdsk_nx`). + :param ny: int, number of vertical points. + Defaults to the value from PLEQUE settings (`io.geqdsk_ny`). :param use_basedata: The original basedata of equilibrium are used instead of interpolation splines. If this option is chosen, the nx and ny parameters are ignored. :param q_positive: Save q value always positive. @@ -2085,7 +2110,7 @@ def is_limter_plasma(self): def _map_midplane2psi(self): from scipy.interpolate import UnivariateSpline - r_mid = np.linspace(0, self.R_max - self._mg_axis[0], 100) + r_mid = np.linspace(0, self.R_max - self._mg_axis[0], get_settings().flux_surfaces.midplane_map_points) psi_mid = self.psi(r_mid + self._mg_axis[0], self._mg_axis[1] * np.ones_like(r_mid), grid=False) if self._psi_axis < self._psi_lcfs: @@ -2098,16 +2123,18 @@ def _map_midplane2psi(self): psi_mid = psi_mid[idxs] r_mid = r_mid[idxs] - self._rmid_spl = UnivariateSpline(psi_mid, r_mid, k=3, s=0) + spl_cfg = get_settings().splines + self._rmid_spl = UnivariateSpline(psi_mid, r_mid, k=spl_cfg.profile_order, s=spl_cfg.profile_smooth) def _init_fluxsurfaces(self, npsi: int | None = None, psi_n_levels: Sequence[float] | None = None): if psi_n_levels and npsi: raise ValueError("npsi and psi_n_levels cannot be used simultaneously.") if psi_n_levels is None: - psi0 = settings.psin0 + fs_cfg = get_settings().flux_surfaces + psi0 = fs_cfg.psi_n_min if npsi is None: - npsi = settings.npsi_grid + npsi = fs_cfg.n_psi psi_n_levels = np.linspace(psi0, 1, npsi) else: npsi = len(psi_n_levels) @@ -2121,7 +2148,8 @@ def _init_fluxsurfaces(self, npsi: int | None = None, psi_n_levels: Sequence[flo self._flux_surfaces = surfs def _init_q(self): - psi_n = np.arange(0.01, 1, 0.005) + fs_cfg = get_settings().flux_surfaces + psi_n = np.arange(fs_cfg.q_psi_n_min, 1, fs_cfg.q_psi_n_step) qs = [] if self._verbose: @@ -2134,6 +2162,7 @@ def _init_q(self): qs.append(c.eval_q) qs = np.array(qs) - self._q_spl = UnivariateSpline(psi_n, qs, s=0, k=3) + spl_cfg = get_settings().splines + self._q_spl = UnivariateSpline(psi_n, qs, s=spl_cfg.profile_smooth, k=spl_cfg.profile_order) self._dq_dpsin_spl = self._q_spl.derivative() self._q_anideriv_spl = self._q_spl.antiderivative() diff --git a/pleque/io/geqdsk.py b/pleque/io/geqdsk.py index e5d3bfb..2983a71 100644 --- a/pleque/io/geqdsk.py +++ b/pleque/io/geqdsk.py @@ -3,10 +3,15 @@ import numpy as np import pleque +from pleque.config.settings import get_settings from ._geqdsk import read_as_equilibrium from ._geqdsk import write as write_geqdsk +# Sentinel distinguishing "argument not given" from an explicit None (None has its own +# meaning for `nbdry` in `write`). +_UNSET = object() + def read(file, cocos=3, first_wall=None): """ @@ -164,18 +169,19 @@ def basedata_to_dict(equilibrium: pleque.Equilibrium, cocos_out=13): return data -def write(equilibrium: pleque.Equilibrium, file, nx=64, ny=128, nbdry=200, label=None, cocos_out=3, +def write(equilibrium: pleque.Equilibrium, file, nx=None, ny=None, nbdry=_UNSET, label=None, cocos_out=3, q_positive=True, use_basedata=False): """ Write a GEQDSK equilibrium file. :param equilibrium: pleque.Equilibrum :param file: str, file name - :param nx: int, R-dimension - :param ny: int, Z-dimension + :param nx: int, R-dimension. Defaults to the value from PLEQUE settings (`io.geqdsk_nx`). + :param ny: int, Z-dimension. Defaults to the value from PLEQUE settings (`io.geqdsk_ny`). :param nbdry: int, None - Maximal number of points used to describe boundary (LCFS). If None, default number obtained by FluxFunc will be used. Note: Actual number number of points can be in (n-bry/2, n-bry) + Defaults to the value from PLEQUE settings (`io.geqdsk_nbdry`). :param label: str, max 11 characters long text added on the beginning of g-file (default is PLEQUE) :param cocos_out: At the moment only perform 2pi normalization (!) (TODO) :param q_positive: always save q positive @@ -202,6 +208,14 @@ def write(equilibrium: pleque.Equilibrium, file, nx=64, ny=128, nbdry=200, label psi 2D array (nx,ny) of poloidal flux """ + io_cfg = get_settings().io + if nx is None: + nx = io_cfg.geqdsk_nx + if ny is None: + ny = io_cfg.geqdsk_ny + if nbdry is _UNSET: + nbdry = io_cfg.geqdsk_nbdry + data = dict.fromkeys(['nx', 'ny', 'rdim', 'zdim', 'rcentr', 'bcentr', 'rleft', 'zmid', 'rmagx', 'zmagx', 'simagx', 'sibdry', 'cpasma', 'F', 'pres', 'q', 'psi']) diff --git a/pleque/io/omas.py b/pleque/io/omas.py index 843c9e8..b2b6279 100644 --- a/pleque/io/omas.py +++ b/pleque/io/omas.py @@ -2,6 +2,7 @@ import omas import xarray as xr +from pleque.config.settings import get_settings from pleque.core import Equilibrium @@ -137,11 +138,12 @@ def write(equilibrium: Equilibrium, grid_1d=None, grid_2d=None, gridtype=1, ods= if ods is None: ods = omas.ODS(cocosio=cocosio) + io_cfg = get_settings().io if grid_1d is None: - grid_1d = equilibrium.coordinates(psi_n=np.linspace(0, 1, 200)) + grid_1d = equilibrium.coordinates(psi_n=np.linspace(0, 1, io_cfg.omas_n_psi)) if grid_2d is None: - grid_2d = equilibrium.grid(resolution=(1e-3, 1e-3), dim="step") + grid_2d = equilibrium.grid(resolution=(io_cfg.omas_grid_step, io_cfg.omas_grid_step), dim="step") shot_time = equilibrium.time if equilibrium.time_unit == "ms": diff --git a/pleque/utils/equi_tools.py b/pleque/utils/equi_tools.py index 3e720c0..bfd16bd 100644 --- a/pleque/utils/equi_tools.py +++ b/pleque/utils/equi_tools.py @@ -12,6 +12,7 @@ import xarray as xa import pleque.utils.surfaces as surf +from pleque.config.settings import get_settings from pleque.utils.surfaces import find_contour, points_inside_curve @@ -23,31 +24,39 @@ def psi_grad_sq(x): return psi_grad_sq -def _get_psi_n_on_q(eq, q, max_psi_n=0.99): +def _get_psi_n_on_q(eq, q, max_psi_n=None): + psi_n_cap = get_settings().flux_surfaces.q_search_psi_n_cap + if max_psi_n is None: + max_psi_n = psi_n_cap + # todo: brentq method is probably not the fastest. if not (np.abs(eq.q(0)) < q < np.abs(eq.q(max_psi_n))): return None - psi_n = brentq(lambda psi_n: np.abs(eq.q(psi_n)) - q, 0, 0.99) + psi_n = brentq(lambda psi_n: np.abs(eq.q(psi_n)) - q, 0, psi_n_cap) return psi_n -def get_psi_n_on_q(eq, q, max_psi_n=0.95): +def get_psi_n_on_q(eq, q, max_psi_n=None): + if max_psi_n is None: + max_psi_n = get_settings().flux_surfaces.q_search_psi_n_max if isinstance(q, Iterable): return [_get_psi_n_on_q(eq, _q, max_psi_n=max_psi_n) for _q in q] return _get_psi_n_on_q(eq, q, max_psi_n=max_psi_n) -def is_monotonic(f, x0, x1, n_test=10): +def is_monotonic(f, x0, x1, n_test=None): """ Test whether line connection of two points is monotonic on f. :param f: 2D spline `f(x[0], x[1])` :param x0: start point (2d) of the line :param x1: end point (2d) of the line - :param n_test: number of points which are tested. + :param n_test: number of points which are tested. Defaults to the value from PLEQUE settings. :return: logic value """ + if n_test is None: + n_test = get_settings().critical_points.monotonicity_test_points rpts = np.linspace(x0[0], x1[0], n_test) zpts = np.linspace(x0[1], x1[1], n_test) psi_test = f(rpts, zpts, grid=False) @@ -68,23 +77,25 @@ def minimize_in_vicinity(point, func, r_lims, z_lims): # minimize in the vicinity: # Study different methods and find the most propriate and fastest! - bounds = ((np.max((r_lims[0], point[0] - 0.1)), - np.min((r_lims[-1], point[0] + 0.1))), - (np.max((z_lims[0], point[1] - 0.1)), - np.min((z_lims[-1], point[1] + 0.1)))) - - res = minimize(func, point, method='Powell', options={'xtol': 1e-7}) + cp_cfg = get_settings().critical_points + vicinity = cp_cfg.vicinity_radius + bounds = ((np.max((r_lims[0], point[0] - vicinity)), + np.min((r_lims[-1], point[0] + vicinity))), + (np.max((z_lims[0], point[1] - vicinity)), + np.min((z_lims[-1], point[1] + vicinity)))) + + res = minimize(func, point, method='Powell', options={'xtol': cp_cfg.minimizer_xtol}) res_point = np.array((res['x'][0], res['x'][1])) # If unbounded Powell algorithm finds wrong minimum, algorithm with bounds is used. - if np.sum(res_point ** 2 - point ** 2) > 1e-2: - res = minimize(func, point, method='TNC', bounds=bounds, options={'xtol': 1e-7}) + if np.sum(res_point ** 2 - point ** 2) > cp_cfg.relocation_threshold: + res = minimize(func, point, method='TNC', bounds=bounds, options={'xtol': cp_cfg.minimizer_xtol}) res_point = np.array((res['x'][0], res['x'][1])) return res_point -def find_extremes(rs, zs, psi_spl, order=20): +def find_extremes(rs, zs, psi_spl, order=None): """ Find the extremes on grid given by rs and zs. x-points: Candidates for x-point @@ -95,11 +106,15 @@ def find_extremes(rs, zs, psi_spl, order=20): :param order: int, order used by scipy argrelmin function (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.argrelmin.html) How many points on each side to use for the comparison to - consider comparator(n, n+x) to be True + consider comparator(n, n+x) to be True. + Defaults to the value from PLEQUE settings. :param psi_spl: :return: tuple(x-points, o-points) of arrays(N, 2) """ + cp_cfg = get_settings().critical_points + if order is None: + order = cp_cfg.find_extremes_order psi_xysq_func = _make_psi_grad_sq(psi_spl) @@ -116,7 +131,7 @@ def psi_2nd_derivatives(r_coord, z_coord): x_points = [] iteration = 0 - while len(o_points) == 0 and order > 0 and iteration < 10: + while len(o_points) == 0 and order > 0 and iteration < cp_cfg.find_extremes_max_iter: iteration += 1 o_points = [] @@ -142,7 +157,7 @@ def psi_2nd_derivatives(r_coord, z_coord): z_ex = zs[az] # XXX Remove bad candidates for the extreme (this is potentional trouble point): - if psi_xysq_func((r_ex, z_ex)) > 1: # 1e3 * dpsidx: + if psi_xysq_func((r_ex, z_ex)) > cp_cfg.gradient_threshold: # 1e3 * dpsidx: continue psi_xx, psi_yy, psi_xy2 = psi_2nd_derivatives(r_ex, z_ex) @@ -196,11 +211,13 @@ def recognize_mg_axis(o_points, psi_spl, r_lims, z_lims, first_wall=None, mg_axi :return: tuple with recognize magnetic axis point and arguments of sorted o-points """ + cp_cfg = get_settings().critical_points + if mg_axis_candidate is None: r_centr = (r_lims[0] + r_lims[-1]) / 2 z_centr = (z_lims[0] + z_lims[-1]) / 2 # vertical distance if favoured - op_dist = 5 * (o_points[:, 0] - r_centr) ** 2 + (o_points[:, 1] - z_centr) ** 2 + op_dist = cp_cfg.axis_vertical_weight * (o_points[:, 0] - r_centr) ** 2 + (o_points[:, 1] - z_centr) ** 2 else: op_dist = (o_points[:, 0] - mg_axis_candidate[0]) ** 2 + (o_points[:, 1] - mg_axis_candidate[1]) ** 2 # normalise the maximal distance to one @@ -217,7 +234,7 @@ def recognize_mg_axis(o_points, psi_spl, r_lims, z_lims, first_wall=None, mg_axi if first_wall is not None and len(first_wall) > 2: mask_in = points_inside_curve(o_points, first_wall) op_in_first_wall[mask_in] = 1 - op_in_first_wall[np.logical_not(mask_in)] = 1e-3 + op_in_first_wall[np.logical_not(mask_in)] = cp_cfg.axis_out_of_wall_penalty sortidx = np.argsort(op_dist * op_psiscale * (1 - op_in_first_wall)) @@ -254,12 +271,15 @@ def recognize_x_points(x_points, mg_axis, psi_axis, psi_spl, r_lims, z_lims, psi len_diff = (x_point_candidates[0] - x_points[:, 0]) ** 2 + (x_point_candidates[1] - x_points[:, 1]) ** 2 len_diff = len_diff / np.max(len_diff) + cp_cfg = get_settings().critical_points + sort_eps = cp_cfg.x_point_sort_epsilon + for i, xpoint in enumerate(x_points): - monotonic[i] = is_monotonic(psi_spl, mg_axis, xpoint, 10) - monotonic[i] = (1 - monotonic[i] * 1) + 1e-3 + monotonic[i] = is_monotonic(psi_spl, mg_axis, xpoint, cp_cfg.monotonicity_test_points) + monotonic[i] = (1 - monotonic[i] * 1) + sort_eps - # The monotonic points are preferred (addition of 1e-3 is to avoid zero difference) - sortidx = np.argsort((psi_diff + 1e-3) * monotonic * len_diff) + # The monotonic points are preferred (addition of the epsilon is to avoid zero difference) + sortidx = np.argsort((psi_diff + sort_eps) * monotonic * len_diff) xp1 = x_points[sortidx[0]] if len(x_points) > 1: @@ -302,7 +322,8 @@ def recognize_plasma_type(x_point, first_wall, mg_axis, psi_axis, psi_spl): i = 0 - while not (i == len(idxs_wall) or is_monotonic(psi_spl, first_wall[idxs_wall[i]], mg_axis, 50)): + n_test = get_settings().critical_points.limiter_monotonicity_test_points + while not (i == len(idxs_wall) or is_monotonic(psi_spl, first_wall[idxs_wall[i]], mg_axis, n_test)): i += 1 if i == len(idxs_wall): @@ -337,7 +358,7 @@ def find_close_lcfs(psi_lcfs, rs, zs, psi_spl, mg_axis, psi_axis=0): :return: """ - new_psi_lcfs = psi_lcfs - 1e-4 * (psi_lcfs - psi_axis) + new_psi_lcfs = psi_lcfs - get_settings().lcfs.initial_psi_offset * (psi_lcfs - psi_axis) contours = find_contour(psi_spl(rs, zs, grid=True).T, new_psi_lcfs, rs, zs) @@ -395,8 +416,9 @@ def find_surface_step(psi_spl, psi_target, flux_surf): psix = psix / (deriv_norm ** 2) psiy = psiy / (deriv_norm ** 2) - flux_surf[:, 0] -= 0.99 * psix * (psi - psi_target) - flux_surf[:, 1] -= 0.99 * psiy * (psi - psi_target) + damping = get_settings().lcfs.surface_step_damping + flux_surf[:, 0] -= damping * psix * (psi - psi_target) + flux_surf[:, 1] -= damping * psiy * (psi - psi_target) return flux_surf diff --git a/pleque/utils/field_line_tracers.py b/pleque/utils/field_line_tracers.py index 371bc0b..0e82cc8 100644 --- a/pleque/utils/field_line_tracers.py +++ b/pleque/utils/field_line_tracers.py @@ -2,6 +2,8 @@ import numpy as np +from pleque.config.settings import get_settings + def dphi_tracer_factory(BR_func, BZ_func, Bphi_func, BR_pert_func=None, BZ_pert_func=None, direction=1): """Factory for function $d[R,Z]/d\\phi=f(\\phi, [R,Z])$ @@ -85,7 +87,7 @@ def ds_func(s, x): return ds_func -def poloidal_angle_stopper_factory(y0, y_center, direction, stop_res=np.pi / 360): +def poloidal_angle_stopper_factory(y0, y_center, direction, stop_res=None): """Factory for function which stops field line tracing close to the original poloidal angle Suitable for the *events* argument of :func:`scipy.integrate.solve_ivp` @@ -99,8 +101,11 @@ def poloidal_angle_stopper_factory(y0, y_center, direction, stop_res=np.pi / 360 sign of dtheta/dphi derivative, depends on plasma current and toroidal mag. field stop_res : float stopping offset to stop before initial position - necessary to prevent stopping at initial position + necessary to prevent stopping at initial position. + Defaults to the value from PLEQUE settings (`field_line_tracing.default_stop_resolution`). """ + if stop_res is None: + stop_res = get_settings().field_line_tracing.default_stop_resolution y_center = np.asarray(y_center) # should be [R0, Z0] def full_arc(y): @@ -186,10 +191,11 @@ def ds_func(s, x): return ds_func -def rz_target_s_min_stopper_factory(y0, s_min, atol=1e-6): +def rz_target_s_min_stopper_factory(y0, s_min, atol=None): """Stopper which returns the (squared) distance after a minimum traced length atol is subtracted from the distance to guarantee 0-crossing + (defaults to the value from PLEQUE settings, `field_line_tracing.target_stopper_atol`) This should prevent early stopping by having the initial distance 0 @@ -197,6 +203,9 @@ def rz_target_s_min_stopper_factory(y0, s_min, atol=1e-6): A good estimate for s_min could be r_mid of y0=[Rt, Zt], because the total flux surface length will be at least pi times larger """ + if atol is None: + atol = get_settings().field_line_tracing.target_stopper_atol + def stopper(s, y): if s < s_min: # only started tracing return 42 diff --git a/pleque/utils/plotting.py b/pleque/utils/plotting.py index 50d13a7..1fd7a0c 100644 --- a/pleque/utils/plotting.py +++ b/pleque/utils/plotting.py @@ -2,6 +2,7 @@ import numpy as np import pleque +from pleque.config.settings import get_settings from pleque.utils.equi_tools import get_psi_n_on_q @@ -26,12 +27,13 @@ def _plot_debug(eq: pleque.Equilibrium, ax: plt.Axes = None, levels=None, colorb if ax is None: ax = plt.gca() - rs = np.linspace(eq.R_min, eq.R_max, 400) - zs = np.linspace(eq.Z_min, eq.Z_max, 600) + plot_cfg = get_settings().plotting + rs = np.linspace(eq.R_min, eq.R_max, plot_cfg.grid_nr) + zs = np.linspace(eq.Z_min, eq.Z_max, plot_cfg.grid_nz) try: if levels is None: - levels = 60 + levels = plot_cfg.debug_contour_levels cl = ax.contour(rs, zs, eq._spl_psi(rs, zs).T, levels) if colorbar: plt.contour(cl) @@ -101,14 +103,15 @@ def plot_rational_surface(eq, ax, q_tuple, linestyles="--", colors="C3"): qarr = np.array(q_tuple) q = qarr[:, 0] / qarr[:, 1] - psi_n = get_psi_n_on_q(eq, q, max_psi_n=0.95) + plot_cfg = get_settings().plotting + psi_n = get_psi_n_on_q(eq, q, max_psi_n=plot_cfg.rational_q_psi_n_max) i_ok = np.nonzero(psi_n)[0] psi_n = np.atleast_1d(psi_n)[i_ok] qarr = qarr[i_ok, :] - coords = eq.grid((400, 600), 'size') + coords = eq.grid((plot_cfg.grid_nr, plot_cfg.grid_nz), 'size') mask_inlcfs = eq.in_lcfs(coords) psi_ns = np.ma.masked_array(coords.psi_n, np.logical_not(mask_inlcfs)) @@ -145,7 +148,8 @@ def plot_lcfs(eq, ax, color="C1", lw=2, ls="--"): def plot_psi_contours(eq, ax, where="in_lcfs", alpha=1, **kwargs): - coords = eq.grid((400, 600), 'size') + plot_cfg = get_settings().plotting + coords = eq.grid((plot_cfg.grid_nr, plot_cfg.grid_nz), 'size') psi = eq.psi(coords) mask_inlcfs = eq.in_lcfs(coords) @@ -155,21 +159,25 @@ def plot_psi_contours(eq, ax, where="in_lcfs", alpha=1, **kwargs): elif where.lower() == "out_lcfs": psi = np.ma.masked_array(psi, mask_inlcfs) - cl = ax.contour(coords.R, coords.Z, psi, 20, alpha=alpha, **kwargs) + cl = ax.contour(coords.R, coords.Z, psi, plot_cfg.psi_contour_levels, alpha=alpha, **kwargs) return cl -def plot_near_sol(eq: pleque.Equilibrium, ax: plt.Axes, colors="C0", dr:float = 2e-3, lw=0.7, ls="solid"): +def plot_near_sol(eq: pleque.Equilibrium, ax: plt.Axes, colors="C0", dr: float | None = None, lw=0.7, ls="solid"): + plot_cfg = get_settings().plotting + if dr is None: + dr = plot_cfg.sol_spacing contour_out = eq.coordinates(r=eq.lcfs.r_mid[0] + dr * np.arange(1, 6), theta=np.zeros(5), grid=False) - coords = eq.grid((400, 600), 'size') + coords = eq.grid((plot_cfg.grid_nr, plot_cfg.grid_nz), 'size') ax.contour(coords.R, coords.Z, coords.psi, np.sort(np.squeeze(contour_out.psi)), colors=colors, linewidths=lw, linestyles=ls) -def plot_equilibrium(eq: pleque.Equilibrium, ax: plt.Axes = None, colorbar=False, dr_sol:float = 2e-3, **kwargs): +def plot_equilibrium(eq: pleque.Equilibrium, ax: plt.Axes = None, colorbar=False, dr_sol: float | None = None, + **kwargs): if ax is None: ax = plt.gca() @@ -206,13 +214,15 @@ def normalize_axis_xylim_by_first_wall(eq, ax): rlim = [np.min(eq.first_wall.R), np.max(eq.first_wall.R)] zlim = [np.min(eq.first_wall.Z), np.max(eq.first_wall.Z)] + margin = get_settings().plotting.axis_margin + size = rlim[1] - rlim[0] - rlim[0] -= size / 12 - rlim[1] += size / 12 + rlim[0] -= size * margin + rlim[1] += size * margin size = zlim[1] - zlim[0] - zlim[0] -= size / 12 - zlim[1] += size / 12 + zlim[0] -= size * margin + zlim[1] += size * margin ax.set_xlim(*rlim) ax.set_ylim(*zlim) diff --git a/pleque/utils/surfaces.py b/pleque/utils/surfaces.py index ee55c6c..9a07710 100644 --- a/pleque/utils/surfaces.py +++ b/pleque/utils/surfaces.py @@ -2,6 +2,8 @@ import shapely.geometry as geo from skimage import measure +from pleque.config.settings import get_settings + def find_contour(array, level, r=None, z=None, fully_connected="low", positive_orientation="low"): """ @@ -108,20 +110,25 @@ def points_inside_curve(points, contour): return measure.points_in_poly(points, contour) -def get_surface(equilibrium, psi, r=100, z=100, norm=True, closed=True, insidelcfs=True): +def get_surface(equilibrium, psi, r=None, z=None, norm=True, closed=True, insidelcfs=True): """ Finds points of surface with given value of psi. :param equilibrium: Equilibrium object :param psi: Value of psi to get the surface for :param r: If number, specifies number of points in the r dimension of the mesh. If numpy array, - gives r grid points. + gives r grid points. Defaults to the value from PLEQUE settings (`flux_surfaces.contour_grid_nr`). :param z: If number, specifies number of points in the z dimension of the mesh. If numpy array, - gives z grid points. + gives z grid points. Defaults to the value from PLEQUE settings (`flux_surfaces.contour_grid_nz`). :param norm: Specifies whether we are working with normalised values of psi :param closed: Are we looking for a closed surface? :param insidelcfs: Are we looking for a closed surface inside lcfs? :return: List of contours with surface coordinates """ + fs_cfg = get_settings().flux_surfaces + if r is None: + r = fs_cfg.contour_grid_nr + if z is None: + z = fs_cfg.contour_grid_nz # if r is integer make r grid if not isinstance(r, np.ndarray): @@ -166,7 +173,7 @@ def curve_is_closed(points): return np.isclose(points[0, 0], points[-1, 0]) and np.isclose(points[0, 1], points[-1, 1]) -def point_inside_fluxsurface(equilibrium, points, psi, r=100, z=100, norm=True, +def point_inside_fluxsurface(equilibrium, points, psi, r=None, z=None, norm=True, insidelcfs=True): """ Checks if a point is inside a flux surface with specified value of psi. @@ -175,9 +182,9 @@ def point_inside_fluxsurface(equilibrium, points, psi, r=100, z=100, norm=True, :param points: 2d numpy array (N, 2) of points coordinates :param psi: value of the psi on the surface :param r: If number, specifies number of points in the r dimension of the mesh. If numpy array, - gives r grid points. + gives r grid points. Defaults to the value from PLEQUE settings (`flux_surfaces.contour_grid_nr`). :param z: If number, specifies number of points in the z dimension of the mesh. If numpy array, - gives z grid points. + gives z grid points. Defaults to the value from PLEQUE settings (`flux_surfaces.contour_grid_nz`). :param norm: Specifies whether we are working with normalised values of psi :param insidelcfs: Are we looking for a closed surface inside lcfs? :return: array of bool @@ -206,7 +213,7 @@ def point_in_first_wall(equilibrium, points): return isinside -def track_plasma_boundary(equilibrium, xp, xp_shift=1e-6, vect_no=0, phi_0=0): +def track_plasma_boundary(equilibrium, xp, xp_shift=None, vect_no=0, phi_0=0): """ Tracing one of two separa trix branches (switched by `vect_no`) which are goes around magnetic axis for `direction = 1` or in the opposite direction for `direction = -1`. @@ -214,6 +221,8 @@ def track_plasma_boundary(equilibrium, xp, xp_shift=1e-6, vect_no=0, phi_0=0): :param equilibrium: :type equilibrium: pleque.Equilibrium :param xp: x-point position + :param xp_shift: Distance by which the tracing start point is shifted from the x-point. + Defaults to the value from PLEQUE settings (`lcfs.x_point_shift`). :param vect_no: (0, 1) Choose one of the eigen vectors of matrix of field line differential equation. :return: """ @@ -221,6 +230,9 @@ def track_plasma_boundary(equilibrium, xp, xp_shift=1e-6, vect_no=0, phi_0=0): from pleque.utils.tools import xp_vecs + if xp_shift is None: + xp_shift = get_settings().lcfs.x_point_shift + evecs, _ = xp_vecs(equilibrium._spl_psi, *xp) mg_axis = equilibrium._mg_axis diff --git a/poetry.lock b/poetry.lock index a3bbc6e..482aff0 100644 --- a/poetry.lock +++ b/poetry.lock @@ -2821,7 +2821,7 @@ description = "Provides an object-oriented python interface to the netCDF versio optional = false python-versions = ">=3.10" groups = ["main"] -markers = "python_version >= \"3.11\" or platform_system != \"Windows\" or platform_machine != \"ARM64\"" +markers = "platform_system != \"Windows\" or platform_machine != \"ARM64\" or python_version >= \"3.11\"" files = [ {file = "netcdf4-1.7.4-cp310-cp310-macosx_13_0_x86_64.whl", hash = "sha256:b1c1a7ea3678db76bf33d14f7e202385d634db38c5e70d8cf4895971023eebb9"}, {file = "netcdf4-1.7.4-cp310-cp310-macosx_14_0_arm64.whl", hash = "sha256:d3f9497873454207f9480847d02b1b19a4bc81ad6e9166e1c17d4e2f8f3555d1"}, @@ -2851,8 +2851,8 @@ files = [ certifi = "*" cftime = "*" numpy = [ - {version = ">=2.3.0", markers = "platform_system == \"Windows\" and platform_machine == \"ARM64\""}, {version = ">=1.21.2", markers = "platform_system != \"Windows\" or platform_machine != \"ARM64\""}, + {version = ">=2.3.0", markers = "platform_system == \"Windows\" and platform_machine == \"ARM64\""}, ] [package.extras] @@ -2911,7 +2911,7 @@ description = "Python package for creating and manipulating graphs and networks" optional = false python-versions = "!=3.14.1,>=3.11" groups = ["main"] -markers = "python_version < \"3.14\" and python_version >= \"3.11\"" +markers = "python_version >= \"3.11\" and python_version < \"3.14\"" files = [ {file = "networkx-3.6.1-py3-none-any.whl", hash = "sha256:d47fbf302e7d9cbbb9e2555a0d267983d2aa476bac30e90dfbe5669bd57f3762"}, {file = "networkx-3.6.1.tar.gz", hash = "sha256:26b7c357accc0c8cde558ad486283728b65b6a95d85ee1cd66bafab4c8168509"}, @@ -3897,6 +3897,7 @@ files = [ [package.dependencies] pydantic = ">=2.7.0" python-dotenv = ">=0.21.0" +tomli = {version = ">=2.0.1", optional = true, markers = "extra == \"toml\""} typing-inspection = ">=0.4.0" [package.extras] @@ -5167,7 +5168,7 @@ description = "Python documentation generator" optional = false python-versions = ">=3.11" groups = ["compass", "dev"] -markers = "python_version == \"3.11\"" +markers = "python_version >= \"3.11\" and python_version < \"3.14\"" files = [ {file = "sphinx-9.0.4-py3-none-any.whl", hash = "sha256:5bebc595a5e943ea248b99c13814c1c5e10b3ece718976824ffa7959ff95fffb"}, {file = "sphinx-9.0.4.tar.gz", hash = "sha256:594ef59d042972abbc581d8baa577404abe4e6c3b04ef61bd7fc2acbd51f3fa3"}, @@ -5199,7 +5200,7 @@ description = "Python documentation generator" optional = false python-versions = ">=3.12" groups = ["compass", "dev"] -markers = "python_version >= \"3.12\"" +markers = "python_version >= \"3.14\"" files = [ {file = "sphinx-9.1.0-py3-none-any.whl", hash = "sha256:c84fdd4e782504495fe4f2c0b3413d6c2bf388589bb352d439b2a3bb99991978"}, {file = "sphinx-9.1.0.tar.gz", hash = "sha256:7741722357dd75f8190766926071fed3bdc211c74dd2d7d4df5404da95930ddb"}, @@ -5377,7 +5378,7 @@ description = "Read and write TIFF files" optional = false python-versions = ">=3.11" groups = ["main"] -markers = "python_version == \"3.11\"" +markers = "python_version >= \"3.11\" and python_version < \"3.14\"" files = [ {file = "tifffile-2026.3.3-py3-none-any.whl", hash = "sha256:e8be15c94273113d31ecb7aa3a39822189dd11c4967e3cc88c178f1ad2fd1170"}, {file = "tifffile-2026.3.3.tar.gz", hash = "sha256:d9a1266bed6f2ee1dd0abde2018a38b4f8b2935cb843df381d70ac4eac5458b7"}, @@ -5401,7 +5402,7 @@ description = "Read and write TIFF files" optional = false python-versions = ">=3.12" groups = ["main"] -markers = "python_version >= \"3.12\"" +markers = "python_version >= \"3.14\"" files = [ {file = "tifffile-2026.6.1-py3-none-any.whl", hash = "sha256:0d7382d2769b855b81ce358528e2b40c16d48aa39031746efa81215205332a8d"}, {file = "tifffile-2026.6.1.tar.gz", hash = "sha256:626c892c0e899d959b9438e7c0e1491dc154a7fead1f1f37a991724a50eceba9"}, @@ -5443,8 +5444,7 @@ version = "2.4.1" description = "A lil' TOML parser" optional = false python-versions = ">=3.8" -groups = ["compass", "dev"] -markers = "python_version == \"3.10\"" +groups = ["main", "compass", "dev"] files = [ {file = "tomli-2.4.1-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:f8f0fc26ec2cc2b965b7a3b87cd19c5c6b8c5e5f436b984e85f486d652285c30"}, {file = "tomli-2.4.1-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:4ab97e64ccda8756376892c53a72bd1f964e519c77236368527f758fbc36a53a"}, @@ -5494,6 +5494,7 @@ files = [ {file = "tomli-2.4.1-py3-none-any.whl", hash = "sha256:0d85819802132122da43cb86656f8d1f8c6587d54ae7dcaf30e90533028b49fe"}, {file = "tomli-2.4.1.tar.gz", hash = "sha256:7c7e1a961a0b2f2472c1ac5b69affa0ae1132c39adcb67aba98568702b9cc23f"}, ] +markers = {compass = "python_version == \"3.10\"", dev = "python_version == \"3.10\""} [[package]] name = "tomlkit" @@ -5827,4 +5828,4 @@ viz = ["matplotlib", "nc-time-axis", "seaborn"] [metadata] lock-version = "2.1" python-versions = "^3.10" -content-hash = "2de640dd7c270dd239793235cf8386487ec1f52e0dc89781fbec7cb7c3cd6315" +content-hash = "d6ff6c345afaa528b1b564eae7c5e2a1f544ab2aeab6399fbf0901ce4ed592e1" diff --git a/pyproject.toml b/pyproject.toml index 5ac705a..86e1796 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,7 @@ netcdf4 = "^1.6.3" h5py = "^3.8.0" scikit-image = ">0.20.0" matplotlib = { version = "^3.7.1", optional = false } -pydantic-settings = "^2.10.1" +pydantic-settings = { version = "^2.10.1", extras = ["toml"] } pydantic = "^2.13.4" [tool.poetry.group.dev.dependencies] diff --git a/tests/test_settings.py b/tests/test_settings.py new file mode 100644 index 0000000..9da5d33 --- /dev/null +++ b/tests/test_settings.py @@ -0,0 +1,196 @@ +"""Tests of the PLEQUE settings module: defaults, file loading precedence, source tracking.""" + +import os + +import pytest +from pydantic import ValidationError + +from pleque.config import settings as settings_module +from pleque.config.settings import ( + CONFIG_FILE_ENV_VAR, + CONFIG_FILE_NAME, + Settings, + get_settings, + load_settings, + reload_settings, + user_config_path, +) + + +@pytest.fixture(autouse=True) +def _isolated_settings(tmp_path, monkeypatch): + """Run each test in an empty cwd with no PLEQUE env vars and a cold settings cache.""" + for key in list(os.environ): + if key.startswith("PLEQUE_"): + monkeypatch.delenv(key) + # point the user config dir somewhere empty so a real ~/.config/pleque does not leak in + monkeypatch.setenv("XDG_CONFIG_HOME", str(tmp_path / "xdg")) + monkeypatch.chdir(tmp_path) + get_settings.cache_clear() + yield + get_settings.cache_clear() + + +def test_defaults_without_any_config(tmp_path): + settings = get_settings() + assert settings.flux_surfaces.n_psi == 200 + assert settings.flux_surfaces.psi_n_min == 0.01 + assert settings.lcfs.search_grid_nr == 700 + assert settings.lcfs.search_grid_nz == 1200 + assert settings.default_cocos == 3 + assert settings.config_source is None + assert tmp_path / CONFIG_FILE_NAME in settings.config_files_consulted + + +def test_load_from_cwd_toml(tmp_path): + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 321\n") + settings = get_settings() + assert settings.flux_surfaces.n_psi == 321 + # untouched values keep their defaults + assert settings.lcfs.search_grid_nr == 700 + assert settings.config_source == tmp_path / CONFIG_FILE_NAME + + +def test_cwd_beats_user_config_dir(tmp_path): + user_file = user_config_path() + user_file.parent.mkdir(parents=True) + user_file.write_text("[flux_surfaces]\nn_psi = 111\n\n[lcfs]\nsearch_grid_nr = 999\n") + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 222\n") + + settings = get_settings() + assert settings.flux_surfaces.n_psi == 222 + assert settings.config_source == tmp_path / CONFIG_FILE_NAME + # files are not merged: the loser file is ignored completely + assert settings.lcfs.search_grid_nr == 700 + + +def test_user_config_dir_posix(tmp_path): + user_file = user_config_path() + assert user_file == tmp_path / "xdg" / "pleque" / CONFIG_FILE_NAME + user_file.parent.mkdir(parents=True) + user_file.write_text("[grid]\ndefault_nr = 555\n") + + settings = get_settings() + assert settings.grid.default_nr == 555 + assert settings.config_source == user_file + + +def test_user_config_dir_windows(tmp_path, monkeypatch): + appdata = tmp_path / "AppData" / "Roaming" + monkeypatch.setenv("APPDATA", str(appdata)) + assert settings_module._windows_config_dir() == appdata / "pleque" + + # without APPDATA the home directory is used as a fallback + monkeypatch.delenv("APPDATA") + monkeypatch.setenv("HOME", str(tmp_path / "home")) + assert settings_module._windows_config_dir() == tmp_path / "home" / "pleque" + + +def test_pyproject_tool_pleque(tmp_path, monkeypatch): + (tmp_path / "pyproject.toml").write_text("[tool.pleque.grid]\ndefault_nr = 555\n") + subdir = tmp_path / "some" / "subdir" + subdir.mkdir(parents=True) + monkeypatch.chdir(subdir) + + settings = get_settings() + assert settings.grid.default_nr == 555 + assert settings.config_source == tmp_path / "pyproject.toml" + + +def test_pyproject_without_table_ignored(tmp_path): + (tmp_path / "pyproject.toml").write_text("[tool.other]\nfoo = 1\n") + settings = get_settings() + assert settings.grid.default_nr == 1000 + assert settings.config_source is None + + +def test_explicit_config_file_env_var(tmp_path, monkeypatch): + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 111\n") + explicit = tmp_path / "elsewhere" / "my_config.toml" + explicit.parent.mkdir() + explicit.write_text("[flux_surfaces]\nn_psi = 333\n") + monkeypatch.setenv(CONFIG_FILE_ENV_VAR, str(explicit)) + + settings = get_settings() + assert settings.flux_surfaces.n_psi == 333 + assert settings.config_source == explicit + assert settings.config_files_consulted == (explicit,) + + +def test_explicit_config_file_env_var_missing_raises(tmp_path, monkeypatch): + monkeypatch.setenv(CONFIG_FILE_ENV_VAR, str(tmp_path / "does_not_exist.toml")) + with pytest.raises(FileNotFoundError): + load_settings() + + +def test_env_overrides_file(tmp_path, monkeypatch): + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 111\n") + monkeypatch.setenv("PLEQUE_FLUX_SURFACES__N_PSI", "999") + + settings = get_settings() + assert settings.flux_surfaces.n_psi == 999 + # the file is still recorded as the source of the remaining values + assert settings.config_source == tmp_path / CONFIG_FILE_NAME + + +def test_nested_env_var_types(monkeypatch): + monkeypatch.setenv("PLEQUE_LCFS__REFINEMENT_TOLERANCE", "1e-8") + settings = get_settings() + assert settings.lcfs.refinement_tolerance == pytest.approx(1e-8) + + monkeypatch.setenv("PLEQUE_LCFS__REFINEMENT_TOLERANCE", "not-a-number") + with pytest.raises(ValidationError): + load_settings() + + +def test_init_kwargs_override_everything(tmp_path): + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 111\n") + settings = Settings(flux_surfaces={"n_psi": 42}) + assert settings.flux_surfaces.n_psi == 42 + + +def test_reload_settings(tmp_path): + assert get_settings().flux_surfaces.n_psi == 200 + (tmp_path / CONFIG_FILE_NAME).write_text("[flux_surfaces]\nn_psi = 321\n") + + settings = reload_settings() + assert settings.flux_surfaces.n_psi == 321 + assert get_settings() is settings + + +def test_get_settings_is_cached(): + assert get_settings() is get_settings() + assert load_settings() is not get_settings() + + +def test_no_import_time_settings_capture(): + """Call sites must read settings at use time, not capture them at import time.""" + import pleque.core.equilibrium as equilibrium_module + + assert not hasattr(equilibrium_module, "settings") + + +def test_settings_used_at_call_time(tmp_path): + """Functions resolving defaults from settings must see freshly reloaded values.""" + from pleque.utils import surfaces + + (tmp_path / CONFIG_FILE_NAME).write_text("[lcfs]\nx_point_shift = 0.125\n") + reload_settings() + assert settings_module.get_settings().lcfs.x_point_shift == 0.125 + # get_settings imported inside pleque.utils.surfaces is the same cached object + assert surfaces.get_settings().lcfs.x_point_shift == 0.125 + + +def test_docs_settings_reference_complete(): + """Every settings section and field must be documented in docs/source/configuration.rst.""" + from pathlib import Path + + docs = Path(__file__).parent.parent / "docs" / "source" / "configuration.rst" + text = docs.read_text() + + for section_name, field in Settings.model_fields.items(): + assert section_name in text, f"settings entry '{section_name}' missing in configuration.rst" + annotation = field.annotation + if hasattr(annotation, "model_fields"): + for field_name in annotation.model_fields: + assert field_name in text, f"setting '{section_name}.{field_name}' missing in configuration.rst"