diff --git a/dumux/decoupled/2p2c/2p2cproblem.hh b/dumux/decoupled/2p2c/2p2cproblem.hh
index c4b43b206b61a73c8e22007cf6fd358945b24f40..da6872f495fd2a1991dabad27f9dcdc2d44cbec6 100644
--- a/dumux/decoupled/2p2c/2p2cproblem.hh
+++ b/dumux/decoupled/2p2c/2p2cproblem.hh
@@ -25,8 +25,7 @@
#define DUMUX_IMPETPROBLEM_2P2C_HH
#include "boundaryconditions2p2c.hh"
-#include <dumux/decoupled/common/impet.hh>
-#include <dumux/decoupled/common/impetproblem.hh>
+#include <dumux/decoupled/2p/impes/impesproblem2p.hh>
#include <dumux/decoupled/2p2c/variableclass2p2c.hh>
#include <dumux/decoupled/2p2c/2p2cproperties.hh>
@@ -45,9 +44,9 @@ namespace Dumux
* in the specific problem.
*/
template<class TypeTag>
-class IMPETProblem2P2C : public IMPETProblem<TypeTag>
+class IMPETProblem2P2C : public IMPESProblem2P<TypeTag>
{
- typedef IMPETProblem<TypeTag> ParentType;
+ typedef IMPESProblem2P<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
@@ -59,7 +58,6 @@ class IMPETProblem2P2C : public IMPETProblem<TypeTag>
typedef typename GridView::Traits::template Codim<0>::Entity Element;
// material properties
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
@@ -78,16 +76,8 @@ public:
* \param verbose Output flag for the time manager.
*/
IMPETProblem2P2C(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView),
- gravity_(0)
- {
- newSpatialParams_ = true;
- spatialParameters_ = new SpatialParameters(gridView);
-
- gravity_ = 0;
- if (GET_PARAM(TypeTag, bool, EnableGravity))
- gravity_[dim - 1] = - 9.81;
- }
+ : ParentType(timeManager, gridView)
+ { }
/*!
* \brief The constructor
*
@@ -96,14 +86,8 @@ public:
* \param verbose Output flag for the time manager.
*/
IMPETProblem2P2C(TimeManager &timeManager, const GridView &gridView, SpatialParameters &spatialParameters)
- : ParentType(timeManager, gridView),
- gravity_(0),spatialParameters_(&spatialParameters)
- {
- newSpatialParams_ = false;
- gravity_ = 0;
- if (GET_PARAM(TypeTag, bool, EnableGravity))
- gravity_[dim - 1] = - 9.81;
- }
+ : ParentType(timeManager, gridView, spatialParameters)
+ { }
/*!
* \brief The constructor
@@ -113,16 +97,8 @@ public:
*/
IMPETProblem2P2C(const GridView &gridView, bool verbose = true)
DUNE_DEPRECATED // use IMPETProblem2P2C(TimeManager&, const GridView &)
- : ParentType(gridView, verbose),
- gravity_(0)
- {
- newSpatialParams_ = true;
- spatialParameters_ = new SpatialParameters(gridView);
-
- gravity_ = 0;
- if (GET_PARAM(TypeTag, bool, EnableGravity))
- gravity_[dim - 1] = - 9.81;
- }
+ : ParentType(gridView, verbose)
+ { }
/*!
* \brief The constructor
*
@@ -132,38 +108,57 @@ public:
*/
IMPETProblem2P2C(const GridView &gridView, SpatialParameters &spatialParameters, bool verbose = true)
DUNE_DEPRECATED // use IMPETProblem2P2C(TimeManager&, const GridView &)
- : ParentType(gridView, verbose),
- gravity_(0),spatialParameters_(&spatialParameters)
- {
- newSpatialParams_ = false;
- gravity_ = 0;
- if (GET_PARAM(TypeTag, bool, EnableGravity))
- gravity_[dim - 1] = - 9.81;
- }
+ : ParentType(gridView, spatialParameters, verbose)
+ { }
virtual ~IMPETProblem2P2C()
- {
- if (newSpatialParams_)
- {
- delete spatialParameters_;
- }
- }
+ { }
/*!
* \name Problem parameters
*/
// \{
-
- /*!
- * \copydoc Dumux::IMPESProblem2P::temperature()
+ //! Saturation initial condition (dimensionless)
+ /*! The problem is initialized with the following saturation. Both
+ * phases are assumed to contain an equilibrium concentration of the
+ * correspondingly other component.
+ * \param element The element.
*/
- Scalar temperature() const
- { return this->asImp_()->temperature(); };
+ Scalar initSat(const Element& element) const
+ {
+ return asImp_().initSatAtPos(element.geometry().center());
+ }
+ //! Saturation initial condition (dimensionless) at given position
+ /*! Has to be provided if initSat() is not used in the specific problem.
+ * \param globalPos The global position.
+ */
+ Scalar initSatAtPos(const GlobalPosition& globalPos) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "please specify initial saturation in the problem"
+ " using an initSatAtPos() method!");
+ return NAN;
+ }
+ //! Concentration initial condition (dimensionless)
+ /*! The problem is initialized with the following concentration.
+ */
+ Scalar initConcentration(const Element& element) const
+ {
+ return asImp_().initConcentrationAtPos(element.geometry().center());
+ }
+ //! Concentration initial condition (dimensionless)
+ /*! Has to be provided if initConcentration() is not used in the specific problem.
+ * \param globalPos The global position.
+ */
+ Scalar initConcentrationAtPos(const GlobalPosition& globalPos) const
+ {
+ DUNE_THROW(Dune::NotImplemented, "please specify initial Concentration in the problem"
+ " using an initConcentrationAtPos() method!");
+ return NAN;
+ }
+ // \}
/*!
- * \copydoc Dumux::IMPESProblem2P::gravity()
+ * \name Deprecated Problem parameters
*/
- const GlobalPosition &gravity() const
- { return gravity_; }
//! Saturation initial condition (dimensionless)
/*! The problem is initialized with the following saturation. Both
@@ -171,6 +166,7 @@ public:
* correspondingly other component.
*/
Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+ DUNE_DEPRECATED // use initSat(const Element& element)
{
DUNE_THROW(Dune::NotImplemented, "please specify initial saturation in the problem!");
return NAN;
@@ -179,30 +175,13 @@ public:
/*! The problem is initialized with the following concentration.
*/
Scalar initConcentration(const GlobalPosition& globalPos, const Element& element) const
+ DUNE_DEPRECATED // use initConcentration(const Element& element)
{
DUNE_THROW(Dune::NotImplemented, "please specify initial Concentration in the problem!");
return NAN;
}
// \}
- /*!
- * \name Access functions
- */
- // \{
- /*!
- * \copydoc Dumux::IMPESProblem2P::spatialParameters()
- */
- SpatialParameters &spatialParameters()
- { return *spatialParameters_; }
-
- /*!
- * \copydoc Dumux::IMPESProblem2P::spatialParameters()
- */
- const SpatialParameters &spatialParameters() const
- { return *spatialParameters_; }
-
- // \}
-
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
@@ -211,12 +190,6 @@ private:
//! \copydoc Dumux::IMPETProblem::asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
-
- GlobalPosition gravity_;
-
- // fluids and material properties
- SpatialParameters* spatialParameters_;
- bool newSpatialParams_;
};
}
diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index 65b44a757401f41274ba73c4f74d1cb2acef6d1e..10ec927fd289fa8db3bd7e389cb6f5e52937be5d 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -85,7 +85,7 @@ NEW_PROP_TAG( EnableMultiPointFluxApproximationOnAdaptiveGrids ); // Two-point f
//////////////////////////////////////////////////////////////////
SET_TYPE_PROP(DecoupledTwoPTwoC, TwoPTwoCIndices,DecoupledTwoPTwoCIndices<TypeTag>);
-SET_INT_PROP(DecoupledTwoPTwoC, NumEq, 2);
+SET_INT_PROP(DecoupledTwoPTwoC, NumEq, 3);
// set fluid/component information
SET_PROP(DecoupledTwoPTwoC, NumPhases) //!< The number of phases in the 2p model is 2
@@ -167,8 +167,22 @@ private:
public:
// Component indices
- static const int wCompIdx = FluidSystem::wCompIdx; //!< Index of the wetting phase in a phase vector
- static const int nCompIdx = FluidSystem::nCompIdx; //!< Index of the non-wetting phase in a phase vector
+ static const int wCompIdx = FluidSystem::wPhaseIdx; //!< Component index equals phase index
+ static const int nCompIdx = FluidSystem::nPhaseIdx; //!< Component index equals phase index
+
+ // Equation indices
+ static const int pressureEqIdx = 0;
+ static const int transportEqOffset = pressureEqIdx + 1; //!< Offset to access transport (mass conservation -) equations
+ static const int contiWEqIdx = transportEqOffset + wCompIdx; //!< Index of the wetting component transport equation
+ static const int contiNEqIdx = transportEqOffset + nCompIdx; //!< Index of the nonwetting component transport equation
+
+ //! Type of value on the Boundary
+ enum BoundaryFormulation
+ {
+ saturation=-1,
+ concentration=-2,
+ };
+
// BoundaryCondition flags
static const int satDependent = 0;
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 4573e546f43d10e8476dbb93d752d41a8d91b7be..0477fe927c407f2385918184938879416d638cbb 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -74,6 +74,7 @@ template<class TypeTag> class FVPressure2P2C
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
@@ -93,7 +94,8 @@ template<class TypeTag> class FVPressure2P2C
enum
{
wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
- wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx
+ wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
+ contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
};
// typedefs to abbreviate several dune classes...
@@ -104,10 +106,10 @@ template<class TypeTag> class FVPressure2P2C
typedef typename GridView::IntersectionIterator IntersectionIterator;
// convenience shortcuts for Vectors/Matrices
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
// the typenames used for the stiffness matrix and solution vector
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
@@ -493,11 +495,12 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
/****************implement sources and sinks************************/
- Dune::FieldVector<Scalar,2> source(problem_.source(globalPos, *eIt));
+ PrimaryVariables source(NAN);
+ problem_.source(source, *eIt);
if(first)
{
- source[wPhaseIdx] /= densityWI;
- source[nPhaseIdx] /= densityNWI;
+ source[contiWEqIdx] /= densityWI;
+ source[contiNEqIdx] /= densityNWI;
}
else
{
@@ -508,10 +511,10 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
problem_.variables().dv(globalIdxI, nPhaseIdx),
problem_.variables().dv_dp(globalIdxI));
- source[wPhaseIdx] *= problem_.variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[nPhaseIdx] *= problem_.variables().dv(globalIdxI, nPhaseIdx);
+ source[contiWEqIdx] *= problem_.variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[contiNEqIdx] *= problem_.variables().dv(globalIdxI, nPhaseIdx);
}
- f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ f_[globalIdxI] = volume * (source[contiWEqIdx] + source[contiNEqIdx]);
/***********************************/
// get absolute permeability
@@ -760,14 +763,15 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
unitDistVec /= dist;
//get boundary condition for boundary face center
- BoundaryConditions::Flags bctype = problem_.bcTypePress(globalPosFace, *isIt);
+ BoundaryTypes bcType;
+ problem_.boundaryTypes(bcType, *isIt);
// prepare pressure boundary condition
PhaseVector pressBC(0.);
Scalar pcBound (0.);
/********** Dirichlet Boundary *************/
- if (bctype == BoundaryConditions::dirichlet)
+ if (bcType.isDirichlet(Indices::pressureEqIdx))
{
//permeability vector at boundary
Dune::FieldVector<Scalar, dim> permeability(0);
@@ -776,13 +780,16 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
// create a fluid state for the boundary
FluidState BCfluidState;
- Scalar temperatureBC = problem_.temperature(globalPosFace, *eIt);
+ Scalar temperatureBC = problem_.temperatureAtPos(globalPosFace);
+ //read boundary values
+ PrimaryVariables primaryVariablesOnBoundary(NAN);
+ problem_.dirichlet(primaryVariablesOnBoundary, *isIt);
if (first)
{
Scalar lambda = lambdaWI+lambdaNWI;
A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
- Scalar pressBC = problem_.dirichletPress(globalPosFace, *isIt);
+ Scalar pressBC = primaryVariablesOnBoundary[Indices::pressureEqIdx];
f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
Scalar rightentry = (fractionalWI * densityWI
+ fractionalNWI * densityNWI)
@@ -863,19 +870,19 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
case pw:
{
potentialW = (unitOuterNormal * distVec) *
- (problem_.variables().pressure()[globalIdxI] - problem_.dirichletPress(globalPosFace, *isIt)) / (dist * dist);
+ (problem_.variables().pressure()[globalIdxI] - pressBC[wPhaseIdx]) / (dist * dist);
potentialNW = (unitOuterNormal * distVec) *
(problem_.variables().pressure()[globalIdxI] + pcI
- - problem_.dirichletPress(globalPosFace, *isIt) - pcBound) / (dist * dist);
+ - pressBC[wPhaseIdx] - pcBound) / (dist * dist);
break;
}
case pn:
{
potentialW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] - pcI -
- problem_.dirichletPress(globalPosFace, *isIt) + pcBound)
+ pressBC[nPhaseIdx] + pcBound)
/ (dist * dist);
potentialNW = (unitOuterNormal * distVec) * (problem_.variables().pressure()[globalIdxI] -
- problem_.dirichletPress(globalPosFace, *isIt)) / (dist * dist);
+ pressBC[nPhaseIdx]) / (dist * dist);
break;
}
}
@@ -956,7 +963,7 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
// set diagonal entry and right hand side entry
A_[globalIdxI][globalIdxI] += entry;
- f_[globalIdxI] += entry * problem_.dirichletPress(globalPosFace, *isIt);
+ f_[globalIdxI] += entry * primaryVariablesOnBoundary[Indices::pressureEqIdx];
f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
} //end of if(first) ... else{...
} // end dirichlet
@@ -964,22 +971,25 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
/**********************************
* set neumann boundary condition
**********************************/
- else
+ else if(bcType.isNeumann(Indices::pressureEqIdx))
{
- Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
+ PrimaryVariables J(NAN);
+ problem_.neumann(J, *isIt);
if (first)
{
- J[wPhaseIdx] /= densityWI;
- J[nPhaseIdx] /= densityNWI;
+ J[contiWEqIdx] /= densityWI;
+ J[contiNEqIdx] /= densityNWI;
}
else
{
- J[wPhaseIdx] *= dv_dC1;
- J[nPhaseIdx] *= dv_dC2;
+ J[contiWEqIdx] *= dv_dC1;
+ J[contiNEqIdx] *= dv_dC2;
}
- f_[globalIdxI] -= (J[wPhaseIdx] + J[nPhaseIdx]) * faceArea;
+ f_[globalIdxI] -= (J[contiWEqIdx] + J[contiNEqIdx]) * faceArea;
}
+ else
+ DUNE_THROW(Dune::NotImplemented, "Boundary Condition neither Dirichlet nor Neumann!");
}
} // end all intersections
// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
@@ -1080,26 +1090,27 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
int globalIdx = problem_.variables().index(*eIt);
// get the temperature
- Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+ Scalar temperature_ = problem_.temperatureAtPos(globalPos);
// initial conditions
PhaseVector pressure(0.);
Scalar sat_0=0.;
- BoundaryConditions2p2c::Flags ictype = problem_.initFormulation(globalPos, *eIt); // get type of initial condition
+ typename Indices::BoundaryFormulation icFormulation;
+ problem_.initialFormulation(icFormulation, *eIt); // get type of initial condition
if(!compositional) //means that we do the first approximate guess without compositions
{
// phase pressures are unknown, so start with an exemplary
- Scalar exemplaryPressure = problem_.referencePressure(globalPos, *eIt);
+ Scalar exemplaryPressure = problem_.referencePressure(*eIt);
pressure[wPhaseIdx] = pressure[nPhaseIdx] = problem_.variables().pressure()[globalIdx] = exemplaryPressure;
problem_.variables().capillaryPressure(globalIdx) = 0.;
- if (ictype == BoundaryConditions2p2c::saturation) // saturation initial condition
+ if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
{
sat_0 = problem_.initSat(globalPos, *eIt);
fluidState.satFlash(sat_0, pressure, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
}
- else if (ictype == BoundaryConditions2p2c::concentration) // concentration initial condition
+ else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
{
Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
fluidState.update(Z1_0, pressure, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
@@ -1107,7 +1118,7 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
}
else if(compositional) //means we regard compositional effects since we know an estimate pressure field
{
- if (ictype == Dumux::BoundaryConditions2p2c::saturation) // saturation initial condition
+ if (icFormulation == Indices::BoundaryFormulation::saturation) // saturation initial condition
{
//get saturation, determine pc
sat_0 = problem_.initSat(globalPos, *eIt);
@@ -1139,7 +1150,7 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
fluidState.satFlash(sat_0, pressure, problem_.spatialParameters().porosity(globalPos, *eIt), temperature_);
}
- else if (ictype == Dumux::BoundaryConditions2p2c::concentration) // concentration initial condition
+ else if (icFormulation == Indices::BoundaryFormulation::concentration) // concentration initial condition
{
Scalar Z1_0 = problem_.initConcentration(globalPos, *eIt);
// If total concentrations are given at the boundary, saturation is unknown.
@@ -1176,7 +1187,7 @@ void FVPressure2P2C<TypeTag>::initialMaterialLaws(bool compositional)
Scalar oldPc = pc;
//update with better pressures
fluidState.update(Z1_0, pressure, problem_.spatialParameters().porosity(globalPos, *eIt),
- problem_.temperature(globalPos, *eIt));
+ problem_.temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
@@ -1254,7 +1265,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
int globalIdx = problem_.variables().index(*eIt);
- Scalar temperature_ = problem_.temperature(globalPos, *eIt);
+ Scalar temperature_ = problem_.temperatureAtPos(globalPos);
// reset volume error
problem_.variables().volErr()[globalIdx] = 0;
@@ -1345,7 +1356,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
oldPc = pc;
//update with better pressures
fluidState.update(Z1, pressure, problem_.spatialParameters().porosity(globalPos, *eIt),
- problem_.temperature(globalPos, *eIt));
+ problem_.temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
@@ -1437,7 +1448,7 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
int globalIdx = problem_.variables().index(*ep);
// get cell temperature
- Scalar temperature_ = problem_.temperature(globalPos, *ep);
+ Scalar temperature_ = problem_.temperatureAtPos(globalPos);
// initialize an Fluid state for the update
FluidState updFluidState;
@@ -1535,7 +1546,8 @@ void FVPressure2P2C<TypeTag>::volumeDerivatives(GlobalPosition globalPos, Elemen
//check routines if derivatives are meaningful
if (isnan(dv_dC1) || isinf(dv_dC1) || isnan(dv_dC2) || isinf(dv_dC2)|| isnan(dv_dp) || isinf(dv_dp))
{
- DUNE_THROW(Dune::MathError, "NAN/inf of dV_dm. If that happens in first timestep, try smaller firstDt!");
+ std::cout << "blubbbber";
+// DUNE_THROW(Dune::MathError, "NAN/inf of dV_dm. If that happens in first timestep, try smaller firstDt!");
}
}
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index 0a8b6f813cf2be3b723ea8f2c7e35e17e3afca7d..e8814a51465a5e553644f1ddda6b5a4a907c5f1e 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -104,10 +104,11 @@ class FVPressure2P2CMultiPhysics : public FVPressure2P2C<TypeTag>
typedef typename GridView::IntersectionIterator IntersectionIterator;
// convenience shortcuts for Vectors/Matrices
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
//! Access functions to the current problem object
Problem& problem()
@@ -248,11 +249,12 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
Scalar densityNWI = problem().variables().densityNonwetting(globalIdxI);
/****************implement sources and sinks************************/
- Dune::FieldVector<Scalar,2> source(problem().source(globalPos, *eIt));
+ PrimaryVariables source(NAN);
+ problem().source(source, *eIt);
if(first || problem().variables().subdomain(globalIdxI) == 1)
{
- source[wPhaseIdx] /= densityWI;
- source[nPhaseIdx] /= densityNWI;
+ source[Indices::contiWEqIdx] /= densityWI;
+ source[Indices::contiNEqIdx] /= densityNWI;
}
else
{
@@ -263,10 +265,10 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
problem().variables().dv(globalIdxI, nPhaseIdx),
problem().variables().dv_dp(globalIdxI));
- source[wPhaseIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
- source[nPhaseIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
+ source[Indices::contiWEqIdx] *= problem().variables().dv(globalIdxI, wPhaseIdx); // note: dV_[i][1] = dv_dC1 = dV/dm1
+ source[Indices::contiNEqIdx] *= problem().variables().dv(globalIdxI, nPhaseIdx);
}
- this->f_[globalIdxI] = volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ this->f_[globalIdxI] = volume * (source[Indices::contiWEqIdx] + source[Indices::contiNEqIdx]);
/********************************************************************/
// get absolute permeability
@@ -584,14 +586,15 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
unitDistVec /= dist;
//get boundary condition for boundary face center
- BoundaryConditions::Flags bctype = problem().bcTypePress(globalPosFace, *isIt);
+ typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) bcType;
+ problem().boundaryTypes(bcType, *isIt);
// prepare pressure boundary condition
PhaseVector pressBC(0.);
Scalar pcBound (0.);
/********** Dirichlet Boundary *************/
- if (bctype == BoundaryConditions::dirichlet)
+ if (bcType.isDirichlet(Indices::pressureEqIdx))
{
//permeability vector at boundary
Dune::FieldVector<Scalar, dim> permeability(0);
@@ -600,13 +603,16 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
// create a fluid state for the boundary
FluidState BCfluidState;
- Scalar temperatureBC = problem().temperature(globalPosFace, *eIt);
+ //read boundary values
+ PrimaryVariables primaryVariablesOnBoundary(NAN);
+ problem().dirichlet(primaryVariablesOnBoundary, *isIt);
+ Scalar temperatureBC = problem().temperatureAtPos(globalPosFace);
if (first)
{
Scalar lambda = lambdaWI+lambdaNWI;
this->A_[globalIdxI][globalIdxI] += lambda * faceArea * (permeability * unitOuterNormal) / (dist);
- Scalar pressBC = problem().dirichletPress(globalPosFace, *isIt);
+ Scalar pressBC = primaryVariablesOnBoundary[Indices::pressureEqIdx];
this->f_[globalIdxI] += lambda * faceArea * pressBC * (permeability * unitOuterNormal) / (dist);
rightEntry = (fractionalWI * densityWI
+ fractionalNWI * densityNWI)
@@ -655,6 +661,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
lambdaNWBound = MaterialLaw::krn(
problem().spatialParameters().materialLawParams(globalPos, *eIt), BCfluidState.saturation(wPhaseIdx))
/ viscosityNWBound;
+ break;
}
}
@@ -681,34 +688,19 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
case pw:
{
potentialW = (unitOuterNormal * distVec) *
- (problem().variables().pressure()[globalIdxI] - problem().dirichletPress(globalPosFace, *isIt)) / (dist * dist);
+ (problem().variables().pressure()[globalIdxI] - pressBC[wPhaseIdx]) / (dist * dist);
potentialNW = (unitOuterNormal * distVec) *
(problem().variables().pressure()[globalIdxI] + pcI
- - problem().dirichletPress(globalPosFace, *isIt) - pcBound) / (dist * dist);
+ - pressBC[wPhaseIdx] - pcBound) / (dist * dist);
break;
}
case pn:
{
potentialW = (unitOuterNormal * distVec) * (problem().variables().pressure()[globalIdxI] - pcI -
- problem().dirichletPress(globalPosFace, *isIt) + pcBound)
+ pressBC[nPhaseIdx] + pcBound)
/ (dist * dist);
potentialNW = (unitOuterNormal * distVec) * (problem().variables().pressure()[globalIdxI] -
- problem().dirichletPress(globalPosFace, *isIt)) / (dist * dist);
- break;
- }
- case pglobal:
- {
- Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNWBound);
- Scalar fractionalNWBound = lambdaNWBound / (lambdaWBound + lambdaNWBound);
- Scalar fMeanW = 0.5 * (fractionalWI + fractionalWBound);
- Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWBound);
-
- potentialW = (unitOuterNormal * distVec)
- * (problem().variables().pressure()[globalIdxI] - problem().dirichletPress(globalPosFace, *isIt)
- - fMeanNW * (pcI - pcBound)) / (dist * dist);
- potentialNW = (unitOuterNormal * distVec)
- * (problem().variables().pressure()[globalIdxI] - problem().dirichletPress(globalPosFace, *isIt)
- + fMeanW * (pcI - pcBound)) / (dist * dist);
+ pressBC[nPhaseIdx]) / (dist * dist);
break;
}
}
@@ -821,7 +813,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
// set diagonal entry and right hand side entry
this->A_[globalIdxI][globalIdxI] += entry;
- this->f_[globalIdxI] += entry * problem().dirichletPress(globalPosFace, *isIt);
+ this->f_[globalIdxI] += entry * primaryVariablesOnBoundary[Indices::pressureEqIdx];
this->f_[globalIdxI] -= rightEntry * (unitOuterNormal * unitDistVec);
} //end of if(first) ... else{...
} // end dirichlet
@@ -831,7 +823,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
**********************************/
else
{
- Dune::FieldVector<Scalar,2> J = problem().neumann(globalPosFace, *isIt);
+ PrimaryVariables J(NAN);
+ problem().neumann(J, *isIt);
if (first || problem().variables().subdomain(globalIdxI)==1)
{
J[wPhaseIdx] /= densityWI;
@@ -868,7 +861,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if(problem().variables().saturation(globalIdxI)==1) //only w-phase
{
Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx,
- problem().temperature(globalPos, *eIt),
+ problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
@@ -878,7 +871,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
{
p_ -= 2*incp;
Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx,
- problem().temperature(globalPos, *eIt),
+ problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
@@ -888,7 +881,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if(problem().variables().saturation(globalIdxI)==0.) //only nw-phase
{
Scalar v_n_ = 1. / FluidSystem::phaseDensity(nPhaseIdx,
- problem().temperature(globalPos, *eIt),
+ problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
@@ -898,7 +891,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
{
p_ -= 2*incp;
Scalar v_n_ = 1. / FluidSystem::phaseDensity(nPhaseIdx,
- problem().temperature(globalPos, *eIt),
+ problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
@@ -985,7 +978,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
int globalIdx = problem().variables().index(*eIt);
- Scalar temperature_ = problem().temperature(globalPos, *eIt);
+ Scalar temperature_ = problem().temperatureAtPos(globalPos);
// reset volErr
problem().variables().volErr()[globalIdx] = 0;
@@ -1086,7 +1079,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
oldPc = pc;
//update with better pressures
fluidState.update(Z1, pressure, problem().spatialParameters().porosity(globalPos, *eIt),
- problem().temperature(globalPos, *eIt));
+ problem().temperatureAtPos(globalPos));
pc = MaterialLaw::pC(problem().spatialParameters().materialLawParams(globalPos, *eIt),
fluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index 6c914ea6df841cd3c9bd3907b3ce82d82396503f..9f9be20e8b85433f204ddee769e5f1d31384cf1f 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -55,6 +55,7 @@ class FVTransport2P2C
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
@@ -88,9 +89,9 @@ class FVTransport2P2C
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
public:
virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes);
@@ -123,27 +124,7 @@ public:
*/
FVTransport2P2C(Problem& problem) :
problem_(problem), switchNormals(false)
- {
- const int velocityType = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
- if (GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility)) && velocityType == vt)
- {
- DUNE_THROW(Dune::NotImplemented,
- "Total velocity - global pressure - model cannot be used with compressible fluids!");
- }
- const int saturationType = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
- if (saturationType != Sw && saturationType != Sn)
- {
- DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
- }
- if (pressureType != pw && pressureType != pn && pressureType != pglobal)
- {
- DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
- }
- if (velocityType != vw && velocityType != vn && velocityType != vt)
- {
- DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
- }
- }
+ { }
~FVTransport2P2C()
{
@@ -393,10 +374,11 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
FluidState BCfluidState;
//get boundary type
- BoundaryConditions::Flags bcTypeTransport_ = problem_.bcTypeTransport(globalPosFace, *isIt);
+ BoundaryTypes bcTypes;
+ problem_.boundaryTypes(bcTypes, *isIt);
/********** Dirichlet Boundary *************/
- if (bcTypeTransport_ == BoundaryConditions::dirichlet)
+ if (bcTypes.isDirichlet(Indices::contiWEqIdx)) // if contiWEq is Dirichlet, so is contiNEq
{
//get dirichlet pressure boundary condition
PhaseVector pressBound(0.);
@@ -415,10 +397,10 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
Scalar densityWBound = BCfluidState.density(wPhaseIdx);
Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
+ problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
+ problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
pcBound = BCfluidState.capillaryPressure();
@@ -436,19 +418,19 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
case pw:
{
potentialW = (K * unitOuterNormal) *
- (pressI - problem_.dirichletPress(globalPosFace, *isIt)) / dist;
+ (pressI - pressBound[wPhaseIdx]) / dist;
potentialNW = (K * unitOuterNormal) *
- (pressI + pcI - problem_.dirichletPress(globalPosFace, *isIt) - pcBound)
+ (pressI + pcI - pressBound[wPhaseIdx] - pcBound)
/ dist;
break;
}
case pn:
{
potentialW = (K * unitOuterNormal) *
- (pressI - pcI - problem_.dirichletPress(globalPosFace, *isIt) + pcBound)
+ (pressI - pcI - pressBound[nPhaseIdx] + pcBound)
/ dist;
potentialNW = (K * unitOuterNormal) *
- (pressI - problem_.dirichletPress(globalPosFace, *isIt)) / dist;
+ (pressI - pressBound[nPhaseIdx]) / dist;
break;
}
}
@@ -505,13 +487,13 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
+ velocityJIn * (1. - Xn1Bound) * densityNWBound
- velocityIJn * (1. - Xn1_I) * densityNWI ;
}//end dirichlet boundary
-
- if (bcTypeTransport_ == BoundaryConditions::neumann)
+ else if (bcTypes.isDirichlet(Indices::contiWEqIdx))
{
// Convention: outflow => positive sign : has to be subtracted from update vec
- Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
- updFactor[wCompIdx] = - J[0] * faceArea / volume;
- updFactor[nCompIdx] = - J[1] * faceArea / volume;
+ PrimaryVariables J(NAN);
+ problem_.neumann(J, *isIt);
+ updFactor[wCompIdx] = - J[Indices::contiWEqIdx] * faceArea / volume;
+ updFactor[nCompIdx] = - J[Indices::contiNEqIdx] * faceArea / volume;
// for timestep control
#define cflIgnoresNeumann
@@ -548,12 +530,11 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt, TransportSolut
}// end all intersections
- /************************************
- * Handle source term
- ***********************************/
- Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
- updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
- updateVec[nCompIdx][globalIdxI] += q[nCompIdx];
+ /*********** Handle source term ***************/
+ PrimaryVariables q(NAN);
+ problem_.source(q, *eIt);
+ updateVec[wCompIdx][globalIdxI] += q[Indices::contiWEqIdx];
+ updateVec[nCompIdx][globalIdxI] += q[Indices::contiNEqIdx];
// account for porosity
sumfactorin = std::max(sumfactorin,sumfactorout)
@@ -591,10 +572,15 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
PhaseVector &pressBound)
{
// read boundary values
- BoundaryConditions2p2c::Flags bctype = problem_.bcFormulation(globalPosFace, *isIt);
- if (bctype == BoundaryConditions2p2c::saturation)
+ PrimaryVariables primaryVariablesOnBoundary(0.);
+ problem_.dirichlet(primaryVariablesOnBoundary, *isIt);
+
+ // read boundary type
+ typename Indices::BoundaryFormulation bcType;
+ problem_.boundaryFormulation(bcType, *isIt);
+ if (bcType == Indices::BoundaryFormulation::saturation)
{
- Scalar satBound = problem_.dirichletTransport(globalPosFace, *isIt);
+ Scalar satBound = primaryVariablesOnBoundary[Indices::contiWEqIdx];
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
{
Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPosFace, *eIt),
@@ -603,33 +589,33 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
{
case pw:
{
- pressBound[wPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
- pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt) + pcBound;
+ pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
+ pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] + pcBound;
break;
}
case pn:
{
- pressBound[wPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
- pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
+ pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] - pcBound;
+ pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
break;
}
}
}
else // capillarity neglected
- pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
+ pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
BCfluidState.satFlash(satBound, pressBound, problem_.spatialParameters().porosity(globalPosFace, *eIt),
- problem_.temperature(globalPosFace, *eIt));
+ problem_.temperatureAtPos(globalPosFace));
}
- if (bctype == BoundaryConditions2p2c::concentration)
+ else if (bcType == Indices::BoundaryFormulation::concentration)
{
// saturation and hence pc and hence corresponding pressure unknown
- pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
- Scalar Z1Bound = problem_.dirichletTransport(globalPosFace, *isIt);
+ pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
+ Scalar Z1Bound = primaryVariablesOnBoundary[Indices::contiWEqIdx];
BCfluidState.update(Z1Bound, pressBound, problem_.spatialParameters().porosity(globalPosFace, *eIt),
- problem_.temperature(globalPosFace, *eIt));
+ problem_.temperatureAtPos(globalPosFace));
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
{
@@ -644,14 +630,14 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
{
case pw:
{
- pressBound[wPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
- pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt) + pcBound;
+ pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
+ pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] + pcBound;
break;
}
case pn:
{
- pressBound[wPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt) - pcBound;
- pressBound[nPhaseIdx] = problem_.dirichletPress(globalPosFace, *isIt);
+ pressBound[wPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx] - pcBound;
+ pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
break;
}
}
@@ -660,7 +646,7 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
Scalar oldPc = pcBound;
//update with better pressures
BCfluidState.update(Z1Bound, pressBound, problem_.spatialParameters().porosity(globalPosFace, *eIt),
- problem_.temperature(globalPosFace, *eIt));
+ problem_.temperatureAtPos(globalPosFace));
pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPosFace, *eIt),
BCfluidState.saturation(wPhaseIdx));
// TODO: get right criterion, do output for evaluation
@@ -670,6 +656,8 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
}
}
}
+ else
+ DUNE_THROW(Dune::NotImplemented, "Boundary Formulation neither Concentration nor Saturation??");
}
}
diff --git a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
index 9bc15cd36130675f832b79118cd0d74375978fcc..6da49d1f7da3da5dd215b6fd558ea50a9f95033e 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
@@ -93,9 +93,9 @@ class FVTransport2P2CMultiPhysics : public FVTransport2P2C<TypeTag>
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, 2> PhaseVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
public:
virtual void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec, bool impes);
@@ -347,9 +347,10 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
FluidState BCfluidState;
//get boundary type
- BoundaryConditions::Flags bcTypeTransport_ = problem_.bcTypeTransport(globalPosFace, *isIt);
+ typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) bcTypes;
+ problem_.boundaryTypes(bcTypes, *isIt);
- if (bcTypeTransport_ == BoundaryConditions::dirichlet)
+ if (bcTypes.isDirichlet(Indices::contiWEqIdx)) // if contiWEq is Dirichlet, so is contiNEq
{
//get dirichlet pressure boundary condition
PhaseVector pressBound(0.);
@@ -368,10 +369,10 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
Scalar densityWBound = BCfluidState.density(wPhaseIdx);
Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
+ problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
- problem_.temperature(globalPosFace, *eIt),
+ problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
pcBound = BCfluidState.capillaryPressure();
@@ -390,19 +391,19 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
case pw:
{
potentialW = (K * unitOuterNormal) *
- (pressI - problem_.dirichletPress(globalPosFace, *isIt)) / dist;
+ (pressI - pressBound[wPhaseIdx]) / dist;
potentialNW = (K * unitOuterNormal) *
- (pressI + pcI - problem_.dirichletPress(globalPosFace, *isIt) - pcBound)
+ (pressI + pcI - pressBound[wPhaseIdx] - pcBound)
/ dist;
break;
}
case pn:
{
potentialW = (K * unitOuterNormal) *
- (pressI - pcI - problem_.dirichletPress(globalPosFace, *isIt) + pcBound)
+ (pressI - pcI - pressBound[nPhaseIdx] + pcBound)
/ dist;
potentialNW = (K * unitOuterNormal) *
- (pressI - problem_.dirichletPress(globalPosFace, *isIt)) / dist;
+ (pressI - pressBound[nPhaseIdx]) / dist;
break;
}
}
@@ -460,13 +461,13 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
+ velocityJIn * (1. - Xn1Bound) * densityNWBound
- velocityIJn * (1. - Xn1_I) * densityNWI ;
}//end dirichlet boundary
-
- if (bcTypeTransport_ == BoundaryConditions::neumann)
+ else if (bcTypes.isDirichlet(Indices::contiWEqIdx))
{
// Convention: outflow => positive sign : has to be subtracted from update vec
- Dune::FieldVector<Scalar,2> J = problem_.neumann(globalPosFace, *isIt);
- updFactor[wCompIdx] = - J[0] * faceArea / volume;
- updFactor[nCompIdx] = - J[1] * faceArea / volume;
+ PrimaryVariables J(NAN);
+ problem_.neumann(J, *isIt);
+ updFactor[wCompIdx] = - J[Indices::contiWEqIdx] * faceArea / volume;
+ updFactor[nCompIdx] = - J[Indices::contiNEqIdx] * faceArea / volume;
// for timestep control
#define cflIgnoresNeumann
@@ -508,9 +509,10 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
/************************************
* Handle source term
***********************************/
- Dune::FieldVector<double,2> q = problem_.source(globalPos, *eIt);
- updateVec[wCompIdx][globalIdxI] += q[wCompIdx];
- updateVec[nCompIdx][globalIdxI] += q[nCompIdx];
+ PrimaryVariables q(NAN);
+ problem_.source(q, *eIt);
+ updateVec[wCompIdx][globalIdxI] += q[Indices::contiWEqIdx];
+ updateVec[nCompIdx][globalIdxI] += q[Indices::contiNEqIdx];
// account for porosity
sumfactorin = std::max(sumfactorin,sumfactorout)
diff --git a/test/decoupled/2p2c/test_dec2p2c.cc b/test/decoupled/2p2c/test_dec2p2c.cc
index d849518d47c5500f0be71f7ff24b76310de25366..99fd5619a9f4872b1e94e76dd17f37ef68744834 100644
--- a/test/decoupled/2p2c/test_dec2p2c.cc
+++ b/test/decoupled/2p2c/test_dec2p2c.cc
@@ -71,7 +71,7 @@ int main(int argc, char** argv)
double restartTime = 0;
// deal with start parameters
double tEnd= 3e3;
- double firstDt = 1e3;
+ double firstDt = 200;
if (argc != 1)
{
// deal with the restart stuff
diff --git a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
index 1175941a3e60e696a7173165cc7328f898f7b022..cd2173dcb3b64b5a5d788fdf76efa65d5d62be03 100644
--- a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
+++ b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
@@ -100,7 +100,7 @@ public:
// entry pressures function
materialLawParams_.setEntryPC(0);
- materialLawParams_.setMaxPC(0);
+ materialLawParams_.setMaxPC(10000);
for(int i = 0; i < dim; i++)
{
diff --git a/test/decoupled/2p2c/test_dec2p2cproblem.hh b/test/decoupled/2p2c/test_dec2p2cproblem.hh
index 25c2af581ad2f2c69eb0fc7535329f16226c7a15..ad9c3ae36d1455b9487a35631b28bdeeac4cf15b 100644
--- a/test/decoupled/2p2c/test_dec2p2cproblem.hh
+++ b/test/decoupled/2p2c/test_dec2p2cproblem.hh
@@ -82,7 +82,6 @@ SET_PROP(TestDecTwoPTwoCProblem, PressureModel)
SET_INT_PROP(TestDecTwoPTwoCProblem, PressureFormulation,
GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::pressureNW);
-
// Select fluid system
SET_PROP(TestDecTwoPTwoCProblem, FluidSystem)
{
@@ -115,6 +114,7 @@ public:
// Enable gravity
SET_BOOL_PROP(TestDecTwoPTwoCProblem, EnableGravity, true);
+SET_BOOL_PROP(TestDecTwoPTwoCProblem, EnableCapillarity, true);
SET_INT_PROP(TestDecTwoPTwoCProblem,
BoundaryMobility,
GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::satDependent);
@@ -146,6 +146,10 @@ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+// boundary typedefs
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
@@ -193,8 +197,9 @@ bool shouldWriteRestartFile() const
//! Returns the temperature within the domain.
/*! This problem assumes a temperature of 10 degrees Celsius.
+ * \param globalPos The global Position
*/
-Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return 273.15 + 10; // -> 10°C
}
@@ -204,94 +209,110 @@ Scalar temperature(const GlobalPosition& globalPos, const Element& element) cons
/*This pressure is used in order to calculate the material properties
* at the beginning of the initialization routine. It should lie within
* a reasonable pressure range for the current problem.
+ * \param globalPos The global Position
*/
-Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
{
return 1e6;
}
-//! Type of pressure boundary condition.
-/*! Defines the type the boundary condition for the pressure equation,
- * either pressure (dirichlet) or flux (neumann).
+//! Type of boundary condition.
+/*! Defines the type the boundary condition for the conservation equation,
+ * e.g. Dirichlet or Neumann. The Pressure equation is acessed via
+ * Indices::pressureEqIdx, while the transport (or mass conservation -)
+ * equations are reached with Indices::contiWEqIdx and Indices::contiNEqIdx
+ *
+ * \param bcTypes The boundary types for the conservation equations
+ * \param globalPos The global Position
*/
-typename BoundaryConditions::Flags bcTypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
{
if (globalPos[0] > 10-1E-6 || globalPos[0] < 1e-6)
- return BoundaryConditions::dirichlet;
- // all other boundaries
- return BoundaryConditions::neumann;
-}
-//! Type of Transport boundary condition.
-/*! Defines the type the boundary condition for the transport equation,
- * either saturation (dirichlet) or flux (neumann).
- */
-typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
-{
- return bcTypePress(globalPos, intersection);
+ bcTypes.setAllDirichlet();
+ else
+ // all other boundaries
+ bcTypes.setAllNeumann();
}
+
//! Flag for the type of Dirichlet conditions
/*! The Dirichlet BCs can be specified by a given concentration (mass of
* a component per total mass inside the control volume) or by means
* of a saturation.
+ *
+ * \param bcFormulation The boundary formulation for the conservation equations.
+ * \param intersection The intersection on the boundary.
*/
-BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
+const void boundaryFormulation(typename Indices::BoundaryFormulation &bcFormulation, const Intersection& intersection) const
{
- return BoundaryConditions2p2c::concentration;
+ bcFormulation = Indices::BoundaryFormulation::concentration;
}
-//! Value for dirichlet pressure boundary condition \f$ [Pa] \f$.
-/*! In case of a dirichlet BC for the pressure equation, the pressure
- * have to be defined on boundaries.
+//! Values for dirichlet boundary condition \f$ [Pa] \f$ for pressure and \f$ \frac{mass}{totalmass} \f$ or \f$ S_{\alpha} \f$ for transport.
+/*! In case of a dirichlet BC, values for all primary variables have to be set. In the sequential 2p2c model, a pressure
+ * is required for the pressure equation and component fractions for the transport equations. Although one BC for the two
+ * transport equations can be deduced by the other, it is seperately defined for consistency reasons with other models.
+ * Depending on the boundary Formulation, either saturation or total mass fractions can be defined.
+ *
+ * \param bcValues Vector holding values for all equations (pressure and transport).
+ * \param globalPos The global Position
*/
-Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+void dirichletAtPos(PrimaryVariables &bcValues ,const GlobalPosition& globalPos) const
{
- const Element& element = *(intersection.inside());
+ Scalar pRef = referencePressureAtPos(globalPos);
+ Scalar temp = temperatureAtPos(globalPos);
- Scalar pRef = referencePressure(globalPos, element);
- Scalar temp = temperature(globalPos, element);
-
- // all other boundaries
- return (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
+ // Dirichlet for pressure equation
+ bcValues[Indices::pressureEqIdx] = (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
: (2e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1]);
+
+ // Dirichlet values for transport equations
+ bcValues[Indices::contiWEqIdx] = 1.;
+ bcValues[Indices::contiNEqIdx] = 1.- bcValues[Indices::contiWEqIdx];
+
}
-//! Value for transport dirichlet boundary condition (dimensionless).
-/*! In case of a dirichlet BC for the transport equation,
- * has to be defined on boundaries.
- */
-Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
-{
- return 1.;
-}
-//! Value for pressure neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] \f$.
+
+//! Value for neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] \f$.
/*! In case of a neumann boundary condition, the flux of matter
- * is returned as a vector.
+ * is returned. Both pressure and transport module regard the neumann-bc values
+ * with the tranport (mass conservation -) equation indices. So the first entry
+ * of the vector is superflous.
+ * An influx into the domain has negative sign.
+ *
+ * \param neumannValues Vector holding values for all equations (pressure and transport).
+ * \param globalPos The global Position
*/
-Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+void neumannAtPos(PrimaryVariables &neumannValues, const GlobalPosition& globalPos) const
{
- Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+ neumannValues[Indices::contiNEqIdx] =0.;
+ neumannValues[Indices::contiWEqIdx] =0.;
// if (globalPos[1] < 15 && globalPos[1]> 5)
// {
-// neumannFlux[nPhaseIdx] = -0.015;
+// neumannValues[Indices::contiNEqIdx] = -0.015;
// }
- return neumannFlux;
}
//! Source of mass \f$ [\frac{kg}{m^3 \cdot s}] \f$
/*! Evaluate the source term for all phases within a given
* volume. The method returns the mass generated (positive) or
* annihilated (negative) per volume unit.
+ * Both pressure and transport module regard the neumann-bc values
+ * with the tranport (mass conservation -) equation indices. So the first entry
+ * of the vector is superflous.
+ *
+ * \param sourceValues Vector holding values for all equations (pressure and transport).
+ * \param globalPos The global Position
*/
-Dune::FieldVector<Scalar,2> source(const GlobalPosition& globalPos, const Element& element)
+void sourceAtPos(PrimaryVariables &sourceValues, const GlobalPosition& globalPos) const
{
- Dune::FieldVector<Scalar,2> q_(0);
- if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1) q_[1] = 0.0001;
- return q_;
+ sourceValues[Indices::contiWEqIdx]=0.;
+ if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1)
+ sourceValues[Indices::contiNEqIdx] = 0.0001;
}
//! Flag for the type of initial conditions
/*! The problem can be initialized by a given concentration (mass of
* a component per total mass inside the control volume) or by means
* of a saturation.
*/
-const BoundaryConditions2p2c::Flags initFormulation (const GlobalPosition& globalPos, const Element& element) const
+const void initialFormulation(typename Indices::BoundaryFormulation &initialFormulation, const Element& element) const
{
- return BoundaryConditions2p2c::concentration;
+ initialFormulation = Indices::BoundaryFormulation::concentration;
}
//! Saturation initial condition (dimensionless)
/*! The problem is initialized with the following saturation. Both
diff --git a/test/decoupled/2p2c/test_multiphysics2p2c.cc b/test/decoupled/2p2c/test_multiphysics2p2c.cc
index fbb815fad8329911feb1d83584e7bd6b9d2ce726..3c047a685b8b26f74717c96ff334c11bf911b585 100644
--- a/test/decoupled/2p2c/test_multiphysics2p2c.cc
+++ b/test/decoupled/2p2c/test_multiphysics2p2c.cc
@@ -70,7 +70,7 @@ int main(int argc, char** argv)
double restartTime = 0;
// deal with start parameters
double tEnd= 3e3;
- double firstDt = 1e3;
+ double firstDt = 200;
if (argc != 1)
{
// deal with the restart stuff
diff --git a/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
index 6b2f71d635a67cdad648198bba2b00c7c7a37b93..b6ff52f51b084dc87804bad46aaeababe403c97b 100644
--- a/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
+++ b/test/decoupled/2p2c/test_multiphysics2p2cproblem.hh
@@ -79,11 +79,8 @@ SET_PROP(TestMultTwoPTwoCProblem, PressureModel)
typedef Dumux::FVPressure2P2CMultiPhysics<TTAG(TestMultTwoPTwoCProblem)> type;
};
-SET_INT_PROP(TestMultTwoPTwoCProblem, VelocityFormulation,
- GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::velocityW);
-
SET_INT_PROP(TestMultTwoPTwoCProblem, PressureFormulation,
- GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::pressureW);
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::pressureNW);
//// Select fluid system
//SET_PROP(TestMultTwoPTwoCProblem, FluidSystem)
@@ -117,6 +114,7 @@ public:
// Enable gravity
SET_BOOL_PROP(TestMultTwoPTwoCProblem, EnableGravity, true);
+SET_BOOL_PROP(TestMultTwoPTwoCProblem, EnableCapillarity, true);
SET_INT_PROP(DecoupledTwoPTwoC,
BoundaryMobility,
GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices))::satDependent);
@@ -136,7 +134,7 @@ SET_SCALAR_PROP(TestMultTwoPTwoCProblem, CFLFactor, 0.8);
* description in the pressure module)
*
* To run the simulation execute the following line in shell:
- * <tt>./test_dec2p2c</tt>
+ * <tt>./test_multiphyiscs2p2c</tt>
* Optionally, simulation endtime and first timestep size can be
* specified by programm arguments.
*/
@@ -151,6 +149,10 @@ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+// boundary typedefs
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
@@ -172,7 +174,9 @@ public:
TestMultTwoPTwoCProblem(TimeManager &timeManager, const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
ParentType(timeManager, gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
{
+ // Specifies how many time-steps are done before output will be written.
this->setOutputInterval(10);
+
// initialize the tables of the fluid system
// FluidSystem::init();
}
@@ -200,91 +204,82 @@ bool shouldWriteRestartFile() const
//! Returns the temperature within the domain.
/*! This problem assumes a temperature of 10 degrees Celsius.
*/
-Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return 273.15 + 10; // -> 10°C
}
// \}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::referencePressure()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::referencePressureAtPos()
*/
-Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
{
return 1e6;
}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::bcTypePress()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::boundaryTypesAtPos()
*/
-typename BoundaryConditions::Flags bcTypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
{
if (globalPos[0] > 10-1E-6 || globalPos[0] < 1e-6)
- return BoundaryConditions::dirichlet;
- // all other boundaries
- return BoundaryConditions::neumann;
-}
-/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::bcTypeTransport()
- */
-typename BoundaryConditions::Flags bcTypeTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
-{
- return bcTypePress(globalPos, intersection);
+ bcTypes.setAllDirichlet();
+ else
+ // all other boundaries
+ bcTypes.setAllNeumann();
}
+
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::bcFormulation()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::boundaryFormulation()
*/
-BoundaryConditions2p2c::Flags bcFormulation(const GlobalPosition& globalPos, const Intersection& intersection) const
+const void boundaryFormulation(typename Indices::BoundaryFormulation &bcFormulation, const Intersection& intersection) const
{
- return BoundaryConditions2p2c::concentration;
+ bcFormulation = Indices::BoundaryFormulation::concentration;
}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::dirichletPress()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::dirichletAtPos()
*/
-Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+void dirichletAtPos(PrimaryVariables &bcValues ,const GlobalPosition& globalPos) const
{
- const Element& element = *(intersection.inside());
+ Scalar pRef = referencePressureAtPos(globalPos);
+ Scalar temp = temperatureAtPos(globalPos);
- Scalar pRef = referencePressure(globalPos, element);
- Scalar temp = temperature(globalPos, element);
-
- // all other boundaries
- return (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
+ // Dirichlet for pressure equation
+ bcValues[Indices::pressureEqIdx] = (globalPos[0] < 1e-6) ? (2.5e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1])
: (2e5 - FluidSystem::H2O::liquidDensity(temp, pRef) * this->gravity()[dim-1]);
+
+ // Dirichlet values for transport equations
+ bcValues[Indices::contiWEqIdx] = 1.;
+ bcValues[Indices::contiNEqIdx] = 1.- bcValues[Indices::contiWEqIdx];
+
}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::dirichletTransport()
- */
-Scalar dirichletTransport(const GlobalPosition& globalPos, const Intersection& intersection) const
-{
- return 1.;
-}
-/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::neumann()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::neumannAtPos()
*/
-Dune::FieldVector<Scalar,2> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
+void neumannAtPos(PrimaryVariables &neumannValues, const GlobalPosition& globalPos) const
{
- Dune::FieldVector<Scalar,2> neumannFlux(0.0);
+ neumannValues[Indices::contiNEqIdx] = 0.;
+ neumannValues[Indices::contiWEqIdx] = 0.;
// if (globalPos[1] < 15 && globalPos[1]> 5)
// {
-// neumannFlux[nPhaseIdx] = -0.015;
+// neumannValues[Indices::contiNEqIdx] = -0.015;
// }
- return neumannFlux;
}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::source()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::sourceAtPos()
*/
-Dune::FieldVector<Scalar,2> source(const GlobalPosition& globalPos, const Element& element)
+void sourceAtPos(PrimaryVariables &sourceValues, const GlobalPosition& globalPos) const
{
- Dune::FieldVector<Scalar,2> q_(0);
- if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1) q_[1] = 0.0001;
- return q_;
+ sourceValues[Indices::contiWEqIdx]=0.;
+ if (fabs(globalPos[0] - 4.5) < 1 && fabs(globalPos[1] - 4.5) < 1)
+ sourceValues[Indices::contiNEqIdx] = 0.0001;
}
/*!
- * \copydoc Dumux::TestDecTwoPTwoCProblem::initFormulation()
+ * \copydoc Dumux::TestDecTwoPTwoCProblem::initialFormulation()
*/
-const BoundaryConditions2p2c::Flags initFormulation (const GlobalPosition& globalPos, const Element& element) const
+const void initialFormulation(typename Indices::BoundaryFormulation &initialFormulation, const Element& element) const
{
- return BoundaryConditions2p2c::concentration;
+ initialFormulation = Indices::BoundaryFormulation::concentration;
}
/*!
* \copydoc Dumux::TestDecTwoPTwoCProblem::initSat()