diff --git a/test/implicit/2p/lensccproblem.hh b/test/implicit/2p/lensccproblem.hh
deleted file mode 100644
index fb2b8e78f3dc2b0e7d38d9e21f061f1b8ab75912..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 "lensspatialparams.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/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/test_2p.cc b/test/implicit/2p/test_2p.cc
index a4dc66c7bf6ac09c33dd0e80fd6edade21519590..afec6a3114977f23ae9458950581d97fc130889f 100644
--- a/test/implicit/2p/test_2p.cc
+++ b/test/implicit/2p/test_2p.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_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);
}