diff --git a/test/porousmediumflow/2p/implicit/lensproblem.hh b/test/porousmediumflow/2p/implicit/lensproblem.hh index a17e557d7bf745db9e73fd8e11492e10fa0187f9..3783e4c1a4d1736d53ea49b91c5789c52016dc3c 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 ff2d532850bbfe72b565b0257529512b31659b00..8d4b53a53e2081fd5de6fe5f169b6eb9c56b5da1 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 0fd096940640e0db0c2d72396b04b0f52423d95c..ff185bdc7ca3c38f1df59a5a962291f21ba25e57 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