diff --git a/dumux/decoupled/common/decoupledproperties.hh b/dumux/decoupled/common/decoupledproperties.hh
index bb8dd62c932fc47157823c14103acd1723f96ad2..30bb05d5d8abf3c519ef81aa4b032124fae7f674 100644
--- a/dumux/decoupled/common/decoupledproperties.hh
+++ b/dumux/decoupled/common/decoupledproperties.hh
@@ -52,6 +52,7 @@ NEW_PROP_TAG( Scalar);
//! This means vectors of primary variables, solution functions on the
//! grid, and elements, and shape functions.
NEW_PROP_TAG( SolutionTypes);
+NEW_PROP_TAG( TransportSolutionType);
NEW_PROP_TAG( Grid); //!< The type of the DUNE grid
NEW_PROP_TAG( GridView); //!< The type of the grid view
@@ -61,6 +62,7 @@ NEW_PROP_TAG( ReferenceElements); //!< DUNE reference elements to be used
NEW_PROP_TAG( Problem); //!< The type of the problem
NEW_PROP_TAG( Model); //!< The type of the discretizations
NEW_PROP_TAG( NumPhases); //!< Number of phases in the system
+NEW_PROP_TAG( NumComponents); //!< Number of components in the system
NEW_PROP_TAG( Variables); //!< The type of the container of global variables
NEW_PROP_TAG( LocalStiffness); //!< The type of communication needed for the mimetic operator
}
@@ -153,7 +155,9 @@ SET_PROP(DecoupledModel, SolutionTypes)
enum
{
- dim = GridView::dimension, numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
+ dim = GridView::dimension,
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents))
};
template<int dim>
@@ -187,12 +191,22 @@ public:
* This defines the primary and secondary variable vectors at each degree of freedom.
*/
typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarSolution;//!<type for vector of scalars
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, numPhases> > PhaseProperty;//!<type for vector of phase properties
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, numPhases> > FluidProperty;//!<type for vector of fluid properties
+ typedef Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numComponents> ComponentProperty;//!<type for vector of phase properties
+ typedef Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numPhases> PhaseProperty;//!<type for vector of phase properties
+ typedef Dune::FieldVector<Dune::BlockVector<Dune::FieldVector<Scalar,1> >, numPhases> FluidProperty;//!<type for vector of fluid properties: Vector[element][phase]
typedef Dune::BlockVector<Dune::FieldVector<Dune::FieldVector<Scalar, numPhases>, 2*dim > > PhasePropertyElemFace;//!<type for vector of vectors (of size 2 x dimension) of scalars
typedef Dune::BlockVector<Dune::FieldVector<Dune::FieldVector<Scalar, dim>, 2*dim > > DimVecElemFace;//!<type for vector of vectors (of size 2 x dimension) of vector (of size dimension) of scalars
};
+SET_PROP_DEFAULT(TransportSolutionType)
+{
+ private:
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionType;
+
+ public:
+ typedef typename SolutionType::ScalarSolution type;//!<type for vector of scalar properties
+};
+
// \}
}
diff --git a/dumux/decoupled/common/impes.hh b/dumux/decoupled/common/impes.hh
index b94ad30adff019a8e413ea73bc3bf3455c967ba2..eefd13b48c29d1fb2a279283fecea5740fa89538 100644
--- a/dumux/decoupled/common/impes.hh
+++ b/dumux/decoupled/common/impes.hh
@@ -60,6 +60,7 @@ template<class TypeTag> class IMPES
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
@@ -82,6 +83,11 @@ public:
problem.pressureModel().initial();
problem.pressureModel().calculateVelocity();
+ //update constitutive functions
+ problem.pressureModel().updateMaterialLaws();
+
+ Dune::dinfo << "IMPES: initialization done." << std::endl;
+
return;
}
@@ -95,19 +101,18 @@ public:
* Calculates the new pressure and velocity and determines the time step size and the saturation update for the explicit time step
* Called from Dumux::Timestep.
*/
- void update(const Scalar t, Scalar& dt, ScalarSolutionType& updateVec)
+ void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
{
+ // the method is valid for any transported quantity. For sake of brevity,
+ // variable names refer to the saturation as an exemplary transported quantity.
int satSize = problem.variables().gridSize();
- ScalarSolutionType saturation(problem.variables().saturation());
- ScalarSolutionType satOldIter(problem.variables().saturation());
- ScalarSolutionType satHelp(satSize);
- ScalarSolutionType satDiff(satSize);
- ScalarSolutionType updateOldIter(satSize);
- ScalarSolutionType updateHelp(satSize);
- ScalarSolutionType updateDiff(satSize);
-
- //update constitutive functions
- problem.pressureModel().updateMaterialLaws();
+ TransportSolutionType saturation(problem.variables().transportedQuantity());
+ TransportSolutionType satOldIter(problem.variables().transportedQuantity());
+ TransportSolutionType satHelp(satSize);
+ TransportSolutionType satDiff(satSize);
+ TransportSolutionType updateOldIter(satSize);
+ TransportSolutionType updateHelp(satSize);
+ TransportSolutionType updateDiff(satSize);
bool converg = false;
int iter = 0;
@@ -124,12 +129,12 @@ public:
problem.pressureModel().calculateVelocity();
//calculate saturation defect
- problem.saturationModel().update(t, dt, updateVec, cFLFactor_, true);
+ problem.saturationModel().update(t, dt, updateVec, true);
if (iterFlag_)
{ // only needed if iteration has to be done
updateHelp = updateVec;
- saturation = problem.variables().saturation();
+ saturation = problem.variables().transportedQuantity();
saturation += (updateHelp *= (dt * cFLFactor_));
saturation *= omega_;
satHelp = satOldIter;
diff --git a/dumux/decoupled/common/impesproblem.hh b/dumux/decoupled/common/impesproblem.hh
index 4f3706ca5bf69112ae0b290e8b1b38aa7971cce3..843fa1b82912bef861de3859d96dfb898d1043ba 100644
--- a/dumux/decoupled/common/impesproblem.hh
+++ b/dumux/decoupled/common/impesproblem.hh
@@ -61,6 +61,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Variables)) Variables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) IMPESModel;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureModel)) PressureModel;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SaturationModel)) SaturationModel;
@@ -155,10 +156,9 @@ public:
void timeIntegration()
{
// allocate temporary vectors for the updates
- typedef typename IMPESModel::SolutionType Solution;
- Solution k1 = (*asImp_()).variables().saturation();
+ typedef TransportSolutionType Solution;
+ Solution k1 = (*asImp_()).variables().transportedQuantity();
- dt_ = 1e100;
Scalar t = timeManager().time();
// obtain the first update and the time step size
@@ -166,10 +166,23 @@ public:
//make sure t_old + dt is not larger than tend
dt_ = std::min(dt_, timeManager().episodeMaxTimeStepSize());
+
+ // check if we are in first TS and an initialDt was assigned
+ if (t==0. && timeManager().timeStepSize()!=0.)
+ {
+ // check if assigned initialDt is in accordance with dt_ from first transport step
+ if (timeManager().timeStepSize() > dt_)
+ Dune::dwarn << "initial timestep of size " << timeManager().timeStepSize()
+ << "is larger then dt_= "<<dt_<<" from transport" << std::endl;
+ // internally assign next tiestep size
+ dt_ = std::min(dt_, timeManager().timeStepSize());
+ }
+
+ //assign next tiestep size
timeManager().setTimeStepSize(dt_);
// explicit Euler: Sat <- Sat + dt*N(Sat)
- (*asImp_()).variables().saturation() += (k1 *= dt_);
+ (*asImp_()).variables().transportedQuantity() += (k1 *= timeManager().timeStepSize());
}
/*!
@@ -181,7 +194,9 @@ public:
* current solution to disk.
*/
void postTimeStep()
- { };
+ {
+ (*asImp_()).pressureModel().updateMaterialLaws();
+ };
/*!
* \brief Returns the current time step size [seconds].