Commit fdb398a5 authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

SpatialParameters -> SpatialParams. To keep 2p and 1p modules, old properties...

SpatialParameters -> SpatialParams. To keep 2p and 1p modules, old properties are still working, but with a warning.
	reviewed by Markus
introduced naming conventions decided on dumux meeting for 2p2c modules.


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8115 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 9092285d
......@@ -63,8 +63,8 @@ NEW_TYPE_TAG(DecoupledTwoP, INHERITS_FROM(DecoupledModel))
//////////////////////////////////////////////////////////////////
// Property tags
//////////////////////////////////////////////////////////////////
NEW_PROP_TAG( SpatialParameters )
NEW_PROP_TAG( SpatialParams ); //!< The type of the soil properties object
NEW_PROP_TAG( SpatialParameters ); //!< DEPRECATED The old type of the soil properties object
; //!< The type of the spatial parameters object
NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
NEW_PROP_TAG(MaterialLawParams); //!< The material law parameters (extracted from the spatial parameters)
......@@ -172,6 +172,13 @@ private:
public:
typedef IsothermalImmiscibleFluidState<Scalar, FluidSystem> type;
};
//! DEPRECATED SpatialParameters property
#warning Please use SpatialParams instead of SpatialParameters
//TODO: next line enables old Models using SpatialParameters to work
// with base class impesproblem2p. If models are adapted,
// l180 should be replaced by l181 to deprecate old problems using SpatialParameters.
SET_TYPE_PROP(DecoupledTwoP, SpatialParams, typename GET_PROP_TYPE(TypeTag, SpatialParameters));
//SET_TYPE_PROP(DecoupledTwoP, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams));
/*!
* \brief Set the property for the material parameters by extracting
......
......@@ -58,7 +58,7 @@ class IMPESProblem2P : public IMPETProblem<TypeTag>
// material properties
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
enum {
......@@ -86,7 +86,7 @@ public:
gravity_(0)
{
newSpatialParams_ = true;
spatialParameters_ = new SpatialParameters(gridView);
spatialParams_ = new SpatialParams(gridView);
gravity_ = 0;
if (GET_PARAM(TypeTag, bool, EnableGravity))
......@@ -99,9 +99,9 @@ public:
* \param gridView The grid view
* \param spatialParameters SpatialParameters instantiation
*/
IMPESProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters)
IMPESProblem2P(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams)
: ParentType(timeManager, gridView),
gravity_(0),spatialParameters_(&spatialParameters)
gravity_(0),spatialParams_(&spatialParams)
{
newSpatialParams_ = false;
gravity_ = 0;
......@@ -115,7 +115,7 @@ public:
{
if (newSpatialParams_)
{
delete spatialParameters_;
delete spatialParams_;
}
}
......@@ -186,15 +186,21 @@ public:
/*!
* \brief Returns the spatial parameters object.
*/
SpatialParameters &spatialParameters()
{ return *spatialParameters_; }
SpatialParams &spatialParams()
{ return *spatialParams_; }
DUMUX_DEPRECATED_MSG("use spatialParams() method instead")
SpatialParams &spatialParameters()
{ return *spatialParams_; }
/*!
* \brief Returns the spatial parameters object.
*/
const SpatialParameters &spatialParameters() const
{ return *spatialParameters_; }
const SpatialParams &spatialParams() const
{ return *spatialParams_; }
DUMUX_DEPRECATED_MSG("use spatialParams() method instead")
const SpatialParams &spatialParameters() const
{ return *spatialParams_; }
// \}
......@@ -210,7 +216,7 @@ private:
GlobalPosition gravity_;
// fluids and material properties
SpatialParameters* spatialParameters_;
SpatialParams* spatialParams_;
bool newSpatialParams_;
};
......
......@@ -57,7 +57,7 @@ class IMPETProblem2P2C : public IMPESProblem2P<TypeTag>
typedef typename GridView::Traits::template Codim<0>::Entity Element;
// material properties
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
enum {
......@@ -85,10 +85,10 @@ public:
*
* \param timeManager The time manager
* \param gridView The grid view
* \param spatialParameters SpatialParameters instantiation
* \param spatialParams SpatialParameters instantiation
*/
IMPETProblem2P2C(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters)
: ParentType(timeManager, gridView, spatialParameters)
IMPETProblem2P2C(TimeManager &timeManager, const GridView &gridView, SpatialParams &spatialParams)
: ParentType(timeManager, gridView, spatialParams)
{ }
virtual ~IMPETProblem2P2C()
......
......@@ -67,7 +67,8 @@ NEW_TYPE_TAG(DecoupledTwoPTwoC, INHERITS_FROM(Pressure, Transport, IMPET));
// Property tags
//////////////////////////////////////////////////////////////////
NEW_PROP_TAG( Indices );
NEW_PROP_TAG( SpatialParameters ); //!< The type of the soil properties object
NEW_PROP_TAG( SpatialParams ); //!< The type of the soil properties object
NEW_PROP_TAG( SpatialParameters ); //!< DEPRECATED The old type of the soil properties object
NEW_PROP_TAG( EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG( PressureFormulation); //!< The formulation of the model
NEW_PROP_TAG( SaturationFormulation); //!< The formulation of the model
......@@ -153,14 +154,15 @@ SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false);
SET_BOOL_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, false);
SET_PROP(DecoupledTwoPTwoC, BoundaryMobility)
{
static const int value = DecoupledTwoPTwoCIndices<TypeTag>::satDependent;
};
{ static const int value = DecoupledTwoPTwoCIndices<TypeTag>::satDependent;};
SET_TYPE_PROP(DecoupledTwoPTwoC, Variables, VariableClass<TypeTag>);
SET_TYPE_PROP(DecoupledTwoPTwoC, CellData, CellData2P2C<TypeTag>);
SET_TYPE_PROP(DecoupledTwoPTwoC, FluidState, DecoupledTwoPTwoCFluidState<TypeTag>);
//! DEPRECATED SpatialParameters property
SET_TYPE_PROP(DecoupledTwoPTwoC, SpatialParameters, typename GET_PROP_TYPE(TypeTag, SpatialParams));
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableMultiPointFluxApproximationOnAdaptiveGrids, false);
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableVolumeIntegral, true);
......
......@@ -124,7 +124,7 @@ template<class TypeTag> class FVPressure2P2C
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Grid Grid;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::template Codim<0>::EntityPointer ElementPtr;
typedef typename GridView::Intersection Intersection;
// convenience shortcuts for Vectors/Matrices
......@@ -354,19 +354,19 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
const Intersection& intersection, const CellData& cellDataI, const bool first)
{
entries = 0.;
ElementPointer elementPointerI = intersection.inside();
ElementPtr elementPtrI = intersection.inside();
// get global coordinate of cell center
const GlobalPosition& globalPos = elementPointerI->geometry().center();
const GlobalPosition& globalPos = elementPtrI->geometry().center();
// cell volume & perimeter, assume linear map here
Scalar volume = elementPointerI->geometry().volume();
Scalar volume = elementPtrI->geometry().volume();
Scalar perimeter = cellDataI.perimeter();
const GlobalPosition& gravity_ = problem().gravity();
// get absolute permeability
FieldMatrix permeabilityI(problem().spatialParameters().intrinsicPermeability(*elementPointerI));
FieldMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
// get mobilities and fractional flow factors
Scalar fractionalWI=0, fractionalNWI=0;
......@@ -385,12 +385,12 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
Scalar faceArea = intersection.geometry().volume();
// access neighbor
ElementPointer neighborPointer = intersection.outside();
int globalIdxJ = problem().variables().index(*neighborPointer);
ElementPtr neighborPtr = intersection.outside();
int globalIdxJ = problem().variables().index(*neighborPtr);
CellData& cellDataJ = problem().variables().cellData(globalIdxJ);
// gemotry info of neighbor
const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
const GlobalPosition& globalPosNeighbor = neighborPtr->geometry().center();
// distance vector between barycenters
GlobalPosition distVec = globalPosNeighbor - globalPos;
......@@ -402,7 +402,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
unitDistVec /= dist;
FieldMatrix permeabilityJ
= problem().spatialParameters().intrinsicPermeability(*neighborPointer);
= problem().spatialParams().intrinsicPermeability(*neighborPtr);
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
......@@ -442,7 +442,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
{
// determine volume derivatives
if (!cellDataJ.hasVolumeDerivatives())
this->volumeDerivatives(globalPosNeighbor, *neighborPointer);
this->volumeDerivatives(globalPosNeighbor, *neighborPtr);
Scalar dv_dC1 = (cellDataJ.dv(wPhaseIdx)
+ cellDataI.dv(wPhaseIdx)) / 2; // dV/dm1= dv/dC^1
......@@ -580,8 +580,8 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
{
entries = 0.;
// get global coordinate of cell center
ElementPointer elementPointerI = intersection.inside();
const GlobalPosition& globalPos = elementPointerI->geometry().center();
ElementPtr elementPtrI = intersection.inside();
const GlobalPosition& globalPos = elementPtrI->geometry().center();
// get normal vector
const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
......@@ -613,7 +613,7 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
if (bcType.isDirichlet(Indices::pressureEqIdx))
{
// get absolute permeability
FieldMatrix permeabilityI(problem().spatialParameters().intrinsicPermeability(*elementPointerI));
FieldMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPtrI));
const GlobalPosition& gravity_ = problem().gravity();
//permeability vector at boundary
......@@ -678,10 +678,10 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
else if(GET_PROP_VALUE(TypeTag, BoundaryMobility) == Indices::permDependent)
{
lambdaWBound
= MaterialLaw::krw(problem().spatialParameters().materialLawParams(*elementPointerI),
= MaterialLaw::krw(problem().spatialParams().materialLawParams(*elementPtrI),
BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
lambdaNWBound
= MaterialLaw::krn(problem().spatialParameters().materialLawParams(*elementPointerI),
= MaterialLaw::krn(problem().spatialParams().materialLawParams(*elementPtrI),
BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
}
// get average density
......@@ -926,13 +926,13 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
}
//complete fluid state
fluidState.update(Z1, pressure, problem().spatialParameters().porosity(elementI), temperature_);
fluidState.update(Z1, pressure, problem().spatialParams().porosity(elementI), temperature_);
// iterations part in case of enabled capillary pressure
Scalar pc(0.), oldPc(0.);
if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
{
pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(elementI),
pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(elementI),
fluidState.saturation(wPhaseIdx));
int maxiter = 5; int iterout = -1;
//start iteration loop
......@@ -958,9 +958,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
//store old pc
oldPc = pc;
//update with better pressures
fluidState.update(Z1, pressure, problem().spatialParameters().porosity(elementI),
fluidState.update(Z1, pressure, problem().spatialParams().porosity(elementI),
problem().temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(elementI),
pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(elementI),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
//converge criterion
......@@ -976,9 +976,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
cellData.setViscosity(nPhaseIdx, FluidSystem::viscosity(fluidState, nPhaseIdx));
// initialize mobilities
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem().spatialParameters().materialLawParams(elementI), fluidState.saturation(wPhaseIdx))
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI), fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(wPhaseIdx));
cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem().spatialParameters().materialLawParams(elementI), fluidState.saturation(wPhaseIdx))
cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI), fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(nPhaseIdx));
// determine volume mismatch between actual fluid volume and pore volume
......@@ -992,7 +992,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
Scalar vol = massw / cellData.density(wPhaseIdx) + massn / cellData.density(nPhaseIdx);
if (problem().timeManager().timeStepSize() != 0)
{
cellData.volumeError()=(vol - problem().spatialParameters().porosity(elementI));
cellData.volumeError()=(vol - problem().spatialParams().porosity(elementI));
if (std::isnan(cellData.volumeError()))
{
......@@ -1000,7 +1000,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
<< "volErr[" << globalIdx << "] isnan: vol = " << vol
<< ", massw = " << massw << ", rho_l = " << cellData.density(wPhaseIdx)
<< ", massn = " << massn << ", rho_g = " << cellData.density(nPhaseIdx)
<< ", poro = " << problem().spatialParameters().porosity(elementI)
<< ", poro = " << problem().spatialParams().porosity(elementI)
<< ", dt = " << problem().timeManager().timeStepSize());
}
}
......
......@@ -117,7 +117,7 @@ class FVPressure2P2CMultiPhysics : public FVPressure2P2C<TypeTag>
// convenience shortcuts for Vectors/Matrices
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
......@@ -441,7 +441,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
// int globalIdxI = problem().variables().index(*elementPointerI);
// get absolute permeability
FieldMatrix permeabilityI(problem().spatialParameters().intrinsicPermeability(*elementPointerI));
DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPointerI));
// get normal vector
const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
......@@ -466,11 +466,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>
GlobalPosition unitDistVec(distVec);
unitDistVec /= dist;
FieldMatrix permeabilityJ
= problem().spatialParameters().intrinsicPermeability(*neighborPointer);
DimMatrix permeabilityJ
= problem().spatialParams().intrinsicPermeability(*neighborPointer);
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
DimMatrix meanPermeability(0);
Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
Dune::FieldVector<Scalar, dim> permeability(0);
......@@ -543,7 +543,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFluxOnBoundary(Dune::FieldVector<
if (bcType.isDirichlet(Indices::pressureEqIdx))
{
// get absolute permeability
FieldMatrix permeabilityI(problem().spatialParameters().intrinsicPermeability(*elementPointerI));
DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPointerI));
// get mobilities and fractional flow factors
Scalar lambdaI = cellDataI.mobility(phaseIdx);
......@@ -585,11 +585,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFluxOnBoundary(Dune::FieldVector<
{
if (phaseIdx == wPhaseIdx)
lambdaBound = MaterialLaw::krw(
problem().spatialParameters().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
problem().spatialParams().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
/ viscosityBound;
else
lambdaBound = MaterialLaw::krn(
problem().spatialParameters().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
problem().spatialParams().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
/ viscosityBound;
break;
}
......@@ -733,7 +733,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
// both phase pressures are necessary for the case 1p domain is assigned for
// the next 2p subdomain
PhaseVector pressure(0.);
Scalar pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(*eIt),
Scalar pc = MaterialLaw::pC(problem().spatialParams().materialLawParams(*eIt),
((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign Sw = 1 if wPhase present, else 0
if(pressureType == wPhaseIdx)
{
......@@ -771,11 +771,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
// initialize mobilities
if(presentPhaseIdx == wPhaseIdx)
cellData.setMobility(wPhaseIdx,
MaterialLaw::krw(problem().spatialParameters().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
MaterialLaw::krw(problem().spatialParams().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
/ cellData.viscosity(wPhaseIdx));
else
cellData.setMobility(nPhaseIdx,
MaterialLaw::krn(problem().spatialParameters().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
MaterialLaw::krn(problem().spatialParams().materialLawParams(*eIt), pseudoFluidState.saturation(wPhaseIdx))
/ cellData.viscosity(nPhaseIdx));
// error term handling
......@@ -783,7 +783,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
vol = sumConc / pseudoFluidState.density(presentPhaseIdx);
if (dt != 0)
cellData.volumeError() = (vol - problem().spatialParameters().porosity(*eIt));
cellData.volumeError() = (vol - problem().spatialParams().porosity(*eIt));
}
......
......@@ -277,8 +277,8 @@ public:
// get position, index
GlobalPosition globalPos = eIt->geometry().center();
int globalIdx = problem_.variables().index(*eIt);
poro_[globalIdx] = problem_.spatialParameters().porosity(*eIt);
perm_[globalIdx] = problem_.spatialParameters().intrinsicPermeability(*eIt)[0][0];
poro_[globalIdx] = problem_.spatialParams().porosity(*eIt);
perm_[globalIdx] = problem_.spatialParams().intrinsicPermeability(*eIt)[0][0];
}
*poroPtr = poro_;
*permPtr = perm_;
......@@ -343,6 +343,7 @@ private:
template<class TypeTag>
void FVPressureCompositional<TypeTag>::initialize(bool solveTwice)
{
// inizialize matrix, RHS, and condition Matrix
ParentType::initialize();
// initialguess: set saturations, determine visco and mobility for initial pressure equation
......@@ -470,12 +471,12 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
if (icFormulation == Indices::saturation) // saturation initial condition
{
sat_0 = problem_.initSat(*eIt);
fluidState.satFlash(sat_0, pressure, problem_.spatialParameters().porosity(*eIt), temperature_);
fluidState.satFlash(sat_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
Scalar Z1_0 = problem_.initConcentration(*eIt);
fluidState.update(Z1_0, pressure, problem_.spatialParameters().porosity(*eIt), temperature_);
fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
}
else if(compositional) //means we regard compositional effects since we know an estimate pressure field
......@@ -487,7 +488,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
Scalar pc=0.;
if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
{
pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt),
pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt),
sat_0);
}
else
......@@ -510,7 +511,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
}
}
fluidState.satFlash(sat_0, pressure, problem_.spatialParameters().porosity(*eIt), temperature_);
fluidState.satFlash(sat_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
......@@ -548,16 +549,16 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
//store old pc
Scalar oldPc = pc;
//update with better pressures
fluidState.update(Z1_0, pressure, problem_.spatialParameters().porosity(*eIt),
fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt),
problem_.temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt),
pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
//converge criterion
if (abs(oldPc-pc)<10)
iter = maxiter;
pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt),
pc = MaterialLaw::pC(problem_.spatialParams().materialLawParams(*eIt),
fluidState.saturation(wPhaseIdx));
}
}
......@@ -565,10 +566,10 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
{
pressure[wPhaseIdx] = pressure[nPhaseIdx]
= this->pressure()[globalIdx];
fluidState.update(Z1_0, pressure, problem_.spatialParameters().porosity(*eIt), temperature_);
fluidState.update(Z1_0, pressure, problem_.spatialParams().porosity(*eIt), temperature_);
}
fluidState.calculateMassConcentration(problem_.spatialParameters().porosity(*eIt));
fluidState.calculateMassConcentration(problem_.spatialParams().porosity(*eIt));
} //end conc initial condition
} //end compositional
......@@ -580,9 +581,9 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
cellData.setViscosity(nPhaseIdx, FluidSystem::viscosity(fluidState, nPhaseIdx));
// initialize mobilities
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), fluidState.saturation(wPhaseIdx))
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(wPhaseIdx));
cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), fluidState.saturation(wPhaseIdx))
cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(nPhaseIdx));
// calculate perimeter used as weighting factor
......@@ -652,7 +653,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
for(int phaseIdx = 0; phaseIdx< numPhases; phaseIdx++)
specificVolume += cellData.phaseMassFraction(phaseIdx) / cellData.density(phaseIdx);
Scalar volalt = mass.one_norm() * specificVolume;
// volalt = cellData.volumeError()+problem_.spatialParameters().porosity(element);
// volalt = cellData.volumeError()+problem_.spatialParams().porosity(element);
// = \sum_{\kappa} C^{\kappa} + \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
/**********************************
......@@ -678,7 +679,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
p_ += pressure;
Scalar Z1 = mass[0] / mass.one_norm();
updFluidState.update(Z1,
p_, problem_.spatialParameters().porosity(element), temperature_);
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
for(int phaseIdx = 0; phaseIdx< numPhases; phaseIdx++)
......@@ -693,7 +694,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
p_ -= 2*incp;
updFluidState.update(Z1,
p_, problem_.spatialParameters().porosity(element), temperature_);
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
for(int phaseIdx = 0; phaseIdx< numPhases; phaseIdx++)
......@@ -708,21 +709,21 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
}
// numerical derivative of fluid volume with respect to mass of components
for (int comp = 0; comp<numComponents; comp++)
for (int compIdx = 0; compIdx<numComponents; compIdx++)
{
mass[comp] += massIncrement[comp];
mass[compIdx] += massIncrement[compIdx];
Z1 = mass[0] / mass.one_norm();
updFluidState.update(Z1, pressure, problem_.spatialParameters().porosity(element), temperature_);
updFluidState.update(Z1, pressure, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
for(int phaseIdx = 0; phaseIdx< numPhases; phaseIdx++)
specificVolume += updFluidState.phaseMassFraction(phaseIdx) / updFluidState.density(phaseIdx);
cellData.dv(comp) = ((mass.one_norm() * specificVolume) - volalt) / massIncrement[comp];
mass[comp] -= massIncrement[comp];
cellData.dv(compIdx) = ((mass.one_norm() * specificVolume) - volalt) / massIncrement[compIdx];
mass[compIdx] -= massIncrement[compIdx];
//check routines if derivatives are meaningful
if (isnan(cellData.dv(comp)) || isinf(cellData.dv(comp)) )
if (isnan(cellData.dv(compIdx)) || isinf(cellData.dv(compIdx)) )
{
DUNE_THROW(Dune::MathError, "NAN/inf of dV_dm. If that happens in first timestep, try smaller firstDt!");
}
......
......@@ -97,6 +97,7 @@ class FVTransport2P2C
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
......@@ -292,7 +293,7 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt,
// account for porosity in fluxes for time-step
sumfactorin = std::max(sumfactorin,sumfactorout)
/ problem().spatialParameters().porosity(*eIt);
/ problem().spatialParams().porosity(*eIt);
if ( 1./sumfactorin < dt)
{
......@@ -348,25 +349,25 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
fluxEntries = 0.;
timestepFlux = 0.;
// cell information
ElementPointer elementI= intersection.inside();
int globalIdxI = problem().variables().index(*elementI);
ElementPointer elementPtrI= intersection.inside();
int globalIdxI = problem().variables().index(*elementPtrI);
// get position
const GlobalPosition globalPos = elementI->geometry().center();
const GlobalPosition globalPos = elementPtrI->geometry().center();
const GlobalPosition& gravity_ = problem().gravity();
// cell volume, assume linear map here
Scalar volume = elementI->geometry().volume();
Scalar volume = elementPtrI->geometry().volume();
// get values of cell I
Scalar pressI = problem().pressureModel().pressure(globalIdxI);
Scalar pcI = cellDataI.capillaryPressure();
Dune::FieldMatrix<Scalar,dim,dim> K_I(problem().spatialParameters().intrinsicPermeability(*elementI));
DimMatrix K_I(problem().spatialParams().intrinsicPermeability(*elementPtrI));
Scalar SwmobI = std::max((cellDataI.saturation(wPhaseIdx)
- problem().spatialParameters().materialLawParams(*elementI).Swr())
- problem().spatialParams().materialLawParams(*elementPtrI).Swr())
, 1e-2);
Scalar SnmobI = std::max((cellDataI.saturation(nPhaseIdx)
- problem().spatialParameters().materialLawParams(*elementI).Snr())
- problem().spatialParams().materialLawParams(*elementPtrI).Snr())
, 1e-2);
Scalar densityWI (0.), densityNWI(0.);
......@@ -386,12 +387,12 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
Scalar potentialW(0.), potentialNW(0.);