diff --git a/dumux/boxmodels/1p/1pfluxvariables.hh b/dumux/boxmodels/1p/1pfluxvariables.hh
index e797b7af56dafca192d3322a5e8c0ea13afcecce..064a860cc5a7def2e52787a6c52fc5a5e5c6ee66 100644
--- a/dumux/boxmodels/1p/1pfluxvariables.hh
+++ b/dumux/boxmodels/1p/1pfluxvariables.hh
@@ -23,6 +23,8 @@
#ifndef DUMUX_1P_FLUX_VARIABLES_HH
#define DUMUX_1P_FLUX_VARIABLES_HH
+#include "1pproperties.hh"
+
#include <dumux/common/math.hh>
namespace Dumux
@@ -61,6 +63,9 @@ class OnePFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePIndices)) Indices;
+ typedef Dune::FieldVector<Scalar, dim> Vector;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
+
public:
OnePFluxVariables(const Problem &problem,
const Element &element,
@@ -71,24 +76,53 @@ public:
{
scvfIdx_ = faceIdx;
- densityAtIP = 0;
- viscosityAtIP = 0;
- pressureGrad = Scalar(0);
-
calculateGradients_(problem, element, elemDat);
- calculateVelocities_(problem, element, elemDat);
+ calculateK_(problem, element);
};
+
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
+ /*!
+ * \brief Return the intrinsic permeability.
+ */
+ const Tensor &intrinsicPermeability() const
+ { return K_; }
+
+ /*!
+ * \brief Return the pressure potential gradient.
+ */
+ const Vector &potentialGrad() const
+ { return potentialGrad_; }
+
+ /*!
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the upstream control volume
+ * for a given phase.
+ */
+ int upstreamIdx(Scalar normalFlux) const
+ { return (normalFlux >= 0)?face().i:face().j; }
+
+ /*!
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the downstream control volume
+ * for a given phase.
+ */
+ int downstreamIdx(Scalar normalFlux) const
+ { return (normalFlux >= 0)?face().j:face().i; }
+
private:
void calculateGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
- // calculate gradients
- GlobalPosition tmp(0.0);
+ potentialGrad_ = 0.0;
+ Scalar densityAtIP = 0.0;
+
+ // calculate potential gradient
for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
@@ -97,72 +131,50 @@ private:
const LocalPosition &feGrad = face().grad[idx];
// the pressure gradient
- tmp = feGrad;
- tmp *= elemDat[idx].pressure;
- pressureGrad += tmp;
-
+ Vector tmp(feGrad);
+ tmp *= elemDat[idx].pressure();
+ potentialGrad_ += tmp;
+
// fluid density
densityAtIP +=
- elemDat[idx].density*face().shapeValue[idx];
- // fluid viscosity
- viscosityAtIP +=
- elemDat[idx].viscosity*face().shapeValue[idx];
+ elemDat[idx].density()*face().shapeValue[idx];
}
// correct the pressure gradients by the hydrostatic
// pressure due to gravity
- tmp = problem.gravity();
+ Vector tmp(problem.gravity());
tmp *= densityAtIP;
- pressureGrad -= tmp;
+ potentialGrad_ -= tmp;
}
- void calculateVelocities_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemDat)
+ void calculateK_(const Problem &problem,
+ const Element &element)
{
const SpatialParameters &spatialParams = problem.spatialParameters();
- typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
- Tensor K;
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().j));
-
- // temporary vector for the Darcy velocity
- GlobalPosition vDarcy;
- K.mv(pressureGrad, vDarcy); // vDarcy = K * grad p
- vDarcyNormal = vDarcy*face().normal;
-
- // set the upstream and downstream vertices
- upstreamIdx = face().i;
- downstreamIdx = face().j;
- if (vDarcyNormal > 0)
- std::swap(upstreamIdx, downstreamIdx);
+ spatialParams.meanK(K_,
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().i),
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().j));
}
-public:
+protected:
const FVElementGeometry &fvElemGeom_;
int scvfIdx_;
// gradients
- GlobalPosition pressureGrad;
-
- // density of the fluid at the integration point
- Scalar densityAtIP;
- // viscosity of the fluid at the integration point
- Scalar viscosityAtIP;
+ Vector potentialGrad_;
- // darcy velocity in direction of the face normal
- Scalar vDarcyNormal;
+ // intrinsic permeability
+ Tensor K_;
// local index of the upwind vertex
- int upstreamIdx;
+ int upstreamIdx_;
// local index of the downwind vertex
- int downstreamIdx;
+ int downstreamIdx_;
};
} // end namepace
diff --git a/dumux/boxmodels/1p/1pindices.hh b/dumux/boxmodels/1p/1pindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..67fb76670200915e258b13973df16b0db34889b8
--- /dev/null
+++ b/dumux/boxmodels/1p/1pindices.hh
@@ -0,0 +1,37 @@
+// $Id: 1pproperties.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Indices for the single phase model.
+ */
+#ifndef DUMUX_1P_INDICES_HH
+#define DUMUX_1P_INDICES_HH
+
+namespace Dumux
+{
+/*!
+ * \brief Indices for the single phase model.
+ */
+struct OnePIndices
+{
+ static const int pressureIdx = 0;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/1p/1plocalresidual.hh b/dumux/boxmodels/1p/1plocalresidual.hh
index 220be3cdda34b0209c0ef6bcbc75e362ae835359..ac7a86762a0d23be10a1cca8f4ef3837d9e3a400 100644
--- a/dumux/boxmodels/1p/1plocalresidual.hh
+++ b/dumux/boxmodels/1p/1plocalresidual.hh
@@ -62,13 +62,13 @@ class OnePLocalResidual : public BoxLocalResidual<TypeTag>
pressureIdx = Indices::pressureIdx,
};
+ static const Scalar upwindWeight = GET_PROP_VALUE(TypeTag, PTAG(UpwindWeight));
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
public:
@@ -92,7 +92,7 @@ public:
const VolumeVariables &volVars = elemVars[scvIdx];
// partial time derivative of the wetting phase mass
- result[pressureIdx] = volVars.density * volVars.porosity;
+ result[pressureIdx] = volVars.density() * volVars.porosity();
}
@@ -102,13 +102,25 @@ public:
*/
void computeFlux(PrimaryVariables &flux, int faceId) const
{
- FluxVariables vars(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
- faceId,
- this->curVolVars_());
-
- flux[pressureIdx] = vars.densityAtIP * vars.vDarcyNormal / vars.viscosityAtIP;
+ FluxVariables fluxVars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceId,
+ this->curVolVars_());
+
+ Vector tmpVec;
+ fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(),
+ tmpVec);
+ Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
+
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
+ flux[pressureIdx] =
+ (( upwindWeight)*(up.density()/up.viscosity())
+ +
+ (1 - upwindWeight)*(dn.density()/dn.viscosity()))
+ *
+ normalFlux;
}
/*!
diff --git a/dumux/boxmodels/1p/1pmodel.hh b/dumux/boxmodels/1p/1pmodel.hh
index 7ff7847337243b75c8137d4bb530ecb44f241f5c..5519535448774964db35f8ebb8b08504e45dfa72 100644
--- a/dumux/boxmodels/1p/1pmodel.hh
+++ b/dumux/boxmodels/1p/1pmodel.hh
@@ -118,7 +118,7 @@ public:
i,
false);
- (*p)[globalIdx] = volVars.pressure;
+ (*p)[globalIdx] = volVars.pressure();
};
}
@@ -128,4 +128,6 @@ public:
};
}
+#include "1ppropertydefaults.hh"
+
#endif
diff --git a/dumux/boxmodels/1p/1pproperties.hh b/dumux/boxmodels/1p/1pproperties.hh
index 29936e46a6042b541f37e45b8293e3f4c599fd41..bc6aa89b7e50ef0d03cf7817e75369d998e77eb0 100644
--- a/dumux/boxmodels/1p/1pproperties.hh
+++ b/dumux/boxmodels/1p/1pproperties.hh
@@ -30,29 +30,6 @@ namespace Dumux
* \addtogroup OnePBoxModel
*/
// \{
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-template<class TypeTag>
-class OnePBoxModel;
-
-template<class TypeTag>
-class OnePLocalResidual;
-
-template <class TypeTag>
-class OnePVolumeVariables;
-
-template <class TypeTag>
-class OnePFluxVariables;
-
-/*!
- * \brief Indices for the single phase model.
- */
-struct OnePIndices
-{
- static const int pressureIdx = 0;
-};
-
///////////////////////////////////////////////////////////////////////////
// properties for the isothermal single phase model
///////////////////////////////////////////////////////////////////////////
@@ -74,30 +51,7 @@ NEW_PROP_TAG(OnePIndices); //!< Enumerations for the 1p models
NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-SET_INT_PROP(BoxOneP, NumEq, 1);
-SET_INT_PROP(BoxOneP, NumPhases, 1);
-
-//! The local residual function
-SET_TYPE_PROP(BoxOneP,
- LocalResidual,
- OnePLocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(BoxOneP, Model, OnePBoxModel<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(BoxOneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
-
-//! the FluxVariables property
-SET_TYPE_PROP(BoxOneP, FluxVariables, OnePFluxVariables<TypeTag>);
-
-//! The indices required by the isothermal single-phase model
-SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
+NEW_PROP_TAG(UpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
// \}
};
diff --git a/dumux/boxmodels/1p/1ppropertydefaults.hh b/dumux/boxmodels/1p/1ppropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0e2499721ed78dce234436c0a82192f19c95588e
--- /dev/null
+++ b/dumux/boxmodels/1p/1ppropertydefaults.hh
@@ -0,0 +1,68 @@
+// $Id: 1pproperties.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the onephase BOX model.
+ */
+#ifndef DUMUX_1P_PROPERTY_DEFAULTS_HH
+#define DUMUX_1P_PROPERTY_DEFAULTS_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+#include "1pmodel.hh"
+#include "1plocalresidual.hh"
+#include "1pvolumevariables.hh"
+#include "1pfluxvariables.hh"
+#include "1pindices.hh"
+
+namespace Dumux
+{
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+SET_INT_PROP(BoxOneP, NumEq, 1);
+SET_INT_PROP(BoxOneP, NumPhases, 1);
+
+//! The local residual function
+SET_TYPE_PROP(BoxOneP,
+ LocalResidual,
+ OnePLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxOneP, Model, OnePBoxModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxOneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxOneP, FluxVariables, OnePFluxVariables<TypeTag>);
+
+//! The indices required by the isothermal single-phase model
+SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
+
+//! The weight of the upwind control volume when calculating
+//! fluxes. Use central differences by default.
+SET_SCALAR_PROP(BoxOneP, UpwindWeight, 0.5);
+
+// \}
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/1p/1pvolumevariables.hh b/dumux/boxmodels/1p/1pvolumevariables.hh
index 3e3b5edf6c4278614bec514ba75a3d7880f3a7be..5e5d4530e13728708063db317f5d90277ae86c8f 100644
--- a/dumux/boxmodels/1p/1pvolumevariables.hh
+++ b/dumux/boxmodels/1p/1pvolumevariables.hh
@@ -23,6 +23,7 @@
#define DUMUX_1P_VOLUME_VARIABLES_HH
#include "1pproperties.hh"
+#include <dumux/boxmodels/common/boxvolumevariables.hh>
namespace Dumux
{
@@ -33,8 +34,10 @@ namespace Dumux
* finite volume in the one-phase model.
*/
template <class TypeTag>
-class OnePVolumeVariables
+class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
{
+ typedef BoxVolumeVariables<TypeTag> ParentType;
+
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
@@ -71,40 +74,74 @@ public:
int scvIdx,
bool isOldSol)
{
- primaryVars_ = priVars;
+ ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
- typedef Indices I;
+ this->asImp_().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
- Scalar temperature = problem.temperature(element, elemGeom, scvIdx);
- pressure = priVars[I::pressureIdx];
- density = Fluid::density(temperature, pressure);
- viscosity = Fluid::viscosity(temperature, pressure);
+ pressure_ = priVars[Indices::pressureIdx];
+ density_ = Fluid::density(temperature(), pressure());
+ viscosity_ = Fluid::viscosity(temperature(), pressure());
// porosity
- porosity = problem.spatialParameters().porosity(element,
- elemGeom,
- scvIdx);
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
};
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperature of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ */
+ Scalar pressure() const
+ { return pressure_; }
/*!
- * \brief Sets the evaluation point used in the by the local jacobian.
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
*/
- void setEvalPoint(const Implementation *ep)
- { }
+ Scalar density() const
+ { return density_; }
/*!
- * \brief Return the vector of primary variables
+ * \brief Returns the dynamic viscosity of the fluid within the
+ * control volume.
*/
- const PrimaryVariables &primaryVars() const
- { return primaryVars_; }
+ Scalar viscosity() const
+ { return viscosity_; }
- Scalar pressure;
- Scalar density;
- Scalar viscosity;
- Scalar porosity;
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
protected:
- PrimaryVariables primaryVars_;
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ { temperature_ = problem.temperature(element, elemGeom, scvIdx); }
+
+ Scalar temperature_;
+ Scalar pressure_;
+ Scalar density_;
+ Scalar viscosity_;
+ Scalar porosity_;
};
}
diff --git a/dumux/boxmodels/1p2c/1p2cfluidstate.hh b/dumux/boxmodels/1p2c/1p2cfluidstate.hh
index a545fafc267938dd93963fc039c03308d99e0c56..81feb9b10d03deeaa0fb3717aadf5f354b275273 100644
--- a/dumux/boxmodels/1p2c/1p2cfluidstate.hh
+++ b/dumux/boxmodels/1p2c/1p2cfluidstate.hh
@@ -37,15 +37,20 @@ class OnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTA
OnePTwoCFluidState<TypeTag> >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
-
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum {
- konti = Indices::konti,
- transport = Indices::transport
+ pressureIdx = Indices::pressureIdx,
+ x1Idx = Indices::x1Idx,
+
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx,
+
+ phaseIndex = GET_PROP_VALUE(TypeTag, PTAG(PhaseIndex)),
+ comp1Index = GET_PROP_VALUE(TypeTag, PTAG(Comp1Index)),
+ comp2Index = GET_PROP_VALUE(TypeTag, PTAG(Comp2Index)),
};
public:
@@ -62,26 +67,106 @@ public:
temperature_ = temperature;
- phasePressure_[0] = primaryVars[konti];
-
- moleFraction_[0][konti] = 1.0 - primaryVars[transport];
- moleFraction_[0][transport] = primaryVars[transport];
+ phasePressure_ = primaryVars[pressureIdx];
+ x1_ = primaryVars[x1Idx];
+ meanMolarMass_ =
+ (1 - x1_)*FluidSystem::molarMass(comp1Index) +
+ (x1_ )*FluidSystem::molarMass(comp2Index);
+
+ density_ = FluidSystem::phaseDensity(phaseIndex, temperature_, phasePressure_, *this);
+
+ Valgrind::CheckDefined(x1_);
+ Valgrind::CheckDefined(phasePressure_);
+ Valgrind::CheckDefined(density_);
+ Valgrind::CheckDefined(meanMolarMass_);
+ Valgrind::CheckDefined(temperature_);
+ Valgrind::CheckDefined(*this);
}
-public:
+ /*!
+ * \brief Returns the saturation of a phase.
+ */
+ Scalar saturation(int phaseIdx) const
+ { return (phaseIndex == phaseIdx)?1.0:0.0; };
+
/*!
* \brief Returns the molar fraction of a component in a fluid phase.
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
+ {
+ // we are a single phase model!
+ if (phaseIdx != phaseIndex) return 0.0;
+
+ if (compIdx==comp1Index)
+ return 1-x1_;
+ else if (compIdx==comp2Index)
+ return x1_;
+ return 0.0;
+ }
+
+ /*!
+ * \brief Returns the total concentration of a phase [mol / m^3].
+ *
+ * This is equivalent to the sum of all component concentrations.
+ */
+ Scalar phaseConcentration(int phaseIdx) const
+ {
+ if (phaseIdx != phaseIndex)
+ return 0;
+ return density_/meanMolarMass_;
+ };
+
+ /*!
+ * \brief Returns the concentration of a component in a phase [mol / m^3].
+ */
+ Scalar concentration(int phaseIdx, int compIdx) const
+ { return phaseConcentration(phaseIdx)*moleFrac(phaseIdx, compIdx); };
+
+
+ /*!
+ * \brief Returns the mass fraction of a component in a phase.
+ */
+ Scalar massFrac(int phaseIdx, int compIdx) const
{
- return moleFraction_[phaseIdx][compIdx];
+ if (phaseIdx != phaseIndex)
+ return 0;
+ return
+ moleFrac(phaseIdx, compIdx)*
+ FluidSystem::molarMass(compIdx)
+ /
+ meanMolarMass_;
}
+ /*!
+ * \brief Returns the density of a phase [kg / m^3].
+ */
+ Scalar density(int phaseIdx) const
+ {
+ if (phaseIdx != phaseIndex)
+ return 0;
+ return density_;
+ }
+
+ /*!
+ * \brief Returns mean molar mass of a phase [kg / mol].
+ *
+ * This is equivalent to the sum of all component molar masses
+ * weighted by their respective mole fraction.
+ */
+ Scalar meanMolarMass(int phaseIdx) const
+ {
+ if (phaseIdx != phaseIndex)
+ return 0;
+ return meanMolarMass_;
+ };
+
/*!
* \brief Returns the pressure of a fluid phase [Pa].
*/
Scalar phasePressure(int phaseIdx) const
- { return phasePressure_[phaseIdx]; }
+ {
+ return phasePressure_;
+ }
/*!
* \brief Returns the temperature of the fluids [K].
@@ -93,8 +178,10 @@ public:
{ return temperature_; };
public:
- Scalar moleFraction_[numPhases][numComponents];
- Scalar phasePressure_[numPhases];
+ Scalar x1_;
+ Scalar phasePressure_;
+ Scalar density_;
+ Scalar meanMolarMass_;
Scalar temperature_;
};
diff --git a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
index f1ce7eb944cb75aa8d26ba9376e18000da6e8f56..52e67cccfdd227480cdb2f3cb4d480a5353def0c 100644
--- a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
@@ -27,6 +27,8 @@
#ifndef DUMUX_1P2C_FLUX_VARIABLES_HH
#define DUMUX_1P2C_FLUX_VARIABLES_HH
+#include "1p2cproperties.hh"
+
#include <dumux/common/math.hh>
namespace Dumux
@@ -45,243 +47,238 @@ class OnePTwoCFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
-
- typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
- enum {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
- };
+ enum { dim = GridView::dimension };
+ enum { dimWorld = GridView::dimensionworld };
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
+ typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolume SCV;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
- typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+ typedef typename GridView::template Codim<0>::Entity Element;
public:
OnePTwoCFluxVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int faceIdx,
- const ElementVolumeVariables &elemDat)
- : fvElemGeom(elemGeom)
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvfIdx,
+ const ElementVolumeVariables &elemDat)
+ : fvElemGeom_(elemGeom)
{
- face = &fvElemGeom.subContVolFace[faceIdx];
-
- int i = face->i;
- int j = face->j;
-
- insideSCV = &fvElemGeom.subContVol[i];
- outsideSCV = &fvElemGeom.subContVol[j];
-
- densityAtIP = 0;
- molarDensityAtIP = 0;
- viscosityAtIP = 0;
- pressureGrad = Scalar(0);
- concentrationGrad = Scalar(0);
+ scvfIdx_ = scvfIdx;
calculateGradients_(problem, element, elemDat);
- calculateVelocities_(problem, element, elemDat);
+ calculateK_(problem, element, elemDat);
calculateDiffCoeffPM_(problem, element, elemDat);
calculateDispersionTensor_(problem, element, elemDat);
};
-private:
+public:
+ const SCVFace &face() const
+ { return fvElemGeom_.subContVolFace[scvfIdx_]; }
+
+ /*!
+ * \brief Return the intrinsic permeability.
+ */
+ const Tensor &intrinsicPermeability() const
+ { return K_; }
+
+ /*!
+ * \brief Return the dispersion tensor.
+ */
+ const Tensor &dispersionTensor() const
+ { return dispersionTensor_; }
+
+ /*!
+ * \brief Return the pressure potential gradient.
+ */
+ const Vector &potentialGrad() const
+ { return potentialGrad_; }
+
+ /*!
+ * \brief Return the concentration gradient.
+ */
+ const Vector &concentrationGrad(int compIdx) const
+ {
+ if (compIdx != 1)
+ { DUNE_THROW(Dune::InvalidStateException,
+ "The 1p2c model is supposed to need "
+ "only the concentration gradient of "
+ "the second component!"); }
+ return concentrationGrad_;
+ };
+
+ Scalar porousDiffCoeff() const
+ {
+ // TODO: tensorial diffusion coefficients
+ return diffCoeffPM_;
+ };
+
+
+ /*!
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the upstream control volume
+ * for a given phase.
+ */
+ int upstreamIdx(Scalar normalFlux) const
+ { return (normalFlux >= 0)?face().i:face().j; }
+
+ /*!
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the downstream control volume
+ * for a given phase.
+ */
+ int downstreamIdx(Scalar normalFlux) const
+ { return (normalFlux > 0)?face().j:face().i; }
+
+protected:
void calculateGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
- GlobalPosition tmp;
- if (!problem.spatialParameters().useTwoPointGradient(element, face->i, face->j)) {
+ const VolumeVariables &vVars_i = elemDat[face().i];
+ const VolumeVariables &vVars_j = elemDat[face().j];
+
+ potentialGrad_ = 0.0;
+ concentrationGrad_ = 0.0;
+
+ Vector tmp;
+ if (!problem.spatialParameters().useTwoPointGradient(element, face().i, face().j)) {
// use finite-element gradients
tmp = 0.0;
- for (int idx = 0;
- idx < fvElemGeom.numVertices;
- idx++) // loop over adjacent vertices
+ int n = element.template count<dim>();
+ for (int idx = 0; idx < n; idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
- const LocalPosition &feGrad = face->grad[idx];
+ const Vector &feGrad = face().grad[idx];
// the pressure gradient
tmp = feGrad;
- tmp *= elemDat[idx].pressure;
- pressureGrad += tmp;
+ tmp *= elemDat[idx].pressure();
+ potentialGrad_ += tmp;
// the concentration gradient
tmp = feGrad;
- tmp *= elemDat[idx].molefraction;
- concentrationGrad += tmp;
-
- // phase density
- densityAtIP +=
- elemDat[idx].density*face->shapeValue[idx];
-
- // molar phase density
- molarDensityAtIP +=
- elemDat[idx].molarDensity*face->shapeValue[idx];
-
- // phase viscosity
- viscosityAtIP +=
- elemDat[idx].viscosity*face->shapeValue[idx];
+ tmp *= elemDat[idx].moleFrac(1);
+ concentrationGrad_ += tmp;
}
}
else {
// use two-point gradients
- tmp = element.geometry().corner(face->i);
- tmp -= element.geometry().corner(face->j);
+ tmp = element.geometry().corner(face().i);
+ tmp -= element.geometry().corner(face().j);
Scalar dist = tmp.two_norm();
- tmp = face->normal;
- tmp /= face->normal.two_norm()*dist;
+ tmp = face().normal;
+ tmp /= face().normal.two_norm()*dist;
- pressureGrad = tmp;
- pressureGrad *= elemDat[face->j].pressure - elemDat[face->i].pressure;
- concentrationGrad = tmp;
- concentrationGrad *= elemDat[face->j].molefraction - elemDat[face->i].molefraction;
- densityAtIP = (elemDat[face->j].density + elemDat[face->i].density)/2;
- molarDensityAtIP = (elemDat[face->j].molarDensity + elemDat[face->i].molarDensity)/2;
- viscosityAtIP = (elemDat[face->j].viscosity + elemDat[face->i].viscosity)/2;
+ potentialGrad_ = tmp;
+ potentialGrad_ *= vVars_j.pressure() - vVars_i.pressure();
+ concentrationGrad_ = tmp;
+ concentrationGrad_ *= vVars_j.moleFrac(1) - vVars_i.moleFrac(1);
}
// correct the pressure by the hydrostatic pressure due to
// gravity
if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity))) {
tmp = problem.gravity();
- tmp *= densityAtIP;
- pressureGrad -= tmp;
+ tmp *= 0.5*(vVars_i.density() + vVars_j.density());
+ potentialGrad_ -= tmp;
}
}
- void calculateVelocities_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemDat)
+ void calculateK_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
{
- Tensor K;
- problem.spatialParameters().meanK(K,
- problem.spatialParameters().intrinsicPermeability(element,
- fvElemGeom,
- face->i),
- problem.spatialParameters().intrinsicPermeability(element,
- fvElemGeom,
- face->j));
-
- // vector for the Darcy velocity
- K.mv(pressureGrad, vDarcy); // vDarcy = K * grad p
- vDarcyNormal = vDarcy*face->normal;
-
- // set the upstream and downstream vertices
- upstreamIdx = face->i;
- downstreamIdx = face->j;
- if (vDarcyNormal > 0)
- std::swap(upstreamIdx, downstreamIdx);
+ const SpatialParameters &sp = problem.spatialParameters();
+ sp.meanK(K_,
+ sp.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().i),
+ sp.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().j));
}
void calculateDiffCoeffPM_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
- const VolumeVariables &vDat_i = elemDat[face->i];
- const VolumeVariables &vDat_j = elemDat[face->j];
+ const VolumeVariables &vDat_i = elemDat[face().i];
+ const VolumeVariables &vDat_j = elemDat[face().j];
// Diffusion coefficient in the porous medium
- diffCoeffPM
- = 1./2*(vDat_i.porosity * vDat_i.tortuosity * vDat_i.diffCoeff +
- vDat_j.porosity * vDat_j.tortuosity * vDat_j.diffCoeff);
+ diffCoeffPM_
+ = 1./2*(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff() +
+ vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff());
}
void calculateDispersionTensor_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemDat)
+ const Element &element,
+ const ElementVolumeVariables &elemDat)
{
- const VolumeVariables &vDat_i = elemDat[face->i];
- const VolumeVariables &vDat_j = elemDat[face->j];
+ const VolumeVariables &vDat_i = elemDat[face().i];
+ const VolumeVariables &vDat_j = elemDat[face().j];
//calculate dispersivity at the interface: [0]: alphaL = longitudinal disp. [m], [1] alphaT = transverse disp. [m]
- Dune::FieldVector<Scalar, 2> dispersivity(0);
- dispersivity[0] = 0.5 * (vDat_i.dispersivity[0] + vDat_j.dispersivity[0]);
- dispersivity[1] = 0.5 * (vDat_i.dispersivity[1] + vDat_j.dispersivity[1]);
+ Scalar dispersivity[2];
+ dispersivity[0] = 0.5 * (vDat_i.dispersivity()[0] + vDat_j.dispersivity()[0]);
+ dispersivity[1] = 0.5 * (vDat_i.dispersivity()[1] + vDat_j.dispersivity()[1]);
//calculate velocity at interface: v = -1/mu * vDarcy = -1/mu * K * grad(p)
- GlobalPosition velocity(vDarcy);
- velocity *= -1/viscosityAtIP;
-
- dispersionTensor = 0;
+ Vector velocity;
+ Valgrind::CheckDefined(potentialGrad());
+ Valgrind::CheckDefined(K_);
+ K_.mv(potentialGrad(), velocity);
+ velocity /= - 0.5 * (vDat_i.viscosity() + vDat_j.viscosity());
//matrix multiplication of the velocity at the interface: vv^T
+ dispersionTensor_ = 0;
for (int i=0; i<dim; i++)
- {
for (int j = 0; j<dim; j++)
- {
- dispersionTensor[i][j]=velocity[i]*velocity[j];
- }
- }
+ dispersionTensor_[i][j]=velocity[i]*velocity[j];
+
//normalize velocity product --> vv^T/||v||, [m/s]
Scalar vNorm = velocity.two_norm();
- dispersionTensor /= vNorm;
+ dispersionTensor_ /= vNorm;
if (vNorm == 0)
- {
- dispersionTensor = 0;
- }
+ dispersionTensor_ = 0;
//multiply with dispersivity difference: vv^T/||v||*(alphaL - alphaT), [m^2/s] --> alphaL = longitudinal disp., alphaT = transverse disp.
- dispersionTensor *= (dispersivity[0] - dispersivity[1]);
+ dispersionTensor_ *= (dispersivity[0] - dispersivity[1]);
//add ||v||*alphaT to the main diagonal:vv^T/||v||*(alphaL - alphaT) + ||v||*alphaT, [m^2/s]
for (int i = 0; i<dim; i++)
- {
- dispersionTensor[i][i]+=vNorm*dispersivity[1];
- }
+ dispersionTensor_[i][i] += vNorm*dispersivity[1];
}
-public:
- const FVElementGeometry &fvElemGeom;
- const SCVFace *face;
- const SCV *insideSCV;
- const SCV *outsideSCV;
+ const FVElementGeometry &fvElemGeom_;
+ int scvfIdx_;
- //! pressure gradient
- GlobalPosition pressureGrad;
+ //! pressure potential gradient
+ Vector potentialGrad_;
//! concentratrion gradient
- GlobalPosition concentrationGrad;
-
- //! density of the fluid at the integration point
- Scalar densityAtIP;
-
- //! molar density of the fluid at the integration point
- Scalar molarDensityAtIP;
-
- //! viscosity of the fluid at the integration point
- Scalar viscosityAtIP;
+ Vector concentrationGrad_;
//! the effective diffusion coefficent in the porous medium
- Scalar diffCoeffPM;
-
+ Scalar diffCoeffPM_;
+
//! the dispersion tensor in the porous medium
- Tensor dispersionTensor;
-
- //! darcy velocity at the face
- GlobalPosition vDarcy;
-
- //! darcy velocity in direction of the face normal
- Scalar vDarcyNormal;
+ Tensor dispersionTensor_;
- //! local index of the upwind vertex for each phase
- int upstreamIdx;
- //! local index of the downwind vertex for each phase
- int downstreamIdx;
+ //! the intrinsic permeability tensor
+ Tensor K_;
};
} // end namepace
diff --git a/dumux/boxmodels/1p2c/1p2cindices.hh b/dumux/boxmodels/1p2c/1p2cindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..28e9bab32cbde45e5496c085c9132e76512ee090
--- /dev/null
+++ b/dumux/boxmodels/1p2c/1p2cindices.hh
@@ -0,0 +1,48 @@
+// $Id: 1p2cproperties.hh 3838 2010-07-15 08:31:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines the primary variable and equation indices used by
+ * the 1p2c model
+ */
+
+#ifndef DUMUX_1P2C_INDICES_HH
+#define DUMUX_1P2C_INDICES_HH
+
+namespace Dumux
+{
+
+/*!
+ * \brief The indices for the isothermal single-phase, two-component model.
+ */
+struct OnePTwoCIndices
+{
+ // Equation indices
+ static const int contiEqIdx = 0; //!< continuity equation index
+ static const int transEqIdx = 1; //!< transport equation index
+
+ // primary variable indices
+ static const int pressureIdx = 0; //!< pressure
+ static const int x1Idx = 1; //!< mole fraction of the second component
+};
+
+}
+
+#endif
+
diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
index cf234883f6c171301366343c8cbf642a005c16ba..0a4ab8bd239f0c015aeae141acf0d6f71bdf8649 100644
--- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh
+++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
@@ -21,9 +21,7 @@
#include <dumux/boxmodels/common/boxmodel.hh>
#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
-
#include <dumux/boxmodels/1p2c/1p2cvolumevariables.hh>
-
#include <dumux/boxmodels/1p2c/1p2cfluxvariables.hh>
#include <dune/common/collectivecommunication.hh>
@@ -44,45 +42,34 @@ protected:
typedef BoxLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
enum
{
- dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
- numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
- numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+ // indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ x1Idx = Indices::x1Idx,
- konti = Indices::konti,
- transport = Indices::transport,
+ // indices of the equations
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx,
};
+ static const Scalar upwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(UpwindAlpha));
+
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
-
- typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
- typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
-
- static const Scalar upwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(UpwindAlpha));
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
+ typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
public:
/*!
@@ -96,14 +83,19 @@ public:
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
- const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
- const VolumeVariables &vertDat = elemDat[scvIdx];
+ const VolumeVariables &volVars =
+ usePrevSol ?
+ this->prevVolVars_(scvIdx) :
+ this->curVolVars_(scvIdx);
// storage term of continuity equation
- result[konti] = 0;
+ result[contiEqIdx] =
+ volVars.density()*volVars.porosity();
// storage term of the transport equation
- result[transport] = vertDat.molarDensity * vertDat.porosity * vertDat.molefraction;
+ result[transEqIdx] =
+ volVars.concentration(1) *
+ volVars.porosity();
}
/*!
@@ -112,44 +104,50 @@ public:
*/
void computeFlux(PrimaryVariables &flux, int faceId) const
{
- FluxVariables vars(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
- faceId,
- this->curVolVars_());
flux = 0;
-
- // data attached to upstream and the downstream vertices
- // of the current phase
- const VolumeVariables &up = this->curVolVars_(vars.upstreamIdx);
- const VolumeVariables &dn = this->curVolVars_(vars.downstreamIdx);
-
- flux[konti] = vars.vDarcyNormal / vars.viscosityAtIP;
-
- // advective flux
- flux[transport] +=
- vars.vDarcyNormal *
- ( upwindAlpha*
- ( up.molarDensity * up.molefraction/up.viscosity )
- +
- (1 - upwindAlpha)*
- ( dn.molarDensity * dn.molefraction/dn.viscosity ) );
-
- Dune::FieldVector<Scalar,dim> unitNormal(vars.face->normal);
- unitNormal/=vars.face->normal.two_norm();
-
- // diffusive flux
- flux[transport] +=
- vars.molarDensityAtIP * vars.diffCoeffPM *
- (vars.concentrationGrad * vars.face->normal);
-
- //multiply hydrodynamic dispersion tensor with the face normal
- Dune::FieldVector<Scalar,dim> normalDisp;
- vars.dispersionTensor.mv(vars.face->normal, normalDisp);
-
- //add dispersive flux
- flux[transport] +=
- vars.molarDensityAtIP * (normalDisp * vars.concentrationGrad);
+ FluxVariables fluxVars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceId,
+ this->curVolVars_());
+
+ Vector tmpVec;
+ fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(), tmpVec);
+
+ // "intrinsic" flux from cell i to cell j
+ Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
+
+ // total mass flux
+ flux[contiEqIdx] =
+ normalFlux *
+ (( upwindAlpha)*up.density()/up.viscosity()
+ +
+ ((1 - upwindAlpha)*dn.density()/dn.viscosity()));
+
+ // advective flux of the second component
+ flux[transEqIdx] +=
+ normalFlux *
+ (( upwindAlpha)*up.concentration(1)/up.viscosity()
+ +
+ (1 - upwindAlpha)*dn.concentration(1)/dn.viscosity());
+
+ // diffusive flux of second component
+ Scalar c = (up.concentration(1) + dn.concentration(1))/2;
+ flux[transEqIdx] +=
+ c * fluxVars.porousDiffCoeff() *
+ (fluxVars.concentrationGrad(1) * fluxVars.face().normal);
+
+ // dispersive flux of second component
+ Vector normalDisp;
+ fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
+ flux[transEqIdx] +=
+ c * (normalDisp * fluxVars.concentrationGrad(1));
+
+ // we need to calculate the flux from i to j, not the other
+ // way round...
+ flux *= -1;
}
/*!
diff --git a/dumux/boxmodels/1p2c/1p2cmodel.hh b/dumux/boxmodels/1p2c/1p2cmodel.hh
index 7641412ad6fb44f9eb82e8ac61639fed28e41af7..6c2c4cbebe8bc012892ec6a886728e83d2105e4b 100644
--- a/dumux/boxmodels/1p2c/1p2cmodel.hh
+++ b/dumux/boxmodels/1p2c/1p2cmodel.hh
@@ -18,9 +18,11 @@
#ifndef DUMUX_ONEP_TWOC_MODEL_HH
#define DUMUX_ONEP_TWOC_MODEL_HH
-#include "1p2clocalresidual.hh"
+#include "1p2cproperties.hh"
#include "1p2cproblem.hh"
+#include <dumux/boxmodels/common/boxmodel.hh>
+
namespace Dumux
{
@@ -105,7 +107,10 @@ public:
// create the required scalar fields
unsigned numVertices = this->problem_().gridView().size(dim);
ScalarField *pressure = writer.template createField<Scalar, 1>(numVertices);
- ScalarField *molefraction = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *moleFrac0 = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *moleFrac1 = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *massFrac0 = writer.template createField<Scalar, 1>(numVertices);
+ ScalarField *massFrac1 = writer.template createField<Scalar, 1>(numVertices);
unsigned numElements = this->gridView_().size(0);
ScalarField *rank =
@@ -137,17 +142,25 @@ public:
false);
- (*pressure)[globalIdx] = volVars.pressure;
- (*molefraction)[globalIdx] = volVars.molefraction;
+ (*pressure)[globalIdx] = volVars.pressure();
+ (*moleFrac0)[globalIdx] = volVars.moleFrac(0);
+ (*moleFrac1)[globalIdx] = volVars.moleFrac(1);
+ (*massFrac0)[globalIdx] = volVars.massFrac(0);
+ (*massFrac1)[globalIdx] = volVars.massFrac(1);
};
}
- writer.addVertexData(pressure, "pressure");
- writer.addVertexData(molefraction, "molefraction");
+ writer.addVertexData(pressure, "p");
+ writer.addVertexData(moleFrac0, "x_0");
+ writer.addVertexData(moleFrac1, "x_1");
+ writer.addVertexData(massFrac0, "X_0");
+ writer.addVertexData(massFrac1, "X_1");
writer.addCellData(rank, "process rank");
}
};
}
+#include "1p2cpropertydefaults.hh"
+
#endif
diff --git a/dumux/boxmodels/1p2c/1p2cproperties.hh b/dumux/boxmodels/1p2c/1p2cproperties.hh
index 3156cc8b0398827802fbdaee40a9607b3f090460..836cef83647a0dcf184670864c077d64d6b2f114 100644
--- a/dumux/boxmodels/1p2c/1p2cproperties.hh
+++ b/dumux/boxmodels/1p2c/1p2cproperties.hh
@@ -29,31 +29,6 @@
namespace Dumux
{
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-template<class TypeTag>
-class OnePTwoCBoxModel;
-
-template<class TypeTag>
-class OnePTwoCLocalResidual;
-
-template <class TypeTag>
-class OnePTwoCVolumeVariables;
-
-template <class TypeTag>
-class OnePTwoCFluxVariables;
-
-/*!
- * \brief The indices for the isothermal single-phase, two-component model.
- */
-struct OnePTwoCIndices
-{
- // Primary variable indices
- static const int konti = 0; //!< pressure in the solution vector
- static const int transport = 1; //!< mole fraction in the solution vector
-};
-
namespace Properties
{
@@ -75,37 +50,9 @@ NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(UpwindAlpha); //!< The default value of the upwind parameter
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-SET_INT_PROP(BoxOnePTwoC, NumEq, 2); //!< set the number of equations to 2
-SET_INT_PROP(BoxOnePTwoC, NumPhases, 1); //!< The number of phases in the 1p2c model is 1
-SET_INT_PROP(BoxOnePTwoC, NumComponents, 2); //!< The number of components in the 1p2c model is 2
-
-//! Use the 1p2c local jacobian operator for the 1p2c model
-SET_TYPE_PROP(BoxOnePTwoC,
- LocalResidual,
- OnePTwoCLocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(BoxOnePTwoC, Model, OnePTwoCBoxModel<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(BoxOnePTwoC, VolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
-
-
-
-
-//! the FluxVariables property
-SET_TYPE_PROP(BoxOnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
-
-//! the default upwind factor. Default 1.0, i.e. fully upwind...
-SET_SCALAR_PROP(BoxOnePTwoC, UpwindAlpha, 1.0);
-
-//! The indices required by the isothermal 1p2c model
-SET_TYPE_PROP(BoxOnePTwoC, OnePTwoCIndices, Dumux::OnePTwoCIndices);
+NEW_PROP_TAG(PhaseIndex); //!< The index of the phase in the fluid system
+NEW_PROP_TAG(Comp1Index); //!< The index of the first component in the fluid system
+NEW_PROP_TAG(Comp2Index); //!< The index of the second component in the fluid system
}
}
diff --git a/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..94c7a67903b44f94bf906b7432e13d9138c18030
--- /dev/null
+++ b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
@@ -0,0 +1,79 @@
+// $Id: 1p2cproperties.hh 3838 2010-07-15 08:31:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines some default values for the properties of the the
+ * single-phase, two-compenent BOX model.
+ */
+
+#ifndef DUMUX_1P2C_PROPERTY_DEFAULTS_HH
+#define DUMUX_1P2C_PROPERTY_DEFAULTS_HH
+
+#include "1p2cproperties.hh"
+
+#include "1p2cmodel.hh"
+#include "1p2clocalresidual.hh"
+#include "1p2cvolumevariables.hh"
+#include "1p2cfluxvariables.hh"
+#include "1p2cindices.hh"
+
+namespace Dumux
+{
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+SET_INT_PROP(BoxOnePTwoC, NumEq, 2); //!< set the number of equations to 2
+SET_INT_PROP(BoxOnePTwoC, NumPhases, 1); //!< The number of phases in the 1p2c model is 1
+SET_INT_PROP(BoxOnePTwoC, NumComponents, 2); //!< The number of components in the 1p2c model is 2
+
+//! Use the 1p2c local residual function for the 1p2c model
+SET_TYPE_PROP(BoxOnePTwoC, LocalResidual, OnePTwoCLocalResidual<TypeTag>);
+
+//! define the model
+SET_TYPE_PROP(BoxOnePTwoC, Model, OnePTwoCBoxModel<TypeTag>);
+
+//! define the VolumeVariables
+SET_TYPE_PROP(BoxOnePTwoC, VolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
+
+//! define the FluxVariables
+SET_TYPE_PROP(BoxOnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
+
+//! set default upwind factor to 1.0, i.e. fully upwind
+SET_SCALAR_PROP(BoxOnePTwoC, UpwindAlpha, 1.0);
+
+//! Set the indices used by the 1p2c model
+SET_TYPE_PROP(BoxOnePTwoC, OnePTwoCIndices, Dumux::OnePTwoCIndices);
+
+//! Set the default phase used by the fluid system to the first one
+SET_INT_PROP(BoxOnePTwoC, PhaseIndex, 0);
+
+//! Set the default for the first component the fluid system's first component
+SET_INT_PROP(BoxOnePTwoC, Comp1Index, 0);
+
+//! Set the default for the second component the fluid system's second component
+SET_INT_PROP(BoxOnePTwoC, Comp2Index, 1);
+}
+
+}
+
+#endif
+
diff --git a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
index 97c01cd0c0f5155a5896e9853b71df805e164170..7fc3ecd05bb3cf3ada19401caaff5bb4183c16bf 100644
--- a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
@@ -26,6 +26,8 @@
#include "1p2cfluidstate.hh"
+#include <dumux/boxmodels/common/boxvolumevariables.hh>
+
namespace Dumux
{
@@ -34,40 +36,37 @@ namespace Dumux
* finite volume in the single-phase, two-component model.
*/
template <class TypeTag>
-class OnePTwoCVolumeVariables
+class OnePTwoCVolumeVariables : public BoxVolumeVariables<TypeTag>
{
+ typedef BoxVolumeVariables<TypeTag> ParentType;
+
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
- typedef OnePTwoCFluidState<TypeTag> FluidState;
-
- typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
- dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
- konti = Indices::konti,
- transport = Indices::transport
+ phaseIndex = GET_PROP_VALUE(TypeTag, PTAG(PhaseIndex)),
+ comp1Index = GET_PROP_VALUE(TypeTag, PTAG(Comp1Index)),
+ comp2Index = GET_PROP_VALUE(TypeTag, PTAG(Comp2Index)),
+
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx
};
-
- typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
- typedef typename RefElemProp::Container ReferenceElements;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
-
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef OnePTwoCFluidState<TypeTag> FluidState;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar,dimWorld> Vector;
public:
/*!
@@ -80,48 +79,121 @@ public:
int scvIdx,
bool isOldSol)
{
- primaryVars_ = priVars;
+ ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
- porosity = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
- tortuosity = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx);
- dispersivity = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx);
+ porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
+ tortuosity_ = problem.spatialParameters().tortuosity(element, elemGeom, scvIdx);
+ dispersivity_ = problem.spatialParameters().dispersivity(element, elemGeom, scvIdx);
Scalar temperature = problem.temperature(element, elemGeom, scvIdx);
fluidState_.update(priVars, temperature);
- pressure = fluidState_.phasePressure(konti);
- molefraction = fluidState_.moleFrac(konti, transport);
- density = FluidSystem::phaseDensity(konti, temperature, pressure, fluidState_);
- molarDensity = FluidSystem::molarDensity(konti, temperature, pressure, fluidState_);
- viscosity = FluidSystem::phaseViscosity(konti, temperature, pressure, fluidState_);
- diffCoeff = FluidSystem::diffCoeff(konti, transport, konti, temperature, pressure, fluidState_);
+ viscosity_ = FluidSystem::phaseViscosity(phaseIndex,
+ temperature,
+ pressure(),
+ fluidState_);
+ diffCoeff_ = FluidSystem::diffCoeff(phaseIndex,
+ comp1Index,
+ comp2Index,
+ temperature,
+ pressure(),
+ *this);
+
+ Valgrind::CheckDefined(porosity_);
+ Valgrind::CheckDefined(viscosity_);
+ Valgrind::CheckDefined(tortuosity_);
+ Valgrind::CheckDefined(dispersivity_);
+ Valgrind::CheckDefined(diffCoeff_);
+ Valgrind::CheckDefined(fluidState_);
+ Valgrind::CheckDefined(*this);
}
+
+ /*!
+ * \brief Return the fluid configuration at the given primary
+ * variables
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
/*!
- * \brief Sets the evaluation point used in the by the local jacobian.
+ * \brief Returns density the of the fluid phase.
*/
- void setEvalPoint(const Implementation *ep)
- { }
+ Scalar density() const
+ { return fluidState_.density(phaseIndex); }
/*!
- * \brief Return the vector of primary variables
+ * \brief Returns mole fraction of a component in the phase
*/
- const PrimaryVariables &primaryVars() const
- { return primaryVars_; }
-
- Scalar porosity;
- Scalar density;
- Scalar viscosity;
- Scalar tortuosity;
- Dune::FieldVector<Scalar,dim> dispersivity;
- Scalar diffCoeff;
- Scalar molefraction;
- Scalar pressure;
- Scalar molarDensity;
- FluidState fluidState_;
+ Scalar moleFrac(int compIdx) const
+ { return fluidState_.moleFrac(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
+
+ /*!
+ * \brief Returns mass fraction of a component in the phase
+ */
+ Scalar massFrac(int compIdx) const
+ { return fluidState_.massFrac(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
+
+ /*!
+ * \brief Returns concentration of a component in the phase
+ */
+ Scalar concentration(int compIdx) const
+ { return fluidState_.concentration(phaseIndex, (compIdx==0)?comp1Index:comp2Index); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ */
+ Scalar pressure() const
+ { return fluidState_.phasePressure(phaseIndex); }
+
+ /*!
+ * \brief Returns the binary diffusion coefficient in the fluid
+ */
+ Scalar diffCoeff() const
+ { return diffCoeff_; }
+
+ /*!
+ * \brief Returns the tortuosity of the streamlines of the fluid.
+ */
+ Scalar tortuosity() const
+ { return tortuosity_; }
+
+ /*!
+ * \brief Returns the dispersivity of the fluid's streamlines.
+ */
+ const Vector &dispersivity() const
+ { return dispersivity_; }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperature of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(); }
+
+ /*!
+ * \brief Returns the dynamic viscosity [Pa s] of a given phase
+ * within the control volume.
+ */
+ Scalar viscosity() const
+ { return viscosity_; }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
protected:
- PrimaryVariables primaryVars_;
+ Scalar porosity_;
+ Scalar viscosity_;
+ Scalar tortuosity_;
+ Vector dispersivity_;
+ Scalar diffCoeff_;
+ FluidState fluidState_;
};
}
diff --git a/dumux/boxmodels/2p/2pfluxvariables.hh b/dumux/boxmodels/2p/2pfluxvariables.hh
index 94665217c4ee77ae1f659f8d5b67f9168e0da1e6..18eda48ac47320745347afe9435695701270c54b 100644
--- a/dumux/boxmodels/2p/2pfluxvariables.hh
+++ b/dumux/boxmodels/2p/2pfluxvariables.hh
@@ -26,6 +26,8 @@
#ifndef DUMUX_2P_FLUX_VARIABLES_HH
#define DUMUX_2P_FLUX_VARIABLES_HH
+#include "2pproperties.hh"
+
#include <dumux/common/math.hh>
namespace Dumux
@@ -47,7 +49,6 @@ class TwoPFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GridView::ctype CoordScalar;
typedef typename GridView::template Codim<0>::Entity Element;
@@ -59,16 +60,13 @@ class TwoPFluxVariables
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
};
- typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename FVElementGeometry::SubControlVolume SCV;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
public:
TwoPFluxVariables(const Problem &problem,
@@ -81,7 +79,6 @@ public:
scvfIdx_ = faceIdx;
for (int phase = 0; phase < numPhases; ++phase) {
- densityAtIP_[phase] = Scalar(0);
potentialGrad_[phase] = Scalar(0);
}
@@ -90,41 +87,35 @@ public:
};
public:
- /*!
- * \brief Return the pressure potential multiplied with the
- * intrinsic permeability which goes from vertex i to
- * vertex j.
- *
- * Note that the length of the face's normal is the area of the
- * phase, so this is not the actual velocity by the integral of
- * the velocity over the face's area. Also note that the phase
- * mobility is not yet included here since this would require a
- * decision on the upwinding approach (which is done in the
- * actual model).
+ /*
+ * \brief Return the intrinsic permeability.
*/
- Scalar KmvpNormal(int phaseIdx) const
- { return KmvpNormal_[phaseIdx]; }
+ const Tensor &intrinsicPermeability() const
+ { return K_; }
/*!
- * \brief Return the local index of the upstream control volume
- * for a given phase.
+ * \brief Return the pressure potential gradient.
*/
- int upstreamIdx(int phaseIdx) const
- { return upstreamIdx_[phaseIdx]; }
+ const Vector &potentialGrad(int phaseIdx) const
+ { return potentialGrad_[phaseIdx]; }
/*!
- * \brief Return the local index of the downstream control volume
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the downstream control volume
* for a given phase.
*/
- int downstreamIdx(int phaseIdx) const
- { return downstreamIdx_[phaseIdx]; }
+ int downstreamIdx(Scalar normalFlux) const
+ { return (normalFlux >= 0)?face().j:face().i; }
/*!
- * \brief Return density [kg/m^3] of a phase at the integration
- * point.
+ * \brief Given the intrinisc permeability times the pressure
+ * potential gradient and SCV face normal for a phase,
+ * return the local index of the upstream control volume
+ * for a given phase.
*/
- Scalar densityAtIP(int phaseIdx) const
- { return densityAtIP_[phaseIdx]; }
+ int upstreamIdx(Scalar normalFlux) const
+ { return (normalFlux > 0)?face().i:face().j; }
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
@@ -135,23 +126,9 @@ protected:
// gradients
Vector potentialGrad_[numPhases];
- Vector concentrationGrad_[numPhases];
-
- // density of each face at the integration point
- Scalar densityAtIP_[numPhases];
-
- // intrinsic permeability times pressure potential gradient
- // projected on the face normal
- Scalar KmvpNormal_[numPhases];
-
- // local index of the upwind vertex for each phase
- int upstreamIdx_[numPhases];
- // local index of the downwind vertex for each phase
- int downstreamIdx_[numPhases];
-
- // the diffusion coefficient for the porous medium
- Scalar porousDiffCoeff_[numPhases];
+ // intrinsic permeability
+ Tensor K_;
private:
void calculateGradients_(const Problem &problem,
@@ -159,7 +136,6 @@ private:
const ElementVolumeVariables &elemDat)
{
// calculate gradients
- Vector tmp(0.0);
for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
@@ -171,26 +147,34 @@ private:
for (int phase = 0; phase < numPhases; phase++)
{
// the pressure gradient
- tmp = feGrad;
+ Vector tmp(feGrad);
tmp *= elemDat[idx].pressure(phase);
potentialGrad_[phase] += tmp;
-
- // phase density
- densityAtIP_[phase]
- +=
- elemDat[idx].density(phase) *
- face().shapeValue[idx];
}
}
// correct the pressure gradients by the hydrostatic
// pressure due to gravity
- for (int phase=0; phase < numPhases; phase++)
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
- tmp = problem.gravity();
- tmp *= densityAtIP_[phase];
-
- potentialGrad_[phase] -= tmp;
+ // calculate the phase density at the integration
+ // point. for this we make sure only existing phases
+ // matter
+ Scalar SI = elemDat[face().i].saturation(phaseIdx);
+ Scalar SJ = elemDat[face().j].saturation(phaseIdx);
+ Scalar rhoI = elemDat[face().i].density(phaseIdx);
+ Scalar rhoJ = elemDat[face().j].density(phaseIdx);
+ Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
+ Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
+ if (fI + fJ == 0)
+ fI = fJ = 0.5; // doesn't matter because no phase is
+ // present in both cells!
+ Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
+
+ Vector tmp(problem.gravity());
+ tmp *= density;
+
+ potentialGrad_[phaseIdx] -= tmp;
}
}
@@ -199,34 +183,14 @@ private:
const ElementVolumeVariables &elemDat)
{
const SpatialParameters &spatialParams = problem.spatialParameters();
- // multiply the pressure potential with the intrinsic
- // permeability
- Tensor K;
- Vector Kmvp;
- for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
- {
- spatialParams.meanK(K,
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().i),
- spatialParams.intrinsicPermeability(element,
- fvElemGeom_,
- face().j));
- K.mv(potentialGrad_[phaseIdx], Kmvp);
- KmvpNormal_[phaseIdx] = - (Kmvp * face().normal);
- }
-
- // set the upstream and downstream vertices
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- upstreamIdx_[phaseIdx] = face().i;
- downstreamIdx_[phaseIdx] = face().j;
-
- if (KmvpNormal_[phaseIdx] < 0) {
- std::swap(upstreamIdx_[phaseIdx],
- downstreamIdx_[phaseIdx]);
- }
- }
+ // calculate the intrinsic permeability
+ spatialParams.meanK(K_,
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().i),
+ spatialParams.intrinsicPermeability(element,
+ fvElemGeom_,
+ face().j));
}
};
diff --git a/dumux/boxmodels/2p/2plocalresidual.hh b/dumux/boxmodels/2p/2plocalresidual.hh
index d0dc353cc64d04f2ca0c042ecea5bd6bc8b2c3e1..3080de8d92d0ff8862f014fb5c5d8ca02973b301 100644
--- a/dumux/boxmodels/2p/2plocalresidual.hh
+++ b/dumux/boxmodels/2p/2plocalresidual.hh
@@ -26,14 +26,8 @@
#include "2pproperties.hh"
-#include "2pvolumevariables.hh"
-#include "2pfluxvariables.hh"
-#include "2pfluidstate.hh"
-#include <dune/common/collectivecommunication.hh>
-#include <vector>
-#include <iostream>
namespace Dumux
{
@@ -53,7 +47,6 @@ protected:
typedef TwoPLocalResidual<TypeTag> ThisType;
typedef BoxLocalResidual<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
@@ -80,23 +73,12 @@ protected:
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(DofMapper)) DofMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
-
- typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
-
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
static const Scalar mobilityUpwindAlpha =
@@ -134,10 +116,10 @@ public:
void computeFlux(PrimaryVariables &flux, int faceIdx) const
{
FluxVariables vars(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
- faceIdx,
- this->curVolVars_());
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_());
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
asImp_()->computeDiffusiveFlux(flux, vars);
@@ -151,25 +133,34 @@ public:
* This method is called by compute flux and is mainly there for
* derived models to ease adding equations selectively.
*/
- void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
////////
// advective fluxes of all components in all phases
////////
+ Vector tmpVec;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
+ // calculate the flux in the normal direction of the
+ // current sub control volume face
+ fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx),
+ tmpVec);
+ Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
+
// data attached to upstream and the downstream vertices
// of the current phase
- const VolumeVariables &up = this->curVolVars_(vars.upstreamIdx(phaseIdx));
- const VolumeVariables &dn = this->curVolVars_(vars.downstreamIdx(phaseIdx));
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
// add advective flux of current component in current
// phase
int eqIdx = (phaseIdx == wPhaseIdx) ? contiWEqIdx : contiNEqIdx;
- flux[eqIdx] += vars.KmvpNormal(phaseIdx) * (mobilityUpwindAlpha * // upstream vertex
- (up.density(phaseIdx) * up.mobility(phaseIdx)) + (1
- - mobilityUpwindAlpha) * // downstream vertex
- (dn.density(phaseIdx) * dn.mobility(phaseIdx)));
+ flux[eqIdx] +=
+ normalFlux
+ *
+ (( mobilityUpwindAlpha)*up.density(phaseIdx)*up.mobility(phaseIdx)
+ +
+ (1 - mobilityUpwindAlpha)*dn.density(phaseIdx)*dn.mobility(phaseIdx));
}
}
diff --git a/dumux/boxmodels/2p/2pmodel.hh b/dumux/boxmodels/2p/2pmodel.hh
index d73bcdf4109ad49880ed260ae01b569607ec3546..1b203fc3885678ddf9c962219b5b2601a2f6f460 100644
--- a/dumux/boxmodels/2p/2pmodel.hh
+++ b/dumux/boxmodels/2p/2pmodel.hh
@@ -22,8 +22,8 @@
#ifndef DUMUX_TWOP_MODEL_HH
#define DUMUX_TWOP_MODEL_HH
+#include "2pproperties.hh"
#include "2plocalresidual.hh"
-#include "2pnewtoncontroller.hh"
#include "2pproblem.hh"
namespace Dumux
@@ -80,18 +80,13 @@ class TwoPModel : public BoxModel<TypeTag>
typedef TwoPModel<TypeTag> ThisType;
typedef BoxModel<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
@@ -102,6 +97,17 @@ class TwoPModel : public BoxModel<TypeTag>
};
public:
+ /*!
+ * \brief Returns the relative weight of a primary variable for
+ * calculating relative errors.
+ */
+ Scalar primaryVarWeight(int vertIdx, int pvIdx) const
+ {
+ if (Indices::pressureIdx == pvIdx)
+ return 1e-7;
+ return 1;
+ }
+
/*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
@@ -184,4 +190,6 @@ public:
};
}
+#include "2ppropertydefaults.hh"
+
#endif
diff --git a/dumux/boxmodels/2p/2pnewtoncontroller.hh b/dumux/boxmodels/2p/2pnewtoncontroller.hh
deleted file mode 100644
index 0e65a90bd131cf7f0b0307ea1a6da20dd1f9e53a..0000000000000000000000000000000000000000
--- a/dumux/boxmodels/2p/2pnewtoncontroller.hh
+++ /dev/null
@@ -1,135 +0,0 @@
-// $Id$
-/*****************************************************************************
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Institute of Hydraulic Engineering *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version, as long as this copyright notice *
- * is included in its original form. *
- * *
- * This program is distributed WITHOUT ANY WARRANTY. *
- *****************************************************************************/
-/*!
- * \file
- * \brief A newton controller for two-phase problems.
- *
- * This controller 'knows' what a 'physically meaningful' solution is
- * which allows the newton method to abort quicker if the solution is
- * way out of bounds.
- */
-#ifndef DUMUX_2P_NEWTON_CONTROLLER_HH
-#define DUMUX_2P_NEWTON_CONTROLLER_HH
-
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-namespace Dumux {
-/*!
- * \ingroup TwoPBoxModel
- * \brief A newton controller for two-phase problems.
- *
- * This controller 'knows' what a 'physically meaningful' solution is
- * which allows the newton method to abort quicker if the solution is
- * way out of bounds.
- */
-template <class TypeTag>
-class TwoPNewtonController : public NewtonController<TypeTag>
-{
- typedef TwoPNewtonController<TypeTag> ThisType;
- typedef NewtonController<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
-
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
- enum {
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
- pressureIdx = Indices::pressureIdx
- };
-
-public:
- TwoPNewtonController()
- {};
-
-#if 0
- /*!
- * \brief Update the error of the solution compared to the
- * previous iteration.
- */
- void newtonUpdateRelError(const SolutionVector &uOld,
- const SolutionVector &deltaU)
- {
- typedef typename SolutionVector::block_type FV;
-
- // get the maximum value of each primary in either the old or
- // the new solution
- Dune::FieldVector<Scalar, numEq> weight(1.0);
-
- // a change in pressure is considered to be less severe by a
- // factor of 10000 than a change saturation.
- weight[pressureIdx] = 1;// 1e-4;
-
- // calculate the relative error as the maximum relative
- // deflection in any degree of freedom.
- this->error_ = 0;
- int offenderVertIdx = -1;
- int offenderPVIdx = -1;
- Scalar offenderDelta = 0;
- for (int i = 0; i < int(uOld.size()); ++i) {
- for (int j = 0; j < FV::size; ++j) {
- // calculate the relative error at the current vertex
- // i and the current primary variable j
-// Scalar curErr = std::abs(deltaU[i][j] * weight[j]);
- Scalar curErr = std::abs(deltaU[i][j]/(1.0 + std::abs(uOld[i][j])));
-#if 0
- // make sure that the specified tolerance is not below
- // machine precision!
- typedef std::numeric_limits<Scalar> limits;
- const Scalar machinePrec =
- limits::epsilon()*100
- * std::max<Scalar>(Scalar(1e10)/limits::max(),
- std::abs(uOld[i][j]));
- if (this->tolerance_ < machinePrec*weight[j])
- {
- std::cerr << "Allowed tolerance ("
- << this->tolerance_
- << ") is below machine precision ("
- << machinePrec*weight[j]
- << ") at vertex " << i
- << ", primary var " << j << "!\n";
- if (curErr < machinePrec*weight[j])
- // if the machine precision is reached, we
- // accept the solution even if we're above the
- // tolerance!
- curErr = this->tolerance_/100;
- };
-#endif
-
- if (this->error_ < curErr) {
- offenderVertIdx = i;
- offenderPVIdx = j;
- offenderDelta = deltaU[i][j];
- this->error_ = curErr;
- };
- }
- };
-
- this->endIterMsg() << ", worst offender: vertex " << offenderVertIdx
- << ", primary var " << offenderPVIdx
- << ", delta: " << offenderDelta;
-
- this->model().gridView().comm().max(this->error_);
- }
-#endif
-};
-}
-
-#endif
diff --git a/dumux/boxmodels/2p/2pproblem.hh b/dumux/boxmodels/2p/2pproblem.hh
index 8e103203ea8de0927da99316b7182e8db4156345..3207088685d8a920bd7060c821b1cda1ce368129 100644
--- a/dumux/boxmodels/2p/2pproblem.hh
+++ b/dumux/boxmodels/2p/2pproblem.hh
@@ -20,8 +20,9 @@
#ifndef DUMUX_2P_PROBLEM_HH
#define DUMUX_2P_PROBLEM_HH
+#include "2pproperties.hh"
+
#include <dumux/boxmodels/common/boxproblem.hh>
-#include <dumux/material/fluidsystems/2p_system.hh>
namespace Dumux
{
@@ -43,7 +44,6 @@ class TwoPProblem : public BoxProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
// material properties
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
enum {
diff --git a/dumux/boxmodels/2p/2pproperties.hh b/dumux/boxmodels/2p/2pproperties.hh
index dd5475a6401183b50dccee3423de2d61aea1cff8..cd89088c00a40e68fe438952950aa89c1a3938c7 100644
--- a/dumux/boxmodels/2p/2pproperties.hh
+++ b/dumux/boxmodels/2p/2pproperties.hh
@@ -20,8 +20,10 @@
* \brief Defines the properties required for the twophase BOX model.
*/
-#ifndef DUMUX_2PPROPERTIES_HH
-#define DUMUX_2PPROPERTIES_HH
+#ifndef DUMUX_2P_PROPERTIES_HH
+#define DUMUX_2P_PROPERTIES_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
namespace Dumux
{
@@ -30,94 +32,6 @@ namespace Dumux
*/
// \{
-////////////////////////////////
-// forward declarations
-////////////////////////////////
-template<class TypeTag>
-class TwoPModel;
-
-template<class TypeTag>
-class TwoPLocalResidual;
-
-template<class TypeTag>
-class TwoPNewtonController;
-
-template <class TypeTag>
-class TwoPProblem;
-
-template <class TypeTag>
-class TwoPVolumeVariables;
-
-template <class TypeTag>
-class TwoPFluxVariables;
-
-template<class TypeTag>
-class FluidSystem2P;
-
-template<class TypeTag>
-class TwoPFluidState;
-
-/*!
- * \brief The common indices for the isothermal two-phase model.
- */
-struct TwoPCommonIndices
-{
- // Formulations
- static const int pwSn = 0; //!< Pw and Sn as primary variables
- static const int pnSw = 1; //!< Pn and Sw as primary variables
-
- // Phase indices
- static const int wPhaseIdx = 0; //!< Index of the wetting phase in a phase vector
- static const int nPhaseIdx = 1; //!< Index of the non-wetting phase in a phase vector
-};
-
-/*!
- * \brief The indices for the \f$p_w-S_n\f$ formulation of the
- * isothermal two-phase model.
- *
- * \tparam PVOffset The first index in a primary variable vector.
- */
-template <int formulation = TwoPCommonIndices::pwSn, int PVOffset = 0>
-struct TwoPIndices : public TwoPCommonIndices
-{
- // Primary variable indices
- static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
- static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
-
- // indices of the primary variables
- static const int pwIdx = PVOffset + 0; //!< Pressure index of the wetting phase
- static const int SnIdx = PVOffset + 1; //!< Saturation index of the wetting phase
-
- // indices of the equations
- static const int contiWEqIdx = PVOffset + 0; //!< Index of the continuity equation of the wetting phase
- static const int contiNEqIdx = PVOffset + 1; //!< Index of the continuity equation of the non-wetting phase
-};
-
-/*!
- * \brief The indices for the \f$p_w-S_n\f$ formulation of the
- * isothermal two-phase model.
- *
- * \tparam PVOffset The first index in a primary variable vector.
- */
-template <int PVOffset>
-struct TwoPIndices<TwoPCommonIndices::pnSw, PVOffset>
- : public TwoPCommonIndices
-{
- // Primary variable indices
- static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
- static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
-
- // indices of the primary variables
- static const int pnIdx = PVOffset + 0; //!< Pressure index of the wetting phase
- static const int SwIdx = PVOffset + 1; //!< Saturation index of the wetting phase
-
- // indices of the equations
- static const int contiNEqIdx = PVOffset + 0; //!< Index of the continuity equation of the non-wetting phase
- static const int contiWEqIdx = PVOffset + 1; //!< Index of the continuity equation of the wetting phase
-};
-
-// \}
-
////////////////////////////////
// properties
////////////////////////////////
@@ -150,83 +64,8 @@ NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extra
NEW_PROP_TAG(MaterialLawParams); //!< The context material law (extracted from the spatial parameters)
NEW_PROP_TAG(WettingPhase); //!< The wetting phase for two-phase models
NEW_PROP_TAG(NonwettingPhase); //!< The non-wetting phase for two-phase models
-NEW_PROP_TAG( FluidSystem ); //!<The fluid systems including the information about the phases
-NEW_PROP_TAG( FluidState ); //!<The phases state
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-SET_INT_PROP(BoxTwoP, NumEq, 2); //!< set the number of equations to 2
-SET_INT_PROP(BoxTwoP, NumPhases, 2); //!< The number of phases in the 2p model is 2
-
-//! Set the default formulation to pWsN
-SET_INT_PROP(BoxTwoP,
- Formulation,
- TwoPCommonIndices::pwSn);
-
-//! Use the 2p local jacobian operator for the 2p model
-SET_TYPE_PROP(BoxTwoP,
- LocalResidual,
- TwoPLocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(BoxTwoP, Model, TwoPModel<TypeTag>);
-
-//! the default newton controller for two-phase problems
-SET_TYPE_PROP(BoxTwoP, NewtonController, TwoPNewtonController<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(BoxTwoP, VolumeVariables, TwoPVolumeVariables<TypeTag>);
-
-//! the FluxVariables property
-SET_TYPE_PROP(BoxTwoP, FluxVariables, TwoPFluxVariables<TypeTag>);
-
-//! the upwind factor for the mobility.
-SET_SCALAR_PROP(BoxTwoP, MobilityUpwindAlpha, 1.0);
-
-//! The indices required by the isothermal 2p model
-SET_PROP(BoxTwoP, TwoPIndices)
-{
- typedef TwoPIndices<GET_PROP_VALUE(TypeTag, PTAG(Formulation)), 0> type;
-};
-
-/*!
- * \brief Set the property for the material law by retrieving it from
- * the spatial parameters.
- */
-SET_PROP(BoxTwoP, MaterialLaw)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
-
-public:
- typedef typename SpatialParameters::MaterialLaw type;
-};
-
-/*!
- * \brief Set the property for the material parameters by extracting
- * it from the material law.
- */
-SET_PROP(BoxTwoP, MaterialLawParams)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
-
-public:
- typedef typename MaterialLaw::Params type;
-};
-
-SET_TYPE_PROP(BoxTwoP, FluidSystem, FluidSystem2P<TypeTag>);
-
-SET_TYPE_PROP(BoxTwoP, FluidState, TwoPFluidState<TypeTag>);
-
-// enable jacobian matrix recycling by default
-SET_BOOL_PROP(BoxTwoP, EnableJacobianRecycling, true);
-// enable partial reassembling by default
-SET_BOOL_PROP(BoxTwoP, EnablePartialReassemble, true);
-// enable time-step ramp up by default
-SET_BOOL_PROP(BoxTwoP, EnableTimeStepRampUp, true);
+NEW_PROP_TAG(FluidSystem); //!<The fluid systems including the information about the phases
+NEW_PROP_TAG(FluidState); //!<The phases state
// \}
}
diff --git a/dumux/boxmodels/2p/2ppropertydefaults.hh b/dumux/boxmodels/2p/2ppropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e769025d6be88511de5403bdb362d8c43e307bd4
--- /dev/null
+++ b/dumux/boxmodels/2p/2ppropertydefaults.hh
@@ -0,0 +1,180 @@
+// $Id: 2pproperties.hh 3357 2010-03-25 13:02:05Z lauser $
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the twophase BOX model.
+ */
+#ifndef DUMUX_2P_PROPERTY_DEFAULTS_HH
+#define DUMUX_2P_PROPERTY_DEFAULTS_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+#include "2pmodel.hh"
+#include "2pproblem.hh"
+#include "2pfluxvariables.hh"
+#include "2pvolumevariables.hh"
+#include "2pfluidstate.hh"
+#include <dumux/material/fluidsystems/2p_system.hh>
+
+namespace Dumux
+{
+/*!
+ * \addtogroup TwoPBoxModel
+ */
+// \{
+
+/*!
+ * \brief The common indices for the isothermal two-phase model.
+ */
+struct TwoPCommonIndices
+{
+ // Formulations
+ static const int pwSn = 0; //!< Pw and Sn as primary variables
+ static const int pnSw = 1; //!< Pn and Sw as primary variables
+
+ // Phase indices
+ static const int wPhaseIdx = 0; //!< Index of the wetting phase in a phase vector
+ static const int nPhaseIdx = 1; //!< Index of the non-wetting phase in a phase vector
+};
+
+/*!
+ * \brief The indices for the \f$p_w-S_n\f$ formulation of the
+ * isothermal two-phase model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <int formulation = TwoPCommonIndices::pwSn, int PVOffset = 0>
+struct TwoPIndices : public TwoPCommonIndices
+{
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
+
+ // indices of the primary variables
+ static const int pwIdx = PVOffset + 0; //!< Pressure index of the wetting phase
+ static const int SnIdx = PVOffset + 1; //!< Saturation index of the wetting phase
+
+ // indices of the equations
+ static const int contiWEqIdx = PVOffset + 0; //!< Index of the continuity equation of the wetting phase
+ static const int contiNEqIdx = PVOffset + 1; //!< Index of the continuity equation of the non-wetting phase
+};
+
+/*!
+ * \brief The indices for the \f$p_w-S_n\f$ formulation of the
+ * isothermal two-phase model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <int PVOffset>
+struct TwoPIndices<TwoPCommonIndices::pnSw, PVOffset>
+ : public TwoPCommonIndices
+{
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
+
+ // indices of the primary variables
+ static const int pnIdx = PVOffset + 0; //!< Pressure index of the wetting phase
+ static const int SwIdx = PVOffset + 1; //!< Saturation index of the wetting phase
+
+ // indices of the equations
+ static const int contiNEqIdx = PVOffset + 0; //!< Index of the continuity equation of the non-wetting phase
+ static const int contiWEqIdx = PVOffset + 1; //!< Index of the continuity equation of the wetting phase
+};
+
+// \}
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property defaults
+//////////////////////////////////////////////////////////////////
+SET_INT_PROP(BoxTwoP, NumEq, 2); //!< set the number of equations to 2
+SET_INT_PROP(BoxTwoP, NumPhases, 2); //!< The number of phases in the 2p model is 2
+
+//! Set the default formulation to pWsN
+SET_INT_PROP(BoxTwoP,
+ Formulation,
+ TwoPCommonIndices::pwSn);
+
+//! Use the 2p local jacobian operator for the 2p model
+SET_TYPE_PROP(BoxTwoP,
+ LocalResidual,
+ TwoPLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoP, Model, TwoPModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoP, VolumeVariables, TwoPVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoP, FluxVariables, TwoPFluxVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(BoxTwoP, MobilityUpwindAlpha, 1.0);
+
+//! The indices required by the isothermal 2p model
+SET_PROP(BoxTwoP, TwoPIndices)
+{
+ typedef TwoPIndices<GET_PROP_VALUE(TypeTag, PTAG(Formulation)), 0> type;
+};
+
+/*!
+ * \brief Set the property for the material law by retrieving it from
+ * the spatial parameters.
+ */
+SET_PROP(BoxTwoP, MaterialLaw)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+
+public:
+ typedef typename SpatialParameters::MaterialLaw type;
+};
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_PROP(BoxTwoP, MaterialLawParams)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+
+public:
+ typedef typename MaterialLaw::Params type;
+};
+
+SET_TYPE_PROP(BoxTwoP, FluidSystem, FluidSystem2P<TypeTag>);
+
+SET_TYPE_PROP(BoxTwoP, FluidState, TwoPFluidState<TypeTag>);
+
+// enable jacobian matrix recycling by default
+SET_BOOL_PROP(BoxTwoP, EnableJacobianRecycling, true);
+// enable partial reassembling by default
+SET_BOOL_PROP(BoxTwoP, EnablePartialReassemble, true);
+// enable time-step ramp up by default
+SET_BOOL_PROP(BoxTwoP, EnableTimeStepRampUp, false);
+
+// \}
+}
+
+}
+
+#endif
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index bc03a0b386c201ac1aa93dbfa42734cd0892e0e6..0d1dcf1607e9cd4676046041091cf37f301f618f 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -23,6 +23,9 @@
#define DUMUX_2P_VOLUME_VARIABLES_HH
#include "2pproperties.hh"
+#include "2pfluidstate.hh"
+
+#include <dune/common/fvector.hh>
namespace Dumux
{
@@ -142,15 +145,6 @@ public:
scvIdx);
}
- void updateTemperature_(const PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int scvIdx,
- const Problem &problem)
- {
- temperature_ = problem.temperature(element, elemGeom, scvIdx);
- }
-
/*!
* \brief Return the vector of primary variables
*/
@@ -220,6 +214,15 @@ public:
{ return porosity_; }
protected:
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ temperature_ = problem.temperature(element, elemGeom, scvIdx);
+ }
+
PrimaryVariables primaryVars_;
FluidState fluidState_;
Scalar porosity_;
diff --git a/dumux/boxmodels/2p2c/2p2cfluidstate.hh b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
index 6e342d0ba6def6cd6163399ff09714bf9776a493..51ee03f639530b865433e92d2f9f53ba692ac064 100644
--- a/dumux/boxmodels/2p2c/2p2cfluidstate.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
@@ -299,8 +299,8 @@ public:
/*!
* \brief Return the fugacity of a component [Pa].
*/
- Scalar fugacity(int phaseIdx, int compIdx) const
- { return moleFrac(phaseIdx, compIdx)*phasePressure(compIdx); };
+ Scalar fugacity(int compIdx) const
+ { return moleFrac(gPhaseIdx, compIdx)*phasePressure(gPhaseIdx); };
/*!
* \brief Returns the capillary pressure [Pa]
diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
index 08815ae866f4cdecae1234f4131c81a4e4ccc565..ecf24523b4cf365dad13b35f0b1f01f346f5df5d 100644
--- a/dumux/boxmodels/common/boxelementboundarytypes.hh
+++ b/dumux/boxmodels/common/boxelementboundarytypes.hh
@@ -18,8 +18,9 @@
#include "boxproperties.hh"
-#include <vector>
-#include <dumux/common/boundarytypes.hh>
+#include <dune/grid/common/geometry.hh>
+
+#include <dumux/common/valgrind.hh>
namespace Dumux
{
@@ -36,7 +37,6 @@ class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTa
typedef std::vector<BoundaryTypes> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
diff --git a/dumux/boxmodels/common/boxelementvolumevariables.hh b/dumux/boxmodels/common/boxelementvolumevariables.hh
index 012028c3ee86229258d8a7d5b6e1d1bc6c7bd54c..5e30507c236e2b6328e343bd8d548073d2926e84 100644
--- a/dumux/boxmodels/common/boxelementvolumevariables.hh
+++ b/dumux/boxmodels/common/boxelementvolumevariables.hh
@@ -18,8 +18,6 @@
#include "boxproperties.hh"
-#include <vector>
-#include <dumux/common/boundarytypes.hh>
namespace Dumux
{
@@ -36,7 +34,6 @@ class BoxElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(Type
typedef std::vector<VolumeVariables> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
diff --git a/dumux/boxmodels/common/boxfvelementgeometry.hh b/dumux/boxmodels/common/boxfvelementgeometry.hh
index 8ac1925afe60518843ffd45d62626d6062ef301c..3225ce0aba7c5d47364ea093432efe74a38a8040 100644
--- a/dumux/boxmodels/common/boxfvelementgeometry.hh
+++ b/dumux/boxmodels/common/boxfvelementgeometry.hh
@@ -17,7 +17,6 @@
#define DUMUX_BOX_FV_ELEMENTGEOMETRY_HH
#include <dune/grid/common/intersectioniterator.hh>
-#include <dune/grid/common/capabilities.hh>
#include <dumux/common/propertysystem.hh>
namespace Dumux
diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index a0e9775d4800a400d45f0e1be5fa14546d53786e..fc77fb39ec9d03a9817bb4797b4a2bbb0b89de83 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -26,24 +26,9 @@
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH
-#include <dune/common/exceptions.hh>
-
-#include <dumux/common/valgrind.hh>
-
-#include <dune/grid/common/genericreferenceelements.hh>
-
-#include <boost/format.hpp>
-
-#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
-#include "boxelementvolumevariables.hh"
-#include "boxfvelementgeometry.hh"
-#include "boxlocalresidual.hh"
-
-#include "boxproperties.hh"
-
-
+#include "boxelementboundarytypes.hh"
namespace Dumux
{
@@ -96,7 +81,6 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index e2d26a40213a918de536bbbed278216227a0b288..13f530a21bb6c0c54a86e54506c475768893fab4 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -26,19 +26,12 @@
#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
#define DUMUX_BOX_LOCAL_RESIDUAL_HH
-#include <dune/common/exceptions.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
-#include <dune/grid/common/genericreferenceelements.hh>
-
-#include <boost/format.hpp>
-
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
#include "boxproperties.hh"
-#include "boxfvelementgeometry.hh"
-#include "boxelementboundarytypes.hh"
namespace Dumux
{
@@ -69,8 +62,8 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GridView::Grid::ctype CoordScalar;
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -84,7 +77,6 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index 4044bcbf9a44a93aa142cf789f3e803b0dafefcb..0099d56858cd947e3a2ba78eefc9a05165624839 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -17,19 +17,13 @@
#ifndef DUMUX_BOX_MODEL_HH
#define DUMUX_BOX_MODEL_HH
-#include <dumux/common/valgrind.hh>
-#include <dune/grid/common/genericreferenceelements.hh>
-
-#include <boost/format.hpp>
-
#include "boxproperties.hh"
+#include "boxpropertydefaults.hh"
#include "boxelementvolumevariables.hh"
#include "boxlocaljacobian.hh"
#include "boxlocalresidual.hh"
-#include "pdelabboxlocaloperator.hh"
-
namespace Dumux
{
@@ -68,7 +62,6 @@ class BoxModel
typedef typename GET_PROP_TYPE(TypeTag, PTAG(DofMapper)) DofMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
enum {
@@ -85,7 +78,6 @@ class BoxModel
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
@@ -96,6 +88,8 @@ class BoxModel
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };
+
// copying a model is not a good idea
BoxModel(const BoxModel &);
@@ -204,13 +198,30 @@ public:
* \brief Reference to the current solution as a block vector.
*/
const SolutionVector &curSol() const
- { return uCur_; }
+ {
+ return uCur_;
+ if (enablePartialReassemble &&
+ jacobianAssembler().inJacobianAssemble())
+ {
+ return jacobianAssembler().evalPoint();
+ }
+ return uCur_;
+ }
/*!
* \brief Reference to the current solution as a block vector.
*/
SolutionVector &curSol()
- { return uCur_; }
+ {
+ return uCur_;
+
+ if (enablePartialReassemble &&
+ jacobianAssembler().inJacobianAssemble())
+ {
+ return jacobianAssembler().evalPoint();
+ }
+ return uCur_;
+ }
/*!
* \brief Reference to the previous solution as a block vector.
@@ -264,6 +275,37 @@ public:
const LocalResidual &localResidual() const
{ return localJacobian().localResidual(); }
+ /*!
+ * \brief Returns the relative weight of a primary variable for
+ * calculating relative errors.
+ */
+ Scalar primaryVarWeight(int vertIdx, int pvIdx) const
+ { return 1.0; }
+
+ /*!
+ * \brief Returns the relative error between two vectors of
+ * primary variables.
+ *
+ * \todo The vertexIdx argument is pretty hacky. it is required by
+ * models with pseudo primary variables (i.e. the primary
+ * variable switching models). the clean solution would be
+ * to access the pseudo primary variables from the primary
+ * variables.
+ */
+ Scalar relativeErrorVertex(int vertexIdx,
+ const PrimaryVariables &pv1,
+ const PrimaryVariables &pv2)
+ {
+ Scalar result = 0.0;
+ for (int j = 0; j < numEq; ++j) {
+ Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
+ Scalar eqErr = std::abs(pv1[j] - pv2[j])*weight;
+
+ result = std::max(result, eqErr);
+ }
+ return result;
+ }
+
/*!
* \brief Try to progress the model to the next timestep.
*/
diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index 79c28fe11b1e5922e830d516717ce26b1e1e977d..fdd496d1d441305e2a1ba556026e2ad12449abd3 100644
--- a/dumux/boxmodels/common/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -25,8 +25,6 @@
#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/io/restart.hh>
-#include <dumux/common/timemanager.hh>
-
namespace Dumux
{
/*!
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index 5517f9cc5ceaa20f3ec6bdb0ea1ad473a630c7c9..f064db2d8ad118bf663837d91f1f9ee0b79fc9f7 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -108,280 +108,6 @@ NEW_PROP_TAG(DofMapper);
}
}
-
-#include <dumux/nonlinear/newtonmethod.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-#include "boxfvelementgeometry.hh"
-#include "pdelabboxassembler.hh"
-#include "pdelabboxistlvectorbackend.hh"
-//#include <dumux/boxmodels/pdelab/boxdirichletconstraints.hh>
-
-#include <dumux/common/boundarytypes.hh>
-#include <dumux/common/timemanager.hh>
-
-namespace Dumux {
-
-template<typename TypeTag>
-class BoxFVElementGeometry;
-
-template<typename TypeTag>
-class BoxElementBoundaryTypes;
-
-template<typename TypeTag>
-class BoxElementVolumeVariables;
-
-template<typename TypeTag>
-class BoxLocalJacobian;
-
-namespace PDELab {
-template<typename TypeTag>
-class BoxLocalOperator;
-
-template<typename TypeTag>
-class BoxAssembler;
-
-template<typename TypeTag>
-class BoxISTLVectorBackend;
-};
-
-namespace Properties {
-//////////////////////////////////////////////////////////////////
-// Some defaults for very fundamental properties
-//////////////////////////////////////////////////////////////////
-
-//! Set the default type for scalar values to double
-SET_PROP_DEFAULT(Scalar)
-{ typedef double type; };
-
-//! Set the default type for the time manager
-SET_PROP_DEFAULT(TimeManager)
-{ typedef Dumux::TimeManager<TypeTag> type; };
-
-/*!
- * \brief Specify the reference elements which we ought to use.
- *
- * We use Dune::ReferenceElements by default (-> old entity
- * numbering).
- *
- * TODO: Some specialization if the grid only supports one kind of
- * cells would be nice. this would be better fixed inside DUNE,
- * though. something like:
- * Dune::GenericReferenceElements<Dune::GeometryType<cube, dim> >
- */
-SET_PROP_DEFAULT(ReferenceElements)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-
- typedef typename Grid::ctype CoordScalar;
- static const int dim = Grid::dimension;
-
-public:
- typedef Dune::GenericReferenceElements<CoordScalar, dim> Container;
- typedef Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-};
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-//! Use the leaf grid view if not defined otherwise
-SET_PROP(BoxModel, GridView)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-
-public:
- typedef typename Grid::LeafGridView type;
-};
-
-//! Set the default for the FVElementGeometry
-SET_PROP(BoxModel, FVElementGeometry)
-{
- typedef Dumux::BoxFVElementGeometry<TypeTag> type;
-};
-
-//! Set the default for the ElementBoundaryTypes
-SET_PROP(BoxModel, ElementBoundaryTypes)
-{ typedef Dumux::BoxElementBoundaryTypes<TypeTag> type; };
-
-//! use the plain newton method for the box scheme by default
-SET_PROP(BoxModel, NewtonMethod)
-{public:
- typedef Dumux::NewtonMethod<TypeTag> type;
-};
-
-//! use the plain newton controller for the box scheme by default
-SET_PROP(BoxModel, NewtonController)
-{public:
- typedef Dumux::NewtonController<TypeTag> type;
-};
-
-//! Mapper for the grid view's vertices.
-SET_PROP(BoxModel, VertexMapper)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- template<int dim>
- struct VertexLayout {
- bool contains (Dune::GeometryType gt) const
- { return gt.dim() == 0; }
- };
-public:
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, VertexLayout> type;
-};
-
-//! Mapper for the grid view's elements.
-SET_PROP(BoxModel, ElementMapper)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- template<int dim>
- struct ElementLayout {
- bool contains (Dune::GeometryType gt) const
- { return gt.dim() == dim; }
- };
-public:
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, ElementLayout> type;
-};
-
-//! Mapper for the degrees of freedoms.
-SET_PROP(BoxModel, DofMapper)
-{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) type; };
-
-//! The local jacobian operator for the box scheme
-SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);
-
-/*!
- * \brief The type of a solution for the whole grid at a fixed time.
- */
-SET_PROP(BoxModel, SolutionVector)
-{ private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
-public:
- typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type type;
-};
-
-/*!
- * \brief The type of a solution for a whole element.
- */
-SET_PROP(BoxModel, ElementSolutionVector)
-{ private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
-public:
- typedef Dune::BlockVector<PrimaryVariables> type;
-};
-
-/*!
- * \brief A vector of primary variables.
- */
-SET_PROP(BoxModel, PrimaryVariables)
-{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector))::block_type type; };
-
-/*!
- * \brief An array of secondary variable containers.
- */
-SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
-
-/*!
- * \brief Boundary types at a single degree of freedom.
- */
-SET_PROP(BoxModel, BoundaryTypes)
-{ private:
- enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
-public:
- typedef Dumux::BoundaryTypes<numEq> type;
-};
-
-/*!
- * \brief Assembler for the global jacobian matrix.
- */
-SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::PDELab::BoxAssembler<TypeTag>);
-
-//! Extract the type of a global jacobian matrix from the solution types
-SET_PROP(BoxModel, JacobianMatrix)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
-public:
- typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type type;
-};
-
-//! use the local FEM space associated with cubes by default
-SET_PROP(BoxModel, LocalFEMSpace)
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum{dim = GridView::dimension};
-
-public:
- typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type;
-};
-
-SET_PROP(BoxModel, GridFunctionSpace)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-public:
- //typedef MyBoxConstraints Constraints;
- typedef Dune::PDELab::NoConstraints Constraints;
-
- typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints, Dumux::PDELab::BoxISTLVectorBackend<TypeTag> >
- ScalarGridFunctionSpace;
-
- typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, Dune::PDELab::GridFunctionSpaceBlockwiseMapper>
- type;
-
- typedef typename type::template ConstraintsContainer<Scalar>::Type
- ConstraintsTrafo;
-};
-
-SET_PROP(BoxModel, ConstraintsTrafo)
-{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ConstraintsTrafo type; };
-SET_PROP(BoxModel, Constraints)
-{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::Constraints type; };
-SET_PROP(BoxModel, ScalarGridFunctionSpace)
-{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ScalarGridFunctionSpace type; };
-
-SET_PROP(BoxModel, GridOperatorSpace)
-{ private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-
-public:
- typedef Dumux::PDELab::BoxLocalOperator<TypeTag> LocalOperator;
-
- typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
- GridFunctionSpace,
- LocalOperator,
- ConstraintsTrafo,
- ConstraintsTrafo,
- Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
- true
- > type;
-};
-
-
-SET_PROP(BoxModel, LocalOperator)
-{ typedef typename GET_PROP(TypeTag, PTAG(GridOperatorSpace))::LocalOperator type; };
-
-// disable jacobian matrix recycling by default
-SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
-// disable partial reassembling by default
-SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);
-// disable time-step ramp up by default
-SET_BOOL_PROP(BoxModel, EnableTimeStepRampUp, false);
-
// \}
-}
-}
-
#endif
diff --git a/dumux/boxmodels/common/boxpropertydefaults.hh b/dumux/boxmodels/common/boxpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2ecbbfafac6333cef1d11a76378e8681d6a96fd1
--- /dev/null
+++ b/dumux/boxmodels/common/boxpropertydefaults.hh
@@ -0,0 +1,277 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2010 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_BOX_PROPERTY_DEFAULTS_HH
+#define DUMUX_BOX_PROPERTY_DEFAULTS_HH
+
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include "pdelabboxassembler.hh"
+#include "pdelabboxistlvectorbackend.hh"
+#include "boxfvelementgeometry.hh"
+#include "boxelementboundarytypes.hh"
+#include "boxlocaljacobian.hh"
+#include "boxelementvolumevariables.hh"
+
+#include <dune/pdelab/finiteelementmap/p1fem.hh>
+#include <dune/pdelab/finiteelementmap/q1fem.hh>
+#include <dune/pdelab/backend/istlmatrixbackend.hh>
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/timemanager.hh>
+
+#include "boxproperties.hh"
+
+namespace Dumux {
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Some defaults for very fundamental properties
+//////////////////////////////////////////////////////////////////
+
+//! Set the default type for scalar values to double
+SET_PROP_DEFAULT(Scalar)
+{ typedef double type; };
+
+//! Set the default type for the time manager
+SET_PROP_DEFAULT(TimeManager)
+{ typedef Dumux::TimeManager<TypeTag> type; };
+
+/*!
+ * \brief Specify the reference elements which we ought to use.
+ *
+ * We use Dune::ReferenceElements by default (-> old entity
+ * numbering).
+ *
+ * TODO: Some specialization if the grid only supports one kind of
+ * cells would be nice. this would be better fixed inside DUNE,
+ * though. something like:
+ * Dune::GenericReferenceElements<Dune::GeometryType<cube, dim> >
+ */
+SET_PROP_DEFAULT(ReferenceElements)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+
+ typedef typename Grid::ctype CoordScalar;
+ static const int dim = Grid::dimension;
+
+public:
+ typedef Dune::GenericReferenceElements<CoordScalar, dim> Container;
+ typedef Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
+};
+
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+//! Use the leaf grid view if not defined otherwise
+SET_PROP(BoxModel, GridView)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+
+public:
+ typedef typename Grid::LeafGridView type;
+};
+
+//! Set the default for the FVElementGeometry
+SET_PROP(BoxModel, FVElementGeometry)
+{
+ typedef Dumux::BoxFVElementGeometry<TypeTag> type;
+};
+
+//! Set the default for the ElementBoundaryTypes
+SET_PROP(BoxModel, ElementBoundaryTypes)
+{ typedef Dumux::BoxElementBoundaryTypes<TypeTag> type; };
+
+//! use the plain newton method for the box scheme by default
+SET_PROP(BoxModel, NewtonMethod)
+{public:
+ typedef Dumux::NewtonMethod<TypeTag> type;
+};
+
+//! use the plain newton controller for the box scheme by default
+SET_PROP(BoxModel, NewtonController)
+{public:
+ typedef Dumux::NewtonController<TypeTag> type;
+};
+
+//! Mapper for the grid view's vertices.
+SET_PROP(BoxModel, VertexMapper)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ template<int dim>
+ struct VertexLayout {
+ bool contains (Dune::GeometryType gt) const
+ { return gt.dim() == 0; }
+ };
+public:
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, VertexLayout> type;
+};
+
+//! Mapper for the grid view's elements.
+SET_PROP(BoxModel, ElementMapper)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ template<int dim>
+ struct ElementLayout {
+ bool contains (Dune::GeometryType gt) const
+ { return gt.dim() == dim; }
+ };
+public:
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, ElementLayout> type;
+};
+
+//! Mapper for the degrees of freedoms.
+SET_PROP(BoxModel, DofMapper)
+{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) type; };
+
+//! The local jacobian operator for the box scheme
+SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);
+
+/*!
+ * \brief The type of a solution for the whole grid at a fixed time.
+ */
+SET_PROP(BoxModel, SolutionVector)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+public:
+ typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type type;
+};
+
+/*!
+ * \brief The type of a solution for a whole element.
+ */
+SET_PROP(BoxModel, ElementSolutionVector)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+public:
+ typedef Dune::BlockVector<PrimaryVariables> type;
+};
+
+/*!
+ * \brief A vector of primary variables.
+ */
+SET_PROP(BoxModel, PrimaryVariables)
+{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector))::block_type type; };
+
+/*!
+ * \brief An array of secondary variable containers.
+ */
+SET_TYPE_PROP(BoxModel, ElementVolumeVariables, Dumux::BoxElementVolumeVariables<TypeTag>);
+
+/*!
+ * \brief Boundary types at a single degree of freedom.
+ */
+SET_PROP(BoxModel, BoundaryTypes)
+{ private:
+ enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
+public:
+ typedef Dumux::BoundaryTypes<numEq> type;
+};
+
+/*!
+ * \brief Assembler for the global jacobian matrix.
+ */
+SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::PDELab::BoxAssembler<TypeTag>);
+
+//! Extract the type of a global jacobian matrix from the solution types
+SET_PROP(BoxModel, JacobianMatrix)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
+public:
+ typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type type;
+};
+
+//! use the local FEM space associated with cubes by default
+SET_PROP(BoxModel, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type;
+};
+
+SET_PROP(BoxModel, GridFunctionSpace)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
+public:
+ //typedef MyBoxConstraints Constraints;
+ typedef Dune::PDELab::NoConstraints Constraints;
+
+ typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints, Dumux::PDELab::BoxISTLVectorBackend<TypeTag> >
+ ScalarGridFunctionSpace;
+
+ typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, Dune::PDELab::GridFunctionSpaceBlockwiseMapper>
+ type;
+
+ typedef typename type::template ConstraintsContainer<Scalar>::Type
+ ConstraintsTrafo;
+};
+
+SET_PROP(BoxModel, ConstraintsTrafo)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ConstraintsTrafo type; };
+SET_PROP(BoxModel, Constraints)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::Constraints type; };
+SET_PROP(BoxModel, ScalarGridFunctionSpace)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ScalarGridFunctionSpace type; };
+
+SET_PROP(BoxModel, GridOperatorSpace)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
+
+public:
+ typedef Dumux::PDELab::BoxLocalOperator<TypeTag> LocalOperator;
+
+ typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
+ GridFunctionSpace,
+ LocalOperator,
+ ConstraintsTrafo,
+ ConstraintsTrafo,
+ Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
+ true
+ > type;
+};
+
+
+SET_PROP(BoxModel, LocalOperator)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridOperatorSpace))::LocalOperator type; };
+
+// disable jacobian matrix recycling by default
+SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
+// disable partial reassembling by default
+SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);
+// disable time-step ramp up by default
+SET_BOOL_PROP(BoxModel, EnableTimeStepRampUp, false);
+
+} // namespace Properties
+} // namespace Dumux
+
+#endif
diff --git a/dumux/boxmodels/common/boxvolumevariables.hh b/dumux/boxmodels/common/boxvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f858ad92f4ba7050194b4eee83921836fdb85852
--- /dev/null
+++ b/dumux/boxmodels/common/boxvolumevariables.hh
@@ -0,0 +1,116 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Base class for the model specific classes which provide
+ * access to all volume averaged quantities.
+ */
+#ifndef DUMUX_BOX_VOLUME_VARIABLES_HH
+#define DUMUX_BOX_VOLUME_VARIABLES_HH
+
+#include "boxproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \brief Base class for the model specific classes which provide
+ * access to all volume averaged quantities.
+ */
+template <class TypeTag>
+class BoxVolumeVariables
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+
+public:
+ // default constructor
+ BoxVolumeVariables()
+ { evalPoint_ = 0; };
+
+ // copy constructor
+ BoxVolumeVariables(const BoxVolumeVariables &v)
+ {
+ primaryVars_ = v.primaryVars_;
+ evalPoint_ = 0;
+ };
+
+ // assignment operator
+ BoxVolumeVariables &operator=(const BoxVolumeVariables &v)
+ {
+ evalPoint_ = 0;
+ return *this;
+ };
+
+ /*!
+ * \brief Sets the evaluation point used in the by the local jacobian.
+ */
+ void setEvalPoint(const Implementation *ep)
+ { evalPoint_ = ep; }
+
+ /*!
+ * \brief Returns the evaluation point used in the by the local jacobian.
+ */
+ const Implementation &evalPoint() const
+ { return (evalPoint_ == 0)?asImp_():evalPoint_; }
+
+ /*!
+ * \brief Update all quantities for a given control volume.
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ bool isOldSol)
+ {
+ primaryVars_ = priVars;
+ }
+
+ /*!
+ * \brief Return the vector of primary variables
+ */
+ const PrimaryVariables &primaryVars() const
+ { return primaryVars_; }
+
+ /*!
+ * \brief Return a component of primary variable vector
+ */
+ Scalar primaryVars(int pvIdx) const
+ { return primaryVars_[pvIdx]; }
+
+protected:
+ const Implementation &asImp_() const
+ { return *static_cast<Implementation*>(this); }
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ PrimaryVariables primaryVars_;
+
+ // the evaluation point of the local jacobian
+ const Implementation *evalPoint_;
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index aee42fe7cff5e3b9f909610efa2e368c4df2c12e..f3f87d80eb9e77473e9a1695c90b98b9794f02fb 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -16,12 +16,6 @@
#ifndef DUMUX_PDELAB_BOX_ASSEMBLER_HH
#define DUMUX_PDELAB_BOX_ASSEMBLER_HH
-#include<dune/pdelab/finiteelementmap/p1fem.hh>
-#include<dune/pdelab/finiteelementmap/q1fem.hh>
-#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
-#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
-#include<dune/pdelab/backend/istlvectorbackend.hh>
-#include<dune/pdelab/backend/istlmatrixbackend.hh>
//#include "pdelabboundarytypes.hh"
#include "pdelabboxlocaloperator.hh"
@@ -35,9 +29,7 @@ class BoxAssembler
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum{dim = GridView::dimension};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
@@ -50,10 +42,11 @@ class BoxAssembler
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalOperator)) LocalOperator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ enum{dim = GridView::dimension};
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -66,7 +59,9 @@ class BoxAssembler
enum {
enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)),
- enableJacobianRecycling = GET_PROP_VALUE(TypeTag, PTAG(EnableJacobianRecycling))
+ enableJacobianRecycling = GET_PROP_VALUE(TypeTag, PTAG(EnableJacobianRecycling)),
+
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))
};
// copying the jacobian assembler is not a good idea
@@ -106,6 +101,7 @@ public:
void init(Problem& problem)
{
+ inJacobianAssemble_ = false;
problemPtr_ = &problem;
fem_ = new FEM();
//cn_ = new Constraints(*problemPtr_);
@@ -145,27 +141,44 @@ public:
int numVerts = gridView_().size(dim);
int numElems = gridView_().size(0);
residual_.resize(numVerts);
-
+
// initialize data needed for partial reassembly
if (enablePartialReassemble) {
+ evalPoint_.resize(numVerts);
vertexColor_.resize(numVerts);
elementColor_.resize(numElems);
}
- std::fill(vertexColor_.begin(),
- vertexColor_.end(),
- Red);
- std::fill(elementColor_.begin(),
- elementColor_.end(),
- Red);
+ reassembleAll();
}
+ SolutionVector &evalPoint()
+ { return evalPoint_; }
+
+ const SolutionVector &evalPoint() const
+ { return evalPoint_; }
+
+ bool inJacobianAssemble() const
+ { return inJacobianAssemble_; }
+
void assemble(SolutionVector &u)
{
// assemble the global jacobian matrix
if (!reuseMatrix_) {
+ if (enablePartialReassemble) {
+ // move the evaluation points of red vertices
+ for (int i = 0; i < vertexColor_.size(); ++i) {
+ if (vertexColor_[i] == Red)
+ evalPoint_[i] = model_().curSol()[i];
+ }
+ }
+
+ reassembleTolerance_ = nextReassembleTolerance_;
+
// we actually need to reassemle!
resetMatrix_();
+ inJacobianAssemble_ = true;
gridOperatorSpace_->jacobian(u, *matrix_);
+ inJacobianAssemble_ = false;
}
reuseMatrix_ = false;
@@ -195,31 +208,58 @@ public:
reuseMatrix_ = yesno;
}
+
void reassembleAll()
{
- std::fill(vertexColor_.begin(),
- vertexColor_.end(),
- Red);
- std::fill(elementColor_.begin(),
- elementColor_.end(),
- Red);
+ nextReassembleTolerance_ = 0.0;
+ if (enablePartialReassemble) {
+ std::fill(vertexColor_.begin(),
+ vertexColor_.end(),
+ Red);
+ std::fill(elementColor_.begin(),
+ elementColor_.end(),
+ Red);
+ }
}
-
+
+ Scalar reassembleTolerance() const
+ { return reassembleTolerance_; }
+
void markVertexRed(int globalVertIdx, bool yesno)
{
if (enablePartialReassemble)
vertexColor_[globalVertIdx] = yesno?Red:Green;
}
- void computeColors()
+ void computeColors(Scalar relTol)
{
-
if (!enablePartialReassemble)
return;
ElementIterator elemIt = gridView_().template begin<0>();
ElementIterator elemEndIt = gridView_().template end<0>();
+ nextReassembleTolerance_ = 0;
+
+ int redVert = 0;
+ // mark the red vertices
+ for (int i = 0; i < vertexColor_.size(); ++i) {
+ const PrimaryVariables &evalPv = evalPoint_[i];
+ const PrimaryVariables &solPv = model_().curSol()[i];
+
+ Scalar vertErr = model_().relativeErrorVertex(i,
+ evalPv,
+ solPv);
+
+ // mark vertex as red or green
+ vertexColor_[i] = (vertErr > relTol)?Red:Green;
+ if (vertErr > relTol)
+ ++redVert;
+ else
+ nextReassembleTolerance_ =
+ std::max(nextReassembleTolerance_, vertErr);
+ };
+
// Mark all red elements
for (; elemIt != elemEndIt; ++elemIt) {
bool needReassemble = false;
@@ -281,12 +321,16 @@ public:
}
int numTot = gridView_().size(0);
+ numTot = gridView_().comm().sum(numTot);
+ numGreen = gridView_().comm().sum(numGreen);
+ elemsReasm = numTot - numGreen;
problem_().newtonController().endIterMsg()
<< ", reassemble "
<< numTot - numGreen << "/" << numTot
<< " (" << 100*Scalar(numTot - numGreen)/numTot << "%) elems";
};
+ int elemsReasm;
int vertexColor(const Element &element, int vertIdx) const
{
if (!enablePartialReassemble)
@@ -388,10 +432,14 @@ private:
Matrix *matrix_;
bool reuseMatrix_;
+ bool inJacobianAssemble_;
+ SolutionVector evalPoint_;
std::vector<EntityColor> vertexColor_;
std::vector<EntityColor> elementColor_;
std::vector<int> ghostIndices_;
+ Scalar nextReassembleTolerance_;
+ Scalar reassembleTolerance_;
SolutionVector residual_;
};
diff --git a/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh b/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
index 20fd2e9375f15d2e112169cf3a45840b6c2f6f69..b0031deb85b1399dfa359506e4805fb811ad7491 100644
--- a/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
+++ b/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
@@ -7,10 +7,10 @@
#ifndef DUMUX_BOX_SOLUTION_VECTOR_HH
#define DUMUX_BOX_SOLUTION_VECTOR_HH
-#include<vector>
+#include "boxproperties.hh"
-#include<dune/common/fvector.hh>
-#include<dune/istl/bvector.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
namespace Dumux {
namespace Properties {
diff --git a/dumux/boxmodels/common/pdelabboxlocaloperator.hh b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
index e60f24c99fb367877d63cabd122b3c1699cc3264..3ae4b10e9e56cae1956c0547aefe60d481ac2adf 100644
--- a/dumux/boxmodels/common/pdelabboxlocaloperator.hh
+++ b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
@@ -16,15 +16,11 @@
#ifndef DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
#define DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
-#include<vector>
-#include<dune/common/fvector.hh>
-#include<dune/common/geometrytype.hh>
-#include<dune/grid/common/quadraturerules.hh>
-#include<dune/pdelab/gridoperatorspace/gridoperatorspace.hh>
-#include<dune/pdelab/gridoperatorspace/gridoperatorspaceutilities.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
+#include "boxproperties.hh"
+
namespace Dumux {
namespace PDELab {
@@ -41,9 +37,6 @@ class BoxLocalOperator
BoxLocalOperator(const BoxLocalOperator &);
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
public:
diff --git a/dumux/common/boundarytypes.hh b/dumux/common/boundarytypes.hh
index af4994c0f61b61645f774f7726be26d84a9d0e76..d4f4762518862c8eb474e039047c14338b89c32b 100644
--- a/dumux/common/boundarytypes.hh
+++ b/dumux/common/boundarytypes.hh
@@ -20,10 +20,7 @@
#ifndef BOUNDARY_TYPES_HH
#define BOUNDARY_TYPES_HH
-#include <boost/format.hpp>
-#include <dune/common/timer.hh>
-#include <dune/common/mpihelper.hh>
#include <dumux/common/valgrind.hh>
diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index 86d71e4c6becc12e289d03e7d1991ee84c71e28b..27949ece11477ab8bb8c39a4cae1bf571e795126 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -20,11 +20,12 @@
#ifndef DUMUX_MATH_HH
#define DUMUX_MATH_HH
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
#include <cmath>
namespace Dumux
{
-
/*!
* \brief Calculate the harmonic mean of two scalar values.
*/
@@ -109,7 +110,7 @@ bool isLarger(const Dune::FieldVector<Scalar, dim> &pos,
*/
template <class Scalar, int dim>
bool isSmaller(const Dune::FieldVector<Scalar, dim> &pos,
- const Dune::FieldVector<Scalar, dim> &largerVec)
+ const Dune::FieldVector<Scalar, dim> &largerVec)
{
for (int i=0; i < dim; i++)
{
diff --git a/dumux/common/pdelabpreconditioner.hh b/dumux/common/pdelabpreconditioner.hh
index 211d7384c505823c9e2d474e6b95cbf8687b947a..ebc69260bca21c514945cbada2d2e2dd9451358c 100644
--- a/dumux/common/pdelabpreconditioner.hh
+++ b/dumux/common/pdelabpreconditioner.hh
@@ -16,25 +16,40 @@
#ifndef DUMUX_PDELAB_PRECONDITIONER_HH
#define DUMUX_PDELAB_PRECONDITIONER_HH
-#include<dune/pdelab/backend/istlsolverbackend.hh>
+#include <dumux/common/pardiso.hh>
+#include <dumux/common/propertysystem.hh>
+#include <dune/pdelab/backend/istlsolverbackend.hh>
namespace Dumux {
-namespace PDELab {
+// forward declaration of property tags
+namespace Properties {
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(Model);
+NEW_PROP_TAG(GridView);
+NEW_PROP_TAG(NumEq);
+NEW_PROP_TAG(Grid);
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(VertexMapper);
+NEW_PROP_TAG(GridOperatorSpace);
+NEW_PROP_TAG(JacobianMatrix);
+NEW_PROP_TAG(ReferenceElements);
+NEW_PROP_TAG(GridFunctionSpace);
+NEW_PROP_TAG(ConstraintsTrafo);
+};
+
+namespace PDELab {
template<class TypeTag>
class Exchanger
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
dim = GridView::dimension
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
typedef typename JacobianMatrix::block_type BlockType;
@@ -296,7 +311,6 @@ template<class TypeTag>
class ISTLBackend_NoOverlap_BCGS_ILU
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef Dune::PDELab::ParallelISTLHelper<GridFunctionSpace> PHELPER;
@@ -384,7 +398,6 @@ template<class TypeTag>
class ISTLBackend_NoOverlap_Loop_Pardiso
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
diff --git a/dumux/common/propertysystem.hh b/dumux/common/propertysystem.hh
index c3a5e7de4e7912aff2953a960fd63543169ee2f6..71c973fb2cfbea53bfd5bc6ee07c59f837b3661b 100644
--- a/dumux/common/propertysystem.hh
+++ b/dumux/common/propertysystem.hh
@@ -39,12 +39,10 @@
#include <boost/type_traits.hpp>
// Integral Constant Expressions
-#include <boost/type_traits/ice.hpp>
// string formating
#include <boost/format.hpp>
-#include <iostream>
namespace Dumux
{
diff --git a/dumux/common/start.hh b/dumux/common/start.hh
index 9f5df7b5a935bef94eb73b485e3bf21a405f3a75..dff01bc06ffed83695da65f767a4b439db33cacd 100644
--- a/dumux/common/start.hh
+++ b/dumux/common/start.hh
@@ -22,19 +22,23 @@
#include <dumux/common/propertysystem.hh>
+#include <dune/grid/io/file/dgfparser.hh>
#include <dune/common/mpihelper.hh>
-
-#include <dune/grid/common/gridinfo.hh>
-
-#include <boost/format.hpp>
-
#include <iostream>
namespace Dumux
{
+// forward declaration of property tags
+namespace Properties
+{
+NEW_PROP_TAG(Grid);
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(TimeManager);
+}
+
void printUsage(const char *progname)
{
- std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd dt\n")%progname;
+ std::cout << "usage: " << progname << " [--restart restartTime] gridFile.dgf tEnd dt\n";
exit(1);
};
@@ -53,7 +57,6 @@ int startFromDGF(int argc, char **argv)
try {
#endif
- typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Grid)) Grid;
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(TimeManager)) TimeManager;
@@ -91,7 +94,6 @@ int startFromDGF(int argc, char **argv)
GridPointer gridPtr(dgfFileName);
if (mpiHelper.size() > 1)
gridPtr->loadBalance();
- Dune::gridinfo(*gridPtr);
// instantiate and run the concrete problem
TimeManager timeManager;
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index ac118299e1186a75ab5cf488340e532d24480616..60304f11a36a78f899f803a08df34c771b7e2c02 100644
--- a/dumux/common/timemanager.hh
+++ b/dumux/common/timemanager.hh
@@ -27,7 +27,6 @@
#include <dune/common/timer.hh>
#include <dune/common/mpihelper.hh>
-#include <dumux/io/restart.hh>
namespace Dumux
{
@@ -162,7 +161,9 @@ public:
* step size won't exceed the episode, though.
*/
void setTimeStepSize(Scalar stepSize)
- { timeStepSize_ = stepSize; }
+ {
+ timeStepSize_ = std::min(stepSize, maxTimeStepSize());
+ }
/*!
* \brief Returns a suggested timestep length so that we don't
@@ -197,6 +198,19 @@ public:
bool willBeFinished() const
{ return finished_ || time() + timeStepSize() >= endTime(); }
+ /*!
+ * \brief Aligns dt to the episode boundary or the end time of the
+ * simulation.
+ */
+ Scalar maxTimeStepSize() const
+ {
+ if (finished())
+ return 0.0;
+
+ return
+ std::min(episodeMaxTimeStepSize(),
+ std::max(0.0, endTime() - time()));
+ };
/*
* @}
@@ -335,9 +349,10 @@ public:
// notify the problem that the timestep is done and ask it
// for a suggestion for the next timestep size
- Scalar nextDt =
- std::min(problem_->nextTimeStepSize(),
- episodeMaxTimeStepSize());
+ // set the time step size for the next step
+ setTimeStepSize(problem_->nextTimeStepSize());
+
+ Scalar nextDt = timeStepSize();
if (verbose_) {
std::cout <<
@@ -345,8 +360,6 @@ public:
%timeStepIndex()%timer.elapsed()%time()%dt%nextDt;
}
- // set the time step size for the next step
- setTimeStepSize(nextDt);
}
if (verbose_)
diff --git a/dumux/material/components/brine.hh b/dumux/material/components/brine.hh
index cad93e679711973e58a570f02fd3125fe2eb3965..8624c071f37527b4d9e90bc81ff3c507e7a0759f 100644
--- a/dumux/material/components/brine.hh
+++ b/dumux/material/components/brine.hh
@@ -21,14 +21,10 @@
#ifndef DUMUX_BRINE_HH
#define DUMUX_BRINE_HH
-#include <dune/common/exceptions.hh>
-
-#include "h2o.hh"
#include "component.hh"
#include <cmath>
-#include <iostream>
namespace Dumux
{
diff --git a/dumux/material/components/ch4.hh b/dumux/material/components/ch4.hh
index f0cf690415feb1170e98621588559eeec49de89a..6cf7beca56907dd012a9bd5eb24cd0d7c07a0999 100644
--- a/dumux/material/components/ch4.hh
+++ b/dumux/material/components/ch4.hh
@@ -21,7 +21,6 @@
#define DUMUX_CH4_HH
#include <dumux/material/idealgas.hh>
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/components/h2.hh b/dumux/material/components/h2.hh
index 30f1d1ea497bc193273162398e7f10c8221ba959..103cf36ee2e9245cae39addc4d39370773a02a56 100644
--- a/dumux/material/components/h2.hh
+++ b/dumux/material/components/h2.hh
@@ -21,7 +21,6 @@
#define DUMUX_H2_HH
#include <dumux/material/idealgas.hh>
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh
index d7af489398b4c245d92ecd73b487c9b9c1d4d5e0..d2fee8c2d3b86234fd94fe04a5683835212d7c51 100644
--- a/dumux/material/components/h2o.hh
+++ b/dumux/material/components/h2o.hh
@@ -33,7 +33,6 @@
#include "iapws/region4.hh"
#include <cmath>
-#include <iostream>
namespace Dumux
diff --git a/dumux/material/components/n2.hh b/dumux/material/components/n2.hh
index 0640d6c78177f65c832258971708cf3a133ca113..46b3c33c1c15f6d7b7ed1b704405aadef4aae5be 100644
--- a/dumux/material/components/n2.hh
+++ b/dumux/material/components/n2.hh
@@ -22,7 +22,6 @@
#define DUMUX_N2_HH
#include <dumux/material/idealgas.hh>
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/components/o2.hh b/dumux/material/components/o2.hh
index 73be4d2a3e5389e448a83518a7a4cb3ab91cabf9..72700adc73b8ebcd4389844dfe6bfd720977a5a4 100644
--- a/dumux/material/components/o2.hh
+++ b/dumux/material/components/o2.hh
@@ -21,7 +21,6 @@
#define DUMUX_O2_HH
#include <dumux/material/idealgas.hh>
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/components/oil.hh b/dumux/material/components/oil.hh
index 5e176b5d9085c7ab46fd88c3b07e98be5188b947..135a44beabcea6b3eaa2438e8a5033a8d6902c2f 100644
--- a/dumux/material/components/oil.hh
+++ b/dumux/material/components/oil.hh
@@ -20,7 +20,6 @@
#ifndef DUMUX_OIL_HH
#define DUMUX_OIL_HH
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/components/simpleco2.hh b/dumux/material/components/simpleco2.hh
index a6e8675bf730de1356cefa04853a1df36c6a8a3b..67e2559871190da2378cb4ee6ea7eaf2095b416a 100644
--- a/dumux/material/components/simpleco2.hh
+++ b/dumux/material/components/simpleco2.hh
@@ -21,13 +21,11 @@
#ifndef DUMUX_SIMPLE_CO2_HH
#define DUMUX_SIMPLE_CO2_HH
-#include <dune/common/exceptions.hh>
#include <dumux/material/idealgas.hh>
#include "component.hh"
#include <cmath>
-#include <iostream>
namespace Dumux
{
diff --git a/dumux/material/components/simplednapl.hh b/dumux/material/components/simplednapl.hh
index 2ec18d3b36b8ea3ce9b14851f1e3d40f02813edf..ba94b1e3c90b3125a0bdb504a1b46765919c4d7f 100644
--- a/dumux/material/components/simplednapl.hh
+++ b/dumux/material/components/simplednapl.hh
@@ -22,12 +22,9 @@
#ifndef DUMUX_SIMPLE_DNAPL_HH
#define DUMUX_SIMPLE_DNAPL_HH
-#include <dune/common/exceptions.hh>
#include "component.hh"
-#include <cmath>
-#include <iostream>
namespace Dumux
{
diff --git a/dumux/material/components/simpleh2o.hh b/dumux/material/components/simpleh2o.hh
index 20df9f76aca1fdde50e6aac89867597ca96afad8..f8f09cf52387790195a2dafe37a20609133d5767 100644
--- a/dumux/material/components/simpleh2o.hh
+++ b/dumux/material/components/simpleh2o.hh
@@ -23,12 +23,10 @@
#define DUMUX_SIMPLE_H2O_HH
#include <dumux/material/idealgas.hh>
-#include <dune/common/exceptions.hh>
#include "component.hh"
#include <cmath>
-#include <iostream>
namespace Dumux
{
diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh
index e754f6cd5019664760ac3ba8c9d6c96e2aeda92d..1d133b5d772c3f6453b9ac514eb2d6d1822a29b9 100644
--- a/dumux/material/components/tabulatedcomponent.hh
+++ b/dumux/material/components/tabulatedcomponent.hh
@@ -25,8 +25,12 @@
#ifndef DUMUX_TABULATED_COMPONENT_HH
#define DUMUX_TABULATED_COMPONENT_HH
+#include <dumux/common/exceptions.hh>
+
#include <boost/math/special_functions/fpclassify.hpp>
+#include <iostream>
+
namespace Dumux
{
diff --git a/dumux/material/components/unit.hh b/dumux/material/components/unit.hh
index 333b290fd6eca77b0b8dc47df2c523b39f318b68..5b071eaad395170dcee7eb8e035c82c347e60915 100644
--- a/dumux/material/components/unit.hh
+++ b/dumux/material/components/unit.hh
@@ -21,7 +21,6 @@
#ifndef DUMUX_UNIT_HH
#define DUMUX_UNIT_HH
-#include <dune/common/exceptions.hh>
#include "component.hh"
diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index 2bf6fa2f0d702ef3ff63447d7ea23f47ba0d9464..7b6eb70abed48a53f038af1e37a6d263cb38ee2e 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -26,8 +26,6 @@
#include <algorithm>
-#include <math.h>
-#include <assert.h>
namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
index 7e96ef2b33c582903c468b80a224e902338dd276..6a8bf280ebdb231e375fc3af13bca704e07922d0 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
@@ -22,6 +22,8 @@
#include "linearmaterialparams.hh"
+#include <algorithm>
+
namespace Dumux
{
/*!
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
index e5b4b664b59216fecfe9d4319ce769e0ea648529..c38a172ffe00019723320b24312d4ec040464828 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
@@ -25,10 +25,7 @@
#include "brookscorey.hh"
#include "regularizedbrookscoreyparams.hh"
-#include <algorithm>
-#include <math.h>
-#include <assert.h>
#include <dumux/common/spline.hh>
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
index c65b2869e355c77754af013a2b6cb6e21efb5ddc..65b520a6dcc6ab39286df8d0388d12e84a7cb792 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
@@ -52,11 +52,22 @@ public:
* \brief Threshold saturation below which the capillary pressure
* is regularized.
*
- * This is just 5%. If you need a different value, overload this
+ * This is just 1%. If you need a different value, overload this
* class.
*/
Scalar thresholdSw() const
- { return 0.05; }
+ {
+ // Some problems are very sensitive to this value
+ // (e.g. makeing it smaller might result in negative
+ // pressures), if you change it here, you will almost
+ // certainly break someone's code!
+ //
+ // If you want to use a different regularization threshold,
+ // overload this class and supply the new class as second
+ // template parameter for the RegularizedVanGenuchten law!
+ return /* PLEASE DO _NOT_ */ 1e-2; /* CHANGE THIS VALUE. READ
+ * COMMENT ABOVE! */
+ }
};
}; // namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
index be33b5fb71194fbb2f769fd5a9d3a9e6024acf48..926579900f556bd2af101e3f3708519f8a33bc05 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
@@ -20,12 +20,8 @@
#define DUMUX_REGULARIZED_LINEAR_MATERIAL_HH
#include "linearmaterial.hh"
-#include "linearmaterialparams.hh"
-#include <algorithm>
-#include <math.h>
-#include <assert.h>
#include <dumux/common/spline.hh>
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index 80233098b23258e086f2e2ce3518f899a1ff92a9..91eccad9aeeb60c8182de4cfa7c8915830d421f1 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -27,8 +27,6 @@
#include <algorithm>
-#include <math.h>
-#include <assert.h>
#include <dumux/common/spline.hh>
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
index 40955fbe6ba8056b451af8ebd19eb8f98f1a0c18..8a3c0091c701b833b4bf0f4cbfbbbfc8ad920fd0 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
@@ -21,7 +21,6 @@
#ifndef REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
#define REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
-#include "vangenuchten.hh"
#include "vangenuchtenparams.hh"
namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
index 88603f05c3c8a9b328e9ef6fa1b7860e14c12e81..d44fb9350f6282cfa921639bb3b8c8f686d214f5 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
@@ -24,8 +24,6 @@
#include <algorithm>
-#include <math.h>
-#include <assert.h>
namespace Dumux
{
diff --git a/dumux/material/fluidstate.hh b/dumux/material/fluidstate.hh
index 3e0dc01f919b29f2ef98a8a6071f7c686a158f18..771de886ec359edce094050de79d8a306b03f51f 100644
--- a/dumux/material/fluidstate.hh
+++ b/dumux/material/fluidstate.hh
@@ -23,6 +23,8 @@
#ifndef DUMUX_FLUID_STATE_HH
#define DUMUX_FLUID_STATE_HH
+#include <dumux/common/exceptions.hh>
+
namespace Dumux
{
/*!
diff --git a/dumux/material/fluidsystems/2p_system.hh b/dumux/material/fluidsystems/2p_system.hh
index eea985f1ed9fdd5443f4199305cf8eada7409742..b9670cd2d9e851ae77aac6222a8189cf0cb7070d 100644
--- a/dumux/material/fluidsystems/2p_system.hh
+++ b/dumux/material/fluidsystems/2p_system.hh
@@ -25,8 +25,18 @@
#include "liquidphase.hh"
#include "gasphase.hh"
-namespace Dumux
-{
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/propertysystem.hh>
+
+namespace Dumux {
+// forward defintions of the property tags
+namespace Properties {
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(WettingPhase);
+NEW_PROP_TAG(NonwettingPhase);
+};
+
/*!
* \brief A compositional fluid with water and molecular nitrogen as
diff --git a/dumux/material/fluidsystems/h2o_n2_system.hh b/dumux/material/fluidsystems/h2o_n2_system.hh
index 83a80a3460380f47e3b420b43b2a26df83de2f9d..a9ec010b5e952a71d903a10ed110aa6043d61834 100644
--- a/dumux/material/fluidsystems/h2o_n2_system.hh
+++ b/dumux/material/fluidsystems/h2o_n2_system.hh
@@ -33,6 +33,12 @@
namespace Dumux
{
+// forward defintions of the property tags
+namespace Properties {
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(Components);
+};
+
/*!
* \brief A compositional fluid with water and molecular nitrogen as
* components in both, the liquid and the gas phase.
diff --git a/dumux/material/fluidsystems/isfluid_trail_system.hh b/dumux/material/fluidsystems/isfluid_trail_system.hh
index 83fe66b74cc68e59c1f5e71482ed52416d88caa4..d9ffb4acde05e7c1476b9961a6795e87b73fca11 100644
--- a/dumux/material/fluidsystems/isfluid_trail_system.hh
+++ b/dumux/material/fluidsystems/isfluid_trail_system.hh
@@ -20,6 +20,8 @@
#ifndef DUMUX_ISFLUID_TRAIL_SYSTEM_HH
#define DUMUX_ISFLUID_TRAIL_SYSTEM_HH
+#include <dune/common/exceptions.hh>
+
#include <dumux/common/propertysystem.hh>
#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
@@ -29,6 +31,7 @@ namespace Dumux
namespace Properties
{
NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(OnePTwoCIndices);
};
/*!
@@ -41,12 +44,13 @@ class ISFluid_Trail_System
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+public:
enum {
- isfluid = Indices::konti,
- trail = Indices::transport
+ // component indices
+ isFluidIdx = 0,
+ trailIdx = 1
};
-public:
static void init()
{}
@@ -57,9 +61,9 @@ public:
{
switch(compIdx)
{
- case isfluid:
+ case isFluidIdx:
return "ISFluid";
- case trail:
+ case trailIdx:
return "Trail";
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
@@ -67,41 +71,39 @@ public:
}
/*!
- * \brief Given all mole fractions in a phase, return the phase
- * density [kg/m^3].
+ * \brief Return the molar mass of a component [kg/mol].
*/
- template <class FluidState>
- static Scalar phaseDensity(int phaseIdx,
- Scalar temperature,
- Scalar pressure,
- const FluidState &fluidState)
+ static Scalar molarMass(int compIdx)
{
- if (phaseIdx == 0)
- return 1.03e-9; // in [kg /microm^3]
- // Umwandlung in Pascal notwendig -> density*10^6
- // will man wieder mit kg/m^3 rechnen muss g auch wieder geändert werden
-
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ switch (compIdx) {
+ case isFluidIdx:
+ // TODO: this is just a rough guess
+ return 22e-3; // [kg/mol]
+ case trailIdx:
+ return 567e-3; // [kg/mol]
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
+
/*!
* \brief Given all mole fractions in a phase, return the phase
* density [kg/m^3].
*/
template <class FluidState>
- static Scalar molarDensity(int phaseIdx,
+ static Scalar phaseDensity(int phaseIdx,
Scalar temperature,
Scalar pressure,
const FluidState &fluidState)
{
if (phaseIdx == 0)
- return 3.035e-16; // in [mol/microm^3] = 303.5 mol/m^3
+ return 1.03e3; // in [kg /m^3]
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
- * \brief Return the viscosity of a phase.
+ * \brief Return the dynamic viscosity of a phase.
*/
template <class FluidState>
static Scalar phaseViscosity(int phaseIdx,
@@ -110,7 +112,7 @@ public:
const FluidState &fluidState)
{
if (phaseIdx == 0)
- return 0.00069152; // in [Pas]
+ return 0.00069152; // in [Pa*s]
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
@@ -127,7 +129,8 @@ public:
Scalar pressure,
const FluidState &fluidState)
{
- return 0.088786695; // in [microm^2/s] = 3.7378e-12 m^2/s
+ // 3.7378e-12
+ return 8.8786695-14; // in [m^2/s]
}
};
diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh
index 421f15fd54e61e43c6b6a849c304709f29f2820d..e9c2a41c0d0ddeb98afa51bbe5591832a59d919a 100644
--- a/dumux/material/spatialparameters/boxspatialparameters.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters.hh
@@ -29,8 +29,12 @@
#include <dune/common/fmatrix.hh>
-namespace Dumux
-{
+namespace Dumux {
+// forward declation of property tags
+namespace Properties {
+NEW_PROP_TAG(SpatialParameters);
+};
+
/**
* \brief The base class for spatial parameters of problems using the
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 7c915f049b4047ddbe889dceb68087764066e780..cee80557a4dd3b3278e908b40739aca779ac69fe 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -25,29 +25,40 @@
#include <dumux/common/exceptions.hh>
-#include <dune/istl/overlappingschwarz.hh>
-//#include <dune/istl/schwarz.hh>
-#include <dune/istl/preconditioners.hh>
-#include <dune/istl/solvers.hh>
-//#include "dune/istl/owneroverlapcopy.hh"
-
-#include <dune/istl/io.hh>
-
-#include <dune/common/mpihelper.hh>
-
-#include <iostream>
-#include <boost/format.hpp>
+#include <queue> // for std::priority_queue
#include <dumux/common/pardiso.hh>
+#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/common/pdelabpreconditioner.hh>
-#include <dune/pdelab/backend/istlsolverbackend.hh>
+
namespace Dumux
{
namespace Properties
{
+//! specifies the implementation of the newton controller
+NEW_PROP_TAG(NewtonController);
+
+//! specifies the type of the actual newton method
+NEW_PROP_TAG(NewtonMethod);
+
+//! specifies the type of a solution
+NEW_PROP_TAG(SolutionVector);
+
+//! specifies the type of a vector of primary variables at an DOF
+NEW_PROP_TAG(PrimaryVariables);
+
+//! specifies the type of a global jacobian matrix
+NEW_PROP_TAG(JacobianMatrix);
+
+//! specifies the type of the jacobian matrix assembler
+NEW_PROP_TAG(JacobianAssembler);
+
+//! specifies the type of the time manager
+NEW_PROP_TAG(TimeManager);
+
//! specifies the verbosity of the linear solver (by default it is 0,
//! i.e. it doesn't print anything)
NEW_PROP_TAG(NewtonLinearSolverVerbosity);
@@ -56,6 +67,14 @@ NEW_PROP_TAG(NewtonLinearSolverVerbosity);
//! gets written out to disk for every newton iteration (default is false)
NEW_PROP_TAG(NewtonWriteConvergence);
+//! specifies whether time step size should be increased during the
+//! newton methods first few iterations
+NEW_PROP_TAG(EnableTimeStepRampUp);
+
+//! specifies whether the jacobian matrix should only be reassembled
+//! if the current solution deviates too much from the evaluation point
+NEW_PROP_TAG(EnablePartialReassemble);
+
//! specifies whether the update should be done using the line search
//! method instead of the "raw" newton method. whether this property
//! has any effect depends on wether the line search method is
@@ -185,13 +204,37 @@ class NewtonController
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef NewtonConvergenceWriter<TypeTag, GET_PROP_VALUE(TypeTag, PTAG(NewtonWriteConvergence))> ConvergenceWriter;
enum { enableTimeStepRampUp = GET_PROP_VALUE(TypeTag, PTAG(EnableTimeStepRampUp)) };
+ enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };
+ enum { Red = JacobianAssembler::Red };
+
+ // class to keep track of the most offending vertices in a way
+ // compatible with std::priority_queue
+ class VertexError
+ {
+ public:
+ VertexError(int idx, Scalar err)
+ {
+ idx_ = idx;
+ err_ = err;
+ }
+
+ int index() const
+ { return idx_; }
+
+ bool operator<(const VertexError &a) const
+ { return a.err_ < err_; }
+
+ private:
+ int idx_;
+ Scalar err_;
+ };
public:
NewtonController()
@@ -199,13 +242,21 @@ public:
convergenceWriter_(asImp_())
{
numSteps_ = 0;
- rampUpSteps_ = 4;
- // maximum acceptable difference of any primary variable
- // between two iterations for convergence
- setRelTolerance(1e-7);
- setTargetSteps(9);
- setMaxSteps(15);
+ this->setRelTolerance(1e-8);
+ this->rampUpSteps_ = 0;
+
+ if (enableTimeStepRampUp) {
+ this->rampUpSteps_ = 9;
+
+ // the ramp-up steps are not counting
+ this->setTargetSteps(10);
+ this->setMaxSteps(12);
+ }
+ else {
+ this->setTargetSteps(10);
+ this->setMaxSteps(18);
+ }
};
/*!
@@ -229,22 +280,35 @@ public:
void setMaxSteps(int maxSteps)
{ maxSteps_ = maxSteps; }
+ /*!
+ * \brief Returns the number of iterations used for the time step
+ * ramp-up.
+ */
+ Scalar rampUpSteps() const
+ { return enableTimeStepRampUp?rampUpSteps_:0; }
+
+ /*!
+ * \brief Returns whether the time-step ramp-up is still happening
+ */
+ bool inRampUp() const
+ { return numSteps_ < rampUpSteps(); }
+
/*!
* \brief Returns true if another iteration should be done.
*/
bool newtonProceed(const SolutionVector &u)
{
- if (numSteps_ < rampUpSteps_)
+ if (numSteps_ < rampUpSteps() + 2)
return true; // we always do at least two iterations
- else if (numSteps_ >= maxSteps_) {
+ else if (asImp_().newtonConverged())
+ return false; // we are below the desired tolerance
+ else if (numSteps_ >= rampUpSteps() + maxSteps_) {
// we have exceeded the allowed number of steps. if the
// relative error was reduced by a factor of at least 4,
// we proceed even if we are above the maximum number of
// steps
- return error_*4.0 < lastError_ && !asImp_().newtonConverged();
+ return error_*4.0 < lastError_;
}
- else if (asImp_().newtonConverged())
- return false; // we are below the desired tolerance
return true;
}
@@ -255,7 +319,9 @@ public:
*/
bool newtonConverged() const
{
- return error_ <= tolerance_;
+ return
+ error_ <= tolerance_ &&
+ model_().jacobianAssembler().reassembleTolerance() <= tolerance_/2;
}
/*!
@@ -267,16 +333,19 @@ public:
method_ = &method;
numSteps_ = 0;
- if (enableTimeStepRampUp) {
- destTimeStepSize_ = timeManager_().timeStepSize();
-
- // reduce initial time step size for ramp-up
- timeManager_().setTimeStepSize(destTimeStepSize_ / rampUpSteps_);
+ model_().jacobianAssembler().reassembleAll();
- // reduce the tolerance for partial reassembly during the
- // ramp up
- finalTolerance_ = tolerance_;
- tolerance_ = tolerance_*1e8;
+ dtInitial_ = timeManager_().timeStepSize();
+ if (enableTimeStepRampUp) {
+ rampUpDelta_ =
+ timeManager_().timeStepSize()
+ /
+ rampUpSteps()
+ *
+ 2;
+
+ // reduce initial time step size for ramp-up.
+ timeManager_().setTimeStepSize(rampUpDelta_);
}
convergenceWriter_.beginTimestep();
@@ -308,34 +377,26 @@ public:
error_ = 0;
int idxI = -1;
- int idxJ = -1;
+ int aboveTol = 0;
for (int i = 0; i < int(uOld.size()); ++i) {
- bool needReassemble = false;
- for (int j = 0; j < FV::size; ++j) {
- Scalar tmp
- =
- std::abs(deltaU[i][j])
- / std::max(std::abs(uOld[i][j]), Scalar(1e-4));
- if (tmp > tolerance_ / 10) {
- needReassemble = true;
- }
- if (tmp > error_)
- {
- idxI = i;
- idxJ = j;
- error_ = tmp;
- }
- }
+ PrimaryVariables uNewI = uOld[i];
+ uNewI -= deltaU[i];
+ Scalar vertErr =
+ model_().relativeErrorVertex(i,
+ uOld[i],
+ uNewI);
- model_().jacobianAssembler().markVertexRed(i,
- needReassemble);
+ if (vertErr > tolerance_)
+ ++aboveTol;
+ if (vertErr > error_) {
+ idxI = i;
+ error_ = vertErr;
+ }
}
- model_().jacobianAssembler().computeColors();
error_ = gridView_().comm().max(error_);
}
-
/*!
* \brief Solve the linear system of equations \f$ \mathbf{A}x - b
* = 0\f$.
@@ -401,11 +462,30 @@ public:
writeConvergence_(uOld, deltaU);
newtonUpdateRelError(uOld, deltaU);
-
+
deltaU *= -1;
deltaU += uOld;
- }
+ // compute the vertex and element colors for partial
+ // reassembly
+ if (enablePartialReassemble) {
+ Scalar maxDelta = 0;
+ for (int i = 0; i < int(uOld.size()); ++i) {
+ const PrimaryVariables &uEval = this->model_().jacobianAssembler().evalPoint()[i];
+ const PrimaryVariables &uSol = this->model_().curSol()[i];
+ Scalar tmp =
+ model_().relativeErrorVertex(i,
+ uEval,
+ uSol);
+ maxDelta = std::max(tmp, maxDelta);
+ }
+
+ Scalar reassembleTol = std::max(maxDelta/10, this->tolerance_/5);
+ if (error_ < 10*tolerance_)
+ reassembleTol = tolerance_/5;
+ this->model_().jacobianAssembler().computeColors(reassembleTol);
+ }
+ }
/*!
* \brief Indicates that one newton iteration was finished.
@@ -413,21 +493,20 @@ public:
void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
{
++numSteps_;
-
- if (enableTimeStepRampUp) {
- if (numSteps_ < rampUpSteps_) {
- // increase time step size
- Scalar dt = destTimeStepSize_ * (numSteps_ + 1) / rampUpSteps_;
- timeManager_().setTimeStepSize(dt);
- }
- else // if we're not in the ramp-up reduce the target
- // tolerance
- tolerance_ = finalTolerance_;
+
+ Scalar realError = error_;
+ if (inRampUp() && error_ < 1.0) {
+ // change time step size
+ Scalar dt = timeManager_().timeStepSize();
+ dt += rampUpDelta_;
+ timeManager_().setTimeStepSize(dt);
+
+ endIterMsg() << ", dt=" << timeManager_().timeStepSize() << ", ddt=" << rampUpDelta_;
}
if (verbose())
std::cout << "\rNewton iteration " << numSteps_ << " done: "
- << "error=" << error_ << endIterMsg().str() << "\n";
+ << "error=" << realError << endIterMsg().str() << "\n";
endIterMsgStream_.str("");
}
@@ -436,8 +515,6 @@ public:
*/
void newtonEnd()
{
- if (enableTimeStepRampUp)
- timeManager_().setTimeStepSize(destTimeStepSize_);
convergenceWriter_.endTimestep();
}
@@ -448,6 +525,7 @@ public:
*/
void newtonFail()
{
+ timeManager_().setTimeStepSize(dtInitial_);
numSteps_ = targetSteps_*2;
}
@@ -468,10 +546,12 @@ public:
*/
Scalar suggestTimeStepSize(Scalar oldTimeStep) const
{
+ if (enableTimeStepRampUp)
+ return oldTimeStep;
+
Scalar n = numSteps_;
- if (enableTimeStepRampUp) {
- n -= rampUpSteps_/2;
- }
+ n -= rampUpSteps();
+
// be agressive reducing the timestep size but
// conservative when increasing it. the rationale is
// that we want to avoid failing in the next newton
@@ -538,6 +618,12 @@ protected:
TimeManager &timeManager_()
{ return problem_().timeManager(); }
+ /*!
+ * \brief Returns a reference to the time manager.
+ */
+ const TimeManager &timeManager_() const
+ { return problem_().timeManager(); }
+
/*!
* \brief Returns a reference to the problem.
*/
@@ -588,7 +674,7 @@ protected:
typedef Dumux::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
Solver solver(problem_(), 500, verbosity);
#else
- typedef Dumux::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
+ typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
Solver solver(500, verbosity);
#endif // HAVE_MPI
#endif // HAVE_PARDISO
@@ -618,15 +704,16 @@ protected:
ConvergenceWriter convergenceWriter_;
Scalar tolerance_;
- Scalar finalTolerance_;
Scalar error_;
Scalar lastError_;
// number of iterations for the time-step ramp-up
Scalar rampUpSteps_;
- // time step size after the ramp-up
- Scalar destTimeStepSize_;
+ // the increase of the time step size during the rampup
+ Scalar rampUpDelta_;
+
+ Scalar dtInitial_; // initial time step size
// optimal number of iterations we want to achive
int targetSteps_;
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index 30e733d2f3e93404a3758982ba6bc4e4fda740d7..11e3bc31ec08d503dbe6229a149b65d491b76b93 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -23,15 +23,26 @@
#ifndef DUMUX_NEWTONMETHOD_HH
#define DUMUX_NEWTONMETHOD_HH
-#include <limits>
#include <dumux/common/exceptions.hh>
+#include <dumux/common/propertysystem.hh>
-#include <dumux/io/vtkmultiwriter.hh>
+#include <dune/istl/istlexception.hh>
-#include <dune/istl/io.hh>
+#include <iostream>
namespace Dumux
{
+// forward declaration of property tags
+namespace Properties
+{
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(Model);
+NEW_PROP_TAG(NewtonController);
+NEW_PROP_TAG(SolutionVector);
+NEW_PROP_TAG(JacobianAssembler);
+};
+
/*!
* \brief The algorithmic part of the multi dimensional newton method.
*
diff --git a/test/boxmodels/1p/1pspatialparameters.hh b/test/boxmodels/1p/1pspatialparameters.hh
index 61a7e4135ebfe43569470cdfae3ab817cca5bd2a..f81b93dae4e5a6c43175fea4f19b0e69b53db1bd 100644
--- a/test/boxmodels/1p/1pspatialparameters.hh
+++ b/test/boxmodels/1p/1pspatialparameters.hh
@@ -22,24 +22,18 @@ namespace Dumux
{
/** \todo Please doc me! */
-
template<class TypeTag>
class OnePSpatialParameters : public BoxSpatialParameters<TypeTag>
{
typedef BoxSpatialParameters<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GridView::ctype CoordScalar;
- enum {
- dimWorld=GridView::dimensionworld,
- };
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
public:
-
OnePSpatialParameters(const GridView& gridView)
: ParentType(gridView)
{}
@@ -56,18 +50,12 @@ public:
Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
- {
- return 1e-10;
- }
+ { return 1e-10; }
Scalar porosity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
- {
- return 0.4;
- }
-
-private:
+ { return 0.4; }
};
} // end namespace
diff --git a/test/boxmodels/1p/1ptestproblem.hh b/test/boxmodels/1p/1ptestproblem.hh
index c72209187b07aacb27e7b0e916f47dca1dc6524e..29bbf8215715b791d24056b01b9e4de0b1360e57 100644
--- a/test/boxmodels/1p/1ptestproblem.hh
+++ b/test/boxmodels/1p/1ptestproblem.hh
@@ -47,12 +47,7 @@ public:
};
// Set the grid type
-#if HAVE_UG
-SET_PROP(OnePTestProblem, Grid) { typedef Dune::ALUCubeGrid<3,3> type; };
-#else
-//SET_PROP(OnePTestProblem, Grid) { typedef Dune::SGrid<3, 3> type; };
SET_TYPE_PROP(OnePTestProblem, Grid, Dune::YaspGrid<3>);
-#endif
SET_PROP(OnePTestProblem, LocalFEMSpace)
{
@@ -65,36 +60,15 @@ public:
// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
};
-SET_PROP(OnePTestProblem, NewtonLinearSolverVerbosity)
-{public:
- static const int value = 1;
-};
+SET_INT_PROP(OnePTestProblem, NewtonLinearSolverVerbosity, 0);
// Set the problem property
SET_PROP(OnePTestProblem, Problem)
-{
- typedef Dumux::OnePTestProblem<TTAG(OnePTestProblem)> type;
-};
-
-// Set the wetting phase
-//SET_TYPE_PROP(OnePTestProblem, Fluid, Dumux::Air);
+{ typedef Dumux::OnePTestProblem<TTAG(OnePTestProblem)> type; };
-//! Use the leaf grid view if not defined otherwise
-//SET_PROP(OnePTestProblem, GridView)
-//{
-//private:
-// typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-//
-//public:
-// typedef typename Grid::LevelGridView type;
-//};
-
-// Set the spatial parameters
// Set the spatial parameters
SET_PROP(OnePTestProblem, SpatialParameters)
-{
- typedef Dumux::OnePSpatialParameters<TypeTag> type;
-};
+{ typedef Dumux::OnePSpatialParameters<TypeTag> type; };
// Enable gravity
SET_BOOL_PROP(OnePTestProblem, EnableGravity, true);
@@ -287,7 +261,7 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ //const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
values[pressureIdx] = 1.0e+5;// + 9.81*1.23*(20-globalPos[dim-1]);
}
diff --git a/test/boxmodels/1p2c/grids/test_1p2c.dgf b/test/boxmodels/1p2c/grids/test_1p2c.dgf
index d7715bbd3388f39243de8e8b1af52843a19ac60b..3d48df44204ee9114553557e23e7853a28d195f4 100644
--- a/test/boxmodels/1p2c/grids/test_1p2c.dgf
+++ b/test/boxmodels/1p2c/grids/test_1p2c.dgf
@@ -1,7 +1,7 @@
DGF
Interval
0 0 % first corner
-22 22 % second corner
+22e-3 22e-3 % second corner
44 44 % x-Achse wird 44x unterteilt und y-Achse
#
diff --git a/test/boxmodels/1p2c/tissue_tumor_problem.hh b/test/boxmodels/1p2c/tissue_tumor_problem.hh
index 0fbac25d02d0a2a2c4d8dd2eaa5dde67b3539927..c0910a13fe3cd60985583dbd13c380ae99199ddb 100644
--- a/test/boxmodels/1p2c/tissue_tumor_problem.hh
+++ b/test/boxmodels/1p2c/tissue_tumor_problem.hh
@@ -121,6 +121,10 @@ class TissueTumorProblem : public OnePTwoCBoxProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
// copy some indices for convenience
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
@@ -130,28 +134,39 @@ class TissueTumorProblem : public OnePTwoCBoxProblem<TypeTag>
dimWorld = GridView::dimensionworld,
// indices of the primary variables
- konti = Indices::konti,
- transport = Indices::transport,
+ pressureIdx = Indices::pressureIdx,
+ x1Idx = Indices::x1Idx,
+
+ // indices of the equations
+ contiEqIdx = Indices::contiEqIdx,
+ transEqIdx = Indices::transEqIdx,
};
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
-
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::Intersection Intersection;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
TissueTumorProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
{
+ // calculate the injection volume
+ totalInjectionVolume_ = 0;
+ FVElementGeometry fvGeom;
+ ElementIterator elemIt = gridView.template begin<0>();
+ const ElementIterator endIt = gridView.template end<0>();
+ for (; elemIt != endIt; ++ elemIt) {
+ fvGeom.update(gridView, *elemIt);
+ for (int i = 0; i < fvGeom.numVertices; ++i) {
+ const GlobalPosition &pos = fvGeom.subContVol[i].global;
+ if (inInjectionVolume_(pos))
+ totalInjectionVolume_ += fvGeom.subContVol[i].volume;
+ };
+ }
}
/*!
@@ -242,10 +257,11 @@ public:
//Scalar lambda = (globalPos[1])/height_;
if (globalPos[0] < eps_ ) {
- values[konti] = -3.8676e-8;
- //values[transport] = -4.35064e-10;
+ values[contiEqIdx] = -3.8676e-2; // [kg/(m^2 * s)]
+
+ //values[transEqIdx] = -4.35064e-4; // [mol/(m^2*s)
//Robin-Boundary
- //values[transport] = (*this->model().curSolFunction())[globalIdx][transport];
+ //values[transEqIdx] = (*this->model().curSolFunction())[globalIdx][transEqIdx];
}
}
@@ -270,13 +286,29 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
-
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
values = Scalar(0.0);
-
- if(globalPos[0]>10 && globalPos[0]<12 && globalPos[1]>10 && globalPos[1]<12)
- values[0]= 1.5e-6;
+ if (inInjectionVolume_(globalPos)) {
+ // total volumetric injection rate in ml/h
+ Scalar injRateVol = 0.1;
+ // convert to m^3/s
+ injRateVol *= 1e-6/3600;
+ // total mass injection rate. assume a density of 1030kg/m^3
+ Scalar injRateMass = injRateVol*1030.0;
+
+ // trail concentration in injected fluid in [mol/ml]
+ Scalar trailInjRate = 1e-5;
+ // convert to mol/m^3
+ trailInjRate *= 1e6;
+ // convert to mol/s
+ trailInjRate *= injRateVol;
+
+ // source term of the total mass
+ values[contiEqIdx] = injRateMass / totalInjectionVolume_; // [kg/(s*m^3)]
+ values[transEqIdx] = trailInjRate / totalInjectionVolume_; // [mol/(s*m^3)]
+ }
}
/*!
@@ -299,16 +331,21 @@ public:
// \}
private:
+ bool inInjectionVolume_(const GlobalPosition &globalPos) const
+ {
+ return
+ 10e-3 < globalPos[0] && globalPos[0] < 12e-3 &&
+ 10e-3 < globalPos[1] && globalPos[1] < 12e-3;
+ };
// the internal method for the initial condition
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
-
- values[konti] = -1067; //initial condition for the pressure
- values[transport] = 1.1249e-8; //initial condition for the molefraction
-
+ values[pressureIdx] = 0; //initial condition for the pressure
+ values[x1Idx] = 0; //initial condition for the trail molefraction
}
+ Scalar totalInjectionVolume_;
static const Scalar eps_ = 1e-6;
}; //end namespace
}
diff --git a/test/boxmodels/1p2c/tissue_tumor_spatialparameters.hh b/test/boxmodels/1p2c/tissue_tumor_spatialparameters.hh
index 6e1f11d82aa71b6e4e24a4f328584804afe9c7d2..9f4f2419421a514f485cbda617b8de49729a3560 100644
--- a/test/boxmodels/1p2c/tissue_tumor_spatialparameters.hh
+++ b/test/boxmodels/1p2c/tissue_tumor_spatialparameters.hh
@@ -169,8 +169,8 @@ public:
private:
bool isTumor_(const Dune::FieldVector<CoordScalar,dim>& globalPos) const
{
- if(globalPos[0] > 0.99 && globalPos[0] < 2.01 &&
- globalPos[1] > 0.99 && globalPos[1] < 3.01)
+ if(10e-3 < globalPos[0] && globalPos[0] < 15e-3 &&
+ 10e-3 < globalPos[1] && globalPos[1] < 15e-3)
return true;
return false;
}
diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh
index d970a81e54fa6df3946e19a4569603bb88a8e7c7..cbb694640f5f87f854e22a3e8d6035adcc2787e5 100644
--- a/test/boxmodels/2p/lensspatialparameters.hh
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -24,6 +24,8 @@
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/boxmodels/2p/2pmodel.hh>
+
/**
* @file
* @brief Class for defining spatial parameters