diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp index 36f57c4..a830faf 100644 --- a/src/solvers/polonsky_keer_rey.cpp +++ b/src/solvers/polonsky_keer_rey.cpp @@ -1,303 +1,305 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-19 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) { const Real target = target_v.back(); Real G = 0, Gold = 1, error = 0, error_norm = 1; UInt n = 0; bool delta = false; *search_direction = 0; - //*dual = 0; + + // 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; // 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) enforceMeanValue(target); error = computeError() / error_norm; 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 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 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) { integral_op->apply(*search_direction, *projected_search_direction); 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(const Real& p, const Real& q, const Real& t) { return (p > 0) ? q * t : 0; }, *primal, *dual, *search_direction); 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) { 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(); Real norm = 1; 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."); 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); } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::enforceMeanValue(Real mean) { *primal *= mean / primal->mean(); } /* -------------------------------------------------------------------------- */ void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) { functional.addFunctionalTerm(func, false); } } // namespace tamaas /* -------------------------------------------------------------------------- */