From 011ca01a9986193ffa9136e4c2df8754e0806276 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sat, 7 Dec 2024 22:02:18 +0200 Subject: [PATCH 01/18] Initial vtkHyperTreeGrid wrapper - works to construct a working HTG and places CellIDs and fileIndices onto the HTG. Need to figure out how to map data from file to VTKs arrays next. --- pyVlsv/vlsvfile.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyVlsv/vlsvfile.py b/pyVlsv/vlsvfile.py index e6c125b6..071795a2 100644 --- a/pyVlsv/vlsvfile.py +++ b/pyVlsv/vlsvfile.py @@ -38,6 +38,7 @@ from vlsvreader import fsDecompositionFromGlobalIds,fsReadGlobalIdsPerRank,fsGlobalIdToGlobalIndex from vlsvwriter import VlsvWriter from vlasiatorreader import VlasiatorReader +from vlsvvtkinterface import vtkVlsvHyperTreeGrid from vlsvparticles import VlsvParticles import reduction From 1c3c814f7ab01dd06f89fb0fda2c4b1ef9f7776c Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sun, 8 Dec 2024 01:22:28 +0200 Subject: [PATCH 02/18] Fileindex-mapped data conversion for SpatialGrid from vlsv to vtkHyperTreeGrid. HTG object (as the extended class) works. examples/write_htg.py added. Starts to look nifty. vtkHyperTreeGridSource availability next? --- examples/write_htg.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 examples/write_htg.py diff --git a/examples/write_htg.py b/examples/write_htg.py new file mode 100644 index 00000000..ad16e5bc --- /dev/null +++ b/examples/write_htg.py @@ -0,0 +1,16 @@ +import vtk +import pytools as pt + +# This initializes a hypertreegrid from the given reader. +reader = pt.vlsvfile.VlsvReader("../../bulk.0002189.vlsv") +htg = pt.vlsvfile.vtkVlsvHyperTreeGrid(reader) + +# These functions grab one SpatialGrid variable and map that to +# the hypertreegrid. Variable vector sizes of 1,2,3,4,9 supported. +htg.addArrayFromVlsv("vg_b_vol") +htg.addArrayFromVlsv("vg_beta_star") + +writer = vtk.vtkXMLHyperTreeGridWriter() +writer.SetFileName("output_EGE.htg") +writer.SetInputData(htg) +writer.Write() From 58c09c83da621bf6acc85b6a6424674f5962d336 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sun, 8 Dec 2024 01:26:13 +0200 Subject: [PATCH 03/18] Adding requirements file --- requirements.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 requirements.txt diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 00000000..b291109e --- /dev/null +++ b/requirements.txt @@ -0,0 +1 @@ +numpy scipy matplotlib scikit-image vtk From a666577a82f3d08bba5fcd0b0e0e835a2ec2d310 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sun, 8 Dec 2024 01:27:53 +0200 Subject: [PATCH 04/18] Linebreaks! --- requirements.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index b291109e..a2682ed3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1,5 @@ -numpy scipy matplotlib scikit-image vtk +numpy +scipy +matplotlib +scikit-image +vtk From 3515be3f6f2a7b32abcecc4ac5af5a60b2a52371 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sun, 8 Dec 2024 01:30:22 +0200 Subject: [PATCH 05/18] Helps if I add the actual HTG file to the index --- pyVlsv/vlsvvtkinterface.py | 298 +++++++++++++++++++++++++++++++++++++ 1 file changed, 298 insertions(+) create mode 100644 pyVlsv/vlsvvtkinterface.py diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py new file mode 100644 index 00000000..4dff49a3 --- /dev/null +++ b/pyVlsv/vlsvvtkinterface.py @@ -0,0 +1,298 @@ +# +# This file is part of Analysator. +# Copyright 2013-2016 Finnish Meteorological Institute +# Copyright 2017-2024 University of Helsinki +# +# For details of usage, see the COPYING file and read the "Rules of the Road" +# at http://www.physics.helsinki.fi/vlasiator/ +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# + +import logging +import struct +import xml.etree.ElementTree as ET +import ast +import numpy as np +import os +import sys +import re +import numbers + +import vtk.numpy_interface +import vtk.numpy_interface.dataset_adapter +from vlsvreader import VlsvReader + +try: + import vtk +except Exception as e: + logging.error("VTK import did not succeed") + raise(e) + + + +'''Extend vtkHyperTreeGrid with a vlsvReader wrapper. + +The extension consists of automatic mapping of a vlsv +SpatialGrid to the HTG structure. Functions are provided +for mapping more SpatialGrid data to the HTG structure as +needed. How well does this function as a VTK source already? +''' +class vtkVlsvHyperTreeGrid(vtk.vtkHyperTreeGrid): + + ''' Set the vlsvReader object for the backend, + otherwise use parent class init, and populate the + hypertreegrid object with CellID data and the file- + index information for easier mapping of vlsv data to + HTG. + ''' + def __init__(self, vlsvFile, *args, **kwargs): + self.vlsvreader = vlsvFile #VlsvReader(vlsvFile) + f = self.vlsvreader + vtk.vtkHyperTreeGrid.__init__(self, args, kwargs) + self.Initialize() + basegridsize = f.get_spatial_mesh_size() + ext = f.get_fsgrid_mesh_extent() + nodeCoordinatesX = f.read(tag="MESH_NODE_CRDS_X", mesh='SpatialGrid') + nodeCoordinatesY = f.read(tag="MESH_NODE_CRDS_Y", mesh='SpatialGrid') + nodeCoordinatesZ = f.read(tag="MESH_NODE_CRDS_Z", mesh='SpatialGrid') + + self.SetDimensions([len(nodeCoordinatesX),len(nodeCoordinatesY),len(nodeCoordinatesZ)]) + self.SetBranchFactor(2) + + xValues = vtk.vtkDoubleArray() + xValues.SetNumberOfValues(int(basegridsize[0]+1)) + # print (basegridsizeNodes, nodeCoordinatesX) + for i,x in enumerate(nodeCoordinatesX): + xValues.SetValue(i, x) + self.SetXCoordinates(xValues) + + yValues = vtk.vtkDoubleArray() + yValues.SetNumberOfValues(int(basegridsize[1]+1)) + for i,x in enumerate(nodeCoordinatesY): + yValues.SetValue(i, x) + self.SetYCoordinates(yValues) + + zValues = vtk.vtkDoubleArray() + zValues.SetNumberOfValues(int(basegridsize[2]+1)) + for i,x in enumerate(nodeCoordinatesZ): + zValues.SetValue(i, x) + self.SetZCoordinates(zValues) + + + # Here we need to actually init the hypertreegrid, + # and at least construct the mapping to leaf nodes + + + cid_offsets = {0: 0} + for i in range(1,f.get_max_refinement_level()+1): + cid_offsets[i] = int(cid_offsets[i-1] + np.prod(np.array([f._VlsvReader__xcells, f._VlsvReader__ycells, f._VlsvReader__zcells]) * 2**(i-1))) # Check next AMR level + + xc = f._VlsvReader__xcells + yc = f._VlsvReader__ycells + zc = f._VlsvReader__zcells + cid_offsets = np.zeros(f.get_max_refinement_level()+1, dtype=np.int64) + isum = 0 + for i in range(0,f.get_max_refinement_level()): + isum = isum + 2**(3*i) * xc * yc * zc + cid_offsets[i+1] = isum + + print(cid_offsets) + + xcells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) + ycells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) + zcells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) + for r in range(f.get_max_refinement_level()+1): + xcells[r] = f._VlsvReader__xcells*2**(r) + ycells[r] = f._VlsvReader__ycells*2**(r) + zcells[r] = f._VlsvReader__zcells*2**(r) + mins = np.array([f._VlsvReader__xmin,f._VlsvReader__ymin,f._VlsvReader__zmin]) + + f._VlsvReader__read_fileindex_for_cellid() + + cidArray = vtk.vtkDoubleArray() + cidArray.SetName('CellID') + cidArray.SetNumberOfValues(0) + self.GetCellData().AddArray(cidArray) + + self.fileIndexArray = vtk.vtkDoubleArray() + self.fileIndexArray.SetName('fileIndex') + self.fileIndexArray.SetNumberOfValues(0) + self.GetCellData().AddArray(self.fileIndexArray) + + f.read_variable_to_cache("CellID") + + cursor = vtk.vtkHyperTreeGridNonOrientedCursor() + offsetIndex = 0 + numvalues = 0 + self.idxToFileIndex = {} + + def refine_cell(cid, level): + # ave = np.array([0.0,0.0,0.0,0.0]) + # print(cid, f.get_cell_indices(cid, reflevels = 0), cid in f._VlsvReader__fileindex_for_cellid.keys()) + cell_lengths = np.array([(f._VlsvReader__xmax - f._VlsvReader__xmin)/(xcells[level]), + (f._VlsvReader__ymax - f._VlsvReader__ymin)/(ycells[level]), + (f._VlsvReader__zmax - f._VlsvReader__zmin)/(zcells[level])]) + # print(f.get_cell_coordinates(cid-1)) + cellind = f.get_cell_indices(cid, reflevels = level) + + # print(mins, cellind, cell_lengths) + cellcoordinates = mins + (cellind + 0.25)*cell_lengths + # print(cellind,cellcoordinates) + + cell_lengths = np.array([(f._VlsvReader__xmax - f._VlsvReader__xmin)/(xcells[level+1]), + (f._VlsvReader__ymax - f._VlsvReader__ymin)/(ycells[level+1]), + (f._VlsvReader__zmax - f._VlsvReader__zmin)/(zcells[level+1])]) + # print(cell_lengths) + cellind = np.array((cellcoordinates - mins)/cell_lengths, dtype = np.int32) + # print(cellind) + # print(xcells[level],xcells[level+1]) + # print(xcells) + # sys.exit() + # delta = [[0,0,0],[0,0,1],[0,1,0],[0,1,1], + # [1,0,0],[1,0,1],[1,1,0],[1,1,1]] + delta = [[0,0,0],[1,0,0],[0,1,0],[1,1,0], + [0,0,1],[1,0,1],[0,1,1],[1,1,1]] + + for c in range(8): + cursor.ToChild(c) # cell 0/0 + # print(cellind, delta[c]) + ci = cellind + delta[c] + # print(ci) + cellid = int(cid_offsets[level+1] + ci[0] + 2**(level+1)*ci[1]*xc + 4**(level+1)*ci[2]*xc*yc + 1) + # print(f.get_cell_indices(cid, reflevels = level)) + if (cellid in f._VlsvReader__fileindex_for_cellid.keys()): + # Assigns an index for the value of the cell pointed to by the current cursor state + idx = cursor.GetGlobalNodeIndex() + self.GetCellData().SetActiveScalars('CellID') + cidArray.InsertTuple1(idx, cellid) + self.GetCellData().SetActiveScalars('fileIndex') + self.fileIndexArray.InsertTuple1(idx, f._VlsvReader__fileindex_for_cellid[cellid]) + # numvalues += 1 + self.idxToFileIndex[idx] = f._VlsvReader__fileindex_for_cellid[cellid] + # ave += np.array([var_x(cellid, 'proton/vg_v'),*var(cellid, 'vg_b_vol')]) + else: + # print(cellid) + idx = cursor.GetGlobalNodeIndex() + cursor.SubdivideLeaf() # cell 0 + # print(cellid, "not found, refining again") + refine_cell(cellid, level+1) + # ave += refval + # htg.GetCellData().SetActiveScalars('vg_v') + #scalarArray.InsertTuple1(idx, refval[0]) + # htg.GetCellData().SetActiveVectors('vg_b_vol') + #vectorArray.InsertTuple3(idx, *refval[1:]) + + cursor.ToParent() # cell 0 + idx = cursor.GetGlobalNodeIndex() + + + return #ave/8 + + + def handle_cell(cid): + # Leaf cell, just add + if (cid in f._VlsvReader__fileindex_for_cellid.keys()): + # Assigns an index for the value of the cell pointed to by the current cursor state + idx = cursor.GetGlobalNodeIndex() + self.GetCellData().SetActiveScalars('CellID') + cidArray.InsertTuple1(idx, cid) + self.GetCellData().SetActiveScalars('fileIndex') + self.fileIndexArray.InsertTuple1(idx, f._VlsvReader__fileindex_for_cellid[cid]) + # numvalues += 1 + self.idxToFileIndex[idx] = f._VlsvReader__fileindex_for_cellid[cid] + # Cell not in list -> must be refined (assume cid is in basegrid) + else: + idx = cursor.GetGlobalNodeIndex() + + # print(cid , "not found, refining") + cursor.SubdivideLeaf() # cell 0 + refine_cell(cid, 0) + # self.GetCellData().SetActiveScalars('vg_v') + #scalarArray.InsertTuple1(idx, ave[0]) + # self.GetCellData().SetActiveVectors('vg_b_vol') + #vectorArray.InsertTuple3(idx, *ave[1:]) + + # pass + # idx = cursor.GetGlobalNodeIndex() + # scalarArray.InsertTuple1(idx, vardummy(cid)) + + nc = 0 + for basegridIdx in range(int(np.prod(basegridsize))): + cid = basegridIdx +1 + # Set cursor on cell #0 and set the start index + self.InitializeNonOrientedCursor(cursor, basegridIdx, True) + cursor.SetGlobalIndexStart(offsetIndex) # cell 0 + # print("handling ",cid) + handle_cell(cid) + + offsetIndex += cursor.GetTree().GetNumberOfVertices() + nc += 1 + if cid > np.prod(basegridsize): + print("tried to go over") + break + + '''This function adds one SpatialGrid variable from the reader object and maps + that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. + ''' + def addArrayFromVlsv(self, varname): + array = vtk.vtkDoubleArray() + # cidArray2.DeepCopy(fileIndexArray) + data = self.vlsvreader.read_variable(varname) + + test_var = data[0] + if np.isscalar(test_var): + varlen = 1 + else: + varlen = len(test_var) + # varlen = data[:,np.newaxis].shape[1] + print(varlen) + array.SetNumberOfComponents(varlen) + array.SetNumberOfTuples(self.fileIndexArray.GetNumberOfTuples()) + array.SetName(varname) + # print(cidArray2) + # newdata = np.array([ciddata[self.idxToFileIndex[idx]] for idx in self.idxToFileIndex.keys()]) + if varlen == 1: + for idx,fileIndex in self.idxToFileIndex.items(): + array.SetTuple1(idx, data[fileIndex]) + elif varlen == 2: + for idx,fileIndex in self.idxToFileIndex.items(): + array.SetTuple2(idx, *data[fileIndex]) + elif varlen == 3: + for idx,fileIndex in self.idxToFileIndex.items(): + array.SetTuple3(idx, *data[fileIndex]) + elif varlen == 4: + for idx,fileIndex in self.idxToFileIndex.items(): + array.SetTuple4(idx, *data[fileIndex]) + elif varlen == 9: + for idx,fileIndex in self.idxToFileIndex.items(): + array.SetTuple9(idx, *data[fileIndex]) + else: + raise RuntimeError("No vtk SetTuple wrapper function for varlen = " + str(varlen)) + + + self.GetCellData().AddArray(array) + + # vtk.numpy_interface.numpy_to_vtk(newdata) + # cidArray2.SetData(newdata, len(self.idxToFileIndex), 1) + + + # ... we do not have intermediate node data stored. + # it could be calculated while constructing, but that + # defeats the purpose (limit file read/memory use to + # a threshold). Also, if we can just map directly to + # the file/a buffer, we do not need to do fancy recursive + # traversal downsampling. + From 180fc675477d7fcb506152d108de3b8a0348fd4b Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Tue, 10 Dec 2024 15:32:34 +0200 Subject: [PATCH 06/18] A functional vtk reader object, but with very few and hardcoded functions (just gets one scalar) --- pyVlsv/vlsvvtkinterface.py | 113 ++++++++++++++++++++++++++++++++++++- 1 file changed, 111 insertions(+), 2 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 4dff49a3..4c68fd85 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -33,8 +33,10 @@ import vtk.numpy_interface import vtk.numpy_interface.dataset_adapter -from vlsvreader import VlsvReader - +if __name__ != "__main__": + from vlsvreader import VlsvReader +else: + import pytools as pt try: import vtk except Exception as e: @@ -296,3 +298,110 @@ def addArrayFromVlsv(self, varname): # the file/a buffer, we do not need to do fancy recursive # traversal downsampling. +from vtk.util.vtkAlgorithm import VTKPythonAlgorithmBase + +class VlsvVtkReader(VTKPythonAlgorithmBase): + def __init__(self): + VTKPythonAlgorithmBase.__init__(self, nInputPorts = 0, nOutputPorts=1, outputType='vtkHyperTreeGrid') + self.__FileName = None + self.__htg = None + self.__reader = None + + def SetFileName(self, filename): + if filename != self.__FileName: + self.Modified() + self.__FileName = filename + self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) + + + def GetFileName(self): + return self.__FileName + + ''' + def FillOutputPortInformation(self, port, info): + pass + + def RequestDataObject(self, request, inInfo, outInfo): + #This is where you can create output data objects if the output DATA_TYPE_NAME() is not a concrete type. + pass + ''' + + def RequestInformation(self, request, inInfo, outInfo): + + info = outInfo.GetInformationObject(0) + print("VlsvVtkReader RequestInformation:") + print(info) + if self.__htg is None: + self.__htg = vtkVlsvHyperTreeGrid(self.__reader) + self.__htg.addArrayFromVlsv("vg_beta_star") + self.__htg.GetCellData().SetActiveScalars("vg_beta_star") + + dims = self.__htg.GetExtent() + info.Set(vtk.vtkStreamingDemandDrivenPipeline.WHOLE_EXTENT(), *dims) + + return 1 + + ''' + def RequestUpdateExtent(self, request, inInfo, outInfo): + pass + ''' + + def RequestData(self, request, inInfo, outInfo): + + print("VlsvVtkReader RequestData:") + print(outInfo) + + if self.__htg is None: + self.__htg = vtkVlsvHyperTreeGrid(self.__reader) + + output = vtk.vtkHyperTreeGrid.GetData(outInfo) + + output.ShallowCopy(self.__htg) + + return 1 + + + +def __main__(): + import pytools as pt + # This initializes a hypertreegrid from the given reader. + reader = VlsvVtkReader() + reader.SetFileName("../../bulk.0002189.vlsv") + + cf = vtk.vtkHyperTreeGridContour() + cf.SetInputConnection(reader.GetOutputPort()) + cf.SetValue(0, 1) + + m = vtk.vtkPolyDataMapper() + m.SetInputConnection(cf.GetOutputPort()) + + a = vtk.vtkActor() + a.SetMapper(m) + + ren = vtk.vtkRenderer() + ren.AddActor(a) + + renWin = vtk.vtkRenderWindow() + renWin.AddRenderer(ren) + renWin.SetSize(600, 600) + + renWin.Render() + import time + time.sleep(10) + htg = reader.GetOutputDataObject(0) + # htg.__vtkVlsvHyperTreeGrid_VlsvAddArray('vg_b_vol') + # htg = reader.GetOutputData() + # print(htg) + # These functions grab one SpatialGrid variable and map that to + # the hypertreegrid. Variable vector sizes of 1,2,3,4,9 supported. + # htg.addArrayFromVlsv("vg_b_vol") + # htg.addArrayFromVlsv("vg_beta_star") + + writer = vtk.vtkXMLHyperTreeGridWriter() + writer.SetFileName("output_EGE.htg") + writer.SetInputData(htg) + writer.Write() + + +if __name__ == "__main__": + __main__() \ No newline at end of file From 1fb49772629c7a4d62cb2987ac74085067dc4079 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Wed, 11 Dec 2024 12:41:18 +0200 Subject: [PATCH 07/18] Paraview reader plugin, and optimizations on the constructor (which should be done in C++ or as a stored descriptor, apparently those are fast to reuse) --- pyVlsv/vlsvVtkParaReader.py | 164 ++++++++++++++++++++++++ pyVlsv/vlsvfile.py | 2 +- pyVlsv/vlsvreader.py | 35 ++++++ pyVlsv/vlsvvtkinterface.py | 241 +++++++++++++++++++++++++----------- 4 files changed, 366 insertions(+), 76 deletions(-) create mode 100644 pyVlsv/vlsvVtkParaReader.py diff --git a/pyVlsv/vlsvVtkParaReader.py b/pyVlsv/vlsvVtkParaReader.py new file mode 100644 index 00000000..49225a2b --- /dev/null +++ b/pyVlsv/vlsvVtkParaReader.py @@ -0,0 +1,164 @@ +"""This module demonstrates various ways of adding +VTKPythonAlgorithmBase subclasses as filters, sources, readers, +and writers in ParaView""" + +# https://github.com/Kitware/ParaView/blob/master/Examples/Plugins/PythonAlgorithm/PythonAlgorithmExamples.py + +# This is module to import. It provides VTKPythonAlgorithmBase, the base class +# for all python-based vtkAlgorithm subclasses in VTK and decorators used to +# 'register' the algorithm with ParaView along with information about UI. +from paraview.util.vtkAlgorithm import * +import pytools as pt + +# Use the analysator reader, but include Paraview decorators via inheriting/chaining +# necessary methods etc + +#------------------------------------------------------------------------------ +# A reader example. +#------------------------------------------------------------------------------ + +# To add a reader, we can use the following decorators +# @smproxy.source(name="PythonCSVReader", label="Python-based CSV Reader") +# @smhint.xml("""""") +# or directly use the "@reader" decorator. +@smproxy.reader(name="PythonVLSVReader", label="Python-based VLSV Reader, outputs htg", + extensions="vlsv", file_description="VLSV files") +class PythonPvVLSVReader(pt.vlsvfile.VlsvVtkReader): + """A reader that wraps an Analysator VLSV file reader. + """ + def __init__(self): + super().__init__() + + + @smproperty.stringvector(name="FileName") + @smdomain.filelist() + @smhint.filechooser(extensions="vlsv", file_description="Vlasiator VLSV files") + def SetFileName(self, filename): + """Specify filename for the file to read.""" + print(filename) + super().SetFileName(filename) + + @smproperty.doublevector(name="TimestepValues", information_only="1", si_class="vtkSITimeStepsProperty") + def GetTimestepValues(self): + return None + + # Array selection API is typical with readers in VTK + # This is intended to allow ability for users to choose which arrays to + # load. To expose that in ParaView, simply use the + # smproperty.dataarrayselection(). + # This method **must** return a `vtkDataArraySelection` instance. + @smproperty.dataarrayselection(name="Arrays") + def GetDataArraySelection(self): + return super().GetDataArraySelection() + + +#------------------------------------------------------------------------------ +# A writer example. +#------------------------------------------------------------------------------ +''' +@smproxy.writer(extensions="npz", file_description="NumPy Compressed Arrays", support_reload=False) +@smproperty.input(name="Input", port_index=0) +@smdomain.datatype(dataTypes=["vtkTable"], composite_data_supported=False) +class NumpyWriter(VTKPythonAlgorithmBase): + def __init__(self): + VTKPythonAlgorithmBase.__init__(self, nInputPorts=1, nOutputPorts=0, inputType='vtkTable') + self._filename = None + + @smproperty.stringvector(name="FileName", panel_visibility="never") + @smdomain.filelist() + def SetFileName(self, fname): + """Specify filename for the file to write.""" + if self._filename != fname: + self._filename = fname + self.Modified() + + def RequestData(self, request, inInfoVec, outInfoVec): + from vtkmodules.vtkCommonDataModel import vtkTable + from vtkmodules.numpy_interface import dataset_adapter as dsa + + table = dsa.WrapDataObject(vtkTable.GetData(inInfoVec[0], 0)) + kwargs = {} + for aname in table.RowData.keys(): + kwargs[aname] = table.RowData[aname] + + import numpy + numpy.savez_compressed(self._filename, **kwargs) + return 1 + + def Write(self): + self.Modified() + self.Update() +''' +#------------------------------------------------------------------------------ +# A filter example. +#------------------------------------------------------------------------------ +''' +@smproxy.filter() +@smproperty.input(name="InputTable", port_index=1) +@smdomain.datatype(dataTypes=["vtkTable"], composite_data_supported=False) +@smproperty.input(name="InputDataset", port_index=0) +@smdomain.datatype(dataTypes=["vtkDataSet"], composite_data_supported=False) +class ExampleTwoInputFilter(VTKPythonAlgorithmBase): + def __init__(self): + VTKPythonAlgorithmBase.__init__(self, nInputPorts=2, nOutputPorts=1, outputType="vtkPolyData") + + def FillInputPortInformation(self, port, info): + if port == 0: + info.Set(self.INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet") + else: + info.Set(self.INPUT_REQUIRED_DATA_TYPE(), "vtkTable") + return 1 + + def RequestData(self, request, inInfoVec, outInfoVec): + from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData + input0 = vtkDataSet.GetData(inInfoVec[0], 0) + input1 = vtkDataSet.GetData(inInfoVec[1], 0) + output = vtkPolyData.GetData(outInfoVec, 0) + # do work + print("Pretend work done!") + return 1 + +@smproxy.filter() +@smproperty.input(name="Input") +@smdomain.datatype(dataTypes=["vtkDataSet"], composite_data_supported=False) +class PreserveInputTypeFilter(VTKPythonAlgorithmBase): + """ + Example filter demonstrating how to write a filter that preserves the input + dataset type. + """ + def __init__(self): + super().__init__(nInputPorts=1, nOutputPorts=1, outputType="vtkDataSet") + + def RequestDataObject(self, request, inInfo, outInfo): + inData = self.GetInputData(inInfo, 0, 0) + outData = self.GetOutputData(outInfo, 0) + assert inData is not None + if outData is None or (not outData.IsA(inData.GetClassName())): + outData = inData.NewInstance() + outInfo.GetInformationObject(0).Set(outData.DATA_OBJECT(), outData) + return super().RequestDataObject(request, inInfo, outInfo) + + def RequestData(self, request, inInfo, outInfo): + inData = self.GetInputData(inInfo, 0, 0) + outData = self.GetOutputData(outInfo, 0) + print("input type =", inData.GetClassName()) + print("output type =", outData.GetClassName()) + assert outData.IsA(inData.GetClassName()) + return 1 +''' + +def test_PythonVLSVReader(fname): + reader = PythonVLSVReader() + reader.SetFileName(fname) + reader.Update() + assert reader.GetOutputDataObject(0).GetNumberOfRows() > 0 + +if __name__ == "__main__": + #test_PythonSuperquadricSource() + test_PythonVLSVReader("/tmp/data.csv") + + from paraview.detail.pythonalgorithm import get_plugin_xmls + from xml.dom.minidom import parseString + for xml in get_plugin_xmls(globals()): + dom = parseString(xml) + print(dom.toprettyxml(" ","\n")) diff --git a/pyVlsv/vlsvfile.py b/pyVlsv/vlsvfile.py index 071795a2..2d645172 100644 --- a/pyVlsv/vlsvfile.py +++ b/pyVlsv/vlsvfile.py @@ -38,7 +38,7 @@ from vlsvreader import fsDecompositionFromGlobalIds,fsReadGlobalIdsPerRank,fsGlobalIdToGlobalIndex from vlsvwriter import VlsvWriter from vlasiatorreader import VlasiatorReader -from vlsvvtkinterface import vtkVlsvHyperTreeGrid +from vlsvvtkinterface import vtkVlsvHyperTreeGrid,VlsvVtkReader from vlsvparticles import VlsvParticles import reduction diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py index df0157d0..b856f230 100644 --- a/pyVlsv/vlsvreader.py +++ b/pyVlsv/vlsvreader.py @@ -611,6 +611,41 @@ def __check_datareducer(self, name, reducer): return reducer_ok + def get_variables(self): + + varlist = [] + + for child in self.__xml_root: + if child.tag == "VARIABLE" and "name" in child.attrib: + name = child.attrib["name"] + varlist.append(name) + + return varlist + + def get_reducers(self): + + varlist = [] + + reducer_max_len = 0 + + for reducer_reg in [datareducers, multipopdatareducers, v5reducers, multipopv5reducers]: + for k in reducer_reg.keys(): + reducer_max_len = max(reducer_max_len, len(k)) + + + + for reducer_reg in [datareducers, multipopdatareducers, v5reducers, multipopv5reducers]: + for name, reducer in reducer_reg.items(): + self.__current_reducer_tree_nodes.clear() + if self.__check_datareducer(name,reducer): + if name[:3] == 'pop': + for pop in self.active_populations: + varlist.append(pop+'/'+name[4:]) + else: + varlist.append(name) + + return varlist + def list(self, parameter=True, variable=True, mesh=False, datareducer=False, operator=False, other=False): ''' Print out a description of the content of the file. Useful diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 4c68fd85..4761081c 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -22,19 +22,14 @@ # import logging -import struct -import xml.etree.ElementTree as ET -import ast import numpy as np -import os -import sys -import re -import numbers +import cProfile import vtk.numpy_interface import vtk.numpy_interface.dataset_adapter if __name__ != "__main__": from vlsvreader import VlsvReader + import pytools as pt else: import pytools as pt try: @@ -44,7 +39,6 @@ raise(e) - '''Extend vtkHyperTreeGrid with a vlsvReader wrapper. The extension consists of automatic mapping of a vlsv @@ -111,15 +105,15 @@ def __init__(self, vlsvFile, *args, **kwargs): isum = isum + 2**(3*i) * xc * yc * zc cid_offsets[i+1] = isum - print(cid_offsets) - - xcells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) - ycells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) - zcells = np.zeros((f.get_max_refinement_level()+1), dtype=np.int64) - for r in range(f.get_max_refinement_level()+1): - xcells[r] = f._VlsvReader__xcells*2**(r) - ycells[r] = f._VlsvReader__ycells*2**(r) - zcells[r] = f._VlsvReader__zcells*2**(r) + # print(cid_offsets) + max_ref_level = f.get_max_refinement_level() + xcells = np.zeros((max_ref_level+1), dtype=np.int64) + ycells = np.zeros((max_ref_level+1), dtype=np.int64) + zcells = np.zeros((max_ref_level+1), dtype=np.int64) + for r in range(max_ref_level+1): + xcells[r] = xc*2**(r) + ycells[r] = yc*2**(r) + zcells[r] = zc*2**(r) mins = np.array([f._VlsvReader__xmin,f._VlsvReader__ymin,f._VlsvReader__zmin]) f._VlsvReader__read_fileindex_for_cellid() @@ -137,53 +131,78 @@ def __init__(self, vlsvFile, *args, **kwargs): f.read_variable_to_cache("CellID") cursor = vtk.vtkHyperTreeGridNonOrientedCursor() - offsetIndex = 0 - numvalues = 0 self.idxToFileIndex = {} + xmin,xmax,ymin,ymax,zmin,zmax = f._VlsvReader__xmin,f._VlsvReader__xmax,f._VlsvReader__ymin,f._VlsvReader__ymax,f._VlsvReader__zmin,f._VlsvReader__zmax, + + fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid + delta = [[0,0,0],[1,0,0],[0,1,0],[1,1,0], + [0,0,1],[1,0,1],[0,1,1],[1,1,1]] + + def get_cell_indices(cellids, reflevels): + ''' Single-cellid specialization with external constants. + Returns a given cell's indices as a numpy array + + :param cellid: The cell's ID, numpy array + :param reflevel: The cell's refinement level in the AMR + :returns: a numpy array with the coordinates + + .. seealso:: :func:`get_cellid` + + .. note:: The cell ids go from 1 .. max not from 0 + ''' + + # # Calculating the index of the first cell at this reflevel + # index_at_reflevel = np.zeros(max_refinement_level+1, dtype=np.int64) + # isum = 0 + # for i in range(0,max_refinement_level): + # isum = isum + 2**(3*i) * xcells * ycells * zcells + # index_at_reflevel[i+1] = isum + + # Get cell indices + # print(cellids, reflevels, cid_offsets) + + cellids = cellids - 1 - cid_offsets[reflevels] + cellindices = np.full(3, -1) + cellindices[0] = (cellids)%(xcells[reflevels]) + cellindices[1] = ((cellids)//(xcells[reflevels]))%(ycells[reflevels]) + cellindices[2] = (cellids)//(xcells[reflevels]*ycells[reflevels]) + + # print(cellindices, f.get_cell_indices(cellids,reflevels)) + + return cellindices + + cell_lengths_levels = [] + for level in range(max_ref_level+1): + cell_lengths_levels.append(np.array([(xmax - xmin)/(xcells[level]), + (ymax - ymin)/(ycells[level]), + (zmax - zmin)/(zcells[level])])) + def refine_cell(cid, level): - # ave = np.array([0.0,0.0,0.0,0.0]) - # print(cid, f.get_cell_indices(cid, reflevels = 0), cid in f._VlsvReader__fileindex_for_cellid.keys()) - cell_lengths = np.array([(f._VlsvReader__xmax - f._VlsvReader__xmin)/(xcells[level]), - (f._VlsvReader__ymax - f._VlsvReader__ymin)/(ycells[level]), - (f._VlsvReader__zmax - f._VlsvReader__zmin)/(zcells[level])]) - # print(f.get_cell_coordinates(cid-1)) - cellind = f.get_cell_indices(cid, reflevels = level) - - # print(mins, cellind, cell_lengths) - cellcoordinates = mins + (cellind + 0.25)*cell_lengths - # print(cellind,cellcoordinates) - - cell_lengths = np.array([(f._VlsvReader__xmax - f._VlsvReader__xmin)/(xcells[level+1]), - (f._VlsvReader__ymax - f._VlsvReader__ymin)/(ycells[level+1]), - (f._VlsvReader__zmax - f._VlsvReader__zmin)/(zcells[level+1])]) - # print(cell_lengths) - cellind = np.array((cellcoordinates - mins)/cell_lengths, dtype = np.int32) - # print(cellind) - # print(xcells[level],xcells[level+1]) - # print(xcells) - # sys.exit() - # delta = [[0,0,0],[0,0,1],[0,1,0],[0,1,1], - # [1,0,0],[1,0,1],[1,1,0],[1,1,1]] - delta = [[0,0,0],[1,0,0],[0,1,0],[1,1,0], - [0,0,1],[1,0,1],[0,1,1],[1,1,1]] + + cellind = get_cell_indices(cid, level) + cellcoordinates = mins + (cellind + 0.25)*cell_lengths_levels[level] + + cellind = np.array(np.divide((cellcoordinates - mins),cell_lengths_levels[level+1]), dtype = np.int32) + for c in range(8): cursor.ToChild(c) # cell 0/0 # print(cellind, delta[c]) ci = cellind + delta[c] # print(ci) - cellid = int(cid_offsets[level+1] + ci[0] + 2**(level+1)*ci[1]*xc + 4**(level+1)*ci[2]*xc*yc + 1) + cellid = int(cid_offsets[level+1] + ci[0] + xcells[level+1]*ci[1] + ci[2]*xcells[level+1]*ycells[level+1] + 1) # print(f.get_cell_indices(cid, reflevels = level)) - if (cellid in f._VlsvReader__fileindex_for_cellid.keys()): + if (cellid in fileindex_for_cellid.keys()): # Assigns an index for the value of the cell pointed to by the current cursor state idx = cursor.GetGlobalNodeIndex() + # print(idx, cellid) self.GetCellData().SetActiveScalars('CellID') cidArray.InsertTuple1(idx, cellid) self.GetCellData().SetActiveScalars('fileIndex') - self.fileIndexArray.InsertTuple1(idx, f._VlsvReader__fileindex_for_cellid[cellid]) + self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) # numvalues += 1 - self.idxToFileIndex[idx] = f._VlsvReader__fileindex_for_cellid[cellid] + self.idxToFileIndex[idx] = fileindex_for_cellid[cellid] # ave += np.array([var_x(cellid, 'proton/vg_v'),*var(cellid, 'vg_b_vol')]) else: # print(cellid) @@ -203,25 +222,25 @@ def refine_cell(cid, level): return #ave/8 - - def handle_cell(cid): + def handle_cell(cellid): # Leaf cell, just add - if (cid in f._VlsvReader__fileindex_for_cellid.keys()): + if (cellid in fileindex_for_cellid.keys()): # Assigns an index for the value of the cell pointed to by the current cursor state idx = cursor.GetGlobalNodeIndex() self.GetCellData().SetActiveScalars('CellID') - cidArray.InsertTuple1(idx, cid) + # print(idx, cid) + cidArray.InsertTuple1(idx, cellid) self.GetCellData().SetActiveScalars('fileIndex') - self.fileIndexArray.InsertTuple1(idx, f._VlsvReader__fileindex_for_cellid[cid]) + self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) # numvalues += 1 - self.idxToFileIndex[idx] = f._VlsvReader__fileindex_for_cellid[cid] + self.idxToFileIndex[idx] = fileindex_for_cellid[cellid] # Cell not in list -> must be refined (assume cid is in basegrid) else: idx = cursor.GetGlobalNodeIndex() # print(cid , "not found, refining") cursor.SubdivideLeaf() # cell 0 - refine_cell(cid, 0) + refine_cell(cellid, 0) # self.GetCellData().SetActiveScalars('vg_v') #scalarArray.InsertTuple1(idx, ave[0]) # self.GetCellData().SetActiveVectors('vg_b_vol') @@ -231,20 +250,45 @@ def handle_cell(cid): # idx = cursor.GetGlobalNodeIndex() # scalarArray.InsertTuple1(idx, vardummy(cid)) - nc = 0 - for basegridIdx in range(int(np.prod(basegridsize))): - cid = basegridIdx +1 - # Set cursor on cell #0 and set the start index - self.InitializeNonOrientedCursor(cursor, basegridIdx, True) - cursor.SetGlobalIndexStart(offsetIndex) # cell 0 - # print("handling ",cid) - handle_cell(cid) - - offsetIndex += cursor.GetTree().GetNumberOfVertices() - nc += 1 - if cid > np.prod(basegridsize): - print("tried to go over") - break + + def loop(): + offsetIndex = 0 + + for basegridIdx in range(int(np.prod(basegridsize))): + cid = basegridIdx +1 + # Set cursor on cell #0 and set the start index + self.InitializeNonOrientedCursor(cursor, basegridIdx, True) + cursor.SetGlobalIndexStart(offsetIndex) # cell 0 + # print("handling ",cid) + handle_cell(cid) + + offsetIndex += cursor.GetTree().GetNumberOfVertices() + if cid > np.prod(basegridsize): + print("tried to go over") + break + + import pstats + from pstats import SortKey + with cProfile.Profile() as pr: + loop() + ps = pstats.Stats(pr).sort_stats().reverse_order() + ps.print_stats() + + # loop() + + + + def findVariablesFromVlsv(self): + + vars = self.vlsvreader.get_variables() + + # for var in vars: + # array = vtk.vtkDoubleArray() + # array.SetName(var) + # self.GetCellData().AddArray(array) + + return vars + '''This function adds one SpatialGrid variable from the reader object and maps that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. @@ -260,7 +304,7 @@ def addArrayFromVlsv(self, varname): else: varlen = len(test_var) # varlen = data[:,np.newaxis].shape[1] - print(varlen) + # print(varlen) array.SetNumberOfComponents(varlen) array.SetNumberOfTuples(self.fileIndexArray.GetNumberOfTuples()) array.SetName(varname) @@ -300,12 +344,29 @@ def addArrayFromVlsv(self, varname): from vtk.util.vtkAlgorithm import VTKPythonAlgorithmBase + +def createModifiedCallback(anobject): + import weakref + weakref_obj = weakref.ref(anobject) + anobject = None + def _markmodified(*args, **kwars): + o = weakref_obj() + if o is not None: + o.Modified() + return _markmodified + +# Should implement these, more or less: +# https://vtk.org/doc/nightly/html/classvtkXMLHyperTreeGridReader.html + class VlsvVtkReader(VTKPythonAlgorithmBase): def __init__(self): VTKPythonAlgorithmBase.__init__(self, nInputPorts = 0, nOutputPorts=1, outputType='vtkHyperTreeGrid') self.__FileName = None self.__htg = None self.__reader = None + self.__cellarrays = [] + self._arrayselection = vtk.vtkDataArraySelection() + self._arrayselection.AddObserver("ModifiedEvent", createModifiedCallback(self)) def SetFileName(self, filename): if filename != self.__FileName: @@ -317,6 +378,9 @@ def SetFileName(self, filename): def GetFileName(self): return self.__FileName + + + ''' def FillOutputPortInformation(self, port, info): pass @@ -333,8 +397,14 @@ def RequestInformation(self, request, inInfo, outInfo): print(info) if self.__htg is None: self.__htg = vtkVlsvHyperTreeGrid(self.__reader) - self.__htg.addArrayFromVlsv("vg_beta_star") - self.__htg.GetCellData().SetActiveScalars("vg_beta_star") + vars = self.__htg.findVariablesFromVlsv() + for name in vars: + if "vg" in name: + # print(name) + self._arrayselection.AddArray(name) + self.__cellarrays.append(name) + # self.__htg.addArrayFromVlsv("vg_beta_star") + # self.__htg.GetCellData().SetActiveScalars("vg_beta_star") dims = self.__htg.GetExtent() info.Set(vtk.vtkStreamingDemandDrivenPipeline.WHOLE_EXTENT(), *dims) @@ -353,26 +423,44 @@ def RequestData(self, request, inInfo, outInfo): if self.__htg is None: self.__htg = vtkVlsvHyperTreeGrid(self.__reader) + + for name in self.__cellarrays: + if self._arrayselection.ArrayIsEnabled(name): + self.__htg.addArrayFromVlsv(name) + output = vtk.vtkHyperTreeGrid.GetData(outInfo) + + output.ShallowCopy(self.__htg) return 1 + + def ReadAllScalarsOn(self): + self._arrayselection.EnableAllArrays() + # for name in self.__cellarrays: + + # self.__htg.addArrayFromVlsv(name) + def GetDataArraySelection(self): + return self._arrayselection + def __main__(): import pytools as pt # This initializes a hypertreegrid from the given reader. reader = VlsvVtkReader() - reader.SetFileName("../../bulk.0002189.vlsv") - + reader.SetFileName("/home/mjalho/Downloads/bulk.0002189.vlsv") + reader.ReadAllScalarsOn() + reader.Update() cf = vtk.vtkHyperTreeGridContour() cf.SetInputConnection(reader.GetOutputPort()) cf.SetValue(0, 1) m = vtk.vtkPolyDataMapper() + m.SetInputConnection(cf.GetOutputPort()) a = vtk.vtkActor() @@ -385,10 +473,12 @@ def __main__(): renWin.AddRenderer(ren) renWin.SetSize(600, 600) - renWin.Render() + # renWin.Render() import time time.sleep(10) + htg = reader.GetOutputDataObject(0) + print(htg) # htg.__vtkVlsvHyperTreeGrid_VlsvAddArray('vg_b_vol') # htg = reader.GetOutputData() # print(htg) @@ -399,6 +489,7 @@ def __main__(): writer = vtk.vtkXMLHyperTreeGridWriter() writer.SetFileName("output_EGE.htg") + writer.SetInputData(htg) writer.Write() From 89976990a3ac916e042556d18b2d41ea5672882e Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Wed, 11 Dec 2024 13:59:28 +0200 Subject: [PATCH 08/18] Datareducer visibility on --- pyVlsv/vlsvvtkinterface.py | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 4761081c..08efe26b 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -278,9 +278,14 @@ def loop(): - def findVariablesFromVlsv(self): + def findVariablesFromVlsv(self, getReducers = False): vars = self.vlsvreader.get_variables() + if getReducers: + reducers = self.vlsvreader.get_reducers() + vars_set = set(vars) + vars_set.update(reducers) + vars = list(vars_set) # for var in vars: # array = vtk.vtkDoubleArray() @@ -299,10 +304,22 @@ def addArrayFromVlsv(self, varname): data = self.vlsvreader.read_variable(varname) test_var = data[0] + if test_var is None: + logging.warning("varname " + varname + "returned None; skipping") + return False if np.isscalar(test_var): varlen = 1 + elif test_var.ndim == 1: + varlen = len(test_var) + elif test_var.ndim > 1: + varlen = np.prod(test_var.shape) else: - varlen = len(test_var) + logging.warning("Weird output from " + varname) + return False + # varlen = 0 + # print("test var is ", test_var, 'for', varname) + + # varlen = data[:,np.newaxis].shape[1] # print(varlen) array.SetNumberOfComponents(varlen) @@ -324,12 +341,13 @@ def addArrayFromVlsv(self, varname): array.SetTuple4(idx, *data[fileIndex]) elif varlen == 9: for idx,fileIndex in self.idxToFileIndex.items(): - array.SetTuple9(idx, *data[fileIndex]) + array.SetTuple9(idx, *np.reshape(data[fileIndex],(9))) else: raise RuntimeError("No vtk SetTuple wrapper function for varlen = " + str(varlen)) self.GetCellData().AddArray(array) + return True # vtk.numpy_interface.numpy_to_vtk(newdata) # cidArray2.SetData(newdata, len(self.idxToFileIndex), 1) @@ -397,7 +415,7 @@ def RequestInformation(self, request, inInfo, outInfo): print(info) if self.__htg is None: self.__htg = vtkVlsvHyperTreeGrid(self.__reader) - vars = self.__htg.findVariablesFromVlsv() + vars = self.__htg.findVariablesFromVlsv(getReducers=True) for name in vars: if "vg" in name: # print(name) @@ -426,7 +444,9 @@ def RequestData(self, request, inInfo, outInfo): for name in self.__cellarrays: if self._arrayselection.ArrayIsEnabled(name): - self.__htg.addArrayFromVlsv(name) + success = self.__htg.addArrayFromVlsv(name) + if not success: + self._arrayselection.RemoveArrayByName(name) output = vtk.vtkHyperTreeGrid.GetData(outInfo) From f8512ace593acfec99abc5b008193e54cfd46d43 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Fri, 13 Dec 2024 10:50:45 +0200 Subject: [PATCH 09/18] A bunch of late-night optimizations to repo --- pyVlsv/vlsvvtkinterface.py | 221 +++++++++++++++++++++++++++++-------- requirements.txt | 3 +- 2 files changed, 179 insertions(+), 45 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 08efe26b..20bd0c50 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -27,18 +27,14 @@ import vtk.numpy_interface import vtk.numpy_interface.dataset_adapter -if __name__ != "__main__": - from vlsvreader import VlsvReader - import pytools as pt -else: - import pytools as pt +import pytools as pt try: import vtk except Exception as e: logging.error("VTK import did not succeed") raise(e) - +from numba import jit '''Extend vtkHyperTreeGrid with a vlsvReader wrapper. The extension consists of automatic mapping of a vlsv @@ -87,6 +83,7 @@ def __init__(self, vlsvFile, *args, **kwargs): zValues.SetValue(i, x) self.SetZCoordinates(zValues) + self.__descriptor = "" # Here we need to actually init the hypertreegrid, # and at least construct the mapping to leaf nodes @@ -118,14 +115,14 @@ def __init__(self, vlsvFile, *args, **kwargs): f._VlsvReader__read_fileindex_for_cellid() - cidArray = vtk.vtkDoubleArray() - cidArray.SetName('CellID') - cidArray.SetNumberOfValues(0) - self.GetCellData().AddArray(cidArray) + # cidArray = vtk.vtkDoubleArray() + # cidArray.SetName('CellID') + # cidArray.SetNumberOfValues(0) + # self.GetCellData().AddArray(cidArray) self.fileIndexArray = vtk.vtkDoubleArray() self.fileIndexArray.SetName('fileIndex') - self.fileIndexArray.SetNumberOfValues(0) + self.GetCellData().AddArray(self.fileIndexArray) f.read_variable_to_cache("CellID") @@ -136,9 +133,10 @@ def __init__(self, vlsvFile, *args, **kwargs): xmin,xmax,ymin,ymax,zmin,zmax = f._VlsvReader__xmin,f._VlsvReader__xmax,f._VlsvReader__ymin,f._VlsvReader__ymax,f._VlsvReader__zmin,f._VlsvReader__zmax, fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid - delta = [[0,0,0],[1,0,0],[0,1,0],[1,1,0], - [0,0,1],[1,0,1],[0,1,1],[1,1,1]] + delta = np.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0], + [0,0,1],[1,0,1],[0,1,1],[1,1,1]],dtype=np.int32) + # @jit(nopython=True) def get_cell_indices(cellids, reflevels): ''' Single-cellid specialization with external constants. Returns a given cell's indices as a numpy array @@ -163,7 +161,7 @@ def get_cell_indices(cellids, reflevels): # print(cellids, reflevels, cid_offsets) cellids = cellids - 1 - cid_offsets[reflevels] - cellindices = np.full(3, -1) + cellindices = np.full(3, -1,dtype=np.int32) cellindices[0] = (cellids)%(xcells[reflevels]) cellindices[1] = ((cellids)//(xcells[reflevels]))%(ycells[reflevels]) cellindices[2] = (cellids)//(xcells[reflevels]*ycells[reflevels]) @@ -176,31 +174,51 @@ def get_cell_indices(cellids, reflevels): for level in range(max_ref_level+1): cell_lengths_levels.append(np.array([(xmax - xmin)/(xcells[level]), (ymax - ymin)/(ycells[level]), - (zmax - zmin)/(zcells[level])])) + (zmax - zmin)/(zcells[level])],np.int32)) + @jit(nopython=True) + def children(cid, level): + # cellind = get_cell_indices(cid, level) # get_cellind here for compilation + cellids = cid - 1 - cid_offsets[level] + cellind = np.full(3, -1,dtype=np.int32) + cellind[0] = (cellids)%(xcells[level]) + cellind[1] = ((cellids)//(xcells[level]))%(ycells[level]) + cellind[2] = (cellids)//(xcells[level]*ycells[level]) + + # cellcoordinates = mins + (cellind + 0.25)*cell_lengths_levels[level] + # cellind = np.int64(np.divide((cellcoordinates - mins),cell_lengths_levels[level+1])) + cellind = cellind*2 + out = np.zeros((8,),dtype=np.int64) + # c = np.array((0,1,2,3,4,5,6,7)) + out[:] = cid_offsets[level+1] + (cellind[0] + delta[:,0]) + xcells[level+1]*(cellind[1] + delta[:,1]) + (cellind[2] + delta[:,2])*xcells[level+1]*ycells[level+1] + 1 + # out2 = [int(cid_offsets[level+1] + (cellind + delta[c])[0] + xcells[level+1]*(cellind + delta[c])[1] + (cellind + delta[c])[2]*xcells[level+1]*ycells[level+1] + 1) for c in range(8)] + # print(out, out2) + return out - def refine_cell(cid, level): - cellind = get_cell_indices(cid, level) - cellcoordinates = mins + (cellind + 0.25)*cell_lengths_levels[level] + def refine_cell(cid, level): - cellind = np.array(np.divide((cellcoordinates - mins),cell_lengths_levels[level+1]), dtype = np.int32) + # cellind = get_cell_indices(cid, level) + # cellcoordinates = mins + (cellind + 0.25)*cell_lengths_levels[level] + # cellind = np.array(np.divide((cellcoordinates - mins),cell_lengths_levels[level+1]), dtype = np.int32) - for c in range(8): + childs = children(cid,level) + for c,cellid in enumerate(childs): + # print(c) cursor.ToChild(c) # cell 0/0 # print(cellind, delta[c]) - ci = cellind + delta[c] + # ci = cellind + delta[c] # print(ci) - cellid = int(cid_offsets[level+1] + ci[0] + xcells[level+1]*ci[1] + ci[2]*xcells[level+1]*ycells[level+1] + 1) + # cellid = int(cid_offsets[level+1] + ci[0] + xcells[level+1]*ci[1] + ci[2]*xcells[level+1]*ycells[level+1] + 1) # print(f.get_cell_indices(cid, reflevels = level)) if (cellid in fileindex_for_cellid.keys()): # Assigns an index for the value of the cell pointed to by the current cursor state idx = cursor.GetGlobalNodeIndex() # print(idx, cellid) - self.GetCellData().SetActiveScalars('CellID') - cidArray.InsertTuple1(idx, cellid) - self.GetCellData().SetActiveScalars('fileIndex') - self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) + # self.GetCellData().SetActiveScalars('CellID') + # cidArray.InsertTuple1(idx, cellid) + # self.GetCellData().SetActiveScalars('fileIndex') + # self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) # numvalues += 1 self.idxToFileIndex[idx] = fileindex_for_cellid[cellid] # ave += np.array([var_x(cellid, 'proton/vg_v'),*var(cellid, 'vg_b_vol')]) @@ -217,7 +235,7 @@ def refine_cell(cid, level): #vectorArray.InsertTuple3(idx, *refval[1:]) cursor.ToParent() # cell 0 - idx = cursor.GetGlobalNodeIndex() + # idx = cursor.GetGlobalNodeIndex() return #ave/8 @@ -227,11 +245,11 @@ def handle_cell(cellid): if (cellid in fileindex_for_cellid.keys()): # Assigns an index for the value of the cell pointed to by the current cursor state idx = cursor.GetGlobalNodeIndex() - self.GetCellData().SetActiveScalars('CellID') + # self.GetCellData().SetActiveScalars('CellID') # print(idx, cid) - cidArray.InsertTuple1(idx, cellid) - self.GetCellData().SetActiveScalars('fileIndex') - self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) + # cidArray.InsertTuple1(idx, cellid) + # self.GetCellData().SetActiveScalars('fileIndex') + # self.fileIndexArray.InsertTuple1(idx, fileindex_for_cellid[cellid]) # numvalues += 1 self.idxToFileIndex[idx] = fileindex_for_cellid[cellid] # Cell not in list -> must be refined (assume cid is in basegrid) @@ -269,13 +287,126 @@ def loop(): import pstats from pstats import SortKey - with cProfile.Profile() as pr: - loop() - ps = pstats.Stats(pr).sort_stats().reverse_order() - ps.print_stats() + print("Walking HTG") + # with cProfile.Profile() as pr: + loop() + # print("loop done") + # writer = vtk.vtkXMLHyperTreeGridWriter() + # writer.SetFileName("output_FID_1.htg") + + # writer.SetInputPort(src.GetOutputPort()) + # self.fileIndexArray.SetNumberOfValues(1) + # self.fileIndexArray.SetNumberOfTuples(len(self.idxToFileIndex)) + for idx,fileIndex in self.idxToFileIndex.items(): + self.fileIndexArray.InsertTuple1(idx, fileIndex) + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() + # ps.print_stats() + # self >> writer + # writer.Write() + - # loop() + + + print("Adding CellID array") + with cProfile.Profile() as pr: + self.addArrayFromVlsv("CellID") + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() + # ps.print_stats() + if False: + + from io import StringIO + descr = StringIO() + print("Building descriptor") + subdivided = [[] for l in range(max_ref_level+1)] + with cProfile.Profile() as pr: + for c in range(1,int(np.prod(basegridsize))+1): + if c in fileindex_for_cellid.keys(): + # self.__descriptor += "." + descr.write(".") + else: + # self.__descriptor += "R" + descr.write("R") + subdivided[0].append(c) + + # self.__descriptor += "|" + descr.write("|") + for l in range(1,max_ref_level+1): + for c in subdivided[l-1]: + for child in children(c, l-1): + if child in fileindex_for_cellid.keys(): + # self.__descriptor += "." + descr.write(".") + else: + # self.__descriptor += "R" + descr.write("R") + subdivided[l].append(child) + if l < max_ref_level: + # self.__descriptor += "|" + descr.write("|") + self.__descriptor = descr.getvalue() + ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() + ps.print_stats() + # self.SetDimensions([len(nodeCoordinatesX),len(nodeCoordinatesY),len(nodeCoordinatesZ)]) + # self.SetBranchFactor(2) + print("Building with vtkHyperTreeGridSource and descriptor, len ", len(self.__descriptor)) + with cProfile.Profile() as pr: + src = vtk.vtkHyperTreeGridSource(max_depth = max_ref_level, + dimensions=(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)), + grid_scale = (f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz), + branch_factor = 2, + descriptor = self.__descriptor) + + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() + # ps.print_stats() + + # ahtg = src.GetOutput() + print("foo") + writer = vtk.vtkXMLHyperTreeGridWriter() + writer.SetFileName("output_FID_2.htg") + + # # writer.SetInputPort(src.GetOutputPort()) + src.Update() + ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() + ps.print_stats() + # src >> writer + # writer.Write() + # print("Stats after write") + + # Hyper tree grid to unstructured grid filter. + # htg2ug = vtk.vtkHyperTreeGridToUnstructuredGrid() + # # data = src.RequestData() + + + # shrink = vtk.vtkShrinkFilter(shrink_factor=0.5) + + # mapper = vtk.vtkDataSetMapper(scalar_visibility=False) + + # # src >> htg2ug >> shrink >> mapper + # src >> htg2ug >> mapper + + # actor = vtk.vtkActor(mapper=mapper) + # # actor.property.diffuse_color = vtk.colors.GetColor3d('Burlywood') + + # renderer = vtk.vtkRenderer()#background=vtk.colors.GetColor3d('SlateGray')) + # render_window = vtk.vtkRenderWindow(size=(640, 480), window_name='HyperTreeGridSource') + # render_window.AddRenderer(renderer) + # interactor = vtk.vtkRenderWindowInteractor() + # interactor.render_window = render_window + + # renderer.AddActor(actor) + # renderer.ResetCamera() + # renderer.active_camera.Azimuth(150) + # renderer.active_camera.Elevation(30) + # renderer.ResetCameraClippingRange() + + # render_window.Render() + # interactor.Start() + # render_window.Close() + + + + # print(ahtg) def findVariablesFromVlsv(self, getReducers = False): @@ -411,8 +542,8 @@ def RequestDataObject(self, request, inInfo, outInfo): def RequestInformation(self, request, inInfo, outInfo): info = outInfo.GetInformationObject(0) - print("VlsvVtkReader RequestInformation:") - print(info) + # print("VlsvVtkReader RequestInformation:") + # print(info) if self.__htg is None: self.__htg = vtkVlsvHyperTreeGrid(self.__reader) vars = self.__htg.findVariablesFromVlsv(getReducers=True) @@ -420,6 +551,7 @@ def RequestInformation(self, request, inInfo, outInfo): if "vg" in name: # print(name) self._arrayselection.AddArray(name) + self._arrayselection.DisableArray(name) self.__cellarrays.append(name) # self.__htg.addArrayFromVlsv("vg_beta_star") # self.__htg.GetCellData().SetActiveScalars("vg_beta_star") @@ -436,8 +568,8 @@ def RequestUpdateExtent(self, request, inInfo, outInfo): def RequestData(self, request, inInfo, outInfo): - print("VlsvVtkReader RequestData:") - print(outInfo) + # print("VlsvVtkReader RequestData:") + # print(outInfo) if self.__htg is None: self.__htg = vtkVlsvHyperTreeGrid(self.__reader) @@ -472,8 +604,9 @@ def __main__(): import pytools as pt # This initializes a hypertreegrid from the given reader. reader = VlsvVtkReader() - reader.SetFileName("/home/mjalho/Downloads/bulk.0002189.vlsv") - reader.ReadAllScalarsOn() + # reader.SetFileName("/home/mjalho/Downloads/bulk.0002189.vlsv") + reader.SetFileName("/home/mjalho/bulk1.0001418.vlsv") + # reader.ReadAllScalarsOn() reader.Update() cf = vtk.vtkHyperTreeGridContour() cf.SetInputConnection(reader.GetOutputPort()) @@ -498,7 +631,7 @@ def __main__(): time.sleep(10) htg = reader.GetOutputDataObject(0) - print(htg) + # print(htg) # htg.__vtkVlsvHyperTreeGrid_VlsvAddArray('vg_b_vol') # htg = reader.GetOutputData() # print(htg) @@ -508,7 +641,7 @@ def __main__(): # htg.addArrayFromVlsv("vg_beta_star") writer = vtk.vtkXMLHyperTreeGridWriter() - writer.SetFileName("output_EGE.htg") + writer.SetFileName("output_FID.htg") writer.SetInputData(htg) writer.Write() diff --git a/requirements.txt b/requirements.txt index a2682ed3..f35d2899 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,5 @@ numpy scipy matplotlib scikit-image -vtk +vtk>=9.2 +numba \ No newline at end of file From e188bceebfa085be630da3b98250c04933d3ea70 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Wed, 25 Dec 2024 15:32:02 +0200 Subject: [PATCH 10/18] No reloading of array --- pyVlsv/vlsvvtkinterface.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 20bd0c50..b3b8f22e 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -289,7 +289,8 @@ def loop(): from pstats import SortKey print("Walking HTG") # with cProfile.Profile() as pr: - loop() + if True: + loop() # print("loop done") # writer = vtk.vtkXMLHyperTreeGridWriter() # writer.SetFileName("output_FID_1.htg") @@ -297,8 +298,8 @@ def loop(): # writer.SetInputPort(src.GetOutputPort()) # self.fileIndexArray.SetNumberOfValues(1) # self.fileIndexArray.SetNumberOfTuples(len(self.idxToFileIndex)) - for idx,fileIndex in self.idxToFileIndex.items(): - self.fileIndexArray.InsertTuple1(idx, fileIndex) + for idx,fileIndex in self.idxToFileIndex.items(): + self.fileIndexArray.InsertTuple1(idx, fileIndex) # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() # self >> writer @@ -307,9 +308,9 @@ def loop(): - print("Adding CellID array") - with cProfile.Profile() as pr: - self.addArrayFromVlsv("CellID") + print("Adding CellID array") + with cProfile.Profile() as pr: + self.addArrayFromVlsv("CellID") # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() @@ -355,7 +356,7 @@ def loop(): dimensions=(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)), grid_scale = (f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz), branch_factor = 2, - descriptor = self.__descriptor) + descriptor = descr.getvalue())#self.__descriptor) # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() @@ -430,6 +431,11 @@ def findVariablesFromVlsv(self, getReducers = False): that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. ''' def addArrayFromVlsv(self, varname): + + # Do not re-add an already existing array + if self.GetCellData().HasArray(varname): + return True + array = vtk.vtkDoubleArray() # cidArray2.DeepCopy(fileIndexArray) data = self.vlsvreader.read_variable(varname) From e239c641581770a35be35bf746106c28e493ef7b Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Thu, 2 Jan 2025 13:47:58 +0200 Subject: [PATCH 11/18] Saving descriptor as a metadata file and fast reconstruction from that --- miscellaneous/rankine.py | 2 +- pyCalculations/gyrophaseangle.py | 6 +- pyCalculations/themis_observation.py | 2 +- pyCalculations/timeevolution.py | 2 +- pyPlots/plot_variables.py | 4 +- pyVlsv/reduction.py | 2 +- pyVlsv/vlsvVtkParaReader.py | 2 + pyVlsv/vlsvvtkinterface.py | 230 +++++++++++++++++++-------- 8 files changed, 173 insertions(+), 77 deletions(-) diff --git a/miscellaneous/rankine.py b/miscellaneous/rankine.py index e1927723..bef3988e 100644 --- a/miscellaneous/rankine.py +++ b/miscellaneous/rankine.py @@ -23,7 +23,7 @@ # Rankine-Hugoniot conditions to determine the right state from the left state of an MHD oblique shock import numpy as np -import pylab as pl +import matplotlib.pyplot as pl import sys from scipy import optimize from cutthrough import cut_through diff --git a/pyCalculations/gyrophaseangle.py b/pyCalculations/gyrophaseangle.py index 698bd3c5..a20de01a 100644 --- a/pyCalculations/gyrophaseangle.py +++ b/pyCalculations/gyrophaseangle.py @@ -22,7 +22,7 @@ # import numpy as np -import pylab as pl +import matplotlib.pyplot as pl from rotation import rotateVectorToVector import logging @@ -39,7 +39,7 @@ def gyrophase_angles_from_file( vlsvReader, cellid): vlsvReader = VlsvReader("fullf.0001.vlsv") result = gyrophase_angles_from_file( vlsvReader=vlsvReader, cellid=1924) # Plot the data - import pylab as pl + import matplotlib.pyplot as pl pl.hist(result[0].data, weights=result[1].data, bins=100, log=False) ''' # Read the velocity cells: @@ -95,7 +95,7 @@ def gyrophase_angles(bulk_velocity, B_unit, velocity_cell_data, velocity_coordin vlsvReader = VlsvReader("fullf.0001.vlsv") result = gyrophase_angles_from_file( vlsvReader=vlsvReader, cellid=1924, cosine=True, plasmaframe=False ) # Plot the data - import pylab as pl + import matplotlib.pyplot as pl pl.hist(result[0].data, weights=result[1].data, bins=100, log=False) ''' diff --git a/pyCalculations/themis_observation.py b/pyCalculations/themis_observation.py index 9b41f454..cf46a46a 100644 --- a/pyCalculations/themis_observation.py +++ b/pyCalculations/themis_observation.py @@ -22,7 +22,7 @@ # import numpy as np -import pylab as pl +import matplotlib.pyplot as pl import matplotlib from rotation import rotateVectorToVector from scipy.interpolate import griddata diff --git a/pyCalculations/timeevolution.py b/pyCalculations/timeevolution.py index 26837ac6..1fb9d905 100644 --- a/pyCalculations/timeevolution.py +++ b/pyCalculations/timeevolution.py @@ -38,7 +38,7 @@ def cell_time_evolution( vlsvReader_list, variables, cellids, units="" ): .. code-block:: python - import pytools as pt; import pylab as pl + import pytools as pt; import matplotlib.pyplot as pl # Example of usage: time_data = pt.calculations.cell_time_evolution( vlsvReader_list=[VlsvReader("bulk.000.vlsv"), VlsvReader("bulk.001.vlsv"), VlsvReader("bulk.002.vlsv")], variables=["rho", "Pressure", "B"], cellids=[2,4], units=["N", "Pascal", "T"] ) diff --git a/pyPlots/plot_variables.py b/pyPlots/plot_variables.py index 7f8421de..a9e81725 100644 --- a/pyPlots/plot_variables.py +++ b/pyPlots/plot_variables.py @@ -29,7 +29,7 @@ ############################### -import pylab as pl +import matplotlib.pyplot as pl import logging import numpy as np from matplotlib.ticker import MaxNLocator @@ -135,7 +135,7 @@ def plot_multiple_variables( variables_x_list, variables_y_list, figure=[], clea length_of_list = len(variables_x_list) if figure != []: - fig = pl.figure + fig = pl.figure() if len(fig.get_axes()) < length_of_list: for i in (np.arange(length_of_list-len(fig.get_axes())) + len(fig.get_axes())): fig.add_subplot(length_of_list,1,i) diff --git a/pyVlsv/reduction.py b/pyVlsv/reduction.py index a02a82f8..5fcd149c 100644 --- a/pyVlsv/reduction.py +++ b/pyVlsv/reduction.py @@ -25,7 +25,7 @@ ''' import logging import numpy as np -import pylab as pl +import matplotlib.pyplot as pl from reducer import DataReducerVariable from rotation import rotateTensorToVector, rotateArrayTensorToVector from gyrophaseangle import gyrophase_angles diff --git a/pyVlsv/vlsvVtkParaReader.py b/pyVlsv/vlsvVtkParaReader.py index 49225a2b..5fe07521 100644 --- a/pyVlsv/vlsvVtkParaReader.py +++ b/pyVlsv/vlsvVtkParaReader.py @@ -35,6 +35,8 @@ def __init__(self): @smhint.filechooser(extensions="vlsv", file_description="Vlasiator VLSV files") def SetFileName(self, filename): """Specify filename for the file to read.""" + if filename is None or filename == "None": + return print(filename) super().SetFileName(filename) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index b3b8f22e..b8fd6e7f 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -24,16 +24,22 @@ import logging import numpy as np import cProfile +import os.path +import pickle import vtk.numpy_interface import vtk.numpy_interface.dataset_adapter +import vtk.vtkCommonColor +import vtkmodules.vtkCommonColor +import vtkmodules.vtkCommonColor +import vtkmodules.vtkCommonColor import pytools as pt try: import vtk except Exception as e: logging.error("VTK import did not succeed") raise(e) - +import time from numba import jit '''Extend vtkHyperTreeGrid with a vlsvReader wrapper. @@ -53,8 +59,10 @@ class vtkVlsvHyperTreeGrid(vtk.vtkHyperTreeGrid): def __init__(self, vlsvFile, *args, **kwargs): self.vlsvreader = vlsvFile #VlsvReader(vlsvFile) f = self.vlsvreader - vtk.vtkHyperTreeGrid.__init__(self, args, kwargs) + vtk.vtkHyperTreeGrid.__init__(self)#, args, kwargs) self.Initialize() + self.metafile = self.vlsvreader.file_name[:-5]+"_meta" + basegridsize = f.get_spatial_mesh_size() ext = f.get_fsgrid_mesh_extent() nodeCoordinatesX = f.read(tag="MESH_NODE_CRDS_X", mesh='SpatialGrid') @@ -175,6 +183,8 @@ def get_cell_indices(cellids, reflevels): cell_lengths_levels.append(np.array([(xmax - xmin)/(xcells[level]), (ymax - ymin)/(ycells[level]), (zmax - zmin)/(zcells[level])],np.int32)) + + @jit(nopython=True) def children(cid, level): # cellind = get_cell_indices(cid, level) # get_cellind here for compilation @@ -287,19 +297,44 @@ def loop(): import pstats from pstats import SortKey - print("Walking HTG") + # with cProfile.Profile() as pr: - if True: + reinit = True + + try: + if reinit: + raise Exception("reinit = True") + with open(self.metafile,'rb') as mfile: + data = pickle.load(mfile) + # print(data) + self.idxToFileIndex = data['idxToFileIndexMap'] + self.__descriptor = data['descr'] + except Exception as e: + print("Re-initializing HTG, no metadata accessible because of ", e) + reinit = True + + + if False: + t0 = time.time() + print("Walking HTG") loop() - # print("loop done") - # writer = vtk.vtkXMLHyperTreeGridWriter() - # writer.SetFileName("output_FID_1.htg") + print("loop done") + writer = vtk.vtkXMLHyperTreeGridWriter() + writer.SetFileName("output_FID_1.htg") + self.fileIndexArray.SetNumberOfValues(1) + self.fileIndexArray.SetNumberOfTuples(len(self.idxToFileIndex)) + for idx,fileIndex in self.idxToFileIndex.items(): + self.fileIndexArray.InsertTuple1(idx, fileIndex) + writer.SetInputData(self) + writer.Write() - # writer.SetInputPort(src.GetOutputPort()) - # self.fileIndexArray.SetNumberOfValues(1) - # self.fileIndexArray.SetNumberOfTuples(len(self.idxToFileIndex)) + + if not reinit: for idx,fileIndex in self.idxToFileIndex.items(): self.fileIndexArray.InsertTuple1(idx, fileIndex) + # else: + + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() # self >> writer @@ -307,56 +342,88 @@ def loop(): - - print("Adding CellID array") - with cProfile.Profile() as pr: - self.addArrayFromVlsv("CellID") + t1 = time.time() + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() + # print("Actual time elapsed:", time.time()-t0, "of which", time.time()-t1, "spent adding an array") - if False: - + print("pre-calling the jit function to pre-compile") + children(1,0) + self.idxToCellID = {} + self.idxToFileIndex = {} + if reinit: + t0 = time.time() from io import StringIO descr = StringIO() print("Building descriptor") subdivided = [[] for l in range(max_ref_level+1)] - with cProfile.Profile() as pr: - for c in range(1,int(np.prod(basegridsize))+1): - if c in fileindex_for_cellid.keys(): - # self.__descriptor += "." - descr.write(".") - else: - # self.__descriptor += "R" - descr.write("R") - subdivided[0].append(c) - - # self.__descriptor += "|" - descr.write("|") - for l in range(1,max_ref_level+1): - for c in subdivided[l-1]: - for child in children(c, l-1): - if child in fileindex_for_cellid.keys(): - # self.__descriptor += "." - descr.write(".") - else: - # self.__descriptor += "R" - descr.write("R") - subdivided[l].append(child) - if l < max_ref_level: - # self.__descriptor += "|" - descr.write("|") - self.__descriptor = descr.getvalue() - ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() - ps.print_stats() + idx = 0 + # with cProfile.Profile() as pr: + for c in range(1,int(np.prod(basegridsize))+1): + if c in fileindex_for_cellid.keys(): + # self.__descriptor += "." + descr.write(".") + self.idxToFileIndex[idx] = fileindex_for_cellid[c] + self.idxToCellID[idx] = c + else: + # self.__descriptor += "R" + descr.write("R") + subdivided[0].append(c) + idx += 1 + + # self.__descriptor += "|" + descr.write("|") + for l in range(1,max_ref_level+1): + for c in subdivided[l-1]: + for child in children(c, l-1): + if child in fileindex_for_cellid.keys(): + # self.__descriptor += "." + descr.write(".") + self.idxToFileIndex[idx] = fileindex_for_cellid[child] + self.idxToCellID[idx] = c + else: + # self.__descriptor += "R" + descr.write("R") + subdivided[l].append(child) + idx += 1 + if l < max_ref_level: + # self.__descriptor += "|" + descr.write("|") + self.__descriptor = descr.getvalue() + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE) + # ps.print_stats() + print("Actual time spent:",time.time()-t0) # self.SetDimensions([len(nodeCoordinatesX),len(nodeCoordinatesY),len(nodeCoordinatesZ)]) # self.SetBranchFactor(2) print("Building with vtkHyperTreeGridSource and descriptor, len ", len(self.__descriptor)) + + if reinit: + with open(self.metafile,'wb') as mfile: + pickle.dump({"descr":self.__descriptor, "idxToFileIndexMap":self.idxToFileIndex}, mfile) + + if True: with cProfile.Profile() as pr: - src = vtk.vtkHyperTreeGridSource(max_depth = max_ref_level, - dimensions=(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)), - grid_scale = (f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz), - branch_factor = 2, - descriptor = descr.getvalue())#self.__descriptor) + # src = vtk.vtkHyperTreeGridSource(max_depth = max_ref_level+1, + # dimensions=(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)), + # grid_scale = (f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz), + # branch_factor = 2,) + src = vtk.vtkHyperTreeGridSource() + src.SetMaxDepth(max_ref_level+1) + src.SetDimensions(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)) + src.SetGridScale(f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz) + src.SetBranchFactor(2) + self.__descriptor = src.ConvertDescriptorStringToBitArray(self.__descriptor) + src.SetDescriptorBits(self.__descriptor) + # src.SetDescriptor(self.__descriptor) + + src.Update() + + htg = src.GetHyperTreeGridOutput() + htg.SetXCoordinates(xValues) + htg.SetYCoordinates(yValues) + htg.SetZCoordinates(zValues) + # ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() # ps.print_stats() @@ -365,48 +432,69 @@ def loop(): print("foo") writer = vtk.vtkXMLHyperTreeGridWriter() writer.SetFileName("output_FID_2.htg") + writer.SetInputData(htg) - # # writer.SetInputPort(src.GetOutputPort()) - src.Update() + print("foois") + # src >> writer + writer.Write() + print("Stats after write") ps = pstats.Stats(pr).sort_stats(SortKey.CUMULATIVE).reverse_order() ps.print_stats() - # src >> writer - # writer.Write() - # print("Stats after write") + print("foo2") + self.CopyStructure(htg) + + - # Hyper tree grid to unstructured grid filter. + # # # Hyper tree grid to unstructured grid filter. # htg2ug = vtk.vtkHyperTreeGridToUnstructuredGrid() - # # data = src.RequestData() + # # # # data = src.RequestData() - # shrink = vtk.vtkShrinkFilter(shrink_factor=0.5) + # # # shrink = vtk.vtkShrinkFilter(shrink_factor=0.5) - # mapper = vtk.vtkDataSetMapper(scalar_visibility=False) + # mapper = vtk.vtkDataSetMapper() - # # src >> htg2ug >> shrink >> mapper + # # # # src >> htg2ug >> shrink >> mapper # src >> htg2ug >> mapper + # # htg2ug.SetInputData(src.GetHyperTreeGridOutput()) + # # htg2ug.SetInputData(self) + # # mapper.SetInputData(htg2ug.GetOutput()) - # actor = vtk.vtkActor(mapper=mapper) - # # actor.property.diffuse_color = vtk.colors.GetColor3d('Burlywood') + # actor = vtk.vtkActor()#(mapper=mapper) + # actor.SetMapper(mapper) + + # #actor.property.diffuse_color = vtk.colors.GetColor3d('Burlywood') + # import vtkmodules.vtkCommonColor as vtkCommonColor + # colors = vtkCommonColor.vtkNamedColors() + + # actor.property.diffuse_color = colors.GetColor3d('Burlywood') # renderer = vtk.vtkRenderer()#background=vtk.colors.GetColor3d('SlateGray')) - # render_window = vtk.vtkRenderWindow(size=(640, 480), window_name='HyperTreeGridSource') + # renderer.background = colors.GetColor3d('SlateGray') + # render_window = vtk.vtkRenderWindow()#(size=(640, 480), window_name='HyperTreeGridSource') + # render_window.SetSize(640,480) # render_window.AddRenderer(renderer) # interactor = vtk.vtkRenderWindowInteractor() # interactor.render_window = render_window # renderer.AddActor(actor) # renderer.ResetCamera() - # renderer.active_camera.Azimuth(150) - # renderer.active_camera.Elevation(30) + # renderer.GetActiveCamera().Azimuth(150) + # renderer.GetActiveCamera().Elevation(30) # renderer.ResetCameraClippingRange() # render_window.Render() # interactor.Start() + # time.sleep(20) # render_window.Close() - + # self. + # if not reinit: + for idx,fileIndex in self.idxToFileIndex.items(): + self.fileIndexArray.InsertTuple1(idx, fileIndex) + # self.addArrayFromVlsv("CellID") + # print(ahtg) @@ -426,6 +514,8 @@ def findVariablesFromVlsv(self, getReducers = False): return vars + def getDecriptor(): + return self.__descriptor '''This function adds one SpatialGrid variable from the reader object and maps that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. @@ -434,6 +524,7 @@ def addArrayFromVlsv(self, varname): # Do not re-add an already existing array if self.GetCellData().HasArray(varname): + print("skipped existing array") return True array = vtk.vtkDoubleArray() @@ -527,7 +618,8 @@ def SetFileName(self, filename): if filename != self.__FileName: self.Modified() self.__FileName = filename - self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) + if self.__FileName is not None: + self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) def GetFileName(self): @@ -551,6 +643,7 @@ def RequestInformation(self, request, inInfo, outInfo): # print("VlsvVtkReader RequestInformation:") # print(info) if self.__htg is None: + print("Creating htg via RequestInformation") self.__htg = vtkVlsvHyperTreeGrid(self.__reader) vars = self.__htg.findVariablesFromVlsv(getReducers=True) for name in vars: @@ -578,6 +671,7 @@ def RequestData(self, request, inInfo, outInfo): # print(outInfo) if self.__htg is None: + print("Creating htg via RequestData") self.__htg = vtkVlsvHyperTreeGrid(self.__reader) for name in self.__cellarrays: @@ -633,8 +727,8 @@ def __main__(): renWin.SetSize(600, 600) # renWin.Render() - import time - time.sleep(10) + + # time.sleep(10) htg = reader.GetOutputDataObject(0) # print(htg) From daa62dd31f5d7b59a8b07bbc1cf33ea95bbe3a0d Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Thu, 2 Jan 2025 13:58:32 +0200 Subject: [PATCH 12/18] Quick logic fixes --- pyVlsv/vlsvvtkinterface.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index b8fd6e7f..4a3e70aa 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -128,10 +128,6 @@ def __init__(self, vlsvFile, *args, **kwargs): # cidArray.SetNumberOfValues(0) # self.GetCellData().AddArray(cidArray) - self.fileIndexArray = vtk.vtkDoubleArray() - self.fileIndexArray.SetName('fileIndex') - - self.GetCellData().AddArray(self.fileIndexArray) f.read_variable_to_cache("CellID") @@ -299,7 +295,7 @@ def loop(): from pstats import SortKey # with cProfile.Profile() as pr: - reinit = True + reinit = False try: if reinit: @@ -329,9 +325,9 @@ def loop(): writer.Write() - if not reinit: - for idx,fileIndex in self.idxToFileIndex.items(): - self.fileIndexArray.InsertTuple1(idx, fileIndex) + # if not reinit: + # for idx,fileIndex in self.idxToFileIndex.items(): + # self.fileIndexArray.InsertTuple1(idx, fileIndex) # else: @@ -350,9 +346,10 @@ def loop(): print("pre-calling the jit function to pre-compile") children(1,0) - self.idxToCellID = {} - self.idxToFileIndex = {} + if reinit: + self.idxToCellID = {} + self.idxToFileIndex = {} t0 = time.time() from io import StringIO descr = StringIO() @@ -491,6 +488,14 @@ def loop(): # self. # if not reinit: + self.fileIndexArray = vtk.vtkDoubleArray() + self.fileIndexArray.SetName('fileIndex') + + self.GetCellData().AddArray(self.fileIndexArray) + + # self.fileIndexArray.SetNumberOfValues(1) + # self.fileIndexArray.SetNumberOfTuples(len(self.idxToFileIndex)) + for idx,fileIndex in self.idxToFileIndex.items(): self.fileIndexArray.InsertTuple1(idx, fileIndex) # self.addArrayFromVlsv("CellID") From 066dd63cd3b568a567cba488b7fe24de7a359dc5 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Sat, 4 Jan 2025 14:31:05 +0200 Subject: [PATCH 13/18] Use the VlsvVtkReader class instead for wrapping the HTG with vtkHyperTreeGridSource, moved the Paraview plugin to a new folder --- .../vlsvParaReader.py | 9 +- pyVlsv/vlsvvtkinterface.py | 235 +++++++++++++++++- 2 files changed, 232 insertions(+), 12 deletions(-) rename pyVlsv/vlsvVtkParaReader.py => plugins/vlsvParaReader.py (95%) diff --git a/pyVlsv/vlsvVtkParaReader.py b/plugins/vlsvParaReader.py similarity index 95% rename from pyVlsv/vlsvVtkParaReader.py rename to plugins/vlsvParaReader.py index 5fe07521..241dd877 100644 --- a/pyVlsv/vlsvVtkParaReader.py +++ b/plugins/vlsvParaReader.py @@ -1,4 +1,9 @@ -"""This module demonstrates various ways of adding +""" A rudimentary wrapper from the vlsvVtkInterface. Add and load from Paraview +plugins menu. + +Built on Kitware example: +https://github.com/Kitware/ParaView/blob/master/Examples/Plugins/PythonAlgorithm/PythonAlgorithmExamples.py +This module demonstrates various ways of adding VTKPythonAlgorithmBase subclasses as filters, sources, readers, and writers in ParaView""" @@ -150,7 +155,7 @@ def RequestData(self, request, inInfo, outInfo): ''' def test_PythonVLSVReader(fname): - reader = PythonVLSVReader() + reader = PythonPvVLSVReader() reader.SetFileName(fname) reader.Update() assert reader.GetOutputDataObject(0).GetNumberOfRows() > 0 diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 4a3e70aa..660168f9 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -41,7 +41,10 @@ raise(e) import time from numba import jit -'''Extend vtkHyperTreeGrid with a vlsvReader wrapper. +'''Extend vtkHyperTreeGrid. +This class is not used right now: makes probably more sense to have +most of the features in the reader class below, esp. since the +vtkHyperTreeGridSource is the fast way of constructing the HTG object. The extension consists of automatic mapping of a vlsv SpatialGrid to the HTG structure. Functions are provided @@ -519,8 +522,7 @@ def findVariablesFromVlsv(self, getReducers = False): return vars - def getDecriptor(): - return self.__descriptor + '''This function adds one SpatialGrid variable from the reader object and maps that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. @@ -609,6 +611,15 @@ def _markmodified(*args, **kwars): # Should implement these, more or less: # https://vtk.org/doc/nightly/html/classvtkXMLHyperTreeGridReader.html +'''Experimental Python wrapper to expose VLSV files as VTK. SpatialGrid supported for now. + +This will create an internal VlsvReader object and cache the HTG descriptor and +a corresponding file layout mapping to disk for repeated accesses: +building the descriptor takes ~30s (with Python; maybe the whole loop could be @jit-ted) + +This class functions like a VTK source. Calling `addArrayFromVlsv` is necessary +to load and map data from the VLSV to HTG. +''' class VlsvVtkReader(VTKPythonAlgorithmBase): def __init__(self): VTKPythonAlgorithmBase.__init__(self, nInputPorts = 0, nOutputPorts=1, outputType='vtkHyperTreeGrid') @@ -625,14 +636,155 @@ def SetFileName(self, filename): self.__FileName = filename if self.__FileName is not None: self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) + self.__metafile = self.__FileName[:-5]+"_meta" + + def GetFileName(self): + return self.__FileName + + def buildDescriptor(self): + f = self.__reader + f._VlsvReader__read_fileindex_for_cellid() + fileindex_for_cellid = f._VlsvReader__fileindex_for_cellid + xc = f._VlsvReader__xcells + yc = f._VlsvReader__ycells + zc = f._VlsvReader__zcells + max_ref_level = f.get_max_refinement_level() + + cid_offsets = np.zeros(max_ref_level+1, dtype=np.int64) + isum = 0 + for i in range(0,max_ref_level): + isum = isum + 2**(3*i) * xc * yc * zc + cid_offsets[i+1] = isum + + + xcells = np.zeros((max_ref_level+1), dtype=np.int64) + ycells = np.zeros((max_ref_level+1), dtype=np.int64) + zcells = np.zeros((max_ref_level+1), dtype=np.int64) + for r in range(max_ref_level+1): + xcells[r] = xc*2**(r) + ycells[r] = yc*2**(r) + zcells[r] = zc*2**(r) + delta = np.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0], + [0,0,1],[1,0,1],[0,1,1],[1,1,1]],dtype=np.int32) + + + idxToFileIndex = {} + from io import StringIO + descr = StringIO() + print("Building descriptor") + subdivided = [[] for l in range(max_ref_level+1)] + idx = 0 + for c in range(1,int(np.prod(f.get_spatial_mesh_size()))+1): + if c in fileindex_for_cellid.keys(): + descr.write(".") + idxToFileIndex[idx] = fileindex_for_cellid[c] + else: + descr.write("R") + subdivided[0].append(c) + idx += 1 + + @jit(nopython=True) + def children(cid, level): + # cellind = get_cell_indices(cid, level) # get_cellind here for compilation + cellids = cid - 1 - cid_offsets[level] + cellind = np.full(3, -1,dtype=np.int32) + cellind[0] = (cellids)%(xcells[level]) + cellind[1] = ((cellids)//(xcells[level]))%(ycells[level]) + cellind[2] = (cellids)//(xcells[level]*ycells[level]) + + cellind = cellind*2 + out = np.zeros((8,),dtype=np.int64) + out[:] = cid_offsets[level+1] + (cellind[0] + delta[:,0]) + xcells[level+1]*(cellind[1] + delta[:,1]) + (cellind[2] + delta[:,2])*xcells[level+1]*ycells[level+1] + 1 + return out + + descr.write("|") + for l in range(1,max_ref_level+1): + for c in subdivided[l-1]: + for child in children(c, l-1): + if child in fileindex_for_cellid.keys(): + descr.write(".") + idxToFileIndex[idx] = fileindex_for_cellid[child] + else: + descr.write("R") + subdivided[l].append(child) + idx += 1 + if l < max_ref_level: + descr.write("|") + + return descr.getvalue(), idxToFileIndex + + + def getDescriptor(self, reinit=False): + + try: + if reinit: + raise Exception("reinit = True") + with open(self.__metafile,'rb') as mfile: + data = pickle.load(mfile) + # print(data) + self.__idxToFileIndex = data['idxToFileIndexMap'] + self.__descriptor = data['descr'] + except Exception as e: + print("Re-initializing HTG, no metadata accessible because of ", e) + self.__descriptor, self.__idxToFileIndex = self.buildDescriptor() + with open(self.__metafile,'wb') as mfile: + pickle.dump({"descr":self.__descriptor, "idxToFileIndexMap":self.idxToFileIndex}, mfile) + + return self.__descriptor, self.__idxToFileIndex + + def getHTG(self): + f = self.__reader + + descr, idxToFileIndex = self.getDescriptor() + basegridsize = f.get_spatial_mesh_size() + + src = vtk.vtkHyperTreeGridSource() + src.SetMaxDepth(f.get_max_refinement_level()+1) + src.SetDimensions(int(basegridsize[0]+1),int(basegridsize[1]+1),int(basegridsize[2]+1)) + src.SetGridScale(f._VlsvReader__dx,f._VlsvReader__dy,f._VlsvReader__dz) + src.SetBranchFactor(2) - def GetFileName(self): - return self.__FileName + #descr = src.ConvertDescriptorStringToBitArray(descr) + #src.SetDescriptorBits(descr) + src.SetDescriptor(descr) + src.Update() + htg = src.GetHyperTreeGridOutput() + xValues = vtk.vtkDoubleArray() + xValues.SetNumberOfValues(int(basegridsize[0]+1)) + nodeCoordinatesX = f.read(tag="MESH_NODE_CRDS_X", mesh='SpatialGrid') + nodeCoordinatesY = f.read(tag="MESH_NODE_CRDS_Y", mesh='SpatialGrid') + nodeCoordinatesZ = f.read(tag="MESH_NODE_CRDS_Z", mesh='SpatialGrid') + for i,x in enumerate(nodeCoordinatesX): + xValues.SetValue(i, x) + htg.SetXCoordinates(xValues) + + yValues = vtk.vtkDoubleArray() + yValues.SetNumberOfValues(int(basegridsize[1]+1)) + for i,x in enumerate(nodeCoordinatesY): + yValues.SetValue(i, x) + htg.SetYCoordinates(yValues) + + zValues = vtk.vtkDoubleArray() + zValues.SetNumberOfValues(int(basegridsize[2]+1)) + for i,x in enumerate(nodeCoordinatesZ): + zValues.SetValue(i, x) + htg.SetZCoordinates(zValues) + + htg.fileIndexArray = vtk.vtkDoubleArray() + htg.fileIndexArray.SetName('fileIndex') + + for idx,fileIndex in self.__idxToFileIndex.items(): + htg.fileIndexArray.InsertTuple1(idx, fileIndex) + htg.GetCellData().AddArray(htg.fileIndexArray) + + + return htg + ''' def FillOutputPortInformation(self, port, info): pass @@ -641,6 +793,16 @@ def RequestDataObject(self, request, inInfo, outInfo): #This is where you can create output data objects if the output DATA_TYPE_NAME() is not a concrete type. pass ''' + def findVariablesFromVlsv(self, getReducers = False): + + vars = self.__reader.get_variables() + if getReducers: + reducers = self.__reader.get_reducers() + vars_set = set(vars) + vars_set.update(reducers) + vars = list(vars_set) + + return vars def RequestInformation(self, request, inInfo, outInfo): @@ -649,8 +811,8 @@ def RequestInformation(self, request, inInfo, outInfo): # print(info) if self.__htg is None: print("Creating htg via RequestInformation") - self.__htg = vtkVlsvHyperTreeGrid(self.__reader) - vars = self.__htg.findVariablesFromVlsv(getReducers=True) + self.__htg = self.getHTG() + vars = self.findVariablesFromVlsv(getReducers=True) for name in vars: if "vg" in name: # print(name) @@ -669,7 +831,60 @@ def RequestInformation(self, request, inInfo, outInfo): def RequestUpdateExtent(self, request, inInfo, outInfo): pass ''' - + + '''This function adds one SpatialGrid variable from the reader object and maps + that to the hypertreegrid object. Variable vector sizes of 1,2,3,4,9 supported. + ''' + def addArrayFromVlsv(self, htg, varname): + + # Do not re-add an already existing array + if htg.GetCellData().HasArray(varname): + print("skipped existing array") + return True + + array = vtk.vtkDoubleArray() + # cidArray2.DeepCopy(fileIndexArray) + data = self.__reader.read_variable(varname) + + test_var = data[0] + if test_var is None: + logging.warning("varname " + varname + "returned None; skipping") + return False + if np.isscalar(test_var): + varlen = 1 + elif test_var.ndim == 1: + varlen = len(test_var) + elif test_var.ndim > 1: + varlen = np.prod(test_var.shape) + else: + logging.warning("Weird output from " + varname) + return False + + array.SetNumberOfComponents(varlen) + array.SetNumberOfTuples(htg.fileIndexArray.GetNumberOfTuples()) + array.SetName(varname) + if varlen == 1: + for idx,fileIndex in self.__idxToFileIndex.items(): + array.SetTuple1(idx, data[fileIndex]) + elif varlen == 2: + for idx,fileIndex in self.__idxToFileIndex.items(): + array.SetTuple2(idx, *data[fileIndex]) + elif varlen == 3: + for idx,fileIndex in self.__idxToFileIndex.items(): + array.SetTuple3(idx, *data[fileIndex]) + elif varlen == 4: + for idx,fileIndex in self.__idxToFileIndex.items(): + array.SetTuple4(idx, *data[fileIndex]) + elif varlen == 9: + for idx,fileIndex in self.__idxToFileIndex.items(): + array.SetTuple9(idx, *np.reshape(data[fileIndex],(9))) + else: + raise RuntimeError("No vtk SetTuple wrapper function for varlen = " + str(varlen)) + + + htg.GetCellData().AddArray(array) + return True + def RequestData(self, request, inInfo, outInfo): # print("VlsvVtkReader RequestData:") @@ -677,11 +892,11 @@ def RequestData(self, request, inInfo, outInfo): if self.__htg is None: print("Creating htg via RequestData") - self.__htg = vtkVlsvHyperTreeGrid(self.__reader) + self.__htg = self.getHTG() for name in self.__cellarrays: if self._arrayselection.ArrayIsEnabled(name): - success = self.__htg.addArrayFromVlsv(name) + success = self.addArrayFromVlsv(self.__htg, name) if not success: self._arrayselection.RemoveArrayByName(name) From 93abb638e41efe0fc833b9a558745beeed1a90de Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Thu, 9 Jan 2025 14:06:24 +0200 Subject: [PATCH 14/18] endlines --- pyVlsv/vlsvvtkinterface.py | 2 +- requirements.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index 660168f9..d8c66612 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -968,4 +968,4 @@ def __main__(): if __name__ == "__main__": - __main__() \ No newline at end of file + __main__() diff --git a/requirements.txt b/requirements.txt index f35d2899..544ce163 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,4 +3,4 @@ scipy matplotlib scikit-image vtk>=9.2 -numba \ No newline at end of file +numba From 6d62b36acfcea8a526a72ec2f3a46fcc6b128899 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Mon, 13 Jan 2025 13:55:35 +0200 Subject: [PATCH 15/18] Some cleanup --- plugins/vlsvParaReader.py | 107 ++----------------------------------- pyVlsv/vlsvfile.py | 2 +- pyVlsv/vlsvreader.py | 4 +- pyVlsv/vlsvvtkinterface.py | 14 ++--- 4 files changed, 14 insertions(+), 113 deletions(-) diff --git a/plugins/vlsvParaReader.py b/plugins/vlsvParaReader.py index 241dd877..e2845579 100644 --- a/plugins/vlsvParaReader.py +++ b/plugins/vlsvParaReader.py @@ -60,112 +60,11 @@ def GetDataArraySelection(self): #------------------------------------------------------------------------------ -# A writer example. +# Todo: some tests #------------------------------------------------------------------------------ -''' -@smproxy.writer(extensions="npz", file_description="NumPy Compressed Arrays", support_reload=False) -@smproperty.input(name="Input", port_index=0) -@smdomain.datatype(dataTypes=["vtkTable"], composite_data_supported=False) -class NumpyWriter(VTKPythonAlgorithmBase): - def __init__(self): - VTKPythonAlgorithmBase.__init__(self, nInputPorts=1, nOutputPorts=0, inputType='vtkTable') - self._filename = None - - @smproperty.stringvector(name="FileName", panel_visibility="never") - @smdomain.filelist() - def SetFileName(self, fname): - """Specify filename for the file to write.""" - if self._filename != fname: - self._filename = fname - self.Modified() - - def RequestData(self, request, inInfoVec, outInfoVec): - from vtkmodules.vtkCommonDataModel import vtkTable - from vtkmodules.numpy_interface import dataset_adapter as dsa - - table = dsa.WrapDataObject(vtkTable.GetData(inInfoVec[0], 0)) - kwargs = {} - for aname in table.RowData.keys(): - kwargs[aname] = table.RowData[aname] - - import numpy - numpy.savez_compressed(self._filename, **kwargs) - return 1 - - def Write(self): - self.Modified() - self.Update() -''' -#------------------------------------------------------------------------------ -# A filter example. -#------------------------------------------------------------------------------ -''' -@smproxy.filter() -@smproperty.input(name="InputTable", port_index=1) -@smdomain.datatype(dataTypes=["vtkTable"], composite_data_supported=False) -@smproperty.input(name="InputDataset", port_index=0) -@smdomain.datatype(dataTypes=["vtkDataSet"], composite_data_supported=False) -class ExampleTwoInputFilter(VTKPythonAlgorithmBase): - def __init__(self): - VTKPythonAlgorithmBase.__init__(self, nInputPorts=2, nOutputPorts=1, outputType="vtkPolyData") - - def FillInputPortInformation(self, port, info): - if port == 0: - info.Set(self.INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet") - else: - info.Set(self.INPUT_REQUIRED_DATA_TYPE(), "vtkTable") - return 1 - - def RequestData(self, request, inInfoVec, outInfoVec): - from vtkmodules.vtkCommonDataModel import vtkTable, vtkDataSet, vtkPolyData - input0 = vtkDataSet.GetData(inInfoVec[0], 0) - input1 = vtkDataSet.GetData(inInfoVec[1], 0) - output = vtkPolyData.GetData(outInfoVec, 0) - # do work - print("Pretend work done!") - return 1 - -@smproxy.filter() -@smproperty.input(name="Input") -@smdomain.datatype(dataTypes=["vtkDataSet"], composite_data_supported=False) -class PreserveInputTypeFilter(VTKPythonAlgorithmBase): - """ - Example filter demonstrating how to write a filter that preserves the input - dataset type. - """ - def __init__(self): - super().__init__(nInputPorts=1, nOutputPorts=1, outputType="vtkDataSet") - - def RequestDataObject(self, request, inInfo, outInfo): - inData = self.GetInputData(inInfo, 0, 0) - outData = self.GetOutputData(outInfo, 0) - assert inData is not None - if outData is None or (not outData.IsA(inData.GetClassName())): - outData = inData.NewInstance() - outInfo.GetInformationObject(0).Set(outData.DATA_OBJECT(), outData) - return super().RequestDataObject(request, inInfo, outInfo) - - def RequestData(self, request, inInfo, outInfo): - inData = self.GetInputData(inInfo, 0, 0) - outData = self.GetOutputData(outInfo, 0) - print("input type =", inData.GetClassName()) - print("output type =", outData.GetClassName()) - assert outData.IsA(inData.GetClassName()) - return 1 -''' def test_PythonVLSVReader(fname): - reader = PythonPvVLSVReader() - reader.SetFileName(fname) - reader.Update() - assert reader.GetOutputDataObject(0).GetNumberOfRows() > 0 + pass if __name__ == "__main__": - #test_PythonSuperquadricSource() - test_PythonVLSVReader("/tmp/data.csv") - - from paraview.detail.pythonalgorithm import get_plugin_xmls - from xml.dom.minidom import parseString - for xml in get_plugin_xmls(globals()): - dom = parseString(xml) - print(dom.toprettyxml(" ","\n")) + pass diff --git a/pyVlsv/vlsvfile.py b/pyVlsv/vlsvfile.py index 2d645172..e8c3d1be 100644 --- a/pyVlsv/vlsvfile.py +++ b/pyVlsv/vlsvfile.py @@ -38,7 +38,7 @@ from vlsvreader import fsDecompositionFromGlobalIds,fsReadGlobalIdsPerRank,fsGlobalIdToGlobalIndex from vlsvwriter import VlsvWriter from vlasiatorreader import VlasiatorReader -from vlsvvtkinterface import vtkVlsvHyperTreeGrid,VlsvVtkReader +from vlsvvtkinterface import VlsvVtkReader from vlsvparticles import VlsvParticles import reduction diff --git a/pyVlsv/vlsvreader.py b/pyVlsv/vlsvreader.py index b856f230..9bdc8843 100644 --- a/pyVlsv/vlsvreader.py +++ b/pyVlsv/vlsvreader.py @@ -2032,9 +2032,9 @@ def get_max_refinement_level(self): if self.__max_spatial_amr_level < 0: # Read the file index for cellid cellids=self.read(mesh="SpatialGrid",name="CellID", tag="VARIABLE") - maxcellid = np.amax([cellids]) + maxcellid = np.int64(np.amax([cellids])) - AMR_count = 0 + AMR_count = np.int64(0) while (maxcellid > 0): maxcellid -= 2**(3*(AMR_count))*(self.__xcells*self.__ycells*self.__zcells) AMR_count += 1 diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index d8c66612..f488c971 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -40,7 +40,7 @@ logging.error("VTK import did not succeed") raise(e) import time -from numba import jit +# from numba import jit '''Extend vtkHyperTreeGrid. This class is not used right now: makes probably more sense to have most of the features in the reader class below, esp. since the @@ -184,7 +184,7 @@ def get_cell_indices(cellids, reflevels): (zmax - zmin)/(zcells[level])],np.int32)) - @jit(nopython=True) + # @jit(nopython=True) def children(cid, level): # cellind = get_cell_indices(cid, level) # get_cellind here for compilation cellids = cid - 1 - cid_offsets[level] @@ -636,7 +636,7 @@ def SetFileName(self, filename): self.__FileName = filename if self.__FileName is not None: self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) - self.__metafile = self.__FileName[:-5]+"_meta" + self.__metafile = os.path.join(os.path.dirname(self.__FileName),"vlsvmeta",self.__FileName[:-5]+"_meta") def GetFileName(self): return self.__FileName @@ -683,7 +683,7 @@ def buildDescriptor(self): subdivided[0].append(c) idx += 1 - @jit(nopython=True) + # @jit(nopython=True) def children(cid, level): # cellind = get_cell_indices(cid, level) # get_cellind here for compilation cellids = cid - 1 - cid_offsets[level] @@ -728,6 +728,8 @@ def getDescriptor(self, reinit=False): except Exception as e: print("Re-initializing HTG, no metadata accessible because of ", e) self.__descriptor, self.__idxToFileIndex = self.buildDescriptor() + if not os.path.isdir(os.path.dirname(self.__metafile)): + os.mkdir(os.path.dirname(self.__metafile)) with open(self.__metafile,'wb') as mfile: pickle.dump({"descr":self.__descriptor, "idxToFileIndexMap":self.idxToFileIndex}, mfile) @@ -924,8 +926,8 @@ def __main__(): import pytools as pt # This initializes a hypertreegrid from the given reader. reader = VlsvVtkReader() - # reader.SetFileName("/home/mjalho/Downloads/bulk.0002189.vlsv") - reader.SetFileName("/home/mjalho/bulk1.0001418.vlsv") + reader.SetFileName("/home/mjalho/Downloads/bulk.0002189.vlsv") + # reader.SetFileName("/home/mjalho/bulk1.0001418.vlsv") # reader.ReadAllScalarsOn() reader.Update() cf = vtk.vtkHyperTreeGridContour() From 0960bab3eba05c7413ab9e16d0eeb672cfe1a69c Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Fri, 17 Jan 2025 14:41:42 +0200 Subject: [PATCH 16/18] new folder for metafile, no numba for now --- pyVlsv/vlsvvtkinterface.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyVlsv/vlsvvtkinterface.py b/pyVlsv/vlsvvtkinterface.py index f488c971..548e952e 100644 --- a/pyVlsv/vlsvvtkinterface.py +++ b/pyVlsv/vlsvvtkinterface.py @@ -636,7 +636,8 @@ def SetFileName(self, filename): self.__FileName = filename if self.__FileName is not None: self.__reader = pt.vlsvfile.VlsvReader(self.__FileName) - self.__metafile = os.path.join(os.path.dirname(self.__FileName),"vlsvmeta",self.__FileName[:-5]+"_meta") + fn = os.path.basename(self.__FileName) + self.__metafile = os.path.join(os.path.dirname(self.__FileName),"vlsvmeta",fn[:-5]+"_meta") def GetFileName(self): return self.__FileName @@ -731,7 +732,7 @@ def getDescriptor(self, reinit=False): if not os.path.isdir(os.path.dirname(self.__metafile)): os.mkdir(os.path.dirname(self.__metafile)) with open(self.__metafile,'wb') as mfile: - pickle.dump({"descr":self.__descriptor, "idxToFileIndexMap":self.idxToFileIndex}, mfile) + pickle.dump({"descr":self.__descriptor, "idxToFileIndexMap":self.__idxToFileIndex}, mfile) return self.__descriptor, self.__idxToFileIndex From 71d8b5beb5f72b91c05438614689961b3b9744e6 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Fri, 17 Jan 2025 14:44:41 +0200 Subject: [PATCH 17/18] Refreshed example --- examples/write_htg.py | 12 +++++++----- requirements.txt | 1 - 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/write_htg.py b/examples/write_htg.py index ad16e5bc..84dc1346 100644 --- a/examples/write_htg.py +++ b/examples/write_htg.py @@ -2,15 +2,17 @@ import pytools as pt # This initializes a hypertreegrid from the given reader. -reader = pt.vlsvfile.VlsvReader("../../bulk.0002189.vlsv") -htg = pt.vlsvfile.vtkVlsvHyperTreeGrid(reader) +reader = VlsvVtkReader() +reader.SetFileName("../../bulk.0002189.vlsv") # These functions grab one SpatialGrid variable and map that to # the hypertreegrid. Variable vector sizes of 1,2,3,4,9 supported. -htg.addArrayFromVlsv("vg_b_vol") -htg.addArrayFromVlsv("vg_beta_star") +reader.addArrayFromVlsv("vg_b_vol") +reader.addArrayFromVlsv("vg_beta_star") + +reader.Update() writer = vtk.vtkXMLHyperTreeGridWriter() writer.SetFileName("output_EGE.htg") -writer.SetInputData(htg) +writer.SetInputData(reader.GetOutputDataObject(0)) writer.Write() diff --git a/requirements.txt b/requirements.txt index 544ce163..5dd55481 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,4 +3,3 @@ scipy matplotlib scikit-image vtk>=9.2 -numba From 926e854dcd66360c088e5753afa5485938840703 Mon Sep 17 00:00:00 2001 From: Alho Markku J Date: Fri, 17 Jan 2025 14:46:50 +0200 Subject: [PATCH 18/18] module path to example --- examples/write_htg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/write_htg.py b/examples/write_htg.py index 84dc1346..ed20e424 100644 --- a/examples/write_htg.py +++ b/examples/write_htg.py @@ -2,7 +2,7 @@ import pytools as pt # This initializes a hypertreegrid from the given reader. -reader = VlsvVtkReader() +reader = pt.vlsvfile.VlsvVtkReader() reader.SetFileName("../../bulk.0002189.vlsv") # These functions grab one SpatialGrid variable and map that to