From 9d50eef20c139e00d6adf51fba43a44e51801b8b Mon Sep 17 00:00:00 2001
From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>
Date: Tue, 27 Sep 2016 15:29:34 +0200
Subject: [PATCH] [TwoPAdaptive] Update of lens test problems
The tests have been adapted such that the total mass is calculated before
and after the adaption prodecure. Furthermore, the material parameters
are set element-wise.
---
.../2p/implicit/lensproblem.hh | 90 +++++++++++++++----
.../2p/implicit/lensspatialparams.hh | 6 +-
.../2p/implicit/test_boxadaptive2p.input | 6 +-
3 files changed, 80 insertions(+), 22 deletions(-)
diff --git a/test/porousmediumflow/2p/implicit/lensproblem.hh b/test/porousmediumflow/2p/implicit/lensproblem.hh
index a17e557d7b..3783e4c1a4 100644
--- a/test/porousmediumflow/2p/implicit/lensproblem.hh
+++ b/test/porousmediumflow/2p/implicit/lensproblem.hh
@@ -62,12 +62,12 @@ SET_TYPE_PROP(LensCCProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(LensBoxProblem, Grid, Dune::YaspGrid<2>);
#endif
-SET_TYPE_PROP(LensBoxAdaptiveProblem, Grid, Dune::UGGrid<2>);
-
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(LensCCAdaptiveProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
#endif
+SET_TYPE_PROP(LensBoxAdaptiveProblem, Grid, Dune::UGGrid<2>);
+
// Set the problem property
SET_TYPE_PROP(LensProblem, Problem, LensProblem<TypeTag>);
@@ -94,18 +94,19 @@ SET_TYPE_PROP(LensCCProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag> );
SET_TYPE_PROP(LensBoxProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag> );
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(LensCCAdaptiveProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag> );
-SET_TYPE_PROP(LensBoxAdaptiveProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag> );
SET_BOOL_PROP(LensCCAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionHelper, TwoPAdaptionHelper<TypeTag>);
+#endif
+
+SET_TYPE_PROP(LensBoxAdaptiveProblem, LinearSolver, ILU0BiCGSTABBackend<TypeTag> );
SET_BOOL_PROP(LensBoxAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
-//SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
+SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionHelper, TwoPAdaptionHelper<TypeTag>);
-#endif
NEW_PROP_TAG(BaseProblem);
SET_TYPE_PROP(LensBoxProblem, BaseProblem, ImplicitPorousMediaProblem<TypeTag>);
@@ -179,11 +180,18 @@ class LensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
// world dimension
- dimWorld = GridView::dimensionworld
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
};
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+ enum { adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid) };
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -225,19 +233,28 @@ public:
}
/*!
- * \brief User defined output after the time integration
- *
- * Will be called diretly after the time integration.
+ * \brief User defined output before the time integration
*/
- void postTimeStep()
+ void preTimeStep()
{
- // Calculate storage terms
- PrimaryVariables storage;
- this->model().globalStorage(storage);
+ if(adaptiveGrid)
+ {
+ PrimaryVariables totalMass = calculateTotalMass();
+ std::cout << "Total mass before grid adaption: " << std::endl;
+ std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << std::endl;
+
+ ParentType::preTimeStep();
- // Write mass balance information for rank 0
- if (this->gridView().comm().rank() == 0) {
- std::cout<<"Storage: " << storage << std::endl;
+ totalMass = calculateTotalMass();
+ std::cout << "Total mass after grid adaption: " << std::endl;
+ std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << std::endl;
+ }
+ else
+ {
+ ParentType::preTimeStep();
+ PrimaryVariables totalMass = calculateTotalMass();
+ std::cout << "Total mass: " << std::endl;
+ std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << std::endl;
}
}
@@ -420,6 +437,47 @@ private:
return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
}
+ PrimaryVariables calculateTotalMass()
+ {
+ PrimaryVariables totalMass(0);
+
+ for (const auto& element : elements(this->gridView()))
+ {
+ if(element.partitionType() == Dune::InteriorEntity)
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(this->gridView(), element);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ // get index
+ int dofIdx = this->model().dofMapper().subIndex(element, scvIdx, dofCodim);
+
+ VolumeVariables volVars;
+ volVars.update(this->model().curSol()[dofIdx],
+ *this,
+ element,
+ fvGeometry,
+ scvIdx,
+ false);
+
+ Scalar volume = fvGeometry.subContVol[scvIdx].volume;
+
+ totalMass[nPhaseIdx] += volume*volVars.density(nPhaseIdx)
+ * volVars.porosity() * volVars.saturation(nPhaseIdx);
+ totalMass[wPhaseIdx] += volume*volVars.density(wPhaseIdx)
+ * volVars.porosity() * volVars.saturation(wPhaseIdx);
+ }
+
+ }
+ }
+
+ if (this->gridView().comm().size() > 1)
+ totalMass = this->gridView().comm().sum(totalMass);
+
+ return totalMass;
+ }
+
Scalar temperature_;
Scalar eps_;
std::string name_;
diff --git a/test/porousmediumflow/2p/implicit/lensspatialparams.hh b/test/porousmediumflow/2p/implicit/lensspatialparams.hh
index ff2d532850..8d4b53a53e 100644
--- a/test/porousmediumflow/2p/implicit/lensspatialparams.hh
+++ b/test/porousmediumflow/2p/implicit/lensspatialparams.hh
@@ -161,10 +161,10 @@ public:
const FVElementGeometry &fvGeometry,
int scvIdx) const
{
- const GlobalPosition& globalPos = fvGeometry.subContVol[scvIdx].global;
+ const GlobalPosition& globalPos = element.geometry().center();
-// if (isInLens_(globalPos))
-// return lensMaterialParams_;
+ if (isInLens_(globalPos))
+ return lensMaterialParams_;
return outerMaterialParams_;
}
diff --git a/test/porousmediumflow/2p/implicit/test_boxadaptive2p.input b/test/porousmediumflow/2p/implicit/test_boxadaptive2p.input
index 0fd0969406..ff185bdc7c 100644
--- a/test/porousmediumflow/2p/implicit/test_boxadaptive2p.input
+++ b/test/porousmediumflow/2p/implicit/test_boxadaptive2p.input
@@ -5,10 +5,10 @@ TEnd = 3000 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 6 4
-#Cells = 48 32
-Cells = 24 16
-Refinement = 2
+Cells = 48 32
+Refinement = 0
CellType = Cube
+ClosureType = Green
[SpatialParams]
LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
--
GitLab