diff --git a/dumux/decoupled/2p/cellData2p_adaptive.hh b/dumux/decoupled/2p/cellData2p_adaptive.hh
index dcb10838a9ab191fec48852decb812de9618969f..f1c342cd402493dfe85ab99b45d7677f121b8a73 100644
--- a/dumux/decoupled/2p/cellData2p_adaptive.hh
+++ b/dumux/decoupled/2p/cellData2p_adaptive.hh
@@ -148,11 +148,11 @@ public:
static void setAdaptionValues(AdaptedValues& adaptedValues, AdaptedValues& adaptedValuesFather, const Problem& problem)
{
- adaptedValues.saturationW += adaptedValuesFather.saturationW / adaptedValuesFather.count;
- adaptedValues.saturationNW += adaptedValuesFather.saturationNW / adaptedValuesFather.count;
- adaptedValues.pressW += adaptedValuesFather.pressW / adaptedValuesFather.count;
- adaptedValues.pressNW += adaptedValuesFather.pressNW / adaptedValuesFather.count;
- adaptedValues.volCorr += adaptedValuesFather.volCorr / adaptedValuesFather.count;
+ adaptedValues.saturationW = adaptedValuesFather.saturationW / adaptedValuesFather.count;
+ adaptedValues.saturationNW = adaptedValuesFather.saturationNW / adaptedValuesFather.count;
+ adaptedValues.pressW = adaptedValuesFather.pressW / adaptedValuesFather.count;
+ adaptedValues.pressNW = adaptedValuesFather.pressNW / adaptedValuesFather.count;
+ adaptedValues.volCorr = adaptedValuesFather.volCorr / adaptedValuesFather.count;
}
};
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index 7d645e3be2df54f064eca60171a2dbfb77bc4d0f..26cfa08a42c6b6f2cca6f768818e542c3b58bb0e 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -901,8 +901,6 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
if (compressibility_)
{
FluidState& fluidState = cellData.fluidState();
- fluidState.setSaturation(wPhaseIdx, satW);
- fluidState.setSaturation(nPhaseIdx, satNW);
fluidState.setTemperature(temperature);
fluidState.setPressure(wPhaseIdx, pressW);
@@ -924,8 +922,6 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
}
else
{
- cellData.setSaturation(wPhaseIdx, satW);
- cellData.setSaturation(nPhaseIdx, satNW);
cellData.setCapillaryPressure(pc);
if (pressureType_ != pglobal)
@@ -938,8 +934,6 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
// initialize mobilities
Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[wPhaseIdx];
Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosity_[nPhaseIdx];
-// std::cout<<"MobilityW: "<<mobilityW <<"\n"
-// "MobilityNW"<< mobilityNW<<"\n";
if (compressibility_)
{
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
index 1c4a6d3e4d46bcf36b9b51c8a5c13f2a7e985bf3..13c211fad30cc61047a3cc964ff276715c76aed1 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
@@ -27,6 +27,7 @@
// dumux environment
#include <dumux/decoupled/2p/2pproperties.hh>
#include <dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh>
+#include <dumux/decoupled/common/fv/fvvelocity.hh>
/**
* @file
@@ -126,15 +127,25 @@ public:
ParentType::update();
+ velocity_.calculateVelocity();
+
return;
}
+ //! \copydoc Dumux::FVPressure1P::addOutputVtkFields(MultiWriter &writer)
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ ParentType::addOutputVtkFields(writer);
+ velocity_.addOutputVtkFields(writer);
+ }
+
//! Constructs a FVPressure2PAdaptive object
/**
* \param problem a problem class object
*/
FVPressure2PAdaptive(Problem& problem) :
- ParentType(problem), problem_(problem), gravity_(problem.gravity())
+ ParentType(problem), problem_(problem), velocity_(problem), gravity_(problem.gravity())
{
if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pglobal)
{
@@ -167,6 +178,7 @@ public:
private:
Problem& problem_;
+ FVVelocity<TypeTag, typename GET_PROP_TYPE(TypeTag, Velocity)> velocity_;
const GlobalPosition& gravity_; //!< vector including the gravity constant
Scalar density_[numPhases];
@@ -423,19 +435,22 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entrie
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
- this->A_[globalIdxI][globalIdxJ] -= 0.5 * entries[matrix];
- this->A_[globalIdxI][globalIdxK] -= 0.5 * entries[matrix];
+ //set off-diagonal
+ this->A_[globalIdxI][globalIdxJ] -= entries[matrix];
- // set off-diagonal entry
+ // set entries for cell J
this->A_[globalIdxJ][globalIdxI] -= entries[matrix];
- this->A_[globalIdxJ][globalIdxJ] += 0.5 * entries[matrix];
- this->A_[globalIdxJ][globalIdxK] += 0.5 * entries[matrix];
+ this->A_[globalIdxJ][globalIdxJ] += entries[matrix];
//set entries to zero -> matrix already written!
entries = 0.;
// std::cout<<"finished hanging node!\n";
}
+ else
+ {
+ entries = 0;
+ }
return;
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index a9d07c86bdcb1ef6473067e16c0a25679bc3d12d..522d19dcdd650553ca0dfd0ef9dbb673a9f7dd13 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -183,9 +183,9 @@ public:
{
int isIndex = isIt->indexInInside();
- fluxW[isIndex] = isIt->geometry().volume()
+ fluxW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
- fluxNW[isIndex] = isIt->geometry().volume()
+ fluxNW[isIndex] += isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
index 5c6d6e160cb033ca8a6adc5341fbe95c8cac615a..e04dd6d47f65280ba085b7bbc97088a2de29f253 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
@@ -154,7 +154,7 @@ public:
bool calculateVelocityInTransport()
{
- return GET_PROP_VALUE(TypeTag, CalculateVelocityInTransport);
+ return false;
}
private:
@@ -376,8 +376,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
//store potentials for further calculations (velocity, saturation, ...)
// these quantities only have correct sign (needed for upwinding)
// potentials are defined slightly different for adaptive scheme
- cellDataI.fluxData().setPotential(wPhaseIdx, isIndexI, potentialWIJ);
- cellDataI.fluxData().setPotential(nPhaseIdx, isIndexI, potentialNWIJ);
+ cellDataI.fluxData().addPotential(wPhaseIdx, isIndexI, potentialWIJ);
+ cellDataI.fluxData().addPotential(nPhaseIdx, isIndexI, potentialNWIJ);
cellDataJ.fluxData().setPotential(wPhaseIdx, isIndexJ, -potentialWIJ);
cellDataJ.fluxData().setPotential(nPhaseIdx, isIndexJ, -potentialNWIJ);
@@ -399,16 +399,6 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
}
- fractionalWIJ = (potentialWIJ > 0.) ? fractionalWI : fractionalWJ;
- fractionalNWIJ = (potentialNWIJ > 0.) ? fractionalNWI : fractionalNWJ;
- fractionalWIK = (potentialWIK > 0.) ? fractionalWI : fractionalWK;
- fractionalNWIK = (potentialNWIK > 0.) ? fractionalNWI : fractionalNWK;
-
- fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
- fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
- fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
- fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-
//calculate velocities and the gravity term
Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
@@ -456,12 +446,13 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNW);
cellDataJ.fluxData().setVelocityMarker(isIndexJ);
+ //times 0.5 because cell face with hanging node is called twice! Do not set marker because it should be called twice!
velocityW *= 0.5;
velocityNW *= 0.5;
- cellDataI.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW);
- cellDataI.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNW);
- cellDataI.fluxData().setVelocityMarker(isIndexI);
+ cellDataI.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW);
+ cellDataI.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNW);
}
+
return;
}
}
diff --git a/dumux/decoupled/2p/fluxData2p.hh b/dumux/decoupled/2p/fluxData2p.hh
index 204b742d7a9140fba72a406019fbdc6b7e14fb84..92fec2590382050210487dbbbcab163b0d4443f7 100644
--- a/dumux/decoupled/2p/fluxData2p.hh
+++ b/dumux/decoupled/2p/fluxData2p.hh
@@ -109,11 +109,16 @@ public:
return velocity_[phaseIdx][indexInInside];
}
- void setVelocity(int phaseIdx, int indexInInside, FieldVector& velocity)
+ void setVelocity(int phaseIdx, int indexInInside, const FieldVector& velocity)
{
velocity_[phaseIdx][indexInInside] = velocity;
}
+ void addVelocity(int phaseIdx, int indexInInside, const FieldVector& velocity)
+ {
+ velocity_[phaseIdx][indexInInside] += velocity;
+ }
+
void resetVelocity()
{
for (int i = 0; i < 2 * dim; i++)
@@ -121,6 +126,7 @@ public:
for (int j = 0; j < numPhases; j++)
{
velocity_[j][i] = 0.;
+ potential_[j][i] = 0.;
}
velocityMarker_[i] = false;
}
@@ -179,7 +185,10 @@ public:
{
potential_[indexInInside][phaseIdx] = pot;
}
-
+ void addPotential(int phaseIdx, int indexInInside, Scalar pot)
+ {
+ potential_[indexInInside][phaseIdx] += pot;
+ }
};
}
#endif
diff --git a/dumux/decoupled/2p/transport/gridadaptionindicator2p.hh b/dumux/decoupled/2p/transport/gridadaptionindicator2p.hh
index a961a4420b5d69f47a4bf6bb50082853113781cc..a3ab09712f0e7aa2a4c045ef2236274bede03e26 100644
--- a/dumux/decoupled/2p/transport/gridadaptionindicator2p.hh
+++ b/dumux/decoupled/2p/transport/gridadaptionindicator2p.hh
@@ -73,8 +73,8 @@ public:
if(indicatorVector_.size() != problem_.variables().cellDataGlobal().size())
{
indicatorVector_.resize(problem_.variables().cellDataGlobal().size());
- indicatorVector_ = -1e100;
};
+ indicatorVector_ = -1e100;
Scalar globalMax = -1e100;
Scalar globalMin = 1e100;
@@ -137,7 +137,7 @@ public:
}
}
- Scalar globaldelta = globalMax- globalMin;
+ Scalar globaldelta = globalMax - globalMin;
refineBound_ = refinetol_*globaldelta;
coarsenBound_ = coarsentol_*globaldelta;
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index d5e2358151017b8aa16c4b5dff895136116e0cde..42c682ce4974b2ebab76ebd1a4b45044ebf1f3c2 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -358,11 +358,13 @@ void FVPressure<TypeTag>::solve()
f_[idxFixPressureAtIndex_] = fixPressureAtIndex_;
}
+// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
+// printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
+
Solver solver(problem_);
solver.solve(A_, pressure_, f_);
-// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
-// printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
-// printvector(std::cout, pressure_, "pressure", "row", 200, 1, 3);
+
+// printvector(std::cout, pressure_, "pressure", "row", 200, 1, 3);
return;
}
diff --git a/dumux/decoupled/common/fv/fvvelocity.hh b/dumux/decoupled/common/fv/fvvelocity.hh
index dfeb98717d8ec32a8961b68a3b60d3deaa7b14cd..6974779e266f3a7e00db26da57755e1bf5c8ad83 100644
--- a/dumux/decoupled/common/fv/fvvelocity.hh
+++ b/dumux/decoupled/common/fv/fvvelocity.hh
@@ -112,13 +112,10 @@ void FVVelocity<TypeTag, Velocity>::calculateVelocity()
/************* handle interior face *****************/
if (isIt->neighbor())
{
- int globalIdxJ = problem_.variables().index(*(isIt->outside()));
+ int isIndex = isIt->indexInInside();
- //calculate only from one side
- if (globalIdxI > globalIdxJ)
- continue;
-
- velocity_.calculateVelocity(*isIt, cellDataI);
+ if (!cellDataI.fluxData().haveVelocity(isIndex))
+ velocity_.calculateVelocity(*isIt, cellDataI);
} // end neighbor