diff --git a/cmake/SpecMiCPOptions.cmake b/cmake/SpecMiCPOptions.cmake
index e5cea45..02b4f3a 100644
--- a/cmake/SpecMiCPOptions.cmake
+++ b/cmake/SpecMiCPOptions.cmake
@@ -1,33 +1,34 @@
 
 # For an explanation of the options see the INSTALL file
 option( SPECMICP_USE_OPENMP          "Use OpenMP for parallelisation"   ON  )
 option( SPECMICP_NO_DEBUG            "Disable SpecMiCP assert"          OFF )
 option( SPECMICP_BUILD_STATIC        "Build static libraries"           OFF )
 option( SPECMICP_BUILD_EXAMPLE       "Build the examples"               ON  )
 option( SPECMICP_BENCHMARK           "Build benchmark"                  OFF )
 option( SPECMICP_TEST                "Enable testing"                   ON  )
 option( SPECMICP_BINARIES_USE_STATIC "Executables use static libraries" OFF )
 
 option( SPECMICP_OPTIMIZE_FOR_NATIVE "Optimize modules for native machine" ON )
 
 #  Linear algebra
 option( SPECMICP_USE_BLAS            "Eigen uses BLAS/LAPACK"           OFF )
 
 # PGO sequence
 option( SPECMICP_PROFILE_GENERATE "Generate profile for PGO optimization" OFF )
 option( SPECMICP_PROFILE_USE      "Use profile for PGO optimization"      OFF )
 
 # LTO optimization
 option( SPECMICP_LD_GOLD "Use GNU gold linker"        ON  )
 option( SPECMICP_LTO     "Use link time optimization" OFF )
 option( SPECMICP_FAT_LTO "Use link time optimization with fat objects" ON )
 
 # the following is only a debug options for developpers
 
 # This options turns on the finite difference jacobian in specmicp system
 option( SPECMICP_DEBUG_EQUATION_FD_JACOBIAN "Use a finite difference jacobian" OFF )
 
 # ReactMiCP systems
 option( REACTMICP_ENABLE_SYSTEM_SATURATED "Enable Saturated ReactMiCP system" ON )
 option( REACTMICP_ENABLE_SYSTEM_UNSATURATED "Enable Unsaturated ReactMiCP system" ON )
 option( REACTMICP_ENABLE_SYSTEM_CHLORIDE "Enable Chloride ReactMiCP system" ON )
+option( REACTMICP_ENABLE_SYSTEM_SINGLE "Enable Single ReactMiCP system" ON )
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index c32fa36..cec5116 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,252 +1,266 @@
 ###### Installation #######################
 
 set(EXAMPLES_INSTALL_DIR ${SHARE_INSTALL_DIR}/examples)
 
 # SpecMiCP
 # ========
 
 set (SPECMICP_EXAMPLES
     specmicp/adimensional/thermocarbo.cpp
     specmicp/adimensional/equilibrium_curve.cpp
     specmicp/adimensional/momas_thermo.cpp
     specmicp/adimensional/carbofe.cpp
 )
 install(FILES ${SPECMICP_EXAMPLES}
     DESTINATION ${EXAMPLES_INSTALL_DIR}/specmicp
 )
 
 # ReactMiCP
 # =========
 
 set (REACTMICP_SATURATED_EXAMPLES
     reactmicp/equilibrium_curve.cpp
 
     reactmicp/saturated_react/react_leaching.cpp
     reactmicp/saturated_react/carbonation.cpp
     reactmicp/saturated_react/carbonation.yaml # the configuration file
     reactmicp/saturated_react/carbonationfe.cpp
     reactmicp/saturated_react/leaching_kinetic.cpp
 
     reactmicp/unsaturated/drying.cpp
     reactmicp/unsaturated/acc_carbo.cpp
 
     reactmicp/chloride/migration.cpp
     reactmicp/chloride/benchmark_rasouli.cpp
 )
 
 install(FILES ${REACTMICP_SATURATED_EXAMPLES}
     DESTINATION ${EXAMPLES_INSTALL_DIR}/reactmicp/saturated
 )
 
 file(COPY reactmicp/saturated_react/carbonation.yaml
      DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
 )
 
 ###### Build (optional) ####################
 
 if (NOT SPECMICP_BUILD_EXAMPLE)
     set(EXAMPLE_IS_OPTIONAL EXCLUDE_FROM_ALL)
 else()
     set(EXAMPLE_IS_OPTIONAL )
 endif()
 
 # SpecMiCP
 # ========
 
 #            thermocarbo
 
 add_executable(ex_adim_thermocarbo
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/thermocarbo.cpp
 )
 target_link_libraries(ex_adim_thermocarbo ${SPECMICP_LIBS})
 
 
 #            carbofe
 
 add_executable(ex_adim_carbofe
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/carbofe.cpp
 )
 target_link_libraries(ex_adim_carbofe ${SPECMICP_LIBS})
 
 #           Momas thermodynamic example
 
 add_executable(ex_adim_momas
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/momas_thermo.cpp
 )
 target_link_libraries(ex_adim_momas ${SPECMICP_LIBS})
 
 
 #            equilibrium curve
 
 add_executable(ex_adim_equilibriumcurve
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/equilibrium_curve.cpp
 )
 target_link_libraries(ex_adim_equilibriumcurve ${SPECMICP_LIBS})
 
 
 
 # ReactMiCP
 # =========
 
 # EquilibriumCurve
 
 add_executable(ex_reactmicp_equilibriumcurve
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/equilibrium_curve.cpp
 )
 target_link_libraries(ex_reactmicp_equilibriumcurve ${REACTMICP_LIBS})
 
 # a simple macro to find if the examples from a reactmicp systems
 # must be included
 #
 # Check the option REACTMICP_ENABLE_SYSTEM_XXXX
 macro(set_exclude_if_no_reactmicp_system exclude_var reactmicp_option)
     if (NOT ${reactmicp_option} )
         set( ${exclude_var} EXCLUDE_FROM_ALL )
     else ()
         set( ${exclude_var} ${EXAMPLE_IS_OPTIONAL})
     endif ()
 endmacro(set_exclude_if_no_reactmicp_system)
 
 
 # Saturated system
 # ----------------
 
 set_exclude_if_no_reactmicp_system(
     DONT_COMPILE_SATURATED_EXAMPLE
     REACTMICP_ENABLE_SYSTEM_SATURATED
 )
 
 # Leaching
 
 add_executable(ex_reactmicp_leaching
     ${DONT_COMPILE_SATURATED_EXAMPLE}
     reactmicp/saturated_react/react_leaching.cpp
 )
 target_link_libraries(ex_reactmicp_leaching ${REACTMICP_LIBS})
 
 # Momas benchmark
 
 add_executable(ex_reactmicp_momas
     ${DONT_COMPILE_SATURATED_EXAMPLE}
     reactmicp/saturated_react/momas_benchmark.cpp
 )
 target_link_libraries(ex_reactmicp_momas ${REACTMICP_LIBS})
 
 # Carbonation
 
 add_custom_target(ex_reactmicp_carbo_conf
     SOURCES reactmicp/saturated_react/carbonation.yaml)
 
 add_executable(ex_reactmicp_carbo
     ${DONT_COMPILE_SATURATED_EXAMPLE}
     reactmicp/saturated_react/carbonation.cpp
 )
 target_link_libraries(ex_reactmicp_carbo ${REACTMICP_LIBS})
 
 add_executable(ex_reactmicp_carbo_static
     EXCLUDE_FROM_ALL
     reactmicp/saturated_react/carbonation.cpp
 )
 target_link_libraries(ex_reactmicp_carbo_static ${REACTMICP_STATIC_LIBS})
 
 
 # Carbonation with Fe
 
 add_executable(ex_reactmicp_carbofe
     ${DONT_COMPILE_SATURATED_EXAMPLE}
     reactmicp/saturated_react/carbonationfe.cpp
 )
 target_link_libraries(ex_reactmicp_carbofe ${REACTMICP_LIBS})
 
 # Leaching with hydration
 
 add_executable(ex_reactmicp_leaching_kinetic
     ${DONT_COMPILE_SATURATED_EXAMPLE}
     reactmicp/saturated_react/leaching_kinetic.cpp
 )
 target_link_libraries(ex_reactmicp_leaching_kinetic ${REACTMICP_STATIC_LIBS})
 
 
 
 # Unsaturated system
 # ==================
 
 set_exclude_if_no_reactmicp_system(
     DONT_COMPILE_UNSATURATED_EXAMPLE
     REACTMICP_ENABLE_SYSTEM_UNSATURATED
 )
 
 # drying
 # ------
 
 add_executable(ex_reactmicp_drying
     ${DONT_COMPILE_UNSATURATED_EXAMPLE}
     reactmicp/unsaturated/drying.cpp
 )
 
 set(CEMDATA_PATH \"../data/cemdata.yaml\")
 set_source_files_properties(reactmicp/unsaturated/drying.cpp PROPERTIES COMPILE_DEFINITIONS
 "CEMDATA_PATH=${CEMDATA_PATH}"
 )
 
 target_link_libraries(ex_reactmicp_drying ${REACTMICP_STATIC_LIBS})
 
 # carbonation
 # -----------
 
 if (HDF5_FOUND)
 add_executable(ex_reactmicp_acc_carbo
     ${DONT_COMPILE_UNSATURATED_EXAMPLE}
     reactmicp/unsaturated/acc_carbo.cpp
 )
 
 set(CEMDATA_PATH \"../data/cemdata.yaml\")
 set_source_files_properties(reactmicp/unsaturated/acc_carbo.cpp PROPERTIES COMPILE_DEFINITIONS
 "CEMDATA_PATH=${CEMDATA_PATH}"
 )
 include_directories(${HDF5_INCLUDE_DIRS})
 set_target_properties(ex_reactmicp_acc_carbo
 	PROPERTIES COMPILE_DEFINITIONS HDF5_DEFINITIONS)
 target_link_libraries(ex_reactmicp_acc_carbo ${REACTMICP_STATIC_LIBS})
 endif()
 
 # Chloride
 # ========
 
 # migration
 # ---------
 
 set_exclude_if_no_reactmicp_system(
     DONT_COMPILE_CHLORIDE_EXAMPLE
     REACTMICP_ENABLE_SYSTEM_CHLORIDE
 )
 
 add_executable(ex_reactmicp_migration
     ${DONT_COMPILE_CHLORIDE_EXAMPLE}
     reactmicp/chloride/migration.cpp
 )
 
 set(CEMDATA_PATH \"../data/cemdata.yaml\")
 set_source_files_properties(reactmicp/chloride/migration.cpp PROPERTIES COMPILE_DEFINITIONS
 "CEMDATA_PATH=${CEMDATA_PATH}"
 )
 
 target_link_libraries(ex_reactmicp_migration ${REACTMICP_STATIC_LIBS})
 
 
 add_executable(ex_reactmicp_benchmark_rasouli
     ${DONT_COMPILE_CHLORIDE_EXAMPLE}
     reactmicp/chloride/benchmark_rasouli.cpp
 )
 
 set(CEMDATA_PATH \"../data/cemdata.yaml\")
 set_source_files_properties(reactmicp/chloride/benchmark_rasouli.cpp PROPERTIES COMPILE_DEFINITIONS
 "CEMDATA_PATH=${CEMDATA_PATH}"
 )
 
 target_link_libraries(ex_reactmicp_benchmark_rasouli ${REACTMICP_STATIC_LIBS})
+
+# Single
+# ======
+
+set_exclude_if_no_reactmicp_system(
+    DONT_COMPILE_SINGLE_EXAMPLE
+    REACTMICP_ENABLE_SYSTEM_SINGLE
+    )
+
+add_executable(ex_reactmicp_single
+    ${DONT_COMPILE_SINGLE_EXAMPLE}
+    reactmicp/single/single.cpp
+)
+target_link_libraries(ex_reactmicp_single ${REACTMICP_STATIC_LIBS})
diff --git a/examples/reactmicp/single/conf_single.yaml b/examples/reactmicp/single/conf_single.yaml
new file mode 100644
index 0000000..664d62f
--- /dev/null
+++ b/examples/reactmicp/single/conf_single.yaml
@@ -0,0 +1,18 @@
+parameters:
+     rate_factor: 1
+     freund_coeff: 0.2
+     freund_exp: 2
+     diffusion_coeff: 1e-5 
+     porosity: 0.2
+     saturation: 1
+     
+mesh:
+    type: uniform_mesh
+    uniform_mesh:
+        dx: 0.05
+        nb_nodes: 200
+        section: 5.0
+        
+runner:
+    save_time: 7200
+    run_time: 5*72000
diff --git a/examples/reactmicp/single/single.cpp b/examples/reactmicp/single/single.cpp
new file mode 100644
index 0000000..82fe86c
--- /dev/null
+++ b/examples/reactmicp/single/single.cpp
@@ -0,0 +1,197 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#include "reactmicp/reactmicp_single.hpp"
+#include "reactmicp/solver/runner.hpp"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+
+#include "specmicp_common/io/safe_config.hpp"
+
+using namespace specmicp;
+using namespace specmicp::io;
+using namespace specmicp::reactmicp::systems::single;
+
+struct SingleSimpleConf {
+    double rate_factor;
+    double diffusion_coeff;
+    double freund_coeff;
+    double freund_exp;
+    double porosity;
+    double saturation;
+};
+
+std::shared_ptr<SingleSimpleConf>
+configure_single_simple(
+        YAMLConfigHandle&& configuration
+        )
+{
+    auto conf = std::make_shared<SingleSimpleConf>();
+    conf->rate_factor = configuration.get_required_attribute<scalar_t>("rate_factor");
+    conf->freund_coeff = configuration.get_required_attribute<scalar_t>("freund_coeff");
+    conf->freund_exp = configuration.get_required_attribute<scalar_t>("freund_exp");
+    conf->diffusion_coeff = configuration.get_required_attribute<scalar_t>("diffusion_coeff");
+    conf->porosity = configuration.get_required_attribute<scalar_t>("porosity");
+    conf->saturation = configuration.get_required_attribute<scalar_t>("saturation");
+    return conf;
+}
+
+
+class RateFunctor
+{
+public:
+    RateFunctor(scalar_t coeff, scalar_t exp, scalar_t rate):
+        m_coeff(coeff), m_exp(exp), m_rate(rate)
+    {}
+
+    scalar_t operator()(index_t node, scalar_t dt, SingleVariables * const var) {
+        auto c_eq = std::pow(var->bound_concentration(node)/m_coeff, m_exp);
+        auto rate =  m_rate * (var->fluid_concentration(node) - c_eq);
+        //if (node == 1) {
+        //    std::cout << "rate : " << rate << " - " << c_eq << std::endl;
+        //}
+        if (rate < 0) {
+            rate = 0;
+        }
+        auto factor = var->porosity(node)*var->saturation(node);
+        auto new_c = var->fluid_concentration(node)-rate/factor*dt;
+        if (new_c < 0) {
+            rate = factor*var->fluid_concentration(node)/dt;
+        }
+        return rate;
+    }
+private:
+    scalar_t m_coeff;
+    scalar_t m_exp;
+    scalar_t m_rate;
+};
+
+void output(scalar_t time, specmicp::reactmicp::solver::VariablesBasePtr variables)
+{
+    auto vars = static_cast<specmicp::reactmicp::systems::single::SingleVariables*>(variables.get());
+
+    auto mesh = vars->get_mesh();
+
+    //Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> toplot(mesh->nb_nodes());
+    std::ofstream ofile;
+    ofile.open("output/"+std::to_string(time)+".csv");
+    for (int i=0;i<mesh->nb_nodes();++i) {
+        ofile << mesh->get_position(i) << ", " << vars->fluid_concentration(i) << ", " << vars->bound_concentration(i) << "\n";
+    }
+    ofile.close();
+}
+
+
+int main(int argc, char *argv[])
+{
+
+    if (argc == 1) {
+        std::cerr << "Error: expected one argument" << std::endl;
+        return 1;
+    }
+
+    std::string config_filepath(argv[1]);
+    auto config = YAMLConfigFile::load(config_filepath);
+
+//    specmicp::mesh::Uniform1DMeshGeometry geom;
+//    geom.nb_nodes = 200;
+//    geom.dx = 0.05;
+//    geom.section = 5.0;
+    specmicp::mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh(config.get_section("mesh"));
+    // specmicp::mesh::uniform_mesh1d(geom);
+    auto params = configure_single_simple(config.get_section("parameters"));
+
+    std::vector<index_t> list_fixed_nodes = {0,};
+
+    Vector init_cond = Vector::Zero(the_mesh->nb_nodes());
+    init_cond(0) = 1.0;
+    Vector init_cond_solid = Vector::Zero(the_mesh->nb_nodes());
+
+
+    auto variables = init_variables(the_mesh, list_fixed_nodes, init_cond, init_cond_solid);
+    for (int i=1;i<the_mesh->nb_nodes();++i) {
+        variables->diffusion_coefficient(i) = params->diffusion_coeff;
+        variables->porosity(i) = params->porosity;
+        variables->saturation(i) = params->saturation;
+    }
+    variables->diffusion_coefficient(0) = 1;
+    variables->porosity(0) = 1;
+    variables->saturation(0) = 1;
+
+
+    // transport stagger
+    auto transport_stagger =
+            std::make_shared<SingleTransportStagger>(variables, list_fixed_nodes);
+    auto& options = transport_stagger->options_solver();
+    options.alpha = 1.0;
+    //options.residuals_tolerance = 1e-10;
+    //options.quasi_newton = 3;
+
+    // chemistry stagger
+    auto functor = RateFunctor(params->freund_coeff, params->freund_exp, params->rate_factor);
+    std::shared_ptr<ChemistryProgram> chem_program = std::make_shared<DefaultChemistryProgram>(functor);
+    auto chemistry_stagger = std::make_shared<ChemistryStagger>(chem_program);
+
+    // upscaling stagger
+    auto upscaling_stagger = std::make_shared<SingleUpscalingStagger>();
+
+    // reactmicp solver
+    auto solver = specmicp::reactmicp::solver::ReactiveTransportSolver(transport_stagger, chemistry_stagger, upscaling_stagger);
+    auto& opts = solver.get_options();
+    opts.maximum_iterations = 51;
+    opts.residuals_tolerance = 1e-5;
+    //opts.
+
+    // runner
+
+    auto conf_runner = config.get_section("runner");
+
+    auto info = specmicp::reactmicp::solver::SimulationInformation("ex_reactmicp_single", conf_runner.get_required_attribute<scalar_t>("save_time"));
+    auto runner = std::make_shared<specmicp::reactmicp::solver::ReactiveTransportRunner>(solver, 0.01, 200, info);
+    runner->set_output_policy(output);
+    runner->run_until(conf_runner.get_required_attribute<scalar_t>("run_time"), variables);
+
+    for (int i=0; i<the_mesh->nb_nodes(); ++i) {
+        std::cout << variables->fluid_concentration(i) << "\n";
+    }
+    std::cout << "----------\n" << std::endl;
+    for (int i=0; i<the_mesh->nb_nodes(); ++i) {
+        std::cout << variables->bound_concentration(i) << "\n";
+    }
+    std::cout << std::endl;
+    return 0;
+}
diff --git a/src/reactmicp/CMakeLists.txt b/src/reactmicp/CMakeLists.txt
index d639bb3..f55ea61 100644
--- a/src/reactmicp/CMakeLists.txt
+++ b/src/reactmicp/CMakeLists.txt
@@ -1,120 +1,121 @@
 #   ReactMiCP
 # =============
 
 
 # var to store the files
 
 set( reactmicp_srcs ""
     CACHE INTERNAL "reactmicp files" FORCE)
 
 set( reactmicp_headers ""
     CACHE INTERNAL "reactmicp headers" FORCE
 )
 
 # macro to add to vars
 
 macro(add_to_reactmicp_srcs_list LIST_NAME)
     set(tmp "")
     foreach (src ${${LIST_NAME}})
         list(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${src})
     endforeach(src)
     set(reactmicp_srcs "${reactmicp_srcs};${tmp}" CACHE INTERNAL
             "reactmicp files" FORCE)
 endmacro(add_to_reactmicp_srcs_list)
 
 
 macro(add_to_reactmicp_headers_list LIST_NAME)
     set( tmp "")
     foreach(header ${${LIST_NAME}})
        LIST(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${header})
     endforeach(header)
     set(reactmicp_headers "${reactmicp_headers};${tmp}"
             CACHE INTERNAL "reactmicp headers" FORCE)
 endmacro(add_to_reactmicp_headers_list)
 
 # gather files
 
 set(reactmicp_main_headers
     reactmicp_core.hpp
     reactmicp.hpp
     reactmicp_saturated.hpp
     reactmicp_unsaturated.hpp
     reactmicp_chloride.hpp
+    reactmicp_single.hpp
 )
 add_to_reactmicp_headers_list(reactmicp_main_headers)
 
 INSTALL(FILES ${reactmicp_main_headers}
     DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp
 )
 
 
 add_subdirectory( equilibrium_curve )
 add_subdirectory( io                )
 add_subdirectory( solver            )
 add_subdirectory( systems           )
 
 # Library
 # ========
 
 spc_set_pgo_flag(${reactmicp_srcs})
 
 add_library(objreactmicp OBJECT
     ${reactmicp_srcs}
     ${reactmicp_headers}
     )
 target_include_directories(objreactmicp PUBLIC
     ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS} ${YAML_INCLUDE_DIRS})
 if(EIGEN3_UNSUPPORTED_FOUND)
     set_property(TARGET objreactmicp APPEND PROPERTY
         COMPILE_DEFINITIONS EIGEN3_UNSUPPORTED_FOUND)
     target_include_directories(objreactmicp PRIVATE
         ${EIGEN3_UNSUPPORTED_INCLUDE_DIR}) # added after
 endif(EIGEN3_UNSUPPORTED_FOUND)
 set_target_properties(objreactmicp PROPERTIES
     POSITION_INDEPENDENT_CODE 1
     CXX_STANDARD ${SPECMICP_CXX_STANDARD}
     CXX_STANDARD_REQUIRED ON
 )
 spc_set_objlib_optimization_flags(objreactmicp)
 spc_set_visibility_flag(objreactmicp)
 
 
 add_library(reactmicp SHARED $<TARGET_OBJECTS:objreactmicp>)
 
 target_link_libraries(reactmicp PUBLIC dfpm specmicp_common specmicp specmicp_database)
 target_link_libraries(reactmicp PUBLIC PkgConfig::YAML)
 spc_target_link_profile(reactmicp)
 spc_set_cxx_version_interface(reactmicp)
 spc_set_target_optimization_flags(reactmicp)
 install(TARGETS reactmicp
     EXPORT ${SPECMICP_TARGET}
     LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR}
 )
 
 
 # static libraries
 # ----------------
 
 if(SPECMICP_BUILD_STATIC)
     add_library(reactmicp_static STATIC $<TARGET_OBJECTS:objreactmicp>)
     install(TARGETS reactmicp_static
         EXPORT ${SPECMICP_TARGET}
         ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR}
     )
 else()
     add_library(reactmicp_static EXCLUDE_FROM_ALL STATIC $<TARGET_OBJECTS:objreactmicp>)
 endif()
 
 set_target_properties(reactmicp_static  PROPERTIES
     OUTPUT_NAME reactmicp
 )
 spc_set_cxx_version_interface(reactmicp_static)
 spc_set_target_optimization_flags(reactmicp_static)
 target_link_libraries(reactmicp_static PUBLIC
     dfpm_static
     specmicp_static
     specmicp_database_static
     specmicp_common_static
     )
 target_link_libraries(reactmicp_static PUBLIC PkgConfig::YAML)
 spc_target_link_profile(reactmicp_static)
diff --git a/src/reactmicp/reactmicp_single.hpp b/src/reactmicp/reactmicp_single.hpp
new file mode 100644
index 0000000..ca5dd27
--- /dev/null
+++ b/src/reactmicp/reactmicp_single.hpp
@@ -0,0 +1,56 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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 reactmicp_single.hpp
+//! \brief Include this header to use the single system of ReactMiCP
+
+#ifndef SPECMICP_REACTMICP_SINGLE_HPP
+#define SPECMICP_REACTMICP_SINGLE_HPP
+
+// Core files
+
+#include "reactmicp_core.hpp"
+
+// systems
+
+#include "reactmicp/systems/single/variables.hpp"
+#include "reactmicp/systems/single/chemistry_program.hpp"
+#include "reactmicp/systems/single/chemistry_stagger.hpp"
+#include "reactmicp/systems/single/transport_stagger.hpp"
+#include "reactmicp/systems/single/init_variables.hpp"
+#include "reactmicp/systems/single/default_upscaling_stagger.hpp"
+
+#endif // SPECMICP_REACTMICP_SINGLE_HPP
diff --git a/src/reactmicp/systems/CMakeLists.txt b/src/reactmicp/systems/CMakeLists.txt
index 84f957b..2833f9e 100644
--- a/src/reactmicp/systems/CMakeLists.txt
+++ b/src/reactmicp/systems/CMakeLists.txt
@@ -1,15 +1,19 @@
 # The different systems
 
 add_subdirectory( common )
 
 if ( REACTMICP_ENABLE_SYSTEM_SATURATED )
 add_subdirectory( saturated_react )
 endif ( REACTMICP_ENABLE_SYSTEM_SATURATED )
 
 if ( REACTMICP_ENABLE_SYSTEM_UNSATURATED )
     add_subdirectory( unsaturated )
 endif ( REACTMICP_ENABLE_SYSTEM_UNSATURATED )
 
 if ( REACTMICP_ENABLE_SYSTEM_CHLORIDE )
     add_subdirectory( chloride )
 endif ( REACTMICP_ENABLE_SYSTEM_CHLORIDE )
+
+if ( REACTMICP_ENABLE_SYSTEM_SINGLE )
+    add_subdirectory( single )
+endif ( REACTMICP_ENABLE_SYSTEM_SINGLE )
diff --git a/src/reactmicp/systems/single/CMakeLists.txt b/src/reactmicp/systems/single/CMakeLists.txt
new file mode 100644
index 0000000..c1b0481
--- /dev/null
+++ b/src/reactmicp/systems/single/CMakeLists.txt
@@ -0,0 +1,30 @@
+set( reactmicp_single_srcs
+
+    default_upscaling_stagger.cpp
+    chemistry_stagger.cpp
+    chemistry_program.cpp
+    init_variables.cpp
+    transport_program.cpp
+    transport_stagger.cpp
+    variables.cpp
+)
+
+set( reactmicp_single_headers
+    variablesfwd.hpp
+
+    default_upscaling_stagger.hpp
+    chemistry_program.hpp
+    chemistry_stagger.hpp
+    init_variables.hpp
+    transport_program.hpp
+    transport_stagger.hpp
+    variables.hpp
+)
+
+add_to_reactmicp_srcs_list(reactmicp_single_srcs)
+add_to_reactmicp_headers_list(reactmicp_single_headers)
+
+INSTALL(FILES ${reactmicp_saturated_headers}
+    DESTINATION ${INCLUDE_INSTALL_DIR}/reactmicp/systems/single
+)
+
diff --git a/src/reactmicp/systems/single/chemistry_program.cpp b/src/reactmicp/systems/single/chemistry_program.cpp
new file mode 100644
index 0000000..9f12d0a
--- /dev/null
+++ b/src/reactmicp/systems/single/chemistry_program.cpp
@@ -0,0 +1,68 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+
+#include "chemistry_program.hpp"
+#include "variables.hpp"
+#include "dfpm/mesh.hpp"
+//#include <iostream>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+ChemistryProgramReturn DefaultChemistryProgram::solve_one_node(index_t node, scalar_t dt, SingleVariables *const var)
+{
+    auto solid_rate = m_rate_function(node, dt, var);
+    auto chem_return = ChemistryProgramReturn();
+
+    chem_return.bound_concentration = var->bound_concentration(node, var->predictor())+dt*solid_rate;
+
+    auto factor = var->porosity(node)*var->saturation(node);
+    // FIXME alpha
+    chem_return.fluid_concentration = var->fluid_concentration(node, var->predictor())+dt/factor*(
+                -solid_rate+var->fluid_concentration(node, var->transport_rate())/var->get_mesh()->get_volume_cell(node));
+
+
+    chem_return.retcode = ReturnCode::SuccessWithoutWarning;
+
+    return chem_return;
+}
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
diff --git a/src/reactmicp/systems/single/chemistry_program.hpp b/src/reactmicp/systems/single/chemistry_program.hpp
new file mode 100644
index 0000000..9c850b4
--- /dev/null
+++ b/src/reactmicp/systems/single/chemistry_program.hpp
@@ -0,0 +1,88 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYPROGRAM_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYPROGRAM_HPP
+
+#include "specmicp_common/types.hpp"
+#include "specmicp_common/return_code.hpp"
+#include "variablesfwd.hpp"
+
+#include <functional>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+struct ChemistryProgramReturn
+{
+    ReturnCode retcode;
+    // only assume that we get the concentration (i.e. case of equilibrium stagger)
+    scalar_t fluid_concentration;
+    scalar_t bound_concentration;
+};
+
+class ChemistryProgram
+{
+public:
+    ChemistryProgram() {}
+
+    virtual ChemistryProgramReturn solve_one_node(index_t node, scalar_t dt, SingleVariables * const var) =0;
+};
+
+using rate_f = std::function<scalar_t (index_t, scalar_t, SingleVariables * const)>;
+
+class DefaultChemistryProgram:
+        public ChemistryProgram
+{
+public:
+    DefaultChemistryProgram(rate_f rate_function):
+        m_rate_function(rate_function)
+    {}
+
+    virtual ChemistryProgramReturn solve_one_node(index_t node, scalar_t dt, SingleVariables * const var) override;
+
+private:
+    rate_f m_rate_function;
+};
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYPROGRAM_HPP
diff --git a/src/reactmicp/systems/single/chemistry_stagger.cpp b/src/reactmicp/systems/single/chemistry_stagger.cpp
new file mode 100644
index 0000000..6f981bd
--- /dev/null
+++ b/src/reactmicp/systems/single/chemistry_stagger.cpp
@@ -0,0 +1,162 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+#include "chemistry_stagger.hpp"
+#include "variables.hpp"
+#include "chemistry_program.hpp"
+
+#include "specmicp_common/log.hpp"
+#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
+
+#include "specmicp_common/config.h"
+#include "specmicp_common/return_code.hpp"
+
+#include <iostream>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+using TrueConstPtr = SingleVariables * const; //!< const Ptr to the variables
+
+//! \brief cast base variables to saturated variables
+inline TrueConstPtr cast_to_var(solver::VariablesBase * const var)
+{
+    return static_cast<TrueConstPtr>(var);
+}
+
+//! \brief Initialize the stagger at the beginning of an iteration
+void ChemistryStagger::initialize_timestep(
+        scalar_t dt,
+        VariablesBase * const var)
+{
+    m_dt = dt;
+
+    TrueConstPtr true_var = cast_to_var(var);
+
+    // Initialize velocity using values from previous timestep
+    for (index_t node=0; node<true_var->nb_nodes(); ++node)
+    {
+        if (true_var->is_fixed_composition(node)) continue;
+        scalar_t alpha = 1.0;
+        alpha = std::max(alpha,
+                       0.9*dt*true_var->bound_concentration(node, true_var->chemistry_rate())
+                       /(true_var->bound_concentration(node, true_var->displacement()))
+                         );
+
+        // auto solid_velocity = true_var->bound_concentration(node, true_var->velocity());
+        auto solid_chemistry_rate = true_var->bound_concentration(node, true_var->chemistry_rate());
+        true_var->bound_concentration(node, true_var->velocity()) = 1/alpha * solid_chemistry_rate;
+    }
+    //true_var->transport_rate().setZero();
+    //true_var->chemistry_rate().setZero();
+}
+
+//! \brief Solve the equation for the timestep
+solver::StaggerReturnCode
+ChemistryStagger::restart_timestep(VariablesBase * const var)
+{
+    TrueConstPtr true_var = cast_to_var(var);
+    int failed_chemistry = 0;
+#ifdef SPECMICP_HAVE_OPENMP
+#pragma omp parallel  default(none) shared(failed_chemistry, true_var)
+    {
+#pragma omp for schedule(dynamic, 5)
+        for (index_t node=0; node<true_var->nb_nodes(); ++node)
+        {
+            // only solve if necessary
+            if (true_var->is_fixed_composition(node) or failed_chemistry > 0) continue;
+            const auto retcode = solve_one_node(node, true_var);
+            if (retcode > 0)
+            {
+                ++failed_chemistry;
+            }
+        }
+    }
+#else
+    {
+        for (index_t node=0; node<true_var->nb_nodes(); ++node)
+        {
+            if (true_var->is_fixed_composition(node)) continue;
+            const auto retcode = solve_one_node(node, true_var);
+            if (retcode > 0)
+            {
+                ++failed_chemistry;
+                break;
+            }
+        }
+    }
+#endif // SPECMICP_HAVE_OPENMP
+    if (failed_chemistry > 0)
+        return solver::StaggerReturnCode::UnknownError;
+    return solver::StaggerReturnCode::ResidualMinimized;
+}
+
+//!
+int ChemistryStagger::solve_one_node(
+        index_t node,
+        SingleVariables * const true_var
+        )
+{
+
+    auto chem_return = m_program->solve_one_node(node, m_dt, true_var);
+    if (chem_return.retcode <= ReturnCode::NotFinishedYet)
+    {
+        ERROR << "Failed to solve chemistry problem at node " << node
+              << ", return code = " << static_cast<int>(chem_return.retcode);
+        return  1;
+    }
+
+    auto b0 = true_var->bound_concentration(node, true_var->predictor());
+    auto b_vel = (chem_return.bound_concentration-b0)/m_dt;
+
+    true_var->bound_concentration(node, true_var->displacement()) = chem_return.bound_concentration;
+    true_var->bound_concentration(node, true_var->velocity()) = b_vel;
+    true_var->bound_concentration(node, true_var->chemistry_rate()) = b_vel;
+
+    auto c0 = true_var->fluid_concentration(node, true_var->predictor());
+    auto c_vel = (chem_return.fluid_concentration-c0)/m_dt;
+
+    true_var->fluid_concentration(node, true_var->displacement()) = chem_return.fluid_concentration;
+    true_var->fluid_concentration(node, true_var->velocity()) = c_vel;
+
+    return 0;
+}
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
diff --git a/src/reactmicp/systems/single/chemistry_stagger.hpp b/src/reactmicp/systems/single/chemistry_stagger.hpp
new file mode 100644
index 0000000..515c68a
--- /dev/null
+++ b/src/reactmicp/systems/single/chemistry_stagger.hpp
@@ -0,0 +1,107 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYSTAGGER_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYSTAGGER_HPP
+
+#include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp"
+#include "variablesfwd.hpp"
+
+#include "specmicp/adimensional/adimensional_system_solver_structs.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+class ChemistryProgram;
+using ChemistryProgramPtr = std::shared_ptr<ChemistryProgram>;
+
+using VariablesBase = solver::VariablesBase; //!< Type of the variables
+
+
+//! \brief Solve the equilibrium problem
+class SPECMICP_DLL_PUBLIC ChemistryStagger:
+        public solver::ChemistryStaggerBase
+{
+public:
+    //! \brief Constructor with uniform chemistry
+    ChemistryStagger(ChemistryProgramPtr& chem_program):
+        m_program(chem_program)
+    {}
+
+    //! \brief Initialize the stagger at the beginning of the computation
+    virtual void initialize(VariablesBase * const var) override {}
+
+    //! \brief Initialize the stagger at the beginning of the computation
+    void initialize(std::shared_ptr<SingleVariables>& var) {}
+
+    //! \brief Initialize the stagger at the beginning of an iteration
+    virtual void initialize_timestep(
+            scalar_t dt,
+            VariablesBase * const var
+            ) override;
+
+    //! \brief Solve the equation for the timestep
+    virtual solver::StaggerReturnCode restart_timestep(
+            VariablesBase * const var
+            ) override;
+
+    //! \brief Solve one node
+    int solve_one_node(
+            index_t node,
+            SingleVariables * const var
+            );
+
+
+    //! \brief Return the solution for one node
+    virtual ReturnCode solve_equilibrium_at_node(
+            index_t node, VariablesBase * const var, AdimensionalSystemSolution& out) override {
+        return ReturnCode::NotFinishedYet;
+    }; // N/A
+
+private:
+    scalar_t m_dt; //!< The timestep
+
+    ChemistryProgramPtr m_program; // The chemistry program
+};
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_CHEMISTRYSTAGGER_HPP
diff --git a/src/reactmicp/systems/single/default_upscaling_stagger.cpp b/src/reactmicp/systems/single/default_upscaling_stagger.cpp
new file mode 100644
index 0000000..e69de29
diff --git a/src/reactmicp/systems/single/default_upscaling_stagger.hpp b/src/reactmicp/systems/single/default_upscaling_stagger.hpp
new file mode 100644
index 0000000..439e04c
--- /dev/null
+++ b/src/reactmicp/systems/single/default_upscaling_stagger.hpp
@@ -0,0 +1,92 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_DEFAULTUPSCALINGPROGRAM_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_DEFAULTUPSCALINGPROGRAM_HPP
+
+#include "reactmicp/solver/staggers_base/variables_base.hpp"
+#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
+#include "variablesfwd.hpp"
+#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+
+#ifndef SPC_DOXYGEN_SHOULD_SKIP_THIS
+using StaggerReturnCode = solver::StaggerReturnCode;
+using VariablesBase = solver::VariablesBase;
+#endif // SPC_DOXYGEN_SHOULD_SKIP_THIS
+
+namespace systems {
+namespace single {
+
+class SPECMICP_DLL_PUBLIC SingleUpscalingStagger:
+        public solver::UpscalingStaggerBase
+{
+public:
+    //! \brief Constructor
+    //!
+    //! By defautl does nothing
+    SingleUpscalingStagger():
+        m_dt(HUGE_VAL)
+    {}
+
+    //! \brief Initialize the stagger at the beginning of the computation
+    //!
+    //! The porosity and transport properties should be initialized here
+    void initialize(VariablesBase * const var) override {}
+
+    //! \brief Initialize the stagger at the beginning of a timestep
+    //!
+    //! Typically do nothing but erase the "dot" variable (velocity such as the vel_porosity)
+    void initialize_timestep(scalar_t dt, VariablesBase * const var) override {
+        m_dt = dt;
+    }
+
+    //! \brief This is the main function called during a timestep
+    StaggerReturnCode restart_timestep(VariablesBase * const var) override {
+        return StaggerReturnCode::ResidualMinimized;
+    }
+private:
+    scalar_t m_dt; //!< The current timestep
+};
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_DEFAULTUPSCALINGPROGRAM_HPP
diff --git a/src/reactmicp/systems/single/init_variables.cpp b/src/reactmicp/systems/single/init_variables.cpp
new file mode 100644
index 0000000..6d953ba
--- /dev/null
+++ b/src/reactmicp/systems/single/init_variables.cpp
@@ -0,0 +1,114 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+
+#include "init_variables.hpp"
+#include "variables.hpp"
+#include "dfpm/meshes/mesh1d.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+SingleVariablesPtr SPECMICP_DLL_PUBLIC init_variables(
+        mesh::Mesh1DPtr the_mesh,
+        const std::vector<index_t>& list_fixed_nodes,
+        const Vector& liquid_concentration,
+        const Vector& solid_concentration
+        )
+{
+    auto factory = SingleVariablesFactory(the_mesh, list_fixed_nodes, liquid_concentration, solid_concentration);
+    return factory.get_variable();
+}
+
+SingleVariablesFactory::SingleVariablesFactory(
+        mesh::Mesh1DPtr the_mesh,
+        const std::vector<index_t>& list_fixed_nodes,
+        const Vector& liquid_concentration,
+        const Vector& solid_concentration
+        ):
+    m_variable(std::make_shared<SingleVariables>(the_mesh)),
+    nb_nodes(the_mesh->nb_nodes())
+{
+    assert(liquid_concentration.size() == nb_nodes);
+    assert(solid_concentration.size() == nb_nodes);
+
+    init_size();
+    set_fixed_nodes(list_fixed_nodes);
+
+    set_concentrations(liquid_concentration, solid_concentration);
+}
+
+void SingleVariablesFactory::init_size()
+{
+    index_t main_ndf = m_variable->ndf()*nb_nodes;
+
+    m_variable->displacement() = Vector::Zero(main_ndf);
+    m_variable->chemistry_rate() = Vector::Zero(main_ndf);
+    m_variable->transport_rate() = Vector::Zero(main_ndf);
+    m_variable->predictor() = Vector::Zero(main_ndf);
+    m_variable->velocity() = Vector::Zero(main_ndf);
+
+    m_variable->m_upscaling = Vector::Zero(m_variable->ndf_upscaling()*nb_nodes);
+}
+
+
+void SingleVariablesFactory::set_fixed_nodes(const std::vector<index_t>& list_fixed_nodes)
+{
+    m_variable->m_is_fixed_composition = std::vector<bool>(nb_nodes, false);
+    for (index_t node: list_fixed_nodes)
+    {
+        m_variable->m_is_fixed_composition[node] = true;
+    }
+}
+
+
+void SingleVariablesFactory::set_concentrations(
+        const Vector& liquid_concentration,
+        const Vector& solid_concentration) {
+
+    for (index_t node=0; node<nb_nodes; ++node)
+    {
+        m_variable->fluid_concentration(node) = liquid_concentration(node);
+        m_variable->bound_concentration(node) = solid_concentration(node);
+    }
+}
+
+
+} // namespace single
+} // namespace systems
+} // namespace reactmicp
+} // namespace specmicp
diff --git a/src/reactmicp/systems/single/init_variables.hpp b/src/reactmicp/systems/single/init_variables.hpp
new file mode 100644
index 0000000..c298d22
--- /dev/null
+++ b/src/reactmicp/systems/single/init_variables.hpp
@@ -0,0 +1,100 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 2021-     F. Georget <georget@ibac.rwth-aachen.de> RWTH Aachen
+ 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. *
+
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_INITVARIABLES_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_INITVARIABLES_HPP
+
+#include "variablesfwd.hpp"
+#include "dfpm/meshes/mesh1dfwd.hpp"
+#include "specmicp_database/database_fwd.hpp"
+#include "specmicp_common/physics/units.hpp"
+
+#include <vector>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+//! \brief Initialize an instance of Saturated Variables
+class SingleVariablesFactory
+{
+public:
+    //! \brief Constructor
+    //!
+    //! \param the_mesh the mesh
+    //! \param the_database SpecMiCP database
+    //! \param the_units units of the variables and used in computation
+    //! \param list_fixed_nodes list of nodes where variables are fixed
+    //! \param list_initial_states a list of the different initial conditions
+    //! \param index_initial_state a map combining nodes (index) to the initial conditions (value)
+    SingleVariablesFactory(
+            mesh::Mesh1DPtr the_mesh,
+            const std::vector<index_t>& list_fixed_nodes,
+            const Vector& liquid_concentration,
+            const Vector& solid_concentration
+            );
+    //! \brief Initialize the main vectors
+    void init_size();
+    //! \brief Initialize the BC
+    void set_fixed_nodes(const std::vector<index_t>& list_fixed_nodes);
+    //! \brief Initialize the initial conditions
+    void set_concentrations(
+            const Vector& liquid_concentration,
+            const Vector& solid_concentration);
+    //! \brief Return the variables
+    SingleVariablesPtr get_variable() {return m_variable;}
+
+private:
+    SingleVariablesPtr m_variable;
+    index_t nb_nodes;
+};
+
+//! \brief Initialise an instance of SaturatedVariables
+SingleVariablesPtr SPECMICP_DLL_PUBLIC init_variables(
+        mesh::Mesh1DPtr the_mesh,
+        const std::vector<index_t>& list_fixed_nodes,
+        const Vector& liquid_concentration,
+        const Vector& solid_concentration
+        );
+
+} // end namespace satdiff
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_INITVARIABLES_HPP
+
diff --git a/src/reactmicp/systems/single/transport_program.cpp b/src/reactmicp/systems/single/transport_program.cpp
new file mode 100644
index 0000000..44050e4
--- /dev/null
+++ b/src/reactmicp/systems/single/transport_program.cpp
@@ -0,0 +1,329 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#include "transport_program.hpp"
+#include "variables.hpp"
+#include "dfpm/meshes/mesh1d.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+//class SingleDiffusion::
+
+
+SingleDiffusion::SingleDiffusion(
+        SingleVariables* variables,
+        std::vector<index_t> list_fixed_nodes
+        ):
+    m_ndf(2),
+    m_tot_ndf(2*variables->get_mesh()->nb_nodes()),
+    m_mesh(variables->get_mesh()),
+    m_variables(variables),
+    m_is_in_residual_computation(false)
+{
+    number_equations(list_fixed_nodes, {});
+}
+
+
+SingleDiffusion::SingleDiffusion(
+        SingleVariables* variables,
+        std::vector<index_t> list_fixed_nodes,
+        std::map<index_t, index_t>  list_slave_nodes
+        ):
+    m_ndf(2),
+    m_tot_ndf(2*variables->get_mesh()->nb_nodes()),
+    m_mesh(variables->get_mesh()),
+    m_variables(variables),
+    m_is_in_residual_computation(false)
+{
+    number_equations(list_fixed_nodes, list_slave_nodes);
+}
+
+void SingleDiffusion::number_equations(
+        std::vector<index_t> list_fixed_nodes,
+        std::map<index_t, index_t>  list_slave_nodes
+        )
+{
+
+    m_ideq.resizeLike(m_variables->displacement());
+    m_ideq.setZero();
+    // flag fixed nodes
+    for (index_t node: list_fixed_nodes)
+    {
+        m_ideq(m_variables->dof_fluid_concentration(node)) = no_equation;
+    }
+    // flag slaves nodes
+    // we flag them by making their ideq more negative than no_equation
+    for (auto slave_pair: list_slave_nodes)
+    {
+        m_ideq(m_variables->dof_fluid_concentration(slave_pair.first)) = no_equation-1;
+    }
+    // set equation numbers
+    index_t neq = 0;
+    for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
+    {
+        const index_t dof = m_variables->dof_fluid_concentration(node);
+        if (m_ideq(dof) > no_equation) // attribute an equation number if it is NOT a slave nor a fixed node
+        {
+            m_ideq(dof) = neq;
+            ++neq;
+        }
+        m_ideq(m_variables->dof_bound_concentration(node)) = no_equation;
+    }
+    // slave nodes
+    // attribute the correct equation number
+    for (auto slave_pair: list_slave_nodes)
+    {
+        m_ideq(m_variables->dof_fluid_concentration(slave_pair.first)) =
+                    m_ideq(m_variables->dof_fluid_concentration(slave_pair.second));
+    }
+    m_neq = neq;
+}
+
+void SingleDiffusion::compute_residuals(
+        const Vector& displacement,
+        const Vector& velocity,
+        Vector& residual,
+        bool use_chemistry_rate
+        )
+{
+    residual = Vector::Zero(get_neq());
+    m_is_in_residual_computation = true;
+    m_variables->transport_rate().setZero(); // important !
+    for (index_t element: m_mesh->range_elements())
+    {
+        Eigen::Vector2d element_residual;
+        element_residual.setZero();
+        residuals_element_component(element, displacement, velocity, element_residual, use_chemistry_rate);
+        for (index_t en=0; en<2; ++en)
+        {
+                const index_t node = m_mesh->get_node(element, en);
+                const index_t id = m_ideq(m_variables->dof_fluid_concentration(node));
+                if (id != no_equation) {residual(id) += element_residual(en);}
+        }
+    }
+    m_is_in_residual_computation = false;
+}
+
+void SingleDiffusion::residuals_element_component(
+        index_t element,
+        const Vector& displacement,
+        const Vector& velocity,
+        Eigen::Vector2d& element_residual,
+        bool use_chemistry_rate
+        )
+{
+    const scalar_t mass_coeff_0 = m_mesh->get_volume_cell_element(element, 0);
+    const scalar_t mass_coeff_1 = m_mesh->get_volume_cell_element(element, 1);
+
+    const index_t node_0 = m_mesh->get_node(element, 0);
+    const index_t node_1 = m_mesh->get_node(element, 1);
+
+    const index_t dof_0 = m_variables->dof_fluid_concentration(node_0);
+    const index_t dof_1 = m_variables->dof_fluid_concentration(node_1);
+
+
+    const scalar_t velocity_0 = mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0)*m_variables->saturation(node_0)
+                                        +m_variables->vel_porosity(node_0)*displacement(dof_0)*m_variables->saturation(node_0)
+                                        +m_variables->vel_saturation(node_0)*displacement(dof_0)*m_variables->porosity(node_0)
+                                              );
+
+    const scalar_t velocity_1 =  mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1)*m_variables->saturation(node_1)
+                                         + m_variables->vel_porosity(node_1)*displacement(dof_1)*m_variables->saturation(node_1)
+                                         + m_variables->vel_saturation(node_1)*displacement(dof_1)*m_variables->porosity(node_1)
+                                               );
+    scalar_t transport_0 = 0.0;
+    scalar_t transport_1 = 0.0;
+
+
+    // transport
+    scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) +
+                               0.5/m_variables->diffusion_coefficient(node_1));
+
+    scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
+                             * diff_coeff
+                           );
+
+
+    // diffusion
+    scalar_t flux_diffusion = flux_coeff*(displacement(dof_0) - displacement(dof_1));
+
+    transport_0 = flux_diffusion;
+    transport_1 = - flux_diffusion;
+
+//    // advection
+
+//    if (m_variables->fluid_velocity(element) != 0.0)
+//    {
+//        scalar_t flux_advection = (m_mesh->get_face_area(element))
+//                *m_variables->fluid_velocity(element);
+//        if (m_variables->fluid_velocity(element) > 0)
+//        {
+//            flux_advection *= (displacement(dof_0) - displacement(dof_1));
+//            transport_1 += flux_advection;
+//        }
+//        else
+//        {
+//            flux_advection *= (displacement(dof_1) - displacement(dof_0));
+//            transport_0 -= flux_advection;
+//        }
+//    }
+
+    if (m_is_in_residual_computation)
+    {
+        m_variables->fluid_concentration(node_0, m_variables->transport_rate()) += transport_0;
+        m_variables->fluid_concentration(node_1, m_variables->transport_rate()) += transport_1;
+    }
+
+    if (use_chemistry_rate)
+    {
+        const scalar_t chemistry_0 = mass_coeff_0*(
+                    -  m_variables->bound_concentration(node_0, m_variables->chemistry_rate())
+                    +  m_variables->fluid_concentration(node_0, m_variables->chemistry_rate())
+                    );
+        const scalar_t chemistry_1 = mass_coeff_1*(
+                    -  m_variables->bound_concentration(node_1, m_variables->chemistry_rate())
+                    +  m_variables->fluid_concentration(node_1, m_variables->chemistry_rate())
+                    );
+        element_residual(0) = + transport_0 + chemistry_0 - velocity_0;
+        element_residual(1) = + transport_1 + chemistry_1 - velocity_1;
+    }
+    else
+    {
+        element_residual(0) = + transport_0 - velocity_0;
+        element_residual(1) = + transport_1 - velocity_1;
+    }
+}
+
+
+void SingleDiffusion::compute_jacobian(
+        Vector& displacement,
+        Vector& velocity,
+        Eigen::SparseMatrix<scalar_t>& jacobian,
+        scalar_t alphadt
+        )
+{
+    dfpm::list_triplet_t jacob;
+    const index_t ncomp = m_variables->nb_component();
+    const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
+    jacob.reserve(estimation);
+    for (index_t element: m_mesh->range_elements())
+    {
+        jacobian_element(element, displacement, velocity, jacob, alphadt);
+    }
+    jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
+    jacobian.setFromTriplets(jacob.begin(), jacob.end());
+}
+
+void SingleDiffusion::jacobian_element(
+        index_t element,
+        Vector& displacement,
+        Vector& velocity,
+        dfpm::list_triplet_t& jacobian,
+        scalar_t alphadt)
+{
+    Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero());
+    // note : jacobian computation do not need to use the chemistry rate
+    residuals_element_component(element, displacement, velocity, element_residual_orig, NO_CHEMISTRY_RATE);
+
+    for (index_t en=0; en<2; ++en)
+    {
+        Eigen::Vector2d element_residual;
+
+        const index_t node = m_mesh->get_node(element, en);
+        const index_t dof = m_variables->dof_fluid_concentration(node);
+        const index_t idc = m_ideq(dof);
+
+        if (idc == no_equation) continue;
+
+        const scalar_t tmp_v = velocity(dof);
+        const scalar_t tmp_d = displacement(dof);
+
+        scalar_t h = eps_jacobian*std::abs(tmp_v);
+        if (h < 1e-4*eps_jacobian) h = eps_jacobian;
+        velocity(dof) = tmp_v + h;
+        h = velocity(dof) - tmp_v;
+
+        displacement(dof) = tmp_d + alphadt*h;
+
+        element_residual.setZero();
+        residuals_element_component(element, displacement, velocity, element_residual, NO_CHEMISTRY_RATE);
+        velocity(dof) = tmp_v;
+        displacement(dof) = tmp_d;
+
+        for (index_t enr=0; enr<2; ++enr)
+        {
+            const index_t noder = m_mesh->get_node(element, enr);
+            const index_t idr = m_ideq(m_variables->dof_fluid_concentration(noder));
+
+            if (idr == no_equation) continue;
+            jacobian.push_back(dfpm::triplet_t(
+                                   idr,
+                                   idc,
+                                   (element_residual(enr) - element_residual_orig(enr))/h
+                                   ));
+        }
+    }
+}
+
+//! \brief Update the solutions
+void SingleDiffusion::update_solution(const Vector& update,
+                     scalar_t lambda,
+                     scalar_t alpha_dt,
+                     Vector& predictor,
+                     Vector& displacement,
+                     Vector& velocity)
+{
+    for (index_t node: m_mesh->range_nodes())
+    {
+        const index_t dof = m_variables->dof_fluid_concentration(node);
+        const index_t id = m_ideq(dof);
+        if (id == no_equation)
+            continue;
+        velocity(dof) += lambda*update(id);
+        //displacement(dof) = predictor(dof) + alpha_dt*velocity(dof);
+
+    }
+    //displacement = m_variables->predictor() + alpha_dt*velocity;
+
+    displacement = predictor + alpha_dt*velocity;
+}
+
+} // end namespace satdiff
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
diff --git a/src/reactmicp/systems/single/transport_program.hpp b/src/reactmicp/systems/single/transport_program.hpp
new file mode 100644
index 0000000..e7de3cf
--- /dev/null
+++ b/src/reactmicp/systems/single/transport_program.hpp
@@ -0,0 +1,187 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTPROGRAM_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTPROGRAM_HPP
+
+//! \name reactmicp/single/transport_program.hpp
+//! \brief The diffusion equations
+
+#include "variablesfwd.hpp"
+
+#include "dfpm/meshes/mesh1dfwd.hpp"
+#include "dfpm/types.hpp"
+#include "dfpm/solver/parabolic_program.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+//! \brief Diffusion equations for the single system
+//!
+//! Solve all equations together.
+class SingleDiffusion:
+        public dfpmsolver::ParabolicProgram<SingleDiffusion>
+{
+public:
+    static constexpr bool NO_CHEMISTRY_RATE {false}; //!< do not include chemistry rates
+    static constexpr bool CHEMISTRY_RATE {true};     //!< include chemistry rates
+
+    //! \brief Fancy constructor
+    SingleDiffusion(SingleVariables* variables,
+                       std::vector<index_t> list_fixed_nodes,
+                       std::map<index_t, index_t>  list_slave_nodes
+                    );
+
+    //! \brief Constructor
+    SingleDiffusion(SingleVariables* variables,
+                       std::vector<index_t> list_fixed_nodes);
+
+    //! \brief Return the number of equations
+    index_t get_neq() const {return m_neq;}
+    //! \brief Return the number of degrees of freedom per node
+    index_t get_ndf() const {return m_ndf;}
+    //! \brief Return the total number of degrees of freedom
+    index_t get_tot_ndf() const  {return m_tot_ndf;}
+
+    //! \brief Method to update the variables
+    void set_variables(SingleVariables* variables) {m_variables = variables;}
+
+    //! \brief Return the id of the equation corresponding to the degree of freedom 'id_dof'
+    //!
+    //! Return 'no_equation' if no equation exist
+    index_t id_equation(index_t id_dof) const {return m_ideq(id_dof);}
+
+    //! \brief Compute the residuals
+    //!
+    //! \param displacement variables
+    //! \param velocity time derivatives of the variables
+    //! \param residual vector containing the residuals
+    //! \param use_chemistry_rate if true compute residuals with chemistry rate
+    void compute_residuals(const Vector& displacement,
+                           const Vector& velocity,
+                           Vector& residual,
+                           bool use_chemistry_rate
+                           );
+    //! \brief Compute the residualts
+    //!
+    //! \param displacement variables
+    //! \param velocity time derivatives of the variables
+    //! \param residual vector containing the residuals
+    void compute_residuals(const Vector& displacement,
+                           const Vector& velocity,
+                           Vector& residual
+                           )
+    {
+        compute_residuals(displacement, velocity, residual, CHEMISTRY_RATE);
+    }
+
+
+    //! \brief Compute the residuals inside 'element' for 'component'
+    void residuals_element_component(
+            index_t element,
+            const Vector& displacement,
+            const Vector& velocity,
+            Eigen::Vector2d& element_residual,
+            bool use_chemistry_rate
+            );
+    //! \brief Compute the residuals inside 'element' for 'component'
+    void residuals_element_component(
+            index_t element,
+            const Vector& displacement,
+            const Vector& velocity,
+            Eigen::Vector2d& element_residual
+            )
+    {
+        residuals_element_component(element, displacement, velocity, element_residual, CHEMISTRY_RATE);
+    }
+
+    //! \brief Compute the jacobian
+    void compute_jacobian(Vector& displacement,
+                          Vector& velocity,
+                          Eigen::SparseMatrix<scalar_t>& jacobian,
+                          scalar_t alphadt
+                          );
+    //! \brief Compute the contribution of 'element' in the jacobian
+    void jacobian_element(
+            index_t element,
+            Vector& displacement,
+            Vector& velocity,
+            dfpm::list_triplet_t& jacobian,
+            scalar_t alphadt);
+
+    //! \brief Update the solutions
+    void update_solution(const Vector& update,
+                         scalar_t lambda,
+                         scalar_t alpha_dt,
+                         Vector& predictor,
+                         Vector& displacement,
+                         Vector& velocity);
+
+    //! \brief Apply boundary conditions to the velocity vector
+    //!
+    //! by default do nothing.
+    void apply_bc(scalar_t dt,
+                  const Vector& displacement,
+                  Vector& velocity)
+    {}
+
+
+private:
+    //! \brief Number the equations
+    void number_equations(std::vector<index_t> list_fixed_nodes,
+                          std::map<index_t, index_t>  list_slave_nodes
+                          );
+
+
+    index_t m_neq;     //!< Number fo equations
+    index_t m_ndf;     //!< Number of degree of freedoms per node
+    index_t m_tot_ndf; //!< Total number of degrees of freedoms
+
+    Eigen::Matrix<index_t, Eigen::Dynamic, 1> m_ideq; //!< equations indexes
+
+    mesh::Mesh1DPtr m_mesh;          //!< The mesh
+    SingleVariables* m_variables; //!< Variables, Don't own it, reset every timestep
+
+    bool m_is_in_residual_computation; //!< True if algo currently computing residuals
+};
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTPROGRAM_HPP
diff --git a/src/reactmicp/systems/single/transport_stagger.cpp b/src/reactmicp/systems/single/transport_stagger.cpp
new file mode 100644
index 0000000..fb81880
--- /dev/null
+++ b/src/reactmicp/systems/single/transport_stagger.cpp
@@ -0,0 +1,223 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#include "transport_stagger.hpp"
+
+#include "variables.hpp"
+#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
+#include "transport_program.hpp"
+
+#include "dfpm/solver/parabolic_driver.hpp"
+#include "reactmicp/solver/staggers_base/parabolic_driver_help.hpp"
+
+#include "specmicp_common/compat.hpp"
+
+#include <iostream>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+// Declaration of the implementation of the stagger
+struct SingleTransportStagger::SingleTransportStaggerImpl
+{
+    scalar_t m_dt {-1};
+    scalar_t m_residual_0 {-1};
+    SingleDiffusion m_program;
+    dfpmsolver::ParabolicDriver<SingleDiffusion> m_solver;
+
+    SingleTransportStaggerImpl(SingleVariables* variables,
+                               std::vector<index_t>& list_fixed_nodes):
+        m_program(variables, list_fixed_nodes),
+        m_solver(m_program)
+    {}
+
+    SingleTransportStaggerImpl(
+            SingleVariables* variables,
+            std::vector<index_t>& list_fixed_nodes,
+            std::map<index_t, index_t>&  list_slave_nodes):
+        m_program(variables,
+                  list_fixed_nodes,
+                  list_slave_nodes),
+        m_solver(m_program)
+    {}
+
+
+    void initialize_timestep(
+            scalar_t dt,
+            SingleVariables * const variables
+            );
+    solver::StaggerReturnCode restart_timestep(
+            SingleVariables * const variables
+            );
+    scalar_t get_residual(
+            SingleVariables * const var
+            );
+    scalar_t get_residual_0() {return m_residual_0;}
+
+    dfpmsolver::ParabolicDriverOptions& get_options() {
+        return m_solver.get_options();
+    }
+};
+
+SingleTransportStagger::SingleTransportStagger(
+        SingleVariablesPtr variables,
+        std::vector<index_t> list_fixed_nodes
+        ):
+    m_impl(make_unique<SingleTransportStaggerImpl>(
+               variables.get(),
+               list_fixed_nodes)
+           )
+{
+}
+
+SingleTransportStagger::SingleTransportStagger(
+        SingleVariablesPtr variables,
+        std::vector<index_t> list_fixed_nodes,
+        std::map<index_t, index_t>  list_slave_nodes):
+    m_impl(make_unique<SingleTransportStaggerImpl>(
+               variables.get(),
+               list_fixed_nodes,
+               list_slave_nodes)
+           )
+{
+}
+
+SingleTransportStagger::~SingleTransportStagger() = default;
+
+using VarConstPtr = SingleVariables * const; //!< Const pointer to the variables
+
+//! \brief Initialize the stagger at the beginning of an iteration
+void SingleTransportStagger::initialize_timestep(
+        scalar_t dt,
+        VariablesBase * const var
+        )
+{
+    VarConstPtr true_var = static_cast<VarConstPtr>(var);
+    m_impl->initialize_timestep(dt, true_var);
+}
+
+//! \brief Solve the equation for the timestep
+solver::StaggerReturnCode
+SingleTransportStagger::restart_timestep(VariablesBase* var)
+{
+    SingleVariables* true_var = static_cast<SingleVariables*>(var);
+    return m_impl->restart_timestep(true_var);
+}
+
+scalar_t SingleTransportStagger::get_update(VariablesBase * const var)
+{
+    return static_cast<VarConstPtr>(var)->velocity().norm();
+}
+
+//! \brief Compute the residuals norm
+scalar_t SingleTransportStagger::get_residual(VariablesBase * const var)
+{
+    VarConstPtr true_var = static_cast<VarConstPtr>(var);
+    auto res = m_impl->get_residual(true_var);
+
+    return res;
+}
+
+scalar_t SingleTransportStagger::get_residual_0(VariablesBase * const _)
+{
+    return m_impl->get_residual_0();
+}
+
+dfpmsolver::ParabolicDriverOptions& SingleTransportStagger::options_solver()
+{
+    return m_impl->get_options();
+}
+
+// ################ //
+//                  //
+//  Implementation  //
+//                  //
+// ###############  //
+
+void
+SingleTransportStagger::SingleTransportStaggerImpl::initialize_timestep(
+        scalar_t dt,
+        SingleVariables * const var
+        )
+{
+    m_dt = dt;
+    // Set the predictor
+    var->predictor() = var->displacement();
+
+    // Reset the values
+    var->velocity().setZero();
+    var->transport_rate().setZero();
+    m_program.set_variables(var);
+    // Initialize the DFPM solver
+    m_solver.initialize_timestep(dt, var->displacement());
+
+    // Compute the initial residual without external rate :
+    Eigen::VectorXd residuals;
+    m_program.compute_residuals(var->displacement(), var->velocity(), residuals, m_program.NO_CHEMISTRY_RATE);
+    m_residual_0 = residuals.norm();
+}
+
+//! \brief Solve the equation for the timestep
+solver::StaggerReturnCode
+SingleTransportStagger::SingleTransportStaggerImpl::restart_timestep(
+        SingleVariables * const var
+        )
+{
+    m_solver.velocity() = var->velocity();
+    dfpmsolver::ParabolicDriverReturnCode retcode = m_solver.restart_timestep(var->displacement());
+    // copy variables if successful
+    if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet)
+    {
+        var->velocity() = m_solver.velocity();
+    }
+    return solver::parabolic2StaggerReturnCode(retcode);
+}
+
+
+scalar_t
+SingleTransportStagger::SingleTransportStaggerImpl::get_residual(
+        SingleVariables * const var
+        )
+{
+    Eigen::VectorXd residuals;
+    m_program.compute_residuals(var->displacement(), var->velocity(), residuals);
+    return residuals.norm();
+}
+
+} // end namespace satdiff
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
diff --git a/src/reactmicp/systems/single/transport_stagger.hpp b/src/reactmicp/systems/single/transport_stagger.hpp
new file mode 100644
index 0000000..d11545d
--- /dev/null
+++ b/src/reactmicp/systems/single/transport_stagger.hpp
@@ -0,0 +1,104 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTSTAGGER_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTSTAGGER_HPP
+
+#include "variablesfwd.hpp"
+#include "reactmicp/solver/staggers_base/transport_stagger_base.hpp"
+#include "dfpm/solver/parabolic_structs.hpp"
+
+#include <vector>
+#include <map>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+using VariablesBase = solver::VariablesBase;
+
+//! \brief The transport stagger for the single system
+class SPECMICP_DLL_PUBLIC SingleTransportStagger: public solver::TransportStaggerBase
+{
+public:
+    //! \brief Simple Constructor
+    SingleTransportStagger(SingleVariablesPtr variables,
+                           std::vector<index_t> list_fixed_nodes);
+
+    //! \brief Constructor
+    //!
+    //! \param variables the variables
+    //! \param list_fixed_nodes list of fixed nodes
+    //! \param list_slave_nodes list of slave nodes (same equation)
+    SingleTransportStagger(SingleVariablesPtr variables,
+                              std::vector<index_t> list_fixed_nodes,
+                              std::map<index_t, index_t>  list_slave_nodes
+                           );
+
+    ~SingleTransportStagger();
+    SingleTransportStagger& operator= (const SingleTransportStagger& other) = delete;
+    SingleTransportStagger(const SingleTransportStagger& other) = delete;
+
+    //! \brief Return the options of the solver
+    dfpmsolver::ParabolicDriverOptions& options_solver();
+
+    //! \brief Initialize the stagger at the beginning of an iteration
+    void initialize_timestep(scalar_t dt, VariablesBase * const var) override;
+
+    //! \brief Solve the equation for the timestep
+    solver::StaggerReturnCode restart_timestep(VariablesBase * const var) override;
+
+    //! \brief Compute the residuals norm
+    scalar_t get_residual(VariablesBase * const var) override;
+    //! \brief Compute the residuals norm
+    scalar_t get_residual_0(VariablesBase * const var) override;
+
+    //! \brief Obtain the update
+    scalar_t get_update(VariablesBase * const var) override;
+
+private:
+    //! \brief Implementation structure
+    struct SingleTransportStaggerImpl;
+    std::unique_ptr<SingleTransportStaggerImpl> m_impl; //!< Implementation
+
+};
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_TRANSPORTSTAGGER_HPP
diff --git a/src/reactmicp/systems/single/variables.cpp b/src/reactmicp/systems/single/variables.cpp
new file mode 100644
index 0000000..b3b8d78
--- /dev/null
+++ b/src/reactmicp/systems/single/variables.cpp
@@ -0,0 +1,61 @@
+/* =============================================================================
+
+ Copyright (c) 2014 - 2017
+ F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#include "reactmicp/systems/single/variables.hpp"
+
+#include "dfpm/mesh.hpp"
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+SingleVariables::SingleVariables(mesh::Mesh1DPtr the_mesh):
+    m_mesh(the_mesh)
+{
+    reset_main_variables();
+}
+
+void SingleVariables::reset_main_variables() {
+    displacement() = predictor();
+    velocity().setZero();
+    transport_rate().setZero();
+    chemistry_rate().setZero();
+}
+
+} // namespace
+} // namespace systems
+} // namespace reactmicp
+} // namespace specmicp
diff --git a/src/reactmicp/systems/single/variables.hpp b/src/reactmicp/systems/single/variables.hpp
new file mode 100644
index 0000000..fdce62a
--- /dev/null
+++ b/src/reactmicp/systems/single/variables.hpp
@@ -0,0 +1,226 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SINGLE_VARIABLES_HPP
+#define SPECMICP_REACTMICP_SINGLE_VARIABLES_HPP
+
+#include <vector>
+
+#include "specmicp_common/types.hpp"
+#include "reactmicp/solver/staggers_base/variables_base.hpp"
+#include "specmicp_database/database_fwd.hpp"
+
+#include "dfpm/meshes/mesh1dfwd.hpp"
+namespace specmicp {
+namespace reactmicp {
+namespace solver {
+    using VariablesBasePtr = std::shared_ptr<VariablesBase>;
+}
+} // end namespace reactmicp
+} // end namespace specmicp
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+
+//! \brief Variables for the saturated reactive transport system
+//!
+//! Contain all the variables that need to be shared between the staggers
+class SPECMICP_DLL_PUBLIC SingleVariables: public solver::VariablesBase
+{
+    // SingleVariablesFactory should be the class to use to inialize
+    // the variables correctly
+    friend class SingleVariablesFactory;
+
+public:
+    //! \brief Constructor
+    SingleVariables(
+            mesh::Mesh1DPtr the_mesh);
+
+    //! \brief Return the mesh
+    mesh::Mesh1DPtr get_mesh() {return m_mesh;}
+
+    //! \brief Return the number of components
+    index_t nb_component() {return 1;}
+    //! \brief Return the number of nodes
+    index_t nb_nodes() {return m_is_fixed_composition.size();}
+    //! \brief Return true if 'node' has a fixed composition
+    index_t is_fixed_composition(index_t node) {return m_is_fixed_composition[node];}
+    // Main variables
+    // ==============
+
+    //! \brief Return the main variable vector
+    Vector& displacement() {return m_displacement;}
+    //! \brief Return the main variable vector at the beginning of the timestep
+    Vector& predictor() {return m_predictor;}
+    //! \brief Return the velocity of the main variables
+    Vector& velocity() {return m_velocity;}
+    //! \brief Return the rate of change of the main variables due to the transport operator
+    Vector& transport_rate() {return m_transport_rate;}
+    //! \brief Return the rate of change of the main variables due to the chemistry operator
+    Vector& chemistry_rate() {return m_chemistry_rate;}
+
+    // Access to main variables
+    // ========================
+
+    //! \brief Return the number of degree of freedom (per node) in the main variables vector
+    index_t ndf() {return 2;} // fluid_conc and bound_conc
+    //! \brief Return the offset of 'node' in the main variables vector
+    index_t offset_node(index_t node) {return node*ndf();}
+    //! \brief Return the offset of the aqueous concentration variables in the main variables vector
+    index_t offset_fluid_concentration() {return 0;}
+    //! \brief Return the offset of the aqueous concentrations variables in the main variables vector
+    index_t offset_fluid_concentration(index_t node) {
+        return offset_fluid_concentration()+offset_node(node);
+    }
+    //! \brief Return the offset of the bound concentration variables in the main variables vector
+    index_t offset_bound_concentration() {return 1;}
+    //! \brief Return the offset of the bound concentrations variables in the main variables vector
+    index_t offset_bound_concentration(index_t node) {
+        return offset_bound_concentration()+offset_node(node);
+    }
+    //! \brief Return the degree of freedom number for the fluid concentration at 'node'
+    index_t dof_fluid_concentration(index_t node) {
+        return offset_fluid_concentration(node);
+    }
+    //! \brief Return the degree of freedom number for the bound concentration at 'node'
+    index_t dof_bound_concentration(index_t node) {
+        return offset_bound_concentration(node);
+    }
+    //! \brief Return the fluid concentration at 'node' in 'var'
+    //!
+    //! 'var' is any of the main variables vector, it may be a velocity vector
+    scalar_t& fluid_concentration(index_t node, Vector& var) {
+        return var(dof_fluid_concentration(node));
+    }
+    //! \brief Return the fluid concentration at 'node' in main variables
+    scalar_t& fluid_concentration(index_t node) {
+        return m_displacement(dof_fluid_concentration(node));
+    }
+    //! \brief Return the bound concentration  at 'node' in 'var'
+    //!
+    //! 'var' is any of the main variables vector, it may be a velocity vector
+    scalar_t& bound_concentration(index_t node, Vector& var){
+        return var(dof_bound_concentration(node));
+    }
+    //! \brief Return the solid concentration at 'node' in main variables
+    //!
+    //! 'var' is any of the main variables vector, it may be a velocity vector
+    scalar_t& bound_concentration(index_t node){
+        return m_displacement(dof_bound_concentration(node));
+    }
+    //! \brief Return a vector containing the total concentrations computed from the main variables
+    //!
+    //! This is to be used to restart the chemistry computation
+    Vector total_concentrations(index_t node);
+
+    // Chemistry variables
+    // ===================
+
+    //! \brief Return the chemistry variables
+    Vector chemistry_variables() {return m_chemistry_variables;}
+
+    // Upscaling
+    // =========
+    //! \brief Return the offset for 'node' in the upscaling variables vector
+    index_t offset_node_upscaling(index_t node) {return ndf_upscaling()*node;}
+
+    //! \brief Return the number fo degree of freedom (per node) for the upscaling vector
+    index_t ndf_upscaling() {return 5;}
+    //! \brief Refturn the dof for the diffusion coefficient at 'node'
+    index_t dof_diffusion_coefficient(index_t node) {return 0 + offset_node_upscaling(node);}
+    //! \brief Return the degree of freedom for the porosity at 'node'
+    index_t dof_porosity(index_t node) {return 1 + offset_node_upscaling(node);}
+    //! \brief Return the degree of freedom of the porosity velocity at 'node'
+    index_t dof_vel_porosity(index_t node) {return 2 + offset_node_upscaling(node);}
+    //! \brief Return the degree of freedom of the saturation at 'node'
+    index_t dof_saturation(index_t node) {return 3 + offset_node_upscaling(node);}
+    //! \brief Return the degree of freedom of the saturation velocity at 'node'
+    index_t dof_vel_saturation(index_t node) {return 4 + offset_node_upscaling(node);}
+
+    //! \brief Return the diffusion coefficient at 'node'
+    scalar_t& diffusion_coefficient(index_t node) {return m_upscaling(dof_diffusion_coefficient(node));}
+    //! \brief Return the porosity at 'node'
+    scalar_t& porosity(index_t node) {return m_upscaling(dof_porosity(node));}
+    //! \brief Return the rate of change of the porosity at 'node'
+    scalar_t& vel_porosity(index_t node) {return m_upscaling(dof_vel_porosity(node));}
+    //! \brief Return the saturation at 'node'
+    scalar_t& saturation(index_t node) {return m_upscaling(dof_saturation(node));}
+    //! \brief Return the saturation velocity at 'node'
+    scalar_t& vel_saturation(index_t node) {return m_upscaling(dof_vel_saturation(node));}
+
+    //! \brief Return the vector of upscaling variables
+    Vector& upscaling_variables() {return m_upscaling;}
+
+
+    //! \brief Reset the main variables
+    void reset_main_variables() override;
+
+private:
+
+    // ############ //
+    //  Attributes  //
+    // ############ //
+
+    mesh::Mesh1DPtr m_mesh;
+    std::vector<bool> m_is_fixed_composition;
+
+    // Main variables
+    // ==============
+
+    Vector m_displacement;
+    Vector m_predictor;
+    Vector m_velocity;
+    Vector m_transport_rate;
+    Vector m_chemistry_rate;
+
+    // Chemistry variables
+    // ===================
+
+    Vector m_chemistry_variables;
+
+    // Upscaling
+    // =========
+
+    Vector m_upscaling;
+};
+
+} // namespace single
+} // namespace systems
+} // namespace reactmicp
+} // namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SINGLE_VARIABLES_HPP
diff --git a/src/reactmicp/systems/single/variablesfwd.hpp b/src/reactmicp/systems/single/variablesfwd.hpp
new file mode 100644
index 0000000..6fee795
--- /dev/null
+++ b/src/reactmicp/systems/single/variablesfwd.hpp
@@ -0,0 +1,55 @@
+/* =============================================================================
+
+ Copyright (c)
+ 2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
+ 2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
+ 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. *
+
+============================================================================= */
+
+#ifndef SPECMICP_REACTMICP_SYSTEMS_SINGLE_VARIABLESFWD_HPP
+#define SPECMICP_REACTMICP_SYSTEMS_SINGLE_VARIABLESFWD_HPP
+
+#include <memory>
+
+namespace specmicp {
+namespace reactmicp {
+namespace systems {
+namespace single {
+
+class SingleVariables;
+using SingleVariablesPtr = std::shared_ptr<SingleVariables>;
+
+} // end namespace single
+} // end namespace systems
+} // end namespace reactmicp
+} // end namespace specmicp
+
+#endif // SPECMICP_REACTMICP_SYSTEMS_SINGLE_VARIABLESFWD_HPP
diff --git a/tests/reactmicp/CMakeLists.txt b/tests/reactmicp/CMakeLists.txt
index 2e2e25d..9d3e693 100644
--- a/tests/reactmicp/CMakeLists.txt
+++ b/tests/reactmicp/CMakeLists.txt
@@ -1,102 +1,126 @@
 
 
 # Reactmicp : Reactive transport solver
 # -------------------------------------
 
 set(test_reactmicp_files
     solver/test_reactive_transport_solver.cpp
     solver/test_coupling.cpp
 )
 add_catch_test(NAME reactmicp SOURCES ${test_reactmicp_files}
  LINK_LIBRARIES ${REACTMICP_STATIC_LIBS})
 
 
 # Saturated diffusion using new reactive transport solver
 # =======================================================
 
 set(test_reactmicp_saturated_files
     systems/saturated_react/test_reactmicp_saturated_react.cpp
     systems/saturated_react/speciation_system.cpp
     systems/saturated_react/variables.cpp
     systems/saturated_react/transport_program.cpp
     systems/saturated_react/equilibrium_stagger.cpp
 )
 if ( REACTMICP_ENABLE_SYSTEM_SATURATED )
     add_catch_test(
         NAME reactmicp_saturated_react
         SOURCES ${test_reactmicp_saturated_files}
         LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}
     )
 endif( REACTMICP_ENABLE_SYSTEM_SATURATED )
 
 # Unsaturated reactive transport system
 # =====================================
 
 set(UNSATURATED_DIR systems/unsaturated)
 
 set(test_reactmicp_unsaturated_files
     ${UNSATURATED_DIR}/test_reactmicp_unsaturated.cpp
 
     ${UNSATURATED_DIR}/aqueous_equation.cpp
     ${UNSATURATED_DIR}/aqueous_pressure_equation.cpp
     ${UNSATURATED_DIR}/boundary_conditions.cpp
     ${UNSATURATED_DIR}/configuration.cpp
     ${UNSATURATED_DIR}/equilibrium_stagger.cpp
     ${UNSATURATED_DIR}/hdf5_unsaturated.cpp
     ${UNSATURATED_DIR}/pressure_equation.cpp
     ${UNSATURATED_DIR}/saturation_equation.cpp
     ${UNSATURATED_DIR}/saturation_pressure_equation.cpp
     ${UNSATURATED_DIR}/transport_stagger.cpp
     ${UNSATURATED_DIR}/variables.cpp
 )
 
 set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\")
 set_source_files_properties(
     ${test_reactmicp_unsaturated_files}
      PROPERTIES COMPILE_DEFINITIONS "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}"
 )
 
 if ( REACTMICP_ENABLE_SYSTEM_UNSATURATED )
 add_catch_test(
     NAME reactmicp_unsaturated
     SOURCES ${test_reactmicp_unsaturated_files}
     LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}
 )
 
 # unsaturated binary
 add_subdirectory( binary/unsaturated )
 
 endif ( REACTMICP_ENABLE_SYSTEM_UNSATURATED )
 
 # Chloride reactive transport system
 # ===================================
 
 set(CHLORIDE_DIR systems/chloride)
 
 set(test_reactmicp_chloride_files
     ${CHLORIDE_DIR}/test_reactmicp_chloride.cpp
 
     ${CHLORIDE_DIR}/boundary_conditions.cpp
     ${CHLORIDE_DIR}/variables.cpp
     ${CHLORIDE_DIR}/transport_program.cpp
     ${CHLORIDE_DIR}/transport_stagger.cpp
     ${CHLORIDE_DIR}/hdf5_chloride.cpp
 )
 
 set(TEST_CEMDATA_PATH \"../../data/cemdata.yaml\")
 set_source_files_properties(
     ${test_reactmicp_chloride_files}
      PROPERTIES COMPILE_DEFINITIONS "TEST_CEMDATA_PATH=${TEST_CEMDATA_PATH}"
 )
 
 if ( REACTMICP_ENABLE_SYSTEM_CHLORIDE )
 add_catch_test(
     NAME reactmicp_chloride
     SOURCES ${test_reactmicp_chloride_files}
     LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}
 )
 
 # chloride binary
 add_subdirectory( binary/chloride )
 
 endif ( REACTMICP_ENABLE_SYSTEM_CHLORIDE )
+
+# Single-reactive system
+# ======================
+
+set(SINGLE_DIR systems/single)
+
+set(test_reactmicp_single_files
+    ${SINGLE_DIR}/test_reactmicp_single.cpp
+
+    ${SINGLE_DIR}/variables.cpp
+    ${SINGLE_DIR}/transport_program.cpp
+    ${SINGLE_DIR}/chemistry.cpp
+)
+
+if( REACTMICP_ENABLE_SYSTEM_SINGLE )
+    add_catch_test(
+        NAME reactmicp_single
+        SOURCES ${test_reactmicp_single_files}
+        LINK_LIBRARIES ${REACTMICP_STATIC_LIBS}
+
+        )
+
+endif ( REACTMICP_ENABLE_SYSTEM_SINGLE )
+
diff --git a/tests/reactmicp/systems/single/chemistry.cpp b/tests/reactmicp/systems/single/chemistry.cpp
new file mode 100644
index 0000000..1861958
--- /dev/null
+++ b/tests/reactmicp/systems/single/chemistry.cpp
@@ -0,0 +1,74 @@
+#include "catch.hpp"
+
+#include "dfpm/mesh.hpp"
+
+#include "reactmicp/systems/single/variables.hpp"
+#include "reactmicp/systems/single/init_variables.hpp"
+#include "reactmicp/systems/single/chemistry_program.hpp"
+#include "reactmicp/systems/single/chemistry_stagger.hpp"
+
+#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
+#include <memory>
+
+using namespace specmicp;
+using namespace specmicp::reactmicp::systems::single;
+
+scalar_t rate_function_ex(index_t node, scalar_t dt, SingleVariables * const var) {
+
+    return 0.2*(var->fluid_concentration(node)-var->bound_concentration(node)/10);
+
+}
+
+
+TEST_CASE("Single default chemistry program", "[SingleReact, chemistry, stagger]") {
+    specmicp::mesh::Uniform1DMeshGeometry geom;
+    geom.nb_nodes = 5;
+    geom.dx = 0.1;
+    geom.section = 5.0;
+    specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(geom);
+
+    std::vector<specmicp::index_t> list_fixed_nodes = {0, };
+
+    specmicp::Vector init(5);
+    init << 10,1,1,1,1;
+
+    specmicp::Vector solid_init(5);
+    solid_init << 0, 0, 0, 0, 0;
+
+    auto variables =
+            specmicp::reactmicp::systems::single::init_variables(
+                the_mesh,
+                list_fixed_nodes,
+                init,
+                solid_init
+                );
+
+    for (specmicp::index_t node: the_mesh->range_nodes())
+    {
+        variables->porosity(node) = 1;
+        variables->saturation(node) = 1;
+        variables->diffusion_coefficient(node) = 2;
+    }
+
+    SECTION("Program") {
+
+        auto program = DefaultChemistryProgram(rate_function_ex);
+
+        auto ret = program.solve_one_node(1, 1, variables.get());
+
+        CHECK(ret.retcode > ReturnCode::NotFinishedYet);
+        CHECK(ret.fluid_concentration == Approx(0.8));
+        CHECK(ret.bound_concentration == Approx(0.2));
+    }
+
+    SECTION("Stagger") {
+
+        std::shared_ptr<ChemistryProgram> program = std::make_shared<DefaultChemistryProgram>(rate_function_ex);
+        auto stagger= ChemistryStagger(program);
+        stagger.initialize_timestep(1, variables.get());
+
+        auto ret = stagger.restart_timestep(variables.get());
+        REQUIRE(ret > specmicp::reactmicp::solver::StaggerReturnCode::NotConvergedYet);
+    }
+
+}
diff --git a/tests/reactmicp/systems/single/test_reactmicp_single.cpp b/tests/reactmicp/systems/single/test_reactmicp_single.cpp
new file mode 100644
index 0000000..178916e
--- /dev/null
+++ b/tests/reactmicp/systems/single/test_reactmicp_single.cpp
@@ -0,0 +1,3 @@
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
diff --git a/tests/reactmicp/systems/single/transport_program.cpp b/tests/reactmicp/systems/single/transport_program.cpp
new file mode 100644
index 0000000..f963999
--- /dev/null
+++ b/tests/reactmicp/systems/single/transport_program.cpp
@@ -0,0 +1,78 @@
+#include "catch.hpp"
+
+#include "dfpm/mesh.hpp"
+#include "dfpm/solver/parabolic_driver.hpp"
+
+#include "reactmicp/systems/single/variables.hpp"
+#include "reactmicp/systems/single/init_variables.hpp"
+#include "reactmicp/systems/single/transport_program.hpp"
+
+#include "reactmicp/systems/single/transport_stagger.hpp"
+
+TEST_CASE("Single transport program", "[SingleReact, transport, program, stagger]") {
+    specmicp::mesh::Uniform1DMeshGeometry geom;
+    geom.nb_nodes = 5;
+    geom.dx = 0.1;
+    geom.section = 5.0;
+    specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(geom);
+
+    std::vector<specmicp::index_t> list_fixed_nodes = {0, };
+
+    specmicp::Vector init(5);
+    init << 1,0,0,0,0;
+
+    specmicp::Vector solid_init(5);
+    solid_init << 0, 0, 0, 0, 0;
+
+    auto variables =
+            specmicp::reactmicp::systems::single::init_variables(
+                the_mesh,
+                list_fixed_nodes,
+                init,
+                solid_init
+                );
+
+    for (specmicp::index_t node: the_mesh->range_nodes())
+    {
+        variables->porosity(node) = 0.2;
+        variables->diffusion_coefficient(node) = 2;
+    }
+
+    SECTION("Initialisation") {
+
+        auto program = specmicp::reactmicp::systems::single::SingleDiffusion(variables.get(), list_fixed_nodes);
+
+        specmicp::Vector residuals;
+
+        REQUIRE(program.get_ndf() == 2); // 2*nb_comp
+        REQUIRE(program.get_tot_ndf() == 10); // 2*nb_comp*nb_nodes
+        REQUIRE(program.get_neq() == 4); // nb_nodes*active_comp
+
+        program.compute_residuals(variables->displacement(), variables->velocity(), residuals);
+
+    }
+
+    SECTION("Solving the problem manually") {
+        specmicp::reactmicp::systems::single::SingleDiffusion program(variables.get(), list_fixed_nodes);
+
+        specmicp::dfpmsolver::ParabolicDriver<specmicp::reactmicp::systems::single::SingleDiffusion> solver(program);
+
+        variables->predictor() = variables->displacement();
+
+         int retcode =  (int) solver.solve_timestep(100, variables->displacement());
+
+         REQUIRE(retcode > 0);
+    }
+
+
+    SECTION("Transport stagger") {
+        specmicp::reactmicp::systems::single::SingleTransportStagger stagger(variables, list_fixed_nodes);
+
+        stagger.initialize_timestep(100, variables.get());
+
+        auto retcode =  stagger.restart_timestep(variables.get());
+
+        REQUIRE((int) retcode > 0);
+
+    }
+}
diff --git a/tests/reactmicp/systems/single/variables.cpp b/tests/reactmicp/systems/single/variables.cpp
new file mode 100644
index 0000000..bc37fba
--- /dev/null
+++ b/tests/reactmicp/systems/single/variables.cpp
@@ -0,0 +1,49 @@
+#include "catch.hpp"
+
+#include "dfpm/mesh.hpp"
+
+#include "reactmicp/systems/single/variables.hpp"
+#include "reactmicp/systems/single/init_variables.hpp"
+
+
+TEST_CASE("Single Variables test", "[Single, variables, initialisation]") {
+
+    SECTION("Creation") {
+
+        specmicp::mesh::Uniform1DMeshGeometry geom;
+        geom.nb_nodes = 5;
+        geom.dx = 0.1;
+        geom.section = 5.0;
+        specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(geom);
+
+
+        auto vars = specmicp::reactmicp::systems::single::SingleVariables(the_mesh);
+
+        REQUIRE( vars.get_mesh()->nb_nodes() == 5);
+    }
+
+    SECTION("Initialisation") {
+        specmicp::mesh::Uniform1DMeshGeometry geom;
+        geom.nb_nodes = 5;
+        geom.dx = 0.1;
+        geom.section = 5.0;
+        specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(geom);
+
+        std::vector<specmicp::index_t> fixed_nodes = {1,};
+
+        specmicp::Vector liquid_conc(5);
+        liquid_conc << 1,2,3,4,5;
+        specmicp::Vector solid_conc(5);
+        solid_conc << 6,7,8,9,10;
+
+        auto vars = specmicp::reactmicp::systems::single::init_variables(
+                    the_mesh, fixed_nodes, liquid_conc, solid_conc);
+
+
+        REQUIRE( vars->get_mesh()->nb_nodes() == 5);
+        REQUIRE( vars->chemistry_rate()(4)== 0);
+        REQUIRE( vars->fluid_concentration(2) == 3);
+        REQUIRE( vars->bound_concentration(4) == 10);
+    }
+
+}