Roll dev forward with master commits

CITATION.cff

@@ -1,12 +1,16 @@
 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: "Pfau-Kempf"
   given-names: "Yann"
   orcid: "https://orcid.org/0000-0001-5793-7070"
@@ -48,6 +52,12 @@ authors:
 - family-names: "Turc"
   given-names: "Lucile"
   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"
@@ -57,6 +67,9 @@ authors:
 - family-names: "Honkonen"
   given-names: "Ilja"
   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"
@@ -67,8 +80,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.8
 doi: 10.5281/zenodo.4462514
 date-released: 2024-11-26
 url: "https://github.com/fmihpc/analysator"

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

Documentation/sphinx/index.rst

@@ -24,6 +24,7 @@ analysator documentation
    analysator
    analysator_exercises
    analysator_supported
+   hpcusage
 
 .. toctree::
    :maxdepth: 2 

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

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)
@@ -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] 

analysator/calculations/spectra.py

@@ -41,7 +41,15 @@ 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)
+   # 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
       if vlsvReader.check_variable("fSaved"):
@@ -124,7 +132,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,density=False)
 
    if (bindifferential): # finish differential flux per d(eV)
       nhist = np.divide(nhist,dE)
@@ -151,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):
@@ -230,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]

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)]

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)")
 
@@ -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:

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
 

containerization/lumi_container.sh

@@ -0,0 +1,56 @@
+#!/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' )
+else
+  PY_VERSION="3.12"
+  PY_PACKAGE="312"
+fi
+
+if [ ! -f ./container.sif ]; 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
+       # 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
+
+        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 --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
+
+  #Build the container
+  module load LUMI systools
+  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
+
+
+

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)
 
@@ -67,7 +74,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
@@ -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
     """
@@ -118,20 +125,59 @@ 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])
+        vg_beta_star =f.read_variable("vg_beta_star")
+        vg_conn = f.read_variable("vg_connection")
+        vg_classifier = vg_beta_star*np.minimum(vg_conn,1) # if closed-closed fieldlines, set var to zero
+        if True:
+
+            vtkreader = pt.vlsvfile.VlsvVtkReader()
+            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 = ["vg_betastar_classifier", "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)
+
+            threshold0 = vtk.vtkThreshold()
+            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])
+            threshold0.Update()
+
+            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)
+                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