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()