seissolxdm/seissolxdmfwriter: bugfixes

seissolxdmf/seissolxdmf/__init__.py

@@ -1,2 +1,2 @@
 from .seissolxdmf import *
-__version__ = '0.1.2'
+__version__ = '0.1.3'

seissolxdmf/seissolxdmf/seissolxdmf.py

@@ -238,7 +238,7 @@ def Read1dData(self, dataName, nElements, isInt=False):
             myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, 0, MemDimension)
         else:
             filename = dataLocation
-            myData = self.ReadSimpleBinaryFile(path + dataLocation, nElements, data_prec, isInt)
+            myData = self.ReadSimpleBinaryFile(path + dataLocation, nElements, data_prec, isInt, 0)
         return myData
 
     def ReadPartition(self):

seissolxdmfwriter/README.md

@@ -36,7 +36,7 @@ sxw.write(
 
 sxw.write_from_seissol_output(
     'test-fault-sx',
-    fn,
+    sx,
     ['SRs', 'SRd','fault-tag', 'partition'],
     [3,4],
     reduce_precision=True,

seissolxdmfwriter/seissolxdmfwriter/__init__.py

@@ -1,2 +1,2 @@
 from .seissolxdmfwriter import *
-__version__ = '0.3.0'
+__version__ = '0.4.0'

seissolxdmfwriter/seissolxdmfwriter/seissol_output_extractor.py

@@ -5,6 +5,7 @@
 import seissolxdmfwriter as sxw
 import numpy as np
 import argparse
+from warnings import warn
 
 
 def generate_new_prefix(prefix, append2prefix):
@@ -29,7 +30,7 @@ def generate_new_prefix(prefix, append2prefix):
     "--add2prefix",
     help="string to append to the prefix in the new file",
     type=str,
-    default="_resampled",
+    default="_extracted",
 )
 parser.add_argument(
     "--variables",
@@ -131,6 +132,12 @@ def ReadData(self, dataName, idt=-1):
         else:
             return super().ReadData(dataName, idt)
 
+    def GetDataLocationPrecisionMemDimension(self, dataName):
+        if dataName == "SR" and "SR" not in self.ReadAvailableDataFields():
+            return super().GetDataLocationPrecisionMemDimension("SRs")
+        else:
+            return super().GetDataLocationPrecisionMemDimension(dataName)
+
 
 def main():
     sx = SeissolxdmfExtended(args.xdmfFilename)
@@ -139,7 +146,7 @@ def main():
     if spatial_filtering:
         xyz = sx.ReadGeometry()
         connect = sx.ReadConnect()
-        print("Warning: spatial filtering significantly slows down this script")
+        warn("spatial filtering significantly slows down this script")
         ids = range(0, sx.nElements)
         xyzc = (
             xyz[connect[:, 0], :] + xyz[connect[:, 1], :] + xyz[connect[:, 2], :]
@@ -182,9 +189,16 @@ def filter_cells(coords, filter_range):
         args.variables = sorted(sx.ReadAvailableDataFields())
         print(f"args.variables was set to all and now contains {args.variables}")
 
+    if args.backend == "hdf5" and args.compression > 0:
+        print(
+            "Writing hdf5 output with compression enabled"
+            f" (compression_level={args.compression}). \n"
+            "Use --compression=0 if you want to speed-up data extraction."
+        )
+
     sxw.write_from_seissol_output(
         prefix_new,
-        args.xdmfFilename,
+        sx,
         args.variables,
         indices,
         reduce_precision=True,

seissolxdmfwriter/seissolxdmfwriter/seissolxdmfwriter.py

@@ -42,6 +42,7 @@ def write_timeseries_xdmf(
     reduce_precision,
     backend,
 ):
+    bn_prefix = os.path.basename(prefix)
     data_format = "HDF" if backend == "hdf5" else "Binary"
     lDataName = list(dictDataTypes.keys())
     lDataTypes = list(dictDataTypes.values())
@@ -50,8 +51,8 @@ def write_timeseries_xdmf(
 <!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
 <Xdmf Version="2.0">
  <Domain>"""
-    geometry_location = dataLocation(prefix, "geometry", backend)
-    connect_location = dataLocation(prefix, "connect", backend)
+    geometry_location = dataLocation(bn_prefix, "geometry", backend)
+    connect_location = dataLocation(bn_prefix, "connect", backend)
 
     for i, ctime in enumerate(timeValues):
         xdmf += f"""
@@ -65,7 +66,7 @@ def write_timeseries_xdmf(
     </Geometry>
     <Time Value="{ctime}"/>"""
         for k, dataName in enumerate(list(dictDataTypes.keys())):
-            data_location = dataLocation(prefix, dataName, backend)
+            data_location = dataLocation(bn_prefix, dataName, backend)
             prec, number_type = dictDataTypes[dataName]
             if dataName in known_1d_arrays:
                 xdmf += f"""
@@ -90,6 +91,8 @@ def write_timeseries_xdmf(
     with open(prefix + ".xdmf", "w") as fid:
         fid.write(xdmf)
     print(f"done writing {prefix}.xdmf")
+    full_path = os.path.abspath(f"{prefix}.xdmf")
+    print(f"full path: {full_path}")
 
 
 def write_mesh_xdmf(
@@ -101,6 +104,7 @@ def write_mesh_xdmf(
     reduce_precision,
     backend,
 ):
+    bn_prefix = os.path.basename(prefix)
     data_format = "HDF" if backend == "hdf5" else "Binary"
     lDataName = list(dictDataTypes.keys())
     lDataTypes = list(dictDataTypes.values())
@@ -109,8 +113,8 @@ def write_mesh_xdmf(
 <!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
 <Xdmf Version="2.0">
  <Domain>"""
-    geometry_location = dataLocation(prefix, "geometry", backend)
-    connect_location = dataLocation(prefix, "connect", backend)
+    geometry_location = dataLocation(bn_prefix, "geometry", backend)
+    connect_location = dataLocation(bn_prefix, "connect", backend)
 
     xdmf += f"""
   <Grid Name="puml mesh" GridType="Uniform">
@@ -121,7 +125,7 @@ def write_mesh_xdmf(
      <DataItem NumberType="Float" Precision="8" Format="{data_format}" Dimensions="{nNodes} 3">{geometry_location}</DataItem>
     </Geometry>"""
     for k, dataName in enumerate(list(dictDataTypes.keys())):
-        data_location = dataLocation(prefix, dataName, backend)
+        data_location = dataLocation(bn_prefix, dataName, backend)
         prec, number_type = dictDataTypes[dataName]
         xdmf += f"""
     <Attribute Name="{dataName}" Center="Cell">
@@ -135,6 +139,8 @@ def write_mesh_xdmf(
     with open(prefix + ".xdmf", "w") as fid:
         fid.write(xdmf)
     print(f"done writing {prefix}.xdmf")
+    full_path = os.path.abspath(f"{prefix}.xdmf")
+    print(f"full path: {full_path}")
 
 
 def output_type(input_array, reduce_precision):
@@ -191,7 +197,7 @@ def read_non_temporal(sx, ar_name, filtered_cells):
         elif ar_name == "connect":
             return sx.ReadConnect()[filtered_cells, :]
         else:
-            return sx.Read1dData(ar_name, sx.nElements, isInt=True)[:, filtered_cells]
+            return sx.Read1dData(ar_name, sx.nElements, isInt=True)[filtered_cells]
 
     nel = infer_n_elements(sx, filtered_cells)
     if backend == "hdf5":
@@ -215,14 +221,14 @@ def read_non_temporal(sx, ar_name, filtered_cells):
                     desc=ar_name,
                     dynamic_ncols=False,
                 ):
+                    my_array = sx.ReadData(ar_name, idt)[filtered_cells]
                     if i == 0:
                         h5f.create_dataset(
                             f"/{ar_name}",
                             (len(time_indices), nel),
                             dtype=str(output_type(my_array, reduce_precision)),
                             **compression_options,
                         )
-                    my_array = sx.ReadData(ar_name, idt)[filtered_cells]
                     if my_array.shape[0] == 0:
                         print(
                             f"time step {idt} of {ar_name} is corrupted, replacing with 0s"
@@ -438,7 +444,7 @@ def write(
 
 def write_from_seissol_output(
     prefix,
-    input_file,
+    sx,
     var_names,
     time_indices,
     reduce_precision=False,
@@ -449,21 +455,17 @@ def write_from_seissol_output(
     """
     Write hdf5/xdmf files output, readable by ParaView from a seissolxdmf object
     prefix: file
-    input_file: filename of the seissol input
+    sx: seissolxdmf object
     var_names: list of variables to extract
     time_indices: list of times indices to extract
     reduce_precision: convert double to float and i64 to i32 if True
     backend: data format ("hdf5" or "raw")
     """
-    import seissolxdmf as sx
-
     if backend not in ("hdf5", "raw"):
         raise ValueError(f"Invalid backend {backend}. Must be 'hdf5' or 'raw'.")
     if compression_level < 0 or compression_level > 9:
         raise ValueError("compression_level has to be in 0-9")
 
-    sx = sx.seissolxdmf(input_file)
-
     non_temporal_array_names = ["geometry", "connect"]
     to_move = [name for name in var_names if name in known_1d_arrays]
     dictDataTypes = {}

seissolxdmfwriter/setup.py

@@ -34,5 +34,5 @@ def get_property(prop, project):
         "Operating System :: OS Independent",
     ],
     python_requires=">=3.6",
-    install_requires=["numpy", "h5py","seissolxdmf>=0.1.2"],
+    install_requires=["numpy", "h5py","seissolxdmf>=0.1.2", "tqdm"],
 )