diff --git a/src/SConscript b/src/SConscript index e082bc2..4ca9821 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,170 +1,167 @@ # -*- mode:python; coding: utf-8 -*- # vim: set ft=python: # @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 . import os Import('main_env') def prepend(path, list): return [os.path.join(path, x) for x in list] env = main_env.Clone() # Core core_list = """ -fft_plan_manager.cpp -fftransform.cpp -fftransform_fftw.cpp fft_engine.cpp grid.cpp grid_hermitian.cpp statistics.cpp surface.cpp tamaas.cpp legacy_types.cpp loop.cpp computes.cpp logger.cpp """.split() if env['legacy_bem']: core_list.append('surface_statistics.cpp') core_list = prepend('core', core_list) info_file = main_env.Substfile('tamaas_info.cpp', 'tamaas_info.cpp.in') core_list.append(info_file) # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_filter.cpp surface_generator_filter_fft.cpp surface_generator_random_phase.cpp isopowerlaw.cpp regularized_powerlaw.cpp """.split() generator_list = prepend('surface', generator_list) # Lib PERCOLATION percolation_list = """ flood_fill.cpp """.split() percolation_list = prepend('percolation', percolation_list) # BEM PERCOLATION bem_list = """ bem_kato.cpp bem_polonski.cpp bem_gigi.cpp bem_gigipol.cpp bem_penalty.cpp bem_uzawa.cpp bem_fft_base.cpp bem_functional.cpp bem_meta_functional.cpp elastic_energy_functional.cpp exponential_adhesion_functional.cpp squared_exponential_adhesion_functional.cpp maugis_adhesion_functional.cpp complimentary_term_functional.cpp bem_grid.cpp bem_grid_polonski.cpp bem_grid_kato.cpp bem_grid_teboulle.cpp bem_grid_condat.cpp """.split() bem_list = prepend('bem', bem_list) # Model model_list = """ model.cpp model_factory.cpp model_type.cpp model_template.cpp be_engine.cpp westergaard.cpp elastic_functional.cpp meta_functional.cpp adhesion_functional.cpp volume_potential.cpp kelvin.cpp mindlin.cpp boussinesq.cpp elasto_plastic/isotropic_hardening.cpp elasto_plastic/residual.cpp integration/element.cpp """.split() model_list = prepend('model', model_list) # Solvers solvers_list = """ contact_solver.cpp polonsky_keer_rey.cpp kato_saturated.cpp kato.cpp beck_teboulle.cpp condat.cpp polonsky_keer_tan.cpp ep_solver.cpp epic.cpp """.split() solvers_list = prepend('solvers', solvers_list) # GPU API gpu_list = """ fftransform_cufft.cpp """.split() gpu_list = prepend('gpu', gpu_list) # Assembling total list rough_contact_list = \ core_list + generator_list + percolation_list + model_list + solvers_list # Adding legacy code if env['legacy_bem']: rough_contact_list += bem_list # Adding GPU if needed if env['backend'] == 'cuda': rough_contact_list += gpu_list # Adding extra warnings for Tamaas base lib env.AppendUnique(CXXFLAGS=['-Wextra']) # Build library libTamaas = env.SharedLibrary('Tamaas', rough_contact_list, SHLIBVERSION=env['version']) # Defining alias for cpp builds main_env.Alias('build-cpp', libTamaas) # Specify install target to install lib lib_prefix = os.path.join(env['prefix'], 'lib') lib_install = env.InstallVersionedLib(target=lib_prefix, source=libTamaas) main_env.Alias('install-lib', lib_install) # Export target for use in python builds Export('libTamaas') diff --git a/src/bem/bem_interface.hh b/src/bem/bem_interface.hh index 2cc620b..fab64dd 100644 --- a/src/bem/bem_interface.hh +++ b/src/bem/bem_interface.hh @@ -1,100 +1,99 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BEM_INTERFACE_H #define BEM_INTERFACE_H #include "fft_engine.hh" -#include "fftransform.hh" #include "surface.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { class BemInterface { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: BemInterface(Surface& surface) : surface(), spectral_surface(surface.size(), surface.getL()) { this->surface.wrap(surface); this->computeSpectralSurface(); }; virtual ~BemInterface(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: //! compute the displacements from the tractions virtual void computeDisplacementsFromTractions(){ SURFACE_FATAL("TO BE IMPLEMENTED")}; //! compute the tractions from the displacements virtual void computeTractionsFromDisplacements(){ SURFACE_FATAL("TO BE IMPLEMENTED")}; public: //! compute the spectral surface virtual void computeSpectralSurface() { engine.forward(surface, spectral_surface); }; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: //! return the actual tractions virtual const Surface& getTractions() const = 0; //! return the actual displacements virtual const Surface& getDisplacements() const = 0; //! return the manipulated surface const Surface& getSurface() const { return this->surface; }; //! return the manipulated surface in fourier space const SurfaceComplex& getSpectralSurface() const { return this->spectral_surface; }; //! return the manipulated surface void setSurface(Surface& surface) { this->surface = surface; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: //! the surface profile we work with Surface surface; //! the surface profile in fourier space SurfaceComplex spectral_surface; //! fft engine FFTEngine engine; }; } // namespace tamaas #endif /* BEM_INTERFACE_H */ diff --git a/src/bem/bem_polonski.cpp b/src/bem/bem_polonski.cpp index 9cbfbea..454e586 100644 --- a/src/bem/bem_polonski.cpp +++ b/src/bem/bem_polonski.cpp @@ -1,466 +1,464 @@ /** * * @author Guillaume Anciaux * * @section LICENSE * * Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "bem_polonski.hh" -#include "fft_plan_manager.hh" -#include "fftransform.hh" #include "surface.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ BemPolonski::BemPolonski(Surface& p) : BemFFTBase(p), surface_t(p.size(), p.getL()), surface_r(p.size(), p.getL()), pold(p.size(), p.getL()), p_sat(-1.) { e_star = 1.; max_iterations = 2000; surface_rms_heights = SurfaceStatistics::computeStdev(this->surface); } /* -------------------------------------------------------------------------- */ BemPolonski::~BemPolonski() {} /* -------------------------------------------------------------------------- */ Real BemPolonski::computeEquilibrium(Real epsilon, Real pressure) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions); std::cout << "current pressure " << current_pressure << std::endl; if (current_pressure <= 0.) surface_tractions = pressure; this->enforcePressureBalance(pressure); convergence_iterations.clear(); nb_iterations = 0; while (f > epsilon && nb_iterations < max_iterations) { fPrevious = f; this->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); this->updateT(G, Gold, delta); Real tau = this->computeTau(); pold = this->surface_tractions; Gold = G; delta = this->updateTractions(tau); // Projection on admissible space this->enforcePressureBalance(pressure); // std::cout << std::scientific << std::setprecision(10) // << SurfaceStatistics::computeAverage(surface_tractions) << // std::fixed << std::endl; // f = this->computeF(); if (nb_iterations % dump_freq == 0) { Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions); std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f - fPrevious << " " << A << std::fixed << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); return f; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeEquilibriuminit(Real epsilon, Real pressure, Surface& init) { this->computeSpectralInfluenceOverDisplacement(); this->surface_t = 0.; this->surface_r = 0.; this->pold = 0.; this->surface_displacements = 0.; Real delta = 0.; Real Gold = 1.; Real f = 1e300; Real fPrevious = 1e300; this->surface_tractions = init; Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions); if (current_pressure <= 0.) surface_tractions = pressure; this->enforcePressureBalance(pressure); convergence_iterations.clear(); nb_iterations = 0; this->computeDisplacementsFromTractions(); f = this->computeF(); std::cout << std::scientific << std::setprecision(10) << " " << f << std::endl; while (f > epsilon && nb_iterations < max_iterations) { fPrevious = f; this->computeDisplacementsFromTractions(); this->computeGaps(); Real gbar = this->computeMeanGapsInContact(); this->gap -= gbar; Real G = this->computeG(); this->updateT(G, Gold, delta); Real tau = this->computeTau(); pold = this->surface_tractions; Gold = G; delta = this->updateTractions(tau); // Projection on admissible space this->enforcePressureBalance(pressure); // std::cout << std::scientific << std::setprecision(10) // << SurfaceStatistics::computeAverage(surface_tractions) << // std::fixed << std::endl; // f = this->computeF(); if (nb_iterations % dump_freq == 0) { Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions); std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f - fPrevious << " " << A << std::fixed << std::endl; } convergence_iterations.push_back(f); ++nb_iterations; } this->computeTrueDisplacements(); return f; } /* -------------------------------------------------------------------------- */ void BemPolonski::computeGaps() { UInt size = surface.size(); #pragma omp parallel for for (UInt i = 0; i < size * size; i++) { gap(i) = surface_displacements(i) - surface(i); } } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeMeanGapsInContact() { UInt n = surface.size(); UInt size = n * n; Real res = 0.; UInt nb_contact = 0; #pragma omp parallel for reduction(+ : nb_contact, res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0 || std::abs(this->surface_tractions(i) - p_sat) < 1.0e-10) continue; ++nb_contact; res += this->gap(i); } res /= nb_contact; return res; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeG() { UInt n = surface.size(); UInt size = n * n; Real res = 0.; #pragma omp parallel for reduction(+ : res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0 || std::abs(this->surface_tractions(i) - p_sat) < 1.0e-10) continue; Real val = this->gap(i); res += val * val; } return res; } /* -------------------------------------------------------------------------- */ void BemPolonski::updateT(Real G, Real Gold, Real delta) { UInt n = surface.size(); UInt size = n * n; Real factor = delta * G / Gold; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.; else { this->surface_t(i) *= factor; this->surface_t(i) += this->gap(i); } } } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeTau() { this->applyInfluenceFunctions(surface_t, surface_r); const UInt size = surface_t.size() * surface_t.size(); Real rbar = 0; UInt nb_contact = 0; #pragma omp parallel for reduction(+ : nb_contact, rbar) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; ++nb_contact; rbar += surface_r(i); } rbar /= nb_contact; surface_r -= rbar; Real tau_sum1 = 0., tau_sum2 = 0.; #pragma omp parallel for reduction(+ : tau_sum1, tau_sum2) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; tau_sum1 += this->gap(i) * surface_t(i); tau_sum2 += surface_r(i) * surface_t(i); } Real tau = tau_sum1 / tau_sum2; return tau; } /* -------------------------------------------------------------------------- */ Real BemPolonski::updateTractions(Real tau) { UInt n = surface.size(); UInt size = n * n; #pragma omp parallel for for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0) continue; this->surface_tractions(i) -= tau * this->surface_t(i); if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0; } // compute number of interpenetration without contact UInt nb_iol = 0; #pragma omp parallel for reduction(+ : nb_iol) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.) { this->surface_tractions(i) -= tau * this->gap(i); ++nb_iol; } } // compute number of interpenetration without contact UInt nb_sl = 0; #pragma omp parallel for reduction(+ : nb_sl) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) == p_sat && 0. < this->gap(i)) { this->surface_tractions(i) -= tau * this->gap(i); ++nb_sl; } } Real delta = 0; if (nb_iol > 0 || nb_sl > 0) delta = 0.; else delta = 1.; return delta; } /* -------------------------------------------------------------------------- */ void BemPolonski::enforcePressureBalance(Real applied_pressure) { UInt n = surface.size(); UInt size = n * n; // If we have no saturation, scale to match applied_pressure if (this->p_sat <= 0. || SurfaceStatistics::computeMaximum( this->surface_tractions) < this->p_sat) { Real pressure = 0.; #pragma omp parallel for reduction(+ : pressure) for (UInt i = 0; i < size; ++i) { pressure += this->surface_tractions(i); } pressure *= 1. / size; this->surface_tractions *= applied_pressure / pressure; } // If we have saturation, use secant method to find zero of saturationFunction else { // Checking if we can find a zero Real limit = this->p_sat * SurfaceStatistics::computeContactArea( this->surface_tractions) - applied_pressure; if (limit < 0) SURFACE_FATAL("Saturation pressure is too small for applied pressure"); // Initial points Real x_n_2 = 0., x_n_1 = 1., x_n = 0.; Real f_n_2 = 0., f_n_1 = 0., f_n = 0.; // Secant loop do { f_n_2 = saturationFunction(x_n_2, applied_pressure); f_n_1 = saturationFunction(x_n_1, applied_pressure); x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2); f_n = saturationFunction(x_n, applied_pressure); x_n_2 = x_n_1; x_n_1 = x_n; } while (std::abs(f_n) > 1e-16); this->surface_tractions *= x_n; // Truncating pressures #pragma omp parallel for for (UInt i = 0; i < size; i++) { if (this->surface_tractions(i) > this->p_sat) this->surface_tractions(i) = this->p_sat; } // Real mu_I_sat=0.; // #pragma omp parallel for reduction(+: mu_I_sat) // for (UInt i = 0; i < size; ++i) { // if ( this->p_sat< this->surface_tractions(i) ){ // mu_I_sat+=1.; // } // } // Real pc=0.; // #pragma omp parallel for reduction(+: pc) // for (UInt i = 0; i < size; ++i) { // if ( this->surface_tractions(i) < this->p_sat ){ // pc+=this->surface_tractions(i);} // } // #pragma omp parallel for // for (UInt i = 0; i < size; ++i) { // if ( this->surface_tractions(i) < this->p_sat ){ // this->surface_tractions(i)= // this->surface_tractions(i)*((size*applied_pressure-this->p_sat*mu_I_sat)/pc);} // } // #pragma omp parallel for // for (UInt i = 0; i < size; ++i) { // if ( this->p_sat< this->surface_tractions(i) ){ // this->surface_tractions(i)=this->p_sat; // } // } } } /* -------------------------------------------------------------------------- */ Real BemPolonski::saturationFunction(Real alpha, Real applied_pressure) { const UInt n = surface.size(); const UInt size = n * n; Real sum = 0.; #pragma omp parallel for reduction(+ : sum) for (UInt i = 0; i < size; i++) { Real alpha_p = alpha * this->surface_tractions(i); Real alpha_p_I = (alpha_p > this->p_sat) ? alpha_p - this->p_sat : 0.; sum += alpha_p - alpha_p_I; } sum /= size; return sum - applied_pressure; } /* -------------------------------------------------------------------------- */ Real BemPolonski::computeF() { UInt n = surface.size(); UInt size = n * n; Real res = 0; Real t_sum = SurfaceStatistics::computeSum(surface_tractions); computeTrueDisplacements(); #pragma omp parallel for reduction(+ : res) for (UInt i = 0; i < size; ++i) { if (this->surface_tractions(i) == this->p_sat) continue; res += surface_tractions(i) * (this->true_displacements(i) - surface(i)); // res += // this->surface_tractions[i].real() // *(surface_displacements[i].real() - surface[i].real()); } return std::abs(res / t_sum / surface_rms_heights / size); } /* -------------------------------------------------------------------------- */ void BemPolonski::computeTrueDisplacements() { this->applyInfluenceFunctions(this->surface_tractions, this->surface_displacements); this->true_displacements = this->surface_displacements; UInt n = surface.size(); UInt size = n * n; Real shift = 1e300; #pragma omp parallel for reduction(min : shift) for (UInt i = 0; i < size; ++i) { if (surface_displacements(i) - shift - surface(i) < 0. && (this->p_sat < 0 || this->surface_tractions(i) < this->p_sat)) { shift = surface_displacements(i) - surface(i); } } #pragma omp parallel for for (UInt i = 0; i < size; ++i) { true_displacements(i) = surface_displacements(i) - shift; } } /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/core/fft_plan_manager.cpp b/src/core/fft_plan_manager.cpp deleted file mode 100644 index f829d7d..0000000 --- a/src/core/fft_plan_manager.cpp +++ /dev/null @@ -1,98 +0,0 @@ -/** - * @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 . - * - */ -/* -------------------------------------------------------------------------- */ -#include "fft_plan_manager.hh" -#include "fftransform_fftw.hh" -#if defined(USE_CUFFT) -#include "fftransform_cufft.hh" -#endif -#include -/* -------------------------------------------------------------------------- */ - -namespace tamaas { -/* -------------------------------------------------------------------------- */ - -FFTPlanManager::~FFTPlanManager() { clean(); } - -/* -------------------------------------------------------------------------- */ - -namespace detail { -/// class purely for singleton instanciation -struct PublicFFTPlanManager : public FFTPlanManager {}; -} // namespace detail - -FFTPlanManager& FFTPlanManager::get() { - if (FFTPlanManager::singleton == nullptr) - FFTPlanManager::singleton = - std::make_unique(); - - return *FFTPlanManager::singleton; -} -/* -------------------------------------------------------------------------- */ - -void FFTPlanManager::clean() { - for (auto pair : plans) { - delete std::get<0>(pair.second); - delete std::get<1>(pair.second); - } - plans.clear(); -} - -/* -------------------------------------------------------------------------- */ - -template <> -void FFTPlanManager::destroyPlan(const Real* some) { - std::list pairs; - for (auto & pair : plans) { - if (pair.first.first == some) { // C++14 use get(pair.first) - delete std::get<0>(pair.second); - delete std::get<1>(pair.second); - pairs.push_back(pair.first); - } - } - - for (auto & pair : pairs) - plans.erase(pair); -} - -/* -------------------------------------------------------------------------- */ - -template <> -void FFTPlanManager::destroyPlan(const Complex* some) { - std::list pairs; - for (auto & pair : plans) { - if (pair.first.second == some) { - delete std::get<0>(pair.second); - delete std::get<1>(pair.second); - pairs.push_back(pair.first); - } - } - - for (auto & pair : pairs) - plans.erase(pair); -} - -/* -------------------------------------------------------------------------- */ -std::unique_ptr FFTPlanManager::singleton = nullptr; -/* -------------------------------------------------------------------------- */ - -} // namespace tamaas diff --git a/src/core/fft_plan_manager.hh b/src/core/fft_plan_manager.hh deleted file mode 100644 index 6256629..0000000 --- a/src/core/fft_plan_manager.hh +++ /dev/null @@ -1,132 +0,0 @@ -/** - * @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 . - * - */ -/* -------------------------------------------------------------------------- */ -#ifndef FFT_PLAN_MANAGER_H -#define FFT_PLAN_MANAGER_H -/* -------------------------------------------------------------------------- */ -#include "fftransform.hh" -#include "fftransform_fftw.hh" -#include -#ifdef USE_CUDA -#include "fftransform_cufft.hh" -#endif -/* -------------------------------------------------------------------------- */ -namespace tamaas { -/* -------------------------------------------------------------------------- */ - -/** - * @brief Singleton class for FFT plan management - */ -class FFTPlanManager { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -protected: - FFTPlanManager() = default; - -public: - ~FFTPlanManager(); - -public: - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - - /// Get singleton instance - static FFTPlanManager& get(); - - /// Create/retrieve a plan from two surfaces - template - FFTransform& createPlan(Grid& input, - GridHermitian& output); - - /// Remove all plans - void clean(); - - /// Destroy a plan from two surfaces - template - void destroyPlan(Grid& input, GridHermitian& output); - - /// Destroy any plan containing a given data pointer - template - void destroyPlan(const T* some); - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -private: - using FFTMap = - std::map, - std::tuple*, FFTransform*>>; - FFTMap plans; - static std::unique_ptr singleton; -}; - -/* -------------------------------------------------------------------------- */ - -template -FFTransform& -FFTPlanManager::createPlan(Grid& input, - GridHermitian& output) { - static_assert(dim <= 2, "Cannot do FFT of dimension higher than 2"); - - auto index = std::make_pair(const_cast(input.getInternalData()), - const_cast(output.getInternalData())); - - auto it = plans.find(index); - auto end = plans.end(); - - if (it == end) { // we need to create a plan - std::get(plans[index]) = -#ifdef USE_CUDA - new FFTransformCUFFT(input, output); -#else - new FFTransformFFTW(input, output); -#endif - } - - return *std::get(plans[index]); -} - -/* -------------------------------------------------------------------------- */ - -template -void FFTPlanManager::destroyPlan(Grid& input, - GridHermitian& output) { - - auto index = std::make_pair(const_cast(input.getInternalData()), - const_cast(output.getInternalData())); - - auto it = plans.find(index); - auto end = plans.end(); - - if (it != end) { - delete std::get(plans[index]); - plans.erase(index); - } -} - -/* -------------------------------------------------------------------------- */ -} // namespace tamaas -/* -------------------------------------------------------------------------- */ - -#endif /* FFT_PLAN_MANAGER_H */ diff --git a/src/core/fftransform.cpp b/src/core/fftransform.cpp deleted file mode 100644 index e559656..0000000 --- a/src/core/fftransform.cpp +++ /dev/null @@ -1,55 +0,0 @@ -/** - * @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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -#include "fftransform.hh" - -namespace tamaas { - -template -FFTransform::FFTransform(Grid& real_, - GridHermitian& spectral_) - : real(), spectral() { - this->real.wrap(real_); - this->spectral.wrap(spectral_); -} - -/* -------------------------------------------------------------------------- */ - -template -FFTransform::~FFTransform() = default; - -/* -------------------------------------------------------------------------- */ - -template -void FFTransform::normalize() { - real *= 1. / real.getNbPoints(); -} - -/* -------------------------------------------------------------------------- */ - -template class FFTransform; -template class FFTransform; -template class FFTransform; - -} // namespace tamaas diff --git a/src/core/fftransform.hh b/src/core/fftransform.hh deleted file mode 100644 index 01d6563..0000000 --- a/src/core/fftransform.hh +++ /dev/null @@ -1,113 +0,0 @@ -/** - * @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 . - * - */ - -/* -------------------------------------------------------------------------- */ -#ifndef FFTRANSFORM_HH -#define FFTRANSFORM_HH -/* -------------------------------------------------------------------------- */ -#include "grid.hh" -#include "grid_hermitian.hh" -#include "tamaas.hh" -#include "static_types.hh" -#include -#include -/* -------------------------------------------------------------------------- */ -namespace tamaas { -/* -------------------------------------------------------------------------- */ - -/** - * @brief Generic class for FFT interface wrapping - */ -template -class FFTransform { -public: - /// Constructor - FFTransform(Grid& real, GridHermitian& spectral); - - /// Destructor - virtual ~FFTransform(); - -public: - /// Perform FFT - virtual void forwardTransform() = 0; - - /// Perform IFFT - virtual void backwardTransform() = 0; - - /// Normalize the real surface after backward transform - virtual void normalize(); - -public: - /// Fill a grid with modevectors: boolean if grid has hermitian dimensions - template - static Grid computeFrequencies(const std::array& sizes); - -protected: - Grid real; - GridHermitian spectral; -}; - -/* -------------------------------------------------------------------------- */ - -template -template -Grid -FFTransform::computeFrequencies(const std::array& sizes) { - // If hermitian is true, we suppose the dimensions of freq are - // reduced based on hermitian symetry and that it has dim components - - auto& n = sizes; - Grid freq(n, dim); - constexpr UInt dmax = dim - static_cast(hermitian); - -#pragma omp parallel for // to get rid of a compilation warning from nvcc - for (UInt i = 0; i < freq.getNbPoints(); ++i) { - UInt index = i; - VectorProxy wavevector(freq(index * freq.getNbComponents())); - std::array tuple{{0}}; - /// Computing tuple from index - for (Int d = dim - 1; d >= 0; d--) { - tuple[d] = index % n[d]; - index -= tuple[d]; - index /= n[d]; - } - - if (hermitian) - wavevector(dim - 1) = tuple[dim - 1]; - - for (UInt d = 0; d < dmax; d++) { - // Type conversion - T td = tuple[d]; - T nd = n[d]; - T q = (tuple[d] < n[d] / 2) ? td : td - nd; - wavevector(d) = q; - } - } - - return freq; -} - -/* -------------------------------------------------------------------------- */ -} // namespace tamaas -/* -------------------------------------------------------------------------- */ - -#endif // FFTRANSFORM_HH diff --git a/src/core/fftransform_fftw.cpp b/src/core/fftransform_fftw.cpp deleted file mode 100644 index 87559f6..0000000 --- a/src/core/fftransform_fftw.cpp +++ /dev/null @@ -1,103 +0,0 @@ -/** - * @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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -#include "fftransform_fftw.hh" - -#define TAMAAS_FFTW_FLAG FFTW_ESTIMATE - -namespace tamaas { - -template -FFTransformFFTW::FFTransformFFTW(Grid& real, - GridHermitian& spectral) - : FFTransform(real, spectral), forward_plan(nullptr), - backward_plan(nullptr) { - TAMAAS_ASSERT(real.getNbComponents() == spectral.getNbComponents(), - "components for FFTW do not match"); - - // Data pointers - Real* in = this->real.getInternalData(); - auto* out = - reinterpret_cast(this->spectral.getInternalData()); - int components = real.getNbComponents(); - - // fftw parameters - int rank = dim; // dimension of fft - int n[dim] = {0}; - std::copy(real.sizes().begin(), real.sizes().end(), - n); // size of individual fft - int howmany = components; // how many fft to compute - int idist = 1, odist = 1; // components are next to each other in memory - int istride = real.getStrides().back() * components, - ostride = spectral.getStrides().back() * components; - int *inembed = nullptr, *onembed = nullptr; // row major - -#if TAMAAS_FFTW_FLAG != FFTW_ESTIMATE - // backup the input - Grid backup(real); -#endif - - forward_plan = - fftw::plan_many_forward(rank, n, howmany, in, inembed, istride, idist, - out, onembed, ostride, odist, TAMAAS_FFTW_FLAG); - backward_plan = - fftw::plan_many_backward(rank, n, howmany, out, onembed, ostride, odist, - in, inembed, istride, idist, TAMAAS_FFTW_FLAG); - -#if TAMAAS_FFTW_FLAG != FFTW_ESTIMATE - // restore the backup - real = backup; -#endif -} - -/* -------------------------------------------------------------------------- */ - -template -FFTransformFFTW::~FFTransformFFTW() { - fftw::destroy(forward_plan); - fftw::destroy(backward_plan); -} - -/* -------------------------------------------------------------------------- */ - -template -void FFTransformFFTW::forwardTransform() { - fftw::execute(forward_plan); -} - -/* -------------------------------------------------------------------------- */ - -template -void FFTransformFFTW::backwardTransform() { - fftw::execute(backward_plan); - this->normalize(); -} - -/* -------------------------------------------------------------------------- */ - -template class FFTransformFFTW; -template class FFTransformFFTW; -template class FFTransformFFTW; - -} // namespace tamaas diff --git a/src/core/fftransform_fftw.hh b/src/core/fftransform_fftw.hh deleted file mode 100644 index be67fd2..0000000 --- a/src/core/fftransform_fftw.hh +++ /dev/null @@ -1,62 +0,0 @@ -/** - * @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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -#ifndef FFTRANSFORM_FFTW_HH -#define FFTRANSFORM_FFTW_HH - -#include "fftransform.hh" -#include "fftw_interface.hh" - -namespace tamaas { - -/** - * @brief Class wrapping fftw's interface - */ -template -class FFTransformFFTW : public FFTransform { - using plan_t = typename fftw::helper::plan; - using complex_t = typename fftw::helper::complex; - -public: - /// Constructor - FFTransformFFTW(Grid& real, GridHermitian& spectral); - - /// Destructor - ~FFTransformFFTW() override; - -public: - /// Perform FFT - void forwardTransform() override; - - /// Perform IFFT - void backwardTransform() override; - -private: - plan_t forward_plan; - plan_t backward_plan; -}; - -} // namespace tamaas - -#endif // FFTRANSFORM_FFTW_HH diff --git a/src/core/tamaas.cpp b/src/core/tamaas.cpp index 8ff2c25..ccd5ffa 100644 --- a/src/core/tamaas.cpp +++ b/src/core/tamaas.cpp @@ -1,55 +1,54 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" -#include "fft_plan_manager.hh" #include "loop.hh" #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { void initialize(UInt num_threads) { if (num_threads) omp_set_num_threads(num_threads); // set user-defined number of threads // fftw_import_wisdom_from_filename("fftransform_fftw_wisdom"); #ifndef USE_CUDA if (!fftw_init_threads()) TAMAAS_EXCEPTION("FFTW could not initialize threads!"); fftw_plan_with_nthreads(omp_get_max_threads()); #endif } /* -------------------------------------------------------------------------- */ void finalize() { // fftw_export_wisdom_to_filename("fftransform_fftw_wisdom"); fftw_cleanup_threads(); } } // namespace tamaas diff --git a/src/solvers/contact_solver.cpp b/src/solvers/contact_solver.cpp index 4812cf8..12ef5e9 100644 --- a/src/solvers/contact_solver.cpp +++ b/src/solvers/contact_solver.cpp @@ -1,65 +1,64 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "contact_solver.hh" -#include "fft_plan_manager.hh" #include "logger.hh" #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { ContactSolver::ContactSolver(Model& model, const GridBase& surface, Real tolerance) : model(model), surface(), functional(), tolerance(tolerance), surface_stddev(std::sqrt(surface.var())) { auto discretization = model.getBoundaryDiscretization(); auto n_points = model.getTraction().getNbPoints(); if (n_points != surface.dataSize()) TAMAAS_EXCEPTION("Model size and surface size do not match!"); this->surface.wrap(surface); _gap = allocateGrid(model.getType(), discretization); model.registerField("gap", _gap); } /* -------------------------------------------------------------------------- */ void ContactSolver::printState(UInt iter, Real cost_f, Real error) { // UInt precision = std::cout.precision(); if (iter % dump_frequency) return; Logger().get(LogLevel::info) << std::setw(5) << iter << " " << std::setw(15) << std::scientific << cost_f << " " << std::setw(15) << error << std::endl << std::fixed; } /* -------------------------------------------------------------------------- */ void ContactSolver::addFunctionalTerm(functional::Functional* func) { functional.addFunctionalTerm(func, false); } } // namespace tamaas /* -------------------------------------------------------------------------- */ diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp index 3d48965..1632e05 100644 --- a/src/solvers/kato.cpp +++ b/src/solvers/kato.cpp @@ -1,502 +1,501 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "kato.hh" #include "elastic_functional.hh" -#include "fft_plan_manager.hh" #include "logger.hh" #include "loop.hh" #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { Kato::Kato(Model& model, const GridBase& surface, Real tolerance, Real mu) : ContactSolver(model, surface, tolerance), engine(model.getBEEngine()), mu(mu) { if (model.getType() != model_type::surface_1d && model.getType() != model_type::surface_2d) { TAMAAS_EXCEPTION("Model type is not compatible with Kato solver"); } gap = this->_gap.get(); // locally allocated pressure = &model.getTraction(); N = pressure->getNbPoints(); if (model.getType() == model_type::surface_1d) { initSurfaceWithComponents(); } else { initSurfaceWithComponents(); } engine.registerNeumann(); } /* -------------------------------------------------------------------------- */ Real Kato::solve(GridBase& p0, UInt proj_iter) { if (p0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean pressure does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveTmpl(p0, proj_iter); break; case model_type::surface_2d: cost = solveTmpl(p0, proj_iter); break; default: break; } return cost; } template Real Kato::solveTmpl(GridBase& p0, UInt proj_iter) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers Logger().get(LogLevel::info) << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; pressure->uniformSetComponents(p0); do { computeGradient(); *pressure -= *gap; enforcePressureConstraints(p0, proj_iter); cost = computeCost(); printState(n, cost, cost); } while (cost > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ Real Kato::solveRelaxed(GridBase& g0) { if (g0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean gap does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveRelaxedTmpl(g0); break; case model_type::surface_2d: cost = solveRelaxedTmpl(g0); break; default: break; } return cost; } template Real Kato::solveRelaxedTmpl(GridBase& g0) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers Logger().get(LogLevel::info) << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; *pressure = 0; do { engine.solveNeumann(*pressure, *gap); addUniform(*gap, g0); *gap -= *surfaceComp; *pressure -= *gap; enforcePressureCoulomb(); cost = computeCost(); printState(n, cost, cost); } while (cost > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ Real Kato::solveRegularized(GridBase& p0, Real r) { if (p0.getNbPoints() != pressure->getNbComponents()) { TAMAAS_EXCEPTION( "Target mean pressure does not have the right number of components"); } Real cost = 0; switch (model.getType()) { case model_type::surface_1d: cost = solveRegularizedTmpl(p0, r); break; case model_type::surface_2d: cost = solveRegularizedTmpl(p0, r); break; default: break; } return cost; } template Real Kato::solveRegularizedTmpl(GridBase& p0, Real r) { constexpr UInt comp = model_type_traits::components; Real cost = 0; UInt n = 0; // Printing column headers Logger().get(LogLevel::info) << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; pressure->uniformSetComponents(p0); do { // enforcePressureMean(p0); engine.solveNeumann(*pressure, *gap); *gap -= *surfaceComp; // Impose zero tangential displacement in non-sliding zone UInt count_static = 0; Vector g_static = Loop::reduce( [&] CUDA_LAMBDA(VectorProxy g, VectorProxy p) -> Vector { VectorProxy p_T(p(0)); Real p_T_norm = p_T.l2norm(); Real p_N = p(comp - 1); if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact count_static++; return g; // to compute mean of g_T } else { return 0; } }, range>(*gap), range>(*pressure)); g_static /= count_static != 0 ? count_static : 1; g_static(comp - 1) = 0; Loop::loop( [this, r, g_static] CUDA_LAMBDA(VectorProxy p, VectorProxy g) { // Add frictional term to gradient of functional g -= g_static; // Vector _g = g; // copy VectorProxy g_T(g(0)); VectorProxy g_N(g(comp - 1)); Real g_T_norm = g_T.l2norm(); // g_N += mu * regularize(g_T_norm, r) * g_T_norm; g_N += mu * g_T_norm; // Update pressure with gradient // _g *= 0.1; p -= g; // Truncate negative normal pressure VectorProxy p_T(p(0)); VectorProxy p_N(p(comp - 1)); if (p_N(0) < 0) p_N = 0; // Set tangential pressure p_T = g_T; if (g_T_norm != 0) p_T *= -mu * p_N(0) * regularize(g_T_norm, r) / g_T_norm; }, range>(*pressure), range>(*gap)); // enforcePressureMean(p0); // enforcePressureCoulomb(); enforcePressureConstraints(p0, 50); cost = computeCost(); printState(n, cost, cost); } while (std::abs(cost) > this->tolerance && n++ < this->max_iterations); computeFinalGap(); return cost; } /* -------------------------------------------------------------------------- */ template void Kato::initSurfaceWithComponents() { constexpr UInt comp = model_type_traits::components; surfaceComp = allocateGrid(type, model.getDiscretization(), comp); *surfaceComp = 0; Loop::loop([] CUDA_LAMBDA(Real & s, VectorProxy sc) { sc(comp - 1) = s; }, surface, range>(*surfaceComp)); } /* -------------------------------------------------------------------------- */ template void Kato::enforcePressureMean(GridBase& p0) { Vector corr = computeMean(*pressure); VectorProxy _p0(p0(0)); corr -= _p0; *pressure -= corr; } /* -------------------------------------------------------------------------- */ // template // void Kato::enforcePressureCoulomb() { // Loop::stridedLoop( // [this] CUDA_LAMBDA(VectorProxy&& p) { // VectorProxy p_T(p(0)); // Real p_N = p(comp - 1); // Real p_T_sqrd= p_T.l2squared(); // // Projection normale au cône de friction // bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N); // bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu); // if (cond2) { // p_T = 0; // p(comp - 1) = 0; // } else if (!cond1) { // Real p_T_norm = std::sqrt(p_T_sqrd); // Real k = (p_N + mu * p_T_norm) / (1 + mu * mu); // p_T *= k * mu / p_T_norm; // p(comp - 1) = k; // } // }, // *pressure); // } /* -------------------------------------------------------------------------- */ template void Kato::enforcePressureTresca() { Loop::loop( [this] CUDA_LAMBDA(VectorProxy p) { VectorProxy p_T(p(0)); Real p_N = p(comp - 1); Real p_T_norm = p_T.l2norm(); if (p_N < 0) { p_T = 0; p(comp - 1) = 0; } else if (p_T_norm > mu) { // TODO: replace name of variable mu p_T *= mu / p_T_norm; // p(comp - 1) += p_T_norm - mu; } }, range>(*pressure)); } template void Kato::enforcePressureTresca<2>(); template void Kato::enforcePressureTresca<3>(); /* -------------------------------------------------------------------------- */ /** * Compute mean of the field taking each component separately. */ // template // Vector Kato::computeMean(GridBase& field) { // Vector mean = Loop::stridedReduce( // [] CUDA_LAMBDA(VectorProxy&& f) -> Vector { // return f; // }, // field); // mean /= N; // return mean; // } /* -------------------------------------------------------------------------- */ // template // void Kato::addUniform(GridBase& field, GridBase& vec) { // VectorProxy _vec(vec(0)); // field += _vec; // } /* -------------------------------------------------------------------------- */ Real Kato::computeCost(bool use_tresca) { UInt N = pressure->getNbPoints(); Grid lambda({N}, 1); Grid eta({N}, 1); Grid p_N({N}, 1); Grid p_C({N}, 1); switch (model.getType()) { case model_type::surface_1d: if (!use_tresca) { computeValuesForCost(lambda, eta, p_N, p_C); } else { computeValuesForCostTresca(lambda, eta, p_N, p_C); } break; case model_type::surface_2d: if (!use_tresca) { computeValuesForCost(lambda, eta, p_N, p_C); } else { computeValuesForCostTresca(lambda, eta, p_N, p_C); } break; default: break; } return p_N.dot(lambda) + p_C.dot(eta); } /* -------------------------------------------------------------------------- */ template void Kato::computeValuesForCost(GridBase& lambda, GridBase& eta, GridBase& p_N, GridBase& p_C) { constexpr UInt comp = model_type_traits::components; Real g_N_min = Loop::reduce( [] CUDA_LAMBDA(VectorProxy g) { return g(comp - 1); }, range>(*gap)); Loop::loop( [this, g_N_min] CUDA_LAMBDA(VectorProxy p, VectorProxy g, Real & lambda_, Real & eta_, Real & p_N_, Real & p_C_) { VectorProxy g_T(g(0)); Real g_N = g(comp - 1); Real g_T_norm = g_T.l2norm(); lambda_ = g_N - g_N_min; eta_ = g_T_norm; VectorProxy p_T(p(0)); Real p_N = p(comp - 1); Real p_T_norm = p_T.l2norm(); p_N_ = p_N; p_C_ = (p_N > 0) ? mu * p_N - p_T_norm : 0; }, range>(*pressure), range>(*gap), lambda, eta, p_N, p_C); } /* -------------------------------------------------------------------------- */ template void Kato::computeValuesForCostTresca(GridBase& lambda, GridBase& eta, GridBase& p_N, GridBase& p_C) { constexpr UInt comp = model_type_traits::components; Real g_N_min = Loop::reduce( [] CUDA_LAMBDA(VectorProxy g) { return g(comp - 1); }, range>(*gap)); Real p_C_min = Loop::reduce( [this] CUDA_LAMBDA(VectorProxy p) { VectorProxy p_T(p(0)); Real p_T_norm = p_T.l2norm(); return mu - p_T_norm; }, range>(*pressure)); Loop::loop( [this, g_N_min, p_C_min] CUDA_LAMBDA( VectorProxy p, VectorProxy g, Real & lambda_, Real & eta_, Real & p_N_, Real & p_C_) { VectorProxy g_T(g(0)); Real g_N = g(comp - 1); Real g_T_norm = g_T.l2norm(); lambda_ = g_N - g_N_min; eta_ = g_T_norm; VectorProxy p_T(p(0)); Real p_N = p(comp - 1); Real p_T_norm = p_T.l2norm(); p_N_ = p_N; p_C_ = (p_N > 0) ? mu - p_T_norm - p_C_min : 0; }, range>(*pressure), range>(*gap), lambda, eta, p_N, p_C); } /* -------------------------------------------------------------------------- */ template void Kato::computeFinalGap() { engine.solveNeumann(*pressure, *gap); *gap -= *surfaceComp; Real g_N_min = Loop::reduce( [] CUDA_LAMBDA(VectorProxy g) { return g(comp - 1); }, range>(*gap)); Grid g_shift({comp}, 1); g_shift = 0; g_shift(comp - 1) = -g_N_min; *gap += *surfaceComp; addUniform(*gap, g_shift); model.getDisplacement() = *gap; } /* -------------------------------------------------------------------------- */ Real Kato::regularize(Real x, Real r) { Real xr = x / r; return xr / (1 + std::abs(xr)); } } // namespace tamaas /* -------------------------------------------------------------------------- */ diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp index 03e3e93..6a9346c 100644 --- a/src/solvers/polonsky_keer_rey.cpp +++ b/src/solvers/polonsky_keer_rey.cpp @@ -1,308 +1,307 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "polonsky_keer_rey.hh" #include "elastic_functional.hh" -#include "fft_plan_manager.hh" #include "logger.hh" #include "loop.hh" #include "model_type.hh" #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase& surface, Real tolerance, type variable_type, type constraint_type) : ContactSolver(model, surface, tolerance), variable_type(variable_type), constraint_type(constraint_type) { #define SET_VIEWS(_, __, type) \ case type: { \ setViews(); \ break; \ } switch (model.getType()) { BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES); } #undef SET_VIEWS search_direction = allocateGrid( operation_type, model.getBoundaryDiscretization()); projected_search_direction = allocateGrid( operation_type, model.getBoundaryDiscretization()); switch (variable_type) { case gap: { model.getBEEngine().registerDirichlet(); primal = gap_view.get(); dual = pressure_view.get(); this->functional.addFunctionalTerm( new functional::ElasticFunctionalGap(*integral_op, this->surface), true); break; } case pressure: { model.getBEEngine().registerNeumann(); primal = pressure_view.get(); dual = gap_view.get(); this->functional.addFunctionalTerm( new functional::ElasticFunctionalPressure(*integral_op, this->surface), true); break; } } } /* -------------------------------------------------------------------------- */ Real PolonskyKeerRey::solve(std::vector target_v) { // Update in case of a surface change between solves this->surface_stddev = std::sqrt(this->surface.var()); // Printing column headers Logger().get(LogLevel::info) << std::setw(5) << "Iter" << " " << std::setw(15) << "Cost_f" << " " << std::setw(15) << "Error" << '\n' << std::fixed; const Real target = target_v.back(); // Setting uniform value if constraint if (constraint_type == variable_type && std::abs(primal->sum()) <= 0) *primal = target; else if (constraint_type == variable_type) *primal *= target / primal->mean(); else if (constraint_type != variable_type) *primal = this->surface_stddev; // Algorithm variables Real Gold = 1, error = 0; UInt n = 0; bool delta = false; *search_direction = 0; do { // Computing functional gradient functional.computeGradF(*primal, *dual); const Real dbar = meanOnUnsaturated(*dual); // Enforcing dual constraint via gradient if (constraint_type != variable_type) { *dual += 2 * target + dbar; } else { // Centering dual on primal > 0 *dual -= dbar; } // Computing gradient norm const Real G = computeSquaredNorm(*dual); // Updating search direction (conjugate gradient) updateSearchDirection(delta * G / Gold); Gold = G; // Computing critical step const Real tau = computeCriticalStep(target); // Update primal variable delta = updatePrimal(tau); // We should scale to constraint if (constraint_type == variable_type) enforceMeanValue(target); error = computeError(); const Real cost_f = functional.computeF(*primal, *dual); printState(n, cost_f, error); } while (error > this->tolerance and n++ < this->max_iterations); // Final update of dual variable functional.computeGradF(*primal, *dual); enforceAdmissibleState(); return error; } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::enforceAdmissibleState() { /// Make dual admissible const Real dual_min = dual->min(); *dual -= dual_min; // Primal is pressure: need to make sure gap is positive if (variable_type == pressure) { *displacement_view = *dual; *displacement_view += this->surface; } // Primal is gap: need to make sure dual is positive (pressure + adhesion) else { *displacement_view = *primal; *displacement_view += this->surface; integral_op->apply(*displacement_view, *pressure_view); *pressure_view -= dual_min; } } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$ */ Real PolonskyKeerRey::meanOnUnsaturated(const GridBase& field) const { // Sum on unsaturated const Real sum = Loop::reduce( [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; }, *primal, field); // Number of unsaturated points const UInt n_unsat = Loop::reduce( [] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal); return sum / n_unsat; } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$ */ Real PolonskyKeerRey::computeSquaredNorm(const GridBase& field) const { return Loop::reduce( [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; }, *primal, field); } /* -------------------------------------------------------------------------- */ /** * Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i' * t_i} } \f$ */ Real PolonskyKeerRey::computeCriticalStep(Real target) { integral_op->apply(*search_direction, *projected_search_direction); const Real rbar = meanOnUnsaturated(*projected_search_direction); if (constraint_type == variable_type) *projected_search_direction -= rbar; else { *projected_search_direction += 2 * target + rbar; } const Real num = Loop::reduce( [] CUDA_LAMBDA(const Real& p, const Real& q, const Real& t) { return (p > 0) ? q * t : 0; }, *primal, *dual, *search_direction); const Real denum = Loop::reduce( [] CUDA_LAMBDA(const Real& p, const Real& r, const Real& t) { return (p > 0) ? r * t : 0; }, *primal, *projected_search_direction, *search_direction); return num / denum; } /* -------------------------------------------------------------------------- */ /** * Update steps: * 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$ * 2. Truncate all \f$p\f$ negative * 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i = * p_i - \tau q_i\f$ */ bool PolonskyKeerRey::updatePrimal(Real step) { const UInt na_num = Loop::reduce( [step] CUDA_LAMBDA(Real & p, const Real& q, const Real& t) -> UInt { p -= step * t; // Updating primal if (p < 0) p = 0; // Truncating neg values if (p == 0 && q < 0) { // If non-admissible state p -= step * q; return 1; } else return 0; }, *primal, *dual, *search_direction); return na_num == 0; } /* -------------------------------------------------------------------------- */ /** * Error is based on \f$ \sum{p_i q_i} \f$ */ Real PolonskyKeerRey::computeError() { /// Making sure dual respects constraint *dual -= dual->min(); /// Complementarity error const Real error = primal->dot(*dual); if (std::isnan(error)) TAMAAS_EXCEPTION("Encountered NaN in complementarity error: this may be " "caused by a contact area of a single node."); const Real norm = [this]() { if (variable_type == pressure) return std::abs(primal->sum() * this->surface_stddev); else return std::abs(dual->sum() * this->surface_stddev); }(); return std::abs(error) / (norm * primal->getNbPoints()); } /* -------------------------------------------------------------------------- */ /** * Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t} * \f$ */ void PolonskyKeerRey::updateSearchDirection(Real factor) { Loop::loop( [factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) { t = (p > 0) ? q + factor * t : 0; }, *primal, *dual, *search_direction); } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::enforceMeanValue(Real mean) { *primal *= mean / primal->mean(); } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) { functional.addFunctionalTerm(func, false); } } // namespace tamaas /* -------------------------------------------------------------------------- */ diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp index a33b0ce..1ea70bc 100644 --- a/tests/test_fft.cpp +++ b/tests/test_fft.cpp @@ -1,417 +1,231 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "fft_engine.hh" -#include "fft_plan_manager.hh" -#include "fftransform.hh" #include "grid.hh" #include "grid_hermitian.hh" #include "grid_view.hh" #include "test.hh" using namespace tamaas; using fft = fftw::helper; /* -------------------------------------------------------------------------- */ template struct span { T* ptr; std::size_t size; ~span() { fftw::free(ptr); } const T* begin() const { return ptr; } const T* end() const { return ptr + size; } T* begin() { return ptr; } T* end() { return ptr + size; } operator T*() { return ptr; } }; -TEST(TestFFTInterface, FFT1D) { - constexpr UInt size = 100000; - - span data{fft::alloc_real(size), size}; - span solution{fft::alloc_complex(size / 2 + 1), size / 2 + 1}; - - std::iota(data.begin(), data.end(), 0); - fftw::plan solution_plan{ - fftw::plan_1d_forward(size, data, solution, FFTW_ESTIMATE)}; - fftw::execute(solution_plan); - - Grid grid({size}, 1); - std::iota(grid.begin(), grid.end(), 0); - GridHermitian result({size / 2 + 1}, 1); - - FFTPlanManager::get().createPlan(grid, result).forwardTransform(); -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - ASSERT_TRUE(compare(result, solution, AreComplexEqual())) - << "1D FFTW transform failed"; - FFTPlanManager::get().destroyPlan(grid, result); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFT1D) { constexpr UInt size = 1000; FFTEngine engine; span data{fft::alloc_real(size), size}; span solution{fft::alloc_complex(size / 2 + 1), size / 2 + 1}; fftw::plan solution_plan{ fftw::plan_1d_forward(size, data, solution, engine.flags())}; std::iota(data.begin(), data.end(), 0); fftw::execute(solution_plan); Grid grid({size}, 1); GridHermitian result({size / 2 + 1}, 1); std::iota(grid.begin(), grid.end(), 0); engine.forward(grid, result); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "1D FFTW transform failed"; } -/* -------------------------------------------------------------------------- */ -TEST(TestFFTInterface, FFT2D) { - constexpr UInt size = 100; - constexpr UInt rsize = size * size; - constexpr UInt csize = size * (size / 2 + 1); - - span data{fft::alloc_real(rsize), rsize}; - span solution{fft::alloc_complex(csize), csize}; - - fftw::plan solution_plan{ - fftw::plan_2d_forward(size, size, data, solution, FFTW_ESTIMATE)}; - - std::iota(data.begin(), data.end(), 0); - fftw::execute(solution_plan); - - Grid grid({size, size}, 1); - GridHermitian result({size, size / 2 + 1}, 1); - - std::iota(grid.begin(), grid.end(), 0); - FFTPlanManager::get().createPlan(grid, result).forwardTransform(); - -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - ASSERT_TRUE(compare(result, solution, AreComplexEqual())) - << "2D FFTW transform failed"; - FFTPlanManager::get().destroyPlan(grid, result); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFT2D) { constexpr UInt size = 100; constexpr UInt rsize = size * size; constexpr UInt csize = size * (size / 2 + 1); FFTEngine engine; span data{fft::alloc_real(rsize), rsize}; span solution{fft::alloc_complex(csize), csize}; fftw::plan solution_plan{ fftw::plan_2d_forward(size, size, data, solution, engine.flags())}; std::iota(data.begin(), data.end(), 0); fftw::execute(solution_plan); Grid grid({size, size}, 1); GridHermitian result({size, size / 2 + 1}, 1); std::iota(grid.begin(), grid.end(), 0); engine.forward(grid, result); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "2D FFTW transform failed"; } -/* -------------------------------------------------------------------------- */ - -TEST(TestFFTInterface, FFT1D2Comp) { - constexpr UInt size = 20; - - /// 1D single component FFT should be working here - Grid grid({size}, 2), data({size}, 1); - std::iota(grid.begin(), grid.end(), 0); - std::iota(data.begin(), data.end(), 0); - GridHermitian result({size / 2 + 1}, 2), solution({size / 2 + 1}, 1); - - FFTPlanManager::get().createPlan(grid, result).forwardTransform(); -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - auto& plan = FFTPlanManager::get().createPlan(data, solution); - std::iota(data.begin(), data.end(), 0); - data *= 2; - plan.forwardTransform(); - - ASSERT_TRUE( - compare(make_component_view(result, 0), solution, AreComplexEqual())) - << "1D FFTW transform with 2 components failed on 1st component"; - - data += 1; - plan.forwardTransform(); - - ASSERT_TRUE( - compare(make_component_view(result, 1), solution, AreComplexEqual())) - << "1D FFTW transform with 2 components failed on 2nd component"; - - FFTPlanManager::get().destroyPlan(grid, result); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFT1D2Comp) { constexpr UInt size = 20; /// 1D single component FFT should be working here Grid grid({size}, 2), data({size}, 1); std::iota(grid.begin(), grid.end(), 0); std::iota(data.begin(), data.end(), 0); GridHermitian result({size / 2 + 1}, 2), solution({size / 2 + 1}, 1); FFTEngine engine; engine.forward(grid, result); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif std::iota(data.begin(), data.end(), 0); data *= 2; engine.forward(data, solution); ASSERT_TRUE( compare(make_component_view(result, 0), solution, AreComplexEqual())) << "1D FFTW transform with 2 components failed on 1st component"; data += 1; engine.forward(data, solution); ASSERT_TRUE( compare(make_component_view(result, 1), solution, AreComplexEqual())) << "1D FFTW transform with 2 components failed on 2nd component"; } -/* -------------------------------------------------------------------------- */ -TEST(TestFFTInterface, FFT2D3Comp) { - constexpr UInt size = 20; - - /// 2D single component FFT should be working here - Grid grid({size, size}, 3), data({size, size}, 1); - std::iota(grid.begin(), grid.end(), 0); - std::iota(data.begin(), data.end(), 0); - data *= 3; - - GridHermitian result({size, size / 2 + 1}, 3), - solution({size, size / 2 + 1}, 1); - - FFTPlanManager::get().createPlan(grid, result).forwardTransform(); - auto& plan = FFTPlanManager::get().createPlan(data, solution); - -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - for (UInt i = 0; i < 3; ++i) { - plan.forwardTransform(); - ASSERT_TRUE( - compare(make_component_view(result, i), solution, AreComplexEqual())) - << "2D FFTW transform with 3 components failed on " << i - << "th component"; - data += 1; - } - - FFTPlanManager::get().destroyPlan(grid, result); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFT2D3Comp) { constexpr UInt size = 20; /// 2D single component FFT should be working here Grid grid({size, size}, 3), data({size, size}, 1); std::iota(grid.begin(), grid.end(), 0); std::iota(data.begin(), data.end(), 0); data *= 3; GridHermitian result({size, size / 2 + 1}, 3), solution({size, size / 2 + 1}, 1); FFTEngine engine; engine.forward(grid, result); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif for (UInt i = 0; i < 3; ++i) { engine.forward(data, solution); ASSERT_TRUE( compare(make_component_view(result, i), solution, AreComplexEqual())) << "2D FFTW transform with 3 components failed on " << i << "th component"; data += 1; } } -/* -------------------------------------------------------------------------- */ -TEST(TestFFTInterface, FFT2DViewTransform) { - constexpr UInt size = 20; - - Grid data({size, size}, 1); - GridHermitian solution({size, size / 2 + 1}, 1); - std::iota(std::begin(data), std::end(data), 0); - FFTPlanManager::get().createPlan(data, solution).forwardTransform(); - - Grid grid({size, size}, 3); - auto view = make_component_view(grid, 1); - std::iota(view.begin(), view.end(), 0); - - GridHermitian result({size, size / 2 + 1}, 1); - FFTPlanManager::get().createPlan(view, result).forwardTransform(); - - ASSERT_TRUE(compare(result, solution, AreComplexEqual())) - << "Fourier transform on component view fail"; -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFT2DViewTransform) { constexpr UInt size = 20; Grid data({size, size}, 1); GridHermitian solution({size, size / 2 + 1}, 1); std::iota(std::begin(data), std::end(data), 0); FFTEngine engine; engine.forward(data, solution); Grid grid({size, size}, 3); auto view = make_component_view(grid, 1); std::iota(view.begin(), view.end(), 0); GridHermitian result({size, size / 2 + 1}, 1); engine.forward(view, result); ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "Fourier transform on component view fail"; } -/* -------------------------------------------------------------------------- */ -TEST(TestFFTInterface, FFTI1D2Comp) { - constexpr UInt size = 20; - - Grid grid({size}, 2); - std::iota(grid.begin(), grid.end(), 0); - GridHermitian grid_hermitian({size / 2 + 1}, 2); - Grid result({size}, 2); - - FFTPlanManager::get().createPlan(grid, grid_hermitian).forwardTransform(); - FFTPlanManager::get().createPlan(result, grid_hermitian).backwardTransform(); -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - ASSERT_TRUE(compare(grid, result, AreFloatEqual())) - << "1D FFTI transform with 2 components failed"; - FFTPlanManager::get().destroyPlan(grid, grid_hermitian); - FFTPlanManager::get().destroyPlan(result, grid_hermitian); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFTI1D2Comp) { constexpr UInt size = 20; Grid grid({size}, 2); std::iota(grid.begin(), grid.end(), 0); GridHermitian grid_hermitian({size / 2 + 1}, 2); Grid result({size}, 2); FFTEngine engine; engine.forward(grid, grid_hermitian); engine.backward(result, grid_hermitian); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif ASSERT_TRUE(compare(grid, result, AreFloatEqual())) << "1D FFTI transform with 2 components failed"; } -/* -------------------------------------------------------------------------- */ -TEST(TestFFTInterface, FFTI2D3Comp) { - constexpr UInt size = 20; - - Grid grid({size, size}, 3); - std::iota(grid.begin(), grid.end(), 0); - GridHermitian grid_hermitian({size, size / 2 + 1}, 3); - Grid result({size, size}, 3); - - FFTPlanManager::get().createPlan(grid, grid_hermitian).forwardTransform(); - FFTPlanManager::get().createPlan(result, grid_hermitian).backwardTransform(); -#ifdef USE_CUDA - cudaDeviceSynchronize(); -#endif - - ASSERT_TRUE(compare(grid, result, AreFloatEqual())) - << "2D FFTI transform with 3 components failed"; - FFTPlanManager::get().destroyPlan(grid, grid_hermitian); - FFTPlanManager::get().destroyPlan(result, grid_hermitian); -} - /* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, FFTI2D3Comp) { constexpr UInt size = 20; Grid grid({size, size}, 3); std::iota(grid.begin(), grid.end(), 0); GridHermitian grid_hermitian({size, size / 2 + 1}, 3); Grid result({size, size}, 3); FFTEngine engine; engine.forward(grid, grid_hermitian); engine.backward(result, grid_hermitian); #ifdef USE_CUDA cudaDeviceSynchronize(); #endif ASSERT_TRUE(compare(grid, result, AreFloatEqual())) << "2D FFTI transform with 3 components failed"; } diff --git a/tests/test_fftfreq.cpp b/tests/test_fftfreq.cpp index 7cd23d6..2749695 100644 --- a/tests/test_fftfreq.cpp +++ b/tests/test_fftfreq.cpp @@ -1,117 +1,76 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "fft_engine.hh" -#include "fftransform.hh" #include "grid.hh" #include "grid_hermitian.hh" #include "grid_view.hh" #include "test.hh" #include #include #include using namespace tamaas; namespace py = pybind11; /* -------------------------------------------------------------------------- */ - -TEST(TestFFTInterface, Frequencies1D) { - std::array sizes = {{10}}; - auto freq = FFTransform::computeFrequencies(sizes); - - py::module fftfreq = py::module::import("fftfreq"); - - std::vector reference(freq.dataSize()); - py::array py_arg(reference.size(), reference.data(), py::none()); - - fftfreq.attr("frequencies1D")(py_arg); - - ASSERT_TRUE(compare(reference, freq, AreFloatEqual())) - << "Non hermitian frequencies are wrong"; - - auto hfreq = FFTransform::computeFrequencies(sizes); - std::iota(reference.begin(), reference.end(), 0); - - ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual())) - << "Hermitian frequencies are wrong"; -} - -TEST(TestFFTInterface, Frequencies2D) { - std::array sizes = {{10, 10}}; - auto freq = FFTransform::computeFrequencies(sizes); - - py::module fftfreq = py::module::import("fftfreq"); - - std::vector reference(freq.dataSize()); - py::array py_arg({10, 10, 2}, reference.data(), py::none()); - - fftfreq.attr("frequencies2D")(py_arg); - ASSERT_TRUE(compare(reference, freq, AreFloatEqual())) - << "Non hermitian frequencies are wrong"; - - auto hfreq = FFTransform::computeFrequencies(sizes); - fftfreq.attr("hfrequencies2D")(py_arg); - ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual())) - << "Hermitian frequencies are wrong"; -} - TEST(TestFFTEngine, Frequencies1D) { std::array sizes = {{10}}; auto freq = FFTEngine::computeFrequencies(sizes); py::module fftfreq = py::module::import("fftfreq"); std::vector reference(freq.dataSize()); py::array py_arg(reference.size(), reference.data(), py::none()); fftfreq.attr("frequencies1D")(py_arg); ASSERT_TRUE(compare(reference, freq, AreFloatEqual())) << "Non hermitian frequencies are wrong"; auto hfreq = FFTEngine::computeFrequencies(sizes); std::iota(reference.begin(), reference.end(), 0); ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual())) << "Hermitian frequencies are wrong"; } +/* -------------------------------------------------------------------------- */ TEST(TestFFTEngine, Frequencies2D) { std::array sizes = {{10, 10}}; auto freq = FFTEngine::computeFrequencies(sizes); py::module fftfreq = py::module::import("fftfreq"); std::vector reference(freq.dataSize()); py::array py_arg({10, 10, 2}, reference.data(), py::none()); fftfreq.attr("frequencies2D")(py_arg); ASSERT_TRUE(compare(reference, freq, AreFloatEqual())) << "Non hermitian frequencies are wrong"; auto hfreq = FFTEngine::computeFrequencies(sizes); fftfreq.attr("hfrequencies2D")(py_arg); ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual())) << "Hermitian frequencies are wrong"; }