add comment for idt=0

seissolxdmf/seissolxdmf/seissolxdmf.py

@@ -2,29 +2,33 @@
 import os
 import xml.etree.ElementTree as ET
 
+
 def find_line_number_endtag_xdmf(alines):
     for n, line in enumerate(alines):
-       if '</Xdmf>' in line:
-          return n + 1
+        if "</Xdmf>" in line:
+            return n + 1
     raise ValueError("</Xdmf> not found")
 
+
 class seissolxdmf:
     def __init__(self, xdmfFilename):
         self.xdmfFilename = xdmfFilename
-        with open (xdmfFilename, "r") as fid:
-            lines=fid.readlines()
+        with open(xdmfFilename, "r") as fid:
+            lines = fid.readlines()
         # Remove potential extra content at the end of the file
         nlines = find_line_number_endtag_xdmf(lines)
-        if nlines!=len(lines):
-            print(f'Warning: extra content at the end of {xdmfFilename} detected')
-        file_txt = ' '.join([line for line in lines[0:nlines]])
+        if nlines != len(lines):
+            print(f"Warning: extra content at the end of {xdmfFilename} detected")
+        file_txt = " ".join([line for line in lines[0:nlines]])
         self.tree = ET.ElementTree(ET.fromstring(file_txt))
         self.ndt = self.ReadNdt()
         self.nElements = self.ReadNElements()
 
-    def ReadHdf5DatasetChunk(self, absolute_path, hdf5var, firstElement, nchunk, idt=-1):
-        """ Read block of data in hdf5 format
-        idt!=-1 loads only one time step """
+    def ReadHdf5DatasetChunk(
+        self, absolute_path, hdf5var, firstElement, nchunk, idt=-1
+    ):
+        """Read block of data in hdf5 format
+        idt!=-1 loads only one time step"""
         import h5py
 
         lastElement = firstElement + nchunk
@@ -53,11 +57,13 @@ def GetDtype(self, data_prec, isInt):
             else:
                 return np.dtype("d")
 
-    def ReadSimpleBinaryFile(self, absolute_path, MemDimension, data_prec, isInt, idt=-1):
+    def ReadSimpleBinaryFile(
+        self, absolute_path, MemDimension, data_prec, isInt, idt=-1
+    ):
         """Read block of data in binary format (posix)
-        idt!=-1 loads only one time step 
+        idt!=-1 loads only one time step
         this function is a special case of ReadSimpleBinaryFileChunk
-        but kept for performance reasons """
+        but kept for performance reasons"""
         oneDtMem = True if idt != -1 else False
         data_type = self.GetDtype(data_prec, isInt)
 
@@ -72,10 +78,19 @@ def ReadSimpleBinaryFile(self, absolute_path, MemDimension, data_prec, isInt, id
         fid.close()
         return myData
 
-    def ReadSimpleBinaryFileChunk(self, absolute_path, MemDimension, data_prec, isInt, firstElement, nchunk, idt=-1):
+    def ReadSimpleBinaryFileChunk(
+        self,
+        absolute_path,
+        MemDimension,
+        data_prec,
+        isInt,
+        firstElement,
+        nchunk,
+        idt=-1,
+    ):
         """Read block of data in binary format (posix)
         same as ReadSimpleBinaryFile: but reads a subset of the second dimension
-        idt!=-1 loads only one time step """
+        idt!=-1 loads only one time step"""
         oneDtMem = True if idt != -1 else False
         data_type = self.GetDtype(data_prec, isInt)
 
@@ -93,7 +108,7 @@ def ReadSimpleBinaryFileChunk(self, absolute_path, MemDimension, data_prec, isIn
         return myData
 
     def GetDataLocationPrecisionNElementsMemDimension(self, attribute):
-        """ Common function called by ReadTopologyOrGeometry """
+        """Common function called by ReadTopologyOrGeometry"""
         root = self.tree.getroot()
         for Property in root.findall(".//%s" % (attribute)):
             nElements = int(Property.get("NumberOfElements"))
@@ -106,7 +121,7 @@ def GetDataLocationPrecisionNElementsMemDimension(self, attribute):
         return [dataLocation, data_prec, nElements, MemDimension]
 
     def GetDataLocationPrecisionMemDimension(self, dataName):
-        """ Common function called by ReadData """
+        """Common function called by ReadData"""
 
         def get(prop):
             dataLocation = prop.text
@@ -131,9 +146,14 @@ def get(prop):
         raise NameError("%s not found in dataset" % (dataName))
 
     def ReadTopologyOrGeometry(self, attribute):
-        """ Common function to read either connect or geometry """
+        """Common function to read either connect or geometry"""
         path = os.path.join(os.path.dirname(self.xdmfFilename), "")
-        dataLocation, data_prec, nElements, MemDimension = self.GetDataLocationPrecisionNElementsMemDimension(attribute)
+        (
+            dataLocation,
+            data_prec,
+            nElements,
+            MemDimension,
+        ) = self.GetDataLocationPrecisionNElementsMemDimension(attribute)
         # 3 for surface, 4 for volume
         dim2 = MemDimension[1]
         splitArgs = dataLocation.split(":")
@@ -143,21 +163,23 @@ def ReadTopologyOrGeometry(self, attribute):
             filename, hdf5var = splitArgs
             myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, 0, dim2)
         else:
-            myData = self.ReadSimpleBinaryFile(path + dataLocation, dim2, data_prec, isInt)
+            myData = self.ReadSimpleBinaryFile(
+                path + dataLocation, dim2, data_prec, isInt
+            )
             # due to zero padding in SeisSol for memory alignement, the array read may be larger than the actual data
             myData = myData[0:nElements, :]
         return myData
 
     def ReadConnect(self):
-        """ Read the connectivity matrice defining the cells """
+        """Read the connectivity matrice defining the cells"""
         return self.ReadTopologyOrGeometry("Topology")
 
     def ReadGeometry(self):
-        """ Read the connectivity matrice defining the cells """
+        """Read the connectivity matrice defining the cells"""
         return self.ReadTopologyOrGeometry("Geometry")
 
     def ReadNdt(self):
-        """ read number of time steps in the file """
+        """read number of time steps in the file"""
         root = self.tree.getroot()
         ndt = 0
         for Property in root.findall(".//Grid"):
@@ -189,33 +211,36 @@ def ReadAttributeValue(self, attribute, list_possible_location):
         raise NameError(f"{attribute} not found in {list_possible_location}")
 
     def ReadNElements(self):
-        """ read number of cell elements of the mesh """
+        """read number of cell elements of the mesh"""
         list_possible_location = ["Domain/Grid/Grid/Topology", "Domain/Grid/Topology"]
-        value  = self.ReadAttributeValue("NumberOfElements", list_possible_location)
+        value = self.ReadAttributeValue("NumberOfElements", list_possible_location)
         return int(value)
 
     def ReadNNodes(self):
-        """ read number of vertex of the mesh """
+        """read number of vertex of the mesh"""
         list_possible_location = ["Domain/Grid/Grid/Geometry", "Domain/Grid/Geometry"]
-        value  = self.ReadAttributeValue("NumberOfElements", list_possible_location)
+        value = self.ReadAttributeValue("NumberOfElements", list_possible_location)
         return int(value)
 
     def ReadNodesPerElement(self):
-        """ read number of nodes per elements of the mesh """
-        list_possible_location = ["Domain/Grid/Grid/Topology/DataItem", "Domain/Grid/Topology/DataItem"]
-        value  = self.ReadAttributeValue("Dimensions", list_possible_location)
+        """read number of nodes per elements of the mesh"""
+        list_possible_location = [
+            "Domain/Grid/Grid/Topology/DataItem",
+            "Domain/Grid/Topology/DataItem",
+        ]
+        value = self.ReadAttributeValue("Dimensions", list_possible_location)
         return int(value.split()[1])
 
     def ReadAvailableDataFields(self):
-        """ read all available data fields, e.g. SRs or P_n """
+        """read all available data fields, e.g. SRs or P_n"""
         root = self.tree.getroot()
         availableDataFields = set()
         for Property in root.findall(".//Attribute"):
             availableDataFields.add(Property.get("Name"))
         return availableDataFields
 
     def ReadTimeStep(self):
-        """ reading the time step (dt) in the xdmf file """
+        """reading the time step (dt) in the xdmf file"""
         root = self.tree.getroot()
         i = 0
         for Property in root.findall("Domain/Grid/Grid/Time"):
@@ -228,52 +253,84 @@ def ReadTimeStep(self):
         raise NameError("time step could not be determined")
 
     def Read1dData(self, dataName, nElements, isInt=False):
-        """ Read 1 dimension array (used by ReadPartition) """
+        """Read 1 dimension array (used by ReadPartition)"""
         path = os.path.join(os.path.dirname(self.xdmfFilename), "")
-        dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName)
+        (
+            dataLocation,
+            data_prec,
+            MemDimension,
+        ) = self.GetDataLocationPrecisionMemDimension(dataName)
         splitArgs = dataLocation.split(":")
         isHdf5 = True if len(splitArgs) == 2 else False
         if isHdf5:
             filename, hdf5var = splitArgs
-            myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, 0, MemDimension)
+            myData = self.ReadHdf5DatasetChunk(
+                path + filename, hdf5var, 0, MemDimension
+            )
         else:
             filename = dataLocation
-            myData = self.ReadSimpleBinaryFile(path + dataLocation, nElements, data_prec, isInt, 0)
+            # The Dimensions of the 1D arrays (e.g. partition) in the xdmf files
+            # are not not accounting for potential 0 padding (bug in xdmfwriter?).
+            # Therefore, we precise that there is only one "time-step" to read at location 0
+            # by setting idt=0
+            myData = self.ReadSimpleBinaryFile(
+                path + dataLocation, nElements, data_prec, isInt, 0
+            )
         return myData
 
     def ReadPartition(self):
-        """ Read partition array """
+        """Read partition array"""
         partition = self.Read1dData("partition", self.nElements, isInt=True)
         return partition
 
     def ReadData(self, dataName, idt=-1):
-        """ Load a data array named 'dataName' (e.g. SRs)
+        """Load a data array named 'dataName' (e.g. SRs)
         if idt!=-1, only the time step idt is loaded
-        else all time steps are loaded   """
+        else all time steps are loaded"""
 
-        return self.ReadDataChunk(dataName, firstElement=0, nchunk=self.nElements, idt=idt)
+        return self.ReadDataChunk(
+            dataName, firstElement=0, nchunk=self.nElements, idt=idt
+        )
 
     def ReadDataChunk(self, dataName, firstElement, nchunk, idt=-1):
-        """ Load a chunk of a data array named 'dataName' (e.g. SRs)
+        """Load a chunk of a data array named 'dataName' (e.g. SRs)
         That is instead of loading 0:nElements, load firstElement:firstElement+nchunk
         This function is used for generating in parallel Ground motion estimate maps
         if idt!=-1, only the time step idt is loaded
-        else all time steps are loaded   """
+        else all time steps are loaded"""
         path = os.path.join(os.path.dirname(self.xdmfFilename), "")
-        dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName)
+        (
+            dataLocation,
+            data_prec,
+            MemDimension,
+        ) = self.GetDataLocationPrecisionMemDimension(dataName)
         splitArgs = dataLocation.split(":")
         isHdf5 = True if len(splitArgs) == 2 else False
         oneDtMem = True if idt != -1 else False
         if isHdf5:
             filename, hdf5var = splitArgs
-            myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, firstElement, nchunk, idt)
+            myData = self.ReadHdf5DatasetChunk(
+                path + filename, hdf5var, firstElement, nchunk, idt
+            )
         else:
-            myData = self.ReadSimpleBinaryFileChunk(path + dataLocation, MemDimension, data_prec, isInt=False, firstElement=firstElement, nchunk=nchunk, idt=idt)
+            myData = self.ReadSimpleBinaryFileChunk(
+                path + dataLocation,
+                MemDimension,
+                data_prec,
+                isInt=False,
+                firstElement=firstElement,
+                nchunk=nchunk,
+                idt=idt,
+            )
         return myData
 
     def LoadData(self, dataName, nElements, idt=0, oneDtMem=False, firstElement=-1):
-        """ Do the same as ReadDataChunk. here for backward compatibility """
-        dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName)
+        """Do the same as ReadDataChunk. here for backward compatibility"""
+        (
+            dataLocation,
+            data_prec,
+            MemDimension,
+        ) = self.GetDataLocationPrecisionMemDimension(dataName)
         if not oneDtMem:
             idt = -1
         if firstElement == -1: