diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel
index ac1ae8e598..809ae95cbc 100644
--- a/ortools/sat/BUILD.bazel
+++ b/ortools/sat/BUILD.bazel
@@ -981,6 +981,7 @@ cc_library(
         "//ortools/port:proto_utils",
         "//ortools/util:logging",
         "//ortools/util:sorted_interval_list",
+        "@com_google_absl//absl/base:log_severity",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/types:span",
     ],
diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc
index df44d8936b..b11a050b13 100644
--- a/ortools/sat/cp_model_checker.cc
+++ b/ortools/sat/cp_model_checker.cc
@@ -636,6 +636,14 @@ std::string ValidateRoutesConstraint(const ConstraintProto& ct) {
         "All nodes in a route constraint must have incident arcs");
   }
 
+  // If the demands field is present, it must be of size nodes.size().
+  if (!ct.routes().demands().empty() &&
+      ct.routes().demands().size() != nodes.size()) {
+    return absl::StrCat(
+        "If the demands fields is set, it must be of size num_nodes:",
+        nodes.size());
+  }
+
   return ValidateGraphInput(/*is_route=*/true, ct.routes());
 }
 
@@ -793,6 +801,15 @@ std::string ValidateReservoirConstraint(const CpModelProto& model,
   }
   for (const LinearExpressionProto& expr : ct.reservoir().time_exprs()) {
     RETURN_IF_NOT_EMPTY(ValidateAffineExpression(model, expr));
+    // We want to be able to safely put time_exprs[i]-time_exprs[j] in a linear.
+    if (MinOfExpression(model, expr) <=
+            -std::numeric_limits<int64_t>::max() / 4 ||
+        MaxOfExpression(model, expr) >=
+            std::numeric_limits<int64_t>::max() / 4) {
+      return absl::StrCat(
+          "Potential integer overflow on time_expr of a reservoir: ",
+          ProtobufShortDebugString(ct));
+    }
   }
   for (const LinearExpressionProto& expr : ct.reservoir().level_changes()) {
     RETURN_IF_NOT_EMPTY(ValidateConstantAffineExpression(model, expr));
diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc
index c706273acb..f4be6c86a8 100644
--- a/ortools/sat/cp_model_expand.cc
+++ b/ortools/sat/cp_model_expand.cc
@@ -56,7 +56,6 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
                                  int64_t sum_of_negative_demand,
                                  ConstraintProto* reservoir_ct,
                                  PresolveContext* context) {
-  SolutionCrush& crush = context->solution_crush();
   const ReservoirConstraintProto& reservoir = reservoir_ct->reservoir();
   const int num_events = reservoir.time_exprs_size();
 
@@ -72,73 +71,6 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
   std::vector<int> level_vars(num_events);
   for (int i = 0; i < num_events; ++i) {
     level_vars[i] = context->NewIntVar(Domain(var_min, var_max));
-    if (crush.HintIsLoaded()) {
-      // The hint of active events is set later.
-      crush.SetNewVariableHint(level_vars[i], 0);
-    }
-  }
-
-  // The hints of the active events, in the order they should appear in the
-  // circuit. The hints are collected first, and sorted later.
-  struct ReservoirEventHint {
-    int index;  // In the reservoir constraint.
-    int64_t time;
-    int64_t level_change;
-  };
-  std::vector<ReservoirEventHint> active_event_hints;
-  bool has_complete_hint = false;
-  if (crush.HintIsLoaded()) {
-    has_complete_hint = true;
-    for (int i = 0; i < num_events && has_complete_hint; ++i) {
-      if (crush.VarHasSolutionHint(PositiveRef(reservoir.active_literals(i)))) {
-        if (crush.LiteralSolutionHint(reservoir.active_literals(i))) {
-          const std::optional<int64_t> time_hint =
-              crush.GetExpressionSolutionHint(reservoir.time_exprs(i));
-          const std::optional<int64_t> change_hint =
-              crush.GetExpressionSolutionHint(reservoir.level_changes(i));
-          if (time_hint.has_value() && change_hint.has_value()) {
-            active_event_hints.push_back(
-                {i, time_hint.value(), change_hint.value()});
-          } else {
-            has_complete_hint = false;
-          }
-        }
-      } else {
-        has_complete_hint = false;
-      }
-    }
-  }
-  // Update the `level_vars` hints by computing the level at each active event.
-  if (has_complete_hint) {
-    std::sort(active_event_hints.begin(), active_event_hints.end(),
-              [](const ReservoirEventHint& a, const ReservoirEventHint& b) {
-                return a.time < b.time;
-              });
-    int64_t current_level = 0;
-    for (int i = 0; i < active_event_hints.size(); ++i) {
-      int j = i;
-      // Adjust the order of the events occurring at the same time, in the
-      // circuit, so that, at each node, the level is between `var_min` and
-      // `var_max`. For instance, if e1 = {t, +1} and e2 = {t, -1}, and if
-      // `current_level` = 0, `var_min` = -1 and `var_max` = 0, then e2 must
-      // occur before e1.
-      while (j < active_event_hints.size() &&
-             active_event_hints[j].time == active_event_hints[i].time &&
-             (current_level + active_event_hints[j].level_change < var_min ||
-              current_level + active_event_hints[j].level_change > var_max)) {
-        ++j;
-      }
-      if (j < active_event_hints.size() &&
-          active_event_hints[j].time == active_event_hints[i].time) {
-        if (i != j) std::swap(active_event_hints[i], active_event_hints[j]);
-        current_level += active_event_hints[i].level_change;
-        crush.UpdateVarSolutionHint(level_vars[active_event_hints[i].index],
-                                    current_level);
-      } else {
-        has_complete_hint = false;
-        break;
-      }
-    }
   }
 
   // For the corner case where all events are absent, we need a potential
@@ -148,17 +80,8 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
     circuit->add_tails(num_events);
     circuit->add_heads(num_events);
     circuit->add_literals(all_inactive);
-    if (has_complete_hint) {
-      crush.SetNewVariableHint(all_inactive, active_event_hints.empty());
-    }
   }
 
-  // The index of each event in `active_event_hints`, or -1 if the event's
-  // "active" hint is false.
-  std::vector<int> active_event_hint_index(num_events, -1);
-  for (int i = 0; i < active_event_hints.size(); ++i) {
-    active_event_hint_index[active_event_hints[i].index] = i;
-  }
   for (int i = 0; i < num_events; ++i) {
     if (!reservoir.active_literals().empty()) {
       // Add self arc to represent absence.
@@ -175,11 +98,6 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
       circuit->add_tails(num_events);
       circuit->add_heads(i);
       circuit->add_literals(start_var);
-      if (has_complete_hint) {
-        crush.SetNewVariableHint(start_var,
-                                 !active_event_hints.empty() &&
-                                     active_event_hints.front().index == i);
-      }
 
       // Add enforced linear for demand.
       {
@@ -200,11 +118,6 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
       circuit->add_tails(i);
       circuit->add_heads(num_events);
       circuit->add_literals(end_var);
-      if (has_complete_hint) {
-        crush.SetNewVariableHint(end_var,
-                                 !active_event_hints.empty() &&
-                                     active_event_hints.back().index == i);
-      }
     }
 
     for (int j = 0; j < num_events; ++j) {
@@ -224,13 +137,6 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
       circuit->add_tails(i);
       circuit->add_heads(j);
       circuit->add_literals(arc_i_j);
-      if (has_complete_hint) {
-        const int hint_i_index = active_event_hint_index[i];
-        const int hint_j_index = active_event_hint_index[j];
-        crush.SetNewVariableHint(arc_i_j, hint_i_index != -1 &&
-                                              hint_j_index != -1 &&
-                                              hint_j_index == hint_i_index + 1);
-      }
 
       // Add enforced linear for time.
       {
@@ -261,6 +167,8 @@ void ExpandReservoirUsingCircuit(int64_t sum_of_positive_demand,
       }
     }
   }
+  context->solution_crush().SetReservoirCircuitVars(reservoir, var_min, var_max,
+                                                    level_vars, *circuit);
 
   reservoir_ct->Clear();
   context->UpdateRuleStats("reservoir: expanded using circuit.");
@@ -337,6 +245,9 @@ void ExpandReservoirUsingPrecedences(bool max_level_is_constraining,
 
     // Canonicalize the newly created constraint.
     context->CanonicalizeLinearConstraint(new_cumul);
+
+    DCHECK(!PossibleIntegerOverflow(*context->working_model, new_linear->vars(),
+                                    new_linear->coeffs()));
   }
 
   reservoir_ct->Clear();
@@ -539,39 +450,6 @@ void ExpandIntMod(ConstraintProto* ct, PresolveContext* context) {
     return;
   }
 
-  bool has_complete_hint = false;
-  bool enforced_hint = true;
-  int64_t expr_hint = 0;
-  int64_t mod_expr_hint = 0;
-  int64_t target_expr_hint = 0;
-  SolutionCrush& crush = context->solution_crush();
-  if (crush.HintIsLoaded()) {
-    has_complete_hint = true;
-    for (const int lit : ct->enforcement_literal()) {
-      if (!crush.VarHasSolutionHint(PositiveRef(lit))) {
-        has_complete_hint = false;
-        break;
-      }
-      enforced_hint = enforced_hint && crush.LiteralSolutionHint(lit);
-    }
-    if (has_complete_hint && enforced_hint) {
-      has_complete_hint = false;
-      std::optional<int64_t> hint = crush.GetExpressionSolutionHint(expr);
-      if (hint.has_value()) {
-        expr_hint = hint.value();
-        hint = crush.GetExpressionSolutionHint(mod_expr);
-        if (hint.has_value()) {
-          mod_expr_hint = hint.value();
-          hint = crush.GetExpressionSolutionHint(target_expr);
-          if (hint.has_value()) {
-            target_expr_hint = hint.value();
-            has_complete_hint = true;
-          }
-        }
-      }
-    }
-  }
-
   // Create a new constraint with the same enforcement as ct.
   auto new_enforced_constraint = [&]() {
     ConstraintProto* new_ct = context->working_model->add_constraints();
@@ -583,13 +461,6 @@ void ExpandIntMod(ConstraintProto* ct, PresolveContext* context) {
   const int div_var = context->NewIntVar(
       context->DomainSuperSetOf(expr).PositiveDivisionBySuperset(
           context->DomainSuperSetOf(mod_expr)));
-  if (has_complete_hint) {
-    if (enforced_hint) {
-      crush.SetNewVariableHint(div_var, expr_hint / mod_expr_hint);
-    } else {
-      crush.SetNewVariableHint(div_var, context->MinOf(div_var));
-    }
-  }
   LinearExpressionProto div_expr;
   div_expr.add_vars(div_var);
   div_expr.add_coeffs(1);
@@ -607,13 +478,6 @@ void ExpandIntMod(ConstraintProto* ct, PresolveContext* context) {
           .IntersectionWith(context->DomainSuperSetOf(expr).AdditionWith(
               context->DomainSuperSetOf(target_expr).Negation()));
   const int prod_var = context->NewIntVar(prod_domain);
-  if (has_complete_hint) {
-    if (enforced_hint) {
-      crush.SetNewVariableHint(prod_var, expr_hint - target_expr_hint);
-    } else {
-      crush.SetNewVariableHint(prod_var, context->MinOf(prod_var));
-    }
-  }
   LinearExpressionProto prod_expr;
   prod_expr.add_vars(prod_var);
   prod_expr.add_coeffs(1);
@@ -633,15 +497,18 @@ void ExpandIntMod(ConstraintProto* ct, PresolveContext* context) {
   AddLinearExpressionToLinearConstraint(prod_expr, -1, lin);
   AddLinearExpressionToLinearConstraint(target_expr, -1, lin);
 
+  context->solution_crush().SetIntModExpandedVars(*ct, div_var, prod_var,
+                                                  context->MinOf(div_var),
+                                                  context->MinOf(prod_var));
   ct->Clear();
   context->UpdateRuleStats("int_mod: expanded");
 }
 
 void ExpandNonBinaryIntProd(ConstraintProto* ct, PresolveContext* context) {
   CHECK_GT(ct->int_prod().exprs_size(), 2);
-  SolutionCrush& crush = context->solution_crush();
   std::deque<LinearExpressionProto> terms(
       {ct->int_prod().exprs().begin(), ct->int_prod().exprs().end()});
+  std::vector<int> new_vars;
   while (terms.size() > 2) {
     const LinearExpressionProto& left = terms[0];
     const LinearExpressionProto& right = terms[1];
@@ -649,16 +516,7 @@ void ExpandNonBinaryIntProd(ConstraintProto* ct, PresolveContext* context) {
         context->DomainSuperSetOf(left).ContinuousMultiplicationBy(
             context->DomainSuperSetOf(right));
     const int new_var = context->NewIntVar(new_domain);
-    if (crush.HintIsLoaded()) {
-      const std::optional<int64_t> left_hint =
-          crush.GetExpressionSolutionHint(left);
-      const std::optional<int64_t> right_hint =
-          crush.GetExpressionSolutionHint(right);
-      if (left_hint.has_value() && right_hint.has_value()) {
-        crush.SetNewVariableHint(new_var,
-                                 left_hint.value() * right_hint.value());
-      }
-    }
+    new_vars.push_back(new_var);
     LinearArgumentProto* const int_prod =
         context->working_model->add_constraints()->mutable_int_prod();
     *int_prod->add_exprs() = left;
@@ -666,8 +524,7 @@ void ExpandNonBinaryIntProd(ConstraintProto* ct, PresolveContext* context) {
     int_prod->mutable_target()->add_vars(new_var);
     int_prod->mutable_target()->add_coeffs(1);
     terms.pop_front();
-    terms.pop_front();
-    terms.push_back(int_prod->target());
+    terms.front() = int_prod->target();
   }
 
   LinearArgumentProto* const final_int_prod =
@@ -676,6 +533,7 @@ void ExpandNonBinaryIntProd(ConstraintProto* ct, PresolveContext* context) {
   *final_int_prod->add_exprs() = terms[1];
   *final_int_prod->mutable_target() = ct->int_prod().target();
 
+  context->solution_crush().SetIntProdExpandedVars(ct->int_prod(), new_vars);
   context->UpdateRuleStats(absl::StrCat(
       "int_prod: expanded int_prod with arity ", ct->int_prod().exprs_size()));
   ct->Clear();
diff --git a/ortools/sat/cp_model_expand_test.cc b/ortools/sat/cp_model_expand_test.cc
index d3038f154b..1ed3d4a668 100644
--- a/ortools/sat/cp_model_expand_test.cc
+++ b/ortools/sat/cp_model_expand_test.cc
@@ -701,7 +701,8 @@ TEST(IntProdExpandTest, TestLargerAffineProd) {
 
 TEST(IntProdExpansionTest, ExpandNonBinaryIntProdPreservesSolutionHint) {
   CpModelProto initial_model = ParseTestProto(R"pb(
-    variables { domain: [ 0, 100 ] }
+    variables { domain: [ 0, 500 ] }
+    variables { domain: [ 0, 10 ] }
     variables { domain: [ 0, 10 ] }
     variables { domain: [ 0, 10 ] }
     variables { domain: [ 0, 10 ] }
@@ -711,6 +712,7 @@ TEST(IntProdExpansionTest, ExpandNonBinaryIntProdPreservesSolutionHint) {
         exprs { vars: 1 coeffs: 1 }
         exprs { vars: 2 coeffs: 1 }
         exprs { vars: 3 coeffs: 1 }
+        exprs { vars: 4 coeffs: 1 }
       }
     }
     objective {
@@ -718,8 +720,8 @@ TEST(IntProdExpansionTest, ExpandNonBinaryIntProdPreservesSolutionHint) {
       coeffs: [ 1 ]
     }
     solution_hint {
-      vars: [ 0, 1, 2, 3 ]
-      values: [ 60, 3, 4, 5 ]
+      vars: [ 0, 1, 2, 3, 4 ]
+      values: [ 360, 3, 4, 5, 6 ]
     }
   )pb");
 
diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc
index 096a74257d..786e43f660 100644
--- a/ortools/sat/cp_model_postsolve.cc
+++ b/ortools/sat/cp_model_postsolve.cc
@@ -18,6 +18,7 @@
 #include <limits>
 #include <vector>
 
+#include "absl/base/log_severity.h"
 #include "absl/log/check.h"
 #include "absl/types/span.h"
 #include "ortools/base/logging.h"
@@ -354,6 +355,30 @@ void PostsolveIntProd(const ConstraintProto& ct, std::vector<Domain>* domains) {
   (*domains)[target.vars(0)] = Domain(target_value);
 }
 
+namespace {
+void CheckPostsolveFixedVariables(const ConstraintProto& ct,
+                                  absl::Span<const Domain> domains) {
+  if (DEBUG_MODE) {
+    for (const int var : UsedVariables(ct)) {
+      DCHECK(domains[PositiveRef(var)].IsFixed())
+          << "Variable " << PositiveRef(var)
+          << " is not fixed after postsolve for constraint: \""
+          << ProtobufShortDebugString(ct)
+          << "\" with domain: " << domains[PositiveRef(var)] << ".";
+    }
+  }
+}
+
+void ArbitrarilyFixVariables(const ConstraintProto& ct,
+                             absl::Span<Domain> domains) {
+  for (const int var : UsedVariables(ct)) {
+    if (!domains[PositiveRef(var)].IsFixed()) {
+      domains[PositiveRef(var)] = Domain(domains[PositiveRef(var)].Min());
+    }
+  }
+}
+}  // namespace
+
 void PostsolveResponse(const int64_t num_variables_in_original_model,
                        const CpModelProto& mapping_proto,
                        absl::Span<const int> postsolve_mapping,
@@ -392,7 +417,11 @@ void PostsolveResponse(const int64_t num_variables_in_original_model,
         break;
       }
     }
-    if (constraint_can_be_ignored) continue;
+    if (constraint_can_be_ignored) {
+      ArbitrarilyFixVariables(ct, absl::MakeSpan(domains));
+      CheckPostsolveFixedVariables(ct, domains);
+      continue;
+    }
 
     switch (ct.constraint_case()) {
       case ConstraintProto::kBoolOr:
@@ -422,6 +451,7 @@ void PostsolveResponse(const int64_t num_variables_in_original_model,
         LOG(FATAL) << "Unsupported constraint: "
                    << ProtobufShortDebugString(ct);
     }
+    CheckPostsolveFixedVariables(ct, domains);
   }
 
   // Fill the response.
diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc
index d93d7f22a6..bbf3d08947 100644
--- a/ortools/sat/cp_model_presolve.cc
+++ b/ortools/sat/cp_model_presolve.cc
@@ -2478,7 +2478,14 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) {
   // instead of adding a reduced version of it each time a new singleton
   // variable appear in the same constraint later. That would work but would
   // also force the postsolve to take search decisions...
-  *context_->NewMappingConstraint(__FILE__, __LINE__) = *ct;
+  if (absl::GetFlag(FLAGS_cp_model_debug_postsolve)) {
+    auto* new_ct = context_->NewMappingConstraint(*ct, __FILE__, __LINE__);
+    const std::string name = new_ct->name();
+    *new_ct = *ct;
+    new_ct->set_name(absl::StrCat(ct->name(), " copy ", name));
+  } else {
+    *context_->NewMappingConstraint(*ct, __FILE__, __LINE__) = *ct;
+  }
 
   int new_size = 0;
   for (int i = 0; i < num_vars; ++i) {
diff --git a/ortools/sat/cp_model_symmetries.cc b/ortools/sat/cp_model_symmetries.cc
index d29bec4bce..4bd47c6fdf 100644
--- a/ortools/sat/cp_model_symmetries.cc
+++ b/ortools/sat/cp_model_symmetries.cc
@@ -966,7 +966,7 @@ void UpdateHintAfterFixingBoolToBreakSymmetry(
   const std::vector<int> generator_idx =
       TracePoint(target_var, schrier_vector, generators);
   for (const int i : generator_idx) {
-    crush.PermuteHintValues(*generators[i]);
+    crush.PermuteVariables(*generators[i]);
   }
 
   DCHECK(crush.VarHasSolutionHint(var));
diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc
index fd9b36b0db..6e41755bce 100644
--- a/ortools/sat/linear_relaxation.cc
+++ b/ortools/sat/linear_relaxation.cc
@@ -1299,9 +1299,10 @@ void AppendLinearConstraintRelaxation(const ConstraintProto& ct,
     // <=> Sum_i (~ei * (rhs_domain_min - min_activity)) + terms >=
     // rhs_domain_min
     LinearConstraintBuilder lc(model, rhs_domain_min, kMaxIntegerValue);
+    const IntegerValue term = CapSubI(rhs_domain_min, min_activity);
+    if (AtMinOrMaxInt64I(term) && !enforcing_literals.empty()) return;
     for (const Literal& literal : enforcing_literals) {
-      CHECK(
-          lc.AddLiteralTerm(literal.Negated(), rhs_domain_min - min_activity));
+      CHECK(lc.AddLiteralTerm(literal.Negated(), term));
     }
     for (int i = 0; i < ct.linear().vars_size(); i++) {
       const int ref = ct.linear().vars(i);
@@ -1319,9 +1320,10 @@ void AppendLinearConstraintRelaxation(const ConstraintProto& ct,
     // <=> Sum_i (~ei * (rhs_domain_max - max_activity)) + terms <=
     // rhs_domain_max
     LinearConstraintBuilder lc(model, kMinIntegerValue, rhs_domain_max);
+    const IntegerValue term = CapSubI(rhs_domain_max, max_activity);
+    if (AtMinOrMaxInt64I(term) && !enforcing_literals.empty()) return;
     for (const Literal& literal : enforcing_literals) {
-      CHECK(
-          lc.AddLiteralTerm(literal.Negated(), rhs_domain_max - max_activity));
+      CHECK(lc.AddLiteralTerm(literal.Negated(), term));
     }
     for (int i = 0; i < ct.linear().vars_size(); i++) {
       const int ref = ct.linear().vars(i);
diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc
index 8edfa6915b..33f1856557 100644
--- a/ortools/sat/presolve_context.cc
+++ b/ortools/sat/presolve_context.cc
@@ -89,9 +89,10 @@ int PresolveContext::NewIntVarWithDefinition(
 
   // Create new linear constraint new_var = definition.
   // TODO(user): When we encounter overflow (rare), we still create a variable.
-  auto* new_linear = append_constraint_to_mapping_model
-                         ? mapping_model->add_constraints()->mutable_linear()
-                         : working_model->add_constraints()->mutable_linear();
+  auto* new_linear =
+      append_constraint_to_mapping_model
+          ? NewMappingConstraint(__FILE__, __LINE__)->mutable_linear()
+          : working_model->add_constraints()->mutable_linear();
   for (const auto [var, coeff] : definition) {
     new_linear->add_vars(var);
     new_linear->add_coeffs(coeff);
@@ -2345,12 +2346,18 @@ int PresolveContext::GetOrCreateReifiedPrecedenceLiteral(
   const auto& it = reified_precedences_cache_.find(key);
   if (it != reified_precedences_cache_.end()) return it->second;
 
-  const int result = NewBoolVar("");
+  const int result = NewBoolVar("precedences");
   reified_precedences_cache_[key] = result;
 
   solution_crush_.SetVarToReifiedPrecedenceLiteral(result, time_i, time_j,
                                                    active_i, active_j);
 
+  if (!IsFixed(time_i) && !IsFixed(time_j)) {
+    DCHECK(!PossibleIntegerOverflow(*working_model,
+                                    {time_i.vars(0), time_j.vars(0)},
+                                    {-time_i.coeffs(0), time_j.coeffs(0)}));
+  }
+
   // result => (time_i <= time_j) && active_i && active_j.
   ConstraintProto* const lesseq = working_model->add_constraints();
   lesseq->add_enforcement_literal(result);
diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py
index e5bb3bdeb4..195535d0d6 100644
--- a/ortools/sat/python/cp_model.py
+++ b/ortools/sat/python/cp_model.py
@@ -420,11 +420,11 @@ class VariableList:
     def __init__(self) -> None:
         self.__var_list: list[IntVar] = []
 
-    def store(self, var: IntVar) -> None:
+    def append(self, var: IntVar) -> None:
         assert var.index == len(self.__var_list)
         self.__var_list.append(var)
 
-    def query(self, index: int) -> IntVar:
+    def get(self, index: int) -> IntVar:
         if index < 0 or index >= len(self.__var_list):
             raise ValueError("Index out of bounds.")
         return self.__var_list[index]
@@ -438,14 +438,14 @@ def rebuild_expr(
         if num_elements == 0:
             return proto.offset
         elif num_elements == 1:
-            var = self.query(proto.vars[0])
+            var = self.get(proto.vars[0])
             return LinearExpr.affine(
                 var, proto.coeffs[0], proto.offset
             )  # pytype: disable=bad-return-type
         else:
             variables = []
             for var_index in range(len(proto.vars)):
-                var = self.query(var_index)
+                var = self.get(var_index)
                 variables.append(var)
             if proto.offset != 0:
                 coeffs = []
@@ -644,13 +644,13 @@ def name(self, name: str):
 
     # Integer variable.
 
-    def _store_int_var(self, var: IntVar) -> IntVar:
+    def _append_int_var(self, var: IntVar) -> IntVar:
         """Appends an integer variable to the list of variables."""
-        self.__var_list.store(var)
+        self.__var_list.append(var)
         return var
 
-    def _query_int_var(self, index: int) -> IntVar:
-        return self.__var_list.query(index)
+    def _get_int_var(self, index: int) -> IntVar:
+        return self.__var_list.get(index)
 
     def rebuild_from_linear_expression_proto(
         self,
@@ -674,7 +674,7 @@ def new_int_var(self, lb: IntegralT, ub: IntegralT, name: str) -> IntVar:
           a variable whose domain is [lb, ub].
         """
         domain_is_boolean = lb >= 0 and ub <= 1
-        return self._store_int_var(
+        return self._append_int_var(
             IntVar(
                 self.__model,
                 sorted_interval_list.Domain(lb, ub),
@@ -700,20 +700,20 @@ def new_int_var_from_domain(
             a variable whose domain is the given domain.
         """
         domain_is_boolean = domain.min() >= 0 and domain.max() <= 1
-        return self._store_int_var(
+        return self._append_int_var(
             IntVar(self.__model, domain, domain_is_boolean, name)
         )
 
     def new_bool_var(self, name: str) -> IntVar:
         """Creates a 0-1 variable with the given name."""
-        return self._store_int_var(
+        return self._append_int_var(
             IntVar(self.__model, sorted_interval_list.Domain(0, 1), True, name)
         )
 
     def new_constant(self, value: IntegralT) -> IntVar:
         """Declares a constant integer."""
         index: int = self.get_or_make_index_from_constant(value)
-        return self._query_int_var(index)
+        return self._get_int_var(index)
 
     def new_int_var_series(
         self,
@@ -768,7 +768,7 @@ def new_int_var_series(
             index=index,
             data=[
                 # pylint: disable=g-complex-comprehension
-                self._store_int_var(
+                self._append_int_var(
                     IntVar(
                         model=self.__model,
                         name=f"{name}[{i}]",
@@ -808,7 +808,7 @@ def new_bool_var_series(
             index=index,
             data=[
                 # pylint: disable=g-complex-comprehension
-                self._store_int_var(
+                self._append_int_var(
                     IntVar(
                         model=self.__model,
                         name=f"{name}[{i}]",
@@ -2063,6 +2063,12 @@ def clone(self) -> "CpModel":
         clone.rebuild_var_and_constant_map()
         return clone
 
+    def __copy__(self):
+        return self.clone()
+
+    def __deepcopy__(self, memo):
+        return self.clone()
+
     def rebuild_var_and_constant_map(self):
         """Internal method used during model cloning."""
         for i, var in enumerate(self.__model.variables):
@@ -2071,11 +2077,11 @@ def rebuild_var_and_constant_map(self):
             is_boolean = (
                 len(var.domain) == 2 and var.domain[0] >= 0 and var.domain[1] <= 1
             )
-            self.__var_list.store(IntVar(self.__model, i, is_boolean, None))
+            self.__var_list.append(IntVar(self.__model, i, is_boolean, None))
 
     def get_bool_var_from_proto_index(self, index: int) -> IntVar:
         """Returns an already created Boolean variable from its index."""
-        result = self._query_int_var(index)
+        result = self._get_int_var(index)
         if not result.is_boolean:
             raise ValueError(
                 f"get_bool_var_from_proto_index: index {index} does not reference a"
@@ -2085,7 +2091,7 @@ def get_bool_var_from_proto_index(self, index: int) -> IntVar:
 
     def get_int_var_from_proto_index(self, index: int) -> IntVar:
         """Returns an already created integer variable from its index."""
-        return self._query_int_var(index)
+        return self._get_int_var(index)
 
     def get_interval_var_from_proto_index(self, index: int) -> IntervalVar:
         """Returns an already created interval variable from its index."""
diff --git a/ortools/sat/python/cp_model_test.py b/ortools/sat/python/cp_model_test.py
index 0daf692e3e..f7d23701ae 100644
--- a/ortools/sat/python/cp_model_test.py
+++ b/ortools/sat/python/cp_model_test.py
@@ -12,6 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+import copy
 import itertools
 import sys
 import time
@@ -1730,16 +1731,38 @@ def testCopyModel(self) -> None:
         lin = model.add(x + y <= 10)
 
         new_model = model.clone()
-        copy_b = new_model.get_bool_var_from_proto_index(b.index)
-        copy_x = new_model.get_int_var_from_proto_index(x.index)
-        copy_y = new_model.get_int_var_from_proto_index(y.index)
-        copy_i = new_model.get_interval_var_from_proto_index(i.index)
+        clone_b = new_model.get_bool_var_from_proto_index(b.index)
+        clone_x = new_model.get_int_var_from_proto_index(x.index)
+        clone_y = new_model.get_int_var_from_proto_index(y.index)
+        clone_i = new_model.get_interval_var_from_proto_index(i.index)
+
+        self.assertEqual(b.index, clone_b.index)
+        self.assertEqual(x.index, clone_x.index)
+        self.assertEqual(y.index, clone_y.index)
+        self.assertEqual(i.index, clone_i.index)
+
+        model_copy = copy.copy(model)
+        copy_b = model_copy.get_bool_var_from_proto_index(b.index)
+        copy_x = model_copy.get_int_var_from_proto_index(x.index)
+        copy_y = model_copy.get_int_var_from_proto_index(y.index)
+        copy_i = model_copy.get_interval_var_from_proto_index(i.index)
 
         self.assertEqual(b.index, copy_b.index)
         self.assertEqual(x.index, copy_x.index)
         self.assertEqual(y.index, copy_y.index)
         self.assertEqual(i.index, copy_i.index)
 
+        model_deepcopy = copy.deepcopy(model)
+        deepcopy_b = model_deepcopy.get_bool_var_from_proto_index(b.index)
+        deepcopy_x = model_deepcopy.get_int_var_from_proto_index(x.index)
+        deepcopy_y = model_deepcopy.get_int_var_from_proto_index(y.index)
+        deepcopy_i = model_deepcopy.get_interval_var_from_proto_index(i.index)
+
+        self.assertEqual(b.index, deepcopy_b.index)
+        self.assertEqual(x.index, deepcopy_x.index)
+        self.assertEqual(y.index, deepcopy_y.index)
+        self.assertEqual(i.index, deepcopy_i.index)
+
         with self.assertRaises(ValueError):
             new_model.get_bool_var_from_proto_index(-1)
 
diff --git a/ortools/sat/python/linear_expr.h b/ortools/sat/python/linear_expr.h
index 9b04174372..87e311c2f2 100644
--- a/ortools/sat/python/linear_expr.h
+++ b/ortools/sat/python/linear_expr.h
@@ -567,6 +567,10 @@ class NotBooleanVariable : public Literal {
   std::string DebugString() const override;
 
  private:
+  // We keep a weak ptr to the base variable to avoid a circular dependency.
+  // The base variable holds a shared pointer to the negated variable.
+  // Any call to a risky method is checked at the pybind11 level to raise a
+  // python exception before the call is made.
   std::weak_ptr<BaseIntVar> var_;
 };
 
diff --git a/ortools/sat/solution_crush.cc b/ortools/sat/solution_crush.cc
index db44664c38..d83c72b064 100644
--- a/ortools/sat/solution_crush.cc
+++ b/ortools/sat/solution_crush.cc
@@ -13,6 +13,7 @@
 
 #include "ortools/sat/solution_crush.h"
 
+#include <algorithm>
 #include <cstdint>
 #include <optional>
 #include <utility>
@@ -30,6 +31,35 @@
 namespace operations_research {
 namespace sat {
 
+void SolutionCrush::Resize(int new_size) {
+  hint_.resize(new_size, 0);
+  hint_has_value_.resize(new_size, false);
+}
+
+void SolutionCrush::LoadSolution(
+    const absl::flat_hash_map<int, int64_t>& solution) {
+  CHECK(!hint_is_loaded_);
+  model_has_hint_ = !solution.empty();
+  hint_is_loaded_ = true;
+  for (const auto [var, value] : solution) {
+    hint_has_value_[var] = true;
+    hint_[var] = value;
+  }
+}
+
+void SolutionCrush::MaybeSetLiteralToValueEncoding(int literal, int var,
+                                                   int64_t value) {
+  if (hint_is_loaded_) {
+    const int bool_var = PositiveRef(literal);
+    DCHECK(RefIsPositive(var));
+    if (!hint_has_value_[bool_var] && hint_has_value_[var]) {
+      const int64_t bool_value = hint_[var] == value ? 1 : 0;
+      hint_has_value_[bool_var] = true;
+      hint_[bool_var] = RefIsPositive(literal) ? bool_value : 1 - bool_value;
+    }
+  }
+}
+
 void SolutionCrush::SetVarToLinearExpression(
     int new_var, absl::Span<const std::pair<int, int64_t>> linear) {
   // We only fill the hint of the new variable if all the variable involved
@@ -94,25 +124,117 @@ void SolutionCrush::SetVarToConjunction(int new_var,
   }
 }
 
+void SolutionCrush::SetVarToValueIfLinearConstraintViolated(
+    int new_var, int64_t value,
+    absl::Span<const std::pair<int, int64_t>> linear, const Domain& domain) {
+  if (hint_is_loaded_) {
+    int64_t linear_value = 0;
+    bool all_have_hint = true;
+    for (const auto [var, coeff] : linear) {
+      if (!hint_has_value_[var]) {
+        all_have_hint = false;
+        break;
+      }
+      linear_value += hint_[var] * coeff;
+    }
+    if (all_have_hint && !domain.Contains(linear_value)) {
+      hint_has_value_[new_var] = true;
+      hint_[new_var] = value;
+    }
+  }
+}
+
+void SolutionCrush::SetLiteralToValueIfLinearConstraintViolated(
+    int literal, bool value, absl::Span<const std::pair<int, int64_t>> linear,
+    const Domain& domain) {
+  SetVarToValueIfLinearConstraintViolated(
+      PositiveRef(literal), RefIsPositive(literal) ? value : !value, linear,
+      domain);
+}
+
+void SolutionCrush::SetVarToValueIf(int var, int64_t value, int condition_lit) {
+  SetVarToValueIfLinearConstraintViolated(
+      var, value, {{PositiveRef(condition_lit), 1}},
+      Domain(RefIsPositive(condition_lit) ? 0 : 1));
+}
+
+void SolutionCrush::SetLiteralToValueIf(int literal, bool value,
+                                        int condition_lit) {
+  SetLiteralToValueIfLinearConstraintViolated(
+      literal, value, {{PositiveRef(condition_lit), 1}},
+      Domain(RefIsPositive(condition_lit) ? 0 : 1));
+}
+
+void SolutionCrush::MakeLiteralsEqual(int lit1, int lit2) {
+  if (!hint_is_loaded_) return;
+  if (hint_has_value_[PositiveRef(lit2)]) {
+    SetLiteralHint(lit1, LiteralSolutionHint(lit2));
+  } else if (hint_has_value_[PositiveRef(lit1)]) {
+    SetLiteralHint(lit2, LiteralSolutionHint(lit1));
+  }
+}
+
 void SolutionCrush::UpdateVarToDomain(int var, const Domain& domain) {
   if (VarHasSolutionHint(var)) {
     UpdateVarSolutionHint(var, domain.ClosestValue(SolutionHint(var)));
   }
 
 #ifdef CHECK_HINT
-  if (model_has_hint_ && HintIsLoaded() &&
-      !domains_[var].Contains(hint_[var])) {
+  if (model_has_hint_ && hint_is_loaded_ && !domain.Contains(hint_[var])) {
     LOG(FATAL) << "Hint with value " << hint_[var]
                << " infeasible when changing domain of " << var << " to "
-               << domains_[var];
+               << domain;
   }
 #endif
 }
 
-void SolutionCrush::PermuteHintValues(const SparsePermutation& perm) {
-  CHECK(hint_is_loaded_);
-  perm.ApplyToDenseCollection(hint_);
-  perm.ApplyToDenseCollection(hint_has_value_);
+void SolutionCrush::UpdateLiteralsToFalseIfDifferent(int lit1, int lit2) {
+  // Set lit1 and lit2 to false if "lit1 - lit2 == 0" is violated.
+  const int sign1 = RefIsPositive(lit1) ? 1 : -1;
+  const int sign2 = RefIsPositive(lit2) ? 1 : -1;
+  const std::vector<std::pair<int, int64_t>> linear = {
+      {PositiveRef(lit1), sign1}, {PositiveRef(lit2), -sign2}};
+  const Domain domain = Domain((sign1 == 1 ? 0 : -1) - (sign2 == 1 ? 0 : -1));
+  SetLiteralToValueIfLinearConstraintViolated(lit1, false, linear, domain);
+  SetLiteralToValueIfLinearConstraintViolated(lit2, false, linear, domain);
+}
+
+void SolutionCrush::UpdateLiteralsWithDominance(int lit, int dominating_lit) {
+  if (LiteralSolutionHintIs(lit, true) &&
+      LiteralSolutionHintIs(dominating_lit, false)) {
+    UpdateLiteralSolutionHint(lit, false);
+    UpdateLiteralSolutionHint(dominating_lit, true);
+  }
+}
+
+void SolutionCrush::UpdateRefsWithDominance(
+    int ref, int64_t min_value, int64_t max_value,
+    absl::Span<const std::pair<int, Domain>> dominating_refs) {
+  const std::optional<int64_t> ref_hint = GetRefSolutionHint(ref);
+  if (!ref_hint.has_value()) return;
+  // This can happen if the solution hint is not initially feasible (in which
+  // case we can't fix it).
+  if (*ref_hint < min_value) return;
+  // If the solution hint is already in the new domain there is nothing to do.
+  if (*ref_hint <= max_value) return;
+  // The quantity to subtract from the solution hint of `ref`.
+  const int64_t ref_hint_delta = *ref_hint - max_value;
+
+  UpdateRefSolutionHint(ref, *ref_hint - ref_hint_delta);
+  int64_t remaining_delta = ref_hint_delta;
+  for (const auto& [dominating_ref, dominating_ref_domain] : dominating_refs) {
+    const std::optional<int64_t> dominating_ref_hint =
+        GetRefSolutionHint(dominating_ref);
+    if (!dominating_ref_hint.has_value()) continue;
+    const int64_t new_dominating_ref_hint =
+        dominating_ref_domain.ValueAtOrBefore(*dominating_ref_hint +
+                                              remaining_delta);
+    // This might happen if the solution hint is not initially feasible.
+    if (!dominating_ref_domain.Contains(new_dominating_ref_hint)) continue;
+    UpdateRefSolutionHint(dominating_ref, new_dominating_ref_hint);
+    remaining_delta -= (new_dominating_ref_hint - *dominating_ref_hint);
+    if (remaining_delta == 0) break;
+  }
 }
 
 bool SolutionCrush::MaybeSetVarToAffineEquationSolution(int var_x, int var_y,
@@ -133,7 +255,7 @@ bool SolutionCrush::MaybeSetVarToAffineEquationSolution(int var_x, int var_y,
 #ifdef CHECK_HINT
   const int64_t vx = hint_[var_x];
   const int64_t vy = hint_[var_y];
-  if (model_has_hint_ && HintIsLoaded() && vx != vy * coeff + offset) {
+  if (model_has_hint_ && hint_is_loaded_ && vx != vy * coeff + offset) {
     LOG(FATAL) << "Affine relation incompatible with hint: " << vx
                << " != " << vy << " * " << coeff << " + " << offset;
   }
@@ -141,31 +263,90 @@ bool SolutionCrush::MaybeSetVarToAffineEquationSolution(int var_x, int var_y,
   return true;
 }
 
-void SolutionCrush::Resize(int new_size) {
-  hint_.resize(new_size, 0);
-  hint_has_value_.resize(new_size, false);
-}
+void SolutionCrush::SetReservoirCircuitVars(
+    const ReservoirConstraintProto& reservoir, int64_t min_level,
+    int64_t max_level, absl::Span<const int> level_vars,
+    const CircuitConstraintProto& circuit) {
+  if (!hint_is_loaded_) return;
+  // The hints of the active events, in the order they should appear in the
+  // circuit. The hints are collected first, and sorted later.
+  struct ReservoirEventHint {
+    int index;  // In the reservoir constraint.
+    int64_t time;
+    int64_t level_change;
+  };
+  const int num_events = reservoir.time_exprs_size();
+  std::vector<ReservoirEventHint> active_event_hints;
+  for (int i = 0; i < num_events; ++i) {
+    if (!VarHasSolutionHint(PositiveRef(reservoir.active_literals(i)))) return;
+    if (LiteralSolutionHint(reservoir.active_literals(i))) {
+      const std::optional<int64_t> time_hint =
+          GetExpressionSolutionHint(reservoir.time_exprs(i));
+      const std::optional<int64_t> change_hint =
+          GetExpressionSolutionHint(reservoir.level_changes(i));
+      if (!time_hint.has_value() || !change_hint.has_value()) return;
+      active_event_hints.push_back({i, time_hint.value(), change_hint.value()});
+    }
+  }
 
-void SolutionCrush::LoadSolution(
-    const absl::flat_hash_map<int, int64_t>& solution) {
-  CHECK(!hint_is_loaded_);
-  model_has_hint_ = !solution.empty();
-  hint_is_loaded_ = true;
-  for (const auto [var, value] : solution) {
-    hint_has_value_[var] = true;
-    hint_[var] = value;
+  // Update the `level_vars` hints by computing the level at each active event.
+  std::sort(active_event_hints.begin(), active_event_hints.end(),
+            [](const ReservoirEventHint& a, const ReservoirEventHint& b) {
+              return a.time < b.time;
+            });
+  int64_t current_level = 0;
+  for (int i = 0; i < active_event_hints.size(); ++i) {
+    int j = i;
+    // Adjust the order of the events occurring at the same time, in the
+    // circuit, so that, at each node, the level is between `var_min` and
+    // `var_max`. For instance, if e1 = {t, +1} and e2 = {t, -1}, and if
+    // `current_level` = 0, `var_min` = -1 and `var_max` = 0, then e2 must occur
+    // before e1.
+    while (j < active_event_hints.size() &&
+           active_event_hints[j].time == active_event_hints[i].time &&
+           (current_level + active_event_hints[j].level_change < min_level ||
+            current_level + active_event_hints[j].level_change > max_level)) {
+      ++j;
+    }
+    if (j < active_event_hints.size() &&
+        active_event_hints[j].time == active_event_hints[i].time) {
+      if (i != j) std::swap(active_event_hints[i], active_event_hints[j]);
+      current_level += active_event_hints[i].level_change;
+      SetVarHint(level_vars[active_event_hints[i].index], current_level);
+    } else {
+      return;
+    }
   }
-}
 
-void SolutionCrush::MaybeSetLiteralToValueEncoding(int literal, int var,
-                                                   int64_t value) {
-  if (hint_is_loaded_) {
-    const int bool_var = PositiveRef(literal);
-    DCHECK(RefIsPositive(var));
-    if (!hint_has_value_[bool_var] && hint_has_value_[var]) {
-      const int64_t bool_value = hint_[var] == value ? 1 : 0;
-      hint_has_value_[bool_var] = true;
-      hint_[bool_var] = RefIsPositive(literal) ? bool_value : 1 - bool_value;
+  // The index of each event in `active_event_hints`, or -1 if the event's
+  // "active" hint is false.
+  std::vector<int> active_event_hint_index(num_events, -1);
+  for (int i = 0; i < active_event_hints.size(); ++i) {
+    active_event_hint_index[active_event_hints[i].index] = i;
+  }
+  for (int i = 0; i < circuit.literals_size(); ++i) {
+    const int head = circuit.heads(i);
+    const int tail = circuit.tails(i);
+    const int literal = circuit.literals(i);
+    if (tail == num_events) {
+      if (head == num_events) {
+        // Self-arc on the start and end node.
+        SetLiteralHint(literal, active_event_hints.empty());
+      } else {
+        // Arc from the start node to an event node.
+        SetLiteralHint(literal, !active_event_hints.empty() &&
+                                    active_event_hints.front().index == head);
+      }
+    } else if (head == num_events) {
+      // Arc from an event node to the end node.
+      SetLiteralHint(literal, !active_event_hints.empty() &&
+                                  active_event_hints.back().index == tail);
+    } else if (tail != head) {
+      // Arc between two different event nodes.
+      const int tail_index = active_event_hint_index[tail];
+      const int head_index = active_event_hint_index[head];
+      SetLiteralHint(literal, tail_index != -1 && tail_index != -1 &&
+                                  head_index == tail_index + 1);
     }
   }
 }
@@ -189,5 +370,59 @@ void SolutionCrush::SetVarToReifiedPrecedenceLiteral(
   }
 }
 
+void SolutionCrush::SetIntModExpandedVars(const ConstraintProto& ct,
+                                          int div_var, int prod_var,
+                                          int64_t default_div_value,
+                                          int64_t default_prod_value) {
+  if (!hint_is_loaded_) return;
+  bool enforced_hint = true;
+  for (const int lit : ct.enforcement_literal()) {
+    if (!VarHasSolutionHint(PositiveRef(lit))) return;
+    enforced_hint = enforced_hint && LiteralSolutionHint(lit);
+  }
+  if (!enforced_hint) {
+    SetVarHint(div_var, default_div_value);
+    SetVarHint(prod_var, default_prod_value);
+    return;
+  }
+  const LinearArgumentProto& int_mod = ct.int_mod();
+  std::optional<int64_t> hint = GetExpressionSolutionHint(int_mod.exprs(0));
+  if (!hint.has_value()) return;
+  const int64_t expr_hint = hint.value();
+
+  hint = GetExpressionSolutionHint(int_mod.exprs(1));
+  if (!hint.has_value()) return;
+  const int64_t mod_expr_hint = hint.value();
+
+  hint = GetExpressionSolutionHint(int_mod.target());
+  if (!hint.has_value()) return;
+  const int64_t target_expr_hint = hint.value();
+
+  // target_expr_hint should be equal to "expr_hint % mod_expr_hint".
+  SetVarHint(div_var, expr_hint / mod_expr_hint);
+  SetVarHint(prod_var, expr_hint - target_expr_hint);
+}
+
+void SolutionCrush::SetIntProdExpandedVars(const LinearArgumentProto& int_prod,
+                                           absl::Span<const int> prod_vars) {
+  DCHECK_EQ(int_prod.exprs_size(), prod_vars.size() + 2);
+  if (!hint_is_loaded_) return;
+  std::optional<int64_t> hint = GetExpressionSolutionHint(int_prod.exprs(0));
+  if (!hint.has_value()) return;
+  int64_t last_prod_hint = hint.value();
+  for (int i = 1; i < int_prod.exprs_size() - 1; ++i) {
+    hint = GetExpressionSolutionHint(int_prod.exprs(i));
+    if (!hint.has_value()) return;
+    last_prod_hint *= hint.value();
+    SetVarHint(prod_vars[i - 1], last_prod_hint);
+  }
+}
+
+void SolutionCrush::PermuteVariables(const SparsePermutation& permutation) {
+  CHECK(hint_is_loaded_);
+  permutation.ApplyToDenseCollection(hint_);
+  permutation.ApplyToDenseCollection(hint_has_value_);
+}
+
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/solution_crush.h b/ortools/sat/solution_crush.h
index dddc37ee6e..cdee5e0e6d 100644
--- a/ortools/sat/solution_crush.h
+++ b/ortools/sat/solution_crush.h
@@ -39,14 +39,19 @@ class SolutionCrush {
   SolutionCrush& operator=(const SolutionCrush&) = delete;
   SolutionCrush& operator=(SolutionCrush&&) = delete;
 
+  // Resizes the solution to contain `new_size` variables. Does not change the
+  // value of existing variables, and does not set any value for the new
+  // variables.
+  // WARNING: the methods below do not automatically resize the solution. To set
+  // the value of a new variable with one of them, call this method first.
+  void Resize(int new_size);
+
   // Sets the given values in the solution. `solution` must be a map from
   // variable indices to variable values. This must be called only once, before
   // any other method (besides `Resize`).
   // TODO(user): revisit this near the end of the refactoring.
   void LoadSolution(const absl::flat_hash_map<int, int64_t>& solution);
 
-  void PermuteHintValues(const SparsePermutation& perm);
-
   // Solution hint accessors.
   bool VarHasSolutionHint(int var) const { return hint_has_value_[var]; }
   int64_t SolutionHint(int var) const { return hint_[var]; }
@@ -118,17 +123,9 @@ class SolutionCrush {
   void SetNewVariableHint(int var, int64_t value) {
     CHECK(hint_is_loaded_);
     CHECK(!hint_has_value_[var]);
-    hint_has_value_[var] = true;
-    hint_[var] = value;
+    SetVarHint(var, value);
   }
 
-  // Resizes the solution to contain `new_size` variables. Does not change the
-  // value of existing variables, and does not set any value for the new
-  // variables.
-  // WARNING: the methods below do not automatically resize the solution. To set
-  // the value of a new variable with one of them, call this method first.
-  void Resize(int new_size);
-
   // Sets the value of `literal` to "`var`'s value == `value`". Does nothing if
   // `literal` already has a value.
   void MaybeSetLiteralToValueEncoding(int literal, int var, int64_t value);
@@ -148,17 +145,73 @@ class SolutionCrush {
   // literal indices.
   void SetVarToConjunction(int var, absl::Span<const int> conjunction);
 
+  // Sets the value of `var` to `value` if the value of the given linear
+  // expression is not in `domain` (or does nothing otherwise). `linear` must be
+  // a list of (variable index, coefficient) pairs.
+  void SetVarToValueIfLinearConstraintViolated(
+      int var, int64_t value, absl::Span<const std::pair<int, int64_t>> linear,
+      const Domain& domain);
+
+  // Sets the value of `literal` to `value` if the value of the given linear
+  // expression is not in `domain` (or does nothing otherwise). `linear` must be
+  // a list of (variable index, coefficient) pairs.
+  void SetLiteralToValueIfLinearConstraintViolated(
+      int literal, bool value, absl::Span<const std::pair<int, int64_t>> linear,
+      const Domain& domain);
+
+  // Sets the value of `var` to `value` if the value of `condition_lit` is true.
+  void SetVarToValueIf(int var, int64_t value, int condition_lit);
+
+  // Sets the value of `literal` to `value` if the value of `condition_lit` is
+  // true.
+  void SetLiteralToValueIf(int literal, bool value, int condition_lit);
+
+  // If one literal does not have a value, and the other one does, sets the
+  // value of the latter to the value of the former. If both literals have a
+  // value, sets the value of `lit1` to the value of `lit2`.
+  void MakeLiteralsEqual(int lit1, int lit2);
+
   // Updates the value of the given variable to be within the given domain. The
   // variable is updated to the closest value within the domain. `var` must
   // already have a value.
   void UpdateVarToDomain(int var, const Domain& domain);
 
+  // Updates the value of the given literals to false if their current values
+  // are different (or does nothing otherwise).
+  void UpdateLiteralsToFalseIfDifferent(int lit1, int lit2);
+
+  // Decrements the value of `lit` and increments the value of `dominating_lit`
+  // if their values are equal to 1 and 0, respectively.
+  void UpdateLiteralsWithDominance(int lit, int dominating_lit);
+
+  // Decrements the value of `ref` by the minimum amount necessary to be in
+  // [`min_value`, `max_value`], and increments the value of one or more
+  // `dominating_refs` by the same total amount (or less if it is not possible
+  // to exactly match this amount), while staying within their respective
+  // domains. The value of a negative reference index r is the opposite of the
+  // value of the variable PositiveRef(r).
+  //
+  // `min_value` must be the minimum value of `ref`'s current domain D, and
+  // `max_value` must be in D.
+  void UpdateRefsWithDominance(
+      int ref, int64_t min_value, int64_t max_value,
+      absl::Span<const std::pair<int, Domain>> dominating_refs);
+
   // Sets the value of `var_y` so that "`var_x`'s value = `var_y`'s value
   // * `coeff` + `offset`". Does nothing if `var_y` already has a value.
   // Returns whether the update was successful.
   bool MaybeSetVarToAffineEquationSolution(int var_x, int var_y, int64_t coeff,
                                            int64_t offset);
 
+  // Sets the value of the variables in `level_vars` and in `circuit` if all the
+  // variables in `reservoir` have a value. This assumes that there is one level
+  // variable and one circuit node per element in `reservoir` (in the same
+  // order) -- plus one last node representing the start and end of the circuit.
+  void SetReservoirCircuitVars(const ReservoirConstraintProto& reservoir,
+                               int64_t min_level, int64_t max_level,
+                               absl::Span<const int> level_vars,
+                               const CircuitConstraintProto& circuit);
+
   // Sets the value of `var` to "`time_i`'s value <= `time_j`'s value &&
   // `active_i`'s value == true && `active_j`'s value == true".
   void SetVarToReifiedPrecedenceLiteral(int var,
@@ -166,7 +219,39 @@ class SolutionCrush {
                                         const LinearExpressionProto& time_j,
                                         int active_i, int active_j);
 
+  // Sets the value of `div_var` and `prod_var` if all the variables in the
+  // IntMod `ct` constraint have a value, assuming that this "target = x % m"
+  // constraint is expanded into "div_var = x / m", "prod_var = div_var * m",
+  // and "target = x - prod_var" constraints. If `ct` is not enforced, sets the
+  // values of `div_var` and `prod_var` to `default_div_value` and
+  // `default_prod_value`, respectively.
+  void SetIntModExpandedVars(const ConstraintProto& ct, int div_var,
+                             int prod_var, int64_t default_div_value,
+                             int64_t default_prod_value);
+
+  // Sets the value of as many variables in `prod_vars` as possible (depending
+  // on how many expressions in `int_prod` have a value), assuming that the
+  // `int_prod` constraint "target = x_0 * x_1 * ... * x_n" is expanded into
+  // "prod_var_1 = x_0 * x1",
+  // "prod_var_2 = prod_var_1 * x_2",
+  //  ...,
+  // "prod_var_(n-1) = prod_var_(n-2) * x_(n-1)",
+  // and "target = prod_var_(n-1) * x_n" constraints.
+  void SetIntProdExpandedVars(const LinearArgumentProto& int_prod,
+                              absl::Span<const int> prod_vars);
+
+  void PermuteVariables(const SparsePermutation& permutation);
+
  private:
+  void SetVarHint(int var, int64_t value) {
+    hint_has_value_[var] = true;
+    hint_[var] = value;
+  }
+
+  void SetLiteralHint(int lit, bool value) {
+    SetVarHint(PositiveRef(lit), RefIsPositive(lit) == value ? 1 : 0);
+  }
+
   // This contains all the hinted value or zero if the hint wasn't specified.
   // We try to maintain this as we create new variable.
   bool model_has_hint_ = false;
diff --git a/ortools/sat/var_domination.cc b/ortools/sat/var_domination.cc
index 1e84abcf70..243b8a08ab 100644
--- a/ortools/sat/var_domination.cc
+++ b/ortools/sat/var_domination.cc
@@ -883,9 +883,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
           // decreasing [`ct` = enf => (var = implied)] does not apply. We can
           // thus set the hint of `positive_ref` to `bound` to preserve the hint
           // feasibility.
-          if (crush.LiteralSolutionHintIs(enf, false)) {
-            crush.UpdateRefSolutionHint(positive_ref, bound);
-          }
+          crush.SetVarToValueIf(positive_ref, bound, NegatedRef(enf));
           if (RefIsPositive(enf)) {
             // positive_ref = enf * implied + (1 - enf) * bound.
             if (!context->StoreAffineRelation(
@@ -914,9 +912,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
           // hint(ref) is 1. In this case the only locking constraint `ct` does
           // not apply and thus does not prevent decreasing the hint of ref in
           // order to preserve the hint feasibility.
-          if (crush.LiteralSolutionHintIs(enf, false)) {
-            crush.UpdateLiteralSolutionHint(ref, false);
-          }
+          crush.SetLiteralToValueIf(ref, false, NegatedRef(enf));
           context->AddImplication(NegatedRef(enf), NegatedRef(ref));
           context->UpdateNewConstraintsVariableUsage();
           continue;
@@ -972,9 +968,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
               // break the hint only if hint(ref) = hint(encoding_lit) = 1. But
               // in this case `ct` is actually not blocking ref from decreasing.
               // We can thus set its hint to 0 to preserve the hint feasibility.
-              if (crush.LiteralSolutionHintIs(encoding_lit, true)) {
-                crush.UpdateLiteralSolutionHint(ref, false);
-              }
+              crush.SetLiteralToValueIf(ref, false, encoding_lit);
               if (!context->StoreBooleanEqualityRelation(encoding_lit,
                                                          NegatedRef(ref))) {
                 return false;
@@ -986,9 +980,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
               // = 1. But in this case `ct` is actually not blocking ref from
               // decreasing. We can thus set its hint to 0 to preserve the hint
               // feasibility.
-              if (crush.LiteralSolutionHintIs(encoding_lit, false)) {
-                crush.UpdateLiteralSolutionHint(ref, false);
-              }
+              crush.SetLiteralToValueIf(ref, false, NegatedRef(encoding_lit));
               if (!context->StoreBooleanEqualityRelation(encoding_lit, ref)) {
                 return false;
               }
@@ -1012,10 +1004,8 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
           // set the hint of `ref` to false. This should be safe since the only
           // constraint blocking `ref` from decreasing is `ct` = not(ref) =>
           // (`var` in `rhs`) -- which does not apply when `ref` is true.
-          const std::optional<int64_t> var_hint = crush.GetRefSolutionHint(var);
-          if (var_hint.has_value() && !complement.Contains(*var_hint)) {
-            crush.UpdateLiteralSolutionHint(ref, false);
-          }
+          crush.SetLiteralToValueIfLinearConstraintViolated(
+              ref, false, {{var, 1}}, complement);
           ConstraintProto* new_ct = context->working_model->add_constraints();
           new_ct->add_enforcement_literal(ref);
           new_ct->mutable_linear()->add_vars(var);
@@ -1095,13 +1085,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
               // not blocking the ref at 1 from decreasing. Hence we can set its
               // hint to false to preserve the hint feasibility despite the new
               // Boolean equality constraint.
-              if (crush.VarHasSolutionHint(PositiveRef(ref)) &&
-                  crush.VarHasSolutionHint(PositiveRef(other_ref)) &&
-                  crush.LiteralSolutionHint(ref) !=
-                      crush.LiteralSolutionHint(other_ref)) {
-                crush.UpdateLiteralSolutionHint(ref, false);
-                crush.UpdateLiteralSolutionHint(other_ref, false);
-              }
+              crush.UpdateLiteralsToFalseIfDifferent(ref, other_ref);
               if (!context->StoreBooleanEqualityRelation(ref, other_ref)) {
                 return false;
               }
@@ -1176,10 +1160,8 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) {
     // If hint(a) is false we can always set it to hint(b) since this can only
     // increase its value. If hint(a) is true then hint(b) must be true as well
     // if the hint is feasible, due to the a => b constraint. Setting hint(a) to
-    // hint(b) is thus always safe.
-    if (crush.VarHasSolutionHint(PositiveRef(b))) {
-      crush.UpdateLiteralSolutionHint(a, crush.LiteralSolutionHint(b));
-    }
+    // hint(b) is thus always safe. The opposite is true as well.
+    crush.MakeLiteralsEqual(a, b);
     if (!context->StoreBooleanEqualityRelation(a, b)) return false;
     context->UpdateRuleStats("dual: enforced equivalence");
   }
@@ -1472,19 +1454,6 @@ void ScanModelForDualBoundStrengthening(
 }
 
 namespace {
-// Decrements the solution hint of `lit` and increments the solution hint of
-// `dominating_lit` if both hint values are present and equal to 1 and 0,
-// respectively.
-void MaybeUpdateLiteralHintFromDominance(PresolveContext& context, int lit,
-                                         int dominating_lit) {
-  SolutionCrush& crush = context.solution_crush();
-  if (crush.LiteralSolutionHintIs(lit, true) &&
-      crush.LiteralSolutionHintIs(dominating_lit, false)) {
-    crush.UpdateLiteralSolutionHint(lit, false);
-    crush.UpdateLiteralSolutionHint(dominating_lit, true);
-  }
-}
-
 // Decrements the solution hint of `ref` by the minimum amount necessary to be
 // in `domain`, and increments the solution hint of one or more
 // `dominating_variables` by the same total amount (or less if it is not
@@ -1494,36 +1463,19 @@ void MaybeUpdateLiteralHintFromDominance(PresolveContext& context, int lit,
 // domain D in `context`, and whose upper bound must be in D.
 void MaybeUpdateRefHintFromDominance(
     PresolveContext& context, int ref, const Domain& domain,
-    const absl::Span<const IntegerVariable> dominating_variables) {
-  SolutionCrush& crush = context.solution_crush();
-  const std::optional<int64_t> ref_hint = crush.GetRefSolutionHint(ref);
-  if (!ref_hint.has_value()) return;
-  // The quantity to subtract from the solution hint of `ref`. If the closest
-  // value of *ref_hint in `domain` is not *ref_hint then it is either the lower
-  // or upper bound of `domain`, which by hypothesis are in `ref`'s current
-  // domain D. Hence, in any case, this closest value is in D.
-  const int64_t ref_hint_delta = *ref_hint - domain.ClosestValue(*ref_hint);
-  // If it is 0 there is nothing to do. It might be negative if the solution
-  // hint is not initially feasible (in which case we can't fix it).
-  if (ref_hint_delta <= 0) return;
-
-  crush.UpdateRefSolutionHint(ref, *ref_hint - ref_hint_delta);
-  int64_t remaining_delta = ref_hint_delta;
-  for (const IntegerVariable ivar : dominating_variables) {
-    const int dominating_ref = VarDomination::IntegerVariableToRef(ivar);
-    const std::optional<int64_t> dominating_ref_hint =
-        crush.GetRefSolutionHint(dominating_ref);
-    if (!dominating_ref_hint.has_value()) continue;
-    const Domain& dominating_ref_domain = context.DomainOf(dominating_ref);
-    const int64_t new_dominating_ref_hint =
-        dominating_ref_domain.ValueAtOrBefore(*dominating_ref_hint +
-                                              remaining_delta);
-    // This might happen if the solution hint is not initially feasible.
-    if (!dominating_ref_domain.Contains(new_dominating_ref_hint)) continue;
-    crush.UpdateRefSolutionHint(dominating_ref, new_dominating_ref_hint);
-    remaining_delta -= (new_dominating_ref_hint - *dominating_ref_hint);
-    if (remaining_delta == 0) break;
+    absl::Span<const IntegerVariable> dominating_variables) {
+  DCHECK_EQ(domain.NumIntervals(), 1);
+  DCHECK_EQ(domain.Min(), context.DomainOf(ref).Min());
+  DCHECK(context.DomainOf(ref).Contains(domain.Max()));
+  std::vector<std::pair<int, Domain>> dominating_refs;
+  dominating_refs.reserve(dominating_variables.size());
+  for (const IntegerVariable var : dominating_variables) {
+    const int dominating_ref = VarDomination::IntegerVariableToRef(var);
+    dominating_refs.push_back(
+        {dominating_ref, context.DomainOf(dominating_ref)});
   }
+  context.solution_crush().UpdateRefsWithDominance(
+      ref, domain.Min(), domain.Max(), dominating_refs);
 }
 
 bool ProcessAtMostOne(
@@ -1550,7 +1502,7 @@ bool ProcessAtMostOne(
       // constraint. Hence the hint feasibility can always be preserved (if the
       // hint value of `ref` is 0 the hint does not need to be updated).
       context->UpdateRuleStats(message);
-      MaybeUpdateLiteralHintFromDominance(*context, ref, iref);
+      context->solution_crush().UpdateLiteralsWithDominance(ref, iref);
       if (!context->SetLiteralToFalse(ref)) return false;
       break;
     }
@@ -1599,6 +1551,7 @@ bool ExploitDominanceRelations(const VarDomination& var_domination,
   util_intops::StrongVector<IntegerVariable, bool> in_constraints(num_vars * 2,
                                                                   false);
 
+  SolutionCrush& crush = context->solution_crush();
   absl::flat_hash_set<std::pair<int, int>> implications;
   const int num_constraints = cp_model.constraints_size();
   for (int c = 0; c < num_constraints; ++c) {
@@ -1622,7 +1575,7 @@ bool ExploitDominanceRelations(const VarDomination& var_domination,
           const int ref = VarDomination::IntegerVariableToRef(ivar);
           if (ref == NegatedRef(b)) {
             context->UpdateRuleStats("domination: in implication");
-            MaybeUpdateLiteralHintFromDominance(*context, a, ref);
+            crush.UpdateLiteralsWithDominance(a, ref);
             if (!context->SetLiteralToFalse(a)) return false;
             break;
           }
@@ -1638,7 +1591,7 @@ bool ExploitDominanceRelations(const VarDomination& var_domination,
           const int ref = VarDomination::IntegerVariableToRef(ivar);
           if (ref == a) {
             context->UpdateRuleStats("domination: in implication");
-            MaybeUpdateLiteralHintFromDominance(*context, NegatedRef(b), ref);
+            crush.UpdateLiteralsWithDominance(NegatedRef(b), ref);
             if (!context->SetLiteralToTrue(b)) return false;
             break;
           }
@@ -1978,7 +1931,7 @@ bool ExploitDominanceRelations(const VarDomination& var_domination,
         // of `ref` is 1 and the hint value of `dom_ref` is 0. In this case the
         // call below fixes it by negating both values. Otherwise it does
         // nothing and thus preserves its feasibility.
-        MaybeUpdateLiteralHintFromDominance(*context, ref, dom_ref);
+        crush.UpdateLiteralsWithDominance(ref, dom_ref);
         context->UpdateNewConstraintsVariableUsage();
         implications.insert({ref, dom_ref});
         implications.insert({NegatedRef(dom_ref), NegatedRef(ref)});