diff --git a/test/implicit/2p/Makefile.am b/test/implicit/2p/Makefile.am
index 890f94bff76fcfb03658b6cf730111768c650b9d..bb63e52dcf69a53e0c961616e415e4bd79b8d24e 100644
--- a/test/implicit/2p/Makefile.am
+++ b/test/implicit/2p/Makefile.am
@@ -1,9 +1,9 @@
-check_PROGRAMS = test_2p test_cc2p
+check_PROGRAMS = test_box2p test_cc2p
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_2p_SOURCES = test_2p.cc
+test_box2p_SOURCES = test_box2p.cc
test_cc2p_SOURCES = test_cc2p.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/2p/lensccproblem.hh b/test/implicit/2p/lensccproblem.hh
deleted file mode 100644
index f8dff793dd051b12be157ecadc5dd4f7452b804f..0000000000000000000000000000000000000000
--- a/test/implicit/2p/lensccproblem.hh
+++ /dev/null
@@ -1,392 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2007-2008 by Klaus Mosthaf *
- * Copyright (C) 2007-2008 by Bernd Flemisch *
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Institute for Modelling Hydraulic and Environmental Systems *
- * 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Soil contamination problem where DNAPL infiltrates a fully
- * water saturated medium.
- */
-
-#ifndef DUMUX_LENSPROBLEM_HH
-#define DUMUX_LENSPROBLEM_HH
-
-#if HAVE_UG
-#include <dune/grid/uggrid.hh>
-#endif
-
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/sgrid.hh>
-
-#include <dumux/implicit/cellcentered/ccpropertydefaults.hh>
-#include <dumux/implicit/cellcentered/porousmediaccproblem.hh>
-#include <dumux/implicit/2p/2pmodel.hh>
-#include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/components/dnapl.hh>
-
-#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
-
-#include "lensccspatialparams.hh"
-
-namespace Dumux
-{
-
-template <class TypeTag>
-class LensProblem;
-
-//////////
-// Specify the properties for the lens problem
-//////////
-namespace Properties
-{
-NEW_TYPE_TAG(LensProblem, INHERITS_FROM(CCTwoP, LensSpatialParams));
-
-// Set the grid type
-#if 0//HAVE_UG
-SET_TYPE_PROP(LensProblem, Grid, Dune::UGGrid<2>);
-#else
-SET_TYPE_PROP(LensProblem, Grid, Dune::YaspGrid<2>);
-#endif
-
-// Set the problem property
-SET_TYPE_PROP(LensProblem, Problem, Dumux::LensProblem<TypeTag>);
-
-// TODO: remove this macro switch
-#if 1
-// Set the wetting phase
-SET_PROP(LensProblem, WettingPhase)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-public:
- typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
-};
-
-// Set the non-wetting phase
-SET_PROP(LensProblem, NonwettingPhase)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-public:
- typedef Dumux::LiquidPhase<Scalar, Dumux::DNAPL<Scalar> > type;
-};
-#else
-// OR: set the fluid system
-SET_TYPE_PROP(LensProblem, FluidSystem, H2ON2FluidSystem<TypeTag>);
-#endif
-
-// Enable partial reassembly of the jacobian matrix?
-SET_BOOL_PROP(LensProblem, ImplicitEnablePartialReassemble, true);
-
-// Enable reuse of jacobian matrices?
-SET_BOOL_PROP(LensProblem, ImplicitEnableJacobianRecycling, true);
-
-// Write the solutions of individual newton iterations?
-SET_BOOL_PROP(LensProblem, NewtonWriteConvergence, false);
-
-// Use forward differences instead of central differences
-SET_INT_PROP(LensProblem, ImplicitNumericDifferenceMethod, +1);
-
-// Linear solver settings
-SET_TYPE_PROP(LensProblem, LinearSolver, Dumux::BoxBiCGStabILU0Solver<TypeTag> );
-SET_INT_PROP(LensProblem, LinearSolverVerbosity, 0);
-SET_INT_PROP(LensProblem, LinearSolverPreconditionerIterations, 1);
-SET_SCALAR_PROP(LensProblem, LinearSolverPreconditionerRelaxation, 1.0);
-
-// Enable gravity
-SET_BOOL_PROP(LensProblem, ProblemEnableGravity, true);
-}
-
-/*!
- * \ingroup TwoPCCModel
- * \ingroup CCTestProblems
- * \brief Soil contamination problem where DNAPL infiltrates a fully
- * water saturated medium.
- *
- * The domain is sized 6m times 4m and features a rectangular lens
- * with low permeablility which spans from (1 m , 2 m) to (4 m, 3 m)
- * and is surrounded by a medium with higher permability. Note that
- * this problem is discretized using only two dimensions, so from the
- * point of view of the two-phase model, the depth of the domain
- * implicitly is 1 m everywhere.
- *
- * On the top and the bottom of the domain neumann boundary conditions
- * are used, while dirichlet conditions apply on the left and right
- * boundaries.
- *
- * DNAPL is injected at the top boundary from 3m to 4m at a rate of
- * 0.04 kg/(s m^2), the remaining neumann boundaries are no-flow
- * boundaries.
- *
- * The dirichlet boundaries on the left boundary is the hydrostatic
- * pressure scaled by a factor of 1.125, while on the right side it is
- * just the hydrostatic pressure. The DNAPL saturation on both sides
- * is zero.
- *
- * This problem uses the \ref TwoPCCModel.
- *
- * This problem should typically be simulated until \f$t_{\text{end}}
- * \approx 20\,000\;s\f$ is reached. A good choice for the initial time step
- * size is \f$t_{\text{inital}} = 250\;s\f$.
- *
- * To run the simulation execute the following line in shell:
- * <tt>./test_2p -parameterFile test_2p.input</tt>
- */
-template <class TypeTag >
-class LensProblem : public PorousMediaCCProblem<TypeTag>
-{
- typedef PorousMediaCCProblem<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
- typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
-
- enum {
-
- // primary variable indices
- pwIdx = Indices::pwIdx,
- SnIdx = Indices::SnIdx,
-
- // equation indices
- contiNEqIdx = Indices::contiNEqIdx,
-
- // phase indices
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
-
-
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
-
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-public:
- /*!
- * \brief The constructor
- *
- * \param timeManager The time manager
- * \param gridView The grid view
- */
- LensProblem(TimeManager &timeManager,
- const GridView &gridView)
- : ParentType(timeManager, gridView)
- {
- eps_ = 3e-6;
- temperature_ = 273.15 + 20; // -> 20°C
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const char *name() const
- { return "lens"; }
-
- /*!
- * \brief Called directly after the time integration.
- */
- void postTimeStep()
- {
- // Calculate storage terms
- PrimaryVariables storage;
- this->model().globalStorage(storage);
-
- // Write mass balance information for rank 0
- if (this->gridView().comm().rank() == 0) {
- std::cout<<"Storage: " << storage << std::endl;
- }
- }
-
- /*!
- * \brief Returns the temperature within the domain.
- *
- * This problem assumes a uniform temperature of 10 degrees Celsius.
- */
- Scalar temperature() const
- { return temperature_; };
-
-
- void sourceAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- values = 0;
- }
-
- // \}
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary control volume.
- *
- * \param values The boundary types for the conservation equations
- * \param globalPos The position of the center of the finite volume
- */
- void boundaryTypesAtPos(BoundaryTypes &values,
- const GlobalPosition &globalPos) const
- {
- if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) {
- values.setAllDirichlet();
- }
- else {
- values.setAllNeumann();
- }
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * control volume.
- *
- * \param values The dirichlet values for the primary variables
- * \param globalPos The center of the finite volume which ought to be set.
- *
- * For this method, the \a values parameter stores primary variables.
- */
- void dirichletAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
- fluidState.setTemperature(temperature_);
- fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
- fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
-
- Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
-
- Scalar height = 4.0;//this->bboxMax()[1] - this->bboxMin()[1];
- Scalar depth = 4.0 - globalPos[1];//this->bboxMax()[1] - globalPos[1];
- Scalar alpha = 1 + 1.5/height;
- Scalar factor = (6.0*alpha + (1.0 - alpha)*globalPos[0])/6.0;
-
- // hydrostatic pressure scaled by alpha
- values[pwIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
- values[SnIdx] = 0.0;
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
- * \param globalPos The position of the integration point of the boundary segment.
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- void neumannAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- values = 0.0;
- if (onInlet_(globalPos)) {
- values[contiNEqIdx] = -0.04; // kg / (m * s)
- }
- }
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The initial values for the primary variables
- * \param globalPos The center of the finite volume which ought to be set.
- *
- * For this method, the \a values parameter stores primary
- * variables.
- */
- void initialAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- Scalar depth = 4.0/*this->bboxMax()[1]*/ - globalPos[1];
-
- typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
- fluidState.setTemperature(temperature_);
- fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
- fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
-
- Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
-
- // hydrostatic pressure
- values[pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
- values[SnIdx] = 0.0;
- }
- // \}
-
-private:
-
- bool onLeftBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[0] < this->bboxMin()[0] + eps_;
- }
-
- bool onRightBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[0] > this->bboxMax()[0] - eps_;
- }
-
- bool onLowerBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[1] < this->bboxMin()[1] + eps_;
- }
-
- bool onUpperBoundary_(const GlobalPosition &globalPos) const
- {
- return globalPos[1] > this->bboxMax()[1] - eps_;
- }
-
- bool onInlet_(const GlobalPosition &globalPos) const
- {
- Scalar width = 6.0;//this->bboxMax()[0] - this->bboxMin()[0];
- Scalar lambda = (6.0/*this->bboxMax()[0]*/ - globalPos[0])/width;
- return (globalPos[1] > 4.0 - eps_) && 0.5 < lambda && lambda < 2.0/3.0;
- }
-
- Scalar temperature_;
- Scalar eps_;
-};
-} //end namespace
-
-#endif
diff --git a/test/implicit/2p/lensccspatialparams.hh b/test/implicit/2p/lensccspatialparams.hh
deleted file mode 100644
index 73eed21cb273c35ec99a01098cc8def55828068f..0000000000000000000000000000000000000000
--- a/test/implicit/2p/lensccspatialparams.hh
+++ /dev/null
@@ -1,220 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010 by Markus Wolff *
- * Copyright (C) 2007-2008 by Klaus Mosthaf *
- * Copyright (C) 2007-2008 by Bernd Flemisch *
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Institute for Modelling Hydraulic and Environmental Systems *
- * 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief The spatial parameters for the LensProblem which uses the
- * twophase box model
- */
-#ifndef DUMUX_LENS_SPATIAL_PARAMS_HH
-#define DUMUX_LENS_SPATIAL_PARAMS_HH
-
-#include <dumux/material/spatialparams/boxspatialparams.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-
-#include <dumux/implicit/2p/2pmodel.hh>
-
-namespace Dumux
-{
-
-//forward declaration
-template<class TypeTag>
-class LensSpatialParams;
-
-namespace Properties
-{
-// The spatial parameters TypeTag
-NEW_TYPE_TAG(LensSpatialParams);
-
-// Set the spatial parameters
-SET_TYPE_PROP(LensSpatialParams, SpatialParams, Dumux::LensSpatialParams<TypeTag>);
-
-// Set the material Law
-SET_PROP(LensSpatialParams, MaterialLaw)
-{
-private:
- // define the material law which is parameterized by effective
- // saturations
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef RegularizedVanGenuchten<Scalar> EffectiveLaw;
-public:
- // define the material law parameterized by absolute saturations
- typedef EffToAbsLaw<EffectiveLaw> type;
-};
-}
-/*!
- * \ingroup TwoPCCModel
- * \ingroup CCTestProblems
- * \brief The spatial parameters for the LensProblem which uses the
- * twophase box model
- */
-template<class TypeTag>
-class LensSpatialParams : public BoxSpatialParams<TypeTag>
-{
- typedef BoxSpatialParams<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename Grid::ctype CoordScalar;
-
- enum {
- dim=GridView::dimension,
- dimWorld=GridView::dimensionworld
- };
-
- typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
-public:
- //get the material law from the property system
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename MaterialLaw::Params MaterialLawParams;
-
- LensSpatialParams(const GridView& gridView)
- : ParentType(gridView)
- {
- try
- {
- lensLowerLeft_[0] = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.LensLowerLeftX);
- lensLowerLeft_[1] = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.LensLowerLeftY);
- lensUpperRight_[0] = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.LensUpperRightX);
- lensUpperRight_[1] = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.LensUpperRightY);
- }
- catch (Dumux::ParameterException &e) {
- std::cerr << e << ". Abort!\n";
- exit(1) ;
- }
- catch (...) {
- std::cerr << "Unknown exception thrown!\n";
- exit(1);
- }
-
- // residual saturations
- lensMaterialParams_.setSwr(0.18);
- lensMaterialParams_.setSnr(0.0);
- outerMaterialParams_.setSwr(0.05);
- outerMaterialParams_.setSnr(0.0);
-
- // parameters for the Van Genuchten law
- // alpha and n
- lensMaterialParams_.setVgAlpha(0.00045);
- lensMaterialParams_.setVgN(7.3);
- outerMaterialParams_.setVgAlpha(0.0037);
- outerMaterialParams_.setVgN(4.7);
-
- // parameters for the linear law
- // minimum and maximum pressures
- // lensMaterialParams_.setEntryPC(0);
-// outerMaterialParams_.setEntryPC(0);
-// lensMaterialParams_.setMaxPC(0);
-// outerMaterialParams_.setMaxPC(0);
-
- lensK_ = 9.05e-12;
- outerK_ = 4.6e-10;
- }
-
- /*!
- * \brief Intrinsic permability
- *
- * \param element The current element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume.
- * \return Intrinsic permeability
- */
- Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- if (scvIdx > 0)
- return intrinsicPermeability(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
-
- const GlobalPosition &globalPos = element.geometry().center();
- if (isInLens_(globalPos))
- return lensK_;
- return outerK_;
- }
-
- /*!
- * \brief Porosity
- *
- * \param element The current element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume.
- * \return Porosity
- */
- Scalar porosity(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- { return 0.4; }
-
- /*!
- * \brief Function for defining the parameters needed by constitutive relationships (kr-Sw, pc-Sw, etc.).
- *
- * \param element The current element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume.
- * \return the material parameters object
- */
- const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
- {
- if (scvIdx > 0)
- return materialLawParams(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
-
- const GlobalPosition &globalPos = element.geometry().center();
-
- if (isInLens_(globalPos))
- return lensMaterialParams_;
- return outerMaterialParams_;
- }
-
-
-private:
- bool isInLens_(const GlobalPosition &pos) const
- {
- for (int i = 0; i < dim; ++i) {
- if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
- return false;
- }
- return true;
- }
-
- GlobalPosition lensLowerLeft_;
- GlobalPosition lensUpperRight_;
-
- Scalar lensK_;
- Scalar outerK_;
- MaterialLawParams lensMaterialParams_;
- MaterialLawParams outerMaterialParams_;
-};
-
-} // end namespace
-#endif
-
diff --git a/test/implicit/2p/lensproblem.hh b/test/implicit/2p/lensproblem.hh
index b4762561527d39c35ce5982733bc90f556ecb9f5..cd2507f91d82c0fe01449788c8ad52d3dbdfe167 100644
--- a/test/implicit/2p/lensproblem.hh
+++ b/test/implicit/2p/lensproblem.hh
@@ -37,6 +37,8 @@
#include <dumux/material/components/dnapl.hh>
#include <dumux/implicit/2p/2pmodel.hh>
#include <dumux/implicit/box/porousmediaboxproblem.hh>
+#include <dumux/implicit/cellcentered/porousmediaccproblem.hh>
+#include <dumux/implicit/cellcentered/ccpropertydefaults.hh>
#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
@@ -53,7 +55,9 @@ class LensProblem;
//////////
namespace Properties
{
-NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP, LensSpatialParams));
+NEW_TYPE_TAG(LensProblem, INHERITS_FROM(TwoP, LensSpatialParams));
+NEW_TYPE_TAG(LensBoxProblem, INHERITS_FROM(BoxModel, LensProblem));
+NEW_TYPE_TAG(LensCCProblem, INHERITS_FROM(CCModel, LensProblem));
// Set the grid type
#if HAVE_UG
@@ -109,6 +113,11 @@ SET_SCALAR_PROP(LensProblem, LinearSolverPreconditionerRelaxation, 1.0);
// Enable gravity
SET_BOOL_PROP(LensProblem, ProblemEnableGravity, true);
+
+NEW_PROP_TAG(BaseProblem);
+SET_TYPE_PROP(LensBoxProblem, BaseProblem, PorousMediaBoxProblem<TypeTag>);
+SET_TYPE_PROP(LensCCProblem, BaseProblem, PorousMediaCCProblem<TypeTag>);
+
}
/*!
@@ -147,9 +156,9 @@ SET_BOOL_PROP(LensProblem, ProblemEnableGravity, true);
* <tt>./test_2p -parameterFile test_2p.input</tt>
*/
template <class TypeTag >
-class LensProblem : public PorousMediaBoxProblem<TypeTag>
+class LensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
{
- typedef PorousMediaBoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
@@ -193,7 +202,8 @@ public:
*/
LensProblem(TimeManager &timeManager,
const GridView &gridView)
- : ParentType(timeManager, gridView)
+ : ParentType(timeManager, gridView),
+ isBox_(GET_PROP_VALUE(TypeTag, ImplicitIsBox))
{
eps_ = 3e-6;
temperature_ = 273.15 + 20; // -> 20°C
@@ -210,7 +220,12 @@ public:
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
- { return "lens"; }
+ {
+ if (isBox_)
+ return "lensbox";
+ else
+ return "lenscc";
+ }
/*!
* \brief Called directly after the time integration.
@@ -283,29 +298,43 @@ public:
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
-
+
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
-
- if (onLeftBoundary_(globalPos))
+
+ if (isBox_)
{
- Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
- Scalar depth = this->bboxMax()[1] - globalPos[1];
- Scalar alpha = (1 + 1.5/height);
-
- // hydrostatic pressure scaled by alpha
- values[pwIdx] = 1e5 - alpha*densityW*this->gravity()[1]*depth;
- values[SnIdx] = 0.0;
+ if (onLeftBoundary_(globalPos))
+ {
+ Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
+ Scalar depth = this->bboxMax()[1] - globalPos[1];
+ Scalar alpha = (1 + 1.5/height);
+
+ // hydrostatic pressure scaled by alpha
+ values[pwIdx] = 1e5 - alpha*densityW*this->gravity()[1]*depth;
+ values[SnIdx] = 0.0;
+ }
+ else if (onRightBoundary_(globalPos))
+ {
+ Scalar depth = this->bboxMax()[1] - globalPos[1];
+
+ // hydrostatic pressure
+ values[pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+ values[SnIdx] = 0.0;
+ }
+ else
+ values = 0.0;
}
- else if (onRightBoundary_(globalPos))
+ else
{
- Scalar depth = this->bboxMax()[1] - globalPos[1];
-
- // hydrostatic pressure
- values[pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+ Scalar height = 4.0;//this->bboxMax()[1] - this->bboxMin()[1];
+ Scalar depth = 4.0 - globalPos[1];//this->bboxMax()[1] - globalPos[1];
+ Scalar alpha = 1 + 1.5/height;
+ Scalar factor = (6.0*alpha + (1.0 - alpha)*globalPos[0])/6.0;
+
+ // hydrostatic pressure scaled by alpha
+ values[pwIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
values[SnIdx] = 0.0;
}
- else
- values = 0.0;
}
/*!
@@ -346,15 +375,19 @@ public:
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
- Scalar depth = this->bboxMax()[1] - globalPos[1];
-
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
-
+
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
+ Scalar depth;
+ if (isBox_)
+ depth = this->bboxMax()[1] - globalPos[1];
+ else
+ depth = 4.0 - globalPos[1];
+
// hydrostatic pressure
values[pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
values[SnIdx] = 0.0;
@@ -385,13 +418,23 @@ private:
bool onInlet_(const GlobalPosition &globalPos) const
{
- Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
- Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
- return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+ if (isBox_)
+ {
+ Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+ }
+ else
+ {
+ Scalar width = 6.0;//this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (6.0/*this->bboxMax()[0]*/ - globalPos[0])/width;
+ return (globalPos[1] > 4.0 - eps_) && 0.5 < lambda && lambda < 2.0/3.0;
+ }
}
Scalar temperature_;
Scalar eps_;
+ const bool isBox_;
};
} //end namespace
diff --git a/test/implicit/2p/lensspatialparams.hh b/test/implicit/2p/lensspatialparams.hh
index 068695f2be2b39b01949abd896815cbf52cf0a6d..50d24f97836a13c45c894b86d82232b1781996f5 100644
--- a/test/implicit/2p/lensspatialparams.hh
+++ b/test/implicit/2p/lensspatialparams.hh
@@ -91,7 +91,8 @@ public:
typedef typename MaterialLaw::Params MaterialLawParams;
LensSpatialParams(const GridView& gridView)
- : ParentType(gridView)
+ : ParentType(gridView),
+ isBox_(GET_PROP_VALUE(TypeTag, ImplicitIsBox))
{
try
{
@@ -145,7 +146,11 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (!isBox_ && scvIdx > 0)
+ return intrinsicPermeability(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
+
+ const GlobalPosition& globalPos = fvElemGeom.subContVol[scvIdx].global;
+
if (isInLens_(globalPos))
return lensK_;
return outerK_;
@@ -176,8 +181,11 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
-
+ if (!isBox_ && scvIdx > 0)
+ return materialLawParams(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
+
+ const GlobalPosition& globalPos = fvElemGeom.subContVol[scvIdx].global;
+
if (isInLens_(globalPos))
return lensMaterialParams_;
return outerMaterialParams_;
@@ -201,6 +209,7 @@ private:
Scalar outerK_;
MaterialLawParams lensMaterialParams_;
MaterialLawParams outerMaterialParams_;
+ const bool isBox_;
};
} // end namespace
diff --git a/test/implicit/2p/test_2p.cc b/test/implicit/2p/test_box2p.cc
similarity index 98%
rename from test/implicit/2p/test_2p.cc
rename to test/implicit/2p/test_box2p.cc
index a4dc66c7bf6ac09c33dd0e80fd6edade21519590..afec6a3114977f23ae9458950581d97fc130889f 100644
--- a/test/implicit/2p/test_2p.cc
+++ b/test/implicit/2p/test_box2p.cc
@@ -61,6 +61,6 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
int main(int argc, char** argv)
{
- typedef TTAG(LensProblem) TypeTag;
+ typedef TTAG(LensBoxProblem) TypeTag;
return Dumux::start<TypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/2p/test_2p.input b/test/implicit/2p/test_box2p.input
similarity index 100%
rename from test/implicit/2p/test_2p.input
rename to test/implicit/2p/test_box2p.input
diff --git a/test/implicit/2p/test_cc2p.cc b/test/implicit/2p/test_cc2p.cc
index c3b3a6d8393c2853ae37594a3614511c77af0f1a..747e9d350a9f0b8926015a89a5da8a5354113b5d 100644
--- a/test/implicit/2p/test_cc2p.cc
+++ b/test/implicit/2p/test_cc2p.cc
@@ -27,7 +27,7 @@
* \brief test for the two-phase box model
*/
#include "config.h"
-#include "lensccproblem.hh"
+#include "lensproblem.hh"
#include <dumux/common/start.hh>
/*!
@@ -66,6 +66,6 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
int main(int argc, char** argv)
{
- typedef TTAG(LensProblem) TypeTag;
+ typedef TTAG(LensCCProblem) TypeTag;
return Dumux::start<TypeTag>(argc, argv, usage);
}