diff --git a/CHANGELOG.md b/CHANGELOG.md
index ef266eccedca1b75f8e8a369c1dabf737c1f757a..04c6261fbf904bc4279b9753c6bd66fd9304f1d2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,9 +11,11 @@ Differences Between DuMu<sup>x</sup> 3.8 and DuMu<sup>x</sup> 3.7
 - __Pore network__: Added the dual network model proposed in Koch et al (2021) https://doi.org/10.1007/s11242-021-01602-5
 - __Pore network__: Added a model for heat conduction in a solid grain network
 - __Nonisothermal__: Added a specialized local residual for incompressible flows where the pressure work term cancels with the pressure-dependent part of the enthalpy so that only the internal energy remains in the advective heat flux.
+- __Tpfa__: The scv/scvf interfaces no longer store the geometry (can be obtained via local view of grid geometry). This means the memory footprint and also construction performance of the scv/scvfs is greatly improved, especially in higher dimensions. (Similar changes for the CVFE methods are already implemented since the previous release.)
 
 ### Immediate interface changes not allowing/requiring a deprecation period:
 
+- __Newton__: The Newton solver no longer supports linear solvers without a `norm` interface when computing the resisual norm is required. The linear solvers available in Dumux all have such an interface.
 
 ### Deprecated properties/classes/functions/files, to be removed after 3.8:
 
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index 32544f8b5aa87f56ae24784741b7939419a92a28..f2340daee9c997fc9e651e5c6ff04e466dcd6cba 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -216,53 +216,6 @@ public:
         });
     }
 
-    //! compute a residual's vector norm (this is a temporary interface introduced during the deprecation period)
-    [[deprecated("Use the linear solver's norm. Will be deleted after 3.7")]]
-    Scalar normOfResidual(ResidualType& residual) const
-    {
-        // issue a warning if the calculation is used in parallel with overlap
-        static bool warningIssued = false;
-
-        if (gridView().comm().size() > 1 && gridView().overlapSize(0) == 0)
-        {
-            if constexpr (isBox)
-            {
-                using DM = typename GridGeometry::VertexMapper;
-                using PVHelper = ParallelVectorHelper<GridView, DM, GridView::dimension>;
-
-                PVHelper vectorHelper(gridView(), gridGeometry_->vertexMapper());
-
-                vectorHelper.makeNonOverlappingConsistent(residual);
-            }
-        }
-        else if (!warningIssued)
-        {
-            if (gridView().comm().size() > 1 && gridView().comm().rank() == 0)
-                std::cout << "\nWarning: norm calculation adds entries corresponding to\n"
-                << "overlapping entities multiple times. Please use the norm\n"
-                << "function provided by a linear solver instead." << std::endl;
-
-            warningIssued = true;
-        }
-
-        // calculate the square norm of the residual
-        Scalar result2 = residual.two_norm2();
-        if (gridView().comm().size() > 1)
-            result2 = gridView().comm().sum(result2);
-
-        using std::sqrt;
-        return sqrt(result2);
-    }
-
-    //! compute the residual and return it's vector norm
-    [[deprecated("Use assembleResidual and the linear solver's norm. Will be deleted after 3.7")]]
-    Scalar residualNorm(const SolutionVector& curSol) const
-    {
-        ResidualType residual(numDofs());
-        assembleResidual(residual, curSol);
-        return normOfResidual(residual);
-    }
-
     /*!
      * \brief Tells the assembler which jacobian and residual to use.
      *        This also resizes the containers to the required sizes and sets the
diff --git a/dumux/discretization/box/elementsolution.hh b/dumux/discretization/box/elementsolution.hh
deleted file mode 100644
index d3e572eb2b5a160082752c64fad09e2988cbb0f0..0000000000000000000000000000000000000000
--- a/dumux/discretization/box/elementsolution.hh
+++ /dev/null
@@ -1,18 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup BoxDiscretization
- * \brief The local element solution class for the box method
- */
-#ifndef DUMUX_BOX_ELEMENT_SOLUTION_HH
-#define DUMUX_BOX_ELEMENT_SOLUTION_HH
-
-#include <dumux/discretization/cvfe/elementsolution.hh>
-#warning "This header is deprecated and will be removed after release 3.7. Include dumux/discretization/cvfe/elementsolution.hh"
-
-#endif
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 1372cc07076d23e469a93df7246420ef9d75067c..b4de9dd31bcb9266a688d47419e2743718b2d948 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -180,15 +180,15 @@ public:
     bool hasBoundaryScvf() const
     { return gridGeometry().hasBoundaryScvf(scvIndices_[0]); }
 
+    //! Create the geometry of a given sub control volume
+    typename Element::Geometry geometry(const SubControlVolume& scv) const
+    { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
+
     // suppress warnings due to current implementation
     // these interfaces should be used!
     #pragma GCC diagnostic push
     #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
 
-    //! Create the geometry of a given sub control volume
-    typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
-    { return scv.geometry(); }
-
     //! Create the geometry of a given sub control volume face
     typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
     { return scvf.geometry(); }
diff --git a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
index dcc2a419c012be5f4a3aa24eea5a54d897940693..d677c55da6d041b04515784f3eb6aeb4c1887dc1 100644
--- a/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh
@@ -185,14 +185,6 @@ public:
     std::size_t corners() const
     { return corners_.size(); }
 
-    //! Returns the corner for a given local index
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! Returns the global position of the vertex the scvf is connected to
     const GlobalPosition& vertexCorner() const
     { return corners_.back(); }
diff --git a/dumux/discretization/cellcentered/subcontrolvolume.hh b/dumux/discretization/cellcentered/subcontrolvolume.hh
index ee61d5ecea537341afda774db45863defb680b04..ced0fddd48ca4630c0dd8acbd88377fb5822eec2 100644
--- a/dumux/discretization/cellcentered/subcontrolvolume.hh
+++ b/dumux/discretization/cellcentered/subcontrolvolume.hh
@@ -12,10 +12,6 @@
 #ifndef DUMUX_DISCRETIZATION_CC_SUBCONTROLVOLUME_HH
 #define DUMUX_DISCRETIZATION_CC_SUBCONTROLVOLUME_HH
 
-#include <memory>
-
-#include <dune/common/fvector.hh>
-
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/subcontrolvolumebase.hh>
 
@@ -51,23 +47,9 @@ class CCSubControlVolume
 {
     using ThisType = CCSubControlVolume<GV, T>;
     using ParentType = SubControlVolumeBase<ThisType, T>;
-    using Geometry = typename T::Geometry;
     using GridIndexType = typename T::GridIndexType;
     using LocalIndexType = typename T::LocalIndexType;
     using Scalar = typename T::Scalar;
-
-    // In the following, the correct parameter type for the geometry passed to
-    // the constructor below is determined. It depends upon whether the
-    // `geometry()` method of the `Element` returns a copy or a reference.
-    // In the first case, the correct type is `Geometry&&`, while it is
-    // `Geometry` for the second case. Although returning by copy is prescribed
-    // by the Dune interface, the grid implementation CpGrid uses a const
-    // reference as of Opm 2018.04. Once this is fixed, the parameter type can
-    // be hardcoded to `Geometry&&` again.
-    using Element = typename GV::template Codim<0>::Entity;
-    using GeometryRT = decltype(std::declval<Element>().geometry());
-    static constexpr bool grtIsReference = std::is_lvalue_reference<GeometryRT>::value;
-    using GeometryParamType = std::conditional_t<grtIsReference, Geometry, Geometry&&>;
 public:
     //! export the type used for global coordinates
     using GlobalPosition = typename T::GlobalPosition;
@@ -76,32 +58,15 @@ public:
 
     CCSubControlVolume() = default;
 
-    // See the explanation above for deriving `GeometryParamType`.
-    CCSubControlVolume(GeometryParamType geometry,
+    template<class Geometry>
+    CCSubControlVolume(Geometry&& geometry,
                        GridIndexType elementIndex)
     : ParentType()
-    , geometry_(std::make_unique<Geometry>(std::move(geometry)))
-    , center_(geometry_->center())
+    , volume_(geometry.volume())
+    , center_(geometry.center())
     , elementIndex_(elementIndex)
     {}
 
-    //! The copy constructor
-    CCSubControlVolume(const CCSubControlVolume& other)
-    { deepCopy_(other); }
-
-    //! The move constructor
-    CCSubControlVolume(CCSubControlVolume&& other) = default;
-
-    //! The copy assignment operator
-    CCSubControlVolume& operator=(const CCSubControlVolume& other)
-    {
-        deepCopy_(other);
-        return *this;
-    }
-
-    //! The move assignment operator
-    CCSubControlVolume& operator=(CCSubControlVolume&& other) = default;
-
     //! The center of the sub control volume
     const GlobalPosition& center() const
     {
@@ -111,15 +76,7 @@ public:
     //! The volume of the sub control volume
     Scalar volume() const
     {
-        return geometry_->volume();
-    }
-
-    //! The geometry of the sub control volume
-    // e.g. for integration
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).")]]
-    const Geometry& geometry() const
-    {
-        return *geometry_;
+        return volume_;
     }
 
     //! The index of the dof this scv is embedded in (the global index of this scv)
@@ -153,26 +110,8 @@ public:
         return elementIndex_;
     }
 
-    //! Return the corner for the given local index
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).corner(i).")]]
-    GlobalPosition corner(LocalIndexType localIdx) const
-    {
-        return geometry_->corner(localIdx);
-    }
-
 private:
-    void deepCopy_(const CCSubControlVolume& other)
-    {
-        if (other.geometry_)
-            geometry_ = std::make_unique<Geometry>(*other.geometry_);
-        else
-            geometry_.reset();
-        center_ = other.center_;
-        elementIndex_ = other.elementIndex_;
-    }
-
-    // Work around the fact that geometry is not default-constructible and not copy-assignable
-    std::unique_ptr<Geometry> geometry_;
+    Scalar volume_;
     GlobalPosition center_;
     GridIndexType elementIndex_;
 };
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
index 8cd28a1a8ced38f21f4021f94058b7cf6af1f3b9..9d5f7e5afc1907946b36af6971b7a34df8763f27 100644
--- a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
@@ -82,7 +82,6 @@ class CCTpfaSubControlVolumeFace
     using Scalar = typename T::Scalar;
     using CornerStorage = typename T::CornerStorage;
     using GridIndexStorage = typename T::GridIndexStorage;
-    using Geometry = typename T::Geometry;
     using BoundaryFlag = typename T::BoundaryFlag;
 
 public:
@@ -110,7 +109,6 @@ public:
                                const GridIndexStorage& scvIndices,
                                bool isBoundary)
     : ParentType()
-    , geomType_(isGeometry.type())
     , area_(isGeometry.volume())
     , center_(isGeometry.center())
     , unitOuterNormal_(is.centerUnitOuterNormal())
@@ -118,11 +116,7 @@ public:
     , scvIndices_(scvIndices)
     , boundary_(isBoundary)
     , boundaryFlag_{is}
-    {
-        corners_.resize(isGeometry.corners());
-        for (int i = 0; i < isGeometry.corners(); ++i)
-            corners_[i] = isGeometry.corner(i);
-    }
+    {}
 
     //! The center of the sub control volume face
     const GlobalPosition& center() const
@@ -180,21 +174,6 @@ public:
         return scvfIndex_;
     }
 
-    //! return the i-th corner of this sub control volume face
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(int i) const
-    {
-        assert(i < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[i];
-    }
-
-    //! The geometry of the sub control volume face
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
-    Geometry geometry() const
-    {
-        return Geometry(geomType_, corners_);
-    }
-
     //! Return the boundary flag
     typename BoundaryFlag::value_type boundaryFlag() const
     {
@@ -202,8 +181,6 @@ public:
     }
 
 private:
-    Dune::GeometryType geomType_;
-    CornerStorage corners_;
     Scalar area_;
     GlobalPosition center_;
     GlobalPosition unitOuterNormal_;
diff --git a/dumux/discretization/facecentered/diamond/elementsolution.hh b/dumux/discretization/facecentered/diamond/elementsolution.hh
deleted file mode 100644
index 1d74cc771be640a67058c0b73c9583e33722b4cb..0000000000000000000000000000000000000000
--- a/dumux/discretization/facecentered/diamond/elementsolution.hh
+++ /dev/null
@@ -1,17 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup DiamondDiscretization
- */
-#ifndef DUMUX_DISCRETIZATION_FACECENTERED_DIAMOND_ELEMENT_SOLUTION_HH
-#define DUMUX_DISCRETIZATION_FACECENTERED_DIAMOND_ELEMENT_SOLUTION_HH
-
-#include <dumux/discretization/cvfe/elementsolution.hh>
-#warning "This header is deprecated and will be removed after release 3.7. Include dumux/discretization/cvfe/elementsolution.hh"
-
-#endif
diff --git a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
index efc1e907e84e94fcfbcad63805c5a6cfb7f94ac9..5d4bc1024f7b6362b2058a816a614a3372d3d0eb 100644
--- a/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/facecentered/staggered/fvelementgeometry.hh
@@ -223,6 +223,7 @@ public:
     { return scvf.geometry(); }
 
     #pragma GCC diagnostic pop
+
 private:
 
     const auto& scvfIndices_() const
diff --git a/dumux/discretization/facecentered/staggered/subcontrolvolume.hh b/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
index 33ba0c3ad78e1324269fc9ad9717cd7d0fb5dd1f..a33bc9656bf6183bce4cbb277939f132c086f0f6 100644
--- a/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
+++ b/dumux/discretization/facecentered/staggered/subcontrolvolume.hh
@@ -139,13 +139,6 @@ public:
     bool boundary() const
     { return boundary_; }
 
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! The geometry of the sub control volume face
     [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).")]]
     Geometry geometry() const
diff --git a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
index 86f71d9d5d289c095dc2db857fbe48b8f369245e..9adb53bb86e59e2504a2b38880d81af5e899c237 100644
--- a/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/facecentered/staggered/subcontrolvolumeface.hh
@@ -203,13 +203,6 @@ public:
     std::int_least8_t directionSign() const
     { return outerNormalSign_; }
 
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! The geometry of the sub control volume face
     [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
     Geometry geometry() const
diff --git a/dumux/discretization/porenetwork/subcontrolvolume.hh b/dumux/discretization/porenetwork/subcontrolvolume.hh
index 5a5ecd6826d0195f4ee9a79033927d640f042ccf..59c3def0bc5a162ad9b4748c339c62b3a24993fe 100644
--- a/dumux/discretization/porenetwork/subcontrolvolume.hh
+++ b/dumux/discretization/porenetwork/subcontrolvolume.hh
@@ -36,7 +36,6 @@ struct PNMDefaultScvGeometryTraits
     using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
     using Scalar = typename Grid::ctype;
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-    using CornerStorage = std::array<GlobalPosition, 2>;
     using Geometry = Dune::AffineGeometry<Scalar, 1, dimWorld>;
 };
 
@@ -56,7 +55,6 @@ class PNMSubControlVolume
     using GridIndexType = typename T::GridIndexType;
     using LocalIndexType = typename T::LocalIndexType;
     using Scalar = typename T::Scalar;
-    using CornerStorage = typename T::CornerStorage;
     using Geometry = typename T::Geometry;
 
 public:
@@ -75,12 +73,12 @@ public:
                         GridIndexType elementIndex,
                         Corners&& corners,
                         const Scalar volume)
-    : center_((corners[0]+corners[1])/2.0),
-      corners_(std::forward<CornerStorage>(corners)),
-      volume_(volume),
-      elementIndex_(elementIndex),
-      localDofIdx_(scvIdx),
-      dofIndex_(dofIndex)
+    : center_(0.5*(corners[0]+corners[1]))
+    , dofPosition_(corners[0])
+    , volume_(volume)
+    , elementIndex_(elementIndex)
+    , localDofIdx_(scvIdx)
+    , dofIndex_(dofIndex)
     {}
 
     //! The center of the sub control volume (return pore center).
@@ -106,30 +104,15 @@ public:
 
     // The position of the dof this scv is embedded in
     const GlobalPosition& dofPosition() const
-    { return corners_[0]; }
+    { return dofPosition_; }
 
     //! The global index of the element this scv is embedded in
     GridIndexType elementIndex() const
     { return elementIndex_; }
 
-    //! Return the corner for the given local index
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).corner(i).")]]
-    const GlobalPosition& corner(LocalIndexType localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
-    //! The geometry of the sub control volume e.g. for integration
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).")]]
-    Geometry geometry() const
-    {
-        return Geometry(Dune::GeometryTypes::simplex(1), corners_);
-    }
-
 private:
     GlobalPosition center_;
-    CornerStorage corners_;
+    GlobalPosition dofPosition_;
     Scalar volume_;
     GridIndexType elementIndex_;
     LocalIndexType localDofIdx_;
diff --git a/dumux/discretization/pq1bubble/elementsolution.hh b/dumux/discretization/pq1bubble/elementsolution.hh
deleted file mode 100644
index 1b3dba9f9944f536d73ec22bb835cc9463586e39..0000000000000000000000000000000000000000
--- a/dumux/discretization/pq1bubble/elementsolution.hh
+++ /dev/null
@@ -1,18 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup PQ1BubbleDiscretization
- * \brief The local element solution class for the pq1bubble method
- */
-#ifndef DUMUX_PQ1BUBBLE_ELEMENT_SOLUTION_HH
-#define DUMUX_PQ1BUBBLE_ELEMENT_SOLUTION_HH
-
-#include <dumux/discretization/cvfe/elementsolution.hh>
-#warning "This header is deprecated and will be removed after release 3.7. Include dumux/discretization/cvfe/elementsolution.hh"
-
-#endif
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index c0cb4154cc80ccb6b61889522a8343abebb52e86..54d558a9f2ecd18d1b8061e33baadb2e3c5063a1 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -222,14 +222,6 @@ public:
         return scvfIndex_;
     }
 
-    //! The positions of the corners
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! The geometry of the sub control volume face
     [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
     const Geometry geometry() const
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
index e8bb50960d0feb850918791ca4f28603291ce634..6adf5ef48044914edb3a56eea382d3f08ec799e9 100644
--- a/dumux/discretization/staggered/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -200,14 +200,6 @@ public:
         return scvfIndex_;
     }
 
-    //! The positions of the corners
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! The geometry of the sub control volume face
     [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
     const Geometry geometry() const
diff --git a/dumux/freeflow/navierstokes/momentum/CMakeLists.txt b/dumux/freeflow/navierstokes/momentum/CMakeLists.txt
index f738353326d3fa26e98e4b57f7a1d27b0736732d..6b1c27a787085208851060aac988cb6876f7f79a 100644
--- a/dumux/freeflow/navierstokes/momentum/CMakeLists.txt
+++ b/dumux/freeflow/navierstokes/momentum/CMakeLists.txt
@@ -2,8 +2,6 @@
 # SPDX-License-Identifier: GPL-3.0-or-later
 
 add_subdirectory(cvfe)
-add_subdirectory(diamond)
-add_subdirectory(pq1bubble)
 
 file(GLOB DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_HEADERS *.hh *.inc)
 install(FILES ${DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_HEADERS}
diff --git a/dumux/freeflow/navierstokes/momentum/diamond/CMakeLists.txt b/dumux/freeflow/navierstokes/momentum/diamond/CMakeLists.txt
deleted file mode 100644
index 1cb11d599a9083428e1b5c6c02ed7d6b0f5a8931..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/diamond/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-# SPDX-License-Identifier: GPL-3.0-or-later
-
-file(GLOB DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_DIAMOND_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_DIAMOND_HEADERS}
-        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/momentum/diamond)
diff --git a/dumux/freeflow/navierstokes/momentum/diamond/indices.hh b/dumux/freeflow/navierstokes/momentum/diamond/indices.hh
deleted file mode 100644
index e2c57c36ade7b84858a2dad55c8f3f7d020d8d5a..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/diamond/indices.hh
+++ /dev/null
@@ -1,26 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- * \copydoc Dumux::NavierStokesIndices
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_INDICES_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_INDICES_HH
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/indices.hh>
-
-namespace Dumux {
-
-template <int dimension>
-using NavierStokesMomentumDiamondIndices [[deprecated("Use NavierStokesMomentumCVFEIndices")]] = NavierStokesMomentumCVFEIndices<dimension>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/diamond/localresidual.hh b/dumux/freeflow/navierstokes/momentum/diamond/localresidual.hh
deleted file mode 100644
index cb99104abea719227caf47b9cbcb85eed777ce7b..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/diamond/localresidual.hh
+++ /dev/null
@@ -1,26 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- * \copydoc Dumux::NavierStokesResidualImpl
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_LOCAL_RESIDUAL_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_LOCAL_RESIDUAL_HH
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/localresidual.hh>
-
-namespace Dumux {
-
-template<class TypeTag>
-using NavierStokesMomentumDiamondResidual [[deprecated("Use NavierStokesMomentumCVFELocalResidual")]] = NavierStokesMomentumCVFELocalResidual<TypeTag>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/diamond/model.hh b/dumux/freeflow/navierstokes/momentum/diamond/model.hh
deleted file mode 100644
index a53a918598b6dc8a904bde23180aeddad92585dc..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/diamond/model.hh
+++ /dev/null
@@ -1,53 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- *
- * \brief A single-phase, isothermal Navier-Stokes model
- *
- * This model implements a single-phase, isothermal Navier-Stokes model, solving the <B> momentum balance equation </B>
- * \f[
- \frac{\partial (\varrho \textbf{v})}{\partial t} + \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}}) = \nabla \cdot (\mu (\nabla \textbf{v} + \nabla \textbf{v}^{\text{T}}))
-   - \nabla p + \varrho \textbf{g} - \textbf{f}
- * \f]
- * By setting the runtime parameter <code>Problem.EnableInertiaTerms</code> to <code>false</code> the Stokes
- * equation can be solved. In this case the term
- * \f[
- *    \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}})
- * \f]
- * is neglected.
- *
- * The <B> mass balance equation </B>
- * \f[
-       \frac{\partial \varrho}{\partial t} + \nabla \cdot (\varrho \textbf{v}) - q = 0
- * \f]
- *
- * closes the system.
- */
-
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_MODEL_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_MODEL_HH
-
-#warning "This file is deprecated and will be removed after 3.7. Use NavierStokesMomentumCVFE type tag."
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/model.hh>
-
-///////////////////////////////////////////////////////////////////////////
-// properties for the single-phase Navier-Stokes model
-///////////////////////////////////////////////////////////////////////////
-namespace Dumux::Properties {
-
-// Create new type tags
-namespace TTag {
-//! The type tag for the single-phase, isothermal Navier-Stokes model
-struct NavierStokesMomentumDiamond { using InheritsFrom = std::tuple<NavierStokesMomentumCVFE>; };
-} // end namespace TTag
-
-} // end namespace Dumux::Properties
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/diamond/volumevariables.hh b/dumux/freeflow/navierstokes/momentum/diamond/volumevariables.hh
deleted file mode 100644
index 6e2b986f5bc75cc8613ddaa2a08f5cad43cd33f4..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/diamond/volumevariables.hh
+++ /dev/null
@@ -1,27 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- *
- * \copydoc Dumux::NavierStokesVolumeVariables
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_VOLUME_VARIABLES_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_DIAMOND_VOLUME_VARIABLES_HH
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/volumevariables.hh>
-
-namespace Dumux {
-
-template<class TypeTag>
-using NavierStokesMomentumDiamondVolumeVariables [[deprecated("Use NavierStokesMomentumCVFEVolumeVariables")]] = NavierStokesMomentumCVFEVolumeVariables<TypeTag>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/pq1bubble/CMakeLists.txt b/dumux/freeflow/navierstokes/momentum/pq1bubble/CMakeLists.txt
deleted file mode 100644
index f78d8b8e9e1cae7fb8a8606638cb3826c99fa063..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/pq1bubble/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-# SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-# SPDX-License-Identifier: GPL-3.0-or-later
-
-file(GLOB DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_HEADERS}
-        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/navierstokes/momentum/pq1bubble)
diff --git a/dumux/freeflow/navierstokes/momentum/pq1bubble/indices.hh b/dumux/freeflow/navierstokes/momentum/pq1bubble/indices.hh
deleted file mode 100644
index 305419e50d3f437d94144cf0e52aec582cb7fdb7..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/pq1bubble/indices.hh
+++ /dev/null
@@ -1,27 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- * \copydoc Dumux::NavierStokesIndices
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_INDICES_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_INDICES_HH
-
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/indices.hh>
-
-namespace Dumux {
-
-template <int dimension>
-using NavierStokesMomentumPQ1BubbleIndices [[deprecated("Use NavierStokesMomentumCVFEIndices")]] = NavierStokesMomentumCVFEIndices<dimension>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/pq1bubble/localresidual.hh b/dumux/freeflow/navierstokes/momentum/pq1bubble/localresidual.hh
deleted file mode 100644
index 0933f81a61fa4422287c9a12dd0e99d29a49011e..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/pq1bubble/localresidual.hh
+++ /dev/null
@@ -1,26 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- * \copydoc Dumux::NavierStokesResidualImpl
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_LOCAL_RESIDUAL_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_LOCAL_RESIDUAL_HH
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/localresidual.hh>
-
-namespace Dumux {
-
-template<class TypeTag>
-using NavierStokesMomentumPQ1BubbleLocalResidual [[deprecated("Use NavierStokesMomentumCVFELocalResidual")]] = NavierStokesMomentumCVFELocalResidual<TypeTag>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/navierstokes/momentum/pq1bubble/model.hh b/dumux/freeflow/navierstokes/momentum/pq1bubble/model.hh
deleted file mode 100644
index c5472cc699044802de1d98414d5261c90f9b3382..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/pq1bubble/model.hh
+++ /dev/null
@@ -1,53 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- *
- * \brief A single-phase, isothermal Navier-Stokes model
- *
- * This model implements a single-phase, isothermal Navier-Stokes model, solving the <B> momentum balance equation </B>
- * \f[
- \frac{\partial (\varrho \textbf{v})}{\partial t} + \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}}) = \nabla \cdot (\mu (\nabla \textbf{v} + \nabla \textbf{v}^{\text{T}}))
-   - \nabla p + \varrho \textbf{g} - \textbf{f}
- * \f]
- * By setting the runtime parameter <code>Problem.EnableInertiaTerms</code> to <code>false</code> the Stokes
- * equation can be solved. In this case the term
- * \f[
- *    \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\text{T}})
- * \f]
- * is neglected.
- *
- * The <B> mass balance equation </B>
- * \f[
-       \frac{\partial \varrho}{\partial t} + \nabla \cdot (\varrho \textbf{v}) - q = 0
- * \f]
- *
- * closes the system.
- */
-
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_MODEL_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_MODEL_HH
-
-#warning "This file is deprecated and will be removed after 3.7. Use NavierStokesMomentumCVFE type tag."
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/model.hh>
-
-///////////////////////////////////////////////////////////////////////////
-// properties for the single-phase Navier-Stokes model
-///////////////////////////////////////////////////////////////////////////
-namespace Dumux::Properties {
-
-// Create new type tags
-namespace TTag {
-//! The type tag for the single-phase, isothermal Navier-Stokes model
-struct NavierStokesMomentumPQ1Bubble { using InheritsFrom = std::tuple<NavierStokesMomentumCVFE>; };
-} // end namespace TTag
-
-} // end namespace Dumux::Properties
-
-#endif
\ No newline at end of file
diff --git a/dumux/freeflow/navierstokes/momentum/pq1bubble/volumevariables.hh b/dumux/freeflow/navierstokes/momentum/pq1bubble/volumevariables.hh
deleted file mode 100644
index 4b5a24afefca007883b27098156bd4d9b5e9aa87..0000000000000000000000000000000000000000
--- a/dumux/freeflow/navierstokes/momentum/pq1bubble/volumevariables.hh
+++ /dev/null
@@ -1,31 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-//
-// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
-// SPDX-License-Identifier: GPL-3.0-or-later
-//
-/*!
- * \file
- * \ingroup NavierStokesModel
- *
- * \copydoc Dumux::NavierStokesVolumeVariables
- */
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_VOLUME_VARIABLES_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_PQ1BUBBLE_VOLUME_VARIABLES_HH
-
-#warning "This header is deprecated and will be removed after 3.7"
-
-#include <dumux/freeflow/navierstokes/momentum/cvfe/volumevariables.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup NavierStokesModel
- * \brief Volume variables for the single-phase Navier-Stokes model.
- */
-template <class Traits>
-using NavierStokesMomentumPQ1BubbleVolumeVariables [[deprecated("Use NavierStokesMomentumCVFEVolumeVariables")]] = NavierStokesMomentumCVFEVolumeVariables<Traits>;
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh
index 916740a31545dab5742e9b844734da1667938484..103961b4b78a1cc5d7557a0e912a07d9f2461c95 100644
--- a/dumux/linear/seqsolverbackend.hh
+++ b/dumux/linear/seqsolverbackend.hh
@@ -18,8 +18,6 @@
 
 #include <dune/istl/preconditioners.hh>
 #include <dune/istl/solvers.hh>
-#include <dune/istl/superlu.hh>
-#include <dune/istl/umfpack.hh>
 #include <dune/istl/io.hh>
 #include <dune/common/indices.hh>
 #include <dune/common/hybridutilities.hh>
@@ -55,6 +53,7 @@ class IterativePreconditionedSolverImpl
 public:
 
     template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
+    [[deprecated("Removed after 3.8. Use solver from istlsolvers.hh")]]
     static bool solve(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                       const std::string& modelParamGroup = "")
     {
@@ -75,6 +74,7 @@ public:
     }
 
     template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
+    [[deprecated("Removed after 3.8. Use solver from istlsolvers.hh")]]
     static bool solveWithGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                                const std::string& modelParamGroup = "")
     {
@@ -98,6 +98,7 @@ public:
     }
 
     template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
+    [[deprecated("Removed after 3.8. Use solver from istlsolvers.hh")]]
     static bool solveWithILU0Prec(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                                   const std::string& modelParamGroup = "")
     {
@@ -118,6 +119,7 @@ public:
 
     // solve with RestartedGMRes (needs restartGMRes as additional argument)
     template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector>
+    [[deprecated("Removed after 3.8. Use solver from istlsolvers.hh")]]
     static bool solveWithILU0PrecGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b,
                                        const std::string& modelParamGroup = "")
     {
@@ -172,565 +174,6 @@ constexpr std::size_t preconditionerBlockLevel() noexcept
     return isMultiTypeBlockMatrix<M>::value ? 2 : 1;
 }
 
-/*!
- * \ingroup Linear
- * \brief Sequential ILU(n)-preconditioned BiCSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. It can be applied to
- * nonsymmetric matrices.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: ILU(n) incomplete LU factorization. The order n can be
- * provided by the parameter LinearSolver.PreconditionerIterations and controls
- * the fill-in. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILUnBiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILUn preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential SOR-preconditioned BiCSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. It can be applied to
- * nonsymmetric matrices.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: SOR successive overrelaxation method. The relaxation is
- * controlled by the parameter LinearSolver.PreconditionerRelaxation. In each
- * preconditioning step, it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SORBiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SOR preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential SSOR-preconditioned BiCGSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. While, it can be
- * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: SSOR symmetric successive overrelaxation method. The
- * relaxation is controlled by the parameter LinearSolver.PreconditionerRelaxation.
- * In each preconditioning step, it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SSORBiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SSOR preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential GS-preconditioned BiCGSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. While, it can be
- * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: GS Gauss-Seidel method. It can be damped by the relaxation
- * parameter LinearSolver.PreconditionerRelaxation. In each preconditioning step,
- * it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] GSBiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SSOR preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential Jacobi-preconditioned BiCSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. While, it can be
- * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: Jacobi method. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation. In each preconditioning step, it is
- * applied as often as given by the parameter LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] JacBiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "Jac preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential ILU(n)-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: ILU(n) incomplete LU factorization. The order n can be
- * provided by the parameter LinearSolver.PreconditionerIterations and controls
- * the fill-in. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILUnCGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILUn preconditioned CG solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential SOR-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: SOR successive overrelaxation method. The relaxation is
- * controlled by the parameter LinearSolver.PreconditionerRelaxation. In each
- * preconditioning step, it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SORCGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SOR preconditioned CG solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential SSOR-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: SSOR symmetric successive overrelaxation method. The
- * relaxation is controlled by the parameter LinearSolver.PreconditionerRelaxation.
- * In each preconditioning step, it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SSORCGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SSOR preconditioned CG solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential GS-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: GS Gauss-Seidel method. It can be damped by the relaxation
- * parameter LinearSolver.PreconditionerRelaxation. In each preconditioning step,
- * it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] GSCGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "GS preconditioned CG solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential Jacobi-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: Jacobi method. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation. In each preconditioning step, it is
- * applied as often as given by the parameter LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] JacCGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "GS preconditioned CG solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential SSOR-preconditioned GMRes solver.
- *
- * Solver: The GMRes (generalized minimal residual) method is an iterative
- * method for the numerical solution of a nonsymmetric system of linear
- * equations.\n
- * See: Saad, Y., Schultz, M. H. (1986). "GMRES: A generalized minimal residual
- * algorithm for solving nonsymmetric linear systems." SIAM J. Sci. and Stat.
- * Comput. 7: 856–869.
- *
- * Preconditioner: SSOR symmetric successive overrelaxation method. The
- * relaxation is controlled by the parameter LinearSolver.PreconditionerRelaxation.
- * In each preconditioning step, it is applied as often as given by the parameter
- * LinearSolver.PreconditionerIterations.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SSORRestartedGMResBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::RestartedGMResSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "SSOR preconditioned GMRes solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential ILU(0)-preconditioned BiCGSTAB solver.
- *
- * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
- * faster and smoother convergence than the original BiCG. It can be applied to
- * nonsymmetric matrices.\n
- * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
- * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
- * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
- *
- * Preconditioner: ILU(0) incomplete LU factorization. The order 0 indicates
- * that no fill-in is allowed. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILU0BiCGSTABBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::BiCGSTABSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILU0 preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential ILU(0)-preconditioned CG solver.
- *
- * Solver: CG (conjugate gradient) is an iterative method for solving linear
- * systems with a symmetric, positive definite matrix.\n
- * See:  Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
- * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
- * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
- *
- * Preconditioner: ILU(0) incomplete LU factorization. The order 0 indicates
- * that no fill-in is allowed. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILU0CGBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::CGSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILU0 preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential ILU0-preconditioned GMRes solver.
- *
- * Solver: The GMRes (generalized minimal residual) method is an iterative
- * method for the numerical solution of a nonsymmetric system of linear
- * equations.\n
- * See: Saad, Y., Schultz, M. H. (1986). "GMRES: A generalized minimal residual
- * algorithm for solving nonsymmetric linear systems." SIAM J. Sci. and Stat.
- * Comput. 7: 856–869.
- *
- * Preconditioner: ILU(0) incomplete LU factorization. The order 0 indicates
- * that no fill-in is allowed. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILU0RestartedGMResBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::RestartedGMResSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solveWithILU0PrecGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILU0 preconditioned BiCGSTAB solver";
-    }
-};
-
-/*!
- * \ingroup Linear
- * \brief Sequential ILU(n)-preconditioned GMRes solver.
- *
- * Solver: The GMRes (generalized minimal residual) method is an iterative
- * method for the numerical solution of a nonsymmetric system of linear
- * equations.\n
- * See: Saad, Y., Schultz, M. H. (1986). "GMRES: A generalized minimal residual
- * algorithm for solving nonsymmetric linear systems." SIAM J. Sci. and Stat.
- * Comput. 7: 856–869.
- *
- * Preconditioner: ILU(n) incomplete LU factorization. The order n can be
- * provided by the parameter LinearSolver.PreconditionerIterations and controls
- * the fill-in. It can be damped by the relaxation parameter
- * LinearSolver.PreconditionerRelaxation.\n
- * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] ILUnRestartedGMResBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        constexpr auto precondBlockLevel = preconditionerBlockLevel<Matrix>();
-        using Preconditioner = Dune::SeqILU<Matrix, Vector, Vector, precondBlockLevel>;
-        using Solver = Dune::RestartedGMResSolver<Vector>;
-
-        return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, this->paramGroup());
-    }
-
-    std::string name() const
-    {
-        return "ILUn preconditioned GMRes solver";
-    }
-};
-
 /*!
  * \ingroup Linear
  * \brief Solver for simple block-diagonal matrices (e.g. from explicit time stepping schemes)
@@ -761,158 +204,6 @@ public:
     }
 };
 
-#if HAVE_SUPERLU
-/*!
- * \ingroup Linear
- * \brief Direct linear solver using the SuperLU library.
- *
- * See: Li, X. S. (2005). "An overview of SuperLU: Algorithms, implementation,
- * and user interface." ACM Transactions on Mathematical Software (TOMS) 31(3): 302-325.
- * http://crd-legacy.lbl.gov/~xiaoye/SuperLU/
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] SuperLUBackend : public LinearSolver
-{
-public:
-    using LinearSolver::LinearSolver;
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        static_assert(isBCRSMatrix<Matrix>::value, "SuperLU only works with BCRS matrices!");
-        using BlockType = typename Matrix::block_type;
-        static_assert(BlockType::rows == BlockType::cols, "Matrix block must be quadratic!");
-        constexpr auto blockSize = BlockType::rows;
-
-        Dune::SuperLU<Matrix> solver(A, this->verbosity() > 0);
-
-        Vector bTmp(b);
-        solver.apply(x, bTmp, result_);
-
-        int size = x.size();
-        for (int i = 0; i < size; i++)
-        {
-            for (int j = 0; j < blockSize; j++)
-            {
-                using std::isnan;
-                using std::isinf;
-                if (isnan(x[i][j]) || isinf(x[i][j]))
-                {
-                    result_.converged = false;
-                    break;
-                }
-            }
-        }
-
-        return result_.converged;
-    }
-
-    std::string name() const
-    {
-        return "SuperLU solver";
-    }
-
-    const Dune::InverseOperatorResult& result() const
-    {
-        return result_;
-    }
-
-private:
-    Dune::InverseOperatorResult result_;
-};
-#endif // HAVE_SUPERLU
-
-#if HAVE_UMFPACK
-/*!
- * \ingroup Linear
- * \brief Direct linear solver using the UMFPack library.
- *
- * See: Davis, Timothy A. (2004). "Algorithm 832". ACM Transactions on
- * Mathematical Software 30 (2): 196–199. doi:10.1145/992200.992206.
- * http://faculty.cse.tamu.edu/davis/suitesparse.html
- *
- * You can choose from one of the following ordering strategies using the input
- * parameter "LinearSolver.UMFPackOrdering":
- * \verbatim
- * 0: UMFPACK_ORDERING_CHOLMOD
- * 1: UMFPACK_ORDERING_AMD (default)
- * 2: UMFPACK_ORDERING_GIVEN
- * 3: UMFPACK_ORDERING_METIS
- * 4: UMFPACK_ORDERING_BEST
- * 5: UMFPACK_ORDERING_NONE
- * 6: UMFPACK_ORDERING_USER
- * \endverbatim
- *
- * See https://fossies.org/linux/SuiteSparse/UMFPACK/Doc/UMFPACK_UserGuide.pdf page 17 for details.
- */
-class [[deprecated("Removed after 3.7. Use solver from istlsolvers.hh")]] UMFPackBackend : public LinearSolver
-{
-public:
-
-    UMFPackBackend(const std::string& paramGroup = "")
-    : LinearSolver(paramGroup)
-    {
-        ordering_ = getParamFromGroup<int>(this->paramGroup(), "LinearSolver.UMFPackOrdering", 1);
-    }
-
-    //! set ordering strategy
-    void setOrdering(int i)
-    { ordering_ = i; }
-
-    //! the ordering strategy
-    int ordering() const
-    { return ordering_; }
-
-    template<class Matrix, class Vector>
-    bool solve(const Matrix& A, Vector& x, const Vector& b)
-    {
-        static_assert(isBCRSMatrix<Matrix>::value, "UMFPack only works with BCRS matrices!");
-        using BlockType = typename Matrix::block_type;
-        static_assert(BlockType::rows == BlockType::cols, "Matrix block must be quadratic!");
-        constexpr auto blockSize = BlockType::rows;
-
-        Dune::UMFPack<Matrix> solver;
-        solver.setVerbosity(this->verbosity() > 0);
-        solver.setOption(UMFPACK_ORDERING, ordering_);
-        solver.setMatrix(A);
-
-        Vector bTmp(b);
-        solver.apply(x, bTmp, result_);
-
-        int size = x.size();
-        for (int i = 0; i < size; i++)
-        {
-            for (int j = 0; j < blockSize; j++)
-            {
-                using std::isnan;
-                using std::isinf;
-                if (isnan(x[i][j]) || isinf(x[i][j]))
-                {
-                    result_.converged = false;
-                    break;
-                }
-            }
-        }
-
-        return result_.converged;
-    }
-
-    std::string name() const
-    {
-        return "UMFPack solver";
-    }
-
-    const Dune::InverseOperatorResult& result() const
-    {
-        return result_;
-    }
-
-private:
-    Dune::InverseOperatorResult result_;
-    int ordering_;
-};
-#endif // HAVE_UMFPACK
-
-
 /*!
  * \name Solver for MultiTypeBlockMatrix's
  */
diff --git a/dumux/multidomain/facet/box/subcontrolvolumeface.hh b/dumux/multidomain/facet/box/subcontrolvolumeface.hh
index 919e399a83b20b515222962cf367e05c9328e3e9..6f33bb60f3bdd5dcce454483d5d52e2c2db0bcec 100644
--- a/dumux/multidomain/facet/box/subcontrolvolumeface.hh
+++ b/dumux/multidomain/facet/box/subcontrolvolumeface.hh
@@ -184,14 +184,6 @@ public:
     typename BoundaryFlag::value_type boundaryFlag() const
     { return boundaryFlag_.get(); }
 
-    //! returns the position of a corner of the face
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
 private:
     // geometrical information
     CornerStorage corners_;
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 347fdc50854a4295f2dc9cc903a08e0347b591ff..de5fd2e06a5a4cc7879db628e92d5ce3068b0a9d 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -268,67 +268,6 @@ public:
         });
     }
 
-    //! compute a residual's vector norm (this is a temporary interface introduced during the deprecation period)
-    [[deprecated("Use the linear solver's norm. Will be deleted after 3.7")]]
-    Scalar normOfResidual(ResidualType& residual) const
-    {
-        // calculate the squared norm of the residual
-        Scalar resultSquared = 0.0;
-
-        // for box communicate the residual with the neighboring processes
-        using namespace Dune::Hybrid;
-        forEach(integralRange(Dune::Hybrid::size(residual)), [&](const auto domainId)
-        {
-            const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
-            const auto& gridView = gridGeometry->gridView();
-
-            if (gridView.comm().size() > 1 && gridView.overlapSize(0) == 0)
-            {
-                if constexpr (GridGeometry<domainId>::discMethod == DiscretizationMethods::box)
-                {
-                    using GV = typename GridGeometry<domainId>::GridView;
-                    using DM = typename GridGeometry<domainId>::VertexMapper;
-                    using PVHelper = ParallelVectorHelper<GV, DM, GV::dimension>;
-
-                    PVHelper vectorHelper(gridView, gridGeometry->vertexMapper());
-
-                    vectorHelper.makeNonOverlappingConsistent(residual[domainId]);
-                }
-            }
-            else if (!warningIssued_)
-            {
-                if (gridView.comm().size() > 1 && gridView.comm().rank() == 0)
-                    std::cout << "\nWarning: norm calculation adds entries corresponding to\n"
-                              << "overlapping entities multiple times. Please use the norm\n"
-                              << "function provided by a linear solver instead." << std::endl;
-
-                warningIssued_ = true;
-            }
-
-            Scalar localNormSquared = residual[domainId].two_norm2();
-
-            if (gridView.comm().size() > 1)
-            {
-                localNormSquared = gridView.comm().sum(localNormSquared);
-            }
-
-            resultSquared += localNormSquared;
-        });
-
-        using std::sqrt;
-        return sqrt(resultSquared);
-    }
-
-    //! compute the residual and return it's vector norm
-    [[deprecated("Use norm(curSol) provided by the linear solver class instead. Will be deleted after 3.7")]]
-    Scalar residualNorm(const SolutionVector& curSol)
-    {
-        ResidualType residual;
-        setResidualSize_(residual);
-        assembleResidual(residual, curSol);
-        return normOfResidual(residual);
-    }
-
     /*!
      * \brief Tells the assembler which jacobian and residual to use.
      *        This also resizes the containers to the required sizes and sets the
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 61098b492c5014819572dcd8ade35e8f9e66a09a..73c70df6075bd16f70c2e82cd89b31cc0bc0d0d6 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -84,14 +84,6 @@ struct supportsPartialReassembly
     {}
 };
 
-// helper struct and function detecting if the linear solver features a norm() function
-template <class LinearSolver, class Residual>
-using NormDetector = decltype(std::declval<LinearSolver>().norm(std::declval<Residual>()));
-
-template<class LinearSolver, class Residual>
-static constexpr bool hasNorm()
-{ return Dune::Std::is_detected<NormDetector, LinearSolver, Residual>::value; }
-
 // helpers to implement max relative shift
 template<class C> using dynamicIndexAccess = decltype(std::declval<C>()[0]);
 template<class C> using staticIndexAccess = decltype(std::declval<C>()[Dune::Indices::_0]);
@@ -173,15 +165,6 @@ void assign(To& to, const From& from)
         DUNE_THROW(Dune::Exception, "Values are not assignable to each other!");
 }
 
-template<class Residual, class LinearSolver, class Assembler>
-typename Assembler::Scalar residualNorm(Residual& residual, const LinearSolver& linearSolver, const Assembler& assembler)
-{
-    if constexpr (Detail::Newton::hasNorm<LinearSolver, Residual>())
-        return linearSolver.norm(residual);
-    else // fallback to deprecated interface
-        return assembler.normOfResidual(residual);
-}
-
 } // end namespace Dumux::Detail::Newton
 
 namespace Dumux {
@@ -497,9 +480,7 @@ public:
         try
         {
             if (numSteps_ == 0)
-                initialResidual_ = Detail::Newton::residualNorm(
-                    this->assembler().residual(), this->linearSolver(), this->assembler()
-                );
+                initialResidual_ = this->linearSolver().norm(this->assembler().residual());
 
             // solve by calling the appropriate implementation depending on whether the linear solver
             // is capable of handling MultiType matrices or not
@@ -866,9 +847,7 @@ protected:
         else
             this->assembler().assembleResidual(vars);
 
-        residualNorm_ = Detail::Newton::residualNorm(
-            this->assembler().residual(), this->linearSolver(), this->assembler()
-        );
+        residualNorm_ = this->linearSolver().norm(this->assembler().residual());
 
         reduction_ = residualNorm_/initialResidual_;
     }
diff --git a/dumux/porousmediumflow/boxdfm/geometryhelper.hh b/dumux/porousmediumflow/boxdfm/geometryhelper.hh
index 6f4464c4036a6244c024ab7f68044020263d946c..fdad1d0df476e8d9f7468c2e839544fd7c703083 100644
--- a/dumux/porousmediumflow/boxdfm/geometryhelper.hh
+++ b/dumux/porousmediumflow/boxdfm/geometryhelper.hh
@@ -135,15 +135,6 @@ public:
                                                             << " type=" << type);
     }
 
-    //! Create the sub control volume face geometries on an intersection marked as fracture
-    [[deprecated("Will be removed after release 3.6. Use other signature.")]]
-    ScvfCornerStorage getFractureScvfCorners(const Intersection& is,
-                                             const typename Intersection::Geometry& isGeom,
-                                             unsigned int edgeIndexInIntersection) const
-    {
-        return getFractureScvfCorners(is.indexInInside(), edgeIndexInIntersection);
-    }
-
     //! get fracture scvf normal vector
     typename ScvfType::Traits::GlobalPosition
     fractureNormal(const ScvfCornerStorage& scvfCorners,
diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
index e33d28065bd25d6e51a0ea9781f26535917711da..cf8504ea5a529aa60849c5398e5a70352a92797a 100644
--- a/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
+++ b/dumux/porousmediumflow/boxdfm/subcontrolvolume.hh
@@ -202,14 +202,6 @@ public:
         return Geometry(Dune::GeometryTypes::cube(dim), corners_);
     }
 
-    //! Return the corner for the given local index
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scv).corner(i).")]]
-    const GlobalPosition& corner(LocalIndexType localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
 private:
     bool isFractureScv_;
     CornerStorage corners_;
diff --git a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
index 71d19e97a9a3008d752caa5ae7b0c05b7ed9d542..85fb5af05d1c4978f5ecff9f785bcf804e986f7c 100644
--- a/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
+++ b/dumux/porousmediumflow/boxdfm/subcontrolvolumeface.hh
@@ -224,14 +224,6 @@ public:
         return static_cast<std::size_t>(!boundary());
     }
 
-    //! Returns a corner of the sub control volume face
-    [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).corner(i).")]]
-    const GlobalPosition& corner(unsigned int localIdx) const
-    {
-        assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
-        return corners_[localIdx];
-    }
-
     //! The geometry of the sub control volume face
     [[deprecated("Will be removed after 3.7. Use fvGeometry.geometry(scvf).")]]
     Geometry geometry() const
diff --git a/test/porousmediumflow/1p/compressible/instationary/assembler.hh b/test/porousmediumflow/1p/compressible/instationary/assembler.hh
index f5618f15c026000b2e7b29283569452db67e4d8c..f677cbad893988906d4a5afe130725c0af448a9d 100644
--- a/test/porousmediumflow/1p/compressible/instationary/assembler.hh
+++ b/test/porousmediumflow/1p/compressible/instationary/assembler.hh
@@ -31,10 +31,6 @@ public:
 
     void assembleResidual(const GridVariables& gridVars)
     { Assembler::assembleResidual(gridVars.dofs()); }
-
-    //! Remove residualNorm once this is fixed in the solvers !2113
-    auto residualNorm(const GridVariables& gridVars)
-    { return Assembler::residualNorm(gridVars.dofs()); }
 };
 
 } // end namespace Dumux::OnePCompressibleTest