diff --git a/tests/reactmicp/binary/unsaturated/acc_carbo.yml b/tests/reactmicp/binary/unsaturated/acc_carbo.yml
index 60f1c78..4a6ef9c 100644
--- a/tests/reactmicp/binary/unsaturated/acc_carbo.yml
+++ b/tests/reactmicp/binary/unsaturated/acc_carbo.yml
@@ -1,196 +1,197 @@
 ---
 plugins:
     to_load:
         - acc_carbo_ups.so
         - acc_carbo_variables.so
 
 logs:
     level:  warning
     output: cerr
 
 configuration_logs:
     output: file
     file: conf_log_acc_carbo
 
 database:
     path: cemdata.yaml
     swap_components:
         - out: HCO3[-]
           in:   CO2
         #- out: H[+]
         #  in: HO[-]
         - out: Al[3+]
           in:  Al(OH)4[-]
     remove_half_cells: true
     remove_sorbeds: true
     list_solids_to_keep: # the solids phases at equilibrium to keep
         - "Portlandite"
         - "CSH,jennite"
         - "CSH,tobermorite"
         - "SiO2(am)"
         - "Calcite"
         - "Al(OH)3(mic)"
         - "C3AH6"
         - "C3ASO_41H5_18"
         - "C3ASO_84H4_32"
         - "C4AH19"
         - "C4AH13"
         - "C2AH7_5"
         - "CAH10"
         - "Monosulfoaluminate"
         - "Tricarboaluminate"
         - "Monocarboaluminate"
         - "Hemicarboaluminate"
         # - "Straetlingite"
         - "Gypsum"
         - "Ettringite"
         # - "Thaumasite"
 
 
 units:
     length: centimeter
     quantity: millimoles
 
 mesh:
     type: uniform_mesh
     uniform_mesh:
         dx: 0.01
         section: 5.0
         nb_nodes: 50
 
 variables:
     type: plugin
     plugin_name: acc_carbo_variables
 
     bc_p_co2: 1000
     bc_p_h2o: 2000
     scaling_sat: 1.0e3/18.3e-3*1e-3
     scaling_aqueous: 1e-6
     scaling_gas: 1e-6
 
 boundary_conditions:
     #gas_nodes: 0
     #fixed_nodes: 0
     constraints:
         default:
             #charge_keeper: "H[+]"
             #charge_keeper: "HO[-]"
             conservation_water: true
             #inert_volume_fraction: 0.1
         others:
             - name: carbonated
               fork_from: default
               nodes: 0
               fixed_fugacity:
                   - label_component: CO2
                     label_gas: CO2(g)
                     amount: 0.01
 
 staggers:
     upscaling:
         type: plugin
         plugin_name: acc_carbo_ups
 
         parameters:
             capillary_a: 37.5479e6
             capillary_b: 2.1684
 
             pv_sat: 3170
 
             exp_r_l_d: 3.00
             exp_i_g_d: 2.74
             exp_r_g_d: 4.20
             exp_r_l_k: 0.4191
 
             k_0: 1.0e-21*1e4
             d_0: 2.3e-13*1e4
             phi_0: 0.25
 
     chemistry:
         type: equilibrium
         equilibrium_options:
             default:
                 residual_tolerance: 1e-10
                 step_tolerance: 1e-12
                 maximum_step_length: 1.0
                 maximum_iterations: 500
                 maximum_step_maximum_iterations: 500
-                non_ideality_tolerance: 1e-10
+                non_ideality_tolerance: 1e-12
                 non_ideality_maximum_iterations: 20
-                cutoff_total_concentration: 1e-8
+                cutoff_total_concentration: 1e-10
                 restart_concentration: -4
                 restart_water_volume_fraction: 0.5
                 under_relaxation: 0.5
                 threshold_cycling_linesearch: 1e-7
                 threshold_condition_check: -1
 
             others:
                 - name: carbonated
                   fork_from: default
                   nodes: 0
-                  residual_tolerance: 1e-12
+                  residual_tolerance: 1e-10
                   cutoff_total_concentration: 1e-14
     transport:
         # general options
         #merge_saturation: true
         #merge_aqueous: true
         merge_saturation: false
         merge_aqueous: false
 
         # dfpm options
         default_options:
             residual_tolerance: 1e-8
             step_tolerance: 1e-12
             absolute_tolerance: 1e-14
             threshold_stationary: 1e-6
-            quasi_newton: 3
+            quasi_newton: 1
             #linesearch: strang
             linesearch: backtracking
-            sparse_solver: LU
+            sparse_solver: QR
 
         saturation_options:
             threshold_stationary: 1e-6
             maximum_step_length: 0.1
             maximum_step_maximum_iterations: 500
             maximum_iterations: 500
 
         aqueous_options:
-            residual_tolerance: 1e-10
+            residual_tolerance: 1e-8
+            step_tolerance: 1e-14
             maximum_step_length: 10.0
             maximum_step_maximum_iterations: 500
             maximum_iterations: 500
 
         gas_options:
             residual_tolerance: 1e-8
-            step_tolerance: 1e-10
-            maximum_step_length: 100.0
+            step_tolerance: 1e-14
+            maximum_step_length: 1000.0
             maximum_step_maximum_iterations: 1000
             maximum_iterations: 1000
 
 reactmicp_options:
-    residual_tolerance: 5e-3
+    residual_tolerance: 1e-3
     good_enough_tolerance: 0.1
-    maximum_iterations: 50
+    maximum_iterations: 100
     step_tolerance: 1e-12
     implicit_upscaling: true
 
 simulation:
     name: acc_carbo
     output_prefix: acc_carbo_
     print_iter_info: true
     output_step: 3600
 
 run:
     run_until: 30*24*3600
 
 timestepper:
-    minimum_dt: 0.001
+    minimum_dt: 1e-1
     maximum_dt: 5.0
-    restart_dt: 0.001
+    restart_dt: 1e-1
     lower_iterations_target: 2.1
 
 
 reactmicp_output:
     database: acc_carbo_db.yaml
     type: hdf5
     file: acc_carbo.hdf5
diff --git a/tests/reactmicp/binary/unsaturated/acc_carbo_ups.cpp b/tests/reactmicp/binary/unsaturated/acc_carbo_ups.cpp
index 1336d33..866c7ae 100644
--- a/tests/reactmicp/binary/unsaturated/acc_carbo_ups.cpp
+++ b/tests/reactmicp/binary/unsaturated/acc_carbo_ups.cpp
@@ -1,338 +1,338 @@
 #include "reactmicp/reactmicp_unsaturated.hpp"
 
 #include "specmicp_common/physics/unsaturated_laws.hpp"
 
 #include "specmicp_common/plugins/plugin_interface.h"
 #include "specmicp_common/plugins/plugin_base.hpp"
 #include "specmicp_common/compat.hpp"
 
 #include "specmicp_common/io/safe_config.hpp"
 
 #include <iostream>
 
 using namespace specmicp;
 using namespace specmicp::reactmicp::solver;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 class SPECMICP_DLL_PUBLIC UpscalingStagger:
         public specmicp::reactmicp::solver::UpscalingStaggerBase
 {
 public:
     UpscalingStagger(
             SimulationData& data,
             io::YAMLConfigHandle&& conf
             );
 
     ~UpscalingStagger() {}
 
     //! \brief Initialize the stagger at the beginning of the computation
     //!
     //! \param var a shared_ptr to the variables
     virtual void initialize(VariablesBase * const var) override;
 
     //! \brief Initialize the stagger at the beginning of an iteration
     //!
     //! This is where the predictor can be saved, the first trivial iteration done, ...
     //!
     //! \param dt the duration of the timestep
     //! \param var a shared_ptr to the variables
     virtual void initialize_timestep(
             scalar_t dt,
             VariablesBase * const var
             ) override;
     //! \param var a shared_ptr to the variables
     virtual StaggerReturnCode restart_timestep(
             VariablesBase * const var) override;
 
     virtual void register_chemistry_stagger(
                 std::shared_ptr<ChemistryStaggerBase> chem_stagger
                 ) override {
         EquilibriumStagger* stagger = static_cast<EquilibriumStagger*>(chem_stagger.get());
         m_chem_opts = stagger->get_options();
     }
 
     // User models
     scalar_t capillary_pressure(index_t node, scalar_t saturation);
 
     scalar_t vapor_pressure(index_t node, scalar_t saturation);
 
     scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation);
 
     scalar_t relative_liquid_permeability(index_t node, scalar_t saturation);
 
 
 //! \brief Solve the equation for the timestep
 //!
 
     scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation);
 
     scalar_t intrinsic_permeability(scalar_t porosity);
 
     scalar_t intrinsic_liquid_diffusivity(scalar_t porosity);
 
     scalar_t resistance_gas_diffusivity(scalar_t porosity);
 
 private:
     scalar_t m_a {37.5479e6};
     scalar_t m_b {2.1684};
     scalar_t m_alpha {7.27};
     scalar_t m_beta {0.440};
     scalar_t m_pv_sat {3170};
     scalar_t m_exp_r_l_d {3.00};
     scalar_t m_exp_i_g_d {2.74};
     scalar_t m_exp_r_g_d {4.20};
     scalar_t m_exp_r_l_k {0.4191};
 
     scalar_t m_k_0 {1e-21*1e4};
     scalar_t m_d_0 {2.3e-13*1e4};
     scalar_t m_phi_0 {0.2};
 
     std::shared_ptr<EquilibriumOptionsVector> m_chem_opts;
     std::shared_ptr<BoundaryConditions> m_bcs;
 };
 
 
 UpscalingStagger::UpscalingStagger(
         SimulationData& data,
         io::YAMLConfigHandle&& conf
         )
 {
     m_bcs = data.bcs; // register boundary conditions
 
     auto param = conf.get_section("parameters");
     m_a = param.get_required_attribute<scalar_t>("capillary_a");
     m_b = param.get_required_attribute<scalar_t>("capillary_b");
 
     m_pv_sat = param.get_required_attribute<scalar_t>("pv_sat");
 
     m_exp_r_l_d = param.get_required_attribute<scalar_t>("exp_r_l_d");
     m_exp_i_g_d = param.get_required_attribute<scalar_t>("exp_i_g_d");
     m_exp_r_g_d = param.get_required_attribute<scalar_t>("exp_r_g_d");
     m_exp_r_l_k = param.get_required_attribute<scalar_t>("exp_r_l_k");
 
     m_k_0 = param.get_required_attribute<scalar_t>("k_0");
     m_d_0 = param.get_required_attribute<scalar_t>("d_0");
     m_phi_0 = param.get_required_attribute<scalar_t>("phi_0");
 }
 
 class SPECMICP_DLL_PUBLIC AccCarboUpscalingFactory: public UpscalingStaggerFactory
 {
 public:
     static AccCarboUpscalingFactory* get() {
         return new AccCarboUpscalingFactory();
     }
 
     virtual std::shared_ptr<UpscalingStaggerBase> get_upscaling_stagger(
             SimulationData& data,
             io::YAMLConfigHandle&& configuration
             ) override
     {
         //auto bcs = data.bcs;
         //bcs->set_fixed_flux_liquid_dof(1, 0, 2.2e-8*1e4);
 
         return std::make_shared<UpscalingStagger>(data, std::move(configuration));
     }
 };
 
 
 class SPECMICP_DLL_PUBLIC AccCarboUpscalingPlugin:
         public specmicp::plugins::PluginBase
 {
 public:
     AccCarboUpscalingPlugin():
         PluginBase()
     {
         set_api_version(0, 1, 0);
     }
 
     bool initialize(const specmicp::plugins::PluginManagerServices& services) override
     {
         auto retcode = services.register_object(
                            "reactmicp/unsaturated/upscaling", "acc_carbo_ups",
                            &AccCarboUpscalingFactory::get
                            );
         if (not retcode) {
             return false;
         }
         return true;
     }
 
 };
 
 SPECMICP_PLUGIN(AccCarboUpscalingPlugin)
 
 
 // implementation
 
 void UpscalingStagger::initialize(VariablesBase * const var)
 {
     UnsaturatedVariables * const true_var = static_cast<UnsaturatedVariables * const>(var);
+    m_bcs->add_fixed_flux_gas_dof(0, 0);
 
     // set relative models
 
     true_var->set_capillary_pressure_model(
                 [this](index_t x, scalar_t sat){
                     return this->capillary_pressure(x, sat);
                 });
     true_var->set_relative_liquid_permeability_model(
                 [this](index_t x, scalar_t sat){
                     return this->relative_liquid_permeability(x, sat);
                 });
     true_var->set_relative_liquid_diffusivity_model(
                 [this](index_t x, scalar_t sat){
                     return this->relative_liquid_diffusivity(x, sat);
                 });
     true_var->set_relative_gas_diffusivity_model(
                 [this](index_t x, scalar_t sat){
                     return this->relative_gas_diffusivity(x, sat);
                 });
     true_var->set_vapor_pressure_model(
                 [this](index_t x, scalar_t sat){
                     return this->vapor_pressure(x, sat);
                 });
 
     // initialization
 
     SecondaryTransientVariable& porosity = true_var->get_porosity();
     SecondaryVariable& permeability = true_var->get_liquid_permeability();
     SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
     SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
 
     //true_var->set_relative_variables();
     //gas_diffusivity(0) = resistance_gas_diffusivity(1.0);
     //true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0);
 
     for (index_t node=0; node<true_var->get_mesh()->nb_nodes(); ++node)
     {
         scalar_t phi = porosity(node);
 
         permeability(node) = intrinsic_permeability(phi);
         liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
         gas_diffusivity(node) = resistance_gas_diffusivity(phi);
 
         true_var->set_relative_variables(node);
     }
 
    // true_var->set_relative_variables(0);
 }
 
 scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation)
 {
     if (saturation == 0.0) return 0.0;
 
     return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b);
 }
 
 scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation)
 {
     scalar_t pv;
     if (saturation == 0.0) pv = 0.0;
     else if (saturation >= 1.0) pv = m_pv_sat;
     else pv = m_pv_sat*laws::invert_kelvin_equation(
                 capillary_pressure(node, saturation));
     //if (node == 1)
         //std::cout << "saturation " << saturation << " - pv " << pv << "\n";
     return pv;
 }
 
 scalar_t UpscalingStagger::relative_liquid_diffusivity(
         index_t node,
         scalar_t saturation
         )
 {
     if (saturation == 0.0) return 0.0;
     return std::pow(saturation, m_exp_r_l_d);
 }
 
 scalar_t UpscalingStagger::relative_liquid_permeability(
         index_t node,
         scalar_t saturation
         )
 {
     if (saturation == 0.0) return 0.0;
 
     return laws::relative_liquid_permeability_Mualem(saturation, 1.0/m_b);
 }
 
 scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation)
 {
     if (saturation == 0.0) return 1.0;
 
     return std::pow(1.0 - saturation, m_exp_r_g_d);
 }
 
 scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity)
 {
     scalar_t tmp_1 = porosity / m_phi_0;
     scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity);
     tmp_1 = std::pow(tmp_1, 3);
     tmp_2 = std::pow(tmp_2, 2);
     return m_k_0 * tmp_1 * tmp_2;
 }
 
 scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity)
 {
     return m_d_0 * std::exp(-9.95*porosity);
 }
 
 scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity)
 {
     return std::pow(porosity, m_exp_i_g_d);
 }
 
 void UpscalingStagger::initialize_timestep(
         scalar_t dt,
         VariablesBase * const var
         )
 {
 
     UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
 
-    m_bcs->add_fixed_flux_gas_dof(0, 0);
 
     scalar_t flux = - 1.22e-7 * true_var->get_rt() * (
                         true_var->get_pressure_main_variables(0).variable(0) - 2000);
 
     m_bcs->set_fixed_flux_gas_dof(0, 0, flux);
 
 //    const auto& the_mesh = true_var->get_mesh();
 //    const auto& raw_data = true_var->get_database();
 
 //    auto& liq_sat = true_var->get_liquid_saturation();
 //    auto& ca_conc = true_var->get_solid_concentration(raw_data->get_id_component("Ca[2+]"));
 //    auto& co2_conc = true_var->get_solid_concentration(raw_data->get_id_component("CO2"));
 
 
 //    for (index_t node: the_mesh->range_nodes()) {
 //        if (not liq_sat(node) > 0) continue;
 
 //        if (co2_conc(node) >= ca_conc(node)) {
 //            WARNING << "Node " << node << "is carbonated";
 //            m_chem_opts->set(node, "carbonated");
 //        }
 //    }
 }
 
 StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBase * const var)
 {
     UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var);
 
     SecondaryTransientVariable& porosity = true_var->get_porosity();
     SecondaryVariable& permeability = true_var->get_liquid_permeability();
     SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity();
     SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity();
 
 
     for (index_t node=0; node<true_var->get_mesh()->nb_nodes(); ++node)
     {
         scalar_t phi = porosity(node);
 
         permeability(node) = intrinsic_permeability(phi);
         liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi);
         gas_diffusivity(node) = resistance_gas_diffusivity(phi);
 
         true_var->set_relative_variables(node);
     }
 
     return StaggerReturnCode::ResidualMinimized;
 }
diff --git a/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp b/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp
index c7119ca..2741a9f 100644
--- a/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp
+++ b/tests/reactmicp/binary/unsaturated/acc_carbo_variables.cpp
@@ -1,301 +1,301 @@
 #include "reactmicp/reactmicp_unsaturated.hpp"
 
 #include "specmicp_common/physics/unsaturated_laws.hpp"
 
 #include "specmicp_common/plugins/plugin_interface.h"
 #include "specmicp_common/plugins/plugin_base.hpp"
 #include "specmicp_common/compat.hpp"
 
 #include "specmicp_common/io/safe_config.hpp"
 
 #include "specmicp/problem_solver/reactant_box.hpp"
 #include "specmicp/problem_solver/smart_solver.hpp"
 
 #include <iostream>
 #include <Eigen/QR>
 
 using namespace specmicp;
 using namespace specmicp::reactmicp::solver;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 
 AdimensionalSystemSolution get_initial_condition(
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         );
 
 AdimensionalSystemSolution carbonate_initial_condition(
         AdimensionalSystemSolution& init_cond,
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         );
 
 
 scalar_t vapor_pressure(index_t node, scalar_t saturation);
 
 class SPECMICP_DLL_PUBLIC AccCarboVariablesFactory: public InitializeVariables
 {
 public:
     static AccCarboVariablesFactory* get() {
         return new AccCarboVariablesFactory();
     }
 
     virtual std::shared_ptr<UnsaturatedVariables> initialize_variables(
             std::shared_ptr<database::DataContainer> the_database,
             std::shared_ptr<mesh::Mesh1D> the_mesh,
             const units::UnitsSet& the_units,
             io::YAMLConfigHandle& configuration
             ) override
     {
 
         AdimensionalSystemSolution init = get_initial_condition(
                     the_database,
                     the_units);
         AdimensionalSystemSolutionExtractor extr(init, the_database, the_units);
 
 
         AdimensionalSystemSolution init_carb = carbonate_initial_condition(
                     init,
                     the_database,
                     the_units);
         AdimensionalSystemSolutionExtractor extr_carb(init_carb, the_database, the_units);
 
 
 
 
         std::cout<< "Hop : " << init.main_variables(0) << std::endl;
         std::cout<< "Hop : " << extr.saturation_water() << std::endl;
         std::cout<< "Hop : " << extr.porosity() << std::endl;
 
         //throw std::runtime_error("that's enough");
 
         index_t id_co2 = the_database->get_id_component("CO2");
 
         VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, the_units);
 
         variable_interface.set_binary_gas_diffusivity(0,      0.240);
         variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);
 
         variable_interface.set_vapor_pressure_model(&vapor_pressure);
 
         variable_interface.initialize_variables(0, extr_carb);
         for (index_t node=1; node<the_mesh->nb_nodes(); ++node)
         {
             variable_interface.initialize_variables(node, extr);
         }
 
         //variable_interface.set_porosity(0, 1.0);
         //variable_interface.set_liquid_saturation(0, 0.0);
 
         const auto p_co2 = configuration.get_required_attribute<scalar_t>("bc_p_co2");
         variable_interface.set_partial_pressure(id_co2, 0, p_co2);
         const auto p_h2o = configuration.get_required_attribute<scalar_t>("bc_p_h2o");
         variable_interface.set_partial_pressure(0, 0, p_h2o);
 
         std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
 
         const auto scal_sat = configuration.get_required_attribute<scalar_t>(
                     "scaling_sat");
         //vars->set_aqueous_scaling(0, scal_sat);
-        vars->set_aqueous_scaling(0, extr.total_aqueous_concentration(0));
+        vars->set_aqueous_scaling(0, extr_carb.total_aqueous_concentration(0));
         const auto scal_aqueous = configuration.get_required_attribute<scalar_t>(
                     "scaling_aqueous");
         for (auto aqc: the_database->range_aqueous_component())
         {
 
             //vars->set_aqueous_scaling(aqc, scal_aqueous);
-            vars->set_aqueous_scaling(aqc, std::max(1e-3, std::abs(extr.total_aqueous_concentration(aqc))));
+            vars->set_aqueous_scaling(aqc, std::max(1e-3, std::abs(extr_carb.total_aqueous_concentration(aqc))));
         }
 
         const auto scal_gas = configuration.get_required_attribute<scalar_t>(
                     "scaling_gas");
         for (auto akk: the_database->range_component())
         {
             vars->set_gaseous_scaling(akk, scal_gas);
             return vars;
         }
     }
 };
 
 
 class SPECMICP_DLL_PUBLIC AccCarboVariablesPlugin: public specmicp::plugins::PluginBase
 {
 public:
     AccCarboVariablesPlugin():
         PluginBase()
     {
         set_api_version(0, 1, 0);
     }
 
     bool initialize(const specmicp::plugins::PluginManagerServices& services) override
     {
         auto retcode = services.register_object(
                            "reactmicp/unsaturated/variables", "acc_carbo_variables",
                            &AccCarboVariablesFactory::get
                            );
         if (not retcode) {
             return false;
         }
 
         return true;
     }
 
 };
 
 SPECMICP_PLUGIN(AccCarboVariablesPlugin)
 
 
 // implementation
 
 scalar_t capillary_pressure(index_t node, scalar_t saturation)
 {
     if (saturation == 0.0) return 0.0;
 
     return laws::capillary_pressure_BaroghelBouny(saturation, 37.5479e6, 2.1684);
 }
 
 scalar_t vapor_pressure(index_t node, scalar_t saturation)
 {
     scalar_t pv;
 
     if (saturation == 0.0) pv = 0.0;
     else if (saturation >= 1.0) pv = 3170;
     else pv = 3170*laws::invert_kelvin_equation(
                 capillary_pressure(node, saturation));
 
     std::cout << "saturation " << saturation << " - pv " << pv << std::endl;
     return pv;
 }
 
 AdimensionalSystemSolution get_initial_condition(
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         )
 {
 
     ReactantBox reactant_box(the_database, units_set);
     reactant_box.set_solution(0.57, "kg/L");
     scalar_t factor = 1.0/100*1.2;
     reactant_box.add_solid_phase("CaO(ox)",   factor*67.87, "kg/L");
     reactant_box.add_solid_phase("SiO2(ox)",  factor*22.66, "kg/L");
     reactant_box.add_solid_phase("Al2O3(ox)", factor* 5.19, "kg/L");
     reactant_box.add_solid_phase("SO3(ox)",   factor* 2.96, "kg/L");
     reactant_box.add_solid_phase("Calcite", factor*1.0, "kg/L");
     //reactant_box.add_component("CO2");
 
     Vector tot_conc = reactant_box.get_total_concentration(true);// also simplify the db
 
     std::cout << "\n ------ \n " << tot_conc << "\n ------ \n" << std::endl;
 
     AdimensionalSystemConstraints constraints;
     constraints.set_total_concentrations(tot_conc);
     //constraints.set_saturated_system();
     constraints.set_inert_volume_fraction(0.0);
 
     constraints.set_water_partial_pressure_model(
                 std::bind(&vapor_pressure,
                           1,
                           std::placeholders::_1)
                 );
     //constraints.charge_keeper = the_database->get_id_component_from_element("H");
 
     //constraints.set_inert_volume_fraction(0.1);
 
     AdimensionalSystemSolverOptions opts;
     opts.units_set = units_set;
     opts.system_options.non_ideality_tolerance = 1e-14;
     opts.system_options.restart_water_volume_fraction = 0.8;
     opts.system_options.cutoff_total_concentration = 1e-14;
     opts.solver_options.set_tolerance(1e-10, 1e-12);
     //micpsolver::MiCPSolverOptions& solver_options = opts.solver_options;
     //solver_options.maxstep = 10.0;
     //solver_options.maxiter_maxstep = 200;
     //solver_options.max_iter = 200;
     //solver_options.set_tolerance(1e-6, 1e-12);
 
 
     SmartAdimSolver smart_solver(the_database, constraints, opts);
     smart_solver.solve();
     if (not smart_solver.is_solved())
     {
         throw std::runtime_error("Error : failed to solve initial conditions");
     }
     AdimensionalSystemSolution sol = smart_solver.get_solution();
 
 
 //    std::cout << "normal solver " << std::endl;
 //    AdimensionalSystemSolver solver(the_database, constraints, opts);
 //    Vector x;
 //    solver.set_options(opts);
 //    solver.initialise_variables(x, 0.5, -6);
 //    micpsolver::MiCPPerformance perf = solver.solve(x, false);
 //    std::cout << perf.current_residual << " - "
 //              << perf.current_update << " - "
 //              << perf.nb_iterations
 //              << std::endl;
 //    if (perf.return_code <= micpsolver::MiCPSolverReturnCode::Success) {
 //        throw std::runtime_error("Error : failed to solve initial conditions");
 //    }
 //    AdimensionalSystemSolution sol = solver.get_raw_solution(x);
 
 
    std::cout << sol.main_variables << std::endl;
    AdimensionalSystemSolutionExtractor extr(sol, the_database, units_set);
 
     for (index_t idc: the_database->range_component())
     {
         std::cout << the_database->get_label_component(idc) << " : "
                   << extr.total_concentration(idc) << " - "
                   << tot_conc(idc) << " - "
                   << (extr.total_concentration(idc) - tot_conc(idc))/tot_conc(idc) << std::endl;
     }
 
     //throw std::runtime_error("That's enough");
 
     return sol;
 }
 
 
 AdimensionalSystemSolution carbonate_initial_condition(
         AdimensionalSystemSolution& init_cond,
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         )
 {
     AdimensionalSystemSolutionExtractor extr(init_cond, the_database, units_set);
     scalar_t init_porosity = extr.porosity();
     AdimensionalSystemConstraints constraints;
 
     Vector tot_conc = extr.total_concentrations();
     index_t id_ca = the_database->get_id_component_from_element("Ca");
     index_t id_co2 = the_database->get_id_component_from_element("C");
 
     tot_conc(id_co2) = tot_conc(id_ca);
 
     constraints.set_total_concentrations(tot_conc);
     constraints.set_saturated_system();
 
     AdimensionalSystemSolverOptions opts;
     opts.units_set = units_set;
     opts.system_options.non_ideality_tolerance = 1e-14;
     opts.system_options.restart_water_volume_fraction = 0.8;
     opts.system_options.cutoff_total_concentration = 1e-14;
     opts.solver_options.set_tolerance(1e-10, 1e-12);
 
     AdimensionalSystemSolver solver(the_database, constraints, init_cond, opts);
 
     Vector x = init_cond.main_variables;
     auto perf = solver.solve(x, false);
     if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
     {
         throw std::runtime_error("Failed to solve carbonated initial condition");
     }
     AdimensionalSystemSolution solution = solver.get_raw_solution(x);
 
     // set the same saturation for both cases
     AdimensionalSystemSolutionExtractor ext2r(solution, the_database, units_set);
     scalar_t porosity = ext2r.porosity();
     solution.main_variables(0) = porosity/init_porosity * init_cond.main_variables(0);
 
     return solution;
 }