Add 3D models to executable

dune/copasi/common/muparser_data_handler.hh

@@ -108,4 +108,4 @@ std::vector<std::shared_ptr<T>> MuParserDataHandler<T>::_functions = {};
 
 } // namespace Dune::Copasi
 
-#endif // DUNE_MUPARSER_DATA_HANDLER_HH
\ No newline at end of file
+#endif // DUNE_MUPARSER_DATA_HANDLER_HH

lib/dune_copasi_md.hh

@@ -11,12 +11,16 @@
 namespace Dune {
 namespace Copasi {
 
-constexpr int dim = 2;
-using HostGrid = Dune::UGGrid<dim>;
-using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
-using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+using HostGrid2D = Dune::UGGrid<2>;
+using MDGTraits2D = Dune::mdgrid::DynamicSubDomainCountTraits<2, 1>;
+using Grid2D = Dune::mdgrid::MultiDomainGrid<HostGrid2D, MDGTraits2D>;
+
+
+using HostGrid3D = Dune::UGGrid<3>;
+using MDGTraits3D = Dune::mdgrid::DynamicSubDomainCountTraits<3, 1>;
+using Grid3D = Dune::mdgrid::MultiDomainGrid<HostGrid3D, MDGTraits3D>;
 
 } // namespace Dorie
 } // namespace Dune
 
-#endif // DUNE_COPASI_MD_HH
\ No newline at end of file
+#endif // DUNE_COPASI_MD_HH

lib/dune_copasi_md_cg.cc

@@ -7,13 +7,22 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits1 =
-  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid, 1>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits1>;
+using ModelTraits2DPk1 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid2D, 1>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk1>;
+
+using ModelTraits2DPk2 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid2D, 2>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk2>;
+
+using ModelTraits3DPk1 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid3D, 1>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk1>;
+
+using ModelTraits3DPk2 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid3D, 2>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk2>;
 
-using ModelTraits2 =
-  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid, 2>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits2>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

lib/dune_copasi_md_fv.cc

@@ -7,9 +7,13 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits0 =
-  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid, 0>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits0>;
+using ModelTraits2DPk0 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid2D, 0>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk0>;
+
+using ModelTraits3DPk0 =
+  Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid3D, 0>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk0>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

lib/dune_copasi_md_fv_cg.cc

@@ -7,13 +7,29 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits01 =
-  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, 1>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits01>;
+using ModelTraits2DPk0 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid2D, 0>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk0>;
 
-using ModelTraits02 =
-  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, 2>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits02>;
+using ModelTraits2DPk1 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid2D, 1>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk1>;
+
+using ModelTraits2DPk2 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid2D, 2>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits2DPk2>;
+
+using ModelTraits3DPk0 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid3D, 0>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk0>;
+
+using ModelTraits3DPk1 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid3D, 1>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk1>;
+
+using ModelTraits3DPk2 =
+  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid3D, 2>;
+template class ModelMultiDomainDiffusionReaction<ModelTraits3DPk2>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

lib/dune_copasi_md_fv_cg.cc.orig

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "dune_copasi_md.hh"
-
-namespace Dune {
-namespace Copasi {
-
-using ModelTraits01 =
-  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, 1>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits01>;
-
-using ModelTraits02 =
-  Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, 2>;
-template class ModelMultiDomainDiffusionReaction<ModelTraits02>;
-
-
-} // namespace Dorie
-} // namespace Dune
\ No newline at end of file

lib/dune_copasi_sd.hh

@@ -11,15 +11,23 @@
 namespace Dune {
 namespace Copasi {
 
-constexpr int dim = 2;
-using HostGrid = Dune::UGGrid<dim>;
-using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
-using MDGrid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+using HostGrid2D = Dune::UGGrid<2>;
+using MDGTraits2D = Dune::mdgrid::DynamicSubDomainCountTraits<2, 1>;
+using MDGrid2D = Dune::mdgrid::MultiDomainGrid<HostGrid2D, MDGTraits2D>;
+
+using Grid2D = typename MDGrid2D::SubDomainGrid;
+using GridView2D = typename Grid2D::Traits::LeafGridView;
+
+
+using HostGrid3D = Dune::UGGrid<3>;
+using MDGTraits3D = Dune::mdgrid::DynamicSubDomainCountTraits<3, 1>;
+using MDGrid3D = Dune::mdgrid::MultiDomainGrid<HostGrid3D, MDGTraits3D>;
+
+using Grid3D = typename MDGrid3D::SubDomainGrid;
+using GridView3D = typename Grid3D::Traits::LeafGridView;
 
-using Grid = typename MDGrid::SubDomainGrid;
-using GridView = typename Grid::Traits::LeafGridView;
 
 } // namespace Dorie
 } // namespace Dune
 
-#endif // DUNE_COPASI_SD_HH
\ No newline at end of file
+#endif // DUNE_COPASI_SD_HH

lib/dune_copasi_sd_cg.cc

@@ -7,13 +7,22 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits1 =
-  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid, GridView, 1>;
-template class ModelDiffusionReaction<ModelTraits1>;
+using ModelTraits2DPk1 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid2D, GridView2D, 1>;
+template class ModelDiffusionReaction<ModelTraits2DPk1>;
+
+using ModelTraits2DPk2 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid2D, GridView2D, 2>;
+template class ModelDiffusionReaction<ModelTraits2DPk2>;
+
+using ModelTraits3DPk1 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid3D, GridView3D, 1>;
+template class ModelDiffusionReaction<ModelTraits3DPk1>;
+
+using ModelTraits3DPk2 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid3D, GridView3D, 2>;
+template class ModelDiffusionReaction<ModelTraits3DPk2>;
 
-using ModelTraits2 =
-  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid, GridView, 2>;
-template class ModelDiffusionReaction<ModelTraits2>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

lib/dune_copasi_sd_fv.cc

@@ -7,9 +7,13 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits0 =
-  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid, GridView, 0>;
-template class ModelDiffusionReaction<ModelTraits0>;
+using ModelTraits2DPk0 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid2D, GridView2D, 0>;
+template class ModelDiffusionReaction<ModelTraits2DPk0>;
+
+using ModelTraits3DPk0 =
+  Dune::Copasi::ModelPkDiffusionReactionTraits<Grid3D, GridView3D, 0>;
+template class ModelDiffusionReaction<ModelTraits3DPk0>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

lib/dune_copasi_sd_fv_cg.cc

@@ -7,13 +7,30 @@
 namespace Dune {
 namespace Copasi {
 
-using ModelTraits01 =
-  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid, GridView, 1>;
-template class ModelDiffusionReaction<ModelTraits01>;
 
-using ModelTraits02 =
-  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid, GridView, 2>;
-template class ModelDiffusionReaction<ModelTraits02>;
+using ModelTraits2DPk0 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid2D, GridView2D, 0>;
+template class ModelDiffusionReaction<ModelTraits2DPk0>;
+
+using ModelTraits2DPk1 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid2D, GridView2D, 1>;
+template class ModelDiffusionReaction<ModelTraits2DPk1>;
+
+using ModelTraits2DPk2 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid2D, GridView2D, 2>;
+template class ModelDiffusionReaction<ModelTraits2DPk2>;
+
+using ModelTraits3DPk0 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid3D, GridView3D, 0>;
+template class ModelDiffusionReaction<ModelTraits3DPk0>;
+
+using ModelTraits3DPk1 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid3D, GridView3D, 1>;
+template class ModelDiffusionReaction<ModelTraits3DPk1>;
+
+using ModelTraits3DPk2 =
+  Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid3D, GridView3D, 2>;
+template class ModelDiffusionReaction<ModelTraits3DPk2>;
 
 } // namespace Dorie
-} // namespace Dune
\ No newline at end of file
+} // namespace Dune

src/dune_copasi_md.cc

@@ -16,6 +16,7 @@
 #include <dune/logging/logging.hh>
 
 #include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/parametertreeparser.hh>
@@ -61,74 +62,80 @@ main(int argc, char** argv)
       log.detail(2, "{}"_fmt, line);
     log.detail(2, "----"_fmt);
 
-    // create a grid
-    constexpr int dim = 2;
-    using HostGrid = Dune::UGGrid<dim>;
-    using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
-    using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+    int order = config.template get<int>("model.order");
+    int dim = config.template get<int>("grid.dimension");
 
-    auto& grid_config = config.sub("grid", true);
-    auto level = grid_config.get<int>("initial_level", 0);
+    if (dim != 2 and dim != 3)
+      DUNE_THROW(Dune::IOError, "Only 2D and 3D grids are alloed!");
 
-    auto grid_file = grid_config.get<std::string>("file");
+    // lambda that instantiates and evolves a multidomain model
+    auto evolve_model = [](const auto& config, auto dim, auto order) {
+      using namespace Dune::Copasi;
 
-    auto [md_grid_ptr, host_grid_ptr] =
-      Dune::Copasi::MultiDomainGmshReader<Grid>::read(grid_file);
+      using HostGrid = Dune::UGGrid<dim>;
+      using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
+      using Grid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
 
-    auto grid_log = Dune::Logging::Logging::componentLogger({}, "grid");
-    grid_log.detail("Applying refinement of level: {}"_fmt, level);
+      auto& grid_config = config.sub("grid", true);
+      auto level = grid_config.template get<int>("initial_level");
 
-    for (int i = 1; i <= level; i++) {
-      Dune::Copasi::mark_stripes(*host_grid_ptr);
-      md_grid_ptr->preAdapt();
-      md_grid_ptr->adapt();
-      md_grid_ptr->postAdapt();
-    }
+      auto grid_file = grid_config.template get<std::string>("file");
 
-    auto& model_config = config.sub("model", true);
-    int order = model_config.get<int>("order");
+      auto [md_grid_ptr, host_grid_ptr] =
+        MultiDomainGmshReader<Grid>::read(grid_file);
 
-    // create time stepper
-    auto timestep_config = model_config.sub("time_stepping", true);
-    auto end_time = timestep_config.template get<double>("end");
-    auto initial_step = timestep_config.template get<double>("initial_step");
-    auto stepper = Dune::Copasi::make_default_stepper(timestep_config);
+      auto grid_log = Dune::Logging::Logging::componentLogger({}, "grid");
+      grid_log.detail("Applying refinement of level: {}"_fmt, level);
 
-    auto file = model_config.get("writer.file_path", "");
-    auto write_output = [=](const auto& state) {
-      if (not file.empty())
-        state.write(file, true);
-    };
+      for (int i = 1; i <= level; i++) {
+        mark_stripes(*host_grid_ptr);
+        md_grid_ptr->preAdapt();
+        md_grid_ptr->adapt();
+        md_grid_ptr->postAdapt();
+      }
 
-    if (order == 0) {
-      constexpr int Order = 0;
-      using ModelTraits =
-        Dune::Copasi::ModelMultiDomainPkDiffusionReactionTraits<Grid, Order>;
-      Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
-        md_grid_ptr, model_config);
-      write_output(model.state()); // write initial condition
-      stepper.evolve(model, initial_step, end_time, write_output);
-    } else if (order == 1) {
-      constexpr int Order = 1;
-      using ModelTraits =
-        Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, Order>;
-      Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
-        md_grid_ptr, model_config);
-      write_output(model.state()); // write initial condition
-      stepper.evolve(model, initial_step, end_time, write_output);
-    } else if (order == 2) {
-      constexpr int Order = 2;
-      using ModelTraits =
-        Dune::Copasi::ModelMultiDomainP0PkDiffusionReactionTraits<Grid, Order>;
-      Dune::Copasi::ModelMultiDomainDiffusionReaction<ModelTraits> model(
-        md_grid_ptr, model_config);
+      auto& model_config = config.sub("model", true);
+      using Traits = ModelMultiDomainP0PkDiffusionReactionTraits<Grid, order>;
+      ModelMultiDomainDiffusionReaction<Traits> model(md_grid_ptr,
+                                                      model_config);
+
+      // setup writer
+      auto file = model_config.get("writer.file_path", "");
+      auto write_output = [=](const auto& state) {
+        if (not file.empty())
+          state.write(file, true);
+      };
       write_output(model.state()); // write initial condition
+
+      // create time stepper
+      auto timestep_config = model_config.sub("time_stepping", true);
+      auto end_time = timestep_config.template get<double>("end");
+      auto initial_step = timestep_config.template get<double>("initial_step");
+      auto stepper = Dune::Copasi::make_default_stepper(timestep_config);
+
       stepper.evolve(model, initial_step, end_time, write_output);
-    } else {
+    };
+
+    // maximum polynomial order instantiated
+    constexpr auto max_order = 2;
+    if (order > max_order)
       DUNE_THROW(Dune::IOError,
                  "Finite element order " << order
                                          << " is not supported by dune-copasi");
-    }
+
+    // static loop that instantiates 2D and 3D polynomial orders but runs the
+    // dynamic ones
+    auto order_range = Dune::range(Dune::index_constant<max_order + 1>{});
+    Dune::Hybrid::forEach(order_range, [&](auto static_order) {
+      // instantiate models with 2D grids
+      if (dim == 2 and order == static_order)
+        evolve_model(config, Dune::Indices::_2, static_order);
+
+      // instantiate models with 3D grids
+      if (dim == 3 and order == static_order)
+        evolve_model(config, Dune::Indices::_3, static_order);
+    });
+
   } catch (Dune::Exception& e) {
     log.error("Dune reported error:"_fmt);
     log.error(2, "{}"_fmt, e.what());

src/dune_copasi_sd.cc

@@ -16,6 +16,7 @@
 #include <dune/logging/logging.hh>
 
 #include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/parametertreeparser.hh>
@@ -61,90 +62,93 @@ main(int argc, char** argv)
       log.detail(2, "{}"_fmt, line);
     log.detail(2, "----"_fmt);
 
-    // Grid setup
-    constexpr int dim = 2;
-    using HostGrid = Dune::UGGrid<dim>;
-    using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
-    using MDGrid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
-
-    auto& grid_config = config.sub("grid", true);
-    auto level = grid_config.get<int>("initial_level", 0);
-
-    auto grid_file = grid_config.get<std::string>("file");
-
-    auto [md_grid_ptr, host_grid_ptr] =
-      Dune::Copasi::MultiDomainGmshReader<MDGrid>::read(grid_file);
-
-    using Grid = typename MDGrid::SubDomainGrid;
-    using GridView = typename Grid::Traits::LeafGridView;
-
-    auto grid_log = Dune::Logging::Logging::componentLogger({}, "grid");
-    grid_log.detail("Applying refinement of level: {}"_fmt, level);
-
-    for (int i = 1; i <= level; i++) {
-      Dune::Copasi::mark_stripes(*host_grid_ptr);
-      md_grid_ptr->preAdapt();
-      md_grid_ptr->adapt();
-      md_grid_ptr->postAdapt();
-    }
-
-    auto& model_config = config.sub("model", true);
-    auto& compartments_map = model_config.sub("compartments", true);
-    int order = model_config.get<int>("order");
-
-    // TODO: Use OS for different domains when is ready
-    if (compartments_map.getValueKeys().size() != 1)
-      DUNE_THROW(
-        Dune::NotImplemented,
-        "Multiple compartments per model are not allowed in this executable");
-
-    // get subdomain grid as a shared pointer
-    auto compartment = compartments_map.getValueKeys().front();
-    int domain = compartments_map.template get<int>(compartment);
-    std::shared_ptr<Grid> grid_ptr =
-      Dune::stackobject_to_shared_ptr(md_grid_ptr->subDomain(domain));
-
-    // create time stepper
-    auto timestep_config = model_config.sub("time_stepping", true);
-    auto end_time = timestep_config.template get<double>("end");
-    auto initial_step = timestep_config.template get<double>("initial_step");
-    auto stepper = Dune::Copasi::make_default_stepper(timestep_config);
-
-    auto file = model_config.get("writer.file_path", "");
-    auto write_output = [=](const auto& state) {
-      if (not file.empty())
-        state.write(file, true);
-    };
-
-    if (order == 0) {
-      constexpr int Order = 0;
-      using ModelTraits =
-        Dune::Copasi::ModelPkDiffusionReactionTraits<Grid, GridView, Order>;
-      Dune::Copasi::ModelDiffusionReaction<ModelTraits> model(grid_ptr,
-                                                              model_config);
-      write_output(model.state()); // write initial condition
-      stepper.evolve(model, initial_step, end_time, write_output);
-    } else if (order == 1) {
-      constexpr int Order = 1;
-      using ModelTraits =
-        Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid, GridView, Order>;
-      Dune::Copasi::ModelDiffusionReaction<ModelTraits> model(grid_ptr,
-                                                              model_config);
-      write_output(model.state()); // write initial condition
-      stepper.evolve(model, initial_step, end_time, write_output);
-    } else if (order == 2) {
-      constexpr int Order = 2;
-      using ModelTraits =
-        Dune::Copasi::ModelP0PkDiffusionReactionTraits<Grid, GridView, Order>;
-      Dune::Copasi::ModelDiffusionReaction<ModelTraits> model(grid_ptr,
-                                                              model_config);
+    int order = config.template get<int>("model.order");
+    int dim = config.template get<int>("grid.dimension");
+
+    if (dim != 2 and dim != 3)
+      DUNE_THROW(Dune::IOError, "Only 2D and 3D grids are alloed!");
+
+    // lambda that instantiates and evolves a multidomain model
+    auto evolve_model = [](const auto& config, auto dim, auto order) {
+      using namespace Dune::Copasi;
+
+      using HostGrid = Dune::UGGrid<dim>;
+      using MDGTraits = Dune::mdgrid::DynamicSubDomainCountTraits<dim, 1>;
+      using MDGrid = Dune::mdgrid::MultiDomainGrid<HostGrid, MDGTraits>;
+
+      auto& grid_config = config.sub("grid", true);
+      auto level = grid_config.template get<int>("initial_level");
+      auto grid_file = grid_config.template get<std::string>("file");
+
+      auto [md_grid_ptr, host_grid_ptr] =
+        MultiDomainGmshReader<MDGrid>::read(grid_file);
+
+      auto grid_log = Dune::Logging::Logging::componentLogger({}, "grid");
+      grid_log.detail("Applying refinement of level: {}"_fmt, level);
+
+      for (int i = 1; i <= level; i++) {
+        mark_stripes(*host_grid_ptr);
+        md_grid_ptr->preAdapt();
+        md_grid_ptr->adapt();
+        md_grid_ptr->postAdapt();
+      }
+
+      auto& model_config = config.sub("model", true);
+      // TODO: Use OS for different domains when is ready
+      auto& compartments_map = model_config.sub("compartments", true);
+      if (compartments_map.getValueKeys().size() != 1)
+        DUNE_THROW(
+          Dune::NotImplemented,
+          "Multiple compartments per model are not allowed in this executable");
+
+      // get subdomain grid as a shared pointer
+      auto compartment = compartments_map.getValueKeys().front();
+      int domain = compartments_map.template get<int>(compartment);
+      auto grid_ptr =
+        Dune::stackobject_to_shared_ptr(md_grid_ptr->subDomain(domain));
+
+      using Grid = typename MDGrid::SubDomainGrid;
+      using GridView = typename Grid::Traits::LeafGridView;
+      using Traits = ModelP0PkDiffusionReactionTraits<Grid, GridView, order>;
+      ModelDiffusionReaction<Traits> model(grid_ptr, model_config);
+
+      // setup writer
+      auto file = model_config.get("writer.file_path", "");
+      auto write_output = [=](const auto& state) {
+        if (not file.empty())
+          state.write(file, true);
+      };
       write_output(model.state()); // write initial condition
+
+      // create time stepper
+      auto timestep_config = model_config.sub("time_stepping", true);
+      auto end_time = timestep_config.template get<double>("end");
+      auto initial_step = timestep_config.template get<double>("initial_step");
+      auto stepper = Dune::Copasi::make_default_stepper(timestep_config);
+
       stepper.evolve(model, initial_step, end_time, write_output);
-    } else {
+    };
+
+    // maximum polynomial order instantiated
+    constexpr auto max_order = 2;
+    if (order > max_order)
       DUNE_THROW(Dune::IOError,
-                 "Finite element order "
-                   << order << " is not supported by dune-copasi(sd)");
-    }
+                 "Finite element order " << order
+                                         << " is not supported by dune-copasi");
+
+    // static loop that instantiates 2D and 3D polynomial orders but runs the
+    // dynamic ones
+    auto order_range = Dune::range(Dune::index_constant<max_order + 1>{});
+    Dune::Hybrid::forEach(order_range, [&](auto static_order) {
+      // instantiate models with 2D grids
+      if (dim == 2 and order == static_order)
+        evolve_model(config, Dune::Indices::_2, static_order);
+
+      // instantiate models with 3D grids
+      if (dim == 3 and order == static_order)
+        evolve_model(config, Dune::Indices::_3, static_order);
+    });
+
   } catch (Dune::Exception& e) {
     log.error("Dune reported error:"_fmt);
     log.error(2, "{}"_fmt, e.what());

test/CMakeLists.txt

@@ -82,6 +82,10 @@ if (dune-testtools_FOUND)
     SCRIPT          dune_vtkcompare.py
     INIFILE         test_cell.mini)
 
+  add_system_test_per_target(
+    TARGET          dune_copasi_md
+    INIFILE         test_in_out.mini)
+
   add_system_test_per_target(
     TARGET          dune_copasi_sd
     INIFILE         test_no_writer.mini)

test/NFKappaB.mini

 __name = NFkappaB
 
 [grid]
+dimension = 2
 file = grids/test_cell.msh
 initial_level = 0
 

test/grids/test_cell_membrane.geo

 // Refinement value
-r = 2e-1;
-mr = 4e-2;
-
+refinement = DefineNumber[ 2, Name "Parameters/refinement" ];
 
 // Cytoplasm radious
 cr = 2;
 // Nucleus radious
 nr = 0.7;
 // Membrane width
-mw = 0.05;
+mw = 0.01;
 
 // Cytoplasm, Nucleus, and Membrane points
-Point(1) = {-cr, 0, 0, r};
-Point(2) = {-nr-mw, 0, 0, mr};
-Point(3) = {-nr, 0, 0, mr};
-Point(4) = {0, 0, 0, r};
-Point(5) = {nr, 0, 0, mr};
-Point(6) = {nr+mw, 0, 0, mr};
-Point(7) = {cr, 0, 0, r};
+Point(1) = {-cr, 0, 0, refinement};
+Point(2) = {-nr-mw, 0, 0, refinement};
+Point(3) = {-nr, 0, 0, refinement};
+Point(4) = {0, 0, 0, refinement};
+Point(5) = {nr, 0, 0, refinement};
+Point(6) = {nr+mw, 0, 0, refinement};
+Point(7) = {cr, 0, 0, refinement};
 
 // outer cytoplasm lines
 Circle(1) = {1,4,7};
@@ -50,4 +48,10 @@ Physical Surface("Membrane") = {2};
 Plane Surface(3) = {3};
 Physical Surface("Nucleus") = {3};
 
-Point{4} In Surface{3};
\ No newline at end of file
+Point{4} In Surface{3};
+
+membrane_size = DefineNumber[ 20, Name "Parameters/membrane_size" ];
+
+Transfinite Curve {5, 6, 4, 3} = membrane_size Using Progression 1;
+
+Recombine Surface {2};

test/test_cell.mini

 __name = test_cell
 
 [grid]
+dimension = 2
 file = grids/test_cell.msh
 initial_level = 0
 

test/test_exp.mini

@@ -4,6 +4,7 @@ _geometry_type = cubes, triangles | expand geometry_type
 __name = test_exp_{_geometry_type}
 
 [grid]
+dimension = 2
 file = grids/square_{_geometry_type}.msh
 initial_level = 0
 

test/test_gauss.mini

@@ -6,6 +6,7 @@ _geometry_type = cubes, triangles | expand geometry_type
 __name = test_gauss_{_geometry_type}
 
 [grid]
+dimension = 2
 file = grids/square_{_geometry_type}.msh
 initial_level = 4
 

test/test_gauss_mixed_geometry.mini

@@ -5,6 +5,7 @@ _t0 = 1.
 _gauss_eq = exp(-((x-0.5)^2+(y-0.5)^2)/(4*t*{_diffusion})) / (4*pi*t*{_diffusion})
 
 [grid]
+dimension = 2
 file = grids/square_mixed_geometry.msh
 initial_level = 4