diff --git a/demos/FiniteVolume/AMR_Burgers_Hat.cpp b/demos/FiniteVolume/AMR_Burgers_Hat.cpp
index 0a975c4e8..2d2310ae0 100644
--- a/demos/FiniteVolume/AMR_Burgers_Hat.cpp
+++ b/demos/FiniteVolume/AMR_Burgers_Hat.cpp
@@ -197,7 +197,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the Burgers equation in 2d using AMR", argc, argv);
 
     constexpr std::size_t dim = 1; // cppcheck-suppress unreadVariable
-    using Config              = samurai::amr::Config<dim>;
 
     // Simulation parameters
     double left_box  = -3;
@@ -208,10 +207,7 @@ int main(int argc, char* argv[])
     std::string restart_file;
 
     // AMR parameters
-    std::size_t start_level = 7;
-    std::size_t min_level   = 2;
-    std::size_t max_level   = 7;
-    bool correction         = false;
+    bool correction = false;
 
     // Output parameters
     fs::path path        = fs::current_path();
@@ -224,9 +220,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--start-level", start_level, "Start level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--min-level", min_level, "Minimum level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--max-level", max_level, "Maximum level of AMR")->capture_default_str()->group("AMR parameters");
     app.add_option("--with-correction", correction, "Apply flux correction at the interface of two refinement levels")
         ->capture_default_str()
         ->group("AMR parameters");
@@ -236,12 +229,13 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box({left_box}, {right_box});
-    samurai::amr::Mesh<Config> mesh;
-    auto phi = samurai::make_scalar_field<double>("phi", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(7).max_stencil_size(2).start_level(7).disable_minimal_ghost_width();
+    auto mesh   = samurai::amr::make_empty_Mesh(config);
+    auto phi    = samurai::make_scalar_field<double>("phi", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, start_level, min_level, max_level};
+        mesh = samurai::amr::make_Mesh(config, box);
         init_solution(phi);
     }
     else
@@ -256,7 +250,7 @@ int main(int argc, char* argv[])
     auto tag = samurai::make_scalar_field<int>("tag", mesh);
     const xt::xtensor_fixed<int, xt::xshape<2, 1>> stencil_grad{{1}, {-1}};
 
-    const double dx      = mesh.cell_length(max_level);
+    const double dx      = mesh.min_cell_length();
     double dt            = 0.99 * dx;
     const double dt_save = Tf / static_cast<double>(nfiles);
 
diff --git a/demos/FiniteVolume/BZ/bz_2d.cpp b/demos/FiniteVolume/BZ/bz_2d.cpp
index 05d7efd64..8e430d573 100644
--- a/demos/FiniteVolume/BZ/bz_2d.cpp
+++ b/demos/FiniteVolume/BZ/bz_2d.cpp
@@ -289,8 +289,8 @@ int main()
     const double regularity = 1.;    // Regularity guess for multiresolution
     const double epsilon_MR = 2.e-4; // Threshold used by multiresolution
 
-    using Config = samurai::MRConfig<dim, 2>;
-    samurai::MRMesh<Config> mesh(box, min_level, max_level); // Constructing mesh from the box
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).max_stencil_radius(2);
+    auto mesh   = samurai::make_MRMesh(config, box); // Constructing mesh from the box
 
     using mesh_id_t = typename decltype(mesh)::mesh_id_t;
 
diff --git a/demos/FiniteVolume/advection_1d.cpp b/demos/FiniteVolume/advection_1d.cpp
index 3df099cc7..8cf3b2003 100644
--- a/demos/FiniteVolume/advection_1d.cpp
+++ b/demos/FiniteVolume/advection_1d.cpp
@@ -62,7 +62,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the advection equation in 1d using multiresolution", argc, argv);
 
     constexpr std::size_t dim = 1; // cppcheck-suppress unreadVariable
-    using Config              = samurai::MRConfig<dim>;
 
     // Simulation parameters
     double left_box  = -2;
@@ -74,10 +73,6 @@ int main(int argc, char* argv[])
     double t         = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 6;
-    std::size_t max_level = 12;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "FV_advection_1d";
@@ -91,8 +86,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -100,12 +93,14 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box({left_box}, {right_box});
-    samurai::MRMesh<Config> mesh;
-    auto u = samurai::make_scalar_field<double>("u", mesh);
+
+    auto config = samurai::mesh_config<dim>().min_level(6).max_level(12).periodic(is_periodic).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh = samurai::make_empty_MRMesh(config);
+    auto u    = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level, std::array<bool, dim>{is_periodic}};
+        mesh = samurai::make_MRMesh(config, box);
         init(u);
     }
     else
@@ -113,7 +108,7 @@ int main(int argc, char* argv[])
         samurai::load(restart_file, mesh, u);
     }
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     if (!is_periodic)
diff --git a/demos/FiniteVolume/advection_2d.cpp b/demos/FiniteVolume/advection_2d.cpp
index 7066e7119..a3cbbdd91 100644
--- a/demos/FiniteVolume/advection_2d.cpp
+++ b/demos/FiniteVolume/advection_2d.cpp
@@ -63,7 +63,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the advection equation in 2d using multiresolution", argc, argv);
 
     constexpr std::size_t dim = 2;
-    using Config              = samurai::MRConfig<dim>;
 
     // Simulation parameters
     xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
@@ -76,10 +75,6 @@ int main(int argc, char* argv[])
     double t   = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 4;
-    std::size_t max_level = 10;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "FV_advection_2d";
@@ -92,8 +87,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -101,12 +94,13 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
-    samurai::MRMesh<Config> mesh;
-    auto u = samurai::make_scalar_field<double>("u", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(10).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
+    auto u      = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         init(u);
     }
     else
@@ -115,7 +109,7 @@ int main(int argc, char* argv[])
     }
     samurai::make_bc<samurai::Dirichlet<1>>(u, 0.);
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
diff --git a/demos/FiniteVolume/advection_2d_user_bc.cpp b/demos/FiniteVolume/advection_2d_user_bc.cpp
index ac4d5e163..301bf3361 100644
--- a/demos/FiniteVolume/advection_2d_user_bc.cpp
+++ b/demos/FiniteVolume/advection_2d_user_bc.cpp
@@ -95,7 +95,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the advection equation in 2d using multiresolution", argc, argv);
 
     constexpr std::size_t dim = 2;
-    using Config              = samurai::MRConfig<dim>;
 
     // Simulation parameters
     xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
@@ -108,10 +107,6 @@ int main(int argc, char* argv[])
     double t   = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 4;
-    std::size_t max_level = 10;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "FV_advection_2d";
@@ -124,8 +119,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -133,12 +126,13 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
-    samurai::MRMesh<Config> mesh;
-    auto u = samurai::make_scalar_field<double>("u", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(10).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
+    auto u      = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         init(u);
     }
     else
@@ -147,7 +141,7 @@ int main(int argc, char* argv[])
     }
     samurai::make_bc<Mybc>(u, 0.);
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
diff --git a/demos/FiniteVolume/burgers.cpp b/demos/FiniteVolume/burgers.cpp
index c965fe9eb..8d5506322 100644
--- a/demos/FiniteVolume/burgers.cpp
+++ b/demos/FiniteVolume/burgers.cpp
@@ -52,7 +52,6 @@ int main_dim(int argc, char* argv[])
 {
     auto& app = samurai::initialize("Finite volume example for the Burgers equation", argc, argv);
 
-    using Config  = samurai::MRConfig<dim, 3>;
     using Box     = samurai::Box<double, dim>;
     using point_t = typename Box::point_t;
 
@@ -91,8 +90,6 @@ int main_dim(int argc, char* argv[])
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -108,7 +105,9 @@ int main_dim(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh;
+
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).max_stencil_size(6);
+    auto mesh   = samurai::make_empty_MRMesh(config);
 
     auto u    = samurai::make_vector_field<n_comp>("u", mesh);
     auto u1   = samurai::make_vector_field<n_comp>("u1", mesh);
@@ -117,7 +116,7 @@ int main_dim(int argc, char* argv[])
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         u.resize();
 
         // Initial solution
@@ -221,7 +220,7 @@ int main_dim(int argc, char* argv[])
     //   Time iteration   //
     //--------------------//
 
-    double dx = mesh.cell_length(max_level);
+    double dx = mesh.min_cell_length();
     dt        = cfl * dx / pow(2, dim);
 
     auto MRadaptation = samurai::make_MRAdapt(u);
diff --git a/demos/FiniteVolume/burgers_mra.cpp b/demos/FiniteVolume/burgers_mra.cpp
index ea6e74dbe..5e92301c4 100644
--- a/demos/FiniteVolume/burgers_mra.cpp
+++ b/demos/FiniteVolume/burgers_mra.cpp
@@ -226,7 +226,6 @@ void run_simulation(Field& u,
 int main(int argc, char* argv[])
 {
     constexpr std::size_t dim = 1;
-    using Config              = samurai::MRConfig<dim, 2, 2>;
     using Box                 = samurai::Box<double, dim>;
 
     auto& app = samurai::initialize("Finite volume example for the Burgers equation in 1d", argc, argv);
@@ -245,10 +244,6 @@ int main(int argc, char* argv[])
     // Time integration
     double cfl = 0.95;
 
-    // Multiresolution parameters
-    std::size_t min_level = 1;
-    std::size_t max_level = 7;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "burgers_mra";
@@ -256,8 +251,6 @@ int main(int argc, char* argv[])
 
     app.add_option("--init-sol", init_sol, "Initial solution: hat/gaussian")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Ouput");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Ouput");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Ouput");
@@ -270,8 +263,10 @@ int main(int argc, char* argv[])
     //--------------------//
 
     Box box({left_box}, {right_box});
-    samurai::MRMesh<Config> mesh{box, min_level, max_level};
-    samurai::MRMesh<Config> max_level_mesh{box, max_level, max_level};
+    auto config           = samurai::mesh_config<dim>().min_level(1).max_level(7).max_stencil_radius(2).graduation_width(2);
+    auto mesh             = samurai::make_MRMesh(config, box);
+    auto max_level_config = samurai::mesh_config<dim>().min_level(mesh.max_level()).max_level(mesh.max_level()).max_stencil_radius(2);
+    auto max_level_mesh   = samurai::make_MRMesh(max_level_config, box);
 
     auto u    = samurai::make_scalar_field<double>("u", mesh);
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
diff --git a/demos/FiniteVolume/burgers_os.cpp b/demos/FiniteVolume/burgers_os.cpp
index e9965ec1c..3e51be9cf 100644
--- a/demos/FiniteVolume/burgers_os.cpp
+++ b/demos/FiniteVolume/burgers_os.cpp
@@ -90,9 +90,9 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the linear convection equation", argc, argv);
 
     static constexpr std::size_t dim = 1;
-    using Config                     = samurai::MRConfig<dim, 2, 2>;
-    using Box                        = samurai::Box<double, dim>;
-    using point_t                    = typename Box::point_t;
+
+    using Box     = samurai::Box<double, dim>;
+    using point_t = typename Box::point_t;
 
     std::cout << "------------------------- Burgers 1D -------------------------" << std::endl;
 
@@ -125,8 +125,6 @@ int main(int argc, char* argv[])
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -134,8 +132,6 @@ int main(int argc, char* argv[])
 
     SAMURAI_PARSE(argc, argv);
 
-    std::cout << "  max_level = " << max_level << "   min_level = " << min_level << std::endl;
-
     //--------------------//
     // Problem definition //
     //--------------------//
@@ -148,13 +144,18 @@ int main(int argc, char* argv[])
     Box box_left(box_corner1, 0.5 * (box_corner1 + box_corner2));
     Box box_right(0.5 * (box_corner1 + box_corner2), box_corner2);
 
-    samurai::LevelCellArray<1> lca_left(max_level, box_left, {-1}, 0.05, 2);
-    samurai::LevelCellArray<1> lca_right(max_level, box_right, {-1}, 0.05, 2);
+    samurai::LevelCellArray<1> lca_left(max_level, box_left, {-1}, -1., 2);
+    samurai::LevelCellArray<1> lca_right(max_level, box_right, {-1}, -1., 2);
 
-    std::array<bool, dim> periodic;
-    periodic.fill(true);
-    samurai::MRMesh<Config> mesh{box, min_level, max_level, periodic, 0.05, 2};
+    auto config = samurai::mesh_config<dim>()
+                      .min_level(min_level)
+                      .max_level(max_level)
+                      .periodic(true)
+                      .scaling_factor(2)
+                      .max_stencil_radius(1)
+                      .graduation_width(2);
 
+    auto mesh = samurai::make_MRMesh(config, box);
     auto u    = samurai::make_scalar_field<double>("u", mesh);
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
 
@@ -187,7 +188,7 @@ int main(int argc, char* argv[])
     //   Time iteration   //
     //--------------------//
 
-    double dx = mesh.cell_length(max_level);
+    double dx = mesh.min_cell_length();
     dt        = cfl * dx / velocity(0);
 
     while (t != Tf)
diff --git a/demos/FiniteVolume/heat.cpp b/demos/FiniteVolume/heat.cpp
index 0062ae1dc..4e81822a9 100644
--- a/demos/FiniteVolume/heat.cpp
+++ b/demos/FiniteVolume/heat.cpp
@@ -49,7 +49,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the heat equation", argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim>;
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
@@ -83,10 +82,6 @@ int main(int argc, char* argv[])
     double t0            = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 3;
-    std::size_t max_level = dim == 1 ? 5 : 6;
-
     // Output parameters
     fs::path path              = fs::current_path();
     std::string filename       = "heat_" + std::to_string(dim) + "D";
@@ -102,8 +97,6 @@ int main(int argc, char* argv[])
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_flag("--save-final-state-only", save_final_state_only, "Save final state only")->group("Output");
@@ -121,13 +114,14 @@ int main(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh;
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(dim == 1 ? 5 : 8).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
 
     auto u = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         u.resize();
         // Initial solution
         if (init_sol == "dirac")
@@ -175,7 +169,7 @@ int main(int argc, char* argv[])
 
     if (explicit_scheme)
     {
-        double dx = mesh.cell_length(max_level);
+        double dx = mesh.min_cell_length();
         dt        = cfl * (dx * dx) / (pow(2, dim) * diff_coeff);
     }
 
diff --git a/demos/FiniteVolume/heat_heterogeneous.cpp b/demos/FiniteVolume/heat_heterogeneous.cpp
index 8b4823da2..82008351c 100644
--- a/demos/FiniteVolume/heat_heterogeneous.cpp
+++ b/demos/FiniteVolume/heat_heterogeneous.cpp
@@ -37,7 +37,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the heat equation", argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim>;
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
@@ -59,10 +58,6 @@ int main(int argc, char* argv[])
     double t             = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 3;
-    std::size_t max_level = 6;
-
     // Output parameters
     fs::path path              = fs::current_path();
     std::string filename       = "heat_heterog_" + std::to_string(dim) + "D";
@@ -74,8 +69,6 @@ int main(int argc, char* argv[])
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_flag("--save-final-state-only", save_final_state_only, "Save final state only")->group("Output");
@@ -93,14 +86,15 @@ int main(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh;
+    auto config = samurai::mesh_config<dim>().min_level(3).max_level(6).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
 
     auto u    = samurai::make_scalar_field<double>("u", mesh);
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         u.resize();
         // Initial solution: crenel
         samurai::for_each_cell(mesh,
@@ -153,7 +147,7 @@ int main(int argc, char* argv[])
     if (explicit_scheme)
     {
         double diff_coeff = 4;
-        double dx         = mesh.cell_length(max_level);
+        double dx         = mesh.min_cell_length();
         dt                = cfl * (dx * dx) / (pow(2, dim) * diff_coeff);
     }
 
diff --git a/demos/FiniteVolume/heat_nonlinear.cpp b/demos/FiniteVolume/heat_nonlinear.cpp
index c9c259e11..b1c27b8b2 100644
--- a/demos/FiniteVolume/heat_nonlinear.cpp
+++ b/demos/FiniteVolume/heat_nonlinear.cpp
@@ -102,7 +102,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the heat equation", argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim>;
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
@@ -134,10 +133,6 @@ int main(int argc, char* argv[])
     double t             = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 4;
-    std::size_t max_level = 4;
-
     // Output parameters
     fs::path path              = fs::current_path();
     std::string filename       = "heat_nonlinear_" + std::to_string(dim) + "D";
@@ -149,8 +144,6 @@ int main(int argc, char* argv[])
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_flag("--save-final-state-only", save_final_state_only, "Save final state only")->group("Output");
@@ -168,12 +161,14 @@ int main(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh;
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(4).max_stencil_size(2).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
+    ;
     auto u = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         u    = samurai::make_scalar_field<double>("u",
                                                mesh,
                                                [&](const auto& coords)
@@ -204,7 +199,7 @@ int main(int argc, char* argv[])
     if (explicit_scheme)
     {
         double diff_coeff = 1;
-        double dx         = mesh.cell_length(max_level);
+        double dx         = mesh.min_cell_length();
         dt                = cfl * (dx * dx) / (pow(2, dim) * diff_coeff);
     }
 
diff --git a/demos/FiniteVolume/level_set.cpp b/demos/FiniteVolume/level_set.cpp
index 49d2b7c60..eaf2ea092 100644
--- a/demos/FiniteVolume/level_set.cpp
+++ b/demos/FiniteVolume/level_set.cpp
@@ -90,7 +90,7 @@ void AMR_criteria(const Field& f, Tag& tag)
     samurai::for_each_cell(mesh[mesh_id_t::cells],
                            [&](auto cell)
                            {
-                               const double dx = mesh.cell_length(max_level);
+                               const double dx = mesh.min_cell_length();
 
                                if (std::abs(f[cell]) < 1.2 * 5 * std::sqrt(2.) * dx)
                                {
@@ -253,7 +253,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example with a level set in 2d using AMR", argc, argv);
 
     constexpr std::size_t dim = 2;
-    using Config              = samurai::amr::Config<dim, 2>;
 
     // Simulation parameters
     xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
@@ -264,10 +263,7 @@ int main(int argc, char* argv[])
     std::string restart_file;
 
     // AMR parameters
-    std::size_t start_level = 8;
-    std::size_t min_level   = 4;
-    std::size_t max_level   = 8;
-    bool correction         = false;
+    bool correction = false;
 
     // Output parameters
     fs::path path        = fs::current_path();
@@ -280,9 +276,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--start-level", start_level, "Start level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--min-level", min_level, "Minimum level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--max-level", max_level, "Maximum level of AMR")->capture_default_str()->group("AMR parameters");
     app.add_option("--with-correction", correction, "Apply flux correction at the interface of two refinement levels")
         ->capture_default_str()
         ->group("AMR parameters");
@@ -292,12 +285,13 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
-    samurai::amr::Mesh<Config> mesh;
-    auto phi = samurai::make_scalar_field<double>("phi", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(8).start_level(8).max_stencil_radius(2);
+    auto mesh   = samurai::amr::make_empty_Mesh(config);
+    auto phi    = samurai::make_scalar_field<double>("phi", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, start_level, min_level, max_level};
+        mesh = samurai::amr::make_Mesh(config, box);
         init_level_set(phi);
     }
     else
@@ -305,7 +299,7 @@ int main(int argc, char* argv[])
         samurai::load(restart_file, mesh, phi);
     }
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     auto u = init_velocity(mesh);
@@ -373,9 +367,9 @@ int main(int argc, char* argv[])
             // TVD-RK2
             samurai::update_ghost(phi);
             phihat.resize();
-            phihat = phi - dt_fict * H_wrap(phi, phi_0, max_level);
+            phihat = phi - dt_fict * H_wrap(phi, phi_0, mesh.max_level());
             samurai::update_ghost(phihat);
-            phinp1 = .5 * phi_0 + .5 * (phihat - dt_fict * H_wrap(phihat, phi_0, max_level));
+            phinp1 = .5 * phi_0 + .5 * (phihat - dt_fict * H_wrap(phihat, phi_0, mesh.max_level()));
 
             std::swap(phi.array(), phinp1.array());
         }
diff --git a/demos/FiniteVolume/level_set_from_scratch.cpp b/demos/FiniteVolume/level_set_from_scratch.cpp
index 2000b1adc..e8c5865e1 100644
--- a/demos/FiniteVolume/level_set_from_scratch.cpp
+++ b/demos/FiniteVolume/level_set_from_scratch.cpp
@@ -9,6 +9,7 @@
 #include <samurai/io/hdf5.hpp>
 #include <samurai/io/restart.hpp>
 #include <samurai/mesh.hpp>
+#include <samurai/mesh_config.hpp>
 #include <samurai/mr/operators.hpp>
 #include <samurai/samurai.hpp>
 
@@ -51,16 +52,10 @@ struct fmt::formatter<SimpleID> : formatter<string_view>
     }
 };
 
-template <std::size_t dim_>
-struct AMRConfig
+template <std::size_t dim_, int prediction_stencil_radius_ = 1, std::size_t max_refinement_level_ = 20, class interval_t_ = samurai::Interval<int>>
+struct AMRConfig : samurai::mesh_config<dim_, prediction_stencil_radius_, max_refinement_level_, interval_t_>
 {
-    static constexpr std::size_t dim                  = dim_;
-    static constexpr std::size_t max_refinement_level = 20;
-    static constexpr int max_stencil_width            = 2;
-    static constexpr int ghost_width                  = 2;
-    static constexpr std::size_t prediction_order     = 1;
-    using interval_t                                  = samurai::Interval<int>;
-    using mesh_id_t                                   = SimpleID;
+    using mesh_id_t = SimpleID;
 };
 
 template <class Config>
@@ -86,35 +81,37 @@ class AMRMesh : public samurai::Mesh_base<AMRMesh<Config>, Config>
     {
     }
 
-    inline AMRMesh(const cl_type& cl, std::size_t min_level, std::size_t max_level)
-        : base_type(cl, min_level, max_level)
+    inline AMRMesh(const samurai::mesh_config<Config::dim>& cfg, const cl_type& cl)
+        : base_type(cfg, cl)
     {
     }
 
-    inline AMRMesh(const ca_type& ca, std::size_t min_level, std::size_t max_level)
-        : base_type(ca, min_level, max_level)
+    inline AMRMesh(const samurai::mesh_config<Config::dim>& cfg, const ca_type& ca)
+        : base_type(cfg, ca)
     {
     }
 
-    inline AMRMesh(const samurai::Box<double, dim>& b, std::size_t start_level, std::size_t min_level, std::size_t max_level)
-        : base_type(b, start_level, min_level, max_level)
+    inline AMRMesh(samurai::mesh_config<Config::dim>& cfg, const samurai::Box<double, dim>& b)
+        : base_type(cfg, b)
     {
     }
 
     inline void update_sub_mesh_impl()
     {
         cl_type cl;
-        for_each_interval(
-            this->cells()[mesh_id_t::cells],
-            [&](std::size_t level, const auto& interval, const auto& index_yz)
-            {
-                samurai::static_nested_loop<dim - 1, -config::ghost_width, config::ghost_width + 1>(
-                    [&](auto stencil)
-                    {
-                        auto index = xt::eval(index_yz + stencil);
-                        cl[level][index].add_interval({interval.start - config::ghost_width, interval.end + config::ghost_width});
-                    });
-            });
+        auto ghost_width = this->cfg().ghost_width();
+        for_each_interval(this->cells()[mesh_id_t::cells],
+                          [&](std::size_t level, const auto& interval, const auto& index_yz)
+                          {
+                              samurai::static_nested_loop<dim - 1>(
+                                  -ghost_width,
+                                  ghost_width + 1,
+                                  [&](auto stencil)
+                                  {
+                                      auto index = xt::eval(index_yz + stencil);
+                                      cl[level][index].add_interval({interval.start - ghost_width, interval.end + ghost_width});
+                                  });
+                          });
         this->cells()[mesh_id_t::cells_and_ghosts] = {cl, false};
     }
 };
@@ -335,7 +332,7 @@ void AMR_criteria(const Field& f, Tag& tag)
     samurai::for_each_cell(mesh[SimpleID::cells],
                            [&](auto cell)
                            {
-                               const double dx = mesh.cell_length(max_level);
+                               const double dx = mesh.min_cell_length();
 
                                if (std::abs(f[cell]) < 1.2 * 5 * std::sqrt(2.) * dx)
                                {
@@ -606,10 +603,7 @@ int main(int argc, char* argv[])
     std::string restart_file;
 
     // AMR parameters
-    std::size_t start_level = 8;
-    std::size_t min_level   = 4;
-    std::size_t max_level   = 8;
-    bool correction         = false;
+    bool correction = false;
 
     // Output parameters
     fs::path path        = fs::current_path();
@@ -622,9 +616,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--start-level", start_level, "Start level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--min-level", min_level, "Minimum level of AMR")->capture_default_str()->group("AMR parameters");
-    app.add_option("--max-level", max_level, "Maximum level of AMR")->capture_default_str()->group("AMR parameters");
     app.add_flag("--with-correction", correction, "Apply flux correction at the interface of two refinement levels")
         ->capture_default_str()
         ->group("AMR parameters");
@@ -634,13 +625,15 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(8).start_level(8).max_stencil_radius(2);
+    config.parse_args();
     AMRMesh<Config> mesh;
 
     auto phi = samurai::make_scalar_field<double>("phi", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, max_level, min_level, max_level};
+        mesh = AMRMesh<Config>(config, box);
         init_level_set(phi);
     }
     else
@@ -648,7 +641,7 @@ int main(int argc, char* argv[])
         samurai::load(restart_file, mesh, phi);
     }
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     // We initialize the level set function
@@ -715,9 +708,9 @@ int main(int argc, char* argv[])
             update_ghosts(phi, u);
             auto phihat = samurai::make_scalar_field<double>("phi", mesh);
             samurai::make_bc<samurai::Neumann<1>>(phihat, 0.);
-            phihat = phi - dt_fict * H_wrap(phi, phi_0, max_level);
+            phihat = phi - dt_fict * H_wrap(phi, phi_0, mesh.max_level());
             update_ghosts(phihat, u);
-            phinp1 = .5 * phi_0 + .5 * (phihat - dt_fict * H_wrap(phihat, phi_0, max_level));
+            phinp1 = .5 * phi_0 + .5 * (phihat - dt_fict * H_wrap(phihat, phi_0, mesh.max_level()));
             std::swap(phi.array(), phinp1.array());
         }
 
diff --git a/demos/FiniteVolume/lid_driven_cavity.cpp b/demos/FiniteVolume/lid_driven_cavity.cpp
index 004cdde63..b6b1007a8 100644
--- a/demos/FiniteVolume/lid_driven_cavity.cpp
+++ b/demos/FiniteVolume/lid_driven_cavity.cpp
@@ -96,12 +96,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Lid-driven cavity", argc, argv);
 
     constexpr std::size_t dim = 2;
-    using Config              = samurai::MRConfig<dim, 2>;
-    using Mesh                = samurai::MRMesh<Config>;
-    using mesh_id_t           = typename Mesh::mesh_id_t;
-
-    using Config2 = samurai::MRConfig<dim, 3>;
-    using Mesh2   = samurai::MRMesh<Config2>;
 
     static constexpr bool is_soa     = false;
     static constexpr bool monolithic = true;
@@ -110,11 +104,9 @@ int main(int argc, char* argv[])
     //   Parameters   //
     //----------------//
 
-    std::size_t min_level = 3;
-    std::size_t max_level = 6;
-    double Tf             = 5.;
-    double dt             = Tf / 100;
-    double cfl            = 0.95;
+    double Tf  = 5.;
+    double dt  = Tf / 100;
+    double cfl = 0.95;
 
     int ink_init = 20;
 
@@ -128,8 +120,6 @@ int main(int argc, char* argv[])
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -158,7 +148,10 @@ int main(int argc, char* argv[])
     // where v = velocity
     //       p = pressure
 
-    auto mesh = Mesh(box, min_level, max_level);
+    auto config = samurai::mesh_config<dim>().min_level(3).max_level(6).max_stencil_radius(2);
+    auto mesh   = samurai::make_MRMesh(config, box);
+
+    using mesh_id_t = typename decltype(mesh)::mesh_id_t;
 
     // Fields for the Navier-Stokes equations
     auto velocity     = samurai::make_vector_field<double, dim, is_soa>("velocity", mesh);
@@ -235,7 +228,11 @@ int main(int argc, char* argv[])
     //               d(i)/dt + conv(i) = 0,       where conv(i) = v.grad(i).
 
     // 2nd mesh
-    auto mesh2 = Mesh2(box, 1, max_level);
+    auto config2 = samurai::mesh_config<dim>().min_level(1).max_level(mesh.max_level()).max_stencil_size(6).disable_args_parse();
+    auto mesh2   = samurai::make_MRMesh(config2, box);
+
+    std::cout << config2.start_level() << " " << config2.min_level() << " " << config2.max_level() << std::endl;
+    std::cout << mesh2.cfg().start_level() << " " << mesh2.cfg().min_level() << " " << mesh2.cfg().max_level() << std::endl;
 
     // Ink data fields
     auto ink     = samurai::make_scalar_field<double>("ink", mesh2);
@@ -325,7 +322,7 @@ int main(int argc, char* argv[])
         std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t, dt);
 
         // Mesh adaptation for Navier-Stokes
-        if (min_level != max_level)
+        if (mesh.min_level() != mesh.max_level())
         {
             MRadaptation(mra_config);
             min_level_np1    = mesh[mesh_id_t::cells].min_level();
diff --git a/demos/FiniteVolume/linear_convection.cpp b/demos/FiniteVolume/linear_convection.cpp
index a7ae59070..b99e792d5 100644
--- a/demos/FiniteVolume/linear_convection.cpp
+++ b/demos/FiniteVolume/linear_convection.cpp
@@ -37,7 +37,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the linear convection equation", argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim, 3>;
     using Box                        = samurai::Box<double, dim>;
     using point_t                    = typename Box::point_t;
 
@@ -59,10 +58,6 @@ int main(int argc, char* argv[])
     double t   = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 1;
-    std::size_t max_level = dim == 1 ? 6 : 4;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "linear_convection_" + std::to_string(dim) + "D";
@@ -76,8 +71,6 @@ int main(int argc, char* argv[])
     app.add_flag("--implicit", implicit_scheme, "Implicit scheme instead of explicit")->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -93,14 +86,13 @@ int main(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    std::array<bool, dim> periodic;
-    periodic.fill(true);
-    samurai::MRMesh<Config> mesh;
-    auto u = samurai::make_scalar_field<double>("u", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(1).max_level(dim == 1 ? 6 : 4).periodic(true).max_stencil_size(6);
+    auto mesh   = samurai::make_empty_MRMesh(config);
+    auto u      = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level, periodic};
+        mesh = samurai::make_MRMesh(config, box);
         // Initial solution
         u = samurai::make_scalar_field<double>("u",
                                                mesh,
@@ -145,7 +137,7 @@ int main(int argc, char* argv[])
 
     if (dt == 0)
     {
-        double dx             = mesh.cell_length(max_level);
+        double dx             = mesh.min_cell_length();
         auto a                = xt::abs(velocity);
         double sum_velocities = xt::sum(xt::abs(velocity))();
         dt                    = cfl * dx / sum_velocities;
diff --git a/demos/FiniteVolume/linear_convection_obstacle.cpp b/demos/FiniteVolume/linear_convection_obstacle.cpp
index 00a048657..922104048 100644
--- a/demos/FiniteVolume/linear_convection_obstacle.cpp
+++ b/demos/FiniteVolume/linear_convection_obstacle.cpp
@@ -36,8 +36,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the linear convection equation", argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim, 3>;
-    using Mesh                       = samurai::MRMesh<Config>;
 
     std::cout << "------------------------- Linear convection -------------------------" << std::endl;
 
@@ -62,8 +60,6 @@ int main(int argc, char* argv[])
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -77,7 +73,8 @@ int main(int argc, char* argv[])
     samurai::DomainBuilder<dim> domain({-1., -1.}, {1., 1.});
     domain.remove({0.0, 0.0}, {0.4, 0.4});
 
-    Mesh mesh(domain, min_level, max_level);
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).max_stencil_size(6);
+    auto mesh   = samurai::make_MRMesh(config, domain);
 
     // Initial solution
     auto u = samurai::make_scalar_field<double>("u",
@@ -117,7 +114,7 @@ int main(int argc, char* argv[])
 
     if (dt == 0)
     {
-        double dx             = mesh.cell_length(max_level);
+        double dx             = mesh.min_cell_length();
         auto a                = xt::abs(constant_velocity);
         double sum_velocities = xt::sum(xt::abs(constant_velocity))();
         dt                    = cfl * dx / sum_velocities;
diff --git a/demos/FiniteVolume/manual_block_matrix_assembly.cpp b/demos/FiniteVolume/manual_block_matrix_assembly.cpp
index 03b86f243..5c165d18e 100644
--- a/demos/FiniteVolume/manual_block_matrix_assembly.cpp
+++ b/demos/FiniteVolume/manual_block_matrix_assembly.cpp
@@ -188,17 +188,14 @@ int main(int argc, char* argv[])
     samurai::initialize(argc, argv);
 
     static constexpr std::size_t dim = 2;
-    using Config                     = samurai::MRConfig<dim>;
     using Box                        = samurai::Box<double, dim>;
 
     std::cout << "------------------------- Begin -------------------------" << std::endl;
 
-    std::size_t min_level = 3;
-    std::size_t max_level = 3;
-
     Box box({0, 0}, {1, 1});
-    samurai::MRMesh<Config> mesh_e{box, min_level, max_level};
-    samurai::MRMesh<Config> mesh_s{box, min_level, max_level};
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(3).max_level(3);
+    auto mesh_e   = samurai::make_MRMesh(mesh_cfg, box);
+    auto mesh_s   = samurai::make_MRMesh(mesh_cfg, box);
 
     //-------------------------------//
     // Fields and auxiliary unknowns //
diff --git a/demos/FiniteVolume/nagumo.cpp b/demos/FiniteVolume/nagumo.cpp
index 5f051d90e..7ad13ada2 100644
--- a/demos/FiniteVolume/nagumo.cpp
+++ b/demos/FiniteVolume/nagumo.cpp
@@ -37,7 +37,6 @@ int main(int argc, char* argv[])
 
     static constexpr std::size_t dim    = 1;
     static constexpr std::size_t n_comp = 1;
-    using Config                        = samurai::MRConfig<dim>;
     using Box                           = samurai::Box<double, dim>;
     using point_t                       = typename Box::point_t;
 
@@ -70,10 +69,6 @@ int main(int argc, char* argv[])
     double t   = 0.;
     std::string restart_file;
 
-    // Multiresolution parameters
-    std::size_t min_level = 4;
-    std::size_t max_level = 8;
-
     // Output parameters
     fs::path path              = fs::current_path();
     std::string filename       = "nagumo";
@@ -90,8 +85,6 @@ int main(int argc, char* argv[])
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
     app.add_flag("--explicit-reaction", explicit_reaction, "Explicit the reaction term")->capture_default_str()->group("Simulation parameters");
     app.add_flag("--explicit-diffusion", explicit_diffusion, "Explicit the diffusion term")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_flag("--save-final-state-only", save_final_state_only, "Save final state only")->group("Output");
@@ -109,9 +102,9 @@ int main(int argc, char* argv[])
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh{box, min_level, max_level};
-
-    auto u = samurai::make_vector_field<double, n_comp>("u", mesh);
+    auto config = samurai::mesh_config<dim>().min_level(4).max_level(8).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_MRMesh(config, box);
+    auto u      = samurai::make_vector_field<double, n_comp>("u", mesh);
 
     double z0 = left_box / 5;    // wave initial position
     double c  = sqrt(k * D / 2); // wave velocity
@@ -129,7 +122,6 @@ int main(int argc, char* argv[])
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
         u.resize();
         // Initial solution
         samurai::for_each_cell(mesh,
@@ -177,7 +169,7 @@ int main(int argc, char* argv[])
         dt = Tf / 100;
         if (explicit_diffusion)
         {
-            double dx = mesh.cell_length(max_level);
+            double dx = mesh.min_cell_length();
             dt        = cfl * (dx * dx) / (pow(2, dim) * D);
         }
     }
diff --git a/demos/FiniteVolume/scalar_burgers_2d.cpp b/demos/FiniteVolume/scalar_burgers_2d.cpp
index 7e1c1a324..c9af54bb4 100644
--- a/demos/FiniteVolume/scalar_burgers_2d.cpp
+++ b/demos/FiniteVolume/scalar_burgers_2d.cpp
@@ -182,7 +182,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Finite volume example for the scalar Burgers equation in 2d using multiresolution", argc, argv);
 
     constexpr size_t dim = 2;
-    using Config         = samurai::MRConfig<dim>;
 
     // Simulation parameters
     xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
@@ -212,8 +211,6 @@ int main(int argc, char* argv[])
     app.add_option("--Ti", t, "Initial time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--restart-file", restart_file, "Restart file")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--with-correction", correction, "Apply flux correction at the interface of two refinement levels")
         ->capture_default_str()
         ->group("Multiresolution");
@@ -224,13 +221,14 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
-    samurai::MRMesh<Config> mesh;
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).disable_minimal_ghost_width();
+    auto mesh   = samurai::make_empty_MRMesh(config);
 
     auto u = samurai::make_scalar_field<double>("u", mesh);
 
     if (restart_file.empty())
     {
-        mesh = {box, min_level, max_level};
+        mesh = samurai::make_MRMesh(config, box);
         init(u);
     }
     else
@@ -240,7 +238,7 @@ int main(int argc, char* argv[])
 
     samurai::make_bc<samurai::Dirichlet<1>>(u, 0.);
 
-    double dt            = cfl * mesh.cell_length(max_level);
+    double dt            = cfl * mesh.min_cell_length();
     const double dt_save = Tf / static_cast<double>(nfiles);
 
     auto unp1 = samurai::make_scalar_field<double>("unp1", mesh);
diff --git a/demos/FiniteVolume/stokes_2d.cpp b/demos/FiniteVolume/stokes_2d.cpp
index 463786eb7..ca96994d9 100644
--- a/demos/FiniteVolume/stokes_2d.cpp
+++ b/demos/FiniteVolume/stokes_2d.cpp
@@ -153,9 +153,6 @@ int main(int argc, char* argv[])
     auto& app = samurai::initialize("Stokes problem", argc, argv);
 
     constexpr std::size_t dim        = 2;
-    using Config                     = samurai::MRConfig<dim, 2>;
-    using Mesh                       = samurai::MRMesh<Config>;
-    using mesh_id_t                  = typename Mesh::mesh_id_t;
     static constexpr bool is_soa     = false;
     static constexpr bool monolithic = true;
 
@@ -165,10 +162,8 @@ int main(int argc, char* argv[])
 
     std::string test_case = "ns";
 
-    std::size_t min_level = 5;
-    std::size_t max_level = 5;
-    double Tf             = 1.;
-    double dt             = Tf / 100;
+    double Tf = 1.;
+    double dt = Tf / 100;
 
     std::size_t nfiles = 50;
 
@@ -180,8 +175,6 @@ int main(int argc, char* argv[])
         ->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -196,8 +189,11 @@ int main(int argc, char* argv[])
     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 
-    auto box  = samurai::Box<double, dim>({0, 0}, {1, 1});
-    auto mesh = Mesh(box, min_level, max_level);
+    auto box    = samurai::Box<double, dim>({0, 0}, {1, 1});
+    auto config = samurai::mesh_config<dim>().min_level(5).max_level(5).max_stencil_size(2);
+    auto mesh   = samurai::make_MRMesh(config, box);
+
+    using mesh_id_t = typename decltype(mesh)::mesh_id_t;
 
     //--------------------//
     // Stationary problem //
@@ -467,7 +463,7 @@ int main(int argc, char* argv[])
             std::cout << fmt::format("iteration {}: t = {:.2f}, dt = {}", nt++, t_np1, dt);
 
             // Mesh adaptation
-            if (min_level != max_level)
+            if (mesh.min_level() != mesh.max_level())
             {
                 MRadaptation(mra_config);
                 min_level_np1    = mesh[mesh_id_t::cells].min_level();
@@ -548,7 +544,7 @@ int main(int argc, char* argv[])
                 samurai::save(path, fmt::format("{}_ite_{}", filename, nsave), velocity.mesh(), velocity, div_velocity);
             }
 
-            if (min_level != max_level)
+            if (mesh.min_level() != mesh.max_level())
             {
                 // Reconstruction on the finest level
                 auto velocity_recons = samurai::reconstruction(velocity);
diff --git a/demos/Weno/VF_level_set_houc5_amr.cpp b/demos/Weno/VF_level_set_houc5_amr.cpp
index 541859416..83915234d 100644
--- a/demos/Weno/VF_level_set_houc5_amr.cpp
+++ b/demos/Weno/VF_level_set_houc5_amr.cpp
@@ -261,19 +261,15 @@ int main()
 {
     samurai::initialize();
 
-    constexpr std::size_t dim         = 2;
-    constexpr std::size_t ghost_width = 3;
-
-    std::size_t start_level = 8;
-    std::size_t min_level   = 2;
-    std::size_t max_level   = 8;
+    constexpr std::size_t dim = 2;
 
     samurai::Box<double, dim> box{
         {0, 0},
         {1, 1}
     };
-    using Config = samurai::amr::Config<dim, ghost_width>;
-    samurai::amr::Mesh<Config> mesh(box, start_level, min_level, max_level);
+    using Config  = samurai::amr::Config<dim>;
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(8).start_level(8).max_stencil_radius(3);
+    samurai::amr::Mesh<Config> mesh(mesh_cfg, box);
 
     auto field = init_field(mesh);
 
@@ -282,7 +278,7 @@ int main()
         update_bc_for_level_(field, level);
     };
 
-    double dt = 0.5 / (1 << max_level);
+    double dt = 0.5 * mesh.min_cell_length();
     double Tf = 2.; // Final time
 
     std::size_t max_ite = Tf / dt;
diff --git a/demos/Weno/VF_level_set_os5_amr.cpp b/demos/Weno/VF_level_set_os5_amr.cpp
index 612435186..fc130c05e 100644
--- a/demos/Weno/VF_level_set_os5_amr.cpp
+++ b/demos/Weno/VF_level_set_os5_amr.cpp
@@ -278,18 +278,15 @@ int main()
 {
     samurai::initialize();
 
-    constexpr std::size_t dim         = 2;
-    constexpr std::size_t ghost_width = 3;
-    std::size_t start_level           = 6;
-    std::size_t min_level             = 2;
-    std::size_t max_level             = 8;
+    constexpr std::size_t dim = 2;
 
     samurai::Box<double, dim> box{
         {0, 0},
         {1, 1}
     };
-    using Config = samurai::amr::Config<dim, ghost_width>;
-    samurai::amr::Mesh<Config> mesh(box, start_level, min_level, max_level);
+    using Config  = samurai::amr::Config<dim>;
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(8).start_level(6).max_stencil_radius(3);
+    samurai::amr::Mesh<Config> mesh(mesh_cfg, box);
 
     auto field = init_field(mesh);
 
diff --git a/demos/Weno/heat_weno5_amr.cpp b/demos/Weno/heat_weno5_amr.cpp
index e3bffc6a9..dee5aa7e2 100644
--- a/demos/Weno/heat_weno5_amr.cpp
+++ b/demos/Weno/heat_weno5_amr.cpp
@@ -720,18 +720,15 @@ int main()
 {
     samurai::initialize();
 
-    constexpr std::size_t dim         = 2;
-    constexpr std::size_t ghost_width = 3;
-    std::size_t start_level           = 7;
-    std::size_t min_level             = 2;
-    std::size_t max_level             = 10;
+    constexpr std::size_t dim = 2;
 
     samurai::Box<double, dim> box{
         {0, 0, 0},
         {1, 1, 1}
     };
-    using Config = samurai::amr::Config<dim, ghost_width>;
-    samurai::amr::Mesh<Config> mesh(box, start_level, min_level, max_level);
+    using Config  = samurai::amr::Config<dim>;
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(10).start_level(7).max_stencil_radius(3);
+    samurai::amr::Mesh<Config> mesh(mesh_cfg, box);
     using mesh_id_t = typename Config::mesh_id_t;
 
     auto field = init_field::call(samurai::Dim<dim>{}, mesh);
diff --git a/demos/Weno/weno5_amr.cpp b/demos/Weno/weno5_amr.cpp
index 5492c6f2a..10eeed861 100644
--- a/demos/Weno/weno5_amr.cpp
+++ b/demos/Weno/weno5_amr.cpp
@@ -437,18 +437,15 @@ int main()
 {
     samurai::initialize();
 
-    constexpr std::size_t dim         = 2;
-    constexpr std::size_t ghost_width = 3;
-    std::size_t start_level           = 7;
-    std::size_t min_level             = 2;
-    std::size_t max_level             = 10;
+    constexpr std::size_t dim = 2;
 
     samurai::Box<double, dim> box{
         {0, 0, 0},
         {1, 1, 1}
     };
-    using Config = samurai::amr::Config<dim, ghost_width>;
-    samurai::amr::Mesh<Config> mesh(box, start_level, min_level, max_level);
+    using Config  = samurai::amr::Config<dim>;
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(10).start_level(7).max_stencil_radius(3);
+    samurai::amr::Mesh<Config> mesh(mesh_cfg, box);
     using mesh_id_t = typename Config::mesh_id_t;
 
     auto field = init_field::call(samurai::Dim<dim>{}, mesh);
diff --git a/demos/from_obj/main.cpp b/demos/from_obj/main.cpp
index 8821a4b69..1b7ada890 100644
--- a/demos/from_obj/main.cpp
+++ b/demos/from_obj/main.cpp
@@ -30,7 +30,7 @@ void save_mesh(const fs::path& path, const std::string& filename, const Mesh& me
 int main(int argc, char** argv)
 {
     constexpr std::size_t dim = 3;
-    std::size_t start_level   = 1;
+    std::size_t initial_level = 1;
     std::size_t max_level     = 8;
     bool keep_inside          = false;
     bool keep_outside         = false;
@@ -41,7 +41,7 @@ int main(int argc, char** argv)
 
     CLI::App app{"Create an adapted mesh from an OBJ file"};
     app.add_option("--input", input_file, "input File")->required()->check(CLI::ExistingFile);
-    app.add_option("--start-level", start_level, "Start level of the output adaptive mesh")->capture_default_str();
+    app.add_option("--init-level", initial_level, "Start level of the output adaptive mesh")->capture_default_str();
     app.add_flag("--keep-inside", keep_inside, "Keep the cells inside the object")->capture_default_str();
     app.add_flag("--keep-outside", keep_outside, "Keep the cells outside the object")->capture_default_str();
     app.add_option("--max-level", max_level, "Maximum level of the output adaptive mesh")->capture_default_str();
@@ -50,7 +50,7 @@ int main(int argc, char** argv)
 
     std::string output_file = fs::path(input_file).stem();
 
-    auto mesh = samurai::from_geometry<dim>(input_file, start_level, max_level, keep_outside, keep_inside);
+    auto mesh = samurai::from_geometry<dim>(input_file, initial_level, max_level, keep_outside, keep_inside);
     make_graduation(mesh);
 
     save_mesh(path, fmt::format("mesh_{}", output_file), mesh);
diff --git a/demos/highorder/main.cpp b/demos/highorder/main.cpp
index 9c68c95a0..56519636b 100644
--- a/demos/highorder/main.cpp
+++ b/demos/highorder/main.cpp
@@ -112,19 +112,15 @@ int main(int argc, char* argv[])
 {
     auto& app = samurai::initialize("Finite volume example for the advection equation in 2d using multiresolution", argc, argv);
 
-    constexpr std::size_t dim              = 2;
-    constexpr std::size_t stencil_width    = 2;
-    constexpr std::size_t graduation_width = 4;
-    constexpr std::size_t prediction_order = 1;
-    using Config                           = samurai::MRConfig<dim, stencil_width, graduation_width, prediction_order>;
+    constexpr std::size_t dim = 2;
+
+    constexpr std::size_t prediction_stencil_radius = 1;
 
     // Simulation parameters
     xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
     xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1.};
 
     // Multiresolution parameters
-    std::size_t min_level  = 4;
-    std::size_t max_level  = 4;
     std::size_t refinement = 5;
     bool correction        = false;
 
@@ -135,8 +131,6 @@ int main(int argc, char* argv[])
 
     app.add_option("--min-corner", min_corner, "The min corner of the box")->capture_default_str()->group("Simulation parameters");
     app.add_option("--max-corner", min_corner, "The max corner of the box")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--refinement", refinement, "Number of refinement")->capture_default_str()->group("Multiresolution");
     app.add_option("--with-correction", correction, "Apply flux correction at the interface of two refinement levels")
         ->capture_default_str()
@@ -147,10 +141,19 @@ int main(int argc, char* argv[])
     SAMURAI_PARSE(argc, argv);
 
     samurai::Box<double, dim> box(min_corner, max_corner);
-    using mesh_t    = samurai::MRMesh<Config>;
+
+    // clang-format off
+    auto config = samurai::mesh_config<dim, prediction_stencil_radius>()
+    .min_level(4)
+    .max_level(4)
+    .graduation_width(4)
+    .max_stencil_size(4);
+    // clang-format on
+    auto init_mesh = samurai::make_MRMesh(config, box);
+
+    using mesh_t    = decltype(init_mesh);
     using mesh_id_t = typename mesh_t::mesh_id_t;
     using cl_type   = typename mesh_t::cl_type;
-    mesh_t init_mesh{box, min_level, max_level};
 
     auto adapt_field = samurai::make_scalar_field<double>("adapt_field",
                                                           init_mesh,
@@ -196,7 +199,9 @@ int main(int argc, char* argv[])
                                                });
                                        });
             // mesh = {cl, mesh.min_level() + 1, mesh.max_level() + 1};
-            mesh = {cl, min_level + i_ref + 1, max_level + i_ref + 1};
+            auto mesh_cfg = mesh.cfg();
+            mesh_cfg.min_level(init_mesh.min_level() + i_ref + 1).max_level(init_mesh.max_level() + i_ref + 1);
+            mesh = samurai::make_MRMesh(mesh_cfg, cl);
         }
         // std::cout << mesh << std::endl;
         // samurai::save("refine_mesh", mesh);
@@ -228,7 +233,7 @@ int main(int argc, char* argv[])
         // set.apply_op(samurai::projection(f));
         // auto f_recons = samurai::make_scalar_field<double>("f_recons", mesh);
         // auto error_f = samurai::make_scalar_field<double>("error", mesh);
-        // set.apply_op(samurai::prediction<prediction_order, true>(f_recons, f));
+        // set.apply_op(samurai::prediction<prediction_stencil_radius, true>(f_recons, f));
         // samurai::for_each_interval(mesh[mesh_id_t::cells], [&](std::size_t level,
         // const auto& i, const auto& index)
         // {
@@ -310,8 +315,8 @@ int main(int argc, char* argv[])
                                });
 
         samurai::save(fmt::format("error_ref_{}", ite), mesh, error_field);
-        samurai::save(fmt::format("solution_{}_{}_ref_{}", min_level, max_level, ite), mesh, u);
-        samurai::save(fmt::format("solution_recons_{}_{}_ref_{}", min_level, max_level, ite), u_recons.mesh(), u_recons);
+        samurai::save(fmt::format("solution_{}_{}_ref_{}", mesh.min_level(), mesh.max_level(), ite), mesh, u);
+        samurai::save(fmt::format("solution_recons_{}_{}_ref_{}", mesh.min_level(), mesh.max_level(), ite), u_recons.mesh(), u_recons);
         h_coarse            = h;
         error_coarse        = error;
         error_recons_coarse = error_recons;
diff --git a/demos/multigrid/main.cpp b/demos/multigrid/main.cpp
index c7eb04e8d..947f92d6a 100644
--- a/demos/multigrid/main.cpp
+++ b/demos/multigrid/main.cpp
@@ -66,21 +66,21 @@ static char help[] = "Solution of the Poisson problem in the domain [0,1]^d.\n"
                      "-log_view -pc_mg_log           Monitors the multigrid performance\n"
                      "\n";
 
-template <class Mesh>
-Mesh create_uniform_mesh(std::size_t level)
+template <std::size_t dim>
+auto create_uniform_mesh(std::size_t level)
 {
-    using Box = samurai::Box<double, Mesh::dim>;
+    using Box = samurai::Box<double, dim>;
 
     Box box;
-    if constexpr (Mesh::dim == 1)
+    if constexpr (dim == 1)
     {
         box = Box({0}, {1});
     }
-    else if constexpr (Mesh::dim == 2)
+    else if constexpr (dim == 2)
     {
         box = Box({0, 0}, {1, 1});
     }
-    else if constexpr (Mesh::dim == 3)
+    else if constexpr (dim == 3)
     {
         box = Box({0, 0, 0}, {1, 1, 1});
     }
@@ -89,44 +89,44 @@ Mesh create_uniform_mesh(std::size_t level)
     min_level   = level;
     max_level   = level;
 
-    return Mesh(box, start_level, min_level, max_level); // amr::Mesh
-    // return Mesh(box, /*start_level,*/ min_level, max_level); // MRMesh
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).start_level(start_level).disable_minimal_ghost_width();
+    return samurai::amr::make_Mesh(config, box); // amr::Mesh
+    // return samurai::make_MRMesh(box, /*start_level,*/ min_level, max_level); // MRMesh
 }
 
-template <class Mesh>
-[[maybe_unused]] Mesh create_refined_mesh(std::size_t level)
+template <std::size_t dim>
+[[maybe_unused]] auto create_refined_mesh(std::size_t level)
 {
-    using cl_type = typename Mesh::cl_type;
-
     std::size_t min_level, max_level;
     min_level = level - 1;
     max_level = level;
 
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).disable_minimal_ghost_width();
+
+    using cl_type = typename decltype(config)::cl_type;
+
     int i = 1 << min_level;
 
     cl_type cl;
-    if constexpr (Mesh::dim == 1)
+    if constexpr (dim == 1)
     {
         cl[min_level][{}].add_interval({0, i / 2});
         cl[max_level][{}].add_interval({i, 2 * i});
     }
-    static_assert(Mesh::dim == 1, "create_refined_mesh() not implemented for this dimension");
+    static_assert(dim == 1, "create_refined_mesh() not implemented for this dimension");
 
-    return Mesh(cl, min_level, max_level); // amr::Mesh
-    // return Mesh(box, /*start_level,*/ min_level, max_level); // MRMesh
+    return samurai::amr::make_Mesh(config, cl); // amr::Mesh
+    // return samurai::make_MRMesh(box, /*start_level,*/ min_level, max_level); // MRMesh
 }
 
 int main(int argc, char* argv[])
 {
     samurai::initialize(argc, argv);
 
-    constexpr std::size_t dim = 2;
-    using Config              = samurai::amr::Config<dim>;
-    using Mesh                = samurai::amr::Mesh<Config>;
-    // using Config = samurai::MRConfig<dim>;
-    // using Mesh = samurai::MRMesh<Config>;
+    constexpr std::size_t dim     = 2;
     constexpr unsigned int n_comp = 1;
     constexpr bool is_soa         = true;
+    using Mesh                    = decltype(create_uniform_mesh<dim>(1));
     using Field                   = samurai::VectorField<Mesh, double, n_comp, is_soa>;
 
     PetscMPIInt size;
@@ -197,7 +197,7 @@ int main(int argc, char* argv[])
     // Mesh creation //
     //---------------//
 
-    Mesh mesh = create_uniform_mesh<Mesh>(static_cast<std::size_t>(level));
+    auto mesh = create_uniform_mesh<dim>(static_cast<std::size_t>(level));
     // print_mesh(mesh);
 
     if (save_mesh)
diff --git a/demos/multigrid/samurai_new/multigrid/LevelCtx.hpp b/demos/multigrid/samurai_new/multigrid/LevelCtx.hpp
index 230d82ef1..8358d07ad 100644
--- a/demos/multigrid/samurai_new/multigrid/LevelCtx.hpp
+++ b/demos/multigrid/samurai_new/multigrid/LevelCtx.hpp
@@ -33,24 +33,24 @@ namespace samurai_new
         LevelContext* finer   = nullptr;
         LevelContext* coarser = nullptr;
         TransferOperators transfer_ops;
-        int prediction_order;
+        int prediction_stencil_radius;
 
         LevelContext(Dsctzr& d, Mesh& m, TransferOperators to, int pred_order)
             : _mesh(m)
             , _discretizer(d)
         {
-            level            = 0;
-            transfer_ops     = to;
-            prediction_order = pred_order;
+            level                     = 0;
+            transfer_ops              = to;
+            prediction_stencil_radius = pred_order;
         }
 
         LevelContext(LevelContext& fine_ctx)
             : _mesh(samurai_new::coarsen(fine_ctx.mesh()))
             , _discretizer(Dsctzr::create_coarse(fine_ctx.assembly(), _mesh))
         {
-            level            = fine_ctx.level + 1;
-            transfer_ops     = fine_ctx.transfer_ops;
-            prediction_order = fine_ctx.prediction_order;
+            level                     = fine_ctx.level + 1;
+            transfer_ops              = fine_ctx.transfer_ops;
+            prediction_stencil_radius = fine_ctx.prediction_stencil_radius;
 
             this->finer      = &fine_ctx;
             fine_ctx.coarser = this;
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp b/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
index 1e2603914..76043a2c2 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/GeometricMultigrid.hpp
@@ -57,8 +57,8 @@ namespace samurai_new
                 PetscOptionsGetInt(NULL, NULL, "--samg_transfer_ops", &transfer_ops_arg, NULL);
                 samurai_new::TransferOperators transfer_ops = static_cast<samurai_new::TransferOperators>(transfer_ops_arg);
 
-                PetscInt prediction_order = 0;
-                PetscOptionsGetInt(NULL, NULL, "--samg_pred_order", &prediction_order, NULL);
+                PetscInt prediction_stencil_radius = 0;
+                PetscOptionsGetInt(NULL, NULL, "--samg_pred_order", &prediction_stencil_radius, NULL);
 
                 PetscBool smoother_is_set = PETSC_FALSE;
                 Smoothers smoother        = SymGaussSeidel;
@@ -121,9 +121,9 @@ namespace samurai_new
                 }
                 std::cout << std::endl;
 
-                std::cout << "    prediction order  : " << prediction_order << std::endl;
+                std::cout << "    prediction order  : " << prediction_stencil_radius << std::endl;
 
-                _samuraiDM = new SamuraiDM<Dsctzr>(PETSC_COMM_SELF, *_discretizer, *_mesh, transfer_ops, prediction_order);
+                _samuraiDM = new SamuraiDM<Dsctzr>(PETSC_COMM_SELF, *_discretizer, *_mesh, transfer_ops, prediction_stencil_radius);
                 KSPSetDM(ksp, _samuraiDM->PetscDM());
 
                 // Default outer solver: CG
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp b/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
index 3c56c1458..8571731bb 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/SamuraiDM.hpp
@@ -21,8 +21,8 @@ namespace samurai_new
             using Mesh  = typename Dsctzr::Mesh;
             using Field = typename Dsctzr::field_t;
 
-            SamuraiDM(MPI_Comm comm, Dsctzr& assembly, Mesh& mesh, TransferOperators to, int prediction_order)
-                : _ctx(assembly, mesh, to, prediction_order)
+            SamuraiDM(MPI_Comm comm, Dsctzr& assembly, Mesh& mesh, TransferOperators to, int prediction_stencil_radius)
+                : _ctx(assembly, mesh, to, prediction_stencil_radius)
             {
                 DMShellCreate(comm, &_dm);
                 DefineShellFunctions(_dm, _ctx);
@@ -155,7 +155,7 @@ namespace samurai_new
                 {
                     Field coarse_field("coarse_field", coarse_ctx->mesh());
                     samurai::petsc::copy(x, coarse_field);
-                    Field fine_field = multigrid::prolong(coarse_field, fine_ctx->mesh(), coarse_ctx->prediction_order);
+                    Field fine_field = multigrid::prolong(coarse_field, fine_ctx->mesh(), coarse_ctx->prediction_stencil_radius);
                     samurai::petsc::copy(fine_field, y);
 
                     // std::cout << "prolongated vector (marche):" << std::endl;
@@ -168,7 +168,7 @@ namespace samurai_new
                     double* yarray;
                     VecGetArrayRead(x, &xarray);
                     VecGetArray(y, &yarray);
-                    multigrid::prolong(coarse_ctx->mesh(), fine_ctx->mesh(), xarray, yarray, coarse_ctx->prediction_order);
+                    multigrid::prolong(coarse_ctx->mesh(), fine_ctx->mesh(), xarray, yarray, coarse_ctx->prediction_stencil_radius);
                     VecRestoreArrayRead(x, &xarray);
                     VecRestoreArray(y, &yarray);
 
@@ -227,7 +227,7 @@ namespace samurai_new
                 MatSetSizes(*P, nf, nc, nf, nc);
                 MatSetFromOptions(*P);
                 MatSeqAIJSetPreallocation(*P, stencil_size, NULL);
-                multigrid::set_prolong_matrix(coarse_ctx->mesh(), fine_ctx->mesh(), *P, coarse_ctx->prediction_order);
+                multigrid::set_prolong_matrix(coarse_ctx->mesh(), fine_ctx->mesh(), *P, coarse_ctx->prediction_stencil_radius);
                 MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
                 MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);
 
diff --git a/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp b/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
index b0d6853a2..7baa18637 100644
--- a/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
+++ b/demos/multigrid/samurai_new/multigrid/petsc/intergrid_operators.hpp
@@ -11,7 +11,7 @@ namespace samurai_new
         {
 
             template <class Mesh>
-            void prolong(const Mesh& coarse_mesh, const Mesh& fine_mesh, const double* carray, double* farray, int prediction_order)
+            void prolong(const Mesh& coarse_mesh, const Mesh& fine_mesh, const double* carray, double* farray, int prediction_stencil_radius)
             {
                 using mesh_id_t                  = typename Mesh::mesh_id_t;
                 static constexpr std::size_t dim = Mesh::dim;
@@ -46,7 +46,7 @@ namespace samurai_new
                             {
                                 auto i_f = fine_mesh.get_index(level, (2 * i).start);
                                 auto i_c = coarse_mesh.get_index(level - 1, i.start);
-                                if (prediction_order == 0)
+                                if (prediction_stencil_radius == 0)
                                 {
                                     for (std::size_t ii = 0; ii < i.size(); ++ii)
                                     {
@@ -55,7 +55,7 @@ namespace samurai_new
                                         farray[i_f + 2 * ii + 1] = carray[i_c + ii];
                                     }
                                 }
-                                else if (prediction_order == 1)
+                                else if (prediction_stencil_radius == 1)
                                 {
                                     for (std::size_t ii = 0; ii < i.size(); ++ii)
                                     {
@@ -74,7 +74,7 @@ namespace samurai_new
                                 auto i_f_2j   = fine_mesh.get_index(level, (2 * i).start, 2 * j);
                                 auto i_f_2jp1 = fine_mesh.get_index(level, (2 * i).start, 2 * j + 1);
 
-                                if (prediction_order == 0)
+                                if (prediction_stencil_radius == 0)
                                 {
                                     for (std::size_t ii = 0; ii < i.size(); ++ii)
                                     {
@@ -85,7 +85,7 @@ namespace samurai_new
                                         farray[i_f_2jp1 + 2 * ii + 1] = carray[i_c_j + ii];
                                     }
                                 }
-                                else if (prediction_order == 1)
+                                else if (prediction_stencil_radius == 1)
                                 {
                                     auto i_c_jm1 = coarse_mesh.get_index(level - 1, i.start, j - 1);
                                     auto i_c_jp1 = coarse_mesh.get_index(level - 1, i.start, j + 1);
@@ -198,7 +198,7 @@ namespace samurai_new
             }
 
             template <class Mesh>
-            void set_prolong_matrix(const Mesh& coarse_mesh, const Mesh& fine_mesh, Mat& P, int prediction_order)
+            void set_prolong_matrix(const Mesh& coarse_mesh, const Mesh& fine_mesh, Mat& P, int prediction_stencil_radius)
             {
                 using mesh_id_t                  = typename Mesh::mesh_id_t;
                 static constexpr std::size_t dim = Mesh::dim;
@@ -235,7 +235,7 @@ namespace samurai_new
                             {
                                 auto i_f = static_cast<int>(fine_mesh.get_index(level, (2 * i).start));
                                 auto i_c = static_cast<int>(coarse_mesh.get_index(level - 1, i.start));
-                                if (prediction_order == 0)
+                                if (prediction_stencil_radius == 0)
                                 {
                                     for (int ii = 0; ii < static_cast<int>(i.size()); ++ii)
                                     {
@@ -244,7 +244,7 @@ namespace samurai_new
                                         MatSetValue(P, i_f + 2 * ii + 1, i_c + ii, 1, INSERT_VALUES);
                                     }
                                 }
-                                else if (prediction_order == 1)
+                                else if (prediction_stencil_radius == 1)
                                 {
                                     for (int ii = 0; ii < static_cast<int>(i.size()); ++ii)
                                     {
@@ -273,7 +273,7 @@ namespace samurai_new
                                 auto i_f_2j   = static_cast<int>(fine_mesh.get_index(level, (2 * i).start, 2 * j));
                                 auto i_f_2jp1 = static_cast<int>(fine_mesh.get_index(level, (2 * i).start, 2 * j + 1));
 
-                                if (prediction_order == 0)
+                                if (prediction_stencil_radius == 0)
                                 {
                                     for (int ii = 0; ii < static_cast<int>(i.size()); ++ii)
                                     {
@@ -291,7 +291,7 @@ namespace samurai_new
                                         MatSetValue(P, fine_top_right, coarse, 1, INSERT_VALUES);
                                     }
                                 }
-                                else if (prediction_order == 1)
+                                else if (prediction_stencil_radius == 1)
                                 {
                                     auto i_c_jm1 = static_cast<int>(coarse_mesh.get_index(level - 1, i.start, j - 1));
                                     auto i_c_jp1 = static_cast<int>(coarse_mesh.get_index(level - 1, i.start, j + 1));
@@ -487,7 +487,7 @@ namespace samurai_new
             }
 
             template <class Field>
-            Field prolong(const Field& coarse_field, typename Field::mesh_t& fine_mesh, int prediction_order)
+            Field prolong(const Field& coarse_field, typename Field::mesh_t& fine_mesh, int prediction_stencil_radius)
             {
                 using mesh_id_t                  = typename Field::mesh_t::mesh_id_t;
                 static constexpr std::size_t dim = Field::mesh_t::dim;
@@ -512,11 +512,11 @@ namespace samurai_new
 
                     // Coarse cells to fine cells: prediction
                     auto others = samurai::intersection(cm[level - 1], fm[level]).on(level - 1);
-                    if (prediction_order == 0)
+                    if (prediction_stencil_radius == 0)
                     {
                         others.apply_op(samurai::prediction<0, true>(fine_field, coarse_field));
                     }
-                    else if (prediction_order == 1)
+                    else if (prediction_stencil_radius == 1)
                     {
                         others.apply_op(samurai::prediction<1, true>(fine_field, coarse_field));
                     }
diff --git a/demos/p4est/simple_2d.cpp b/demos/p4est/simple_2d.cpp
index e7fea6a03..c016a5236 100644
--- a/demos/p4est/simple_2d.cpp
+++ b/demos/p4est/simple_2d.cpp
@@ -187,7 +187,7 @@ void refine_2(mesh_t& mesh, std::size_t max_level)
                                        }
                                    });
 
-        mesh = {cl, mesh.min_level(), mesh.max_level()};
+        mesh = samurai::make_MRMesh(mesh.cfg(), cl);
     }
 }
 
@@ -201,17 +201,15 @@ int main(int argc, char* argv[])
 
     constexpr size_t dim = 2;
 
-    // Adaptation parameters
-    std::size_t max_level = 9;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "simple_2d";
 
-    app.add_option("--max-level", max_level, "Maximum level of the adaptation")->capture_default_str()->group("Adaptation parameters");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     SAMURAI_PARSE(argc, argv);
+    auto config = samurai::mesh_config<dim>().min_level(1).max_level(9);
+    config.parse_args();
 
     if (!fs::exists(path))
     {
@@ -232,19 +230,17 @@ int main(int argc, char* argv[])
     std::cout << "nb_cells: " << mesh_1.nb_cells() << "\n";
 
     tic();
-    refine_1(mesh_1, max_level);
+    refine_1(mesh_1, config.max_level());
     auto duration = toc();
     std::cout << "Version 1: " << duration << "s" << std::endl;
     toc();
     std::cout << "nb_cells: " << mesh_1.nb_cells() << "\n";
 
-    // samurai::CellArray<dim> mesh_2(cl);
-    using Config    = samurai::MRConfig<dim>;
-    using mesh_id_t = typename samurai::MRMesh<Config>::mesh_id_t;
-    samurai::MRMesh<Config> mesh_2(cl, 1, max_level);
+    auto mesh_2     = samurai::make_MRMesh(config, cl);
+    using mesh_id_t = typename decltype(mesh_2)::mesh_id_t;
 
     tic();
-    refine_2(mesh_2, max_level);
+    refine_2(mesh_2, config.max_level());
     duration = toc();
     std::cout << "Version 2: " << duration << "s" << std::endl;
     std::cout << "nb_cells: " << mesh_2.nb_cells(mesh_id_t::cells) << "\n";
diff --git a/demos/pablo/bubble_2d.cpp b/demos/pablo/bubble_2d.cpp
index 3d97eb37b..10bfcea1b 100644
--- a/demos/pablo/bubble_2d.cpp
+++ b/demos/pablo/bubble_2d.cpp
@@ -244,9 +244,9 @@ int main(int argc, char* argv[])
     double dt                        = 0.5;
 
     // Adaptation parameters
-    std::size_t start_level = 4;
-    std::size_t min_level   = 1;
-    std::size_t max_level   = 9;
+    std::size_t init_level = 4;
+    std::size_t min_level  = 1;
+    std::size_t max_level  = 9;
 
     // Output parameters
     fs::path path        = fs::current_path();
@@ -258,9 +258,9 @@ int main(int argc, char* argv[])
     app.add_option("--nb-bubbles", nb_bubbles, "Number of bubbles")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
     app.add_option("--dt", dt, "Time step")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--start-level", start_level, "Start level of AMR")->capture_default_str()->group("Adaptation parameters");
-    app.add_option("--min-level", min_level, "Minimum level of AMR")->capture_default_str()->group("Adaptation parameters");
-    app.add_option("--max-level", max_level, "Maximum level of AMR")->capture_default_str()->group("Adaptation parameters");
+    app.add_option("--initial-level", init_level, "Initial level of AMR")->capture_default_str()->group("Adaptation parameters");
+    app.add_option("--minimum-level", min_level, "Minimum level of AMR")->capture_default_str()->group("Adaptation parameters");
+    app.add_option("--maximum-level", max_level, "Maximum level of AMR")->capture_default_str()->group("Adaptation parameters");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -289,7 +289,7 @@ int main(int argc, char* argv[])
     samurai::CellArray<dim> mesh;
 
     const samurai::Box<double, dim> box(min_corner, max_corner);
-    mesh[start_level] = {start_level, box};
+    mesh[init_level] = {init_level, box};
 
     const double dt_save = Tf / static_cast<double>(nfiles);
     double t             = 0.;
diff --git a/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
index 920d331aa..6c8f464b3 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_3/main.cpp
@@ -56,8 +56,8 @@ int main(int argc, char* argv[])
     const std::size_t max_level   = 6; // <--------------------------------
 
     const samurai::Box<double, dim> box({-3}, {3});
-    Mesh<MeshConfig<dim>> mesh(box, start_level, min_level,
-                               max_level); // <--------------------------------
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).start_level(start_level);
+    Mesh<MeshConfig<dim>> mesh(config, box); // <--------------------------------
 
     auto phi = init_sol(mesh);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_3/mesh.hpp b/demos/tutorial/AMR_1D_Burgers/step_3/mesh.hpp
index 3f8cef37d..e30c6f293 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_3/mesh.hpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_3/mesh.hpp
@@ -64,14 +64,14 @@ class Mesh : public samurai::Mesh_base<Mesh<Config>, Config>
     Mesh() = default;
 
     // Constructor starting from a cell list
-    inline Mesh(const cl_type& cl, std::size_t min_level, std::size_t max_level)
-        : base_type(cl, min_level, max_level)
+    inline Mesh(const samurai::mesh_config<Config::dim>& cfg, const cl_type& cl)
+        : base_type(cfg, cl)
     {
     }
 
     // Constructor from a given box (domain)
-    inline Mesh(const samurai::Box<double, dim>& b, std::size_t start_level, std::size_t min_level, std::size_t max_level)
-        : base_type(b, start_level, min_level, max_level, 0, 1)
+    inline Mesh(samurai::mesh_config<Config::dim>& cfg, const samurai::Box<double, dim>& b)
+        : base_type(cfg.approx_box_tol(0).scaling_factor(1), b)
     {
     }
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_4/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_4/main.cpp
index 5bdbe5353..5fb129533 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_4/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_4/main.cpp
@@ -40,9 +40,6 @@ int main(int argc, char* argv[])
     fs::path path        = fs::current_path();
     std::string filename = "amr_1d_burgers_step_4";
 
-    app.add_option("--start-level", start_level, "Start level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--min-level", min_level, "Minimum level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--max-level", max_level, "Maximum level of the AMR")->capture_default_str()->group("AMR parameter");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     SAMURAI_PARSE(argc, argv);
@@ -55,7 +52,8 @@ int main(int argc, char* argv[])
     constexpr std::size_t dim = 1;
 
     const samurai::Box<double, dim> box({-3}, {3});
-    Mesh<MeshConfig<dim>> mesh(box, start_level, min_level, max_level);
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).start_level(start_level);
+    Mesh<MeshConfig<dim>> mesh(config, box);
 
     auto phi = init_sol(mesh);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp b/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
index 49a13c983..5d09621df 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_4/update_mesh.hpp
@@ -62,7 +62,7 @@ bool update_mesh(Field& f, const Tag& tag)
                                    }
                                });
 
-    mesh_t new_mesh(cell_list, mesh.min_level(), mesh.max_level());
+    mesh_t new_mesh(mesh.cfg(), cell_list);
 
     if (new_mesh == mesh)
     {
diff --git a/demos/tutorial/AMR_1D_Burgers/step_5/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_5/main.cpp
index ea3348292..8a9935b86 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_5/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_5/main.cpp
@@ -42,9 +42,6 @@ int main(int argc, char* argv[])
     fs::path path        = fs::current_path();
     std::string filename = "amr_1d_burgers_step_5";
 
-    app.add_option("--start-level", start_level, "Start level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--min-level", min_level, "Minimum level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--max-level", max_level, "Maximum level of the AMR")->capture_default_str()->group("AMR parameter");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     SAMURAI_PARSE(argc, argv);
@@ -57,7 +54,8 @@ int main(int argc, char* argv[])
     constexpr std::size_t dim = 1;
 
     const samurai::Box<double, dim> box({-3}, {3});
-    Mesh<MeshConfig<dim>> mesh(box, start_level, min_level, max_level);
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).start_level(start_level);
+    Mesh<MeshConfig<dim>> mesh(config, box);
 
     auto phi = init_sol(mesh);
 
diff --git a/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp b/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
index 48d4573fd..3a6b7715d 100644
--- a/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
+++ b/demos/tutorial/AMR_1D_Burgers/step_6/main.cpp
@@ -53,9 +53,6 @@ int main(int argc, char* argv[])
 
     app.add_option("--cfl", cfl, "The CFL")->capture_default_str()->group("Simulation parameters");
     app.add_option("--Tf", Tf, "Final time")->capture_default_str()->group("Simulation parameters");
-    app.add_option("--start-level", start_level, "Start level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--min-level", min_level, "Minimum level of the AMR")->capture_default_str()->group("AMR parameter");
-    app.add_option("--max-level", max_level, "Maximum level of the AMR")->capture_default_str()->group("AMR parameter");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
     app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output");
@@ -69,7 +66,8 @@ int main(int argc, char* argv[])
     constexpr std::size_t dim = 1;
 
     const samurai::Box<double, dim> box({-3}, {3});
-    Mesh<MeshConfig<dim>> mesh(box, start_level, min_level, max_level);
+    auto config = samurai::mesh_config<dim>().min_level(min_level).max_level(max_level).start_level(start_level);
+    Mesh<MeshConfig<dim>> mesh(config, box);
 
     auto phi = init_sol(mesh);
 
diff --git a/demos/tutorial/graduation_case_1.cpp b/demos/tutorial/graduation_case_1.cpp
index aa9997765..994467d9b 100644
--- a/demos/tutorial/graduation_case_1.cpp
+++ b/demos/tutorial/graduation_case_1.cpp
@@ -69,7 +69,7 @@ int main(int argc, char* argv[])
     fs::path path        = fs::current_path();
     std::string filename = "graduation_case_1";
 
-    app.add_option("--start-level", start_level, "where to start the mesh generator")->capture_default_str();
+    app.add_option("--starting-level", start_level, "where to start the mesh generator")->capture_default_str();
     app.add_option("--max-refinement-level", max_refinement_level, "Maximum level of the mesh generator")->capture_default_str();
     app.add_flag("--with-corner", with_corner, "Make the graduation including the diagonal")->capture_default_str();
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
diff --git a/demos/tutorial/graduation_case_2.cpp b/demos/tutorial/graduation_case_2.cpp
index 46c24b83e..da336852e 100644
--- a/demos/tutorial/graduation_case_2.cpp
+++ b/demos/tutorial/graduation_case_2.cpp
@@ -50,8 +50,8 @@ int main(int argc, char* argv[])
     fs::path path        = fs::current_path();
     std::string filename = "graduation_case_2";
 
-    app.add_option("--min-level", min_level, "Minimum level of the mesh generator")->capture_default_str();
-    app.add_option("--max-level", max_level, "Maximum level of the mesh generator")->capture_default_str();
+    app.add_option("--minimum-level", min_level, "Minimum level of the mesh generator")->capture_default_str();
+    app.add_option("--maximum-level", max_level, "Maximum level of the mesh generator")->capture_default_str();
     app.add_flag("--with-corner", with_corner, "Make the graduation including the diagonal")->capture_default_str();
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
diff --git a/demos/tutorial/graduation_case_3.cpp b/demos/tutorial/graduation_case_3.cpp
index ef4b76c49..761ddc9be 100644
--- a/demos/tutorial/graduation_case_3.cpp
+++ b/demos/tutorial/graduation_case_3.cpp
@@ -42,8 +42,8 @@ int main(int argc, char* argv[])
     fs::path path        = fs::current_path();
     std::string filename = "graduation_case_3";
 
-    app.add_option("--start-level", start_level, "where to start the mesh generator")->capture_default_str();
-    app.add_option("--max-level", max_level, "Maximum level of the mesh generator")->capture_default_str();
+    app.add_option("--starting-level", start_level, "where to start the mesh generator")->capture_default_str();
+    app.add_option("--maximum-level", max_level, "Maximum level of the mesh generator")->capture_default_str();
     app.add_flag("--with-graduation", with_graduation, "Make the mesh graduated")->capture_default_str();
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
diff --git a/demos/tutorial/proj_on_mesh.cpp b/demos/tutorial/proj_on_mesh.cpp
index 36e91fe18..965509115 100644
--- a/demos/tutorial/proj_on_mesh.cpp
+++ b/demos/tutorial/proj_on_mesh.cpp
@@ -54,15 +54,13 @@ int main()
     samurai::initialize();
 
     constexpr std::size_t dim = 3;
-    using Config              = samurai::MRConfig<dim>;
     auto box                  = samurai::Box<double, dim>({0., 0., 0.}, {1., 1., 1.});
 
-    const std::size_t min_level = 2;
-    const std::size_t max_level = 8;
-    auto mesh1                  = samurai::MRMesh<Config>(box, min_level, max_level);
-    auto f1                     = init_field(mesh1, 0);
+    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(8);
+    auto mesh1    = samurai::make_MRMesh(mesh_cfg, box);
+    auto f1       = init_field(mesh1, 0);
 
-    auto mesh2 = samurai::MRMesh<Config>(box, min_level, max_level);
+    auto mesh2 = samurai::make_MRMesh(mesh_cfg, box);
     auto f2    = init_field(mesh2, 0.1);
 
     auto adapt_1 = samurai::make_MRAdapt(f1);
diff --git a/demos/tutorial/reconstruction_1d.cpp b/demos/tutorial/reconstruction_1d.cpp
index 3e0946dc1..bb67bf86e 100644
--- a/demos/tutorial/reconstruction_1d.cpp
+++ b/demos/tutorial/reconstruction_1d.cpp
@@ -58,12 +58,7 @@ int main(int argc, char* argv[])
 {
     auto& app = samurai::initialize("1d reconstruction of an adapted solution using multiresolution", argc, argv);
 
-    constexpr size_t dim                        = 1;
-    constexpr std::size_t max_stencil_width_    = 2;
-    constexpr std::size_t graduation_width_     = 2;
-    constexpr std::size_t max_refinement_level_ = samurai::default_config::max_level;
-    constexpr std::size_t prediction_order_     = 1;
-    using MRConfig = samurai::MRConfig<dim, max_stencil_width_, graduation_width_, prediction_order_, max_refinement_level_>;
+    constexpr size_t dim = 1;
 
     Case test_case{Case::abs};
     const std::map<std::string, Case> map{
@@ -72,17 +67,11 @@ int main(int argc, char* argv[])
         {"tanh", Case::tanh}
     };
 
-    // Adaptation parameters
-    std::size_t min_level = 3;
-    std::size_t max_level = 8;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "reconstruction_1d";
 
     app.add_option("--case", test_case, "Test case")->capture_default_str()->transform(CLI::CheckedTransformer(map, CLI::ignore_case));
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
 
@@ -93,15 +82,23 @@ int main(int argc, char* argv[])
         fs::create_directory(path);
     }
 
-    using MRMesh      = samurai::MRMesh<MRConfig>;
-    using mrmesh_id_t = typename MRMesh::mesh_id_t;
-
     using UConfig = samurai::UniformConfig<dim>;
     using UMesh   = samurai::UniformMesh<UConfig>;
 
     const samurai::Box<double, dim> box({-1}, {1});
-    MRMesh mrmesh{box, min_level, max_level, 0, 1};
-    UMesh umesh{box, max_level, 0, 1};
+    // clang-format off
+    auto config = samurai::mesh_config<dim>()
+    .min_level(3)
+    .max_level(8)
+    .scaling_factor(1)
+    .graduation_width(2)
+    .max_stencil_radius(2);
+    // clang-format on
+    auto mrmesh = samurai::make_MRMesh(config, box);
+    UMesh umesh{box, mrmesh.max_level(), 0, 1};
+
+    using mrmesh_id_t = typename decltype(mrmesh)::mesh_id_t;
+
     auto u       = init(mrmesh, test_case);
     auto u_exact = init(umesh, test_case);
 
diff --git a/demos/tutorial/reconstruction_2d.cpp b/demos/tutorial/reconstruction_2d.cpp
index 8dbd50d88..93e24a2e6 100644
--- a/demos/tutorial/reconstruction_2d.cpp
+++ b/demos/tutorial/reconstruction_2d.cpp
@@ -86,12 +86,7 @@ int main(int argc, char* argv[])
 {
     auto& app = samurai::initialize("2d reconstruction of an adapted solution using multiresolution", argc, argv);
 
-    constexpr size_t dim                        = 2;
-    constexpr std::size_t max_stencil_width_    = 2;
-    constexpr std::size_t graduation_width_     = 2;
-    constexpr std::size_t max_refinement_level_ = samurai::default_config::max_level;
-    constexpr std::size_t prediction_order_     = 1;
-    using MRConfig = samurai::MRConfig<dim, max_stencil_width_, graduation_width_, prediction_order_, max_refinement_level_>;
+    constexpr size_t dim = 2;
 
     Case test_case{Case::abs};
     const std::map<std::string, Case> map{
@@ -100,17 +95,11 @@ int main(int argc, char* argv[])
         {"tanh", Case::tanh}
     };
 
-    // Adaptation parameters
-    std::size_t min_level = 3;
-    std::size_t max_level = 8;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "reconstruction_2d";
 
     app.add_option("--case", test_case, "Test case")->capture_default_str()->transform(CLI::CheckedTransformer(map, CLI::ignore_case));
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
 
@@ -121,18 +110,25 @@ int main(int argc, char* argv[])
         fs::create_directory(path);
     }
 
-    using MRMesh      = samurai::MRMesh<MRConfig>;
-    using mrmesh_id_t = typename MRMesh::mesh_id_t;
-
     using UConfig = samurai::UniformConfig<dim>;
     using UMesh   = samurai::UniformMesh<UConfig>;
 
     const samurai::Box<double, dim> box({-1, -1}, {1, 1});
-    MRMesh mrmesh{box, min_level, max_level, 0, 1};
-    UMesh umesh{box, max_level, 0, 1};
+    // clang-format off
+    auto config = samurai::mesh_config<dim>()
+    .min_level(3)
+    .max_level(8)
+    .scaling_factor(1)
+    .graduation_width(2)
+    .max_stencil_radius(2);
+    // clang-format on
+    auto mrmesh = samurai::make_MRMesh(config, box);
+    UMesh umesh{box, mrmesh.max_level(), 0, 1};
     auto u       = init(mrmesh, test_case);
     auto u_exact = init(umesh, test_case);
 
+    using mrmesh_id_t = typename decltype(mrmesh)::mesh_id_t;
+
     auto MRadaptation = samurai::make_MRAdapt(u);
     auto mra_config   = samurai::mra_config().regularity(2);
     MRadaptation(mra_config);
diff --git a/demos/tutorial/reconstruction_3d.cpp b/demos/tutorial/reconstruction_3d.cpp
index b2def1b5b..2a8ba255e 100644
--- a/demos/tutorial/reconstruction_3d.cpp
+++ b/demos/tutorial/reconstruction_3d.cpp
@@ -90,14 +90,7 @@ int main(int argc, char* argv[])
 {
     auto& app = samurai::initialize("3d reconstruction of an adapted solution using multiresolution", argc, argv);
 
-    constexpr size_t dim                     = 3;
-    constexpr std::size_t max_stencil_width_ = 2;
-    // constexpr std::size_t graduation_width_ =
-    // samurai::default_config::graduation_width;
-    constexpr std::size_t graduation_width_     = 2;
-    constexpr std::size_t max_refinement_level_ = samurai::default_config::max_level;
-    constexpr std::size_t prediction_order_     = 1;
-    using MRConfig = samurai::MRConfig<dim, max_stencil_width_, graduation_width_, prediction_order_, max_refinement_level_>;
+    constexpr size_t dim = 3;
 
     Case test_case{Case::abs};
     const std::map<std::string, Case> map{
@@ -106,17 +99,11 @@ int main(int argc, char* argv[])
         {"tanh", Case::tanh}
     };
 
-    // Adaptation parameters
-    std::size_t min_level = 2;
-    std::size_t max_level = 5;
-
     // Output parameters
     fs::path path        = fs::current_path();
     std::string filename = "reconstruction_3d";
 
     app.add_option("--case", test_case, "Test case")->capture_default_str()->transform(CLI::CheckedTransformer(map, CLI::ignore_case));
-    app.add_option("--min-level", min_level, "Minimum level of the multiresolution")->capture_default_str()->group("Multiresolution");
-    app.add_option("--max-level", max_level, "Maximum level of the multiresolution")->capture_default_str()->group("Multiresolution");
     app.add_option("--path", path, "Output path")->capture_default_str()->group("Output");
     app.add_option("--filename", filename, "File name prefix")->capture_default_str()->group("Output");
 
@@ -127,18 +114,25 @@ int main(int argc, char* argv[])
         fs::create_directory(path);
     }
 
-    using MRMesh      = samurai::MRMesh<MRConfig>;
-    using mrmesh_id_t = typename MRMesh::mesh_id_t;
-
     using UConfig = samurai::UniformConfig<dim>;
     using UMesh   = samurai::UniformMesh<UConfig>;
 
     const samurai::Box<double, dim> box({-1, -1, -1}, {1, 1, 1});
-    MRMesh mrmesh{box, min_level, max_level, 0, 1};
-    UMesh umesh{box, max_level, 0, 1};
+    // clang-format off
+    auto config = samurai::mesh_config<dim>()
+    .min_level(2)
+    .max_level(5)
+    .scaling_factor(1)
+    .graduation_width(2)
+    .max_stencil_radius(2);
+    // clang-format on
+    auto mrmesh = samurai::make_MRMesh(config, box);
+    UMesh umesh{box, mrmesh.max_level(), 0, 1};
     auto u       = init(mrmesh, test_case);
     auto u_exact = init(umesh, test_case);
 
+    using mrmesh_id_t = typename decltype(mrmesh)::mesh_id_t;
+
     auto MRadaptation = samurai::make_MRAdapt(u);
     auto mra_config   = samurai::mra_config().regularity(2);
     MRadaptation(mra_config);
diff --git a/docs/source/howto/snippet/field/scalar_field.cpp b/docs/source/howto/snippet/field/scalar_field.cpp
index 0a1794235..85c805e7c 100644
--- a/docs/source/howto/snippet/field/scalar_field.cpp
+++ b/docs/source/howto/snippet/field/scalar_field.cpp
@@ -5,10 +5,11 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field = samurai::make_scalar_field<double>("u", mesh);
 
diff --git a/docs/source/howto/snippet/field/vector_field.cpp b/docs/source/howto/snippet/field/vector_field.cpp
index 2627bd985..b1acd32f8 100644
--- a/docs/source/howto/snippet/field/vector_field.cpp
+++ b/docs/source/howto/snippet/field/vector_field.cpp
@@ -5,10 +5,11 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field = samurai::make_vector_field<double, 3>("v", mesh); // 3 components
 
diff --git a/docs/source/howto/snippet/loop/for_each_cell_field.cpp b/docs/source/howto/snippet/loop/for_each_cell_field.cpp
index da2a00e72..aea7c8a9b 100644
--- a/docs/source/howto/snippet/loop/for_each_cell_field.cpp
+++ b/docs/source/howto/snippet/loop/for_each_cell_field.cpp
@@ -6,10 +6,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field = samurai::make_scalar_field<double>("u", mesh);
 
diff --git a/docs/source/howto/snippet/loop/for_each_cell_mesh.cpp b/docs/source/howto/snippet/loop/for_each_cell_mesh.cpp
index 6213dbf40..f977bdf20 100644
--- a/docs/source/howto/snippet/loop/for_each_cell_mesh.cpp
+++ b/docs/source/howto/snippet/loop/for_each_cell_mesh.cpp
@@ -5,10 +5,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     samurai::for_each_cell(mesh,
                            [&](const auto& cell)
diff --git a/docs/source/howto/snippet/loop/for_each_interval_field.cpp b/docs/source/howto/snippet/loop/for_each_interval_field.cpp
index 41c21f086..0f0f117f7 100644
--- a/docs/source/howto/snippet/loop/for_each_interval_field.cpp
+++ b/docs/source/howto/snippet/loop/for_each_interval_field.cpp
@@ -6,10 +6,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({-1.0, -1.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 0, 2); // min level 0, max level 2
+    auto config = samurai::mesh_config<dim>().min_level(0).max_level(2);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field = samurai::make_scalar_field<double>("u", mesh);
 
diff --git a/docs/source/howto/snippet/loop/for_each_interval_mesh.cpp b/docs/source/howto/snippet/loop/for_each_interval_mesh.cpp
index 131945bb4..fb3875a92 100644
--- a/docs/source/howto/snippet/loop/for_each_interval_mesh.cpp
+++ b/docs/source/howto/snippet/loop/for_each_interval_mesh.cpp
@@ -5,10 +5,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     samurai::for_each_interval(mesh,
                                [&](std::size_t level, const auto& interval, const auto& index)
diff --git a/docs/source/howto/snippet/mesh/amrmesh.cpp b/docs/source/howto/snippet/mesh/amrmesh.cpp
index 5c377011e..e5e837f40 100644
--- a/docs/source/howto/snippet/mesh/amrmesh.cpp
+++ b/docs/source/howto/snippet/mesh/amrmesh.cpp
@@ -4,10 +4,11 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::amr::Config<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::amr::Mesh<config_t> mesh(box, 4, 2, 5); // start level 4, min level 2, max level 5
+
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5).start_level(4);
+    auto mesh   = samurai::amr::make_Mesh(config, box);
 
     return 0;
 }
diff --git a/docs/source/howto/snippet/mesh/mrmesh.cpp b/docs/source/howto/snippet/mesh/mrmesh.cpp
index 7ab70e068..dd503914b 100644
--- a/docs/source/howto/snippet/mesh/mrmesh.cpp
+++ b/docs/source/howto/snippet/mesh/mrmesh.cpp
@@ -4,10 +4,11 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     return 0;
 }
diff --git a/docs/source/howto/snippet/save/dump_mesh.cpp b/docs/source/howto/snippet/save/dump_mesh.cpp
index c5473bd0b..df2ef81e5 100644
--- a/docs/source/howto/snippet/save/dump_mesh.cpp
+++ b/docs/source/howto/snippet/save/dump_mesh.cpp
@@ -9,10 +9,10 @@ namespace fs = std::filesystem;
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field_1 = samurai::make_scalar_field<double>("u", mesh);
     auto field_2 = samurai::make_vector_field<double, 3>("v", mesh);
diff --git a/docs/source/howto/snippet/save/save_all_submeshes.cpp b/docs/source/howto/snippet/save/save_all_submeshes.cpp
index 45ac0601b..0aef00012 100644
--- a/docs/source/howto/snippet/save/save_all_submeshes.cpp
+++ b/docs/source/howto/snippet/save/save_all_submeshes.cpp
@@ -9,10 +9,10 @@ namespace fs = std::filesystem;
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field_1 = samurai::make_scalar_field<double>("u", mesh);
     auto field_2 = samurai::make_vector_field<double, 3>("v", mesh);
diff --git a/docs/source/howto/snippet/save/save_field.cpp b/docs/source/howto/snippet/save/save_field.cpp
index 19fe42689..9af2aa13e 100644
--- a/docs/source/howto/snippet/save/save_field.cpp
+++ b/docs/source/howto/snippet/save/save_field.cpp
@@ -6,10 +6,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field_1 = samurai::make_scalar_field<double>("u", mesh);
     auto field_2 = samurai::make_vector_field<double, 3>("v", mesh);
diff --git a/docs/source/howto/snippet/save/save_mesh.cpp b/docs/source/howto/snippet/save/save_mesh.cpp
index fdb3c15c2..eed964c30 100644
--- a/docs/source/howto/snippet/save/save_mesh.cpp
+++ b/docs/source/howto/snippet/save/save_mesh.cpp
@@ -5,10 +5,10 @@
 int main()
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     samurai::save("output_path", "mesh_filename", mesh);
     // or
diff --git a/docs/source/howto/snippet/timers/custom_timer.cpp b/docs/source/howto/snippet/timers/custom_timer.cpp
index 048d2a755..305b102cd 100644
--- a/docs/source/howto/snippet/timers/custom_timer.cpp
+++ b/docs/source/howto/snippet/timers/custom_timer.cpp
@@ -8,12 +8,12 @@
 int main(int argc, char** argv)
 {
     static constexpr std::size_t dim = 2;
-    using config_t                   = samurai::MRConfig<dim>;
 
     samurai::initialize("Custom timer example", argc, argv);
 
     samurai::Box<double, dim> box({0.0, 0.0}, {1.0, 1.0});
-    samurai::MRMesh<config_t> mesh(box, 2, 5); // min level 2, max level 5
+    auto config = samurai::mesh_config<dim>().min_level(2).max_level(5);
+    auto mesh   = samurai::make_MRMesh(config, box);
 
     auto field = samurai::make_scalar_field<double>("u", mesh);
 
diff --git a/docs/source/tutorial/reaction_diffusion.rst b/docs/source/tutorial/reaction_diffusion.rst
index facc13951..49e04b7f2 100644
--- a/docs/source/tutorial/reaction_diffusion.rst
+++ b/docs/source/tutorial/reaction_diffusion.rst
@@ -45,18 +45,16 @@ The mesh is created by the following code:
     double right_box = 10;
 
     // Mesh creation
-    using Config  = samurai::MRConfig<dim>;
     using Box     = samurai::Box<double, dim>;
     using point_t = typename Box::point_t;
 
-    std::size_t min_level = 0;
-    std::size_t max_level = 4;
+    auto config = samurai::mesh_config<dim>().min_level(0).max_level(4);
 
     point_t box_corner1, box_corner2;
     box_corner1.fill(left_box);
     box_corner2.fill(right_box);
     Box box(box_corner1, box_corner2);
-    samurai::MRMesh<Config> mesh{box, min_level, max_level};
+    auto mesh = samurai::make_MRMesh(config, box);
 
 Here, we have used the mesh type :code:`MRMesh`, but any other type of mesh can be chosen.
 Then, solution fields at two time steps (:math:`u_n` and :math:`u_{n+1}`) are declared, and boundary conditions are attached:
diff --git a/include/samurai/algorithm/graduation.hpp b/include/samurai/algorithm/graduation.hpp
index 1827820f1..0f2b60cc3 100644
--- a/include/samurai/algorithm/graduation.hpp
+++ b/include/samurai/algorithm/graduation.hpp
@@ -135,8 +135,6 @@ namespace samurai
 
         std::size_t max_level = mesh.max_level();
 
-        constexpr int ghost_width = mesh_t::config::graduation_width; // cppcheck-suppress unreadVariable
-
         for (std::size_t level = max_level; level > 0; --level)
         {
             /**
@@ -159,7 +157,18 @@ namespace samurai
 
             auto subset_2 = intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]);
 
-            subset_2.apply_op(tag_to_keep<ghost_width>(tag, CellFlag::refine));
+            auto ghost_width = mesh.cfg().graduation_width();
+            assert(ghost_width < 10 && "Graduation not implemented for ghost_width higher than 10");
+            // maximum ghost width is set to 9
+            static_for<1, 10>::apply(
+                [&](auto static_ghost_width_)
+                {
+                    static constexpr int static_ghost_width = static_cast<int>(static_ghost_width_());
+                    if (ghost_width == static_ghost_width)
+                    {
+                        subset_2.apply_op(tag_to_keep<static_ghost_width>(tag, CellFlag::refine));
+                    }
+                });
 
             /**
              *      K     C                          K     K
@@ -321,13 +330,13 @@ namespace samurai
 
     template <size_t dim, typename TInterval, size_t max_size, typename TCoord>
     void list_interval_to_refine_for_contiguous_boundary_cells(
-        const size_t half_stencil_width,
+        const int max_stencil_radius,
         const CellArray<dim, TInterval, max_size>& ca,
         const LevelCellArray<dim, TInterval>& domain,
         const std::array<bool, dim>& is_periodic,
         std::array<ArrayOfIntervalAndPoint<TInterval, TCoord>, CellArray<dim, TInterval, max_size>::max_size>& out)
     {
-        if (half_stencil_width == 1)
+        if (max_stencil_radius == 1)
         {
             return;
         }
@@ -348,13 +357,13 @@ namespace samurai
                     // Case where the boundary is at level L and the jump is going down to L-1:
                     //     We want to have enough contiguous boundary cells to ensure that the stencil at the lower level
                     //     won't go outside the domain.
-                    //     To ensure half_stencil_width at L-1, we need 2*half_stencil_width at level L.
+                    //     To ensure max_stencil_radius at L-1, we need 2*max_stencil_radius at level L.
                     //     However, since we project the B.C. in the first outside ghost at level L-1, we can reduce the number of
-                    //     contiguous cells by 1 at level L-1. This makes, at level L, 2*(half_stencil_width - 2) contiguous cells.
-                    //     (One cell is a real cell, the other is a ghost cell outside of the domain, which makes half_stencil_width - 2
+                    //     contiguous cells by 1 at level L-1. This makes, at level L, 2*(max_stencil_radius - 2) contiguous cells.
+                    //     (One cell is a real cell, the other is a ghost cell outside of the domain, which makes max_stencil_radius - 2
                     //     ghosts cells inside the domain).
 
-                    int n_contiguous_boundary_cells = std::max(int(half_stencil_width), 2 * (int(half_stencil_width) - 2));
+                    int n_contiguous_boundary_cells = std::max(max_stencil_radius, 2 * (max_stencil_radius - 2));
 
                     if (n_contiguous_boundary_cells > 1)
                     {
@@ -382,16 +391,16 @@ namespace samurai
 
                     // 2. Jump level --> level+1
                     // Case where the boundary is at level L and jump is going up:
-                    //    If the number of boundary contiguous cells is >= ceil(half_stencil_width/2), then there is nothing to do,
-                    //    since the half stencil at L+1 will not go out of the domain. Here, we just test if half_stencil_width > 2 by
+                    //    If the number of boundary contiguous cells is >= ceil(max_stencil_radius/2), then there is nothing to do,
+                    //    since the half stencil at L+1 will not go out of the domain. Here, we just test if max_stencil_radius > 2 by
                     //    simplicity, but at some point it would be nice to implement the real test. Otherwise, ensuring
-                    //    half_stencil_width contiguous cells at level L+1 is enough.
-                    if (half_stencil_width > 2)
+                    //    max_stencil_radius contiguous cells at level L+1 is enough.
+                    if (max_stencil_radius > 2)
                     {
                         for (size_t level = max_level - 1; level != min_level - 1; --level)
                         {
                             auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation));
-                            for (size_t i = 1; i != half_stencil_width; ++i)
+                            for (int i = 1; i != max_stencil_radius; ++i)
                             {
                                 auto refine_subset = translate(
                                                          intersection(translate(boundaryCells, -i * translation), ca[level + 1]).on(level),
@@ -410,8 +419,8 @@ namespace samurai
     }
 
     template <size_t dim, typename TInterval, typename MeshType, size_t max_size, typename TCoord>
-    void list_intervals_to_refine(const size_t grad_width,
-                                  const size_t half_stencil_width,
+    void list_intervals_to_refine(const std::size_t grad_width,
+                                  const int max_stencil_radius,
                                   const CellArray<dim, TInterval, max_size>& ca,
                                   const LevelCellArray<dim, TInterval>& domain,
                                   [[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
@@ -422,7 +431,7 @@ namespace samurai
         list_interval_to_refine_for_graduation(grad_width, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, out);
         if (!domain.empty())
         {
-            list_interval_to_refine_for_contiguous_boundary_cells(half_stencil_width, ca, domain, is_periodic, out);
+            list_interval_to_refine_for_contiguous_boundary_cells(max_stencil_radius, ca, domain, is_periodic, out);
         }
     }
 
@@ -486,8 +495,8 @@ namespace samurai
                            const LevelCellArray<dim, TInterval>& domain,
                            [[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
                            const std::array<bool, dim>& is_periodic,
-                           const size_t grad_width         = 1,
-                           const size_t half_stencil_width = 1 // half of width of the numerical scheme's stencil.
+                           const size_t grad_width      = 1,
+                           const int max_stencil_radius = 1 // half of width of the numerical scheme's stencil.
     )
     {
         using ca_type    = CellArray<dim, TInterval, max_size>;
@@ -537,7 +546,7 @@ namespace samurai
             // Then, if the non-graduated is not tagged as keep, we coarsen it
             ca_add_p.clear();
             ca_remove_p.clear();
-            list_intervals_to_refine(grad_width, half_stencil_width, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all);
+            list_intervals_to_refine(grad_width, max_stencil_radius, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all);
 
             add_p_interval.clear();
             add_p_inner_stencil.clear();
diff --git a/include/samurai/algorithm/update.hpp b/include/samurai/algorithm/update.hpp
index bd9c95823..f7cf35ace 100644
--- a/include/samurai/algorithm/update.hpp
+++ b/include/samurai/algorithm/update.hpp
@@ -4,6 +4,7 @@
 #pragma once
 
 #include <algorithm>
+#include <type_traits>
 
 #include <xtensor/containers/xfixed.hpp>
 
@@ -36,7 +37,7 @@ namespace samurai
     void update_ghost(Field& field, Fields&... fields)
     {
         using mesh_id_t                  = typename Field::mesh_t::mesh_id_t;
-        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_order;
+        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_stencil_radius;
 
         auto& mesh            = field.mesh();
         std::size_t max_level = mesh.max_level();
@@ -61,7 +62,7 @@ namespace samurai
     void update_ghost_mro(Field& field)
     {
         using mesh_id_t                  = typename Field::mesh_t::mesh_id_t;
-        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_order;
+        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_stencil_radius;
         auto& mesh                       = field.mesh();
 
         std::size_t max_level = mesh.max_level();
@@ -102,10 +103,10 @@ namespace samurai
     {
         using mesh_id_t = typename Field::mesh_t::mesh_id_t;
 
-        assert(layer > 0 && layer <= Field::mesh_t::config::max_stencil_width);
-
         auto& mesh = field.mesh();
 
+        assert(layer > 0 && layer <= mesh.max_stencil_radius());
+
         auto domain = self(mesh.domain()).on(proj_level);
 
         auto& inner = mesh.get_union()[proj_level];
@@ -125,7 +126,7 @@ namespace samurai
             // We don't want to fill by projection the ghosts that have been/will be filled by the B.C. in other directions.
             // This can happen when there is a hole in the domain.
 
-            std::size_t n_bc_ghosts = Field::mesh_t::config::max_stencil_width;
+            int n_bc_ghosts = mesh.max_stencil_radius();
             if (!field.get_bc().empty())
             {
                 n_bc_ghosts = field.get_bc().front()->stencil_size() / 2;
@@ -278,7 +279,7 @@ namespace samurai
 
         auto& mesh = field.mesh();
 
-        std::size_t n_bc_ghosts = Field::mesh_t::config::max_stencil_width;
+        int n_bc_ghosts = mesh.max_stencil_radius();
         if (!field.get_bc().empty())
         {
             n_bc_ghosts = field.get_bc().front()->stencil_size() / 2;
@@ -392,7 +393,7 @@ namespace samurai
     template <class Field>
     void update_outer_ghosts(std::size_t level, Field& field)
     {
-        static_assert(Field::mesh_t::config::prediction_order <= 1);
+        static_assert(Field::mesh_t::config::prediction_stencil_radius <= 1);
 
         constexpr std::size_t dim = Field::dim;
 
@@ -432,12 +433,12 @@ namespace samurai
                         // - the B.C. is projected onto the lower ghosts
                         // - those lower ghosts are projected onto the even lower ghosts
 
-                        std::size_t n_bc_ghosts = Field::mesh_t::config::max_stencil_width;
+                        int n_bc_ghosts = mesh.max_stencil_radius();
                         if (!field.get_bc().empty())
                         {
-                            n_bc_ghosts = field.get_bc().front()->stencil_size() / 2;
+                            n_bc_ghosts = static_cast<int>(field.get_bc().front()->stencil_size()) / 2;
                         }
-                        int max_coarse_layer = static_cast<int>(n_bc_ghosts % 2 == 0 ? n_bc_ghosts / 2 : (n_bc_ghosts + 1) / 2);
+                        int max_coarse_layer = n_bc_ghosts % 2 == 0 ? n_bc_ghosts / 2 : (n_bc_ghosts + 1) / 2;
                         for (int layer = 1; layer <= max_coarse_layer; ++layer)
                         {
                             project_bc(level, direction, layer, field); // project from level+1 to level
@@ -525,7 +526,7 @@ namespace samurai
     void update_ghost_mr(Field& field, Fields&... other_fields)
     {
         using mesh_id_t                  = typename Field::mesh_t::mesh_id_t;
-        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_order;
+        constexpr std::size_t pred_order = Field::mesh_t::config::prediction_stencil_radius;
 
         times::timers.start("ghost update");
 
@@ -602,7 +603,7 @@ namespace samurai
         using interval_t = typename Field::mesh_t::interval_t;
         using coord_t    = typename lca_t::coord_type;
 
-        static constexpr std::size_t ghost_width = Field::mesh_t::config::ghost_width;
+        int ghost_width = field.mesh().ghost_width();
 
         ArrayOfIntervalAndPoint<interval_t, coord_t> interval_list;
 
@@ -954,7 +955,6 @@ namespace samurai
         using field_value_t = typename Field::value_type;
 #endif
         using mesh_id_t        = typename Field::mesh_t::mesh_id_t;
-        using config           = typename Field::mesh_t::config;
         using lca_type         = typename Field::mesh_t::lca_type;
         using interval_value_t = typename Field::interval_t::value_t;
         using box_t            = Box<interval_value_t, Field::dim>;
@@ -977,8 +977,8 @@ namespace samurai
 
         for (std::size_t d = 0; d < dim; ++d)
         {
-            min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
-            max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+            min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
+            max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
             shift[d]      = 0;
         }
 #ifdef SAMURAI_WITH_MPI
@@ -997,23 +997,23 @@ namespace samurai
                 const auto shift_interval = shift[0];
                 const auto shift_index    = xt::view(shift, xt::range(1, _));
 
-                min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
+                min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
                 max_corner[d] = (min_indices[d] >> delta_l);
 
                 lca_type lca_min_m(level, box_t(min_corner, max_corner));
 
-                min_corner[d] = (max_indices[d] >> delta_l) - config::ghost_width;
+                min_corner[d] = (max_indices[d] >> delta_l) - mesh.ghost_width();
                 max_corner[d] = (max_indices[d] >> delta_l);
 
                 lca_type lca_max_m(level, box_t(min_corner, max_corner));
 
                 min_corner[d] = (min_indices[d] >> delta_l);
-                max_corner[d] = (min_indices[d] >> delta_l) + config::ghost_width;
+                max_corner[d] = (min_indices[d] >> delta_l) + mesh.ghost_width();
 
                 lca_type lca_min_p(level, box_t(min_corner, max_corner));
 
                 min_corner[d] = (max_indices[d] >> delta_l);
-                max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+                max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
 
                 lca_type lca_max_p(level, box_t(min_corner, max_corner));
 
@@ -1084,8 +1084,8 @@ namespace samurai
 #endif // SAMURAI_WITH_MPI
                 /* reset variables for next iterations. */
                 shift[d]      = 0;
-                min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
-                max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+                min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
+                max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
             }
         }
     }
@@ -1125,7 +1125,6 @@ namespace samurai
         using tag_value_type = typename Tag::value_type;
 #endif
         using mesh_id_t           = typename Tag::mesh_t::mesh_id_t;
-        using config              = typename Tag::mesh_t::config;
         using lca_type            = typename Tag::mesh_t::lca_type;
         using interval_value_t    = typename Tag::interval_t::value_t;
         using box_t               = Box<interval_value_t, Tag::dim>;
@@ -1147,8 +1146,8 @@ namespace samurai
         for (std::size_t d = 0; d < dim; ++d)
         {
             shift[d]      = 0;
-            min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
-            max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+            min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
+            max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
         }
 #ifdef SAMURAI_WITH_MPI
         using tag_value_type = typename Tag::value_type;
@@ -1167,23 +1166,23 @@ namespace samurai
                 const auto shift_interval = shift[0];
                 const auto shift_index    = xt::view(shift, xt::range(1, _));
 
-                min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
+                min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
                 max_corner[d] = (min_indices[d] >> delta_l);
 
                 lca_type lca_min_m(level, box_t(min_corner, max_corner));
 
-                min_corner[d] = (max_indices[d] >> delta_l) - config::ghost_width;
+                min_corner[d] = (max_indices[d] >> delta_l) - mesh.ghost_width();
                 max_corner[d] = (max_indices[d] >> delta_l);
 
                 lca_type lca_max_m(level, box_t(min_corner, max_corner));
 
                 min_corner[d] = (min_indices[d] >> delta_l);
-                max_corner[d] = (min_indices[d] >> delta_l) + config::ghost_width;
+                max_corner[d] = (min_indices[d] >> delta_l) + mesh.ghost_width();
 
                 lca_type lca_min_p(level, box_t(min_corner, max_corner));
 
                 min_corner[d] = (max_indices[d] >> delta_l);
-                max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+                max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
 
                 lca_type lca_max_p(level, box_t(min_corner, max_corner));
 
@@ -1313,8 +1312,8 @@ namespace samurai
 #endif // SAMURAI_WITH_MPI
                 /* reset variables for next iterations. */
                 shift[d]      = 0;
-                min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
-                max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+                min_corner[d] = (min_indices[d] >> delta_l) - mesh.ghost_width();
+                max_corner[d] = (max_indices[d] >> delta_l) + mesh.ghost_width();
             }
         }
     }
diff --git a/include/samurai/amr/mesh.hpp b/include/samurai/amr/mesh.hpp
index e5b7dbd5d..5b171598e 100644
--- a/include/samurai/amr/mesh.hpp
+++ b/include/samurai/amr/mesh.hpp
@@ -25,25 +25,6 @@ namespace samurai::amr
         // reference = cells_and_ghosts
     };
 
-    template <std::size_t dim_,
-              std::size_t max_stencil_width_    = default_config::ghost_width,
-              std::size_t graduation_width_     = default_config::graduation_width,
-              std::size_t max_refinement_level_ = default_config::max_level,
-              std::size_t prediction_order_     = default_config::prediction_order,
-              class TInterval                   = default_config::interval_t>
-    struct Config
-    {
-        static constexpr std::size_t dim                  = dim_;
-        static constexpr std::size_t max_refinement_level = max_refinement_level_;
-        static constexpr int max_stencil_width            = max_stencil_width_;
-        static constexpr int prediction_order             = prediction_order_;
-        static constexpr int ghost_width      = std::max(static_cast<int>(max_stencil_width), static_cast<int>(prediction_order));
-        static constexpr int graduation_width = graduation_width_;
-
-        using interval_t = TInterval;
-        using mesh_id_t  = AMR_Id;
-    };
-
     /////////////////////////
     // AMR mesh definition //
     /////////////////////////
@@ -69,11 +50,13 @@ namespace samurai::amr
         Mesh() = default;
         Mesh(const ca_type& ca, const self_type& ref_mesh);
         Mesh(const cl_type& cl, const self_type& ref_mesh);
-        Mesh(const cl_type& cl, std::size_t min_level, std::size_t max_level);
-        Mesh(const ca_type& ca, std::size_t min_level, std::size_t max_level);
-        Mesh(const Box<double, dim>& b, std::size_t start_level, std::size_t min_level, std::size_t max_level);
+        Mesh(const mesh_config<Config::dim>& config, const cl_type& cl);
+        Mesh(const mesh_config<Config::dim>& config, const ca_type& ca);
+        Mesh(const mesh_config<Config::dim>& config, const Box<double, dim>& b);
 
         void update_sub_mesh_impl();
+
+        using base_type::cfg;
     };
 
     /////////////////////////////
@@ -93,20 +76,20 @@ namespace samurai::amr
     }
 
     template <class Config>
-    inline Mesh<Config>::Mesh(const cl_type& cl, std::size_t min_level, std::size_t max_level)
-        : base_type(cl, min_level, max_level)
+    inline Mesh<Config>::Mesh(const mesh_config<Config::dim>& config, const cl_type& cl)
+        : base_type(config, cl)
     {
     }
 
     template <class Config>
-    inline Mesh<Config>::Mesh(const ca_type& ca, std::size_t min_level, std::size_t max_level)
-        : base_type(ca, min_level, max_level)
+    inline Mesh<Config>::Mesh(const mesh_config<Config::dim>& config, const ca_type& ca)
+        : base_type(config, ca)
     {
     }
 
     template <class Config>
-    inline Mesh<Config>::Mesh(const Box<double, dim>& b, std::size_t start_level, std::size_t min_level, std::size_t max_level)
-        : base_type(b, start_level, min_level, max_level)
+    inline Mesh<Config>::Mesh(const mesh_config<Config::dim>& config, const Box<double, dim>& b)
+        : base_type(config, b)
     {
     }
 
@@ -114,15 +97,18 @@ namespace samurai::amr
     inline void Mesh<Config>::update_sub_mesh_impl()
     {
         cl_type cl;
+        auto ghost_width = cfg().ghost_width();
         for_each_interval(this->cells()[mesh_id_t::cells],
                           [&](std::size_t level, const auto& interval, const auto& index_yz)
                           {
                               lcl_type& lcl = cl[level];
-                              static_nested_loop<dim - 1, -config::ghost_width, config::ghost_width + 1>(
+                              static_nested_loop<dim - 1>(
+                                  -ghost_width,
+                                  ghost_width + 1,
                                   [&](auto stencil)
                                   {
                                       auto index = xt::eval(index_yz + stencil);
-                                      lcl[index].add_interval({interval.start - config::ghost_width, interval.end + config::ghost_width});
+                                      lcl[index].add_interval({interval.start - ghost_width, interval.end + ghost_width});
                                   });
                           });
         this->cells()[mesh_id_t::cells_and_ghosts] = {cl, false};
@@ -176,11 +162,12 @@ namespace samurai::amr
                 [&](const auto& interval, const auto& index_yz)
                 {
                     // add ghosts for the prediction
-                    static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                    static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>(
                         [&](auto stencil)
                         {
                             auto index = xt::eval(index_yz + stencil);
-                            lcl[index].add_interval({interval.start - config::prediction_order, interval.end + config::prediction_order});
+                            lcl[index].add_interval(
+                                {interval.start - config::prediction_stencil_radius, interval.end + config::prediction_stencil_radius});
                         });
                 });
         }
@@ -199,7 +186,42 @@ namespace samurai::amr
             this->cells()[mesh_id_t::all_cells][level] = {lcl};
         }
     }
-}
+
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, AMR_Id>>
+    auto make_empty_Mesh(const mesh_config_t&)
+    {
+        return Mesh<complete_mesh_config_t>();
+    }
+
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, AMR_Id>>
+    auto make_Mesh(const mesh_config_t& cfg, const typename Mesh<complete_mesh_config_t>::cl_type& cl)
+    {
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+
+        return Mesh<complete_mesh_config_t>(mesh_cfg, cl);
+    }
+
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, AMR_Id>>
+    auto make_Mesh(const mesh_config_t& cfg, const typename Mesh<complete_mesh_config_t>::ca_type& ca)
+    {
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+
+        return Mesh<complete_mesh_config_t>(mesh_cfg, ca);
+    }
+
+    template <class mesh_config_t>
+    auto make_Mesh(const mesh_config_t& cfg, const samurai::Box<double, mesh_config_t::dim>& b)
+    {
+        using complete_cfg_t = complete_mesh_config<mesh_config_t, AMR_Id>;
+
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+
+        return Mesh<complete_cfg_t>(mesh_cfg, b);
+    }
+} // namespace samurai::amr
 
 template <>
 struct fmt::formatter<samurai::amr::AMR_Id> : formatter<string_view>
diff --git a/include/samurai/arguments.hpp b/include/samurai/arguments.hpp
index 0b00c62c3..dea8ebe4d 100644
--- a/include/samurai/arguments.hpp
+++ b/include/samurai/arguments.hpp
@@ -8,6 +8,13 @@ namespace samurai
 {
     namespace args
     {
+        // Mesh arguments
+        static std::size_t min_level        = std::numeric_limits<std::size_t>::max();
+        static std::size_t max_level        = std::numeric_limits<std::size_t>::max();
+        static std::size_t start_level      = std::numeric_limits<std::size_t>::max();
+        static std::size_t graduation_width = std::numeric_limits<std::size_t>::max();
+        static int max_stencil_radius       = std::numeric_limits<int>::max();
+
         static bool timers = false;
 #ifdef SAMURAI_WITH_MPI
         static bool dont_redirect_output = false;
@@ -24,6 +31,12 @@ namespace samurai
 
     inline void read_samurai_arguments(CLI::App& app, int& argc, char**& argv)
     {
+        app.add_option("--min-level", args::min_level, "The minimum level of the mesh")->group("SAMURAI");
+        app.add_option("--max-level", args::max_level, "The maximum level of the mesh")->group("SAMURAI");
+        app.add_option("--start-level", args::start_level, "Start level of AMR")->group("SAMURAI");
+        app.add_option("--graduation-width", args::graduation_width, "The graduation width of the mesh")->group("SAMURAI");
+        app.add_option("--max-stencil-radius", args::max_stencil_radius, "The maximum number of neighbour in each direction")->group("SAMURAI");
+
 #ifdef SAMURAI_WITH_MPI
         app.add_flag("--dont-redirect-output", args::dont_redirect_output, "Redirect the output for all ranks different of 0")
             ->capture_default_str()
diff --git a/include/samurai/bc/apply_field_bc.hpp b/include/samurai/bc/apply_field_bc.hpp
index 37d9689eb..c75eb7f12 100644
--- a/include/samurai/bc/apply_field_bc.hpp
+++ b/include/samurai/bc/apply_field_bc.hpp
@@ -366,25 +366,25 @@ namespace samurai
     template <class Field>
     void update_further_ghosts_by_polynomial_extrapolation(std::size_t level, const DirectionVector<Field::dim>& direction, Field& field)
     {
-        static constexpr std::size_t ghost_width                     = Field::mesh_t::config::ghost_width;
+        int ghost_width                                              = field.mesh().ghost_width();
         static constexpr std::size_t max_stencil_size_implemented_PE = PolynomialExtrapolation<Field, 2>::max_stencil_size_implemented_PE;
 
         // 1. We fill the ghosts that are further than those filled by the B.C. (where there are boundary cells)
 
-        std::size_t ghost_layers_filled_by_bc = 0;
+        int ghost_layers_filled_by_bc = 0;
         for (auto& bc : field.get_bc())
         {
             ghost_layers_filled_by_bc = std::max(ghost_layers_filled_by_bc, bc->stencil_size() / 2);
         }
 
         // We populate the ghosts sequentially from the closest to the farthest.
-        for (std::size_t ghost_layer = ghost_layers_filled_by_bc + 1; ghost_layer <= ghost_width; ++ghost_layer)
+        for (int ghost_layer = ghost_layers_filled_by_bc + 1; ghost_layer <= ghost_width; ++ghost_layer)
         {
-            std::size_t stencil_s = 2 * ghost_layer;
-            static_for<2, std::min(max_stencil_size_implemented_PE, 2 * ghost_width) + 1>::apply(
-                [&](auto integral_constant_i)
+            int stencil_s = 2 * ghost_layer;
+            static_for<2, max_stencil_size_implemented_PE + 1>::apply(
+                [&](auto stencil_size_)
                 {
-                    static constexpr std::size_t stencil_size = integral_constant_i();
+                    static constexpr int stencil_size = static_cast<int>(stencil_size_());
 
                     if constexpr (stencil_size % 2 == 0) // (because PolynomialExtrapolation is only implemented for even stencil_size)
                     {
@@ -405,13 +405,13 @@ namespace samurai
 
         const std::size_t ghost_layers_filled_by_projection_bc = 1;
 
-        for (std::size_t ghost_layer = ghost_layers_filled_by_projection_bc + 1; ghost_layer <= ghost_width; ++ghost_layer)
+        for (int ghost_layer = ghost_layers_filled_by_projection_bc + 1; ghost_layer <= ghost_width; ++ghost_layer)
         {
-            std::size_t stencil_s = 2 * ghost_layer;
-            static_for<2, std::min(max_stencil_size_implemented_PE, 2 * ghost_width) + 1>::apply(
-                [&](auto integral_constant_i)
+            int stencil_s = 2 * ghost_layer;
+            static_for<2, max_stencil_size_implemented_PE + 1>::apply(
+                [&](auto stencil_size_)
                 {
-                    static constexpr std::size_t stencil_size = integral_constant_i();
+                    static constexpr int stencil_size = static_cast<int>(stencil_size_());
 
                     if constexpr (stencil_size % 2 == 0) // (because PolynomialExtrapolation is only implemented for even stencil_size)
                     {
diff --git a/include/samurai/bc/bc.hpp b/include/samurai/bc/bc.hpp
index 7eb9d56f0..97247aeeb 100644
--- a/include/samurai/bc/bc.hpp
+++ b/include/samurai/bc/bc.hpp
@@ -31,32 +31,31 @@
         return line_stencil<dim, 0, STENCIL_SIZE>();                                                                                      \
     }
 
-#define INIT_BC(NAME, STENCIL_SIZE)                                                                                                \
-    using base_t      = samurai::Bc<Field>;                                                                                        \
-    using cell_t      = typename base_t::cell_t;                                                                                   \
-    using value_t     = typename base_t::value_t;                                                                                  \
-    using direction_t = typename base_t::direction_t;                                                                              \
-    using base_t::base_t;                                                                                                          \
-    using base_t::dim;                                                                                                             \
-    using base_t::get_apply_function;                                                                                              \
-    using base_t::get_stencil;                                                                                                     \
-                                                                                                                                   \
-    using stencil_t               = samurai::Stencil<STENCIL_SIZE, dim>;                                                           \
-    using constant_stencil_size_t = std::integral_constant<std::size_t, STENCIL_SIZE>;                                             \
-    using stencil_cells_t         = std::array<cell_t, STENCIL_SIZE>;                                                              \
-    using apply_function_t        = std::function<void(Field&, const std::array<cell_t, STENCIL_SIZE>&, const value_t&)>;          \
-                                                                                                                                   \
-    static_assert(STENCIL_SIZE <= base_t::max_stencil_size_implemented, "The stencil size is too large.");                         \
-    static_assert(Field::mesh_t::config::ghost_width >= STENCIL_SIZE / 2, "Not enough ghost layers for this boundary condition."); \
-                                                                                                                                   \
-    std::unique_ptr<base_t> clone() const override                                                                                 \
-    {                                                                                                                              \
-        return std::make_unique<NAME>(*this);                                                                                      \
-    }                                                                                                                              \
-                                                                                                                                   \
-    std::size_t stencil_size() const override                                                                                      \
-    {                                                                                                                              \
-        return STENCIL_SIZE;                                                                                                       \
+#define INIT_BC(NAME, STENCIL_SIZE)                                                                                       \
+    using base_t      = samurai::Bc<Field>;                                                                               \
+    using cell_t      = typename base_t::cell_t;                                                                          \
+    using value_t     = typename base_t::value_t;                                                                         \
+    using direction_t = typename base_t::direction_t;                                                                     \
+    using base_t::base_t;                                                                                                 \
+    using base_t::dim;                                                                                                    \
+    using base_t::get_apply_function;                                                                                     \
+    using base_t::get_stencil;                                                                                            \
+                                                                                                                          \
+    using stencil_t               = samurai::Stencil<STENCIL_SIZE, dim>;                                                  \
+    using constant_stencil_size_t = std::integral_constant<std::size_t, STENCIL_SIZE>;                                    \
+    using stencil_cells_t         = std::array<cell_t, STENCIL_SIZE>;                                                     \
+    using apply_function_t        = std::function<void(Field&, const std::array<cell_t, STENCIL_SIZE>&, const value_t&)>; \
+                                                                                                                          \
+    static_assert(STENCIL_SIZE <= base_t::max_stencil_size_implemented, "The stencil size is too large.");                \
+                                                                                                                          \
+    std::unique_ptr<base_t> clone() const override                                                                        \
+    {                                                                                                                     \
+        return std::make_unique<NAME>(*this);                                                                             \
+    }                                                                                                                     \
+                                                                                                                          \
+    int stencil_size() const override                                                                                     \
+    {                                                                                                                     \
+        return STENCIL_SIZE;                                                                                              \
     }
 
 namespace samurai
@@ -592,9 +591,9 @@ namespace samurai
         Bc& operator=(Bc&& bc) noexcept = default;
 
         virtual std::unique_ptr<Bc> clone() const = 0;
-        virtual std::size_t stencil_size() const  = 0;
+        virtual int stencil_size() const          = 0;
 
-        static constexpr std::size_t max_stencil_size_implemented = 10; // cppcheck-suppress unusedStructMember
+        static constexpr int max_stencil_size_implemented = 10; // cppcheck-suppress unusedStructMember
         APPLY_AND_STENCIL_FUNCTIONS(1)
         APPLY_AND_STENCIL_FUNCTIONS(2)
         APPLY_AND_STENCIL_FUNCTIONS(3)
diff --git a/include/samurai/boundary.hpp b/include/samurai/boundary.hpp
index ac505f65e..8977c27ed 100644
--- a/include/samurai/boundary.hpp
+++ b/include/samurai/boundary.hpp
@@ -49,7 +49,7 @@ namespace samurai
     }
 
     template <class Mesh>
-    auto domain_boundary_outer_layer(const Mesh& mesh, std::size_t level, std::size_t layer_width)
+    auto domain_boundary_outer_layer(const Mesh& mesh, std::size_t level, int layer_width)
     {
         using lca_t = typename Mesh::lca_type;
         using lcl_t = typename Mesh::lcl_type;
@@ -62,7 +62,7 @@ namespace samurai
             [&](const auto& direction)
             {
                 auto inner_boundary = domain_boundary(mesh, level, direction);
-                for (std::size_t layer = 1; layer <= layer_width; ++layer)
+                for (int layer = 1; layer <= layer_width; ++layer)
                 {
                     auto outer_layer = difference(translate(inner_boundary, layer * direction), domain);
                     outer_layer(
@@ -76,8 +76,7 @@ namespace samurai
     }
 
     template <class Mesh>
-    auto
-    domain_boundary_outer_layer(const Mesh& mesh, std::size_t level, const DirectionVector<Mesh::dim>& direction, std::size_t layer_width)
+    auto domain_boundary_outer_layer(const Mesh& mesh, std::size_t level, const DirectionVector<Mesh::dim>& direction, int layer_width)
     {
         using lca_t = typename Mesh::lca_type;
         using lcl_t = typename Mesh::lcl_type;
@@ -87,7 +86,7 @@ namespace samurai
         lcl_t outer_boundary_lcl(level, mesh.origin_point(), mesh.scaling_factor());
 
         auto inner_boundary = domain_boundary(mesh, level, direction);
-        for (std::size_t layer = 1; layer <= layer_width; ++layer)
+        for (int layer = 1; layer <= layer_width; ++layer)
         {
             auto outer_layer = difference(translate(inner_boundary, layer * direction), domain);
             outer_layer(
diff --git a/include/samurai/domain_builder.hpp b/include/samurai/domain_builder.hpp
index 8c9ebba6b..60452f42b 100644
--- a/include/samurai/domain_builder.hpp
+++ b/include/samurai/domain_builder.hpp
@@ -3,6 +3,8 @@
 
 #pragma once
 
+#include <numeric>
+
 #include "box.hpp"
 
 namespace samurai
@@ -81,15 +83,20 @@ namespace samurai
         {
             double largest_subdivision = m_added_boxes[0].min_length();
 
-            // The largest subdivision must be smaller than the smallest legnth of all boxes
-            for (const auto& box : m_added_boxes)
-            {
-                largest_subdivision = gcd_float(largest_subdivision, box.min_length());
-            }
-            for (const auto& box : m_removed_boxes)
+            auto gcd_loop = [](double largest_subdiv, const auto& boxes)
             {
-                largest_subdivision = gcd_float(largest_subdivision, box.min_length());
-            }
+                return std::accumulate(boxes.begin(),
+                                       boxes.end(),
+                                       largest_subdiv,
+                                       [](double l, const auto& box)
+                                       {
+                                           return gcd_float(l, box.min_length());
+                                       });
+            };
+
+            // The largest subdivision must be smaller than the smallest length of all boxes
+            largest_subdivision = gcd_loop(largest_subdivision, m_added_boxes);
+            largest_subdivision = gcd_loop(largest_subdivision, m_removed_boxes);
 
             // The largest subdivision must be smaller than the smallest length of all differences
             for (const auto& box : m_added_boxes)
@@ -99,10 +106,7 @@ namespace samurai
                     if (rbox.intersects(box))
                     {
                         std::vector<box_t> diff = box.difference(rbox);
-                        for (const auto& dbox : diff)
-                        {
-                            largest_subdivision = gcd_float(largest_subdivision, dbox.min_length());
-                        }
+                        largest_subdivision     = gcd_loop(largest_subdivision, diff);
                     }
                 }
             }
diff --git a/include/samurai/field.hpp b/include/samurai/field.hpp
index 1f037057f..340b8dc48 100644
--- a/include/samurai/field.hpp
+++ b/include/samurai/field.hpp
@@ -702,6 +702,14 @@ namespace samurai
     template <class Bc_derived>
     inline auto VectorField<mesh_t, value_t, n_comp_, SOA>::attach_bc(const Bc_derived& bc)
     {
+        if (bc.stencil_size() > this->mesh().cfg().max_stencil_size())
+        {
+            std::cerr << "The stencil size required by this boundary condition (" << bc.stencil_size()
+                      << ") is larger than the max_stencil_size parameter of the mesh (" << this->mesh().cfg().max_stencil_size()
+                      << ").\nYou can set it with mesh_config.max_stencil_radius(" << bc.stencil_size() / 2
+                      << ") or mesh_config.max_stencil_size(" << bc.stencil_size() << ")." << std::endl;
+            exit(EXIT_FAILURE);
+        }
         p_bc.push_back(bc.clone());
         return p_bc.back().get();
     }
@@ -1331,6 +1339,14 @@ namespace samurai
     template <class Bc_derived>
     inline auto ScalarField<mesh_t, value_t>::attach_bc(const Bc_derived& bc)
     {
+        if (bc.stencil_size() > this->mesh().cfg().max_stencil_size())
+        {
+            std::cerr << "The stencil size required by this boundary condition (" << bc.stencil_size()
+                      << ") is larger than the max_stencil_size parameter of the mesh (" << this->mesh().cfg().max_stencil_size()
+                      << ").\nYou can set it with mesh_config.max_stencil_radius(" << bc.stencil_size() / 2
+                      << ") or mesh_config.max_stencil_size(" << bc.stencil_size() << ")." << std::endl;
+            exit(EXIT_FAILURE);
+        }
         p_bc.push_back(bc.clone());
         return p_bc.back().get();
     }
diff --git a/include/samurai/io/restart.hpp b/include/samurai/io/restart.hpp
index a4a1259d8..3343acc58 100644
--- a/include/samurai/io/restart.hpp
+++ b/include/samurai/io/restart.hpp
@@ -415,7 +415,8 @@ namespace samurai
 
         ca_type ca;
         load(file, ca);
-        Mesh new_mesh{ca, min_level, max_level};
+        auto mesh_cfg = mesh_config<Mesh::dim>().min_level(min_level).max_level(max_level).disable_args_parse();
+        Mesh new_mesh{mesh_cfg, ca};
         std::swap(mesh, new_mesh);
         load_fields(file, mesh, fields...);
     }
diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp
index 1fa2192f7..8a39c6626 100644
--- a/include/samurai/mesh.hpp
+++ b/include/samurai/mesh.hpp
@@ -11,6 +11,7 @@
 #include "cell_array.hpp"
 #include "cell_list.hpp"
 #include "domain_builder.hpp"
+#include "mesh_config.hpp"
 #include "static_algorithm.hpp"
 #include "stencil.hpp"
 #include "subset/node.hpp"
@@ -98,12 +99,18 @@ namespace samurai
         std::size_t min_level() const;
         std::size_t& min_level();
 
+        std::size_t graduation_width() const;
+        int ghost_width() const;
+        int max_stencil_radius() const;
+        auto& cfg() const;
+
         auto& origin_point() const;
         void set_origin_point(const coords_t& origin_point);
         double scaling_factor() const;
         void set_scaling_factor(double scaling_factor);
         void scale_domain(double domain_scaling_factor);
         double cell_length(std::size_t level) const;
+        double min_cell_length() const;
         const lca_type& domain() const;
         const lca_type& subdomain() const;
         const ca_type& get_union() const;
@@ -155,27 +162,34 @@ namespace samurai
         Mesh_base() = default; // cppcheck-suppress uninitMemberVar
         Mesh_base(const ca_type& ca, const self_type& ref_mesh);
         Mesh_base(const cl_type& cl, const self_type& ref_mesh);
-        Mesh_base(const cl_type& cl, std::size_t min_level, std::size_t max_level);
-        Mesh_base(const ca_type& ca, std::size_t min_level, std::size_t max_level);
-        Mesh_base(const samurai::Box<double, dim>& b,
-                  std::size_t start_level,
-                  std::size_t min_level,
-                  std::size_t max_level,
-                  double approx_box_tol = lca_type::default_approx_box_tol,
-                  double scaling_factor = 0);
+        Mesh_base(const mesh_config<Config::dim>& config, const cl_type& cl);
+        Mesh_base(const mesh_config<Config::dim>& config, const ca_type& ca);
+        Mesh_base(const mesh_config<Config::dim>& config, const samurai::Box<double, dim>& b);
+
+        Mesh_base(const mesh_config<Config::dim>& config, const samurai::DomainBuilder<dim>& domain_builder);
+
+        // cppcheck-suppress uninitMemberVar
+        Mesh_base(const samurai::Box<double, dim>&, std::size_t, std::size_t, std::size_t, double, double)
+        {
+            std::cerr << "Delete min_level and max_level from CLI11 options and use mesh_config object and the make_Mesh function"
+                      << std::endl;
+            exit(EXIT_FAILURE);
+        }
+
         Mesh_base(const samurai::DomainBuilder<dim>& domain_builder,
                   std::size_t start_level,
                   std::size_t min_level,
                   std::size_t max_level,
                   double approx_box_tol = lca_type::default_approx_box_tol,
                   double scaling_factor = 0);
-        Mesh_base(const samurai::Box<double, dim>& b,
-                  std::size_t start_level,
-                  std::size_t min_level,
-                  std::size_t max_level,
-                  const std::array<bool, dim>& periodic,
-                  double approx_box_tol = lca_type::default_approx_box_tol,
-                  double scaling_factor = 0);
+
+        // cppcheck-suppress uninitMemberVar
+        Mesh_base(const samurai::Box<double, dim>&, std::size_t, std::size_t, std::size_t, const std::array<bool, dim>&, double, double)
+        {
+            std::cerr << "Delete min_level and max_level from CLI11 options and use mesh_config object and the make_Mesh function"
+                      << std::endl;
+            exit(EXIT_FAILURE);
+        }
 
         derived_type& derived_cast() & noexcept;
         const derived_type& derived_cast() const& noexcept;
@@ -201,15 +215,14 @@ namespace samurai
 
         lca_type m_domain;
         lca_type m_subdomain;
-        std::size_t m_min_level;
-        std::size_t m_max_level;
-        std::array<bool, dim> m_periodic;
         mesh_t m_cells;
         ca_type m_union;
         std::vector<lca_type> m_corners;
         // std::vector<int> m_neighbouring_ranks;
         std::vector<mpi_subdomain_t> m_mpi_neighbourhood;
 
+        mesh_config<dim> m_config;
+
 #ifdef SAMURAI_WITH_MPI
         friend class boost::serialization::access;
 
@@ -223,8 +236,7 @@ namespace samurai
             ar & m_domain;
             ar & m_subdomain;
             ar & m_union;
-            ar & m_min_level;
-            ar & m_max_level;
+            ar & m_config;
         }
 #endif
     };
@@ -248,24 +260,18 @@ namespace samurai
     }
 
     template <class D, class Config>
-    inline Mesh_base<D, Config>::Mesh_base(const samurai::Box<double, dim>& b,
-                                           std::size_t start_level,
-                                           std::size_t min_level,
-                                           std::size_t max_level,
-                                           double approx_box_tol,
-                                           double scaling_factor_)
-        : m_domain{start_level, b, approx_box_tol, scaling_factor_}
-        , m_min_level{min_level}
-        , m_max_level{max_level}
+    inline Mesh_base<D, Config>::Mesh_base(const mesh_config<Config::dim>& config, const samurai::Box<double, dim>& b)
+        : m_domain{config.start_level(), b, config.approx_box_tol(), config.scaling_factor()}
+        , m_config(config)
     {
-        assert(min_level <= max_level);
-        m_periodic.fill(false);
-
 #ifdef SAMURAI_WITH_MPI
-        partition_mesh(start_level, b);
+        partition_mesh(m_config.start_level(), b);
         // load_balancing();
 #else
-        this->m_cells[mesh_id_t::cells][start_level] = {start_level, b, approx_box_tol, scaling_factor_};
+        this->m_cells[mesh_id_t::cells][m_config.start_level()] = {m_config.start_level(),
+                                                                   b,
+                                                                   m_config.approx_box_tol(),
+                                                                   m_config.scaling_factor()};
 #endif
         construct_subdomain();
         construct_union();
@@ -279,28 +285,36 @@ namespace samurai
     }
 
     template <class D, class Config>
-    Mesh_base<D, Config>::Mesh_base(const samurai::DomainBuilder<dim>& domain_builder,
-                                    [[maybe_unused]] std::size_t start_level,
-                                    std::size_t min_level,
-                                    std::size_t max_level,
-                                    [[maybe_unused]] double approx_box_tol,
-                                    double scaling_factor_)
-        : m_min_level{min_level}
-        , m_max_level{max_level}
+    Mesh_base<D, Config>::Mesh_base(const mesh_config<Config::dim>& config, const samurai::DomainBuilder<dim>& domain_builder)
+        : m_config(config)
     {
-        assert(min_level <= max_level);
-        m_periodic.fill(false);
+        if (std::any_of(config.periodic().begin(),
+                        config.periodic().end(),
+                        [](bool b)
+                        {
+                            return b;
+                        }))
+        {
+            std::cerr << "Periodicity is not implemented with DomainBuilder." << std::endl;
+            exit(EXIT_FAILURE);
+        }
 
 #ifdef SAMURAI_WITH_MPI
         std::cerr << "MPI is not implemented with DomainBuilder." << std::endl;
-        std::exit(1);
-        // partition_mesh(start_level, b);
+        std::exit(EXIT_FAILURE);
+        // partition_mesh(m_config.start_level(), b);
         //  load_balancing();
+
 #else
+        double scaling_factor_ = m_config.scaling_factor();
         compute_scaling_factor(domain_builder, scaling_factor_);
+        m_config.scaling_factor() = scaling_factor_;
 
         // Build the domain by adding and removing boxes
-        cl_type domain_cl = construct_initial_mesh(domain_builder, start_level, approx_box_tol, scaling_factor_);
+        cl_type domain_cl = construct_initial_mesh(domain_builder,
+                                                   m_config.start_level(),
+                                                   m_config.approx_box_tol(),
+                                                   m_config.scaling_factor());
 
         this->m_cells[mesh_id_t::cells] = {domain_cl, false};
 #endif
@@ -313,50 +327,13 @@ namespace samurai
         update_mesh_neighbour();
 
         set_origin_point(domain_builder.origin_point());
-        set_scaling_factor(scaling_factor_);
-    }
-
-    template <class D, class Config>
-    inline Mesh_base<D, Config>::Mesh_base(const samurai::Box<double, dim>& b,
-                                           std::size_t start_level,
-                                           std::size_t min_level,
-                                           std::size_t max_level,
-                                           const std::array<bool, dim>& periodic,
-                                           double approx_box_tol,
-                                           double scaling_factor_)
-        : m_domain{start_level, b, approx_box_tol, scaling_factor_}
-        , m_min_level{min_level}
-        , m_max_level{max_level}
-        , m_periodic{periodic}
-    {
-        assert(min_level <= max_level);
-
-#ifdef SAMURAI_WITH_MPI
-        partition_mesh(start_level, b);
-        // load_balancing();
-#else
-        this->m_cells[mesh_id_t::cells][start_level] = {start_level, b, approx_box_tol, scaling_factor_};
-#endif
-
-        construct_subdomain();
-        construct_union();
-        update_sub_mesh();
-        construct_corners();
-        renumbering();
-        update_mesh_neighbour();
-
-        set_origin_point(origin_point());
-        set_scaling_factor(scaling_factor());
+        set_scaling_factor(m_config.scaling_factor());
     }
 
     template <class D, class Config>
-    inline Mesh_base<D, Config>::Mesh_base(const cl_type& cl, std::size_t min_level, std::size_t max_level)
-        : m_min_level{min_level}
-        , m_max_level{max_level}
+    inline Mesh_base<D, Config>::Mesh_base(const mesh_config<Config::dim>& config, const cl_type& cl)
+        : m_config(config)
     {
-        m_periodic.fill(false);
-        assert(min_level <= max_level);
-
         this->m_cells[mesh_id_t::cells] = {cl};
 
         construct_subdomain();
@@ -372,13 +349,9 @@ namespace samurai
     }
 
     template <class D, class Config>
-    inline Mesh_base<D, Config>::Mesh_base(const ca_type& ca, std::size_t min_level, std::size_t max_level)
-        : m_min_level{min_level}
-        , m_max_level{max_level}
+    inline Mesh_base<D, Config>::Mesh_base(const mesh_config<Config::dim>& config, const ca_type& ca)
+        : m_config(config)
     {
-        m_periodic.fill(false);
-        assert(min_level <= max_level);
-
         this->m_cells[mesh_id_t::cells] = ca;
 
         construct_subdomain();
@@ -408,11 +381,8 @@ namespace samurai
     template <class D, class Config>
     inline Mesh_base<D, Config>::Mesh_base(const ca_type& ca, const self_type& ref_mesh)
         : m_domain(ref_mesh.m_domain)
-        , m_min_level(ref_mesh.m_min_level)
-        , m_max_level(ref_mesh.m_max_level)
-        , m_periodic(ref_mesh.m_periodic)
         , m_mpi_neighbourhood(ref_mesh.m_mpi_neighbourhood)
-
+        , m_config(ref_mesh.m_config)
     {
         m_cells[mesh_id_t::cells] = ca;
 
@@ -430,11 +400,8 @@ namespace samurai
     template <class D, class Config>
     inline Mesh_base<D, Config>::Mesh_base(const cl_type& cl, const self_type& ref_mesh)
         : m_domain(ref_mesh.m_domain)
-        , m_min_level(ref_mesh.m_min_level)
-        , m_max_level(ref_mesh.m_max_level)
-        , m_periodic(ref_mesh.m_periodic)
         , m_mpi_neighbourhood(ref_mesh.m_mpi_neighbourhood)
-
+        , m_config(ref_mesh.m_config)
     {
         m_cells[mesh_id_t::cells] = {cl, false};
 
@@ -460,7 +427,7 @@ namespace samurai
         // negartive directions, we actually need 2 times the stencil width inside the hole.
 
         // min_level where the BC can be applied
-        std::size_t min_level_bc = m_min_level;
+        std::size_t min_level_bc = min_level();
         if (scaling_factor <= 0)
         {
             scaling_factor = domain_builder.largest_subdivision();
@@ -468,7 +435,7 @@ namespace samurai
             auto largest_cell_length = samurai::cell_length(scaling_factor, min_level_bc);
             for (const auto& box : domain_builder.removed_boxes())
             {
-                while (box.min_length() < 2 * largest_cell_length * config::max_stencil_width)
+                while (box.min_length() < 2 * largest_cell_length * max_stencil_radius())
                 {
                     scaling_factor /= 2;
                     largest_cell_length /= 2;
@@ -480,10 +447,10 @@ namespace samurai
             auto largest_cell_length = samurai::cell_length(scaling_factor, min_level_bc);
             for (const auto& box : domain_builder.removed_boxes())
             {
-                if (box.min_length() < 2 * largest_cell_length * config::max_stencil_width)
+                if (box.min_length() < 2 * largest_cell_length * max_stencil_radius())
                 {
                     std::cerr << "The hole " << box << " is too small to apply the BC at level " << min_level_bc
-                              << " with the given scaling factor. We need to be able to construct " << (2 * config::max_stencil_width)
+                              << " with the given scaling factor. We need to be able to construct " << (2 * max_stencil_radius())
                               << " ghosts in each direction inside the hole." << std::endl;
                     std::cerr << "Please choose a smaller scaling factor or enlarge the hole." << std::endl;
                     std::exit(1);
@@ -576,25 +543,49 @@ namespace samurai
     template <class D, class Config>
     inline std::size_t Mesh_base<D, Config>::max_level() const
     {
-        return m_max_level;
+        return m_config.max_level();
     }
 
     template <class D, class Config>
     inline std::size_t& Mesh_base<D, Config>::max_level()
     {
-        return m_max_level;
+        return m_config.max_level();
     }
 
     template <class D, class Config>
     inline std::size_t Mesh_base<D, Config>::min_level() const
     {
-        return m_min_level;
+        return m_config.min_level();
     }
 
     template <class D, class Config>
     inline std::size_t& Mesh_base<D, Config>::min_level()
     {
-        return m_min_level;
+        return m_config.min_level();
+    }
+
+    template <class D, class Config>
+    inline std::size_t Mesh_base<D, Config>::graduation_width() const
+    {
+        return m_config.graduation_width();
+    }
+
+    template <class D, class Config>
+    inline int Mesh_base<D, Config>::ghost_width() const
+    {
+        return m_config.ghost_width();
+    }
+
+    template <class D, class Config>
+    inline int Mesh_base<D, Config>::max_stencil_radius() const
+    {
+        return m_config.max_stencil_radius();
+    }
+
+    template <class D, class Config>
+    inline auto& Mesh_base<D, Config>::cfg() const
+    {
+        return m_config;
     }
 
     template <class D, class Config>
@@ -645,6 +636,12 @@ namespace samurai
         return samurai::cell_length(scaling_factor(), level);
     }
 
+    template <class D, class Config>
+    inline double Mesh_base<D, Config>::min_cell_length() const
+    {
+        return cell_length(max_level());
+    }
+
     template <class D, class Config>
     inline auto Mesh_base<D, Config>::domain() const -> const lca_type&
     {
@@ -731,8 +728,8 @@ namespace samurai
     template <class D, class Config>
     inline bool Mesh_base<D, Config>::is_periodic() const
     {
-        return std::any_of(m_periodic.cbegin(),
-                           m_periodic.cend(),
+        return std::any_of(m_config.periodic().cbegin(),
+                           m_config.periodic().cend(),
                            [](bool v)
                            {
                                return v;
@@ -742,13 +739,13 @@ namespace samurai
     template <class D, class Config>
     inline bool Mesh_base<D, Config>::is_periodic(std::size_t d) const
     {
-        return m_periodic[d];
+        return m_config.periodic(d);
     }
 
     template <class D, class Config>
     inline auto Mesh_base<D, Config>::periodicity() const -> const std::array<bool, dim>&
     {
-        return m_periodic;
+        return m_config.periodic();
     }
 
     template <class D, class Config>
@@ -772,8 +769,7 @@ namespace samurai
         swap(m_subdomain, mesh.m_subdomain);
         swap(m_mpi_neighbourhood, mesh.m_mpi_neighbourhood);
         swap(m_union, mesh.m_union);
-        swap(m_max_level, mesh.m_max_level);
-        swap(m_min_level, mesh.m_min_level);
+        swap(m_config, mesh.m_config);
     }
 
     template <class D, class Config>
@@ -941,7 +937,7 @@ namespace samurai
     inline void Mesh_base<D, Config>::construct_domain()
     {
 #ifdef SAMURAI_WITH_MPI
-        lcl_type lcl = {m_max_level};
+        lcl_type lcl = {max_level()};
         mpi::communicator world;
         std::vector<lca_type> all_subdomains(static_cast<std::size_t>(world.size()));
         mpi::all_gather(world, m_subdomain, all_subdomains);
@@ -965,12 +961,12 @@ namespace samurai
     inline void Mesh_base<D, Config>::construct_subdomain()
     {
         // lcl_type lcl = {m_cells[mesh_id_t::cells].max_level()};
-        lcl_type lcl = {m_max_level};
+        lcl_type lcl = {max_level()};
 
         for_each_interval(m_cells[mesh_id_t::cells],
                           [&](std::size_t level, const auto& i, const auto& index)
                           {
-                              std::size_t shift = m_max_level - level;
+                              std::size_t shift = max_level() - level;
                               interval_t to_add = i << shift;
                               auto shift_index  = index << shift;
                               static_nested_loop<dim - 1>(0,
@@ -995,7 +991,7 @@ namespace samurai
     template <class D, class Config>
     inline void Mesh_base<D, Config>::construct_union()
     {
-        std::size_t max_lvl = m_max_level;
+        std::size_t max_lvl = max_level();
 
         // Construction of union cells
         // ===========================
@@ -1046,7 +1042,7 @@ namespace samurai
                 }
                 for (std::size_t d = 0; d < dim; ++d)
                 {
-                    if (m_periodic[d])
+                    if (m_config.periodic(d))
                     {
                         auto shift             = get_periodic_shift(m_domain, m_subdomain.level(), d);
                         auto periodic_set_left = intersection(nestedExpand(m_subdomain, 1), translate(neighbours[i], -shift));
diff --git a/include/samurai/mesh_config.hpp b/include/samurai/mesh_config.hpp
new file mode 100644
index 000000000..7da91bca5
--- /dev/null
+++ b/include/samurai/mesh_config.hpp
@@ -0,0 +1,461 @@
+// Copyright 2018-2025 the samurai's authors
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#pragma once
+
+#include "arguments.hpp"
+#include "cell_array.hpp"
+#include "samurai_config.hpp"
+
+#include <array>
+
+namespace samurai
+{
+
+    template <std::size_t dim_,
+              int prediction_stencil_radius_    = default_config::prediction_stencil_radius,
+              std::size_t max_refinement_level_ = default_config::max_level,
+              class interval_t_                 = default_config::interval_t>
+    class mesh_config
+    {
+      public:
+
+        static constexpr std::size_t dim                  = dim_;
+        static constexpr int prediction_stencil_radius    = prediction_stencil_radius_;
+        static constexpr std::size_t max_refinement_level = max_refinement_level_;
+
+        using interval_t = interval_t_;
+
+        using cell_t   = Cell<dim, interval_t>;
+        using cl_type  = CellList<dim, interval_t, max_refinement_level>;
+        using lcl_type = typename cl_type::lcl_type;
+
+        using ca_type  = CellArray<dim, interval_t, max_refinement_level>;
+        using lca_type = typename ca_type::lca_type;
+
+      private:
+
+        int m_max_stencil_radius       = 1;
+        std::size_t m_graduation_width = default_config::graduation_width;
+        int m_ghost_width              = default_config::ghost_width;
+
+        std::size_t m_min_level   = 0;
+        std::size_t m_max_level   = 6;
+        std::size_t m_start_level = 6;
+
+        double m_approx_box_tol = 0.05;
+        double m_scaling_factor = 0;
+
+        std::array<bool, dim> m_periodic;
+
+        bool m_disable_args_parse          = false;
+        bool m_disable_minimal_ghost_width = false;
+
+      public:
+
+        mesh_config()
+        {
+            m_periodic.fill(false);
+        }
+
+        // m_max_stencil_radius ---------------------------
+
+        /**
+         * @brief set max stencil radius in chained config
+         *
+         * @param stencil_radius
+         * @return auto& returns this object
+         */
+        auto& max_stencil_radius(int stencil_radius)
+        {
+            m_max_stencil_radius = stencil_radius;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on max stencil radius
+         */
+        auto& max_stencil_radius()
+        {
+            return m_max_stencil_radius;
+        }
+
+        /**
+         * @brief get a reference on max stencil radius
+         */
+        const auto& max_stencil_radius() const
+        {
+            return m_max_stencil_radius;
+        }
+
+        /**
+         * @brief set stencil radius from a size (size is twice the size of the radius)
+         *
+         * @param stencil_size
+         * @return auto& returns this object
+         */
+        auto& max_stencil_size(int stencil_size)
+        {
+            m_max_stencil_radius = stencil_size / 2;
+            if (stencil_size % 2 == 1)
+            {
+                m_max_stencil_radius += 1;
+            }
+            return *this;
+        }
+
+        /**
+         * @brief get value of max stencil size
+         */
+        auto max_stencil_size() const
+        {
+            return m_max_stencil_radius * 2;
+        }
+
+        // m_graduation_width -----------------------------
+
+        /**
+         * @brief set graduation width in chained config
+         *
+         * @param grad_width
+         * @return auto& returns this object
+         */
+        auto& graduation_width(std::size_t grad_width)
+        {
+            m_graduation_width = grad_width;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on graduation width
+         */
+        auto& graduation_width()
+        {
+            return m_graduation_width;
+        }
+
+        /**
+         * @brief get a reference on graduation width
+         */
+        const auto& graduation_width() const
+        {
+            return m_graduation_width;
+        }
+
+        // m_ghost_width ----------------------------------
+
+        /**
+         * @brief get a reference on ghost width
+         */
+        const auto& ghost_width() const
+        {
+            return m_ghost_width;
+        }
+
+        // m_min_level ------------------------------------
+
+        /**
+         * @brief set min level in chained config
+         *
+         * @param level
+         * @return auto& returns this object
+         */
+        auto& min_level(std::size_t level)
+        {
+            m_min_level = level;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on min level
+         */
+        auto& min_level()
+        {
+            return m_min_level;
+        }
+
+        /**
+         * @brief get a reference on min level
+         */
+        const auto& min_level() const
+        {
+            return m_min_level;
+        }
+
+        // m_max_level ------------------------------------
+
+        /**
+         * @brief set max level in chained config
+         *
+         * @param level
+         * @return auto& returns this object
+         */
+        auto& max_level(std::size_t level)
+        {
+            m_max_level = level;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on max level
+         */
+        auto& max_level()
+        {
+            return m_max_level;
+        }
+
+        /**
+         * @brief get a reference on max level
+         */
+        const auto& max_level() const
+        {
+            return m_max_level;
+        }
+
+        // m_start_level ------------------------------------
+
+        /**
+         * @brief set start level in chained config
+         *
+         * @param level
+         * @return auto& returns this object
+         */
+        auto& start_level(std::size_t level)
+        {
+            m_start_level = level;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on start level
+         */
+        auto& start_level()
+        {
+            return m_start_level;
+        }
+
+        /**
+         * @brief get a reference on start level
+         */
+        const auto& start_level() const
+        {
+            return m_start_level;
+        }
+
+        // m_approx_box_tol -------------------------------
+
+        /**
+         * @brief set approximation box tolerance in chained config
+         *
+         * @param tol
+         * @return auto& returns this object
+         */
+        auto& approx_box_tol(double tol)
+        {
+            m_approx_box_tol = tol;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on approximation box tolerance
+         */
+        auto& approx_box_tol()
+        {
+            return m_approx_box_tol;
+        }
+
+        /**
+         * @brief get a reference on approximation box tolerance
+         */
+        const auto& approx_box_tol() const
+        {
+            return m_approx_box_tol;
+        }
+
+        // m_scaling_factor -------------------------------
+
+        /**
+         * @brief set scaling factor in chained config
+         *
+         * @param factor
+         * @return auto& returns this object
+         */
+        auto& scaling_factor(double factor)
+        {
+            m_scaling_factor = factor;
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on scaling factor
+         */
+        auto& scaling_factor()
+        {
+            return m_scaling_factor;
+        }
+
+        /**
+         * @brief get a reference on scaling factor
+         */
+        const auto& scaling_factor() const
+        {
+            return m_scaling_factor;
+        }
+
+        // m_periodic -------------------------------------
+
+        /**
+         * @brief set periodicity in chained config
+         *
+         * @param periodicity
+         * @return auto& returns this object
+         */
+        auto& periodic(std::array<bool, dim> const& periodicity)
+        {
+            m_periodic = periodicity;
+            return *this;
+        }
+
+        /**
+         * @brief set periodicity in chained config with a value to fill in each direction
+         *
+         * @param periodicity
+         * @return auto& returns this object
+         */
+        auto& periodic(bool periodicity)
+        {
+            m_periodic.fill(periodicity);
+            return *this;
+        }
+
+        /**
+         * @brief get a reference on periodicity array
+         */
+        auto& periodic()
+        {
+            return m_periodic;
+        }
+
+        /**
+         * @brief get a reference on periodicity array
+         */
+        const auto& periodic() const
+        {
+            return m_periodic;
+        }
+
+        /**
+         * @brief get a reference on periodicity in direction i
+         */
+        auto& periodic(std::size_t i)
+        {
+            return m_periodic[i];
+        }
+
+        /**
+         * @brief get a reference on periodicity in direction i
+         */
+        const auto& periodic(std::size_t i) const
+        {
+            return m_periodic[i];
+        }
+
+        // m_disable_args_parse ---------------------------
+
+        /**
+         * @brief disable argument parse
+         */
+        auto& disable_args_parse()
+        {
+            m_disable_args_parse = true;
+            return *this;
+        }
+
+        // disable_minimal_ghost_width --------------------
+        auto& disable_minimal_ghost_width()
+        {
+            m_disable_minimal_ghost_width = true;
+            return *this;
+        }
+
+        /**
+         * @brief parse arguments and set value to default samurai config value if needed
+         */
+        void parse_args()
+        {
+            if (!m_disable_args_parse)
+            {
+                if (args::max_stencil_radius != std::numeric_limits<int>::max())
+                {
+                    m_max_stencil_radius = args::max_stencil_radius;
+                }
+                if (args::graduation_width != std::numeric_limits<std::size_t>::max())
+                {
+                    m_graduation_width = args::graduation_width;
+                }
+                if (args::min_level != std::numeric_limits<std::size_t>::max())
+                {
+                    m_min_level = args::min_level;
+                }
+                if (args::max_level != std::numeric_limits<std::size_t>::max())
+                {
+                    m_max_level = args::max_level;
+                }
+                if (args::start_level != std::numeric_limits<std::size_t>::max())
+                {
+                    m_start_level = args::start_level;
+                }
+                // if (args::approx_box_tol != std::numeric_limits<double>::infinity())
+                // {
+                //     m_approx_box_tol = args::approx_box_tol;
+                // }
+                // if (args::scaling_factor != std::numeric_limits<double>::infinity())
+                // {
+                //     m_scaling_factor = args::scaling_factor;
+                // }
+                if (m_max_level < m_min_level)
+                {
+                    std::cerr << "Max level must be greater than min level." << std::endl;
+                    exit(EXIT_FAILURE);
+                }
+            }
+
+            m_ghost_width = std::max(m_max_stencil_radius, static_cast<int>(prediction_stencil_radius));
+
+            if (!m_disable_minimal_ghost_width)
+            {
+                // 2 is because prediction_stencil_radius=1, if >1 we don't know what to do...
+                // The idea is to have enough ghosts at the boundary for the reconstruction and the transfer to work.
+                m_max_stencil_radius = std::max(m_max_stencil_radius, 2);
+            }
+        }
+
+      private:
+
+#ifdef SAMURAI_WITH_MPI
+        friend class boost::serialization::access;
+
+        template <class Archive>
+        void serialize(Archive& ar, const unsigned long)
+        {
+            // ar & m_max_stencil_radius;
+            // ar & m_graduation_width;
+            // ar & m_ghost_width;
+            ar & m_min_level;
+            ar & m_max_level;
+            // ar & m_approx_box_tol;
+            // ar & m_scaling_factor;
+            // ar & m_disable_args_parse;
+        }
+#endif
+    };
+
+    template <class mesh_cfg_t, class mesh_id_t_>
+    class complete_mesh_config
+        : public mesh_config<mesh_cfg_t::dim, mesh_cfg_t::prediction_stencil_radius, mesh_cfg_t::max_refinement_level, typename mesh_cfg_t::interval_t>
+    {
+      public:
+
+        using mesh_id_t = mesh_id_t_;
+    };
+} // namespace samurai
diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp
index 958c14a7f..b8af6fd37 100644
--- a/include/samurai/mr/adapt.hpp
+++ b/include/samurai/mr/adapt.hpp
@@ -220,7 +220,7 @@ namespace samurai
     {
         // Since the adaptation process starts at max_level, we just need to flag to `keep` the boundary cells at max_level only.
         // There will never be boundary cells at lower levels.
-        auto bdry = domain_boundary_layer(mesh, mesh.max_level(), direction, Mesh::config::max_stencil_width);
+        auto bdry = domain_boundary_layer(mesh, mesh.max_level(), direction, static_cast<std::size_t>(mesh.max_stencil_radius()));
         for_each_cell(mesh,
                       bdry,
                       [&](auto& cell)
@@ -380,12 +380,7 @@ namespace samurai
         using ca_type = typename mesh_t::ca_type;
 
         ca_type new_ca = update_cell_array_from_tag(mesh[mesh_id_t::cells], m_tag);
-        make_graduation(new_ca,
-                        mesh.domain(),
-                        mesh.mpi_neighbourhood(),
-                        mesh.periodicity(),
-                        mesh_t::config::graduation_width,
-                        mesh_t::config::max_stencil_width);
+        make_graduation(new_ca, mesh.domain(), mesh.mpi_neighbourhood(), mesh.periodicity(), mesh.graduation_width(), mesh.max_stencil_radius());
         mesh_t new_mesh{new_ca, mesh};
 #ifdef SAMURAI_WITH_MPI
         mpi::communicator world;
diff --git a/include/samurai/mr/mesh.hpp b/include/samurai/mr/mesh.hpp
index 10912c566..6556974c7 100644
--- a/include/samurai/mr/mesh.hpp
+++ b/include/samurai/mr/mesh.hpp
@@ -30,29 +30,6 @@ namespace samurai
         reference        = all_cells
     };
 
-    template <std::size_t dim_,
-              std::size_t max_stencil_width_    = default_config::ghost_width,
-              std::size_t graduation_width_     = default_config::graduation_width,
-              std::size_t prediction_order_     = default_config::prediction_order,
-              std::size_t max_refinement_level_ = default_config::max_level,
-              class TInterval                   = default_config::interval_t>
-    struct MRConfig
-    {
-        static constexpr std::size_t dim                  = dim_;
-        static constexpr std::size_t max_refinement_level = max_refinement_level_;
-        static constexpr int max_stencil_width            = max_stencil_width_;
-        static constexpr std::size_t graduation_width     = graduation_width_;
-        static constexpr int prediction_order             = prediction_order_;
-
-        // static constexpr int ghost_width = std::max(std::max(2 *
-        // static_cast<int>(graduation_width) - 1,
-        //                                                      static_cast<int>(max_stencil_width)),
-        //                                             static_cast<int>(prediction_order));
-        static constexpr int ghost_width = std::max(static_cast<int>(max_stencil_width), static_cast<int>(prediction_order));
-        using interval_t                 = TInterval;
-        using mesh_id_t                  = MRMeshId;
-    };
-
     template <class Config>
     class MRMesh : public samurai::Mesh_base<MRMesh<Config>, Config>
     {
@@ -72,11 +49,18 @@ namespace samurai
         using ca_type  = typename base_type::ca_type;
         using lca_type = typename base_type::lca_type;
 
+        using base_type::ghost_width;
+        using base_type::max_stencil_radius;
+
         MRMesh() = default;
         MRMesh(const ca_type& ca, const self_type& ref_mesh);
         MRMesh(const cl_type& cl, const self_type& ref_mesh);
-        MRMesh(const cl_type& cl, std::size_t min_level, std::size_t max_level);
-        MRMesh(const ca_type& ca, std::size_t min_level, std::size_t max_level);
+        MRMesh(const mesh_config<Config::dim>& config, const cl_type& cl);
+        MRMesh(const mesh_config<Config::dim>& config, const ca_type& ca);
+        MRMesh(const mesh_config<Config::dim>& config, const samurai::Box<double, dim>& b);
+        MRMesh(const mesh_config<Config::dim>& config, const samurai::DomainBuilder<dim>& domain_builder);
+
+        // deprecated constructors
         MRMesh(const samurai::Box<double, dim>& b,
                std::size_t min_level,
                std::size_t max_level,
@@ -113,14 +97,26 @@ namespace samurai
     }
 
     template <class Config>
-    inline MRMesh<Config>::MRMesh(const cl_type& cl, std::size_t min_level, std::size_t max_level)
-        : base_type(cl, min_level, max_level)
+    inline MRMesh<Config>::MRMesh(const mesh_config<Config::dim>& config, const cl_type& cl)
+        : base_type(config, cl)
     {
     }
 
     template <class Config>
-    inline MRMesh<Config>::MRMesh(const ca_type& ca, std::size_t min_level, std::size_t max_level)
-        : base_type(ca, min_level, max_level)
+    inline MRMesh<Config>::MRMesh(const mesh_config<Config::dim>& config, const ca_type& ca)
+        : base_type(config, ca)
+    {
+    }
+
+    template <class Config>
+    inline MRMesh<Config>::MRMesh(const mesh_config<Config::dim>& config, const samurai::Box<double, dim>& b)
+        : base_type(config, b)
+    {
+    }
+
+    template <class Config>
+    inline MRMesh<Config>::MRMesh(const mesh_config<Config::dim>& config, const samurai::DomainBuilder<dim>& domain_builder)
+        : base_type(config, domain_builder)
     {
     }
 
@@ -186,18 +182,19 @@ namespace samurai
         //
         // level 0 |.......|-------|.......|       |.......|-------|.......|
         //
-        for_each_interval(
-            this->cells()[mesh_id_t::cells],
-            [&](std::size_t level, const auto& interval, const auto& index_yz)
-            {
-                lcl_type& lcl = cell_list[level];
-                static_nested_loop<dim - 1, -config::max_stencil_width, config::max_stencil_width + 1>(
-                    [&](auto stencil)
-                    {
-                        auto index = xt::eval(index_yz + stencil);
-                        lcl[index].add_interval({interval.start - config::max_stencil_width, interval.end + config::max_stencil_width});
-                    });
-            });
+        for_each_interval(this->cells()[mesh_id_t::cells],
+                          [&](std::size_t level, const auto& interval, const auto& index_yz)
+                          {
+                              lcl_type& lcl = cell_list[level];
+                              static_nested_loop<dim - 1>(
+                                  -max_stencil_radius(),
+                                  max_stencil_radius() + 1,
+                                  [&](auto stencil)
+                                  {
+                                      auto index = xt::eval(index_yz + stencil);
+                                      lcl[index].add_interval({interval.start - max_stencil_radius(), interval.end + max_stencil_radius()});
+                                  });
+                          });
         this->cells()[mesh_id_t::cells_and_ghosts] = {cell_list, false};
 
         // Add cells for the MRA
@@ -213,12 +210,12 @@ namespace samurai
                     {
                         lcl_type& lcl = cell_list[level - 1];
 
-                        static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                        static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>(
                             [&](auto stencil)
                             {
                                 auto new_interval = interval >> 1;
-                                lcl[(index_yz >> 1) + stencil].add_interval(
-                                    {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
+                                lcl[(index_yz >> 1) + stencil].add_interval({new_interval.start - config::prediction_stencil_radius,
+                                                                             new_interval.end + config::prediction_stencil_radius});
                             });
                     });
 
@@ -231,12 +228,12 @@ namespace samurai
                         {
                             lcl_type& lcl = cell_list[level - 2];
 
-                            static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                            static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>(
                                 [&](auto stencil)
                                 {
                                     auto new_interval = interval >> 2;
-                                    lcl[(index_yz >> 2) + stencil].add_interval(
-                                        {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
+                                    lcl[(index_yz >> 2) + stencil].add_interval({new_interval.start - config::prediction_stencil_radius,
+                                                                                 new_interval.end + config::prediction_stencil_radius});
                                 });
                         }
                     });
@@ -261,12 +258,12 @@ namespace samurai
                             {
                                 lcl_type& lclm1 = cell_list[level - 1];
 
-                                static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                                static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>(
                                     [&](auto stencil)
                                     {
                                         auto new_interval = interval >> 1;
-                                        lclm1[(index_yz >> 1) + stencil].add_interval(
-                                            {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
+                                        lclm1[(index_yz >> 1) + stencil].add_interval({new_interval.start - config::prediction_stencil_radius,
+                                                                                       new_interval.end + config::prediction_stencil_radius});
                                     });
                             }
                         });
@@ -313,8 +310,8 @@ namespace samurai
                         box_t box;
                         for (std::size_t d = 0; d < dim; ++d)
                         {
-                            box.min_corner()[d] = (neighbor_min_indices[d] >> delta_l) - config::ghost_width;
-                            box.max_corner()[d] = (neighbor_max_indices[d] >> delta_l) + config::ghost_width;
+                            box.min_corner()[d] = (neighbor_min_indices[d] >> delta_l) - ghost_width();
+                            box.max_corner()[d] = (neighbor_max_indices[d] >> delta_l) + ghost_width();
                         }
                         neighbourhood_extended_subdomain[neighbor_id][level] = {level, box};
                     }
@@ -329,8 +326,8 @@ namespace samurai
 
                 for (std::size_t d = 0; d < dim; ++d)
                 {
-                    min_corner[d] = (subdomain_min_indices[d] >> delta_l) - config::ghost_width;
-                    max_corner[d] = (subdomain_max_indices[d] >> delta_l) + config::ghost_width;
+                    min_corner[d] = (subdomain_min_indices[d] >> delta_l) - ghost_width();
+                    max_corner[d] = (subdomain_max_indices[d] >> delta_l) + ghost_width();
                 }
                 lca_type lca_extended_subdomain(level, box_t(min_corner, max_corner));
                 for (std::size_t d = 0; d < dim; ++d)
@@ -339,13 +336,13 @@ namespace samurai
                     {
                         shift[d] = (domain_max_indices[d] - domain_min_indices[d]) >> delta_l;
 
-                        min_corner[d] = (domain_min_indices[d] >> delta_l) - config::ghost_width;
-                        max_corner[d] = (domain_min_indices[d] >> delta_l) + config::ghost_width;
+                        min_corner[d] = (domain_min_indices[d] >> delta_l) - ghost_width();
+                        max_corner[d] = (domain_min_indices[d] >> delta_l) + ghost_width();
 
                         lca_type lca_min(level, box_t(min_corner, max_corner));
 
-                        min_corner[d] = (domain_max_indices[d] >> delta_l) - config::ghost_width;
-                        max_corner[d] = (domain_max_indices[d] >> delta_l) + config::ghost_width;
+                        min_corner[d] = (domain_max_indices[d] >> delta_l) - ghost_width();
+                        max_corner[d] = (domain_max_indices[d] >> delta_l) + ghost_width();
 
                         lca_type lca_max(level, box_t(min_corner, max_corner));
 
@@ -390,8 +387,8 @@ namespace samurai
 
                         /* reset variables for next iterations. */
                         shift[d]      = 0;
-                        min_corner[d] = (subdomain_min_indices[d] >> delta_l) - config::ghost_width;
-                        max_corner[d] = (subdomain_max_indices[d] >> delta_l) + config::ghost_width;
+                        min_corner[d] = (subdomain_min_indices[d] >> delta_l) - ghost_width();
+                        max_corner[d] = (subdomain_max_indices[d] >> delta_l) + ghost_width();
                     }
                 }
             }
@@ -421,7 +418,7 @@ namespace samurai
             {
                 for (std::size_t level = 0; level <= this->max_level(); ++level)
                 {
-                    auto expr = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width),
+                    auto expr = intersection(nestedExpand(self(this->subdomain()).on(level), ghost_width()),
                                              neighbour.mesh[mesh_id_t::reference][level]);
                     expr(
                         [&](const auto& interval, const auto& index_yz)
@@ -434,7 +431,7 @@ namespace samurai
                         if (this->is_periodic(d))
                         {
                             auto domain_shift = get_periodic_shift(this->domain(), level, d);
-                            auto expr_left    = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width),
+                            auto expr_left    = intersection(nestedExpand(self(this->subdomain()).on(level), ghost_width()),
                                                           translate(neighbour.mesh[mesh_id_t::reference][level], -domain_shift));
                             expr_left(
                                 [&](const auto& interval, const auto& index_yz)
@@ -443,7 +440,7 @@ namespace samurai
                                     lcl[index_yz].add_interval(interval);
                                 });
 
-                            auto expr_right = intersection(nestedExpand(self(this->subdomain()).on(level), config::ghost_width),
+                            auto expr_right = intersection(nestedExpand(self(this->subdomain()).on(level), ghost_width()),
                                                            translate(neighbour.mesh[mesh_id_t::reference][level], domain_shift));
                             expr_right(
                                 [&](const auto& interval, const auto& index_yz)
@@ -493,6 +490,57 @@ namespace samurai
             return out;
         }
     }
+
+    // create an empty mesh
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, MRMeshId>>
+    auto make_empty_MRMesh(const mesh_config_t&)
+    {
+        return MRMesh<complete_mesh_config_t>();
+    }
+
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, MRMeshId>>
+    auto make_MRMesh(const mesh_config_t& cfg, const typename MRMesh<complete_mesh_config_t>::cl_type& cl)
+    {
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+        mesh_cfg.start_level() = mesh_cfg.max_level(); // cppcheck-suppress unreadVariable
+
+        return MRMesh<complete_mesh_config_t>(cfg, cl);
+    }
+
+    template <class mesh_config_t, class complete_mesh_config_t = complete_mesh_config<mesh_config_t, MRMeshId>>
+    auto make_MRMesh(const mesh_config_t& cfg, const typename MRMesh<complete_mesh_config_t>::ca_type& ca)
+    {
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+        mesh_cfg.start_level() = mesh_cfg.max_level(); // cppcheck-suppress unreadVariable
+
+        return MRMesh<complete_mesh_config_t>(cfg, ca);
+    }
+
+    template <class mesh_config_t>
+    auto make_MRMesh(const mesh_config_t& cfg, const samurai::Box<double, mesh_config_t::dim>& b)
+    {
+        using complete_cfg_t = complete_mesh_config<mesh_config_t, MRMeshId>;
+
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+        mesh_cfg.start_level() = mesh_cfg.max_level(); // cppcheck-suppress unreadVariable
+
+        return MRMesh<complete_cfg_t>(mesh_cfg, b);
+    }
+
+    template <class mesh_config_t>
+    auto make_MRMesh(const mesh_config_t& cfg, const samurai::DomainBuilder<mesh_config_t::dim>& domain_builder)
+    {
+        using complete_cfg_t = complete_mesh_config<mesh_config_t, MRMeshId>;
+
+        auto mesh_cfg = cfg;
+        mesh_cfg.parse_args();
+        mesh_cfg.start_level() = mesh_cfg.max_level(); // cppcheck-suppress unreadVariable
+
+        return MRMesh<complete_cfg_t>(mesh_cfg, domain_builder);
+    }
 } // namespace samurai
 
 template <>
diff --git a/include/samurai/mr/mesh_with_overleaves.hpp b/include/samurai/mr/mesh_with_overleaves.hpp
index f94ee3e15..dd7de08fd 100644
--- a/include/samurai/mr/mesh_with_overleaves.hpp
+++ b/include/samurai/mr/mesh_with_overleaves.hpp
@@ -28,21 +28,21 @@ namespace samurai
     };
 
     template <std::size_t dim_,
-              std::size_t max_stencil_width_    = default_config::ghost_width,
-              std::size_t graduation_width_     = default_config::graduation_width,
-              std::size_t max_refinement_level_ = default_config::max_level,
-              std::size_t prediction_order_     = default_config::prediction_order,
-              class TInterval                   = default_config::interval_t>
+              std::size_t max_stencil_width_         = default_config::ghost_width,
+              std::size_t graduation_width_          = default_config::graduation_width,
+              std::size_t max_refinement_level_      = default_config::max_level,
+              std::size_t prediction_stencil_radius_ = default_config::prediction_stencil_radius,
+              class TInterval                        = default_config::interval_t>
     struct MROConfig
     {
-        static constexpr std::size_t dim                  = dim_;
-        static constexpr std::size_t max_refinement_level = max_refinement_level_;
-        static constexpr std::size_t max_stencil_width    = max_stencil_width_;
-        static constexpr std::size_t graduation_width     = graduation_width_;
-        static constexpr std::size_t prediction_order     = prediction_order_;
+        static constexpr std::size_t dim                       = dim_;
+        static constexpr std::size_t max_refinement_level      = max_refinement_level_;
+        static constexpr std::size_t max_stencil_width         = max_stencil_width_;
+        static constexpr std::size_t graduation_width          = graduation_width_;
+        static constexpr std::size_t prediction_stencil_radius = prediction_stencil_radius_;
 
         static constexpr int ghost_width = std::max(std::max(2 * static_cast<int>(graduation_width) - 1, static_cast<int>(max_stencil_width)),
-                                                    static_cast<int>(prediction_order));
+                                                    static_cast<int>(prediction_stencil_radius));
         using interval_t = TInterval;
         using mesh_id_t  = MROMeshId;
     };
diff --git a/include/samurai/mr/operators.hpp b/include/samurai/mr/operators.hpp
index 256ddedda..01f8ecdf8 100644
--- a/include/samurai/mr/operators.hpp
+++ b/include/samurai/mr/operators.hpp
@@ -248,7 +248,7 @@ namespace samurai
 
         INIT_OPERATOR(compute_detail_op)
 
-        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_order>
+        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_stencil_radius>
         inline void operator()(Dim<1>, T1& detail, const T2& field) const
         {
             if constexpr (order == 0)
@@ -265,7 +265,7 @@ namespace samurai
             }
         }
 
-        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_order>
+        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_stencil_radius>
         inline void operator()(Dim<2>, T1& detail, const T2& field) const
         {
             if constexpr (order == 0)
@@ -311,7 +311,7 @@ namespace samurai
             }
         }
 
-        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_order>
+        template <class T1, class T2, std::size_t order = T2::mesh_t::config::prediction_stencil_radius>
         inline void operator()(Dim<3>, T1& detail, const T2& field) const
         {
             if constexpr (order == 0)
diff --git a/include/samurai/petsc/fv/FV_scheme_assembly.hpp b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
index abe9ab2b9..f68d1ad87 100644
--- a/include/samurai/petsc/fv/FV_scheme_assembly.hpp
+++ b/include/samurai/petsc/fv/FV_scheme_assembly.hpp
@@ -28,25 +28,25 @@ namespace samurai
 
           public:
 
-            using cfg_t                                           = typename Scheme::cfg_t;
-            using bdry_cfg_t                                      = typename Scheme::bdry_cfg;
-            using field_t                                         = typename Scheme::field_t;
-            using output_field_t                                  = typename cfg_t::output_field_t;
-            using mesh_t                                          = typename field_t::mesh_t;
-            using mesh_id_t                                       = typename mesh_t::mesh_id_t;
-            using interval_t                                      = typename mesh_t::interval_t;
-            using mesh_interval_t                                 = typename mesh_t::mesh_interval_t;
-            using field_value_type                                = typename field_t::value_type; // double
-            using coord_index_t                                   = typename interval_t::coord_index_t;
-            using index_t                                         = typename interval_t::index_t;
-            static constexpr std::size_t dim                      = field_t::dim;
-            static constexpr std::size_t input_n_comp             = field_t::n_comp;
-            static constexpr std::size_t output_n_comp            = output_field_t::n_comp;
-            static constexpr std::size_t prediction_order         = mesh_t::config::prediction_order;
-            static constexpr std::size_t bdry_neighbourhood_width = bdry_cfg_t::neighbourhood_width;
-            static constexpr std::size_t bdry_stencil_size        = bdry_cfg_t::stencil_size;
-            static constexpr std::size_t nb_bdry_ghosts           = bdry_cfg_t::nb_ghosts;
-            using cell_t                                          = Cell<dim, interval_t>;
+            using cfg_t                                            = typename Scheme::cfg_t;
+            using bdry_cfg_t                                       = typename Scheme::bdry_cfg;
+            using field_t                                          = typename Scheme::field_t;
+            using output_field_t                                   = typename cfg_t::output_field_t;
+            using mesh_t                                           = typename field_t::mesh_t;
+            using mesh_id_t                                        = typename mesh_t::mesh_id_t;
+            using interval_t                                       = typename mesh_t::interval_t;
+            using mesh_interval_t                                  = typename mesh_t::mesh_interval_t;
+            using field_value_type                                 = typename field_t::value_type; // double
+            using coord_index_t                                    = typename interval_t::coord_index_t;
+            using index_t                                          = typename interval_t::index_t;
+            static constexpr std::size_t dim                       = field_t::dim;
+            static constexpr std::size_t input_n_comp              = field_t::n_comp;
+            static constexpr std::size_t output_n_comp             = output_field_t::n_comp;
+            static constexpr std::size_t prediction_stencil_radius = mesh_t::config::prediction_stencil_radius;
+            static constexpr std::size_t bdry_neighbourhood_width  = bdry_cfg_t::neighbourhood_width;
+            static constexpr std::size_t bdry_stencil_size         = bdry_cfg_t::stencil_size;
+            static constexpr std::size_t nb_bdry_ghosts            = bdry_cfg_t::nb_ghosts;
+            using cell_t                                           = Cell<dim, interval_t>;
 
             using dirichlet_t = DirichletImpl<nb_bdry_ghosts, field_t>;
             using neumann_t   = NeumannImpl<nb_bdry_ghosts, field_t>;
@@ -695,7 +695,7 @@ namespace samurai
                                 //                                                  1 + 9=10 in 2D
                                 // Order 2: cell + hypercube of 5 coarser cells --> 1 + 5= 6 in 1D
                                 //                                                  1 +25=21 in 2D
-                                static constexpr std::size_t pred_stencil_size = 1 + ce_pow(2 * prediction_order + 1, dim);
+                                static constexpr std::size_t pred_stencil_size = 1 + ce_pow(2 * prediction_stencil_radius + 1, dim);
 
                                 nnz[static_cast<std::size_t>(row_index(ghost, field_i))] = pred_stencil_size;
                             }
@@ -828,19 +828,20 @@ namespace samurai
                             auto ig      = ii >> 1;
                             double isign = (ii & 1) ? -1 : 1;
 
-                            auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
+                            auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
 
                             auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig)), field_i);
                             MatSetValue(A, ghost_index, parent_index, -scaling, current_insert_mode());
 
                             for (std::size_t ci = 0; ci < interpx.size(); ++ci)
                             {
-                                if (ci != prediction_order)
+                                if (ci != prediction_stencil_radius)
                                 {
                                     double value           = -interpx[ci];
                                     auto coarse_cell_index = this->col_index(
                                         static_cast<PetscInt>(
-                                            this->mesh().get_index(ghost.level - 1, ig + static_cast<coord_index_t>(ci - prediction_order))),
+                                            this->mesh().get_index(ghost.level - 1,
+                                                                   ig + static_cast<coord_index_t>(ci - prediction_stencil_radius))),
                                         field_i);
                                     MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
                                 }
@@ -859,17 +860,17 @@ namespace samurai
                 auto ig      = ii >> 1;
                 double isign = (ii & 1) ? -1 : 1;
 
-                auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
+                auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
 
                 auto parent_index = this->mesh().get_index(ghost.level - 1, ig);
                 linear_comb.emplace_back(parent_index, 1.);
 
                 for (std::size_t ci = 0; ci < interpx.size(); ++ci)
                 {
-                    if (ci != prediction_order)
+                    if (ci != prediction_stencil_radius)
                     {
                         auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
-                                                                        ig + static_cast<coord_index_t>(ci - prediction_order));
+                                                                        ig + static_cast<coord_index_t>(ci - prediction_stencil_radius));
                         linear_comb.emplace_back(coarse_cell_index, interpx[ci]);
                     }
                 }
@@ -896,8 +897,8 @@ namespace samurai
                             double isign = (ii & 1) ? -1 : 1;
                             double jsign = (j & 1) ? -1 : 1;
 
-                            auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
-                            auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign);
+                            auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
+                            auto interpy = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(jsign);
 
                             auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig, jg)),
                                                                 field_i);
@@ -907,14 +908,15 @@ namespace samurai
                             {
                                 for (std::size_t cj = 0; cj < interpy.size(); ++cj)
                                 {
-                                    if (ci != prediction_order || cj != prediction_order)
+                                    if (ci != prediction_stencil_radius || cj != prediction_stencil_radius)
                                     {
                                         double value           = -interpx[ci] * interpy[cj];
-                                        auto coarse_cell_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(
-                                                                                     ghost.level - 1,
-                                                                                     ig + static_cast<coord_index_t>(ci - prediction_order),
-                                                                                     jg + static_cast<coord_index_t>(cj - prediction_order))),
-                                                                                 field_i);
+                                        auto coarse_cell_index = this->col_index(
+                                            static_cast<PetscInt>(
+                                                this->mesh().get_index(ghost.level - 1,
+                                                                       ig + static_cast<coord_index_t>(ci - prediction_stencil_radius),
+                                                                       jg + static_cast<coord_index_t>(cj - prediction_stencil_radius))),
+                                            field_i);
                                         MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
                                     }
                                 }
@@ -936,8 +938,8 @@ namespace samurai
                 double isign = (ii & 1) ? -1 : 1;
                 double jsign = (j & 1) ? -1 : 1;
 
-                auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
-                auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign);
+                auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
+                auto interpy = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(jsign);
 
                 auto parent_index = this->mesh().get_index(ghost.level - 1, ig, jg);
                 linear_comb.emplace_back(parent_index, 1.);
@@ -946,11 +948,11 @@ namespace samurai
                 {
                     for (std::size_t cj = 0; cj < interpy.size(); ++cj)
                     {
-                        if (ci != prediction_order || cj != prediction_order)
+                        if (ci != prediction_stencil_radius || cj != prediction_stencil_radius)
                         {
                             auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
-                                                                            ig + static_cast<coord_index_t>(ci - prediction_order),
-                                                                            jg + static_cast<coord_index_t>(cj - prediction_order));
+                                                                            ig + static_cast<coord_index_t>(ci - prediction_stencil_radius),
+                                                                            jg + static_cast<coord_index_t>(cj - prediction_stencil_radius));
                             linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj]);
                         }
                     }
@@ -981,9 +983,9 @@ namespace samurai
                             double jsign = (j & 1) ? -1 : 1;
                             double ksign = (k & 1) ? -1 : 1;
 
-                            auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
-                            auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign);
-                            auto interpz = samurai::interp_coeffs<2 * prediction_order + 1>(ksign);
+                            auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
+                            auto interpy = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(jsign);
+                            auto interpz = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(ksign);
 
                             auto parent_index = this->col_index(static_cast<PetscInt>(this->mesh().get_index(ghost.level - 1, ig, jg, kg)),
                                                                 field_i);
@@ -995,15 +997,16 @@ namespace samurai
                                 {
                                     for (std::size_t ck = 0; ck < interpz.size(); ++ck)
                                     {
-                                        if (ci != prediction_order || cj != prediction_order || ck != prediction_order)
+                                        if (ci != prediction_stencil_radius || cj != prediction_stencil_radius
+                                            || ck != prediction_stencil_radius)
                                         {
                                             double value           = -interpx[ci] * interpy[cj] * interpz[ck];
                                             auto coarse_cell_index = this->col_index(
-                                                static_cast<PetscInt>(
-                                                    this->mesh().get_index(ghost.level - 1,
-                                                                           ig + static_cast<coord_index_t>(ci - prediction_order),
-                                                                           jg + static_cast<coord_index_t>(cj - prediction_order),
-                                                                           kg + static_cast<coord_index_t>(ck - prediction_order))),
+                                                static_cast<PetscInt>(this->mesh().get_index(
+                                                    ghost.level - 1,
+                                                    ig + static_cast<coord_index_t>(ci - prediction_stencil_radius),
+                                                    jg + static_cast<coord_index_t>(cj - prediction_stencil_radius),
+                                                    kg + static_cast<coord_index_t>(ck - prediction_stencil_radius))),
                                                 field_i);
                                             MatSetValue(A, ghost_index, coarse_cell_index, scaling * value, current_insert_mode());
                                         }
@@ -1030,9 +1033,9 @@ namespace samurai
                 double jsign = (j & 1) ? -1 : 1;
                 double ksign = (k & 1) ? -1 : 1;
 
-                auto interpx = samurai::interp_coeffs<2 * prediction_order + 1>(isign);
-                auto interpy = samurai::interp_coeffs<2 * prediction_order + 1>(jsign);
-                auto interpz = samurai::interp_coeffs<2 * prediction_order + 1>(ksign);
+                auto interpx = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(isign);
+                auto interpy = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(jsign);
+                auto interpz = samurai::interp_coeffs<2 * prediction_stencil_radius + 1>(ksign);
 
                 auto parent_index = this->mesh().get_index(ghost.level - 1, ig, jg, kg);
                 linear_comb.emplace_back(parent_index, 1.);
@@ -1043,12 +1046,13 @@ namespace samurai
                     {
                         for (std::size_t ck = 0; ck < interpz.size(); ++ck)
                         {
-                            if (ci != prediction_order || cj != prediction_order || ck != prediction_order)
+                            if (ci != prediction_stencil_radius || cj != prediction_stencil_radius || ck != prediction_stencil_radius)
                             {
-                                auto coarse_cell_index = this->mesh().get_index(ghost.level - 1,
-                                                                                ig + static_cast<coord_index_t>(ci - prediction_order),
-                                                                                jg + static_cast<coord_index_t>(cj - prediction_order),
-                                                                                kg + static_cast<coord_index_t>(ck - prediction_order));
+                                auto coarse_cell_index = this->mesh().get_index(
+                                    ghost.level - 1,
+                                    ig + static_cast<coord_index_t>(ci - prediction_stencil_radius),
+                                    jg + static_cast<coord_index_t>(cj - prediction_stencil_radius),
+                                    kg + static_cast<coord_index_t>(ck - prediction_stencil_radius));
                                 linear_comb.emplace_back(coarse_cell_index, interpx[ci] * interpy[cj] * interpz[ck]);
                             }
                         }
diff --git a/include/samurai/reconstruction.hpp b/include/samurai/reconstruction.hpp
index f692c1d21..00af7d166 100644
--- a/include/samurai/reconstruction.hpp
+++ b/include/samurai/reconstruction.hpp
@@ -373,7 +373,7 @@ namespace samurai
         // inline void operator()(Dim<d>, std::size_t& reconstruct_level, T1& dest, const T2& src) const
         // {
         //     using index_t                          = typename T2::interval_t::value_t;
-        //     constexpr std::size_t prediction_order = T2::mesh_t::config::prediction_order;
+        //     constexpr std::size_t prediction_stencil_radius = T2::mesh_t::config::prediction_stencil_radius;
 
         //     std::size_t delta_l = reconstruct_level - level;
         //     if (delta_l == 0)
@@ -386,7 +386,7 @@ namespace samurai
         //         detail::multi_dim_loop(interp_coeff_values,
         //                                  [&](auto&... out_indices)
         //                                  {
-        //                                      auto& pred = prediction<prediction_order, index_t>(delta_l, out_indices...);
+        //                                      auto& pred = prediction<prediction_stencil_radius, index_t>(delta_l, out_indices...);
         //                                      for (const auto& kv : pred.coeff)
         //                                      {
         //                                          std::apply(
@@ -404,8 +404,8 @@ namespace samurai
         template <class T1, class T2>
         inline void operator()(Dim<1>, std::size_t& reconstruct_level, T1& dest, const T2& src) const
         {
-            using index_t                          = typename T2::interval_t::value_t;
-            constexpr std::size_t prediction_order = T2::mesh_t::config::prediction_order;
+            using index_t                                   = typename T2::interval_t::value_t;
+            constexpr std::size_t prediction_stencil_radius = T2::mesh_t::config::prediction_stencil_radius;
 
             std::size_t delta_l = reconstruct_level - level;
             if (delta_l == 0)
@@ -417,7 +417,7 @@ namespace samurai
                 index_t nb_cells = 1 << delta_l;
                 for (index_t ii = 0; ii < nb_cells; ++ii)
                 {
-                    auto pred = prediction<prediction_order, index_t>(delta_l, ii);
+                    auto pred = prediction<prediction_stencil_radius, index_t>(delta_l, ii);
                     for (const auto& kv : pred.coeff)
                     {
                         auto i_f = (i << delta_l) + ii;
@@ -432,8 +432,8 @@ namespace samurai
         inline void operator()(Dim<2>, std::size_t& reconstruct_level, T1& dest, const T2& src) const
 
         {
-            using index_t                          = typename T2::interval_t::value_t;
-            constexpr std::size_t prediction_order = T2::mesh_t::config::prediction_order;
+            using index_t                                   = typename T2::interval_t::value_t;
+            constexpr std::size_t prediction_stencil_radius = T2::mesh_t::config::prediction_stencil_radius;
 
             std::size_t delta_l = reconstruct_level - level;
             if (delta_l == 0)
@@ -448,7 +448,7 @@ namespace samurai
                     auto j_f = (j << delta_l) + jj;
                     for (index_t ii = 0; ii < nb_cells; ++ii)
                     {
-                        const auto& pred = prediction<prediction_order, index_t>(delta_l, ii, jj);
+                        const auto& pred = prediction<prediction_stencil_radius, index_t>(delta_l, ii, jj);
                         auto i_f         = (i << delta_l) + ii;
                         i_f.step         = nb_cells;
 
@@ -464,8 +464,8 @@ namespace samurai
         template <class T1, class T2>
         inline void operator()(Dim<3>, std::size_t& reconstruct_level, T1& dest, const T2& src) const
         {
-            using index_t                          = typename T2::interval_t::value_t;
-            constexpr std::size_t prediction_order = T2::mesh_t::config::prediction_order;
+            using index_t                                   = typename T2::interval_t::value_t;
+            constexpr std::size_t prediction_stencil_radius = T2::mesh_t::config::prediction_stencil_radius;
 
             std::size_t delta_l = reconstruct_level - level;
             if (delta_l == 0)
@@ -483,7 +483,7 @@ namespace samurai
                         auto j_f = (j << delta_l) + jj;
                         for (index_t ii = 0; ii < nb_cells; ++ii)
                         {
-                            const auto& pred = prediction<prediction_order, index_t>(delta_l, ii, jj, kk);
+                            const auto& pred = prediction<prediction_stencil_radius, index_t>(delta_l, ii, jj, kk);
                             auto i_f         = (i << delta_l) + ii;
                             i_f.step         = nb_cells;
 
@@ -514,6 +514,13 @@ namespace samurai
         using mesh_id_t = typename mesh_t::mesh_id_t;
         using ca_type   = typename mesh_t::ca_type;
 
+        if (field.mesh().max_stencil_radius() < 2)
+        {
+            std::cerr << "The reconstruction function requires at least 2 ghosts on the boundary.\nTo fix this issue, remove mesh_config.disable_minimal_ghost_width()."
+                      << std::endl;
+            exit(EXIT_FAILURE);
+        }
+
         update_ghost_mr_if_needed(field);
 
         auto make_field_like = [](std::string const& name, auto& mesh)
@@ -552,7 +559,7 @@ namespace samurai
 
     namespace detail
     {
-        template <std::size_t prediction_order, class Field, class... index_t>
+        template <std::size_t prediction_stencil_radius, class Field, class... index_t>
             requires(Field::dim == sizeof...(index_t) + 1 && (std::same_as<typename Field::interval_t, index_t> && ...))
         decltype(auto) get_prediction(std::size_t level, std::size_t delta_l, const std::tuple<index_t...>& ii)
         {
@@ -564,7 +571,7 @@ namespace samurai
             auto iter = std::apply(
                 [&](auto&... index)
                 {
-                    return values.find({prediction_order, level, index...});
+                    return values.find({prediction_stencil_radius, level, index...});
                 },
                 ii);
 
@@ -577,7 +584,9 @@ namespace samurai
                         std::apply(
                             [&](auto&... index)
                             {
-                                values[{prediction_order, level, index...}] += prediction<prediction_order, value_t>(delta_l, ii_...);
+                                values[{prediction_stencil_radius, level, index...}] += prediction<prediction_stencil_radius, value_t>(
+                                    delta_l,
+                                    ii_...);
                             },
                             ii);
                     });
@@ -586,12 +595,12 @@ namespace samurai
             return std::apply(
                 [&](auto&... index) -> auto&
                 {
-                    return values[{prediction_order, level, index...}];
+                    return values[{prediction_stencil_radius, level, index...}];
                 },
                 ii);
         }
 
-        template <std::size_t prediction_order, class Field, class... index_t>
+        template <std::size_t prediction_stencil_radius, class Field, class... index_t>
             requires((std::same_as<typename Field::interval_t::value_t, index_t> && ...))
         decltype(auto) get_prediction(std::size_t, std::size_t delta_l, const std::tuple<index_t...>& ii)
         {
@@ -599,12 +608,12 @@ namespace samurai
             return std::apply(
                 [delta_l](const auto&... index) -> auto&
                 {
-                    return prediction<prediction_order, value_t>(delta_l, index...);
+                    return prediction<prediction_stencil_radius, value_t>(delta_l, index...);
                 },
                 ii);
         }
 
-        template <std::size_t prediction_order, class Field, class Func, class... index_t, class... cell_index_t>
+        template <std::size_t prediction_stencil_radius, class Field, class Func, class... index_t, class... cell_index_t>
             requires(Field::dim == sizeof...(index_t) + 1 && Field::dim == sizeof...(cell_index_t)
                      && ((std::same_as<typename Field::interval_t, cell_index_t> && ...)
                          || (std::same_as<typename Field::interval_t::value_t, cell_index_t> && ...))
@@ -618,7 +627,7 @@ namespace samurai
         {
             using result_t = std::decay_t<decltype(result)>;
 
-            const auto& pred = get_prediction<prediction_order, Field>(level, delta_l, ii);
+            const auto& pred = get_prediction<prediction_stencil_radius, Field>(level, delta_l, ii);
 
             if constexpr (std::is_same_v<result_t, double>)
             {
@@ -661,7 +670,7 @@ namespace samurai
         {
             return f(element, level, indices...);
         };
-        detail::portion_impl<Field::mesh_t::config::prediction_order, Field>(result, f, get_f, level, delta_l, i, ii);
+        detail::portion_impl<Field::mesh_t::config::prediction_stencil_radius, Field>(result, f, get_f, level, delta_l, i, ii);
     }
 
     template <class Field, class... index_t, class... cell_index_t>
@@ -682,7 +691,7 @@ namespace samurai
         return result;
     }
 
-    template <std::size_t prediction_order, class Field, class... index_t, class... cell_index_t>
+    template <std::size_t prediction_stencil_radius, class Field, class... index_t, class... cell_index_t>
     void portion(auto& result,
                  const Field& f,
                  std::size_t level,
@@ -695,7 +704,7 @@ namespace samurai
             return f(level, indices...);
         };
 
-        detail::portion_impl<prediction_order, Field>(result, get_f, level, delta_l, i, ii);
+        detail::portion_impl<prediction_stencil_radius, Field>(result, get_f, level, delta_l, i, ii);
     }
 
     template <class Field, class... index_t, class... cell_index_t>
@@ -706,10 +715,10 @@ namespace samurai
                  const std::tuple<typename Field::interval_t, index_t...>& i,
                  const std::tuple<cell_index_t...>& ii)
     {
-        portion<Field::mesh_t::config::prediction_order>(result, f, level, delta_l, i, ii);
+        portion<Field::mesh_t::config::prediction_stencil_radius>(result, f, level, delta_l, i, ii);
     }
 
-    template <std::size_t prediction_order, class Field, class... index_t, class... cell_index_t>
+    template <std::size_t prediction_stencil_radius, class Field, class... index_t, class... cell_index_t>
     auto portion(const Field& f,
                  std::size_t level,
                  std::size_t delta_l,
@@ -722,7 +731,7 @@ namespace samurai
                 return zeros_like(f(level, indices...));
             },
             i);
-        portion<prediction_order>(result, f, level, delta_l, i, ii);
+        portion<prediction_stencil_radius>(result, f, level, delta_l, i, ii);
         return result;
     }
 
@@ -733,7 +742,7 @@ namespace samurai
                  const std::tuple<typename Field::interval_t, index_t...>& i,
                  const std::tuple<cell_index_t...>& ii)
     {
-        return portion<Field::mesh_t::config::prediction_order>(f, level, delta_l, i, ii);
+        return portion<Field::mesh_t::config::prediction_stencil_radius>(f, level, delta_l, i, ii);
     }
 
     namespace detail
@@ -794,7 +803,7 @@ namespace samurai
         auto src_tuple = detail::extract_src_tuple<dim, interval_t>(src_indices);
         auto dst_tuple = detail::extract_dst_tuple<dim>(delta_l, dst_indices);
 
-        detail::portion_impl<Field::mesh_t::config::prediction_order, Field>(result, get_f, level, delta_l, src_tuple, dst_tuple);
+        detail::portion_impl<Field::mesh_t::config::prediction_stencil_radius, Field>(result, get_f, level, delta_l, src_tuple, dst_tuple);
     }
 
     template <class Field_src, class Field_dst>
@@ -808,6 +817,13 @@ namespace samurai
         auto& mesh_src                   = field_src.mesh();
         auto& mesh_dst                   = field_dst.mesh();
 
+        if (field_src.mesh().max_stencil_radius() < 2)
+        {
+            std::cerr << "The transfer function requires at least 2 ghosts on the boundary.\nTo fix this issue, remove mesh_config.disable_minimal_ghost_width()."
+                      << std::endl;
+            exit(EXIT_FAILURE);
+        }
+
         update_ghost_mr_if_needed(field_src);
 
         field_dst.fill(0.);
diff --git a/include/samurai/samurai_config.hpp b/include/samurai/samurai_config.hpp
index b4394f61a..ab4782295 100644
--- a/include/samurai/samurai_config.hpp
+++ b/include/samurai/samurai_config.hpp
@@ -15,10 +15,10 @@ namespace samurai
 
     namespace default_config
     {
-        static constexpr std::size_t max_level        = 20;
-        static constexpr std::size_t ghost_width      = 1;
-        static constexpr std::size_t graduation_width = 1;
-        static constexpr std::size_t prediction_order = 1;
+        static constexpr std::size_t max_level         = 20;
+        static constexpr int ghost_width               = 1;
+        static constexpr int graduation_width          = 1;
+        static constexpr int prediction_stencil_radius = 1;
 
         using index_t    = signed long long int;
         using value_t    = int;
@@ -26,7 +26,7 @@ namespace samurai
 
         inline auto default_prediction_fn = [](auto& new_field, const auto& old_field) // cppcheck-suppress constParameterReference
         {
-            constexpr std::size_t pred_order = std::decay_t<decltype(new_field)>::mesh_t::config::prediction_order;
+            constexpr std::size_t pred_order = std::decay_t<decltype(new_field)>::mesh_t::config::prediction_stencil_radius;
             return prediction<pred_order, true>(new_field, old_field);
         };
     }
diff --git a/include/samurai/schemes/fv/explicit_FV_scheme.hpp b/include/samurai/schemes/fv/explicit_FV_scheme.hpp
index 5b6b096f3..c61a69662 100644
--- a/include/samurai/schemes/fv/explicit_FV_scheme.hpp
+++ b/include/samurai/schemes/fv/explicit_FV_scheme.hpp
@@ -75,6 +75,18 @@ namespace samurai
 
         virtual void apply(output_field_t& output_field, input_field_t& input_field)
         {
+            static constexpr int scheme_stencil_size = static_cast<int>(scheme_t::cfg_t::stencil_size);
+            int mesh_stencil_size                    = input_field.mesh().cfg().max_stencil_size();
+
+            if (scheme_stencil_size > mesh_stencil_size)
+            {
+                std::cerr << "The stencil size required by the scheme '" << scheme().name() << "' (" << scheme_stencil_size
+                          << ") is larger than the max_stencil_size parameter of the mesh (" << mesh_stencil_size
+                          << ").\nYou can set it with mesh_config.max_stencil_radius(" << scheme_stencil_size / 2
+                          << ") or mesh_config.max_stencil_size(" << scheme_stencil_size << ")." << std::endl;
+                exit(EXIT_FAILURE);
+            }
+
             for (std::size_t d = 0; d < dim; ++d)
             {
                 apply(d, output_field, input_field);
diff --git a/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp b/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
index c8c5eced7..18aaf8b63 100644
--- a/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
+++ b/include/samurai/schemes/fv/flux_based/flux_based_scheme__nonlin.hpp
@@ -201,7 +201,7 @@ namespace samurai
             {
                 copy_stencil_values(field, cells, stencil_values_list[0]);
             }
-            else if constexpr (mesh_t::config::prediction_order == 0 && stencil_size <= 4)
+            else if constexpr (mesh_t::config::prediction_stencil_radius == 0 && stencil_size <= 4)
             {
                 for (std::size_t fine_flux_index = 0; fine_flux_index < flux_params.n_fine_fluxes; ++fine_flux_index)
                 {
diff --git a/include/samurai/schemes/fv/operators/convection_lin.hpp b/include/samurai/schemes/fv/operators/convection_lin.hpp
index 8a414411f..67f9505ed 100644
--- a/include/samurai/schemes/fv/operators/convection_lin.hpp
+++ b/include/samurai/schemes/fv/operators/convection_lin.hpp
@@ -89,8 +89,6 @@ namespace samurai
     template <class Field>
     auto make_convection_weno5(const VelocityVector<Field::dim>& velocity)
     {
-        static_assert(Field::mesh_t::config::ghost_width >= 3, "WENO5 requires at least 3 ghosts.");
-
         static constexpr std::size_t dim          = Field::dim;
         static constexpr bool is_soa              = detail::is_soa_v<Field>;
         static constexpr std::size_t stencil_size = 6;
@@ -259,8 +257,6 @@ namespace samurai
         requires IsField<VelocityField>
     auto make_convection_weno5(VelocityField& velocity_field)
     {
-        static_assert(Field::mesh_t::config::ghost_width >= 3, "WENO5 requires at least 3 ghosts.");
-
         static constexpr std::size_t dim = Field::dim;
         static constexpr bool is_soa     = detail::is_soa_v<Field>;
 
diff --git a/include/samurai/schemes/fv/operators/convection_nonlin.hpp b/include/samurai/schemes/fv/operators/convection_nonlin.hpp
index ca93c7a83..bb3c4b880 100644
--- a/include/samurai/schemes/fv/operators/convection_nonlin.hpp
+++ b/include/samurai/schemes/fv/operators/convection_nonlin.hpp
@@ -84,8 +84,6 @@ namespace samurai
     {
         using field_value_t = typename Field::value_type;
 
-        static_assert(Field::mesh_t::config::ghost_width >= 3, "WENO5 requires at least 3 ghosts.");
-
         static constexpr std::size_t dim = Field::dim;
 
         static constexpr std::size_t stencil_size = 6;
diff --git a/include/samurai/static_algorithm.hpp b/include/samurai/static_algorithm.hpp
index d91a4f745..c948c1e31 100644
--- a/include/samurai/static_algorithm.hpp
+++ b/include/samurai/static_algorithm.hpp
@@ -188,6 +188,12 @@ namespace samurai
         detail::static_nested_loop_impl<nloops>(start, end, step, std::forward<Func>(f), index, std::integral_constant<std::size_t, 0>{});
     }
 
+    template <std::size_t nloops, class Func>
+    inline void static_nested_loop(int start, int end, Func&& f)
+    {
+        static_nested_loop<nloops>(start, end, 1, std::forward<Func>(f));
+    }
+
     /**
      * constexpr power function
      */
diff --git a/include/samurai/uniform_mesh.hpp b/include/samurai/uniform_mesh.hpp
index f16e928bf..7fc1694fe 100644
--- a/include/samurai/uniform_mesh.hpp
+++ b/include/samurai/uniform_mesh.hpp
@@ -9,6 +9,7 @@
 #include "level_cell_array.hpp"
 #include "level_cell_list.hpp"
 #include "mesh.hpp"
+#include "mesh_config.hpp"
 #include "samurai_config.hpp"
 
 namespace samurai
@@ -77,6 +78,8 @@ namespace samurai
         void set_scaling_factor(double scaling_factor);
         double cell_length(std::size_t level) const;
 
+        auto cfg() const;
+
       private:
 
         void update_sub_mesh();
@@ -241,6 +244,12 @@ namespace samurai
         return samurai::cell_length(scaling_factor(), level);
     }
 
+    template <class Config>
+    inline auto UniformMesh<Config>::cfg() const
+    {
+        return mesh_config<dim>();
+    }
+
     template <class Config>
     inline bool operator==(const UniformMesh<Config>& mesh1, const UniformMesh<Config>& mesh2)
     {
diff --git a/tests/test_adapt.cpp b/tests/test_adapt.cpp
index c4b86c39d..a87407a9f 100644
--- a/tests/test_adapt.cpp
+++ b/tests/test_adapt.cpp
@@ -30,12 +30,13 @@ namespace samurai
         ::samurai::initialize();
 
         static constexpr std::size_t dim = TypeParam::value;
-        using config                     = MRConfig<dim>;
-        using box_t                      = Box<double, dim>;
-        auto mesh                        = MRMesh<config>(box_t{xt::zeros<double>({dim}), xt::ones<double>({dim})}, 2, 4);
-        auto u_1                         = make_scalar_field<double>("u_1", mesh);
-        auto u_2                         = make_vector_field<double, 3, true>("u_2", mesh);
-        auto u_3                         = make_vector_field<double, 2>("u_3", mesh);
+
+        using box_t   = Box<double, dim>;
+        auto mesh_cfg = mesh_config<dim>().min_level(2).max_level(4);
+        auto mesh     = make_MRMesh(mesh_cfg, box_t{xt::zeros<double>({dim}), xt::ones<double>({dim})});
+        auto u_1      = make_scalar_field<double>("u_1", mesh);
+        auto u_2      = make_vector_field<double, 3, true>("u_2", mesh);
+        auto u_3      = make_vector_field<double, 2>("u_3", mesh);
 
         auto adapt      = make_MRAdapt(u_1, u_2, u_3);
         auto mra_config = samurai::mra_config().regularity(2.);
diff --git a/tests/test_bc.cpp b/tests/test_bc.cpp
index e0e428c21..8b4fb7aac 100644
--- a/tests/test_bc.cpp
+++ b/tests/test_bc.cpp
@@ -57,10 +57,10 @@ namespace samurai
     TEST(bc, scalar_function)
     {
         static constexpr std::size_t dim = 1;
-        using config                     = MRConfig<dim>;
 
         Box<double, dim> box = {{0}, {1}};
-        auto mesh            = MRMesh<config>(box, 2, 4);
+        auto mesh_cfg        = mesh_config<dim>().min_level(2).max_level(4);
+        auto mesh            = make_MRMesh(mesh_cfg, box);
         auto u               = make_scalar_field<double>("u", mesh);
 
         make_bc<Dirichlet<1>>(u,
diff --git a/tests/test_corner_projection.cpp b/tests/test_corner_projection.cpp
index a7e78e3f3..3f5da7a4c 100644
--- a/tests/test_corner_projection.cpp
+++ b/tests/test_corner_projection.cpp
@@ -9,19 +9,17 @@ namespace samurai
     TEST(subset, corner_projection)
     {
         static constexpr std::size_t dim = 2;
-        using Config                     = MRConfig<dim>;
         using Box                        = Box<double, dim>;
-        using Mesh                       = MRMesh<Config>;
-        using interval_t                 = typename Mesh::interval_t;
-        using mesh_id_t                  = typename Mesh::mesh_id_t;
         using direction_t                = DirectionVector<dim>;
 
         Box box({-1., -1.}, {1., 1.});
 
-        std::size_t min_level = 2;
-        std::size_t max_level = 6;
+        auto mesh_cfg = mesh_config<dim>().min_level(2).max_level(6);
+        auto mesh     = make_MRMesh(mesh_cfg, box);
 
-        Mesh mesh{box, min_level, max_level};
+        using Mesh       = decltype(mesh);
+        using interval_t = typename Mesh::interval_t;
+        using mesh_id_t  = typename Mesh::mesh_id_t;
 
         direction_t direction = {-1, -1}; // corner direction
         std::size_t level     = 6;
diff --git a/tests/test_domain_with_hole.cpp b/tests/test_domain_with_hole.cpp
index a69527962..20e06dc11 100644
--- a/tests/test_domain_with_hole.cpp
+++ b/tests/test_domain_with_hole.cpp
@@ -9,15 +9,15 @@ namespace samurai
     {
         static constexpr std::size_t dim = 2;
 
-        using Config    = samurai::MRConfig<dim>;
-        using Mesh      = samurai::MRMesh<Config>;
+        std::size_t level = 3;
+        auto mesh_cfg     = mesh_config<dim>().min_level(level).max_level(level);
+
+        using Mesh      = decltype(make_empty_MRMesh(mesh_cfg));
         using mesh_id_t = typename Mesh::mesh_id_t;
         using cl_t      = typename Mesh::cl_type;
         using lca_t     = typename Mesh::lca_type;
         using Box       = samurai::Box<double, dim>;
 
-        std::size_t level = 3;
-
         const Box domain_box({-1., -1.}, {1., 1.});
         const Box hole_box({0.0, 0.0}, {0.2, 0.2});
 
@@ -36,7 +36,7 @@ namespace samurai
                 domain_with_hole_cl[level][index_y].add_interval({interval});
             });
 
-        samurai::MRMesh<Config> mesh{domain_with_hole_cl, level, level};
+        auto mesh = samurai::make_MRMesh(mesh_cfg, domain_with_hole_cl);
 
         EXPECT_EQ(mesh.nb_cells(mesh_id_t::cells), domain_lca.nb_cells() - hole_lca.nb_cells());
     }
diff --git a/tests/test_field.cpp b/tests/test_field.cpp
index f8f6fb5d1..9c0ce447d 100644
--- a/tests/test_field.cpp
+++ b/tests/test_field.cpp
@@ -12,7 +12,6 @@ namespace samurai
     TEST(field, from_expr)
     {
         Box<double, 1> box{{0}, {1}};
-        // using Config = MRConfig<1>;
         // auto mesh    = MRMesh<Config>(box, 3, 3);
 
         using Config = UniformConfig<1>;
@@ -98,14 +97,14 @@ namespace samurai
 
     TEST(field, iterator)
     {
-        using config = MRConfig<2>;
         CellList<2> cl;
         cl[1][{0}].add_interval({0, 2});
         cl[1][{0}].add_interval({4, 6});
         cl[2][{0}].add_interval({4, 8});
 
-        auto mesh  = MRMesh<config>(cl, 1, 2);
-        auto field = make_scalar_field<std::size_t>("u", mesh);
+        auto mesh_cfg = mesh_config<2>().min_level(1).max_level(2);
+        auto mesh     = make_MRMesh(mesh_cfg, cl);
+        auto field    = make_scalar_field<std::size_t>("u", mesh);
 
         std::size_t index = 0;
         for_each_cell(mesh,
diff --git a/tests/test_find.cpp b/tests/test_find.cpp
index 8e1efd1ac..ec1b8ff49 100644
--- a/tests/test_find.cpp
+++ b/tests/test_find.cpp
@@ -10,17 +10,11 @@ namespace samurai
     TEST(find, test)
     {
         static constexpr std::size_t dim = 2;
-        using Config                     = samurai::MRConfig<dim>;
         using Box                        = samurai::Box<double, dim>;
-        using Mesh                       = samurai::MRMesh<Config>;
-        using coords_t                   = typename Mesh::cell_t::coords_t;
-
         Box box({-1., -1.}, {1., 1.});
 
-        std::size_t min_level = 2;
-        std::size_t max_level = 6;
-
-        Mesh mesh{box, min_level, max_level};
+        auto mesh_cfg = mesh_config<dim>().min_level(2).max_level(6);
+        auto mesh     = make_MRMesh(mesh_cfg, box);
 
         auto u = samurai::make_scalar_field<double>("u",
                                                     mesh,
@@ -35,6 +29,7 @@ namespace samurai
         auto mra_config   = samurai::mra_config().epsilon(1e-3);
         MRadaptation(mra_config);
 
+        using coords_t  = typename decltype(mesh)::cell_t::coords_t;
         coords_t coords = {0.4, 0.8};
         auto cell       = samurai::find_cell(mesh, coords);
 
diff --git a/tests/test_for_each.cpp b/tests/test_for_each.cpp
index d207b960c..3cf5e1052 100644
--- a/tests/test_for_each.cpp
+++ b/tests/test_for_each.cpp
@@ -5,8 +5,8 @@ namespace samurai
 {
     auto create_meshes(std::size_t level)
     {
-        using Config  = amr::Config<1>;
-        using Mesh    = amr::Mesh<Config>;
+        using Config  = mesh_config<1>;
+        using Mesh    = decltype(amr::make_empty_Mesh(std::declval<Config>()));
         using cl_type = typename Mesh::cl_type;
 
         cl_type cl1;
@@ -14,22 +14,23 @@ namespace samurai
         cl_type cl2;
         cl2[level][{}].add_interval({0, 3});
 
-        auto m1 = Mesh(cl1, level, level);
-        auto m2 = Mesh(cl2, level, level);
+        auto mesh_cfg = mesh_config<1>().min_level(level).max_level(level).disable_minimal_ghost_width();
+        auto m1       = amr::make_Mesh(mesh_cfg, cl1);
+        auto m2       = amr::make_Mesh(mesh_cfg, cl2);
         return std::make_tuple(m1, m2);
     }
 
     TEST(set, for_each_interval)
     {
-        using Config    = amr::Config<1>;
-        using Mesh      = amr::Mesh<Config>;
-        using mesh_id_t = typename Mesh::mesh_id_t;
-
         std::size_t level = 1;
         auto meshes       = create_meshes(level);
         auto& m1          = std::get<0>(meshes);
         auto& m2          = std::get<1>(meshes);
-        auto set          = intersection(m1[mesh_id_t::cells][level], m2[mesh_id_t::cells][level]);
+
+        using Mesh      = std::tuple_element_t<0, decltype(meshes)>;
+        using mesh_id_t = Mesh::mesh_id_t;
+
+        auto set = intersection(m1[mesh_id_t::cells][level], m2[mesh_id_t::cells][level]);
 
         int nb_intervals = 0;
         for_each_interval(set,
@@ -45,15 +46,15 @@ namespace samurai
 
     TEST(set, for_each_cell)
     {
-        using Config    = amr::Config<1>;
-        using Mesh      = amr::Mesh<Config>;
-        using mesh_id_t = typename Mesh::mesh_id_t;
-
         std::size_t level = 1;
         auto meshes       = create_meshes(level);
         auto& m1          = std::get<0>(meshes);
         auto& m2          = std::get<1>(meshes);
-        auto set          = intersection(m1[mesh_id_t::cells][level], m2[mesh_id_t::cells][level]);
+
+        using Mesh      = std::tuple_element_t<0, decltype(meshes)>;
+        using mesh_id_t = Mesh::mesh_id_t;
+
+        auto set = intersection(m1[mesh_id_t::cells][level], m2[mesh_id_t::cells][level]);
 
         /** Cell indices in m1:
          *
diff --git a/tests/test_hdf5.cpp b/tests/test_hdf5.cpp
index f3567a707..d0d21de87 100644
--- a/tests/test_hdf5.cpp
+++ b/tests/test_hdf5.cpp
@@ -87,14 +87,14 @@ namespace samurai
     TYPED_TEST(hdf5_test, mr_mesh)
     {
         static constexpr std::size_t dim = TypeParam::value;
-        using Config                     = MRConfig<dim>;
-        using Mesh                       = MRMesh<Config>;
+
         xt::xtensor_fixed<double, xt::xshape<dim>> min_corner;
         xt::xtensor_fixed<double, xt::xshape<dim>> max_corner;
         min_corner.fill(-1);
         max_corner.fill(1);
         Box<double, dim> box(min_corner, max_corner);
-        Mesh mesh(box, 4, 4);
+        auto mesh_cfg = mesh_config<dim>().min_level(4).max_level(4);
+        auto mesh     = make_MRMesh(mesh_cfg, box);
         test_save(mesh);
         args::save_debug_fields = true;
         test_save(mesh);
diff --git a/tests/test_mra.cpp b/tests/test_mra.cpp
index 0eeb51cde..b5005b4a0 100644
--- a/tests/test_mra.cpp
+++ b/tests/test_mra.cpp
@@ -11,12 +11,12 @@ namespace samurai
     auto init_mesh()
     {
         constexpr std::size_t dim   = 2;
-        using Config                = samurai::MRConfig<dim>;
         auto box                    = samurai::Box<double, dim>({0., 0.}, {1., 1.});
         const std::size_t min_level = 2;
         const std::size_t max_level = 10;
 
-        return samurai::MRMesh<Config>(box, min_level, max_level);
+        auto mesh_cfg = mesh_config<dim>().min_level(min_level).max_level(max_level);
+        return make_MRMesh(mesh_cfg, box);
     }
 
     void init_field(auto& u, double factor = 1.0)
diff --git a/tests/test_periodic.cpp b/tests/test_periodic.cpp
index 9405df454..85b2c0031 100644
--- a/tests/test_periodic.cpp
+++ b/tests/test_periodic.cpp
@@ -45,7 +45,6 @@ namespace samurai
     TYPED_TEST(dim_test, periodic)
     {
         static constexpr std::size_t dim = TypeParam::value;
-        using Config                     = MRConfig<dim, 1>;
 
         // Simulation parameters
         xt::xtensor_fixed<double, xt::xshape<dim>> min_corner, max_corner;
@@ -53,16 +52,13 @@ namespace samurai
         max_corner.fill(1);
 
         // Multiresolution parameters
-        std::size_t min_level = 3, max_level = 6;
-
         Box<double, dim> box(min_corner, max_corner);
-        std::array<bool, dim> bc;
-        bc.fill(true);
-        MRMesh<Config> mesh{box, min_level, max_level, bc};
-        using mesh_id_t = typename MRMesh<Config>::mesh_id_t;
+        auto mesh_cfg   = samurai::mesh_config<dim>().min_level(3).max_level(6).periodic(true).graduation_width(1);
+        auto mesh       = samurai::make_MRMesh(mesh_cfg, box);
+        using mesh_id_t = typename decltype(mesh)::mesh_id_t;
 
         double dt = 1;
-        double Tf = 2 / mesh.cell_length(max_level);
+        double Tf = 2 / mesh.cell_length(mesh.max_level());
         double t  = 0.;
 
         auto u    = init(mesh);
diff --git a/tests/test_restart.cpp b/tests/test_restart.cpp
index 094c87df6..b6b64e4da 100644
--- a/tests/test_restart.cpp
+++ b/tests/test_restart.cpp
@@ -10,8 +10,6 @@ namespace samurai
     template <std::size_t dim>
     auto create_mesh(double box_boundary)
     {
-        using Config  = samurai::MRConfig<dim>;
-        using Mesh    = samurai::MRMesh<Config>;
         using Box     = samurai::Box<double, dim>;
         using point_t = typename Box::point_t;
 
@@ -19,8 +17,9 @@ namespace samurai
         box_corner1.fill(0);
         box_corner2.fill(box_boundary);
         Box box(box_corner1, box_corner2);
+        auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(5).disable_minimal_ghost_width();
 
-        return Mesh(box, 2, 5);
+        return make_MRMesh(mesh_cfg, box);
     }
 
     TEST(restart, restart_lca)
diff --git a/tests/test_scaling.cpp b/tests/test_scaling.cpp
index bd35f93e2..f5fd05d43 100644
--- a/tests/test_scaling.cpp
+++ b/tests/test_scaling.cpp
@@ -8,8 +8,6 @@ namespace samurai
     template <std::size_t dim>
     auto create_mesh(double box_boundary, std::size_t level)
     {
-        using Config  = samurai::MRConfig<dim>;
-        using Mesh    = samurai::MRMesh<Config>;
         using Box     = samurai::Box<double, dim>;
         using point_t = typename Box::point_t;
 
@@ -17,8 +15,9 @@ namespace samurai
         box_corner1.fill(0);
         box_corner2.fill(box_boundary);
         Box box(box_corner1, box_corner2);
+        auto mesh_cfg = samurai::mesh_config<dim>().min_level(level).max_level(level);
 
-        return Mesh(box, level, level);
+        return samurai::make_MRMesh(mesh_cfg, box);
     }
 
     /**
diff --git a/tests/test_subset.cpp b/tests/test_subset.cpp
index 71e28144e..8406304d1 100644
--- a/tests/test_subset.cpp
+++ b/tests/test_subset.cpp
@@ -483,10 +483,12 @@ namespace samurai
         }
 
         {
-            using Config = MRConfig<2>;
-            const Box<double, 2> box({0, 0}, {1, 1});
-            MRMesh<Config> mesh{box, 0, 3};
-            auto& domain                              = mesh.domain();
+            static constexpr std::size_t dim = 2;
+            const Box<double, dim> box({0, 0}, {1, 1});
+            auto mesh_cfg = samurai::mesh_config<dim>().min_level(0).max_level(3);
+            auto mesh     = samurai::make_MRMesh(mesh_cfg, box);
+            auto& domain  = mesh.domain();
+
             xt::xtensor_fixed<int, xt::xshape<2>> dir = {0, 1 << (3 - 1)};
 
             auto expected = expected_t{