diff --git a/dumux/implicit/adaptive/gridadapt.hh b/dumux/implicit/adaptive/gridadapt.hh
index 8dca1c52f50df12b5bb01a7c5baf9c1967435e5e..015369d71731d6ed1a5ecda031f4fb5942c6d83b 100644
--- a/dumux/implicit/adaptive/gridadapt.hh
+++ b/dumux/implicit/adaptive/gridadapt.hh
@@ -147,14 +147,6 @@ public:
/*!
* @brief Returns true if grid cells have been marked for adaption
*/
- bool wasAdapted()
- {
- int sumMarked = problem_.grid().comm().sum(marked_);
- int sumCoarsened = problem_.grid().comm().sum(coarsened_);
-
- return (sumMarked != 0 || sumCoarsened != 0);
- }
-
bool wasAdapted() const
{
int sumMarked = problem_.grid().comm().sum(marked_);
@@ -438,10 +430,8 @@ public:
{}
void adaptGrid()
{}
- bool wasAdapted()
- {
- return false;
- }
+ bool wasAdapted() const
+ { return false; }
void setLevels(int, int)
{}
void setTolerance(int, int)
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index 6124493f88251dd3a750c93ed34007a104ab9281..590a908bc579dfeebd04b67657bd0ab6218f2153 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -423,23 +423,32 @@ public:
*/
void updateBegin()
{
- if(GET_PROP_VALUE(TypeTag, AdaptiveGrid) && problem_().gridAdapt().wasAdapted())
+ if(GET_PROP_VALUE(TypeTag, AdaptiveGrid) && problem_().gridChanged())
{
- uPrev_ = uCur_;
- prevGlobalVolVars_ = curGlobalVolVars_;
-
+ // update boundary indices
updateBoundaryIndices_();
- int numDofs = asImp_().numDofs();
-
if (isBox)
- boxVolume_.resize(numDofs);
+ boxVolume_.resize(asImp_().numDofs());
- jacAsm_->init(problem_());
- }
+ // update the fv geometry
+ globalFvGeometryPtr_->update(problem_());
- }
+ // resize and update the volVars with the initial solution
+ curGlobalVolVars_.update(problem_(), curSol());
+ // resize the matrix and assembly map if necessary
+ localJacobian().init(problem_());
+ jacobianAssembler().init(problem_());
+
+ // update the flux variables caches
+ globalfluxVarsCache_.update(problem_());
+
+ // update the previous solution and volvars
+ uPrev_ = uCur_;
+ prevGlobalVolVars_ = curGlobalVolVars_;
+ }
+ }
/*!
* \brief Called by the update() method if it was
diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh
index 2389194950576e02d5bfd0f69d827caf41d34bc0..f40db5ef92e1bb6b578c9f3f11241979f1e13233 100644
--- a/dumux/implicit/problem.hh
+++ b/dumux/implicit/problem.hh
@@ -534,7 +534,7 @@ public:
this->gridAdapt().adaptGrid();
// if the grid changed recompute the source map and the bounding box tree
- if (this->gridAdapt().wasAdapted())
+ if (asImp_().gridChanged())
{
// update bounding box tree if it exists
if (boundingBoxTree_)
@@ -918,6 +918,17 @@ public:
return GridCreator::grid();
}
+ /*!
+ * \brief Returns whether the grid has changed
+ */
+ bool gridChanged() const
+ {
+ if (adaptiveGrid)
+ return asImp_().gridAdapt().wasAdapted();
+ else
+ return false;
+ }
+
/*!
* \brief Returns adaptivity model used for the problem.
*/
diff --git a/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanager.hh b/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanager.hh
index 2bd7df5c31b1a6e44f5d0502d44716444846f3cc..1a64a13cb626f878c79dcd152e68c548eae94cb2 100644
--- a/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanager.hh
+++ b/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanager.hh
@@ -131,7 +131,19 @@ public:
* initializing the subproblems / models
*/
void postInit()
- {}
+ {
+ glue_ = std::make_shared<CCMultiDimensionGlue<TypeTag>>(bulkProblem(), lowDimProblem());
+ asImp_().computeLowDimVolumeFractions();
+ }
+
+ /*!
+ * \brief Update after the grid has changed
+ */
+ void update()
+ {
+ asImp_().computePointSourceData();
+ asImp_().computeLowDimVolumeFractions();
+ }
/*!
* \brief The bulk coupling stencil, i.e. which low-dimensional dofs
@@ -344,6 +356,31 @@ public:
std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
}
+ //! Compute the low dim volume fraction in the bulk domain cells
+ void computeLowDimVolumeFractions()
+ {
+ // intersect the bounding box trees
+ glue_->build();
+
+ // compute the low dim volume fractions
+ for (const auto& is : intersections(*glue_))
+ {
+ // all inside elements are identical...
+ const auto& inside = is.inside(0);
+ const auto intersectionGeometry = is.geometry();
+ const unsigned int lowDimElementIdx = lowDimProblem().elementMapper().index(inside);
+
+ // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
+ const auto radius = lowDimProblem().spatialParams().radius(lowDimElementIdx);
+ for (unsigned int outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+ {
+ const auto& outside = is.outside(outsideIdx);
+ const unsigned int bulkElementIdx = bulkProblem().elementMapper().index(outside);
+ lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
+ }
+ }
+ }
+
/*!
* \brief Methods to be accessed by the subproblems
*/
@@ -597,6 +634,9 @@ private:
//! id generator for point sources
unsigned int idCounter_;
+
+ //! The glue object
+ std::shared_ptr<CCMultiDimensionGlue<TypeTag>> glue_;
};
} // end namespace Dumux
diff --git a/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh b/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh
index f48d229325282b38388430f7a7b569a1bf269d4f..611578d116a56914082be9e51d6ce15d4852d7a8 100644
--- a/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh
+++ b/dumux/mixeddimension/embedded/cellcentered/bboxtreecouplingmanagersimple.hh
@@ -123,6 +123,7 @@ public:
*/
void preInit()
{
+ glue_ = std::make_shared<CCMultiDimensionGlue<TypeTag>>(bulkProblem(), lowDimProblem());
asImp_().computePointSourceData();
}
@@ -131,7 +132,18 @@ public:
* initializing the subproblems / models
*/
void postInit()
- {}
+ {
+ asImp_().computeLowDimVolumeFractions();
+ }
+
+ /*!
+ * \brief Update after the grid has changed
+ */
+ void update()
+ {
+ asImp_().computePointSourceData();
+ asImp_().computeLowDimVolumeFractions();
+ }
/*!
* \brief The bulk coupling stencil, i.e. which low-dimensional dofs
@@ -217,10 +229,9 @@ public:
lowDimVolumeInBulkElement_.resize(this->bulkGridView().size(0));
// intersect the bounding box trees
- Dumux::CCMultiDimensionGlue<TypeTag> glue(bulkProblem(), lowDimProblem());
- glue.build();
+ glue_->build();
- for (const auto& is : intersections(glue))
+ for (const auto& is : intersections(*glue_))
{
// all inside elements are identical...
const auto& inside = is.inside(0);
@@ -231,15 +242,6 @@ public:
const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
const unsigned int lowDimElementIdx = lowDimProblem().elementMapper().index(inside);
- // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
- const auto radius = lowDimProblem().spatialParams().radius(lowDimElementIdx);
- for (unsigned int outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
- {
- const auto& outside = is.outside(outsideIdx);
- const unsigned int bulkElementIdx = bulkProblem().elementMapper().index(outside);
- lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
- }
-
// iterate over all quadrature points
for (auto&& qp : quad)
{
@@ -292,6 +294,27 @@ public:
std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
}
+ //! Compute the low dim volume fraction in the bulk domain cells
+ void computeLowDimVolumeFractions()
+ {
+ for (const auto& is : intersections(*glue_))
+ {
+ // all inside elements are identical...
+ const auto& inside = is.inside(0);
+ const auto intersectionGeometry = is.geometry();
+ const unsigned int lowDimElementIdx = lowDimProblem().elementMapper().index(inside);
+
+ // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
+ const auto radius = lowDimProblem().spatialParams().radius(lowDimElementIdx);
+ for (unsigned int outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+ {
+ const auto& outside = is.outside(outsideIdx);
+ const unsigned int bulkElementIdx = bulkProblem().elementMapper().index(outside);
+ lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
+ }
+ }
+ }
+
/*!
* \brief Methods to be accessed by the subproblems
*/
@@ -493,6 +516,9 @@ private:
//! id generator for point sources
unsigned int idCounter_;
+
+ //! The glue object
+ std::shared_ptr<CCMultiDimensionGlue<TypeTag>> glue_;
};
} // end namespace Dumux
diff --git a/dumux/mixeddimension/model.hh b/dumux/mixeddimension/model.hh
index a5e41326c2947cd1a4c085e758d605758d95b007..1027f727aef46220b67471d3374080043e9a8b79 100644
--- a/dumux/mixeddimension/model.hh
+++ b/dumux/mixeddimension/model.hh
@@ -108,6 +108,7 @@ public:
*
* \param problem The object representing the problem which needs to
* be simulated.
+ * \note at this point the sub problems are already initialized
*/
void init(Problem &problem)
{
@@ -122,6 +123,8 @@ public:
uPrev_[lowDimIdx].resize(lowDimNumDofs);
// initialize local jocabian and jacobian assembler
+ // \note these local jacobians are not the localjacobian objects of the submodels
+ // \todo generalize so that there is only one unique local jac object
bulkLocalJacobian_.init(problem_());
lowDimLocalJacobian_.init(problem_());
jacAsm_ = std::make_shared<JacobianAssembler>();
@@ -363,17 +366,17 @@ public:
*/
void updateBegin()
{
- if((GET_PROP_VALUE(BulkProblemTypeTag, AdaptiveGrid) && problem_().bulkProblem().gridAdapt().wasAdapted()) ||
- (GET_PROP_VALUE(LowDimProblemTypeTag, AdaptiveGrid) && problem_().lowDimProblem().gridAdapt().wasAdapted()))
+ problem_().bulkProblem().model().updateBegin();
+ problem_().lowDimProblem().model().updateBegin();
+
+ if (problem_().bulkProblem().gridChanged() || problem_().lowDimProblem().gridChanged())
{
+ // resize matrix adjust prev sol
uPrev_ = uCur_;
+ //! Recompute the assembly map in the localjacobian
+ bulkLocalJacobian_.init(problem_());
+ lowDimLocalJacobian_.init(problem_());
jacAsm_->init(problem_());
-
- if (GET_PROP_VALUE(BulkProblemTypeTag, AdaptiveGrid) && problem_().bulkProblem().gridAdapt().wasAdapted())
- problem_().bulkProblem().model().updateBegin();
-
- if (GET_PROP_VALUE(LowDimProblemTypeTag, AdaptiveGrid) && problem_().lowDimProblem().gridAdapt().wasAdapted())
- problem_().lowDimProblem().model().updateBegin();
}
}
@@ -473,7 +476,6 @@ public:
*/
void resetJacobianAssembler ()
{
- jacAsm_.template reset<JacobianAssembler>(0);
jacAsm_ = std::make_shared<JacobianAssembler>();
jacAsm_->init(problem_());
}
diff --git a/dumux/porousmediumflow/richards/implicit/model.hh b/dumux/porousmediumflow/richards/implicit/model.hh
index d71307fc9b43289df4919d2bf9130f07390569bf..9e6c70d3a8636f60c049de0e24095d881359ad86 100644
--- a/dumux/porousmediumflow/richards/implicit/model.hh
+++ b/dumux/porousmediumflow/richards/implicit/model.hh
@@ -139,7 +139,7 @@ public:
vtkOutputModule.addSecondaryVariable("mobility", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
vtkOutputModule.addSecondaryVariable("kr", [](const VolumeVariables& v){ return v.relativePermeability(wPhaseIdx); });
vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
- if(GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity));
+ if(GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
vtkOutputModule.addSecondaryVariable("pressure head", [](const VolumeVariables& v){ return v.pressureHead(wPhaseIdx); });
vtkOutputModule.addSecondaryVariable("water content", [](const VolumeVariables& v){ return v.waterContent(wPhaseIdx); });
diff --git a/dumux/porousmediumflow/richards/implicit/volumevariables.hh b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
index bf993368f1dc85369b13887a77a33488c0db03a6..82518cdb9a3cb0ec209b0966ea11b4c38a4f2ed1 100644
--- a/dumux/porousmediumflow/richards/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
@@ -170,7 +170,7 @@ public:
*
* \param phaseIdx The index of the fluid phase
*/
- Scalar saturation(const int phaseIdx) const
+ Scalar saturation(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? fluidState_.saturation(wPhaseIdx) : 1.0-fluidState_.saturation(wPhaseIdx); }
/*!
@@ -179,7 +179,7 @@ public:
*
* \param phaseIdx The index of the fluid phase
*/
- Scalar density(const int phaseIdx) const
+ Scalar density(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? fluidState_.density(phaseIdx) : 0.0; }
/*!
@@ -193,7 +193,7 @@ public:
*
* \param phaseIdx The index of the fluid phase
*/
- Scalar pressure(const int phaseIdx) const
+ Scalar pressure(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? fluidState_.pressure(phaseIdx) : pn_; }
/*!
@@ -207,7 +207,7 @@ public:
*
* \param phaseIdx The index of the fluid phase
*/
- Scalar mobility(const int phaseIdx) const
+ Scalar mobility(const int phaseIdx = wPhaseIdx) const
{ return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); }
/*!
@@ -217,7 +217,7 @@ public:
* \param phaseIdx The index of the fluid phase
* \note The non-wetting phase is infinitely mobile
*/
- Scalar viscosity(const int phaseIdx) const
+ Scalar viscosity(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? fluidState_.viscosity(wPhaseIdx) : 0.0; }
/*!
@@ -226,7 +226,7 @@ public:
*
* \param phaseIdx The index of the fluid phase
*/
- Scalar relativePermeability(const int phaseIdx) const
+ Scalar relativePermeability(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? relativePermeabilityWetting_ : 1.0; }
/*!
@@ -260,7 +260,7 @@ public:
* manually do a conversion. It is not correct if the density is not constant
* or the gravity different
*/
- Scalar pressureHead(const int phaseIdx) const
+ Scalar pressureHead(const int phaseIdx = wPhaseIdx) const
{ return 100.0 *(pressure(phaseIdx) - pn_)/density(phaseIdx)/9.81; }
/*!
@@ -274,7 +274,7 @@ public:
* \note this function is here as a convenience to the user to not have to
* manually do a conversion.
*/
- Scalar waterContent(const int phaseIdx) const
+ Scalar waterContent(const int phaseIdx = wPhaseIdx) const
{ return saturation(phaseIdx) * porosity_; }
protected:
diff --git a/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
index 77ad87ab84f7fe3ae828e3eac68b6a0abfbb92a4..e36986b5898e3df2bad70b5900a6498396fe01c2 100644
--- a/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
@@ -139,7 +139,7 @@ public:
*
* We always forward to the fluid state with the phaseIdx property (see class description).
*/
- Scalar molarDensity(const int phaseIdx) const
+ Scalar molarDensity(const int phaseIdx = wPhaseIdx) const
{ return phaseIdx == wPhaseIdx ? this->fluidState_.molarDensity(phaseIdx) : 0.0; }
/*!
diff --git a/test/porousmediumflow/2p/implicit/fractureproblem.hh b/test/porousmediumflow/2p/implicit/fractureproblem.hh
index 7c2baf6c507ca914ce134b51a978ffa75613825d..f60b6bea32a330f9134a53b86e427ec050e4699c 100644
--- a/test/porousmediumflow/2p/implicit/fractureproblem.hh
+++ b/test/porousmediumflow/2p/implicit/fractureproblem.hh
@@ -30,6 +30,7 @@
#include <dumux/porousmediumflow/2p/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include <dumux/implicit/cellcentered/propertydefaults.hh>
#include "fracturespatialparams.hh"
diff --git a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
index d21b87e1d0aba005f0611fcacf359f5130817943..2b94ffc26ff469142398931edd6cb02d0d644ebb 100644
--- a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
@@ -25,6 +25,7 @@
#define DUMUX_INJECTION_PROBLEM_HH
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
#include <dumux/porousmediumflow/2p2c/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>