diff --git a/dumux/adaptive/CMakeLists.txt b/dumux/adaptive/CMakeLists.txt
index a9fe39cd98c6232082fe6da169e535ef117768ff..5f56a25371676e4a1b0700f6431e71e20e7e93e6 100644
--- a/dumux/adaptive/CMakeLists.txt
+++ b/dumux/adaptive/CMakeLists.txt
@@ -2,5 +2,6 @@
install(FILES
adapt.hh
griddatatransfer.hh
+initializationindicator.hh
markelements.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/adaptive)
diff --git a/dumux/adaptive/initializationindicator.hh b/dumux/adaptive/initializationindicator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3571568892ef3f86f68a48099467ddaba622980b
--- /dev/null
+++ b/dumux/adaptive/initializationindicator.hh
@@ -0,0 +1,289 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/**
+* \file
+* \brief Class defining an initialization indicator for grid adaption
+*/
+#ifndef DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
+#define DUMUX_GRIDADAPTINITIALIZATIONINDICATOR_HH
+
+#include <dune/geometry/type.hh>
+
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+
+/*!\ingroup GridAdapt
+ * \brief Class defining an initialization indicator for grid adaption.
+ * Accounts for sources and boundaries. Only for grid initialization!
+ *
+ * \tparam TypeTag The problem TypeTag
+ */
+template<class TypeTag>
+class GridAdaptInitializationIndicator
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::Traits::template Codim<0>::Entity;
+
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+
+public:
+
+ /*! \brief Constructor
+ *
+ * \param problem The problem object
+ * \param fvGridGeometry The finite volume geometry of the grid
+ * \param gridVariables The secondary variables on the grid
+ */
+ GridAdaptInitializationIndicator(std::shared_ptr<const Problem> problem,
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<const GridVariables> gridVariables)
+ : problem_(problem)
+ , fvGridGeometry_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , minLevel_(getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.MinLevel"))
+ , maxLevel_(getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.MaxLevel"))
+ , refineAtDirichletBC_(getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.RefineAtDirichletBC", true))
+ , refineAtFluxBC_(getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.RefineAtFluxBC", true))
+ , refineAtSource_(getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.RefineAtSource", true))
+ , eps_(getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.BCRefinementThreshold", 1e-10))
+ {}
+
+ /*!
+ * \brief Function to set the minimum allowed level.
+ *
+ */
+ void setMinLevel(std::size_t minLevel)
+ {
+ minLevel_ = minLevel;
+ }
+
+ /*!
+ * \brief Function to set the maximum allowed level.
+ *
+ */
+ void setMaxLevel(std::size_t maxLevel)
+ {
+ maxLevel_ = maxLevel;
+ }
+
+ /*!
+ * \brief Function to set the minumum/maximum allowed levels.
+ *
+ */
+ void setLevels(std::size_t minLevel, std::size_t maxLevel)
+ {
+ minLevel_ = minLevel;
+ maxLevel_ = maxLevel;
+ }
+
+ /*!
+ * \brief Function to set if refinement at Dirichlet boundaries is to be used.
+ *
+ */
+ void setRefinementAtDirichletBC(bool refine)
+ {
+ refineAtDirichletBC_ = refine;
+ }
+
+ /*!
+ * \brief Function to set if refinement at sources is to be used.
+ *
+ */
+ void setRefinementAtSources(bool refine)
+ {
+ refineAtSource_ = refine;
+ }
+
+ /*!
+ * \brief Function to set if refinement at sources is to be used.
+ *
+ */
+ void setRefinementAtFluxBoundaries(bool refine)
+ {
+ refineAtFluxBC_ = refine;
+ }
+
+ /*!
+ * \brief Calculates the indicator used for refinement/coarsening for each grid cell.
+ *
+ * \param sol The solution for which indicator is to be calculated
+ */
+ template<class SolutionVector>
+ void calculate(const SolutionVector& sol)
+ {
+ //! prepare an indicator for refinement
+ indicatorVector_.assign(fvGridGeometry_->gridView().size(0), false);
+
+ for (const auto& element : elements(fvGridGeometry_->gridView()))
+ {
+ const auto eIdx = fvGridGeometry_->elementMapper().index(element);
+
+ //! refine any element being below the minimum level
+ if (element.level() < minLevel_)
+ {
+ indicatorVector_[eIdx] = true;
+ continue; //! proceed to the next element
+ }
+
+ //! If refinement at sources/BCs etc is deactivated, skip the rest
+ if (!refineAtSource_ && !refineAtFluxBC_ && !refineAtDirichletBC_)
+ continue;
+
+ //! if the element is already on the maximum permissive level, skip rest
+ if (element.level() == maxLevel_)
+ continue;
+
+ //! get the fvGeometry and elementVolVars needed for the bc and source interfaces
+ auto fvGeometry = localView(*fvGridGeometry_);
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables_->curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, sol);
+
+ //! Check if we have to refine around a source term
+ if (refineAtSource_)
+ {
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ auto source = problem_->source(element, fvGeometry, elemVolVars, scv);
+ if (source.infinity_norm() > eps_)
+ {
+ indicatorVector_[eIdx] = true;
+ break; //! element is marked, escape scv loop
+ }
+ }
+ }
+
+ //! Check if we have to refine at the boundary
+ if (!indicatorVector_[eIdx] //! proceed if element is not already marked
+ && element.hasBoundaryIntersections() //! proceed if element is on boundary
+ && (refineAtDirichletBC_ || refineAtFluxBC_)) //! proceed if boundary refinement is active
+ {
+ //! cell-centered schemes
+ if (!isBox)
+ {
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ //! skip non-boundary scvfs
+ if (!scvf.boundary())
+ continue;
+
+ const auto bcTypes = problem_->boundaryTypes(element, scvf);
+ //! We assume pure BCs, mixed boundary types are not allowed anyway!
+ if(bcTypes.hasOnlyDirichlet() && refineAtDirichletBC_)
+ {
+ indicatorVector_[eIdx] = true;
+ break; //! element is marked, escape scvf loop
+ }
+
+ //! we are on a pure Neumann boundary
+ else if(refineAtFluxBC_)
+ {
+ const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, scvf);
+ if (fluxes.infinity_norm() > eps_)
+ {
+ indicatorVector_[eIdx] = true;
+ break; //! element is marked, escape scvf loop
+ }
+ }
+ }
+ }
+ //! box-scheme
+ else
+ {
+ //! container to store bcTypes
+ std::vector<BoundaryTypes> bcTypes(fvGeometry.numScv());
+
+ //! Get bcTypes and maybe mark for refinement on Dirichlet boundaries
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ bcTypes[scv.indexInElement()] = problem_->boundaryTypes(element, scv);
+ if (refineAtDirichletBC_ && bcTypes[scv.indexInElement()].hasDirichlet())
+ {
+ indicatorVector_[eIdx] = true;
+ break; //! element is marked, escape scv loop
+ }
+ }
+
+ //! If element hasn't been marked because of Dirichlet BCS, check Neumann BCs
+ if (!indicatorVector_[eIdx] && refineAtFluxBC_)
+ {
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ //! check if scvf is on Neumann boundary
+ if (scvf.boundary() && bcTypes[scvf.insideScvIdx()].hasNeumann())
+ {
+ const auto fluxes = problem_->neumann(element, fvGeometry, elemVolVars, scvf);
+ if (fluxes.infinity_norm() > eps_)
+ {
+ indicatorVector_[eIdx] = true;
+ break; //! element is marked, escape scvf loop
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /*! \brief function call operator to return mark
+ *
+ * \return 1 if an element should be refined
+ * -1 if an element should be coarsened
+ * 0 otherwise
+ *
+ * \note In this initialization indicator implementation
+ * element coarsening is not considered.
+ *
+ * \param element A grid element
+ */
+ int operator() (const Element& element) const
+ {
+ if (indicatorVector_[fvGridGeometry_->elementMapper().index(element)])
+ return 1;
+ return 0;
+ }
+
+private:
+ std::shared_ptr<const Problem> problem_; //! The problem to be solved
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry_; //! The finite volume grid geometry
+ std::shared_ptr<const GridVariables> gridVariables_; //! The secondary variables on the grid
+ std::vector<bool> indicatorVector_; //! Indicator for BCs/sources
+
+ int minLevel_; //! The minimum allowed level
+ int maxLevel_; //! The maximum allowed level
+ bool refineAtDirichletBC_; //! Specifies if it should be refined at Dirichlet BCs
+ bool refineAtFluxBC_; //! Specifies if it should be refined at non-zero Neumann BCs
+ bool refineAtSource_; //! Specifies if it should be refined at sources
+ Scalar eps_; //! Threshold for refinement at sources/BCS
+};
+
+}
+#endif
diff --git a/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh b/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
index 23e970d71d8aaffbc706fe782ab6fda4f64bbd29..bb2f8f945f83b19dce51385fd83b51166d810235 100644
--- a/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
+++ b/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
@@ -65,24 +65,27 @@ public:
, maxLevel_(getParamFromGroup<std::size_t>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.MaxLevel", 0))
{}
- /*! \brief Function to set the minumum allowed levels.
- *.
+ /*!
+ * \brief Function to set the minimum allowed level.
+ *
*/
void setMinLevel(std::size_t minLevel)
{
minLevel_ = minLevel;
}
- /*! \brief Function to set the maximum allowed levels.
- *.
+ /*!
+ * \brief Function to set the maximum allowed level.
+ *
*/
void setMaxLevel(std::size_t maxLevel)
{
maxLevel_ = maxLevel;
}
- /*! \brief Function to set the minumum/maximum allowed levels.
- *.
+ /*!
+ * \brief Function to set the minumum/maximum allowed levels.
+ *
*/
void setLevels(std::size_t minLevel, std::size_t maxLevel)
{
@@ -90,7 +93,8 @@ public:
maxLevel_ = maxLevel;
}
- /*! \brief Calculates the indicator used for refinement/coarsening for each grid cell.
+ /*!
+ * \brief Calculates the indicator used for refinement/coarsening for each grid cell.
*
* This standard two-phase indicator is based on the saturation gradient.
*/
diff --git a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
index 5e9aa92d80ce192ec095d283290814e0a977df87..296ed5b9e4448153454202e517be99e11879c8e2 100644
--- a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
+++ b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
@@ -16,8 +16,7 @@ LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
Name = 2padaptive # name passed to the output routines
[Adaptive]
-EnableInitializationIndicator = 1
-RefineAtDirichletBC = 1
+RefineAtDirichletBC = 0
RefineAtFluxBC = 1
MinLevel = 0
MaxLevel = 2
diff --git a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
index 8c62d57878cc12829ca6a5e0f458690b5e7cc795..e6a0bd4bd735a29c3e2b04f25a550008235b02bd 100644
--- a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
+++ b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
@@ -55,6 +55,7 @@
#include <dumux/adaptive/adapt.hh>
#include <dumux/adaptive/markelements.hh>
+#include <dumux/adaptive/initializationindicator.hh>
#include <dumux/porousmediumflow/2p/implicit/griddatatransfer.hh>
#include <dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh>
@@ -129,6 +130,41 @@ int main(int argc, char** argv) try
TwoPGridAdaptIndicator<TypeTag> indicator(fvGridGeometry);
TwoPGridDataTransfer<TypeTag> dataTransfer(problem, fvGridGeometry, gridVariables, x);
+ // Do initial refinement around sources/BCs
+ GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, fvGridGeometry, gridVariables);
+
+ // refine up to the maximum level
+ const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
+ for (std::size_t i = 0; i < maxLevel; ++i)
+ {
+ initIndicator.calculate(x);
+
+ bool wasAdapted = false;
+ if (markElements(GridCreator::grid(), initIndicator))
+ wasAdapted = adapt(GridCreator::grid(), dataTransfer);
+
+ // update grid data after adaption
+ if (wasAdapted)
+ {
+ xOld = x; //! Overwrite the old solution with the new (resized & interpolated) one
+ gridVariables->init(x, xOld); //! Initialize the secondary variables to the new (and "new old") solution
+ }
+ }
+
+ // Do refinement for the initial conditions using our indicator
+ indicator.calculate(x, refineTol, coarsenTol);
+
+ bool wasAdapted = false;
+ if (markElements(GridCreator::grid(), indicator))
+ wasAdapted = adapt(GridCreator::grid(), dataTransfer);
+
+ // update grid data after adaption
+ if (wasAdapted)
+ {
+ xOld = x; //! Overwrite the old solution with the new (resized & interpolated) one
+ gridVariables->init(x, xOld); //! Initialize the secondary variables to the new (and "new old") solution
+ }
+
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");