diff --git a/examples/reactmicp/unsaturated/acc_carbo.cpp b/examples/reactmicp/unsaturated/acc_carbo.cpp
index d2a57be..995e114 100644
--- a/examples/reactmicp/unsaturated/acc_carbo.cpp
+++ b/examples/reactmicp/unsaturated/acc_carbo.cpp
@@ -1,527 +1,527 @@
 /*------------------------------------------------------------------------------
 
 Copyright (c)2015 F. Georget <fabieng@princeton.edu>, Princeton University
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
 1. Redistributions of source code must retain the above copyright notice, this
 list of conditions and the following disclaimer.
 
 2. Redistributions in binary form must reproduce the above copyright notice,
 this list of conditions and the following disclaimer in the documentation and/or
 other materials provided with the distribution.
 
 3. Neither the name of the copyright holder nor the names of its contributors
 may be used to endorse or promote products derived from this software without
 specific prior written permission.
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 ------------------------------------------------------------------------------*/
 
 //! \file acc_carbo.cpp
 //! \brief Accelerated carbonation example
 //!
 //! This is an example to show the features of the unsaturated
 //! reactive transport solver
 //! It is a simple study of the accelerated carbonation of cement paste
 //! (with drying)
 //!
 //!
 
 #include "reactmicp_unsaturated.hpp"
 
 #include "reactmicp/systems/unsaturated/variables.hpp"
 
 #include <H5Cpp.h>
 #include "reactmicp/io/hdf5_unsaturated.hpp"
 #include "utils/io/specmicp_hdf5.hpp"
 #include "dfpm/io/hdf5_mesh.hpp"
 #include "database/io/hdf5_database.hpp"
 
 #include "dfpmsolver/parabolic_structs.hpp"
 
 #include "physics/unsaturated_laws.hpp"
 
 #include "utils/cli/parser.hpp"
 
 #include <iostream>
 
 using namespace specmicp;
 using namespace specmicp::reactmicp;
 using namespace specmicp::reactmicp::systems::unsaturated;
 
 using VariablesBase = solver::VariablesBase;
 using VariablesBasePtr = solver::VariablesBasePtr;
 using StaggerReturnCode = solver::StaggerReturnCode;
 
 
 class UpscalingStagger;
 
 database::RawDatabasePtr get_database();
 
 AdimensionalSystemSolution get_initial_condition(
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         );
 
 // Upscaling stagger -> user models
 
 class UpscalingStagger: public solver::UpscalingStaggerBase
 {
 public:
     UpscalingStagger() {}
 
     ~UpscalingStagger() {}
 
     //! \brief Initialize the stagger at the beginning of the computation
     //!
     //! \param var a shared_ptr to the variables
     void initialize(VariablesBasePtr var);
 
     //! \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
     void initialize_timestep(scalar_t dt, VariablesBasePtr var);
 
     //! \brief Solve the equation for the timestep
     //!
     //! \param var a shared_ptr to the variables
     StaggerReturnCode restart_timestep(VariablesBasePtr var);
 
     // 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);
 
     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};
 };
 
 // implementations
 
 int main(int argc, char* argv[])
 {
     cli::CommandLineParser cli_parser;
     cli_parser.add_option('n', "min_dt", 0.01, "minimum timestep (s)");
     cli_parser.add_option('m', "max_dt", 10.0, "maximum timestep (s)");
     cli_parser.add_option('d', "dx", 0.05, "space step (cm)");
-    cli_parser.add_option('r', "run_until", 30.0*24.0*3600.0, "run simulation until <value> s");
+    cli_parser.add_option('r', "run_until", 10.0*24.0*3600.0, "run simulation until <value> s");
     cli_parser.add_option('p', "p_h20", 1800.0, "vapor pressure at the boundary");
 
     cli_parser.set_help_message("Atmospheric carbonation");
 
     cli_parser.parse(argc, argv);
 
     init_logger(&std::cout, logger::Warning);
 
 
 
     index_t nb_nodes = 50;
     scalar_t dx = cli_parser.get_option<cli::ValueType::floating>("dx"); // cm
     scalar_t cross_section = 5; // cm^2
 
     units::UnitsSet units_set;
     units_set.length = units::LengthUnit::centimeter;
 
     specmicp::RawDatabasePtr the_database = get_database();
 
     AdimensionalSystemSolution init = get_initial_condition(
                 the_database,
                 units_set);
 
     AdimensionalSystemSolutionExtractor extr(init, the_database, units_set);
 
     std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl;
 
     //return EXIT_SUCCESS;
 
     mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section});
 
     index_t id_co2 = the_database->get_id_component("CO2");
 
     VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set);
 
     variable_interface.set_binary_gas_diffusivity(0, 0.240);
     variable_interface.set_binary_gas_diffusivity(id_co2, 0.160);
 
     std::shared_ptr<UpscalingStagger> upscaling_stagger =
             std::make_shared<UpscalingStagger>();
 
     auto* call_ptr = upscaling_stagger.get();
     variable_interface.set_vapor_pressure_model(
                 [call_ptr](index_t node, scalar_t sat) -> scalar_t {
                     return call_ptr->vapor_pressure(node, sat);
                 });
 
     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);
 
     variable_interface.set_partial_pressure(id_co2, 0, 100.0);
     variable_interface.set_partial_pressure(0, 0,
             cli_parser.get_option<cli::ValueType::floating>("p_h20"));
 
     std::shared_ptr<UnsaturatedVariables> vars = variable_interface.get_variables();
 
 
     vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6);
     for (auto aqc: the_database->range_aqueous_component())
     {
         vars->set_aqueous_scaling(aqc, 100e-6);
     }
     //vars->set_gaseous_scaling(0, 1.0/vars->get_rt());
     //vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt());
 
     std::cout << "Youpla : " << vars->get_rt() << std::endl;
 
     std::shared_ptr<EquilibriumConstraints>  eq_constraints = equilibrium_constraints(the_mesh);
     eq_constraints->add_fixed_node(0);
     eq_constraints->set_constraint_for_all(eq_constraints->new_constraints());
     eq_constraints->get_options().units_set = units_set;
 
     micpsolver::MiCPSolverOptions& copts = eq_constraints->get_options().solver_options;
     copts.set_maximum_iterations(200);
     copts.set_maximum_step_length(10, 200);
     copts.set_tolerance(1e-8, 1e-14);
     eq_constraints->get_options().system_options.non_ideality_tolerance = 1e-12;
 
     std::shared_ptr<EquilibriumStagger> chemistry_stagger =
             equilibrium_stagger(vars, eq_constraints);
 
 
 
     upscaling_stagger->initialize(vars);
     vars->set_relative_variables(0);
     chemistry_stagger->initialize(vars);
 
     TransportConstraints transport_constraints;
     transport_constraints.add_gas_node(0);
     transport_constraints.add_fixed_node(0);
 
     std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
             unsaturated_transport_stagger(vars, transport_constraints, true, true);
 
     transport_stagger->initialize(vars);
 
     ((*eq_constraints)[0]).set_water_partial_pressure_model(
                 std::bind(&UpscalingStagger::vapor_pressure,
                           upscaling_stagger.get(),
                           1,
                           std::placeholders::_1)
                 );
 
     auto check = variable_interface.check_variables();
     if (check >= VariablesValidity::error) {
         throw std::runtime_error("Invalid variables");
     }
 
     solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger);
     solver::ReactiveTransportOptions& opts = react_solver.get_options();
     opts.residuals_tolerance = 1e-2;
     opts.good_enough_tolerance = 0.99;
-    opts.maximum_iterations = 500;
-    transport_stagger->get_saturation_options()->maximum_iterations = 5000;
+    opts.maximum_iterations = 100;
+    transport_stagger->get_saturation_options()->maximum_iterations = 500;
     transport_stagger->get_saturation_options()->step_tolerance = 1e-14;
     transport_stagger->get_saturation_options()->residuals_tolerance = 1e-8;
     //transport_stagger->get_gas_options(3)->maximum_iterations = 5000;
     opts.step_tolerance = 1e-10;
 
     io::HDF5File save_file("out_acc_carbo.hdf5", io::HDF5_OpenMode::CreateTruncate);
     io::save_database_labels(save_file, "database", "/", the_database);
     io::save_mesh_coordinates(save_file, "mesh", "/", the_mesh);
     io::UnsaturatedHDF5Saver saver(vars);
 
     saver.save_timepoint(save_file, "0", "/");
     save_file.close();
 
     solver::SimulationInformation simul_info("acc_carbo", 900);
 
     auto out_pol = [&saver](scalar_t timestep, VariablesBasePtr vars) {
         io::HDF5File file("out_acc_carbo.hdf5", io::HDF5_OpenMode::OpenReadWrite);
         saver.save_timepoint(file, std::to_string(timestep), "/");
     };
 
     //solver::ReactiveTransportReturnCode retcode =  react_solver.solve_timestep(0.01, vars);
     //retcode =  react_solver.solve_timestep(1.0, vars);
 
     //std::cout << (int) retcode << std::endl;
 //    retcode =  react_solver.solve_timestep(5.0, vars);
 //    std::cout << (int) retcode << std::endl;
 //    retcode =  react_solver.solve_timestep(10.0, vars);
 //    std::cout << (int) retcode << std::endl;
 //    retcode =  react_solver.solve_timestep(10.0, vars);
 //    std::cout << (int) retcode << std::endl;
 
     solver::ReactiveTransportRunner runner(
                 react_solver,
                 cli_parser.get_option<cli::ValueType::floating>("min_dt"),
                 cli_parser.get_option<cli::ValueType::floating>("max_dt"),
                 simul_info);
     runner.set_output_policy(out_pol);
 
 
     runner.run_until(cli_parser.get_option<cli::ValueType::floating>("run_until"),
                      vars);
 
     std::cout << vars->get_liquid_saturation().variable << std::endl;
     std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl;
 
     std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl;
 }
 
 void UpscalingStagger::initialize(VariablesBasePtr var)
 {
     UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var.get());
 
     // 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=1; 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);
     }
 }
 
 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));
     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/m_b);
+    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);
+    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, VariablesBasePtr var)
 {
 }
 
 StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBasePtr var)
 {
     UnsaturatedVariables* true_var = static_cast<UnsaturatedVariables*>(var.get());
 
     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=1; 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;
 }
 
 specmicp::RawDatabasePtr get_database()
 {
     specmicp::database::Database thedatabase(CEMDATA_PATH);
     std::map<std::string, std::string> swap = {{"HCO3[-]", "CO2"},};
     thedatabase.swap_components(swap);
     database::RawDatabasePtr raw_data = thedatabase.get_database();
     return raw_data;
 }
 
 AdimensionalSystemSolution get_initial_condition(
         database::RawDatabasePtr the_database,
         const units::UnitsSet& units_set
         )
 {
     Formulation initial_formulation;
     
     scalar_t scaling = 0.9e-3;
     initial_formulation.set_mass_solution(0.6*scaling);
     initial_formulation.add_mineral("C3S", 5.5*scaling);
     initial_formulation.add_mineral("C2S", 1.25*scaling);
     initial_formulation.add_mineral("Calcite", 0.05*scaling);
-    
     Dissolver dissolver(the_database);
 
     UpscalingStagger upsstag_init; // just for vapor pressure model
 
     AdimensionalSystemConstraints constraints;
     constraints.set_total_concentrations(
                 dissolver.dissolve(initial_formulation, true)
                 );
+    //constraints.set_inert_volume_fraction(0.1);
 
     AdimensionalSystemSolver the_solver(the_database, constraints);
     the_solver.get_options().units_set = units_set;
     the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
     micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options;
     //solver_options.maxstep = 10.0;
     solver_options->maxiter_maxstep = 200;
     solver_options->set_tolerance(1e-9);
 
     Vector x;
     the_solver.initialise_variables(x, 0.8, -4);
 
     micpsolver::MiCPPerformance perf = the_solver.solve(x);
 
     if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
     {
         throw std::runtime_error("Error : failed to solve initial conditions");
     }
 
     constraints.set_water_partial_pressure_model(
                 std::bind(&UpscalingStagger::vapor_pressure,
                           &upsstag_init,
                           1,
                           std::placeholders::_1)
                 );
     the_solver = AdimensionalSystemSolver (the_database, constraints);
     the_solver.get_options().units_set = units_set;
     the_solver.get_options().system_options.non_ideality_tolerance = 1e-12;
     solver_options = &the_solver.get_options().solver_options;
     //solver_options.maxstep = 10.0;
     solver_options->maxiter_maxstep = 200;
     solver_options->set_tolerance(1e-9);
     perf = the_solver.solve(x);
 
     if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
     {
         throw std::runtime_error("Error : failed to solve initial conditions");
     }
 
     return the_solver.get_raw_solution(x);
 }