diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index 91a5144b9a9454db3b5fb91f3f8261390348d37f..4b82b2cdb7bae10abdb64499e4fba12c662f6f50 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -19,7 +19,6 @@ add_subdirectory(nonequilibrium)
add_subdirectory(nonisothermal)
add_subdirectory(richards)
add_subdirectory(richardsnc)
-add_subdirectory(sequential)
add_subdirectory(solidenergy)
add_subdirectory(tracer)
diff --git a/dumux/porousmediumflow/sequential/CMakeLists.txt b/dumux/porousmediumflow/sequential/CMakeLists.txt
deleted file mode 100644
index 0c41110ad0faca6bb45ed8173a399a208c5a4688..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-add_subdirectory(cellcentered)
-add_subdirectory(mimetic)
-
-file(GLOB DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_HEADERS}
- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/sequential)
diff --git a/dumux/porousmediumflow/sequential/boundarytypes.hh b/dumux/porousmediumflow/sequential/boundarytypes.hh
deleted file mode 100644
index 24b43fe4ebcac0d2839fa4e1fb73e1cb5eaea8f1..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/boundarytypes.hh
+++ /dev/null
@@ -1,121 +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 3 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 Sequential
- * \brief Class to specify the type of a boundary.
- */
-#ifndef DUMUX_SEQUENTIAL_BOUNDARY_TYPES_HH
-#define DUMUX_SEQUENTIAL_BOUNDARY_TYPES_HH
-
-#include <algorithm>
-#include <array>
-
-#include <dumux/common/boundarytypes.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Sequential
- * \brief Class to specify the type of a boundary.
- */
-template <int numEq>
-class SequentialBoundaryTypes : public BoundaryTypes<numEq>
-{
-public:
- SequentialBoundaryTypes()
- : BoundaryTypes<numEq>()
- { reset(); }
-
- /*!
- * \brief Reset the boundary types for one equation.
- */
- void resetEq(int eqIdx)
- {
- BoundaryTypes<numEq>::resetEq(eqIdx);
- seqBoundaryInfo_[eqIdx].isOutflow = false;
- }
-
- /*!
- * \brief Reset the boundary types for all equations.
- *
- * After this method no equations will be disabled and neither
- * Neumann nor Dirichlet conditions will be evaluated. This
- * corresponds to a Neumann zero boundary.
- */
- void reset()
- {
- for (int eqIdx=0; eqIdx < numEq; ++eqIdx)
- resetEq(eqIdx);
- }
-
- /*!
- * \brief Set all boundary conditions to Neumann.
- */
- void setAllOutflow()
- {
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- setOutflow(eqIdx);
- }
-
- /*!
- * \brief Set a Neumann boundary condition for a single equation.
- * \param eqIdx The index of the equation on which the outflow
- * condition applies.
- */
- void setOutflow(int eqIdx)
- {
- resetEq(eqIdx);
- this->boundaryInfo_[eqIdx].visited = true;
- seqBoundaryInfo_[eqIdx].isOutflow = true;
- }
-
- /*!
- * \brief Returns true if an equation is used to specify an
- * outflow condition.
- *
- * \param eqIdx The index of the equation
- */
- bool isOutflow(unsigned eqIdx) const
- { return seqBoundaryInfo_[eqIdx].isOutflow; }
-
- /*!
- * \brief Returns true if some equation is used to specify an
- * outflow condition.
- */
- bool hasOutflow() const
- {
- return std::any_of(seqBoundaryInfo_.begin(),
- seqBoundaryInfo_.end(),
- [](const BoundaryInfo& b){ return b.isOutflow; }
- );
- }
-
-protected:
- //! use bitfields to minimize the size
- struct BoundaryInfo {
- bool isOutflow : 1;
- };
-
- std::array<BoundaryInfo, numEq> seqBoundaryInfo_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/CMakeLists.txt b/dumux/porousmediumflow/sequential/cellcentered/CMakeLists.txt
deleted file mode 100644
index f862e0b7456c7de55ea65408574158edd5c340d2..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/CMakeLists.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-add_subdirectory(mpfa)
-
-file(GLOB DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_CELLCENTERED_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_CELLCENTERED_HEADERS}
- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/sequential/cellcentered)
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/CMakeLists.txt b/dumux/porousmediumflow/sequential/cellcentered/mpfa/CMakeLists.txt
deleted file mode 100644
index f4ac89c9732d48baeae83b22d5cab21ca30e17c7..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/CMakeLists.txt
+++ /dev/null
@@ -1,3 +0,0 @@
-file(GLOB DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_CELLCENTERED_MPFA_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_CELLCENTERED_MPFA_HEADERS}
- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/sequential/cellcentered/mpfa)
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh
deleted file mode 100644
index d9808c44931a8bd24bb26c454b4a2d35d71a5795..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh
+++ /dev/null
@@ -1,430 +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 3 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/>. *
- *****************************************************************************/
-
-#ifndef DUMUX_FVMPFALINTERACTIONVOLUME_HH
-#define DUMUX_FVMPFALINTERACTIONVOLUME_HH
-
-/**
- * @file
- * @brief Class including the information of an interaction volume of a MPFA L-method that does not change with time.
- */
-
-#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
-
-namespace Dumux
-{
-
-//! \ingroup IMPET
-/*! \brief Class including the information of an interaction volume of a MPFA L-method that does not change with time.
- *
- * Includes information needed to calculate the transmissibility matrix of an L-interaction-volume.
- *
- */
-template<class TypeTag>
-class FVMPFALInteractionVolume
-{
-private:
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- enum
- {
- dim = GridView::dimension
- };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementSeed = typename Grid::template Codim<0>::EntitySeed;
-
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
-
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using FieldVectorVector = Dune::FieldVector<DimVector, dim>;
- using IndexVector = Dune::FieldVector<int, dim>;
- using BCTypeVector = std::vector<BoundaryTypes>;
- using BCVector = std::vector<PrimaryVariables>;
-
-public:
- enum FaceTypes
- {
- inside = 1,
- boundary = 0,
- outside = -1
- };
-
- //! Constructs a FVMPFALInteractionVolume object
- FVMPFALInteractionVolume(const Grid& grid)
- : grid_(&grid),
- stored_(false), normal_(FieldVectorVector(DimVector(0.0))),
- facePos_(FieldVectorVector(DimVector(0.0))),
- faceArea_(DimVector(0.0)),
- faceType_(2*dim, inside),
- indexOnElement_(IndexVector(0.0)),
- elements_(2*dim),
- centerVertexPos_(0),
- elementNum_(0)
- {
- faceIndexOnSubVolume_[0][0] = 0;
- faceIndexOnSubVolume_[0][1] = 3;
- faceIndexOnSubVolume_[1][0] = 1;
- faceIndexOnSubVolume_[1][1] = 0;
- faceIndexOnSubVolume_[2][0] = 2;
- faceIndexOnSubVolume_[2][1] = 1;
- faceIndexOnSubVolume_[3][0] = 3;
- faceIndexOnSubVolume_[3][1] = 2;
- }
-
- //! Delete stored information
- void reset()
- {
- stored_ = false;
- elements_.clear();
- elements_.resize(2*dim);
- centerVertexPos_ = 0;
- indexOnElement_ = IndexVector(0.0);
- faceType_.clear();
- faceType_.resize(2*dim, inside);
- faceArea_ = DimVector(0.0);
- facePos_ = FieldVectorVector(DimVector(0.0));
- normal_ = FieldVectorVector(DimVector(0.0));
- elementNum_ = 0;
- boundaryTypes_.clear();
- neumannValues_.clear();
- dirichletValues_.clear();
- }
-
- //! Mark storage as completed
- void setStored()
- {
- stored_ = true;
- }
-
- //! Returns true if information has already been stored
- bool isStored() const
- {
- return stored_;
- }
-
- //! Store position of the center vertex of the interaction volume
- void setCenterPosition(DimVector ¢erVertexPos)
- {
- centerVertexPos_ = centerVertexPos;
- }
-
- //! Store an element of the interaction volume
- /*!
- * \param element The element
- * \param subVolumeIdx The local element index in the interaction volume
- */
- void setSubVolumeElement(const Element& element, int subVolumeIdx)
- {
- elements_[subVolumeIdx].push_back(element.seed());
- elementNum_++;
- }
-
- //! Store position of the center of a flux face
- /*!
- * \param pos Position of face center
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setFacePosition(const DimVector& pos, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- facePos_[subVolumeIdx][subVolumeFaceIdxInInside] = pos;
- }
-
- //! Store a flux face area
- /*!
- * \param faceArea The flux face area
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setFaceArea(Scalar& faceArea, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- faceArea_[subVolumeIdx][subVolumeFaceIdxInInside] = faceArea;
- }
-
- //! Store a flux face normal
- /*!
- * \param normal The normal vector
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- normal_[subVolumeIdx][subVolumeFaceIdxInInside] = normal;
- }
-
- //! Store boundary condtion types for a flux face
- /*!
- * \param boundaryTypes BoundaryTypes object
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
- {
- if (boundaryTypes_.size() == 0)
- {
- boundaryTypes_.resize(2 * dim);
- }
- boundaryTypes_[subVolumeFaceIdx] = boundaryTypes;
- faceType_[subVolumeFaceIdx] = boundary;
- }
-
- //! Mark a flux face to be outside the model domain
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setOutsideFace(int subVolumeFaceIdx)
- {
- faceType_[subVolumeFaceIdx] = outside;
- }
-
- //! Store Dirichlet boundary condtions for a flux face
- /*!
- * \param condition Vector of primary variables
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (dirichletValues_.size() == 0)
- {
- dirichletValues_.resize(2 * dim);
- }
- dirichletValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store Neumann boundary condtions for a flux face
- /*!
- * \param condition Vector phase fluxes
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (neumannValues_.size() == 0)
- {
- neumannValues_.resize(2 * dim);
- }
- neumannValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store map from local interaction volume numbering to numbering of the Dune reference element.
- /*!
- * \param indexInInside Face index of the Dune reference element
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- indexOnElement_[subVolumeIdx][subVolumeFaceIdxInInside] = indexInInside;
- }
-
- //! Get position vector of central vertex
- DimVector& getCenterPosition()
- {
- return centerVertexPos_;
- }
-
- //! Get number of stored elements
- int getElementNumber()
- {
- return elementNum_;
- }
-
- //! Map from local interaction volume numbering to numbering of the Dune reference element.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Face index of the Dune reference element
- */
- int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! Map from local interaction volume numbering on element to numbering on interaction volume.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Local face index int the interaction volume
- */
- int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! Get an element of the interaction volume.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- *
- * \return The interaction volume sub-element.
- */
- Element getSubVolumeElement(int subVolumeIdx)
- {
- return grid_->entity(elements_[subVolumeIdx][0]);
- }
-
- //! Get boundary condtion types for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Object containing information about the boundary types.
- */
- BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
- {
- return boundaryTypes_[subVolumeFaceIdx];
- }
-
- //! Returns true if an interaction volume flux face is outside the model domain.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isOutsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == outside);
- }
-
- //! Returns true if an interaction volume flux face is inside the model domain and no boundary face.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isInsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == inside);
- }
-
- //! Returns true if the interaction volume is completely inside the model domain.
- bool isInnerVolume()
- {
- for (unsigned int i = 0; i < faceType_.size(); i++)
- {
- if (isOutsideFace(i) || isBoundaryFace(i))
- return false;
- }
- return true;
- }
-
- //! Returns true if an interaction volume flux face is a boundary face.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isBoundaryFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == boundary);
- }
-
- //! Get the Dirichlet boundary condtions for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of primary variables.
- */
- PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
- {
- return dirichletValues_[subVolumeFaceIdx];
- }
-
- //! Get the Neumann boundary condtions for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of phase fluxes.
- */
- PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
- {
- return neumannValues_[subVolumeFaceIdx];
- }
-
- //! Get a flux face normal
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return Normal vector
- */
- DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return normal_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Get position of the center of a flux face
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return Position of face center
- */
- DimVector& getFacePosition(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return facePos_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Get a flux face area
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return Area of the flux face
- */
- Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return faceArea_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Outputs stored information
- void printInteractionVolumeInfo()
- {
- std::cout<<"\nNumber of stored elements: "<<elementNum_<<"\n";
- std::cout<<" center position: "<<centerVertexPos_<<"\n";
- for (int i = 0; i < 2*dim; i++)
- {
- if (elements_[i].size() > 0)
- {
- std::cout<<"element "<<i<<":\n";
- std::cout<<"element position: "<<elements_[i][0]->geometry().center()<<"\n";
- std::cout<<"face indices on element: "<<indexOnElement_[i]<<"\n";
- std::cout<<"face normals on element: "<<normal_[i]<<"\n";
- std::cout<<"face areas on element: "<<faceArea_[i]<<"\n";
- std::cout<<"face position on element: "<<facePos_[i]<<"\n";
- std::cout<<"face type: "<<faceType_[i]<<"\n";
- }
- }
- }
-
-private:
- const Grid* grid_;
- bool stored_;
- Dune::FieldVector<FieldVectorVector, 2*dim> normal_;
- Dune::FieldVector<FieldVectorVector, 2*dim> facePos_;
- Dune::FieldVector<DimVector, 2*dim> faceArea_;
- BCTypeVector boundaryTypes_;
- std::vector<int> faceType_;
- Dune::FieldVector<IndexVector, 2*dim> indexOnElement_;
- Dune::FieldVector<IndexVector, 2*dim> faceIndexOnSubVolume_;
- std::vector<std::vector<ElementSeed> > elements_;
- BCVector neumannValues_;
- BCVector dirichletValues_;
- DimVector centerVertexPos_;
- int elementNum_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3d.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3d.hh
deleted file mode 100644
index 605e3d84ba0b8b0934a67e8b8f8253d9f4e00a4e..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3d.hh
+++ /dev/null
@@ -1,704 +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 3 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/>. *
- *****************************************************************************/
-
-#ifndef DUMUX_FVMPFAL3DINTERACTIONVOLUME_HH
-#define DUMUX_FVMPFAL3DINTERACTIONVOLUME_HH
-
-#include "properties.hh"
-
-/**
- * @file
- * @brief Class including the information of an interaction volume of a MPFA 3D method that does not change with time
- * @author Markus Wolff
- */
-
-namespace Dumux
-{
-
-//! \cond \private
-// Mapper for local interaction volume indices (see doc/docextra/3dmpfa).
-class IndexTranslator
-{
-public:
- enum
- {
- subVolumeTotalNum = 8,
- fluxFacesTotalNum = 12,
- fluxFacesNumOnSubVolume = 3,
- fluxEdgesTotalNum = 6,
- edgesNumOnFluxFace = 2
- };
-
- static int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- static int getEdgeIndexFromSubVolumeFace(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
- {
- return edgeIndexOnSubVolumeFace_[subVolumeIdx][subVolumeFaceIdx][subVolumeEdgeIdx];
- }
-
- static int getFaceIndexFromElements(int elem1Idx, int elem2Idx)
- {
- return faceIndexFromElements_[elem1Idx][elem2Idx];
- }
-
-
-private:
- static const int faceIndexOnSubVolume_[subVolumeTotalNum][fluxFacesNumOnSubVolume];
- static const int edgeIndexOnSubVolumeFace_[subVolumeTotalNum][fluxFacesNumOnSubVolume][edgesNumOnFluxFace];
- static const int faceIndexFromElements_[subVolumeTotalNum][subVolumeTotalNum];
-};
-
-const int IndexTranslator::faceIndexOnSubVolume_[subVolumeTotalNum][fluxFacesNumOnSubVolume] =
-{
- {0, 3, 8},
- {1, 0, 9},
- {3, 2, 11},
- {2, 1, 10},
- {8, 4, 7},
- {9, 5, 4},
- {11, 7, 6},
- {10, 6, 5}
-};
-
-const int IndexTranslator::edgeIndexOnSubVolumeFace_[subVolumeTotalNum][fluxFacesNumOnSubVolume][edgesNumOnFluxFace] =
-{
- {{2, 0},{0, 5},{5, 2}},
- {{3, 0},{0, 2},{2, 3}},
- {{5, 0},{0, 4},{4, 5}},
- {{4, 0},{0, 3},{3, 4}},
- {{2, 5},{1, 2},{5, 1}},
- {{3, 2},{1, 3},{2, 1}},
- {{5, 4},{1, 5},{4, 1}},
- {{4, 3},{1, 4},{3, 1}}
-};
-
-const int IndexTranslator::faceIndexFromElements_[subVolumeTotalNum][subVolumeTotalNum] =
-{
- {-1, 0, 3, -1, 8, -1, -1, -1},
- {0, -1, -1, 1, -1, 9, -1, -1},
- {3, -1, -1, 2, -1, -1, 11, -1},
- {-1, 1, 2, -1, -1, -1, -1, 10},
- {8, -1, -1, -1, -1, 4, 7, -1},
- {-1, 9, -1, -1, 4, -1, -1, 5},
- {-1, -1, 11, -1, 7, -1, -1, 6},
- {-1, -1, -1, 10, -1, 5, 6, -1}
-};
-//! \endcond
-
-//! \ingroup IMPET mpfa
-/*! \brief Class including the information of a 3d interaction volume of a MPFA L-method that does not change with time.
- *
- * Includes information needed to calculate the transmissibility matrices of an L-interaction-volume.
- *
- */
-template<class TypeTag>
-class FvMpfaL3dInteractionVolume
-{
-private:
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
-
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementSeed = typename Grid::template Codim<0>::EntitySeed;
-
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
-
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using FieldVectorVector = Dune::FieldVector<DimVector, dim>;
- using FieldVectorVector2 = Dune::FieldVector<DimVector, 2>;
- using FieldVectorVectorVector = Dune::FieldVector<FieldVectorVector2, dim>;
- using IndexVector = Dune::FieldVector<int, dim>;
- using BCTypeVector = std::vector<BoundaryTypes>;
- using BCVector = std::vector<PrimaryVariables>;
-
-public:
- //An interaction volume around a vertex includes in general 12 flux faces (see doc/docextra/3dmpfa).
- //If the vertex is on the boundary these flux faces can be inside the domain on the boundary or outside the domain.
- enum FaceTypes
- {
- inside = 1,//!< Flux face is inside the model domain
- boundary = 0,//!< Flux face is a boundary face
- outside = -1,//!< Flux face is outside the model domain
- };
-
- enum
- {
- subVolumeTotalNum = IndexTranslator::subVolumeTotalNum,//!< Number of sub-volumes in the interaction volume
- fluxFacesTotalNum = IndexTranslator::fluxFacesTotalNum,//!< Number of flux faces in the interaction volume
- fluxEdgesTotalNum = IndexTranslator::fluxEdgesTotalNum//!< Number of edges in the interaction volume
- };
-
- //! Constructs a FvMpfaL3dInteractionVolume object
- /**
- */
- FvMpfaL3dInteractionVolume(const Grid& grid)
- : grid_(&grid)
- , normal_(FieldVectorVector(DimVector(0.0)))
- , facePos_(DimVector(0.0))
- , edgePos_((DimVector(0.0)))
- , faceArea_(0.0)
- , faceType_(fluxFacesTotalNum, inside)
- , indexOnElement_(IndexVector(0.0))
- , elements_(subVolumeTotalNum)
- , centerVertexPos_(0)
- , elementNum_(0)
- {}
-
- //! Reset the interaction volume (deletes stored data)
- void reset()
- {
- elements_.clear();
- elements_.resize(subVolumeTotalNum);
- centerVertexPos_ = 0;
- indexOnElement_ = IndexVector(0.0);
- faceType_.clear();
- faceType_.resize(fluxFacesTotalNum, inside);
- faceArea_ = 0.0;
- facePos_ = DimVector(0.0);
- edgePos_ = DimVector(0.0);
- normal_ = FieldVectorVector(DimVector(0.0));
- elementNum_ = 0;
- boundaryTypes_.clear();
- neumannValues_.clear();
- dirichletValues_.clear();
- }
-
- //! Store position of the central vertex
- /*
- * \param centerVertexPos Position of the central vertex
- */
- void setCenterPosition(const DimVector ¢erVertexPos)
- {
- centerVertexPos_ = centerVertexPos;
- }
- //! Store a dune element as a sub volume element
- /*!
- * \param element The element
- * \param subVolumeIdx The local element index in the interaction volume
- */
- void setSubVolumeElement(const Element& element, int subVolumeIdx)
- {
- if (!hasSubVolumeElement(subVolumeIdx))
- {
- elements_[subVolumeIdx].push_back(element.seed());
- elementNum_++;
- }
- else
- {
- elements_[subVolumeIdx].insert(elements_[subVolumeIdx].begin(), element.seed());
- }
- }
-
- //! Store the position of a flux face
- /*!
- * \param pos Position of face center
- * \param fIdx The interaction volume face index
- */
- void setFacePosition(const DimVector& pos, int fIdx)
- {
- facePos_[fIdx] = pos;
- }
-
- //! Store the center of the edges
- /*!
- * \param pos Position of face center
- * \param edgeIdx The interaction volume edge index
- */
- void setEdgePosition(const DimVector& pos, int edgeIdx)
- {
- edgePos_[edgeIdx] = pos;
- }
-
- //! Store the flux face areas
- /*!
- * \param faceArea The flux face area
- * \param fIdx The interaction volume face index
- */
- void setFaceArea(Scalar faceArea, int fIdx)
- {
- faceArea_[fIdx] = faceArea;
- }
-
- //! Store the normals
- /*!
- * \param normal The normal vector
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- */
- void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdx)
- {
- normal_[subVolumeIdx][subVolumeFaceIdx] = normal;
- }
-
- //! Store the types of boundary conditions
- /*!
- * \param boundaryTypes BoundaryTypes object
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
- {
- if (boundaryTypes_.size() == 0)
- {
- boundaryTypes_.resize(fluxFacesTotalNum);
- }
- boundaryTypes_[subVolumeFaceIdx] = boundaryTypes;
- faceType_[subVolumeFaceIdx] = boundary;
- }
-
- //! Define a flux face to be outside the model domain (for boundary vertices)
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setOutsideFace(int subVolumeFaceIdx)
- {
- faceType_[subVolumeFaceIdx] = outside;
- }
-
- //! Store Dirichlet boundary conditions
- /*!
- * \param condition Vector of primary variables
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (dirichletValues_.size() == 0)
- {
- dirichletValues_.resize(fluxFacesTotalNum, PrimaryVariables(0.0));
- }
- dirichletValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store Neumann boundary conditions
- /*!
- * \param condition Vector phase fluxes
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (neumannValues_.size() == 0)
- {
- neumannValues_.resize(fluxFacesTotalNum, PrimaryVariables(0.0));
- }
- neumannValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store the local dune face index dependent on the local interaction volume indices
- /*!
- * \param indexInInside Face index of the Dune reference element
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- */
- void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdx)
- {
- indexOnElement_[subVolumeIdx][subVolumeFaceIdx] = indexInInside;
- }
-
- //! The position of the central vertex
- DimVector& getCenterPosition()
- {
- return centerVertexPos_;
- }
-
- //! The local interaction volume face index dependent on the local sub volume indices
- /*
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return The interaction volume face index
- */
- int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx);
- }
-
- //! The local interaction volume edge index dependent on the local sub volume indices
- /*
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- * \param subVolumeEdgeIdx The local edge index in the interaction volume element
- *
- * \return The interaction volume edge index
- */
- int getEdgeIndexFromSubVolumeFace(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
- {
- return IndexTranslator::getEdgeIndexFromSubVolumeFace(subVolumeIdx, subVolumeFaceIdx, subVolumeEdgeIdx);
- }
-
- //! Number of stored dune elements
- int getElementNumber()
- {
- return elementNum_;
- }
-
- //! The local dune face index dependent on the local interaction volume indices
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Face index of the Dune reference element
- */
- int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! The dune element of the interaction volume sub-volume
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- *
- * \return The interaction volume sub-element.
- */
- Element getSubVolumeElement(int subVolumeIdx)
- {
- if (hasSubVolumeElement(subVolumeIdx))
- return grid_->entity(elements_[subVolumeIdx][0]);
- else
- {
- std::cout<<"Problems when calling getSubVolumeElement("<<subVolumeIdx<<")\n";
- printInteractionVolumeInfo();
- DUNE_THROW(Dune::RangeError, "element not in interaction volume!");
- }
- }
-
- //! Check if a dune element is stored for an interaction volume sub-volume
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- *
- * \return returns <tt>true</tt> if an element pointer is stored at position <tt>subVolumeIdx</tt>.
- */
- bool hasSubVolumeElement(int subVolumeIdx)
- {
- return (elements_[subVolumeIdx].size() > 0);
- }
-
- //! The boundary types
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Object containing information about the boundary types.
- */
- BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
- {
- return boundaryTypes_[subVolumeFaceIdx];
- }
-
- //! Check if an interaction volume face is outside the model domain (for boundary vertices)
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return returns <tt>true</tt> if the flux face is outside the model domain.
- */
- bool isOutsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == outside);
- }
-
- //! Check if an interaction volume face is inside the model domain (for boundary vertices)
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return returns <tt>true</tt> if the flux face is inside the model domain.
- */
- bool isInsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == inside);
- }
-
- //! Check if an interaction volume face is a model domain boundary (for boundary vertices)
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return returns <tt>true</tt> if the flux face is a boundary face.
- */
- bool isBoundaryFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == boundary);
- }
-
- //! Check if an interaction volume is not located around a boundary vertex
- /*!
- * \return returns <tt>true</tt> if an interaction volume is not located around a boundary vertex
- */
- bool isInnerVolume()
- {
- for (unsigned int i = 0; i < faceType_.size(); i++)
- {
- if (isOutsideFace(i) || isBoundaryFace(i))
- return false;
- }
- return true;
- }
-
- //! Check if an interaction volume is located around a boundary vertex
- /*!
- * \return returns <tt>true</tt> if an interaction volume located around a boundary vertex
- */
- bool isBoundaryInteractionVolume()
- {
- return (boundaryTypes_.size() > 0);
- }
-
- //! The Dirichlet boundary values
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of primary variables.
- */
- PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
- {
- return dirichletValues_[subVolumeFaceIdx];
- }
-
- //! The Neumann boundary values
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of phase fluxes.
- */
- PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
- {
- return neumannValues_[subVolumeFaceIdx];
- }
-
- //! The face normal
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Normal vector
- */
- DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return normal_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! The position of the face center
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Position of face center
- */
- DimVector& getFacePosition(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return facePos_[IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx)];
- }
-
- //! The position of the edge center
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- * \param subVolumeEdgeIdx The local edge index in the interaction volume element
- *
- * \return Position of the edge center
- */
- DimVector& getEdgePosition(int subVolumeIdx, int subVolumeFaceIdx, int subVolumeEdgeIdx)
- {
- return edgePos_[IndexTranslator::getEdgeIndexFromSubVolumeFace(subVolumeIdx, subVolumeFaceIdx, subVolumeEdgeIdx)];
- }
-
- //! The interaction volume flux face area
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Area of the flux face
- */
- Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return faceArea_[IndexTranslator::getFaceIndexFromSubVolume(subVolumeIdx, subVolumeFaceIdx)];
- }
-
- //! The position of the face center
- /*!
- * \param fIdx The local face index in the interaction volume
- *
- * \return Position of face center
- */
- DimVector& getFacePosition(int fIdx)
- {
- return facePos_[fIdx];
- }
-
- //! The position of the edge center
- /*!
- * \param edgeIdx The local edge index in the interaction volume
- *
- * \return Position of the edge center
- */
- DimVector& getEdgePosition(int edgeIdx)
- {
- return edgePos_[edgeIdx];
- }
-
- //! The interaction volume flux face area
- /*!
- * \param fIdx The local face index in the interaction volume
- *
- * \return Area of the flux face
- */
- Scalar& getFaceArea(int fIdx)
- {
- return faceArea_[fIdx];
- }
-
- //! Print the stored interaction volume data
- void printInteractionVolumeInfo()
- {
- std::cout<<"\nNumber of stored elements: "<<elementNum_<<"\n";
- std::cout<<" center position: "<<centerVertexPos_<<"\n";
- for (int i = 0; i < subVolumeTotalNum; i++)
- {
- if (elements_[i].size() > 0)
- {
- std::cout<<"element "<<i<<":\n";
- std::cout<<"element level: "<<grid_->entity(elements_[i][0]).level()<<"\n";
- std::cout<<"element position: "<<grid_->entity(elements_[i][0]).geometry().center()<<"\n";
- std::cout<<"element volume: "<<grid_->entity(elements_[i][0]).geometry().volume()<<"\n";
- std::cout<<"face indices on element: "<<indexOnElement_[i]<<"\n";
- std::cout<<"face normals on element: "<<normal_[i]<<"\n";
- std::cout<<"face areas on element: ";
- for (int j = 0; j < 3; j++)
- {
- std::cout<<getFaceArea(i, j)<<" ";
- }
- std::cout<<"\n";
- std::cout<<"face position on element: ";
- for (int j = 0; j < 3; j++)
- {
- std::cout<<getFacePosition(i, j)<<" ";
- }
- std::cout<<"\n";
- std::cout<<"edgePos:";
- for (int j = 0; j < 3; j++)
- {
- std::cout<<" ";
- for (int k = 0; k < 2; k++)
- {
- std::cout<<getEdgePosition(i, j, k) <<" ";
- }
- }
- std::cout<<"\n";
- }
- }
- }
-
-// void printInteractionVolumeInfoToFile(std::ofstream& dataFile)
-// {
-// dataFile<<"\nNumber of stored elements: "<<elementNum_<<"\n";
-// dataFile<<" center position: "<<centerVertexPos_<<"\n";
-// for (int i = 0; i < subVolumeTotalNum; i++)
-// {
-// if (elements_[i].size() > 0)
-// {
-// dataFile<<"element "<<i<<":\n";
-// dataFile<<"element level: "<<elements_[i][0]->level()<<"\n";
-// dataFile<<"element position: "<<elements_[i][0]->geometry().center()<<"\n";
-// dataFile<<"element volume: "<<elements_[i][0]->geometry().volume()<<"\n";
-// dataFile<<"face indices on element: "<<indexOnElement_[i]<<"\n";
-// dataFile<<"face normals on element: "<<normal_[i]<<"\n";
-// dataFile<<"face areas on element: ";
-// for (int j = 0; j < 3; j++)
-// {
-// dataFile<<getFaceArea(i, j)<<" ";
-// }
-// dataFile<<"\n";
-// dataFile<<"face position on element: ";
-// for (int j = 0; j < 3; j++)
-// {
-// dataFile<<getFacePosition(i, j)<<" ";
-// }
-// dataFile<<"\n";
-// dataFile<<"edgePos:";
-// for (int j = 0; j < 3; j++)
-// {
-// dataFile<<" ";
-// for (int k = 0; k < 2; k++)
-// {
-// dataFile<<getEdgePosition(i, j, k) <<" ";
-// }
-// }
-// dataFile<<"\n";
-//
-// }
-// dataFile<<"face types: ";
-// for (int i = 0; i < fluxFacesTotalNum; i++)
-// {
-// dataFile<<faceType_[i]<<" ";
-// }
-// dataFile<<"\n";
-// if (isBoundaryInteractionVolume())
-// {
-// dataFile<<"Boundaries:\n";
-// dataFile<<"dirichlet: ";
-// for (int i = 0; i < fluxFacesTotalNum; i++)
-// {
-// dataFile<<boundaryTypes_[i].isDirichlet(0)<<" ";
-// dataFile<<boundaryTypes_[i].isDirichlet(1)<<" ";
-// }
-// dataFile<<"\n";
-// dataFile<<"neumann: ";
-// for (int i = 0; i < fluxFacesTotalNum; i++)
-// {
-// dataFile<<boundaryTypes_[i].isNeumann(0)<<" ";
-// dataFile<<boundaryTypes_[i].isNeumann(1)<<" ";
-// }
-// dataFile<<"\n";
-// dataFile<<"values: ";
-// for (int i = 0; i < fluxFacesTotalNum; i++)
-// {
-// if (dirichletValues_.size() > 0)
-// dataFile<<dirichletValues_[i]<<" ";
-// if (neumannValues_.size() > 0)
-// dataFile<<neumannValues_[i]<<" ";
-// }
-// dataFile<<"\n";
-// }
-// }
-
-private:
- const Grid* grid_;
- Dune::FieldVector<FieldVectorVector, subVolumeTotalNum> normal_;
- Dune::FieldVector<DimVector, fluxFacesTotalNum> facePos_;
- Dune::FieldVector<DimVector, fluxEdgesTotalNum> edgePos_;
- Dune::FieldVector<Scalar, fluxFacesTotalNum> faceArea_;
- BCTypeVector boundaryTypes_;
- std::vector<int> faceType_;
- Dune::FieldVector<IndexVector, subVolumeTotalNum> indexOnElement_;
- std::vector<std::vector<ElementSeed> > elements_;
- BCVector neumannValues_;
- BCVector dirichletValues_;
- DimVector centerVertexPos_;
- int elementNum_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3dadaptive.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3dadaptive.hh
deleted file mode 100644
index ee320f9dcb7a48b0d95424791e88c3fc01028462..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3dadaptive.hh
+++ /dev/null
@@ -1,346 +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 3 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/>. *
- *****************************************************************************/
-
-#ifndef DUMUX_FVMPFAL3DINTERACTIONVOLUME_ADAPTIVE_HH
-#define DUMUX_FVMPFAL3DINTERACTIONVOLUME_ADAPTIVE_HH
-
-/**
- * @file
- * @brief Class including the information of an interaction volume of a MPFA 3D method that does not change with time
- * @author Markus Wolff
- */
-
-#include <set>
-
-#include "linteractionvolume3d.hh"
-
-namespace Dumux
-{
-//! \cond \private
-// Mapper for local interaction volume indices (see doc/docextra/3dmpfa).
-class IndexTranslatorAdaptive: public IndexTranslator
-{
-public:
- enum
- {
- subVolumeTotalNum = 8,
- fluxFacesTotalNum = 12,
- fluxFacesNumOnSubVolume = 3,
- fluxEdgesTotalNum = 6,
- edgesNumOnFluxFace = 2
- };
-
- static int getOldElemIdxFromNewFaceIdxto0(int zeroFaceIdx, int elementIdx)
- {
- return oldElemIdxFromNewFaceIdxto0_[zeroFaceIdx][elementIdx];
- }
-
- static int getNewElemIdxFromOldFaceIdxto0(int zeroFaceIdx, int elementIdx)
- {
- return newElemIdxFromOldFaceIdxto0_[zeroFaceIdx][elementIdx];
- }
-
- static int getOldFaceIdxFromNewIdxto0(int zeroFaceIdx, int fIdx)
- {
- return oldFaceIdxFromNewIdxto0_[zeroFaceIdx][fIdx];
- }
-
- static int getNewFaceIdxFromOldIdxto0(int zeroFaceIdx, int fIdx)
- {
- return newFaceIdxFromOldIdxto0_[zeroFaceIdx][fIdx];
- }
-
- static int getOldEdgeIdxFromNewFaceIdxto0(int zeroFaceIdx, int edgeIdx)
- {
- return oldEdgeIdxFromNewFaceIdxto0_[zeroFaceIdx][edgeIdx];
- }
-
- static int getNewEdgeIdxFromOldFaceIdxto0(int zeroFaceIdx, int edgeIdx)
- {
- return newEdgeIdxFromOldFaceIdxto0_[zeroFaceIdx][edgeIdx];
- }
-
-
-private:
- static const int oldElemIdxFromNewFaceIdxto0_[fluxFacesTotalNum][subVolumeTotalNum];
- static const int newElemIdxFromOldFaceIdxto0_[fluxFacesTotalNum][subVolumeTotalNum];
- static const int oldFaceIdxFromNewIdxto0_[fluxFacesTotalNum][fluxFacesTotalNum];
- static const int newFaceIdxFromOldIdxto0_[fluxFacesTotalNum][fluxFacesTotalNum];
- static const int oldEdgeIdxFromNewFaceIdxto0_[fluxFacesTotalNum][fluxEdgesTotalNum];
- static const int newEdgeIdxFromOldFaceIdxto0_[fluxFacesTotalNum][fluxEdgesTotalNum];
-};
-
-const int IndexTranslatorAdaptive::oldElemIdxFromNewFaceIdxto0_[fluxFacesTotalNum][subVolumeTotalNum] =
-{
- {0, 1, 2, 3, 4, 5, 6, 7},
- {1, 3, 0, 2, 5, 7, 4, 6},
- {3, 2, 1, 0, 7, 6, 5, 4},
- {2, 0, 3, 1, 6, 4, 7, 5},
- {5, 4, 7, 6, 1, 0, 3, 2},
- {7, 5, 6, 4, 3, 1, 2, 0},
- {6, 7, 4, 5, 2, 3, 0, 1},
- {4, 6, 5, 7, 0, 2, 1, 3},
- {0, 4, 1, 5, 2, 6, 3, 7},
- {1, 5, 3, 7, 0, 4, 2, 6},
- {3, 7, 2, 6, 1, 5, 0, 4},
- {2, 6, 0, 4, 3, 7, 1, 5}
-};
-
-const int IndexTranslatorAdaptive::newElemIdxFromOldFaceIdxto0_[fluxFacesTotalNum][subVolumeTotalNum] =
-{
- {0, 1, 2, 3, 4, 5, 6, 7},
- {2, 0, 3, 1, 6, 4, 7, 5},
- {3, 2, 1, 0, 7, 6, 5, 4},
- {1, 3, 0, 2, 5, 7, 4, 6},
- {5, 4, 7, 6, 1, 0, 3, 2},
- {7, 5, 6, 4, 3, 1, 2, 0},
- {6, 7, 4, 5, 2, 3, 0, 1},
- {4, 6, 5, 7, 0, 2, 1, 3},
- {0, 2, 4, 6, 1, 3, 5, 7},
- {4, 0, 6, 2, 5, 1, 7, 3},
- {6, 4, 2, 0, 7, 5, 3, 1},
- {2, 6, 0, 4, 3, 7, 1, 5}
-};
-
-const int IndexTranslatorAdaptive::oldFaceIdxFromNewIdxto0_[fluxFacesTotalNum][fluxFacesTotalNum] =
-{
- {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
- {1, 2, 3, 0, 5, 6, 7, 4, 9, 10, 11, 8},
- {2, 3, 0, 1, 6, 7, 4, 5, 10, 11, 8, 9},
- {3, 0, 1, 2, 7, 4, 5, 6, 11, 8, 9, 10},
- {4, 7, 6, 5, 0, 3, 2, 1, 9, 8, 11, 10},
- {5, 4, 7, 6, 1, 0, 3, 2, 10, 9, 8, 11},
- {6, 5, 4, 7, 2, 1, 0, 3, 11, 10, 9, 8},
- {7, 6, 5, 4, 3, 2, 1, 0, 8, 11, 10, 9},
- {8, 4, 9, 0, 11, 6, 10, 2, 3, 7, 5, 1},
- {9, 5, 10, 1, 8, 7, 11, 3, 0, 4, 6, 2},
- {10, 6, 11, 2, 9, 4, 8, 0, 1, 5, 7, 3},
- {11, 7, 8, 3, 10, 5, 9, 1, 2, 6, 4, 0}
-};
-
-const int IndexTranslatorAdaptive::newFaceIdxFromOldIdxto0_[fluxFacesTotalNum][fluxFacesTotalNum] =
-{
- {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
- {3, 0, 1, 2, 7, 4, 5, 6, 11, 8, 9, 10},
- {2, 3, 0, 1, 6, 7, 4, 5, 10, 11, 8, 9},
- {1, 2, 3, 0, 5, 6, 7, 4, 9, 10, 11, 8},
- {4, 7, 6, 5, 0, 3, 2, 1, 9, 8, 11, 10},
- {5, 4, 7, 6, 1, 0, 3, 2, 10, 9, 8, 11},
- {6, 5, 4, 7, 2, 1, 0, 3, 11, 10, 9, 8},
- {7, 6, 5, 4, 3, 2, 1, 0, 8, 11, 10, 9},
- {3, 11, 7, 8, 1, 10, 5, 9, 0, 2, 6, 4},
- {8, 3, 11, 7, 9, 1, 10, 5, 4, 0, 2, 6},
- {7, 8, 3, 11, 5, 9, 1, 10, 6, 4, 0, 2},
- {11, 7, 8, 3, 10, 5, 9, 1, 2, 6, 4, 0}
-};
-
-const int IndexTranslatorAdaptive::oldEdgeIdxFromNewFaceIdxto0_[fluxFacesTotalNum][fluxEdgesTotalNum] =
-{
- {0, 1, 2, 3, 4, 5},
- {0, 1, 3, 4, 5, 2},
- {0, 1, 4, 5, 2, 3},
- {0, 1, 5, 2, 3, 4},
- {1, 0, 2, 5, 4, 3},
- {1, 0, 3, 2, 5, 4},
- {1, 0, 4, 3, 2, 5},
- {1, 0, 5, 4, 3, 2},
- {2, 4, 5, 1, 3, 0},
- {3, 5, 2, 1, 4, 0},
- {4, 2, 3, 1, 5, 0},
- {5, 3, 4, 1, 2, 0}
-};
-
-const int IndexTranslatorAdaptive::newEdgeIdxFromOldFaceIdxto0_[fluxFacesTotalNum][fluxEdgesTotalNum] =
-{
- {0, 1, 2, 3, 4, 5},
- {0, 1, 5, 2, 3, 4},
- {0, 1, 4, 5, 2, 3},
- {0, 1, 3, 4, 5, 2},
- {1, 0, 2, 5, 4, 3},
- {1, 0, 3, 2, 5, 4},
- {1, 0, 4, 3, 2, 5},
- {1, 0, 5, 4, 3, 2},
- {5, 3, 0, 4, 1, 2},
- {5, 3, 2, 0, 4, 1},
- {5, 3, 1, 2, 0, 4},
- {5, 3, 4, 1, 2, 0}
-};
-//! \endcond
-
-//! \ingroup IMPET mpfa
-/*! \brief Class including the information of a 3d interaction volume of an adaptive MPFA L-method that does not change with time.
- *
- * Includes information needed to calculate the transmissibility matrices of an L-interaction-volume.
- *
- */
-template<class TypeTag>
-class FvMpfaL3dInteractionVolumeAdaptive:public FvMpfaL3dInteractionVolume<TypeTag>
-{
-private:
- using ParentType = FvMpfaL3dInteractionVolume<TypeTag>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
-
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementSeed = typename Grid::template Codim<0>::EntitySeed;
-
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
-
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using FieldVectorVector = Dune::FieldVector<DimVector, dim>;
- using FieldVectorVector2 = Dune::FieldVector<DimVector, 2>;
- using FieldVectorVectorVector = Dune::FieldVector<FieldVectorVector2, dim>;
- using IndexVector = Dune::FieldVector<int, dim>;
- using BCTypeVector = std::vector<BoundaryTypes>;
- using BCVector = std::vector<PrimaryVariables>;
-
-public:
-
- //! \cond \private
- enum FaceTypes
- {
- inside = 1,
- boundary = 0,
- outside = -1,
- };
-
- //!\copydoc
- enum
- {
- subVolumeTotalNum = IndexTranslatorAdaptive::subVolumeTotalNum,
- fluxFacesTotalNum = IndexTranslatorAdaptive::fluxFacesTotalNum,
- fluxEdgesTotalNum = IndexTranslatorAdaptive::fluxEdgesTotalNum
- };
- //! \endcond
-
- //! The different hanging node interaction volume types (see dissertation M. Wolff, http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/)
- enum HangingNodeTypes
- {
- noHangingNode = -1,//!< regular interaction volume
- twoSmallCells = 0,//!< hanging-node interaction volume of type 5 or 7
- fourSmallCellsFace = 1,//!< hanging-node interaction volume of type 1
- fourSmallCellsEdge = 2,//!< hanging-node interaction volume of type 3
- fourSmallCellsDiag = 3,//!< hanging-node interaction volume of type 4
- sixSmallCells = 4//!< hanging-node interaction volume of type 2 or 6
- };
-
- //! Constructs a FvMpfaL3dInteractionVolumeAdaptive object
- /**
- */
- FvMpfaL3dInteractionVolumeAdaptive(const Grid& grid)
- : ParentType(grid)
- , hangingNodeType_(noHangingNode)
- {}
-
- //!\copydoc FvMpfaL3dInteractionVolume::reset()
- void reset()
- {
- hangingNodeType_ = noHangingNode;
- existingLevel_.clear();
- }
-
- //!\copydoc FvMpfaL3dInteractionVolume::setSubVolumeElement()
- void setSubVolumeElement(const Element& element, int subVolumeIdx)
- {
- ParentType::setSubVolumeElement(element, subVolumeIdx);
- existingLevel_.insert(element.level());
- }
-
- //! Store the type of hanging-node-interaction volume
- /*!
- * \param hNType the type of hanging-node-interaction volume of type FvMpfaL3dInteractionVolumeAdaptive<TypeTag>::HangingNodeTypes
- */
- void setHangingNodeType(int hNType)
- {
- hangingNodeType_ = hNType;
- }
-
- //! Check if elements in the interaction volume are of the same grid level
- /*!
- * \return <tt>true</tt> if all elements are on the same grid level
- */
- bool sameLevel()
- {
- if (!isHangingNodeVolume())
- return existingLevel_.size() < 2;
- else
- {
- return existingLevel_.size() < 3;
- }
- }
-
- //! Check if an element of a certain grid level is stored
- /*!
- * \param level the dune grid level
- *
- * \return <tt>true</tt> if an element of a certain grid level is stored
- */
- bool hasLevel(int level)
- {
- return existingLevel_.find(level) != existingLevel_.end();
- }
-
- //! Check whether the interaction volume is a hanging-node volume
- /*!
- * \return <tt>true</tt> if the interaction volume is a hanging-node volume
- */
- bool isHangingNodeVolume()
- {
- return hangingNodeType_ != noHangingNode;
- }
-
- //! The type of the interaction volume as type of FvMpfaL3dInteractionVolumeAdaptive<TypeTag>::HangingNodeTypes
- int getHangingNodeType()
- {
- return hangingNodeType_;
- }
-
- //!\copydoc FvMpfaL3dInteractionVolume::printInteractionVolumeInfo()
- void printInteractionVolumeInfo()
- {
- ParentType::printInteractionVolumeInfo();
-
- if (isHangingNodeVolume())
- std::cout<<"hanging node type: "<<hangingNodeType_<<"\n";
- }
-
-// void printInteractionVolumeInfoToFile(std::ofstream& dataFile)
-// {
-// ParentType::printInteractionVolumeInfoToFile(dataFile);
-//
-// if (isHangingNodeVolume())
-// dataFile<<"hanging node type: "<<hangingNodeType_<<"\n";
-// }
-
-private:
- int hangingNodeType_;
- std::set<int> existingLevel_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh
deleted file mode 100644
index 5e7d0128075abab3eaadfe9fb6ec35b6939ec2c8..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh
+++ /dev/null
@@ -1,458 +0,0 @@
-/*****************************************************************************
- * 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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_FVMPFAOINTERACTIONVOLUME_HH
-#define DUMUX_FVMPFAOINTERACTIONVOLUME_HH
-
-/**
- * @file
- * @brief Class including the information of an interaction volume of a MPFA O-method that does not change with time.
- */
-
-#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
-
-namespace Dumux
-{
-
-//! \ingroup IMPET
-/*! \brief Class including the information of an interaction volume of a MPFA O-method that does not change with time.
- *
- * Includes information needed to calculate the transmissibility matrix of an O-interaction-volume.
- *
- */
-template<class TypeTag>
-class FVMPFAOInteractionVolume
-{
-private:
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- enum
- {
- dim = GridView::dimension
- };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementSeed = typename Grid::template Codim<0>::EntitySeed;
-
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
-
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
- using FieldVectorVector = Dune::FieldVector<DimVector, dim>;
- using IndexVector = Dune::FieldVector<int, dim>;
- using BCTypeVector = std::vector<BoundaryTypes>;
- using BCVector = std::vector<PrimaryVariables>;
-
-public:
- enum FaceTypes
- {
- inside = 1,
- boundary = 0,
- outside = -1
- };
-
- //! Constructs a FVMPFAOInteractionVolume object
- FVMPFAOInteractionVolume(const Grid& grid)
- : grid_(&grid),
- stored_(false), permTimesNu_(FieldVectorVector(DimVector(0.0))),
- nu_(FieldVectorVector(DimVector(0.0))), normal_(FieldVectorVector(DimVector(0.0))),
- faceArea_(DimVector(0.0)), dF_(0.0),
- faceType_(2*dim, inside),
- indexOnElement_(IndexVector(0.0)),
- elements_(2*dim)
- {
- faceIndexOnSubVolume_[0][0] = 0;
- faceIndexOnSubVolume_[0][1] = 3;
- faceIndexOnSubVolume_[1][0] = 1;
- faceIndexOnSubVolume_[1][1] = 0;
- faceIndexOnSubVolume_[2][0] = 2;
- faceIndexOnSubVolume_[2][1] = 1;
- faceIndexOnSubVolume_[3][0] = 3;
- faceIndexOnSubVolume_[3][1] = 2;
- }
-
- //! Mark storage as completed
- void setStored()
- {
- stored_ = true;
- }
-
- //! Returns true if information has already been stored
- bool isStored() const
- {
- return stored_;
- }
-
- //! Store an element of the interaction volume
- /*!
- * \param element The element
- * \param subVolumeIdx The local element index in the interaction volume
- */
- void setSubVolumeElement(const Element& element, int subVolumeIdx)
- {
- elements_[subVolumeIdx].push_back(element.seed());
- }
-
- //! Store the \f$ dF \f$ for the transmissiblity calculation
- /*!
- * \param dF Value of \f$ dF \f$
- * \param subVolumeIdx The local element index in the interaction volume
- */
- void setDF(Scalar dF, int subVolumeIdx)
- {
- dF_[subVolumeIdx] = dF;
- }
-
- //! Store a flux face area
- /*!
- * \param faceArea The flux face area
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setFaceArea(Scalar& faceArea, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- faceArea_[subVolumeIdx][subVolumeFaceIdxInInside] = faceArea;
- }
-
- //! Store \f$ \boldsymbol K \boldsymbol \nu\f$ for the transmissiblity calculation
- /*!
- * \param nu \f$ \boldsymbol \nu \f$ vector
- * \param perm Permeability matrix
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setPermTimesNu(DimVector& nu, DimMatrix& perm, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInside] = 0;
-
- perm.mv(nu, permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInside]);
-
- nu_[subVolumeIdx][subVolumeFaceIdxInInside] = nu;
- }
-
- //! Store a flux face normal
- /*!
- * \param normal The normal vector
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setNormal(DimVector& normal, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- normal_[subVolumeIdx][subVolumeFaceIdxInInside] = normal;
- }
-
- //! Store boundary condtion types for a flux face
- /*!
- * \param boundaryTypes BoundaryTypes object
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setBoundary(BoundaryTypes& boundaryTypes, int subVolumeFaceIdx)
- {
-// std::cout<<"old BCType = "<<boundaryType_[subVolumeFaceIdx]<<"\n";
-// std::cout<<"new BCType = "<<boundaryType<<"\n";
- if (boundaryTypes_.size() == 0)
- {
- boundaryTypes_.resize(2 * dim);
- }
- boundaryTypes_[subVolumeFaceIdx] = boundaryTypes;
- faceType_[subVolumeFaceIdx] = boundary;
- }
-
- //! Mark a flux face to be outside the model domain
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setOutsideFace(int subVolumeFaceIdx)
- {
- faceType_[subVolumeFaceIdx] = outside;
- }
-
- //! Store Dirichlet boundary condtions for a flux face
- /*!
- * \param condition Vector of primary variables
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setDirichletCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (dirichletValues_.size() == 0)
- {
- dirichletValues_.resize(2 * dim);
- }
- dirichletValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store Neumann boundary condtions for a flux face
- /*!
- * \param condition Vector phase fluxes
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- void setNeumannCondition(PrimaryVariables& condition, int subVolumeFaceIdx)
- {
- if (neumannValues_.size() == 0)
- {
- neumannValues_.resize(2 * dim);
- }
- neumannValues_[subVolumeFaceIdx] = condition;
- }
-
- //! Store map from local interaction volume numbering to numbering of the Dune reference element.
- /*!
- * \param indexInInside Face index of the Dune reference element
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- */
- void setIndexOnElement(int indexInInside, int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- indexOnElement_[subVolumeIdx][subVolumeFaceIdxInInside] = indexInInside;
- }
-
- //! Map from local interaction volume numbering to numbering of the Dune reference element.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Face index of the Dune reference element
- */
- int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return indexOnElement_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! Map from local interaction volume numbering on element to numbering on interaction volume.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdx The local face index in the interaction volume element
- *
- * \return Local face index int the interaction volume
- */
- int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
- {
- return faceIndexOnSubVolume_[subVolumeIdx][subVolumeFaceIdx];
- }
-
- //! Get an element of the interaction volume.
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- *
- * \return The interaction volume sub-element.
- */
- Element getSubVolumeElement(int subVolumeIdx)
- {
- return grid_->entity(elements_[subVolumeIdx][0]);
- }
-
- //! Get boundary condtion types for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Object containing information about the boundary types.
- */
- BoundaryTypes& getBoundaryType(int subVolumeFaceIdx)
- {
- return boundaryTypes_[subVolumeFaceIdx];
- }
-
- //! Returns true if an interaction volume flux face is outside the model domain.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isOutsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == outside);
- }
-
- //! Returns true if an interaction volume flux face is inside the model domain and no boundary face.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isInsideFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == inside);
- }
-
- //! Returns true if the interaction volume is completely inside the model domain.
- bool isInnerVolume()
- {
- for (unsigned int i = 0; i < faceType_.size(); i++)
- {
- if (isOutsideFace(i) || isBoundaryFace(i))
- return false;
- }
- return true;
- }
-
- //! Returns true if an interaction volume flux face is a boundary face.
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- */
- bool isBoundaryFace(int subVolumeFaceIdx)
- {
- return (faceType_[subVolumeFaceIdx] == boundary);
- }
-
- //! Get the Dirichlet boundary condtions for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of primary variables.
- */
- PrimaryVariables& getDirichletValues(int subVolumeFaceIdx)
- {
- return dirichletValues_[subVolumeFaceIdx];
- }
-
- //! Get the Neumann boundary condtions for a flux face
- /*!
- * \param subVolumeFaceIdx The local face index in the interaction volume
- *
- * \return Vector of phase fluxes.
- */
- PrimaryVariables& getNeumannValues(int subVolumeFaceIdx)
- {
- return neumannValues_[subVolumeFaceIdx];
- }
-
- //! Get a flux face normal
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return Normal vector
- */
- DimVector& getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return normal_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Get a flux face area
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return Area of the flux face
- */
- Scalar& getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return faceArea_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Get \f$ \boldsymbol \nu \f$ vector used for the transmissiblity calculation
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInside The local face index in the interaction volume element
- *
- * \return \f$ \boldsymbol \nu \f$ vector
- */
- DimVector& getNu(int subVolumeIdx, int subVolumeFaceIdxInInside)
- {
- return nu_[subVolumeIdx][subVolumeFaceIdxInInside];
- }
-
- //! Get \f$ \boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \f$ for the transmissiblity calculation
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element
- * for the normal \f$ \boldsymbol n \f$
- * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element
- * for the vector \f$ \boldsymbol \nu \f$
- *
- * \return \f$ \boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \f$
- */
- Scalar getNtkNu(int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
- {
- return normal_[subVolumeIdx][subVolumeFaceIdxInInsideN] * permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInsideNu];
- }
-
- //! Get \f$ \boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu\f$ for the transmissiblity calculation
- /*!
- * \param relPerm relative permeability value (\f$ \boldsymbol n^\text{T} k_{r\alpha} \f$)
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element
- * for the normal \f$ \boldsymbol n \f$
- * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element
- * for the vector \f$ \boldsymbol \nu \f$
- *
- * \return \f$ \boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \f$
- */
- Scalar getNtkrkNu(Scalar& relPerm, int subVolumeIdx, int subVolumeFaceIdxInInsideN,
- int subVolumeFaceIdxInInsideNu) const
- {
- DimVector krKNu(permTimesNu_[subVolumeIdx][subVolumeFaceIdxInInsideNu]);
-
- krKNu *= relPerm;
-
- return normal_[subVolumeIdx][subVolumeFaceIdxInInsideN] * krKNu;
- }
-
- //! Get \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \right) \f$
- //! for the transmissiblity calculation
- /*!
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element
- * for the normal \f$ \boldsymbol n \f$
- * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element
- * for the vector \f$ \boldsymbol \nu \f$
- *
- * \return \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} \boldsymbol K \boldsymbol \nu \right) \f$
- */
- Scalar getNtkNu_df(int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
- {
- return faceArea_[subVolumeIdx][subVolumeFaceIdxInInsideN]
- * getNtkNu(subVolumeIdx, subVolumeFaceIdxInInsideN, subVolumeFaceIdxInInsideNu)
- / dF_[subVolumeIdx];
- }
-
- //! Get \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \right) \f$
- //! for the transmissiblity calculation
- /*!
- * \param relPerm relative permeability value (\f$ \boldsymbol n^\text{T} k_{r\alpha} \f$)
- * \param subVolumeIdx The local element index in the interaction volume
- * \param subVolumeFaceIdxInInsideN The local face index in the interaction volume element
- * for the normal \f$ \boldsymbol n \f$
- * \param subVolumeFaceIdxInInsideNu The local face index in the interaction volume element
- * for the vector \f$ \boldsymbol \nu \f$
- *
- * \return \f$ \frac{1}{dF} \left(\boldsymbol n^\text{T} k_{r\alpha} \boldsymbol K \boldsymbol \nu \right) \f$
- */
- Scalar getNtkrkNu_df(Scalar& relPerm, int subVolumeIdx, int subVolumeFaceIdxInInsideN, int subVolumeFaceIdxInInsideNu) const
- {
- return faceArea_[subVolumeIdx][subVolumeFaceIdxInInsideN]
- * getNtkrkNu(relPerm, subVolumeIdx, subVolumeFaceIdxInInsideN, subVolumeFaceIdxInInsideNu)
- / dF_[subVolumeIdx];
- }
-
-private:
- const Grid* grid_;
- bool stored_;
- Dune::FieldVector<FieldVectorVector, 2*dim> permTimesNu_;
- Dune::FieldVector<FieldVectorVector, 2*dim> nu_;
- Dune::FieldVector<FieldVectorVector, 2*dim> normal_;
- Dune::FieldVector<DimVector, 2*dim> faceArea_;
- Dune::FieldVector<Scalar, 2*dim> dF_;
- BCTypeVector boundaryTypes_;
- std::vector<int> faceType_;
- Dune::FieldVector<IndexVector, 2*dim> indexOnElement_;
- Dune::FieldVector<IndexVector, 2*dim> faceIndexOnSubVolume_;
- std::vector<std::vector<ElementSeed> > elements_;
- BCVector neumannValues_;
- BCVector dirichletValues_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh
deleted file mode 100644
index db50d17519fca541996dd80c8226735864ce532d..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh
+++ /dev/null
@@ -1,131 +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 3 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/>. *
- *****************************************************************************/
-
-/*!
- * \ingroup IMPET
- * \ingroup IMPETProperties mpfa
- * \file
- *
- * \brief Properties for a MPFA method.
- */
-#ifndef DUMUX_FVMPFAPROPERTIES_HH
-#define DUMUX_FVMPFAPROPERTIES_HH
-
-// dumux environment
-#include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-#include <dune/grid/yaspgrid.hh>
-#if HAVE_DUNE_ALUGRID
-#include <dune/alugrid/grid.hh>
-#endif
-#if HAVE_DUNE_UGGRID
-#include <dune/grid/uggrid.hh>
-#endif
-
-namespace Dumux
-{
-/*!
- *
- *
- * \brief Indices denoting the different grid types.
- */
-struct GridTypes
-{
-public:
- //
- static const int any = 0;
- //YaspGrid
- static const int yaspGrid = 2;
- //UGGrid
- static const int ugGrid = 3;
- //ALUGrid
- static const int aluGrid = 4;
-};
-//! \cond \private
-template<class Grid, int dim>
-struct GridImp
-{
- static const int imp = GridTypes::any;
-};
-
-template<int dim>
-struct GridImp<Dune::YaspGrid<dim>, dim>
-{
- static const int imp = GridTypes::yaspGrid;
-};
-
-#if HAVE_DUNE_ALUGRID
-template<int dim>
-struct GridImp<Dune::ALUGrid<dim, dim, Dune::cube, Dune::nonconforming>, dim>
-{
- static const int imp = GridTypes::aluGrid;
-};
-#endif
-
-#if HAVE_DUNE_UGGRID
-template<int dim>
-struct GridImp<Dune::UGGrid<dim>, dim>
-{
- static const int imp = GridTypes::ugGrid;
-};
-#endif
-//! \endcond
-
-namespace Properties
-{
-//! Basic Type tag for MFPA models
-namespace TTag {
-struct MPFAProperties {};
-}
-
-template<class TypeTag, class MyTypeTag>
-struct GridTypeIndices { using type = UndefinedProperty; };//!< The grid type indices to decide which grid is used
-template<class TypeTag, class MyTypeTag>
-struct GridImplementation { using type = UndefinedProperty; }; //!< Gives kind of grid implementation in form of a GridType
-template<class TypeTag, class MyTypeTag>
-struct MPFAInteractionVolume { using type = UndefinedProperty; }; //!< Type of the data container for one interaction volume
-template<class TypeTag, class MyTypeTag>
-struct MPFAInteractionVolumeContainer { using type = UndefinedProperty; }; //!< Type of the data container for one interaction volume
-template<class TypeTag, class MyTypeTag>
-struct MPFATransmissibilityCalculator { using type = UndefinedProperty; }; //!< Type defining the transmissibility calculation
-}
-}
-
-namespace Dumux
-{
-namespace Properties
-{
-
-//! \cond \private
-template<class TypeTag>
-struct GridImplementation<TypeTag, TTag::MPFAProperties>
-{
-private:
- using Grid = GetPropType<TypeTag, Properties::Grid>;
-public:
- static const int value = GridImp<Grid, Grid::dimension>::imp;
-};
-//! \endcond
-
-//! Set grid type indices
-template<class TypeTag>
-struct GridTypeIndices<TypeTag, TTag::MPFAProperties> { using type = GridTypes; };
-}
-}// end of Dune namespace
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/mpfa/velocityintransport.hh b/dumux/porousmediumflow/sequential/cellcentered/mpfa/velocityintransport.hh
deleted file mode 100644
index e0929fdc52915d10b199d3649e9531de67aa9d10..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/mpfa/velocityintransport.hh
+++ /dev/null
@@ -1,98 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_FVMPFAVELOCITYINTRANSPORT_HH
-#define DUMUX_FVMPFAVELOCITYINTRANSPORT_HH
-
-/**
- * @file
- * @brief Implementation of a interface velocity class for MPFA models.
- */
-
-#include <dune/grid/common/gridenums.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-#include "properties.hh"
-
-namespace Dumux
-{
-//! \ingroup IMPET
-/*! \brief Implementation of a interface velocity class for MPFA models.
- *
- * Allows to calculate the MPFA-velocity in the transport model. (Less efficient!)
- *
- * \tparam TypeTag The Type Tag
- */
-template<class TypeTag>
-class FvMpfaVelocityInTransport
-{
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- using CellData = GetPropType<TypeTag, Properties::CellData>;
- using Intersection = typename GridView::Intersection;
-
-
-public:
- //! Constructs a FvMpfaVelocityInTransport object
- /*!
- * \param problem A problem class object
- */
- FvMpfaVelocityInTransport(Problem& problem):
- problem_(problem)
- {
- calcVelocityInTransport_ = getParam<bool>("MPFA.CalcVelocityInTransport");
- }
-
- //! For initialization
- void initialize()
- {}
-
- //! Local velocity calculation
- void calculateVelocity(const Intersection& intersection, CellData& cellData)
- {
- problem_.pressureModel().calculateVelocity(intersection, cellData);
- }
-
- //! Local velocity calculation
- void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
- {
- problem_.pressureModel().calculateVelocityOnBoundary(intersection, cellData);
- }
-
- //! \brief Indicates if velocity is reconstructed the transport step
- bool calculateVelocityInTransport()
- {
- return calcVelocityInTransport_;
- }
-
- /*! \brief Adds velocity output to the output file
- *
- * \tparam MultiWriter Class defining the output writer
- * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
- *
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {}
-
-private:
- Problem& problem_;
- bool calcVelocityInTransport_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/pressure.hh b/dumux/porousmediumflow/sequential/cellcentered/pressure.hh
deleted file mode 100644
index 5addb034d688a697fc3bbc924417e1b4a3e609fc..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/pressure.hh
+++ /dev/null
@@ -1,560 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_FVPRESSURE_HH
-#define DUMUX_FVPRESSURE_HH
-
-// dumux environment
-#include <type_traits>
-#include <dumux/common/math.hh>
-#include <dumux/common/exceptions.hh>
-#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
-#include <map>
-/**
- * @file
- * @brief Finite Volume Diffusion Model
- */
-
-namespace Dumux
-{
-//! The finite volume base class for the solution of a pressure equation
-/*! \ingroup IMPET
- * Base class for finite volume (FV) implementations of a diffusion-like pressure equation.
- * The class provides a methods for assembling of the global matrix and right hand side (RHS)
- * as well as for solving the system of equations.
- * Additionally, it contains the global matrix, the RHS-vector as well as the solution vector.
- * A certain pressure equation defined in the implementation of this base class must be splitted
- * into a storage term, a flux term and a source term.
- * Corresponding functions (<tt>getSource()</tt>, <tt>getStorage()</tt>, <tt>getFlux()</tt> and
- * <tt>getFluxOnBoundary()</tt>) have to be defined in the implementation.
- *
- * \tparam TypeTag The Type Tag
- */
-template<class TypeTag> class FVPressure
-{
- //the model implementation
- using Implementation = GetPropType<TypeTag, Properties::PressureModel>;
-
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using CellData = GetPropType<TypeTag, Properties::CellData>;
-
- // using declarations to abbreviate several dune classes...
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- // the typenames used for the stiffness matrix and solution vector
- using Matrix = GetPropType<TypeTag, Properties::PressureCoefficientMatrix>;
- using RHSVector = GetPropType<TypeTag, Properties::PressureRHSVector>;
- using PressureSolution = GetPropType<TypeTag, Properties::PressureSolutionVector>;
-
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
-
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
-
-protected:
-
- /*! Type of the vector of entries
- *
- * Contains the return values of the get*()-functions (matrix or right-hand side entry).
- */
- using EntryType = Dune::FieldVector<Scalar, 2>;
-
- //! Indices of matrix and rhs entries
- /**
- * During the assembling of the global system of equations get-functions are called (getSource(),
- * getFlux(), etc.), which return global matrix or right hand side entries in a vector.
- * These can be accessed using following indices:
- */
- enum
- {
- rhs = 1,//!<index for the right hand side entry
- matrix = 0//!<index for the global matrix entry
-
- };
-
- enum
- {
- pressEqIdx = Indices::pressureEqIdx,
- };
-
- //!Initialize the global matrix of the system of equations to solve
- void initializeMatrix();
- void initializeMatrixRowSize();
- void initializeMatrixIndices();
-
-
- /*!\brief Function which assembles the system of equations to be solved
- *
- * This function assembles the Matrix and the right hand side (RHS) vector to solve for
- * a pressure field with a Finite-Volume (FV) discretization.
- * Implementations of this base class have to provide the methods <tt>getSource()</tt>,
- * <tt>getStorage()</tt>, <tt>getFlux()</tt> and <tt>getFluxOnBoundary()</tt> if the assemble() method is called!
- *
- * \param first Indicates if function is called at the initialization step or during the simulation
- * (If <tt>first</tt> is <tt>true</tt>, no pressure field of previous iterations is required)
- */
- void assemble(bool first);
-
- //!Solves the global system of equations to get the spatial distribution of the pressure
- void solve();
-
- //!Returns the vector containing the pressure solution
- PressureSolution& pressure()
- { return pressure_;}
-
- //!Returns the vector containing the pressure solution
- const PressureSolution& pressure() const
- { return pressure_;}
-
- //!Initialization of the pressure solution vector: Initialization with meaningful values may
- //result in better convergence of the linear solver!
- void initializePressure()
- {
- for (const auto& element : elements(problem_.gridView()))
- {
- PrimaryVariables initValues;
- problem_.initial(initValues, element);
-
- int eIdxGlobal = problem_.variables().index(element);
-
- pressure_[eIdxGlobal] = initValues[pressEqIdx];
- }
- }
-
-public:
- /*! \brief Function which calculates the source entry
- *
- * Function computes the source term and writes it to the corresponding entry of the entry vector
- *
- * \param entry Vector containing return values of the function
- * \param element Grid element
- * \param cellData Object containing all model relevant cell data
- * \param first Indicates if function is called in the initialization step or during the simulation
- */
- void getSource(EntryType& entry, const Element& element, const CellData& cellData, const bool first);
-
- /*! \brief Function which calculates the storage entry
- *
- * Function computes the storage term and writes it to the corresponding entry of the entry vector
- *
- * \param entry Vector containing return values of the function
- * \param element Grid element
- * \param cellData Object containing all model relevant cell data
- * \param first Indicates if function is called in the initialization step or during the simulation
- */
- void getStorage(EntryType& entry, const Element& element, const CellData& cellData, const bool first);
-
- /*! \brief Function which calculates the flux entry
- *
- * Function computes the inter-cell flux term and writes it to the corresponding entry of the entry vector
- *
- * \param entry Vector containing return values of the function
- * \param intersection Intersection of two grid elements
- * \param cellData Object containing all model relevant cell data
- * \param first Indicates if function is called in the initialization step or during the simulation
- */
- void getFlux(EntryType& entry, const Intersection& intersection, const CellData& cellData, const bool first);
-
- /*! \brief Function which calculates the boundary flux entry
- *
- * Function computes the boundary-flux term and writes it to the corresponding entry of the entry vector
- *
- * \param entry Vector containing return values of the function
- * \param intersection Intersection of two grid elements
- * \param cellData Object containing all model relevant cell data
- * \param first Indicates if function is called in the initialization step or during the simulation
- */
- void getFluxOnBoundary(EntryType& entry,
- const Intersection& intersection, const CellData& cellData, const bool first);
-
- /*! \brief Public access function for the primary pressure variable
- *
- * Function returns the cell pressure value at index <tt>eIdxGlobal</tt>
- *
- * \param eIdxGlobal Global index of a grid cell
- */
- const Scalar pressure(const int eIdxGlobal) const
- { return pressure_[eIdxGlobal];}
-
- //!Returns the global matrix of the last pressure solution step.
- const Matrix& globalMatrix()
- {
- return A_;
- }
-
- //!Returns the right hand side vector of the last pressure solution step.
- const RHSVector& rightHandSide()
- {
- return f_;
- }
-
- /*! \brief Initialize pressure model
- *
- * Function initializes the sparse matrix to solve the global system of equations and sets/calculates the initial pressure
- */
- void initialize()
- {
- int size = problem_.gridView().size(0);//resize to make sure the final grid size (after the problem was completely built) is used!
- A_.setSize(size, size);
- A_.setBuildMode(Matrix::random);
- f_.resize(size);
- pressure_.resize(size);
- initializePressure();
- asImp_().initializeMatrix();// initialize sparse matrix
- }
-
- /*! \brief Pressure update
- *
- * Function reassembles the system of equations and solves for a new pressure solution.
- */
- void update()
- {
- asImp_().assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
- solve();
-
- return;
- }
-
- void calculateVelocity()
- {
- DUNE_THROW(Dune::NotImplemented,"Velocity calculation not implemented in pressure model!");
- }
-
- void updateVelocity()
- {} //empty function for the case the velocity is calculated in the transport model
-
- /*! \brief Function for serialization of the pressure field.
- *
- * Function needed for restart option. Writes the pressure of a grid element to a restart file.
- *
- * \param outstream Stream into the restart file.
- * \param element Grid element
- */
- void serializeEntity(std::ostream &outstream, const Element &element)
- {
- int eIdxGlobal = problem_.variables().index(element);
- outstream << pressure_[eIdxGlobal][0];
- }
-
- /*! \brief Function for deserialization of the pressure field.
- *
- * Function needed for restart option. Reads the pressure of a grid element from a restart file.
- *
- * \param instream Stream from the restart file.
- * \param element Grid element
- */
- void deserializeEntity(std::istream &instream, const Element &element)
- {
- int eIdxGlobal = problem_.variables().index(element);
- instream >> pressure_[eIdxGlobal][0];
- }
-
- /*! \brief Set a pressure to be fixed at a certain cell.
- *
- *Allows to fix a pressure somewhere (at one certain cell) in the domain.
- *This can be necessary e.g. if only Neumann boundary conditions are defined.
- *The pressure is fixed until the <tt>unsetFixPressureAtIndex()</tt> function is called
- *
- * \param pressure Pressure value at eIdxGlobal
- * \param eIdxGlobal Global index of a grid cell
- */
- void setFixPressureAtIndex(Scalar pressure, int eIdxGlobal)
- {
- fixPressure_.insert(std::make_pair(eIdxGlobal, pressure));
- }
-
- /*! \brief Reset the fixed pressure state
- *
- * No pressure is fixed inside the domain until <tt>setFixPressureAtIndex()</tt> function is called again.
- *
- * \param eIdxGlobal Global index of a grid cell
- */
- void unsetFixPressureAtIndex(int eIdxGlobal)
- {
- fixPressure_.erase(eIdxGlobal);
- }
-
- void resetFixPressureAtIndex()
- {
- fixPressure_.clear();
- }
-
- /*! \brief Constructs a FVPressure object
- *
- * \param problem A problem class object
- */
- FVPressure(Problem& problem) :
- problem_(problem)
- {}
-
-private:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this);}
-
- //! \copydoc IMPETProblem::asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this);}
-
- Problem& problem_;
-
- PressureSolution pressure_;
-
- std::string solverName_;
- std::string preconditionerName_;
-protected:
- Matrix A_;//!<Global stiffnes matrix (sparse matrix which is build by the <tt> initializeMatrix()</tt> function)
- RHSVector f_;//!<Right hand side vector
-private:
- std::map<int, Scalar> fixPressure_;
-};
-
-//!Initialize the global matrix of the system of equations to solve
-template<class TypeTag>
-void FVPressure<TypeTag>::initializeMatrix()
-{
- initializeMatrixRowSize();
- A_.endrowsizes();
- initializeMatrixIndices();
- A_.endindices();
-}
-
-//!Initialize the global matrix of the system of equations to solve
-template<class TypeTag>
-void FVPressure<TypeTag>::initializeMatrixRowSize()
-{
- // determine matrix row sizes
- for (const auto& element : elements(problem_.gridView()))
- {
- // cell index
- int eIdxGlobalI = problem_.variables().index(element);
-
- // initialize row size
- int rowSize = 1;
-
- // run through all intersections with neighbors
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- if (intersection.neighbor())
- rowSize++;
- }
- A_.setrowsize(eIdxGlobalI, rowSize);
- }
-}
-
-//!Initialize the global matrix of the system of equations to solve
-template<class TypeTag>
-void FVPressure<TypeTag>::initializeMatrixIndices()
-{
- // determine position of matrix entries
- for (const auto& element : elements(problem_.gridView()))
- {
- // cell index
- int eIdxGlobalI = problem_.variables().index(element);
-
- // add diagonal index
- A_.addindex(eIdxGlobalI, eIdxGlobalI);
-
- // run through all intersections with neighbors
- for (const auto& intersection : intersections(problem_.gridView(), element))
- if (intersection.neighbor())
- {
- // access neighbor
- int eIdxGlobalJ = problem_.variables().index(intersection.outside());
-
- // add off diagonal index
- A_.addindex(eIdxGlobalI, eIdxGlobalJ);
- }
- }
-}
-
-
-/*!\brief Function which assembles the system of equations to be solved
- *
- * This function assembles the Matrix and the right hand side (RHS) vector to solve for
- * a pressure field with a Finite-Volume (FV) discretization.
- * Implementations of this base class have to provide the methods <tt>getSource()</tt>,
- * <tt>getStorage()</tt>, <tt>getFlux()</tt> and <tt>getFluxOnBoundary()</tt> if the assemble() method is called!
- *
- * \param first Indicates if function is called at the initialization step or during the simulation
- * (If <tt>first</tt> is <tt>true</tt>, no pressure field of previous iterations is required)
- */
-template<class TypeTag>
-void FVPressure<TypeTag>::assemble(bool first)
-{
- // initialization: set matrix A_ to zero
- A_ = 0;
- f_ = 0;
-
- for (const auto& element : elements(problem_.gridView()))
- {
- // get the global index of the cell
- int eIdxGlobalI = problem_.variables().index(element);
-
- // assemble interior element contributions
- if (element.partitionType() == Dune::InteriorEntity)
- {
- // get the cell data
- CellData& cellDataI = problem_.variables().cellData(eIdxGlobalI);
-
- EntryType entries(0.);
-
- /***** source term ***********/
- asImp_().getSource(entries, element, cellDataI, first);
- f_[eIdxGlobalI] += entries[rhs];
-
- /***** flux term ***********/
- // iterate over all faces of the cell
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- /************* handle interior face *****************/
- if (intersection.neighbor())
- {
- auto elementNeighbor = intersection.outside();
-
- int eIdxGlobalJ = problem_.variables().index(elementNeighbor);
-
- // check for hanging nodes
- bool haveSameLevel = (element.level() == elementNeighbor.level());
- // calculate only from one side (except for hanging nodes), but add matrix entries for both sides
- // the last condition is needed to properly assemble in the presence
- // of ghost elements
- if (getPropValue<TypeTag, Properties::VisitFacesOnlyOnce>()
- && (eIdxGlobalI > eIdxGlobalJ) && haveSameLevel
- && elementNeighbor.partitionType() == Dune::InteriorEntity)
- continue;
-
- entries = 0;
- asImp_().getFlux(entries, intersection, cellDataI, first);
-
- //set right hand side
- f_[eIdxGlobalI] -= entries[rhs];
-
- // set diagonal entry
- A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
-
- // set off-diagonal entry
- A_[eIdxGlobalI][eIdxGlobalJ] -= entries[matrix];
-
- // The second condition is needed to not spoil the ghost element entries
- if (getPropValue<TypeTag, Properties::VisitFacesOnlyOnce>()
- && elementNeighbor.partitionType() == Dune::InteriorEntity)
- {
- f_[eIdxGlobalJ] += entries[rhs];
- A_[eIdxGlobalJ][eIdxGlobalJ] += entries[matrix];
- A_[eIdxGlobalJ][eIdxGlobalI] -= entries[matrix];
- }
-
- } // end neighbor
-
- /************* boundary face ************************/
- else
- {
- entries = 0;
- asImp_().getFluxOnBoundary(entries, intersection, cellDataI, first);
-
- //set right hand side
- f_[eIdxGlobalI] += entries[rhs];
- // set diagonal entry
- A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
- }
- } //end interfaces loop
- // printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
-
- /***** storage term ***********/
- entries = 0;
- asImp_().getStorage(entries, element, cellDataI, first);
- f_[eIdxGlobalI] += entries[rhs];
- // set diagonal entry
- A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
- }
- // assemble overlap and ghost element contributions
- else
- {
- A_[eIdxGlobalI] = 0.0;
- A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
- f_[eIdxGlobalI] = pressure_[eIdxGlobalI];
- }
- } // end grid traversal
-// printmatrix(std::cout, A_, "global stiffness matrix after assempling", "row", 11,3);
-// printvector(std::cout, f_, "right hand side", "row", 10);
-}
-
-// forward declaration
-template<class T>
-class AMGBiCGSTABBackend;
-
-namespace Detail {
-template<class T> struct isParallelAMGBackend : public std::false_type {};
-template<class T>
-struct isParallelAMGBackend<Dumux::AMGBiCGSTABBackend<T>> : public std::true_type {};
-} // end namespace Detail
-
-template<class Solver, class Problem>
-typename std::enable_if_t<!Detail::isParallelAMGBackend<Solver>::value, Solver>
-getSolver(const Problem& problem)
-{
- return Solver();
-}
-
-template<class Solver, class Problem>
-typename std::enable_if_t<Detail::isParallelAMGBackend<Solver>::value, Solver>
-getSolver(const Problem& problem)
-{
- return Solver(problem.gridView(), problem.model().dofMapper());
-}
-
-//!Solves the global system of equations to get the spatial distribution of the pressure
-template<class TypeTag>
-void FVPressure<TypeTag>::solve()
-{
- using Solver = GetPropType<TypeTag, Properties::LinearSolver>;
-
- int verboseLevelSolver = getParam<int>("LinearSolver.Verbosity", 0);
-
- if (verboseLevelSolver)
- std::cout << __FILE__ << ": solve for pressure" << std::endl;
-
- //set a fixed pressure for a certain cell
- if (fixPressure_.size() > 0)
- {
- for (auto it = fixPressure_.begin(); it != fixPressure_.end(); ++it)
- {
- A_[it->first] = 0;
- A_[it->first][it->first] = 1;
- f_[it->first] = it->second;
- }
- }
-
-// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
-// printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
-
- auto solver = getSolver<Solver>(problem_);
- bool converged = solver.solve(A_, pressure_, f_);
-
- if (!converged)
- DUNE_THROW(Dumux::NumericalProblem, "Pressure solver did not converge!");
-
-// printvector(std::cout, pressure_, "pressure", "row", 200, 1, 3);
-}
-
-} //end namespace Dumux
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/transport.hh b/dumux/porousmediumflow/sequential/cellcentered/transport.hh
deleted file mode 100644
index fe564c699c737614b938cb0f535a799c26e6708c..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/transport.hh
+++ /dev/null
@@ -1,649 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_FVTRANSPORT_HH
-#define DUMUX_FVTRANSPORT_HH
-
-#include <dune/grid/common/gridenums.hh>
-#include <dumux/porousmediumflow/sequential/transportproperties.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-#include <dumux/parallel/vectorcommdatahandle.hh>
-#include <unordered_map>
-
-/**
- * @file
- * @brief Finite Volume discretization of a transport equation
- */
-
-namespace Dumux
-{
-//! \ingroup IMPET
-/*!\brief The finite volume discretization of a transport equation
- *
- * Base class for finite volume (FV) implementations of an explicitly treated transport equation.
- * The class provides a method to calculate the explicit update to get a new solution of the transported quantity:
- * \f[
- * u_{new} = u_{old} + \Delta t \Delta u_{update}
- * \f]
- * A certain transport equation defined in a implementation of this base class must be splitted
- * into a flux term and a source term.
- * Corresponding functions (<tt>getSource()</tt>, <tt>getFlux()</tt> and <tt>getFluxOnBoundary()</tt>)
- * have to be defined in the implementation.
- *
- * \tparam TypeTag The Type Tag
- */
-template<class TypeTag>
-class FVTransport
-{
- using Implementation = GetPropType<TypeTag, Properties::TransportModel>;
-
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- enum
- {
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
- };
-
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
-
- using TransportSolutionType = GetPropType<TypeTag, Properties::TransportSolutionType>;
- using CellData = GetPropType<TypeTag, Properties::CellData>;
-
- using EvalCflFluxFunction = GetPropType<TypeTag, Properties::EvalCflFluxFunction>;
-
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-
- struct LocalTimesteppingData
- {
- Dune::FieldVector<Scalar, 2*dim> faceFluxes;
- Dune::FieldVector<Scalar, 2*dim> faceTargetDt;
- Scalar dt;
- LocalTimesteppingData():faceFluxes(0.0), faceTargetDt(0.0), dt(0)
- {}
- };
-
-protected:
- //! \cond \private
- EvalCflFluxFunction& evalCflFluxFunction()
- {
- return *evalCflFluxFunction_;
- }
-
- const EvalCflFluxFunction& evalCflFluxFunction() const
- {
- return *evalCflFluxFunction_;
- }
- //! \endcond
-
- void innerUpdate(TransportSolutionType& updateVec);
-
-public:
-
- // Calculate the update vector.
- void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet = false);
-
- void updateTransport(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
- {
- asImp_().updateMaterialLaws();
- asImp_().update(t, dt, updateVec);
- }
-
- /*! \brief Function which calculates the flux update
- *
- * Function computes the inter-cell flux term and adds it to the update.
- *
- * \param update The cell update
- * \param intersection Intersection of two grid elements
- * \param cellDataI Object containing all model relevant cell data
- */
- void getFlux(Scalar& update, const Intersection& intersection, CellData& cellDataI);
-
- /*! \brief Function which calculates the boundary flux update
- *
- * Function computes the boundary-flux term and adds it to the update.
- *
- * \param update The cell update
- * \param intersection Intersection of two grid elements
- * \param cellDataI Object containing all model relevant cell data
- */
- void getFluxOnBoundary(Scalar& update, const Intersection& intersection, CellData& cellDataI);
-
- /*! \brief Function which calculates the source update
- *
- * Function computes the source term and adds it to the update.
- *
- * \param update The cell update
- * \param element Grid element
- * \param cellDataI Object containing all model relevant cell data
- */
- void getSource(Scalar& update, const Element& element, CellData& cellDataI);
-
- //! Sets the initial solution \f$ S_0 \f$.
- void initialize()
- {
- evalCflFluxFunction_->initialize();
- }
-
- /*! \brief Updates constitutive relations and stores them in the variable class*/
- void updateMaterialLaws();
-
- /*! \brief Writes the current values of the primary transport variable into the
- * <tt>transportedQuantity</tt>-vector (comes as function argument)
- *
- * \param transportedQuantity Vector of the size of global numbers of degrees of freedom
- * of the primary transport variable.
- */
- void getTransportedQuantity(TransportSolutionType& transportedQuantity);
-
- template<class DataEntry>
- bool inPhysicalRange(DataEntry& entry);
-
- /*! \brief Writes the current values of the primary transport variable into the variable container
- *
- * \param transportedQuantity Vector of the size of global numbers of degrees of freedom of the primary transport variable.
- */
- void setTransportedQuantity(TransportSolutionType& transportedQuantity);
-
- /*! \brief Updates the primary transport variable.
- *
- * \param updateVec Vector containing the global update.
- */
- void updateTransportedQuantity(TransportSolutionType& updateVec);
-
- /*! \brief Adds transport output to the output file
- *
- * \tparam MultiWriter Class defining the output writer
- * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
- *
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {}
-
- /*! \brief Function for serialization of the primary transport variable.
- *
- * Function needed for restart option. Writes the primary transport variable of a grid element to a restart file.
- *
- * \param outstream Stream into the restart file.
- * \param element Grid element
- */
- void serializeEntity(std::ostream &outstream, const Element &element)
- {}
-
- /*! \brief Function for deserialization of the primary transport variable.
- *
- * Function needed for restart option. Reads the the primary transport variable of a grid element from a restart file.
- *
- * \param instream Stream from the restart file.
- * \param element Grid element
- */
- void deserializeEntity(std::istream &instream, const Element &element)
- {}
-
- bool enableLocalTimeStepping()
- {
- return localTimeStepping_;
- }
-
-
- //! Constructs a FVTransport object
- /**
-
- * \param problem A problem class object
- */
- FVTransport(Problem& problem) :
- problem_(problem), switchNormals_(getParam<bool>("Impet.SwitchNormals")),
- subCFLFactor_(1.0), accumulatedDt_(0), dtThreshold_(1e-6)
- {
- evalCflFluxFunction_ = std::make_shared<EvalCflFluxFunction>(problem);
-
- Scalar cFLFactor = getParam<Scalar>("Impet.CFLFactor");
- using std::min;
- subCFLFactor_ = min(getParam<Scalar>("Impet.SubCFLFactor"), cFLFactor);
- verbosity_ = getParam<int>("TimeManager.SubTimestepVerbosity");
-
- localTimeStepping_ = subCFLFactor_/cFLFactor < 1.0 - dtThreshold_;
-
- if (localTimeStepping_)
- std::cout<<"max CFL-Number of "<<cFLFactor<<", max Sub-CFL-Number of "<<subCFLFactor_<<": Enable local time-stepping!\n";
- }
-
-private:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \copydoc IMPETProblem::asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
-
- void updatedTargetDt_(Scalar &dt);
-
- void resetTimeStepData_()
- {
- timeStepData_.clear();
- accumulatedDt_ = 0;
- }
-
- Problem& problem_;
- bool switchNormals_;
-
- std::shared_ptr<EvalCflFluxFunction> evalCflFluxFunction_;
- std::vector<LocalTimesteppingData> timeStepData_;
- bool localTimeStepping_;
- Scalar subCFLFactor_;
- Scalar accumulatedDt_;
- const Scalar dtThreshold_;
- int verbosity_;
-};
-
-
-/*! \brief Calculate the update vector.
- * \param t current time
- * \param dt time step size
- * \param updateVec vector containing the update values
- * \param impet variable should be true if an impet algorithm is used and false if the transport part is solved independently
- *
- * Additionally to the \a update vector, the recommended time step size \a dt is calculated
- * employing a CFL condition.
- */
-template<class TypeTag>
-void FVTransport<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impet)
-{
- if (!impet)
- {
- asImp_().updateMaterialLaws();
- }
-
- unsigned int size = problem_.gridView().size(0);
- if (localTimeStepping_)
- {
- if (timeStepData_.size() != size)
- timeStepData_.resize(size);
- }
- // initialize dt very large
- dt = std::numeric_limits<Scalar>::max();
-
- // resize update vector and set to zero
- updateVec.resize(size);
- updateVec = 0.0;
-
- // compute update vector
- for (const auto& element : elements(problem_.gridView()))
- {
-#if HAVE_MPI
- if (element.partitionType() != Dune::InteriorEntity)
- {
- continue;
- }
-#endif
-
- // cell index
- int globalIdxI = problem_.variables().index(element);
-
- CellData& cellDataI = problem_.variables().cellData(globalIdxI);
-
- Scalar update = 0;
- evalCflFluxFunction().reset();
-
- if (localTimeStepping_)
- {
- LocalTimesteppingData& localData = timeStepData_[globalIdxI];
- for (int i = 0; i < 2*dim; i++)
- {
- if (localData.faceTargetDt[i] < accumulatedDt_ + dtThreshold_)
- {
- localData.faceFluxes[i] = 0.0;
- }
- }
- }
-
- // run through all intersections with neighbors and boundary
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
- if (switchNormals_)
- unitOuterNormal *= -1.0;
-
- int indexInInside = intersection.indexInInside();
-
- // handle interior face
- if (intersection.neighbor())
- {
- if (localTimeStepping_)
- {
- LocalTimesteppingData& localData = timeStepData_[globalIdxI];
-
- if (localData.faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_)
- {
- asImp_().getFlux(localData.faceFluxes[indexInInside], intersection, cellDataI);
- }
- else
- {
- asImp_().getFlux(update, intersection, cellDataI);//only for time-stepping
- }
- }
- else
- {
- //add flux to update
- asImp_().getFlux(update, intersection, cellDataI);
- }
- } //end intersection with neighbor element
- // handle boundary face
- else if (intersection.boundary())
- {
- if (localTimeStepping_)
- {
- LocalTimesteppingData& localData = timeStepData_[globalIdxI];
- if (localData.faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_)
- {
- asImp_().getFluxOnBoundary(localData.faceFluxes[indexInInside], intersection, cellDataI);
- }
- else
- {
- asImp_().getFluxOnBoundary(update, intersection, cellDataI);//only for time-stepping
- }
- }
- else
- {
- //add boundary flux to update
- asImp_().getFluxOnBoundary(update, intersection, cellDataI);
- }
- } //end boundary
- } // end all intersections
-
- if (localTimeStepping_)
- {
- LocalTimesteppingData& localData = timeStepData_[globalIdxI];
- for (int i=0; i < 2*dim; i++)
- {
- updateVec[globalIdxI] += localData.faceFluxes[i];
- }
- }
- else
- {
- //add flux update to global update vector
- updateVec[globalIdxI] += update;
- }
-
- // std::cout<<"updateVec at "<<element.geometry().center()<<" : "<<updateVec[globalIdxI]<<"\n";
-
- //add source to global update vector
- Scalar source = 0.;
- asImp_().getSource(source,element, cellDataI);
- updateVec[globalIdxI] += source;
-
- //calculate time step
- using std::min;
- if (localTimeStepping_)
- {
- Scalar dtCfl = evalCflFluxFunction().getDt(element);
-
- timeStepData_[globalIdxI].dt = dtCfl;
- dt = min(dt, dtCfl);
- }
- else
- {
- //calculate time step
- dt = min(dt, evalCflFluxFunction().getDt(element));
- }
-
- //store update
- cellDataI.setUpdate(updateVec[globalIdxI]);
- } // end grid traversal
-
-
-#if HAVE_MPI
- // communicate updated values
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using ElementMapper = typename SolutionTypes::ElementMapper;
- using DataHandle = VectorCommDataHandleEqual<ElementMapper, Dune::BlockVector<Dune::FieldVector<Scalar, 1> >, 0/*elementCodim*/>;
- DataHandle dataHandle(problem_.elementMapper(), updateVec);
- problem_.gridView().template communicate<DataHandle>(dataHandle,
- Dune::InteriorBorder_All_Interface,
- Dune::ForwardCommunication);
-
- if (localTimeStepping_)
- {
- using TimeDataHandle = VectorCommDataHandleEqual<ElementMapper, std::vector<LocalTimesteppingData>, 0/*elementCodim*/>;
-
- TimeDataHandle timeDataHandle(problem_.elementMapper(), timeStepData_);
- problem_.gridView().template communicate<TimeDataHandle>(timeDataHandle,
- Dune::InteriorBorder_All_Interface,
- Dune::ForwardCommunication);
- }
-
- dt = problem_.gridView().comm().min(dt);
-#endif
-}
-
-template<class TypeTag>
-void FVTransport<TypeTag>::updatedTargetDt_(Scalar &dt)
-{
- dt = std::numeric_limits<Scalar>::max();
-
- // update target time-step-sizes
- for (const auto& element : elements(problem_.gridView()))
- {
-#if HAVE_MPI
- if (element.partitionType() != Dune::InteriorEntity)
- {
- continue;
- }
-#endif
-
- // cell index
- int globalIdxI = problem_.variables().index(element);
-
- LocalTimesteppingData& localDataI = timeStepData_[globalIdxI];
-
-
- using FaceDt = std::unordered_map<int, Scalar>;
- FaceDt faceDt;
-
- // run through all intersections with neighbors and boundary
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- int indexInInside = intersection.indexInInside();
- using std::min;
- if (intersection.neighbor())
- {
- auto neighbor = intersection.outside();
- int globalIdxJ = problem_.variables().index(neighbor);
-
- int levelI = element.level();
- int levelJ = neighbor.level();
-
- if (globalIdxI < globalIdxJ && levelI <= levelJ)
- {
- LocalTimesteppingData& localDataJ = timeStepData_[globalIdxJ];
-
- int indexInOutside = intersection.indexInOutside();
-
- if (localDataI.faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_
- || localDataJ.faceTargetDt[indexInOutside] < accumulatedDt_ + dtThreshold_)
- {
- Scalar timeStep = min(localDataI.dt, localDataJ.dt);
-
- if (levelI < levelJ)
- {
- typename FaceDt::iterator it = faceDt.find(indexInInside);
- if (it != faceDt.end())
- {
- it->second = min(it->second, timeStep);
- }
- else
- {
- faceDt.insert(std::make_pair(indexInInside, timeStep));
- }
- }
- else
- {
- localDataI.faceTargetDt[indexInInside] += subCFLFactor_ * timeStep;
- localDataJ.faceTargetDt[indexInOutside] += subCFLFactor_ * timeStep;
- }
-
- dt = min(dt, timeStep);
- }
- }
- }
- else if (intersection.boundary())
- {
- if (localDataI.faceTargetDt[indexInInside] < accumulatedDt_ + dtThreshold_)
- {
- localDataI.faceTargetDt[indexInInside] += subCFLFactor_ * localDataI.dt;
- dt =min(dt, subCFLFactor_ * localDataI.dt);
- }
- }
- }
- if (faceDt.size() > 0)
- {
- typename FaceDt::iterator it = faceDt.begin();
- for (;it != faceDt.end();++it)
- {
- localDataI.faceTargetDt[it->first] += subCFLFactor_ * it->second;
- }
-
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- if (intersection.neighbor())
- {
- int indexInInside = intersection.indexInInside();
-
- it = faceDt.find(indexInInside);
- if (it != faceDt.end())
- {
- int globalIdxJ = problem_.variables().index(intersection.outside());
-
- LocalTimesteppingData& localDataJ = timeStepData_[globalIdxJ];
-
- int indexInOutside = intersection.indexInOutside();
-
- localDataJ.faceTargetDt[indexInOutside] += subCFLFactor_ * it->second;
- }
- }
- }
- }
- }
-
-#if HAVE_MPI
- // communicate updated values
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using ElementMapper = typename SolutionTypes::ElementMapper;
- using TimeDataHandle = VectorCommDataHandleEqual<ElementMapper, std::vector<LocalTimesteppingData>, 0/*elementCodim*/>;
-
- TimeDataHandle timeDataHandle(problem_.elementMapper(), timeStepData_);
- problem_.gridView().template communicate<TimeDataHandle>(timeDataHandle,
- Dune::InteriorBorder_All_Interface,
- Dune::ForwardCommunication);
-
- dt = problem_.gridView().comm().min(dt);
-#endif
-}
-
-template<class TypeTag>
-void FVTransport<TypeTag>::innerUpdate(TransportSolutionType& updateVec)
-{
- if (localTimeStepping_)
- {
- Scalar realDt = problem_.timeManager().timeStepSize();
-
- Scalar subDt = realDt;
-
- updatedTargetDt_(subDt);
-
- Scalar accumulatedDtOld = accumulatedDt_;
- accumulatedDt_ += subDt;
-
- Scalar t = problem_.timeManager().time();
-
- if (accumulatedDt_ < realDt)
- {
- using std::min;
- while(true)
- {
- Scalar dtCorrection = min(0.0, realDt - accumulatedDt_);
- subDt += dtCorrection;
-
- if (verbosity_ > 0)
- std::cout<<" Sub-time-step size: "<<subDt<<"\n";
-
- TransportSolutionType transportedQuantity;
- asImp_().getTransportedQuantity(transportedQuantity);
-
- bool stopTimeStep = false;
- int size = transportedQuantity.size();
- for (int i = 0; i < size; i++)
- {
- Scalar newVal = transportedQuantity[i] += updateVec[i] * subDt;
- if (!asImp_().inPhysicalRange(newVal))
- {
- stopTimeStep = true;
-
- break;
- }
- }
-
-#if HAVE_MPI
- int rank = 0;
- if (stopTimeStep)
- rank = problem_.gridView().comm().rank();
-
- rank = problem_.gridView().comm().max(rank);
- problem_.gridView().comm().broadcast(&stopTimeStep,1,rank);
-#endif
-
-
- if (stopTimeStep && accumulatedDtOld > dtThreshold_)
- {
- problem_.timeManager().setTimeStepSize(accumulatedDtOld);
- break;
- }
- else
- {
- asImp_().setTransportedQuantity(transportedQuantity);
- }
-
-
- if (accumulatedDt_ >= realDt)
- {
- break;
- }
-
- problem_.model().updateTransport(t, subDt, updateVec);
-
- updatedTargetDt_(subDt);
-
- accumulatedDtOld = accumulatedDt_;
- accumulatedDt_ += subDt;
- }
- }
- else
- {
- asImp_().updateTransportedQuantity(updateVec, realDt);
- }
-
- resetTimeStepData_();
- }
-}
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/velocity.hh b/dumux/porousmediumflow/sequential/cellcentered/velocity.hh
deleted file mode 100644
index 38b999d84b0d7f3d6c6fad12f8031558f27e5b22..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/velocity.hh
+++ /dev/null
@@ -1,129 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_FVVELOCITY_HH
-#define DUMUX_FVVELOCITY_HH
-
-// dumux environment
-#include <dumux/common/math.hh>
-#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
-#include "velocitydefault.hh"
-
-/**
- * @file
- * @brief Finite volume velocity reconstruction
- */
-
-namespace Dumux
-{
-
-/*! \ingroup IMPET
- *
- * \brief Base class for finite volume velocity reconstruction
- *
- * Provides a basic frame for calculating a global velocity field.
- * The definition of the local velocity calculation as well as the storage or other postprocessing
- * has to be provided by the local velocity implementation.
- * This local implementation has to have the form of VelocityDefault.
- *
- * \tparam TypeTag The Type Tag
- * \tparam Velocity The implementation of the local velocity calculation
- */
-template<class TypeTag, class Velocity> class FVVelocity
-{
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
-
- using CellData = GetPropType<TypeTag, Properties::CellData>;
-public:
-
- //!Initialize velocity implementation
- void initialize()
- {
- velocity_.initialize();
- }
-
- //function which iterates through the grid and calculates the global velocity field
- void calculateVelocity();
-
- /*! \brief Adds velocity output to the output file
- *
- * \tparam MultiWriter Class defining the output writer
- * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
- *
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {
- velocity_.addOutputVtkFields(writer);
- }
-
- //! Constructs a FVVelocity object
- /**
- * \param problem A problem class object
- */
- FVVelocity(Problem& problem) :
- problem_(problem), velocity_(problem)
- {}
-
-private:
- Problem& problem_;
- Velocity velocity_;
-};
-
-
-/*! \brief Function which reconstructs a global velocity field
- *
- * Iterates through the grid and calls the local calculateVelocity(...) or calculateVelocityOnBoundary(...)
- * functions which have to be provided by the local velocity implementation (see e.g. VelocityDefault )
- */
-template<class TypeTag, class Velocity>
-void FVVelocity<TypeTag, Velocity>::calculateVelocity()
-{
- for (const auto& element : elements(problem_.gridView()))
- {
- // cell information
- int globalIdxI = problem_.variables().index(element);
- CellData& cellDataI = problem_.variables().cellData(globalIdxI);
-
- /***** flux term ***********/
- // iterate over all faces of the cell
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- /************* handle interior face *****************/
- if (intersection.neighbor())
- {
- int isIndex = intersection.indexInInside();
-
- if (!cellDataI.fluxData().haveVelocity(isIndex))
- velocity_.calculateVelocity(intersection, cellDataI);
- } // end neighbor
-
- /************* boundary face ************************/
- else
- {
- velocity_.calculateVelocityOnBoundary(intersection, cellDataI);
- }
- } //end interfaces loop
- } // end grid traversal
-
- return;
-}
-
-}//end namespace Dumux
-#endif
diff --git a/dumux/porousmediumflow/sequential/cellcentered/velocitydefault.hh b/dumux/porousmediumflow/sequential/cellcentered/velocitydefault.hh
deleted file mode 100644
index c8b790d3696e2363d5f2647394c7f90ad78b5c9e..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/cellcentered/velocitydefault.hh
+++ /dev/null
@@ -1,89 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_VELOCITYDEFAULT_HH
-#define DUMUX_VELOCITYDEFAULT_HH
-
-/**
- * @file
- * @brief Default implementation of velocity class.
- */
-
-#include <dune/grid/common/gridenums.hh>
-
-#include <dumux/porousmediumflow/sequential/properties.hh>
-
-namespace Dumux
-{
-//! \ingroup IMPET
-/*! \brief Default implementation of a velocity class.
- *
- * If the velocity is reconstructed in the pressure model this default implementation is used in the transport model.
- *
- * \tparam TypeTag The Type Tag
- */
-template<class TypeTag>
-class FVVelocityDefault
-{
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- using CellData = GetPropType<TypeTag, Properties::CellData>;
- using Intersection = typename GridView::Intersection;
-
-
-public:
- //! Constructs a FVVelocityDefault object
- /*!
- * \param problem A problem class object
- */
- FVVelocityDefault(Problem& problem)
- {}
-
- //! For initialization
- void initialize()
- {}
-
- //! Local velocity calculation
- void calculateVelocity(const Intersection& intersection, CellData& cellData)
- {}
-
- //! Local velocity calculation
- void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
- {}
-
- //! \brief Indicates if velocity is reconstructed the transport step
- bool calculateVelocityInTransport()
- {
- return false;
- }
-
- /*! \brief Adds velocity output to the output file
- *
- * \tparam MultiWriter Class defining the output writer
- * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
- *
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {}
-
-
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/gridadapt.hh b/dumux/porousmediumflow/sequential/gridadapt.hh
deleted file mode 100644
index 3818baa59510476c756e11b831b3b9b88a473758..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/gridadapt.hh
+++ /dev/null
@@ -1,492 +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 3 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
- * \brief Base class for h-adaptive sequential models.
- */
-#ifndef DUMUX_GRIDADAPT_HH
-#define DUMUX_GRIDADAPT_HH
-
-#include <unordered_map>
-#include <dune/grid/common/partitionset.hh>
-#include <dumux/common/deprecated.hh>
-
-#include "properties.hh"
-#include "gridadaptproperties.hh"
-
-namespace Dumux
-{
-
-/*!\ingroup IMPET
- * @brief Standard Module for h-adaptive simulations
- *
- * This class is created by the problem class with the template
- * parameters <TypeTag, true> and provides basic functionality
- * for adaptive methods:
- *
- * A standard implementation adaptGrid() will prepare everything
- * to calculate the next pressure field on the new grid.
- */
-template<class TypeTag, bool adaptive>
-class GridAdapt
-{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-
- using Grid = typename GridView::Grid;
- using LeafGridView = typename Grid::LeafGridView;
- using Element = typename Grid::template Codim<0>::Entity;
-
- using CellData = GetPropType<TypeTag, Properties::CellData>;
- using AdaptionIndicator = GetPropType<TypeTag, Properties::AdaptionIndicator>;
- using AdaptionInitializationIndicator = GetPropType<TypeTag, Properties::AdaptionInitializationIndicator>;
-
- using IdSet = typename Grid::Traits::LocalIdSet;
- using IdType = typename IdSet::IdType;
-
-public:
- /*!
- * Constructor for h-adaptive simulations (adaptive grids)
- * @param problem The problem
- */
- GridAdapt (Problem& problem)
- : problem_(problem), adaptionIndicator_(problem), marked_(0), coarsened_(0)
- {
- levelMin_ = getParam<int>("GridAdapt.MinLevel");
- levelMax_ = getParam<int>("GridAdapt.MaxLevel");
- adaptationInterval_ = getParam<int>("GridAdapt.AdaptionInterval", 1);
-
- if (levelMin_ < 0)
- Dune::dgrave << __FILE__<< ":" <<__LINE__
- << " : Dune cannot coarsen to gridlevels smaller 0! "<< std::endl;
- }
-
- /*!
- * @brief Initalization method of the h-adaptive module
- *
- * Prepares the grid for simulation after the initialization of the
- * problem. The applied indicator is selectable via the property
- * AdaptionInitializationIndicator
- */
- void init()
- {
- adaptionIndicator_.init();
-
- if (!getParam<bool>("GridAdapt.EnableInitializationIndicator"))
- return;
-
- AdaptionInitializationIndicator adaptionInitIndicator(problem_, adaptionIndicator_);
-
- int maxIter = 2*levelMax_;
- int iter = 0;
- while (iter <= maxIter)
- {
- adaptGrid(adaptionInitIndicator);
-
- if (!wasAdapted())
- {
- break;
- }
-
- int shouldInitialize = adaptionInitIndicator.initializeModel();
- if (problem_.grid().comm().max(shouldInitialize))
- problem_.model().initialize();
-
- iter++;
- }
- }
-
- /*!
- * @brief Standard method to adapt the grid
- *
- * This method is called from IMPETProblem::preTimeStep() if
- * adaptive grids are used in the simulation. It uses the standard
- * indicator (selected by the property AdaptionIndicator) and forwards to
- * with it to the ultimate method adaptGrid(indicator), which
- * uses a standard procedure for adaptivity:
- * 1) Determine the refinement indicator
- * 2) Mark the elements
- * 3) Store primary variables in a map
- * 4) Adapt the grid, adapt variables sizes, update mappers
- * 5) Reconstruct primary variables, regain secondary variables
- */
- void adaptGrid()
- {
- adaptGrid(adaptionIndicator_) ;
- }
-
- /*!
- * @brief Method to adapt the grid with individual indicator vector
- *
- * @param indicator The refinement indicator that is applied
- *
- * This method is called by an user-defined preTimeStep() of
- * the applied problem and takes a given vector with indicator
- * values.
- *
- * It uses a standard procedure for adaptivity:
- * 1) Determine the refinement indicator
- * 2) Mark the elements
- * 3) Store primary variables in a map
- * 4) Adapt the grid, adapt variables sizes, update mappers
- * 5) Reconstruct primary variables, regain secondary variables
- */
- template<class Indicator>
- void adaptGrid(Indicator& indicator)
- {
- // reset internal counter for marked elements
- marked_ = coarsened_ = 0;
-
- // check for adaption interval: Adapt only at certain time step indices
- if (problem_.timeManager().timeStepIndex() % adaptationInterval_ != 0)
- return;
-
- /**** 1) determine refining parameter if standard is used ***/
- // if not, the indicatorVector and refinement Bounds have to
- // specified by the problem through setIndicator()
- indicator.calculateIndicator();
-
- /**** 2) mark elements according to indicator *********/
- markElements(indicator);
-
- // abort if nothing in grid is marked
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
- if (sumMarked == 0 && sumCoarsened == 0)
- return;
- else
- Dune::dinfo << marked_ << " cells have been marked_ to be refined, "
- << coarsened_ << " to be coarsened." << std::endl;
-
- /**** 2b) Do pre-adaption step *****/
- problem_.grid().preAdapt();
- problem_.preAdapt();
-
- /**** 3) Put primary variables in a map *********/
- problem_.variables().storePrimVars(problem_);
-
- /**** 4) Adapt Grid and size of variable vectors *****/
- problem_.grid().adapt();
-
- // forceRefineRatio(1);
-
- // update mapper to new cell indices
- using ElementMapper = std::decay_t<decltype(problem_.variables().elementMapper())>;
- if constexpr (Deprecated::hasUpdateGridView<ElementMapper, GridView>())
- problem_.variables().elementMapper().update(problem_.gridView());
- else
- Deprecated::update(problem_.variables().elementMapper());
-
- // adapt size of vectors
- problem_.variables().adaptVariableSize(problem_.variables().elementMapper().size());
-
- /**** 5) (Re-)construct primary variables to new grid **/
- problem_.variables().reconstructPrimVars(problem_);
-
- // delete markers in grid
- problem_.grid().postAdapt();
-
- // call method in problem for potential output etc.
- problem_.postAdapt();
-
- return;
- }
-
- /*!
- * Mark Elements for grid refinement according to applied Indicator
- * @return Total ammount of marked cells
- */
- template<class Indicator>
- void markElements(Indicator& indicator)
- {
- using CoarsenMarkerType = std::unordered_map<IdType, IdType>;
- CoarsenMarkerType coarsenMarker;
- const IdSet& idSet(problem_.grid().localIdSet());
-
- for (const auto& element : elements(problem_.gridView()))
- {
- // only mark non-ghost elements
- if (element.partitionType() == Dune::GhostEntity)
- continue;
-
- // refine?
- if (indicator.refine(element) && element.level() < levelMax_)
- {
- problem_.grid().mark( 1, element);
- ++marked_;
-
- // this also refines the neighbor elements
- checkNeighborsRefine_(element);
- }
- if (indicator.coarsen(element) && element.hasFather())
- {
- IdType idx = idSet.id(element.father());
- auto it = coarsenMarker.find(idx);
- if (it != coarsenMarker.end())
- {
- ++it->second;
- }
- else
- {
- coarsenMarker[idx] = 1;
- }
- }
- }
- // coarsen
- for (const auto& element : elements(problem_.gridView()))
- {
- // only mark non-ghost elements
- if (element.partitionType() == Dune::GhostEntity)
- continue;
-
- if (indicator.coarsen(element) && element.level() > levelMin_)
- {
- IdType idx = idSet.id(element.father());
- auto it = coarsenMarker.find(idx);
- if (it != coarsenMarker.end())
- {
- if (problem_.grid().getMark(element) == 0
- && it->second == element.geometry().corners())
- {
- // check if coarsening is possible
- bool coarsenPossible = true;
- for(const auto& intersection : intersections(problem_.gridView(), element))
- {
- if(intersection.neighbor())
- {
- auto outside = intersection.outside();
- if ((problem_.grid().getMark(outside) > 0)
- || outside.level() > element.level())
- {
- coarsenPossible = false;
- }
- }
- }
-
- if(coarsenPossible)
- {
- problem_.grid().mark( -1, element );
- ++coarsened_;
- }
- }
- }
- }
- }
- }
-
- /*!
- * @brief Returns true if grid cells have been marked for adaptation
- */
- bool wasAdapted()
- {
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
-
- return (sumMarked != 0 || sumCoarsened != 0);
- }
-
- /*!
- * Sets minimum and maximum refinement levels
- *
- * @param levMin minimum level for coarsening
- * @param levMax maximum level for refinement
- */
- void setLevels(int levMin, int levMax)
- {
- if (levMin < 0)
- Dune::dgrave << __FILE__<< ":" <<__LINE__
- << " : Dune cannot coarsen to gridlevels smaller 0! "<< std::endl;
- levelMin_ = levMin;
- levelMax_ = levMax;
- }
-
- /*!
- * @brief Returns maximum refinement level
- *
- * The value is the assign maximum possible level,
- * not the actual maximum level of the grid.
- * @return levelMax_ maximum level for refinement
- */
- int getMaxLevel() const
- {
- return levelMax_;
- }
- /*!
- * @brief Returns minimum refinement level
- *
- * The value is the assign minimum possible level,
- * not the actual minimum level of the grid.
- * @return levelMin_ minimum level for coarsening
- */
- int getMinLevel() const
- {
- return levelMin_;
- }
-
- AdaptionIndicator& adaptionIndicator()
- {
- return adaptionIndicator_;
- }
-
- AdaptionIndicator& adaptionIndicator() const
- {
- return adaptionIndicator_;
- }
-
-private:
- /*!
- * @brief Method ensuring the refinement ratio of 2:1
- *
- * For any given entity, a loop over the neighbors checks weather the
- * entities refinement would require that any of the neighbors has
- * to be refined, too.
- * This is done recursively over all levels of the grid.
- *
- * @param entity Element of interest that is to be refined
- * @param level level of the refined entity: it is at least 1
- * @return true if everything was successful
- */
- bool checkNeighborsRefine_(const Element &entity, int level = 1)
- {
- // this also refines the neighbor elements
- for(const auto& intersection : intersections(problem_.gridView(), entity))
- {
- if(!intersection.neighbor())
- continue;
-
- auto outside = intersection.outside();
-
- // only mark non-ghost elements
- if (outside.partitionType() == Dune::GhostEntity)
- continue;
-
- if ((outside.level() < levelMax_)
- && (outside.level() < entity.level()))
- {
- problem_.grid().mark(1, outside);
- ++marked_;
-
- if(level != levelMax_)
- checkNeighborsRefine_(outside, ++level);
- }
- }
- return true;
- }
-
-
- /*!
- * \brief Enforces a given refine ratio after grid was adapted
- *
- * If the refine ratio is not taken into consideration during
- * marking, then this method ensures a certain ratio.
- *
- * @param maxLevelDelta The maximum level difference (refine ratio)
- * between neighbors.
- */
- void forceRefineRatio(int maxLevelDelta = 1)
- {
- LeafGridView leafGridView = problem_.gridView();
- // delete all existing marks
- problem_.grid().postAdapt();
- bool done;
- do
- {
- // run through all cells
- done=true;
- for (const auto& element : elements(problem_.gridView()))
- {
- // only mark non-ghost elements
- if (element.partitionType() == Dune::GhostEntity)
- continue;
-
- // run through all neighbor-cells (intersections)
- for (const auto& intersection : intersections(leafGridView, element))
- {
- if(!intersection.neighbor())
- continue;
-
- if (element.level() + maxLevelDelta < intersection.outside().level())
- {
- problem_.grid().mark( 1, element );
- done=false;
- }
- }
- }
- if (done==false)
- {
- // adapt the grid
- problem_.grid().adapt();
- // delete marks
- problem_.grid().postAdapt();
- }
- }
- while (done!=true);
- }
-
- // private Variables
- Problem& problem_;
- AdaptionIndicator adaptionIndicator_;
-
- int marked_;
- int coarsened_;
-
- int levelMin_; //!< minimum allowed level
- int levelMax_; //!< maximum allowed level
-
- int adaptationInterval_; //!< time step interval for adaption
-};
-
-/*!
- * @brief Class for NON-adaptive simulations
- *
- * This class provides empty methods for non-adaptive simulations
- * for compilation reasons. If adaptivity is desired, create the
- * class with template arguments <TypeTag, true> instead.
- */
-template<class TypeTag>
-class GridAdapt<TypeTag, false>
-{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
-
-public:
- void init()
- {}
- void adaptGrid()
- {}
- bool wasAdapted()
- {
- return false;
- }
- void setLevels(int, int)
- {}
- void setTolerance(int, int)
- {}
- void setIndicator(const ScalarSolutionType&,
- const Scalar&, const Scalar&)
- {}
- GridAdapt (Problem& problem)
- {}
-};
-
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/gridadaptinitializationindicator.hh b/dumux/porousmediumflow/sequential/gridadaptinitializationindicator.hh
deleted file mode 100644
index 82af602cfb298535f6bc15b0153661dad82eede2..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/gridadaptinitializationindicator.hh
+++ /dev/null
@@ -1,401 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
-#define DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
-
-#include "properties.hh"
-
-#include <dune/common/dynvector.hh>
-/**
- * @file
- * @brief Class defining a start indicator for grid adaption
- */
-namespace Dumux
-{
-/*!\ingroup IMPES
- * @brief Class defining a start indicator for grid adaption
- *
- * Uses the defined grid adaptation indicator and further accounts for sources and boundaries.
- * Only for grid initialization!
- *
- * \tparam TypeTag The problem TypeTag
- */
-template<class TypeTag>
-class GridAdaptInitializationIndicator
-{
-private:
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using AdaptionIndicator = GetPropType<TypeTag, Properties::AdaptionIndicator>;
-
- using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
-
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- numEq = getPropValue<TypeTag, Properties::NumEq>(),
- numPhases = getPropValue<TypeTag, Properties::NumPhases>()
- };
-
- enum
- {
- refineCell = 1,
- coarsenCell = -1
- };
-
- using LocalPosition = Dune::FieldVector<Scalar, dim>;
- using LocalPositionFace = Dune::FieldVector<Scalar, dim-1>;
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-
- void virtualHierarchicSourceSearch_(PrimaryVariables &source, const Element& element)
- {
- int level = element.level();
-
- if (level == maxAllowedLevel_)
- {
- GlobalPosition globalPos = element.geometry().center();
- problem_.sourceAtPos(source, globalPos);
-
- return;
- }
-
- unsigned int numRefine = maxAllowedLevel_ - level;
- int numCheckCoords = pow(2, numRefine);
-
- LocalPosition localPos(0.0);
- GlobalPosition globalPosCheck(0.0);
- Scalar halfInterval = (1.0/double(numCheckCoords))/2.;
-
- PrimaryVariables sourceCheck(0.0);
-
- using std::abs;
- for (int i = 1; i <= numCheckCoords; i++)
- {
- for (int j = 1; j <= numCheckCoords; j++)
- {
- localPos[0] = double(i)/double(numCheckCoords) - halfInterval;
- localPos[1] = double(j)/double(numCheckCoords) - halfInterval;
- if (dim == 2)
- {
- globalPosCheck = element.geometry().global(localPos);
- problem_.sourceAtPos(sourceCheck, globalPosCheck);
-
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- if (abs(sourceCheck[eqIdx]) > abs(source[eqIdx]))
- {
- source[eqIdx] = sourceCheck[eqIdx];
- }
- }
- }
- else if (dim == 3)
- {
- for (int k = 1; k <= numCheckCoords; k++)
- {
- localPos[2] = double(k)/double(numCheckCoords) - halfInterval;
- globalPosCheck = element.geometry().global(localPos);
- problem_.sourceAtPos(sourceCheck, globalPosCheck);
-
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- if (abs(sourceCheck[eqIdx]) > abs(source[eqIdx]))
- {
- source[eqIdx] = sourceCheck[eqIdx];
- }
- }
- }
- }
- }
- }
- }
-
- void virtualHierarchicBCSearch_(BoundaryTypes &bcTypes, PrimaryVariables &values, const Element& element, const Intersection& intersection)
- {
- int level = element.level();
-
- if (level == maxAllowedLevel_)
- {
- GlobalPosition globalPos = intersection.geometry().center();
- problem_.boundaryTypesAtPos(bcTypes, globalPos);
-
- if (refineAtFluxBC_)
- {
- for (int i = 0; i < numEq; i++)
- {
- if (bcTypes.isNeumann(i))
- {
- PrimaryVariables fluxes;
- problem_.neumannAtPos(fluxes, globalPos);
-
- values += fluxes;
- }
- }
- }
- return;
- }
-
- unsigned int numRefine = maxAllowedLevel_ - level;
- int numCheckCoords = pow(2, numRefine);
-
- LocalPositionFace localPos(0.0);
- GlobalPosition globalPosCheck(0.0);
- Scalar halfInterval = (1.0/double(numCheckCoords))/2.;
-
- PrimaryVariables fluxCheck(0.0);
-
- for (int i = 1; i <= numCheckCoords; i++)
- {
- localPos[0] = double(i)/double(numCheckCoords) - halfInterval;
- if (dim == 2)
- {
- globalPosCheck = intersection.geometry().global(localPos);
- problem_.boundaryTypesAtPos(bcTypes, globalPosCheck);
-
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- if (refineAtDirichletBC_ && bcTypes.isDirichlet(eqIdx))
- {
- return;
- }
- else if (refineAtFluxBC_ && bcTypes.isNeumann(eqIdx))
- {
- problem_.neumannAtPos(fluxCheck, globalPosCheck);
-
- values += fluxCheck;
- }
- }
- }
- else if (dim == 3)
- {
- for (int k = 1; k <= numCheckCoords; k++)
- {
- localPos[1] = double(k)/double(numCheckCoords) - halfInterval;
- globalPosCheck = intersection.geometry().global(localPos);
-
- problem_.boundaryTypesAtPos(bcTypes, globalPosCheck);
-
-
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- if (refineAtDirichletBC_ && bcTypes.isDirichlet(eqIdx))
- {
- return;
- }
- else if (refineAtFluxBC_ && bcTypes.isNeumann(eqIdx))
- {
- problem_.neumannAtPos(fluxCheck, globalPosCheck);
-
- values += fluxCheck;
- }
-
- }
- }
- }
- }
- }
-
-
-public:
- /*! \brief Calculates the indicator used for refinement/coarsening for each grid cell.
- *
- */
- void calculateIndicator()
- {
- //First Adapt for boundary conditions and sources to get a correct pressure solution
- if (nextMaxLevel_ == maxAllowedLevel_)
- adaptionIndicator_.calculateIndicator();
-
- // prepare an indicator for refinement
- indicatorVector_.resize(problem_.variables().cellDataGlobal().size());
-
- indicatorVector_ = coarsenCell;
-
- if (!enableInitializationIndicator_)
- return;
-
- // 1) calculate Indicator -> min, maxvalues
- // Schleife über alle Leaf-Elemente
- using std::abs;
- using std::max;
- using std::min;
- for (const auto& element : elements(problem_.gridView()))
- {
- int globalIdxI = problem_.variables().index(element);
-
- int level = element.level();
- maxLevel_ = max(level, maxLevel_);
-
- if (level < minAllowedLevel_)
- {
- nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
- indicatorVector_[globalIdxI] = refineCell;
- continue;
- }
-
- if (refineAtSource_)
- {
- PrimaryVariables source(0.0);
- virtualHierarchicSourceSearch_(source, element);
- for (int i = 0; i < numEq; i++)
- {
- if (abs(source[i]) > 1e-10)
- {
- nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
- indicatorVector_[globalIdxI] = refineCell;
- break;
- }
- }
- }
-
- if (indicatorVector_[globalIdxI] != refineCell && (refineAtDirichletBC_ || refineAtFluxBC_))
- {
- // Berechne Verfeinerungsindikator an allen Zellen
- for (const auto& intersection : intersections(problem_.gridView(), element))
- {
- if (intersection.boundary() && indicatorVector_[globalIdxI] != refineCell)
- {
- BoundaryTypes bcTypes;
- PrimaryVariables values(0.0);
-
- virtualHierarchicBCSearch_(bcTypes, values, element, intersection);
-
-
- for (int i = 0; i < numEq; i++)
- {
- if (bcTypes.isDirichlet(i) && refineAtDirichletBC_)
- {
- nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
- indicatorVector_[globalIdxI] = refineCell;
- break;
- }
- }
- for (int j = 0; j < numPhases; j++)
- {
- if (abs(values[j]) > 1e-10)
- {
- nextMaxLevel_ = min(max(level + 1, nextMaxLevel_), maxAllowedLevel_);
- indicatorVector_[globalIdxI] = refineCell;
- break;
- }
-
- }
- }
- }
- }
- }
- }
-
- /*! \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);
-
- if (indicatorVector_[idx] == refineCell)
- return true;
- else if (maxLevel_ == maxAllowedLevel_)
- return adaptionIndicator_.refine(element);
- else
- return false;
- }
-
- /*! \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);
-
- if (indicatorVector_[idx] == coarsenCell && maxLevel_ < maxAllowedLevel_)
- return true;
- else if (indicatorVector_[idx] == coarsenCell && !adaptionIndicator_.refine(element))
- return true;
- else
- return false;
- }
-
- int maxLevel()
- {
- return maxLevel_;
- }
-
- /*! \brief Initializes the adaption indicator class */
- void init()
- {}
-
- bool initializeModel()
- {
- return nextMaxLevel_ == maxAllowedLevel_;
- }
-
- /*! \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
- * \param adaptionIndicator Indicator whether a be adapted
- */
- GridAdaptInitializationIndicator(Problem& problem, AdaptionIndicator& adaptionIndicator):
- problem_(problem), adaptionIndicator_(adaptionIndicator), maxLevel_(0), nextMaxLevel_(0)
- {
- minAllowedLevel_ = getParam<int>("GridAdapt.MinLevel");
- maxAllowedLevel_ = getParam<int>("GridAdapt.MaxLevel");
- enableInitializationIndicator_ = getParam<bool>("GridAdapt.EnableInitializationIndicator");
- refineAtDirichletBC_ = getParam<bool>("GridAdapt.RefineAtDirichletBC");
- refineAtFluxBC_ = getParam<bool>("GridAdapt.RefineAtFluxBC");
- refineAtSource_ = getParam<bool>("GridAdapt.RefineAtSource");
-
- if (!refineAtDirichletBC_ && !refineAtFluxBC_ && !refineAtSource_)
- {
- nextMaxLevel_ = maxAllowedLevel_;
- maxLevel_ = maxAllowedLevel_;
- }
- }
-
-private:
- Problem& problem_;
- AdaptionIndicator& adaptionIndicator_;
- Dune::DynamicVector<int> indicatorVector_;
- int maxLevel_;
- int nextMaxLevel_;
- int minAllowedLevel_; //!< minimum allowed level
- int maxAllowedLevel_; //!< maximum allowed level
- bool enableInitializationIndicator_;
- bool refineAtDirichletBC_;
- bool refineAtFluxBC_;
- bool refineAtSource_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/gridadaptinitializationindicatordefault.hh b/dumux/porousmediumflow/sequential/gridadaptinitializationindicatordefault.hh
deleted file mode 100644
index 263077fabe98160c8753189bfa4e126a0f5d3dd3..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/gridadaptinitializationindicatordefault.hh
+++ /dev/null
@@ -1,111 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_GRIDADAPTINITIALIZATIONINDICATORDEFAULT_HH
-#define DUMUX_GRIDADAPTINITIALIZATIONINDICATORDEFAULT_HH
-
-#include "properties.hh"
-#include <dumux/common/properties.hh>
-
-#include <dune/common/dynvector.hh>
-
-/**
- * @file
- * @brief Class defining a start indicator for grid adaption
- */
-namespace Dumux
-{
-/*!\ingroup IMPES
- * @brief Class defining a start indicator for grid adaption
- *
- *Default implementation
- *
- * \tparam TypeTag The problem TypeTag
- */
-template<class TypeTag>
-class GridAdaptInitializationIndicatorDefault
-{
-private:
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using AdaptionIndicator = GetPropType<TypeTag, Properties::AdaptionIndicator>;
-
-public:
- /*! \brief Calculates the indicator used for refinement/coarsening for each grid cell.
- *
- */
- void calculateIndicator()
- {}
-
- /*! \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)
- {
- return false;
- }
-
- /*! \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)
- {
- return false;
- }
-
- int maxLevel()
- {
- return maxLevel_;
- }
-
-
- /*! \brief Initializes the adaption indicator class*/
- void init()
- {}
-
- /*! \brief Returns true if the IMPET-Model needs to be initialized*/
- bool initializeModel()
- {
- return false;
- }
-
- /*! \brief Constructs a GridAdaptionIndicator for initialization of an adaptive grid
- *
- * Default implementation
- *
- * \param problem The problem object
- * \param adaptionIndicator Indicator whether a be adapted
- */
- GridAdaptInitializationIndicatorDefault(Problem& problem, AdaptionIndicator& adaptionIndicator)
- {
- maxLevel_ = getParam<int>("GridAdapt.MaxLevel");
- }
-
-private:
- int maxLevel_; //!< maximum allowed level
-};
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/gridadaptproperties.hh b/dumux/porousmediumflow/sequential/gridadaptproperties.hh
deleted file mode 100644
index e2d98fc6b04061b7c8b86d5c9806f70535bb16ab..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/gridadaptproperties.hh
+++ /dev/null
@@ -1,94 +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 3 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/>. *
- *****************************************************************************/
-/*!
- * \ingroup IMPETProperties
- * \ingroup IMPET
- * \file
- *
- * \brief Defines a type tag and some fundamental properties for
- * linear solvers
- */
-#ifndef DUMUX_GRIDADAPT_PROPERTIES_HH
-#define DUMUX_GRIDADAPT_PROPERTIES_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-
-namespace Dumux
-{
-// forward declarations
-template<class TypeTag>
-class GridAdaptInitializationIndicatorDefault;
-
-template<class TypeTag, bool>
-class GridAdapt;
-
-
-namespace Properties
-{
-//! Grid adaption type tag for all sequential models.
-namespace TTag {
-struct GridAdapt {};
-}
-
-//! Defines if the grid is h-adaptive
-template<class TypeTag, class MyTypeTag>
-struct AdaptiveGrid { using type = UndefinedProperty; };
-
-//! The type of grid adaptation
-template<class TypeTag, class MyTypeTag>
-struct GridAdaptModel { using type = UndefinedProperty; };
-
-//! Class defining the refinement/coarsening indicator
-template<class TypeTag, class MyTypeTag>
-struct AdaptionIndicator { using type = UndefinedProperty; };
-
-//! Class defining the refinement/coarsening indicator for grid initialization
-template<class TypeTag, class MyTypeTag>
-struct AdaptionInitializationIndicator { using type = UndefinedProperty; };
-
-//! Tolerance for refinement
-template<class TypeTag, class MyTypeTag>
-struct GridAdaptRefineThreshold { using type = UndefinedProperty; };
-
-//! Tolerance for coarsening
-template<class TypeTag, class MyTypeTag>
-struct GridAdaptCoarsenThreshold { using type = UndefinedProperty; };
-
-//no adaptive grid
-template<class TypeTag>
-struct AdaptiveGrid<TypeTag, TTag::GridAdapt> { static constexpr bool value = false; };
-
-//Set default class for adaptation initialization indicator
-template<class TypeTag>
-struct AdaptionInitializationIndicator<TypeTag, TTag::GridAdapt> { using type = GridAdaptInitializationIndicatorDefault<TypeTag>; };
-//Set default class for adaptation
-template<class TypeTag>
-struct GridAdaptModel<TypeTag, TTag::GridAdapt> { using type = GridAdapt<TypeTag, getPropValue<TypeTag, Properties::AdaptiveGrid>()>; };
-
-//standard setting
-template<class TypeTag>
-struct GridAdaptRefineThreshold<TypeTag, TTag::GridAdapt> { static constexpr auto value = 0.0; };
-template<class TypeTag>
-struct GridAdaptCoarsenThreshold<TypeTag, TTag::GridAdapt> { static constexpr auto value = 0.0; };
-} // namespace Properties
-} // namespace Dumux
-
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/impet.hh b/dumux/porousmediumflow/sequential/impet.hh
deleted file mode 100644
index 78eb12d5931b40b993d1f18331b9c50cbc840838..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/impet.hh
+++ /dev/null
@@ -1,228 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_IMPET_HH
-#define DUMUX_IMPET_HH
-
-/**
- * @file
- * @brief IMPET scheme
- */
-
-#include "impetproperties.hh"
-
-namespace Dumux
-{
-/**
- * \ingroup IMPET
- * \brief IMplicit Pressure Explicit Transport (IMPET) scheme for the solution of weakly coupled diffusion-transport formulations.
- *
- * The model implements the sequential equations of two-phase flow.
- * These equations can be derived from the two-phase flow equations shown for the two-phase box model (TwoPBoxModel).
- * The first equation to solve is a pressure equation of elliptic character.
- * The second one is a transport equation (e.g. for saturation, concentration,...), which can be hyperbolic or parabolic.
- *
- * As the equations are only weakly coupled they do not have to be solved simultaneously
- * but can be solved sequentially. First the pressure equation is solved implicitly,
- * second the transport equation can be solved explicitly. This solution procedure is called IMPES algorithm
- * (IMplicit Pressure Explicit Saturation) for immiscible flow or IMPEC algorithm
- * (IMplicit Pressure Explicit Concentration) for miscible flow.
- *
- * \tparam TypeTag The Type Tag
- */
-template<class TypeTag> class IMPET
-{
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using ElementMapper = typename SolutionTypes::ElementMapper;
- using TransportSolutionType = GetPropType<TypeTag, Properties::TransportSolutionType>;
-
- enum IterationType
- {
- noIter,
- iterToNumIter,
- iterToConverged,
- };
-
-public:
- using SolutionType = typename SolutionTypes::ScalarSolution;
-
- //! Set initial solution and initialize parameters
- void initialize()
- {
- //initial values of transported quantity
- problem_.transportModel().initialize();
- //call function with true to get a first initialisation of the pressure field
- problem_.pressureModel().initialize();
-
- //update constitutive functions
- problem_.pressureModel().updateMaterialLaws();
-
- Dune::dinfo << "IMPET: initialization done." << std::endl;
-
- return;
- }
-
- //! Calculate the update.
- /*!
- * \param t time
- * \param dt time step size
- * \param updateVec vector for the update values
- *
- * Calculates the new pressure and velocity and determines the time step size
- * and the update of the transported quantity for the explicit time step.
- */
- void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
- {
- if (iterFlag_ == noIter)
- {
- //update pressure
- problem_.pressureModel().update();
-
- //calculate defect of transported quantity
- problem_.transportModel().update(t, dt, updateVec, true);
-
- dt *= cFLFactor_;
- }
- else if (iterFlag_ == iterToNumIter || iterFlag_ == iterToConverged)
- {
- bool converg = false;
- int iter = 0;
- int iterTot = 0;
-
- // the method is valid for any transported quantity.
- TransportSolutionType transValueOldIter;
- problem_.transportModel().getTransportedQuantity(transValueOldIter);
- TransportSolutionType updateOldIter(transValueOldIter);
- updateOldIter = 0;
- TransportSolutionType transportedQuantity(transValueOldIter);
- TransportSolutionType updateHelp(transValueOldIter);
- TransportSolutionType updateDiff(transValueOldIter);
-
- while (!converg)
- {
- iter++;
- iterTot++;
-
- problem_.pressureModel().update();
-
- //calculate defect of transported quantity
- problem_.transportModel().update(t, dt, updateVec, true);
-
- updateHelp = updateVec;
- problem_.transportModel().getTransportedQuantity(transportedQuantity);
- transportedQuantity += (updateHelp *= (dt * cFLFactor_));
- transportedQuantity *= omega_;
- transValueOldIter *= (1 - omega_);
- transportedQuantity += transValueOldIter;
- updateDiff = updateVec;
- updateDiff -= updateOldIter;
- transValueOldIter = transportedQuantity;
- updateOldIter = updateVec;
- Dune::dinfo << " defect = " << dt * updateDiff.two_norm() / transportedQuantity.two_norm();
- // break criteria for iteration loop
- if (iterFlag_ == iterToConverged && dt * updateDiff.two_norm() / transportedQuantity.two_norm() <= maxDefect_)
- {
- converg = true;
- }
- else if (iterFlag_ == iterToNumIter && iter > nIter_)
- {
- converg = true;
- }
-
- if (iterFlag_ == iterToConverged && transportedQuantity.infinity_norm() > (1 + maxDefect_))
- {
- converg = false;
- }
- if (!converg && iter > nIter_)
- {
- converg = true;
- std::cout << "Nonlinear loop in IMPET.update exceeded nIter = " << nIter_ << " iterations."
- << std::endl;
- std::cout << transportedQuantity.infinity_norm() << std::endl;
- }
- }
- // outputs
- if (iterFlag_ == iterToConverged)
- std::cout << "Iteration steps: " << iterTot << std::endl;
- std::cout.setf(std::ios::scientific, std::ios::floatfield);
-
- dt *= cFLFactor_;
- }
- else
- {
- DUNE_THROW(Dune::NotImplemented,"IMPET: Iterationtype not implemented!");
- }
-
- }
-
- void updateTransport(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
- {
- problem_.pressureModel().updateVelocity();
-
- problem_.transportModel().update(t, dt, updateVec, true);
- }
-
- //! \brief Write data files
- /*!
- * calls the output methods of both models, pressure and transport.
- * \param writer the current VTKwriter
- */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {
- problem_.pressureModel().addOutputVtkFields(writer);
- problem_.transportModel().addOutputVtkFields(writer);
-
- return;
- }
-
- /*!
- * \brief Mapper for the entities where degrees of freedoms are defined.
- */
- const ElementMapper &dofMapper() const
- {
- return problem_.elementMapper();
- }
-
- //! Constructs an IMPET object
- /**
- * \param problem Problem
- */
- IMPET(Problem& problem) :
- problem_(problem)
- {
- cFLFactor_ = getParam<Scalar>("Impet.CFLFactor");
- iterFlag_ = getParam<int>("Impet.IterationFlag", 0);
- nIter_ = getParam<int>("Impet.IterationNumber", 2);
- maxDefect_ = getParam<Scalar>("Impet.MaximumDefect", 1e-5);
- omega_ = getParam<Scalar>("Impet.RelaxationFactor", 1.0);
- }
-
-private:
- Problem& problem_; //!< object of type Problem including problem
- Scalar cFLFactor_;
- int iterFlag_; //!< flag to switch the iteration type of the IMPET scheme
- int nIter_; //!< number of iterations if IMPET iterations are enabled by IterationFlag
- Scalar maxDefect_; //!< maximum defect if IMPET iterations are enabled by IterationFlag
- Scalar omega_; //!< 1 = new solution is new solution, 0 = old solution is new solution
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/impetproblem.hh b/dumux/porousmediumflow/sequential/impetproblem.hh
deleted file mode 100644
index fc0f7276bafe985649eefe9eb15518d456e54b06..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/impetproblem.hh
+++ /dev/null
@@ -1,869 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_IMPETPROBLEM_HH
-#define DUMUX_IMPETPROBLEM_HH
-
-#include <dune/common/float_cmp.hh>
-#include "impetproperties.hh"
-#include <dumux/io/vtkmultiwriter.hh>
-#include <dumux/io/restart.hh>
-#include <dumux/io/adaptivegridrestart.hh>
-
-#include <dumux/porousmediumflow/sequential/gridadapt.hh>
-
-/**
- * @file
- * @brief Base class for defining an instance of the diffusion problem
- */
-
-namespace Dumux
-{
-/*!
- * \ingroup IMPET
- * @brief base class for problems using a sequential implicit-explicit strategy
- *
- * \tparam TypeTag problem TypeTag
- * \tparam Implementation problem implementation
- */
-template<class TypeTag>
-class IMPETProblem
-{
-private:
- using Implementation = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Grid = GetPropType<TypeTag, Properties::Grid>;
-
- using TimeManager = GetPropType<TypeTag, Properties::TimeManager>;
-
- using VtkMultiWriter = Dumux::VtkMultiWriter<GridView>;
-
- using Variables = GetPropType<TypeTag, Properties::Variables>;
-
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using VertexMapper = typename SolutionTypes::VertexMapper;
- using ElementMapper = typename SolutionTypes::ElementMapper;
-
- using IMPETModel = GetPropType<TypeTag, Properties::Model>;
- using TransportSolutionType = GetPropType<TypeTag, Properties::TransportSolutionType>;
- using PressureModel = GetPropType<TypeTag, Properties::PressureModel>;
- using TransportModel = GetPropType<TypeTag, Properties::TransportModel>;
-
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
- enum
- {
- wetting = 0, nonwetting = 1,
- adaptiveGrid = getPropValue<TypeTag, Properties::AdaptiveGrid>()
- };
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
- // The module to adapt grid. If adaptiveGrid is false, this model does nothing.
- using GridAdaptModel = GetPropType<TypeTag, Properties::GridAdaptModel>;
-
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
-
- //private!! copy constructor
- IMPETProblem(const IMPETProblem&)
- {}
-
-public:
- /*! \brief Constructs an object of type IMPETProblemProblem
- * \param timeManager The time manager
- * \param grid The grid
- */
- IMPETProblem(TimeManager &timeManager, Grid& grid)
- : IMPETProblem(timeManager, grid, grid.leafGridView())
- {}
-
- IMPETProblem(TimeManager &timeManager, Grid& grid, const GridView& gridView)
- : grid_(&grid),
- gridView_(gridView),
- bBoxMin_(std::numeric_limits<double>::max()),
- bBoxMax_(-std::numeric_limits<double>::max()),
- timeManager_(&timeManager),
- outputInterval_(1),
- outputTimeInterval_(0.0),
- vtkOutputLevel_(-1)
- {
- // calculate the bounding box of the grid view
- using std::min;
- using std::max;
- for (const auto& vertex : vertices(this->gridView())) {
- for (int i=0; i<dim; i++) {
- bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
- bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
- }
- }
-
- // communicate to get the bounding box of the whole domain
- if (this->gridView().comm().size() > 1)
- for (int i = 0; i < dim; ++i) {
- bBoxMin_[i] = this->gridView().comm().min(bBoxMin_[i]);
- bBoxMax_[i] = this->gridView().comm().max(bBoxMax_[i]);
- }
-
- variables_ = std::make_shared<Variables>(this->gridView());
- pressModel_ = std::make_shared<PressureModel>(asImp_());
-
- transportModel_ = std::make_shared<TransportModel>(asImp_());
- model_ = std::make_shared<IMPETModel>(asImp_()) ;
-
- // create an Object to handle adaptive grids
- if (adaptiveGrid)
- gridAdapt_ = std::make_shared<GridAdaptModel>(asImp_());
-
- vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
- dtVariationRestrictionFactor_ = getParam<Scalar>("Impet.DtVariationRestrictionFactor", std::numeric_limits<Scalar>::max());
- maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
- }
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param bcTypes The boundary types for the conservation equations
- * \param intersection The intersection for which the boundary type is set
- */
- void boundaryTypes(BoundaryTypes &bcTypes,
- const Intersection &intersection) const
- {
- // forward it to the method which only takes the global coordinate
- asImp_().boundaryTypesAtPos(bcTypes, intersection.geometry().center());
- }
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param bcTypes The boundary types for the conservation equations
- * \param globalPos The position of the center of the boundary intersection
- */
- void boundaryTypesAtPos(BoundaryTypes &bcTypes,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Dirichlet conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem does not provide "
- "a boundaryTypesAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param intersection The boundary intersection
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichlet(PrimaryVariables &values,
- const Intersection &intersection) const
- {
- // forward it to the method which only takes the global coordinate
- asImp_().dirichletAtPos(values, intersection.geometry().center());
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position of the center of the boundary intersection
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Dirichlet conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem specifies that some boundary "
- "segments are dirichlet, but does not provide "
- "a dirichletAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
- * \param intersection The boundary intersection
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumann(PrimaryVariables &values,
- const Intersection &intersection) const
- {
- // forward it to the interface with only the global position
- asImp_().neumannAtPos(values, intersection.geometry().center());
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
- * \param globalPos The position of the center of the boundary intersection
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumannAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Neumann conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem specifies that some boundary "
- "segments are neumann, but does not provide "
- "a neumannAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the source term
- *
- * \param values The source and sink values for the conservation equations
- * \param element The element
- *
- * For this method, the \a values parameter stores the rate mass
- * generated or annihilate per volume unit. Positive values mean
- * that mass is created, negative ones mean that it vanishes.
- */
- void source(PrimaryVariables &values,
- const Element &element) const
- {
- // forward to generic interface
- asImp_().sourceAtPos(values, element.geometry().center());
- }
-
- /*!
- * \brief Evaluate the source term for all phases within a given
- * sub-control-volume.
- *
- * \param values The source and sink values for the conservation equations
- * \param globalPos The position of the center of the finite volume
- * for which the source term ought to be
- * specified in global coordinates
- *
- * For this method, the \a values parameter stores the rate mass
- * generated or annihilate per volume unit. Positive values mean
- * that mass is created, negative ones mean that it vanishes.
- */
- void sourceAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- { // Throw an exception (there is no initial condition)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem does not provide "
- "a sourceAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The initial values for the primary variables
- * \param element The element
- *
- * For this method, the \a values parameter stores primary
- * variables.
- */
- void initial(PrimaryVariables &values,
- const Element &element) const
- {
- // forward to generic interface
- asImp_().initialAtPos(values, element.geometry().center());
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position of the center of the finite volume
- * for which the initial values ought to be
- * set (in global coordinates)
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // initialize with 0 by default
- values = 0;
- }
-
- /*!
- * \brief Called by the TimeManager in order to
- * initialize the problem.
- */
- void init()
- {
- // set the initial condition of the model
- variables_->initialize();
- model().initialize();
- if (adaptiveGrid)
- gridAdapt().init();
- }
-
- /*!
- * \brief Called by TimeManager just before the time
- * integration.
- */
- void preTimeStep()
- {
- // if adaptivity is used, this method adapts the grid.
- // if it is not used, this method does nothing.
- if (adaptiveGrid && timeManager().timeStepIndex() > 0)
- this->gridAdapt().adaptGrid();
- asImp_().pressureModel().updateMaterialLaws();
- }
-
- /*!
- * \brief Called by TimeManager in order to do a time
- * integration on the model.
- *
- * \note \a timeStepSize and \a nextStepSize are references and may
- * be modified by the timeIntegration(). On exit of this
- * function \a timeStepSize must contain the step size
- * actually used by the time integration for the current
- * steo, and \a nextStepSize must contain a suggestion for the
- * next time step size.
- */
- void timeIntegration()
- {
- // allocate temporary vectors for the updates
- TransportSolutionType updateVector;//(this->transportModel().solutionSize());
-
- Scalar t = timeManager().time();
- Scalar dt = std::numeric_limits<Scalar>::max();
-
- // obtain the first update and the time step size
- model().update(t, dt, updateVector);
-
- using std::max;
- using std::min;
- if (t > 0 && timeManager().timeStepSize()> 0.)
- {
- Scalar oldDt = timeManager().timeStepSize();
- Scalar minDt = max(oldDt - dtVariationRestrictionFactor_ * oldDt, 0.0);
- Scalar maxDt = oldDt + dtVariationRestrictionFactor_ * oldDt;
- dt = max(min(dt, maxDt), minDt);
- }
-
- //make sure t_old + dt is not larger than tend
- dt = min(dt, timeManager().episodeMaxTimeStepSize());
-
- // check if we are in first TS and an initialDt was assigned
- if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(t, 0.0, 1.0e-30)
- && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(timeManager().timeStepSize(), 0.0, 1.0e-30))
- {
- if (gridView().comm().size() > 1)
- dt = gridView().comm().min(dt);
-
- // check if assigned initialDt is in accordance with dt from first transport step
- if (timeManager().timeStepSize() > dt
- && gridView().comm().rank() == 0)
- Dune::dwarn << "Initial timestep of size " << timeManager().timeStepSize()
- << " is larger then dt=" << dt <<" from transport" << std::endl;
- // internally assign next timestep size
- dt = min(dt, timeManager().timeStepSize());
- }
-
- //make sure the right time-step is used by all processes in the parallel case
- if (gridView().comm().size() > 1)
- dt = gridView().comm().min(dt);
-
- // check maximum allowed time step size
- timeManager().setTimeStepSize(dt);
-
- // explicit Euler: Sat <- Sat + dt*N(Sat)
- transportModel().updateTransportedQuantity(updateVector);
- }
-
- /*!
- * \brief Called by TimeManager whenever a solution for a
- * timestep has been computed and the simulation time has
- * been updated.
- *
- * This is used to do some janitorial tasks like writing the
- * current solution to disk.
- */
- void postTimeStep()
- {}
-
- /*!
- * \brief Called by the time manager after everything which can be
- * done about the current time step is finished and the
- * model should be prepared to do the next time integration.
- */
- void advanceTimeLevel()
- {}
-
- /*!
- * \brief Returns the current time step size [seconds].
- */
- Scalar timeStepSize() const
- { return timeManager().timeStepSize(); }
-
- /*!
- * \brief Sets the current time step size [seconds].
- */
- void setTimeStepSize(const Scalar dt)
- { timeManager().setTimeStepSize(dt); }
-
- /*!
- * \brief Called by TimeManager whenever a solution for a
- * timestep has been computed and the simulation time has
- * been updated.
- */
- Scalar nextTimeStepSize(const Scalar dt)
- { return timeManager().timeStepSize();}
-
- /*!
- * \brief Returns the maximum allowed time step size [s]
- *
- * By default this the time step size is unrestricted.
- */
- Scalar maxTimeStepSize() const
- { return maxTimeStepSize_; }
-
- /*!
- * \brief Returns true if a restart file should be written to
- * disk.
- *
- * The default behaviour is to write one restart file every 5 time
- * steps. This file is intented to be overwritten by the
- * implementation.
- */
- bool shouldWriteRestartFile() const
- {
- if (outputInterval_ > 0)
- {
- return
- timeManager().timeStepIndex() > 0 &&
- (timeManager().timeStepIndex() % int(100*outputInterval_) == 0);
- }
- else
- return
- shouldWriteOutput();
- }
-
- /*!
- * \brief Sets a time interval for Output
- *
- * The default is 0.0 -> Output determined by output number interval (<tt>setOutputInterval(int)</tt>)
- */
- void setOutputTimeInterval(const Scalar timeInterval)
- {
- outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
- timeManager().startNextEpisode(outputTimeInterval_);
- }
-
- /*!
- * \brief Sets the interval for Output
- *
- * The default is 1 -> Output every time step
- */
- void setOutputInterval(const int interval)
- {
- using std::max;
- outputInterval_ = max(interval, 0);
- }
-
- /*!
- * \brief Returns true if the current solution should be written to
- * disk (i.e. as a VTK file)
- *
- * The default behaviour is to write out every the solution for
- * very time step. This file is intented to be overwritten by the
- * implementation.
- */
- bool shouldWriteOutput() const
- {
- if (outputInterval_ > 0)
- {
- if (timeManager().timeStepIndex() % outputInterval_ == 0
- || timeManager().willBeFinished()
- || timeManager().episodeWillBeFinished())
- {
- return true;
- }
- }
- else if (timeManager().willBeFinished()
- || timeManager().episodeWillBeFinished() || timeManager().timeStepIndex() == 0)
- {
- return true;
- }
- return false;
- }
-
- /*!
- * \brief Called when the end of an simulation episode is reached.
- */
- void episodeEnd()
- {
- if (outputTimeInterval_ > 0.0 && !timeManager().finished())
- {
- timeManager().startNextEpisode(outputTimeInterval_);
- }
- else if (!timeManager().finished())
- {
- std::cerr << "The end of an episode is reached, but the problem "
- << "does not override the episodeEnd() method. "
- << "Doing nothing!\n";
- }
- }
-
- // \}
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- * It could be either overwritten by the problem files, or simply
- * declared over the setName() function in the application file.
- */
- const std::string& name() const
- {
- return simname_;
- }
-
- /*!
- * \brief Set the problem name.
- *
- * This function sets the simulation name, which should be called before
- * the application porblem is declared! If not, the default name "sim"
- * will be used.
- *
- * \param newName The problem's name
- */
- void setName(std::string newName)
- {
- simname_ = newName;
- }
-
- /*!
- * \brief The GridView which used by the problem.
- */
- const GridView& gridView() const
- {
- return gridView_;
- }
-
- /*!
- * \brief Returns the current grid which used by the problem.
- */
- Grid &grid()
- {
- if (!grid_)
- {
- DUNE_THROW(Dune::InvalidStateException, "Grid was called in problemclass, "
- << "although it is not specified. Do so by using setGrid() method!");
- }
- return *grid_;
- }
- /*!
- * \brief Specifies the grid from outside the problem.
- * \param grid The grid used by the problem.
- */
- void setGrid(Grid &grid)
- {
- grid_ = &grid;
- }
-
- /*!
- * \brief Returns adaptivity model used for the problem.
- */
- GridAdaptModel& gridAdapt()
- {
- return *gridAdapt_;
- }
-
- /*!
- * \brief Capability to introduce problem-specific routines at the
- * beginning of the grid adaptation
- *
- * Function is called at the beginning of the standard grid
- * modification routine, GridAdapt::adaptGrid() .
- */
- void preAdapt()
- {
- }
-
- /*!
- * \brief Capability to introduce problem-specific routines after grid adaptation
- *
- * Function is called at the end of the standard grid
- * modification routine, GridAdapt::adaptGrid() , to allow
- * for problem-specific output etc.
- */
- void postAdapt()
- {
- // write out new grid
- // Dune::VTKWriter<LeafGridView> vtkwriter(problem_.gridView());
- // vtkwriter.write("latestgrid",Dune::VTK::appendedraw);
- }
-
- /*!
- * \brief Returns the mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return variables_->vertexMapper(); }
-
- /*!
- * \brief Returns the mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return variables_->elementMapper(); }
-
- /*!
- * \brief The coordinate of the corner of the GridView's bounding
- * box with the smallest values.
- */
- const GlobalPosition &bBoxMin() const
- { return bBoxMin_; }
-
- /*!
- * \brief The coordinate of the corner of the GridView's bounding
- * box with the largest values.
- */
- const GlobalPosition &bBoxMax() const
- { return bBoxMax_; }
-
- //! \name Access functions
- //@{
- /*!
- * \brief Returns TimeManager object used by the simulation
- */
- TimeManager &timeManager()
- { return *timeManager_; }
-
- //! \copydoc IMPETProblem::timeManager()
- const TimeManager &timeManager() const
- { return *timeManager_; }
-
- /*!
- * \brief Returns variables container
- *
- * This provides access to the important variables that are used in the
- * simulation process, such as pressure, saturation etc.
- */
- Variables& variables ()
- { return *variables_; }
-
- //! \copydoc IMPETProblem::variables ()
- const Variables& variables () const
- { return *variables_; }
-
- /*!
- * \brief Returns numerical model used for the problem.
- */
- IMPETModel &model()
- { return *model_; }
-
- //! \copydoc IMPETProblem::model()
- const IMPETModel &model() const
- { return *model_; }
-
- /*!
- * \brief Returns the pressure model used for the problem.
- */
- PressureModel &pressureModel()
- { return *pressModel_; }
-
- //! \copydoc IMPETProblem::pressureModel()
- const PressureModel &pressureModel() const
- { return *pressModel_; }
-
- /*!
- * \brief Returns transport model used for the problem.
- */
- TransportModel &transportModel()
- { return *transportModel_; }
-
- //! \copydoc IMPETProblem::transportModel()
- const TransportModel &transportModel() const
- { return *transportModel_; }
- // \}
-
- /*!
- * \name Restart mechanism
- */
- // \{
-
- /*!
- * \brief This method writes the complete state of the problem
- * to the harddisk.
- *
- * The file will start with the prefix returned by the name()
- * method, has the current time of the simulation clock in it's
- * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
- * file.) See Restart for details.
- */
- void serialize()
- {
- using Restarter = Restart;
-
- Restarter res;
- res.serializeBegin(asImp_());
- std::cout << "Serialize to file " << res.fileName() << "\n";
-
- timeManager().serialize(res);
- resultWriter().serialize(res);
-
- // do the actual serialization process: write primary variables
- res.template serializeEntities<0> (*pressModel_, gridView());
- res.template serializeEntities<0> (*transportModel_, gridView());
-
- res.serializeEnd();
-
- if (adaptiveGrid)
- {
- AdaptiveGridRestart<Grid>::serializeGrid(asImp_());
- }
- }
-
- /*!
- * \brief This method restores the complete state of the problem
- * from disk.
- *
- * It is the inverse of the serialize() method.
- * @param tRestart Restart time
- */
- void restart(const double tRestart)
- {
- if (adaptiveGrid)
- {
- AdaptiveGridRestart<Grid>::restartGrid(asImp_());
- variables().initialize();
- model().initialize();
- }
-
- using Restarter = Restart;
-
- Restarter res;
- res.deserializeBegin(asImp_(), tRestart);
- std::cout << "Deserialize from file " << res.fileName() << "\n";
-
- timeManager().deserialize(res);
-
- resultWriter().deserialize(res);
-
- // do the actual serialization process: get primary variables
- res.template deserializeEntities<0> (*pressModel_, gridView());
- res.template deserializeEntities<0> (*transportModel_, gridView());
- pressureModel().updateMaterialLaws();
-
- res.deserializeEnd();
- }
-
- void addOutputVtkFields()
- {
- }
- /*!
- * \brief Returns the vtk output verbosity level
- *
- * Level is set by property or input file.
- */
- int vtkOutputLevel() const
- {
- return vtkOutputLevel_;
- }
-
- //! Write the fields current solution into an VTK output file.
- void writeOutput(bool verbose = true)
- {
- if (verbose && gridView().comm().rank() == 0)
- std::cout << "Writing result file for current time step\n";
-
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
- if (adaptiveGrid)
- resultWriter_->gridChanged();
- resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
- model().addOutputVtkFields(*resultWriter_);
- asImp_().addOutputVtkFields();
- resultWriter_->endWrite();
- }
-
- // \}
-
-protected:
- //! Returns the applied VTK-writer for the output
- VtkMultiWriter& resultWriter()
- {
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
- return *resultWriter_;
- }
- //! \copydoc IMPETProblem::resultWriter()
- VtkMultiWriter& resultWriter() const
- {
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
- return *resultWriter_;
- }
-
-private:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \copydoc IMPETProblem::asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
-
- std::string simname_; // a string for the name of the current simulation,
- // which could be set by means of a program argument, for example.
-
- // pointer to a possibly adaptive grid.
- Grid *grid_;
- GridView gridView_;
-
- GlobalPosition bBoxMin_;
- GlobalPosition bBoxMax_;
-
- TimeManager *timeManager_;
- Scalar maxTimeStepSize_;
-
- std::shared_ptr<Variables> variables_;
-
- std::shared_ptr<PressureModel> pressModel_;//!< object including the pressure model
- std::shared_ptr<TransportModel> transportModel_;//!< object including the saturation model
- std::shared_ptr<IMPETModel> model_;
-
- std::shared_ptr<VtkMultiWriter> resultWriter_;
- int outputInterval_;
- Scalar outputTimeInterval_;
- int vtkOutputLevel_;
- std::shared_ptr<GridAdaptModel> gridAdapt_;
-
- Scalar dtVariationRestrictionFactor_;
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/impetproperties.hh b/dumux/porousmediumflow/sequential/impetproperties.hh
deleted file mode 100644
index 699d030dd3b08263e5cff9c3d03d17c3a8df4aed..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/impetproperties.hh
+++ /dev/null
@@ -1,74 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_IMPET_PROPERTIES_HH
-#define DUMUX_IMPET_PROPERTIES_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
-#include <dumux/porousmediumflow/sequential/transportproperties.hh>
-
-/*!
- * \ingroup IMPET
- * \ingroup IMPETProperties
- */
-/*!
- * \file
- * \brief Base file for properties related to sequential IMPET algorithms
- */
-namespace Dumux
-{
-
-template<class TypeTag>
-class IMPET;
-
-namespace Properties
-{
-/*!
- *
- * \brief General properties for sequential IMPET algorithms
- *
- * This class holds properties necessary for the sequential IMPET solution.
- */
-
-//////////////////////////////////////////////////////////////////
-// Type tags tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for models based on the diffusion-scheme
-// Create new type tags
-namespace TTag {
-struct IMPET { using InheritsFrom = std::tuple<SequentialModel>; };
-} // end namespace TTag
-}
-}
-
-#include <dumux/porousmediumflow/sequential/impet.hh>
-
-namespace Dumux
-{
-namespace Properties
-{
-//set impet model
-template<class TypeTag>
-struct Model<TypeTag, TTag::IMPET> { using type = IMPET<TypeTag>; };
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/mimetic/CMakeLists.txt b/dumux/porousmediumflow/sequential/mimetic/CMakeLists.txt
deleted file mode 100644
index c6a6dc215be96c6b8e39e0c712d2ed379b5f04e3..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/mimetic/CMakeLists.txt
+++ /dev/null
@@ -1,3 +0,0 @@
-file(GLOB DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_MIMETIC_HEADERS *.hh *.inc)
-install(FILES ${DUMUX_POROUSMEDIUMFLOW_SEQUENTIAL_MIMETIC_HEADERS}
- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/sequential/mimetic)
diff --git a/dumux/porousmediumflow/sequential/mimetic/properties.hh b/dumux/porousmediumflow/sequential/mimetic/properties.hh
deleted file mode 100644
index 1620faa647906623eb077ba5e37df1bd67b57ac5..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/mimetic/properties.hh
+++ /dev/null
@@ -1,67 +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 3 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/>. *
- *****************************************************************************/
-/*!
- * \ingroup IMPET
- * \ingroup IMPETProperties
- */
-/*!
- * \file
- *
- * \brief Defines the basic properties required for a mimetic method.
- */
-
-#ifndef DUMUX_MIMETICPROPERTIES_SEQUENTIAL_HH
-#define DUMUX_MIMETICPROPERTIES_SEQUENTIAL_HH
-
-//Dumux-includes
-#include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-namespace Dumux
-{
-
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-
-
-////////////////////////////////
-// properties
-////////////////////////////////
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for models using a mimetic method
-namespace TTag {
-struct Mimetic {};
-}
-
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-template<class TypeTag, class MyTypeTag>
-struct LocalStiffness { using type = UndefinedProperty; };
-
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/onemodelproblem.hh b/dumux/porousmediumflow/sequential/onemodelproblem.hh
deleted file mode 100644
index 840329ec37fca0f9e157f6ceb413de4044ea24b2..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/onemodelproblem.hh
+++ /dev/null
@@ -1,683 +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 3 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/>. *
- *****************************************************************************/
-
-#ifndef DUMUX_ONE_MODEL_PROBLEM_HH
-#define DUMUX_ONE_MODEL_PROBLEM_HH
-
-#include <dune/common/shared_ptr.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/porousmediumflow/sequential/properties.hh>
-#include <dumux/io/vtkmultiwriter.hh>
-#include <dumux/io/restart.hh>
-
-
-/**
- * @file
- * @brief Base class for definition of an sequential diffusion (pressure) or transport problem
- */
-
-namespace Dumux
-{
-
-/*! \ingroup IMPET
- *
- * @brief Base class for definition of an sequential diffusion (pressure) or transport problem
- *
- * @tparam TypeTag The Type Tag
- */
-template<class TypeTag>
-class OneModelProblem
-{
-private:
- using Implementation = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Grid = typename GridView::Grid;
-
- using TimeManager = GetPropType<TypeTag, Properties::TimeManager>;
-
- using VtkMultiWriter = Dumux::VtkMultiWriter<GridView>;
-
- using Variables = GetPropType<TypeTag, Properties::Variables>;
-
- using Model = GetPropType<TypeTag, Properties::Model>;
-
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using VertexMapper = typename SolutionTypes::VertexMapper;
- using ElementMapper = typename SolutionTypes::ElementMapper;
-
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
- enum
- {
- wetting = 0, nonwetting = 1
- };
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
- using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
-
- // private!! copy constructor
- OneModelProblem(const OneModelProblem&)
- {}
-
-public:
-
- //! Constructs an object of type OneModelProblemProblem
- /*!
- * \tparam TypeTag The TypeTag
- * \tparam verbose Output level for TimeManager
- */
- OneModelProblem(Grid& grid, bool verbose = true)
- : gridView_(grid.leafGridView()),
- bBoxMin_(std::numeric_limits<double>::max()),
- bBoxMax_(-std::numeric_limits<double>::max()),
- variables_(grid.leafGridView()),
- outputInterval_(1),
- outputTimeInterval_(0)
- {
- // calculate the bounding box of the grid view
- using std::max;
- using std::min;
- for (const auto& vertex : vertices(grid.leafGridView())) {
- for (int i=0; i<dim; i++) {
- bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
- bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
- }
- }
-
- timeManager_ = std::make_shared<TimeManager>(verbose);
-
- model_ = std::make_shared<Model>(asImp_()) ;
- maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
- }
-
- //! Constructs an object of type OneModelProblemProblem
- /*!
- * \tparam TypeTag The TypeTag
- * \tparam verbose Output level for TimeManager
- */
- OneModelProblem(TimeManager& timeManager, Grid& grid)
- : gridView_(grid.leafGridView()),
- bBoxMin_(std::numeric_limits<double>::max()),
- bBoxMax_(-std::numeric_limits<double>::max()),
- variables_(grid.leafGridView()),
- outputInterval_(1),
- outputTimeInterval_(0)
- {
- // calculate the bounding box of the grid view
- using std::max;
- using std::min;
- for (const auto& vertex : vertices(grid.leafGridView())) {
- for (int i=0; i<dim; i++) {
- bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
- bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
- }
- }
-
- timeManager_ = Dune::stackobject_to_shared_ptr<TimeManager>(timeManager);
-
- model_ = std::make_shared<Model>(asImp_()) ;
- maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
- }
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param bcTypes The boundary types for the conservation equations
- * \param intersection The intersection for which the boundary type is set
- */
- void boundaryTypes(BoundaryTypes &bcTypes,
- const Intersection &intersection) const
- {
- // forward it to the method which only takes the global coordinate
- asImp_().boundaryTypesAtPos(bcTypes, intersection.geometry().center());
- }
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param bcTypes The boundary types for the conservation equations
- * \param globalPos The position of the center of the boundary intersection
- */
- void boundaryTypesAtPos(BoundaryTypes &bcTypes,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Dirichlet conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem does not provide "
- "a boundaryTypesAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param intersection The boundary intersection
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichlet(PrimaryVariables &values,
- const Intersection &intersection) const
- {
- // forward it to the method which only takes the global coordinate
- asImp_().dirichletAtPos(values, intersection.geometry().center());
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position of the center of the boundary intersection
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Dirichlet conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem specifies that some boundary "
- "segments are dirichlet, but does not provide "
- "a dirichletAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
- * \param intersection The boundary intersection
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumann(PrimaryVariables &values,
- const Intersection &intersection) const
- {
- // forward it to the interface with only the global position
- asImp_().neumannAtPos(values, intersection.geometry().center());
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
- * \param globalPos The position of the center of the boundary intersection
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumannAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // Throw an exception (there is no reasonable default value
- // for Neumann conditions)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem specifies that some boundary "
- "segments are neumann, but does not provide "
- "a neumannAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the source term
- *
- * \param values The source and sink values for the conservation equations
- * \param element The element
- *
- * For this method, the \a values parameter stores the rate mass
- * generated or annihilate per volume unit. Positive values mean
- * that mass is created, negative ones mean that it vanishes.
- */
- void source(PrimaryVariables &values,
- const Element &element) const
- {
- // forward to generic interface
- asImp_().sourceAtPos(values, element.geometry().center());
- }
-
- /*!
- * \brief Evaluate the source term for all phases within a given
- * sub-control-volume.
- *
- * \param values The source and sink values for the conservation equations
- * \param globalPos The position of the center of the finite volume
- * for which the source term ought to be
- * specified in global coordinates
- *
- * For this method, the \a values parameter stores the rate mass
- * generated or annihilate per volume unit. Positive values mean
- * that mass is created, negative ones mean that it vanishes.
- */
- void sourceAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- { // Throw an exception (there is no initial condition)
- DUNE_THROW(Dune::InvalidStateException,
- "The problem does not provide "
- "a sourceAtPos() method.");
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The initial values for the primary variables
- * \param element The element
- *
- * For this method, the \a values parameter stores primary
- * variables.
- */
- void initial(PrimaryVariables &values,
- const Element &element) const
- {
- // forward to generic interface
- asImp_().initialAtPos(values, element.geometry().center());
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position of the center of the finite volume
- * for which the initial values ought to be
- * set (in global coordinates)
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- // initialize with 0 by default
- values = 0;
- }
-
- /*!
- * \brief Called by the TimeManager in order to
- * initialize the problem.
- */
- void init()
- {
- // set the initial condition of the model
- variables_.initialize();
- model().initialize();
- }
-
- /*!
- * \brief Called by TimeManager just before the time
- * integration.
- */
- void preTimeStep()
- {}
-
- /*!
- * \brief Called by TimeManager in order to do a time
- * integration on the model.
- */
- void timeIntegration()
- {}
-
- /*!
- * \brief Called by TimeManager whenever a solution for a
- * timestep has been computed and the simulation time has
- * been updated.
- *
- * This is used to do some janitorial tasks like writing the
- * current solution to disk.
- */
- void postTimeStep()
- {}
-
- /*!
- * \brief Called by the time manager after everything which can be
- * done about the current time step is finished and the
- * model should be prepared to do the next time integration.
- */
- void advanceTimeLevel()
- {}
-
- /*!
- * \brief Returns the user specified maximum time step size
- *
- * Overload in problem for custom needs.
- */
- Scalar maxTimeStepSize() const
- { return maxTimeStepSize_; }
-
- /*!
- * \brief Returns the current time step size [seconds].
- */
- Scalar timeStepSize() const
- { return timeManager().timeStepSize(); }
-
- /*!
- * \brief Sets the current time step size [seconds].
- */
- void setTimeStepSize(Scalar dt)
- { timeManager().setTimeStepSize(dt); }
-
- /*!
- * \brief Called by TimeManager whenever a solution for a
- * timestep has been computed and the simulation time has
- * been updated.
- */
- Scalar nextTimeStepSize(Scalar dt)
- { return timeManager().timeStepSize();}
-
- /*!
- * \brief Returns true if a restart file should be written to
- * disk.
- *
- * The default behaviour is to write one restart file every 5 time
- * steps. This file is intented to be overwritten by the
- * implementation.
- */
- bool shouldWriteRestartFile() const
- {
- return
- timeManager().timeStepIndex() > 0 &&
- (timeManager().timeStepIndex() % 5 == 0);
- }
-
- /*!
- * \brief Sets a time interval for Output
- *
- * The default is 0.0 -> Output determined by output number interval (<tt>setOutputInterval(int)</tt>)
- */
- void setOutputTimeInterval(const Scalar timeInterval)
- {
- outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
- timeManager().startNextEpisode(outputTimeInterval_);
- }
-
- /*!
- * \brief Sets the interval for Output
- *
- * The default is 1 -> Output every time step
- */
- void setOutputInterval(int interval)
- { outputInterval_ = interval; }
-
- /*!
- * \brief Returns true if the current solution should be written to
- * disk (i.e. as a VTK file)
- *
- * The default behaviour is to write out every the solution for
- * very time step. This file is intented to be overwritten by the
- * implementation.
- */
-
- bool shouldWriteOutput() const
- {
- if (outputInterval_ > 0)
- {
- if (timeManager().timeStepIndex() % outputInterval_ == 0
- || timeManager().willBeFinished()
- || timeManager().episodeWillBeFinished())
- {
- return true;
- }
- }
- else if (timeManager().willBeFinished()
- || timeManager().episodeWillBeFinished() || timeManager().timeStepIndex() == 0)
- {
- return true;
- }
- return false;
- }
-
- void addOutputVtkFields()
- {}
-
- //! Write the fields current solution into an VTK output file.
- void writeOutput(bool verbose = true)
- {
- if (verbose && gridView().comm().rank() == 0)
- std::cout << "Writing result file for current time step\n";
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
- resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
- model().addOutputVtkFields(*resultWriter_);
- asImp_().addOutputVtkFields();
- resultWriter_->endWrite();
- }
-
- /*!
- * \brief Called when the end of an simulation episode is reached.
- */
- void episodeEnd()
- {
- if (outputTimeInterval_ > 0.0 && !timeManager().finished())
- {
- timeManager().startNextEpisode(outputTimeInterval_);
- }
- else if (!timeManager().finished())
- {
- std::cerr << "The end of an episode is reached, but the problem "
- << "does not override the episodeEnd() method. "
- << "Doing nothing!\n";
- }
- }
-
- // \}
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- * It could be either overwritten by the problem files, or simply
- * declared over the setName() function in the application file.
- */
- const std::string& name() const
- {
- return simname_;
- }
-
- /*!
- * \brief Set the problem name.
- *
- * This function sets the simulation name, which should be called before
- * the application problem is declared! If not, the default name "sim"
- * will be used.
- */
- void setName(const std::string& newName)
- {
- simname_ = newName;
- }
-
- /*!
- * \brief The GridView which used by the problem.
- */
- const GridView &gridView() const
- { return gridView_; }
-
- /*!
- * \brief Returns the mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return variables_.vertexMapper(); }
-
- /*!
- * \brief Returns the mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return variables_.elementMapper(); }
-
- /*!
- * \brief The coordinate of the corner of the GridView's bounding
- * box with the smallest values.
- */
- const GlobalPosition &bBoxMin() const
- { return bBoxMin_; }
-
- /*!
- * \brief The coordinate of the corner of the GridView's bounding
- * box with the largest values.
- */
- const GlobalPosition &bBoxMax() const
- { return bBoxMax_; }
-
- /*!
- * \brief Returns TimeManager object used by the simulation
- */
- TimeManager &timeManager()
- { return *timeManager_; }
-
- /*!
- * \brief \copybrief OneModelProblem::timeManager()
- */
- const TimeManager &timeManager() const
- { return *timeManager_; }
-
- /*!
- * \brief Returns variables object.
- */
- Variables& variables ()
- { return variables_; }
-
- /*!
- * \brief \copybrief OneModelProblem::variables()
- */
- const Variables& variables () const
- { return variables_; }
-
- /*!
- * \brief Returns numerical model used for the problem.
- */
- Model &model()
- { return *model_; }
-
- /*!
- * \brief \copybrief OneModelProblem::model()
- */
- const Model &model() const
- { return *model_; }
- // \}
-
-
- /*!
- * \name Restart mechanism
- */
- // \{
-
- /*!
- * \brief This method writes the complete state of the problem
- * to the harddisk.
- *
- * The file will start with the prefix returned by the name()
- * method, has the current time of the simulation clock in it's
- * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
- * file.) See Restart for details.
- */
- void serialize()
- {
- using Restarter = Restart;
-
- Restarter res;
- res.serializeBegin(asImp_());
- std::cout << "Serialize to file " << res.fileName() << "\n";
-
- timeManager().serialize(res);
- resultWriter().serialize(res);
- res.template deserializeEntities<0> (model(), gridView_);
-
- res.serializeEnd();
- }
-
- /*!
- * \brief This method restores the complete state of the problem
- * from disk.
- *
- * It is the inverse of the serialize() method.
- */
- void restart(double tRestart)
- {
- using Restarter = Restart;
-
- Restarter res;
- res.deserializeBegin(asImp_(), tRestart);
- std::cout << "Deserialize from file " << res.fileName() << "\n";
-
- timeManager().deserialize(res);
- resultWriter().deserialize(res);
- res.template deserializeEntities<0> (model(), gridView_);
-
- res.deserializeEnd();
- }
-
- // \}
-
-protected:
- VtkMultiWriter& resultWriter()
- {
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView_, asImp_().name());
- return *resultWriter_;
- }
-
- VtkMultiWriter& resultWriter() const
- {
- if (!resultWriter_)
- resultWriter_ = std::make_shared<VtkMultiWriter>(gridView_, asImp_().name());
- return *resultWriter_;
- }
-
-private:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation &asImp_()
- { return *static_cast<Implementation *>(this); }
-
- //! \brief \copybrief OneModelProblem::asImp_()
- const Implementation &asImp_() const
- { return *static_cast<const Implementation *>(this); }
-
- std::string simname_; // a string for the name of the current simulation,
- // which could be set by means of an program argument,
- // for example.
- const GridView gridView_;
-
- GlobalPosition bBoxMin_;
- GlobalPosition bBoxMax_;
-
- std::shared_ptr<TimeManager> timeManager_;
- Scalar maxTimeStepSize_;
-
- Variables variables_;
-
- std::shared_ptr<Model> model_;
-
- std::shared_ptr<VtkMultiWriter> resultWriter_;
- int outputInterval_;
- Scalar outputTimeInterval_;
-};
-
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/pressureproperties.hh b/dumux/porousmediumflow/sequential/pressureproperties.hh
deleted file mode 100644
index ef76f33577192e90efc582ff564b0cf18bf28d98..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/pressureproperties.hh
+++ /dev/null
@@ -1,116 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_PRESSURE_PROPERTIES_HH
-#define DUMUX_PRESSURE_PROPERTIES_HH
-
-//Dune-includes
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include "properties.hh"
-#include <dumux/linear/seqsolverbackend.hh>
-
-
-/*!
- * \ingroup Sequential
- * \ingroup IMPETProperties
- */
-/*!
- * \file
- * \brief Base file for properties related to sequential IMPET algorithms
- */
-namespace Dumux
-{
-namespace Properties
-{
-/*!
- *
- * \brief General properties for sequential IMPET algorithms
- *
- * This class holds properties necessary for the sequential IMPET solution.
- */
-
-//////////////////////////////////////////////////////////////////
-// Type tags tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for models based on the diffusion-scheme
-// Create new type tags
-namespace TTag {
-struct Pressure { using InheritsFrom = std::tuple<SequentialModel>; };
-} // end namespace TTag
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-//Properties for linear solvers
-template<class TypeTag, class MyTypeTag>
-struct PressureRHSVector { using type = UndefinedProperty; };//!< Type of the right hand side vector given to the linear solver
-template<class TypeTag, class MyTypeTag>
-struct PressureSolutionVector { using type = UndefinedProperty; };//!Type of solution vector or pressure system
-template<class TypeTag, class MyTypeTag>
-struct VisitFacesOnlyOnce { using type = UndefinedProperty; }; //!< Indicates if faces are only regarded from one side
-}
-}
-
-#include <dumux/porousmediumflow/sequential/cellcentered/velocitydefault.hh>
-
-namespace Dumux
-{
-namespace Properties
-{
-//! Faces are only regarded from one side and not from both cells
-template<class TypeTag>
-struct VisitFacesOnlyOnce<TypeTag, TTag::Pressure> { static constexpr bool value = false; };
-
-//Set defaults
-template<class TypeTag>
-struct PressureCoefficientMatrix<TypeTag, TTag::Pressure>
-{
-private:
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using MB = Dune::FieldMatrix<Scalar, 1, 1>;
-
-public:
- using type = Dune::BCRSMatrix<MB>;
-};
-template<class TypeTag>
-struct PressureRHSVector<TypeTag, TTag::Pressure>
-{
-private:
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-
-public:
- using type = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
-};
-
-template<class TypeTag>
-struct PressureSolutionVector<TypeTag, TTag::Pressure> { using type = typename GetProp<TypeTag, SolutionTypes>::ScalarSolution; };
-
-// use the stabilized BiCG solver preconditioned by the ILU-0 by default
-template<class TypeTag>
-struct LinearSolver<TypeTag, TTag::Pressure> { using type = ILU0BiCGSTABBackend ; };
-
-template<class TypeTag>
-struct Velocity<TypeTag, TTag:: Pressure> { using type = FVVelocityDefault<TypeTag>; };
-
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/properties.hh b/dumux/porousmediumflow/sequential/properties.hh
deleted file mode 100644
index 8fc137e275ab469561ff5b70d02df7505348bba2..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/properties.hh
+++ /dev/null
@@ -1,281 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_SEQUENTIAL_PROPERTIES_HH
-#define DUMUX_SEQUENTIAL_PROPERTIES_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/properties/model.hh>
-#include <dumux/common/properties/grid.hh>
-#include <dumux/common/defaultmappertraits.hh>
-#include <dumux/discretization/method.hh>
-#include <dumux/porousmediumflow/sequential/boundarytypes.hh>
-#include <dumux/porousmediumflow/sequential/gridadaptproperties.hh>
-#include <dumux/porousmediumflow/sequential/gridadaptinitializationindicatordefault.hh>
-
-
-/*!
- * \ingroup Sequential
- * \ingroup IMPETProperties
- */
-/*!
- * \file
- * \brief Base file for properties related to sequential models
- */
-namespace Dumux
-{
-namespace Properties
-{
-
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! Create a type tag for all sequential models
-// Create new type tags
-namespace TTag {
-struct SequentialModel { using InheritsFrom = std::tuple<GridAdapt, GridProperties, ModelProperties>; };
-} // end namespace TTag
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-//! Property tag for types associated with the solution of the PDE.
-//! This means vectors of primary variables, solution functions on the
-//! grid, and elements, and shape functions.
-template<class TypeTag, class MyTypeTag>
-struct SolutionTypes { using type = UndefinedProperty; };
-template<class TypeTag, class MyTypeTag>
-struct Indices { using type = UndefinedProperty; };
-template<class TypeTag, class MyTypeTag>
-struct SequentialBoundaryTypes { using type = UndefinedProperty; };
-
-// Some properties that have been removed from numeric model
-template<class TypeTag, class MyTypeTag>
-struct Model { using type = UndefinedProperty; }; //!< The type of the mode
-template<class TypeTag, class MyTypeTag>
-struct DiscretizationMethod { using type = UndefinedProperty; }; //!< The type of discretization method
-
-template<class TypeTag, class MyTypeTag>
-struct PressureModel { using type = UndefinedProperty; }; //!< The type of the discretization of a pressure model
-template<class TypeTag, class MyTypeTag>
-struct TransportModel { using type = UndefinedProperty; }; //!< The type of the discretization of a transport model
-template<class TypeTag, class MyTypeTag>
-struct Velocity { using type = UndefinedProperty; }; //!< The type velocity reconstruction
-template<class TypeTag, class MyTypeTag>
-struct NumEq { using type = UndefinedProperty; }; //!< Number of equations in the system of PDEs
-template<class TypeTag, class MyTypeTag>
-struct NumPhases { using type = UndefinedProperty; }; //!< Number of phases in the system
-template<class TypeTag, class MyTypeTag>
-struct NumComponents { using type = UndefinedProperty; }; //!< Number of components in the system
-template<class TypeTag, class MyTypeTag>
-struct Variables { using type = UndefinedProperty; }; //!< The type of the container of global variables
-template<class TypeTag, class MyTypeTag>
-struct CellData { using type = UndefinedProperty; };//!< Defines data object to be stored
-template<class TypeTag, class MyTypeTag>
-struct MaxIntersections { using type = UndefinedProperty; }; //!< Gives maximum number of intersections of an element and neighboring elements
-template<class TypeTag, class MyTypeTag>
-struct PressureCoefficientMatrix { using type = UndefinedProperty; }; //!< Gives maximum number of intersections of an element and neighboring elements
-}
-}
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dumux/common/timemanager.hh>
-#include <dumux/common/boundarytypes.hh>
-#include<dumux/common/boundaryconditions.hh>
-#include<dumux/discretization/method.hh>
-
-namespace Dumux
-{
-
-template<class TypeTag>
-class VariableClass;
-
-namespace Properties
-{
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-//! Type of the jacobian matrix needed for compatibility with implicit models for the amg backend
-template<class TypeTag>
-struct JacobianMatrix<TypeTag, TTag::SequentialModel> { using type = GetPropType<TypeTag, Properties::PressureCoefficientMatrix>; };
-
-//! Dummy model traits for compatibility with the rest of dumux
-//! until the sequential models are incorporated into the general framework
-template<class TypeTag>
-struct ModelTraits<TypeTag, TTag::SequentialModel>
-{
-private:
- struct DummyTraits
- {
- using Indices = GetPropType<TypeTag, Properties::Indices>;
- static constexpr int numEq() { return getPropValue<TypeTag, Properties::NumEq>(); }
- };
-public:
- using type = DummyTraits;
-};
-
-//! Default number of intersections for quadrilaterals
-template<class TypeTag>
-struct MaxIntersections<TypeTag, TTag::SequentialModel>
-{
-private:
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
-public:
- static constexpr int value = 2*GridView::dimension;
-};
-
-//! A simplified grid geometry for compatibility with new style models
-template<class TypeTag>
-struct GridGeometry<TypeTag, TTag::SequentialModel>
-{
- using GV = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
- struct MockFVGridGeometry
- : public DefaultMapperTraits<GV>
- {
- static constexpr Dumux::DiscretizationMethod discMethod = Dumux::DiscretizationMethod::cctpfa;
- using GridView = GV;
- };
-public:
- using type = MockFVGridGeometry;
-};
-
-//! For compatibility with new style models we need a solution vector type
-template<class TypeTag>
-struct SolutionVector<TypeTag, TTag::SequentialModel>
-{
-public:
- using type = typename GetProp<TypeTag, SolutionTypes>::ScalarSolution;
-};
-
-/*!
- * \brief Specifies the types which are assoicated with a solution.
- *
- * This means shape functions, solution vectors, etc.
- */
-template<class TypeTag>
-struct SolutionTypes<TypeTag, TTag::SequentialModel>
-{
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using Variables = GetPropType<TypeTag, Properties::Variables>;
-
- enum
- {
- dim = GridView::dimension,
- numEq = getPropValue<TypeTag, Properties::NumEq>(),
- numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
- numComponents = getPropValue<TypeTag, Properties::NumComponents>(),
- maxIntersections = getPropValue<TypeTag, Properties::MaxIntersections>()
- };
-
-public:
- /*!
- * \brief Mapper for the grid view's vertices.
- */
- using VertexMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-
- /*!
- * \brief Mapper for the grid view's elements.
- */
- using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-
- /*!
- * \brief The type of a solution at a fixed time.
- *
- * This defines the primary and secondary variable vectors at each degree of freedom.
- */
- using PrimaryVariables = Dune::FieldVector<Scalar, numEq>;
- //! type for vector of scalars
- using ScalarSolution = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
- //! type for vector of phase properties
- using ComponentProperty = Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numComponents>;
- //! type for vector of phase properties
- using PhaseProperty = Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numPhases>;
- //! type for vector of fluid properties: Vector[element][phase]
- using FluidProperty = Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numPhases>;
- //! type for vector of vectors (of size 2 x dimension) of scalars
- using PhasePropertyElemFace = Dune::BlockVector<Dune::FieldVector<Dune::FieldVector<Scalar, numPhases>, maxIntersections> >;
- //! type for vector of vectors (of size 2 x dimension) of vector (of size dimension) of scalars
- using DimVecElemFace = Dune::BlockVector<Dune::FieldVector<Dune::FieldVector<Scalar, dim>, maxIntersections> >;
-};
-
-template<class TypeTag>
-struct Variables<TypeTag, TTag::SequentialModel> { using type = VariableClass<TypeTag>; };
-
-template<class TypeTag>
-struct PrimaryVariables<TypeTag, TTag::SequentialModel> { using type = typename GetProp<TypeTag, SolutionTypes>::PrimaryVariables; };
-
-//! Set the default type for the time manager
-template<class TypeTag>
-struct TimeManager<TypeTag, TTag::SequentialModel> { using type = Dumux::TimeManager<TypeTag>; };
-
-/*!
- * \brief Boundary types at a single degree of freedom.
- */
-template<class TypeTag>
-struct SequentialBoundaryTypes<TypeTag, TTag::SequentialModel>
-{ private:
- enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
-public:
- using type = Dumux::SequentialBoundaryTypes<numEq>;
-};
-
-//! do not specific any model-specific default parameters here
-template<class TypeTag>
-struct ModelDefaultParameters<TypeTag, TTag::SequentialModel>
-{
- static void defaultParams(Dune::ParameterTree& params, const std::string& group = "")
- {
- params["GridAdapt.CoarsenTolerance"] = "0.001"; //!< tolerance for coarsening
- params["GridAdapt.EnableInitializationIndicator"] = "false"; //!< switch for using initial grid adaption
- params["GridAdapt.EnableMultiPointFluxApproximation"] = "true"; //!< apply an mpfa method around hanging nodes
- params["GridAdapt.MaxLevel"] = "1"; //!< maximum allowed level
- params["GridAdapt.MaxInteractionVolumes"] = "4"; //!< use up to four interaction regions
- params["GridAdapt.MinLevel"] = "0"; //!< minimum allowed level
- params["GridAdapt.RefineAtDirichletBC"] = "false"; //!< switch for refinement at Dirichlet BCs
- params["GridAdapt.RefineAtFluxBC"] = "false"; //!< switch for refinement at Neumann BCs
- params["GridAdapt.RefineAtSource"] = "false"; //!< switch for refinement at sources
- params["GridAdapt.RefineTolerance"] = "0.05"; //!< tolerance for refinement
-
- params["Impet.CFLFactor"] = "1.0"; //!< scalar factor for additional scaling of the time step
- params["Impet.EnableVolumeIntegral"] = "true"; //!< regard volume integral in pressure equation
- params["Impet.ErrorTermFactor"] = "0.5"; //!< scaling factor for the error term
- params["Impet.ErrorTermLowerBound"] = "0.1"; //!< lower threshold used for the error term evaluation
- params["Impet.ErrorTermUpperBound"] = "0.9"; //!< upper threshold used for the error term evaluation
- params["Impet.PorosityThreshold"] = "1e-6"; //!< porosity will be set to max(given value, threshold)
- params["Impet.SubCFLFactor"] = "1.0"; //!< scalar factor for scaling of local sub-time-step
- params["Impet.SwitchNormals"] = "false"; //!< don't switch direction of face normal vectors
-
- params["MPFA.CalcVelocityInTransport"] = "false"; //!< disable facewise velocity calculation
-
- params["TimeManager.SubTimestepVerbosity"] = "0"; //!< don't be verbose in local sub-time-steps
-
- params["Vtk.OutputLevel"] = "0"; //!< VTK output contains only the basic values
- }
-};
-
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/transportproperties.hh b/dumux/porousmediumflow/sequential/transportproperties.hh
deleted file mode 100644
index 2ca2edc99160f3737c01370c1e564bd9c5d13ef8..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/transportproperties.hh
+++ /dev/null
@@ -1,76 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_TRANSPORT_PROPERTIES_HH
-#define DUMUX_TRANSPORT_PROPERTIES_HH
-
-#include "properties.hh"
-
-/*!
- * \ingroup Sequential
- * \ingroup IMPETProperties
- */
-/*!
- * \file
- * \brief Specifies the properties for immiscible 2p transport
- */
-namespace Dumux
-{
-
-namespace Properties
-{
-// \{
-
-//////////////////////////////////////////////////////////////////
-// Type tags tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for models based on the diffusion-scheme
-// Create new type tags
-namespace TTag {
-struct Transport { using InheritsFrom = std::tuple<SequentialModel>; };
-} // end namespace TTag
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-template<class TypeTag, class MyTypeTag>
-struct TransportSolutionType { using type = UndefinedProperty; };
-template<class TypeTag, class MyTypeTag>
-struct EvalCflFluxFunction { using type = UndefinedProperty; }; //!< Type of the evaluation of the CFL-condition
-
-/*!
- * \brief Default implementation for the Vector of the transportet quantity
- *
- * This type defines the data type of the transportet quantity. In case of a
- * immiscible 2p system, this would represent a vector holding the saturation
- * of one phase.
- */
-template<class TypeTag>
-struct TransportSolutionType<TypeTag, TTag::Transport>
-{
- private:
- using SolutionType = GetProp<TypeTag, Properties::SolutionTypes>;
-
- public:
- using type = typename SolutionType::ScalarSolution;//!<type for vector of scalar properties
-};
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/sequential/variableclass.hh b/dumux/porousmediumflow/sequential/variableclass.hh
deleted file mode 100644
index ec3e242424fe1c044a1c3bbd3175c24d4a225d88..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/variableclass.hh
+++ /dev/null
@@ -1,191 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_VARIABLECLASS_HH
-#define DUMUX_VARIABLECLASS_HH
-
-#include "properties.hh"
-#include <dumux/common/deprecated.hh>
-
-// for parallelization
-//#include <dumux/parallel/elementhandles.hh>
-
-/**
- * @file
- * @brief Base class holding the variables for sequential models.
- */
-
-namespace Dumux
-{
-
-/*!
- * \ingroup Sequential
- */
-//! Base class holding the variables and discretized data for sequential models.
-/*!
- * Stores global information and variables that are common for all sequential models and also functions needed to access these variables.
- * Can be directly used for a single phase model.
- *
- * @tparam TypeTag The Type Tag
- *
- */
-template<class TypeTag>
-class VariableClass
-{
-private:
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using CellData = GetPropType<TypeTag, Properties::CellData>;
-
- enum
- {
- dim = GridView::dimension,
- };
-
- using Element = typename GridView::Traits::template Codim<0>::Entity;
- using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
-
- using VertexMapper = typename SolutionTypes::VertexMapper;
- using ElementMapper = typename SolutionTypes::ElementMapper;
-
-public:
- using CellDataVector = typename std::vector<CellData>;
-
-private:
- const GridView gridView_;
- ElementMapper elementMapper_;
- VertexMapper vertexMapper_;
- CellDataVector cellDataVector_;
-
-public:
- //! Constructs a VariableClass object
- /**
- * @param gridView a DUNE gridview object corresponding to diffusion and transport equation
- */
- VariableClass(const GridView& gridView) :
- gridView_(gridView),
- elementMapper_(gridView, Dune::mcmgElementLayout()),
- vertexMapper_(gridView, Dune::mcmgVertexLayout())
- {}
-
- //! Initializes the variable class
- /*! Method initializes the cellData vector.
- * Should be called from problem init()
- */
- void initialize()
- {
- if constexpr (Deprecated::hasUpdateGridView<ElementMapper, GridView>())
- elementMapper_.update(gridView_);
- else
- Deprecated::update(elementMapper_);
-
- if constexpr (Deprecated::hasUpdateGridView<VertexMapper, GridView>())
- vertexMapper_.update(gridView_);
- else
- Deprecated::update(vertexMapper_);
-
- cellDataVector_.resize(gridView_.size(0));
- }
-
- //! Resizes sequential variable vectors
- /*! Method that change the size of the vectors for h-adaptive simulations.
- *
- *\param size Size of the current (refined and coarsened) grid
- */
- void adaptVariableSize(const int size)
- {
- cellDataVector_.resize(size);
- }
-
- //! Return the vector holding all cell data
- CellDataVector& cellDataGlobal()
- {
- return cellDataVector_;
- }
-
- const CellDataVector& cellDataGlobal() const
- {
- assert(cellDataVector_.size() == gridView_.size(0));
-
- return cellDataVector_;
- }
-
- //! Return the cell data of a specific cell
- CellData& cellData(const int idx)
- {
- assert(cellDataVector_.size() == gridView_.size(0));
-
- return cellDataVector_[idx];
- }
-
- const CellData& cellData(const int idx) const
- {
- assert(cellDataVector_.size() == gridView_.size(0));
-
- return cellDataVector_[idx];
- }
-
- //! Get index of element (codim 0 entity)
- /*! Get index of element (codim 0 entity).
- * @param element codim 0 entity
- * \return element index
- */
- int index(const Element& element) const
- {
- return elementMapper_.index(element);
- }
-
- //! Get index of vertex (codim dim entity)
- /*! Get index of vertex (codim dim entity).
- * @param vertex codim dim entity
- * \return vertex index
- */
- int index(const Vertex& vertex) const
- {
- return vertexMapper_.index(vertex);
- }
-
- //!Return gridView
- const GridView& gridView() const
- {
- return gridView_;
- }
-
- //! Return mapper for elements (for adaptive grids)
- ElementMapper& elementMapper()
- {
- return elementMapper_;
- }
- //! Return mapper for elements (for static grids)
- const ElementMapper& elementMapper() const
- {
- return elementMapper_;
- }
- //! Return mapper for vertices (for adaptive grids)
- VertexMapper& vertexMapper()
- {
- return vertexMapper_;
- }
- //! Return mapper for vertices (for static grids)
- const VertexMapper& vertexMapper() const
- {
- return vertexMapper_;
- }
-};
-}
-#endif
diff --git a/dumux/porousmediumflow/sequential/variableclassadaptive.hh b/dumux/porousmediumflow/sequential/variableclassadaptive.hh
deleted file mode 100644
index cc0f65aa9af8b698cd86b22b5b7ef712726bd0e2..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/sequential/variableclassadaptive.hh
+++ /dev/null
@@ -1,213 +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 3 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/>. *
- *****************************************************************************/
-#ifndef DUMUX_VARIABLECLASS_ADAPTIVE_HH
-#define DUMUX_VARIABLECLASS_ADAPTIVE_HH
-
-#include <dune/grid/common/partitionset.hh>
-#include <dune/grid/utility/persistentcontainer.hh>
-#include <dumux/parallel/vectorcommdatahandle.hh>
-#include "variableclass.hh"
-
-/**
- * @file
- * @brief Base class holding the variables for sequential models.
- */
-
-namespace Dumux
-{
-/*!
- * \ingroup IMPET
- */
-//! Base class holding the variables and discretized data for sequential models.
-/*!
- * Stores global information and variables that are common for all sequential models and also functions needed to access these variables.
- * Can be directly used for a single phase model.
- *
- * @tparam TypeTag The Type Tag
- *
- */
-template<class TypeTag>
-class VariableClassAdaptive: public VariableClass<TypeTag>
-{
-private:
- using ParentType = VariableClass<TypeTag>;
-
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
- using CellData = GetPropType<TypeTag, Properties::CellData>;
- using AdaptedValues = typename CellData::AdaptedValues;
-
- using Grid = typename GridView::Grid;
- using LevelGridView = typename Grid::LevelGridView;
- using PersistentContainer = Dune::PersistentContainer<Grid, AdaptedValues>;
-
-private:
- const Grid& grid_;
- PersistentContainer adaptationMap_;
-
-public:
- //! Constructs an adaptive VariableClass object
- /**
- * In addition to providing a storage object for cell-centered Methods, this class provides
- * mapping functionality to adapt the grid.
- *
- * @param gridView a DUNE gridview object corresponding to diffusion and transport equation
- */
- VariableClassAdaptive(const GridView& gridView) :
- ParentType(gridView), grid_(gridView.grid()), adaptationMap_(grid_, 0)
- {}
-
-
- /*!
- * Store primary variables
- *
- * To reconstruct the solution in father elements, problem properties might
- * need to be accessed.
- * From upper level on downwards, the old solution is stored into an container
- * object, before the grid is adapted. Father elements hold averaged information
- * from the son cells for the case of the sons being coarsened.
- *
- * @param problem The current problem
- */
- void storePrimVars(const Problem& problem)
- {
- adaptationMap_.resize();
-
- // loop over all levels of the grid
- for (int level = grid_.maxLevel(); level >= 0; level--)
- {
- //get grid view on level grid
- LevelGridView levelView = grid_.levelGridView(level);
-
- for (const auto& element : elements(levelView))
- {
- //get your map entry
- AdaptedValues &adaptedValues = adaptationMap_[element];
-
- // put your value in the map
- if (element.isLeaf())
- {
- // get index
- int indexI = this->index(element);
-
- CellData& cellData = this->cellData(indexI);
-
- cellData.storeAdaptionValues(adaptedValues, element);
-
- adaptedValues.count = 1;
- }
- //Average in father
- if (element.level() > 0)
- {
- auto father = element.father();
- AdaptedValues& adaptedValuesFather = adaptationMap_[father];
- adaptedValuesFather.count += 1;
- CellData::storeAdaptionValues(adaptedValues, adaptedValuesFather, father);
- }
- }
- }
- }
-
- /*!
- * Reconstruct missing primary variables (where elements are created/deleted)
- *
- * To reconstruct the solution in father elements, problem properties might
- * need to be accessed.
- * Starting from the lowest level, the old solution is mapped on the new grid:
- * Where coarsened, new cells get information from old father element.
- * Where refined, a new solution is reconstructed from the old father cell,
- * and then a new son is created. That is then stored into the general data
- * structure (CellData).
- *
- * @param problem The current problem
- */
- void reconstructPrimVars(const Problem& problem)
- {
- adaptationMap_.resize();
-
- for (int level = 0; level <= grid_.maxLevel(); level++)
- {
- LevelGridView levelView = grid_.levelGridView(level);
-
- for (const auto& element : elements(levelView))
- {
- // only treat non-ghosts, ghost data is communicated afterwards
- if (element.partitionType() == Dune::GhostEntity)
- continue;
-
- if (!element.isNew())
- {
- //entry is in map, write in leaf
- if (element.isLeaf())
- {
- AdaptedValues &adaptedValues = adaptationMap_[element];
- int newIdxI = this->index(element);
-
- CellData& cellData = this->cellData(newIdxI);
-
- cellData.setAdaptionValues(adaptedValues, element);
- }
- }
- else
- {
- // value is not in map, interpolate from father element
- if (element.level() > 0)
- {
- // create new entry: reconstruct from adaptationMap_[father] to a new
- // adaptationMap_[son]
- CellData::reconstructAdaptionValues(adaptationMap_, element.father(), element, problem);
-
- // access new son
- AdaptedValues& adaptedValues = adaptationMap_[element];
- adaptedValues.count = 1;
-
- // if we are on leaf, store reconstructed values of son in CellData object
- if (element.isLeaf())
- {
- // acess new CellData object
- int newIdxI = this->index(element);
- CellData& cellData = this->cellData(newIdxI);
-
- cellData.setAdaptionValues(adaptedValues, element);
- }
- }
- }
- }
-
- }
- // reset entries in restrictionmap
- adaptationMap_.resize( typename PersistentContainer::Value() );
- adaptationMap_.shrinkToFit();
- adaptationMap_.fill( typename PersistentContainer::Value() );
-
-#if HAVE_MPI
- // communicate ghost data
- using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
- using ElementMapper = typename SolutionTypes::ElementMapper;
- using DataHandle = VectorCommDataHandleEqual<ElementMapper, std::vector<CellData>, 0/*elementCodim*/>;
- DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
- problem.gridView().template communicate<DataHandle>(dataHandle,
- Dune::InteriorBorder_All_Interface,
- Dune::ForwardCommunication);
-#endif
- }
-
-};
-}
-#endif