diff --git a/site_scons/site_init.py b/site_scons/site_init.py index 3536579..0c345de 100644 --- a/site_scons/site_init.py +++ b/site_scons/site_init.py @@ -1,123 +1,123 @@ from SCons.Script import * from os.path import * from subprocess import check_output, STDOUT import os def fftw(env): """A Tool to search for fftw headers and libraries""" if env.GetOption('clean'): return env.SetDefault(FFTW_VERSION='3') env.SetDefault(FFTW_LIBRARY_WISH=[]) if 'FFTW_INCLUDE_DIR' in env: env.AppendUnique(CPPPATH=env['FFTW_INCLUDE_DIR']) if 'FFTW_LIBRARY_DIR' in env: env.AppendUnique(LIBPATH=env['FFTW_LIBRARY_DIR']) version = env['FFTW_VERSION'] if version == "2": lib_names = {'main': 'fftw'} inc_names = ['fftw.h'] else: lib_names = {'main': 'fftw3', 'thread': 'fftw3_threads', 'omp': 'fftw3_omp'} inc_names = ['fftw3.h'] try: libs = [lib_names[i] for i in env['FFTW_LIBRARY_WISH']] except: raise SCons.Errors.StopError( 'Incompatible wishlist {0} from version {1}'.format( env['FFTW_LIBRARY_WISH'], env['FFTW_VERSION'])) env.AppendUnique(LIBS=libs) if version == "2": env.Append(LIBS='m') conf = Configure(env) if not conf.CheckLibWithHeader(libs, inc_names, 'c++'): raise SCons.Errors.StopError( 'Failed to find libraries {0} or ' 'headers {1}.'.format(str(lib_names), str(inc_names))) env = conf.Finish() def boost(env): """A tool to confugure for boost""" if env.GetOption('clean'): return if 'BOOST_INCLUDE_DIR' in env: env.AppendUnique(CPPPATH=env['BOOST_INCLUDE_DIR']) conf = Configure(env) if not conf.CheckCXXHeader('boost/preprocessor/seq.hpp'): raise SCons.Errors.StopError( 'Failed to find Boost header') env = conf.Finish() def pybind11(env): """A tool to configure pybind11""" if env.GetOption('clean'): return # Run code below to get versions from user preference and not scons versions_script = """ from __future__ import print_function import distutils.sysconfig as dsys from numpy.distutils.misc_util import get_numpy_include_dirs as numpy_dirs print(dsys.get_python_inc()) for d in numpy_dirs(): print(d)""" includes = check_output([env['py_exec'], "-c", versions_script], universal_newlines=True).split('\n') includes += [Dir('#third-party/pybind11/include')] # Extension of shared library for python extension = check_output([env['py_exec'] + '-config', "--extension-suffix"], universal_newlines=True).split('\n')[0] - env['SHLIBSUFFIX'] = extension + '.so' + env['SHLIBSUFFIX'] = extension env.AppendUnique(CPPPATH=includes) conf = Configure(env) if not conf.CheckCXXHeader('pybind11/pybind11.h'): raise SCons.Errors.StopError( 'Failed to find pybind11 header\n' + "Run 'git submodule update --init --recursive third-party/pybind11'") env = conf.Finish() def gtest(env): """A tool to configure GoogleTest""" if env.GetOption('clean'): return conf = Configure(env) if not conf.CheckCXXHeader('gtest/gtest.h'): raise SCons.Errors.StopError( 'Failed to find GoogleTest header\n' + "Run 'git submodule update --init --recursive third-party/googletest'") env = conf.Finish() def thrust(env): """A tool to confugure for thrust""" if env.GetOption('clean'): return if 'THRUST_INCLUDE_DIR' in env: env.AppendUnique(CPPPATH=env['THRUST_INCLUDE_DIR']) conf = Configure(env) if not conf.CheckCXXHeader('thrust/version.h'): raise SCons.Errors.StopError( 'Failed to find Thrust library') env = conf.Finish() diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp index e1f139d..337c533 100644 --- a/src/solvers/polonsky_keer_rey.cpp +++ b/src/solvers/polonsky_keer_rey.cpp @@ -1,308 +1,312 @@ /** * @file * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2016-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 "polonsky_keer_rey.hh" #include "elastic_functional.hh" #include "loop.hh" #include "model_type.hh" +#include "fft_plan_manager.hh" #include #include /* -------------------------------------------------------------------------- */ __BEGIN_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), functional(model.getBEEngine()) { if (model.getType() != model_type::basic_1d && model.getType() != model_type::basic_2d) { TAMAAS_EXCEPTION( "Model type is not compatible with PolonskyKeerRey 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 PolonskyKeerRey::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 PolonskyKeerRey::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 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) { 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 PolonskyKeerRey::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 PolonskyKeerRey::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 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); } /* -------------------------------------------------------------------------- */ PolonskyKeerRey::~PolonskyKeerRey() { 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 PolonskyKeerRey::addFunctionalTerm(functional::Functional * func) { functional.addFunctionalTerm(func, false); } __END_TAMAAS__ /* -------------------------------------------------------------------------- */