diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index 8addc554443d474c48b85f7377b9512f369fb8fd..71a968c2b864435d892f96bad2f8b2a052baf1dc 100644
--- a/dumux/common/CMakeLists.txt
+++ b/dumux/common/CMakeLists.txt
@@ -11,7 +11,6 @@ boundingboxtree.hh
defaultusagemessage.hh
dimensionlessnumbers.hh
dumuxmessage.hh
-eigenvalues.hh
entitymap.hh
exceptions.hh
fixedlengthspline_.hh
@@ -34,7 +33,6 @@ staggeredfvproblem.hh
start.hh
tabulated2dfunction.hh
timemanager.hh
-timesteppingscheme.hh
valgrind.hh
variablelengthspline_.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common)
diff --git a/dumux/common/eigenvalues.hh b/dumux/common/eigenvalues.hh
deleted file mode 100644
index 11f4de1806bded8bf98a87c53d4ef0adae49ee1c..0000000000000000000000000000000000000000
--- a/dumux/common/eigenvalues.hh
+++ /dev/null
@@ -1,219 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Common
- * \brief Defines a functions to calculate eigenvalues and eigenvectors of n x n matrices. For n > 2 a cyclic jacobi method is used.
- * This implementation is not efficient for larg matrices!
- */
-#ifndef DUMUX_EIGENVALUES_HH
-#define DUMUX_EIGENVALUES_HH
-
-#include <algorithm>
-#include <cmath>
-
-#include "math.hh"
-
-namespace Dumux {
-
-template<int dim, class Matrix>
-void identityMatrix(Matrix& matrix)
-{
- for (int i = 0; i < dim; i++)
- for (int j = 0; j < dim; j++)
- matrix[i][j] = (i == j) ? 1.0 : 0.0;
-}
-
-template<int dim, class Matrix>
-double calcOffDiagonalNorm(Matrix& matrix)
-{
- double norm = 0;
- for (int i = 0; i < dim; i++)
- for (int j = 0; j < dim; j++)
- if (i != j)
- {
- norm += matrix[i][j] * matrix[i][j];
- }
-
- using std::sqrt;
- return sqrt(norm);
-}
-
-/*!
- * \brief Function to calculate eigenvalues of n x n matrices
- *
- * \param eigVel Vector for storing the eigenvalues
- * \param matrix n x n matrices for which eigenvalues have to be calculated
- * \param relativeTolerance tolerance for the relative convergence criterion (default: 0.01)
- */
-template<int dim, class EVVectorType, class MatrixType>
-bool calculateEigenValues(EVVectorType &eigVel, MatrixType& matrix, double relativeTolerance = 0.01)
-{
- if (dim == 2)
- {
- eigVel = 0;
-
- double b = -(matrix[0][0] + matrix[1][1]);
- double c = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
-
- using std::sqrt;
- eigVel[0] = (-b + sqrt(b * b - 4.0 * c)) / 2.0;
- eigVel[1] = (-b - sqrt(b * b - 4.0 * c)) / 2.0;
-
- using std::isnan;
- using std::isinf;
- if (isnan(eigVel[0]) || isinf(eigVel[0]))
- return false;
-
- if (isnan(eigVel[1]) || isinf(eigVel[1]))
- return false;
-
- return true;
- }
- else if (dim > 2)
- {
- int maxIter = 100;
- int iter = 0;
- double matrixNorm = matrix.frobenius_norm();
- double offDiagonalNorm = calcOffDiagonalNorm<dim>(matrix);
- MatrixType rotationMatrix(0.0);
- MatrixType evMatrix(matrix);
-
- while (iter < maxIter && offDiagonalNorm > relativeTolerance * matrixNorm)
- {
- for (int i = 0; i < dim - 1; i++)
- {
- for (int j = i + 1; j < dim; j++)
- {
- identityMatrix<dim>(rotationMatrix);
-
- double theta = (evMatrix[i][i] - evMatrix[j][j])
- / (2 * evMatrix[i][j]);
- using std::abs;
- using std::sqrt;
- double t = sign(theta) / (abs(theta) + sqrt(1 + theta * theta));
- double c = 1 / sqrt(1 + t * t);
- double s = c * t;
-
- rotationMatrix[i][i] = c;
- rotationMatrix[j][j] = c;
- rotationMatrix[i][j] = s;
- rotationMatrix[j][i] = -s;
-
- evMatrix.leftmultiply(rotationMatrix);
-
- rotationMatrix[i][j] = -s;
- rotationMatrix[j][i] = s;
-
- evMatrix.rightmultiply(rotationMatrix);
- }
- }
- matrixNorm = evMatrix.frobenius_norm();
- offDiagonalNorm = calcOffDiagonalNorm<dim>(evMatrix);
- iter++;
- }
-
- for (int i = 0; i < dim; i++)
- {
- eigVel[i] = evMatrix[i][i];
- using std::isinf;
- using std::isnan;
- if (isnan(eigVel[i]) || isinf(eigVel[i]))
- return false;
- }
-
- return true;
- }
-
- return false;
-}
-
-//! Function to calculate eigenvalues and eigenvectors of n x n matrices
-/*
- * \param eigVel Vector for storing the eigenvalues
- * \param eigVec n x n for storing the eigenvectors
- * \param matrix n x n matrices for which eigenvalues have to be calculated
- * \param relativeTolerance tolerance for the relative convergence criterion (default: 0.01)
- */
-template<int dim, class EVVectorType, class MatrixType>
-bool calculateEigenValues(EVVectorType &eigVel, MatrixType& eigVec, MatrixType& matrix, double relativeTolerance = 0.01)
-{
- int maxIter = 100;
- int iter = 0;
- double matrixNorm = matrix.frobenius_norm();
- double offDiagonalNorm = calcOffDiagonalNorm<dim>(matrix);
- MatrixType rotationMatrix(0.0);
- MatrixType evMatrix(matrix);
- MatrixType evecMatrix(matrix);
- identityMatrix<dim>(evecMatrix);
-
- while (iter < maxIter && offDiagonalNorm > relativeTolerance * matrixNorm)
- {
- for (int i = 0; i < dim - 1; i++)
- {
- for (int j = i + 1; j < dim; j++)
- {
- identityMatrix<dim>(rotationMatrix);
-
- double theta = (evMatrix[i][i] - evMatrix[j][j])
- / (2 * evMatrix[i][j]);
- using std::abs;
- using std::sqrt;
- double t = sign(theta) / (abs(theta) + sqrt(1 + theta * theta));
- double c = 1 / sqrt(1 + t * t);
- double s = c * t;
-
- rotationMatrix[i][i] = c;
- rotationMatrix[j][j] = c;
- rotationMatrix[i][j] = s;
- rotationMatrix[j][i] = -s;
-
- evMatrix.leftmultiply(rotationMatrix);
-
- rotationMatrix[i][j] = -s;
- rotationMatrix[j][i] = s;
-
- evMatrix.rightmultiply(rotationMatrix);
- evecMatrix.rightmultiply(rotationMatrix);
- }
- }
- matrixNorm = evMatrix.frobenius_norm();
- offDiagonalNorm = calcOffDiagonalNorm<dim>(evMatrix);
- iter++;
- }
-
- for (int i = 0; i < dim; i++)
- {
- eigVel[i] = evMatrix[i][i];
- using std::isinf;
- using std::isnan;
- if (isnan(eigVel[i]) || isinf(eigVel[i]))
- return false;
- for (int j = 0; j < dim; j++)
- {
- eigVec[i][j] = evecMatrix[j][i];
- }
- }
-
- return true;
-}
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/material/fluidmatrixinteractions/CMakeLists.txt b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
index 3b4e7d43340bcb73cea04253db035a45bb7aeae0..72b92a74513e4ca92a22a44a448032e5558bc3ce 100644
--- a/dumux/material/fluidmatrixinteractions/CMakeLists.txt
+++ b/dumux/material/fluidmatrixinteractions/CMakeLists.txt
@@ -10,6 +10,5 @@ install(FILES
diffusivityconstanttortuosity.hh
diffusivitymillingtonquirk.hh
permeabilitykozenycarman.hh
-permeabilityrutqvisttsang.hh
porosityprecipitation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/material/fluidmatrixinteractions)
diff --git a/dumux/material/fluidmatrixinteractions/permeabilityrutqvisttsang.hh b/dumux/material/fluidmatrixinteractions/permeabilityrutqvisttsang.hh
deleted file mode 100644
index e98e7a9904703dfc9ee9c56d1000e8099f067536..0000000000000000000000000000000000000000
--- a/dumux/material/fluidmatrixinteractions/permeabilityrutqvisttsang.hh
+++ /dev/null
@@ -1,93 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Fluidmatrixinteractions
- * \brief Relation for the effective permeability
- */
-#ifndef PERMEABILITYRUTQVISTSTANG_HH
-#define PERMEABILITYRUTQVISTSTANG_HH
-
-#include <algorithm>
-#include <cmath>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup Fluidmatrixinteractions
- * \brief Relation for the effective permeability
- *
- * After Rutqvist and Tsang (2002) \cite rutqvist2002, the effective permeability can be
- * calculated from the intrinsic permeability \f$ k_{\text{0}} \f$, the initial porosity
- * \f$ \phi_{\text{0}} \f$ and the effective porosity \f$ \phi_{\text{eff}} \f$ as used
- * in Darcis, M.: Coupling Models of Different Complexity for the Simulation of CO2
- * Storage in Deep Saline Aquifers, PhD thesis \cite darcis2013.
- * The calculation is shown here for one diagonal entry of the permeability tensor:
- *
- * \f$ k_\text{eff} = k_{\text{0}} \cdot \text{exp} \left[22.2\left(\phi_\text{eff}/ \phi_\text{0} -1 \right) \right] \f$
- *
- */
-template<class Scalar>
-class PermeabilityRutqvistTsang
-{
-public:
- /*!
- * \brief effective permeability tensor \f$\mathrm{[m^{2})]}\f$ after Rutqvist and Tsang (2002) \cite rutqvist2002 <BR>
- *
- * \param volVars volume variables
- * \param spatialParams spatial parameters
- * \param element element (to be passed to spatialParams)
- * \param fvGeometry fvGeometry (to be passed to spatialParams)
- * \param scvIdx control volumne
- *
- * \return effective permeability tensor \f$\mathrm{[m^{2})]}\f$ after Rutqvist and Tsang (2002) \cite rutqvist2002 <BR>
- *
- * This calculates the entries of effective permeability tensor from the intrinsic
- * permeability, the initial and the effective porosity as used in Darcis, M.:
- * Coupling Models of Different Complexity for the Simulation of CO2 Storage in
- * Deep Saline Aquifers, PhD thesis \cite darcis2013 .
- */
- template<class VolumeVariables, class SpatialParams, class Element, class FVGeometry>
- static auto effectivePermeability(const VolumeVariables& volVars,
- const SpatialParams& spatialParams,
- const Element& element,
- const FVGeometry& fvGeometry,
- int scvIdx)
- {
- Scalar effPorosity = volVars.effPorosity;
- Scalar initialPorosity = spatialParams.porosity(element, fvGeometry, scvIdx);
-
- // evaluate effective permeabilities for nodes i and j based on the
- // effective porosities, the initial porosities and the initial permeabilities
- Scalar exponent =
- 22.2
- * (effPorosity / initialPorosity - 1);
-
- auto Keff
- = spatialParams.intrinsicPermeability(element, fvGeometry, scvIdx);
- using std::exp;
- Keff *= exp(exponent);
-
- return Keff;
- }
-
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/CMakeLists.txt b/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/CMakeLists.txt
index 1d7b40451fa4a19e8e6c918d5c828f51c3b880fd..e341057acd1012d640e52f2326d02d0b84f0a6bf 100644
--- a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/CMakeLists.txt
+++ b/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/CMakeLists.txt
@@ -4,6 +4,5 @@ install(FILES
pressure.hh
pressureproperties.hh
pressurevelocity.hh
-pressurevelocityproperties.hh
velocity.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered)
diff --git a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh b/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh
deleted file mode 100644
index 2c65f64449e60cde49153729acac96acec553e76..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh
+++ /dev/null
@@ -1,80 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-
-/*!
- * \file
- * \ingroup SequentialOnePModel
- * \brief Defines the properties required for finite volume pressure models
- */
-
-#ifndef DUMUX_FVPRESSUREVELOCITYPORPERTIES1P_SEQUENTIAL_HH
-#define DUMUX_FVPRESSUREVELOCITYPORPERTIES1P_SEQUENTIAL_HH
-
-//Dumux-includes
-#include <dumux/porousmediumflow/1p/sequential/diffusion/properties.hh>
-
-namespace Dumux
-{
-
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-
-
-////////////////////////////////
-// properties
-////////////////////////////////
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for the one-phase problems using a standard finite volume model
-NEW_TYPE_TAG(FVPressureVelocityOneP, INHERITS_FROM(PressureOneP));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-}
-}
-
-#include "velocity.hh"
-#include "pressurevelocity.hh"
-
-namespace Dumux
-{
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-//! Set velocity reconstruction implementation standard cell centered finite volume schemes as default
-SET_TYPE_PROP( FVPressureVelocityOneP, Velocity, FVVelocity1P<TypeTag> );
-//! Set finite volume implementation of the one-phase pressure equation as default pressure model
-SET_TYPE_PROP(FVPressureVelocityOneP, PressureModel, FVPressureVelocity1P<TypeTag>);
-//! Allow assembling algorithm for the pressure matrix to assemble only from one side of a cell-cell interface
-SET_BOOL_PROP(FVPressureVelocityOneP, VisitFacesOnlyOnce, true);
-// \}
-}
-
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/CMakeLists.txt b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/CMakeLists.txt
index 49859e051379944558a149c17a682b893415817c..7f5a27a05d2a629832239bb9eac751e9e37c3502 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/CMakeLists.txt
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/CMakeLists.txt
@@ -6,7 +6,6 @@ pressure.hh
pressurepropertiesadaptive.hh
pressureproperties.hh
pressurevelocity.hh
-pressurevelocityproperties.hh
velocityadaptive.hh
velocity.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered)
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh
deleted file mode 100644
index 8df6a145d3d19bd740cb0a67f4feab8ca99195b5..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressurevelocityproperties.hh
+++ /dev/null
@@ -1,79 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup SequentialTwoPModel
- * \brief Defines the properties required for finite volume pressure models
- */
-
-#ifndef DUMUX_FVPRESSUREVELOCITYPORPERTIES2P_SEQUENTIAL_HH
-#define DUMUX_FVPRESSUREVELOCITYPORPERTIES2P_SEQUENTIAL_HH
-
-//Dumux-includes
-#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
-
-namespace Dumux
-{
-
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-
-
-////////////////////////////////
-// properties
-////////////////////////////////
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for the one-phase problems using a standard finite volume model
-NEW_TYPE_TAG(FVPressureVelocityTwoP, INHERITS_FROM(PressureTwoP));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-}
-}
-
-#include "velocity.hh"
-#include "pressurevelocity.hh"
-
-namespace Dumux
-{
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-//! Set velocity reconstruction implementation standard cell centered finite volume schemes as default
-SET_TYPE_PROP( FVPressureVelocityTwoP, Velocity, FVVelocity2P<TypeTag> );
-//! Set finite volume implementation of the one-phase pressure equation as default pressure model
-SET_TYPE_PROP(FVPressureVelocityTwoP, PressureModel, FVPressureVelocity2P<TypeTag>);
-//! Allow assembling algorithm for the pressure matrix to assemble only from one side of a cell-cell interface
-SET_BOOL_PROP(FVPressureVelocityTwoP, VisitFacesOnlyOnce, true);
-// \}
-}
-
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p/sequential/impes/CMakeLists.txt b/dumux/porousmediumflow/2p/sequential/impes/CMakeLists.txt
index a146a7a8a49a6a19e64c3f4994be57496cdd9f97..40662e9bf206198788f794dd31ca69dd7a39c287 100644
--- a/dumux/porousmediumflow/2p/sequential/impes/CMakeLists.txt
+++ b/dumux/porousmediumflow/2p/sequential/impes/CMakeLists.txt
@@ -2,7 +2,6 @@
#install headers
install(FILES
gridadaptionindicator.hh
-gridadaptionindicatorlocalflux.hh
gridadaptionindicatorlocal.hh
problem.hh
propertiesadaptive.hh
diff --git a/dumux/porousmediumflow/2p/sequential/impes/gridadaptionindicatorlocalflux.hh b/dumux/porousmediumflow/2p/sequential/impes/gridadaptionindicatorlocalflux.hh
deleted file mode 100644
index 2557018bb9af1f63d703a7c549bc1ae54c4d42a1..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p/sequential/impes/gridadaptionindicatorlocalflux.hh
+++ /dev/null
@@ -1,693 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup SequentialTwoPModel
- * \brief Class defining a standard, saturation dependent indicator for grid adaption
- */
-#ifndef DUMUX_GRIDADAPTIONINDICATOR2PLOCALFLUX_HH
-#define DUMUX_GRIDADAPTIONINDICATOR2PLOCALFLUX_HH
-
-#include <dune/common/version.hh>
-
-#include <dumux/porousmediumflow/sequential/impetproperties.hh>
-#include <dumux/porousmediumflow/2p/sequential/properties.hh>
-
-namespace Dumux
-{
-
-/*!
- * \brief Class defining a standard, saturation dependent indicator for grid adaption
- * \ingroup SequentialTwoPModel
- *
- * \tparam TypeTag The problem TypeTag
- */
-template<class TypeTag>
-class GridAdaptionIndicator2PLocalFlux
-{
-private:
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Element = typename GridView::Traits::template Codim<0>::Entity;
-
- using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- using IndexSet = typename Grid::LevelGridView::IndexSet;
-
- using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
- using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
-
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
-
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
- enum
- {
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
- };
-
- enum
- {
- Sw = Indices::saturationW,
- Sn = Indices::saturationNw,
- eqIdxSat = Indices::satEqIdx
- };
- enum
- {
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
- };
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
- using SetVector = Dune::FieldVector<Scalar, 2>;
-
- struct SetField {
- Scalar indicator;
- Scalar volume;
- int index;
- SetField()
- :indicator(0),volume(0), index(0)
- {}
- };
-
- struct Comparison {
- bool operator() (const SetField& lhs, const SetField& rhs) const
- {return lhs.indicator<rhs.indicator;}
- };
-
- using RangeSet = std::set<SetField, Comparison>;
-
-public:
- /*!
- * \brief Calculates the indicator used for refinement/coarsening for each grid cell.
- *
- * This standard indicator is based on the saturation gradient.
- */
- void calculateIndicator()
- {
-
- int size = problem_.variables().cellDataGlobal().size();
-
- // prepare an indicator for refinement
- indicatorVector_.resize(size);
-
- if (useSatInd_ || usePercentileSat_)
- indicatorVectorSat_.resize(size);
- if (useFluxInd_ || usePercentileFlux_)
- indicatorVectorFlux_.resize(size);
-
-
- RangeSet satRange;
- RangeSet fluxRange;
-
- SetField satVec;
- SetField fluxVec;
-
- Scalar totalVolume = 0;
- Scalar totalVolumeSat = 0;
-
- // 1) calculate Indicator -> min, maxvalues
- // Schleife über alle Leaf-Elemente
- for (const auto& element : elements(problem_.gridView()))
- {
- // Bestimme maximale und minimale Sättigung
- // Index des aktuellen Leaf-Elements
- int globalIdxI = problem_.variables().index(element);
-
- indicatorVector_[globalIdxI] = 0.5 * (refineBound_ + coarsenBound_);
- if (useSatInd_ || usePercentileSat_)
- indicatorVectorSat_[globalIdxI] = -1e100;
- if (useFluxInd_ || usePercentileFlux_)
- indicatorVectorFlux_[globalIdxI] = -1e100;
-
- const CellData& cellDataI = problem_.variables().cellData(globalIdxI);
-
- bool isSpecialCell = false;
-
- if (!checkPhysicalRange_(cellDataI))
- {
- indicatorVector_[globalIdxI] = refineBound_ + 1.0;
- isSpecialCell = true;
- }
-
- Scalar volume = element.geometry().volume();
- totalVolume += volume;
-
- if (refineAtSource_)
- {
- PrimaryVariables source(0.0);
- problem_.sourceAtPos(source, element.geometry().center());
- for (int i = 0; i < 2; i++)
- {
- using std::abs;
- if (abs(source[i]) > 1e-10)
- {
- indicatorVector_[globalIdxI] = refineBound_ + 1.0;
- isSpecialCell = true;
- break;
- }
- }
- }
-
- if (isSpecialCell)
- {
- continue;
- }
-
- Scalar satI = 0.0;
- Scalar satW = 0.0;
- switch (saturationType_)
- {
- case Sw:
- satI = cellDataI.saturation(wPhaseIdx);
- satW = satI;
- break;
- case Sn:
- satI = cellDataI.saturation(nPhaseIdx);
- satW = 1.0-satI;
- break;
- }
-
- const typename Element::Geometry& geometry = element.geometry();
- // get corresponding reference element
- using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
- const auto refElement = ReferenceElements::general(geometry.type());
-#else
- const auto& refElement = ReferenceElements::general(geometry.type());
-#endif
- const int numberOfFaces=refElement.size(1);
-
- std::vector<Scalar> flux(numberOfFaces,0);
-
- // Calculate refinement indicator on all cells
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- if (isSpecialCell)
- {
- break;
- }
-
- if (intersection.neighbor())
- {
- const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(intersection.outside()));
- if (!checkPhysicalRange_(cellDataJ))
- {
- indicatorVector_[globalIdxI] = refineBound_ + 1.0;
- isSpecialCell = true;
- break;
- }
- }
-
- int idxInInside = intersection.indexInInside();
-
- if (useFluxInd_ || usePercentileFlux_)
- {
- int isIndex = intersection.indexInInside();
- flux[isIndex] += (intersection.centerUnitOuterNormal()
- * cellDataI.fluxData().velocityTotal(idxInInside)) * intersection.geometry().volume();
-
- //Scalar velNorm = cellDataI.fluxData().velocityTotal(idxInInside).two_norm();
- //using std::max;
- //indicatorVectorFlux_[globalIdxI] = max(velNorm, indicatorVectorFlux_[globalIdxI]);
- }
-
- // exit, if it is not a neighbor
- if (intersection.boundary())
- {
- BoundaryTypes bcTypes;
- problem_.boundaryTypes(bcTypes, intersection);
-
- for (int i = 0; i < 2; i++)
- {
- if (bcTypes.isNeumann(i))
- {
- PrimaryVariables flux(0.0);
- problem_.neumann(flux, intersection);
-
- bool fluxBound = false;
- for (int j = 0; j < 2; j++)
- {
- using std::abs;
- if (abs(flux[j]) > 1e-10)
- {
- if (refineAtFluxBC_)
- {
- indicatorVector_[globalIdxI] = refineBound_ + 1.0;
- isSpecialCell = true;
- fluxBound = true;
- break;
- }
- }
- }
- if (fluxBound)
- break;
- }
- else if (bcTypes.isDirichlet(i))
- {
- if (refineAtDirichletBC_)
- {
- indicatorVector_[globalIdxI] = refineBound_ + 1.0;
- isSpecialCell = true;
- break;
- }
- }
- }
- if (useSatInd_ || usePercentileSat_)
- {
- Scalar satJ = satI;
- if (bcTypes.isDirichlet(eqIdxSat))
- {
- PrimaryVariables sat(0.0);
- problem_.dirichlet(sat, intersection);
- satJ = sat[eqIdxSat];
- }
- using std::abs;
- Scalar localdelta = abs(satI - satJ);
- using std::max;
- indicatorVectorSat_[globalIdxI] = max(indicatorVectorSat_[globalIdxI], localdelta);
- }
- }
- else
- {
- if (useSatInd_ || usePercentileSat_)
- {
- // get neighbors
- int globalIdxJ = problem_.variables().index(intersection.outside());
-
- Scalar satJ = 0.;
- switch (saturationType_)
- {
- case Sw:
- satJ = problem_.variables().cellData(globalIdxJ).saturation(wPhaseIdx);
- break;
- case Sn:
- satJ = problem_.variables().cellData(globalIdxJ).saturation(nPhaseIdx);
- break;
- }
-
- using std::abs;
- Scalar localdelta = abs(satI - satJ);
- using std::max;
- indicatorVectorSat_[globalIdxI] = max(indicatorVectorSat_[globalIdxI], localdelta);
- }
- }
- }
-
- if (isSpecialCell)
- {
- continue;
- }
-
- if (useFluxInd_ || usePercentileFlux_)
- {
- // calculate velocity on reference element as the Raviart-Thomas-0
- // interpolant of the fluxes
- Dune::FieldVector<Scalar, dim> refVelocity;
- // simplices
- if (refElement.type().isSimplex()) {
- for (int dimIdx = 0; dimIdx < dim; dimIdx++)
- {
- refVelocity[dimIdx] = -flux[dim - 1 - dimIdx];
- for (int fIdx = 0; fIdx < dim + 1; fIdx++)
- {
- refVelocity[dimIdx] += flux[fIdx]/(dim + 1);
- }
- }
- }
- // cubes
- else if (refElement.type().isCube()){
- for (int i = 0; i < dim; i++)
- refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
- }
- // 3D prism and pyramids
- else {
- DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
- }
-
- const DimVector& localPos = refElement.position(0, 0);
-
- // get the transposed Jacobian of the element mapping
- const DimMatrix jacobianT = geometry.jacobianTransposed(localPos);
-
- // calculate the element velocity by the Piola transformation
- DimVector elementVelocity(0);
- jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= geometry.integrationElement(localPos);
-
- Scalar velNorm = elementVelocity.two_norm();
- using std::max;
- indicatorVectorFlux_[globalIdxI] = max(velNorm, indicatorVectorFlux_[globalIdxI]);
- }
-
-
- if (useFluxInd_ || usePercentileFlux_)
- {
- fluxVec.indicator = indicatorVectorFlux_[globalIdxI];
- fluxVec.volume = volume;
- fluxVec.index = globalIdxI;
- fluxRange.insert(fluxVec);
- }
-
- if (useSatInd_ || usePercentileSat_)
- {
- satVec.indicator = indicatorVectorSat_[globalIdxI];
- if (satVec.indicator > 1e-8)
- {
- satVec.volume = volume;
- satVec.index = globalIdxI;
- satRange.insert(satVec);
-
- totalVolumeSat += volume;
- }
- }
- }
-
- if (useSatInd_ || usePercentileSat_)
- {
- int range = satRange.size();
- if (range > 0)
- {
- typename RangeSet::iterator it = satRange.begin();
- Scalar minLocalDelta = (*it).indicator;
-
- typename RangeSet::reverse_iterator rIt = satRange.rbegin();
- Scalar maxLocalDelta = (*rIt).indicator;
-
- if (maxLocalDelta > 0.)
- {
- Scalar globalDeltaDelta = maxLocalDelta - minLocalDelta;
- for (int i = 0; i < size; i++)
- {
- indicatorVectorSat_[i] -= minLocalDelta;
- indicatorVectorSat_[i] /= globalDeltaDelta;
- }
- }
-
- if (usePercentileSat_)
- {
- using std::max;
- Scalar lowerBound = max(0.0,coarsenPercentileSat_ * totalVolumeSat);
-
- Scalar accumulatedVolume = 0;
- while (accumulatedVolume <= lowerBound && it != satRange.end())
- {
- indicatorVector_[(*it).index] = coarsenBound_-1;
- accumulatedVolume += (*it).volume;
- ++it;
- }
- }
- }
- }
-
- if (useFluxInd_ || usePercentileFlux_)
- {
- int range = fluxRange.size();
- if (range > 0)
- {
- typename RangeSet::iterator it = fluxRange.begin();
- Scalar minFlux = (*it).indicator;
-
- typename RangeSet::reverse_iterator rIt = fluxRange.rbegin();
- Scalar maxFlux = (*rIt).indicator;
-
- if (maxFlux > 0.)
- {
- Scalar globalDeltaDelta = maxFlux - minFlux;
- for (int i = 0; i < size; i++)
- {
- indicatorVectorFlux_[i] -= minFlux;
- indicatorVectorFlux_[i] /= globalDeltaDelta;
- }
- }
-
- if (usePercentileFlux_)
- {
- using std::max;
- Scalar lowerBound = max(0.0,coarsenPercentileFlux_ * totalVolume);
- Scalar accumulatedVolume = 0;
- while (accumulatedVolume <= lowerBound && it != fluxRange.end())
- {
- indicatorVector_[(*it).index] = coarsenBound_-1;
- accumulatedVolume += (*it).volume;
- ++it;
- }
- }
- }
- }
-
- if (useSatInd_ || usePercentileSat_)
- {
- int range = satRange.size();
- if (range > 0)
- {
- if (usePercentileSat_)
- {
- typename RangeSet::reverse_iterator rIt = satRange.rbegin();
-
- using std::max;
- Scalar upperBound = max(0.0, refinePercentileSat_ * totalVolumeSat);
-
- Scalar accumulatedVolume = 0;
- while (accumulatedVolume <= upperBound && (*rIt).indicator > 1e-6 && rIt != satRange.rend())
- {
- indicatorVector_[(*rIt).index] = refineBound_+1;
- accumulatedVolume += (*rIt).volume;
- ++rIt;
- }
- }
-
-
- }
- }
-
- if (useFluxInd_ || usePercentileFlux_)
- {
- int range = fluxRange.size();
- if (range > 0)
- {
- typename RangeSet::reverse_iterator rIt = fluxRange.rbegin();
-
- if (usePercentileFlux_)
- {
- using std::min;
- Scalar upperBound = min(totalVolume, refinePercentileFlux_ * totalVolume);
- Scalar accumulatedVolume = 0;
- while (accumulatedVolume <= upperBound && rIt != fluxRange.rend())
- {
- indicatorVector_[(*rIt).index] = refineBound_+1;
- accumulatedVolume += (*rIt).volume;
- ++rIt;
- }
- }
- }
- }
-
- for (int idx = 0; idx < size; idx++)
- {
- if (useSatInd_ && indicatorVectorSat_[idx] > refineThresholdSat_)
- {
- indicatorVector_[idx] = refineBound_+1;
- }
- else if (useFluxInd_ && indicatorVectorFlux_[idx] > refineThresholdFlux_)
- {
- indicatorVector_[idx] = refineBound_+1;
- }
- else if (useSatInd_ && indicatorVectorSat_[idx] < coarsenThresholdSat_ && indicatorVector_[idx] < refineBound_)
- {
- indicatorVector_[idx] = coarsenBound_-1;
- }
- else if (useFluxInd_ && indicatorVectorFlux_[idx] < coarsenThresholdFlux_ && indicatorVector_[idx] < refineBound_)
- {
- indicatorVector_[idx] = coarsenBound_-1;
- }
- }
- }
-
- /*!
- * \brief Check if pressure data is in some physical range
- *
- * Returns true if the cell pressure is in some range
- *
- * \param cellData A grid element
- */
- bool checkPhysicalRange_(const CellData& cellData)
- {
- for (int j = 0; j < 2; j++)
- {
- if (cellData.pressure(j) < lowerPressureBound_ || cellData.pressure(j) > upperPressureBound_)
- return false;
- }
- return true;
- }
-
- /*!
- * \brief Set a lower and upper pressure constraint
- *
- * \param lowerPressureBound lower pressure value
- * \param upperPressureBound upper pressure value
- */
- void setPressureBounds(Scalar lowerPressureBound, Scalar upperPressureBound)
- {
- lowerPressureBound_ = lowerPressureBound;
- upperPressureBound_ = upperPressureBound;
- }
-
- /*!
- * \brief Indicator function for marking of grid cells for refinement
- *
- * Returns true if an element should be refined.
- *
- * \param element A grid element
- */
- bool refine(const Element& element)
- {
- int idx = problem_.elementMapper().index(element);
- return (indicatorVector_[idx] > refineBound_);
- }
-
- /*!
- * \brief Indicator function for marking of grid cells for coarsening
- *
- * Returns true if an element should be coarsened.
- *
- * \param element A grid element
- */
- bool coarsen(const Element& element)
- {
- int idx = problem_.elementMapper().index(element);
- return (indicatorVector_[idx] < coarsenBound_);
- }
-
- /*!
- * \brief Initializes the adaption indicator class
- */
- void init()
- {
- int size = problem_.gridView().size(0);
-
- indicatorVector_.resize(size, -1e100);
- if (useSatInd_)
- indicatorVectorSat_.resize(size, -1e100);
- if (useFluxInd_)
- indicatorVectorFlux_.resize(size, -1e100);
- }
-
- /*!
- * \brief Function for changing the indicatorVector values for refinement
- *
- * \param idx Index of cell which will be refined
- */
- void setIndicatorToRefine(int idx)
- {
- indicatorVector_[idx] = refineBound_+1;
- }
-
- /*!
- * \brief Function for changing the indicatorVector values for coarsening
- *
- * \param idx Index of cell which will be coarsen
- */
- void setIndicatorToCoarse(int idx)
- {
- indicatorVector_[idx] = coarsenBound_ - 1;
- }
-
- /*!
- * \brief Constructs a GridAdaptionIndicator instance
- *
- * This standard indicator is based on the saturation gradient.
- * It checks the local gradient compared to the maximum global gradient.
- * The indicator is compared locally to a refinement/coarsening threshold to decide whether
- * a cell should be marked for refinement or coarsening or should not be adapted.
- *
- * \param problem The problem object
- */
- GridAdaptionIndicator2PLocalFlux (Problem& problem):
- problem_(problem), lowerPressureBound_(0.0), upperPressureBound_(std::numeric_limits<Scalar>::max())
- {
- refineThresholdSat_ = getParam<Scalar>("GridAdapt.RefineThresholdSat", 0.8);
- coarsenThresholdSat_ = getParam<Scalar>("GridAdapt.CoarsenThresholdSat", 0.2);
- refineBound_ = 2.0;
- coarsenBound_ = 1.0;
- refineThresholdFlux_ = getParam<Scalar>("GridAdapt.RefineThresholdFlux", 0.8);
- coarsenThresholdFlux_ = getParam<Scalar>("GridAdapt.CoarsenThresholdFlux", 0.2);
- refinePercentileFlux_ = getParam<Scalar>("GridAdapt.RefinePercentileFlux", 0.8);
- coarsenPercentileFlux_ = getParam<Scalar>("GridAdapt.CoarsenPercentileFlux", 0.2);
- refinePercentileSat_ = getParam<Scalar>("GridAdapt.RefinePercentileSat", 0.8);
- coarsenPercentileSat_ = getParam<Scalar>("GridAdapt.CoarsenPercentileSat", 0.2);
- refineAtDirichletBC_ = getParam<bool>("GridAdapt.RefineAtDirichletBC");
- refineAtFluxBC_ = getParam<bool>("GridAdapt.RefineAtFluxBC");
- refineAtSource_ = getParam<bool>("GridAdapt.RefineAtSource");
-
- useSatInd_ = (refineThresholdSat_ < 1.0 || coarsenThresholdSat_ > 0.0);
- useFluxInd_ = (refineThresholdFlux_ < 1.0 || coarsenThresholdFlux_ > 0.0);
- usePercentileSat_ = (refinePercentileSat_ > 0.0 || coarsenPercentileSat_ > 0.0);
- usePercentileFlux_ = (refinePercentileFlux_ > 0.0 || coarsenPercentileFlux_ > 0.0);
-
- std::cout<<"--------------------------------------------\n";
- std::cout<<"Used adaption indicators:\n";
- if (useSatInd_)
- std::cout<<" - delta S\n";
- if (useFluxInd_)
- std::cout<<" - total flux\n";
- std::cout<<"\n Use percentiles:\n";
- if (usePercentileSat_)
- std::cout<<" - delta S\n";
- if (usePercentileFlux_)
- std::cout<<" - total flux\n";
- std::cout<<"--------------------------------------------\n";
- }
-
-private:
- Problem& problem_;
- Scalar lowerPressureBound_;
- Scalar upperPressureBound_;
- bool useSatInd_;
- bool useFluxInd_;
- bool usePercentileSat_;
- bool usePercentileFlux_;
- Scalar refineThresholdSat_;
- Scalar coarsenThresholdSat_;
- Scalar refineBound_;
- Scalar coarsenBound_;
- Scalar refineThresholdFlux_;
- Scalar coarsenThresholdFlux_;
- Scalar refinePercentileFlux_;
- Scalar coarsenPercentileFlux_;
- Scalar refinePercentileSat_;
- Scalar coarsenPercentileSat_;
- bool refineAtDirichletBC_;
- bool refineAtFluxBC_;
- bool refineAtSource_;
- std::vector<Scalar> indicatorVector_;
- std::vector<Scalar> indicatorVectorSat_;
- std::vector<Scalar> indicatorVectorFlux_;
- static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
- static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
-
-};
-}
-
-#endif