diff --git a/src/solvers/kato_saturated.cpp b/src/solvers/kato_saturated.cpp
index 2503943..20b451b 100644
--- a/src/solvers/kato_saturated.cpp
+++ b/src/solvers/kato_saturated.cpp
@@ -1,240 +1,239 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2024 Lucas Frérot
  *
  *  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 <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "kato_saturated.hh"
 #include "logger.hh"
 #include <iomanip>
 #include <limits>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 KatoSaturated::KatoSaturated(Model& model, const GridBase<Real>& surface,
                              Real tolerance, Real pmax)
     : PolonskyKeerRey(model, surface, tolerance, PolonskyKeerRey::pressure,
                       PolonskyKeerRey::pressure),
       pmax(pmax) {
   model.request<true, Real>("KatoSaturated::residual_displacement",
                             model.getType(), model.getBoundaryDiscretization(),
                             1);
 }
 
 /* -------------------------------------------------------------------------- */
 Real KatoSaturated::solve(std::vector<Real> load) {
   GridBase<Real> initial_surface = surface;
   GridBase<Real>& residual_disp =
       model.field<Real>("KatoSaturated::residual_displacement");
   residual_disp = 0;
 
   const auto norm = surface.var();
 
   auto negative = [this, norm](const GridBase<Real>& f) {
     const auto neg_norm = Loop::reduce<operation::plus>(
         [](const Real& f) { return (f < 0) * f * f; }, f);
     return neg_norm / (norm * f.getGlobalNbPoints()) > this->tolerance;
   };
 
   UInt n = 0;
 
   do {
     surface = initial_surface;
     surface -= residual_disp;
 
     PolonskyKeerRey::solve(load);
 
     // Update the rough surface
     Loop::loop([] CUDA_LAMBDA(Real & h_pl, Real g) { h_pl -= g * (g < 0); },
                residual_disp, *this->dual);
   } while (negative(*this->dual) and n++ < this->max_iterations);
 
   surface = initial_surface;
   *this->displacement_view += residual_disp;
 
   return Real(n >= this->max_iterations);
 }
 
 /* -------------------------------------------------------------------------- */
 Real KatoSaturated::meanOnUnsaturated(const GridBase<Real>& /*field*/) const {
   return 0;
 }
 
 Real KatoSaturated::computeSquaredNorm(const GridBase<Real>& /*field*/) const {
   return 1.;
 }
 
 void KatoSaturated::updateSearchDirection(Real /*factor*/) {
   *this->search_direction = *this->dual;
 }
 
 Real KatoSaturated::computeCriticalStep(Real /*target*/) {
   // integral_op->apply(*search_direction, *projected_search_direction);
   // Real num = search_direction->dot(*search_direction);
   // Real denum = projected_search_direction->dot(*search_direction);
   // return 0.1 * num / denum;
   return 1;
 }
 
 bool KatoSaturated::updatePrimal(Real step) {
   UInt na_num = Loop::reduce<operation::plus>(
       [step] CUDA_LAMBDA(Real & p, const Real& /*q*/, const Real& t) -> UInt {
         p -= step * t;  // Updating primal
         return 0;
         // 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;
 }
 
 /* -------------------------------------------------------------------------- */
 Real KatoSaturated::computeError() const {
   // We shift the gap by the minimum on unsaturated area
   const auto pmax = this->pmax;
   const Real shift = Loop::reduce<operation::min>(
       [pmax] CUDA_LAMBDA(Real p, Real d) {
         return (p < pmax) ? d : std::numeric_limits<Real>::max();
       },
       *primal, *dual);
 
   // Ignore points that are saturated
   const Real error = Loop::reduce<operation::plus>(
       [pmax, shift] CUDA_LAMBDA(Real p, Real d) {
-        return (p < pmax) ? p * (d - shift) : 0;
+        return (p < pmax) * p * (d - shift);
       },
       *primal, *dual);
 
   if (std::isnan(error))
     throw nan_error{
         TAMAAS_MSG("Encountered NaN in complementarity error: this may be ",
                    "caused by a contact area of a single node.")};
 
   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->getGlobalNbPoints();
   return std::abs(error) / norm;
 }
 
 /* -------------------------------------------------------------------------- */
 void KatoSaturated::enforceMeanValue(Real mean) {
   // We want to cancel the difference between saturated alpha + t and the
   // applied pressure
   const auto pmax = this->pmax;
   *primal -= primal->mean();
   auto f = [&](Real scale) {
     Real sum = Loop::reduce<operation::plus>(
         [pmax, scale] CUDA_LAMBDA(Real t) -> Real {
           t += scale;
 
           if (t > pmax)
             return pmax;
           if (t < 0)
             return 0;
           return t;
         },
         *this->primal);
     sum /= this->primal->getGlobalNbPoints();
     return sum - mean;
   };
 
   if (pmax < mean)
     throw std::runtime_error{TAMAAS_MSG("cannot find equilibrium")};
 
   // Dichotomy + Secant Newton on f
 
   // Initial points
   Real x_n_2 = -primal->max(), x_n_1 = -primal->min(), x_n = 0.;
   Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
 
   // Find a not-too-large search interval
   Logger().get(LogLevel::debug) << TAMAAS_MSG("Searching equilibrium interval");
   while (std::signbit(f(x_n_1)) == std::signbit(f(x_n_2))) {
     x_n_1 *= 10;
     x_n_2 *= 10;
   }
 
   Logger().get(LogLevel::debug) << TAMAAS_MSG(
       "Reducing interval [abs(x1/x2) = ", std::abs(x_n_1 / x_n_2), ']');
   UInt n_dic = 10;
   f_n_1 = f(x_n_1);
   for (UInt i = 0; i < n_dic; ++i) {
     x_n = 0.5 * (x_n_1 + x_n_2);
     f_n = f(x_n);
     if (std::signbit(f_n) == std::signbit(f_n_1)) {
       x_n_1 = x_n;
       f_n_1 = f_n;
     } else
       x_n_2 = x_n;
   }
 
   Logger().get(LogLevel::debug) << TAMAAS_MSG(
       "Starting Newton secant [abs(x1/x2) = ", std::abs(x_n_1 / x_n_2), ']');
 
   // Secant loop
   do {
     f_n_2 = f(x_n_2);
     f_n_1 = f(x_n_1);
     if (f_n_1 == f_n_2)
       break;  // Avoid nans
     x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
     f_n = f(x_n);
     x_n_2 = x_n_1;
     x_n_1 = x_n;
   } while (std::abs(f_n / mean) > 1e-14);
 
   // Pressure update
   Loop::loop(
       [pmax, x_n] CUDA_LAMBDA(Real & t) {
         t += x_n;
         if (t > pmax)
           t = pmax;
         else if (t < 0)
           t = 0.;
       },
       *this->primal);
 }
 
 /* -------------------------------------------------------------------------- */
 void KatoSaturated::enforceAdmissibleState() {
   /// Make dual admissible
   const auto pmax = this->pmax;
   Real shift = Loop::reduce<operation::min>(
       [pmax] CUDA_LAMBDA(Real p, Real d) {
         return (p < pmax) ? d : std::numeric_limits<Real>::max();
       },
       *primal, *dual);
   *dual -= shift;
 
   *displacement_view = *dual;
   *displacement_view += this->surface;
 }
 
 }  // namespace tamaas
diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp
index 85a4352..7e8e3ee 100644
--- a/src/solvers/polonsky_keer_rey.cpp
+++ b/src/solvers/polonsky_keer_rey.cpp
@@ -1,332 +1,332 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2024 Lucas Frérot
  *
  *  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 <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "polonsky_keer_rey.hh"
 #include "elastic_functional.hh"
 #include "logger.hh"
 #include "loop.hh"
 #include "model_type.hh"
 #include <iomanip>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase<Real>& surface,
                                  Real tolerance, type variable_type,
                                  type constraint_type)
     : ContactSolver(model, surface, tolerance), variable_type(variable_type),
       constraint_type(constraint_type) {
   model_type_dispatch(
       [this](auto&& _) {
         constexpr auto type = std::decay_t<decltype(_)>::value;
         setViews<type>();
       },
       model.getType());
 
   search_direction = allocateGrid<true, Real>(
       operation_type, model.getBoundaryDiscretization());
 
   projected_search_direction = allocateGrid<true, Real>(
       operation_type, model.getBoundaryDiscretization());
 
   setupFunctional();
 }
 
 /* -------------------------------------------------------------------------- */
 void PolonskyKeerRey::setIntegralOperator(std::string name) {
   auto op = model.getIntegralOperator(name);
   Logger().get(LogLevel::debug) << TAMAAS_MSG("setting operator ", name);
 
   if ((variable_type == gap and op->getKind() != IntegralOperator::dirichlet) or
       (variable_type == pressure and
        op->getKind() != IntegralOperator::neumann))
     throw model_type_error{TAMAAS_MSG("operator kind ", op->getKind(),
                                       " is not compatible with variable type ",
                                       variable_type)};
 
   if (op->getType() != model_type::basic_1d and
       op->getType() != model_type::basic_2d)
     throw model_type_error{TAMAAS_MSG("operator type ", op->getType(),
                                       " should be a 'basic' model type")};
 
   integral_op = std::move(op);
   setupFunctional();
 }
 
 /* --------------------------------------------------------------------------
  */
 void PolonskyKeerRey::setupFunctional() {
   Logger().get(LogLevel::debug) << TAMAAS_MSG("functional setup");
   this->functional.clear();
 
   if (not integral_op) {
     model_type_dispatch(
         [this](auto&& _) {
           constexpr auto type = std::decay_t<decltype(_)>::value;
           defaultOperator<type>();
         },
         model.getType());
   }
 
   switch (variable_type) {
   case gap: {
     model.getBEEngine().registerDirichlet();
     primal = gap_view.get();
     dual = pressure_view.get();
     this->functional.addFunctionalTerm(
         std::make_shared<functional::ElasticFunctionalGap>(*integral_op,
                                                            this->surface));
     break;
   }
   case pressure: {
     model.getBEEngine().registerNeumann();
     primal = pressure_view.get();
     dual = gap_view.get();
     this->functional.addFunctionalTerm(
         std::make_shared<functional::ElasticFunctionalPressure>(*integral_op,
                                                                 this->surface));
     break;
   }
   }
 }
 
 /* --------------------------------------------------------------------------
  */
 Real PolonskyKeerRey::solve(std::vector<Real> 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);
 
   // Log final iteration
   logIteration(n, functional.computeF(*primal, *dual), error);
 
   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<Real>& field) const {
   // Sum on unsaturated
   const auto red = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(const Real& p, const Real& f) {
         return Vector<Real, 2>({(p > 0) * f, static_cast<Real>(p > 0)});
       },
       *primal, field);
   return red(0) / red(1);
 }
 
 /* --------------------------------------------------------------------------
  */
 
 /**
  * Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
  */
 Real PolonskyKeerRey::computeSquaredNorm(const GridBase<Real>& field) const {
   return Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) * f * f; },
       *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);
   const Real shift = [rbar, target, this]() {
     if (constraint_type == variable_type)
       return rbar;
     else
       return -(2 * target + rbar);
   }();
 
   const auto red = Loop::reduce<operation::plus>(
       [shift] CUDA_LAMBDA(const Real& p, const Real& q, const Real& r,
                           const Real& t) {
         return Vector<Real, 2>({(p > 0) * q * t, (p > 0) * (r - shift) * t});
       },
       *primal, *dual, *projected_search_direction, *search_direction);
 
   return red(0) / red(1);
 }
 
 /* --------------------------------------------------------------------------
  */
 
 /**
  * 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<operation::plus>(
       [step] CUDA_LAMBDA(Real & p, const Real& q, const Real& t) -> UInt {
         p -= step * t;  // Updating primal
         p *= p >= 0;    // Truncating neg values
         const auto na = p == 0 and q < 0;
         p -= na * step * q;  // Adding point to active set
         return na;
       },
       *primal, *dual, *search_direction);
 
   return na_num == 0;
 }
 
 /* --------------------------------------------------------------------------
  */
 
 /**
  * Error is based on \f$ \sum{p_i q_i} \f$
  */
 Real PolonskyKeerRey::computeError() const {
   /// Complementarity error
   const Real primal_sum = primal->sum();
   const Real error = primal->dot(*dual) - primal_sum * dual->min();
 
   if (std::isnan(error))
     throw nan_error{
         TAMAAS_MSG("Encountered NaN in complementarity error: this may be "
                    "caused by a contact area of a single node.")};
 
   const Real norm = [this, primal_sum]() {
     if (variable_type == pressure)
       return std::abs(primal_sum * this->surface_stddev);
     return std::abs(dual->sum() * this->surface_stddev);
   }();
 
-  return std::abs(error) / (norm * primal->getGlobalNbPoints());
+  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); },
              *primal, *dual, *search_direction);
 }
 
 /* --------------------------------------------------------------------------
  */
 void PolonskyKeerRey::enforceMeanValue(Real mean) {
   *primal *= mean / primal->mean();
 }
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */