Closed s-lunowa closed 2 years ago
jgostick
commented
2 years ago What do you mean by 'unreadable'? XDMF files open fine for me. I forget the nuances of this io function, but it seems to me that the reversed dimensions were done deliberately.
s-lunowa
commented
2 years ago In the following, I provide a minimal working example to show the issue. The code creates a network consisting of 2 pores and one connecting throat and saves that into a XDMF file:
import openmnm
network = openpnm.network.Cubic(shape=[2, 1, 1])
openpnm.io.XDMF.export_data(network=network, filename='wrong')ParaView version 5.7.0 does open the file. However, the data for the pore coord(inate)s and throat conn(ection)s is missing (only zeros), see the figure:
The same happens for any saved data that has more than one dimension due to the reversed order (see my above comment).
When using the following code instead, this is fixed:
import openpnm
# imports for XDMF
from flatdict import FlatDict
import xml.etree.cElementTree as ET
from openpnm.io import Dict, GenericIO
from openpnm.utils import logging
logger = logging.getLogger(__name__)
# addition functions defined in openpnm.io.XDMF
from openpnm.io.XDMF import create_root, create_domain, create_geometry, create_topology, create_attribute, create_time, create_grid
def create_data_item(value, Dimensions, **attribs):
r"""Correction of misspelled `DataType` in openpnm.io.XDMF.create_data_item"""
element = ET.Element('DataItem')
element.attrib.update({'ItemType': "Uniform",
'Format': "XML",
'DataType': "Float",
'Precision': "4",
'Rank': "1",
'Dimensions': Dimensions,
'Function': None})
element.attrib.update(attribs)
if element.attrib['Function'] is None:
del element.attrib['Function']
element.text = value
return element
class XDMF(openpnm.io.XDMF):
r"""Corretion of the export_data method."""
@classmethod
def export_data(cls, network, phases=[], filename=''):
r"""
Saves (transient/steady-state) data from the given objects into the
specified file.
Parameters
----------
network : OpenPNM Network Object
The network containing the desired data
phases : list of OpenPNM Phase Objects (optional, default is none)
A list of phase objects whose data are to be included
Notes
-----
This method only saves the data, not any of the pore-scale models or
other attributes.
"""
import h5py
project, network, phases = cls._parse_args(network=network,
phases=phases)
network = network[0]
# Check if any of the phases has time series
transient = GenericIO._is_transient(phases=phases)
if filename == '':
filename = project.name
path = cls._parse_filename(filename=filename, ext='xmf')
# Path is a pathlib object, so slice it up as needed
fname_xdf = path.name
d = Dict.to_dict(network, phases=phases, interleave=True,
flatten=False, categorize_by=['element', 'data'])
D = FlatDict(d, delimiter='/')
# Identify time steps
t_steps = []
if transient:
for key in D.keys():
if '@' in key:
t_steps.append(key.split('@')[1])
t_steps = list(set(t_steps))
t_grid = create_grid(Name="TimeSeries", GridType="Collection",
CollectionType="Temporal")
# If steady-state, define '0' time step
if not transient:
t_steps.append('0')
# Setup xdmf file
root = create_root('Xdmf')
domain = create_domain()
# Iterate over time steps present
for i, t_step in enumerate(t_steps):
# Define the hdf file
if not transient:
fname_hdf = path.stem+".hdf"
else:
fname_hdf = path.stem+'@'+t_step+".hdf"
path_p = path.parent
f = h5py.File(path_p.joinpath(fname_hdf), "w")
# Add coordinate and connection information to top of HDF5 file
f["coordinates"] = network["pore.coords"]
f["connections"] = network["throat.conns"]
# geometry coordinates
row, col = f["coordinates"].shape
dims = ' '.join((str(row), str(col)))
hdf_loc = fname_hdf + ":coordinates"
geo_data = create_data_item(value=hdf_loc, Dimensions=dims,
Format='HDF', DataType="Float")
geo = create_geometry(GeometryType="XYZ")
geo.append(geo_data)
# topolgy connections
row, col = f["connections"].shape
dims = ' '.join((str(row), str(col)))
hdf_loc = fname_hdf + ":connections"
topo_data = create_data_item(value=hdf_loc, Dimensions=dims,
Format="HDF", DataType="Int")
topo = create_topology(TopologyType="Polyline",
NodesPerElement=str(2),
NumberOfElements=str(row))
topo.append(topo_data)
# Make HDF5 file with all datasets, and no groups
for item in D.keys():
if D[item].dtype == 'O':
logger.warning(item + ' has dtype object,'
+ ' will not write to file')
del D[item]
elif 'U' in str(D[item][0].dtype):
pass
elif ('@' in item and t_step == item.split('@')[1]):
f.create_dataset(name='/'+item.split('@')[0]+'@t',
shape=D[item].shape,
dtype=D[item].dtype,
data=D[item],
compression="gzip")
elif ('@' not in item and i == 0):
f.create_dataset(name='/'+item, shape=D[item].shape,
dtype=D[item].dtype, data=D[item],
compression="gzip")
# Create a grid
grid = create_grid(Name=t_step, GridType="Uniform")
time = create_time(mode='Single', Value=t_step)
grid.append(time)
# Add pore and throat properties
for item in D.keys():
if item not in ['coordinates', 'connections']:
if ("@" in item and t_step == item.split("@")[1]) or (
"@" not in item
):
attr_type = 'Scalar'
shape = D[item].shape
dims = (' '.join([str(i) for i in shape]))
if '@' in item:
item = item.split('@')[0]+'@t'
hdf_loc = fname_hdf + ":" + item
elif ('@' not in item and i == 0):
hdf_loc = fname_hdf + ":" + item
elif ('@' not in item and i > 0):
hdf_loc = path.stem + '@' + t_steps[0] + ".hdf" + ":" + item
attr = create_data_item(value=hdf_loc,
Dimensions=dims,
Format='HDF',
Precision='8',
DataType='Float')
name = item.replace('/', ' | ')
if 'throat' in item:
Center = "Cell"
else:
Center = "Node"
el_attr = create_attribute(Name=name, Center=Center,
AttributeType=attr_type)
el_attr.append(attr)
grid.append(el_attr)
else:
pass
grid.append(topo)
grid.append(geo)
t_grid.append(grid)
# Close the HDF5 file
f.close()
domain.append(t_grid)
root.append(domain)
with open(path_p.joinpath(fname_xdf), 'w') as file:
file.write(cls._header)
file.write(ET.tostring(root).decode("utf-8"))
# save network with 2 pores and 1 throat
network = openpnm.network.Cubic(shape=[2, 1, 1])
XDMF.export_data(network=network, filename='correct')Note that the produced HDF5 files are the same, only the XMF files differ. Opening this in ParaView shows the following:
I do not know why the dimensions were reversed in the first place, but it seems strange since the HDF5 data has the original dimensions, as can be seen using e.g.
$ h5dump -d /net_01/properties/pore/coords wrong.hdf
HDF5 "correct.hdf" {
DATASET "/net_01/properties/pore/coords" {
DATATYPE H5T_IEEE_F64LE
DATASPACE SIMPLE { ( 2, 3 ) / ( 2, 3 ) }
DATA {
(0,0): 0.5, 0.5, 0.5,
(1,0): 1.5, 0.5, 0.5
}
}
}
Hey @Ni2M , can you take a look into this?
Ni2M
commented
2 years ago Hey @Ni2M , can you take a look into this?
Sure, I'll investigate the issue. Thanks.
ma-sadeghi
commented
2 years ago @s-lunowa Could your pull the latest dev branch and see if this issue has been fixed? Thanks!
s-lunowa
commented
2 years ago @ma-sadeghi Thank you, this bug (and the typo) are fixed by #2270
The XDMF output creates HDF5 files with dimensions given by
ndarray.shape: https://github.com/PMEAL/OpenPNM/blob/f176ea2cd76f7879edc7eabe13115a6ba9a8b18e/openpnm/io/XDMF.py#L109-L118 However, in the XDMF file, the dimensions are reversed: https://github.com/PMEAL/OpenPNM/blob/f176ea2cd76f7879edc7eabe13115a6ba9a8b18e/openpnm/io/XDMF.py#L130-L131 https://github.com/PMEAL/OpenPNM/blob/f176ea2cd76f7879edc7eabe13115a6ba9a8b18e/openpnm/io/XDMF.py#L139-L143 This results in unreadable data for ParaView, which replaces it by zeros. This is fixed, when replacing Line 131 byFurthermore, there is a small typo: https://github.com/PMEAL/OpenPNM/blob/f176ea2cd76f7879edc7eabe13115a6ba9a8b18e/openpnm/io/XDMF.py#L223 This should be: