diff --git a/CHANGELOG.md b/CHANGELOG.md
index ef30afb..ae6418f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,345 +1,349 @@
# 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) for final versions and [PEP440](https://www.python.org/dev/peps/pep-0440/) in case intermediate versions need to be released (e.g. development version `2.2.3.dev1` or release candidates `2.2.3rc1`), or individual commits are packaged.
## Unreleased
### Added
- Added `matvec` interface for Hooke and Kelvin operators
### Changed
- Resizing a grid no longer sets values to zero
- Wheels built for `manylinux2014` instead of `manylinux2010`
- Cleaned up `Array` interface
+### Fixed
+
+- Fixed flood fill algorithm to properly account for PBCs
+
## v2.5.0 -- 2022-07-05
### Added
- Added `tamaas.utils.publications` to print a list of publications relevant to
the parts of Tamaas used in a script. Call at the end of a script for an
exhaustive list.
- Added `tamaas.mpi.gather` and `tamaas.mpi.scatter` to help handling of 2D data
in MPI contexts.
- Added `getBoundaryFields()/boundary_fields` function/property to `Model`.
- Added Python constructor to `Model`, which is an alias to
`ModelFactory.createModel`.
- Added convenience `tamaas.dumpers._helper.netCDFtoParaview` function to
convert model sequence to Paraview data file.
- Added dump of surface fields as `FieldData` in VTK files
- Added `scipy.sparse.linalg.LinearOperator` interface to `Westergaard`. A
`LinearOperator` can be created by calling
`scipy.sparse.linalg.aslinearoperator` on the operator object.
- Allowed sub-classing of `ContactSolver` in Python.
- Exposed elastic contact functionals in Python to help with custom solver
implementation.
- Added an example of custom solver with penalty method using Scipy.
- Added `graphArea` function in statistics to compute area of function graph
- Added `tamaas.utils.radial_average` to radialy average a 2D field
- Added option to not make a dump directory (e.g. `netcdf`, `hdf5` directories
for dump files). This allows to specify an absolute path for dump files, and
is useful if simulation workflow is managed with file targets, as do systems
like GNU Make and Snakemake.
### Changed
- Filter for boundary fields in dumper should be smarter.
- `read` functions in dumpers are now `@classmethod`
- Changed sign of plastic correction in KatoSaturated solver to match the
EPICSolver sign.
- Plastic correction of KatoSaturated is registered in the model as
`KatoSaturated::residual_displacement`.
- Changed to modern Python packaging (PEP517) with `setup.cfg` and
`pyproject.toml`. A `setup.py` script is still provided for editable
installations, see [gh#7953](https://github.com/pypa/pip/issues/7953). Once
setuptools has stable support for
[PEP621](https://peps.python.org/pep-0621/), `setup.cfg` will be merged into
`pyproject.toml`. Note that source distributions that can be built with
[`build`](https://pypi.org/project/build/) are not true source distributions,
since they contain the compiled Python module of Tamaas.
- Updated copyright to reflect Tamaas' development history since leaving EPFL
### Fixed
- Fixed `tamaas.dumpers.NetCDFDumper` to dump in MPI context.
### Deprecated
- SCons versions < 3 and Python < 3.6 are explicitely no longer supported.
## v2.4.0 -- 2022-03-22
### Added
- Added a `tamaas.utils` module.
- Added `tamaas.utils.load_path` generator, which yields model objects for a
sequence of applied loads.
- Added `tamaas.utils.seeded_surfaces` generator, which yields surfaces for a
sequence of random seeds.
- Added `tamaas.utils.hertz_surface` which generates a parabolic (Hertzian)
surface.
- Added `tamaas.utils.log_context` context manager which locally sets the log
level for Tamaas' logger.
- Added `tamaas.compute.from_voigt` to convert Voigt/Mendel fields to dense
tensor fields.
- Added a deep-copy function to `ModelFactory` to copy `Model` objects. Use
`model_copy = copy.deepcopy(model)` in Python to make a copy. Currently only
copies registered fields, same as dumpers/readers.
- Added a `read_sequence` method for dumpers to read model frames.
- Automatic draft release on Zenodo.
### Changed
- `*_ROOT` variables for vendored dependencies GoogleTest and pybind11 now
default to empty strings, so that the vendored trees in `third-party/` are not
selected by default. This is so system packages are given priority. Vendored
dependency submodules will eventually be deprecated.
### Fixed
- Fixed an issue with `scons -c` when GTest was missing.
## v2.3.1 -- 2021-11-08
### Added
- Now using `clang-format`, `clang-tidy` and `flake8` for linting
- Pre-built Tamaas images can be pulled with `docker pull
registry.gitlab.com/tamaas/tamaas`
- Added a `--version` flag to the `tamaas` command
### Changed
- The root `Dockerfile` now compiles Tamaas, so using Tamaas in Docker is easier.
- The model attributes dumped to Numpy files are now written in a JSON-formatted
string to avoid unsafe loading/unpickling of objects.
- Removed the `build_doc` build option: now the doc targets are automatically
added if the dependencies are met, and built if `scons doc` is called.
- Removed the `use_googletest` build option: if tests are built and gtest is
present, the corresponding tests will be built
### Fixed
- The command `tamaas plot` gracefully fails if Matplotlib is missing
- Better heuristic to guess which fields are defined on boundary in dumpers
(still not perfect and may give false negatives)
## v2.3.0 -- 2021-06-15
### 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
- Added a JSON encoder class for models and a JSON dumper
- CUDA compatibility is re-established, but has not been tested
- Docstrings in the Python bindings for many classes/methods
### Changed
- Tamaas version numbers are now managed by
[versioneer](https://github.com/python-versioneer/python-versioneer). This
means that Git tags prefixed with `v` (e.g. `v2.2.3`) carry meaning and
determine the version. When no tag is set, versioneer uses the last tag,
specifies the commit short hash and the distance to the last tag (e.g.
`2.2.2+33.ge314b0e`). This version string is used in the compiled library, the
`setup.py` script and the `__version__` variable in the python module.
- Tamaas migrated to [GitLab](https://gitlab.com/tamaas/tamaas)
- Continuous delivery has been implemented:
- the `master` branch will now automatically build and publish Python wheels
to `https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple`. These
"nightly" builds can be installed with:
pip install \
--extra-index-url https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple \
tamaas
- version tags pushed to `master` will automatically publish the wheels to
[PyPI](https://pypi.org/project/tamaas/)
### Deprecated
- The `finalize()` function is now deprecated, since it is automatically called
when the process terminates
- Python versions 3.5 and below are not supported anymore
### Fixed
- Fixed a host of dump read/write issues when model type was not `volume_*d`.
Dumper tests are now streamlined and systematic.
- Fixed a bug where `Model::solveDirichlet` would not compute correctly
- Fixed a bug where `Statistics::contact` would not normalize by the global
number of surface points
## 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 c4507b4..3db7e28 100644
--- a/src/percolation/flood_fill.cpp
+++ b/src/percolation/flood_fill.cpp
@@ -1,270 +1,276 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* 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 {
// Both i and n should be signed, otherwise weird things happen
Int wrap_pbc(Int i, Int n) { return ((i % n) + n) % n; }
template
std::array wrap_pbc(const std::array& t,
const std::array& n) {
std::array index;
for (UInt i = 0; i < dim; ++i)
index[i] = static_cast(wrap_pbc(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] = static_cast(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& /*p*/) {
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
typename Cluster::BBox Cluster::boundingBox() const {
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
std::array Cluster::extent() const {
const auto bb = boundingBox();
std::array extent;
std::transform(bb.first.begin(), bb.first.end(), bb.second.begin(),
extent.begin(),
[](auto min, auto max) { return max - min + 1; });
return extent;
}
template
Cluster::Cluster(Point start, const Grid& map,
Grid& visited, bool diagonal) {
// Visiting queue: a queue here prioritizes visiting local neighbors
// instead of growing a stack with each new visited point it gets rid of older
// neighbors first this prevents the points positions to be shifted
// indefinitely across periodic boundaries
std::queue visiting;
visiting.push(start);
// Contact sizes
const auto& n = cast_int(map.sizes());
while (not visiting.empty()) {
auto p = visiting.front();
auto p_wrapped = wrap_pbc(p, n);
if (visited(p_wrapped)) {
visiting.pop();
continue;
}
visited(p_wrapped) = true;
points.push_back(visiting.front());
visiting.pop();
auto process = [&](const std::vector& neighbors,
const bool do_perimeter) {
for (const auto& p : neighbors) {
const auto index = wrap_pbc(p, n);
if (not visited(index) and map(index)) {
visiting.push(p);
}
perimeter += do_perimeter and (not map(index));
}
};
process(getNextNeighbors(p), true);
if (diagonal) {
process(getDiagonalNeighbors(p), false);
}
}
}
typename FloodFill::List<1> FloodFill::getSegments(const Grid& map) {
auto n = map.sizes();
Grid visited(n, 1);
visited = false;
List<1> clusters;
for (UInt i = 0; i < n[0]; ++i) {
if (map(i) && !visited(i)) {
clusters.emplace_back(std::array{(Int)i}, map, visited, false);
}
}
return clusters;
}
typename FloodFill::List<2> FloodFill::getClusters(const Grid& map,
bool diagonal) {
auto global_contact = Partitioner<2>::gather(map);
auto n = global_contact.sizes();
Grid visited(n, 1);
visited = false;
List<2> clusters;
- if (mpi::rank() == 0)
+ if (mpi::rank() == 0) {
+ Logger().get(LogLevel::debug) << "Flood fill start\n";
+
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);
}
}
}
+ Logger().get(LogLevel::debug)
+ << "Flood fill end, found " << clusters.size() << " clusters\n";
+ }
+
return clusters;
}
typename FloodFill::List<3> FloodFill::getVolumes(const Grid& map,
bool diagonal) {
auto n = map.sizes();
Grid visited(n, 1);
visited = false;
List<3> 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 (map(i, j, k) && !visited(i, j, k)) {
clusters.emplace_back(std::array{(Int)i, (Int)j, (Int)k}, map,
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 c08cc15..427e3cc 100644
--- a/tests/mpi_routines.py
+++ b/tests/mpi_routines.py
@@ -1,307 +1,307 @@
# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
# Copyright (©) 2020-2022 Lucas Frérot
#
# 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 .
import os
import tempfile
import tamaas as tm
import numpy as np
from mpi4py import MPI
from numpy.linalg import norm
from conftest import HertzFixture
from tamaas.dumpers import MPIIncompatibilityError
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 = 1
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_surface = tm.mpi.scatter(surface)
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 = map(tm.mpi.gather, solve())
if comm.rank == 0:
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)
seq_area = tm.Statistics2D.contact(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']
area = tm.Statistics2D.contact(tractions)
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
assert np.abs(seq_area - area) / seq_area < 1e-13
def model_for_dump():
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
return model
def mpi_h5dump():
try:
from tamaas.dumpers import H5Dumper as Dumper
import h5py
except ImportError as e:
print(e)
return
os.chdir(MPI.COMM_WORLD.bcast(tempfile.mkdtemp(), root=0))
model = model_for_dump()
dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction')
try:
dumper_helper << model
except MPIIncompatibilityError as e:
print(e)
return
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)
def mpi_netcdfdump():
try:
from tamaas.dumpers import NetCDFDumper as Dumper
import netCDF4 as nc
except ImportError as e:
print(e)
return
os.chdir(MPI.COMM_WORLD.bcast(tempfile.mkdtemp(), root=0))
model = model_for_dump()
dumper_helper = Dumper('test_mpi_dump', 'displacement', 'traction')
try:
dumper_helper << model
except MPIIncompatibilityError as e:
print(e)
return
rmodel = dumper_helper.read('netcdf/test_mpi_dump.nc')
assert np.all(rmodel.displacement.astype(int) == MPI.COMM_WORLD.rank)
assert np.all(rmodel.traction.astype(int) == MPI.COMM_WORLD.rank)
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"
+ assert len(clusters) == 2, "Wrong number of clusters"
+ assert [c.getArea() for c in clusters] == [2, 2], "Wrong areas"
if __name__ == '__main__':
mpi_polonsky_keer_compare()