Commit 97f28e7c authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

upwind and calculation of velocity consistent

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5286 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent fa2b9c19
......@@ -606,15 +606,6 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
Scalar graddv_dC1 = (dV_[wPhaseIdx][globalIdxJ] - dV_[wPhaseIdx][globalIdxI]) / dist;
Scalar graddv_dC2 = (dV_[nPhaseIdx][globalIdxJ] - dV_[nPhaseIdx][globalIdxI]) / dist;
// potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
// potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
//
// densityW = (potentialW > 0.) ? densityWI : densityWJ;
// densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
//
// densityW = (potentialW == 0.) ? rhoMeanW : densityW;
// densityNW = (potentialNW == 0.) ? rhoMeanNW : densityNW;
//jochen: central weighting for gravity term
densityW = rhoMeanW; densityNW = rhoMeanNW;
......@@ -980,10 +971,10 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
Scalar hifac = 0.;
hifac /= fac;
if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxErr))
if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxErr) && (!problem_.timeManager().willBeFinished()))
{
if (erri <= x_mi * maxErr)
f_[globalIdxI] += errorCorrection[globalIdxI] = fac* (1-x_mi*(lofac/fac-1)/(x_lo-x_mi) + (lofac/fac-1)/(x_lo-x_mi)*erri/maxErr)
f_[globalIdxI] += errorCorrection[globalIdxI] = fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxErr)
* problem_.variables().volErr()[globalIdxI] * volume;
else
f_[globalIdxI] += errorCorrection[globalIdxI] = fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxErr)
......
......@@ -320,9 +320,9 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
meanK_.umv(unitDistVec,K);
// velocities
double potentialW = (unitOuterNormal * distVec) * (pressI - pressJ) / (dist * dist);
double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
potentialW += densityW_mean * (unitDistVec * gravity_);
double potentialW = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
double potentialNW = potentialW + ((K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
potentialW += ((K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
// upwind mobility
double lambdaW, lambdaN;
......@@ -336,8 +336,8 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
else
lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
double velocityW = lambdaW * potentialW;
double velocityN = lambdaN * potentialW;
// standardized velocity
double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
......@@ -395,21 +395,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
if (bcTypeTransport_ == BoundaryConditions::dirichlet)
{
//get dirichlet pressure boundary condition
Scalar pressBound = 0.;
Scalar pcBound = problem_.variables().capillaryPressure(globalIdxI);
switch (pressureType)
{
case pw:
{
pressBound = problem_.dirichletPress(globalPosFace, *isIt);
break;
}
case pn:
{
pressBound = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
break;
}
}
Scalar pressBound = problem_.dirichletPress(globalPosFace, *isIt);
// read boundary values
BoundaryConditions2p2c::Flags bctype = problem_.bcFormulation(globalPosFace, *isIt);
......@@ -439,16 +425,20 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
pressBound, BCfluidState);
Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
problem_.temperature(globalPosFace, *eIt),
pressBound+pcBound, BCfluidState);
pressBound, BCfluidState);
// average
double densityW_mean = (densityWI + densityWBound) / 2;
double densityNW_mean = (densityNWI + densityNWBound) / 2;
// prepare K
Dune::FieldVector<Scalar,dim> K(0);
K_I.umv(unitDistVec,K);
// velocities
double potentialW = (unitOuterNormal * distVec) * (pressI - pressBound) / (dist * dist);
double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
potentialW += densityW_mean * (unitDistVec * gravity_);
double potentialW = (K * unitOuterNormal) * (pressI - pressBound) / (dist);
double potentialNW = potentialW + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;
potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
// do upwinding for lambdas
double lambdaW, lambdaN;
......@@ -474,20 +464,11 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
problem_.spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
/ viscosityNWBound;
}
// prepare K
Dune::FieldVector<Scalar,dim> K(0);
K_I.umv(unitDistVec,K);
double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressBound)
/ (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressBound)
/ (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
// standardized velocity
double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
double velocityIJn = std::max( velocityN* faceArea / volume, 0.0);
double velocityJIw = std::max(-lambdaW * potentialW * faceArea / volume, 0.0);
double velocityIJw = std::max( lambdaW * potentialW * faceArea / volume, 0.0);
double velocityJIn = std::max(-lambdaN * potentialNW * faceArea / volume, 0.0);
double velocityIJn = std::max( lambdaN * potentialNW* faceArea / volume, 0.0);
// for timestep control
factor[0] = velocityJIw + velocityJIn;
......@@ -539,8 +520,8 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}//end boundary
// correct update Factor by volume error
#ifdef errorInTransport
updFactor[wCompIdx] *= problem_.variables().volErr(globalIdxI);
updFactor[nCompIdx] *= problem_.variables().volErr(globalIdxI);
updFactor[wCompIdx] *= (1+problem_.variables().volErr(globalIdxI));
updFactor[nCompIdx] *= (1+problem_.variables().volErr(globalIdxI));
#endif
// add to update vector
updateVec[wCompIdx][globalIdxI] += updFactor[wCompIdx];
......
......@@ -303,12 +303,12 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
meanK_.umv(unitDistVec,K);
// velocities
double potentialW = (unitOuterNormal * distVec) * (pressI - pressJ) / (dist * dist);
double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
potentialW += densityW_mean * (unitDistVec * gravity_);
double potentialW = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
double potentialNW = potentialW + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;
potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
// lambda = 2 * lambdaI * lambdaJ / (lambdaI + lambdaJ);
double lambdaW, lambdaN;
double lambdaW, lambdaNW;
if (potentialW >= 0.)
lambdaW = problem_.variables().mobilityWetting(globalIdxI);
......@@ -316,18 +316,15 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
if (potentialNW >= 0.)
lambdaN = problem_.variables().mobilityNonwetting(globalIdxI);
lambdaNW = problem_.variables().mobilityNonwetting(globalIdxI);
else
lambdaN = problem_.variables().mobilityNonwetting(globalIdxJ);
lambdaNW = problem_.variables().mobilityNonwetting(globalIdxJ);
double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressJ) / (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
// standardized velocity
double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
double velocityIJn = std::max( velocityN * faceArea / volume, 0.0);
// calculate and standardized velocity
double velocityJIw = std::max(-lambdaW * potentialW * faceArea / volume, 0.0);
double velocityIJw = std::max( lambdaW * potentialW * faceArea / volume, 0.0);
double velocityJIn = std::max(-lambdaNW * potentialNW * faceArea / volume, 0.0);
double velocityIJn = std::max( lambdaNW * potentialNW * faceArea / volume, 0.0);
// for timestep control
factor[0] = velocityJIw + velocityJIn;
......@@ -378,21 +375,7 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
if (bcTypeTransport_ == BoundaryConditions::dirichlet)
{
//get dirichlet pressure boundary condition
Scalar pressBound = 0.;
Scalar pcBound = problem_.variables().capillaryPressure(globalIdxI);
switch (pressureType)
{
case pw:
{
pressBound = problem_.dirichletPress(globalPosFace, *isIt);
break;
}
case pn:
{
pressBound = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
break;
}
}
Scalar pressBound = problem_.dirichletPress(globalPosFace, *isIt);
// read boundary values
BoundaryConditions2p2c::Flags bctype = problem_.bcFormulation(globalPosFace, *isIt);
......@@ -419,22 +402,26 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
pressBound, fluidState);
Scalar densityNWBound = FluidSystem::phaseDensity(nPhaseIdx,
problem_.temperature(globalPosFace, *eIt),
pressBound+pcBound, fluidState);
pressBound, fluidState);
Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
problem_.temperature(globalPosFace, *eIt),
pressBound, fluidState);
Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
problem_.temperature(globalPosFace, *eIt),
pressBound+pcBound, fluidState);
pressBound, fluidState);
// average
double densityW_mean = (densityWI + densityWBound) / 2;
double densityNW_mean = (densityNWI + densityNWBound) / 2;
// prepare K
Dune::FieldVector<Scalar,dim> K(0);
K_I.umv(unitDistVec,K);
// velocities
double potentialW = (unitOuterNormal * distVec) * (pressI - pressBound) / (dist * dist);
double potentialNW = potentialW + densityNW_mean * (unitDistVec * gravity_);
potentialW += densityW_mean * (unitDistVec * gravity_);
double potentialW = (K * unitOuterNormal) * (pressI - pressBound) / (dist);
double potentialNW = potentialW + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean;
potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
double lambdaW, lambdaN;
......@@ -460,21 +447,11 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
problem_.spatialParameters().materialLawParams(globalPos, *eIt), fluidState.saturation(wPhaseIdx))
/ viscosityNWBound;
}
// prepare K
Dune::FieldVector<Scalar,dim> K(0);
K_I.umv(unitDistVec,K);
double velocityW = lambdaW * ((K * unitOuterNormal) * (pressI - pressBound)
/ (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean);
double velocityN = lambdaN * ((K * unitOuterNormal) * (pressI - pressBound)
/ (dist) + (K * gravity_) * (unitOuterNormal * unitDistVec) * densityNW_mean);
// standardized velocity
double velocityJIw = std::max(-velocityW * faceArea / volume, 0.0);
double velocityIJw = std::max( velocityW * faceArea / volume, 0.0);
double velocityJIn = std::max(-velocityN * faceArea / volume, 0.0);
double velocityIJn = std::max( velocityN* faceArea / volume, 0.0);
double velocityJIw = std::max(-lambdaW * potentialW * faceArea / volume, 0.0);
double velocityIJw = std::max( lambdaW * potentialW * faceArea / volume, 0.0);
double velocityJIn = std::max(-lambdaN * potentialNW * faceArea / volume, 0.0);
double velocityIJn = std::max( lambdaN * potentialNW* faceArea / volume, 0.0);
// for timestep control
factor[0] = velocityJIw + velocityJIn;
......
......@@ -147,14 +147,15 @@ public:
void serializeEntity(std::ostream &outstream, const Element &element)
{
Dune::dwarn << "here, rather total concentrations should be used!" << std::endl;
int globalIdx = this->elementMapper().map(element);
outstream << this->pressure()[globalIdx] << " " << saturation_[globalIdx];
outstream << this->pressure()[globalIdx] << " " << totalConcentration(globalIdx, 0)
<< " " << totalConcentration(globalIdx,1);
}
void deserializeEntity(std::istream &instream, const Element &element)
{
int globalIdx = this->elementMapper().map(element);
instream >> this->pressure()[globalIdx] >> saturation_[globalIdx];
instream >> this->pressure()[globalIdx] >> totalConcentration(globalIdx, 0)
>> totalConcentration(globalIdx,1);
}
//@}
......@@ -220,9 +221,9 @@ public:
*pressure = this->pressure();
*saturation = saturation_;
if (GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation)) == GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureW)
writer.addCellData(pressure, "pressure w-phase");
writer.addCellData(pressure, "pressure w_phase");
else
writer.addCellData(pressure, "pressure nw-phase");
writer.addCellData(pressure, "pressure nw_phase");
writer.addCellData(saturation, "water saturation");
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
{
......@@ -236,8 +237,8 @@ public:
ScalarSolutionType *totalConcentration2 = writer.template createField<Scalar, 1> (size_);
*totalConcentration1 = totalConcentration_[0];
*totalConcentration2 = totalConcentration_[1];
writer.addCellData(totalConcentration1, "totalConcentration w-Comp");
writer.addCellData(totalConcentration2, "totalConcentration n-Comp");
writer.addCellData(totalConcentration1, "totalConcentration w_Comp");
writer.addCellData(totalConcentration2, "totalConcentration n_Comp");
// numerical stuff
ScalarSolutionType *volErr = writer.template createField<Scalar, 1> (size_);
......@@ -249,8 +250,8 @@ public:
ScalarSolutionType *massfraction1NW = writer.template createField<Scalar, 1> (size_);
*massfraction1W = massfrac_[0];
*massfraction1NW = massfrac_[1];
writer.addCellData(massfraction1W, "massfraction1 in w-phase");
writer.addCellData(massfraction1NW, "massfraction1NW nw-phase");
writer.addCellData(massfraction1W, "massfraction1 in w_phase");
writer.addCellData(massfraction1NW, "massfraction1NW nw_phase");
// phase properties
ScalarSolutionType *mobilityW = writer.template createField<Scalar, 1> (size_);
......@@ -265,12 +266,12 @@ public:
*densityNW = density_[1];
*numdensityW = density_[0];
*numdensityNW = density_[1];
writer.addCellData(mobilityW, "mobility w-phase");
writer.addCellData(mobilityNW, "mobility nw-phase");
writer.addCellData(densityW, "density w-phase");
writer.addCellData(densityNW, "density nw-phase");
writer.addCellData(numdensityW, "numerical density (mass/volume) w-phase");
writer.addCellData(numdensityNW, "numerical density (mass/volume) nw-phase");
writer.addCellData(mobilityW, "mobility w_phase");
writer.addCellData(mobilityNW, "mobility nw_phase");
writer.addCellData(densityW, "density w_phase");
writer.addCellData(densityNW, "density nw_phase");
writer.addCellData(numdensityW, "numerical density (mass/volume) w_phase");
writer.addCellData(numdensityNW, "numerical density (mass/volume) nw_phase");
}
if (codim_ == dim)
{
......@@ -299,7 +300,7 @@ public:
}
//! Returs a reference to the total concentration vector
const Scalar& totalConcentration() const
const TransportSolutionType& totalConcentration() const
{
return totalConcentration_;
}
......@@ -358,6 +359,11 @@ public:
}
//! Return capillary pressure vector
const ScalarSolutionType& capillaryPressure() const
{
return capillaryPressure_;
}
//! Return capillary pressure of specific Element
/*! \param Idx Element index*/
Scalar& capillaryPressure(int Idx)
{
......
......@@ -203,7 +203,7 @@ public:
{
problem.variables().deserialize<Restarter> (res);
//update constitutive functions
problem.pressureModel().initialize();
problem.pressureModel().updateMaterialLaws();
}
//! Constructs an IMPET object
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment