Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F92350227
statistics.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Nov 19, 14:50
Size
7 KB
Mime Type
text/x-c++
Expires
Thu, Nov 21, 14:50 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
22425207
Attached To
rTAMAAS tamaas
statistics.cpp
View Options
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "statistics.hh"
#include "fftw_engine.hh"
#include "loop.hh"
#include "static_types.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
template <UInt dim>
Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
return std::sqrt(surface.var());
}
template <UInt dim>
Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
const auto h_size =
GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
auto wavevectors =
FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
wavevectors *= 2 * M_PI; // need q for slopes
const auto psd = computePowerSpectrum(surface);
const Real rms_slope_mean = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(auto q, const auto& 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();
else
return 2 * q.l2squared() * psd_val.real();
},
range<VectorProxy<const Real, dim>>(wavevectors), psd);
return std::sqrt(rms_slope_mean);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
GridHermitian<Real, dim>
Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
const auto h_size =
GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
GridHermitian<Real, dim> 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 <UInt dim>
Grid<Real, dim>
Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
auto psd = computePowerSpectrum(surface);
FFTEngine::makeEngine()->backward(acf, psd);
acf *= acf.getGlobalNbPoints();
return acf;
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real Statistics<dim>::contact(const Grid<Real, dim>& tractions,
UInt perimeter) {
Real points = 0;
UInt nc = tractions.getNbComponents();
switch (nc) {
case 1:
points =
Loop::reduce([](const Real& t) -> Real { return t > 0; }, tractions);
break;
case 2:
points = Loop::reduce([](auto t) -> Real { return t.back() > 0; },
range<VectorProxy<const Real, 2>>(tractions));
break;
case 3:
points = Loop::reduce([](auto t) -> Real { return t.back() > 0; },
range<VectorProxy<const Real, 3>>(tractions));
break;
default:
TAMAAS_EXCEPTION("Invalid number of components in traction");
}
auto area = points / tractions.getNbPoints();
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.getNbPoints()) * perimeter;
}
/* -------------------------------------------------------------------------- */
namespace {
template <UInt dim>
class moment_helper {
public:
moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
CUDA_LAMBDA Complex operator()(VectorProxy<Real, dim> 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;
else
return 2 * mul * phi;
}
private:
std::array<UInt, dim> exponent;
};
} // namespace
template <>
std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
constexpr UInt dim = 1;
std::vector<Real> moments(3);
const auto psd = computePowerSpectrum(surface);
auto wavevectors =
FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
// we don't multiply by 2 pi because moments are computed with k
moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0}}},
range<PVector>(wavevectors), psd)
.real();
moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
range<PVector>(wavevectors), psd)
.real();
moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
range<PVector>(wavevectors), psd)
.real();
return moments;
}
template <>
std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
constexpr UInt dim = 2;
std::vector<Real> moments(3);
const auto psd = computePowerSpectrum(surface);
auto wavevectors =
FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
// we don't multiply by 2 pi because moments are computed with k
moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 0}}},
range<PVector>(wavevectors), psd)
.real();
auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
range<PVector>(wavevectors), psd)
.real();
auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
range<PVector>(wavevectors), psd)
.real();
moments[1] = 0.5 * (m02 + m20);
auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
range<PVector>(wavevectors), psd)
.real();
auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
range<PVector>(wavevectors), psd)
.real();
auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
range<PVector>(wavevectors), psd)
.real();
moments[2] = (3 * m22 + m40 + m04) / 3.;
return moments;
}
template struct Statistics<1>;
template struct Statistics<2>;
} // namespace tamaas
Event Timeline
Log In to Comment