diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp index d3c4483..fe78186 100644 --- a/python/wrap/solvers.cpp +++ b/python/wrap/solvers.cpp @@ -1,89 +1,97 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 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 "wrap.hh" #include "contact_solver.hh" #include "polonsky_keer_rey.hh" #include "kato.hh" #include "beck_teboulle.hh" #include "condat.hh" +#include "polonsky_keer_tan.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { /* -------------------------------------------------------------------------- */ using namespace py::literals; /* -------------------------------------------------------------------------- */ void wrapSolvers(py::module& mod) { py::class_(mod, "ContactSolver") .def(py::init&, Real>(), "model"_a, "surface"_a, "tolerance"_a) .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a) .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a) .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm); py::class_ pkr(mod, "PolonskyKeerRey"); pkr.def(py::init&, Real, PolonskyKeerRey::type, PolonskyKeerRey::type>(), "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a, "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &PolonskyKeerRey::solve, "target"_a) .def("computeError", &PolonskyKeerRey::computeError); py::enum_(pkr, "type") .value("gap", PolonskyKeerRey::gap) .value("pressure", PolonskyKeerRey::pressure) .export_values(); py::class_ kato(mod, "Kato"); kato.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a=50) .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a) .def("computeCost", &Kato::computeCost); py::class_ bt(mod, "BeckTeboulle"); bt.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &BeckTeboulle::solve, "p0"_a) .def("computeCost", &BeckTeboulle::computeCost); py::class_ cd(mod, "Condat"); cd.def(py::init&, Real, Real>(), "model"_a, "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) .def("solve", &Condat::solve, "p0"_a, "grad_step"_a=0.9) .def("computeCost", &Condat::computeCost); + + py::class_ pkt(mod, "PolonskyKeerTan"); + pkt.def(py::init&, Real, Real>(), + "model"_a, "surface"_a, "tolerance"_a, "mu"_a, + py::keep_alive<1, 2>(), py::keep_alive<1, 3>()) + .def("solve", &PolonskyKeerTan::solve, "p0"_a) + .def("computeCost", &PolonskyKeerTan::computeCost); } } // namespace wrap __END_TAMAAS__ diff --git a/src/SConscript b/src/SConscript index ef6d985..43fc020 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,117 +1,118 @@ import os Import('main_env') def prepend(path, list): return [os.path.join(path, x) for x in list] env = main_env # Core core_list = """ fft_plan_manager.cpp fftransform.cpp fftransform_fftw.cpp grid.cpp grid_hermitian.cpp statistics.cpp surface.cpp tamaas.cpp legacy_types.cpp loop.cpp """.split() core_list = prepend('core', core_list) core_list.append('tamaas_info.cpp') # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_filter.cpp surface_generator_filter_fft.cpp isopowerlaw.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 westergaard.cpp elastic_functional.cpp meta_functional.cpp adhesion_functional.cpp elasto_plastic_model.cpp ve_engine.cpp kelvin.cpp mindlin.cpp """.split() model_list = prepend('model', model_list) # Solvers solvers_list = """ contact_solver.cpp polonsky_keer_rey.cpp kato.cpp beck_teboulle.cpp condat.cpp +polonsky_keer_tan.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 + bem_list # Adding GPU if needed if env['backend'] == 'cuda': rough_contact_list += gpu_list # Generating dependency files # env.AppendUnique(CXXFLAGS=['-MMD']) # for file_name in rough_contact_list: # obj = file_name.replace('.cpp', '.os') # dep_file_name = file_name.replace('.cpp', '.d') # env.SideEffect(dep_file_name, obj) # env.ParseDepends(dep_file_name) libTamaas = env.SharedLibrary('Tamaas', rough_contact_list) Export('libTamaas') diff --git a/src/solvers/polonsky_keer_tan.cpp b/src/solvers/polonsky_keer_tan.cpp new file mode 100644 index 0000000..5ce344f --- /dev/null +++ b/src/solvers/polonsky_keer_tan.cpp @@ -0,0 +1,323 @@ +/** + * @file + * + * @author Son Pham-Ba + * + * @section LICENSE + * + * Copyright (©) 2016-2018 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 "polonsky_keer_tan.hh" +#include +/* -------------------------------------------------------------------------- */ + +__BEGIN_TAMAAS__ + +PolonskyKeerTan::PolonskyKeerTan(Model& model, const GridBase& surface, + Real tolerance, Real mu) + : Kato(model, surface, tolerance, mu) { + search_direction = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); + + search_direction_backup = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); + + projected_search_direction = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); +} + +/* -------------------------------------------------------------------------- */ + +Real PolonskyKeerTan::solve(GridBase& p0) { + 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); + break; + case model_type::surface_2d: + cost = solveTmpl(p0); + break; + default: + break; + } + + return cost; +} + +template +Real PolonskyKeerTan::solveTmpl(GridBase& p0) { + // Printing column headers + std::cout << std::setw(5) << "Iter" + << " " << std::setw(15) << "Cost_f" + << " " << std::setw(15) << "Error" << '\n' + << std::fixed; + + constexpr UInt comp = model_type_traits::components; + Real cost = 0.0; + UInt n = 0; + + pressure->uniformSetComponents(p0); + Real R_old = 1.0; + *search_direction = 0.0; + + Vector gap_mean; + + do { + enforcePressureMean(p0); + + // Compute functional gradient + engine.solveNeumann(*pressure, *gap); + *gap -= surfaceComp; + computeMean(*gap, gap_mean, true); + *gap -= gap_mean; + + // Compute search direction + Real R = computeSquaredNorm(*gap); + *search_direction *= R / R_old; + *search_direction += *gap; + R_old = R; + truncateSearchDirection(true); + + // Compute step size + Real tau = computeStepSize(true); + + // Update pressure + *search_direction *= tau; + *pressure -= *search_direction; + + // Enforce constraints + enforcePressureCoulomb(); + // *search_direction_backup = *gap; + Real alpha = tau; + while(enforcePressureConstraints(p0, alpha) && + (cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations) { + cost = computeCost(); + printState(n-1, cost, cost); + } + // Kato::enforcePressureConstraints(p0, 50); + // std::cout << "I got out at n = " << n << std::endl; + + cost = computeCost(); + printState(n, cost, cost); + } while ((cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations); + + computeFinalGap(); + + return cost; +} + +/* -------------------------------------------------------------------------- */ + +template +bool PolonskyKeerTan::enforcePressureConstraints( + GridBase& p0, Real& alpha) { + + *search_direction_backup = *gap; + engine.solveNeumann(*pressure, *gap); + *gap -= surfaceComp; + Vector gap_mean; + computeMean(*gap, gap_mean, true); + *gap -= gap_mean; + // computeMean(*search_direction_backup, gap_mean, true); + // *search_direction_backup -= gap_mean; + + UInt na_count = Loop::stridedReduce( + [alpha] CUDA_LAMBDA(VectorProxy&& p, + Vector&& q, + VectorProxy&& g) { + if (p(comp - 1) == 0.0 && g(comp - 1) < 0.0) { + q *= alpha; + p -= q; + return 1; + } else { + return 0; + } + }, + *pressure, *search_direction_backup, *gap); + // std::cout << na_count << std::endl; + if (na_count == 0) return false; + + // Enforce mean value + enforcePressureMean(p0); + + // Compute constraint gradient + computeGradientConstraints(); + + // Compute search direction + *search_direction = *gap; + truncateSearchDirection(true); + + // Compute step size + alpha = computeStepSize(true); + + // Update pressure + *search_direction *= alpha; + *pressure -= *search_direction; + enforcePressureCoulomb(); + + return true; +} + +/* -------------------------------------------------------------------------- */ + +template +void PolonskyKeerTan::enforcePressureMean(GridBase& p0) { + Vector pressure_mean; + computeMean(*pressure, pressure_mean, false); + *pressure -= pressure_mean; + addUniform(*pressure, p0); + + // *pressure /= pressure_mean; // what if == 0 ? + // VectorProxy _p0(p0(0)); + // *pressure *= _p0; + + // Loop::stridedLoop( + // [pressure_mean, p0] CUDA_LAMBDA(VectorProxy&& p) { + // for (UInt i = 0; i < comp - 1; i++) + // p(i) += p0(i) - pressure_mean(i); + // p(comp - 1) *= p0(comp - 1) / pressure_mean(comp - 1); + // }, + // *pressure); +} + +/* -------------------------------------------------------------------------- */ + +template +void PolonskyKeerTan::computeMean(GridBase& field, + Vector& mean, + bool on_c) { + UInt count = 0; + mean = Loop::stridedReduce( + [&] CUDA_LAMBDA(VectorProxy&& f, + VectorProxy&& p) -> Vector { + if ((!on_c) || p(comp - 1) > 0.0) { + count ++; + return f; + } else { + return 0; + } + }, + field, *pressure); + + mean /= count; +} + +/* -------------------------------------------------------------------------- */ + +Real PolonskyKeerTan::computeSquaredNorm(GridBase& field) { + return Loop::reduce( + [] CUDA_LAMBDA(Real& f) { + return f * f; + }, + field); +} + +/* -------------------------------------------------------------------------- */ + +template +void PolonskyKeerTan::truncateSearchDirection(bool on_c) { + if (!on_c) return; + Loop::stridedLoop( + [] CUDA_LAMBDA(VectorProxy&& t, + VectorProxy&& p) { + if (p(comp - 1) == 0.0) + t = 0.0; + }, + *search_direction, *pressure); +} + +/* -------------------------------------------------------------------------- */ + +template +Real PolonskyKeerTan::computeStepSize(bool on_c) { + Vector r_mean; + + engine.solveNeumann(*search_direction, *projected_search_direction); + computeMean(*projected_search_direction, r_mean, on_c); + *projected_search_direction -= r_mean; + + Real num = Loop::stridedReduce( + [] CUDA_LAMBDA(VectorProxy&& q, VectorProxy&& t) { + Real dot = 0.0; + for (UInt i = 0; i < comp; i++) dot += q(i) * t(i); + return dot; + }, + *gap, *search_direction); + + Real denum = Loop::stridedReduce( + [] CUDA_LAMBDA(VectorProxy&& r, VectorProxy&& t) { + Real dot = 0.0; + for (UInt i = 0; i < comp; i++) dot += r(i) * t(i); + return dot; + }, + *projected_search_direction, *search_direction); + + return num / denum; +} + +/* -------------------------------------------------------------------------- */ + +template +void PolonskyKeerTan::computeGradientConstraints() { + engine.solveNeumann(*pressure, *gap); + *gap -= surfaceComp; + + // *search_direction_backup = *gap; + // Vector gap_mean; + // computeMean(*search_direction_backup, gap_mean, true); + // *search_direction_backup -= gap_mean; + + Loop::stridedLoop( + [this] CUDA_LAMBDA(VectorProxy&& g, + VectorProxy&& p) { + VectorProxy g_T(g(0)); + VectorProxy g_N = g(comp - 1); + Real g_T_norm = g_T.l2norm(); + + VectorProxy p_T(p(0)); + Real p_T_norm = p_T.l2norm(); + + // g_N += mu * g_T_norm; + // if (p_T_norm == 0.0) { + // g_T *= -1.0; + // } else { + // g_T = p_T; + // g_T *= - g_T_norm / p_T_norm; + // } + }, + *gap, *pressure); + + Vector gap_mean; + computeMean(*gap, gap_mean, true); + *gap -= gap_mean; +} + +__END_TAMAAS__ +/* -------------------------------------------------------------------------- */ \ No newline at end of file diff --git a/src/solvers/polonsky_keer_tan.hh b/src/solvers/polonsky_keer_tan.hh new file mode 100644 index 0000000..b691640 --- /dev/null +++ b/src/solvers/polonsky_keer_tan.hh @@ -0,0 +1,77 @@ +/** + * @file + * + * @author Son Pham-Ba + * + * @section LICENSE + * + * Copyright (©) 2016-2018 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 __POLONSKY_KEER_TAN_HH__ +#define __POLONSKY_KEER_TAN_HH__ +/* -------------------------------------------------------------------------- */ +#include "kato.hh" +/* -------------------------------------------------------------------------- */ +__BEGIN_TAMAAS__ + +class PolonskyKeerTan : public Kato { +public: + /// Constructor + PolonskyKeerTan(Model& model, const GridBase& surface, Real tolerance, + Real mu); + +public: + /// Solve + Real solve(GridBase& p0); + +private: + /// Template for solve function + template + Real solveTmpl(GridBase& p0); + /// Inner loop to minimize orthogonality + template + bool enforcePressureConstraints(GridBase& p0, Real& alpha); + /// Enforce pressure mean + template + void enforcePressureMean(GridBase& p0); + /// Compute mean of field (only on I_c) + template + void computeMean(GridBase& field, Vector& mean, bool on_c); + /// Compute squared norm + Real computeSquaredNorm(GridBase& field); + /// Restrict search direction on I_c + template + void truncateSearchDirection(bool on_c); + /// Compute optimal step size (only on I_c) + template + Real computeStepSize(bool on_c); + /// Compute gradient for constraints minimization + template + void computeGradientConstraints(); + +private: + std::unique_ptr> search_direction = nullptr; + std::unique_ptr> search_direction_backup = nullptr; + std::unique_ptr> projected_search_direction = nullptr; +}; + +__END_TAMAAS__ + +#endif // __POLONSKY_KEER_TAN_HH__