diff --git a/CHANGELOG.md b/CHANGELOG.md index 38d3cfc..1877b5e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,173 +1,179 @@ # Changelog All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## Unreleased + +### Added + +- Added `read()` method to dumpers to create a model from a dump file + ## v2.2.2 -- 2021-04-02 ### Added - Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try executing `tamaas surface -h` from the command-line! ### Changed - `CXXFLAGS` are now passed to the linker - Added this changelog - Using absolute paths for environmental variables when running `scons test` - Reorganized documentation layout - Gave the build system a facelift (docs are now generated directly with SCons instead of a Makefile) ### Deprecated - Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with a Python 2 wheel. - The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command ### Fixed - `UVWDumper` no longer imports `mpi4py` in sequential - Compiling with different Thrust/FFTW backends ## v2.2.1 -- 2021-03-02 ### Added - Output registered fields and dumpers in `print(model)` - Added `operator[]` to the C++ model class (for fields) - Added `traction` and `displacement` properties to Python model bindings - Added `operators` property to Python model bindings, which provides a dict-like access to registered operators - Added `shape` and `spectrum` to properties to Python surface generator bindings - Surface generator constructor accepts surface global shape as argument - Choice of FFTW thread model ### Changed - Tests use `/tmp` for temporary files - Updated dependency versions (Thrust, Pybind11) ### Deprecated - Most `get___()` and `set___()` in Python bindings have been deprecated. They will generate a `DeprecationWarning`. ### Removed - All legacy code ## v2.2.0 -- 2020-12-31 ### Added - More accurate function for computation of contact area - Function to compute deviatoric of tensor fields - MPI implementation - Convenience `hdf5toVTK` function - Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give shape information ### Changed - Preprocessor defined macros are prefixed with `TAMAAS_` - Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt` ### Fixed - Warning about deprecated constructors with recent GCC versions - Wrong computation of grid strides - Wrong computation of grid sizes in views ## v2.1.4 -- 2020-08-07 ### Added - Possibility to generate a static `libTamaas` - C++ implementation of DFSANE solver - Allowing compilation without OpenMP ### Changed - NetCDF dumper writes frames to a single file ### Fixed - Compatibility with SCons+Python 3 ## v2.1.3 -- 2020-07-27 ### Added - Version number to `TamaasInfo` ### Changed - Prepending root directory when generating archive ## v2.1.2 -- 2020-07-24 This release changes some core internals related to discrete Fourier transforms for future MPI support. ### Added - Caching `CXXFLAGS` in SCons build - SCons shortcut to create code archive - Test of the elastic-plastic contact solver - Paraview data dumper (`.pvd` files) - Compression for UVW dumper - `__contains__` and `__iter__` Python bindings of model - Warning message of possible overflow in Kelvin ### Changed - Simplified `tamaas_info.cpp`, particularly the diff part - Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are re-used as much as possible with different data with the same shape. This is in view of future MPI developments - Redirecting I/O streams in solve functions so they can be used from Python (e.g. in Jupyter notebooks) - Calling `initialize()` and `finalize()` is no longer necessary ### Fixed - Convergence issue with non-linear solvers - Memory error in volume potentials ## v2.1.1 -- 2020-04-22 ### Added - SCons shortcut to run tests ### Fixed - Correct `RPATH` for shared libraries - Issues with SCons commands introduced in v2.1.0 - Tests with Python 2.7 ## v2.1.0 -- 2020-04-17 ### Added - SCons shortcuts to build/install Tamaas and its components - Selection of integration method for Kelvin operator - Compilation option to remove the legacy part of Tamaas - NetCDF dumper ### Fixed - Link bug with clang - NaNs in Kato saturated solver ## v2.0.0 -- 2019-11-11 First public release. Contains relatively mature elastic-plastic contact code. diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py index 6eda25e..b80f42f 100644 --- a/python/tamaas/dumpers/__init__.py +++ b/python/tamaas/dumpers/__init__.py @@ -1,375 +1,393 @@ # -*- mode:python; coding: utf-8 -*- # @file # LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 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 Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . """ Dumpers for the class tamaas.Model """ from __future__ import print_function from sys import stderr from os import makedirs import os.path import numpy as np from .. import ModelDumper, model_type, mpi, type_traits, ModelFactory from ._helper import step_dump, directory_dump, local_slice, _is_surface_field _reverse_trait_map = { 'model_type.' + t.__name__: mtype for mtype, t in type_traits.items() } def _get_attributes(model): "Get model attributes" return { 'model_type': str(model.type), 'system_size': model.system_size, 'discretization': model.global_shape, } def _create_model(attrs): mtype = _reverse_trait_map[attrs['model_type']] return ModelFactory.createModel(mtype, attrs['system_size'], attrs['discretization']) class FieldDumper(ModelDumper): """Abstract dumper for python classes using fields""" postfix = '' extension = '' name_format = "{basename}{postfix}.{extension}" def __init__(self, basename, *fields, **kwargs): """Construct with desired fields""" super(FieldDumper, self).__init__() self.basename = basename self.fields = list(fields) self.all_fields = kwargs.get('all_fields', False) def add_field(self, field): """Add another field to the dump""" if field not in self.fields: self.fields.append(field) def dump_to_file(self, file_descriptor, model): """Dump to a file (name or handle)""" def get_fields(self, model): """Get the desired fields""" if not self.all_fields: requested_fields = self.fields else: requested_fields = list(model) return {field: model[field] for field in requested_fields} def dump(self, model): "Dump model" self.dump_to_file(self.file_path, model) def read(self, file_descriptor): raise RuntimeError('read() method not implemented in {}' .format(type(self).__name__)) @property def file_path(self): """Get the default filename""" return self.name_format.format(basename=self.basename, postfix=self.postfix, extension=self.extension) @directory_dump('numpys') @step_dump class NumpyDumper(FieldDumper): """Dumper to compressed numpy files""" extension = 'npz' def dump_to_file(self, file_descriptor, model): """Saving to compressed multi-field Numpy format""" if mpi.size() > 1: raise RuntimeError("NumpyDumper does not function " "at all in parallel") np.savez_compressed(file_descriptor, attrs=_get_attributes(model), **self.get_fields(model)) def read(self, file_descriptor): + "Create model from Numpy file" data = np.load(file_descriptor, mmap_mode='r') model = _create_model(data['attrs'].item()) for k, v in filter(lambda k: k[0] != 'attrs', data.items()): model[k] = v return model try: import h5py @directory_dump('hdf5') @step_dump class H5Dumper(FieldDumper): """Dumper to HDF5 file format""" extension = 'h5' def dump_to_file(self, file_descriptor, model): """Saving to HDF5 with metadata about the model""" # Setup for MPI if not h5py.get_config().mpi and mpi.size() > 1: raise RuntimeError("HDF5 does not have MPI support") if mpi.size() > 1: from mpi4py import MPI # noqa mpi_args = dict(driver='mpio', comm=MPI.COMM_WORLD) comp_args = {} # compression does not work in parallel else: mpi_args = {} comp_args = dict(compression='gzip', compression_opts=7) with h5py.File(file_descriptor, 'w', **mpi_args) as handle: # Writing data for name, field in self.get_fields(model).items(): shape = list(field.shape) if mpi.size() > 1: xdim = 0 if _is_surface_field(field, model) else 1 shape[xdim] = MPI.COMM_WORLD.allreduce(shape[xdim]) dset = handle.create_dataset(name, shape, field.dtype, **comp_args) dset[local_slice(field, model)] = field # Writing metadata for name, attr in _get_attributes(model).items(): handle.attrs[name] = attr def read(self, file_descriptor): + "Create model from HDF5 file" if mpi.size() > 1: raise RuntimeError("Cannot read in MPI context yet") with h5py.File(file_descriptor, 'r') as handle: model = _create_model(handle.attrs) for k, v in handle.items(): model[k] = v[:] return model except ImportError: pass try: import uvw # noqa @directory_dump('paraview') @step_dump class UVWDumper(FieldDumper): """Dumper to VTK files for elasto-plastic calculations""" extension = 'vtr' forbidden_fields = ['traction', 'gap'] def dump_to_file(self, file_descriptor, model): """Dump displacements, plastic deformations and stresses""" if mpi.size() > 1: raise RuntimeError("UVWDumper does not function " "properly in parallel") bdim = len(model.boundary_shape) # Local MPI size lsize = model.shape gsize = mpi.global_shape(model.boundary_shape) gshape = gsize if len(lsize) > bdim: gshape = [model.shape[0]] + gshape # Space coordinates coordinates = [np.linspace(0, L, N, endpoint=False) for L, N in zip(model.system_size, gshape)] # If model has subsurfce domain, z-coordinate is always first dimension_indices = np.arange(bdim) if len(lsize) > bdim: dimension_indices += 1 dimension_indices = np.concatenate((dimension_indices, [0])) coordinates[0] = \ np.linspace(0, model.system_size[0], gshape[0]) offset = np.zeros_like(dimension_indices) offset[0] = mpi.local_offset(gsize) rectgrid = uvw.RectilinearGrid if mpi.size() == 1 \ else uvw.parallel.PRectilinearGrid # Creating rectilinear grid with correct order for components coordlist = [coordinates[i][o:o+lsize[i]] for i, o in zip(dimension_indices, offset)] grid = rectgrid( file_descriptor, coordlist, compression=True, offsets=offset, ) # Iterator over fields we want to dump fields_it = filter(lambda t: t[0] not in self.forbidden_fields, self.get_fields(model).items()) # We make fields periodic for visualization for name, field in fields_it: array = uvw.DataArray(field, dimension_indices, name) grid.addPointData(array) grid.write() @directory_dump('paraview') class UVWGroupDumper(FieldDumper): "Dumper to ParaViewData files" extension = 'pvd' def __init__(self, basename, *fields, **kwargs): """Construct with desired fields""" super(UVWGroupDumper, self).__init__(basename, *fields, **kwargs) subdir = os.path.join('paraview', basename + '-VTR') if not os.path.exists(subdir): makedirs(subdir) self.uvw_dumper = UVWDumper( os.path.join(basename + '-VTR', basename), *fields, **kwargs ) self.group = uvw.ParaViewData(self.file_path, compression=True) def dump_to_file(self, file_descriptor, model): self.group.addFile( self.uvw_dumper.file_path.replace('paraview/', ''), timestep=self.uvw_dumper.count ) self.group.write() self.uvw_dumper.dump(model) except ImportError: pass try: from netCDF4 import Dataset @directory_dump('netcdf') class NetCDFDumper(FieldDumper): """Dumper to netCDF4 files""" extension = "nc" boundary_fields = ['traction', 'gap'] def _file_setup(self, grp, model): grp.createDimension('frame', None) + # Attibutes + for k, v in _get_attributes(model).items(): + setattr(grp, k, v) + # Local dimensions model_dim = len(model.shape) voigt_dim = type_traits[model.type].voigt self._vec = grp.createDimension('spatial', model_dim) self._tens = grp.createDimension('Voigt', voigt_dim) self.model_info = model.global_shape, model.type global_boundary_shape = mpi.global_shape(model.boundary_shape) # Create boundary dimensions for label, size, length in zip( "xy", global_boundary_shape, model.boundary_system_size ): grp.createDimension(label, size) coord = grp.createVariable(label, 'f8', (label,)) coord[:] = np.linspace(0, length, size, endpoint=False) self._create_variables( grp, model, lambda f: _is_surface_field(f[1], model), global_boundary_shape, "xy" ) # Create volume dimension if model.type in {model_type.volume_1d, model_type.volume_2d}: size = model.shape[0] grp.createDimension("z", size) coord = grp.createVariable("z", 'f8', ("z",)) coord[:] = np.linspace(0, model.system_size[0], size) self._create_variables( grp, model, lambda f: not _is_surface_field(f[1], model), model.global_shape, "zxy" ) self.has_setup = True def dump_to_file(self, file_descriptor, model): if mpi.size() > 1: raise RuntimeError("NetCDFDumper does not function " "properly in parallel") mode = 'a' if os.path.isfile(file_descriptor) \ and getattr(self, 'has_setup', False) else 'w' with Dataset(file_descriptor, mode, format='NETCDF4_CLASSIC', parallel=mpi.size() > 1) as rootgrp: if rootgrp.dimensions == {}: self._file_setup(rootgrp, model) if self.model_info != (model.global_shape, model.type): raise Exception("Unexpected model {}".format(model)) self._dump_generic(rootgrp, model) def _create_variables(self, grp, model, predicate, shape, dimensions): field_dim = len(shape) fields = list(filter(predicate, self.get_fields(model).items())) dim_labels = list(dimensions[:field_dim]) for label, data in fields: local_dim = [] # If we have an extra component if data.ndim > field_dim: if data.shape[-1] == self._tens.size: local_dim = [self._tens.name] elif data.shape[-1] == self._vec.size: local_dim = [self._vec.name] grp.createVariable(label, 'f8', ['frame'] + dim_labels + local_dim, zlib=True) def _dump_generic(self, grp, model): fields = self.get_fields(model).items() new_frame = grp.dimensions['frame'].size for label, data in fields: var = grp[label] slice_in_global = (new_frame,) + local_slice(data, model) var[slice_in_global] = np.array(data, dtype=np.double) + def read(self, file_descriptor): + "Create model with last frame" + with Dataset(file_descriptor, 'r', + format='NETCDF4_CLASSIC') as rootgrp: + attrs = {k: getattr(rootgrp, k) for k in rootgrp.ncattrs()} + model = _create_model(attrs) + dims = rootgrp.dimensions.keys() + for k, v in filter(lambda k: k[0] not in dims, + rootgrp.variables.items()): + model[k] = v[-1, :] + return model + except ImportError: pass diff --git a/tests/test_dumper.py b/tests/test_dumper.py index 625d251..7f7fbd1 100644 --- a/tests/test_dumper.py +++ b/tests/test_dumper.py @@ -1,122 +1,123 @@ # -*- coding: utf-8 -*- # @file # LICENSE # # Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published # by the Free Software Foundation, either version 3 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 Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . from __future__ import division import os import tamaas as tm import numpy as np import pytest from tamaas.dumpers import NumpyDumper class Dumper(tm.ModelDumper): """Simple numpy dumper""" def __init__(self): tm.ModelDumper.__init__(self) def dump(self, model): np.savetxt('tractions.txt', np.ravel(model.traction)) np.savetxt('displacement.txt', np.ravel(model.displacement)) def test_dumpers(tamaas_fixture, tmp_path): os.chdir(tmp_path) model = tm.ModelFactory.createModel(tm.model_type.volume_2d, [1., 1., 1.], [16, 4, 8]) dumper = Dumper() np_dumper = NumpyDumper('test_dump', 'traction', 'displacement') model.addDumper(np_dumper) model.dump() model.dump() dumper << model tractions = np.loadtxt('tractions.txt') displacements = np.loadtxt('displacement.txt') assert np.all(tractions.reshape(model.traction.shape) == model.traction) assert np.all(displacements.reshape(model.displacement.shape) == model.displacement) rmodel = np_dumper.read('numpys/test_dump_0000.npz') assert np.all(model.traction == rmodel.traction) assert np.all(model.displacement == rmodel.displacement) assert model.type == rmodel.type assert os.path.isfile('numpys/test_dump_0001.npz') # Protecting test try: from tamaas.dumpers import H5Dumper import h5py def test_h5dumper(tamaas_fixture, tmp_path): os.chdir(tmp_path) model = tm.ModelFactory.createModel(tm.model_type.volume_2d, [1., 1., 1.], [16, 4, 8]) model['displacement'][...] = 3.1415 dumper = H5Dumper('test_hdf5', 'traction', 'displacement') dumper << model assert os.path.isfile('hdf5/test_hdf5_0000.h5') rmodel = dumper.read('hdf5/test_hdf5_0000.h5') assert np.all(np.abs(rmodel.displacement - 3.1415) < 1e-15) except ImportError: pass try: from tamaas.dumpers import NetCDFDumper from netCDF4 import Dataset def test_netcdfdumper(tamaas_fixture, tmp_path): os.chdir(tmp_path) model = tm.ModelFactory.createModel(tm.model_type.volume_2d, [1., 1., 1.], [16, 4, 8]) model['displacement'][...] = 3.1415 dumper = NetCDFDumper('test_netcdf', 'traction', 'displacement') # Dumping two frames dumper << model dumper << model - with pytest.raises(RuntimeError): - dumper.read('') + rmodel = dumper.read('netcdf/test_netcdf.nc') + assert np.all(np.abs(rmodel.displacement - 3.1415) < 1e-15) assert os.path.isfile('netcdf/test_netcdf.nc') + # Testing all snapshots with Dataset('netcdf/test_netcdf.nc', 'r') as netcdf_file: disp = netcdf_file['displacement'][:] assert np.all(np.abs(disp - 3.1415) < 1e-15) assert disp.shape[1:] == (16, 4, 8, 3) assert disp.shape[0] == 2 with pytest.raises(Exception): model = tm.ModelFactory.createModel(tm.model_type.volume_2d, [1., 1., 1.], [15, 4, 8]) dumper << model except ImportError: pass