diff --git a/CITATION.cff b/CITATION.cff index d66eab2a..3dd47299 100644 --- a/CITATION.cff +++ b/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" diff --git a/Documentation/sphinx/hpcusage.rst b/Documentation/sphinx/hpcusage.rst new file mode 100644 index 00000000..e6447d17 --- /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 fca78f67..3ca5b28c 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/analysator/calculations/__init__.py b/analysator/calculations/__init__.py index 1aea6eaa..011cae3c 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 diff --git a/analysator/calculations/magnetopause_sw_streamline_3d.py b/analysator/calculations/magnetopause_sw_streamline_3d.py index 1cab1e05..8c46162e 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) @@ -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] diff --git a/analysator/calculations/spectra.py b/analysator/calculations/spectra.py index 93f144d4..9145cb9f 100644 --- a/analysator/calculations/spectra.py +++ b/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] diff --git a/analysator/calculations/timeevolution.py b/analysator/calculations/timeevolution.py index 3742c5bd..62c8b989 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)] diff --git a/analysator/vlsvfile/vlsvreader.py b/analysator/vlsvfile/vlsvreader.py index 915bbc18..37d3ecb5 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)") @@ -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: diff --git a/analysator/vlsvfile/vlsvvtkinterface.py b/analysator/vlsvfile/vlsvvtkinterface.py index e0e52953..2e2d010b 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/containerization/lumi_container.sh b/containerization/lumi_container.sh new file mode 100755 index 00000000..729488da --- /dev/null +++ b/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 + + + diff --git a/scripts/magnetopause.py b/scripts/magnetopause.py index 9f8dd66e..1c77f8b2 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) @@ -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 diff --git a/scripts/obliqueshock.py b/scripts/obliqueshock.py index 60cf0b9c..dbe7c152 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 diff --git a/scripts/regions.py b/scripts/regions.py index 7b7895f7..d6da1fd4 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,9 +368,9 @@ 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 - write_flags(writer, magnetopause_SDF, 'SDF_magnetopause') - write_flags(writer, np.where(np.abs(magnetopause_SDF) < 5e6, 1, 0), "flag_magnetopause") + __, 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, '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]) @@ -361,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? @@ -371,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) @@ -381,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]) @@ -392,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') @@ -423,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 @@ -453,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? @@ -477,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 @@ -495,10 +516,16 @@ 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" + 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/magnetopause_sdfs/FID_mpause_{:07d}.vlsv".format(fileid) - RegionFlags(datafile, outfilen, regions=["all"]) + RegionFlags(datafile, outfilen, regions=["magnetopause"]) if __name__ == "__main__": diff --git a/testpackage/run_compare.sh b/testpackage/run_compare.sh index 574a86cb..c0ce387e 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 c8e05a51..f2f6e51c 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 f5d42e73..a05cd73c 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