Skip to content
Snippets Groups Projects
Commit 892e61d0 authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

remove volume-derivatives functions of the old model which was totally accidentially comitted.

This will also remove the strange and unnecessary warning.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7823 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent dd43e872
No related branches found
No related tags found
No related merge requests found
......@@ -605,8 +605,6 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
* \param globalPos The global position of the current element
* \param ep A pointer to the current element
*/
#if 1
#warning decide which one!!!
template<class TypeTag>
void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& globalPos, const Element& element)
{
......@@ -717,125 +715,6 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
}
cellData.confirmVolumeDerivatives();
}
#else
template<class TypeTag>
void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& globalPos, const Element& element)
{
// cell index
int globalIdx = problem_.variables().index(element);
CellData& cellData = problem_.variables().cellData(globalIdx);
// get cell temperature
Scalar temperature_ = problem_.temperatureAtPos(globalPos);
// initialize an Fluid state for the update
FluidState updFluidState;
/**********************************
* a) get necessary variables
**********************************/
//determine phase pressures from primary pressure variable
PhaseVector pressure(0.);
switch (pressureType)
{
//TODO: use pressure from cellData here
case pw:
{
pressure[wPhaseIdx] = this->pressure()[globalIdx];
pressure[nPhaseIdx] = this->pressure()[globalIdx]
+ cellData.capillaryPressure();
break;
}
case pn:
{
pressure[wPhaseIdx] = this->pressure()[globalIdx]
- cellData.capillaryPressure();
pressure[nPhaseIdx] = this->pressure()[globalIdx];
break;
}
}
// mass of components inside the cell
ComponentVector mass(0.);
mass[0] = cellData.massConcentration(wCompIdx);
mass[1] = cellData.massConcentration(nCompIdx);
// shortcuts for density
Scalar densityW = cellData.density(wPhaseIdx);
Scalar densityNW = cellData.density(nPhaseIdx);
// actual fluid volume
Scalar volalt = (mass[0]+mass[1])
* (cellData.phaseMassFraction(wPhaseIdx) / densityW
+ cellData.phaseMassFraction(nPhaseIdx) / densityNW);
/**********************************
* b) define increments
**********************************/
// increments for numerical derivatives
ComponentVector massIncrement(0.);
massIncrement[0] = updateEstimate_[wCompIdx][globalIdx];
massIncrement[1] = updateEstimate_[nCompIdx][globalIdx];
if(fabs(massIncrement[0]) < 1e-8 * densityW)
massIncrement[0] = 1e-8* densityW;
if(fabs(massIncrement[1]) < 1e-8 * densityNW)
massIncrement[1] = 1e-8 * densityNW;
Scalar incp = 1e-2;
/**********************************
* c) Secant method for derivatives
**********************************/
// numerical derivative of fluid volume with respect to pressure
PhaseVector p_(incp);
p_ += pressure;
Scalar Z1 = mass[0] / (mass[0] + mass[1]);
updFluidState.update(Z1,
p_, problem_.spatialParameters().porosity(globalPos, element), temperature_);
cellData.dv_dp() = (((mass[0]+mass[1]) * (updFluidState.phaseMassFraction(wPhaseIdx) /updFluidState.density(wPhaseIdx)
+ updFluidState.phaseMassFraction(nPhaseIdx) /updFluidState.density(nPhaseIdx))) - volalt) /incp;
if (cellData.dv_dp()>0)
{
// dV_dp > 0 is unphysical: Try inverse increment for secant
Dune::dinfo << "dv_dp larger 0 at Idx " << globalIdx << " , try and invert secant"<< std::endl;
p_ -= 2*incp;
updFluidState.update(Z1,
p_, problem_.spatialParameters().porosity(globalPos, element), temperature_);
cellData.dv_dp() = (((mass[0]+mass[1]) * (updFluidState.phaseMassFraction(wPhaseIdx) /updFluidState.density(wPhaseIdx)
+ updFluidState.phaseMassFraction(nPhaseIdx) /updFluidState.density(nPhaseIdx))) - volalt) /incp;
// dV_dp > 0 is unphysical: Try inverse increment for secant
if (cellData.dv_dp()>0)
{
Dune::dinfo << "dv_dp still larger 0 after inverting secant"<< std::endl;
}
}
// numerical derivative of fluid volume with respect to mass of components
for (int comp = 0; comp<numComponents; comp++)
{
mass[comp] += massIncrement[comp];
Z1 = mass[0] / (mass[0] + mass[1]);
updFluidState.update(Z1, pressure, problem_.spatialParameters().porosity(globalPos, element), temperature_);
cellData.dv(comp) = ((mass[0]+mass[1])
* (updFluidState.phaseMassFraction(wPhaseIdx) / updFluidState.density(wPhaseIdx)
+ updFluidState.phaseMassFraction(nPhaseIdx) / updFluidState.density(nPhaseIdx)) - volalt)
/ massIncrement[comp];
mass[comp] -= massIncrement[comp];
//check routines if derivatives are meaningful
if (isnan(cellData.dv(comp)) || isinf(cellData.dv(comp)) )
{
DUNE_THROW(Dune::MathError, "NAN/inf of dV_dm. If that happens in first timestep, try smaller firstDt!");
}
}
cellData.confirmVolumeDerivatives();
}
#endif
}//end namespace Dumux
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment