Repository: PMEAL/OpenPNM
Pull request #1242 (2019-04-16 15:00:59)

Support for transient data
Export transient data in XDMF format.
1- Supports transient and steady-state data
2- Steady-state data (such as geometrical properties) are not duplicated over transient files.
3- On Paraview, the Steady-state data are present at all time steps by pointing the "xmf" file to the initial time step file "hdf" which contains all the steady-state data.
4- Removed the option to export VTK files with time series as it doesn't properly handle non-uniform time steps AND generate large files...
5- When times series are present in phases, proj.export_data method will ALWAYS export in XDMF format an print a warning if the user requested to export into another format.


Patch:
diff --git a/example_transient_Fickian_diffusion.py b/example_transient_Fickian_diffusion.py
new file mode 100644
index 000000000..ff8435413
--- /dev/null
+++ b/example_transient_Fickian_diffusion.py
@@ -0,0 +1,38 @@
+import openpnm as op
+import numpy as np
+
+# work space and project
+ws = op.Workspace()
+ws.settings["loglevel"] = 30
+proj = ws.new_project()
+
+# network
+np.random.seed(7)
+net = op.network.Cubic(shape=[51, 19, 1], spacing=1e-4, project=proj)
+
+# geometry
+geo = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
+
+# phase
+phase = op.phases.Water(network=net)
+
+# physics
+phys = op.physics.GenericPhysics(network=net, phase=phase, geometry=geo)
+phase['pore.diffusivity'] = 2e-09
+phase['throat.diffusivity'] = 2e-09
+
+mod = op.models.physics.diffusive_conductance.ordinary_diffusion
+phys.add_model(propname='throat.diffusive_conductance',
+               model=mod, regen_mode='normal')
+
+# algorithms: Fickian diffusion
+fd = op.algorithms.TransientFickianDiffusion(network=net, phase=phase)
+fd.set_value_BC(pores=net.pores('front'), values=0.5)
+fd.set_value_BC(pores=net.pores('back'), values=0.1)
+fd.setup(t_final=100, t_output=10, t_step=0.5)
+fd.run()
+phase.update(fd.results())
+
+# output results to a vtk file
+phase.update(fd.results())
+proj.export_data(phases=[phase], filename='OUT', filetype='vtk')
diff --git a/openpnm/io/GenericIO.py b/openpnm/io/GenericIO.py
index 882a1ff1b..3ba0adcf7 100644
--- a/openpnm/io/GenericIO.py
+++ b/openpnm/io/GenericIO.py
@@ -90,3 +90,16 @@ class GenericIO():
         else:
             project = network[0].project
         return (project, network, phases)
+
+    @classmethod
+    def is_transient(cls, phases):
+        # Check if any of the phases has time series
+        transient = False
+        if type(phases) == str:
+            transient = True in ['@' in k for k in phases.keys()]
+        elif type(phases) == list:
+            for phase in phases:
+                transient = True in ['@' in k for k in phase.keys()]
+                if transient:
+                    break
+        return (transient)
diff --git a/openpnm/io/VTK.py b/openpnm/io/VTK.py
index ef6910d1e..1c4ec5aa8 100644
--- a/openpnm/io/VTK.py
+++ b/openpnm/io/VTK.py
@@ -3,7 +3,6 @@ from flatdict import FlatDict
 import numpy as np
 from openpnm.io import GenericIO, Dict
 from openpnm.utils import logging, Workspace
-from openpnm.network import GenericNetwork
 logger = logging.getLogger(__name__)
 ws = Workspace()
 
@@ -12,7 +11,6 @@ class VTK(GenericIO):
     r"""
     The Visualization Toolkit (VTK) format defined by Kitware and used by
     Paraview
-
     Because OpenPNM data is unstructured, the actual output format is VTP,
     not VTK.
     """
@@ -40,32 +38,26 @@ class VTK(GenericIO):
         r"""
         Save network and phase data to a single vtp file for visualizing in
         Paraview
-
         Parameters
         ----------
         network : OpenPNM Network Object
             The Network containing the data to be written
-
         phases : list, optional
             A list containing OpenPNM Phase object(s) containing data to be
             written
-
         filename : string, optional
             Filename to write data.  If no name is given the file is named
             after the network
-
         delim : string
             Specify which character is used to delimit the data names.  The
             default is ' | ' which creates a nice clean output in the Paraview
             pipeline viewer (e.g. net | property | pore | diameter)
-
         fill_nans : scalar
             The value to use to replace NaNs with.  The VTK file format does
             not work with NaNs, so they must be dealt with.  The default is
             `None` which means property arrays with NaNs are not written to the
             file.  Other useful options might be 0 or -1, but the user must
             be aware that these are not real values, only place holders.
-
         """
         project, network, phases = cls._parse_args(network=network,
                                                    phases=phases)
@@ -115,14 +107,6 @@ class VTK(GenericIO):
                         continue
                     else:
                         array[np.isnan(array)] = fill_nans
-                if '@' in key:
-                    new_net = GenericNetwork(coords=network['pore.coords'],
-                                             conns=network['throat.conns'])
-                    prop = key.split(delim)[-1]
-                    new_net[prop.split('@')[0]] = array
-                    file, ext = filename.name.split('.')
-                    fname_temp = file + '_' + key.split('@')[1] + '.' + ext
-                    cls.save(network=new_net, filename=fname_temp, delim=delim)
                 element = VTK._array_to_element(key, array)
                 if (array.size == num_points):
                     point_data_node.append(element)
@@ -143,17 +127,14 @@ class VTK(GenericIO):
     def load(cls, filename, project=None, delim=' | '):
         r"""
         Read in pore and throat data from a saved VTK file.
-
         Parameters
         ----------
         filename : string (optional)
             The name of the file containing the data to import.  The formatting
             of this file is outlined below.
-
         project : OpenPNM Project object
             A GenericNetwork is created and added to the specified Project.
             If no Project is supplied then one will be created and returned.
-
         """
         net = {}
 
diff --git a/openpnm/io/XDMF.py b/openpnm/io/XDMF.py
index 860bf5850..1b87afdd7 100644
--- a/openpnm/io/XDMF.py
+++ b/openpnm/io/XDMF.py
@@ -22,7 +22,8 @@ class XDMF(GenericIO):
     @classmethod
     def save(cls, network, phases=[], filename=''):
         r"""
-        Saves data from the given objects into the specified file.
+        Saves (transient/steady-state) data from the given objects into the
+        specified file.
 
         Parameters
         ----------
@@ -42,94 +43,117 @@ class XDMF(GenericIO):
         project, network, phases = cls._parse_args(network=network,
                                                    phases=phases)
         network = network[0]
-
         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
-        fname_hdf = path.stem+".hdf"
-        path = path.parent
-        f = h5py.File(path.joinpath(fname_hdf), "w")
-
         d = Dict.to_dict(network, phases=phases, interleave=True,
                          flatten=False, categorize_by=['element', 'data'])
-
-        # Make HDF5 file with all datasets, and no groups
         D = FlatDict(d, delimiter='/')
-        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
-            else:
-                f.create_dataset(name='/'+item, shape=D[item].shape,
-                                 dtype=D[item].dtype, data=D[item])
-        # Add coordinate and connection information to top of HDF5 file
-        f["coordinates"] = network["pore.coords"]
-        f["connections"] = network["throat.conns"]
-
-        # setup xdmf file
+        # Identify time steps
+        t_steps = []
+        for key in D.keys():
+            if '@' in key:
+                t_steps.append(key.split('@')[1])
+        t_grid = create_grid(Name="TimeSeries", GridType="Collection",
+                             CollectionType="Temporal")
+        # If steady-state, define '0' time step
+        if t_steps == []:
+            t_steps.append('0')
+        # Setup xdmf file
         root = create_root('Xdmf')
         domain = create_domain()
-        grid = create_grid(Name="Structure", GridType="Uniform")
-
-        # 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  # col first then row
-        dims = ' '.join((str(row), str(col)))
-        hdf_loc = fname_hdf + ":connections"
-        topo_data = create_data_item(value=hdf_loc, Dimensions=dims,
-                                     Format="HDF", NumberType="Int")
-        topo = create_topology(TopologyType="Polyline",
-                               NodesPerElement=str(2),
-                               NumberOfElements=str(row))
-        topo.append(topo_data)
-
-        # Add pore and throat properties
-        for item in D.keys():
-            if item not in ['coordinates', 'connections']:
-                attr_type = 'Scalar'
-                shape = f[item].shape
-                dims = ''.join([str(i) + ' ' for i in list(shape)[::-1]])
-                hdf_loc = fname_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)
-
-        grid.append(topo)
-        grid.append(geo)
-        domain.append(grid)
+        # Iterate over time steps present
+        for t in range(len(t_steps)):
+            # Define the hdf file
+            fname_hdf = path.stem+'_'+t_steps[t]+".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  # col first then row
+            dims = ' '.join((str(row), str(col)))
+            hdf_loc = fname_hdf + ":connections"
+            topo_data = create_data_item(value=hdf_loc, Dimensions=dims,
+                                         Format="HDF", NumberType="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 ('@' and '@'+t_steps[t] in item):
+                    f.create_dataset(name='/'+item.split('@')[0]+'_t',
+                                     shape=D[item].shape,
+                                     dtype=D[item].dtype,
+                                     data=D[item])
+                elif ('@' not in item and t == 0):
+                    f.create_dataset(name='/'+item, shape=D[item].shape,
+                                     dtype=D[item].dtype, data=D[item])
+            # Create a grid
+            grid = create_grid(Name=t_steps[t], GridType="Uniform")
+            time = create_time(type='Single', Value=t_steps[t])
+            grid.append(time)
+            # Add pore and throat properties
+            for item in D.keys():
+                if item not in ['coordinates', 'connections']:
+                    if ('@' and '@'+t_steps[t] in item) or ('@' not in item):
+                        attr_type = 'Scalar'
+                        shape = D[item].shape
+                        dims = (''.join([str(i) +
+                                         ' ' for i in list(shape)[::-1]]))
+                        if '@' in item:
+                            item = item.split('@')[0]+'_t'
+                            hdf_loc = fname_hdf + ":" + item
+                        elif ('@' not in item and t == 0):
+                            hdf_loc = fname_hdf + ":" + item
+                        elif ('@' not in item and t > 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.joinpath(fname_xdf), 'w') as file:
+        with open(path_p.joinpath(fname_xdf), 'w') as file:
             file.write(cls._header)
             file.write(ET.tostring(root).decode("utf-8"))
 
-        # CLose the HDF5 file
-        f.close()
-
 
 def create_root(Name):
     return ET.Element(Name)
@@ -165,7 +189,7 @@ def create_attribute(Name, **attribs):
 
 
 def create_time(type='Single', Value=None):
-    element = ET.Element('Attribute')
+    element = ET.Element('Time')
     if type == 'Single' and Value:
         element.attrib['Value'] = Value
     return element
diff --git a/openpnm/utils/Project.py b/openpnm/utils/Project.py
index 29a07336c..a7bf705af 100644
--- a/openpnm/utils/Project.py
+++ b/openpnm/utils/Project.py
@@ -5,8 +5,10 @@ import numpy as np
 import openpnm
 from copy import deepcopy
 from pathlib import Path
-from openpnm.utils import SettingsDict, HealthDict, PrintableList, Workspace
+from openpnm.utils import (SettingsDict, HealthDict, PrintableList, Workspace,
+                           logging)
 ws = Workspace()
+logger = logging.getLogger(__name__)
 
 
 class Project(list):
@@ -579,28 +581,37 @@ class Project(list):
         network = self.network
         if filename is None:
             filename = project.name + '_' + time.strftime('%Y%b%d_%H%M%p')
-        # Infer filetype from extension on file name...if given
+        # Check if any of the phases has time series
+        transient = openpnm.io.GenericIO.is_transient(phases=phases)
+        # Infer filetype from extension on file name..., if given
         if '.' in filename:
             exts = ['vtk', 'vtp', 'vtu', 'csv', 'xmf', 'xdmf', 'hdf', 'hdf5',
                     'h5', 'mat']
             if filename.split('.')[-1] in exts:
                 filename, filetype = filename.rsplit('.', 1)
-        if filetype.lower() in ['vtk', 'vtp', 'vtu']:
-            openpnm.io.VTK.save(network=network, phases=phases,
-                                filename=filename)
-        if filetype.lower() == 'csv':
-            openpnm.io.CSV.save(network=network, phases=phases,
-                                filename=filename)
-        if filetype.lower() in ['xmf', 'xdmf']:
+        if transient:
+            if filetype.lower() not in ['xmf', 'xdmf']:
+                logger.warning(filetype.lower() + ' filetype not supported ' +
+                               'for transient data, xdmf is used instead')
             openpnm.io.XDMF.save(network=network, phases=phases,
                                  filename=filename)
-        if filetype.lower() in ['hdf5', 'hdf', 'h5']:
-            f = openpnm.io.HDF5.to_hdf5(network=network, phases=phases,
-                                        filename=filename)
-            f.close()
-        if filetype.lower() == 'mat':
-            openpnm.io.MAT.save(network=network, phases=phases,
-                                filename=filename)
+        else:
+            if filetype.lower() in ['vtk', 'vtp', 'vtu']:
+                openpnm.io.VTK.save(network=network, phases=phases,
+                                    filename=filename)
+            if filetype.lower() == 'csv':
+                openpnm.io.CSV.save(network=network, phases=phases,
+                                    filename=filename)
+            if filetype.lower() in ['xmf', 'xdmf']:
+                openpnm.io.XDMF.save(network=network, phases=phases,
+                                     filename=filename)
+            if filetype.lower() in ['hdf5', 'hdf', 'h5']:
+                f = openpnm.io.HDF5.to_hdf5(network=network, phases=phases,
+                                            filename=filename)
+                f.close()
+            if filetype.lower() == 'mat':
+                openpnm.io.MAT.save(network=network, phases=phases,
+                                    filename=filename)
 
     def _dump_data(self, mode=['props']):
         r"""


Test patch:
diff --git a/tests/unit/io/XDMFTest.py b/tests/unit/io/XDMFTest.py
index 0acffee52..7165722ec 100644
--- a/tests/unit/io/XDMFTest.py
+++ b/tests/unit/io/XDMFTest.py
@@ -65,7 +65,7 @@ class XDMFTest:
     def test_save(self, tmpdir):
         fname = tmpdir.join('test_file')
         op.io.XDMF.save(network=self.net, phases=self.phase_1, filename=fname)
-        os.remove(tmpdir.join('test_file.hdf'))
+        os.remove(tmpdir.join('test_file_0.hdf'))
         os.remove(tmpdir.join('test_file.xmf'))
 
 
diff --git a/tests/unit/utils/ProjectTest.py b/tests/unit/utils/ProjectTest.py
index 7090f7086..93c6fd86a 100644
--- a/tests/unit/utils/ProjectTest.py
+++ b/tests/unit/utils/ProjectTest.py
@@ -365,7 +365,7 @@ class ProjectTest:
         self.proj.export_data(phases=self.phase1, filename=fname,
                               filetype='xmf')
         os.remove(fname+'.xmf')
-        os.remove(fname+'.hdf')
+        os.remove(fname+'_0.hdf')
         self.proj.export_data(phases=self.phase1, filename=fname,
                               filetype='hdf')
         os.remove(fname+'.hdf')