diff --git a/dumux/geomechanics/poroelastic/couplingmanager.hh b/dumux/geomechanics/poroelastic/couplingmanager.hh
index ad3bef3d551153fc6dc1f44f555db9608c96ea23..2253756a6c83c801f4917cea2b6aa46080e214de 100644
--- a/dumux/geomechanics/poroelastic/couplingmanager.hh
+++ b/dumux/geomechanics/poroelastic/couplingmanager.hh
@@ -66,6 +66,7 @@ class PoroMechanicsCouplingManager : public virtual CouplingManager< MDTraits >
     template<std::size_t id> using PrimaryVariables = typename GridVariables<id>::PrimaryVariables;
     template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables;
     template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView;
+    template<std::size_t id> using VolumeVariables = typename GridVolumeVariables<id>::VolumeVariables;
     template<std::size_t id> using FVGridGeometry = typename GridVariables<id>::GridGeometry;
     template<std::size_t id> using FVElementGeometry = typename FVGridGeometry<id>::LocalView;
     template<std::size_t id> using GridView = typename FVGridGeometry<id>::GridView;
@@ -138,8 +139,10 @@ public:
               std::shared_ptr< Problem<PoroMechId> > poroMechanicalProblem,
               const SolutionVector& curSol)
     {
-        // set up tuple containing the sub-problems
-        problemTuple_ = std::make_tuple(pmFlowProblem, poroMechanicalProblem);
+        // set the sub problems
+        this->setSubProblem(pmFlowProblem, pmFlowId);
+        this->setSubProblem(poroMechanicalProblem, poroMechId);
+
         // copy the solution vector
         ParentType::updateSolution(curSol);
         // set up the coupling map pmfow -> poromechanics
@@ -153,7 +156,7 @@ public:
                                                          const Element<PMFlowId>& element,
                                                          Dune::index_constant<PoroMechId> poroMechDomainId) const
     {
-        return pmFlowCouplingMap_[ problem<PMFlowId>().fvGridGeometry().elementMapper().index(element) ];
+        return pmFlowCouplingMap_[ this->problem(pmFlowId).fvGridGeometry().elementMapper().index(element) ];
     }
 
     /*!
@@ -163,7 +166,7 @@ public:
                                                           const Element<PoroMechId>& element,
                                                           Dune::index_constant<PMFlowId> pmFlowDomainId) const
     {
-        const auto eIdx = problem<PMFlowId>().fvGridGeometry().elementMapper().index(element);
+        const auto eIdx = this->problem(pmFlowId).fvGridGeometry().elementMapper().index(element);
         return CouplingStencilType<PoroMechId>{ {eIdx} };
     }
 
@@ -185,11 +188,12 @@ public:
 
         // prepare the fvGeometry and the element volume variables
         // these quantities will be used later to obtain the effective pressure
-        auto fvGeometry = localView( problem<PMFlowId>().fvGridGeometry() );
+        auto fvGeometry = localView( this->problem(pmFlowId).fvGridGeometry() );
         auto elemVolVars = localView( assembler.gridVariables(Dune::index_constant<PMFlowId>()).curGridVolVars() );
 
         fvGeometry.bindElement(element);
         elemVolVars.bindElement(element, fvGeometry, this->curSol()[Dune::index_constant<PMFlowId>()]);
+
         poroMechCouplingContext_.pmFlowElement = element;
         poroMechCouplingContext_.pmFlowFvGeometry = std::make_unique< FVElementGeometry<PMFlowId> >(fvGeometry);
         poroMechCouplingContext_.pmFlowElemVolVars = std::make_unique< ElementVolumeVariables<PMFlowId> >(elemVolVars);
@@ -346,19 +350,13 @@ public:
                                                                         poroMechLocalAssembler.elemBcTypes());
     }
 
-    //! Return a const reference to one of the problems
-    template<std::size_t id, std::enable_if_t<(id == PMFlowId || id == PoroMechId), int> = 0>
-    const Problem<id>& problem() const
-    { return *std::get<(id == PMFlowId ? 0 : 1)>(problemTuple_); }
-
-    //! Return reference to one of the problems
-    template<std::size_t id, std::enable_if_t<(id == PMFlowId || id == PoroMechId), int> = 0>
-    Problem<id>& problem()
-    { return *std::get<(id == PMFlowId ? 0 : 1)>(problemTuple_); }
-
-    //! Return the coupling context (used in mechanical sub-problem to compute effective pressure)
-    const PoroMechanicsCouplingContext& poroMechanicsCouplingContext() const
-    { return poroMechCouplingContext_; }
+    //! Return the porous medium flow variables an scv of the poromech domain couples to
+    const VolumeVariables<PMFlowId>& getPMFlowVolVars(const Element<PoroMechId>& element) const
+    {
+        //! If we do not yet have the queried object, build it first
+        const auto eIdx = this->problem(poroMechId).fvGridGeometry().elementMapper().index(element);
+        return (*poroMechCouplingContext_.pmFlowElemVolVars)[eIdx];
+    }
 
     /*!
      * \brief the solution vector of the coupled problem
@@ -376,8 +374,8 @@ private:
     void initializeCouplingMap_()
     {
         // some references for convenience
-        const auto& pmFlowGridGeom = problem<PMFlowId>().fvGridGeometry();
-        const auto& poroMechGridGeom = problem<PoroMechId>().fvGridGeometry();
+        const auto& pmFlowGridGeom = this->problem(pmFlowId).fvGridGeometry();
+        const auto& poroMechGridGeom = this->problem(poroMechId).fvGridGeometry();
 
         // make sure the two grids are really the same. Note that if the two grids
         // happen to have equal number of elements by chance, we don't detect this source of error.
@@ -414,11 +412,6 @@ private:
         }
     }
 
-    // tuple for storing pointers to the sub-problems
-    using PMFlowProblemPtr = std::shared_ptr< Problem<PMFlowId> >;
-    using PoroMechanicalProblemPtr = std::shared_ptr< Problem<PoroMechId> >;
-    std::tuple<PMFlowProblemPtr, PoroMechanicalProblemPtr> problemTuple_;
-
     // Container for storing the coupling element stencils for the pm flow domain
     std::vector< CouplingStencilType<PMFlowId> > pmFlowCouplingMap_;
 
diff --git a/test/multidomain/poromechanics/el1p/problem_poroelastic.hh b/test/multidomain/poromechanics/el1p/problem_poroelastic.hh
index 1a3859fc11d322801ac019904ea6327cdb7acb80..2194719b9ede5ef33b9b0a7cd57ecd884b3bf3b3 100644
--- a/test/multidomain/poromechanics/el1p/problem_poroelastic.hh
+++ b/test/multidomain/poromechanics/el1p/problem_poroelastic.hh
@@ -138,11 +138,9 @@ public:
     Scalar effectiveFluidDensity(const Element& element,
                                  const SubControlVolume& scv) const
     {
-        // get context from coupling manager
-        // here, we know that the flow problem uses cell-centered finite volumes,
-        // thus, we simply take the volume variables of the element and return the density
-        const auto& context = couplingManager().poroMechanicsCouplingContext();
-        return (*context.pmFlowElemVolVars)[scv.elementIndex()].density();
+        // get porous medium flow volume variables from coupling manager
+        const auto pmFlowVolVars = couplingManager().getPMFlowVolVars(element);
+        return pmFlowVolVars.density();
     }
 
     /*!
@@ -154,12 +152,9 @@ public:
                                  const ElementVolumeVariables& elemVolVars,
                                  const FluxVarsCache& fluxVarsCache) const
     {
-        // get context from coupling manager
-        // here, we know that the flow problem uses cell-centered finite volumes,
-        // thus, we simply take the volume variables of the element and return the pressure
-        const auto& context = couplingManager().poroMechanicsCouplingContext();
-        const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
-        return (*context.pmFlowElemVolVars)[eIdx].pressure();
+        // get porous medium flow volume variables from coupling manager
+        const auto pmFlowVolVars = couplingManager().getPMFlowVolVars(element);
+        return pmFlowVolVars.pressure();
     }
 
     /*!
diff --git a/test/multidomain/poromechanics/el1p/spatialparams_1p.hh b/test/multidomain/poromechanics/el1p/spatialparams_1p.hh
index cef186006820aa20ee9b5bcdf39673c80b2a450b..206f8503997a375345f3072998e3ae1b45739830 100644
--- a/test/multidomain/poromechanics/el1p/spatialparams_1p.hh
+++ b/test/multidomain/poromechanics/el1p/spatialparams_1p.hh
@@ -75,7 +75,7 @@ public:
     {
         static constexpr auto poroMechId = CouplingManager::poroMechId;
 
-        const auto& poroMechGridGeom = couplingManagerPtr_->template problem<poroMechId>().fvGridGeometry();
+        const auto& poroMechGridGeom = couplingManagerPtr_->problem(poroMechId).fvGridGeometry();
         const auto poroMechElemSol = elementSolution(element, couplingManagerPtr_->curSol()[poroMechId], poroMechGridGeom);
 
         // evaluate the deformation-dependent porosity at the scv center
diff --git a/test/multidomain/poromechanics/el2p/problem_poroelastic.hh b/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
index d4027b1500bf4b470bd98545dabc88527cc8715d..4034d3ea08b3a80cb050b4005996ee30a653f0fd 100644
--- a/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
+++ b/test/multidomain/poromechanics/el2p/problem_poroelastic.hh
@@ -139,16 +139,13 @@ public:
      */
     Scalar effectiveFluidDensity(const Element& element, const SubControlVolume& scv) const
     {
-        // get context from coupling manager
-        const auto& context = couplingManager().poroMechanicsCouplingContext();
-
-        // here, we know that the flow problem uses cell-centered finite volumes, thus,
-        // we simply take the volume variables of the scv (i.e. element) to obtain fluid properties
-        const auto& facetVolVars = (*context.pmFlowElemVolVars)[scv.elementIndex()];
-        Scalar wPhaseDensity = facetVolVars.density(FluidSystem::phase0Idx);
-        Scalar nPhaseDensity = facetVolVars.density(FluidSystem::phase1Idx);
-        Scalar Sw = facetVolVars.saturation(FluidSystem::phase0Idx);
-        Scalar Sn = facetVolVars.saturation(FluidSystem::phase1Idx);
+        // get porous medium flow volume variables from coupling manager
+        const auto pmFlowVolVars = couplingManager().getPMFlowVolVars(element);
+
+        Scalar wPhaseDensity = pmFlowVolVars.density(FluidSystem::phase0Idx);
+        Scalar nPhaseDensity = pmFlowVolVars.density(FluidSystem::phase1Idx);
+        Scalar Sw = pmFlowVolVars.saturation(FluidSystem::phase0Idx);
+        Scalar Sn = pmFlowVolVars.saturation(FluidSystem::phase1Idx);
         return (wPhaseDensity * Sw + nPhaseDensity * Sn);
     }
 
@@ -161,17 +158,13 @@ public:
                                  const ElementVolumeVariables& elemVolVars,
                                  const FluxVarsCache& fluxVarsCache) const
     {
-        // get context from coupling manager
-        const auto& context = couplingManager().poroMechanicsCouplingContext();
-
-        // here, we know that the flow problem uses cell-centered finite volumes, thus,
-        // we simply take the volume variables of the element to obtain fluid properties
-        const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
-        const auto& facetVolVars = (*context.pmFlowElemVolVars)[eIdx];
-        Scalar pw = facetVolVars.pressure(FluidSystem::phase0Idx);
-        Scalar pn = facetVolVars.pressure(FluidSystem::phase1Idx);
-        Scalar Sw = facetVolVars.saturation(FluidSystem::phase0Idx);
-        Scalar Sn = facetVolVars.saturation(FluidSystem::phase1Idx);
+        // get porous medium flow volume variables from coupling manager
+        const auto pmFlowVolVars = couplingManager().getPMFlowVolVars(element);
+
+        Scalar pw = pmFlowVolVars.pressure(FluidSystem::phase0Idx);
+        Scalar pn = pmFlowVolVars.pressure(FluidSystem::phase1Idx);
+        Scalar Sw = pmFlowVolVars.saturation(FluidSystem::phase0Idx);
+        Scalar Sn = pmFlowVolVars.saturation(FluidSystem::phase1Idx);
         return (pw * Sw + pn * Sn);
     }
 
diff --git a/test/multidomain/poromechanics/el2p/spatialparams_2p.hh b/test/multidomain/poromechanics/el2p/spatialparams_2p.hh
index ca1bb8355911a7c01eab283405ddf60c5e0b5a4a..4cba433c4cd04e0cb5c5bfc6374282705ade48d2 100644
--- a/test/multidomain/poromechanics/el2p/spatialparams_2p.hh
+++ b/test/multidomain/poromechanics/el2p/spatialparams_2p.hh
@@ -89,7 +89,7 @@ public:
     {
         static constexpr auto poroMechId = CouplingManager::poroMechId;
 
-        const auto& poroMechGridGeom = couplingManagerPtr_->template problem<poroMechId>().fvGridGeometry();
+        const auto& poroMechGridGeom = couplingManagerPtr_->problem(poroMechId).fvGridGeometry();
         const auto poroMechElemSol = elementSolution(element, couplingManagerPtr_->curSol()[poroMechId], poroMechGridGeom);
 
         // evaluate the deformation-dependent porosity at the scv center