From 64ef94c348e2116ab5f0cbfdccc0080ad96534e6 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Fri, 27 Feb 2026 13:16:43 +0200 Subject: [PATCH 01/26] Added lumicontainer builder script --- lumi_container.sh | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 lumi_container.sh diff --git a/lumi_container.sh b/lumi_container.sh new file mode 100644 index 000000000..ebf428b09 --- /dev/null +++ b/lumi_container.sh @@ -0,0 +1,38 @@ +#!/bin/bash +#Make container.def +if [ -f ./container.def ]; then + cat << EOF > container.def + Bootstrap: docker + From: docker.io/opensuse/leap:15.6 + + %post + # Continue to install software into the container normally. + # Building: + # module load LUMI systools + # singularity build container.sif container.def + # Running: + # For example launching python (one can also launch bash for example) + # singularity exec container.sif python + + zypper -n --no-gpg-checks install python312 python312-pip python312-setuptools python312-devel git + + update-alternatives --install /usr/bin/python python /usr/bin/python3.12 3 + update-alternatives --install /usr/bin/python3 python3 /usr/bin/python3.12 3 + update-alternatives --install /usr/bin/pip pip /usr/bin/pip3.12 3 + update-alternatives --install /usr/bin/pip3 pip3 /usr/bin/pip3.12 3 + + export PYTHONNOUSERSITE=1 + pip install git+https://github.com/alhom/analysator-backends/releases/tag/v0.0.2 + pip install --editable git+https://github.com/fmihpc/analysator#egg=analysator" + EOF + + #Build the container + module load LUMI systools + singularity build container.sif container.def +else + export PYTHONNOUSERSITE=1 + singulairty exec container.def bash +fi + + + From 815de6301fd507539a654d97b25974484d957343 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Thu, 12 Mar 2026 18:49:53 +0200 Subject: [PATCH 02/26] Added way to define python version with PY_VERSION env var, defaults to 3.12 if not givne --- lumi_container.sh | 46 +++++++++++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/lumi_container.sh b/lumi_container.sh index ebf428b09..d7f72dddb 100644 --- a/lumi_container.sh +++ b/lumi_container.sh @@ -1,37 +1,45 @@ #!/bin/bash #Make container.def -if [ -f ./container.def ]; then + +if [[ "$PY_VERSION" ]]; then + PY_PACKAGE=$(echo "$PY_VERSION" | sed 's/\.//g' ) +else + PY_VERSION="3.12" + PY_PACKAGE="312" +fi + +if [ ! -f ./container.def ]; then cat << EOF > container.def - Bootstrap: docker - From: docker.io/opensuse/leap:15.6 +Bootstrap: docker +From: docker.io/opensuse/leap:15.6 - %post - # Continue to install software into the container normally. - # Building: - # module load LUMI systools - # singularity build container.sif container.def - # Running: - # For example launching python (one can also launch bash for example) - # singularity exec container.sif python +%post + # Continue to install software into the container normally. + # Building: + # module load LUMI systools + # singularity build container.sif container.def + # Running: + # For example launching python (one can also launch bash for example) + # singularity exec container.sif python - zypper -n --no-gpg-checks install python312 python312-pip python312-setuptools python312-devel git + zypper -n --no-gpg-checks install python$PY_PACKAGE python$PY_PACKAGE-pip python$PY_PACKAGE-setuptools python$PY_PACKAGE-devel git - update-alternatives --install /usr/bin/python python /usr/bin/python3.12 3 - update-alternatives --install /usr/bin/python3 python3 /usr/bin/python3.12 3 - update-alternatives --install /usr/bin/pip pip /usr/bin/pip3.12 3 - update-alternatives --install /usr/bin/pip3 pip3 /usr/bin/pip3.12 3 + update-alternatives --install /usr/bin/python python /usr/bin/python$PY_VERSION 3 + update-alternatives --install /usr/bin/python3 python3 /usr/bin/python$PY_VERSION 3 + update-alternatives --install /usr/bin/pip pip /usr/bin/pip$PY_VERSION 3 + update-alternatives --install /usr/bin/pip3 pip3 /usr/bin/pip$PY_VERSION 3 export PYTHONNOUSERSITE=1 - pip install git+https://github.com/alhom/analysator-backends/releases/tag/v0.0.2 + pip install --extra-index-url https://version.helsinki.fi/api/v4/projects/5244/packages/pypi/simple vlsvrs pip install --editable git+https://github.com/fmihpc/analysator#egg=analysator" - EOF +EOF #Build the container module load LUMI systools singularity build container.sif container.def else export PYTHONNOUSERSITE=1 - singulairty exec container.def bash + singularity exec container.def bash fi From 47abf57b5dfffbe3820ed758231da2e89e7e733a Mon Sep 17 00:00:00 2001 From: lassjsc Date: Mon, 23 Mar 2026 17:39:29 +0200 Subject: [PATCH 03/26] typo in singularity exec --- lumi_container.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100644 => 100755 lumi_container.sh diff --git a/lumi_container.sh b/lumi_container.sh old mode 100644 new mode 100755 index d7f72dddb..0cc5a9e88 --- a/lumi_container.sh +++ b/lumi_container.sh @@ -39,7 +39,7 @@ EOF singularity build container.sif container.def else export PYTHONNOUSERSITE=1 - singularity exec container.def bash + singularity exec container.sif bash fi From 0a536299507b0d02e880eaa61846f93beb7144d8 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Wed, 22 Apr 2026 11:17:21 +0300 Subject: [PATCH 04/26] Added missing " --- lumi_container.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lumi_container.sh b/lumi_container.sh index 0cc5a9e88..f6855cd20 100755 --- a/lumi_container.sh +++ b/lumi_container.sh @@ -31,7 +31,7 @@ From: docker.io/opensuse/leap:15.6 export PYTHONNOUSERSITE=1 pip install --extra-index-url https://version.helsinki.fi/api/v4/projects/5244/packages/pypi/simple vlsvrs - pip install --editable git+https://github.com/fmihpc/analysator#egg=analysator" + pip install --editable git+"https://github.com/fmihpc/analysator#egg=analysator" EOF #Build the container From d6122c751ceb56a48741ff08b2d4197d671b8a49 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Wed, 22 Apr 2026 11:19:08 +0300 Subject: [PATCH 05/26] Changed the file check such that it makes the container if the .sif (container itself) is not found, previously was checking against .def --- lumi_container.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lumi_container.sh b/lumi_container.sh index f6855cd20..9078468e4 100755 --- a/lumi_container.sh +++ b/lumi_container.sh @@ -8,7 +8,7 @@ else PY_PACKAGE="312" fi -if [ ! -f ./container.def ]; then +if [ ! -f ./container.sif ]; then cat << EOF > container.def Bootstrap: docker From: docker.io/opensuse/leap:15.6 From 9303ba5bb1357f231940f11f8d5bbaf4ff42f2c7 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 23 Apr 2026 14:43:00 +0300 Subject: [PATCH 06/26] CITATION.cff update --- CITATION.cff | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index d66eab2a3..0760688ea 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,12 +1,18 @@ cff-version: 1.2.0 message: "If you use this software, please cite it as below." authors: -- family-names: "Battarbee" - given-names: "Markus" - orcid: "https://orcid.org/0000-0001-7055-551X" - family-names: "Alho" given-names: "Markku" orcid: "https://orcid.org/0000-0001-9762-6795" +- family-names: "Battarbee" + given-names: "Markus" + orcid: "https://orcid.org/0000-0001-7055-551X" +- family-names: "Hannuksela" + given-names: "Otto" +- family-names: "Sihvonen" + given-names: "Lasse" +- family-names: "Koskela" + given-names: "Otto" - family-names: "Pfau-Kempf" given-names: "Yann" orcid: "https://orcid.org/0000-0001-5793-7070" @@ -48,6 +54,10 @@ authors: - family-names: "Turc" given-names: "Lucile" orcid: "https://orcid.org/0000-0002-7576-3251" +- family-names: "Ala-Lahti" + given-names: "Matti" +- family-names: "Kilian" + given-names: "Patrick" - family-names: "Tesema" given-names: "Fasil Kebede" orcid: "https://orcid.org/0000-0002-5978-6955" @@ -57,6 +67,8 @@ authors: - family-names: "Honkonen" given-names: "Ilja" orcid: "https://orcid.org/0000-0002-9542-5866" +- family-names: "Hoilijoki" + given-names: "Sanni" - family-names: "Brito" given-names: "Thiago" orcid: "https://orcid.org/0000-0002-2531-5848" @@ -67,8 +79,10 @@ authors: given-names: "Sameli" - family-names: "Reimi" given-names: "Johanna" +- family-names: "Greguš" + given-names: "Jakub" title: "Analysator" -version: 0.9.6 +version: 0.9.9 doi: 10.5281/zenodo.4462514 date-released: 2024-11-26 url: "https://github.com/fmihpc/analysator" From 4c0eda38482b8f5a6ffded4ebdcf631ea0a29888 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 23 Apr 2026 14:47:16 +0300 Subject: [PATCH 07/26] CITATION.cff update with orcids --- CITATION.cff | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 0760688ea..3dd47299e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -11,8 +11,6 @@ authors: given-names: "Otto" - family-names: "Sihvonen" given-names: "Lasse" -- family-names: "Koskela" - given-names: "Otto" - family-names: "Pfau-Kempf" given-names: "Yann" orcid: "https://orcid.org/0000-0001-5793-7070" @@ -56,8 +54,10 @@ authors: orcid: "https://orcid.org/0000-0002-7576-3251" - family-names: "Ala-Lahti" given-names: "Matti" + orcid: "https://orcid.org/0000-0001-9574-339X" - family-names: "Kilian" given-names: "Patrick" + orcid: "https://orcid.org/0000-0002-8906-7783" - family-names: "Tesema" given-names: "Fasil Kebede" orcid: "https://orcid.org/0000-0002-5978-6955" @@ -69,6 +69,7 @@ authors: orcid: "https://orcid.org/0000-0002-9542-5866" - family-names: "Hoilijoki" given-names: "Sanni" + orcid: "https://orcid.org/0000-0001-5709-2346" - family-names: "Brito" given-names: "Thiago" orcid: "https://orcid.org/0000-0002-2531-5848" @@ -82,7 +83,7 @@ authors: - family-names: "Greguš" given-names: "Jakub" title: "Analysator" -version: 0.9.9 +version: 0.9.8 doi: 10.5281/zenodo.4462514 date-released: 2024-11-26 url: "https://github.com/fmihpc/analysator" From 59413f8df5143528b9bd00666f5f981aaa6b910e Mon Sep 17 00:00:00 2001 From: lassejsc Date: Tue, 26 May 2026 14:53:39 +0300 Subject: [PATCH 08/26] Moving lumi_onctainer to a new folder "containerization", perhaps more can be added later --- lumi_container.sh => containerization/lumi_container.sh | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename lumi_container.sh => containerization/lumi_container.sh (100%) diff --git a/lumi_container.sh b/containerization/lumi_container.sh similarity index 100% rename from lumi_container.sh rename to containerization/lumi_container.sh From af55983f21faa7632e5ba3bc560e4ce486d6e2ac Mon Sep 17 00:00:00 2001 From: lassejsc Date: Thu, 28 May 2026 10:06:39 +0300 Subject: [PATCH 09/26] added documentation about binding folders with a friendly warning for the launch script --- containerization/lumi_container.sh | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/containerization/lumi_container.sh b/containerization/lumi_container.sh index 9078468e4..361013bc6 100755 --- a/containerization/lumi_container.sh +++ b/containerization/lumi_container.sh @@ -21,6 +21,9 @@ From: docker.io/opensuse/leap:15.6 # Running: # For example launching python (one can also launch bash for example) # singularity exec container.sif python + # to bind folders once can use --bind /opt,/data:/mnt + # this will bind /opt to /opt and /data to /mnt, alternatively one can use environment variables + # export SINGULARITY_BIND="/opt,/data:/mnt" zypper -n --no-gpg-checks install python$PY_PACKAGE python$PY_PACKAGE-pip python$PY_PACKAGE-setuptools python$PY_PACKAGE-devel git @@ -39,6 +42,9 @@ EOF singularity build container.sif container.def else export PYTHONNOUSERSITE=1 + if [[ ! $SINGULARITY_BIND ]]; then + echo "Friendly reminder: no folders binded to the container, please use SINGULARITY_BIND environment variable" + fi singularity exec container.sif bash fi From d990d354f3e4328066db7a6f5c76af8835e60da2 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Thu, 28 May 2026 10:39:40 +0300 Subject: [PATCH 10/26] Added documentation for the lumi container --- Documentation/sphinx/hpcusage.rst | 14 ++++++++++++++ Documentation/sphinx/index.rst | 1 + containerization/lumi_container.sh | 4 ++++ 3 files changed, 19 insertions(+) create mode 100644 Documentation/sphinx/hpcusage.rst diff --git a/Documentation/sphinx/hpcusage.rst b/Documentation/sphinx/hpcusage.rst new file mode 100644 index 000000000..b097889ca --- /dev/null +++ b/Documentation/sphinx/hpcusage.rst @@ -0,0 +1,14 @@ +Usage in HPC environment +================= + + +LUMI Container +------------------ +Script for building and running a container on LUMI can be found in ``containerization/lumi_container.sh``. + +Simply execute the bash script to build the container, afterwards the same script can be used for launching the container. + +Before building the container set the python version for the container export the environment variable ``PY_VERSION``, which defaults to ``PY_VERSION="3.12"``. + +For binding folders into the container, set the environment variable ``SINGULARITY_BIND`` before running the container. For example, ``export SINGULARITY_BIND="/opt,/data:/mnt"``, this will bind /opt to /opt and /data to /mnt in the container. +For more information on singularity usage, see https://docs.sylabs.io/guides/3.0/user-guide/index.html diff --git a/Documentation/sphinx/index.rst b/Documentation/sphinx/index.rst index fca78f67d..3ca5b28c5 100644 --- a/Documentation/sphinx/index.rst +++ b/Documentation/sphinx/index.rst @@ -24,6 +24,7 @@ analysator documentation analysator analysator_exercises analysator_supported + hpcusage .. toctree:: :maxdepth: 2 diff --git a/containerization/lumi_container.sh b/containerization/lumi_container.sh index 361013bc6..729488daf 100755 --- a/containerization/lumi_container.sh +++ b/containerization/lumi_container.sh @@ -1,5 +1,9 @@ #!/bin/bash #Make container.def +#Usage: +# Simply run this bash script to build the container, be sure to export PY_VERSION to set the python version desired, for example PY_VERSION=3.12 +# After building one can run the script again to launch the container, alternatively call the .sif file directly with singularity +# For binding folders export SINGULARITY_BIND see for example below in # Running: if [[ "$PY_VERSION" ]]; then PY_PACKAGE=$(echo "$PY_VERSION" | sed 's/\.//g' ) From b27d10648454c35c614f4fb9b68c0bdce001a777 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Thu, 28 May 2026 10:48:40 +0300 Subject: [PATCH 11/26] Fixed the whatever sphinx complained about the underline --- Documentation/sphinx/hpcusage.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Documentation/sphinx/hpcusage.rst b/Documentation/sphinx/hpcusage.rst index b097889ca..e6447d175 100644 --- a/Documentation/sphinx/hpcusage.rst +++ b/Documentation/sphinx/hpcusage.rst @@ -1,9 +1,9 @@ Usage in HPC environment -================= +======================== LUMI Container ------------------- +-------------- Script for building and running a container on LUMI can be found in ``containerization/lumi_container.sh``. Simply execute the bash script to build the container, afterwards the same script can be used for launching the container. From 731dd1d5be932e83fca1c7d7b04edd22aab2cb92 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 28 May 2026 11:44:18 +0300 Subject: [PATCH 12/26] logging config for obliqueshock.py --- scripts/obliqueshock.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/scripts/obliqueshock.py b/scripts/obliqueshock.py index 60cf0b9c3..dbe7c152b 100644 --- a/scripts/obliqueshock.py +++ b/scripts/obliqueshock.py @@ -7,6 +7,10 @@ import math import scipy.optimize import logging +logging.basicConfig( + level=logging.INFO, # Increase level to disable messages + format="%(message)s" +) mu0 = 4*math.pi*1.e-7 mp = 1.67e-27 From 0517a48214abc9babf45af98d12d0b94c42fecc2 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 4 Sep 2025 16:17:22 +0300 Subject: [PATCH 13/26] prototype to ditch Delaunay alpha-hulls when we actually have usable connectivity information --- scripts/magnetopause.py | 87 +++++++++++++++++++++++++++++++++++------ scripts/regions.py | 30 ++++++++++++-- 2 files changed, 100 insertions(+), 17 deletions(-) diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index 9f8dd66ea..1ea51a7c9 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -67,7 +67,7 @@ def write_SDF_to_file(SDF, datafilen, outfilen): -def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.4, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method +def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds=None, return_surface=True, return_SDF=True, SDF_points=None, Delaunay_alpha=None, beta_star_range=[0.0, 0.5], method_args={}): # TODO: separate streamline suface and vtkDelaunay3d surface in streamline method """Finds the magnetopause using the specified method. Surface is constructed using vtk's Delaunay3d triangulation which results in a convex hull if no Delaunay_alpha is given. Returns vtk.vtkDataSetSurfaceFilter object and/or signed distances (negative -> inside magnetopause) (=SDF) to all cells Note that using alpha for Delaunay might make SDF different from expected inside the magnetosphere, especially if surface is constructed with points not everywhere in the magnetosphere (e.g. beta* 0.4-0.5) or if simulation grid size is larger than alpha @@ -119,19 +119,80 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= elif method == "beta_star_with_connectivity": # magnetopause from beta_star, with connectivity if possible betastar_region = regions.treshold_mask(f.read_variable("vg_beta_star"), beta_star_range) - try: - connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines - magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) - contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1]) - np.save("pointcloud.npy", contour_coords) - except: - logging.warning("using field line connectivity for magnetosphere did not work, using only beta*") - #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause - mpause_flags = np.where(betastar_region==1, 1, 0) - contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0]) + if True: + vtkreader = pt.vlsvfile.VlsvVtkReader() + vtkreader.SetFileName(f.file_name) + vtkreader.Update() + # vars =vtkreader.findVariablesFromVlsv(getReducers=False) + # add here more/other datareducer outputs for downstream use if needed + vars = ["proton/vg_beta_star", "vg_connection","cellid"] + + for var in [v for v in vars if ("vg_" in v.lower()) or (v.lower() == "cellid")]: + vtkreader.addArrayFromVlsv(var) + + vtkreader.Modified() + vtkreader.Update() + dataport = vtkreader.GetOutputPort() + + dual = vtk.vtkHyperTreeGridToDualGrid() + dual.SetInputConnection(dataport) + vtkreader.Update() + dual.Update() + + data = dual.GetOutputDataObject(0) + renamer = vtk.vtkArrayRename() + renamer.SetPointArrayName("proton/vg_beta_star", "vg_beta_star") + renamer.SetInputData(data) + renamer.Update() + + data = renamer.GetOutputDataObject(0) + # print(data) + + newarr = vtk.vtkDoubleArray() + newarr.SetName("BL") + newarr.SetNumberOfComponents(1) # This will be a scalar results + newarr.SetNumberOfTuples(data.GetPointData().GetNumberOfTuples()) # Set how many values to allocate + data.GetPointData().AddArray(newarr) # add the array to the dataset + + + # Setting up the actual calulcation + calcBL = vtk.vtkArrayCalculator() + calcBL.SetInputData(data) # This data we ingest + # We need to specify which arrays we will use to derive the value + calcBL.AddScalarArrayName('vg_beta_star') + calcBL.AddScalarArrayName('vg_connection') + calcBL.SetReplaceInvalidValues(True) + calcBL.SetReplacementValue(np.nan) + # The actual function to calculate + calcBL.SetFunction("vg_beta_star*min(vg_connection,1)") + # Set result array name + calcBL.SetResultArrayName("BL") + calcBL.Update() + data = calcBL.GetOutputDataObject(0) + + threshold0 = vtk.vtkThreshold() + threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "BL") + threshold0.SetInputData(data) + threshold0.SetLowerThreshold(beta_star_range[0]) + threshold0.SetUpperThreshold(beta_star_range[1]) + threshold0.Update() + + vtkSurface, SDF = regions.vtkSDF(query_points, threshold0.GetOutputDataObject(0)) + + else: + try: + connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines + magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) + contour_coords = f.get_cell_coordinates(cellids[magnetosphere_proper==1]) + np.save("pointcloud.npy", contour_coords) + except: + logging.warning("using field line connectivity for magnetosphere did not work, using only beta*") + #condition_dict = {"beta_star": [0.5, 0.6]} # FIC: [0.4, 0.5]) # EGE: [0.9, 1.0]) # max 0.6 in FHA to not take flyaways from outside magnetopause + mpause_flags = np.where(betastar_region==1, 1, 0) + contour_coords = f.get_cell_coordinates(cellids[mpause_flags!=0]) - # make a convex hull surface with vtk's Delaunay - vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha) + # make a convex hull surface with vtk's Delaunay + vtkSurface, SDF = regions.vtkDelaunay3d_SDF(query_points, contour_coords, Delaunay_alpha) #elif method == "beta_star_with_fieldlines": # either incredibly slow or does not work, don't use without fixing #TODO proprer measure of actual field line backwall point # # magnetopause from beta_star, with field lines connecting to back wall if possible diff --git a/scripts/regions.py b/scripts/regions.py index 7b7895f77..807da4569 100644 --- a/scripts/regions.py +++ b/scripts/regions.py @@ -67,6 +67,27 @@ def vtkDelaunay3d_SDF(query_points, coordinates, alpha=None): return surface, convexhull_sdf +def vtkSDF(query_points, dualgrid): + ''' Obtain the SDF in relation to a waterproof volumetric grid "dualgrid" + ''' + + # print(dualgrid) + surface = vtk.vtkDataSetSurfaceFilter() + surface.SetInputData(dualgrid) + surface.Update() + + # print(surface.GetOutput()) + + + implicitPolyDataDistance = vtk.vtkImplicitPolyDataDistance() + implicitPolyDataDistance.SetInput(surface.GetOutput()) + import sys + convexhull_sdf = np.zeros(len(query_points)) + for i,coord in enumerate(query_points): + convexhull_sdf[i] = implicitPolyDataDistance.EvaluateFunction(coord) + # sys.exit() + + return surface.GetOutputDataObject(0), convexhull_sdf def treshold_mask(data_array, value): @@ -347,7 +368,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo if magnetopause_kwargs: __, magnetopause_SDF = magnetopause.magnetopause(datafile, **magnetopause_kwargs) else: - __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=None) # default magnetopause: beta*+ B connectivity convex hull + __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=2e6, ) # default magnetopause: beta*+ B connectivity convex hull write_flags(writer, magnetopause_SDF, 'SDF_magnetopause') write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause") @@ -495,10 +516,11 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo def main(): - datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/EGE/bulk/bulk.0002000.vlsv" - outfilen = "EGE_regions_t2000.vlsv" + fileid = 1000 + datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.{:07d}.vlsv".format(fileid) + outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/prototyping/FID_mpause_{:07d}.vlsv".format(fileid) - RegionFlags(datafile, outfilen, regions=["all"]) + RegionFlags(datafile, outfilen, regions=["magnetopause"]) if __name__ == "__main__": From e844c9ef6c18f9341153bcb0ea6dd15585f24f60 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 4 Sep 2025 17:12:29 +0300 Subject: [PATCH 14/26] vg_-prefix to outputs --- scripts/regions.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/scripts/regions.py b/scripts/regions.py index 807da4569..d6da1fd4e 100644 --- a/scripts/regions.py +++ b/scripts/regions.py @@ -369,8 +369,8 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo __, magnetopause_SDF = magnetopause.magnetopause(datafile, **magnetopause_kwargs) else: __, magnetopause_SDF = magnetopause.magnetopause(datafile, method="beta_star_with_connectivity", Delaunay_alpha=2e6, ) # default magnetopause: beta*+ B connectivity convex hull - write_flags(writer, magnetopause_SDF, 'SDF_magnetopause') - write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause") + write_flags(writer, magnetopause_SDF, 'vg_SDF_magnetopause') + write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "vg_flag_magnetopause") # save some magnetopause values for later magnetopause_density = np.mean(variables["density"][np.abs(magnetopause_SDF) < 5e6]) @@ -382,7 +382,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo # magnetosphere from magnetopause SDF if "magnetosphere" in regions: magnetosphere = np.where(magnetopause_SDF<0, 1, 0) - write_flags(writer, magnetosphere, 'flag_magnetosphere') + write_flags(writer, magnetosphere, 'vg_flag_magnetosphere') ## BOW SHOCK ## #TODO: similar kwargs system as magnetopause? @@ -392,8 +392,8 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo bowshock = bowshock_SDF(f, variables, all_points, own_condition_dict=region_conditions["bowshock"]) else: bowshock = bowshock_SDF(f, variables, all_points) # default upstream rho method, might fail with foreshock - write_flags(writer, bowshock, 'SDF_bowshock') - write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "flag_bowshock") + write_flags(writer, bowshock, 'vg_SDF_bowshock') + write_flags(writer, np.where(np.abs(bowshock) < 5e6, 1, 0), "vg_flag_bowshock") # magnetosphere+magnetosheath -area inside_bowshock = np.where(bowshock<0, 1, 0) @@ -402,7 +402,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo if "magnetosheath" in regions: # magnetosheath from bow shock-magnetosphere difference magnetosheath_flags = np.where((inside_bowshock & 1-magnetosphere), 1, 0) - write_flags(writer, magnetosheath_flags, 'flag_magnetosheath') + write_flags(writer, magnetosheath_flags, 'vg_flag_magnetosheath') # save magnetosheath density and temperature for further use #magnetosheath_density = np.mean(variables["density"][magnetosheath_flags == 1]) @@ -413,7 +413,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo ## UPSTREAM ## if "upstream" in regions: # upstream from !bowshock - write_flags(writer, 1-inside_bowshock, 'flag_upstream') + write_flags(writer, 1-inside_bowshock, 'vg_flag_upstream') #write_flags(writer, inside_bowshock, 'flag_inside_bowshock') @@ -444,7 +444,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } cusp_flags = make_region_flags(variables, cusp_conditions, flag_type=region_flag_type, mask=mask_inMagnetosphere) - write_flags(writer, cusp_flags, 'flag_cusps', mask_inMagnetosphere) + write_flags(writer, cusp_flags, 'vg_flag_cusps', mask_inMagnetosphere) # magnetotail lobes @@ -474,11 +474,11 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } lobe_S_flags = make_region_flags(variables, lobe_S_conditions, flag_type=region_flag_type) - write_flags(writer, lobe_N_flags, 'flag_lobe_N') - write_flags(writer, lobe_S_flags, 'flag_lobe_S') + write_flags(writer, lobe_N_flags, 'vg_flag_lobe_N') + write_flags(writer, lobe_S_flags, 'vg_flag_lobe_S') lobes_flags = make_region_flags(variables, lobes_conditions, flag_type=region_flag_type) - write_flags(writer, lobes_flags, 'flag_lobes') + write_flags(writer, lobes_flags, 'vg_flag_lobes') # lobe density from median densities? @@ -498,7 +498,7 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo } central_plasma_sheet_flags = make_region_flags(variables, central_plasma_sheet_conditions,flag_type=region_flag_type, mask=mask_inMagnetosphere) - write_flags(writer, central_plasma_sheet_flags, 'flag_central_plasma_sheet', mask_inMagnetosphere) + write_flags(writer, central_plasma_sheet_flags, 'vg_flag_central_plasma_sheet', mask_inMagnetosphere) ## Other boundary layers, PSBL sometimes works @@ -516,9 +516,14 @@ def errormsg(varstr): logging.warning("{} could not be read, will be ignored".fo def main(): - fileid = 1000 + import sys + + try: + fileid = int(sys.argv[1]) + except Exception as e: + print("Need fileid") datafile = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/bulk1/bulk1.{:07d}.vlsv".format(fileid) - outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/prototyping/FID_mpause_{:07d}.vlsv".format(fileid) + outfilen = "/wrk-vakka/group/spacephysics/vlasiator/3D/FID/postprocessing/magnetopause_sdfs/FID_mpause_{:07d}.vlsv".format(fileid) RegionFlags(datafile, outfilen, regions=["magnetopause"]) From 43bb0126bc6215955abde4fda78bc0bf0e0e2c7c Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Wed, 17 Sep 2025 15:32:31 +0300 Subject: [PATCH 15/26] Add variable cache injection function --- analysator/vlsvfile/vlsvreader.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/analysator/vlsvfile/vlsvreader.py b/analysator/vlsvfile/vlsvreader.py index 915bbc185..d0543edd8 100644 --- a/analysator/vlsvfile/vlsvreader.py +++ b/analysator/vlsvfile/vlsvreader.py @@ -2491,6 +2491,19 @@ def read_variable_to_cache(self, name, operator="pass"): # self.__read_fileindex_for_cellid() return self.__variable_cache[(name,operator)] + def add_cached_variable(self, data, name, operator="pass"): + ''' Add a variable to cache, for the whole grid and after applying + operator. Works by injecting a (variable,operator) key and the data to cache. + + :param data: data array - make sure it is of the correct shape and size. + :param name: Name of the variable (or datareducer) + :param operator: Datareduction operator. "pass" does no operation on data. + + ''' + + # add data to dict, use a tuple of (name,operator) as the key [tuples are immutable and hashable] + self.__variable_cache[(name,operator)] = data + def read_variable(self, name, cellids=-1,operator="pass"): ''' Read variables from the open vlsv file. Arguments: From 238f4013faff1e6481a0213a080df78f8346d880 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 18 Sep 2025 13:46:13 +0300 Subject: [PATCH 16/26] Region VTK pipeline cleanup, temp variable cache injection for vlsvreader and reader-based init for VlsvVtkInterface --- analysator/vlsvfile/vlsvvtkinterface.py | 18 +++++++--- scripts/magnetopause.py | 44 +++++++------------------ 2 files changed, 24 insertions(+), 38 deletions(-) diff --git a/analysator/vlsvfile/vlsvvtkinterface.py b/analysator/vlsvfile/vlsvvtkinterface.py index e0e529534..2e2d010b4 100644 --- a/analysator/vlsvfile/vlsvvtkinterface.py +++ b/analysator/vlsvfile/vlsvvtkinterface.py @@ -78,13 +78,21 @@ def __init__(self): def SetFileName(self, filename): if filename != self.__FileName: - self.Modified() - self.__FileName = filename - if self.__FileName is not None: - self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) - fn = os.path.basename(self.__FileName) + if filename is not None: + if self.__reader is None: + self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) + else: + raise ValueError("Tried to change an existing reader ("+self.__reader.file_name+") to " + filename) + self.Modified() + self.__FileName = filename self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl") + def SetReader(self, reader): + self.__reader = reader + self.Modified() + self.__FileName = reader.file_name + self.__metafile = os.path.join(self.__reader.get_cache_folder(),"vlsvvtkcache.pkl") + def GetFileName(self): return self.__FileName diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index 1ea51a7c9..b466b53af 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -118,14 +118,19 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= elif method == "beta_star_with_connectivity": # magnetopause from beta_star, with connectivity if possible - betastar_region = regions.treshold_mask(f.read_variable("vg_beta_star"), beta_star_range) + vg_beta_star =f.read_variable("vg_beta_star") + vg_conn = f.read_variable("vg_connection") + vg_classifier = vg_beta_star*np.min(vg_conn,1) # if closed-closed fieldlines, set var to zero if True: + vtkreader = pt.vlsvfile.VlsvVtkReader() - vtkreader.SetFileName(f.file_name) + vtkreader.SetReader(f) + f.add_cached_variable(vg_classifier, "vg_betastar_classifier") vtkreader.Update() + # vars =vtkreader.findVariablesFromVlsv(getReducers=False) # add here more/other datareducer outputs for downstream use if needed - vars = ["proton/vg_beta_star", "vg_connection","cellid"] + vars = ["vg_betastar_classifier", "cellid"] for var in [v for v in vars if ("vg_" in v.lower()) or (v.lower() == "cellid")]: vtkreader.addArrayFromVlsv(var) @@ -140,38 +145,9 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= dual.Update() data = dual.GetOutputDataObject(0) - renamer = vtk.vtkArrayRename() - renamer.SetPointArrayName("proton/vg_beta_star", "vg_beta_star") - renamer.SetInputData(data) - renamer.Update() - - data = renamer.GetOutputDataObject(0) - # print(data) - - newarr = vtk.vtkDoubleArray() - newarr.SetName("BL") - newarr.SetNumberOfComponents(1) # This will be a scalar results - newarr.SetNumberOfTuples(data.GetPointData().GetNumberOfTuples()) # Set how many values to allocate - data.GetPointData().AddArray(newarr) # add the array to the dataset - - - # Setting up the actual calulcation - calcBL = vtk.vtkArrayCalculator() - calcBL.SetInputData(data) # This data we ingest - # We need to specify which arrays we will use to derive the value - calcBL.AddScalarArrayName('vg_beta_star') - calcBL.AddScalarArrayName('vg_connection') - calcBL.SetReplaceInvalidValues(True) - calcBL.SetReplacementValue(np.nan) - # The actual function to calculate - calcBL.SetFunction("vg_beta_star*min(vg_connection,1)") - # Set result array name - calcBL.SetResultArrayName("BL") - calcBL.Update() - data = calcBL.GetOutputDataObject(0) threshold0 = vtk.vtkThreshold() - threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "BL") + threshold0.SetInputArrayToProcess(0,0,0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "vg_betastar_classifier") threshold0.SetInputData(data) threshold0.SetLowerThreshold(beta_star_range[0]) threshold0.SetUpperThreshold(beta_star_range[1]) @@ -180,6 +156,8 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= vtkSurface, SDF = regions.vtkSDF(query_points, threshold0.GetOutputDataObject(0)) else: + betastar_region = regions.treshold_mask(vg_classifier, beta_star_range) + try: connectivity_region = regions.treshold_mask(f.read_variable("vg_connection"), 0) # closed-closed magnetic field lines magnetosphere_proper = np.where((connectivity_region==1) | (betastar_region==1), 1, 0) From 5869fb2cdf1c0519302af3e414716da2c39da59d Mon Sep 17 00:00:00 2001 From: Leevi Loikkanen Date: Tue, 2 Jun 2026 17:34:55 +0300 Subject: [PATCH 17/26] corrected VlsvTInterpolator cached variable name (#481) --- analysator/calculations/timeevolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/timeevolution.py b/analysator/calculations/timeevolution.py index 3742c5bd1..62c8b9892 100644 --- a/analysator/calculations/timeevolution.py +++ b/analysator/calculations/timeevolution.py @@ -296,7 +296,7 @@ def __call__(self, t, coordinates, variable, operator='pass', method="linear", e self.activeReaders['low'] = self.readers[max(rti-1,0)] if self.activeReaders['hi']: - if lower_t == self.activeReaders['hi'].read_parameter('time'): + if upper_t == self.activeReaders['hi'].read_parameter('time'): pass else: self.activeReaders['hi'] = self.readers[min(rti,len(self.ts)-1)] From 67cf19ee7d991d08b8eb315620ffd2e750a5a2e6 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 4 Jun 2026 10:39:26 +0300 Subject: [PATCH 18/26] Strip filename string before creating linked readers --- analysator/vlsvfile/vlsvreader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/vlsvfile/vlsvreader.py b/analysator/vlsvfile/vlsvreader.py index 915bbc185..d7a9e6c6a 100644 --- a/analysator/vlsvfile/vlsvreader.py +++ b/analysator/vlsvfile/vlsvreader.py @@ -549,7 +549,7 @@ def add_linked_reader(self, fname): for reader in self.__linked_readers: if fname == reader.file_name: return - self.__linked_readers.add(VlsvReader(fname)) + self.__linked_readers.add(VlsvReader(fname.strip())) else: logging.warning("Could not link "+fname+" (path does not exist)") From ac845c2a5626f8ccac0fcd5d993f4d69046e45e8 Mon Sep 17 00:00:00 2001 From: jreimi Date: Mon, 8 Jun 2026 10:44:37 +0300 Subject: [PATCH 19/26] Added a quick fix for vtk to file writer func to accept both connection and dataobject, bug and typo fix --- scripts/magnetopause.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index b466b53af..1c77f8b2e 100644 --- a/scripts/magnetopause.py +++ b/scripts/magnetopause.py @@ -51,9 +51,16 @@ def write_vtk_surface_to_file(vtkSurface, outfilen): - writer = vtk.vtkXMLPolyDataWriter() - writer.SetInputConnection(vtkSurface.GetOutputPort()) + + writer = vtk.vtkXMLPolyDataWriter() writer.SetFileName(outfilen) + + try: + writer.SetInputConnection(vtkSurface.GetOutputPort()) + except: + writer.SetInputData(vtkSurface) + + writer.Write() logging.info("wrote ", outfilen) @@ -79,7 +86,7 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= :kwarg return_SDF: True/False, return array of distances in m to SDF_points in point input order, negative distance inside the surface :kwarg SDF_points: optionally give array of own points to calculate signed distances to. If not given, distances will be to cell centres in the order of f.read_variable("CellID") output :kwarg Delaunay_alpha: alpha (float) to give to vtkDelaunay3d, None -> convex hull, alpha=__: surface egdes longer than __ will be excluded - :kwarg beta_star_range: [min, max] treshold rage to use with methods "beta_star" and "beta_star_with_connectivity" + :kwarg beta_star_range: [min, max] treshold range to use with methods "beta_star" and "beta_star_with_connectivity" :kwarg method_args: dict of keyword arguments to be passed down to external functions (for streamlines and shue) :returns: vtkDataSetSurfaceFilter object of convex hull or alpha shape if return_surface=True, signed distance field of convex hull or alpha shape of magnetopause if return_SDF=True """ @@ -120,7 +127,7 @@ def magnetopause(datafilen, method="beta_star_with_connectivity", own_tresholds= # magnetopause from beta_star, with connectivity if possible vg_beta_star =f.read_variable("vg_beta_star") vg_conn = f.read_variable("vg_connection") - vg_classifier = vg_beta_star*np.min(vg_conn,1) # if closed-closed fieldlines, set var to zero + vg_classifier = vg_beta_star*np.minimum(vg_conn,1) # if closed-closed fieldlines, set var to zero if True: vtkreader = pt.vlsvfile.VlsvVtkReader() From 2ef557f10c7a93efb6d0e87286bb944fbca9ab6b Mon Sep 17 00:00:00 2001 From: JonasSuni Date: Mon, 15 Jun 2026 15:18:59 +0300 Subject: [PATCH 20/26] Import mp finding functions to calculations --- analysator/calculations/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/analysator/calculations/__init__.py b/analysator/calculations/__init__.py index 1aea6eaac..011cae3cc 100644 --- a/analysator/calculations/__init__.py +++ b/analysator/calculations/__init__.py @@ -64,3 +64,5 @@ from .null_lines import LMN_null_lines_FOTE from .interpolator_amr import AMRInterpolator, supported_amr_interpolators from .virtual_observations import spacecraft_to_simulation_frame,simulation_to_spacecraft_frame,simulation_to_observation_frame +from .magnetopause_sw_streamline_2d import find_magnetopause_sw_streamline_2d +from .magnetopause_sw_streamline_3d import find_magnetopause_sw_streamline_3d From 6e4c0b4c474998279309f2d24eb3f993eb874515 Mon Sep 17 00:00:00 2001 From: jreimi Date: Mon, 15 Jun 2026 17:14:19 +0300 Subject: [PATCH 21/26] magnetopause streamline creation condition fix --- analysator/calculations/magnetopause_sw_streamline_3d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/magnetopause_sw_streamline_3d.py b/analysator/calculations/magnetopause_sw_streamline_3d.py index 1cab1e051..8483db09a 100644 --- a/analysator/calculations/magnetopause_sw_streamline_3d.py +++ b/analysator/calculations/magnetopause_sw_streamline_3d.py @@ -245,7 +245,7 @@ def make_streamlines(vlsvfile, streamline_seeds=None, seeds_n=25, seeds_x0=20*63 f = pt.vlsvfile.VlsvReader(file_name=vlsvfile) # Create streamline starting points if needed - if streamline_seeds == None: + if not np.any(streamline_seeds): streamline_seeds = np.zeros([seeds_n**2, 3]) t = np.linspace(seeds_range[0], seeds_range[1], seeds_n) From 6d2622892db7ee82fcac5bf9f7d95f63c8fac6b6 Mon Sep 17 00:00:00 2001 From: lassejsc Date: Tue, 16 Jun 2026 11:38:05 +0300 Subject: [PATCH 22/26] Changed to --constraint="amd" to eliminate slow nodes --- testpackage/run_compare.sh | 2 +- testpackage/run_testpackage_generate_verf_set.sh | 2 +- testpackage/run_testpackage_workflow.sh | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/testpackage/run_compare.sh b/testpackage/run_compare.sh index 574a86cba..c0ce387e8 100755 --- a/testpackage/run_compare.sh +++ b/testpackage/run_compare.sh @@ -1,7 +1,7 @@ #!/bin/bash -l #SBATCH -t 00:30:00 #SBATCH -J analysator_testpackage_compare -#SBATCH --constraint="ukko|carrington" +#SBATCH --constraint="amd" #SBATCH -p short #SBATCH -n 1 #SBATCH --array=1-10 diff --git a/testpackage/run_testpackage_generate_verf_set.sh b/testpackage/run_testpackage_generate_verf_set.sh index c8e05a51f..f2f6e51c9 100644 --- a/testpackage/run_testpackage_generate_verf_set.sh +++ b/testpackage/run_testpackage_generate_verf_set.sh @@ -4,7 +4,7 @@ #SBATCH -p short #SBATCH -n 1 #SBATCH --array=0-10 -#SBATCH --constraint="carrington|ukko" +#SBATCH --constraint="amd" #SBATCH --no-requeue #SBATCH --mem-per-cpu=16000 diff --git a/testpackage/run_testpackage_workflow.sh b/testpackage/run_testpackage_workflow.sh index f5d42e737..a05cd73cd 100755 --- a/testpackage/run_testpackage_workflow.sh +++ b/testpackage/run_testpackage_workflow.sh @@ -1,7 +1,7 @@ #!/bin/bash -l #SBATCH -t 01:00:00 #SBATCH -J analysator_testpackage -#SBATCH --constraint="carrington|ukko" +#SBATCH --constraint="amd" #SBATCH -p short #SBATCH -n 1 #SBATCH --array=1-14 From a2542fb52c747163e6ea9cf42fb95bbd7e313c27 Mon Sep 17 00:00:00 2001 From: Jonas Suni <38424612+JonasSuni@users.noreply.github.com> Date: Wed, 17 Jun 2026 14:58:11 +0300 Subject: [PATCH 23/26] Fixed logic for dayside and nightside x_points --- analysator/calculations/magnetopause_sw_streamline_3d.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/analysator/calculations/magnetopause_sw_streamline_3d.py b/analysator/calculations/magnetopause_sw_streamline_3d.py index 8483db09a..8c46162ec 100644 --- a/analysator/calculations/magnetopause_sw_streamline_3d.py +++ b/analysator/calculations/magnetopause_sw_streamline_3d.py @@ -301,7 +301,7 @@ def make_magnetopause(streams, end_x=-15*6371000, x_point_n=50, sector_n=36, ign subsolar_x = np.partition(x_axis_points[:,0], ignore)[ignore] # take the nth point as subsolar point # divide the x point numbers between x > 0 (radial) an x < 0 (yz-planes) by ratio - dayside_x_point_n = int((subsolar_x/np.abs(end_x))*x_point_n) + dayside_x_point_n = int((subsolar_x/np.abs(subsolar_x - end_x))*x_point_n) ### dayside magnetopause ### # for x > 0, look for magnetopause radially @@ -355,11 +355,13 @@ def grid_mid_point(theta_idx, phi_idx): dayside_magnetopause[ring_idx] = ring_points - ### x < 0 magnetopause ### # rest: look for magnetopause in yz-planes ## define points in the x axis where to find magnetopause points on the yz-plane - x_points = np.linspace(0.0, end_x, x_point_n-dayside_x_point_n) + if end_x < 0: + x_points = np.linspace(0.0, end_x, int((np.abs(end_x)/np.abs(subsolar_x - end_x))*x_point_n)) + else: + x_points = np.array([],dtype=float) ## interpolate more exact points for streamlines at exery x_point new_streampoints = np.zeros((len(x_points), len(streams), 2)) # new array for keeping interpolated streamlines in form new_streampoints[x_point, streamline, y and z -coordinates] From edd8213ad343bf290bab293483bb12a075477e2b Mon Sep 17 00:00:00 2001 From: Jakub Gregus Date: Sat, 20 Jun 2026 12:55:02 +0300 Subject: [PATCH 24/26] vlsvReader arg takes str filename and already initialized obj, updated np.histogram arguments --- analysator/calculations/spectra.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/analysator/calculations/spectra.py b/analysator/calculations/spectra.py index 93f144d4b..1189e2fc4 100644 --- a/analysator/calculations/spectra.py +++ b/analysator/calculations/spectra.py @@ -41,7 +41,12 @@ def get_spectrum_energy(vlsvReader, EkinBinEdges = np.logspace(np.log10(EMin),np.log10(EMax),nBins+1, endpoint=True) dE = EkinBinEdges[1:] - EkinBinEdges[:-1] - vlsvReader = analysator.vlsvfile.VlsvReader(vlsvReader) + + # initialize a new vlsvReader if passed a str filename of a .vlsv + if isinstance(vlsvReader,str): + vlsvReader = analysator.vlsvfile.VlsvReader(vlsvReader) + # else keep the already initialized vlsvReader if obj supplied + # check if velocity space exists in this cell if not vlsvReader.check_variable("moments"): # restart files have VDFs everywhere if vlsvReader.check_variable("fSaved"): @@ -124,7 +129,7 @@ def get_spectrum_energy(vlsvReader, #Ekin[Ekin > max(EkinBinEdges)] = max(EkinBinEdges) # compute histogram - (nhist,edges) = np.histogram(Ekin,bins=EkinBinEdges,weights=fw,normed=0) + (nhist,edges) = np.histogram(Ekin,bins=EkinBinEdges,weights=fw) if (bindifferential): # finish differential flux per d(eV) nhist = np.divide(nhist,dE) From eb60f78c8ba48adaff8f32d173b3ca2bf4a4ee1e Mon Sep 17 00:00:00 2001 From: Jakub Gregus Date: Sat, 20 Jun 2026 12:59:32 +0300 Subject: [PATCH 25/26] updated np.histogram arguments -> unnormalized hist --- analysator/calculations/spectra.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysator/calculations/spectra.py b/analysator/calculations/spectra.py index 1189e2fc4..70215285f 100644 --- a/analysator/calculations/spectra.py +++ b/analysator/calculations/spectra.py @@ -129,7 +129,7 @@ def get_spectrum_energy(vlsvReader, #Ekin[Ekin > max(EkinBinEdges)] = max(EkinBinEdges) # compute histogram - (nhist,edges) = np.histogram(Ekin,bins=EkinBinEdges,weights=fw) + (nhist,edges) = np.histogram(Ekin,bins=EkinBinEdges,weights=fw,density=False) if (bindifferential): # finish differential flux per d(eV) nhist = np.divide(nhist,dE) From 882de524df5cf32eecbc3930b7a310499236a9bc Mon Sep 17 00:00:00 2001 From: Jakub Gregus Date: Tue, 23 Jun 2026 14:44:30 +0300 Subject: [PATCH 26/26] enforce str or VlsvReader --- analysator/calculations/spectra.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/analysator/calculations/spectra.py b/analysator/calculations/spectra.py index 70215285f..9145cb9f9 100644 --- a/analysator/calculations/spectra.py +++ b/analysator/calculations/spectra.py @@ -45,7 +45,10 @@ def get_spectrum_energy(vlsvReader, # initialize a new vlsvReader if passed a str filename of a .vlsv if isinstance(vlsvReader,str): vlsvReader = analysator.vlsvfile.VlsvReader(vlsvReader) - # else keep the already initialized vlsvReader if obj supplied + # ensure it is a VlsvReader obj instance otherwise + elif not isinstance(vlsvReader,analysator.vlsvfile.VlsvReader): + raise TypeError("Supplied vlsvReader must be either a VlsvReader obj instance or a filename str") + # else: keep the already initialized VlsvReader obj # check if velocity space exists in this cell if not vlsvReader.check_variable("moments"): # restart files have VDFs everywhere @@ -156,7 +159,13 @@ def get_spectrum_alongaxis_vel(vlsvReader, bindifferential=False, # weigh by d(velocity) restart=True): - vlsvReader = analysator.vlsvfile.VlsvReader(vlsvReader) + # initialize a new vlsvReader if passed a str filename of a .vlsv + if isinstance(vlsvReader,str): + vlsvReader = analysator.vlsvfile.VlsvReader(vlsvReader) + # ensure it is a VlsvReader obj instance otherwise + elif not isinstance(vlsvReader,analysator.vlsvfile.VlsvReader): + raise TypeError("Supplied vlsvReader must be either a VlsvReader obj instance or a filename str") + # else: keep the already initialized VlsvReader obj if vectorVar is not None and vector is None: if vlsvReader.check_variable(vectorVar): @@ -235,7 +244,7 @@ def get_spectrum_alongaxis_vel(vlsvReader, latex=r'$f(\vec{r},v)$' weight = 'particles' - (nhist,edges) = np.histogram(Vproj,bins=VBinEdges,weights=fw,normed=0) + (nhist,edges) = np.histogram(Vproj,bins=VBinEdges,weights=fw,density=False) # normalization dv = abs(VBinEdges[1:] - VBinEdges[:-1]) if (differential): # differential flux per [m/s]