diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index 1085ad501533d53e70d075e3d10a7da3c0fee2e7..76ad38a2c0229c7cddcf9661672f925a0066ec95 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -27,6 +27,24 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 #endif
 
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianProximalNewtonSolver<Basis, TargetSpace, Assembler>::
+    setup(const GridType& grid,
+        const Assembler* assembler,
+        const SolutionType& x,
+        const Dune::BitSetVector<blocksize>& dirichletNodes,
+        const Dune::ParameterTree& parameterSet)
+{
+    setup(grid,
+          assembler,
+          x,
+          dirichletNodes,
+          parameterSet.get<double>("tolerance"),
+          parameterSet.get<int>("maxProximalNewtonSteps "),
+          parameterSet.get<double>("initialRegularization"),
+          parameterSet.get<bool>("instrumented", 0));
+}
+
 template <class Basis, class TargetSpace, class Assembler>
 void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
@@ -52,7 +70,7 @@ setup(const GridType& grid,
     //////////////////////////////////////////////////////////////////
     //  Create global numbering for matrix and vector transfer
     //////////////////////////////////////////////////////////////////
-
+    
     globalMapper_ = std::make_unique<GlobalMapper>(grid_->leafGridView());
     // Transfer all Dirichlet data to the master processor
     VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, Dune::BitSetVector<blocksize> > vectorComm(*globalMapper_,
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index e09ed7764ba57fcdcfcd851e407e374ece56375c..915de346baec9100880cd1a1fe31c64904c86dbb 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -4,6 +4,7 @@
 #include <vector>
 
 #include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
@@ -81,6 +82,13 @@ public:
                double initialRegularization,
                bool instrumented);
 
+    /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+    void setup(const GridType& grid,
+               const Assembler* assembler,
+               const SolutionType& x,
+               const Dune::BitSetVector<blocksize>& dirichletNodes,
+               const Dune::ParameterTree& parameterSet);
+
     void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
     {
       scaling_ = scaling;
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index d5bdce74ab503275b860496b21f22fe1c7a8746f..7678a442dfc09472ce7c813e68b987e23b394bf1 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -28,6 +28,29 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 #endif
 
+template <class Basis, class TargetSpace, class Assembler>
+void RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
+    setup(const GridType& grid,
+        const Assembler* assembler,
+        const SolutionType& x,
+        const Dune::BitSetVector<blocksize>& dirichletNodes,
+        const Dune::ParameterTree& parameterSet)
+{
+    setup(grid,
+          assembler,
+          x,
+          dirichletNodes,
+          parameterSet.get<double>("tolerance"),
+          parameterSet.get<int>("maxTrustRegionSteps"),
+          parameterSet.get<double>("initialTrustRegionRadius"),
+          parameterSet.get<int>("numIt"),
+          parameterSet.get<double>("mgTolerance"),
+          parameterSet.get<int>("mu"), parameterSet.get<int>("nu1"), parameterSet.get<int>("nu2"),
+          parameterSet.get<int>("baseIt"),
+          parameterSet.get<double>("baseTolerance"),
+          parameterSet.get<bool>("instrumented", 0));
+}
+
 template <class Basis, class TargetSpace, class Assembler>
 void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
 setup(const GridType& grid,
@@ -692,7 +715,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
                 fwrite(&x_[j], sizeof(TargetSpace), 1, fpIterate);
 
             fclose(fpIterate);
-
+            
         }
 
         if (rank==0)
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 6e88432f166b1c2ff9fae4e07d46a82f4b05ae33..61df18e603de00f93cd808428ee37ab81c08910c 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -4,6 +4,7 @@
 #include <vector>
 
 #include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
@@ -119,6 +120,13 @@ public:
                double baseTolerance,
                bool instrumented);
 
+    /** \brief Set up the solver using a monotone multigrid method as the inner solver */
+    void setup(const GridType& grid,
+               const Assembler* assembler,
+               const SolutionType& x,
+               const Dune::BitSetVector<blocksize>& dirichletNodes,
+               const Dune::ParameterTree& parameterSet);
+
     void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
     {
       scaling_ = scaling;