diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS new file mode 100644 index 00000000..4a9b590c --- /dev/null +++ b/.github/CODEOWNERS @@ -0,0 +1,8 @@ +# CODEOWNERS: Controls automatic assignment of reviewers for pull requests. +# ----------------------------------------------------------------------------- +# When a pull request modifies files matching a pattern, the listed code owners will be requested as reviewers. +# Each line is a file pattern followed by one or more code owners (users or teams) who must review changes to that file. +# Do not edit unless specifically mandated to do so. + +# Global/fallback and technical modifications. +* @SimonPinches @olivhoenen @prasad-sawantdesai diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 00000000..001b22a5 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,15 @@ +version: 2 +updates: + # Enable version updates for Github Actions + - package-ecosystem: "github-actions" + # Look for `/.github/workflows` and `/action.yml` or `.yaml` + directory: "/" + # Check for updates once a week + schedule: + interval: "weekly" + # Group actions version bumps into a single PR + groups: + actions-deps: + patterns: + - "*" + diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml index 392b9602..beb1fa74 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -15,10 +15,10 @@ jobs: steps: - name: Checkout IDStools sources - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v7 with: # until saxonche is available in 3.13 # https://saxonica.plan.io/issues/6561 diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index ac978740..4198a096 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -10,11 +10,11 @@ jobs: name: Build distribution runs-on: ubuntu-22.04 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 with: fetch-depth: 0 - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v7 with: # until saxonche is available in 3.13 # https://saxonica.plan.io/issues/6561 @@ -25,7 +25,7 @@ jobs: - name: Build a binary wheel and a source tarball run: python3 -m build . - name: Store the distribution packages - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: python-package-distributions path: dist/ @@ -43,7 +43,7 @@ jobs: id-token: write # IMPORTANT: mandatory for trusted publishing steps: - name: Download all the dists - uses: actions/download-artifact@v4 + uses: actions/download-artifact@v8 with: name: python-package-distributions path: dist/ @@ -63,7 +63,7 @@ jobs: id-token: write # IMPORTANT: mandatory for trusted publishing steps: - name: Download all the dists - uses: actions/download-artifact@v4 + uses: actions/download-artifact@v8 with: name: python-package-distributions path: dist/ diff --git a/.github/workflows/security.yml b/.github/workflows/security.yml index 63f082c4..ddc45e04 100644 --- a/.github/workflows/security.yml +++ b/.github/workflows/security.yml @@ -10,10 +10,10 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v7 with: python-version: '3.11' cache: 'pip' diff --git a/.github/workflows/test_with_pytest.yml b/.github/workflows/test_with_pytest.yml index ad96ceeb..048cded3 100644 --- a/.github/workflows/test_with_pytest.yml +++ b/.github/workflows/test_with_pytest.yml @@ -17,11 +17,11 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v7 with: python-version: ${{ matrix.python-version }} - name: Display Python version @@ -40,7 +40,7 @@ jobs: TEST_PROFILE: ci - name: Upload coverage report ${{ matrix.python-version }} - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: coverage-report-${{ matrix.python-version }} path: htmlcov @@ -48,7 +48,7 @@ jobs: - name: Upload coverage XML ${{ matrix.python-version }} if: always() - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: coverage-xml-${{ matrix.python-version }} path: coverage.xml @@ -56,7 +56,7 @@ jobs: - name: Upload test report ${{ matrix.python-version }} if: always() - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: test-report-${{ matrix.python-version }} path: junit.xml @@ -64,7 +64,7 @@ jobs: - name: Upload test output logs ${{ matrix.python-version }} if: always() - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: test-output-${{ matrix.python-version }} path: test_output.log diff --git a/.github/workflows/verify_with_sphinx.yml b/.github/workflows/verify_with_sphinx.yml index dc89a842..f5a79d4c 100644 --- a/.github/workflows/verify_with_sphinx.yml +++ b/.github/workflows/verify_with_sphinx.yml @@ -14,10 +14,10 @@ jobs: steps: - name: Checkout IDStools sources - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v7 with: # until saxonche is available in 3.13 # https://saxonica.plan.io/issues/6561 @@ -50,7 +50,7 @@ jobs: - name: Upload Sphinx documentation as artifact if: always() - uses: actions/upload-artifact@v4 + uses: actions/upload-artifact@v7 with: name: sphinx-documentation path: docs/build/html/ diff --git a/README.md b/README.md index ed99be6d..ccf62dca 100644 --- a/README.md +++ b/README.md @@ -60,11 +60,20 @@ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=IT ## Requirements -- Python ≥ 3.8 -- IMAS Python Access Layer (`imas-python`) +- Python ≥ 3.10 + +### Installed automatically via pip - NumPy, Matplotlib, Pandas - Rich (for enhanced terminal output) +### Requires separate installation +- **Tkinter** — usually bundled with Python but may require system packages: + - Linux (Debian/Ubuntu): `sudo apt install python3-tk` + - Linux (RHEL/CentOS/Rocky): `sudo dnf install python3-tkinter` + - Windows: included in the [python.org](https://www.python.org/downloads/) installer ("tcl/tk and IDLE" component, enabled by default) + - macOS (python.org installer): included by default + - macOS (Homebrew): `brew install python-tk` (or `brew install python-tk@3.x` for a specific version) + ## Documentation Full documentation is available at the project repository. Each tool includes built-in help: diff --git a/ci-sdcc/st00-header.sh b/ci-sdcc/st00-header.sh index 5fa6ce35..72e2e8db 100644 --- a/ci-sdcc/st00-header.sh +++ b/ci-sdcc/st00-header.sh @@ -13,7 +13,7 @@ shopt -s expand_aliases #print hostname hostname -f -IMAS_EXISTS=$(module -r -t list 2>&1 | grep -E "IMAS-AL-Core|IMAS-Core" | head -n 1) +IMAS_EXISTS=$(module -r -t list 2>&1 | grep -E "IMAS-Core" | head -n 1) if [ -n "$IMAS_EXISTS" ]; then echo "> Found already loaded IMAS Module : $IMAS_EXISTS" ACCESS_LAYER_VERSION=$(echo "$AL_VERSION" | cut -d '.' -f 1) diff --git a/ci-sdcc/utils.sh b/ci-sdcc/utils.sh index 18708f25..101272f1 100644 --- a/ci-sdcc/utils.sh +++ b/ci-sdcc/utils.sh @@ -19,7 +19,7 @@ getIMASCoreModuleName() { DD_VERSION="$3" fi #Semantic versioning - IMASVERSIONSLIST=$(module -t avail IMAS-AL-Core/ 2>&1 | grep -E "^IMAS-AL-Core/$ACCESS_LAYER_VERSION\.[0-9]+\.[0-9]+-$TOOLCHAIN_VERSION") + IMASVERSIONSLIST=$(module -t avail IMAS-Core/ 2>&1 | grep -E "^IMAS-Core/$ACCESS_LAYER_VERSION\.[0-9]+\.[0-9]+-$TOOLCHAIN_VERSION") if [[ $TOOLCHAIN_VERSION == *"intel"* ]]; then CORE_MODULE_VERSION=$(echo "$IMASVERSIONSLIST" | grep "intel" | sort -rV | head -n 1) diff --git a/docs/source/cheatsheet.rst b/docs/source/cheatsheet.rst index c4557ff5..d3bf0f87 100644 --- a/docs/source/cheatsheet.rst +++ b/docs/source/cheatsheet.rst @@ -73,7 +73,7 @@ This cheat sheet provides quick reference of commonly used commands in IDStools. | | | | | $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134173;run=2326;database=TEST;version=3" | | | -md "imas:mdsplus?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active" | -| | "imas:hdf5?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" --rho | +| | "imas:hdf5?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" --phi | | | | +----------------------------+------------------------------------------------------------------------------------------------------------------------------+ | *plotequicomp* | Shows plasma equilibrium and quantities related with it | diff --git a/docs/source/idsprint.rst b/docs/source/idsprint.rst index a8c44544..d80a27a8 100644 --- a/docs/source/idsprint.rst +++ b/docs/source/idsprint.rst @@ -183,16 +183,18 @@ It can also save generated figures to a specified directory or the default locat $ idsprint --uri "imas:mdsplus?user=public;pulse=122481;run=1;database=ITER;version=3#edge_profiles/ggd[:]/electrons/density[1].values" -p - $ idsprint --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3#core_profiles/profiles_1d(:)/electrons/temperature" -p + $ idsprint --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3#core_profiles/profiles_1d(:)/electrons/temperature" -p - .. image:: _static/images/idsprint_1.png - :alt: image not found - :align: center + $ idsprint --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3#summary/time" - .. image:: _static/images/idsprint_2.png - :alt: image not found - :align: center +.. image:: _static/images/idsprint_1.png + :alt: image not found + :align: center - .. image:: _static/images/idsprint_3.png - :alt: image not found - :align: center \ No newline at end of file +.. image:: _static/images/idsprint_2.png + :alt: image not found + :align: center + +.. image:: _static/images/idsprint_3.png + :alt: image not found + :align: center diff --git a/docs/source/jupyter.rst b/docs/source/jupyter.rst new file mode 100644 index 00000000..9c48301e --- /dev/null +++ b/docs/source/jupyter.rst @@ -0,0 +1,51 @@ +####################### + Jupyter notebook usage +####################### + +IDStools command-line tools can also be used directly from Jupyter notebooks. +Import ``idstools`` once in the notebook kernel to register the IDStools +commands as IPython line magics: + +.. code-block:: python + + import idstools + +After that, commands such as ``idsprint`` and ``plotequilibrium`` can be +called with ``%``: + +.. code-block:: python + + %idsprint -u "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/134174/117#core_profiles/profiles_1d[0]/electrons/temperature" -p + +.. code-block:: python + + %plotequilibrium -u "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/100507/5" + + +****************************** + Interactive Matplotlib plots +****************************** + +For interactive Matplotlib figures in Jupyter, use the ``ipympl`` backend. If +``ipympl`` is installed in the same Python environment as IDStools, the backend +can be selected with Matplotlib's notebook magic before plotting: + +.. code-block:: python + + %matplotlib widget + import idstools + + %idsprint -u "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/134174/117#core_profiles/profiles_1d[0]/electrons/temperature" -p + +Alternatively, the backend can be selected through the IDStools ``--rc`` +option before Matplotlib has been imported in the current kernel: + +.. code-block:: python + + import idstools + + %idsprint -u "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/134174/117#core_profiles/profiles_1d[0]/electrons/temperature" -p --rc "backend='module://ipympl.backend_nbagg'" + +.. code-block:: python + + %plotequilibrium -u "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/100507/5" --rc "backend='module://ipympl.backend_nbagg'" diff --git a/docs/source/plotequilibrium.rst b/docs/source/plotequilibrium.rst index a5f3325b..e2bf5851 100644 --- a/docs/source/plotequilibrium.rst +++ b/docs/source/plotequilibrium.rst @@ -3,9 +3,17 @@ ################# *plotequilibrium* script shows plasma equilibrium. Optionally it also -shows pf coils position and toroidal flux. +shows PF coil positions and a toroidal-flux contour overlay. `refer data dictionary `_. +``--phi`` overlays contours of the toroidal magnetic flux Φ(R,Z) [Wb] read directly from +``equilibrium/time_slice/profiles_2d/phi``: + +.. math:: + + \Phi(R,Z) \quad [\mathrm{Wb}] + + ************************ Syntax plotequilibrium ************************ @@ -19,10 +27,11 @@ shows pf coils position and toroidal flux. .. code-block:: bash - $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3" --rho -md pf_active wall --plots - $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3" --rho -md "imas:mdsplus?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active" "imas:mdsplus?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" - $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134173;run=2326;database=TEST;version=3" --rho --md "imas:mdsplus?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active" "imas:hdf5?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" - + $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3" --phi -md pf_active wall --profiles + $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134174;run=117;database=ITER;version=3" --phi -md "imas:mdsplus?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active" "imas:mdsplus?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" + $ plotequilibrium --uri "imas:mdsplus?user=public;pulse=134173;run=2326;database=TEST;version=3" --phi --md "imas:mdsplus?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active" "imas:hdf5?user=public;pulse=116000;run=4;database=ITER_MD;version=3#wall" + $ plotequilibrium --uri "imas:hdf5?path=/work/imas/shared/imasdb/ITER/3/100507/5" --md "imas:hdf5?path=/work/imas/shared/imasdb/ITER_MD/3/116000/5#wall" --profiles --no-provenance + .. image:: _static/images/plotequilibrium.png :alt: image not found :align: center @@ -33,4 +42,4 @@ shows pf coils position and toroidal flux. .. image:: _static/images/plotequilibrium3.png :alt: image not found - :align: center \ No newline at end of file + :align: center diff --git a/docs/source/tools.rst b/docs/source/tools.rst index a90cc39d..1b5a1685 100644 --- a/docs/source/tools.rst +++ b/docs/source/tools.rst @@ -12,3 +12,4 @@ Following are the different command line tools available in the ids_manipulation_tools database_tools scenariodb_tools + jupyter diff --git a/idstools/compute/core_sources.py b/idstools/compute/core_sources.py index 534900c6..6dc3c9b5 100644 --- a/idstools/compute/core_sources.py +++ b/idstools/compute/core_sources.py @@ -365,8 +365,8 @@ def get_single_and_total_electrons_ions_waveforms(self, time_slice): # single_power_waveform = dict() # waveform # single_particles_waveform = dict() # waveform time_length = len(self.ids.time) - total_power_waveform = np.zeros(time_length) - total_particles_waveform = np.zeros(time_length) + total_ion_power_waveform = np.zeros(time_length) + total_ion_particles_waveform = np.zeros(time_length) single_power_waveform = {} single_particles_waveform = {} dict_single_and_total_electrons_waveforms = self.get_single_and_total_electrons_waveforms(time_slice) @@ -398,16 +398,14 @@ def get_single_and_total_electrons_ions_waveforms(self, time_slice): if total_ion_power < 0: total_ion_power = 0.0 - total_power_waveform[time_index] = ( - total_electron_power_waveform[time_index] + electrons_power + total_ion_power - ) - total_particles_waveform[time_index] = ( - total_electron_particles_waveform[time_index] + electrons_particles - ) + total_ion_particles + total_ion_power_waveform[time_index] += total_ion_power + total_ion_particles_waveform[time_index] += total_ion_particles single_power_waveform[source_index].append(electrons_power + total_ion_power) single_particles_waveform[source_index].append(electrons_particles + total_ion_particles) single_power_waveform[source_index] = np.array(single_power_waveform[source_index]) single_particles_waveform[source_index] = np.array(single_particles_waveform[source_index]) + total_power_waveform = total_electron_power_waveform + total_ion_power_waveform + total_particles_waveform = total_electron_particles_waveform + total_ion_particles_waveform return { "total_power_waveform": total_power_waveform, "total_particles_waveform": total_particles_waveform, diff --git a/idstools/compute/edge_profiles.py b/idstools/compute/edge_profiles.py index 715c3677..a82d82a5 100644 --- a/idstools/compute/edge_profiles.py +++ b/idstools/compute/edge_profiles.py @@ -952,6 +952,242 @@ def get_neutral_density(self, time_slice, x, y): n_neutral_edge = interpolate.griddata((r_edge, z_edge), temp, (x, y)) return n_neutral_edge + def get_cell_polygons_and_values(self, time_slice, values): + """ + Build a list of polygon vertex arrays and corresponding per-cell scalar values for + direct rendering with ``matplotlib.collections.PolyCollection``. + + Args: + time_slice: time index + values: 1-D array of scalar values, either per-cell or per-node. + + Returns: + tuple: ``(polygons, cell_values)`` where *polygons* is a list of ``(n_corners, 2)`` + float arrays (R, Z columns) and *cell_values* is a 1-D float array of the + same length. Returns ``(None, None)`` when the cells grid subset is unavailable. + """ + r_nodes, z_nodes = self.get_rz(time_slice) + values = np.asarray(values, dtype=float) + + cells_subset = None + for grid_subset in self.ids.grid_ggd[time_slice].grid_subset: + if grid_subset.identifier.name.lower() == "cells": + cells_subset = grid_subset + break + + if cells_subset is None or len(cells_subset.element) == 0: + logger.warning("edge_profiles: cells grid subset not found, cannot build PolyCollection") + return None, None + + n_cells = len(cells_subset.element) + per_cell_mode = len(values) == n_cells # True → direct 1:1, False → average from nodes + + polygons = [] + cell_values = [] + n_nodes = len(r_nodes) + i = 0 + for element in cells_subset.element: + for obj in element.object: + space_index = obj.space - 1 + dimension_index = obj.dimension - 1 + object_index = obj.index - 1 + nodes = ( + self.ids.grid_ggd[time_slice] + .space[space_index] + .objects_per_dimension[dimension_index] + .object[object_index] + .nodes + ) + nodes = np.asarray(nodes, dtype=int) - 1 # convert to 0-based + if np.any(nodes >= n_nodes) or np.any(nodes < 0): + logger.warning("edge_profiles: cell node index out of range, skipping cell") + i += 1 + continue + verts = np.column_stack((r_nodes[nodes], z_nodes[nodes])) + polygons.append(verts) + if per_cell_mode: + cell_values.append(float(values[i])) + else: + cell_values.append(float(np.mean(values[nodes]))) + i += 1 + + return polygons, np.array(cell_values, dtype=float) + + def get_electron_density_on_nodes(self, time_slice): + """ + Return raw electron density values at mesh nodes (no rectangular resampling). + + Args: + time_slice: time index + + Returns: + numpy array of node-centred electron density values, or ``None`` if unavailable. + """ + for electrons_density in self.ids.ggd[time_slice].electrons.density: + if electrons_density.grid_subset_index == 1: # nodes + return np.asarray(electrons_density.values) + logger.warning("edge_profiles: electron density values not found for nodes grid_subset") + return None + + def get_electron_density_on_cells(self, time_slice): + """ + Return electron density values. + + cell-centred values (``grid_subset_index == 5``, one value per cell. Falls + back to node values (``grid_subset_index == 1``) which will be averaged per cell further + + Args: + time_slice: time index + + Returns: + numpy array of density values, or ``None`` if unavailable. + """ + for ed in self.ids.ggd[time_slice].electrons.density: + if ed.grid_subset_index == 5: # cell-centred + return np.asarray(ed.values) + return self.get_electron_density_on_nodes(time_slice) + + def get_ion_density_on_nodes(self, time_slice): + """ + Return raw ion density values at mesh nodes (no rectangular resampling). + + Args: + time_slice: time index + + Returns: + numpy array of node-centred ion density values, or ``None`` if unavailable. + """ + for ion_density in self.ids.ggd[time_slice].ion[0].density: + if ion_density.grid_subset_index == 1: # nodes + return np.asarray(ion_density.values) + logger.warning("edge_profiles: ion density values not found for nodes grid_subset") + return None + + def get_ion_density_on_cells(self, time_slice): + """ + Return ion density values. + cell-centred values (``grid_subset_index == 5``), falls back to node values. + """ + for id_ in self.ids.ggd[time_slice].ion[0].density: + if id_.grid_subset_index == 5: + return np.asarray(id_.values) + return self.get_ion_density_on_nodes(time_slice) + + def get_neutral_density_on_nodes(self, time_slice): + """ + Return raw neutral density values at mesh nodes (no rectangular resampling). + + Args: + time_slice: time index + + Returns: + numpy array of node-centred neutral density values, or ``None`` if unavailable. + """ + for neutral_density in self.ids.ggd[time_slice].neutral[0].density: + if neutral_density.grid_subset_index == 1: # nodes + return np.asarray(neutral_density.values) + logger.warning("edge_profiles: neutral density values not found for nodes grid_subset") + return None + + def get_neutral_density_on_cells(self, time_slice): + """ + Return neutral density values. + cell-centred values (``grid_subset_index == 5``), falls back to node values. + """ + for nd in self.ids.ggd[time_slice].neutral[0].density: + if nd.grid_subset_index == 5: + return np.asarray(nd.values) + return self.get_neutral_density_on_nodes(time_slice) + + def get_rz_and_ne_on_subset(self, time_slice, subset_identifier_index): + """ + Return R, Z coordinates and electron density for every element in a named GGD grid subset. + + For 0-D elements (nodes) the node geometry is used directly. For higher-dimension + elements (edges, faces, cells) the centroid of the element's nodes is used. + + Args: + time_slice: time index + subset_identifier_index: GGD ``identifier.index`` of the target subset, + e.g. 11 = outer midplane, 13 = inner divertor target, 14 = outer divertor target. + + Returns: + ``(r, z, ne)`` numpy float arrays in grid-subset element order, + or ``(None, None, None)`` when the subset or its density data is unavailable. + """ + subset_array_index = None + for iset, gs in enumerate(self.ids.grid_ggd[time_slice].grid_subset): + if gs.identifier.index == subset_identifier_index: + subset_array_index = iset + break + + if subset_array_index is None: + logger.warning(f"edge_profiles: grid subset with identifier.index={subset_identifier_index} not found") + return None, None, None + + grid_subset = self.ids.grid_ggd[time_slice].grid_subset[subset_array_index] + if len(grid_subset.element) == 0: + logger.warning(f"edge_profiles: grid subset {subset_identifier_index} has no elements") + return None, None, None + + # grid_subset_index in ggd arrays stores identifier.index (semantic index), not array position + ne_entry = None + ne_all_nodes = None + for dens in self.ids.ggd[time_slice].electrons.density: + if dens.grid_subset_index == subset_identifier_index: + ne_entry = dens + break + if dens.grid_subset_index == 1: # nodes — keep as fallback + ne_all_nodes = np.asarray(dens.values, dtype=float) + + spaces = self.ids.grid_ggd[time_slice].space + r_list, z_list, ne_list = [], [], [] + + for element in grid_subset.element: + for obj in element.object: + space_idx = obj.space - 1 # space + dim_idx = obj.dimension - 1 # 0=nodes, 1=edges, 2=faces, 3=cells + obj_idx = obj.index - 1 # convert to 0-based + space = spaces[space_idx] + n_nodes = len(space.objects_per_dimension[0].object) + ggd_obj = space.objects_per_dimension[dim_idx].object[obj_idx] + if dim_idx == 0: + g = ggd_obj.geometry + r_list.append(float(g[0])) + z_list.append(float(g[1])) + if ne_entry is None and ne_all_nodes is not None: + ne_list.append(float(ne_all_nodes[obj_idx]) if obj_idx < n_nodes else np.nan) + else: + nodes = np.asarray(ggd_obj.nodes, dtype=int) - 1 # 0-based node indices + valid = nodes[(nodes >= 0) & (nodes < n_nodes)] + geom = [space.objects_per_dimension[0].object[n].geometry for n in valid] + r_list.append(float(np.mean([g[0] for g in geom]))) + z_list.append(float(np.mean([g[1] for g in geom]))) + if ne_entry is None and ne_all_nodes is not None: + valid_ne = valid[valid < len(ne_all_nodes)] + ne_list.append(float(np.mean(ne_all_nodes[valid_ne])) if len(valid_ne) else np.nan) + + r = np.array(r_list, dtype=float) + z = np.array(z_list, dtype=float) + + if ne_entry is not None: + ne_values = np.asarray(ne_entry.values, dtype=float) + if len(ne_values) != len(r): + logger.warning( + f"edge_profiles: per-subset density count {len(ne_values)} != element count {len(r)} " + f"for subset {subset_identifier_index}, falling back to node interpolation" + ) + if len(ne_list) == len(r): + return r, z, np.array(ne_list, dtype=float) + return None, None, None + return r, z, ne_values + + if len(ne_list) == len(r): + return r, z, np.array(ne_list, dtype=float) + + logger.warning(f"edge_profiles: electron density not found for grid subset {subset_identifier_index}") + return None, None, None + def get_outer_midplane_array_index(self, time_slice): """ This function `get_outer_midplane_array_index` searches for a specific grid subset with an diff --git a/idstools/compute/equilibrium.py b/idstools/compute/equilibrium.py index d103121b..5f84d4cb 100644 --- a/idstools/compute/equilibrium.py +++ b/idstools/compute/equilibrium.py @@ -16,6 +16,8 @@ from idstools.database import DBMaster +_IDS_VALID_THRESHOLD = abs(imas.ids_defs.EMPTY_FLOAT) + logger = logging.getLogger("module") @@ -41,90 +43,72 @@ def get2d_cartesian_grid(self, time_slice: int, profiles2d_index: int = 0) -> Un This function returns a dictionary containing 2D Cartesian grid coordinates and psi values from an equilibrium IDS object. - Args: time_slice (int): The time slice index of the equilibrium data to be used for generating the - 2D Cartesian grid. Defaults to 0 - profiles2d_index (int): `profiles2d_index` is an integer parameter that represents the index of the - ``profile_2d`` to be used in the calculation. It is used to access the specific 2D profile from the - list of profiles in the `time_slice` object. Defaults to 0 + 2D Cartesian grid. Defaults to 0 + profiles2d_index (int): An integer parameter that represents the index of the + ``profiles_2d`` to be used in the calculation. It is used to access the specific 2D profile from the + list of profiles in the `time_slice` object. Defaults to 0 Returns: A dictionary containing the 2D Cartesian grid coordinates (r2d and z2d) and the corresponding psi - values (psi2d). + values (psi2d), or None if the data is unavailable or invalid. Example: .. code-block:: python - import imas - connection = imas.DBEntry("imas:mdsplus?user=public;pulse=134173;run=106;database=ITER;version=3","r") - idsObj = connection.get('equilibrium') - computeObj = EquilibriumCompute(idsObj) - result = computeObj.get2d_cartesian_grid(time_slice=0) + import imas + connection = imas.DBEntry("imas:mdsplus?user=public;pulse=134173;run=106;database=ITER;version=3","r") + idsObj = connection.get('equilibrium') + computeObj = EquilibriumCompute(idsObj) + result = computeObj.get2d_cartesian_grid(time_slice=0) - {'psi2d': array([[]]), - 'r2d': array([[]]), - 'z2d': array([[]])} + {'r2d': array([...]), 'z2d': array([...]), 'psi2d': array([...])} """ - profiles2d = None + profiles2d = r1d = z1d = None try: - profiles2d = self.ids.time_slice[time_slice].profiles_2d[ - profiles2d_index - ] # using https://docs.python.org/2/glossary.html#term-eafp style + profiles2d = self.ids.time_slice[time_slice].profiles_2d[profiles2d_index] except IndexError: logger.error(f"equilibrium.time_slice[{time_slice}].profiles_2d[{profiles2d_index}] is not available") return None profiles2d = self.ids.time_slice[time_slice].profiles_2d[profiles2d_index] - r2d = profiles2d.r - z2d = profiles2d.z - psi2d = profiles2d.psi - if profiles2d.grid_type.index == 1 and np.size(r2d) == 0: - logger.warning( - f"profiles_2d[{profiles2d_index}].r is not available and grid type is 1.. Calculating from grid" + if profiles2d.grid_type.index == 1 and profiles2d.grid.dim1 is not None and profiles2d.grid.dim2 is not None: + logger.info( + f"Using equilibrium.time_slice[{time_slice}]" + f".profiles_2d[{profiles2d_index}].grid.dim1/dim2 for the 2D grid" ) r1d = profiles2d.grid.dim1 z1d = profiles2d.grid.dim2 - nr = len(r1d) - nz = len(z1d) - r2d = np.empty(shape=(nr, nz)) - z2d = np.empty(shape=(nr, nz)) - for iz in range(nz): - r2d[:, iz] = r1d - for ir in range(nr): - z2d[ir, :] = z1d + else: + logger.error("Only rectangular cylindrical grid (grid_type=1) is supported for now") + return None + + psi2d = profiles2d.psi if np.all(psi2d == 0.0): logger.error( "All values of psi2d are 0. No contour levels were found within the data range, Can not plot contour" ) return None - if np.size(r2d) != np.size(z2d) or np.size(r2d) != np.size(psi2d): - logger.error( - f"r, z and psi have not the same dimension in \ - equilibrium.time_slice[{time_slice}].profiles_2d[{profiles2d_index}]" - ) - return None - return {"r2d": r2d, "z2d": z2d, "psi2d": psi2d} + return {"r2d": r1d, "z2d": z1d, "psi2d": psi2d} - def get_rho2d(self, time_slice: int, profiles2d_index: int = 0) -> Union[np.ndarray, None]: + def get_phi2d(self, time_slice: int, profiles2d_index: int = 0) -> Union[np.ndarray, None]: """ - This function calculates rho(R,Z) using toroidal flux and returns a dictionary containing the result. + Returns the toroidal magnetic flux Φ(R,Z) on the 2D grid. + + Reads ``equilibrium.time_slice[i].profiles_2d[j].phi`` directly from the IDS. Args: - time_slice (int): The time slice is an integer value that represents the index of the time slice in - the equilibrium ids. It is used to select a specific time slice for the calculation of rho(R,Z). - Defaults to 0 - profiles2d_index (int): `profiles2d_index` is an integer parameter that represents the index of the - ``profiles_2d`` to be used for the calculation of rho(R,Z). It is used to access the `profiles_2d` - list in the `time_slice` object. Defaults to 0 + time_slice (int): Index of the time slice in the equilibrium IDS. Defaults to 0. + profiles2d_index (int): Index into ``profiles_2d`` from which + ``phi`` (toroidal flux, Wb) is read. Defaults to 0. Returns: - a value containing the square root of the toroidal flux values divided by the maximum toroidal - flux value, if the length of toroidal flux is greater than 0. If the length of toroidal flux is - less than 1, it returns None. + np.ndarray or None: 2-D array of toroidal flux Φ [Wb] with the same shape as + the ``profiles_2d`` grid, or None if ``phi`` is unavailable or all-NaN. Examples: .. code-block:: python @@ -133,7 +117,7 @@ def get_rho2d(self, time_slice: int, profiles2d_index: int = 0) -> Union[np.ndar connection = imas.DBEntry("imas:mdsplus?user=public;pulse=134173;run=106;database=ITER;version=3", "r") idsObj = connection.get('equilibrium') computeObj = EquilibriumCompute(idsObj) - result = computeObj.get_rho2d(time_slice=0) + result = computeObj.get_phi2d(time_slice=0) """ phi = None @@ -150,7 +134,7 @@ def get_rho2d(self, time_slice: int, profiles2d_index: int = 0) -> Union[np.ndar f"all values are nan for equilibrium.time_slice[{time_slice}].profiles_2d[{profiles2d_index}].phi " ) return None - return np.sqrt(phi / np.amax(phi)) + return phi def get_b_total(self, time_slice: int) -> tuple: """ @@ -259,8 +243,8 @@ def get_flux_surfaces(self, time_slice: int) -> dict: Returns: a dictionary containing information about flux surfaces at a specific time slice. The dictionary includes - a 2D Cartesian grid, a 2D profile index, and a 2D array of rho values. If no profiles are found, - the function returns None. + a 2D Cartesian grid, a 2D profile index, and a 2D array of rho_tor_norm [-] values (dimensionless, + range [0, 1]). If no profiles are found, the function returns None. """ GRID_TYPE_RECTANGULAR = 1 list_of_profiles = self.get2d_profiles_indices(time_slice, GRID_TYPE_RECTANGULAR) @@ -271,10 +255,10 @@ def get_flux_surfaces(self, time_slice: int) -> dict: profile2d_index = list_of_profiles[0] result_dict = self.get2d_cartesian_grid(time_slice, profile2d_index) - rho2d = self.get_rho2d(time_slice, profile2d_index) - if rho2d is None: - rho2d = [] - result_dict["rho2d"] = rho2d + phi2d = self.get_phi2d(time_slice, profile2d_index) + if phi2d is None: + phi2d = [] + result_dict["phi2d"] = phi2d return result_dict def get_ip(self) -> list: @@ -301,6 +285,408 @@ def get_ip(self) -> list: for time_index in range(len(self.ids.time_slice)) ] + def get_boundary_data(self, time_slice: int) -> dict: + """Return boundary data for a given time slice. + + Reads ``boundary/outline``, ``boundary_separatrix`` (DD3), or + ``contour_tree`` (DD4) for the separatrix outline, X-points, and + strike-points. If the separatrix is still missing, falls back to + ``boundary/outline`` for diverted plasmas (``type==1``) or + ``boundary/lcfs`` for limiter/unknown. + + Returns a dict with keys ``bnd_r``, ``bnd_z``, ``bnd_type``, + ``bnd_psi_norm``, ``bnd_geom_r``, ``bnd_geom_z``, ``sep_r``, + ``sep_z``, ``sep_xpoints``, ``sep_strikepoints``. + """ + + def _valid_arr(arr): + a = np.asarray(arr, dtype=float) + return a.size > 0 and np.any(np.isfinite(a) & (np.abs(a) < _IDS_VALID_THRESHOLD)) + + def _valid_scalar(val): + try: + v = float(val) + return np.isfinite(v) and abs(v) < _IDS_VALID_THRESHOLD + except Exception as exc: + logger.debug(f"get_boundary_data: invalid scalar {val!r} ({exc})") + return False + + def _clean(arr): + a = np.array(arr, dtype=float, copy=True) + a[(~np.isfinite(a)) | (np.abs(a) >= _IDS_VALID_THRESHOLD)] = np.nan + return a + + def _read_outline(node): + try: + r = np.asarray(node.outline.r, dtype=float) + z = np.asarray(node.outline.z, dtype=float) + except Exception as exc: + logger.debug(f"get_boundary_data: could not read outline from {node!r}: {exc}") + return None, None + if not (_valid_arr(r) and _valid_arr(z)): + logger.debug("get_boundary_data: outline has no valid data " f"(r.size={r.size}, z.size={z.size})") + return None, None + r, z = _clean(r), _clean(z) + # Insert NaN at large jumps so disconnected arcs are not joined + dist = np.sqrt(np.diff(r) ** 2 + np.diff(z) ** 2) + median_dist = np.nanmedian(dist) + if median_dist > 0: + breaks = np.where(dist > 20.0 * median_dist)[0] + 1 + if len(breaks): + r = np.insert(r, breaks, np.nan) + z = np.insert(z, breaks, np.nan) + return r, z + + def _read_points(node, attr, ids_path): + pts = [] + try: + arr = getattr(node, attr) + except AttributeError: + logger.debug(f"get_boundary_data: {ids_path}/{attr} is not available") + return pts + except Exception as exc: + logger.debug(f"get_boundary_data: could not access {ids_path}/{attr}: {exc}") + return pts + try: + n_points = len(arr) + except Exception as exc: + logger.debug(f"get_boundary_data: could not get length of {ids_path}/{attr}: {exc}") + n_points = None + for pt_index, pt in enumerate(arr): + try: + r, z = float(pt.r), float(pt.z) + except Exception as exc: + logger.debug(f"get_boundary_data: could not read {ids_path}/{attr}[{pt_index}]/r|z: {exc}") + continue + if _valid_scalar(r) and _valid_scalar(z): + pts.append((r, z)) + else: + logger.debug(f"get_boundary_data: {ids_path}/{attr}[{pt_index}]/r|z invalid ({r}, {z})") + logger.debug(f"get_boundary_data: {ids_path}/{attr} — read {len(pts)} valid points out of {n_points}") + return pts + + def _read_contour_tree(ts_node): + """Read separatrix/X-point data from ``time_slice.contour_tree.node``. + + * ``node.critical_type == 1`` for X-points (saddle points) + * first valid X-point ``node.levelset.r/z`` as separatrix contour + """ + sep_r = sep_z = None + xpoints = [] + + try: + nodes = ts_node.contour_tree.node + except Exception as exc: + logger.debug(f"get_boundary_data: could not access contour_tree.node: {exc}") + return sep_r, sep_z, xpoints + + try: + n_nodes = len(nodes) + except Exception as exc: + logger.debug(f"get_boundary_data: could not get length of contour_tree.node: {exc}") + n_nodes = None + + n_saddles = 0 + for node_index, node in enumerate(nodes): + try: + critical_type = int(node.critical_type) + except Exception as exc: + logger.debug( + f"get_boundary_data: could not read contour_tree.node[{node_index}].critical_type: {exc}" + ) + continue + + if critical_type != 1: # 1 = saddle / X-point + continue + n_saddles += 1 + + try: + xr = float(node.r) + xz = float(node.z) + except Exception as exc: + logger.debug(f"get_boundary_data: could not read contour_tree.node[{node_index}].r/z: {exc}") + xr = xz = None + + if xr is not None and _valid_scalar(xr) and xz is not None and _valid_scalar(xz): + xpoints.append((xr, xz)) + else: + logger.debug( + f"get_boundary_data: contour_tree.node[{node_index}] saddle has invalid r/z " f"({xr}, {xz})" + ) + + if sep_r is not None and sep_z is not None: + continue + + try: + r = np.asarray(node.levelset.r, dtype=float) + z = np.asarray(node.levelset.z, dtype=float) + except Exception as exc: + logger.debug( + f"get_boundary_data: could not read contour_tree.node[{node_index}].levelset.r/z: {exc}" + ) + continue + + if not (_valid_arr(r) and _valid_arr(z)): + logger.debug( + f"get_boundary_data: contour_tree.node[{node_index}].levelset has no valid data " + f"(r.size={r.size}, z.size={z.size})" + ) + continue + + sep_r = _clean(r) + sep_z = _clean(z) + + logger.debug( + "get_boundary_data: contour_tree summary " + f"(nodes={n_nodes}, saddles={n_saddles}, xpoints={len(xpoints)}, " + f"has_separatrix={sep_r is not None and sep_z is not None})" + ) + return sep_r, sep_z, xpoints + + result = { + "bnd_r": None, + "bnd_z": None, + "bnd_type": None, + "bnd_psi_norm": None, + "bnd_geom_r": None, + "bnd_geom_z": None, + "sep_r": None, + "sep_z": None, + "sep_xpoints": [], + "sep_strikepoints": [], + } + + try: + ts = self.ids.time_slice[time_slice] + except Exception as exc: + logger.debug(f"get_boundary_data: could not access time_slice[{time_slice}]: {exc}") + return result + + # boundary + try: + bnd = ts.boundary + result["bnd_r"], result["bnd_z"] = _read_outline(bnd) + result["sep_xpoints"] = _read_points(bnd, "x_point", f"time_slice[{time_slice}]/boundary") + result["sep_strikepoints"] = _read_points(bnd, "strike_point", f"time_slice[{time_slice}]/boundary") + logger.debug( + f"get_boundary_data: time_slice[{time_slice}]/boundary summary " + f"(has_outline={result['bnd_r'] is not None and result['bnd_z'] is not None}, " + f"xpoints={len(result['sep_xpoints'])}, strikepoints={len(result['sep_strikepoints'])})" + ) + + bnd_type = int(bnd.type) + if _valid_scalar(bnd_type): + result["bnd_type"] = bnd_type + except Exception as exc: + logger.debug(f"get_boundary_data: could not read time_slice[{time_slice}]/boundary: {exc}") + + try: + psi_norm = float(ts.boundary.psi_norm) + if _valid_scalar(psi_norm): + result["bnd_psi_norm"] = psi_norm + except Exception as exc: + logger.debug(f"get_boundary_data: could not read time_slice[{time_slice}]/boundary/psi_norm: {exc}") + + try: + gax_r = float(ts.boundary.geometric_axis.r) + gax_z = float(ts.boundary.geometric_axis.z) + if _valid_scalar(gax_r) and _valid_scalar(gax_z): + result["bnd_geom_r"] = gax_r + result["bnd_geom_z"] = gax_z + except Exception as exc: + logger.debug( + f"get_boundary_data: could not read time_slice[{time_slice}]/boundary/geometric_axis/r|z: {exc}" + ) + + # boundary_separatrix (DD3 ) + if hasattr(ts, "boundary_separatrix"): + sep = ts.boundary_separatrix + try: + result["sep_r"], result["sep_z"] = _read_outline(sep) + sep_xpoints = _read_points(sep, "x_point", f"time_slice[{time_slice}]/boundary_separatrix") + sep_strikepoints = _read_points(sep, "strike_point", f"time_slice[{time_slice}]/boundary_separatrix") + if sep_xpoints: + result["sep_xpoints"] = sep_xpoints + if sep_strikepoints: + result["sep_strikepoints"] = sep_strikepoints + logger.debug( + f"get_boundary_data: time_slice[{time_slice}]/boundary_separatrix summary " + f"(has_outline={result['sep_r'] is not None and result['sep_z'] is not None}, " + f"xpoints={len(sep_xpoints)}, strikepoints={len(sep_strikepoints)})" + ) + except Exception as exc: + logger.debug(f"get_boundary_data: could not read time_slice[{time_slice}]/boundary_separatrix: {exc}") + + # contour_tree.node (DD4) + if hasattr(ts, "contour_tree") and hasattr(ts.contour_tree, "node"): + contour_sep_r, contour_sep_z, contour_xpoints = _read_contour_tree(ts) + + if ( + (result["sep_r"] is None or result["sep_z"] is None) + and contour_sep_r is not None + and contour_sep_z is not None + ): + result["sep_r"] = contour_sep_r + result["sep_z"] = contour_sep_z + + if not result["sep_xpoints"] and contour_xpoints: + result["sep_xpoints"] = contour_xpoints + + # Separatrix fallback when boundary_separatrix / contour_tree provided nothing. + if result["sep_r"] is None or result["sep_z"] is None: + if result["bnd_type"] == 1: + # type=1 (diverted): boundary/outline IS the separatrix — reuse directly. + if result["bnd_r"] is not None and result["bnd_z"] is not None: + result["sep_r"] = result["bnd_r"] + result["sep_z"] = result["bnd_z"] + logger.debug( + f"get_boundary_data: time_slice[{time_slice}]/boundary/outline/r|z " + f"— sep outline reused (type=1 diverted, {result['sep_r'].size} pts)" + ) + else: + # type=0 (limiter) or unknown: outline is the limiter contour, not the LCFS. + # Fall back to boundary/lcfs + try: + r_raw = np.asarray(ts.boundary.lcfs.r, dtype=float) + z_raw = np.asarray(ts.boundary.lcfs.z, dtype=float) + mask = r_raw > 0 + r_raw, z_raw = _clean(r_raw[mask]), _clean(z_raw[mask]) + if r_raw.size > 0: + result["sep_r"] = r_raw + result["sep_z"] = z_raw + logger.debug( + f"get_boundary_data: time_slice[{time_slice}]/boundary/lcfs/r|z " + f"— sep outline filled ({r_raw.size} pts)" + ) + except Exception as exc: + logger.debug(f"get_boundary_data: could not read time_slice[{time_slice}]/boundary/lcfs/r|z: {exc}") + + logger.debug( + "get_boundary_data: final summary " + f"(has_boundary={result['bnd_r'] is not None and result['bnd_z'] is not None}, " + f"has_separatrix={result['sep_r'] is not None and result['sep_z'] is not None}, " + f"xpoints={len(result['sep_xpoints'])}, strikepoints={len(result['sep_strikepoints'])})" + ) + + return result + + def get_magnetic_axis(self, time_slice: int) -> Union[dict, None]: + """Return the magnetic axis position for a given time slice. + + Reads ``global_quantities.magnetic_axis.r/z`` and validates the + scalar values. + + Args: + time_slice (int): Index into ``time_slice``. + + Returns: + dict with scalar keys ``"r"`` and ``"z"`` (floats), or + ``None`` if the data are absent or invalid. + """ + try: + mag_ax = self.ids.time_slice[time_slice].global_quantities.magnetic_axis + r = float(mag_ax.r) + z = float(mag_ax.z) + except Exception as exc: + logger.debug(f"get_magnetic_axis: could not read magnetic_axis – {exc}") + return None + + def _valid(val): + return np.isfinite(val) and abs(val) < _IDS_VALID_THRESHOLD + + if not (_valid(r) and _valid(z)): + logger.debug("get_magnetic_axis: magnetic_axis contains no valid data") + return None + + return {"r": r, "z": z} + + def get_current_centre(self, time_slice: int) -> Union[dict, None]: + """Return the current centroid position for a given time slice. + + Reads ``global_quantities.current_centre.r/z`` and validates the + scalar values. + + Args: + time_slice (int): Index into ``time_slice``. + + Returns: + dict with scalar keys ``"r"`` and ``"z"`` (floats), or + ``None`` if the data are absent or invalid. + """ + try: + cc = self.ids.time_slice[time_slice].global_quantities.current_centre + r = float(cc.r) + z = float(cc.z) + except Exception as exc: + path = f"time_slice[{time_slice}]/global_quantities/current_centre/r|z" + logger.debug(f"get_current_centre: could not read {path} – {exc}") + return None + + def _valid(val): + return np.isfinite(val) and abs(val) < _IDS_VALID_THRESHOLD + + if not (_valid(r) and _valid(z)): + path = f"time_slice[{time_slice}]/global_quantities/current_centre/r|z" + logger.debug(f"get_current_centre: {path} contains no valid data") + return None + + return {"r": r, "z": z} + + def get_scalar_annotation_quantities(self, time_slice: int) -> list: + """Return validated scalar global/boundary quantities for annotation display. + + Reads a fixed set of scalar fields from ``global_quantities`` and + ``boundary``, validates each value (finite and below the IDS fill + value threshold), and returns + only those with valid data. + + Args: + time_slice (int): Index into ``time_slice``. + + Returns: + list of dicts, each with ``"label"`` (LaTeX str) and ``"text"`` + (formatted value + unit str). Empty list if nothing is valid. + """ + + def _valid(val): + try: + v = float(val) + return np.isfinite(v) and abs(v) < _IDS_VALID_THRESHOLD + except Exception: + return False + + items = [] + ts = self.ids.time_slice[time_slice] + gq = ts.global_quantities + bnd = ts.boundary + + _specs = [ + (lambda: float(gq.ip), lambda v: {"label": "$I_p$", "text": f"{v / 1e6:.3f} MA"}), + ( + lambda: float( + getattr( + gq.magnetic_axis, "b_field_phi" if hasattr(gq.magnetic_axis, "b_field_phi") else "b_field_tor" + ) + ), + lambda v: {"label": r"$B_\phi$(axis)", "text": f"{v:.3f} T"}, + ), + (lambda: float(gq.psi_axis), lambda v: {"label": r"$\psi_{\rm axis}$", "text": f"{v:.4g} Wb"}), + (lambda: float(gq.psi_boundary), lambda v: {"label": r"$\psi_{\rm bnd}$", "text": f"{v:.4g} Wb"}), + (lambda: float(gq.q_axis), lambda v: {"label": "$q_0$", "text": f"{v:.3f}"}), + (lambda: float(gq.q_95), lambda v: {"label": "$q_{95}$", "text": f"{v:.3f}"}), + (lambda: float(bnd.minor_radius), lambda v: {"label": "$a$", "text": f"{v:.3f} m"}), + (lambda: float(bnd.elongation), lambda v: {"label": r"$\kappa$", "text": f"{v:.3f}"}), + (lambda: float(bnd.triangularity), lambda v: {"label": r"$\delta$", "text": f"{v:.3f}"}), + ] + for getter, formatter in _specs: + try: + val = getter() + if _valid(val): + items.append(formatter(val)) + except Exception: + pass + return items + def get_top_view(self, time_slice: int) -> dict: """ The function returns data for plotting the top view of a 2D shape. @@ -998,7 +1384,7 @@ def get_profiles_1d_quantities(self, time_slice, attributes=None): if ids_field.has_value: quantities[attribute] = eval(f"self.ids.time_slice[{time_slice}].profiles_1d.{attribute}") else: - logger.error(f"self.ids.time_slice[{time_slice}].profiles_1d.{attribute} not found") + logger.warning(f"self.ids.time_slice[{time_slice}].profiles_1d.{attribute} not found") return quantities def get_global_quantities(self, time_slice=None, attributes=None): @@ -1040,7 +1426,7 @@ def get_global_quantities(self, time_slice=None, attributes=None): node = eval(f"self.ids.time_slice[{ti}].global_quantities.{attribute}") if info_flag: quantities[attribute]["unit"] = node.metadata.units - quantities[attribute]["coordinate_unit"] = "t" + quantities[attribute]["coordinate_unit"] = self.ids.time.metadata.units or "s" quantities[attribute]["name"] = node.metadata.name quantities[attribute]["coordinate_name"] = "time" @@ -1155,6 +1541,9 @@ def get_equilibria(self, selection=None): name = self.ids.code.name if homogeneous_time == 1: time = self.ids.time + else: + # Extract time from time_slice array + time = np.array([ts.time for ts in self.ids.time_slice]) if len(self.ids.time_slice) > 0 else np.array([]) nt = time.size data = {} @@ -1378,10 +1767,10 @@ def get_equilibria(self, selection=None): z = time_slice.profiles_2d[0].grid.dim2 # Initialize boundary arrays - each time slice can have different size - if need_boundaries and n3 > 0: - rb = [] if "rb" in selection else None - zb = [] if "zb" in selection else None + rb = [] if "rb" in selection else None + zb = [] if "zb" in selection else None + if need_boundaries: for i, time_slice in enumerate(self.ids.time_slice): if time_slice.boundary.outline.r.size > 0 and rb is not None: rb.append(time_slice.boundary.outline.r) diff --git a/idstools/compute/summary.py b/idstools/compute/summary.py index 0847002b..f45277c4 100644 --- a/idstools/compute/summary.py +++ b/idstools/compute/summary.py @@ -90,7 +90,7 @@ def get_summary( # The `energy_mhd` variable in the `get_summary` method is used to retrieve the "energy_mhd" # waveform values from the `global_quantities` of the summary ids object. If the length of the # retrieved values is less than 1, it logs a critical error message indicating that the data - # could not be read and assigns an array of NaN values with the same length as the time vector + # could not be read and assigns an array of nan values with the same length as the time vector # to the `energy_mhd` variable. This ensures that the `energy_mhd` waveform data is handled # properly even in the case of missing or unreadable data. energy_mhd = self.ids.global_quantities.energy_mhd.value diff --git a/idstools/database.py b/idstools/database.py index bb96ae12..a882f05f 100644 --- a/idstools/database.py +++ b/idstools/database.py @@ -10,6 +10,8 @@ import imas import yaml +from idstools.utils.utility_functions import add_query_to_uri + logger = logging.getLogger(f"module.{__name__}") @@ -701,11 +703,16 @@ def create_connection(cls, imasargs, target_dd_version=None): def get_connection(cls, imasargs): connection = None if imasargs.uri != "" and imasargs.uri is not None: + uri = add_query_to_uri( + imasargs.uri, + backend="uda", + query="cache_mode=none", + ) if "mode" in imasargs.__dict__: - connection = imas.DBEntry(imasargs.uri, imasargs.mode) + connection = imas.DBEntry(uri, imasargs.mode) else: try: - connection = imas.DBEntry(imasargs.uri, "r") + connection = imas.DBEntry(uri, "r") except Exception as e: print(e) return connection @@ -954,14 +961,14 @@ def read_scenario_with_args( if test_mode: ids = connection.get_slice(ids_name, test_args_list) else: - ids = connection.get(ids_name) + ids = connection.get(ids_name, autoconvert=False) in_ids_dict[ids_name] = ids for ids_name in out_ids_list: if test_mode: ids = connection.get_slice(ids_name, test_args_list) else: - ids = connection.get(ids_name) + ids = connection.get(ids_name, autoconvert=False) out_ids_dict[ids_name] = ids connection.close() diff --git a/idstools/domain/ecstray.py b/idstools/domain/ecstray.py index 4824f9e4..0a55e2e5 100644 --- a/idstools/domain/ecstray.py +++ b/idstools/domain/ecstray.py @@ -21,7 +21,7 @@ def __init__(self, equilibrium_ids: object, core_profiles_ids: object, waves_ids # self.coreProfilesCompute = coreProfilesIds self.waves_compute = WavesCompute(waves_ids) - def get_resonance_layer(self, coherent_wave_index, time_slice, n_harm=None): + def get_resonance_layer(self, coherent_wave_index, time_slice, n_harm=None, equilibrium_time_slice=None): """This function calculates and returns a dictionary (Resonance Layer) containing r and z values corresponding to the resonance points based on the provided nHarm values, b_resonance, and b_total arrays. @@ -56,11 +56,14 @@ def get_resonance_layer(self, coherent_wave_index, time_slice, n_harm=None): """ if n_harm is None: n_harm = [1, 2, 3, 4] + if equilibrium_time_slice is None: + equilibrium_time_slice = time_slice + b_resonance = self.waves_compute.get_b_resonance(coherent_wave_index, time_slice, harmonic_frequencies=n_harm) - profile2d_index, b_total = self.equilibrium_compute.get_b_total(time_slice) + profile2d_index, b_total = self.equilibrium_compute.get_b_total(equilibrium_time_slice) if profile2d_index != -99: - r = self.equilibrium_compute.ids.time_slice[time_slice].profiles_2d[profile2d_index].grid.dim1 - z = self.equilibrium_compute.ids.time_slice[time_slice].profiles_2d[profile2d_index].grid.dim2 + r = self.equilibrium_compute.ids.time_slice[equilibrium_time_slice].profiles_2d[profile2d_index].grid.dim1 + z = self.equilibrium_compute.ids.time_slice[equilibrium_time_slice].profiles_2d[profile2d_index].grid.dim2 [nr, nz] = np.shape(b_total) b_err = 10 / nr @@ -177,8 +180,8 @@ def get_cutoff_layer(self, coherent_wave_index, time_slice): ) + constants.e * b_total[ir, iz] / (2 * constants.m_e) except Exception as e: # Not defined outside LCFS logger.debug(f"{e}") - ne2d_eq[ir, iz] = -1 # np.NaN - omega_r[ir, iz] = -1 # np.NaN + ne2d_eq[ir, iz] = -1 # np.nan + omega_r[ir, iz] = -1 # np.nan # Find (R,Z) where omega_R = omega_EC (within the tolerance omega_err) [nr, nz] = np.shape(omega_r) diff --git a/idstools/domain/kineticprofiles.py b/idstools/domain/kineticprofiles.py index 80081501..6d39b90a 100644 --- a/idstools/domain/kineticprofiles.py +++ b/idstools/domain/kineticprofiles.py @@ -787,7 +787,7 @@ def get_zeff_profile(self): if len(_zeff) != self.nrho: logger.warning("core_profiles.profiles_1d[:].zeff could not be read.") logger.warning(f"Size mismatch: rho_tor_norm = {self.nrho}, " f"zeff = {len(_zeff)}") - _zeff = np.asarray([np.NaN] * self.nrho) + _zeff = np.asarray([np.nan] * self.nrho) for i in range(self.nrho): zeff[i] = _zeff[i] if self.is_edge_profiles_present: @@ -795,14 +795,14 @@ def get_zeff_profile(self): _zeff = self.edge_profiles.profiles_1d[self.time_index_edge_profiles].zeff.value if len(_zeff) < 1: logger.warning("edge_profiles.profiles_1d[:].zeff could not be read.") - _zeff = np.asarray([np.NaN] * self.erho) + _zeff = np.asarray([np.nan] * self.erho) for i in range(self.erho): zeff[self.nrho + i] = _zeff[i] else: _gset = self.edge_profiles.ggd[self.time_index_edge_profiles].zeff[self.gset].values if len(_gset) < 1: logger.warning("edge_profiles.ggd[:].zeff could not be read.") - _gset = np.asarray([np.NaN] * self.erho) + _gset = np.asarray([np.nan] * self.erho) for i in range(self.erho): zeff[self.nrho + i] = _gset[i] return zeff @@ -814,7 +814,7 @@ def getne_profile(self): Returns: The function `getne_profile` returns a list of electron density values. The electron density values are collected from different sources based on the conditions specified in the function. - The returned list contains electron density values for both core and edge profiles, with NaN + The returned list contains electron density values for both core and edge profiles, with nan values filled in case of missing data. """ electron_density = [0] * (self.nrho + self.erho) @@ -823,7 +823,7 @@ def getne_profile(self): if len(_density) != self.nrho: logger.warning("core_profiles.profiles_1d[:].electrons.density could not be read.") logger.warning(f"Size mismatch: rho_tor_norm = {self.nrho}, electrons.density =" f"{len(_density)}") - _density = np.asarray([np.NaN] * self.nrho) + _density = np.asarray([np.nan] * self.nrho) for i in range(self.nrho): electron_density[i] = _density[i] if self.is_edge_profiles_present: @@ -831,14 +831,14 @@ def getne_profile(self): _density = self.core_profiles.profiles_1d[self.time_index_edge_profiles].electrons.density.value if len(_density) < 1: logger.warning("edge_profiles.profiles_1d[:].electrons.density could not be read.") - _density = np.asarray([np.NaN] * self.erho) + _density = np.asarray([np.nan] * self.erho) for i in range(self.erho): electron_density[self.nrho + i] = _density[i] else: _gset = self.edge_profiles.ggd[self.time_index_edge_profiles].electrons.density[self.gset].values if len() < 1: logger.warning("edge_profiles.ggd[:].electrons.density could not be read.") - _gset = np.asarray([np.NaN] * self.erho) + _gset = np.asarray([np.nan] * self.erho) for i in range(self.erho): electron_density[self.nrho + i] = _gset[i] return electron_density @@ -862,7 +862,7 @@ def gette_profile(self): logger.warning( f"Size mismatch: rho_tor_norm = {self.nrho}, electrons.temperature = " f"{len(_temperature)}" ) - _temperature = np.asarray([np.NaN] * self.nrho) + _temperature = np.asarray([np.nan] * self.nrho) for i in range(self.nrho): electron_temperature[i] = _temperature[i] * 1.0e-3 if self.is_edge_profiles_present: @@ -870,14 +870,14 @@ def gette_profile(self): _temperature = self.edge_profiles.profiles_1d[self.time_index_edge_profiles].electrons.temperature if len(_temperature) < 1: logger.warning("edge_profiles.profiles_1d[:].electrons.temperature could not be read.") - _temperature = np.asarray([np.NaN] * self.erho) + _temperature = np.asarray([np.nan] * self.erho) for i in range(self.erho): electron_temperature[self.nrho + i] = _temperature[i] * 1.0e-3 else: _gset = self.edge_profiles.ggd[self.time_index_edge_profiles].electrons.temperature[self.gset].values if len() < 1: logger.warning("edge_profiles.ggd[:].electrons.temperature could not be read.") - _gset = np.asarray([np.NaN] * self.erho) + _gset = np.asarray([np.nan] * self.erho) for i in range(self.erho): electron_temperature[self.nrho + i] = _gset[i] * 1.0e-3 return electron_temperature @@ -897,7 +897,7 @@ def getti_flag(self): if len(t_i_average) != self.nrho: logger.warning("core_profiles.profiles_1d[:].t_i_average could not be read.") logger.warning(f"Size mismatch: rho_tor_norm = {self.nrho}, t_i_average = " f"{len(t_i_average)}") - t_i_average = np.asarray([np.NaN] * self.nrho) + t_i_average = np.asarray([np.nan] * self.nrho) else: ti_flag = 1 ti_e_flag = 0 @@ -907,14 +907,14 @@ def getti_flag(self): t_i_average = self.edge_profiles.profiles_1d[self.time_index_edge_profiles].t_i_average.value if len() < 1: logger.warning("edge_profiles.profiles_1d[:].t_i_average could not be read.") - t_i_average = np.asarray([np.NaN] * self.erho) + t_i_average = np.asarray([np.nan] * self.erho) else: ti_e_flag = 1 else: t_i_average = self.edge_profiles.ggd[self.time_index_edge_profiles].t_i_average[self.gset].values if len() < 1: logger.warning("edge_profiles.ggd[:].t_i_average could not be read.") - t_i_average = np.asarray([np.NaN] * self.erho) + t_i_average = np.asarray([np.nan] * self.erho) else: ti_e_flag = 1 @@ -927,7 +927,7 @@ def getti_flag(self): logger.warning( f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].temperature = " f"{len(temp)}" ) - temp = np.asarray([np.NaN] * self.nrho) + temp = np.asarray([np.nan] * self.nrho) else: ti_flag = 2 if self.is_edge_profiles_present and ti_e_flag == 0: @@ -944,7 +944,7 @@ def getti_flag(self): logger.warning( f"edge_profiles.profiles_1d[:].ion[{jspecies}].temperature could not be read." ) - temperature = np.asarray([np.NaN] * self.erho) + temperature = np.asarray([np.nan] * self.erho) else: ti_e_flag = 2 else: @@ -954,7 +954,7 @@ def getti_flag(self): if len(temperature) < 1: if ti_e_flag != 1: logger.warning("edge_profiles.ggd[:].ion[:].temperature could not be read.") - temperature[self.gset].values = np.asarray([np.NaN] * self.erho) + temperature[self.gset].values = np.asarray([np.nan] * self.erho) else: ti_e_flag = 2 @@ -1032,7 +1032,7 @@ def get_ion_density(self): logger.warning( f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].density = " f"{len(density)}" ) - density = np.asarray([np.NaN] * self.nrho) + density = np.asarray([np.nan] * self.nrho) for i in range(self.nrho): ion_density[ispecies][i] = density[i] if self.is_edge_profiles_present: @@ -1055,7 +1055,7 @@ def get_ion_density(self): logger.warning( f"edge_profiles.profiles_1d[:].ion[{jspecies}].density could not be read." ) - _density = np.asarray([np.NaN] * self.erho) + _density = np.asarray([np.nan] * self.erho) for i in range(self.erho): ion_density[ispecies][self.nrho + i] = _density[i] else: @@ -1091,7 +1091,7 @@ def get_ion_density(self): ) if len() < 1: logger.warning(f"edge_profiles.ggd[:].ion[{jspecies}.density could not be read.") - _density = np.asarray([np.NaN] * self.erho) + _density = np.asarray([np.nan] * self.erho) for i in range(self.erho): ion_density[ispecies][self.nrho + i] = _density[i] else: @@ -1140,7 +1140,7 @@ def get_v_phi_profile(self): f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].velocity.toroidal = " f"{len(vphioid)}" ) - vphioid = np.asarray([np.NaN] * self.nrho) + vphioid = np.asarray([np.nan] * self.nrho) else: vphi_flag = 1 for i in range(self.nrho): @@ -1152,7 +1152,7 @@ def get_v_phi_profile(self): logger.warning( f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].velocity_tor = " f"{len(vphi)}" ) - vphi = np.asarray([np.NaN] * self.nrho) + vphi = np.asarray([np.nan] * self.nrho) else: if vphi_flag == 0: vphi_flag = 2 @@ -1346,7 +1346,7 @@ def get_vpol_profile(self): f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].velocity.poloidal =" f"{len(vpoloidal)}" ) - vpoloidal = np.asarray([np.NaN] * self.nrho) + vpoloidal = np.asarray([np.nan] * self.nrho) else: vpol_flag = 1 for i in range(self.nrho): @@ -1358,7 +1358,7 @@ def get_vpol_profile(self): logger.warning( f"Size mismatch: rho_tor_norm = {self.nrho}, ion[{ispecies}].velocity_pol = " f"{len(vpol)}" ) - vpol = np.asarray([np.NaN] * self.nrho) + vpol = np.asarray([np.nan] * self.nrho) else: if vpol_flag == 0: vpol_flag = 2 @@ -1920,7 +1920,7 @@ def get_waveform(self): # ) except Exception as e: logger.debug(f"{e}") - waveform["ti"]["central"] = [np.NaN] * len(self.time_array_core_profiles) + waveform["ti"]["central"] = [np.nan] * len(self.time_array_core_profiles) electrons_density = np.array([]) for i, _ in enumerate(self.time_array_core_profiles): @@ -2021,13 +2021,13 @@ def get_waveform(self): except Exception as e: logger.debug(f"{e}") waveform["n_species"][self.species[ispecies]]["density"]["central"] = [ - np.NaN + np.nan ] * self.common_time_length waveform["n_species"][self.species[ispecies]]["vpol"]["central"] = [ - np.NaN + np.nan ] * self.common_time_length waveform["n_species"][self.species[ispecies]]["vphi"]["central"] = [ - np.NaN + np.nan ] * self.common_time_length for itime in range(self.common_time_length): diff --git a/idstools/scripts/__init__.py b/idstools/scripts/__init__.py index a742386d..55d1dfd1 100644 --- a/idstools/scripts/__init__.py +++ b/idstools/scripts/__init__.py @@ -89,6 +89,11 @@ def plotcoresources(): _load_script_main("plotcoresources")() +def plotcoretransport(): + """Plot core transport data.""" + _load_script_main("plotcoretransport")() + + def plotkineticprofiles(): """Plot kinetic profiles.""" _load_script_main("plotkineticprofiles")() @@ -284,6 +289,7 @@ def md_summary(): "plotrotation", "plotequicomp", "plotcoresources", + "plotcoretransport", "plotkineticprofiles", "plothcdwaves", "plotecray", diff --git a/idstools/scripts/bin/dblist b/idstools/scripts/bin/dblist index 57d9e72a..d30a9e75 100644 --- a/idstools/scripts/bin/dblist +++ b/idstools/scripts/bin/dblist @@ -234,7 +234,15 @@ def print_times(dbs, args, print_times=False, pulse_number=None, run_number=None print(TAB * 4 + " Run: " + extended(str(run), RUN_STR_LEN)) available_ids_and_times = get_available_ids_and_times(connection) for idsname, times in available_ids_and_times: - if len(times) == 1 and np.isnan(times[0]): + if times is None: + print( + TAB * 5 + + extended(idsname, IDSNAME_STR_LEN) + + ": " + + extended("N/A", SLICENUM_STR_LEN) + + " slices ( time unavailable )" + ) + elif len(times) == 1 and np.isnan(times[0]): print( TAB * 5 + extended(idsname, IDSNAME_STR_LEN) @@ -311,7 +319,15 @@ def print_times_with_folder(dbs, print_times=False, pulse_number=None, run_numbe continue available_ids_and_times = get_available_ids_and_times(connection) for idsname, times in available_ids_and_times: - if times is not None: + if times is None: + print( + TAB * 15 + + extended(idsname, IDSNAME_STR_LEN) + + ": " + + extended("N/A", SLICENUM_STR_LEN) + + " slices ( time unavailable )" + ) + else: if len(times) == 1 and np.isnan(times[0]): print( TAB * 15 diff --git a/idstools/scripts/bin/dbperf b/idstools/scripts/bin/dbperf index 32849cb7..da5d8265 100644 --- a/idstools/scripts/bin/dbperf +++ b/idstools/scripts/bin/dbperf @@ -195,13 +195,39 @@ if __name__ == "__main__": pulse = entry[0] run = entry[1] + _dummy_hdfde = imas.DBEntry( + imas.ids_defs.HDF5_BACKEND, args.database, pulse, run, args.user, data_version=args.version + ) + _dummy_hdfde.open() + availables_ids = get_available_ids_and_occurrences(_dummy_hdfde) + for idsname, occurrence in availables_ids: + try: + _dummy = _dummy_hdfde.get(idsname, lazy=True, autoconvert=False) + src_dd_version = _dummy.ids_properties.version_put.data_dictionary.value + break + except Exception as e: + pass + _dummy_hdfde.close() + mdsde = imas.DBEntry( - imas.ids_defs.MDSPLUS_BACKEND, args.database, pulse, run, args.user, data_version=args.version + imas.ids_defs.MDSPLUS_BACKEND, + args.database, + pulse, + run, + args.user, + data_version=args.version, + dd_version=src_dd_version, ) mdsde.open() hdfde = imas.DBEntry( - imas.ids_defs.HDF5_BACKEND, args.database, pulse, run, args.user, data_version=args.version + imas.ids_defs.HDF5_BACKEND, + args.database, + pulse, + run, + args.user, + data_version=args.version, + dd_version=src_dd_version, ) hdfde.open() @@ -253,7 +279,6 @@ if __name__ == "__main__": mdsde.close() hdfde.close() - if args.verbose: print("MDSPLUS database closed") print("HDF5 database closed") diff --git a/idstools/scripts/bin/idsdiff b/idstools/scripts/bin/idsdiff index a7cb192b..6cec2dcb 100644 --- a/idstools/scripts/bin/idsdiff +++ b/idstools/scripts/bin/idsdiff @@ -10,6 +10,7 @@ import sys import time from io import BytesIO +from idstools.view.common import PlotCanvas import matplotlib.pyplot as plt import numpy as np import rich @@ -43,8 +44,6 @@ slicing_methods = { def view_plot(ax, field, coordinate, field_name="", coordinate_name="", field_unit="", coordinate_unit="", **kwargs): - from idstools.view.common import PlotCanvas - if not isinstance(field, (imas.ids_primitive.IDSNumericArray, np.ndarray)): print("Not a numeric array, Please select ids path") return @@ -506,7 +505,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -702,8 +704,6 @@ if __name__ == "__main__": ] = compare_data if args.plot and not args.html: - from idstools.view.common import PlotCanvas - canvas = PlotCanvas(1, 1) canvas.update_style(args.rc) ax = canvas.add_axes(title="", xlabel="", row=0, col=0) diff --git a/idstools/scripts/bin/idslist b/idstools/scripts/bin/idslist index 63578322..c5cf7456 100644 --- a/idstools/scripts/bin/idslist +++ b/idstools/scripts/bin/idslist @@ -184,8 +184,8 @@ if __name__ == "__main__": type = time_array elif time_array is None: value = f"{not_applicable}" - type = "unknown" - table.add_row(ids_name, value, Pretty(type)) + type = not_applicable + table.add_row(ids_name, value, type if isinstance(type, str) else Pretty(type)) if args.fullarray is True: with np.printoptions(threshold=sys.maxsize, linewidth=1024, precision=4): console.print(table) diff --git a/idstools/scripts/bin/idsprint b/idstools/scripts/bin/idsprint index f3a05934..e85df44e 100644 --- a/idstools/scripts/bin/idsprint +++ b/idstools/scripts/bin/idsprint @@ -3,16 +3,18 @@ import argparse import logging import os +import re +import shutil import sys try: import imaspy as imas except ImportError: import imas -from itertools import tee import numpy as np import rich +from rich.panel import Panel from rich.columns import Columns from rich.console import Console from rich.pretty import Pretty @@ -23,6 +25,7 @@ from rich_argparse import RichHelpFormatter from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, get_title, @@ -30,10 +33,161 @@ from idstools.utils.clihelper import ( ) from idstools.utils.idshelper import get_available_ids_and_times, parse_uri, partial_get from idstools.utils.idslogger import setup_logger -from idstools.view.common import PlotCanvas logger = setup_logger("module", stdout_level=logging.INFO) +TERMINAL_PLOT_WIDTH = 120 +TERMINAL_PLOT_HEIGHT = 30 +TERMINAL_PLOT_TICKS = 5 +TERMINAL_PLOT_SCIENTIFIC_THRESHOLD = 1e4 + + +def _node_value(node): + if isinstance(node, imas.ids_primitive.IDSPrimitive): + if not node.has_value: + return None + return node.value + return node + + +def _coordinate_1d_array(coordinate, size): + coordinate = _node_value(coordinate) + if coordinate is None or isinstance(coordinate, int): + return np.arange(size, dtype=float) + try: + coordinate = np.asarray(coordinate) + except Exception: + return np.arange(size, dtype=float) + if coordinate.ndim != 1 or coordinate.size != size or not np.issubdtype(coordinate.dtype, np.number): + return np.arange(size, dtype=float) + return coordinate.astype(float, copy=False) + + +def _scientific_ticks(values, count=TERMINAL_PLOT_TICKS): + values = np.asarray(values, dtype=float) + values = values[np.isfinite(values)] + if values.size == 0: + return [], [] + lower = float(values.min()) + upper = float(values.max()) + if lower == upper: + ticks = [lower] + else: + ticks = np.linspace(lower, upper, min(count, values.size)).tolist() + use_scientific = np.max(np.abs(values)) >= TERMINAL_PLOT_SCIENTIFIC_THRESHOLD + labels = [f"{tick:.2e}" if use_scientific else f"{tick:.4g}" for tick in ticks] + return ticks, labels + + +def _series_labels(field_name, count): + """Expand the first path slice into one label per plotted series.""" + match = re.search(r"[\[(]([^\])]*:[^\])]*?)[\])]", field_name) + if match: + parts = match.group(1).split(":") + if len(parts) <= 3: + try: + start = int(parts[0]) if parts[0] else None + step = int(parts[2]) if len(parts) == 3 and parts[2] else 1 + except ValueError: + start = None + step = None + if step: + if start is None: + start = count - 1 if step < 0 else 0 + labels = [] + for series_index in range(count): + selected_index = start + series_index * step + labels.append(f"{field_name[: match.start()]}[{selected_index}]{field_name[match.end() :]}") + return labels + return [f"{field_name} [{series_index}]" for series_index in range(count)] + + +def print_terminal_plot( + field, + coordinate, + field_name="", + coordinate_name="Index", + field_unit="", + coordinate_unit="", +): + raw = _node_value(field) + try: + data = np.asarray(raw).astype(float) + except (ValueError, TypeError): + data = None + + if data is None or data.size == 0 or data.ndim != 1: + logger.warning("Terminal plot supports only 1D numeric arrays") + return False + + # Extract units and coordinate from IDS metadata if not provided + if isinstance(field, imas.ids_primitive.IDSPrimitive): + if not field_unit: + field_unit = str(field.metadata.units or "") + if coordinate is None and isinstance(field, imas.ids_primitive.IDSNumericArray): + for candidate in field.coordinates: + if isinstance(candidate, imas.ids_primitive.IDSPrimitive) and candidate.has_value: + coordinate = candidate + coordinate_name = candidate.metadata.name or candidate.metadata.path or coordinate_name + if not coordinate_unit: + coordinate_unit = str(candidate.metadata.units or "") + break + if isinstance(coordinate, imas.ids_primitive.IDSPrimitive) and not coordinate_unit: + coordinate_unit = str(coordinate.metadata.units or "") + + import plotext as ptx + from rich.align import Align + from rich.text import Text + + console = Console() + + y_label = f"{field_name} [{field_unit}]" if field_unit else field_name + x_label = f"{coordinate_name} [{coordinate_unit}]" if coordinate_unit else (coordinate_name or "Index") + + xdata = _coordinate_1d_array(coordinate, data.size) + finite = np.isfinite(data) & np.isfinite(xdata) + if not np.any(finite): + logger.warning("Terminal plot supports only arrays with finite numeric values") + return False + data, xdata = data[finite], xdata[finite] + + ptx.clear_figure() + ptx.canvas_color("default") + ptx.axes_color("default") + ptx.ticks_color("default") + terminal_width, terminal_height = ptx.terminal_size() + plot_width = min(TERMINAL_PLOT_WIDTH, max(40, terminal_width - 4)) + plot_height = min(TERMINAL_PLOT_HEIGHT, max(12, terminal_height - 8)) + ptx.plot_size(plot_width, plot_height) + xticks, xlabels = _scientific_ticks(xdata) + yticks, ylabels = _scientific_ticks(data) + ptx.xticks(xticks, xlabels) + ptx.yticks(yticks, ylabels) + ptx.plot(xdata.tolist(), data.tolist(), color="cyan", marker="braille") + + stats = Text(justify="center") + stats.append("min: ", style="bright_black") + stats.append(f"{data.min():.2e}", style="bold cyan") + stats.append(" max: ", style="bright_black") + stats.append(f"{data.max():.2e}", style="bold cyan") + stats.append(" mean: ", style="bright_black") + stats.append(f"{data.mean():.2e}", style="bold cyan") + console.print( + Align.center(Panel(Align.center(stats), border_style="bright_blue", padding=(0, 1), width=plot_width)) + ) + + plot = Text.from_ansi(ptx.build()) + console.print(Align.center(plot)) + + axis_labels = Table.grid(padding=(0, 1)) + axis_labels.add_column(justify="right", style="bright_black", no_wrap=True) + axis_labels.add_column(style="cyan") + axis_labels.add_row("x:", x_label) + if y_label: + axis_labels.add_row("y:", y_label) + console.print(Align.center(Panel(axis_labels, border_style="bright_black", padding=(0, 1), width=plot_width))) + return True + def view_plot( ax, @@ -98,11 +252,13 @@ def view_plot( ax.plot(coordinate, data, label=field_name, marker="o", color="red") elif len(data.shape) == 2: + series_labels = _series_labels(field_name, data.shape[1]) + legend_indices = set(range(data.shape[1])) + if data.shape[1] > 20: + legend_indices = set(range(5)) | set(range(data.shape[1] - 5, data.shape[1])) for j in range(data.shape[1]): - if j == 0: - ax.plot(coordinate, data[:, j], label=f"{field_name}") - else: - ax.plot(coordinate, data[:, j]) + label = series_labels[j] if j in legend_indices else "_nolegend_" + ax.plot(coordinate, data[:, j], label=label) else: ax.plot(coordinate, data, label=field_name) @@ -117,12 +273,25 @@ def view_plot( def print_tree(structure, hide_empty_nodes, compact, full_array, depth=None): + console = Console() if full_array: with np.printoptions(threshold=sys.maxsize, linewidth=1024, precision=4): - rich.print(_make_tree(structure, hide_empty_nodes, compact, depth=depth)) + console.print( + Panel( + _make_tree(structure, hide_empty_nodes, compact, depth=depth), + border_style="bright_blue", + padding=(0, 1), + ) + ) else: with np.printoptions(threshold=5, linewidth=1024, precision=4): - rich.print(_make_tree(structure, hide_empty_nodes, compact, depth=depth)) + console.print( + Panel( + _make_tree(structure, hide_empty_nodes, compact, depth=depth), + border_style="bright_blue", + padding=(0, 1), + ) + ) def _make_tree(structure, hide_empty_nodes, compact, *, tree=None, depth=None, current_depth=0): @@ -173,7 +342,19 @@ def _make_tree(structure, hide_empty_nodes, compact, *, tree=None, depth=None, c group = Columns([txt]) tree.add(group) return tree - elif isinstance(structure, np.float64) or isinstance(structure, imas.ids_primitive.IDSPrimitive): + elif isinstance(structure, imas.ids_primitive.IDSPrimitive): + tree = Tree(f"[magenta]{type(structure)}") + if not structure.has_value: + tree.add(f"[grey62]{structure.metadata.name}[/]") + else: + value = Pretty(structure.value) + txt = f"[bright_yellow]{structure.metadata.name}[/]:" + group = Columns([txt, value]) + if compact: + group = Columns([f"[bright_yellow]{structure.metadata.name}[/]"]) + tree.add(group) + return tree + elif isinstance(structure, np.float64): return structure else: tree = Tree(f"[magenta]{structure.metadata.name}") @@ -212,9 +393,9 @@ def _make_tree(structure, hide_empty_nodes, compact, *, tree=None, depth=None, c if isinstance(child, imas.ids_structure.IDSStructure): txt = f"[magenta]{child._path}[/]" # check if structure is not empty - iterator, iterator_copy = tee(child.iter_nonempty_(accept_lazy=True)) + _iter = child.iter_nonempty_(accept_lazy=True) try: - next(iterator_copy) + next(_iter) except StopIteration: pass else: @@ -269,7 +450,7 @@ def _make_dict_tree(structure, hide_empty_nodes, compact, *, tree=None): for counter, ids_structure in enumerate(child): ntree = {} nlist.append(ntree) - nlist.append(_make_dict_tree(ids_structure, hide_empty_nodes, compact, tree=ntree)) + _make_dict_tree(ids_structure, hide_empty_nodes, compact, tree=ntree) return tree @@ -332,7 +513,7 @@ if __name__ == "__main__": "-p", "--plot", action="store_true", - help="plot 1d arrays from leaf nodes", + help="Plot arrays from leaf nodes with Matplotlib", ) parser.add_argument( "--coordinate", @@ -341,7 +522,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -350,6 +534,9 @@ if __name__ == "__main__": default=None, ) args = parser.parse_args() + if args.plot: + from idstools.view.common import PlotCanvas + time_slice = 0.0 if args.time: time_slice = float(args.time) @@ -375,7 +562,10 @@ if __name__ == "__main__": table.add_column("SLICES", style="green") table.add_column("TIME", style="green") for ids_name, time_array in available_ids_and_times: - if len(time_array) == 1 and np.isnan(time_array[0]): + if time_array is None: + value = "N/A" + type = "unknown" + elif len(time_array) == 1 and np.isnan(time_array[0]): value = f"{question_string}" type = "heterogeneous IDS" elif len(time_array) == 1 and time_array[0] == -np.inf: @@ -583,18 +773,25 @@ if __name__ == "__main__": canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) canvas.fig.suptitle(get_title(args, f"{ids_name}/{original_ids_path}", None)) - if args.save: - fname = get_file_name(args, os.path.basename(__file__), None) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__), None) else: print_tree( node, hide_empty_nodes=not args.show_empty, compact=args.compact, full_array=args.full, depth=args.depth ) + arr = _node_value(node) + try: + arr = np.asarray(arr) + except Exception: + arr = None + if arr is not None and arr.ndim == 1 and arr.size > 0 and np.issubdtype(arr.dtype, np.number): + print_terminal_plot( + node, + coordinate, + field_name=f"{ids_name}/{original_ids_path}", + coordinate_name=coordinate_name or "Index", + field_unit=node_unit, + coordinate_unit=coordinate_unit, + ) + connection.close() diff --git a/idstools/scripts/bin/idsquery b/idstools/scripts/bin/idsquery index e23047b9..9d346101 100644 --- a/idstools/scripts/bin/idsquery +++ b/idstools/scripts/bin/idsquery @@ -41,6 +41,7 @@ from idstools.utils.idshelper import ( get_ids_values, ) from idstools.utils.idslogger import setup_logger +from idstools.utils.utility_functions import add_query_to_uri logger = setup_logger("module", stdout_level=logging.INFO) @@ -108,7 +109,12 @@ Example: valpath = path[1 + len(idsname) :] paths_info.append((path, idsname, valpath.replace("(", "[").replace(")", "]").replace("/", "."))) - ids_values = get_ids_values(args.uri, paths_info, dd_update=args.dd_update, verbose=args.verbose) + ids_values = get_ids_values( + add_query_to_uri(args.uri, backend="uda", query="cache_mode=none"), + paths_info, + dd_update=args.dd_update, + verbose=args.verbose, + ) if ids_values: for _path, _value in ids_values.items(): if _value is None: diff --git a/idstools/scripts/bin/idsresample b/idstools/scripts/bin/idsresample index d0c4be01..31de01f0 100644 --- a/idstools/scripts/bin/idsresample +++ b/idstools/scripts/bin/idsresample @@ -117,7 +117,7 @@ if __name__ == "__main__": except Exception as e: logger.error(f"Exception occurred detailed description : {e}") - tmp = imas.DBEntry("imas:memory?path=tmp/1/1", "w") + tmp = imas.DBEntry("imas:memory?path=tmp/1/1", "w", dd_version=src_dd_version) if tmp is None: print("Error creating temporary data-entry in memory!", file=sys.stderr) diff --git a/idstools/scripts/bin/plotcoresources b/idstools/scripts/bin/plotcoresources index c8f7a4b2..6388fc40 100644 --- a/idstools/scripts/bin/plotcoresources +++ b/idstools/scripts/bin/plotcoresources @@ -15,9 +15,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -35,7 +35,10 @@ if __name__ == "__main__": parser.add_argument("-t", "--time", help="time", required=False, type=float, default=-99.0) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -98,17 +101,7 @@ if __name__ == "__main__": ax_torque_waveform = canvas.add_axes(row=1, col=3) ret = core_source_view.view_torque_waveform(ax_torque_waveform, time_slice) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust(top=0.916, bottom=0.09, left=0.044, right=0.953, hspace=0.287, wspace=0.2) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Core sources", time_value)) - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_core_sources", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + show_plot(canvas, args, os.path.basename(__file__) + "_core_sources", time_value) diff --git a/idstools/scripts/bin/plotcoretransport b/idstools/scripts/bin/plotcoretransport index a66560b1..da0a4c96 100644 --- a/idstools/scripts/bin/plotcoretransport +++ b/idstools/scripts/bin/plotcoretransport @@ -22,9 +22,9 @@ from rich_argparse import RichHelpFormatter from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -64,7 +64,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -175,17 +178,8 @@ if __name__ == "__main__": model_index, logscale=args.logscale, ) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.suptitle(get_title(args, "Core transport", time_value)) + canvas.set_sup_title(get_database_path(args, time_value=time_value)) canvas.fig.subplots_adjust(top=0.9, bottom=0.094, left=0.035, right=0.948, hspace=0.417, wspace=0.117) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) canvas.remove_empty_axes() - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_core_transport", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__) + "_core_transport", time_value) diff --git a/idstools/scripts/bin/ploteccomposition b/idstools/scripts/bin/ploteccomposition index 26684b76..adf983df 100644 --- a/idstools/scripts/bin/ploteccomposition +++ b/idstools/scripts/bin/ploteccomposition @@ -14,9 +14,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -52,7 +52,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -113,20 +116,11 @@ if __name__ == "__main__": waves_view.plot_ecrh_waveform(ax3, time_slice) waves_view.plot_e_c_c_d_waveform(ax4, time_slice) - canvas.set_text(text=f"{get_database_path(args, time_value)}") - canvas.fig.suptitle(get_title(args, "EC Composition", time_value)) + canvas.set_sup_title(get_database_path(args, time_value=time_value)) canvas.fig.subplots_adjust(top=0.941, bottom=0.122, left=0.052, right=0.925, hspace=0.2, wspace=0.2) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_EC_Composition", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__) + "_EC_Composition", time_value) conn.close() # def StrToList(string): # if string[0] == "[": diff --git a/idstools/scripts/bin/plotecray b/idstools/scripts/bin/plotecray index 7dfeb4a3..c858fa90 100644 --- a/idstools/scripts/bin/plotecray +++ b/idstools/scripts/bin/plotecray @@ -19,6 +19,7 @@ from idstools.compute.waves import WavesCompute from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, get_title, @@ -67,7 +68,10 @@ if __name__ == "__main__": parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -141,9 +145,8 @@ if __name__ == "__main__": exit(1) # Search for adequate time slice for display - time_array = ids_waves.time - ntime = len(ids_waves.time) - time_slice, time_value = get_nearest_time(time_array, args.time) + time_index_waves, time_value = get_nearest_time(ids_waves.time, args.time) + time_index_equilibrium, _ = get_nearest_time(ids_equilibrium.time, time_value) if len(ids_waves.code.name) > 0: logger.info(f"Code name = {ids_waves.code.name.upper()}") @@ -159,13 +162,18 @@ if __name__ == "__main__": wave_view = WavesView(ids_waves) wave_compute = WavesCompute(ids_waves) - beam_tracing_dict = wave_compute.get_beam_tracing(time_slice) + beam_tracing_dict = wave_compute.get_beam_tracing(time_index_waves) logger.info( f"There are {beam_tracing_dict['active_beams_count']} active beam(s)" f"and each beam has {beam_tracing_dict['max_total_beams']} ray(s)" ) - ecstra_view.plot_poloidal_view(ax_polview, coherent_wave_index=0, time_slice=time_slice) + ecstra_view.plot_poloidal_view( + ax_polview, + coherent_wave_index=0, + time_slice=time_index_waves, + equilibrium_time_slice=time_index_equilibrium, + ) if args.md is True: args.md = ["wall", "pf_active"] @@ -242,32 +250,22 @@ if __name__ == "__main__": wave_view.plot_pol_view_traces( ax_polview, - time_slice, + time_index_waves, color=color, style=style, ) - equi_view.plot_topplotequilibrium(ax_topview, time_slice) - wave_view.plot_top_view_traces(ax_topview, time_slice, color=color, style=style, label=label_code) - - wave_view.plot_electron_power(ax_powview, time_slice, color=color, style=style) - wave_view.plot_power_flow_normal(ax_powparview, time_slice, color=color, style=style) + equi_view.plot_topplotequilibrium(ax_topview, time_index_equilibrium) + wave_view.plot_top_view_traces(ax_topview, time_index_waves, color=color, style=style, label=label_code) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") + wave_view.plot_electron_power(ax_powview, time_index_waves, color=color, style=style) + wave_view.plot_power_flow_normal(ax_powparview, time_index_waves, color=color, style=style) canvas.fig.subplots_adjust(top=0.95, bottom=0.097, left=0, right=0.948, hspace=0.2, wspace=0.108) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "EC rays", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_EC_rays", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_EC_rays", time_value) # def StrToList(string): # if string[0] == "[": diff --git a/idstools/scripts/bin/plotecstrayradiation b/idstools/scripts/bin/plotecstrayradiation index 07080d87..7604bfd3 100644 --- a/idstools/scripts/bin/plotecstrayradiation +++ b/idstools/scripts/bin/plotecstrayradiation @@ -4,6 +4,7 @@ import argparse import logging import os +import sys from rich_argparse import RichHelpFormatter @@ -17,11 +18,11 @@ from idstools.input_processing import ( read_wall, ) from idstools.utils.clihelper import ( + dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, - dbentry_parser, ) from idstools.utils.idslogger import setup_logger from idstools.view.common import PlotCanvas @@ -32,6 +33,15 @@ from idstools.view.polygon import PolygonView from idstools.view.waves import WavesView logger = setup_logger("module", stdout_level=logging.INFO) + + +def _first_existing_path(*paths): + for path in paths: + if os.path.exists(path): + return path + return paths[-1] + + if __name__ == "__main__": parser = argparse.ArgumentParser( description="---- Shows electron cyclotron stray radiation information by showing different plots", @@ -46,7 +56,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -61,20 +74,28 @@ if __name__ == "__main__": time_index_waves = 0 current_file_path = os.path.dirname(os.path.abspath(__file__)) + source_tree_root = os.path.abspath(os.path.join(current_file_path, "../../..")) - scenario_file = os.path.join(current_file_path, "../resources/input/scenario.yaml") - wallfile = os.path.join(current_file_path, "../resources/input/wall2d.txt") - filelaunchers = os.path.join(current_file_path, "../resources/input/ec_waveforms.yaml") - path_result = os.path.join(current_file_path, "../resources/results/") - - if not os.path.exists(scenario_file): - scenario_file = os.path.join(current_file_path, "input/scenario.yaml") - if not os.path.exists(wallfile): - wallfile = os.path.join(current_file_path, "input/wall2d.txt") - if not os.path.exists(filelaunchers): - filelaunchers = os.path.join(current_file_path, "input/ec_waveforms.yaml") - if not os.path.exists(path_result): - path_result = os.path.join(current_file_path, "results/") + scenario_file = _first_existing_path( + os.path.join(source_tree_root, "resources/input/scenario.yaml"), + os.path.join(sys.prefix, "bin/input/scenario.yaml"), + os.path.join(current_file_path, "input/scenario.yaml"), + ) + wallfile = _first_existing_path( + os.path.join(source_tree_root, "resources/input/wall2d.txt"), + os.path.join(sys.prefix, "bin/input/wall2d.txt"), + os.path.join(current_file_path, "input/wall2d.txt"), + ) + filelaunchers = _first_existing_path( + os.path.join(source_tree_root, "resources/input/ec_waveforms.yaml"), + os.path.join(sys.prefix, "bin/input/ec_waveforms.yaml"), + os.path.join(current_file_path, "input/ec_waveforms.yaml"), + ) + path_result = _first_existing_path( + os.path.join(source_tree_root, "resources/results"), + os.path.join(sys.prefix, "bin/results"), + os.path.join(current_file_path, "results"), + ) wall2d = read_wall(wallfile) @@ -166,7 +187,6 @@ if __name__ == "__main__": ax_polygon, wall2d, beam_wall, coherent_wave_index, time_index_waves, time_index_waves ) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value_equilibrium)}") canvas.fig.subplots_adjust( top=0.88, bottom=0.11, @@ -175,14 +195,6 @@ if __name__ == "__main__": hspace=0.458, wspace=0.234, ) - canvas.fig.suptitle(get_title(args, "EC Stray Radiation", time_value_equilibrium)) + canvas.set_sup_title(get_database_path(args, time_value=time_value_equilibrium)) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - if args.save: - fname = get_file_name(args, f"{os.path.basename(__file__)}_Equilibrium", time_value_equilibrium) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, f"{os.path.basename(__file__)}_Equilibrium", time_value_equilibrium) diff --git a/idstools/scripts/bin/plotedgeprofiles b/idstools/scripts/bin/plotedgeprofiles index 4d484430..f554f9de 100644 --- a/idstools/scripts/bin/plotedgeprofiles +++ b/idstools/scripts/bin/plotedgeprofiles @@ -18,9 +18,9 @@ from idstools.compute.wall import WallCompute from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -59,7 +59,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -186,19 +189,9 @@ if __name__ == "__main__": edge_profiles_view.view_equatorial_plane_and_diverter_density(ax4, time_slice, logscale=args.logscale) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust(top=0.93, bottom=0.067, left=0.026, right=0.953, hspace=0.287, wspace=0.12) - canvas.fig.suptitle(get_title(args, "Edge Profiles", time_value)) + canvas.set_sup_title(get_database_path(args, time_value=time_value)) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_edge_profiles", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__) + "_edge_profiles", time_value) diff --git a/idstools/scripts/bin/plotequicomp b/idstools/scripts/bin/plotequicomp index f6716148..0f114532 100644 --- a/idstools/scripts/bin/plotequicomp +++ b/idstools/scripts/bin/plotequicomp @@ -10,16 +10,21 @@ try: except ImportError: import imas import numpy as np + +# PlotCanvas configures a --rc backend before Matplotlib is imported. +from idstools.view.common import PROVENANCE_TITLE_STYLE, PlotCanvas from matplotlib.animation import FuncAnimation from matplotlib.widgets import Slider from rich_argparse import RichHelpFormatter from idstools.compute.common import get_nearest_time from idstools.utils.clihelper import ( + show_plot, rcparam_parser, ) from idstools.utils.idslogger import setup_logger -from idstools.view.common import PlotCanvas +from idstools.utils.utility_functions import add_query_to_uri +from idstools.view.common import PROVENANCE_TITLE_STYLE, PlotCanvas from idstools.view.equilibrium import EquilibriumView from idstools.view.wall import WallView @@ -57,7 +62,10 @@ if __name__ == "__main__": parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -67,6 +75,11 @@ if __name__ == "__main__": ) args = parser.parse_args() + def entry_title(uri, time_value): + if args.no_provenance: + return "" + return f"{uri}\n#time:{time_value:.2f}" + log_level = getattr(logging, args.log_level) logger = setup_logger("module", stdout_level=log_level) equilibrium1 = None @@ -75,7 +88,10 @@ if __name__ == "__main__": wall2 = None if len(args.uri) >= 1: - connection1 = imas.DBEntry(args.uri[0], "r") + connection1 = imas.DBEntry( + add_query_to_uri(args.uri[0], backend="uda", query="cache_mode=none"), + "r", + ) if args.dd_update: equilibrium1 = connection1.get("equilibrium", autoconvert=False) equilibrium1 = imas.convert_ids(equilibrium1, connection1.factory.version) @@ -99,7 +115,10 @@ if __name__ == "__main__": logger.critical("----> Could not open first data entry. Aborted.") exit(1) if len(args.uri) == 2: - connection2 = imas.DBEntry(args.uri[1], "r") + connection2 = imas.DBEntry( + add_query_to_uri(args.uri[1], backend="uda", query="cache_mode=none"), + "r", + ) if args.dd_update: equilibrium2 = connection2.get("equilibrium", autoconvert=False) equilibrium2 = imas.convert_ids(equilibrium2, connection2.factory.version) @@ -126,7 +145,7 @@ if __name__ == "__main__": canvas = PlotCanvas(1, 4) canvas.update_style(args.rc) - canvas.fig.subplots_adjust(top=0.90, bottom=0.12, left=0.044, right=0.946, hspace=0.216, wspace=0.240) + canvas.fig.subplots_adjust(top=0.90, bottom=0.16, left=0.044, right=0.946, hspace=0.216, wspace=0.240) ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0) ax2 = canvas.add_axes(title="", xlabel="", row=0, col=1) @@ -137,7 +156,7 @@ if __name__ == "__main__": current_time_slice1 = [time_slice1] - ax_slider = canvas.fig.add_axes([0.15, 0.02, 0.65, 0.03]) + ax_slider = canvas.fig.add_axes((0.15, 0.055, 0.65, 0.025)) time_slider = Slider( ax_slider, @@ -153,7 +172,7 @@ if __name__ == "__main__": time_text = canvas.fig.text( 0.82, - 0.02, + 0.055, f"Index: {time_slice1} | Time: {time_array1[time_slice1]:.2f} s", ha="left", va="bottom", @@ -212,34 +231,35 @@ if __name__ == "__main__": text_artists["title"] = None if len(args.uri) == 2: - y = 0.99 text_artists["text1"] = canvas.fig.text( 0.5, - y, - args.uri[0] + f"(time {time_val1:.2f})", + 0.985, + entry_title(args.uri[0], time_val1), ha="center", va="top", color="tab:green", transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) text_artists["text2"] = canvas.fig.text( 0.5, - y - 0.03, - args.uri[1] + f"(time {time_val2:.2f})", + 0.94, + entry_title(args.uri[1], time_val2), ha="center", va="top", color="tab:blue", transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) else: text_artists["title"] = canvas.fig.text( 0.5, - 0.99, - f"{args.uri[0]} time {time_val1:.2f}", + 0.985, + entry_title(args.uri[0], time_val1), ha="center", va="top", - fontsize=12, transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) canvas.fig.canvas.draw() @@ -362,44 +382,37 @@ if __name__ == "__main__": # Set initial title using the text_artists dictionary to avoid overdrawing if len(args.uri) == 2: - y = 0.99 text_artists["text1"] = canvas.fig.text( 0.5, - y, - args.uri[0] + f"(time {time_value1:.2f})", + 0.985, + entry_title(args.uri[0], time_value1), ha="center", va="top", color="tab:green", transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) text_artists["text2"] = canvas.fig.text( 0.5, - y - 0.03, - args.uri[1] + f"(time {time_value2:.2f})", + 0.94, + entry_title(args.uri[1], time_value2), ha="center", va="top", color="tab:blue", transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) else: text_artists["title"] = canvas.fig.text( 0.5, - 0.99, - f"{args.uri[0]} time {time_value1:.2f}", + 0.985, + entry_title(args.uri[0], time_value1), ha="center", va="top", - fontsize=12, transform=canvas.fig.transFigure, + **PROVENANCE_TITLE_STYLE, ) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - if args.save: - fname = f"plotequicomp_time_{time_value1}.png" - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, fname=f"plotequicomp_time_{time_value1}.png") diff --git a/idstools/scripts/bin/plotequilibrium b/idstools/scripts/bin/plotequilibrium index f4a3f03e..1da27c78 100644 --- a/idstools/scripts/bin/plotequilibrium +++ b/idstools/scripts/bin/plotequilibrium @@ -1,11 +1,10 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -# python scripts/plotequilibrium -p 134174 -r 117 +# python scripts/plotequilibrium --uri --profiles # -md "imas:mdsplus?user=public;shot=116000;run=2;database=ITER_MD;version=3#wall" # "imas:mdsplus?user=public;shot=111001;run=102;database=ITER_MD;version=3#pf_active" import argparse -import copy import logging import os @@ -13,6 +12,9 @@ try: import imaspy as imas except ImportError: import imas +# PlotCanvas configures a --rc backend before mpl_toolkits imports Matplotlib. +from idstools.view.common import PlotCanvas +from mpl_toolkits.axes_grid1.inset_locator import inset_axes from rich_argparse import RichHelpFormatter from idstools.compute.common import get_nearest_time @@ -21,13 +23,12 @@ from idstools.database import DBMaster from idstools.machinedescription import get_md_data from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger -from idstools.view.common import PlotCanvas from idstools.view.domain.mdplot import plot_machine_description from idstools.view.equilibrium import EquilibriumView @@ -49,14 +50,19 @@ if __name__ == "__main__": ) parser.add_argument("-t", "--time", help="Time (default=middle)", type=float, default=-99.0) parser.add_argument( - "--rho", - help="Show rho overlay on the plot", + "--phi", + help=("Overlay contours of toroidal magnetic flux Φ(R,Z) [Wb] " "from equilibrium profiles_2d/phi"), + action="store_true", + ) + parser.add_argument( + "--no-overlay", + dest="no_overlay", + help="Hide equilibrium overlays", action="store_true", ) parser.add_argument( - "-p", - "--plots", - help="Plots available quantities along with equilibrium", + "--profiles", + help="Plot available 1D profiles and time traces alongside the equilibrium", action="store_true", ) parser.add_argument( @@ -75,13 +81,16 @@ if __name__ == "__main__": """, ) parser.add_argument( - "--show-labels", - help="Show labels", + "--debug", + help="Show diagnostic logging", action="store_true", ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -104,7 +113,7 @@ if __name__ == "__main__": splitted_ids_path_fragment = ids_path_fragment.split("/", 1) occurrence = int(splitted_ids_path_fragment[0]) - logger = setup_logger("module", stdout_level=logging.INFO) + logger = setup_logger("module", stdout_level=logging.DEBUG if args.debug else logging.INFO) connection = DBMaster.get_connection(args) if connection is None: @@ -125,9 +134,7 @@ if __name__ == "__main__": time_slice, time_value = get_nearest_time(ids_obj_equilibrium.time, args.time) view_object = EquilibriumView(ids_obj_equilibrium) - title = f"2D Equilibrium at time {time_value:.3f}" - database_text = "" - if args.plots: + if args.profiles: compute_obj = EquilibriumCompute(ids_obj_equilibrium) profiles_1d_quantities = compute_obj.get_profiles_1d_quantities(time_slice, ["pressure", "q", "beta_pol"]) p1dcounter = sum(1 for value in profiles_1d_quantities.values() if value.has_value) @@ -136,7 +143,6 @@ if __name__ == "__main__": time_slice, ["q_min.value", "q_95", "li_3", "beta_tor", "energy_mhd"] ) gcounter = sum(1 for value in global_quantities.values() if value["has_value"]) - total_plots = p1dcounter + gcounter if total_plots % 2 == 1: @@ -145,6 +151,7 @@ if __name__ == "__main__": col_size = int(total_plots / 2) col_size = col_size + 1 + canvas = PlotCanvas(2, col_size) ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0, rowspan=2) axes_list1 = [] @@ -163,7 +170,11 @@ if __name__ == "__main__": ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0) canvas.update_style(args.rc) - if args.md: + md_overlay_state = {"created": False, "provenance_lines": []} + + def plot_md_overlay(): + if not args.md or md_overlay_state["created"]: + return md_overlay_state["provenance_lines"] idses = "" mduris = [] for mduri in args.md: @@ -188,44 +199,104 @@ if __name__ == "__main__": ids_data = get_md_data(mduris, args.dd_update, idses=idses) else: ids_data = get_md_data(mduris, args.dd_update) - plot_machine_description(ax1, ids_data) + md_provenance = plot_machine_description(ax1, ids_data, main_uri=get_database_path(args).strip()) + if md_provenance: + md_overlay_state["provenance_lines"] = md_provenance.splitlines() + ax1.set_title("") + md_overlay_state["created"] = True + return True + + c_psi, c_rho = view_object.view_magnetic_poloidal_flux( + ax1, + time_slice, + plot_magnetic_axis=False, + plot_current_centre=False, + plot_boundary_data=False, + plot_phi=args.phi, + plot_annotations=False, + ) + legend = ax1.get_legend() + if legend is not None: + legend.set_visible(False) - c_psi, c_rho = view_object.view_magnetic_poloidal_flux(ax1, time_slice, plot_rho=args.rho) if c_psi: - cbar_psi = canvas.fig.colorbar(c_psi, ax=ax1, orientation="horizontal", pad=0.08, fraction=0.03) - cbar_psi.set_label(r"$\psi$ [Wb]") + cax_psi = inset_axes( + ax1, + width=0.24, + height="88%", + loc="lower left", + bbox_to_anchor=(-0.24, 0.05, 1, 1), + bbox_transform=ax1.transAxes, + borderpad=0, + ) + cbar_psi = canvas.fig.colorbar(c_psi, cax=cax_psi) + cbar_psi.set_label(r"$\psi$ [Wb]", fontsize=7, rotation=90, labelpad=4) + cbar_psi.ax.yaxis.set_label_position("left") + cbar_psi.ax.tick_params(labelsize=7) + if c_rho: - cbar_rho = canvas.fig.colorbar(c_rho, ax=ax1, orientation="horizontal", pad=0.08, fraction=0.03) - cbar_rho.set_label(r"$\rho$ [Wb]") - ax1.set_title(title) - - xmin, xmax = ax1.get_xlim() - ymin, ymax = ax1.get_ylim() - ax1.text( - xmax + 0.01 * abs(xmax), - ymin + 0.5 * abs(ymax - ymin), - f"{get_database_path(args, time_value=time_value)}\n{database_text}", - horizontalalignment="left", - verticalalignment="center", - rotation="vertical", - fontsize=7, - ) - if args.plots: + cax_rho = inset_axes( + ax1, + width=0.24, + height="88%", + loc="lower left", + bbox_to_anchor=(-0.58, 0.05, 1, 1), + bbox_transform=ax1.transAxes, + borderpad=0, + ) + cbar_rho = canvas.fig.colorbar(c_rho, cax=cax_rho) + cbar_rho.set_label(r"$\Phi$ [Wb]", fontsize=7, rotation=90, labelpad=4) + cbar_rho.ax.yaxis.set_label_position("left") + cbar_rho.ax.tick_params(labelsize=7) + + plot_md_overlay() + + if args.profiles: view_object.plot_profiles_1d_quantities(axes_list1, time_slice) view_object.plot_global_quantities(axes_list2, time_value) - canvas.fig.suptitle(get_title(args, "Equilibrium", time_value)) - canvas.fig.set_size_inches(14, 8) - canvas.fig.subplots_adjust(top=0.933, bottom=0.05, left=0.024, right=0.988, hspace=0.221, wspace=0.25) - canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - if args.save: - fname = get_file_name(args, f"{os.path.basename(__file__)}_Equilibrium", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) + if args.profiles: + canvas.fig.set_size_inches(10 + col_size * 1.6, 8) + canvas.fig.subplots_adjust(top=0.933, bottom=0.100, left=0.05, right=0.955, hspace=0.221, wspace=0.20) + else: + canvas.fig.set_size_inches(14, 8) + canvas.fig.subplots_adjust(top=0.933, bottom=0.100, left=0.05, right=0.955, hspace=0.221, wspace=0.25) + + def create_overlays(show_legend=True): + plot_md_overlay() + view_object.view_magnetic_poloidal_flux( + ax1, + time_slice, + plot_magnetic_axis=True, + plot_current_centre=True, + plot_boundary_data=True, + plot_boundary_outline=True, + plot_phi=False, + plot_annotations=not args.profiles, + plot_psi=False, + ) + legend = ax1.get_legend() + if legend is not None: + legend.set_visible(show_legend) + + if args.no_overlay: + legend = ax1.get_legend() + if legend is not None: + legend.set_visible(not args.profiles) else: - canvas.show() + create_overlays(show_legend=not args.profiles) + + provenance_parts = [] + if not args.no_provenance: + provenance_parts.append(get_database_path(args).strip()) + if md_overlay_state["provenance_lines"]: + provenance_parts.extend(md_overlay_state["provenance_lines"]) + title = " | ".join(provenance_parts) + if title and time_value is not None: + title += f"\n#time:{time_value:.3f}" + canvas.set_sup_title(title, fontsize=8, y=0.985) + + canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) + show_plot(canvas, args, f"{os.path.basename(__file__)}_Equilibrium", time_value) else: logger.warning("Can not produce plot, equilibrium/time is None") diff --git a/idstools/scripts/bin/plothcd b/idstools/scripts/bin/plothcd index 3dd6fd51..d9001d3f 100644 --- a/idstools/scripts/bin/plothcd +++ b/idstools/scripts/bin/plothcd @@ -14,9 +14,9 @@ from idstools.compute.common import get_nearest_time from idstools.compute.waves import WavesCompute from idstools.database import DBMaster from idstools.utils.clihelper import ( + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -28,6 +28,9 @@ logger = setup_logger("module", stdout_level=logging.INFO) def _show_waves_plots(connargs, args, hold=False, dd_update=False, rc=""): + connargs.no_provenance = args.get("no_provenance", False) + connargs.save = args["save"] + connargs.directory = args["directory"] conn = DBMaster.get_connection(connargs) if conn is None: logger.critical(f"data entry not found : {connargs}") @@ -67,7 +70,6 @@ def _show_waves_plots(connargs, args, hold=False, dd_update=False, rc=""): canvas = PlotCanvas(rows, 2) canvas.update_style(rc) # canvas.setStyle(style="retro") - canvas.set_sup_title(f"HCD Waves Plot {connargs.uri} Time : {time_value:.3f}") ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0) ax2 = canvas.add_axes(title="", xlabel="", row=0, col=1) @@ -96,8 +98,6 @@ def _show_waves_plots(connargs, args, hold=False, dd_update=False, rc=""): else: ax4.get_legend().remove() - canvas.set_text(text=f"{get_database_path(connargs, time_value=time_value)}") - canvas.fig.subplots_adjust( top=0.92, bottom=0.122, @@ -107,12 +107,8 @@ def _show_waves_plots(connargs, args, hold=False, dd_update=False, rc=""): wspace=0.13, ) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(connargs, "HCD Waves Plot", time_value)) - if args["save"]: - fname = get_file_name(connargs, "hcd_waves_plot", time_value) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(connargs, time_value=time_value)) + show_plot(canvas, connargs, "hcd_waves_plot", time_value) conn.close() @@ -145,6 +141,9 @@ def show_waves_plots(args, hold=False, dd_update=False, rc=""): def _show_distribution_plots(connargs, args, hold=False, dd_update=False, rc=""): + connargs.no_provenance = args.get("no_provenance", False) + connargs.save = args["save"] + connargs.directory = args["directory"] conn = DBMaster.get_connection(connargs) if conn is None: logger.critical(f"data entry not found : {connargs}") @@ -168,7 +167,6 @@ def _show_distribution_plots(connargs, args, hold=False, dd_update=False, rc="") canvas = PlotCanvas(3, 2) canvas.update_style(rc) # canvas.setStyle(style="retro") - canvas.set_sup_title(f"HCD Distributions Plot {connargs.uri} Time : {time_value:.3f}") if ntime == 1: logger.info("Only one time slice --> Power and CD waveforms not displayed") ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0) @@ -184,8 +182,6 @@ def _show_distribution_plots(connargs, args, hold=False, dd_update=False, rc="") distributions_view.plot_nbi_fus_power_and_cd_waveforms(ax4, time_slice) distributions_view.plot_cd_waveform(ax5, time_slice) - canvas.set_text(text=f"{get_database_path(connargs, time_value=time_value)}") - canvas.fig.subplots_adjust( top=0.92, bottom=0.122, @@ -195,12 +191,8 @@ def _show_distribution_plots(connargs, args, hold=False, dd_update=False, rc="") wspace=0.328, ) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(connargs, "HCD Distributions Plot", time_value)) - if args["save"]: - fname = get_file_name(connargs, os.path.basename(__file__) + "_Distributions_profile_time", time_value) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(connargs, time_value=time_value)) + show_plot(canvas, connargs, os.path.basename(__file__) + "_Distributions_profile_time", time_value) conn.close() @@ -280,7 +272,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -451,22 +446,26 @@ if __name__ == "__main__": ec_args["force_psi"] = None ec_args["hide_legend"] = False ec_args["save"] = args["save"] + ec_args["directory"] = args["directory"] show_waves_plots(ec_args, hold=holder["ech_hold"], dd_update=args["dd_update"], rc=args["rc"]) # NUCLEAR REACTIONS if hcd_sources["fus"]["run"]: fus_args = hcd_sources["fus"] fus_args["save"] = args["save"] + fus_args["directory"] = args["directory"] show_distribution_plots(fus_args, hold=holder["icrh_hold"], dd_update=args["dd_update"], rc=args["rc"]) # NBI if hcd_sources["nbi"]["run"]: nbi_args = hcd_sources["nbi"] nbi_args["save"] = args["save"] + nbi_args["directory"] = args["directory"] show_distribution_plots(nbi_args, hold=holder["nbi_hold"], dd_update=args["dd_update"], rc=args["rc"]) # ICRH if hcd_sources["icrh"]["run"]: icrh_args = hcd_sources["icrh"] icrh_args["save"] = args["save"] + icrh_args["directory"] = args["directory"] show_distribution_plots(icrh_args, hold=holder["fus_hold"], dd_update=args["dd_update"], rc=args["rc"]) diff --git a/idstools/scripts/bin/plothcddistributions b/idstools/scripts/bin/plothcddistributions index 1c41665b..a03d0c2e 100644 --- a/idstools/scripts/bin/plothcddistributions +++ b/idstools/scripts/bin/plothcddistributions @@ -14,9 +14,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -67,8 +67,6 @@ def show_plots(args): distributions_view.plot_nbi_fus_power_and_cd_waveforms(ax4, time_slice) distributions_view.plot_cd_waveform(ax5, time_slice) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust( top=0.92, bottom=0.122, @@ -78,17 +76,9 @@ def show_plots(args): wspace=0.328, ) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Distributions profile", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_Distributions_profile", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_Distributions_profile", time_value) conn.close() @@ -102,7 +92,10 @@ if __name__ == "__main__": parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( diff --git a/idstools/scripts/bin/plothcdwaves b/idstools/scripts/bin/plothcdwaves index f9c28758..462123e4 100644 --- a/idstools/scripts/bin/plothcdwaves +++ b/idstools/scripts/bin/plothcdwaves @@ -15,9 +15,9 @@ from idstools.compute.waves import WavesCompute from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -96,8 +96,6 @@ def show_plots(args): else: ax4.get_legend().remove() - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust( top=0.92, bottom=0.122, @@ -107,17 +105,9 @@ def show_plots(args): wspace=0.13, ) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "HCD Waves Plot", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_heating_profiles_time", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_heating_profiles_time", time_value) conn.close() @@ -137,7 +127,10 @@ if __name__ == "__main__": parser.add_argument("-l", "--hide_legend", help="remove the legend from graphs", action="store_true") parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( diff --git a/idstools/scripts/bin/plotkineticprofiles b/idstools/scripts/bin/plotkineticprofiles index dd3d556c..078cae1b 100644 --- a/idstools/scripts/bin/plotkineticprofiles +++ b/idstools/scripts/bin/plotkineticprofiles @@ -9,9 +9,9 @@ from rich_argparse import RichHelpFormatter from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -36,7 +36,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -80,21 +83,9 @@ if __name__ == "__main__": kp_view.view_density_profiles(ax7, logscale=args.logscale) # Density profiles kp_view.view_vphi_profile(ax8, logscale=args.logscale) # Vtol profiles - canvas.set_text(text=f"{get_database_path(args, time_value=kp_view.k_profiles.time_value_core_profiles)}") - canvas.fig.subplots_adjust(top=0.928, bottom=0.11, left=0.033, right=0.91, hspace=0.435, wspace=0.518) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Plasma kinetic profiles ", kp_view.k_profiles.time_value_core_profiles)) + canvas.set_sup_title(get_database_path(args, time_value=kp_view.k_profiles.time_value_core_profiles)) - if args.save: - fname = get_file_name( - args, os.path.basename(__file__) + "_Kinetic", kp_view.k_profiles.time_value_core_profiles - ) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__) + "_Kinetic", kp_view.k_profiles.time_value_core_profiles) diff --git a/idstools/scripts/bin/plotmachinedescription b/idstools/scripts/bin/plotmachinedescription index cc5c3d2b..cad8e9eb 100644 --- a/idstools/scripts/bin/plotmachinedescription +++ b/idstools/scripts/bin/plotmachinedescription @@ -9,7 +9,7 @@ import os from rich_argparse import RichHelpFormatter from idstools.machinedescription import get_md_data -from idstools.utils.clihelper import rcparam_parser +from idstools.utils.clihelper import show_plot, rcparam_parser from idstools.utils.idslogger import setup_logger from idstools.view.common import PlotCanvas from idstools.view.domain.mdplot import plot_machine_description @@ -52,7 +52,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -68,16 +71,9 @@ if __name__ == "__main__": mdcanvas = PlotCanvas(1, 1, figsize=(10, 10)) mdcanvas.update_style(args.rc) ax = mdcanvas.add_axes(title="", xlabel="R (m)", ylabel="Z (m)", row=0, col=0) - plot_machine_description(ax, ids_data) + md_provenance = plot_machine_description(ax, ids_data) + ax.set_title("") mdcanvas.fig.subplots_adjust(top=0.916, bottom=0.09, left=0.044, right=0.953, hspace=0.287, wspace=0.2) mdcanvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - mdcanvas.fig.suptitle("Machine Description") - if args.save: - fname = os.path.basename(__file__) + "_machine_description.png" - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - mdcanvas.save(fname) - else: - mdcanvas.show() + mdcanvas.set_sup_title("" if args.no_provenance else md_provenance, fontsize=8) + show_plot(mdcanvas, args, fname=os.path.basename(__file__) + "_machine_description.png") diff --git a/idstools/scripts/bin/plotneutron b/idstools/scripts/bin/plotneutron index f9f40ce5..efc29de4 100644 --- a/idstools/scripts/bin/plotneutron +++ b/idstools/scripts/bin/plotneutron @@ -17,9 +17,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -38,7 +38,10 @@ if __name__ == "__main__": parser.add_argument("-t", "--time", help="Time", required=False, type=float, default=-99.0) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -125,18 +128,8 @@ if __name__ == "__main__": distribution_sources_view.view_neutrons(ax, time_slice) distribution_sources_view.view_time(ax, time_value) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust(top=0.916, bottom=0.09, left=0.044, right=0.953, hspace=0.287, wspace=0.2) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Neutrons profiles", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_Neutrons", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_Neutrons", time_value) diff --git a/idstools/scripts/bin/plotpressure b/idstools/scripts/bin/plotpressure index c263cda2..6751a4ef 100644 --- a/idstools/scripts/bin/plotpressure +++ b/idstools/scripts/bin/plotpressure @@ -17,9 +17,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -40,7 +40,10 @@ if __name__ == "__main__": parser.add_argument("-t", "--time", type=float, help="Time", default=-99) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -69,11 +72,9 @@ if __name__ == "__main__": time_array = ids_core_profiles.time time_slice, time_value = get_nearest_time(time_array, args.time) - title = "Profiles displayed for t = " + "%.1f" % time_value + " s" canvas = PlotCanvas(3, 1) canvas.update_style(args.rc) - canvas.fig.suptitle(title) ax1 = canvas.add_axes(title="", xlabel="", row=0, col=0, colspan=1) ax2 = canvas.add_axes(title="", xlabel="", row=1, col=0, colspan=1) @@ -85,18 +86,8 @@ if __name__ == "__main__": coreprofiles_view.plot_ion_pressure_properties(ax2, time_slice) coreprofiles_view.plot_electron_pressure_properties(ax3, time_slice) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust(top=0.916, bottom=0.09, left=0.044, right=0.953, hspace=0.287, wspace=0.2) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Pressure", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_Pressure", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_Pressure", time_value) diff --git a/idstools/scripts/bin/plotrotation b/idstools/scripts/bin/plotrotation index 2dc5a291..43fcbd46 100644 --- a/idstools/scripts/bin/plotrotation +++ b/idstools/scripts/bin/plotrotation @@ -16,9 +16,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -36,7 +36,10 @@ if __name__ == "__main__": parser.add_argument("-t", "--time", help="Time", required=False, type=float, default=99.0) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -87,18 +90,8 @@ if __name__ == "__main__": ax1.sharex(ax3) ax2.sharex(ax4) - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas.fig.subplots_adjust(top=0.916, bottom=0.09, left=0.044, right=0.953, hspace=0.174, wspace=0.117) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) - canvas.fig.suptitle(get_title(args, "Kinetic profiles", time_value)) - - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_Kinetic_profiles", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) - else: - canvas.show() + canvas.set_sup_title(get_database_path(args, time_value=time_value)) + + show_plot(canvas, args, os.path.basename(__file__) + "_Kinetic_profiles", time_value) diff --git a/idstools/scripts/bin/plotscenario b/idstools/scripts/bin/plotscenario index ffc2e978..1b2a2ad2 100644 --- a/idstools/scripts/bin/plotscenario +++ b/idstools/scripts/bin/plotscenario @@ -16,9 +16,9 @@ from idstools.compute.common import get_nearest_time from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idshelper import get_available_ids_and_occurrences @@ -54,7 +54,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -205,11 +208,10 @@ if __name__ == "__main__": plotequilibrium = EquilibriumView(ids_equilibrium) plotequilibrium.plotequilibrium(ax5, time_slice) - title = get_title(args, "Scenario") if not args.no_profiles: - title = get_title(args, "Scenario", time_value) + title = get_database_path(args, time_value=time_value) else: - title = get_title(args, "Scenario") + title = get_database_path(args) if args.info: title += ( f"\nprovider={ids_summary.ids_properties.provider}, " @@ -218,23 +220,11 @@ if __name__ == "__main__": f"access_layer={ids_summary.ids_properties.version_put.access_layer}" ) - if not args.no_profiles: - canvas.set_text(text=f"{get_database_path(args, time_value=time_value)}") - else: - canvas.set_text(text=f"{get_database_path(args)}") canvas.set_sup_title(title) canvas.get_current_fig_manager().set_window_title(os.path.basename(__file__)) canvas.fig.subplots_adjust(top=0.914, bottom=0.099, left=0.042, right=0.9, hspace=0.113, wspace=0.43) - if args.save: - if not args.no_profiles: - fname = get_file_name(args, os.path.basename(__file__) + "_summary", time_value) - else: - fname = get_file_name(args, os.path.basename(__file__) + "_Scenario_shot") - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas.save(fname) + if not args.no_profiles: + show_plot(canvas, args, os.path.basename(__file__) + "_summary", time_value) else: - canvas.show() + show_plot(canvas, args, os.path.basename(__file__) + "_Scenario_shot") diff --git a/idstools/scripts/bin/plotspectrometry b/idstools/scripts/bin/plotspectrometry index 384e783d..45ff434e 100644 --- a/idstools/scripts/bin/plotspectrometry +++ b/idstools/scripts/bin/plotspectrometry @@ -17,9 +17,9 @@ from idstools.compute.spectrometer_visible import SpectrometerVisibleCompute from idstools.database import DBMaster from idstools.utils.clihelper import ( dbentry_parser, + show_plot, get_database_path, get_file_name, - get_title, rcparam_parser, ) from idstools.utils.idslogger import setup_logger @@ -42,7 +42,10 @@ if __name__ == "__main__": ) parser.add_argument( "--save", - help="Save figure at default location", + help=( + "Save the plot to a file instead of opening a window. " + "If no plot window is available, it is saved automatically." + ), action="store_true", ) parser.add_argument( @@ -97,19 +100,10 @@ if __name__ == "__main__": ax.get_legend().remove() canvas_radiance.fig.subplots_adjust(top=0.88, bottom=0.11, left=0.065, right=0.893, hspace=0.497, wspace=0.243) - canvas_radiance.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas_radiance.fig.suptitle(get_title(args, "Spectrum (Radiance) from spectrometer_visible", time_value)) + canvas_radiance.set_sup_title(get_database_path(args, time_value=time_value)) canvas_radiance.get_current_fig_manager().set_window_title(os.path.basename(__file__) + "-radiance") - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "-radiance", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas_radiance.save(fname) - else: - canvas_radiance.show(block=False) + show_plot(canvas_radiance, args, os.path.basename(__file__) + "-radiance", time_value, show_kwargs={"block": False}) canvas_intensity = PlotCanvas(rows, columns) canvas_intensity.update_style(args.rc) row_counter = 0 @@ -124,18 +118,11 @@ if __name__ == "__main__": if column_counter != 0: ax.get_legend().remove() - canvas_intensity.set_text(text=f"{get_database_path(args, time_value=time_value)}") - canvas_intensity.fig.suptitle(get_title(args, "Spectrum (Intensity) from spectrometer_visible", time_value)) + canvas_intensity.set_sup_title(get_database_path(args, time_value=time_value)) canvas_intensity.fig.subplots_adjust(top=0.88, bottom=0.113, left=0.033, right=0.891, hspace=0.497, wspace=0.18) canvas_intensity.get_current_fig_manager().set_window_title(os.path.basename(__file__) + "-intensity") - if args.save: - fname = get_file_name(args, os.path.basename(__file__) + "_intensity", time_value) - if args.directory: - if not os.path.exists(args.directory): - os.makedirs(args.directory) - fname = os.path.join(args.directory, fname) - canvas_intensity.save(fname) - else: - canvas_intensity.show(block=True) + show_plot( + canvas_intensity, args, os.path.basename(__file__) + "_intensity", time_value, show_kwargs={"block": True} + ) connection.close() diff --git a/idstools/utils/clihelper.py b/idstools/utils/clihelper.py index c7b4cc19..9df05215 100644 --- a/idstools/utils/clihelper.py +++ b/idstools/utils/clihelper.py @@ -1,7 +1,17 @@ import argparse +import logging import os import re -import socket +from datetime import datetime + +from idstools.utils.matplotlib_backend import ( + _configure_backend_from_cli_rc, + _is_jupyter, +) + +_configure_backend_from_cli_rc() + +logger = logging.getLogger("module") try: import imaspy as imas @@ -68,6 +78,11 @@ default="", help="Semicolon-separated rcParams string (e.g., 'lines.linewidth=2;axes.titlesize=14')", ) +rcparam_parser.add_argument( + "--no-provenance", + action="store_true", + help="Hide URI provenance information from plot titles", +) def get_backend_id(name): @@ -192,21 +207,46 @@ def get_title(imasargs, title="", time_value=None): def get_file_name(imasargs, title="", time_value=None): - _file_name = "" - if title: - _file_name += f"{title}_" - if "uri" in imasargs.__dict__ and imasargs.uri: - param = get_details_from_uri(imasargs.uri) - if param["pathPresent"]: - _file_name += f"PATH_{param['path'].replace('/', '_')}_" - else: - _file_name += f"PULSE_{param['pulse']}_RUN_{param['run']}_" - else: - _file_name += f"PULSE_{imasargs.pulse}_RUN_{imasargs.run}_" - if time_value: - _file_name += f"TIME_{time_value:.3f}" - _file_name += ".png" - return _file_name + tool_name = title or "plot" + timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + return f"{tool_name}_{timestamp}.png" + + +def _is_interactive_backend(canvas): + """Return True for GUI and notebook widget backends.""" + import matplotlib + + backend = matplotlib.get_backend().lower() + if backend in ("widget", "ipympl", "module://ipympl.backend_nbagg"): + return True + return getattr(canvas.fig.canvas, "required_interactive_framework", None) is not None + + +def show_plot(canvas, imasargs, title="", time_value=None, fname=None, show_kwargs=None): + """Display interactive plots and save non-interactive plots automatically.""" + noninteractive = not _is_interactive_backend(canvas) + if not imasargs.save and (_is_jupyter() or not noninteractive): + canvas.show(**(show_kwargs or {})) + return + + automatic_save = noninteractive and not imasargs.save + if fname is None: + fname = get_file_name(imasargs, title, time_value) + if automatic_save: + extension = canvas.fig.canvas.get_default_filetype() + fname = f"{os.path.splitext(fname)[0]}.{extension}" + if imasargs.directory: + os.makedirs(imasargs.directory, exist_ok=True) + fname = os.path.join(imasargs.directory, fname) + if automatic_save: + import matplotlib + + logger.info( + "Non-interactive Matplotlib backend '%s' detected; saving figure to %s", + matplotlib.get_backend(), + fname, + ) + canvas.save(fname) def get_database_path(imasargs, time_value=None) -> str: @@ -219,6 +259,9 @@ def get_database_path(imasargs, time_value=None) -> str: Returns: the absolute path of the database. """ + if getattr(imasargs, "no_provenance", False): + return "" + pulse_info = "" database_absolute_path = "" if "uri" in imasargs.__dict__ and imasargs.uri: @@ -241,11 +284,10 @@ def get_database_path(imasargs, time_value=None) -> str: database_absolute_path = database_absolute_path[:-2] time_string = "" if time_value: - time_string = f"time:{time_value:.3f})" - hostdir = f"{socket.gethostname()}:{database_absolute_path} " + time_string = f"time:{time_value:.3f}" + hostdir = f"{database_absolute_path} " if pulse_info: hostdir += f"({pulse_info})" if time_string: hostdir += f"#{time_string}" - # return hostdir diff --git a/idstools/utils/idshelper.py b/idstools/utils/idshelper.py index aed61e8a..e2c41956 100644 --- a/idstools/utils/idshelper.py +++ b/idstools/utils/idshelper.py @@ -439,7 +439,7 @@ def get_available_ids_and_times(db_entry_object) -> list: if homogeneous_time == imas.ids_defs.IDS_TIME_MODE_UNKNOWN: time_array = [] if homogeneous_time == imas.ids_defs.IDS_TIME_MODE_HETEROGENEOUS: - time_array = [np.NaN] + time_array = [np.nan] if homogeneous_time == imas.ids_defs.IDS_TIME_MODE_HOMOGENEOUS: if getattr(ids_object, "time", None): time_array = ids_object.time.value @@ -447,8 +447,12 @@ def get_available_ids_and_times(db_entry_object) -> list: time_array = [-np.inf] except Exception as e: logger.debug(f"{e}") - time_array = [] - logger.info(f"ERROR! IDS {_ids_name} : Reading time array fails due to following problem : {e}") + time_array = None + logger.warning( + "Unable to read the time array for IDS %s: %s", + _ids_name, + e, + ) if occurrence != 0: result.append((f"{_ids_name}/{occurrence}", time_array)) else: diff --git a/idstools/utils/matplotlib_backend.py b/idstools/utils/matplotlib_backend.py new file mode 100644 index 00000000..0b726b1a --- /dev/null +++ b/idstools/utils/matplotlib_backend.py @@ -0,0 +1,47 @@ +"""Configure Matplotlib before importing Matplotlib itself.""" + +import ast +import os +import sys + + +def _backend_from_cli_rc(argv): + """Return a backend requested by a ``--rc backend=...`` option.""" + for index, argument in enumerate(argv): + if argument == "--rc" and index + 1 < len(argv): + rc_string = argv[index + 1] + elif argument.startswith("--rc="): + rc_string = argument.split("=", 1)[1] + else: + continue + for item in rc_string.split(";"): + key, separator, value = item.partition("=") + if separator and key.strip() == "backend": + value = value.strip() + try: + value = ast.literal_eval(value) + except (SyntaxError, ValueError): + pass + return str(value) + return None + + +def _configure_backend_from_cli_rc(argv=None): + """Set ``MPLBACKEND`` from command-line rcParams before Matplotlib loads.""" + requested_backend = _backend_from_cli_rc(sys.argv if argv is None else argv) + if requested_backend: + os.environ["MPLBACKEND"] = requested_backend + return requested_backend + + +def _is_jupyter(): + """Return True if running inside a Jupyter notebook/lab/Colab kernel.""" + try: + from IPython import get_ipython + + shell = get_ipython() + if shell is None: + return False + return shell.__class__.__name__ in ("ZMQInteractiveShell", "Shell") + except ImportError: + return False diff --git a/idstools/utils/utility_functions.py b/idstools/utils/utility_functions.py index d65458b6..d8544b2f 100644 --- a/idstools/utils/utility_functions.py +++ b/idstools/utils/utility_functions.py @@ -1,8 +1,33 @@ import logging +from urllib.parse import urlsplit, urlunsplit logger = logging.getLogger(f"module.{__name__}") +def add_query_to_uri(uri: str, *, query: str, backend=None): + """Add a query to the URI. + + If a backend is given, add the query only when it matches the URI's + backend. Otherwise, return the original URI unchanged. When backend is + None, add the query to any URI. + """ + uri_parts = urlsplit(uri) + uri_backend = uri_parts.path.rsplit("/", 1)[-1] + if backend is not None and uri_backend != backend: + return uri + + updated_query = f"{uri_parts.query};{query}" + return urlunsplit( + ( + uri_parts.scheme, + uri_parts.netloc, + uri_parts.path, + updated_query, + uri_parts.fragment, + ) + ) + + def get_slice_from_array(arr, slice_str): if ":" not in slice_str: index = int(slice_str) diff --git a/idstools/view/common.py b/idstools/view/common.py index e8ba5c3c..67ad60f9 100644 --- a/idstools/view/common.py +++ b/idstools/view/common.py @@ -2,49 +2,21 @@ import os import sys -import matplotlib +from idstools.utils.matplotlib_backend import ( + _configure_backend_from_cli_rc, + _is_jupyter, +) +_configure_backend_from_cli_rc() -def _is_jupyter() -> bool: - """Return True if running inside a Jupyter notebook/lab/Colab kernel.""" - try: - from IPython import get_ipython - - shell = get_ipython() - if shell is None: - return False - shell_class = shell.__class__.__name__ - # ZMQInteractiveShell: Jupyter Notebook/Lab - # Shell: Google Colab - return shell_class in ("ZMQInteractiveShell", "Shell") - except ImportError: - return False - - -# Select the appropriate matplotlib backend -if _is_jupyter(): - if "matplotlib.pyplot" not in sys.modules: - try: - import ipympl # noqa: F401 - imported to check availability - - matplotlib.use("widget") - except ImportError: - matplotlib.use("agg") -elif sys.platform.startswith("win") or "DISPLAY" in os.environ: - - try: - import tkinter # noqa: F401 - imported to check availability - - matplotlib.use("TkAgg") - except (ImportError, ModuleNotFoundError): - matplotlib.use("agg") -else: - matplotlib.use("agg") +import matplotlib # noqa: E402 - backend env must be set before importing matplotlib import matplotlib.pyplot as plt # noqa: E402 logger = logging.getLogger("module") +PROVENANCE_TITLE_STYLE = {"fontsize": 8, "fontweight": "normal"} + current_directory = os.path.abspath(os.path.dirname(__file__)) # reach to `share` directory (sys.prefix won't work if using --prefix option) share_directory = os.path.abspath(os.path.join(current_directory, "../../../../../")) @@ -217,6 +189,9 @@ def set_sup_title(self, text="", *args, **kwargs): >>> canvas = PlotCanvas(nrows=2, ncols=2) >>> canvas.set_sup_title("Main Figure Title", fontsize=16, fontweight='bold') """ + for key, value in PROVENANCE_TITLE_STYLE.items(): + kwargs.setdefault(key, value) + kwargs.setdefault("y", 0.985) plt.suptitle(text, *args, **kwargs) def show(self, *args, **kwargs): @@ -234,8 +209,8 @@ def show(self, *args, **kwargs): None Notes: - Uses the TkAgg backend for window resizing when available. - Other backends (agg, Qt) may not support window maximization. + Window maximization depends on the backend selected by Matplotlib. + Some backends may not support it. Examples: >>> canvas = PlotCanvas() @@ -249,14 +224,14 @@ def show(self, *args, **kwargs): from IPython.display import display display(self.fig) - if backend != "module://matplotlib_ipympl.backend_nbagg": - plt.close("all") + if backend not in ("widget", "ipympl", "module://ipympl.backend_nbagg"): + plt.close(self.fig) except ImportError: pass return wm = self.get_current_fig_manager() try: - # Try to maximize the window (only works with TkAgg backend) + # Try to maximize the window when the active backend exposes one. window = wm.window screen_y = window.winfo_screenheight() screen_x = window.winfo_screenwidth() diff --git a/idstools/view/core_sources.py b/idstools/view/core_sources.py index 247c5878..67ea777e 100644 --- a/idstools/view/core_sources.py +++ b/idstools/view/core_sources.py @@ -209,7 +209,6 @@ def view_power_and_particle_waveforms(self, ax, time_slice, *args, **kwargs): """ The function `view_power_and_particle_waveforms` plots power waveforms for different sources and particles over time. - Args: ax: The parameter `ax` is an instance of the `Axes` class from the `matplotlib.pyplot` module. It represents the axes on which the waveforms will be plotted. @@ -231,37 +230,71 @@ def view_power_and_particle_waveforms(self, ax, time_slice, *args, **kwargs): ) source_names = self.core_sources_compute.get_source_names(time_slice) ax.set_title("POWER AND PARTICLE WAVEFORMS") + + # --- Marker spacing: ~8 markers per curve --- + markevery = max(1, ntime // 8) + + # --- Per-source markers cycle --- + markers = ["o", "s", "D", "^", "v", "P", "X", "*"] + ax.plot( time_array, single_and_total_electrons_ions_waveforms["total_power_waveform"] * 1.0e-6, label=r"Total electrons+ions", + linestyle="-", + linewidth=3.5, + alpha=0.5, + color="black", ) ax.plot( time_array, single_and_total_electrons_waveforms["total_electron_power_waveform"] * 1.0e-6, label=r"Total electrons", + linestyle="--", + linewidth=3.0, + alpha=0.5, + color="tab:blue", ) ax.plot( time_array, single_and_total_ions_waveforms["total_ion_power_waveform"] * 1.0e-6, label=r"Total ions", + linestyle=":", + linewidth=3.0, + alpha=0.5, + color="tab:red", ) - for isource, name in source_names.items(): + marker = markers[isource % len(markers)] ax.plot( time_array, single_and_total_electrons_ions_waveforms["single_power_waveform"][isource] * 1.0e-6, label=name + " [" + str(isource) + "]" + " electrons+ions", + linestyle="-", + linewidth=1.8, + marker=marker, + markersize=5, + markevery=markevery, ) ax.plot( time_array, single_and_total_electrons_waveforms["single_electron_power_waveform"][isource] * 1.0e-6, label=name + " [" + str(isource) + "]" + " electrons", + linestyle="--", + linewidth=1.5, + marker=marker, + markersize=4, + markevery=markevery, ) ax.plot( time_array, single_and_total_ions_waveforms["single_ion_power_waveform"][isource] * 1.0e-6, label=name + " [" + str(isource) + "]" + " ions", + linestyle=":", + linewidth=1.5, + marker=marker, + markersize=4, + markevery=markevery, ) ax.set_ylabel("Power waveforms $\\mathrm{[MW]}$") ax.set_xlabel("Time (s)") diff --git a/idstools/view/domain/ecstray.py b/idstools/view/domain/ecstray.py index 66bfacec..591951b3 100644 --- a/idstools/view/domain/ecstray.py +++ b/idstools/view/domain/ecstray.py @@ -84,27 +84,34 @@ def plot_resonance_layer(self, ax, coherent_wave_index, time_slice, init=1, verb else: ax.set_data(res_layer[i_harm]["r"], res_layer[i_harm]["z"]) - def plot_poloidal_view(self, ax, coherent_wave_index, time_slice): + def plot_poloidal_view(self, ax, coherent_wave_index, time_slice, equilibrium_time_slice=None): n_harm = [1, 2, 3, 4] - - resonance_data = self.ecstray_object.get_resonance_layer(coherent_wave_index, time_slice, n_harm=n_harm) + if equilibrium_time_slice is None: + equilibrium_time_slice = time_slice + + resonance_data = self.ecstray_object.get_resonance_layer( + coherent_wave_index, + time_slice, + n_harm=n_harm, + equilibrium_time_slice=equilibrium_time_slice, + ) profile2d_index = resonance_data["profile2d_index"] resonance_layer = resonance_data["resonance_layer"] - grid_data = self.equilibrium_compute.get2d_cartesian_grid(time_slice, profile2d_index) + grid_data = self.equilibrium_compute.get2d_cartesian_grid(equilibrium_time_slice, profile2d_index) r2d = grid_data["r2d"] z2d = grid_data["z2d"] psi2d = grid_data["psi2d"] - rho2d = self.equilibrium_compute.get_rho2d(time_slice, profile2d_index) + phi2d = self.equilibrium_compute.get_phi2d(equilibrium_time_slice, profile2d_index) # Poloidal view plot - contour_lines = ax.contour(r2d, z2d, psi2d, 50, cmap="summer") + contour_lines = ax.contour(r2d, z2d, psi2d.T, 50, cmap="summer") cbar_psi = plt.colorbar(contour_lines, ax=ax, orientation="horizontal", pad=0.08, fraction=0.03) cbar_psi.set_label(r"$\psi$ [Wb]") - if rho2d is not None and len(rho2d) > 0: - contour_lines_rho = ax.contour(r2d, z2d, rho2d, 50, cmap="YlOrBr") + if phi2d is not None and len(phi2d) > 0: + contour_lines_rho = ax.contour(r2d, z2d, phi2d.T, 50, cmap="YlOrBr") cbar_rho = plt.colorbar(contour_lines_rho, ax=ax, orientation="horizontal", pad=0.08, fraction=0.03) - cbar_rho.set_label(r"$\rho$ [Wb]") + cbar_rho.set_label(r"$\Phi$ [Wb]") # ax_polview.set_xlim(r2d.min(),r2d.max()) ax.set_title("Poloidal view (R,Z)") ax.set_xlabel("R [m]", labelpad=lpad) diff --git a/idstools/view/domain/mdplot.py b/idstools/view/domain/mdplot.py index e93fc1ed..7e6a338a 100644 --- a/idstools/view/domain/mdplot.py +++ b/idstools/view/domain/mdplot.py @@ -36,14 +36,26 @@ def update_labels(ax): ax.figure.canvas.draw_idle() -def plot_machine_description(ax, ids_data): +def plot_machine_description(ax, ids_data, main_uri=None): """ The `plotMachineDescription` method is responsible for plotting the machine description based on the provided pulse list. + Args: + main_uri: If provided, machine description entries with the same URI are omitted + from the provenance text (to avoid duplication when MD and main data share a URI). """ - database_path = "" + provenance_groups = {} + + def _add_provenance(ids_name_label, connection_args): + """Group subsystem labels that originate from the same data-entry URI.""" + uri = get_database_path(connection_args).strip() + if main_uri is not None and uri == main_uri.strip(): + return + labels = provenance_groups.setdefault(uri, []) + if ids_name_label not in labels: + labels.append(ids_name_label) mdlegends = [] mdlabels = [] @@ -70,7 +82,7 @@ def plot_machine_description(ax, ids_data): if _legend: mdlegends.append(_legend) mdlabels.append(f"pf_active:{idsocc}/coil[{select}]") - database_path += "pf_active = " + get_database_path(ids_data_and_config["connectionArgs"]) + "\n" + _add_provenance("pf_active", ids_data_and_config["connectionArgs"]) elif ids_name == "tf": select2 = ":" if len(matches) == 2: @@ -81,7 +93,7 @@ def plot_machine_description(ax, ids_data): if _legend: mdlegends.append(_legend) mdlabels.append(f"tf:{idsocc}/coil[{select}]/conductor[{select2}]") - database_path += "tf = " + get_database_path(ids_data_and_config["connectionArgs"]) + "\n" + _add_provenance("tf", ids_data_and_config["connectionArgs"]) elif ids_name == "pf_passive": pfpassiveview = PFPassiveView(ids_data_and_config["idsData"]) if "loop" in idsfield or idsfield == "": @@ -90,7 +102,7 @@ def plot_machine_description(ax, ids_data): mdlegends.append(_legend) mdlabels.append(f"pf_passive:{idsocc}/loop[{select}]") - database_path += "pf_passive = " + get_database_path(ids_data_and_config["connectionArgs"]) + "\n" + _add_provenance("pf_passive", ids_data_and_config["connectionArgs"]) elif ids_name == "wall": wallview = WallView(ids_data_and_config["idsData"]) select2 = ":" @@ -100,7 +112,7 @@ def plot_machine_description(ax, ids_data): wallview.view_wall_vessel(ax, select_description2d=select, select_unit=select2) if "limiter" in idsfield or idsfield == "": wallview.view_wall_limiter(ax, select_description2d=select, select_unit=select2) - database_path += "wall = " + get_database_path(ids_data_and_config["connectionArgs"]) + "\n" + _add_provenance("wall", ids_data_and_config["connectionArgs"]) elif ids_name == "magnetics": magnetics_view = MagneticsView(ids_data_and_config["idsData"]) if "b_field_phi_probe" in idsfield or idsfield == "": @@ -128,11 +140,9 @@ def plot_machine_description(ax, ids_data): if _legend: mdlegends.append(_legend) mdlabels.append(f"magnetics:{idsocc}/shunt[{select}]") - database_path += "magnetics = " + get_database_path(ids_data_and_config["connectionArgs"]) + "\n" + _add_provenance("magnetics", ids_data_and_config["connectionArgs"]) else: - database_path += ( - f"{ids_name} = " + get_database_path(ids_data_and_config["connectionArgs"]) + "No visualization yet\n" - ) + _add_provenance(f"{ids_name} (No visualization yet)", ids_data_and_config["connectionArgs"]) logger.info(f"Visualization is not implemented yet for machine description {ids_name}") handles, labels = ax.get_legend_handles_labels() @@ -151,14 +161,4 @@ def plot_machine_description(ax, ids_data): # ax.callbacks.connect("ylim_changed", update_labels) ax.plot() - xmin, xmax = ax.get_xlim() - ymin, ymax = ax.get_ylim() - ax.text( - xmax + 0.01 * abs(xmax), - ymin + 0.5 * abs(ymax - ymin), - f"{database_path}", - horizontalalignment="left", - verticalalignment="center", - rotation="vertical", - fontsize=7, - ) + return " | ".join(f"{', '.join(labels)} = {uri}" for uri, labels in provenance_groups.items()) diff --git a/idstools/view/edge_profiles.py b/idstools/view/edge_profiles.py index 225e9fe0..13e12b4a 100644 --- a/idstools/view/edge_profiles.py +++ b/idstools/view/edge_profiles.py @@ -1,6 +1,7 @@ import logging import numpy as np +from matplotlib.collections import PolyCollection from idstools.compute.edge_profiles import EdgeProfilesCompute @@ -109,8 +110,8 @@ def _print_specis_concentration(self, composition_data): def view_electrons_density(self, ax, time_slice, show_separatrix=False): """ - The function `view_electrons_density` plots the electron density on a rectangular grid and adds a - separatrix line. + The function `view_electrons_density` plots the electron density on the original GGD mesh + using tripcolor (Gouraud shading) and adds a separatrix line. Args: ax: The parameter "ax" is an instance of the matplotlib Axes class. It represents the axes on @@ -119,25 +120,31 @@ def view_electrons_density(self, ax, time_slice, show_separatrix=False): plotted. It is an optional parameter with a default value of 0. Defaults to 0 Returns: - the pcolormesh object 'c'. + the tripcolor object 'c'. """ - x, y = self.edge_profiles_compute.get_rectangular_grid(500) + ne_values = self.edge_profiles_compute.get_electron_density_on_cells(time_slice) - ne_edge = self.edge_profiles_compute.get_electron_density(time_slice, x, y) - if ne_edge is not None: + if ne_values is not None: + polygons, cell_values = self.edge_profiles_compute.get_cell_polygons_and_values(time_slice, ne_values) + all_r = np.concatenate([v[:, 0] for v in polygons]) + all_z = np.concatenate([v[:, 1] for v in polygons]) ax.grid(False) - c = ax.pcolormesh(x, y, ne_edge, vmin=0, vmax=5e19, shading="auto") + pc = PolyCollection(polygons, array=cell_values, cmap="coolwarm", linewidths=0) + pc.set_clim(float(cell_values.min()), float(cell_values.max())) + ax.add_collection(pc) + ax.set_xlim(all_r.min(), all_r.max()) + ax.set_ylim(all_z.min(), all_z.max()) core_boundry = self.edge_profiles_compute.get_core_boundry(time_slice) - ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="r", linewidth=0) + ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="none") if show_separatrix: separatrix = self.edge_profiles_compute.get_separatrix(time_slice) if separatrix is not None: - ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x") + ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x", s=10) ax.set_aspect("equal", adjustable="box") ax.set_xlabel("R,m") ax.set_ylabel("Z,m") ax.set_title("Electron density") - return c + return pc else: xmin, xmax = ax.get_xlim() ymin, ymax = ax.get_ylim() @@ -152,7 +159,8 @@ def view_electrons_density(self, ax, time_slice, show_separatrix=False): def view_ion_density(self, ax, time_slice, show_separatrix=False): """ - The function `view_ion_density` plots the ion density on a rectangular grid and adds a separatrix line. + The function `view_ion_density` plots the ion density on the original GGD mesh + using tripcolor (Gouraud shading) and adds a separatrix line. Args: ax: The parameter "ax" is an instance of the matplotlib Axes class. It represents the axes on @@ -161,26 +169,31 @@ def view_ion_density(self, ax, time_slice, show_separatrix=False): be plotted. It is an optional parameter with a default value of 0. Defaults to 0 Returns: - the pcolormesh object 'c'. + the tripcolor object 'c'. """ - x, y = self.edge_profiles_compute.get_rectangular_grid(500) + ni_values = self.edge_profiles_compute.get_ion_density_on_cells(time_slice) - ni_edge = self.edge_profiles_compute.get_ion_density(time_slice, x, y) - if ni_edge is not None: + if ni_values is not None: + polygons, cell_values = self.edge_profiles_compute.get_cell_polygons_and_values(time_slice, ni_values) + all_r = np.concatenate([v[:, 0] for v in polygons]) + all_z = np.concatenate([v[:, 1] for v in polygons]) ax.grid(False) - c = ax.pcolormesh(x, y, ni_edge, vmin=0, vmax=5e19, shading="auto") + pc = PolyCollection(polygons, array=cell_values, cmap="coolwarm", linewidths=0) + pc.set_clim(float(cell_values.min()), float(cell_values.max())) + ax.add_collection(pc) + ax.set_xlim(all_r.min(), all_r.max()) + ax.set_ylim(all_z.min(), all_z.max()) core_boundry = self.edge_profiles_compute.get_core_boundry(time_slice) - ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="r", linewidth=0) + ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="none") if show_separatrix: separatrix = self.edge_profiles_compute.get_separatrix(time_slice) if separatrix is not None: - ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x") - + ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x", s=10) ax.set_aspect("equal", adjustable="box") ax.set_xlabel("R,m") ax.set_ylabel("Z,m") ax.set_title("Ion density") - return c + return pc else: xmin, xmax = ax.get_xlim() ymin, ymax = ax.get_ylim() @@ -195,8 +208,8 @@ def view_ion_density(self, ax, time_slice, show_separatrix=False): def view_neutral_density(self, ax, time_slice, show_separatrix=False): """ - The function `view_neutral_density` plots the neutral density on a rectangular grid and adds a - separatrix line. + The function `view_neutral_density` plots the neutral density on the original GGD mesh + using tripcolor (Gouraud shading) and adds a separatrix line. Args: ax: The parameter "ax" is an instance of the matplotlib Axes class. It represents the axes @@ -205,26 +218,33 @@ def view_neutral_density(self, ax, time_slice, show_separatrix=False): will be plotted. It is an optional parameter with a default value of 0. Defaults to 0 Returns: - the pcolormesh object 'c'. + the tripcolor object 'c'. """ - x, y = self.edge_profiles_compute.get_rectangular_grid(500) - - n_neutral_edge = self.edge_profiles_compute.get_neutral_density(time_slice, x, y) + n_neutral_values = self.edge_profiles_compute.get_neutral_density_on_cells(time_slice) - if n_neutral_edge is not None: + if n_neutral_values is not None: + polygons, cell_values = self.edge_profiles_compute.get_cell_polygons_and_values( + time_slice, n_neutral_values + ) + all_r = np.concatenate([v[:, 0] for v in polygons]) + all_z = np.concatenate([v[:, 1] for v in polygons]) ax.grid(False) - c = ax.pcolormesh(x, y, n_neutral_edge, vmin=0, vmax=5e19, shading="auto") + pc = PolyCollection(polygons, array=cell_values, cmap="coolwarm", linewidths=0) + pc.set_clim(float(cell_values.min()), float(cell_values.max())) + ax.add_collection(pc) + ax.set_xlim(all_r.min(), all_r.max()) + ax.set_ylim(all_z.min(), all_z.max()) core_boundry = self.edge_profiles_compute.get_core_boundry(time_slice) - ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="r", linewidth=0) + ax.fill(core_boundry[:, 0], core_boundry[:, 1], facecolor="w", edgecolor="none") if show_separatrix: separatrix = self.edge_profiles_compute.get_separatrix(time_slice) if separatrix is not None: - ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x") + ax.scatter(separatrix[:, 0], separatrix[:, 1], color="#FF1493", marker="x", s=10) ax.set_aspect("equal", adjustable="box") ax.set_xlabel("R,m") ax.set_ylabel("Z,m") ax.set_title("Neutral density") - return c + return pc else: xmin, xmax = ax.get_xlim() ymin, ymax = ax.get_ylim() @@ -237,23 +257,34 @@ def view_neutral_density(self, ax, time_slice, show_separatrix=False): ) return None - def view_equatorial_plane_and_diverter_density(self, ax, time_slice, logscale=False): - x, y = self.edge_profiles_compute.get_rectangular_grid(500) - ne_edge = self.edge_profiles_compute.get_electron_density(time_slice, x, y) - if ne_edge is not None: - # choose Z position for a radial profile: - z0 = 0.0 - ind = np.argmin(abs(y[:, 0] - z0)) - ax.plot(x[ind, :], ne_edge[ind, :], label="Equatorial plane") + def view_equatorial_plane_and_diverter_density(self, ax, time_slice, logscale=True): + OUTER_MIDPLANE_INDEX = 11 + OUTER_TARGET_INDEX = 13 + INNER_TARGET_INDEX = 14 - z0 = -4.0 - ind = np.argmin(abs(y[:, 0] - z0)) - ax.plot(x[ind, :], ne_edge[ind, :], label="Divertor") + any_data = False + + for subset_idx, label in ( + (OUTER_MIDPLANE_INDEX, "Outer midplane"), + (OUTER_TARGET_INDEX, "Outer divertor target"), + (INNER_TARGET_INDEX, "Inner divertor target"), + ): + r_t, z_t, ne_t = self.edge_profiles_compute.get_rz_and_ne_on_subset(time_slice, subset_idx) + if r_t is None or len(r_t) == 0: + continue + # Arc-length as x-axis + coords = np.column_stack([r_t, z_t]) + seg_len = np.sqrt(np.sum(np.diff(coords, axis=0) ** 2, axis=1)) + arc = np.concatenate([[0.0], np.cumsum(seg_len)]) + ax.plot(arc, ne_t, label=label) + any_data = True + + if any_data: if logscale: ax.set_yscale("log") ax.set_title("Electron density") - ax.set_xlabel("R,m") - # ax.set_ylim([0, 1.5e21]) + ax.set_xlabel("Distance along profile (m)") + ax.set_ylabel("n$_e$ (m$^{-3}$)") ax.legend() else: xmin, xmax = ax.get_xlim() diff --git a/idstools/view/equilibrium.py b/idstools/view/equilibrium.py index 2416967e..590c31bd 100644 --- a/idstools/view/equilibrium.py +++ b/idstools/view/equilibrium.py @@ -13,6 +13,7 @@ except ImportError: import imas import matplotlib.pyplot as plt +from matplotlib.lines import Line2D as ProxyLine import numpy as np from idstools.compute.equilibrium import EquilibriumCompute @@ -40,7 +41,13 @@ def view_magnetic_poloidal_flux( ax: plt.axes, time_slice: int, profiles2d_index: int = 0, - plot_rho: bool = False, + plot_magnetic_axis: bool = True, + plot_current_centre: bool = True, + plot_boundary_data: bool = True, + plot_phi: bool = False, + plot_annotations: bool = True, + plot_psi: bool = True, + plot_boundary_outline: bool = False, ): """ This function plots the magnetic poloidal flux contours on a 2D Cartesian grid. @@ -80,12 +87,22 @@ def view_magnetic_poloidal_flux( :meth:`plotIP` """ contour_lines_psi = contour_lines_rho = None - cartestion_grid = self.compute_obj.get2d_cartesian_grid(time_slice, profiles2d_index) - if cartestion_grid is not None: - levels = 50 - + levels = 50 + cartestion_grid = None + if plot_psi or plot_phi: + cartestion_grid = self.compute_obj.get2d_cartesian_grid(time_slice, profiles2d_index) + if cartestion_grid is not None and plot_psi: + + # As per IMAS data dictionary psi is stored as [R, Z] with shape (N_R, N_Z). + # Check this reference : + # https://imas-data-dictionary.readthedocs.io/en/latest/generated/identifier/poloidal_plane_coordinates_identifier.html + + # For matplotlib contour as per the documentation: + # https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.contour.html + # matplotlib.contour(r, z, Zdata) expects Zdata as [rows=z, cols=r], so the shape must be (N_Z, N_R). + # Therefore we transpose psi before plotting. contour_lines_psi = ax.contour( - cartestion_grid["r2d"], cartestion_grid["z2d"], cartestion_grid["psi2d"], levels, cmap="summer" + cartestion_grid["r2d"], cartestion_grid["z2d"], cartestion_grid["psi2d"].T, levels, cmap="summer" ) # ax.clabel( # contour_lines_psi, @@ -95,23 +112,341 @@ def view_magnetic_poloidal_flux( # # fmt="%.2e", # inline_spacing=1, # ) - if plot_rho: - rho2d = self.compute_obj.get_rho2d(time_slice) - if rho2d is not None: + + # phi (toroidal flux) overlay + if plot_phi: + phi2d = self.compute_obj.get_phi2d(time_slice) + if phi2d is not None: contour_lines_rho = ax.contour( - cartestion_grid["r2d"], cartestion_grid["z2d"], rho2d, levels=levels, cmap="YlOrBr" + cartestion_grid["r2d"], cartestion_grid["z2d"], phi2d.T, levels=levels, cmap="YlOrBr" ) ax.set_aspect("equal", adjustable="box") ax.set_xlabel("$R$ [m]") ax.set_ylabel("$Z$ [m]") - # ax.set_xlim(3.4, cartestionGrid["r2d"].max()) - # ax.set_ylim(cartestionGrid["z2d"].min() * 0.7, cartestionGrid["z2d"].max() * 0.7) + + # Get any handles already in the axes legend (e.g. from machine description). + _existing_legend = ax.get_legend() + if _existing_legend is not None: + _md_handles = list(_existing_legend.legend_handles) + _md_labels = [t.get_text() for t in _existing_legend.get_texts()] + else: + _md_handles, _md_labels = [], [] + overlay_entries = [] + + if plot_magnetic_axis: + mag_ax = self.compute_obj.get_magnetic_axis(time_slice) + if mag_ax is not None: + (marker,) = ax.plot( + mag_ax["r"], + mag_ax["z"], + marker="+", + color="saddlebrown", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + zorder=6, + ) + proxy_mag = ProxyLine( + [0], + [0], + color="saddlebrown", + marker="+", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + label="magnetic axis", + ) + overlay_entries.append((proxy_mag, [marker])) + + if plot_current_centre: + cc = self.compute_obj.get_current_centre(time_slice) + if cc is not None: + (marker,) = ax.plot( + cc["r"], + cc["z"], + marker="+", + color="deeppink", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + zorder=6, + ) + proxy_cc = ProxyLine( + [0], + [0], + color="deeppink", + marker="+", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + label="current centre", + ) + overlay_entries.append((proxy_cc, [marker])) + + if plot_boundary_outline: + bd = self.compute_obj.get_boundary_data(time_slice) + if bd["bnd_r"] is not None and bd["bnd_z"] is not None: + (boundary_line,) = ax.plot( + bd["bnd_r"], + bd["bnd_z"], + color="royalblue", + linewidth=2.0, + linestyle="-", + zorder=5, + ) + proxy_boundary = ProxyLine( + [0], [0], color="royalblue", linewidth=2.0, linestyle="-", label="boundary/outline" + ) + overlay_entries.append((proxy_boundary, [boundary_line])) + + if plot_boundary_data: + + bd = self.compute_obj.get_boundary_data(time_slice) + + separatrix_is_boundary = ( + bd["bnd_r"] is not None + and bd["bnd_z"] is not None + and bd["sep_r"] is not None + and bd["sep_z"] is not None + and np.array_equal(np.asarray(bd["sep_r"]), np.asarray(bd["bnd_r"])) + and np.array_equal(np.asarray(bd["sep_z"]), np.asarray(bd["bnd_z"])) + ) + + # boundary_separatrix outline + if bd["sep_r"] is not None and bd["sep_z"] is not None and not separatrix_is_boundary: + (sep_line,) = ax.plot( + bd["sep_r"], + bd["sep_z"], + color="firebrick", + linewidth=2.0, + linestyle="--", + zorder=4, + ) + proxy_sep_bnd = ProxyLine( + [0], [0], color="firebrick", linewidth=2.0, linestyle="--", label="separatrix" + ) + overlay_entries.append((proxy_sep_bnd, [sep_line])) + + # geometric axis + if bd["bnd_geom_r"] is not None and bd["bnd_geom_z"] is not None: + (gax_marker,) = ax.plot( + bd["bnd_geom_r"], + bd["bnd_geom_z"], + marker="x", + color="darkcyan", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + zorder=6, + ) + proxy_gax = ProxyLine( + [0], + [0], + color="darkcyan", + marker="x", + markersize=6, + markeredgewidth=1.4, + linestyle="None", + label="Geometric axis", + ) + overlay_entries.append((proxy_gax, [gax_marker])) + + # x-points (boundary_separatrix) + point_marker_size = 7 + point_marker_edgewidth = 2.0 + point_label_fontsize = 8 + + _xp_groups = [ + (bd["sep_xpoints"], "red", "x_point"), + ] + for xp_list, xp_color, xp_label in _xp_groups: + _xp_artists = [] + for xp_idx, (xr, xz) in enumerate(xp_list): + (mk,) = ax.plot( + xr, + xz, + marker="x", + color=xp_color, + markersize=point_marker_size, + markeredgewidth=point_marker_edgewidth, + linestyle="None", + zorder=7, + ) + ann = ax.annotate( + f"X{xp_idx}", + xy=(xr, xz), + xytext=(-6, 6), + textcoords="offset points", + fontsize=point_label_fontsize, + ha="right", + color=xp_color, + fontweight="bold", + zorder=8, + ) + _xp_artists.append(mk) + _xp_artists.append(ann) + if _xp_artists: + proxy_xp = ProxyLine( + [0], + [0], + color=xp_color, + marker="x", + markersize=point_marker_size, + markeredgewidth=point_marker_edgewidth, + linestyle="None", + label=xp_label, + ) + overlay_entries.append((proxy_xp, _xp_artists)) + + # strike-points (boundary_separatrix) + _sp_groups = [ + (bd["sep_strikepoints"], "red", "strike_point"), + ] + for sp_list, sp_color, sp_label in _sp_groups: + _sp_artists = [] + for sp_idx, (sr, sz) in enumerate(sp_list): + (mk,) = ax.plot( + sr, + sz, + marker="+", + color=sp_color, + markersize=point_marker_size, + markeredgewidth=point_marker_edgewidth, + linestyle="None", + zorder=7, + ) + ann = ax.annotate( + f"S{sp_idx}", + xy=(sr, sz), + xytext=(-6, 6), + textcoords="offset points", + fontsize=point_label_fontsize, + ha="right", + color=sp_color, + fontweight="bold", + zorder=8, + ) + _sp_artists.append(mk) + _sp_artists.append(ann) + if _sp_artists: + proxy_sp = ProxyLine( + [0], + [0], + color=sp_color, + marker="+", + markersize=point_marker_size, + markeredgewidth=point_marker_edgewidth, + linestyle="None", + label=sp_label, + ) + overlay_entries.append((proxy_sp, _sp_artists)) + + if plot_annotations: + self.view_global_quantities_annotation(ax, time_slice) + + # --- clickable legend + if overlay_entries or _md_handles: + overlay_proxies = [proxy for proxy, _ in overlay_entries] + + all_handles = _md_handles + overlay_proxies + all_labels = _md_labels + [p.get_label() for p in overlay_proxies] + + legend = ax.legend( + handles=all_handles, + labels=all_labels, + loc="upper left", + bbox_to_anchor=(1.15, 1), + fancybox=True, + frameon=False, + framealpha=1.0, + facecolor="white", + fontsize=10, + labelspacing=1.2, + ) + legend.set_zorder(1000) + for text in legend.get_texts(): + text.set_ha("left") + + leg_map = {} + legend_texts = legend.get_texts() + n_md = len(_md_handles) + for i, orig_artist in enumerate(_md_handles): + leg_h = legend.legend_handles[i] + leg_h.set_picker(8) + leg_map[leg_h] = [orig_artist] + legend_texts[i].set_picker(True) + leg_map[legend_texts[i]] = [orig_artist] + + for i, (_, artists) in enumerate(overlay_entries): + leg_h = legend.legend_handles[n_md + i] + leg_text = legend_texts[n_md + i] + leg_h.set_picker(8) + leg_map[leg_h] = artists + leg_text.set_picker(True) + leg_map[leg_text] = artists + if artists and not artists[0].get_visible(): + leg_h.set_alpha(0.3) + leg_text.set_alpha(0.3) + + def on_legend_click(event): + legline = event.artist + if legline not in leg_map: + return + artists = leg_map[legline] + if not artists: + return + visible = not artists[0].get_visible() + for a in artists: + a.set_visible(visible) + legline.set_alpha(1.0 if visible else 0.3) + if legline in legend.legend_handles: + leg_index = legend.legend_handles.index(legline) + legend_texts[leg_index].set_alpha(1.0 if visible else 0.3) + elif legline in legend_texts: + leg_index = legend_texts.index(legline) + legend.legend_handles[leg_index].set_alpha(1.0 if visible else 0.3) + ax.figure.canvas.draw_idle() + + ax.figure.canvas.mpl_connect("pick_event", on_legend_click) + return contour_lines_psi, contour_lines_rho def view_pulse_info(self, ax: plt.axes, title: str, hostdir: str, shot: int, run: int, t: float): self.database_info(ax, title, hostdir, shot, run, t) + def view_global_quantities_annotation(self, ax: plt.axes, time_slice: int): + """Draw a scalar global-quantities text box below the axes. + + Reads validated scalars via + :meth:`idstools.compute.equilibrium.EquilibriumCompute.get_scalar_annotation_quantities` + and renders them as a styled text box just below the axes. + + Args: + ax: matplotlib axes. + time_slice (int): time-slice index. + + Returns: + matplotlib ``Text`` artist, or ``None`` if no valid data. + """ + items = self.compute_obj.get_scalar_annotation_quantities(time_slice) + if not items: + return None + + textstr = "\n".join(f"{d['label']} = {d['text']}" for d in items) + txt = ax.text( + 1.15, + 0.0, + textstr, + transform=ax.transAxes, + fontsize=9, + horizontalalignment="left", + verticalalignment="bottom", + clip_on=False, + bbox=dict(boxstyle="round,pad=0.5", facecolor="white", alpha=1.0, edgecolor="none"), + ) + return txt + def plot_ip(self, ax): """ This function plots the plasma current over time on a given axis. @@ -152,9 +487,9 @@ def plot_poloidal_equilibrium(self, ax, time_slice: int): data = self.compute_obj.get_flux_surfaces(time_slice) r2d = data["r2d"] z2d = data["z2d"] - # rho2d = data["rho2d"] + # phi2d = data["phi2d"] psi2d = data["psi2d"] - cntr = ax.contour(r2d, z2d, psi2d, 50, cmap="summer") + cntr = ax.contour(r2d, z2d, psi2d.T, 50, cmap="summer") cbar = plt.colorbar(cntr, ax=ax, pad=0.08, fraction=0.03) cbar.set_label(r"$\psi$ [Wb]") # if len(rho2d)>0: @@ -207,7 +542,7 @@ def plotequilibrium(self, ax, time_slice): ax.xaxis.tick_top() ax.xaxis.set_label_position("top") - contour_lines = ax.contour(r2d, z2d, psi2d, levels=50, cmap="summer") # ,label=r'$\Psi_{pol}$') + contour_lines = ax.contour(r2d, z2d, psi2d.T, levels=50, cmap="summer") # ,label=r'$\Psi_{pol}$') cbar = plt.colorbar(contour_lines, ax=ax, pad=0.08, fraction=0.03) cbar.set_label(r"$\psi$ [Wb]") ax.set_xlim(r2d.min(), r2d.max()) @@ -247,13 +582,12 @@ def plot_profiles_1d_quantities(self, axes_list, time_slice, attributes=None): coordinate_normalized = (psi - psi_first) / (psi_last - psi_first) axes_list[counter].plot( - coordinate_normalized, copied_field, label=f"{field.metadata.name} ({field.metadata.units})" + coordinate_normalized, copied_field, label=f"{field.metadata.name} [{field.metadata.units}]" ) if coordinate.metadata.name == "psi": - axes_list[counter].set_xlabel(f"{coordinate.metadata.name} (normalized)") + axes_list[counter].set_xlabel(f"{coordinate.metadata.name}_norm [1]") else: - axes_list[counter].set_xlabel(f"{coordinate.metadata.name} ({coordinate.metadata.units})") - axes_list[counter].set_ylabel(name) + axes_list[counter].set_xlabel(f"{coordinate.metadata.name} [{coordinate.metadata.units}]") axes_list[counter].legend(loc="upper right") counter = counter + 1 @@ -267,11 +601,10 @@ def plot_global_quantities(self, axes_list, time_slice, attributes=None): field["node"][field["node"] == imas.ids_defs.EMPTY_FLOAT] = np.nan if field["has_value"]: if len(field["node"]) < 5: - axes_list[counter].scatter(field["coordinate"], field["node"], label=f"{name} ({field['unit']})") + axes_list[counter].scatter(field["coordinate"], field["node"], label=f"{name} [{field['unit']}]") else: - axes_list[counter].plot(field["coordinate"], field["node"], label=f"{name} ({field['unit']})") - axes_list[counter].set_xlabel(f"{field['coordinate_name']} ({field['coordinate_unit']})") - axes_list[counter].set_ylabel(name) + axes_list[counter].plot(field["coordinate"], field["node"], label=f"{name} [{field['unit']}]") + axes_list[counter].set_xlabel(f"{field['coordinate_name']} [{field['coordinate_unit']}]") self.view_time_line(axes_list[counter], time_slice) axes_list[counter].legend(loc="upper right") counter = counter + 1 @@ -424,7 +757,9 @@ def view_profile_plot(self, ax, time_index1, equilibrium2_ids=None, time_index2= logger.warning(f"Equilibrium1 {name}: No valid r or jtor data available") # Set ylim with check for identical values to avoid matplotlib warning - if abs(y_max - y_min) < 1e-10: + if not np.isfinite(y_min) or not np.isfinite(y_max): + y_min, y_max = -1.0, 1.0 + elif abs(y_max - y_min) < 1e-10: # If y_min and y_max are essentially equal, create a small range around the value if abs(y_min) < 1e-10: # If both are near zero, use a default range @@ -619,14 +954,16 @@ def view_current_plot(self, ax, time_index1, equilibrium2_ids=None): line31.set_xdata(time) line31.set_ydata(ip / 1e6) xlims = np.array([min(time), max(time)]) - ylims = np.array([min(ip), max(ip)]) / 1e6 + ip_valid = ip[np.isfinite(ip)] + ylims = np.array([ip_valid.min(), ip_valid.max()]) / 1e6 if ip_valid.size > 0 else np.array([-1.0, 1.0]) if data2 is not None: line32.set_xdata(timeE) line32.set_ydata(ipE / 1e6) + ipE_valid = ipE[np.isfinite(ipE)] xlims[0] = np.minimum(xlims[0], np.min(timeE)) - ylims[0] = np.minimum(ylims[0], np.min(ipE / 1e6)) + ylims[0] = np.minimum(ylims[0], np.min(ipE_valid / 1e6)) if ipE_valid.size > 0 else ylims[0] xlims[1] = np.maximum(xlims[1], np.max(timeE)) - ylims[1] = np.maximum(ylims[1], np.max(ipE / 1e6)) + ylims[1] = np.maximum(ylims[1], np.max(ipE_valid / 1e6)) if ipE_valid.size > 0 else ylims[1] dy = ylims[1] - ylims[0] if dy > 0: ylims[0] = ylims[0] - 0.01 * dy diff --git a/pyproject.toml b/pyproject.toml index 611d08cc..1f3a64f9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ dependencies = [ "numpy", "packaging", "pandas", + "plotext>=5.0", "pyparsing", "python-dateutil", "pyyaml", @@ -78,6 +79,7 @@ plotpressure = "idstools.scripts:plotpressure" plotrotation = "idstools.scripts:plotrotation" plotequicomp = "idstools.scripts:plotequicomp" plotcoresources = "idstools.scripts:plotcoresources" +plotcoretransport = "idstools.scripts:plotcoretransport" plotkineticprofiles = "idstools.scripts:plotkineticprofiles" plothcdwaves = "idstools.scripts:plothcdwaves" plotecray = "idstools.scripts:plotecray" diff --git a/resources/input/ec_waveforms.yaml b/resources/input/ec_waveforms.yaml index b1292dd3..e8a45c3c 100644 --- a/resources/input/ec_waveforms.yaml +++ b/resources/input/ec_waveforms.yaml @@ -68,17 +68,17 @@ dynamic_variables: beam[9].frequency.data: [1.70000e+11] beam[10].frequency.data: [1.70000e+11] - beam[0].mode: -1 # (1) O-mode, (-1) X-mode - beam[1].mode: -1 - beam[2].mode: -1 - beam[3].mode: -1 - beam[4].mode: -1 - beam[5].mode: -1 - beam[6].mode: -1 - beam[7].mode: -1 - beam[8].mode: -1 - beam[9].mode: -1 - beam[10].mode: -1 + beam[0].polarization.o_mode_fraction: [0.0] # 1.0 = O-mode, 0.0 = X-mode + beam[1].polarization.o_mode_fraction: [0.0] + beam[2].polarization.o_mode_fraction: [0.0] + beam[3].polarization.o_mode_fraction: [0.0] + beam[4].polarization.o_mode_fraction: [0.0] + beam[5].polarization.o_mode_fraction: [0.0] + beam[6].polarization.o_mode_fraction: [0.0] + beam[7].polarization.o_mode_fraction: [0.0] + beam[8].polarization.o_mode_fraction: [0.0] + beam[9].polarization.o_mode_fraction: [0.0] + beam[10].polarization.o_mode_fraction: [0.0] beam[0].steering_angle_pol: [0.00178] beam[1].steering_angle_pol: [0.30776] diff --git a/tests/st03_test_analysis_scripts_with_uri.sh b/tests/st03_test_analysis_scripts_with_uri.sh index f2b7234a..d4fd805e 100755 --- a/tests/st03_test_analysis_scripts_with_uri.sh +++ b/tests/st03_test_analysis_scripts_with_uri.sh @@ -33,7 +33,7 @@ SCRIPTS=( "plotecray --uri \"imas:mdsplus?user=public;pulse=134173;run=2326;database=TEST;version=3\" --save --directory $LOG_DIR" "plotecray --uri \"imas:mdsplus?user=public;pulse=134173;run=2326;database=TEST;version=3\" -md wall --save --directory $LOG_DIR" "plotedgeprofiles --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --wall --time 60 --save --directory $LOG_DIR" - "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --rho -md \"imas:hdf5?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active\" \"imas:hdf5?user=public;pulse=116000;run=5;database=ITER_MD;version=3#wall\" --save --directory $LOG_DIR" + "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --phi -md \"imas:hdf5?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active\" \"imas:hdf5?user=public;pulse=116000;run=5;database=ITER_MD;version=3#wall\" --save --directory $LOG_DIR" "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" -md pf_active wall --save --directory $LOG_DIR" "printfluxes --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" -m CLOSEST" "plothcd -ech 134173/101/public/MDSPLUS/TEST/3 -nbi 130012/115/public/MDSPLUS/TEST/3 -fus 130012/115/public/MDSPLUS/TEST/3 -icrh 130012/115/public/MDSPLUS/TEST/3 --save --directory $LOG_DIR" @@ -64,7 +64,7 @@ return $? # "plotecray --uri \"imas:hdf5?user=public;pulse=134173;run=2326;database=TEST;version=3\" --dd-update --save --directory $LOG_DIR" # "plotecray --uri \"imas:hdf5?user=public;pulse=134173;run=2326;database=TEST;version=3\" -md wall --dd-update --save --directory $LOG_DIR" # "plotedgeprofiles --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --wall --time 60 --dd-update --save --directory $LOG_DIR" -# "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --rho -md \"imas:hdf5?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active\" \"imas:hdf5?user=public;pulse=116000;run=5;database=ITER_MD;version=3#wall\" --dd-update --save --directory $LOG_DIR" +# "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" --phi -md \"imas:hdf5?user=public;pulse=111001;run=103;database=ITER_MD;version=3#pf_active\" \"imas:hdf5?user=public;pulse=116000;run=5;database=ITER_MD;version=3#wall\" --dd-update --save --directory $LOG_DIR" # "plotequilibrium --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" -md pf_active wall --dd-update --dd-update --save --directory $LOG_DIR" # "printfluxes --uri \"imas:hdf5?user=public;pulse=134174;run=117;database=ITER;version=3\" -m CLOSEST --dd-update" # "plothcd -ech 134173/101/public/MDSPLUS/TEST/3 -nbi 130012/115/public/MDSPLUS/TEST/3 -fus 130012/115/public/MDSPLUS/TEST/3 -icrh 130012/115/public/MDSPLUS/TEST/3 --dd-update --save --directory $LOG_DIR"