From 234fe48aa952dee790ae55b9be89a8bcff1c7125 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 6 Dec 2017 18:19:05 +0100
Subject: [PATCH] [2p] introduce twop adaptive grid framework and test
---
.../2p/implicit/adaptionhelper.hh | 454 ------------------
.../2p/implicit/gridadaptindicator.hh | 319 ++++++------
.../2p/implicit/griddatatransfer.hh | 419 ++++++++++++++++
.../2p/implicit/CMakeLists.txt | 1 +
.../2p/implicit/adaptive/CMakeLists.txt | 36 ++
.../implicit/adaptive/test_2p_adaptive.input | 25 +
.../implicit/adaptive/test_2p_adaptive_fv.cc | 265 ++++++++++
7 files changed, 912 insertions(+), 607 deletions(-)
delete mode 100644 dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
create mode 100644 dumux/porousmediumflow/2p/implicit/griddatatransfer.hh
create mode 100644 test/porousmediumflow/2p/implicit/adaptive/CMakeLists.txt
create mode 100644 test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
create mode 100644 test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
diff --git a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh b/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
deleted file mode 100644
index 181ba0af5d..0000000000
--- a/dumux/porousmediumflow/2p/implicit/adaptionhelper.hh
+++ /dev/null
@@ -1,454 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-#ifndef DUMUX_TWOP_ADAPTIONHELPER_HH
-#define DUMUX_TWOP_ADAPTIONHELPER_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/implicit/adaptive/adaptionhelper.hh>
-
-namespace Dumux {
-
-/*!
- * \brief Base class holding the variables for implicit models.
- */
-template<class TypeTag>
-class TwoPAdaptionHelper : public ImplicitAdaptionHelper<TypeTag>
-{
-private:
- using ParentType = ImplicitAdaptionHelper<TypeTag>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-
- enum {
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
-
- //indices
- pressureIdx = Indices::pressureIdx,
- saturationIdx = Indices::saturationIdx,
- pwsn = Indices::pwsn,
- pnsw = Indices::pnsw,
- formulation = GET_PROP_VALUE(TypeTag, Formulation),
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx
- };
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
-
- using Grid = typename GridView::Grid;
- using Element = typename GridView::Traits::template Codim<0>::Entity;
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- using CoordScalar = typename GridView::ctype;
- using LocalPosition = Dune::FieldVector<CoordScalar, dim>;
-
- using LocalFiniteElementCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
- using LocalFiniteElement = typename LocalFiniteElementCache::FiniteElementType;
-
- struct AdaptedValues
- {
- ElementSolutionVector u;
- int count;
- PrimaryVariables associatedMass;
- AdaptedValues(): count(0), associatedMass(0.0) {}
- };
-
- using PersistentContainer = Dune::PersistentContainer<Grid, AdaptedValues>;
- PersistentContainer adaptionMap_;
-
-public:
- //! Constructs an adaption helper object
- /**
- * @param gridView a DUNE gridview object
- */
- TwoPAdaptionHelper(Problem& problem) : ParentType(problem), adaptionMap_(problem.grid(), 0)
- {
- if(FluidSystem::isCompressible(wPhaseIdx) || FluidSystem::isCompressible(nPhaseIdx))
- DUNE_THROW(Dune::InvalidStateException, "This adaption helper is only mass conservative for incompressible fluids!");
- }
-
- /*!
- * 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(Problem& problem)
- {
- adaptionMap_.resize();
-
- // loop over all levels of the grid
- for (int level = problem.grid().maxLevel(); level >= 0; level--)
- {
- // get grid view on level grid
- auto levelView = problem.grid().levelGridView(level);
-
- for (const auto& element : elements(levelView))
- {
- //get your map entry
- auto& adaptedValues = adaptionMap_[element];
-
- // put values in the map
- if (element.isLeaf())
- {
- auto fvGeometry = localView(problem.model().globalFvGeometry());
- fvGeometry.bindElement(element);
-
- auto elemVolVars = localView(problem.model().curGlobalVolVars());
- elemVolVars.bindElement(element, fvGeometry, problem.model().curSol());
-
- for (auto&& scv : scvs(fvGeometry))
- {
- adaptedValues.u = problem.model().elementSolution(element, problem.model().curSol());
-
- VolumeVariables volVars;
- volVars.update(adaptedValues.u, problem, element, scv);
-
- adaptedValues.associatedMass[nPhaseIdx] += scv.volume() * volVars.density(nPhaseIdx)
- * volVars.porosity() * volVars.saturation(nPhaseIdx);
- adaptedValues.associatedMass[wPhaseIdx] += scv.volume() * volVars.density(wPhaseIdx)
- * volVars.porosity() * volVars.saturation(wPhaseIdx);
- }
- adaptedValues.count = 1;
- }
- // Average in father
- if (element.level() > 0)
- {
- auto& adaptedValuesFather = adaptionMap_[element.father()];
- // For some grids the father element is identical to the son element.
- // For that case averaging is not necessary.
- if(&adaptedValues != &adaptedValuesFather)
- {
- adaptedValuesFather.count += 1;
- storeAdaptionValues(adaptedValues, adaptedValuesFather);
- }
- }
-
- if(isBox && !element.isLeaf())
- {
- adaptedValues.u = problem.model().elementSolution(element, problem.model().curSol());
- }
- }
- }
- }
-
- /*!
- * 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(Problem& problem)
- {
- adaptionMap_.resize();
- // vectors storing the mass associated with each vertex, when using the box method
- std::vector<Scalar> massCoeff;
- std::vector<Scalar> associatedMass;
-
- if(isBox)
- {
- massCoeff.resize(problem.model().numDofs(), 0.0);
- associatedMass.resize(problem.model().numDofs(), 0.0);
- }
-
- for (int level = 0; level <= problem.grid().maxLevel(); level++)
- {
- auto levelView = problem.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() || element.level() == 0)
- {
- // entry is in map, write in leaf
- if (element.isLeaf())
- {
- auto& adaptedValues = adaptionMap_[element];
-
- auto fvGeometry = localView(problem.model().globalFvGeometry());
- fvGeometry.bindElement(element);
-
- auto elementVolume = element.geometry().volume();
-
- for (auto&& scv : scvs(fvGeometry))
- {
- VolumeVariables volVars;
- volVars.update(adaptedValues.u, problem, element, scv);
-
- problem.model().curSol()[scv.dofIndex()] = getPriVars(adaptedValues, scv);
-
- if (formulation == pwsn)
- {
- problem.model().curSol()[scv.dofIndex()][saturationIdx] = adaptedValues.associatedMass[nPhaseIdx];
- problem.model().curSol()[scv.dofIndex()][saturationIdx] /= elementVolume * volVars.density(nPhaseIdx) * volVars.porosity();
- }
- else if (formulation == pnsw)
- {
- problem.model().curSol()[scv.dofIndex()][saturationIdx] = adaptedValues.associatedMass[wPhaseIdx];
- problem.model().curSol()[scv.dofIndex()][saturationIdx] /= elementVolume * volVars.density(wPhaseIdx) * volVars.porosity();
- }
-
- if(isBox)
- {
- if (int(formulation) == pwsn)
- {
- massCoeff[scv.dofIndex()] += scv.volume() * volVars.density(nPhaseIdx) * volVars.porosity();
- associatedMass[scv.dofIndex()] += scv.volume() / elementVolume * adaptedValues.associatedMass[nPhaseIdx];
- }
- else if (int(formulation) == pnsw)
- {
- massCoeff[scv.dofIndex()] += scv.volume() * volVars.density(wPhaseIdx) * volVars.porosity();
- associatedMass[scv.dofIndex()] += scv.volume() / elementVolume * adaptedValues.associatedMass[wPhaseIdx];
- }
- }
-
- }
-
- }
- }
- else
- {
- // value is not in map, interpolate from father element
- if (element.hasFather())
- {
- auto fatherElement = element.father();
- while(fatherElement.isNew() && fatherElement.level() > 0)
- fatherElement = fatherElement.father();
-
- Scalar massFather = 0.0;
-
- if(!isBox)
- {
- auto& adaptedValuesFather = adaptionMap_[fatherElement];
-
- if (formulation == pwsn)
- massFather = adaptedValuesFather.associatedMass[nPhaseIdx];
-
- else if (formulation == pnsw)
- massFather = adaptedValuesFather.associatedMass[wPhaseIdx];
-
- // access new son
- auto& adaptedValues = adaptionMap_[element];
- adaptedValues.count = 1;
-
- auto fvGeometry = localView(problem.model().globalFvGeometry());
- fvGeometry.bindElement(element);
-
- adaptedValues.u = adaptedValuesFather.u;
-
- for (auto&& scv : scvs(fvGeometry))
- {
- VolumeVariables volVars;
- volVars.update(adaptedValues.u, problem, element, scv);
-
- Scalar massCoeffSon = 0.0;
- if (int(formulation) == pwsn)
- massCoeffSon = scv.volume() * volVars.density(nPhaseIdx) * volVars.porosity();
-
- else if (int(formulation) == pnsw)
- massCoeffSon = scv.volume() * volVars.density(wPhaseIdx) * volVars.porosity();
-
- auto fatherElementVolume = fatherElement.geometry().volume();
- adaptedValues.u[0][saturationIdx] = (scv.volume() / fatherElementVolume * massFather)/massCoeffSon;
- }
-
- // if we are on leaf, store reconstructed values of son in CellData object
- if (element.isLeaf())
- {
- // access new CellData object
- auto newIdxI = problem.elementMapper().index(element);
- problem.model().curSol()[newIdxI] = adaptedValues.u[0];
- }
- }
- else
- {
- // auto& adaptedValuesFather = adaptionMap_[fatherElement];
- // // access new son
- // auto& adaptedValues = adaptionMap_[element];
- // adaptedValues.u.clear();
- // adaptedValues.count = 1;
- //
- // const auto geometryI = element.geometry();
- //
- // FVElementGeometry fvGeometry;
- // fvGeometry.update(problem.gridView(), element);
- //
- // for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- // {
- // auto subEntity = element.template subEntity <dofCodim>(scvIdx);
- //
- // LocalPosition dofCenterPos = geometryI.local(subEntity.geometry().center());
- // const LocalFiniteElementCache feCache;
- // Dune::GeometryType geomType = fatherElement.geometry().type();
- //
- // // Interpolate values from father element by using ansatz functions
- // const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
- // std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
- // localFiniteElement.localBasis().evaluateFunction(dofCenterPos, shapeVal);
- // PrimaryVariables u(0);
- // for (int j = 0; j < shapeVal.size(); ++j)
- // {
- // u.axpy(shapeVal[j], adaptedValuesFather.u[j]);
- // }
- //
- // adaptedValues.u.push_back(u);
- // adaptedValues.count = 1;
- //
- // if (element.isLeaf())
- // {
- // VolumeVariables volVars;
- // volVars.update(adaptedValues.u[scvIdx],
- // problem,
- // element,
- // fvGeometry,
- // scvIdx,
- // false);
- //
- // Scalar volume = fvGeometry.subContVol[scvIdx].volume;
- // Scalar fatherElementVolume = fatherElement.geometry().volume();
- //
- // int dofIdxGlobal = this->dofIndex(problem, element, scvIdx);
- // if (int(formulation) == pwsn)
- // {
- // massCoeff[dofIdxGlobal] += volume * volVars.density(nPhaseIdx) * volVars.porosity();
- // associatedMass[dofIdxGlobal] += volume/fatherElementVolume*adaptedValuesFather.associatedMass[nPhaseIdx];
- // }
- // else if (int(formulation) == pnsw)
- // {
- // massCoeff[dofIdxGlobal] += volume * volVars.density(wPhaseIdx) * volVars.porosity();
- // associatedMass[dofIdxGlobal] += volume/fatherElementVolume*adaptedValuesFather.associatedMass[wPhaseIdx];
- // }
- //
- // this->setAdaptionValues(adaptedValues, problem.model().curSol()[dofIdxGlobal],scvIdx);
- // }
- //
- // }
- }
- }
- else
- {
- DUNE_THROW(Dune::InvalidStateException, "Element is new but has no father element!");
- }
-
- }
- }
-
- }
-
- if(isBox)
- {
- for(int dofIdxGlobal = 0; dofIdxGlobal < problem.model().numDofs(); dofIdxGlobal++)
- problem.model().curSol()[dofIdxGlobal][saturationIdx] = associatedMass[dofIdxGlobal] / massCoeff[dofIdxGlobal];
- }
-
- // reset entries in restrictionmap
- adaptionMap_.resize( typename PersistentContainer::Value() );
- adaptionMap_.shrinkToFit();
- adaptionMap_.fill( typename PersistentContainer::Value() );
-
-//#if HAVE_MPI
-// // communicate ghost data
-// typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
-// typedef typename SolutionTypes::ElementMapper ElementMapper;
-// typedef VectorExchange<ElementMapper, std::vector<CellData> > DataHandle;
-// DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
-// problem.gridView().template communicate<DataHandle>(dataHandle,
-// Dune::InteriorBorder_All_Interface,
-// Dune::ForwardCommunication);
-//#endif
- }
-
- //! Stores sons entries into father element for averaging
- /**
- * Sum up the adaptedValues (sons values) into father element. We store from leaf
- * upwards, so sons are stored first, then cells on the next leaf (=fathers)
- * can be averaged.
- *
- * \param adaptedValues Container for model-specific values to be adapted
- * \param adaptedValuesFather Values to be adapted of father cell
- */
- static void storeAdaptionValues(AdaptedValues& adaptedValues,
- AdaptedValues& adaptedValuesFather)
- {
- if(!isBox)
- {
- if(adaptedValuesFather.u.size() == 0)
- adaptedValuesFather.u.resize(1);
-
- adaptedValuesFather.u[0] += adaptedValues.u[0];
- adaptedValuesFather.u[0] /= adaptedValues.count;
- adaptedValuesFather.associatedMass += adaptedValues.associatedMass;
- }
- else
- {
- adaptedValuesFather.associatedMass += adaptedValues.associatedMass;
- }
- }
- //! Set adapted values in CellData
- /**
- * This methods stores reconstructed values into the cellData object, by
- * this setting a newly mapped solution to the storage container of the
- * sequential models.
- *
- * \param adaptedValues Container for model-specific values to be adapted
- * \param u The variables to be stored
- */
- static PrimaryVariables getPriVars(const AdaptedValues& adaptedValues, const SubControlVolume& scv)
- {
- if(!isBox)
- {
- auto uNew = adaptedValues.u[0];
- uNew /= adaptedValues.count;
- return uNew;
- }
- else
- {
- return adaptedValues.u[scv.index()];
- }
- }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh b/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
index 5d7ac8b703..23e970d71d 100644
--- a/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
+++ b/dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh
@@ -16,170 +16,161 @@
* 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_IMPLICIT_GRIDADAPTINDICATOR2P_HH
-#define DUMUX_IMPLICIT_GRIDADAPTINDICATOR2P_HH
+#ifndef DUMUX_TWOP_ADAPTION_INDICATOR_HH
+#define DUMUX_TWOP_ADAPTION_INDICATOR_HH
+
+#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-//#include <dumux/linear/vectorexchange.hh>
+#include <dumux/discretization/evalsolution.hh>
/**
- * @file
- * @brief Class defining a standard, saturation dependent indicator for grid adaptation
+ * \file
+ * \brief Class defining a standard, saturation dependent indicator for grid adaptation
*/
namespace Dumux
{
-/*!\ingroup IMPES
- * @brief Class defining a standard, saturation dependent indicator for grid adaptation
- *
- * \tparam TypeTag The problem TypeTag
+
+/*!\ingroup TwoPModel
+ * \brief Class defining a standard, saturation dependent indicator for grid adaptation
*/
template<class TypeTag>
-class TwoPImplicitGridAdaptIndicator
+class TwoPGridAdaptIndicator
{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum { saturationIdx = Indices::saturationIdx };
- enum
- {
- saturationIdx = Indices::saturationIdx,
- pressureIdx = Indices::pressureIdx
- };
- enum
+public:
+ /*! \brief The Constructor
+ *
+ * \param fvGridGeometry The finite volume grid geometry
+ *
+ * Note: refineBound_, coarsenBound_ & maxSaturationDelta_ are chosen
+ * in a way such that the indicator returns false for all elements
+ * before having been calculated.
+ */
+ TwoPGridAdaptIndicator(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : fvGridGeometry_(fvGridGeometry)
+ , refineBound_(std::numeric_limits<Scalar>::max())
+ , coarsenBound_(std::numeric_limits<Scalar>::lowest())
+ , maxSaturationDelta_(fvGridGeometry_->gridView().size(0), 0.0)
+ , minLevel_(getParamFromGroup<std::size_t>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.MinLevel", 0))
+ , maxLevel_(getParamFromGroup<std::size_t>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Adaptive.MaxLevel", 0))
+ {}
+
+ /*! \brief Function to set the minumum allowed levels.
+ *.
+ */
+ void setMinLevel(std::size_t minLevel)
{
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx
- };
-
- enum {
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
+ minLevel_ = minLevel;
+ }
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef typename GridView::ctype CoordScalar;
+ /*! \brief Function to set the maximum allowed levels.
+ *.
+ */
+ void setMaxLevel(std::size_t maxLevel)
+ {
+ maxLevel_ = maxLevel;
+ }
- typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
- typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
+ /*! \brief Function to set the minumum/maximum allowed levels.
+ *.
+ */
+ void setLevels(std::size_t minLevel, std::size_t maxLevel)
+ {
+ minLevel_ = minLevel;
+ maxLevel_ = maxLevel;
+ }
-public:
/*! \brief Calculates the indicator used for refinement/coarsening for each grid cell.
*
- * This standard indicator is based on the saturation gradient.
+ * This standard two-phase indicator is based on the saturation gradient.
*/
- void calculateIndicator()
+ void calculate(const SolutionVector& sol,
+ Scalar refineTol = 0.05,
+ Scalar coarsenTol = 0.001)
{
- // prepare an indicator for refinement
- if(indicatorVector_.size() != problem_.gridView().size(0))
- {
- indicatorVector_.resize(problem_.gridView().size(0));
- }
- indicatorVector_ = -1e100;
+ //! reset the indicator to a state that returns false for all elements
+ refineBound_ = std::numeric_limits<Scalar>::max();
+ coarsenBound_ = std::numeric_limits<Scalar>::lowest();
+ maxSaturationDelta_.assign(fvGridGeometry_->gridView().size(0), 0.0);
- Scalar globalMax = -1e100;
- Scalar globalMin = 1e100;
+ //! maxLevel_ must be higher than minLevel_ to allow for refinement
+ if (minLevel_ >= maxLevel_)
+ return;
- // 1) calculate Indicator -> min, maxvalues
- // loop over all leaf-elements
- for (const auto& element : elements(problem_.gridView()))
- {
- // calculate minimum and maximum saturation
- // index of the current leaf-elements
- int globalIdxI = problem_.elementMapper().index(element);
+ //! check for inadmissible tolerance combination
+ if (coarsenTol > refineTol)
+ DUNE_THROW(Dune::InvalidStateException, "Refine tolerance must be higher than coarsen tolerance");
- Scalar satI = 0.0;
+ //! variables to hold the max/mon saturation values on the leaf
+ Scalar globalMax = std::numeric_limits<Scalar>::lowest();
+ Scalar globalMin = std::numeric_limits<Scalar>::max();
- if(!isBox)
- satI = problem_.model().curSol()[globalIdxI][saturationIdx];
- else
- {
- const LocalFiniteElementCache feCache;
- const auto geometryI = element.geometry();
- Dune::GeometryType geomType = geometryI.type();
-
- GlobalPosition centerI = geometryI.local(geometryI.center());
- const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
- std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
- localFiniteElement.localBasis().evaluateFunction(centerI, shapeVal);
-
- for (int i = 0; i < shapeVal.size(); ++i)
- {
- int dofIdxGlobal = problem_.model().dofMapper().subIndex(element, i, dofCodim);
- satI += shapeVal[i]*problem_.model().curSol()[dofIdxGlobal][saturationIdx];
- }
- }
+ //! Calculate minimum and maximum saturation
+ for (const auto& element : elements(fvGridGeometry_->gridView()))
+ {
+ //! index of the current leaf-element
+ const auto globalIdxI = fvGridGeometry_->elementMapper().index(element);
+
+ //! obtain the saturation at the center of the element
+ const auto geometry = element.geometry();
+ const ElementSolution elemSol(element, sol, *fvGridGeometry_);
+ const Scalar satI = evalSolution(element, geometry, *fvGridGeometry_, elemSol, geometry.center())[saturationIdx];
+ //! maybe update the global minimum/maximum
using std::min;
using std::max;
globalMin = min(satI, globalMin);
globalMax = max(satI, globalMax);
- // calculate refinement indicator in all cells
- for (const auto& intersection : intersections(problem_.gridView(), element))
+ //! calculate maximum delta in saturation for this cell
+ for (const auto& intersection : intersections(fvGridGeometry_->gridView(), element))
{
- // Only consider internal intersections
+ //! Only consider internal intersections
if (intersection.neighbor())
{
- // Access neighbor
- auto outside = intersection.outside();
- int globalIdxJ = problem_.elementMapper().index(outside);
+ //! Access neighbor
+ const auto outside = intersection.outside();
+ const auto globalIdxJ = fvGridGeometry_->elementMapper().index(outside);
- // Visit intersection only once
+ //! Visit intersection only once
if (element.level() > outside.level() || (element.level() == outside.level() && globalIdxI < globalIdxJ))
{
- Scalar satJ = 0.0;
-
- if(!isBox)
- satJ = problem_.model().curSol()[globalIdxJ][saturationIdx];
- else
- {
- const LocalFiniteElementCache feCache;
- const auto geometryJ = outside.geometry();
- Dune::GeometryType geomType = geometryJ.type();
-
- GlobalPosition centerJ = geometryJ.local(geometryJ.center());
- const LocalFiniteElement &localFiniteElement = feCache.get(geomType);
- std::vector<Dune::FieldVector<Scalar, 1> > shapeVal;
- localFiniteElement.localBasis().evaluateFunction(centerJ, shapeVal);
-
- for (int i = 0; i < shapeVal.size(); ++i)
- {
- int dofIdxGlobal = problem_.model().dofMapper().subIndex(outside, i, dofCodim);
-
- satJ += shapeVal[i]*problem_.model().curSol()[dofIdxGlobal][saturationIdx];
- }
- }
-
-
+ //! obtain saturation in the neighbor
+ const auto outsideGeometry = outside.geometry();
+ const ElementSolution elemSolJ(outside, sol, *fvGridGeometry_);
+ const Scalar satJ = evalSolution(outside, outsideGeometry, *fvGridGeometry_, elemSolJ, outsideGeometry.center())[saturationIdx];
using std::abs;
Scalar localdelta = abs(satI - satJ);
- using std::max;
- indicatorVector_[globalIdxI][0] = max(indicatorVector_[globalIdxI][0], localdelta);
- indicatorVector_[globalIdxJ][0] = max(indicatorVector_[globalIdxJ][0], localdelta);
+ maxSaturationDelta_[globalIdxI] = max(maxSaturationDelta_[globalIdxI], localdelta);
+ maxSaturationDelta_[globalIdxJ] = max(maxSaturationDelta_[globalIdxJ], localdelta);
}
}
}
}
- Scalar globaldelta = globalMax - globalMin;
+ //! compute the maximum delta in saturation
+ const auto globalDelta = globalMax - globalMin;
- refineBound_ = refinetol_*globaldelta;
- coarsenBound_ = coarsentol_*globaldelta;
+ //! compute the refinement/coarsening bounds
+ refineBound_ = refineTol*globalDelta;
+ coarsenBound_ = coarsenTol*globalDelta;
+// TODO: fix adaptive simulations in parallel
//#if HAVE_MPI
// // communicate updated values
// typedef VectorExchange<ElementMapper, ScalarSolutionType> DataHandle;
-// DataHandle dataHandle(problem_.elementMapper(), indicatorVector_);
+// DataHandle dataHandle(problem_.elementMapper(), maxSaturationDelta_);
// problem_.gridView().template communicate<DataHandle>(dataHandle,
// Dune::InteriorBorder_All_Interface,
// Dune::ForwardCommunication);
@@ -189,62 +180,84 @@ public:
// coarsenBound_ = problem_.gridView().comm().max(coarsenBound_);
//
//#endif
+
+ //! check if neighbors have to be refined too
+ for (const auto& element : elements(fvGridGeometry_->gridView(), Dune::Partitions::interior))
+ if (this->operator()(element) > 0)
+ checkNeighborsRefine_(element);
}
- /*! \brief Indicator function for marking of grid cells for refinement
+ /*! \brief function call operator to return mark
*
- * Returns true if an element should be refined.
+ * \return 1 if an element should be refined
+ * -1 if an element should be coarsened
+ * 0 otherwise
*
* \param element A grid element
*/
- bool refine(const Element& element)
+ int operator() (const Element& element) const
{
- return (indicatorVector_[problem_.elementMapper().index(element)] > refineBound_);
+ if (element.hasFather()
+ && maxSaturationDelta_[fvGridGeometry_->elementMapper().index(element)] < coarsenBound_)
+ {
+ return -1;
+ }
+ else if (element.level() < maxLevel_
+ && maxSaturationDelta_[fvGridGeometry_->elementMapper().index(element)] > refineBound_)
+ {
+ return 1;
+ }
+ else
+ return 0;
}
- /*! \brief Indicator function for marking of grid cells for coarsening
+private:
+ /*!
+ * \brief Method ensuring the refinement ratio of 2:1
*
- * Returns true if an element should be coarsened.
+ * For any given element, a loop over the neighbors checks if 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 element A grid element
+ * \param element Element of interest that is to be refined
+ * \param level level of the refined element: it is at least 1
+ * \return true if everything was successful
*/
- bool coarsen(const Element& element)
+ bool checkNeighborsRefine_(const Element &element, std::size_t level = 1)
{
- return (indicatorVector_[problem_.elementMapper().index(element)] < coarsenBound_);
- }
+ for(const auto& intersection : intersections(fvGridGeometry_->gridView(), element))
+ {
+ if(!intersection.neighbor())
+ continue;
- /*! \brief Initializes the adaptation indicator class*/
- void init()
- {
- refineBound_ = 0.;
- coarsenBound_ = 0.;
- indicatorVector_.resize(problem_.gridView().size(0));
- };
+ // obtain outside element
+ const auto outside = intersection.outside();
- /*! @brief Constructs a GridAdaptIndicator 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
- */
- TwoPImplicitGridAdaptIndicator(Problem& problem):
- problem_(problem)
- {
- refinetol_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, GridAdapt, RefineTolerance);
- coarsentol_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, GridAdapt, CoarsenTolerance);
+ // only mark non-ghost elements
+ if (outside.partitionType() == Dune::GhostEntity)
+ continue;
+
+ if (outside.level() < maxLevel_ && outside.level() < element.level())
+ {
+ // ensure refinement for outside element
+ maxSaturationDelta_[fvGridGeometry_->elementMapper().index(outside)] = std::numeric_limits<Scalar>::max();
+ if(level < maxLevel_)
+ checkNeighborsRefine_(outside, ++level);
+ }
+ }
+
+ return true;
}
-protected:
- Problem& problem_;
- Scalar refinetol_;
- Scalar coarsentol_;
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
+
Scalar refineBound_;
Scalar coarsenBound_;
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > indicatorVector_;
+ std::vector< Scalar > maxSaturationDelta_;
+ std::size_t minLevel_;
+ std::size_t maxLevel_;
};
-}
-#endif
+} // end namespace Dumux
+
+#endif /* DUMUX_TWOP_ADAPTION_INDICATOR_HH */
diff --git a/dumux/porousmediumflow/2p/implicit/griddatatransfer.hh b/dumux/porousmediumflow/2p/implicit/griddatatransfer.hh
new file mode 100644
index 0000000000..0d253be92a
--- /dev/null
+++ b/dumux/porousmediumflow/2p/implicit/griddatatransfer.hh
@@ -0,0 +1,419 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_TWOP_GRIDDATA_TRANSFER_HH
+#define DUMUX_TWOP_GRIDDATA_TRANSFER_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/adaptive/griddatatransfer.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Class performing the transfer of data on a grid from before to after adaptation.
+ */
+template<class TypeTag>
+class TwoPGridDataTransfer : public GridDataTransfer
+{
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ struct AdaptedValues
+ {
+ ElementSolution u;
+ int count = 0;
+ PrimaryVariables associatedMass{0.0};
+ bool wasLeaf = false;
+ };
+
+ using PersistentContainer = Dune::PersistentContainer<Grid, AdaptedValues>;
+
+ static constexpr int dim = Grid::dimension;
+ static constexpr int dimWorld = Grid::dimensionworld;
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+
+ //! export some indices
+ enum {
+ // index of saturation in primary variables
+ saturationIdx = Indices::saturationIdx,
+
+ // phase indices
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ // formulations
+ pwsn = Indices::pwsn,
+ pnsw = Indices::pnsw,
+
+ // the formulation that is actually used
+ formulation = GET_PROP_VALUE(TypeTag, Formulation)
+ };
+
+ //! This won't work (mass conservative) for compressible fluids
+ static_assert(!FluidSystem::isCompressible(wPhaseIdx)
+ && !FluidSystem::isCompressible(nPhaseIdx),
+ "This adaption helper is only mass conservative for incompressible fluids!");
+
+ //! check if the used formulation is implemented here
+ static_assert(formulation == pwsn || formulation == pnsw, "Chosen formulation not known to the TwoPGridDataTransfer");
+
+public:
+ /*! \brief Constructor
+ *
+ * \param problem The DuMuX problem to be solved
+ * \param fvGridGeometry The finite volume grid geometry
+ * \param gridVariables The secondary variables on the grid
+ * \param sol The solution (primary variables) on the grid
+ */
+ TwoPGridDataTransfer(std::shared_ptr<const Problem> problem,
+ std::shared_ptr<FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<const GridVariables> gridVariables,
+ SolutionVector& sol)
+ : GridDataTransfer()
+ , problem_(problem)
+ , fvGridGeometry_(fvGridGeometry)
+ , gridVariables_(gridVariables)
+ , sol_(sol)
+ , adaptionMap_(fvGridGeometry->gridView().grid(), 0)
+ {}
+
+ /*! \brief Stores primary variables and additional data
+ *
+ * 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 a container object, before the grid is adapted. Father elements hold averaged
+ * information from the son cells for the case of the sons being coarsened.
+ */
+ void store() override
+ {
+ adaptionMap_.resize();
+
+ const auto& grid = fvGridGeometry_->gridView().grid();
+ for (auto level = grid.maxLevel(); level >= 0; level--)
+ {
+ for (const auto& element : elements(grid.levelGridView(level)))
+ {
+ //! get map entry
+ auto& adaptedValues = adaptionMap_[element];
+
+ //! put values in the map for leaf elements
+ if (element.isLeaf())
+ {
+ auto fvGeometry = localView(*fvGridGeometry_);
+ fvGeometry.bindElement(element);
+
+ //! store current element solution
+ adaptedValues.u = ElementSolution(element, sol_, *fvGridGeometry_);
+
+ //! compute mass in the scvs
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(adaptedValues.u, *problem_, element, scv);
+
+ const auto poreVolume = scv.volume()*volVars.porosity();
+ adaptedValues.associatedMass[nPhaseIdx] += poreVolume * volVars.density(nPhaseIdx) * volVars.saturation(nPhaseIdx);
+ adaptedValues.associatedMass[wPhaseIdx] += poreVolume * volVars.density(wPhaseIdx) * volVars.saturation(wPhaseIdx);
+ }
+
+ //! leaf elements always start with count = 1
+ adaptedValues.count = 1;
+ adaptedValues.wasLeaf = true;
+ }
+ //! Average in father elements
+ if (element.level() > 0)
+ {
+ auto& adaptedValuesFather = adaptionMap_[element.father()];
+ // For some grids the father element is identical to the son element.
+ // In that case averaging is not necessary.
+ if(&adaptedValues != &adaptedValuesFather)
+ storeAdaptionValues(adaptedValues, adaptedValuesFather);
+ }
+
+ //! The vertices of the non-leaf elements exist on the leaf as well
+ //! This element solution constructor uses the vertex mapper to obtain
+ //! the privars at the vertices, thus, this works for non-leaf elements!
+ if(isBox && !element.isLeaf())
+ adaptedValues.u = ElementSolution(element, sol_, *fvGridGeometry_);
+ }
+ }
+ }
+
+ /*! \brief 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 (AdaptedValues).
+ */
+ void reconstruct() override
+ {
+ //! resize stuff (grid might have changed)
+ adaptionMap_.resize();
+ fvGridGeometry_->update();
+ sol_.resize(fvGridGeometry_->numDofs());
+
+ //! vectors storing the mass associated with each vertex, when using the box method
+ std::vector<Scalar> massCoeff;
+ std::vector<Scalar> associatedMass;
+
+ if(isBox)
+ {
+ massCoeff.resize(fvGridGeometry_->numDofs(), 0.0);
+ associatedMass.resize(fvGridGeometry_->numDofs(), 0.0);
+ }
+
+ //! iterate over leaf and reconstruct the solution
+ for (const auto& element : elements(fvGridGeometry_->gridView().grid().leafGridView(), Dune::Partitions::interior))
+ {
+ if (!element.isNew())
+ {
+ const auto& adaptedValues = adaptionMap_[element];
+
+ auto fvGeometry = localView(*fvGridGeometry_);
+ fvGeometry.bindElement(element);
+
+ //! obtain element solution from map (divide by count!)
+ auto elemSol = adaptedValues.u;
+ elemSol /= adaptedValues.count;
+
+ const auto elementVolume = element.geometry().volume();
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(elemSol, *problem_, element, scv);
+
+ //! write solution at dof in current solution vector
+ sol_[scv.dofIndex()] = elemSol[scv.indexInElement()];
+
+ const auto dofIdxGlobal = scv.dofIndex();
+ //! For cc schemes, overwrite the saturation by a mass conservative one here
+ if (!isBox)
+ {
+ //! only recalculate the saturations if element hasn't been leaf before adaptation
+ if (!adaptedValues.wasLeaf)
+ {
+ if (formulation == pwsn)
+ {
+ sol_[dofIdxGlobal][saturationIdx] = adaptedValues.associatedMass[nPhaseIdx];
+ sol_[dofIdxGlobal][saturationIdx] /= elementVolume * volVars.density(nPhaseIdx) * volVars.porosity();
+ }
+ else if (formulation == pnsw)
+ {
+ sol_[dofIdxGlobal][saturationIdx] = adaptedValues.associatedMass[wPhaseIdx];
+ sol_[dofIdxGlobal][saturationIdx] /= elementVolume * volVars.density(wPhaseIdx) * volVars.porosity();
+ }
+ }
+ }
+
+ //! For the box scheme, add mass & mass coefficient to container (saturations are recalculated at the end)
+ else
+ {
+ const auto scvVolume = scv.volume();
+ if (formulation == pwsn)
+ {
+ massCoeff[dofIdxGlobal] += scvVolume * volVars.density(nPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdxGlobal] += scvVolume / elementVolume * adaptedValues.associatedMass[nPhaseIdx];
+ }
+ else if (formulation == pnsw)
+ {
+ massCoeff[dofIdxGlobal] += scvVolume * volVars.density(wPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdxGlobal] += scvVolume / elementVolume * adaptedValues.associatedMass[wPhaseIdx];
+ }
+ }
+ }
+ }
+ else
+ {
+ //! value is not in map, interpolate from father element
+ assert(element.hasFather() && "new element does not have a father element!");
+
+ //! find the ancestor element that existed on the old grid already
+ auto fatherElement = element.father();
+ while(fatherElement.isNew() && fatherElement.level() > 0)
+ fatherElement = fatherElement.father();
+
+ if(!isBox)
+ {
+ const auto& adaptedValuesFather = adaptionMap_[fatherElement];
+
+ //! obtain the mass contained in father
+ Scalar massFather = 0.0;
+ if (formulation == pwsn)
+ massFather = adaptedValuesFather.associatedMass[nPhaseIdx];
+ else if (formulation == pnsw)
+ massFather = adaptedValuesFather.associatedMass[wPhaseIdx];
+
+ // obtain the element solution through the father
+ auto elemSolSon = adaptedValuesFather.u;
+ elemSolSon /= adaptedValuesFather.count;
+
+ auto fvGeometry = localView(*fvGridGeometry_);
+ fvGeometry.bindElement(element);
+
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(elemSolSon, *problem_, element, scv);
+
+ //! store constructed values of son in the current solution
+ sol_[scv.dofIndex()] = elemSolSon[0];
+
+ //! overwrite the saturation by a mass conservative one here
+ Scalar massCoeffSon = 0.0;
+ if (formulation == pwsn)
+ massCoeffSon = scv.volume() * volVars.density(nPhaseIdx) * volVars.porosity();
+ else if (formulation == pnsw)
+ massCoeffSon = scv.volume() * volVars.density(wPhaseIdx) * volVars.porosity();
+ sol_[scv.dofIndex()][saturationIdx] = (scv.volume() / fatherElement.geometry().volume() * massFather)/massCoeffSon;
+ }
+ }
+ else
+ {
+ auto& adaptedValuesFather = adaptionMap_[fatherElement];
+
+ auto fvGeometry = localView(*fvGridGeometry_);
+ fvGeometry.bindElement(element);
+
+ //! interpolate solution in the father to the vertices of the new son
+ ElementSolution elemSolSon(element, sol_, *fvGridGeometry_);
+ const auto fatherGeometry = fatherElement.geometry();
+ for (const auto& scv : scvs(fvGeometry))
+ elemSolSon[scv.indexInElement()] = evalSolution(fatherElement,
+ fatherGeometry,
+ adaptedValuesFather.u,
+ scv.dofPosition());
+
+ //! compute mass & mass coeffients for the scvs (saturations are recalculated at the end)
+ const auto fatherElementVolume = fatherGeometry.volume();
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(elemSolSon, *problem_, element, scv);
+
+ const auto dofIdxGlobal = scv.dofIndex();
+ const auto scvVolume = scv.volume(); //std::cout << "ratio = " << scvVolume / fatherElementVolume << std::endl;
+ if (int(formulation) == pwsn)
+ {
+ massCoeff[dofIdxGlobal] += scvVolume * volVars.density(nPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdxGlobal] += scvVolume / fatherElementVolume * adaptedValuesFather.associatedMass[nPhaseIdx];
+ }
+ else if (int(formulation) == pnsw)
+ {
+ massCoeff[dofIdxGlobal] += scvVolume * volVars.density(wPhaseIdx) * volVars.porosity();
+ associatedMass[dofIdxGlobal] += scvVolume / fatherElementVolume * adaptedValuesFather.associatedMass[wPhaseIdx];
+ }
+
+ //! store constructed (pressure) values of son in the current solution (saturation comes later)
+ sol_[dofIdxGlobal] = elemSolSon[scv.indexInElement()];
+ }
+ }
+ }
+ }
+
+ if(isBox)
+ {
+ for(std::size_t dofIdxGlobal = 0; dofIdxGlobal < fvGridGeometry_->numDofs(); dofIdxGlobal++)
+ sol_[dofIdxGlobal][saturationIdx] = associatedMass[dofIdxGlobal] / massCoeff[dofIdxGlobal];
+ }
+
+ //! reset entries in adaptation map
+ adaptionMap_.resize( typename PersistentContainer::Value() );
+ adaptionMap_.shrinkToFit();
+ adaptionMap_.fill( typename PersistentContainer::Value() );
+
+//! TODO: fix adaptive simulations in parallel
+//#if HAVE_MPI
+// // communicate ghost data
+// typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+// typedef typename SolutionTypes::ElementMapper ElementMapper;
+// typedef VectorExchange<ElementMapper, std::vector<CellData> > DataHandle;
+// DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal());
+// problem.gridView().template communicate<DataHandle>(dataHandle,
+// Dune::InteriorBorder_All_Interface,
+// Dune::ForwardCommunication);
+//#endif
+ }
+
+ private:
+
+ /*! \brief Stores sons entries into father element for averaging
+ *
+ * Sum up the adaptedValues (sons values) into father element. We store from leaf
+ * upwards, so sons are stored first, then cells on the next leaf (=fathers)
+ * can be averaged.
+ *
+ * \param adaptedValues Container for model-specific values to be adapted
+ * \param adaptedValuesFather Values to be adapted of father cell
+ */
+ static void storeAdaptionValues(AdaptedValues& adaptedValues,
+ AdaptedValues& adaptedValuesFather)
+ {
+ //! Add associated mass of the child to the one of the father
+ adaptedValuesFather.associatedMass += adaptedValues.associatedMass;
+
+ if(!isBox)
+ {
+ //! add the child's primary variables to the ones of father
+ //! we have to divide the child's ones in case it was composed
+ //! of several children as well!
+ auto values = adaptedValues.u[0];
+ values /= adaptedValues.count;
+ adaptedValuesFather.u[0] += values;
+
+ //! keep track of the number of children that composed this father
+ adaptedValuesFather.count += 1;
+
+ //! A father element is never leaf
+ adaptedValuesFather.wasLeaf = false;
+ }
+ else
+ {
+ //! For the box scheme, scaling of primary variables by count is obsolete
+ //! Thus, we always want count = 1
+ adaptedValuesFather.count = 1;
+
+ //! A father element is never leaf
+ adaptedValuesFather.wasLeaf = false;
+ }
+ }
+
+ std::shared_ptr<const Problem> problem_;
+ std::shared_ptr<FVGridGeometry> fvGridGeometry_;
+ std::shared_ptr<const GridVariables> gridVariables_;
+ SolutionVector& sol_;
+ PersistentContainer adaptionMap_;
+};
+
+} // end namespace Dumux
+
+#endif /* DUMUX_TWOP_GRIDDATA_TRANSFER_HH */
diff --git a/test/porousmediumflow/2p/implicit/CMakeLists.txt b/test/porousmediumflow/2p/implicit/CMakeLists.txt
index 1bb24dedd4..aafe2e5b07 100644
--- a/test/porousmediumflow/2p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2p/implicit/CMakeLists.txt
@@ -1,4 +1,5 @@
#add_subdirectory(pointsources)
+add_subdirectory(adaptive)
add_subdirectory(fracture)
add_subdirectory(incompressible)
add_subdirectory(nonisothermal)
diff --git a/test/porousmediumflow/2p/implicit/adaptive/CMakeLists.txt b/test/porousmediumflow/2p/implicit/adaptive/CMakeLists.txt
new file mode 100644
index 0000000000..2fc2d6438d
--- /dev/null
+++ b/test/porousmediumflow/2p/implicit/adaptive/CMakeLists.txt
@@ -0,0 +1,36 @@
+dune_symlink_to_source_files(FILES "test_2p_adaptive.input")
+
+# using tpfa
+dune_add_test(NAME test_2p_adaptive_tpfa
+ SOURCES test_2p_adaptive_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleAdaptiveTpfa
+ CMAKE_GUARD dune-alugrid_FOUND
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/lensccadaptive-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/2p_adaptive_tpfa-00008.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_adaptive_tpfa test_2p_adaptive.input -Problem.Name 2p_adaptive_tpfa")
+
+# using mpfa
+dune_add_test(NAME test_2p_adaptive_mpfa
+ SOURCES test_2p_adaptive_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleAdaptiveMpfa
+ CMAKE_GUARD dune-alugrid_FOUND
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/lensccadaptive-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/2p_adaptive_mpfa-00008.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_adaptive_mpfa test_2p_adaptive.input -Problem.Name 2p_adaptive_mpfa")
+
+# using box
+dune_add_test(NAME test_2p_adaptive_box
+ SOURCES test_2p_adaptive_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleAdaptiveBox
+ CMAKE_GUARD dune-uggrid_FOUND
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/lensboxadaptive-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/2p_adaptive_box-00008.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_adaptive_box test_2p_adaptive.input -Problem.Name 2p_adaptive_box")
+
+set(CMAKE_BUILD_TYPE Release)
diff --git a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
new file mode 100644
index 0000000000..5e9aa92d80
--- /dev/null
+++ b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive.input
@@ -0,0 +1,25 @@
+[TimeLoop]
+DtInitial = 250 # [s]
+TEnd = 3000 # [s]
+
+[Grid]
+UpperRight = 6 4
+Cells = 48 32
+CellType = Cube
+ClosureType = Green
+
+[SpatialParams]
+LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
+LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
+
+[Problem]
+Name = 2padaptive # name passed to the output routines
+
+[Adaptive]
+EnableInitializationIndicator = 1
+RefineAtDirichletBC = 1
+RefineAtFluxBC = 1
+MinLevel = 0
+MaxLevel = 2
+CoarsenTolerance = 1e-4
+RefineTolerance = 1e-4
diff --git a/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
new file mode 100644
index 0000000000..8c62d57878
--- /dev/null
+++ b/test/porousmediumflow/2p/implicit/adaptive/test_2p_adaptive_fv.cc
@@ -0,0 +1,265 @@
+// -*- 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 test for the two-phase porousmedium flow model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+//! Use the incompressible problem for this adaptive test
+#include <test/porousmediumflow/2p/implicit/incompressible/problem.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+#include <dumux/nonlinear/newtonconvergencewriter.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+#include <dumux/adaptive/adapt.hh>
+#include <dumux/adaptive/markelements.hh>
+#include <dumux/porousmediumflow/2p/implicit/griddatatransfer.hh>
+#include <dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh>
+
+//! type tags for the adaptive versions of the two-phase incompressible problem
+namespace Dumux {
+ namespace Properties {
+ //! Type Tags for the adaptive tests
+ NEW_TYPE_TAG(TwoPIncompressibleAdaptiveTpfa, INHERITS_FROM(TwoPIncompressibleTpfa));
+ NEW_TYPE_TAG(TwoPIncompressibleAdaptiveMpfa, INHERITS_FROM(TwoPIncompressibleMpfa));
+ NEW_TYPE_TAG(TwoPIncompressibleAdaptiveBox, INHERITS_FROM(TwoPIncompressibleBox));
+
+ //! Use non-conforming refinement in the cell-centered tests, conforming for box
+ SET_TYPE_PROP(TwoPIncompressibleAdaptiveTpfa, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+ SET_TYPE_PROP(TwoPIncompressibleAdaptiveMpfa, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
+ SET_TYPE_PROP(TwoPIncompressibleAdaptiveBox, Grid, Dune::UGGrid<2>);
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TYPETAG);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // instantiate indicator & data transfer, read parameters for indicator
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
+ const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
+ TwoPGridAdaptIndicator<TypeTag> indicator(fvGridGeometry);
+ TwoPGridDataTransfer<TypeTag> dataTransfer(problem, fvGridGeometry, gridVariables, x);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (haveParam("Restart") || haveParam("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonController = Dumux::NewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // compute refinement indicator
+ indicator.calculate(x, refineTol, coarsenTol);
+
+ // mark elements and maybe adapt grid
+ bool wasAdapted = false;
+ if (markElements(GridCreator::grid(), indicator))
+ wasAdapted = adapt(GridCreator::grid(), dataTransfer);
+
+ if (wasAdapted)
+ {
+ //! Note that if we were using point sources, we would have to update the map here as well
+ xOld = x; //! Overwrite the old solution with the new (resized & interpolated) one
+ assembler->setJacobianPattern(); //! Tell the assembler to resize the matrix and set pattern
+ assembler->setResidualSize(); //! Tell the assembler to resize the residual
+ gridVariables->init(x, xOld); //! Initialize the secondary variables to the new (and "new old") solution
+ }
+
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
--
GitLab