From 758166e135b8bf47040b840c85bb7fe056bb18ce Mon Sep 17 00:00:00 2001 From: Markus Wolff <markus.wolff@twt-gmbh.de> Date: Mon, 8 Aug 2011 08:05:28 +0000 Subject: [PATCH] Adapted test problems to changes in the models git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6406 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- test/boxmodels/1p/1ptestproblem.hh | 6 + test/boxmodels/1p2c/tissue_tumor_problem.hh | 4 +- test/boxmodels/2p/lensproblem.hh | 33 +- test/boxmodels/2p/lensspatialparameters.hh | 33 +- test/boxmodels/2p2c/injectionproblem.hh | 14 +- .../2p2c/injectionspatialparameters.hh | 30 +- test/boxmodels/2p2cni/waterairproblem.hh | 14 +- .../2p2cni/waterairspatialparameters.hh | 29 +- test/boxmodels/2pni/injectionproblem2pni.hh | 16 +- .../boxmodels/richards/richardslensproblem.hh | 24 +- .../richards/richardslensspatialparameters.hh | 38 +- test/common/Makefile.am | 2 +- test/common/generalproblem/CMakeLists.txt | 18 + test/common/generalproblem/Makefile.am | 8 + .../generalproblem/generallensproblem.hh | 471 ++++++++++++++++++ .../generallensspatialparameters.hh | 191 +++++++ .../generalproblem/test_generalproblem_2p.cc | 166 ++++++ test/decoupled/1p/benchmarkresult.hh | 2 +- test/decoupled/1p/test_1p_problem.hh | 11 +- test/decoupled/1p/test_diffusion_problem.hh | 47 +- .../1p/test_diffusion_spatialparams.hh | 45 +- test/decoupled/2p/test_impes_problem.hh | 71 +-- test/decoupled/2p/test_impes_spatialparams.hh | 65 +-- test/decoupled/2p/test_transport_problem.hh | 17 +- .../2p/test_transport_spatialparams.hh | 59 ++- 25 files changed, 1168 insertions(+), 246 deletions(-) create mode 100644 test/common/generalproblem/CMakeLists.txt create mode 100644 test/common/generalproblem/Makefile.am create mode 100644 test/common/generalproblem/generallensproblem.hh create mode 100644 test/common/generalproblem/generallensspatialparameters.hh create mode 100644 test/common/generalproblem/test_generalproblem_2p.cc diff --git a/test/boxmodels/1p/1ptestproblem.hh b/test/boxmodels/1p/1ptestproblem.hh index 93002cc041..24ca13e739 100644 --- a/test/boxmodels/1p/1ptestproblem.hh +++ b/test/boxmodels/1p/1ptestproblem.hh @@ -161,6 +161,12 @@ public: { return 273.15 + 10; } // 10°C // \} + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + /*! * \name Boundary conditions */ diff --git a/test/boxmodels/1p2c/tissue_tumor_problem.hh b/test/boxmodels/1p2c/tissue_tumor_problem.hh index 692b6b3fac..79f3abc63d 100644 --- a/test/boxmodels/1p2c/tissue_tumor_problem.hh +++ b/test/boxmodels/1p2c/tissue_tumor_problem.hh @@ -229,7 +229,6 @@ public: * in normal direction of each component. Negative values mean * influx. */ - using ParentType::neumann; void neumann(PrimaryVariables &values, const Element &element, const FVElementGeometry &fvElemGeom, @@ -268,8 +267,7 @@ public: * unit. Positive values mean that mass is created, negative ones * mean that it vanishes. */ - using ParentType::source; - void source(PrimaryVariables &values, + void sourceAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { values = Scalar(0.0); diff --git a/test/boxmodels/2p/lensproblem.hh b/test/boxmodels/2p/lensproblem.hh index 0f890c1e85..5fbfa04428 100644 --- a/test/boxmodels/2p/lensproblem.hh +++ b/test/boxmodels/2p/lensproblem.hh @@ -55,7 +55,7 @@ class LensProblem; ////////// namespace Properties { -NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP)); +NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP, LensSpatialParameters)); // Set the grid type SET_PROP(LensProblem, Grid) @@ -91,12 +91,6 @@ public: typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleDNAPL<Scalar> > type; }; -// Set the spatial parameters -SET_PROP(LensProblem, SpatialParameters) -{ - typedef Dumux::LensSpatialParameters<TypeTag> type; -}; - // Enable the time step ramp up inside the newton method? SET_BOOL_PROP(LensProblem, EnableTimeStepRampUp, false); @@ -261,6 +255,13 @@ public: Scalar temperature() const { return temperature_; }; + + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + // \} /*! @@ -273,14 +274,11 @@ public: * used for which equation on a given boundary control volume. * * \param values The boundary types for the conservation equations - * \param vertex The vertex on the boundary for which the - * conditions needs to be specified + * \param globalPos The position of the center of the finite volume */ - void boundaryTypes(BoundaryTypes &values, - const Vertex &vertex) const + void boundaryTypesAtPos(BoundaryTypes &values, + const GlobalPosition &globalPos) const { - const GlobalPosition globalPos = vertex.geometry().center(); - if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) { values.setAllDirichlet(); } @@ -298,8 +296,7 @@ public: * * For this method, the \a values parameter stores primary variables. */ - using ParentType::dirichlet; - void dirichlet(PrimaryVariables &values, + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { Scalar densityW = FluidSystem::componentDensity(wPhaseIdx, @@ -339,8 +336,7 @@ public: * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. */ - using ParentType::neumann; - void neumann(PrimaryVariables &values, + void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { values = 0.0; @@ -365,8 +361,7 @@ public: * For this method, the \a values parameter stores primary * variables. */ - using ParentType::initial; - void initial(PrimaryVariables &values, + void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { Scalar depth = this->bboxMax()[1] - globalPos[1]; diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh index 21200a1e50..3188c22e20 100644 --- a/test/boxmodels/2p/lensspatialparameters.hh +++ b/test/boxmodels/2p/lensspatialparameters.hh @@ -39,6 +39,31 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class LensSpatialParameters; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(LensSpatialParameters); + +// Set the spatial parameters +SET_TYPE_PROP(LensSpatialParameters, SpatialParameters, Dumux::LensSpatialParameters<TypeTag>); + +// Set the material Law +SET_PROP(LensSpatialParameters, MaterialLaw) +{ +private: + // define the material law which is parameterized by effective + // saturations + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; +public: + // define the material law parameterized by absolute saturations + typedef EffToAbsLaw<EffectiveLaw> type; +}; +} /*! * \ingroup TwoPBoxModel * @@ -69,13 +94,9 @@ class LensSpatialParameters : public BoxSpatialParameters<TypeTag> typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry; - // define the material law which is parameterized by effective - // saturations - typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; - public: - // define the material law parameterized by absolute saturations - typedef EffToAbsLaw<EffectiveLaw> MaterialLaw; + //get the material law from the property system + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; LensSpatialParameters(const GridView& gridView) diff --git a/test/boxmodels/2p2c/injectionproblem.hh b/test/boxmodels/2p2c/injectionproblem.hh index f775b7bde8..7c059317ac 100644 --- a/test/boxmodels/2p2c/injectionproblem.hh +++ b/test/boxmodels/2p2c/injectionproblem.hh @@ -50,7 +50,7 @@ class InjectionProblem; namespace Properties { -NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(BoxTwoPTwoC)); +NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(BoxTwoPTwoC, InjectionSpatialParameters)); // Set the grid type SET_PROP(InjectionProblem, Grid) @@ -72,11 +72,6 @@ SET_PROP(InjectionProblem, typedef Dumux::H2O_N2_System<TypeTag> type; }; -// Set the spatial parameters -SET_TYPE_PROP(InjectionProblem, - SpatialParameters, - Dumux::InjectionSpatialParameters<TypeTag>); - // Enable gravity SET_BOOL_PROP(InjectionProblem, EnableGravity, true); @@ -203,6 +198,12 @@ public: Scalar temperature() const { return temperature_; }; + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + // \} /*! @@ -257,7 +258,6 @@ public: * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. */ - using ParentType::neumann; void neumann(PrimaryVariables &values, const Element &element, const FVElementGeometry &fvElemGeom, diff --git a/test/boxmodels/2p2c/injectionspatialparameters.hh b/test/boxmodels/2p2c/injectionspatialparameters.hh index eb34e147c5..820ca1a329 100644 --- a/test/boxmodels/2p2c/injectionspatialparameters.hh +++ b/test/boxmodels/2p2c/injectionspatialparameters.hh @@ -36,6 +36,32 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class InjectionSpatialParameters; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(InjectionSpatialParameters); + +// Set the spatial parameters +SET_TYPE_PROP(InjectionSpatialParameters, SpatialParameters, Dumux::InjectionSpatialParameters<TypeTag>); + +// Set the material Law +SET_PROP(InjectionSpatialParameters, MaterialLaw) +{ +private: + // define the material law which is parameterized by effective + // saturations + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw; +public: + // define the material law parameterized by absolute saturations + typedef EffToAbsLaw<EffMaterialLaw> type; +}; +} + /*! * \ingroup TwoPTwoCModel * @@ -75,10 +101,8 @@ class InjectionSpatialParameters : public BoxSpatialParameters<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry; typedef typename GridView::template Codim<0>::Entity Element; - typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw; - //typedef LinearMaterial<Scalar> EffMaterialLaw; public: - typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; /*! diff --git a/test/boxmodels/2p2cni/waterairproblem.hh b/test/boxmodels/2p2cni/waterairproblem.hh index 6a007a5a4c..3ed5227226 100644 --- a/test/boxmodels/2p2cni/waterairproblem.hh +++ b/test/boxmodels/2p2cni/waterairproblem.hh @@ -46,7 +46,7 @@ class WaterAirProblem; namespace Properties { -NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI)); +NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI, WaterAirSpatialParameters)); // Set the grid type SET_PROP(WaterAirProblem, Grid) @@ -63,11 +63,6 @@ SET_PROP(WaterAirProblem, Problem) // Set the wetting phase SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::H2O_N2_System<TypeTag>); -// Set the spatial parameters -SET_TYPE_PROP(WaterAirProblem, - SpatialParameters, - Dumux::WaterAirSpatialParameters<TypeTag>); - // Enable gravity SET_BOOL_PROP(WaterAirProblem, EnableGravity, true); @@ -204,6 +199,12 @@ public: }; #endif + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + // \} /*! @@ -262,7 +263,6 @@ public: * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. */ - using ParentType::neumann; void neumann(PrimaryVariables &values, const Element &element, const FVElementGeometry &fvElemGeom, diff --git a/test/boxmodels/2p2cni/waterairspatialparameters.hh b/test/boxmodels/2p2cni/waterairspatialparameters.hh index 326bb0838d..0c0a2cd19f 100644 --- a/test/boxmodels/2p2cni/waterairspatialparameters.hh +++ b/test/boxmodels/2p2cni/waterairspatialparameters.hh @@ -34,6 +34,32 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class WaterAirSpatialParameters; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(WaterAirSpatialParameters); + +// Set the spatial parameters +SET_TYPE_PROP(WaterAirSpatialParameters, SpatialParameters, Dumux::WaterAirSpatialParameters<TypeTag>); + +// Set the material Law +SET_PROP(WaterAirSpatialParameters, MaterialLaw) +{ +private: + // define the material law which is parameterized by effective + // saturations + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw; +public: + // define the material law parameterized by absolute saturations + typedef EffToAbsLaw<EffMaterialLaw> type; +}; +} + /*! * \ingroup TwoPTwoCNIModel * @@ -75,9 +101,8 @@ class WaterAirSpatialParameters : public BoxSpatialParameters<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem; - typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw; public: - typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; /*! diff --git a/test/boxmodels/2pni/injectionproblem2pni.hh b/test/boxmodels/2pni/injectionproblem2pni.hh index a418d3996c..6f8b92f894 100644 --- a/test/boxmodels/2pni/injectionproblem2pni.hh +++ b/test/boxmodels/2pni/injectionproblem2pni.hh @@ -51,9 +51,9 @@ class InjectionProblem2PNI; namespace Properties { #if !ISOTHERMAL -NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoPNI)); +NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoPNI, InjectionSpatialParameters)); #else -NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoP)); +NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoP, InjectionSpatialParameters)); #endif // Set the grid type @@ -73,11 +73,6 @@ SET_PROP(InjectionProblem2PNI, Problem) typedef Dumux::InjectionProblem2PNI<TypeTag> type; }; -// Set the spatial parameters. we use the same spatial parameters as the -// 2p2c injection problem -SET_PROP(InjectionProblem2PNI, SpatialParameters) -{ typedef InjectionSpatialParameters<TypeTag> type; }; - #if 1 // Use the same fluid system as the 2p2c injection problem SET_PROP(InjectionProblem2PNI, FluidSystem) @@ -237,6 +232,12 @@ public: }; #endif + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + // \} /*! @@ -302,7 +303,6 @@ public: * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. */ - using ParentType::neumann; void neumann(PrimaryVariables &values, const Element &element, const FVElementGeometry &fvElemGeom, diff --git a/test/boxmodels/richards/richardslensproblem.hh b/test/boxmodels/richards/richardslensproblem.hh index 55418a0bfd..e59a32e2e5 100644 --- a/test/boxmodels/richards/richardslensproblem.hh +++ b/test/boxmodels/richards/richardslensproblem.hh @@ -49,7 +49,7 @@ class RichardsLensProblem; ////////// namespace Properties { -NEW_TYPE_TAG(RichardsLensProblem, INHERITS_FROM(BoxRichards)); +NEW_TYPE_TAG(RichardsLensProblem, INHERITS_FROM(BoxRichards, RichardsLensSpatialParameters)); // Set the grid type. Use UG if available, else SGrid #if HAVE_UG @@ -72,10 +72,6 @@ public: typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type; }; -// Set the spatial parameters -SET_PROP(RichardsLensProblem, SpatialParameters) -{ typedef Dumux::RichardsLensSpatialParameters<TypeTag> type; }; - // Enable gravity? SET_BOOL_PROP(RichardsLensProblem, EnableGravity, true); @@ -208,6 +204,12 @@ public: int scvIdx) const { return pnRef_; }; + void sourceAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + values = 0; + } + // \} /*! @@ -222,8 +224,7 @@ public: * \param values The boundary types for the conservation equations * \param vertex The vertex for which the boundary type is set */ - using ParentType::boundaryTypes; - void boundaryTypes(BoundaryTypes &values, + void boundaryTypesAtPos(BoundaryTypes &values, const GlobalPosition &globalPos) const { if (onLeftBoundary_(globalPos) || @@ -244,8 +245,7 @@ public: * * For this method, the \a values parameter stores primary variables. */ - using ParentType::dirichlet; - void dirichlet(PrimaryVariables &values, + void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { // use initial values as boundary conditions @@ -269,8 +269,7 @@ public: * \param boundaryFaceIdx The index of the boundary face of the * finite volume geometry */ - using ParentType::neumann; - void neumann(PrimaryVariables &values, + void neumannAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const { values = 0.0; @@ -298,8 +297,7 @@ public: * \param scvIdx The sub control volume index of the finite * volume geometry */ - using ParentType::initial; - void initial(PrimaryVariables &values, + void initialAtPos(PrimaryVariables &values, const GlobalPosition &pos) const { initial_(values, pos); }; diff --git a/test/boxmodels/richards/richardslensspatialparameters.hh b/test/boxmodels/richards/richardslensspatialparameters.hh index f145bd4547..7e0ee64885 100644 --- a/test/boxmodels/richards/richardslensspatialparameters.hh +++ b/test/boxmodels/richards/richardslensspatialparameters.hh @@ -36,6 +36,32 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class RichardsLensSpatialParameters; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(RichardsLensSpatialParameters); + +// Set the spatial parameters +SET_TYPE_PROP(RichardsLensSpatialParameters, SpatialParameters, Dumux::RichardsLensSpatialParameters<TypeTag>); + +// Set the material Law +SET_PROP(RichardsLensSpatialParameters, MaterialLaw) +{ +private: + // define the material law which is parameterized by effective + // saturations + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; +public: + // define the material law parameterized by absolute saturations + typedef EffToAbsLaw<EffectiveLaw> type; +}; +} + /*! * \ingroup RichardsModel * \brief The spatial parameters for the RichardsLensProblem @@ -59,17 +85,9 @@ class RichardsLensSpatialParameters : public BoxSpatialParameters<TypeTag> typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry; - // define the material law which is parameterized by effective - // saturations - typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; - public: - /*! - * \brief The material law to be used. - * - * This problem uses the RegularizedVanGenuchten material law. - */ - typedef EffToAbsLaw<EffectiveLaw> MaterialLaw; + //get the material law from the property system + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; //! The parameters of the material law to be used typedef typename MaterialLaw::Params MaterialLawParams; diff --git a/test/common/Makefile.am b/test/common/Makefile.am index ff73a4e69f..2f543318c6 100644 --- a/test/common/Makefile.am +++ b/test/common/Makefile.am @@ -1,4 +1,4 @@ -SUBDIRS = pardiso propertysystem spline . +SUBDIRS = generalproblem pardiso propertysystem spline . EXTRA_DIST = CMakeLists.txt diff --git a/test/common/generalproblem/CMakeLists.txt b/test/common/generalproblem/CMakeLists.txt new file mode 100644 index 0000000000..16728c1f53 --- /dev/null +++ b/test/common/generalproblem/CMakeLists.txt @@ -0,0 +1,18 @@ +add_definitions(-DYASPGRID -DGRIDDIM=2 -DENABLE_UG) + +# build target for the simple twophase lens problem in a general definition for both box and decoupled model +ADD_EXECUTABLE("test_generalproblem_2p" test_generalproblem_2p.cc) +TARGET_LINK_LIBRARIES("test_2p" ${DumuxLinkLibraries}) + +# add required libraries and includes to the build flags +LINK_DIRECTORIES(${DumuxLinkDirectories}) +INCLUDE_DIRECTORIES(${DumuxIncludeDirectories}) + +# make sure the grids are present in the build directory +add_custom_command(TARGET "test_generalproblem_2p" + POST_BUILD + COMMAND ${CMAKE_COMMAND} -E + copy_directory + "${CMAKE_CURRENT_SOURCE_DIR}/grids" + "${CMAKE_CURRENT_BINARY_DIR}/grids") + diff --git a/test/common/generalproblem/Makefile.am b/test/common/generalproblem/Makefile.am new file mode 100644 index 0000000000..9f5220ae4f --- /dev/null +++ b/test/common/generalproblem/Makefile.am @@ -0,0 +1,8 @@ +check_PROGRAMS = test_generalproblem_2p + +noinst_HEADERS = *.hh +EXTRA_DIST = CMakeLists.txt + +test_generalproblem_2p_SOURCES = test_generalproblem_2p.cc + +include $(top_srcdir)/am/global-rules diff --git a/test/common/generalproblem/generallensproblem.hh b/test/common/generalproblem/generallensproblem.hh new file mode 100644 index 0000000000..a3eafc5577 --- /dev/null +++ b/test/common/generalproblem/generallensproblem.hh @@ -0,0 +1,471 @@ +/***************************************************************************** + * Copyright (C) 2011 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 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. * + * * + * 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_GENERALLENSPROBLEM_HH +#define DUMUX_GENERALLENSPROBLEM_HH + +//common includes +#if HAVE_UG +#include <dune/grid/uggrid.hh> +#endif +#include <dune/grid/yaspgrid.hh> + +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/components/simplednapl.hh> +#include <dumux/material/fluidsystems/liquidphase.hh> + +//box model +#include <dumux/boxmodels/2p/2pmodel.hh> + +//decoupled model +#include <dumux/decoupled/2p/impes/impesproblem2p.hh> +#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh> +#include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh> +#include<dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh> + +#include "generallensspatialparameters.hh" + +namespace Dumux +{ + +template <class TypeTag> +class GeneralLensProblem; + +////////// +// Specify the properties for the lens problem +////////// +namespace Properties +{ +//Set the general problem TypeTag which does not depend on the model +NEW_TYPE_TAG(GeneralLensProblem, INHERITS_FROM(GeneralLensSpatialParameters)); + +// Property for defining the model specific problem base class +NEW_PROP_TAG(ProblemBaseClass); + +// Set the grid type +SET_PROP(GeneralLensProblem, Grid) +{ +// typedef Dune::UGGrid<2> type; + typedef Dune::YaspGrid<2> type; +}; + +// Set the problem property +SET_PROP(GeneralLensProblem, Problem) +{ + typedef Dumux::GeneralLensProblem<TypeTag> type; +}; + +// Set the wetting phase +SET_PROP(GeneralLensProblem, WettingPhase) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; +public: + typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type; +}; + +// Set the non-wetting phase +SET_PROP(GeneralLensProblem, NonwettingPhase) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; +public: + typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleDNAPL<Scalar> > type; +}; + +// Enable gravity +SET_BOOL_PROP(GeneralLensProblem, EnableGravity, true); + +/////////////////////////////////////////////////// +// Box model TypeTag +////////////////////////////////////////////////// + +NEW_TYPE_TAG(BoxGeneralLensProblem, INHERITS_FROM(GeneralLensProblem, BoxTwoP)); + + +// Set the problem property +SET_PROP(BoxGeneralLensProblem, ProblemBaseClass) +{ + typedef Dumux::TwoPProblem<TypeTag> type; +}; + +// Set the problem property +SET_PROP(BoxGeneralLensProblem, SpatialParamsBaseClass) +{ + typedef Dumux::BoxSpatialParameters<TypeTag> type; +}; + + +/////////////////////////////////////////////////// +// Deoupled model TypeTag +////////////////////////////////////////////////// + +NEW_TYPE_TAG(DecoupledGeneralLensProblem, INHERITS_FROM(DecoupledTwoP, Transport, GeneralLensProblem)); + + +// Set the problem property +SET_PROP(DecoupledGeneralLensProblem, ProblemBaseClass) +{ + typedef Dumux::IMPESProblem2P<TypeTag> type; +}; + +// Set the problem property +SET_PROP(DecoupledGeneralLensProblem, SpatialParamsBaseClass) +{ + typedef Dumux::FVSpatialParameters<TypeTag> type; +}; + +// Set the model properties +SET_TYPE_PROP(DecoupledGeneralLensProblem, TransportModel, Dumux::FVSaturation2P<TTAG(DecoupledGeneralLensProblem)>); + +SET_PROP(DecoupledGeneralLensProblem, PressureModel) +{ + typedef Dumux::FVVelocity2P<TTAG(DecoupledGeneralLensProblem)> type; +}; + +SET_INT_PROP(DecoupledGeneralLensProblem, Formulation, + DecoupledTwoPCommonIndices::pwSn); + +SET_TYPE_PROP(DecoupledGeneralLensProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>); + +SET_SCALAR_PROP(DecoupledGeneralLensProblem, CFLFactor, 1.0); +} + +/*! + * \ingroup TwoPBoxModel + * \ingroup IMPETtests + * + * \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 TwoPBoxModel. + * + * 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 20000 250</tt> + */ +template <class TypeTag > +class GeneralLensProblem : public GET_PROP_TYPE(TypeTag, PTAG(ProblemBaseClass)) +{ + typedef GeneralLensProblem<TypeTag> ThisType; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(ProblemBaseClass)) ParentType; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem; + typedef TwoPFluidState<TypeTag> FluidState; + + enum { + numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)), + + // primary variable indices + pressureIdx = Indices::pressureIdx, + saturationIdx = Indices::saturationIdx, + pwIdx = Indices::pwIdx, + SnIdx = Indices::SnIdx, + + // equation indices + contiWEqIdx = Indices::contiWEqIdx, + 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, 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<dim>::Entity Vertex; + typedef typename GridView::Intersection Intersection; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef Dune::FieldVector<Scalar, dim> LocalPosition; + typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition; + +public: + /*! + * \brief The constructor + * + * \param timeManager The time manager + * \param gridView The grid view + * \param lensLowerLeft Global position of the lenses lower left corner + * \param lensUpperRight Global position of the lenses upper right corner + */ + GeneralLensProblem(TimeManager &timeManager, + const GridView &gridView, + const GlobalPosition &lensLowerLeft, + const GlobalPosition &lensUpperRight) + : ParentType(timeManager, gridView) + { + temperature_ = 273.15 + 20; // -> 20°C + this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight); + this->timeManager().startNextEpisode(500); + } + + /*! + * \name Problem parameters + */ + // \{ + + + bool shouldWriteRestartFile() + { + return false; + } + + bool shouldWriteOutput() const + { + if (this->timeManager().time() < eps_ || this->timeManager().willBeFinished() || this->timeManager().episodeWillBeOver()) + { + return true; + } + return false; + } + + void episodeEnd() + { + if (!this->timeManager().willBeFinished()) + this->timeManager().startNextEpisode(500); + }; + + /*! + * \brief Returns the temperature within the domain. + * + * \param globalPos The global coordinates + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperatureAtPos(const GlobalPosition &globalPos) const + { + return temperature_; + }; + + //! Returns the reference pressure for evaluation of constitutive relations + /*! + * \param globalPos The global coordinates + * + * Only for decoupled model -> incompressible. + */ + Scalar referencePressureAtPos(const GlobalPosition& globalPos) const + { + return 1e5; // -> 10°C + } + + //! Returns the source term + /*! + * \param globalPos The global coordinates + */ + 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 global coordinates of the boundary + */ + 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 + * boundary. + * + * \param values The dirichlet values for the primary variables + * \param globalPos The global coordinates of the boundary + * + * For this method, the \a values parameter stores primary variables. + */ + void dirichletAtPos(PrimaryVariables &values, + const GlobalPosition &globalPos) const + { + Scalar densityW = FluidSystem::componentDensity(wPhaseIdx, + wPhaseIdx, + temperature_, + 1e5); + + 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; + } + + /*! + * \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 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[nPhaseIdx] = -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 = this->bboxMax()[1] - globalPos[1]; + Scalar densityW = FluidSystem::componentDensity(wPhaseIdx, + wPhaseIdx, + temperature_, + 1e5); + + // 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 = 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; + } + + Scalar temperature_; + static constexpr Scalar eps_ = 3e-6; +}; +} //end namespace + +#endif diff --git a/test/common/generalproblem/generallensspatialparameters.hh b/test/common/generalproblem/generallensspatialparameters.hh new file mode 100644 index 0000000000..0f23fa4c6d --- /dev/null +++ b/test/common/generalproblem/generallensspatialparameters.hh @@ -0,0 +1,191 @@ +/***************************************************************************** + * Copyright (C) 2011 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 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. * + * * + * 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 GeneralLensProblem which uses the + * twophase box model or twophase decoupled model + */ +#ifndef DUMUX_GENERALLENSSPATIALPARAMETERS_HH +#define DUMUX_GENERALLENSSPATIALPARAMETERS_HH + +#include <dumux/material/spatialparameters/boxspatialparameters.hh> +#include <dumux/material/spatialparameters/fvspatialparameters.hh> +#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh> +#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> + +namespace Dumux +{ +//forward declaration +template<class TypeTag> +class GeneralLensSpatialParameters; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(GeneralLensSpatialParameters); + +// Property to define the spatial parameters base class -> allows switch with model switch! +NEW_PROP_TAG(SpatialParamsBaseClass); + +// Set the spatial parameters +SET_TYPE_PROP(GeneralLensSpatialParameters, SpatialParameters, Dumux::GeneralLensSpatialParameters<TypeTag>); + +// Set the material Law +SET_PROP(GeneralLensSpatialParameters, MaterialLaw) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedVanGenuchten<Scalar> EffectiveLaw; +public: + typedef EffToAbsLaw<EffectiveLaw> type; +}; +} + +/*! + * \ingroup TwoPBoxModel + * \ingroup IMPETtests + * + * \brief The spatial parameters for the LensProblem which uses the + * twophase box model or twophase decoupled model + */ +template<class TypeTag> +class GeneralLensSpatialParameters : public GET_PROP_TYPE(TypeTag, PTAG(SpatialParamsBaseClass)) +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParamsBaseClass)) ParentType; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef typename Grid::ctype CoordScalar; + + typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices; + + enum { + dim=GridView::dimension, + dimWorld=GridView::dimensionworld, + + wPhaseIdx = Indices::wPhaseIdx, + nPhaseIdx = Indices::nPhaseIdx + }; + + typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition; + typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; + + typedef typename GridView::template Codim<0>::Entity Element; + + +public: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + + GeneralLensSpatialParameters(const GridView& gridView) + : ParentType(gridView) + { + // 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); + + for (int i=0; i < dim; i++) + { + lensK_[i][i] = 9.05e-12; + outerK_[i][i] = 4.6e-10; + } + } + + /*! + * \brief Get the intrinsic permeability tensor + * + * \param globalPos The global coordinates of the finite volume + */ + const FieldMatrix& intrinsicPermeabilityAtPos( + const GlobalPosition &globalPos) const + { + if (isInLens_(globalPos)) + return lensK_; + return outerK_; + } + + /*! + * \brief Get the porosity + * + * \param globalPos The global coordinates of the finite volume + */ + Scalar porosityAtPos(const GlobalPosition &globalPos) const + { return 0.4; } + + + /*! + * \brief Get the material law parameters + * + * \param globalPos The global coordinates of the finite volume + * + * \return the parameter object for the material law which depends on the position + */ + const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition &globalPos) const + { + if (isInLens_(globalPos)) + return lensMaterialParams_; + return outerMaterialParams_; + } + + + //! Set the bounding box of the fine-sand lens + void setLensCoords(const GlobalPosition& lensLowerLeft, + const GlobalPosition& lensUpperRight) + { + lensLowerLeft_ = lensLowerLeft; + lensUpperRight_ = lensUpperRight; + } + +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_; + + FieldMatrix lensK_; + FieldMatrix outerK_; + + MaterialLawParams lensMaterialParams_; + MaterialLawParams outerMaterialParams_; +}; + +} // end namespace +#endif + diff --git a/test/common/generalproblem/test_generalproblem_2p.cc b/test/common/generalproblem/test_generalproblem_2p.cc new file mode 100644 index 0000000000..0538a32324 --- /dev/null +++ b/test/common/generalproblem/test_generalproblem_2p.cc @@ -0,0 +1,166 @@ +/***************************************************************************** + * Copyright (C) 2007-2008 by Klaus Mosthaf * + * Copyright (C) 2007-2008 by Bernd Flemisch * + * 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. * + * * + * 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 test for the two-phase box model + */ +#include "config.h" + +#include "generallensproblem.hh" + +#include <dune/grid/common/gridinfo.hh> +#include <dune/grid/utility/structuredgridfactory.hh> + +#include <dune/common/exceptions.hh> +#include <dune/common/mpihelper.hh> + +#include <iostream> +#include <boost/format.hpp> + + +//////////////////////// +// the main function +//////////////////////// +void usage(const char *progname) +{ + std::cout << boost::format("usage: %s --box/--decoupled tEnd dt [refineLevel]\n")%progname; + exit(1); +} + +int main(int argc, char** argv) +{ +#ifdef NDEBUG + try { +#endif + //TypeTag which chooses the box model + typedef TTAG(BoxGeneralLensProblem) BoxTypeTag; + //TypeTag which chooses the decoupled model + typedef TTAG(DecoupledGeneralLensProblem) DecoupledTypeTag; + typedef GET_PROP_TYPE(BoxTypeTag, PTAG(Scalar)) Scalar; + typedef GET_PROP_TYPE(BoxTypeTag, PTAG(Grid)) Grid; +// typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef GET_PROP_TYPE(BoxTypeTag, PTAG(TimeManager)) BoxTimeManager; + typedef GET_PROP_TYPE(DecoupledTypeTag, PTAG(TimeManager)) DecoupledTimeManager; + typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition; + + static const int dim = Grid::dimension; + + // initialize MPI, finalize is done automatically on exit + Dune::MPIHelper::instance(argc, argv); + + //////////////////////////////////////////////////////////// + // parse the command line arguments + //////////////////////////////////////////////////////////// + if (argc < 4) + usage(argv[0]); + + // deal with the restart stuff + int argPos = 1; + bool useBoxModel = false; + bool useDecoupledModel = false; + if (std::string(argv[argPos]) == "--box") + { + useBoxModel = true; + } + else if (std::string(argv[argPos]) == "--decoupled") + { + useDecoupledModel = true; + } + ++argPos; + + // read the initial time step and the end time + double tEnd, dt; + std::istringstream(argv[argPos++]) >> tEnd; + std::istringstream(argv[argPos++]) >> dt; + + if (argc - argPos > 1) { + usage(argv[0]); + } + + double refineLevel = 0; + if (argc - argPos == 1) + std::istringstream(argv[argPos]) >> refineLevel; + + //////////////////////////////////////////////////////////// + // create the grid + //////////////////////////////////////////////////////////// + GlobalPosition lowerLeftCorner(0.0); + GlobalPosition domainSize; + Dune::array< unsigned int, dim > numberOfCells; // cell resolution + domainSize[0] = 6.0; + domainSize[1] = 4.0; + + numberOfCells[0] = 48; + numberOfCells[1] = 32; + + Dune::shared_ptr<Grid> grid(Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeftCorner, domainSize, numberOfCells)); + grid->globalRefine(refineLevel); + + //////////////////////////////////////////////////////////// + // instantiate and run the concrete problem + //////////////////////////////////////////////////////////// + + // specify dimensions of the low-permeable lens + GlobalPosition lowerLeftLens, upperRightLens; + lowerLeftLens[0] = 1.0; + lowerLeftLens[1] = 2.0; + upperRightLens[0] = 4.0; + upperRightLens[1] = 3.0; + + // instantiate and run the concrete problem + if (useBoxModel) + { + BoxTimeManager timeManager; + Dumux::GeneralLensProblem<BoxTypeTag> problem(timeManager, grid->leafView(), lowerLeftLens, upperRightLens); + problem.setName("generallens_box"); + timeManager.init(problem, 0, dt, tEnd, true); + timeManager.run(); + return 0; + } + else if (useDecoupledModel) + { + DecoupledTimeManager timeManager; + Dumux::GeneralLensProblem<DecoupledTypeTag> problem(timeManager, grid->leafView(), lowerLeftLens, upperRightLens); + problem.setName("generallens_decoupled"); + timeManager.init(problem, 0, dt, tEnd, true); + timeManager.run(); + return 0; + } + else + { + std::cout<<"No valid model chosen!\n"; + usage(argv[0]); + } +#ifdef NDEBUG + } + catch (Dune::Exception &e) { + std::cerr << "Dune reported error: " << e << std::endl; + } + catch (...) { + std::cerr << "Unknown exception thrown!\n"; + throw; + } +#endif // NDEBUG + + return 3; +} diff --git a/test/decoupled/1p/benchmarkresult.hh b/test/decoupled/1p/benchmarkresult.hh index d09a89ea21..99754c4c7b 100644 --- a/test/decoupled/1p/benchmarkresult.hh +++ b/test/decoupled/1p/benchmarkresult.hh @@ -436,7 +436,7 @@ public: sumf += volume*sourceVec[0]; // get the absolute permeability - Dune::FieldMatrix<double,dim,dim> K = problem.spatialParameters().intrinsicPermeability(global, element); + Dune::FieldMatrix<double,dim,dim> K = problem.spatialParameters().intrinsicPermeability(element); int isIdx = -1; Dune::FieldVector<Scalar,2*dim> fluxVector; diff --git a/test/decoupled/1p/test_1p_problem.hh b/test/decoupled/1p/test_1p_problem.hh index ec12033ca2..796b552d51 100644 --- a/test/decoupled/1p/test_1p_problem.hh +++ b/test/decoupled/1p/test_1p_problem.hh @@ -38,7 +38,7 @@ #include <dumux/decoupled/1p/diffusion/diffusionproblem1p.hh> #include <dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh> -#include "test_diffusion_spatialparams.hh" +#include "test_1p_spatialparams.hh" namespace Dumux { @@ -77,7 +77,7 @@ private: typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; public: - typedef Dumux::TestDiffusionSpatialParams<TypeTag> type; + typedef Dumux::TestOnePSpatialParams<TypeTag> type; }; // Enable gravity @@ -89,6 +89,8 @@ SET_TYPE_PROP(TestProblemOneP, Model, Dumux::FVVelocity1P<TypeTag>); //Set the problem SET_TYPE_PROP(TestProblemOneP, Problem, Dumux::TestProblemOneP<TTAG(TestProblemOneP)>); + +SET_INT_PROP(TestProblemOneP, LSVerbosity, 1); } /*! @@ -97,10 +99,9 @@ SET_TYPE_PROP(TestProblemOneP, Problem, Dumux::TestProblemOneP<TTAG(TestProblemO * \brief test problem for the decoupled one-phase model. */ template<class TypeTag = TTAG(TestProblemOneP)> -class TestProblemOneP: public DiffusionProblem1P<TypeTag, TestProblemOneP<TypeTag> > +class TestProblemOneP: public DiffusionProblem1P<TypeTag > { - typedef TestProblemOneP<TypeTag> ThisType; - typedef DiffusionProblem1P<TypeTag, ThisType> ParentType; + typedef DiffusionProblem1P<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid; diff --git a/test/decoupled/1p/test_diffusion_problem.hh b/test/decoupled/1p/test_diffusion_problem.hh index eeb08d1546..4311b58f06 100644 --- a/test/decoupled/1p/test_diffusion_problem.hh +++ b/test/decoupled/1p/test_diffusion_problem.hh @@ -57,8 +57,7 @@ class TestDiffusionProblem; ////////// namespace Properties { -NEW_TYPE_TAG(DiffusionTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties)) - ; +NEW_TYPE_TAG(DiffusionTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties, TestDiffusionSpatialParams)); // Set the grid type SET_PROP(DiffusionTestProblem, Grid) @@ -66,20 +65,8 @@ SET_PROP(DiffusionTestProblem, Grid) typedef Dune::SGrid<2, 2> type; }; -SET_PROP(DiffusionTestProblem, PressurePreconditioner) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix; - typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) Vector; -public: - typedef Dune::SeqILUn<Matrix, Vector, Vector> type; -}; - -//SET_INT_PROP(DiffusionTestProblem, VelocityFormulation, -// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW); - -SET_INT_PROP(DiffusionTestProblem, PressureFormulation, - GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal); +SET_INT_PROP(DiffusionTestProblem, Formulation, + DecoupledTwoPCommonIndices::pGlobalSw); // Set the wetting phase SET_PROP(DiffusionTestProblem, WettingPhase) @@ -99,35 +86,23 @@ public: typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type; }; -// Set the spatial parameters -SET_PROP(DiffusionTestProblem, SpatialParameters) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; - -public: - typedef Dumux::TestDiffusionSpatialParams<TypeTag> type; -}; - // Enable gravity SET_BOOL_PROP(DiffusionTestProblem, EnableGravity, false); // set the types for the 2PFA FV method NEW_TYPE_TAG(FVVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem)); -SET_INT_PROP(FVVelocity2PTestProblem, IterationNumberPreconditioner, 0); SET_TYPE_PROP(FVVelocity2PTestProblem, Model, Dumux::FVVelocity2P<TTAG(FVVelocity2PTestProblem)>); SET_TYPE_PROP(FVVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVVelocity2PTestProblem)>); // set the types for the MPFA-O FV method NEW_TYPE_TAG(FVMPFAOVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem)); -SET_INT_PROP(FVMPFAOVelocity2PTestProblem, IterationNumberPreconditioner, 1); +SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, LinearSolver, Dumux::ILUnBiCGSTABBackend<TypeTag>); +SET_INT_PROP(FVMPFAOVelocity2PTestProblem, PreconditionerIterations, 1); SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Model, Dumux::FVMPFAOVelocity2P<TypeTag>); SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVMPFAOVelocity2PTestProblem)>); // set the types for the mimetic FD method NEW_TYPE_TAG(MimeticPressure2PTestProblem, INHERITS_FROM(DiffusionTestProblem, Mimetic)); -SET_INT_PROP(MimeticPressure2PTestProblem, IterationNumberPreconditioner, 0); SET_TYPE_PROP(MimeticPressure2PTestProblem, Model, Dumux::MimeticPressure2P<TTAG(MimeticPressure2PTestProblem)>); SET_TYPE_PROP(MimeticPressure2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(MimeticPressure2PTestProblem)>); } @@ -157,9 +132,11 @@ class TestDiffusionProblem: public DiffusionProblem2P<TypeTag> enum { - wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, - eqIdxPress = Indices::pressureEq, - eqIdxSat = Indices::saturationEq + wPhaseIdx = Indices::wPhaseIdx, + nPhaseIdx = Indices::nPhaseIdx, + pGlobalIdx = Indices::pGlobalIdx, + SwIdx = Indices::SwIdx, + pressEqIdx = Indices::pressEqIdx, }; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; @@ -236,8 +213,8 @@ public: //! set dirichlet condition (saturation [-]) void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const { - values[eqIdxPress] = exact(globalPos); - values[eqIdxSat] = 1.0; + values[pGlobalIdx] = exact(globalPos); + values[SwIdx] = 1.0; } //! set neumann condition for phases (flux, [kg/(m^2 s)]) diff --git a/test/decoupled/1p/test_diffusion_spatialparams.hh b/test/decoupled/1p/test_diffusion_spatialparams.hh index 770e3f2ca1..a43b184cd4 100644 --- a/test/decoupled/1p/test_diffusion_spatialparams.hh +++ b/test/decoupled/1p/test_diffusion_spatialparams.hh @@ -25,20 +25,44 @@ #ifndef TEST_DIFFUSION_SPATIALPARAMETERS_HH #define TEST_DIFFUSION_SPATIALPARAMETERS_HH - +#include <dumux/material/spatialparameters/fvspatialparameters.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> namespace Dumux { +//forward declaration +template<class TypeTag> +class TestDiffusionSpatialParams; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(TestDiffusionSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(TestDiffusionSpatialParams, SpatialParameters, Dumux::TestDiffusionSpatialParams<TypeTag>); + +// Set the material law +SET_PROP(TestDiffusionSpatialParams, MaterialLaw) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef LinearMaterial<Scalar> RawMaterialLaw; +public: + typedef EffToAbsLaw<RawMaterialLaw> type; +}; +} + /*! * \ingroup IMPETtests * \brief spatial parameters for the test problem for diffusion models. */ template<class TypeTag> -class TestDiffusionSpatialParams +class TestDiffusionSpatialParams: public FVSpatialParameters<TypeTag> { + typedef FVSpatialParameters<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; @@ -49,20 +73,13 @@ class TestDiffusionSpatialParams typedef typename Grid::Traits::template Codim<0>::Entity Element; typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; - typedef Dune::FieldVector<CoordScalar, dim> LocalPosition; typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; - typedef LinearMaterial<Scalar> RawMaterialLaw; public: - typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; - void update (Scalar saturationW, const Element& element) - { - - } - - const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const + const FieldMatrix& intrinsicPermeabilityAtPos (const GlobalPosition& globalPos) const { double rt = globalPos[0]*globalPos[0]+globalPos[1]*globalPos[1]; permeability_[0][0] = (delta_*globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])/rt; @@ -72,14 +89,14 @@ public: return permeability_; } - double porosity(const GlobalPosition& globalPos, const Element& element) const + double porosity(const Element& element) const { return 0.2; } // return the parameter object for the Brooks-Corey material law which depends on the position - const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const + const MaterialLawParams& materialLawParams(const Element &element) const { return materialLawParams_; } @@ -90,7 +107,7 @@ public: } TestDiffusionSpatialParams(const GridView& gridView) - : permeability_(0) + : ParentType(gridView), permeability_(0) { // residual saturations materialLawParams_.setSwr(0.0); diff --git a/test/decoupled/2p/test_impes_problem.hh b/test/decoupled/2p/test_impes_problem.hh index 245484a4aa..342645def5 100644 --- a/test/decoupled/2p/test_impes_problem.hh +++ b/test/decoupled/2p/test_impes_problem.hh @@ -37,18 +37,16 @@ #include <dumux/material/fluidsystems/liquidphase.hh> #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/components/oil.hh> -#include <dumux/linear/impetbicgstabilu0solver.hh> #include <dumux/decoupled/2p/impes/impesproblem2p.hh> #include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh> -#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh> #include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh> #include <dumux/decoupled/2p/transport/fv/capillarydiffusion.hh> #include <dumux/decoupled/2p/transport/fv/gravitypart.hh> #include "test_impes_spatialparams.hh" -//#include<dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh> +#include<dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh> namespace Dumux { @@ -61,7 +59,7 @@ class TestIMPESProblem; ////////// namespace Properties { -NEW_TYPE_TAG(IMPESTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties, Transport)); +NEW_TYPE_TAG(IMPESTestProblem, INHERITS_FROM(DecoupledTwoP, Transport, TestIMPESSpatialParams)); // Set the grid type SET_PROP(IMPESTestProblem, Grid) @@ -71,32 +69,22 @@ SET_PROP(IMPESTestProblem, Grid) }; // Set the problem property -SET_PROP(IMPESTestProblem, Problem) -{ -public: - typedef Dumux::TestIMPESProblem<TTAG(IMPESTestProblem)> type; -}; +SET_TYPE_PROP(IMPESTestProblem, Problem, Dumux::TestIMPESProblem<TTAG(IMPESTestProblem)>); // Set the model properties -SET_PROP(IMPESTestProblem, TransportModel) -{ - typedef Dumux::FVSaturation2P<TTAG(IMPESTestProblem)> type; -}; +SET_TYPE_PROP(IMPESTestProblem, TransportModel, Dumux::FVSaturation2P<TTAG(IMPESTestProblem)>); + SET_TYPE_PROP(IMPESTestProblem, DiffusivePart, Dumux::CapillaryDiffusion<TypeTag>); SET_TYPE_PROP(IMPESTestProblem, ConvectivePart, Dumux::GravityPart<TypeTag>); -SET_TYPE_PROP(IMPESTestProblem, LinearSolver, Dumux::IMPETBiCGStabILU0Solver<TypeTag>); SET_PROP(IMPESTestProblem, PressureModel) { typedef Dumux::FVVelocity2P<TTAG(IMPESTestProblem)> type; -// typedef Dumux::FVMPFAOVelocity2P<TTAG(IMPESTestProblem)> type; }; -//SET_INT_PROP(IMPESTestProblem, VelocityFormulation, -// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW); +//SET_INT_PROP(IMPESTestProblem, Formulation, +// DecoupledTwoPCommonIndices::pnSn); -//SET_INT_PROP(IMPESTestProblem, PressureFormulation, -// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal); // Set the wetting phase SET_PROP(IMPESTestProblem, WettingPhase) @@ -116,21 +104,10 @@ public: typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type; }; -// Set the spatial parameters -SET_PROP(IMPESTestProblem, SpatialParameters) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; - -public: - typedef Dumux::TestIMPESSpatialParams<TypeTag> type; -}; - // Enable gravity SET_BOOL_PROP(IMPESTestProblem, EnableGravity, false); -//SET_TYPE_PROP(IMPESTestProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>); +SET_TYPE_PROP(IMPESTestProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>); SET_SCALAR_PROP(IMPESTestProblem, CFLFactor, 0.95); } @@ -169,8 +146,10 @@ enum enum { wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, - eqIdxPress = Indices::pressureEq, - eqIdxSat = Indices::saturationEq + pWIdx = Indices::pwIdx, + SwIdx = Indices::SwIdx, + eqIdxPress = Indices::pressEqIdx, + eqIdxSat = Indices::satEqIdx }; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; @@ -183,7 +162,7 @@ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes; typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes))::PrimaryVariables PrimaryVariables; public: -TestIMPESProblem(TimeManager& timeManager, const GridView &gridView) : +TestIMPESProblem(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView) { } @@ -213,7 +192,7 @@ bool shouldWriteRestartFile() const * * This problem assumes a temperature of 10 degrees Celsius. */ -Scalar temperatureAtPos(const GlobalPosition& globalPos) const +Scalar temperature(const Element& element) const { return 273.15 + 10; // -> 10°C } @@ -221,12 +200,12 @@ Scalar temperatureAtPos(const GlobalPosition& globalPos) const // \} //! Returns the reference pressure for evaluation of constitutive relations -Scalar referencePressureAtPos(const GlobalPosition& globalPos) const +Scalar referencePressure(const Element& element) const { return 1e5; // -> 10°C } -void sourceAtPos(PrimaryVariables &values,const GlobalPosition& globalPos) const +void source(PrimaryVariables &values,const Element& element) const { values = 0; } @@ -270,18 +249,18 @@ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) c FluidState fluidState; fluidState.update(sat, pRef, pRef, temp); - values[eqIdxPress] = (2e5 + (this->bboxMax()[dim-1] - globalPos[dim-1]) * FluidSystem::phaseDensity(wPhaseIdx, temp, pRef, fluidState) * this->gravity().two_norm()); + values[pWIdx] = (2e5 + (this->bboxMax()[dim-1] - globalPos[dim-1]) * FluidSystem::phaseDensity(wPhaseIdx, temp, pRef, fluidState) * this->gravity().two_norm()); } else { - values[eqIdxPress] = 2e5; + values[pWIdx] = 2e5; } - values[eqIdxSat] = 0.8; + values[SwIdx] = 0.8; } else { - values[eqIdxPress] = 2e5; - values[eqIdxSat] = 0.2; + values[pWIdx] = 2e5; + values[SwIdx] = 0.2; } } @@ -295,11 +274,11 @@ void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) con } } //! return initial solution -> only saturation values have to be given! -void initialAtPos(PrimaryVariables &values, - const GlobalPosition &globalPos) const +void initial(PrimaryVariables &values, + const Element& element) const { - values[eqIdxPress] = 0; - values[eqIdxSat] = 0.2; + values[pWIdx] = 0; + values[SwIdx] = 0.2; } private: diff --git a/test/decoupled/2p/test_impes_spatialparams.hh b/test/decoupled/2p/test_impes_spatialparams.hh index 1ab82602a0..65430df4ab 100644 --- a/test/decoupled/2p/test_impes_spatialparams.hh +++ b/test/decoupled/2p/test_impes_spatialparams.hh @@ -25,7 +25,7 @@ #ifndef TEST_IMPES_SPATIALPARAMETERS_HH #define TEST_IMPES_SPATIALPARAMETERS_HH - +#include <dumux/material/spatialparameters/fvspatialparameters.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> @@ -33,17 +33,41 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class TestIMPESSpatialParams; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(TestIMPESSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(TestIMPESSpatialParams, SpatialParameters, Dumux::TestIMPESSpatialParams<TypeTag>); + +// Set the material law +SET_PROP(TestIMPESSpatialParams, MaterialLaw) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; +public: + typedef EffToAbsLaw<RawMaterialLaw> type; +}; +} + /*! * * \ingroup IMPETtests * \brief spatial parameters for the sequential 2p test */ template<class TypeTag> -class TestIMPESSpatialParams +class TestIMPESSpatialParams: public FVSpatialParameters<TypeTag> { typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef FVSpatialParameters<TypeTag> ParentType; typedef typename Grid::ctype CoordScalar; enum @@ -51,40 +75,34 @@ class TestIMPESSpatialParams typedef typename Grid::Traits::template Codim<0>::Entity Element; typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; - typedef Dune::FieldVector<CoordScalar, dim> LocalPosition; typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; -// typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; - typedef LinearMaterial<Scalar> RawMaterialLaw; public: - typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; - void update (Scalar saturationW, const Element& element) - { - } - - const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const + Scalar intrinsicPermeability (const Element& element) const { - return constPermeability_; + return 1.0e-7; } - double porosity(const GlobalPosition& globalPos, const Element& element) const + double porosity(const Element& element) const { return 0.2; } // return the parameter object for the Brooks-Corey material law which depends on the position - const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const +// const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + const MaterialLawParams& materialLawParams(const Element& element) const { return materialLawParams_; } TestIMPESSpatialParams(const GridView& gridView) - : constPermeability_(0) + : ParentType(gridView) { // residual saturations materialLawParams_.setSwr(0.2); @@ -92,27 +110,18 @@ public: // // parameters for the Brooks-Corey Law // // entry pressures -// materialLawParams_.setPe(0); -// + materialLawParams_.setPe(0); // // Brooks-Corey shape parameters -// materialLawParams_.setLambda(2); + materialLawParams_.setLambda(2); // parameters for the linear // entry pressures function - materialLawParams_.setEntryPC(0); - materialLawParams_.setMaxPC(0); - - for(int i = 0; i < dim; i++) - { - constPermeability_[i][i] = 1e-7; - } - +// materialLawParams_.setEntryPC(0); +// materialLawParams_.setMaxPC(0); } private: MaterialLawParams materialLawParams_; - FieldMatrix constPermeability_; - }; } // end namespace diff --git a/test/decoupled/2p/test_transport_problem.hh b/test/decoupled/2p/test_transport_problem.hh index 54587edf89..a0836bff70 100644 --- a/test/decoupled/2p/test_transport_problem.hh +++ b/test/decoupled/2p/test_transport_problem.hh @@ -49,7 +49,7 @@ class TestTransportProblem; ////////// namespace Properties { -NEW_TYPE_TAG(TransportTestProblem, INHERITS_FROM(DecoupledTwoP, Transport)); +NEW_TYPE_TAG(TransportTestProblem, INHERITS_FROM(DecoupledTwoP, Transport, TestTransportSpatialParams)); // Set the grid type @@ -89,16 +89,8 @@ public: typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type; }; -// Set the spatial parameters -SET_PROP(TransportTestProblem, SpatialParameters) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; - typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; - -public: - typedef Dumux::TestTransportSpatialParams<TypeTag> type; -}; +SET_INT_PROP(TransportTestProblem, VelocityFormulation, + DecoupledTwoPCommonIndices::velocityTotal); // Disable gravity SET_BOOL_PROP(TransportTestProblem, EnableGravity, false); @@ -144,7 +136,8 @@ class TestTransportProblem: public TransportProblem2P<TypeTag> enum { - wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, eqIdxSat = Indices::saturationEq + wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, + satEqIdx = Indices::satEqIdx }; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; diff --git a/test/decoupled/2p/test_transport_spatialparams.hh b/test/decoupled/2p/test_transport_spatialparams.hh index 25454ff8c2..7e394f57a9 100644 --- a/test/decoupled/2p/test_transport_spatialparams.hh +++ b/test/decoupled/2p/test_transport_spatialparams.hh @@ -25,7 +25,7 @@ #ifndef TEST_TRANSPORT_SPATIALPARAMETERS_HH #define TEST_TRANSPORT_SPATIALPARAMETERS_HH - +#include <dumux/material/spatialparameters/fvspatialparameters.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> @@ -33,55 +33,70 @@ namespace Dumux { +//forward declaration +template<class TypeTag> +class TestTransportSpatialParams; + +namespace Properties +{ +// The spatial parameters TypeTag +NEW_TYPE_TAG(TestTransportSpatialParams); + +// Set the spatial parameters +SET_TYPE_PROP(TestTransportSpatialParams, SpatialParameters, Dumux::TestTransportSpatialParams<TypeTag>); + +// Set the material law +SET_PROP(TestTransportSpatialParams, MaterialLaw) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; + typedef LinearMaterial<Scalar> RawMaterialLaw; +public: + typedef EffToAbsLaw<RawMaterialLaw> type; +}; +} + /*! * \ingroup IMPETtests * \brief spatial parameters for the explicit transport test */ template<class TypeTag> -class TestTransportSpatialParams +class TestTransportSpatialParams: public FVSpatialParameters<TypeTag> { + typedef FVSpatialParameters<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid; typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView; typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar; typedef typename Grid::ctype CoordScalar; enum - {dim=Grid::dimension, dimWorld=Grid::dimensionworld, numEq=1}; + {dim=Grid::dimension, dimWorld=Grid::dimensionworld}; typedef typename Grid::Traits::template Codim<0>::Entity Element; - typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; - typedef Dune::FieldVector<CoordScalar, dim> LocalPosition; - typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix; - -// typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; - typedef LinearMaterial<Scalar> RawMaterialLaw; public: - typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams; - void update (Scalar saturationW, const Element& element) - {} - - const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const + Scalar intrinsicPermeability (const Element& element) const { - return constPermeability_; + return 1e-5; } - double porosity(const GlobalPosition& globalPos, const Element& element) const + double porosity(const Element& element) const { return 0.2; } // return the parameter object for the Brooks-Corey material law which depends on the position - const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const + const MaterialLawParams& materialLawParams(const Element &element) const { return materialLawParams_; } TestTransportSpatialParams(const GridView& gridView) - : constPermeability_(0) + : ParentType(gridView) { // residual saturations materialLawParams_.setSwr(0.0); @@ -90,18 +105,10 @@ public: // parameters for the linear entry pressures function materialLawParams_.setEntryPC(0); materialLawParams_.setMaxPC(0); - - for(int i = 0; i < dim; i++) - { - constPermeability_[i][i] = 1e-5; - } - } private: MaterialLawParams materialLawParams_; - FieldMatrix constPermeability_; - }; } // end namespace -- GitLab