diff --git a/CHANGELOG.md b/CHANGELOG.md index 48120d6..63d6b9f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,463 +1,467 @@ # 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. ## Unrelease ### Added - Nonperiodic example with true interference computation +- Added a PETSc non-linear solver class (`tm.nonlinear_solvers.PETScSolver`). + Activate with `scons use_petsc=true`. The class accepts PETSc "command-line" + options, so everything from the linear solver to the iterative scheme can be + tuned at runtime. ### Fixed - Removed some warnings on recent compilers - Compilation issues with Thrust + libcudacxx ### Changed - Moved CI to Debian Bookworm (12) - Changed default backends to host when macros are undefined ## v2.7.1 -- 2023-09-20 ### Fixed - Fixed the elastic energy computation for `Isopowerlaw` when H = 0.5 - Fixed compilation with Thrust versions >= 2 ## v2.7.0 -- 2023-07-03 ### Added - Added a `Material` class tree to decouple material behavior from residuals. The `Material` class can be derived in Python for custom material laws. - `Model` can set the integration method for volume operators with the `setIntegrationMethod` method. - `EPICSolver` can be used in `utils.load_path` - `residual.material` gives a reference to the material object used in the residual - Added a DC-FFT/Boussinesq operator for non-periodic elastic solutions - Can now choose integral operator in `PolonskyKeerRey`, so non-periodic contact is possible! - Python sub-classes of `IntegralOperator` are possible - `TAMAAS_MSG` variadic macro which concatenates arguments with file, line no and function info - `assertion_error`, `not_implemented_error`, `nan_error` and `model_type_error` classes for exceptions, which allow for finer exceptions for debugging - SCons targets `clang-format`, `clang-tidy`, `flake8`, `mypy` for linting, and target alias `lint` to run all available lint tools - Added Anderson mixing solver (see [10.1145/321296.321305](https://doi.org/10.1145/321296.321305) and [10.1006/jcph.1996.0059](https://doi.org/10.1006/jcph.1996.0059) for more details). - `max_iter` property on `EPICSolver` and `DFSANESolver` - Sequential `NewtonRaphsonSolver` using BiCGStab for linear solves - `elasticEnergy` to `Isopowerlaw` to compute full contact energy from PSD ### Changed - Removed `model_type` template for `Residual`: they are only defined on the `model_type::volume_2d` type. - Simplified interface for `Residual`: it now takes a `Material` object at construction which handles the computation of eigenstresses involved in the residual equation - Added function name to `TAMAAS_MSG` - `TAMAAS_ASSERT` raises `assertion_error` and is not disabled in release builds: it should be used to check class invariants - `TAMAAS_EXCEPTION` calls have been replaced by either `throw ` or `TAMAAS_ASSERT` calls, so that specific more error types can be used instead of the catch-all `tamaas::Exception` type - `DFSANESolver` does not throw if non-converged - `scons dev` no longer forces installation in `~/.local`. Run in virtual environment to recover expected behavior. See [FAQ](https://tamaas.readthedocs.io/en/latest/faq.html) for more details. ### Fixed - Random phases set to zero where the C2R transform expects real inputs - Issue where large HDF5 files could not be read in parallel - Copy constructor for `Grid` ### Deprecated - `getStress()`, `setIntegrationMethod()` and `getVector()` in `Residual` are deprecated, use `model["stress"]`, `model.setIntegrationMethod()` and `residual.vector` respectively instead. - `ModelFactor.createResidual` is deprecated. Instead, instanciate a `Material` object and use the `Residual` constructor to create the residual object. ### Removed - `tamaas::Exception` class - `TAMAAS_EXCEPTION`, `TAMAAS_DEBUG_MSG`, `TAMAAS_DEBUG_EXCEPTION` macros - `setup.py` is no longer a viable way of installing the Tamaas python package ## v2.6.0 -- 2022-12-23 ### Added - Added `matvec` interface for the following operators: - `Hooke` - `Kelvin` - `Westergaard` - Documented error handling with MPI - Added `span` class - `FieldContainer` class that defines an interface to access fields registered with an object - Operators now expose their fields, which can be modified from Python - Added `HookeField` operator for heterogeneous elastic properties - Documented proper unit choice - Added `Statistics::computeFDRMSSlope()`, finite difference approximation of RMS of surface slopes - Kato saturated uses L2 norm on residual displacement instead of L∞ - Request system to avoid erasing model data on restart - Boundary field names are dumped ### Changed - Resizing a grid no longer sets values to zero - Wheels built for `manylinux2014` instead of `manylinux2010` - Cleaned up `Array` interface - Shared-memory compile options disabled by default - Unknown build variables raise a warning - Removed conditional branches from PKR - NetCDF dumps in `NETCDF4` format instead of `NETCDF4_CLASSIC` - VTK no longer writes FieldData, this was causing issues reading converted files ### 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/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py index 28c8b42..f11095e 100644 --- a/python/tamaas/nonlinear_solvers/__init__.py +++ b/python/tamaas/nonlinear_solvers/__init__.py @@ -1,234 +1,239 @@ # -*- mode:python; coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 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 . """Nonlinear solvers for plasticity problems. Solvers in this module use :py:mod:`scipy.optimize` to solve the implicit non-linear equation for plastic deformations with fixed contact pressures. """ import numpy as np from functools import wraps from scipy.optimize import newton_krylov, root from scipy.sparse.linalg import bicgstab, LinearOperator from .._tamaas import ( EPSolver, Logger, LogLevel, mpi, _tolerance_manager, _DFSANESolver as DFSANECXXSolver, + TamaasInfo ) __all__ = ['NLNoConvergence', 'DFSANESolver', 'DFSANECXXSolver', 'NewtonKrylovSolver', 'ToleranceManager'] +if TamaasInfo.has_petsc: + from .._tamaas import _PETScSolver as PETScSolver # noqa + __all__.append('PETScSolver') + class NLNoConvergence(RuntimeError): """Convergence not reached exception.""" class NewtonRaphsonSolver(EPSolver): """Basic Newton-Raphson scheme""" def __init__(self, residual): super(NewtonRaphsonSolver, self).__init__(residual) self._x = self.getStrainIncrement() self._residual = self.getResidual() self.max_iter = 1000 def solve(self): EPSolver.beforeSolve(self) # assemble RHS self._residual.computeResidual(self._x) rhs = -self._residual.vector.ravel() norm_0 = np.linalg.norm(rhs) / rhs.size if norm_0 < self.tolerance: Logger().get(LogLevel.info) << ( "NewtonRaphson: successful convergence (0 iterations," f" {self.tolerance}, {norm_0})\n") self._residual.computeResidualDisplacement(self._x) return # define linear operator def residual_tangent(x): v = np.zeros_like(x) self._residual.applyTangent(v, x, self._x) return v.ravel() T = LinearOperator(matvec=residual_tangent, shape=[self._x.size] * 2) error = norm_0 nit = 0 while error > self.tolerance and nit < self.max_iter: self._residual.computeResidual(self._x) rhs = -self._residual.vector.ravel() x, error_code = bicgstab(T, rhs, atol=1e-10) self._x += x.reshape(self._x.shape) error = np.linalg.norm(rhs) / rhs.size nit += 1 if error > self.tolerance: msg = "did not converge" else: msg = "successful convergence" Logger().get(LogLevel.info) << ( f"NewtonRaphson: {msg} ({nit} iterations," f" {self.tolerance}, {error:.3e})\n") self._residual.computeResidualDisplacement(self._x) class ScipySolver(EPSolver): """Base class for solvers wrapping SciPy routines.""" def __init__(self, residual, model=None, callback=None): """Construct nonlinear solver with residual. :param residual: plasticity residual object :param model: Deprecated :param callback: passed on to the Scipy solver """ super(ScipySolver, self).__init__(residual) if mpi.size() > 1: raise RuntimeError("Scipy solvers cannot be used with MPI; " "DFSANECXXSolver can be used instead") self.callback = callback self._x = self.getStrainIncrement() self._residual = self.getResidual() self.options = {'ftol': 0, 'fatol': 1e-9} def solve(self): """Solve R(δε) = 0 using a Scipy function.""" # For initial guess, compute the strain due to boundary tractions # self._residual.computeResidual(self._x) # self._x[...] = self._residual.getVector() EPSolver.beforeSolve(self) # Scipy root callback def compute_residual(vec): self._residual.computeResidual(vec) return self._residual.vector.copy() # Solve Logger().get(LogLevel.debug) << \ "Entering non-linear solve\n" self._x[...] = self._scipy_solve(compute_residual) Logger().get(LogLevel.debug) << \ "Non-linear solve returned" # Computing displacements self._residual.computeResidualDisplacement(self._x) def reset(self): """Set solution vector to zero.""" self._x[...] = 0 def _scipy_solve(self, compute_residual): """Actually call the scipy solver. :param compute_residual: function returning residual for given variable """ raise NotImplementedError() class NewtonKrylovSolver(ScipySolver): """Solve using a finite-difference Newton-Krylov method.""" def _scipy_solve(self, compute_residual): try: return newton_krylov(compute_residual, self._x, f_tol=self.tolerance, verbose=False, callback=self.callback) except Exception: raise NLNoConvergence("Newton-Krylov did not converge") class DFSANESolver(ScipySolver): """Solve using a spectral residual jacobianless method. See :doi:`10.1090/S0025-5718-06-01840-0` for details on method and the relevant Scipy `documentation `_ for details on parameters. """ def _scipy_solve(self, compute_residual): solution = root(compute_residual, self._x, method='df-sane', options={'ftol': 0, 'fatol': self.tolerance}, callback=self.callback) Logger().get(LogLevel.info) << \ "DF-SANE/Scipy: {} ({} iterations, {})".format( solution.message, solution.nit, self.tolerance) if not solution.success: raise NLNoConvergence("DF-SANE/Scipy did not converge") return solution.x.copy() def ToleranceManager(start, end, rate): """Decorate solver to manage tolerance of non-linear solve step.""" def actual_decorator(cls): orig_init = cls.__init__ orig_solve = cls.solve orig_update_state = cls.updateState @wraps(cls.__init__) def __init__(obj, *args, **kwargs): orig_init(obj, *args, **kwargs) obj.setToleranceManager(_tolerance_manager(start, end, rate)) @wraps(cls.solve) def new_solve(obj, *args, **kwargs): ftol = obj.tolerance ftol *= rate obj.tolerance = max(ftol, end) return orig_solve(obj, *args, **kwargs) @wraps(cls.updateState) def updateState(obj, *args, **kwargs): obj.tolerance = start return orig_update_state(obj, *args, **kwargs) cls.__init__ = __init__ # cls.solve = new_solve # cls.updateState = updateState return cls return actual_decorator diff --git a/python/tamaas_module.cpp b/python/tamaas_module.cpp index 865a32e..86e62bd 100644 --- a/python/tamaas_module.cpp +++ b/python/tamaas_module.cpp @@ -1,99 +1,100 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 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 "tamaas.hh" #include "tamaas_info.hh" #include "wrap.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace tamaas { namespace py = pybind11; namespace detail { template struct dtype_helper { static const py::dtype dtype; }; template <> const py::dtype dtype_helper::dtype("=f8"); template <> const py::dtype dtype_helper::dtype("=f16"); } // namespace detail /// Creating the tamaas python module PYBIND11_MODULE(_tamaas, mod) { mod.doc() = "Compiled component of Tamaas"; // Wrapping the base methods mod.def("initialize", &initialize, py::arg("num_threads") = 0, "Initialize tamaas with desired number of threads. Automatically " "called upon import of the tamaas module, but can be manually called " "to set the desired number of threads."); mod.def("finalize", []() { PyErr_WarnEx(PyExc_DeprecationWarning, "finalize() is deprecated, it is now automatically called " "at the end of the program", 1); }); // Default dtype of numpy arrays mod.attr("dtype") = detail::dtype_helper::dtype; // Wrapping release information auto info = py::class_(mod, "TamaasInfo"); info.attr("version") = TamaasInfo::version; info.attr("build_type") = TamaasInfo::build_type; info.attr("branch") = TamaasInfo::branch; info.attr("commit") = TamaasInfo::commit; info.attr("diff") = TamaasInfo::diff; info.attr("remotes") = TamaasInfo::remotes; info.attr("has_mpi") = TamaasInfo::has_mpi; + info.attr("has_petsc") = TamaasInfo::has_petsc; info.attr("backend") = TamaasInfo::backend; // Wrapping tamaas components wrap::wrapCore(mod); wrap::wrapPercolation(mod); wrap::wrapSurface(mod); wrap::wrapModel(mod); wrap::wrapSolvers(mod); wrap::wrapCompute(mod); wrap::wrapMPI(mod); wrap::wrapMaterials(mod); /// Wrapping test features wrap::wrapTestFeatures(mod); // Translate exceptions to python py::register_exception(mod, "AssertionError", PyExc_AssertionError); py::register_exception(mod, "NotImplementedError", PyExc_NotImplementedError); py::register_exception(mod, "ModelTypeError", PyExc_TypeError); py::register_exception(mod, "FloatingPointError", PyExc_FloatingPointError); } } // namespace tamaas diff --git a/src/core/mpi_interface.cpp b/src/core/mpi_interface.cpp index fc17dbb..36d172c 100644 --- a/src/core/mpi_interface.cpp +++ b/src/core/mpi_interface.cpp @@ -1,71 +1,75 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 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 "mpi_interface.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { namespace mpi_dummy { comm& comm::world() { +#if defined(TAMAAS_USE_PETSC) && !defined(TAMAAS_USE_MPI) + static comm _world{MPI_COMM_WORLD}; +#else static comm _world; +#endif return _world; } } // namespace mpi_dummy #ifdef TAMAAS_USE_MPI namespace mpi_impl { comm sequential::old_comm{MPI_COMM_NULL}; comm& comm::world() { static comm _world{MPI_COMM_WORLD}; return _world; } // Define type traits for MPI data types #define TYPE(t, mpi_t) \ const MPI_Datatype type_trait::value { mpi_t } TYPE(float, MPI_FLOAT); TYPE(double, MPI_DOUBLE); TYPE(int, MPI_INT); TYPE(unsigned int, MPI_UNSIGNED); TYPE(long double, MPI_LONG_DOUBLE); TYPE(long, MPI_LONG); TYPE(unsigned long, MPI_UNSIGNED_LONG); TYPE(::tamaas::complex, MPI_CXX_FLOAT_COMPLEX); TYPE(::tamaas::complex, MPI_CXX_DOUBLE_COMPLEX); TYPE(::tamaas::complex, MPI_CXX_LONG_DOUBLE_COMPLEX); TYPE(bool, MPI_CXX_BOOL); #undef TYPE // Define type traits for MPI operations #define OPERATION(op, mpi_op) \ const MPI_Op operation_trait::value { mpi_op } OPERATION(plus, MPI_SUM); OPERATION(min, MPI_MIN); OPERATION(max, MPI_MAX); OPERATION(times, MPI_PROD); #undef OPERATION } // namespace mpi_impl #endif } // namespace tamaas diff --git a/src/core/tamaas.cpp b/src/core/tamaas.cpp index 885bb24..8ab9ef7 100644 --- a/src/core/tamaas.cpp +++ b/src/core/tamaas.cpp @@ -1,125 +1,119 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 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 "tamaas.hh" #include "errors.hh" #include "fftw/interface.hh" #include "logger.hh" #include "mpi_interface.hh" #if TAMAAS_LOOP_BACKEND == TAMAAS_LOOP_BACKEND_OMP #include #endif #ifdef TAMAAS_USE_PETSC #include #endif /* -------------------------------------------------------------------------- */ namespace tamaas { void initialize(UInt num_threads) { static bool has_warned = false; mpi::thread provided = mpi::thread::single; if (not mpi::initialized()) { mpi::init_thread(nullptr, nullptr, mpi::thread::multiple, &provided); } bool should_init_threads = (provided > mpi::thread::single); #if TAMAAS_LOOP_BACKEND == TAMAAS_LOOP_BACKEND_OMP if (num_threads) omp_set_num_threads(num_threads); // set user-defined number of threads else num_threads = omp_get_max_threads(); #else if (num_threads != 0) num_threads = 1; #endif Logger().get(LogLevel::debug) << "requested " << num_threads << " threads\n"; #if TAMAAS_FFTW_BACKEND != TAMAAS_FFTW_BACKEND_NONE if (should_init_threads and (not fftw::init_threads())) { throw std::runtime_error(TAMAAS_MSG("FFTW could not initialize threads!")); } else if (not should_init_threads) Logger().get(LogLevel::debug) << "not initializing FFTW threads\n"; #endif if (mpi::initialized()) { if (not has_warned) { Logger().get(LogLevel::warning) << "experimental MPI support\n"; has_warned = true; } fftw::mpi::init(); } if (should_init_threads) { #if TAMAAS_FFTW_BACKEND != TAMAAS_FFTW_BACKEND_NONE Logger().get(LogLevel::debug) << "initializing FFTW with " << num_threads << " threads\n"; fftw::plan_with_nthreads(num_threads); #endif } #ifdef TAMAAS_USE_CUDA Logger().get(LogLevel::warning) << "experimental Cuda support\n"; #endif #ifdef TAMAAS_USE_PETSC - // if no mpi the macro MPI_COMM_WORLD is defined by petsc - PetscCallAbort(MPI_COMM_WORLD, PetscInitialize(nullptr, nullptr, "", "")); + PetscCallAbort(mpi::comm::world(), PetscInitialize(nullptr, nullptr, "", "")); #endif } /* -------------------------------------------------------------------------- */ void finalize() { if (not mpi::finalized()) { #if TAMAAS_BACKEND != TAMAAS_BACKEND_CPP fftw::cleanup_threads(); #endif fftw::mpi::cleanup(); -#ifdef TAMAAS_USE_PETSC - // if no mpi the macro MPI_COMM_WORLD is defined by petsc - PetscCallAbort(MPI_COMM_WORLD, PetscFinalize()); -#endif - mpi::finalize(); } } namespace { /// Manager for initialize + finalize struct entry_exit_points { entry_exit_points() { initialize(); } ~entry_exit_points() { finalize(); } static const entry_exit_points singleton; }; const entry_exit_points entry_exit_points::singleton; } // namespace } // namespace tamaas diff --git a/src/tamaas_info.cpp.in b/src/tamaas_info.cpp.in index 6fac144..652ed97 100644 --- a/src/tamaas_info.cpp.in +++ b/src/tamaas_info.cpp.in @@ -1,22 +1,28 @@ #include "tamaas_info.hh" // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::version = "@build_version@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::build_type = "@build_type@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::branch = "@branch@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::commit = "@commit@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::diff = "@diff@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::remotes = "@remotes@"; // NOLINTNEXTLINE(readability-redundant-string-init) const std::string tamaas::TamaasInfo::backend = "@backend@"; #if defined(TAMAAS_USE_MPI) const bool tamaas::TamaasInfo::has_mpi = true; #else const bool tamaas::TamaasInfo::has_mpi = false; #endif + +#if defined(TAMAAS_USE_PETSC) +const bool tamaas::TamaasInfo::has_petsc = true; +#else +const bool tamaas::TamaasInfo::has_petsc = false; +#endif diff --git a/src/tamaas_info.hh b/src/tamaas_info.hh index bb8f326..df8d026 100644 --- a/src/tamaas_info.hh +++ b/src/tamaas_info.hh @@ -1,43 +1,44 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef TAMAAS_INFO_HH #define TAMAAS_INFO_HH /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace tamaas { struct TamaasInfo { static const std::string version; static const std::string build_type; static const std::string branch; static const std::string commit; static const std::string remotes; static const std::string diff; static const std::string backend; static const bool has_mpi; + static const bool has_petsc; }; -} +} // namespace tamaas /* -------------------------------------------------------------------------- */ -#endif // TAMAAS_INFO_HH +#endif // TAMAAS_INFO_HH diff --git a/tests/test_patch_plasticity.py b/tests/test_patch_plasticity.py index ff73cdc..aa9a799 100644 --- a/tests/test_patch_plasticity.py +++ b/tests/test_patch_plasticity.py @@ -1,55 +1,63 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-20 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 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 pytest import tamaas from numpy.testing import assert_allclose from tamaas.nonlinear_solvers import \ DFSANESolver as PySolver, \ DFSANECXXSolver as CXXSolver, \ NewtonKrylovSolver, \ NewtonRaphsonSolver -@pytest.mark.parametrize('solvers', [CXXSolver, - NewtonKrylovSolver, - PySolver, - NewtonRaphsonSolver]) +solvers = [CXXSolver, + NewtonKrylovSolver, + PySolver, + NewtonRaphsonSolver] + + +if tamaas.TamaasInfo.has_petsc: + from tamaas.nonlinear_solvers import PETScSolver + solvers.append(PETScSolver) + + +@pytest.mark.parametrize('solvers', solvers) def test_patch_plasticity(patch_isotropic_plasticity, solvers): "Test analyitical solution of 1d plasticity" tamaas.set_log_level(tamaas.LogLevel.info) Solver = solvers 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: assert_allclose((model[key] - solution[key]) / normal[key], 0, atol=3e-15)