diff --git a/CHANGELOG.md b/CHANGELOG.md index 1877b5e..c95df7b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,179 +1,180 @@ # 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 +- `getClusters()` can be called in MPI contact with partial contact maps ## 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/src/percolation/flood_fill.cpp b/src/percolation/flood_fill.cpp index 22bda18..447e5ba 100644 --- a/src/percolation/flood_fill.cpp +++ b/src/percolation/flood_fill.cpp @@ -1,255 +1,258 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "flood_fill.hh" +#include "partitioner.hh" #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { Int unsigned_modulo(Int i, UInt n) { return (i % n + n) % n; } template std::array unsigned_modulo(const std::array& t, const std::array& n) { std::array index; for (UInt i = 0; i < dim; ++i) index[i] = unsigned_modulo(t[i], n[i]); return index; } template std::array cast_int(const std::array& a) { std::array r; for (UInt i = 0; i < dim; ++i) r[i] = (Int)a[i]; return r; } template Cluster::Cluster(const Cluster& other) : points(other.points), perimeter(other.perimeter) {} template <> auto Cluster<1>::getNextNeighbors(const std::array& p) { std::vector neighbors(2); Int i = p[0]; neighbors[0] = {i - 1}; neighbors[1] = {i + 1}; return neighbors; } template <> auto Cluster<1>::getDiagonalNeighbors(const std::array&) { return std::vector(); } template <> auto Cluster<2>::getNextNeighbors(const std::array& p) { std::vector neighbors(4); Int i = p[0], j = p[1]; neighbors[0] = {i + 1, j}; neighbors[1] = {i - 1, j}; neighbors[2] = {i, j - 1}; neighbors[3] = {i, j + 1}; return neighbors; } template <> auto Cluster<2>::getDiagonalNeighbors(const std::array& p) { std::vector neighbors(4); Int i = p[0], j = p[1]; neighbors[0] = {i + 1, j + 1}; neighbors[1] = {i - 1, j - 1}; neighbors[2] = {i - 1, j + 1}; neighbors[3] = {i + 1, j - 1}; return neighbors; } template <> auto Cluster<3>::getNextNeighbors(const std::array& p) { std::vector neighbors(6); Int i = p[0], j = p[1], k = p[2]; neighbors[0] = {i + 1, j, k}; neighbors[1] = {i - 1, j, k}; neighbors[2] = {i, j - 1, k}; neighbors[3] = {i, j + 1, k}; neighbors[4] = {i, j, k - 1}; neighbors[5] = {i, j, k + 1}; return neighbors; } template <> auto Cluster<3>::getDiagonalNeighbors(const std::array& p) { std::vector neighbors(20); Int i = p[0], j = p[1], k = p[2]; // 8 corners neighbors[0] = {i + 1, j + 1, k + 1}; neighbors[1] = {i + 1, j + 1, k - 1}; neighbors[2] = {i + 1, j - 1, k + 1}; neighbors[3] = {i + 1, j - 1, k - 1}; neighbors[4] = {i - 1, j + 1, k + 1}; neighbors[5] = {i - 1, j + 1, k - 1}; neighbors[6] = {i - 1, j - 1, k + 1}; neighbors[7] = {i - 1, j - 1, k - 1}; // 4 diagonals in i = 0 neighbors[8] = {i, j + 1, k + 1}; neighbors[9] = {i, j + 1, k - 1}; neighbors[10] = {i, j - 1, k + 1}; neighbors[11] = {i, j - 1, k - 1}; // 4 diagonals in j = 0 neighbors[12] = {i + 1, j, k + 1}; neighbors[13] = {i + 1, j, k - 1}; neighbors[14] = {i - 1, j, k + 1}; neighbors[15] = {i - 1, j, k - 1}; // 4 diagonals in k = 0 neighbors[16] = {i + 1, j + 1, k}; neighbors[17] = {i + 1, j - 1, k}; neighbors[18] = {i - 1, j + 1, k}; neighbors[19] = {i - 1, j - 1, k}; return neighbors; } template std::pair, std::array> Cluster::boundingBox() const { // TODO make sure periodic boundaries are correctly unwrapped std::array mins, maxs; mins.fill(std::numeric_limits::max()); maxs.fill(std::numeric_limits::min()); for (const auto& p : points) for (UInt i = 0; i < dim; ++i) { mins[i] = std::min(mins[i], p[i]); maxs[i] = std::max(maxs[i], p[i]); } return std::make_pair(mins, maxs); } template Cluster::Cluster(Point start, const Grid& contact, Grid& visited, bool diagonal) { // Visiting stack std::stack> visiting; visiting.push(start); // Contact sizes const auto& n = contact.sizes(); while (not visiting.empty()) { auto p = unsigned_modulo(visiting.top(), n); if (visited(p)) { visiting.pop(); continue; } visited(p) = true; points.push_back(visiting.top()); visiting.pop(); auto process = [&](const std::vector& neighbors, const bool do_perimeter) { for (auto& p : neighbors) { auto index = unsigned_modulo(p, n); if (not visited(index) and contact(index)) { visiting.push(cast_int(index)); } if (do_perimeter and (not contact(index))) ++perimeter; } }; process(getNextNeighbors(p), true); if (diagonal) { process(getDiagonalNeighbors(p), false); } } } std::list> FloodFill::getSegments(const Grid& contact) { auto n = contact.sizes(); Grid visited(n, 1); visited = false; std::list> clusters; for (UInt i = 0; i < n[0]; ++i) { if (contact(i) && !visited(i)) { clusters.emplace_back(std::array{(Int)i}, contact, visited, false); } } return clusters; } std::list> FloodFill::getClusters(const Grid& contact, bool diagonal) { - auto n = contact.sizes(); + auto global_contact = Partitioner<2>::gather(contact); + auto n = global_contact.sizes(); Grid visited(n, 1); visited = false; std::list> clusters; - for (UInt i = 0; i < n[0]; ++i) { - for (UInt j = 0; j < n[1]; ++j) { - if (contact(i, j) && !visited(i, j)) { - clusters.emplace_back(std::array{(Int)i, (Int)j}, contact, - visited, diagonal); + if (mpi::rank() == 0) + for (UInt i = 0; i < n[0]; ++i) { + for (UInt j = 0; j < n[1]; ++j) { + if (global_contact(i, j) && !visited(i, j)) { + clusters.emplace_back(std::array{(Int)i, (Int)j}, + global_contact, visited, diagonal); + } } } - } return clusters; } std::list> FloodFill::getVolumes(const Grid& contact, bool diagonal) { auto n = contact.sizes(); Grid visited(n, 1); visited = false; std::list> clusters; for (UInt i = 0; i < n[0]; ++i) { for (UInt j = 0; j < n[1]; ++j) { for (UInt k = 0; k < n[2]; ++k) { if (contact(i, j, k) && !visited(i, j, k)) { clusters.emplace_back(std::array{(Int)i, (Int)j, (Int)k}, contact, visited, diagonal); } } } } return clusters; } template class Cluster<1>; template class Cluster<2>; template class Cluster<3>; } // namespace tamaas diff --git a/tests/mpi_routines.py b/tests/mpi_routines.py index 057f1f6..9a9e5b1 100644 --- a/tests/mpi_routines.py +++ b/tests/mpi_routines.py @@ -1,268 +1,280 @@ # -*- 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 print_function import shutil import tamaas as tm import numpy as np from mpi4py import MPI from numpy.linalg import norm from conftest import HertzFixture def make_surface(N): spectrum = tm.Isopowerlaw2D() spectrum.q0 = 2 spectrum.q1 = 2 spectrum.q2 = 16 spectrum.hurst = 0.8 generator = tm.SurfaceGeneratorRandomPhase2D(N) generator.spectrum = spectrum generator.random_seed = 0 return generator.buildSurface() def mpi_surface_generator(): N = [128, 128] tm.set_log_level(tm.LogLevel.debug) seq_surface = None comm = MPI.COMM_WORLD print('[{}] {}'.format(comm.rank, comm.size)) with tm.mpi.sequential(): if comm.rank == 0: seq_surface = make_surface(N) surface = make_surface(N) print('[{}] {}'.format(comm.rank, surface.shape)) recv = comm.gather(surface, root=0) if comm.rank == 0: gsurface = np.concatenate(recv) if False: import matplotlib.pyplot as plt plt.imshow(seq_surface) plt.colorbar() plt.figure() plt.imshow(gsurface) plt.colorbar() plt.show() # I think the only reason this assert works is because the # spectrum is compact in the Fourier domain -> surface is regular assert np.all(seq_surface == gsurface) def mpi_model_creation(): N = [20, 50, 50] S = [1., 1., 1.] comm = MPI.COMM_WORLD def get_discretizations(model): return model.shape, model.boundary_shape model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S[1:], N[1:]) n, bn = get_discretizations(model) n[0] = comm.allreduce(n[0]) bn[0] = comm.allreduce(bn[0]) assert n == N[1:] and bn == N[1:] model = tm.ModelFactory.createModel(tm.model_type.volume_2d, S, N) n, bn = get_discretizations(model) n[1] = comm.allreduce(n[1]) bn[0] = comm.allreduce(bn[0]) assert n == N and bn == N[1:] def mpi_polonsky_keer(): N = [1024, 1024] S = [1., 1.] load = 0.00001 comm = MPI.COMM_WORLD hertz = HertzFixture(N[0], load) surface = hertz.surface tm.set_log_level(tm.LogLevel.debug) def solve(): model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = 1, 0 local_n, local_offset = tm.mpi.local_shape(N), tm.mpi.local_offset(N) local_surface = np.copy( surface[local_offset:local_offset+local_n[0], :]) solver = tm.PolonskyKeerRey(model, local_surface, 1e-15) print('Created solver') solver.solve(load) return np.copy(model['traction']), np.copy(model['displacement']) tractions, displacements = solve() precv = comm.gather(tractions, root=0) drecv = comm.gather(displacements, root=0) if comm.rank == 0: tractions = np.concatenate(precv) displacements = np.concatenate(drecv) perror = norm(tractions - hertz.pressure) / norm(hertz.pressure) derror = norm(displacements - hertz.displacement)\ / norm(hertz.displacement) if False: print(perror, derror) import matplotlib.pyplot as plt plt.imshow(tractions - hertz.pressure) plt.colorbar() plt.show() assert perror < 5e-3 assert derror < 8e-3 def mpi_polonsky_keer_compare(): N = [273, 273] S = [1., 1.] E, nu = 0.2, 0.3 load = 0.1 comm = MPI.COMM_WORLD seq_tractions = None rms = 0 with tm.mpi.sequential(): if comm.rank == 0: model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = E, nu surface = make_surface(N) rms = tm.Statistics2D.computeSpectralRMSSlope(surface) surface /= rms solver = tm.PolonskyKeerRey(model, surface, 1e-15) solver.solve(load) seq_tractions = np.copy(model['traction']) rms = comm.bcast(rms, root=0) model = tm.ModelFactory.createModel(tm.model_type.basic_2d, S, N) model.E, model.nu = E, nu surface = make_surface(N) / rms solver = tm.PolonskyKeerRey(model, surface, 1e-15) solver.solve(load) tractions = model['traction'] recv = comm.gather(tractions, root=0) if comm.rank == 0: tractions = np.concatenate(recv) error = np.linalg.norm(seq_tractions - tractions) / seq_tractions.size if False: print(error) import matplotlib.pyplot as plt plt.imshow(seq_tractions - tractions) plt.colorbar() plt.show() assert error < 1e-13 def mpi_h5dump(): try: from tamaas.dumpers import H5Dumper as Dumper import h5py except ImportError: return if not h5py.get_config().mpi: print('Did not test H5Dumper with MPI') return model_type = tm.model_type.volume_2d discretization = [10, 2, 10] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model['displacement'][:] = MPI.COMM_WORLD.rank model['traction'][:] = MPI.COMM_WORLD.rank dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction') dumper_helper << model with h5py.File('hdf5/test_mpi_dump_0000.h5', 'r', driver='mpio', comm=MPI.COMM_WORLD) as handle: assert np.all(handle['displacement'][:, 0, :] == 0) assert np.all(handle['displacement'][:, 1, :] == 1) assert np.all(handle['traction'][0, :] == 0) assert np.all(handle['traction'][1, :] == 1) rmodel = dumper_helper.read('hdf5/test_mpi_dump_0000.h5') assert np.all(rmodel.displacement == MPI.COMM_WORLD.rank) assert np.all(rmodel.traction == MPI.COMM_WORLD.rank) if tm.mpi.rank() == 0: shutil.rmtree('hdf5') def mpi_plastic_solve(): from conftest import UniformPlasticity from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver patch_isotropic_plasticity = UniformPlasticity(tm.model_type.volume_2d, [1.] * 3, [4] * 3) model = patch_isotropic_plasticity.model residual = patch_isotropic_plasticity.residual applied_pressure = 0.1 solver = Solver(residual) solver.tolerance = 1e-15 pressure = model['traction'][..., 2] pressure[:] = applied_pressure solver.solve() solver.updateState() solution, normal = patch_isotropic_plasticity.solution(applied_pressure) for key in solution: error = norm(model[key] - solution[key]) / normal[key] assert error < 2e-15 +def mpi_flood_fill(): + contact = np.zeros([10, 1], dtype=np.bool) + contact[:1, :] = True + contact[-1:, :] = True + + clusters = tm.FloodFill.getClusters(contact, False) + + if tm.mpi.rank() == 0: + assert len(clusters) == 3, "Wrong number of clusters" + assert [c.getArea() for c in clusters] == [1, 2, 1], "Wrong areas" + + if __name__ == '__main__': mpi_polonsky_keer_compare()