diff --git a/src/SConscript b/src/SConscript index 59cb0b8..b4e6fec 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,114 +1,115 @@ 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 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 """.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/kato.cpp b/src/solvers/kato.cpp new file mode 100644 index 0000000..224650e --- /dev/null +++ b/src/solvers/kato.cpp @@ -0,0 +1,312 @@ +/** + * @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 "kato.hh" +#include "elastic_functional.hh" +#include "loop.hh" +#include "model_type.hh" +#include "fft_plan_manager.hh" +#include +#include +/* -------------------------------------------------------------------------- */ + +__BEGIN_TAMAAS__ + +Kato::Kato(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), functional(model.getBEEngine()) { + if (model.getType() != model_type::basic_1d && + model.getType() != model_type::basic_2d) { + TAMAAS_EXCEPTION( + "Model type is not compatible with Kato solver"); + } + + /// Gap is locally allocated + auto gap_ = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); + auto pressure_ = &model.getTraction(); + + switch (variable_type) { + case gap: { + primal = gap_; + dual = pressure_; + functional.addFunctionalTerm(new functional::ElasticFunctionalGap( + model.getBEEngine(), this->surface), + true); + break; + } + case pressure: { + primal = pressure_; + dual = gap_; + functional.addFunctionalTerm(new functional::ElasticFunctionalPressure( + model.getBEEngine(), this->surface), + true); + break; + } + } + + search_direction = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); + + projected_search_direction = + allocateGrid(model.getType(), model.getDiscretization(), + model.getTraction().getNbComponents()); +} + +/* -------------------------------------------------------------------------- */ + +Real Kato::solve(Real target) { + Real G = 0, Gold = 1, error = 0, error_norm = 1; + UInt n = 0; + bool delta = false; + *search_direction = 0; + //*dual = 0; + + // Printing column headers + std::cout << std::setw(5) << "Iter" + << " " << std::setw(15) << "Cost_f" + << " " << std::setw(15) << "Error" << '\n' + << std::fixed; + + // 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; + + do { + // Computing functional gradient + functional.computeGradF(*primal, *dual); + + 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 + G = computeSquaredNorm(*dual); + + // Updating search direction (conjugate gradient) + updateSearchDirection(delta * G / Gold); + Gold = G; + + // Computing critical step + Real tau = computeCriticalStep(target); + + // Update primal variable + delta = updatePrimal(tau); + + // We should scale to constraint + if (constraint_type == variable_type) + *primal *= target / primal->mean(); + + error = computeError() / error_norm; + Real cost_f = functional.computeF(*primal, *dual); + printState(n, cost_f, error); + } while (error > this->tolerance && n++ < this->max_iterations); + + // Final update of dual variable + functional.computeGradF(*primal, *dual); + + Real dual_min = dual->min(); + *dual -= dual_min; + + // Primal is pressure: need to make sure gap is positive + if (variable_type == pressure) { + model.getDisplacement() = *dual; + model.getDisplacement() += this->surface; + } + // Primal is gap: need to make sure dual is positive (pressure + adhesion) + else { + model.getDisplacement() = *primal; + model.getDisplacement() += this->surface; + this->model.solveDirichlet(); + model.getTraction() -= dual_min; + } + + return error; +} + +/* -------------------------------------------------------------------------- */ + +/** + * Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$ + */ +Real Kato::meanOnUnsaturated(const GridBase& field) const { + // Sum on unsaturated + Real sum = Loop::reduce( + [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; }, + *primal, field); + // Number of unsaturated points + 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 Kato::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 Kato::computeCriticalStep(Real target) { + switch (variable_type) { + case gap: + model.getBEEngine().solveDirichlet(*search_direction, + *projected_search_direction); + break; + case pressure: + model.getBEEngine().solveNeumann(*search_direction, + *projected_search_direction); + break; + } + + Real rbar = meanOnUnsaturated(*projected_search_direction); + if (constraint_type == variable_type) + *projected_search_direction -= rbar; + else { + *projected_search_direction += 2 * target + rbar; + } + + Real num = Loop::reduce( + [] CUDA_LAMBDA(Real & p, Real & q, Real & t) { + return (p > 0) ? q * t : 0; + }, + *primal, *dual, *search_direction); + Real denum = Loop::reduce( + [] CUDA_LAMBDA(Real & p, Real & r, 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 Kato::updatePrimal(Real step) { + UInt na_num = Loop::reduce( + [step] CUDA_LAMBDA(Real & p, Real & q, 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 Kato::computeError() { + /// Making sure dual respects constraint + *dual -= dual->min(); + Real error = primal->dot(*dual); + Real norm = 1; + + if (variable_type == pressure) + norm = std::abs(primal->sum() * this->surface_stddev); + else + norm = std::abs(dual->sum() * this->surface_stddev); + + norm *= primal->getNbPoints(); + return std::abs(error) / norm; +} + +/* -------------------------------------------------------------------------- */ + +/** + * Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t} + * \f$ + */ +void Kato::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); +} + +/* -------------------------------------------------------------------------- */ +Kato::~Kato() { + switch (variable_type) { + case gap: + delete primal; + break; + case pressure: + delete dual; + break; + } + + delete search_direction; + delete projected_search_direction; + FFTPlanManager::get().clean(); // TODO find a smarter way to deal with plans +} + +/* -------------------------------------------------------------------------- */ +void Kato::addFunctionalTerm(functional::Functional * func) { + functional.addFunctionalTerm(func, false); +} + +__END_TAMAAS__ +/* -------------------------------------------------------------------------- */ diff --git a/src/solvers/kato.hh b/src/solvers/kato.hh new file mode 100644 index 0000000..bc7196d --- /dev/null +++ b/src/solvers/kato.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 __KATO__ +#define __KATO__ +/* -------------------------------------------------------------------------- */ +#include "contact_solver.hh" +#include "meta_functional.hh" +#include "tamaas.hh" +/* -------------------------------------------------------------------------- */ +__BEGIN_TAMAAS__ + +class Kato : public ContactSolver { +public: + /// Types of algorithm (primal/dual) or constraint + enum type { gap, pressure }; + +public: + /// Constructor + Kato(Model& model, const GridBase& surface, Real tolerance, + type variable_type, type constraint_type); + /// Destructor + ~Kato() override; + +public: + /// Solve + Real solve(Real target); + /// Mean on unsaturated constraint zone + Real meanOnUnsaturated(const GridBase& field) const; + /// Compute squared norm + Real computeSquaredNorm(const GridBase& field) const; + /// Update search direction + void updateSearchDirection(Real factor); + /// Compute critical step + Real computeCriticalStep(Real target = 0); + /// Update primal and check non-admissible state + bool updatePrimal(Real step); + /// Compute error/stopping criterion + Real computeError(); + /// Add term to functional (fuck swig) + void addFunctionalTerm(::tamaas::functional::Functional * func); + +protected: + type variable_type, constraint_type; + GridBase* primal = nullptr; + GridBase* dual = nullptr; + GridBase* search_direction = nullptr; + GridBase* projected_search_direction = nullptr; + functional::MetaFunctional functional; +}; + +__END_TAMAAS__ + +#endif // __KATO__