diff --git a/CHANGELOG.md b/CHANGELOG.md index b5bf3b0..7c6570a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,361 +1,364 @@ # 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 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 ### 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 ### 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/examples/scipy_penalty.py b/examples/scipy_penalty.py index 9d85626..efded2e 100644 --- a/examples/scipy_penalty.py +++ b/examples/scipy_penalty.py @@ -1,103 +1,119 @@ -import operator -import functools +#!/usr/bin/env python3 +# +# Copyright (©) 2016-2022 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 . import tamaas as tm import numpy as np import matplotlib.pyplot as plt import scipy.optimize from tamaas.utils import hertz_surface class GapPenaltyFunctional(tm.Functional): """Quadratic penalty on negative gap functional.""" def __init__(self, penalty): super().__init__() self.penalty = penalty @staticmethod def clip(gap): """Clip negative gap.""" return np.clip(gap, -np.inf, 0) def computeF(self, gap, dual): """Compute penalty.""" g = self.clip(gap) return 0.5 / self.penalty * np.vdot(g, g) / g.size def computeGradF(self, gap, gradient): """Compute gradient.""" gradient[:] += (1 / self.penalty) * self.clip(gap) class ScipyContactSolver(tm.ContactSolver): """Gap-based penalty solver using https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html""" def __init__(self, model, surface, tolerance, penalty): super().__init__(model, surface, tolerance) if model.type != tm.model_type.basic_2d: raise TypeError("model type is not basic_2d") - self.N = functools.reduce(operator.mul, model.shape) self.shape = model.boundary_shape + self.N = np.prod(self.shape) + model['gap'] = np.zeros(list(self.shape) + [1]) model.be_engine.registerDirichlet() self.elastic = tm.ElasticFunctionalGap( self.model.operators['Westergaard::dirichlet'], self.surface) self.penalty = GapPenaltyFunctional(penalty) self.addFunctionalTerm(self.elastic) self.addFunctionalTerm(self.penalty) def solve(self, mean_gap): """Solve with mean gap constraint.""" def cost(gap): gap = gap.reshape(self.shape) grad = np.zeros_like(gap).reshape(self.shape) self.functional.computeGradF(gap, grad) return self.functional.computeF(gap, grad) def jac(gap): gap = gap.reshape(self.shape) grad = np.zeros_like(gap).reshape(self.shape) self.functional.computeGradF(gap, grad) return np.ravel(grad) opt = scipy.optimize.minimize( cost, np.full(self.N, mean_gap), jac=jac, constraints=scipy.optimize.LinearConstraint( np.full((1, self.N), 1 / self.N), mean_gap, mean_gap), method='trust-constr', tol=self.tolerance * mean_gap) print(opt) # Setup solved model self.model['gap'][:] = opt.x.reshape(self.shape) self.model.displacement[:] = self.model['gap'] + self.surface.reshape( self.shape) self.model.solveDirichlet() self.model.traction[:] -= self.model.traction.min() def main(): """Run a simple Hertzian contact with penalty.""" N = 64 model = tm.Model(tm.model_type.basic_2d, [1, 1], [N] * 2) surface = hertz_surface(model.system_size, model.shape, 1) solver = ScipyContactSolver(model, surface, 1e-10, 1) solver.solve(0.04) plt.imshow(model.traction) plt.colorbar() plt.show() if __name__ == "__main__": main() diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp index 05a5a10..7a9f50f 100644 --- a/python/wrap/core.cpp +++ b/python/wrap/core.cpp @@ -1,121 +1,124 @@ /* * 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 "logger.hh" #include "statistics.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ namespace wrap { using namespace py::literals; /* -------------------------------------------------------------------------- */ template void wrapStatistics(py::module& mod) { auto name = makeDimensionName("Statistics", dim); py::class_>(mod, name.c_str()) .def_static("computePowerSpectrum", &Statistics::computePowerSpectrum, py::return_value_policy::move, "Compute PSD of surface") .def_static( "computeAutocorrelation", &Statistics::computeAutocorrelation, py::return_value_policy::move, "Compute autocorrelation of surface") .def_static("computeMoments", &Statistics::computeMoments, "Compute spectral moments") .def_static("computeSpectralRMSSlope", &Statistics::computeSpectralRMSSlope, "Compute hrms' in Fourier space") + .def_static("computeFDRMSSlope", + &Statistics::computeFDRMSSlope, + "Compute hrms' with finite differences") .def_static("computeRMSHeights", &Statistics::computeRMSHeights, "Compute hrms") .def_static("contact", &Statistics::contact, "tractions"_a, "perimeter"_a = 0, "Compute the (corrected) contact area. Permieter is the " "total contact perimeter in number of segments.") .def_static("graphArea", &Statistics::graphArea, "zdisplacement"_a, "Compute area defined by function graph."); } /* -------------------------------------------------------------------------- */ void wrapCore(py::module& mod) { // Exposing logging facility py::enum_(mod, "LogLevel") .value("debug", LogLevel::debug) .value("info", LogLevel::info) .value("warning", LogLevel::warning) .value("error", LogLevel::error); mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); }); mod.def("get_log_level", Logger::getCurrentLevel); py::class_(mod, "Logger") .def(py::init<>()) .def( "get", [](Logger& logger, LogLevel level) -> Logger& { logger.get(level); return logger; }, "Get a logger object for a log level") .def( "__lshift__", [](Logger& logger, std::string msg) -> Logger& { auto level = logger.getWishLevel(); if (level >= Logger::getCurrentLevel() and not(mpi::rank() != 0 and level == LogLevel::info)) { py::print(msg, "file"_a = py::module::import("sys").attr("stderr")); } return logger; }, "Log a message"); wrapStatistics<1>(mod); wrapStatistics<2>(mod); mod.def( "to_voigt", [](const Grid& field) { TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()"); if (field.getNbComponents() == 9) { Grid voigt(field.sizes(), 6); Loop::loop([] CUDA_LAMBDA( MatrixProxy in, SymMatrixProxy out) { out.symmetrize(in); }, range>(field), range>(voigt)); return voigt; } TAMAAS_EXCEPTION("Wrong number of components to symmetrize"); }, "Convert a 3D tensor field to voigt notation", py::return_value_policy::copy); } } // namespace wrap /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index db5c3f3..cb93d0c 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -1,250 +1,275 @@ /* * 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 "statistics.hh" #include "fft_engine.hh" #include "loop.hh" #include "static_types.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { template Real Statistics::computeRMSHeights(Grid& surface) { return std::sqrt(surface.var()); } +template +template +Real Statistics::rmsSlopesFromPSD(const GridHermitian& psd, + const Grid& diff) { + return std::sqrt(Loop::reduce( + [] CUDA_LAMBDA(VectorProxy q, const Complex& psd_val) { + // Checking if we're in the zone that does not have hermitian symmetry + Complex qback{q.back()}, qsq{q.l2squared()}; + auto factor = 2. - static_cast(thrust::abs(qback) < 1e-15); + return factor * thrust::abs(qsq) * psd_val.real(); + }, + range>(diff), psd)); +} + template Real Statistics::computeSpectralRMSSlope(Grid& surface) { const auto h_size = GridHermitian::hermitianDimensions(surface.sizes()); auto wavevectors = FFTEngine::template computeFrequencies(h_size); + + // Diff operator for gradient norm is (qx^2 + qy^2), no need to change wavevectors *= 2 * M_PI; // need q for slopes - const auto psd = computePowerSpectrum(surface); - const Real rms_slope_mean = Loop::reduce( - [] CUDA_LAMBDA(VectorProxy q, const Complex& psd_val) { - // Checking if we're in the zone that does not have hermitian symmetry - if (std::abs(q.back()) < 1e-15) - return q.l2squared() * psd_val.real(); - return 2 * q.l2squared() * psd_val.real(); + return rmsSlopesFromPSD(computePowerSpectrum(surface), wavevectors); +} + +template +Real Statistics::computeFDRMSSlope(Grid& surface) { + const auto h_size = + GridHermitian::hermitianDimensions(surface.sizes()); + auto wavevectors = + FFTEngine::template computeFrequencies(h_size); + + // Diff operator is (exp(iqΔx) - 1)/Δx in each direction + const auto N = Partitioner::global_size(surface); + Loop::loop( + [&N](VectorProxy q) { + for (UInt i = 0; i < dim; ++i) + q(i) = (thrust::polar(1., 2 * M_PI * q(i).real() / N[i]) - 1) * N[i]; }, - range>(wavevectors), psd); + range>(wavevectors)); - return std::sqrt(rms_slope_mean); + return rmsSlopesFromPSD(computePowerSpectrum(surface), wavevectors); } /* -------------------------------------------------------------------------- */ template GridHermitian Statistics::computePowerSpectrum(Grid& surface) { const auto h_size = GridHermitian::hermitianDimensions(surface.sizes()); GridHermitian psd(h_size, surface.getNbComponents()); FFTEngine::makeEngine()->forward(surface, psd); Real factor = 1. / surface.getGlobalNbPoints(); // Squaring the fourier transform of surface and normalizing Loop::loop( [factor] CUDA_LAMBDA(Complex & c) { c *= factor; c *= conj(c); }, psd); return psd; } /* -------------------------------------------------------------------------- */ template Grid Statistics::computeAutocorrelation(Grid& surface) { Grid acf(surface.sizes(), surface.getNbComponents()); auto psd = computePowerSpectrum(surface); FFTEngine::makeEngine()->backward(acf, psd); acf *= acf.getGlobalNbPoints(); return acf; } /* -------------------------------------------------------------------------- */ template Real Statistics::contact(const Grid& tractions, UInt perimeter) { Real points = 0; UInt nc = tractions.getNbComponents(); switch (nc) { case 1: points = Loop::reduce( [] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions); break; case 2: points = Loop::reduce( [] CUDA_LAMBDA(VectorProxy t) -> Real { return t.back() > 0; }, range>(tractions)); break; case 3: points = Loop::reduce( [] CUDA_LAMBDA(VectorProxy t) -> Real { return t.back() > 0; }, range>(tractions)); break; default: TAMAAS_EXCEPTION("Invalid number of components in traction"); } auto area = points / tractions.getGlobalNbPoints(); if (dim == 1) perimeter = 0; // Correction from Yastrebov et al. (Trib. Intl., 2017) // 10.1016/j.triboint.2017.04.023 return area - (M_PI - 1 + std::log(2)) / (24. * tractions.getGlobalNbPoints()) * perimeter; } /* -------------------------------------------------------------------------- */ namespace { template class moment_helper { public: moment_helper(const std::array& exp) : exponent(exp) {} CUDA_LAMBDA Complex operator()(VectorProxy q, const Complex& phi) const { Real mul = 1; for (UInt i = 0; i < dim; ++i) mul *= std::pow(q(i), exponent[i]); // Do not duplicate everything from hermitian symmetry if (std::abs(q.back()) < 1e-15) return mul * phi; return 2 * mul * phi; } private: std::array exponent; }; } // namespace /* -------------------------------------------------------------------------- */ template <> std::vector Statistics<1>::computeMoments(Grid& surface) { constexpr UInt dim = 1; std::vector moments(3); const auto psd = computePowerSpectrum(surface); auto wavevectors = FFTEngine::template computeFrequencies(psd.sizes()); // we don't multiply by 2 pi because moments are computed with k moments[0] = Loop::reduce(moment_helper{{{0}}}, range(wavevectors), psd) .real(); moments[1] = Loop::reduce(moment_helper{{{2}}}, range(wavevectors), psd) .real(); moments[2] = Loop::reduce(moment_helper{{{4}}}, range(wavevectors), psd) .real(); return moments; } template <> std::vector Statistics<2>::computeMoments(Grid& surface) { constexpr UInt dim = 2; std::vector moments(3); const auto psd = computePowerSpectrum(surface); auto wavevectors = FFTEngine::template computeFrequencies(psd.sizes()); // we don't multiply by 2 pi because moments are computed with k moments[0] = Loop::reduce(moment_helper{{{0, 0}}}, range(wavevectors), psd) .real(); auto m02 = Loop::reduce(moment_helper{{{0, 2}}}, range(wavevectors), psd) .real(); auto m20 = Loop::reduce(moment_helper{{{2, 0}}}, range(wavevectors), psd) .real(); moments[1] = (m02 + m20) / 2; auto m22 = Loop::reduce(moment_helper{{{2, 2}}}, range(wavevectors), psd) .real(); auto m40 = Loop::reduce(moment_helper{{{4, 0}}}, range(wavevectors), psd) .real(); auto m04 = Loop::reduce(moment_helper{{{0, 4}}}, range(wavevectors), psd) .real(); moments[2] = (3 * m22 + m40 + m04) / 3; return moments; } /* -------------------------------------------------------------------------- */ template Real Statistics::graphArea(const Grid& zdisplacement) { const auto h_size = GridHermitian::hermitianDimensions(zdisplacement.sizes()); GridHermitian fourier_disp(h_size, 1), fourier_gradient(h_size, dim); Grid gradient(zdisplacement.sizes(), dim); // Compute gradient in Fourier domain FFTEngine::makeEngine()->forward(zdisplacement, fourier_disp); auto wavevectors = FFTEngine::template computeFrequencies(h_size); wavevectors *= 2 * M_PI; Loop::loop( [] CUDA_LAMBDA(VectorProxy q, VectorProxy grad, const Complex& f) { grad = q * f; grad *= Complex{0, 1}; }, range>(wavevectors), range>(fourier_gradient), fourier_disp); FFTEngine::makeEngine()->backward(gradient, fourier_gradient); // Integrate area of graph formula Real factor = Loop::reduce( [] CUDA_LAMBDA(VectorProxy grad) { return std::sqrt(1 + grad.l2squared()); }, range>(gradient)); return factor / gradient.getGlobalNbPoints(); // return 2 * M_PI * factor / gradient.getGlobalNbPoints(); } template struct Statistics<1>; template struct Statistics<2>; } // namespace tamaas diff --git a/src/core/statistics.hh b/src/core/statistics.hh index 90f08a7..e407bef 100644 --- a/src/core/statistics.hh +++ b/src/core/statistics.hh @@ -1,58 +1,65 @@ /* * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef STATISTICS_HH #define STATISTICS_HH /* -------------------------------------------------------------------------- */ #include "fft_engine.hh" #include "grid.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ /** * @brief Suitcase class for all statistics related functions */ template struct Statistics { /// Compute hrms static Real computeRMSHeights(Grid& surface); /// Compute hrms' in fourier space static Real computeSpectralRMSSlope(Grid& surface); + /// Compute hrms' with finite differences + static Real computeFDRMSSlope(Grid& surface); /// Compute PSD of surface static GridHermitian computePowerSpectrum(Grid& surface); /// Compute autocorrelation static Grid computeAutocorrelation(Grid& surface); /// Compute spectral moments static std::vector computeMoments(Grid& surface); /// Compute (corrected) contact area fraction static Real contact(const Grid& tractions, UInt perimeter = 0); /// Compute the area scaling factor of a periodic graph static Real graphArea(const Grid& displacement); using PVector = VectorProxy; + +private: + template + static Real rmsSlopesFromPSD(const GridHermitian& psd, + const Grid& diff); }; /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif // STATISTICS_HH diff --git a/tests/test_surface.py b/tests/test_surface.py index 58fde32..0a23c4c 100644 --- a/tests/test_surface.py +++ b/tests/test_surface.py @@ -1,141 +1,151 @@ # -*- 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 . from __future__ import print_function, division import pytest import numpy as np import tamaas as tm from numpy import pi integrate = pytest.importorskip('scipy.integrate') special = pytest.importorskip('scipy.special') jv = special.jv E = special.ellipe def spectrum(r, isopowerlaw): kl, kr, ks = isopowerlaw.q0, isopowerlaw.q1, isopowerlaw.q2 H = isopowerlaw.hurst C = 1. if r > ks: return 0 elif r > kr: return C * (r / float(kr))**(-2. * (H + 1)) elif r > kl: return C else: return 0 def integrate_bessel(x, spectrum, n): return 2 * pi * integrate.quad( lambda y: y * jv(0, y * x) * spectrum(n * y), 0, np.inf, limit=200)[0] integrate_bessel = np.vectorize(integrate_bessel, otypes=[tm.dtype]) @pytest.fixture( scope="module", params=[tm.SurfaceGeneratorFilter2D, tm.SurfaceGeneratorRandomPhase2D]) def stats(request): iso = tm.Isopowerlaw2D() iso.q0 = 16 iso.q1 = 32 iso.q2 = 64 iso.hurst = 0.8 N = 243 sg = request.param() sg.spectrum = iso sg.shape = [N, N] analytical_stats = { "rms_heights": iso.rmsHeights(), "rms_slopes_spectral": iso.rmsSlopes(), + "rms_slopes_fd": iso.rmsSlopes(), "moments": iso.moments(), "alpha": iso.alpha() } stats = {k: [] for k in analytical_stats} acf = np.zeros(sg.shape, dtype=tm.dtype) realizations = 1000 for i in range(realizations): sg.random_seed = i s = sg.buildSurface() stats['rms_heights'].append(tm.Statistics2D.computeRMSHeights(s)) stats['rms_slopes_spectral'].append( tm.Statistics2D.computeSpectralRMSSlope(s)) + stats['rms_slopes_fd'].append(tm.Statistics2D.computeFDRMSSlope(s)) stats['moments'].append(tm.Statistics2D.computeMoments(s)) acf += tm.Statistics2D.computeAutocorrelation(s) / realizations L = 1 * sg.spectrum.q0 n = sg.shape[0] // 2 x = np.linspace(0, L / 2, n) true_acf = (L / (2 * pi))**2 \ * integrate_bessel(x, lambda x: spectrum(x, sg.spectrum), L / (2 * pi)) return analytical_stats, stats, acf, true_acf -@pytest.mark.parametrize('key', ['rms_heights', 'rms_slopes_spectral']) +quantities = [ + ('rms_heights', 2e-3), + ('rms_slopes_spectral', 3e-3), + ('rms_slopes_fd', 8e-2), +] + + +@pytest.mark.parametrize('key', quantities, ids=lambda k: k[0]) def test_statistics(stats, key): + key, tol = key analytical_stats, stats, _, _ = stats mu = np.mean(stats[key]) error = np.abs(mu - analytical_stats[key]) / analytical_stats[key] - assert error < 3e-3 + assert error < tol def test_autocorrelation(stats): from tamaas.utils import radial_average _, _, acf, true_acf = stats x, y = (np.linspace(-1, 1, acf.shape[i]) for i in range(2)) r = np.linspace(0, 1, true_acf.shape[0]) theta = np.linspace(0, 2 * np.pi, endpoint=False) r_acf = radial_average(x, y, np.fft.fftshift(acf), r, theta, endpoint=False) acf_error = np.sum((r_acf - true_acf)**2) / np.sum(true_acf**2) assert acf_error < 2e-3 def test_graph_area(): N = 100 x = np.linspace(0, 1, N, endpoint=False) f = np.sin(2 * np.pi * x) df = 2 * np.pi * np.cos(2 * np.pi * x) num_area = np.sqrt(1 + df**2).mean() area = 2 * np.sqrt(1 + 4 * np.pi**2) / np.pi * E(4 * np.pi**2 / (1 + 4 * np.pi**2)) fourier_area = tm.Statistics1D.graphArea(f) assert np.abs(num_area - fourier_area) / num_area < 1e-10 assert np.abs(area - fourier_area) / area < 2e-10