diff --git a/include/samurai/petsc/block_assembly.hpp b/include/samurai/petsc/block_assembly.hpp
index 90a4afc00..995c1492a 100644
--- a/include/samurai/petsc/block_assembly.hpp
+++ b/include/samurai/petsc/block_assembly.hpp
@@ -26,8 +26,8 @@ namespace samurai
 
             explicit BlockAssemblyBase(const block_operator_t& block_op)
                 : m_block_operator(block_op)
-                , m_assembly_ops(transform(block_op.operators(),
-                                           [](const auto& op)
+                , m_assembly_ops(transform(m_block_operator.operators(),
+                                           [](auto& op)
                                            {
                                                return make_assembly(op);
                                            }))
diff --git a/include/samurai/petsc/fv/FV_scheme_assembly.hpp b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
index 958fa696c..ca9862f3a 100644
--- a/include/samurai/petsc/fv/FV_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
@@ -74,7 +74,12 @@ namespace samurai
                 this->set_name(scheme.name());
             }
 
-            auto& scheme() const
+            auto& scheme()
+            {
+                return m_scheme;
+            }
+
+            const auto& scheme() const
             {
                 return m_scheme;
             }
diff --git a/include/samurai/petsc/fv/operator_sum_assembly.hpp b/include/samurai/petsc/fv/operator_sum_assembly.hpp
index 39b593c86..deae7a037 100644
--- a/include/samurai/petsc/fv/operator_sum_assembly.hpp
+++ b/include/samurai/petsc/fv/operator_sum_assembly.hpp
@@ -16,16 +16,16 @@ namespace samurai
 
           private:
 
-            const scheme_t* m_sum_scheme;
+            scheme_t* m_sum_scheme;
 
             std::tuple<Assembly<Operators>...> m_assembly_ops;
 
           public:
 
-            explicit Assembly(const scheme_t& sum_scheme)
+            explicit Assembly(scheme_t& sum_scheme)
                 : m_sum_scheme(&sum_scheme)
                 , m_assembly_ops(transform(sum_scheme.operators(),
-                                           [](const auto& op)
+                                           [](auto& op)
                                            {
                                                return make_assembly(op);
                                            }))
@@ -67,6 +67,11 @@ namespace samurai
                          });
             }
 
+            auto& scheme()
+            {
+                return *m_sum_scheme;
+            }
+
             auto& scheme() const
             {
                 return *m_sum_scheme;
diff --git a/include/samurai/petsc/linear_block_solver.hpp b/include/samurai/petsc/linear_block_solver.hpp
index bcbe16f0c..a63be3fbc 100644
--- a/include/samurai/petsc/linear_block_solver.hpp
+++ b/include/samurai/petsc/linear_block_solver.hpp
@@ -33,7 +33,7 @@ namespace samurai
 
             static constexpr bool is_monolithic = monolithic;
 
-            explicit LinearBlockSolver(const block_operator_t& block_op)
+            explicit LinearBlockSolver(block_operator_t& block_op)
                 : base_class(block_op)
             {
                 _configure_solver();
diff --git a/include/samurai/petsc/linear_solver.hpp b/include/samurai/petsc/linear_solver.hpp
index f1959c524..9749a03bb 100644
--- a/include/samurai/petsc/linear_solver.hpp
+++ b/include/samurai/petsc/linear_solver.hpp
@@ -28,7 +28,7 @@ namespace samurai
 
           public:
 
-            explicit LinearSolverBase(const scheme_t& scheme)
+            explicit LinearSolverBase(scheme_t& scheme)
                 : m_assembly(scheme)
             {
                 _configure_solver();
@@ -252,7 +252,7 @@ namespace samurai
 
           public:
 
-            explicit LinearSolver(const scheme_t& scheme)
+            explicit LinearSolver(scheme_t& scheme)
                 : base_class(scheme)
             {
                 _configure_solver();
diff --git a/include/samurai/petsc/matrix_assembly.hpp b/include/samurai/petsc/matrix_assembly.hpp
index 2366da65d..4664148f5 100644
--- a/include/samurai/petsc/matrix_assembly.hpp
+++ b/include/samurai/petsc/matrix_assembly.hpp
@@ -352,11 +352,11 @@ namespace samurai
         };
 
         template <class Scheme>
-        auto make_assembly(const Scheme& s)
+        auto make_assembly(Scheme& s)
         {
             if constexpr (std::is_base_of_v<MatrixAssembly, Scheme>)
             {
-                return *reinterpret_cast<const Assembly<Scheme>*>(&s);
+                return *reinterpret_cast<Assembly<Scheme>*>(&s);
             }
             else
             {
diff --git a/include/samurai/petsc/nonlinear_solver.hpp b/include/samurai/petsc/nonlinear_solver.hpp
index 215678db7..06aeeb876 100644
--- a/include/samurai/petsc/nonlinear_solver.hpp
+++ b/include/samurai/petsc/nonlinear_solver.hpp
@@ -24,7 +24,7 @@ namespace samurai
 
           public:
 
-            explicit NonLinearSolverBase(const scheme_t& scheme)
+            explicit NonLinearSolverBase(scheme_t& scheme)
                 : m_assembly(scheme)
             {
                 _configure_solver();
@@ -308,7 +308,7 @@ namespace samurai
             using base_class::m_is_set_up;
             using base_class::m_J;
 
-            explicit NonLinearSolver(const scheme_t& scheme)
+            explicit NonLinearSolver(scheme_t& scheme)
                 : base_class(scheme)
             {
             }
diff --git a/include/samurai/petsc/solver_helpers.hpp b/include/samurai/petsc/solver_helpers.hpp
index 7d8a62260..c0352655f 100644
--- a/include/samurai/petsc/solver_helpers.hpp
+++ b/include/samurai/petsc/solver_helpers.hpp
@@ -14,21 +14,21 @@ namespace samurai
 
         // Linear solver
         template <class Scheme, std::enable_if_t<Scheme::cfg_t::scheme_type != SchemeType::NonLinear, bool> = true>
-        auto make_solver(const Scheme& scheme)
+        auto make_solver(Scheme& scheme)
         {
             return LinearSolver<Scheme>(scheme);
         }
 
         // Linear block solver (choice monolithic or not)
         template <bool monolithic, std::size_t rows, std::size_t cols, class... Operators>
-        auto make_solver(const BlockOperator<rows, cols, Operators...>& block_operator)
+        auto make_solver(BlockOperator<rows, cols, Operators...>& block_operator)
         {
             return LinearBlockSolver<monolithic, rows, cols, Operators...>(block_operator);
         }
 
         // Linear block solver (monolithic)
         template <std::size_t rows, std::size_t cols, class... Operators>
-        auto make_solver(const BlockOperator<rows, cols, Operators...>& block_operator)
+        auto make_solver(BlockOperator<rows, cols, Operators...>& block_operator)
         {
             static constexpr bool default_monolithic = true;
             return make_solver<default_monolithic, rows, cols, Operators...>(block_operator);
@@ -36,14 +36,14 @@ namespace samurai
 
         // Non-linear solver
         template <class Scheme, std::enable_if_t<Scheme::cfg_t::scheme_type == SchemeType::NonLinear, bool> = true>
-        auto make_solver(const Scheme& scheme)
+        auto make_solver(Scheme& scheme)
         {
             return NonLinearSolver<Scheme>(scheme);
         }
 
         // Non-linear local solvers
         template <class cfg, class bdry_cfg, std::enable_if_t<cfg::scheme_type == SchemeType::NonLinear && cfg::stencil_size == 1, bool> = true>
-        auto make_solver(const CellBasedScheme<cfg, bdry_cfg>& scheme)
+        auto make_solver(CellBasedScheme<cfg, bdry_cfg>& scheme)
         {
             return NonLinearLocalSolvers<CellBasedScheme<cfg, bdry_cfg>>(scheme);
         }
@@ -53,18 +53,26 @@ namespace samurai
          */
 
         template <class Scheme>
-        void solve(const Scheme& scheme, typename Scheme::field_t& unknown, typename Scheme::field_t& rhs)
+        void solve(Scheme& scheme, typename Scheme::field_t& unknown, typename Scheme::field_t& rhs)
         {
             auto solver = make_solver(scheme);
             solver.solve(unknown, rhs);
         }
 
         template <class Scheme, class E>
-        void solve(const Scheme& scheme, typename Scheme::field_t& unknown, const field_expression<E>& rhs_expression)
+        void solve(Scheme& scheme, typename Scheme::field_t& unknown, const field_expression<E>& rhs_expression)
         {
             typename Scheme::field_t rhs = rhs_expression;
             solve(scheme, unknown, rhs);
         }
 
+        template <class Scheme>
+        void solve(Scheme&& scheme, typename Scheme::field_t& unknown, typename Scheme::field_t& rhs)
+        {
+            Scheme scheme_ = std::move(scheme);
+            auto solver    = make_solver(scheme_);
+            solver.solve(unknown, rhs);
+        }
+
     } // end namespace petsc
 } // end namespace samurai
diff --git a/include/samurai/schemes/block_operator.hpp b/include/samurai/schemes/block_operator.hpp
index be8ca68c0..48e08b557 100644
--- a/include/samurai/schemes/block_operator.hpp
+++ b/include/samurai/schemes/block_operator.hpp
@@ -24,7 +24,12 @@ namespace samurai
             static_assert(n_operators == rows * cols, "The number of operators must correspond to rows*cols.");
         }
 
-        auto& operators() const
+        const auto& operators() const
+        {
+            return m_operators;
+        }
+
+        auto& operators()
         {
             return m_operators;
         }
diff --git a/include/samurai/schemes/fv/explicit_operator_sum.hpp b/include/samurai/schemes/fv/explicit_operator_sum.hpp
index 6daffe6e2..76631a672 100644
--- a/include/samurai/schemes/fv/explicit_operator_sum.hpp
+++ b/include/samurai/schemes/fv/explicit_operator_sum.hpp
@@ -17,24 +17,24 @@ namespace samurai
         using size_type      = typename base_class::size_type;
         using base_class::scheme;
 
-        explicit Explicit(const scheme_t& sum_scheme)
+        explicit Explicit(scheme_t& sum_scheme)
             : base_class(sum_scheme)
         {
         }
 
-        void apply(output_field_t& output_field, input_field_t& input_field) const override
+        void apply(output_field_t& output_field, input_field_t& input_field) override
         {
             for_each(scheme().operators(),
-                     [&](const auto& op)
+                     [&](auto& op)
                      {
                          op.apply(output_field, input_field);
                      });
         }
 
-        void apply(std::size_t d, output_field_t& output_field, input_field_t& input_field) const override
+        void apply(std::size_t d, output_field_t& output_field, input_field_t& input_field) override
         {
             for_each(scheme().operators(),
-                     [&](const auto& op)
+                     [&](auto& op)
                      {
                          op.apply(d, output_field, input_field);
                      });
diff --git a/include/samurai/schemes/fv/scheme_operators.hpp b/include/samurai/schemes/fv/scheme_operators.hpp
index 1fb76a7b7..124a5c49a 100644
--- a/include/samurai/schemes/fv/scheme_operators.hpp
+++ b/include/samurai/schemes/fv/scheme_operators.hpp
@@ -130,19 +130,19 @@ namespace samurai
             return ss.str();
         }
 
-        auto operator()(input_field_t& input_field) const
+        auto operator()(input_field_t& input_field)
         {
             auto explicit_scheme = make_explicit(*this);
             return explicit_scheme.apply_to(input_field);
         }
 
-        auto operator()(std::size_t d, input_field_t& input_field) const
+        auto operator()(std::size_t d, input_field_t& input_field)
         {
             auto explicit_scheme = make_explicit(*this);
             return explicit_scheme.apply_to(d, input_field);
         }
 
-        void apply(output_field_t& output_field, input_field_t& input_field) const
+        void apply(output_field_t& output_field, input_field_t& input_field)
         {
             auto explicit_scheme = make_explicit(*this);
             explicit_scheme.apply(output_field, input_field);
diff --git a/include/samurai/static_algorithm.hpp b/include/samurai/static_algorithm.hpp
index cf418a9e0..2474bcbc1 100644
--- a/include/samurai/static_algorithm.hpp
+++ b/include/samurai/static_algorithm.hpp
@@ -151,7 +151,7 @@ namespace samurai
     template <typename... Ts, typename Func, size_t... Is>
     auto transform_impl(const std::tuple<Ts...>& t, Func&& f, std::index_sequence<Is...>)
     {
-        return std::tuple<std::invoke_result_t<Func, Ts>...>{f(std::get<Is>(t))...};
+        return std::make_tuple(f(std::get<Is>(t))...);
     }
 
     template <typename... Ts, typename Func>
@@ -160,6 +160,18 @@ namespace samurai
         return transform_impl(t, std::forward<Func>(f), std::make_index_sequence<sizeof...(Ts)>{});
     }
 
+    template <typename... Ts, typename Func, size_t... Is>
+    auto transform_impl(std::tuple<Ts...>& t, Func&& f, std::index_sequence<Is...>)
+    {
+        return std::make_tuple(f(std::get<Is>(t))...);
+    }
+
+    template <typename... Ts, typename Func>
+    auto transform(std::tuple<Ts...>& t, Func&& f)
+    {
+        return transform_impl(t, std::forward<Func>(f), std::make_index_sequence<sizeof...(Ts)>{});
+    }
+
     /**
      * Static for loop
      */