From 0606de03755c805dfc8f0d7af510680b2107aad1 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Thu, 26 Aug 2010 15:49:56 +0000
Subject: [PATCH] various cleanups
- make all headers in dumux/common and dumux/boxmodels/common self sustained,
i.e. that they can be included without any preconditions
- remove most irrelevant includes in the common code
- bring the 1p and 1p2c box model up to stable standard (accessor functions in
the volume and flux variable classes, make the 1p2c model work for
compressible fluids, various fixes the 1p2c test problem)
- change the isfluid-trail fluid system to standard SI units
- play around with timestep control if timestep ramp-up is enabled
(problem is still not solved, though)
- split the *properties.hh headers into the declaration of the properties and
definition of defaults for the common part of the box models as well as for
the 1p and 1p2c models
- put the indices structure for the 1p and 1p2c models into separate files
- introduce a base class for volume variables of the box models
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4180 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/1p/1pfluxvariables.hh | 114 ++++---
dumux/boxmodels/1p/1pindices.hh | 37 +++
dumux/boxmodels/1p/1plocalresidual.hh | 32 +-
dumux/boxmodels/1p/1pmodel.hh | 4 +-
dumux/boxmodels/1p/1pproperties.hh | 48 +--
dumux/boxmodels/1p/1ppropertydefaults.hh | 68 ++++
dumux/boxmodels/1p/1pvolumevariables.hh | 79 +++--
dumux/boxmodels/1p2c/1p2cfluidstate.hh | 115 ++++++-
dumux/boxmodels/1p2c/1p2cfluxvariables.hh | 297 +++++++++---------
dumux/boxmodels/1p2c/1p2cindices.hh | 48 +++
dumux/boxmodels/1p2c/1p2clocalresidual.hh | 136 ++++----
dumux/boxmodels/1p2c/1p2cmodel.hh | 25 +-
dumux/boxmodels/1p2c/1p2cproperties.hh | 59 +---
dumux/boxmodels/1p2c/1p2cpropertydefaults.hh | 79 +++++
dumux/boxmodels/1p2c/1p2cvolumevariables.hh | 164 +++++++---
dumux/boxmodels/2p/2pfluxvariables.hh | 138 +++-----
dumux/boxmodels/2p/2plocalresidual.hh | 53 ++--
dumux/boxmodels/2p/2pmodel.hh | 20 +-
dumux/boxmodels/2p/2pnewtoncontroller.hh | 135 --------
dumux/boxmodels/2p/2pproblem.hh | 4 +-
dumux/boxmodels/2p/2pproperties.hh | 173 +---------
dumux/boxmodels/2p/2ppropertydefaults.hh | 180 +++++++++++
dumux/boxmodels/2p/2pvolumevariables.hh | 21 +-
dumux/boxmodels/2p2c/2p2cfluidstate.hh | 4 +-
.../common/boxelementboundarytypes.hh | 6 +-
.../common/boxelementvolumevariables.hh | 3 -
.../boxmodels/common/boxfvelementgeometry.hh | 1 -
dumux/boxmodels/common/boxlocaljacobian.hh | 18 +-
dumux/boxmodels/common/boxlocalresidual.hh | 16 +-
dumux/boxmodels/common/boxmodel.hh | 64 +++-
dumux/boxmodels/common/boxproblem.hh | 2 -
dumux/boxmodels/common/boxproperties.hh | 274 ----------------
dumux/boxmodels/common/boxpropertydefaults.hh | 277 ++++++++++++++++
dumux/boxmodels/common/boxvolumevariables.hh | 116 +++++++
dumux/boxmodels/common/pdelabboxassembler.hh | 100 ++++--
.../common/pdelabboxistlvectorbackend.hh | 6 +-
.../common/pdelabboxlocaloperator.hh | 11 +-
dumux/common/boundarytypes.hh | 3 -
dumux/common/math.hh | 5 +-
dumux/common/pdelabpreconditioner.hh | 27 +-
dumux/common/propertysystem.hh | 2 -
dumux/common/start.hh | 18 +-
dumux/common/timemanager.hh | 27 +-
dumux/material/components/brine.hh | 4 -
dumux/material/components/ch4.hh | 1 -
dumux/material/components/h2.hh | 1 -
dumux/material/components/h2o.hh | 1 -
dumux/material/components/n2.hh | 1 -
dumux/material/components/o2.hh | 1 -
dumux/material/components/oil.hh | 1 -
dumux/material/components/simpleco2.hh | 2 -
dumux/material/components/simplednapl.hh | 3 -
dumux/material/components/simpleh2o.hh | 2 -
.../material/components/tabulatedcomponent.hh | 4 +
dumux/material/components/unit.hh | 1 -
.../fluidmatrixinteractions/2p/brookscorey.hh | 2 -
.../2p/linearmaterial.hh | 2 +
.../2p/regularizedbrookscorey.hh | 3 -
.../2p/regularizedbrookscoreyparams.hh | 15 +-
.../2p/regularizedlinearmaterial.hh | 4 -
.../2p/regularizedvangenuchten.hh | 2 -
.../2p/regularizedvangenuchtenparams.hh | 1 -
.../2p/vangenuchten.hh | 2 -
dumux/material/fluidstate.hh | 2 +
dumux/material/fluidsystems/2p_system.hh | 14 +-
dumux/material/fluidsystems/h2o_n2_system.hh | 6 +
.../fluidsystems/isfluid_trail_system.hh | 49 +--
.../spatialparameters/boxspatialparameters.hh | 8 +-
dumux/nonlinear/newtoncontroller.hh | 245 ++++++++++-----
dumux/nonlinear/newtonmethod.hh | 17 +-
test/boxmodels/1p/1pspatialparameters.hh | 22 +-
test/boxmodels/1p/1ptestproblem.hh | 34 +-
test/boxmodels/1p2c/grids/test_1p2c.dgf | 2 +-
test/boxmodels/1p2c/tissue_tumor_problem.hh | 81 +++--
.../1p2c/tissue_tumor_spatialparameters.hh | 4 +-
test/boxmodels/2p/lensspatialparameters.hh | 2 +
76 files changed, 2033 insertions(+), 1515 deletions(-)
create mode 100644 dumux/boxmodels/1p/1pindices.hh
create mode 100644 dumux/boxmodels/1p/1ppropertydefaults.hh
create mode 100644 dumux/boxmodels/1p2c/1p2cindices.hh
create mode 100644 dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
delete mode 100644 dumux/boxmodels/2p/2pnewtoncontroller.hh
create mode 100644 dumux/boxmodels/2p/2ppropertydefaults.hh
create mode 100644 dumux/boxmodels/common/boxpropertydefaults.hh
create mode 100644 dumux/boxmodels/common/boxvolumevariables.hh
diff --git a/dumux/boxmodels/1p/1pfluxvariables.hh b/dumux/boxmodels/1p/1pfluxvariables.hh
index e797b7af56..064a860cc5 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 0000000000..67fb766702
--- /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 220be3cdda..ac7a86762a 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 7ff7847337..5519535448 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 29936e46a6..bc6aa89b7e 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 0000000000..0e2499721e
--- /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 3e3b5edf6c..5e5d4530e1 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 a545fafc26..81feb9b10d 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 f1ce7eb944..52e67cccfd 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 0000000000..28e9bab32c
--- /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 cf234883f6..0a4ab8bd23 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 7641412ad6..6c2c4cbebe 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 3156cc8b03..836cef8364 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 0000000000..94c7a67903
--- /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 97c01cd0c0..7fc3ecd05b 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 94665217c4..18eda48ac4 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 d0dc353cc6..3080de8d92 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 d73bcdf410..1b203fc388 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 0e65a90bd1..0000000000
--- 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 8e103203ea..3207088685 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 dd5475a640..cd89088c00 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 0000000000..e769025d6b
--- /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 bc03a0b386..0d1dcf1607 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 6e342d0ba6..51ee03f639 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 08815ae866..ecf24523b4 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 012028c3ee..5e30507c23 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 8ac1925afe..3225ce0aba 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 a0e9775d48..fc77fb39ec 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 e2d26a4021..13f530a21b 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 4044bcbf9a..0099d56858 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 79c28fe11b..fdd496d1d4 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 5517f9cc5c..f064db2d8a 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 0000000000..2ecbbfafac
--- /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 0000000000..f858ad92f4
--- /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 aee42fe7cf..f3f87d80eb 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 20fd2e9375..b0031deb85 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 e60f24c99f..3ae4b10e9e 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 af4994c0f6..d4f4762518 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 86d71e4c6b..27949ece11 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 211d7384c5..ebc69260bc 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 c3a5e7de4e..71c973fb2c 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 9f5df7b5a9..dff01bc06f 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 ac118299e1..60304f11a3 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 cad93e6797..8624c071f3 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 f0cf690415..6cf7beca56 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 30f1d1ea49..103cf36ee2 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 d7af489398..d2fee8c2d3 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 0640d6c781..46b3c33c1c 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 73be4d2a3e..72700adc73 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 5e176b5d90..135a44beab 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 a6e8675bf7..67e2559871 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 2ec18d3b36..ba94b1e3c9 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 20df9f76ac..f8f09cf523 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 e754f6cd50..1d133b5d77 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 333b290fd6..5b071eaad3 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 2bf6fa2f0d..7b6eb70abe 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 7e96ef2b33..6a8bf280eb 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 e5b4b664b5..c38a172ffe 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 c65b2869e3..65b520a6dc 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 be33b5fb71..926579900f 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 80233098b2..91eccad9ae 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 40955fbe6b..8a3c0091c7 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 88603f05c3..d44fb9350f 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 3e0dc01f91..771de886ec 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 eea985f1ed..b9670cd2d9 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 83a80a3460..a9ec010b5e 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 83fe66b74c..d9ffb4acde 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 421f15fd54..e9c2a41c0d 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 7c915f049b..cee80557a4 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 30e733d2f3..11e3bc31ec 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 61a7e4135e..f81b93dae4 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 c72209187b..29bbf82157 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 d7715bbd33..3d48df4420 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 0fbac25d02..c0910a13fe 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 6e1f11d82a..9f4f241942 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 d970a81e54..cbb694640f 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
--
GitLab