From 03d95e3cff5f70e8a1a627147ce2fe5f7e068f81 Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Wed, 15 Aug 2012 16:54:56 +0000
Subject: [PATCH] added new test for decoupled 2p mpfa models + adapted
existing mpfa tests to changed structure
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8884 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
test/decoupled/1p/test_diffusionproblem.hh | 3 +-
test/decoupled/2p/Makefile.am | 4 +-
test/decoupled/2p/test_mpfa2p.cc | 147 ++++++++
test/decoupled/2p/test_mpfa2p.input | 72 ++++
test/decoupled/2p/test_mpfa2pproblem.hh | 351 ++++++++++++++++++
test/decoupled/2p/test_mpfa2pspatialparams.hh | 283 ++++++++++++++
6 files changed, 858 insertions(+), 2 deletions(-)
create mode 100644 test/decoupled/2p/test_mpfa2p.cc
create mode 100644 test/decoupled/2p/test_mpfa2p.input
create mode 100644 test/decoupled/2p/test_mpfa2pproblem.hh
create mode 100644 test/decoupled/2p/test_mpfa2pspatialparams.hh
diff --git a/test/decoupled/1p/test_diffusionproblem.hh b/test/decoupled/1p/test_diffusionproblem.hh
index 263478c270..c5a20ba4b9 100644
--- a/test/decoupled/1p/test_diffusionproblem.hh
+++ b/test/decoupled/1p/test_diffusionproblem.hh
@@ -39,7 +39,7 @@
#include <dumux/material/components/unit.hh>
#include <dumux/decoupled/2p/diffusion/fv/fvpressureproperties2p.hh>
-#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressureproperties2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressurepropertiessimple2p.hh>
#include <dumux/decoupled/2p/diffusion/mimetic/mimeticpressureproperties2p.hh>
#include <dumux/decoupled/2p/diffusion/diffusionproblem2p.hh>
#include <dumux/decoupled/common/fv/fvvelocity.hh>
@@ -226,6 +226,7 @@ public:
this->variables().cellData(i).setSaturation(wPhaseIdx, 1.0);
}
this->model().initialize();
+ velocity_.initialize();
}
/*!
diff --git a/test/decoupled/2p/Makefile.am b/test/decoupled/2p/Makefile.am
index f1ea957ebc..a74e4c9961 100644
--- a/test/decoupled/2p/Makefile.am
+++ b/test/decoupled/2p/Makefile.am
@@ -1,4 +1,4 @@
-check_PROGRAMS = test_impes test_transport
+check_PROGRAMS = test_impes test_transport test_mpfa2p
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
@@ -7,4 +7,6 @@ test_impes_SOURCES = test_impes.cc
test_transport_SOURCES = test_transport.cc
+test_mpfa2p_SOURCES = test_mpfa2p.cc
+
include $(top_srcdir)/am/global-rules
diff --git a/test/decoupled/2p/test_mpfa2p.cc b/test/decoupled/2p/test_mpfa2p.cc
new file mode 100644
index 0000000000..26b6a00b2d
--- /dev/null
+++ b/test/decoupled/2p/test_mpfa2p.cc
@@ -0,0 +1,147 @@
+/*****************************************************************************
+ * Copyright (C) 2012 by Markus Wolff *
+ * 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/>. *
+ *****************************************************************************/
+
+#define STRUCTUREDGRID 1
+
+#include "config.h"
+
+#include "test_mpfa2pproblem.hh"
+#include <dumux/common/start.hh>
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0)
+ {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe List of Mandatory arguments for this program is:\n"
+ "\t-TEnd The end of the simulation. [s] \n"
+ "\t-DtInitial The initial timestep size. [s] \n"
+ "\t-Grid.NumberOfCellsX Resolution in x-direction [-]\n"
+ "\t-Grid.NumberOfCellsY Resolution in y-direction [-]\n"
+ "\t-Grid.UpperRightX Dimension of the grid [m]\n"
+ "\t-Grid.UpperRightY Dimension of the grid [m]\n";
+ errorMessageOut += "\n\nThe Optional command line argument:\n"
+ "\t-ModelType Can be: FV (standard finite volume), FVAdaptive (adaptive finite volume),\n"
+ "\t MPFAO (MPFA o-method), MPFAL (MPFA l-method), MPFALAdaptive (adaptive MPFA l-method)\n";
+ std::cout << errorMessageOut << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ Dune::ParameterTree paramTree;
+ std::string s(Dumux::readOptions_(argc, argv, paramTree));
+ if (s.empty())
+ {
+ if (paramTree.hasKey("ModelType"))
+ {
+ std::string modelType(paramTree.get<std::string>("ModelType"));
+
+ if (modelType == "FV")
+ {
+ typedef TTAG(FVTwoPTestProblem) ProblemTypeTag;
+ typedef GET_PROP(ProblemTypeTag, ParameterTree) ParamTree;
+ Dune::ParameterTree &rt = ParamTree::runTimeParams();
+ rt["ModelType"]=modelType;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Used standard finite volume model\n";
+ return startReturn;
+ }
+ else if (modelType == "FVAdaptive")
+ {
+ typedef TTAG(FVAdaptiveTwoPTestProblem) ProblemTypeTag;
+ typedef GET_PROP(ProblemTypeTag, ParameterTree) ParamTree;
+ Dune::ParameterTree &rt = ParamTree::runTimeParams();
+ rt["ModelType"]=modelType;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Used adaptive finite volume model\n";
+ return startReturn;
+ }
+ else if (modelType == "MPFAO")
+ {
+ typedef TTAG(MPFAOTwoPTestProblem) ProblemTypeTag;
+ typedef GET_PROP(ProblemTypeTag, ParameterTree) ParamTree;
+ Dune::ParameterTree &rt = ParamTree::runTimeParams();
+ rt["ModelType"]=modelType;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Used finite volume MPFA o-method model\n";
+ return startReturn;
+ }
+ else if (modelType == "MPFAL")
+ {
+ typedef TTAG(MPFALTwoPTestProblem) ProblemTypeTag;
+ typedef GET_PROP(ProblemTypeTag, ParameterTree) ParamTree;
+ Dune::ParameterTree &rt = ParamTree::runTimeParams();
+ rt["ModelType"]=modelType;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Used finite volume MPFA l-method model\n";
+ return startReturn;
+ }
+ else if (modelType == "MPFALAdaptive")
+ {
+ typedef TTAG(MPFALAdaptiveTwoPTestProblem) ProblemTypeTag;
+ typedef GET_PROP(ProblemTypeTag, ParameterTree) ParamTree;
+ Dune::ParameterTree &rt = ParamTree::runTimeParams();
+ rt["ModelType"]=modelType;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Used adaptive finite volume MPFA l-method model\n";
+ return startReturn;
+ }
+ else
+ {
+ typedef TTAG(MPFAOTwoPTestProblem) ProblemTypeTag;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"Unknwon model type "<<modelType<<" specified\n";
+ std::cout<<"Default to finite volume MPFA o-method model\n";
+ return startReturn;
+ }
+ }
+ else
+ {
+ typedef TTAG(MPFAOTwoPTestProblem) ProblemTypeTag;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<"No model type specified\n";
+ std::cout<<"Default to finite volume MPFA o-method model\n";
+ return startReturn;
+ }
+ }
+ else
+ {
+ typedef TTAG(MPFAOTwoPTestProblem) ProblemTypeTag;
+ int startReturn = Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ std::cout<<"######################################################\n";
+ std::cout<<s<<" is not a valid model type specification!\n";
+ std::cout<<"Default to finite volume MPFA o-method model\n";
+ return startReturn;
+ }
+}
diff --git a/test/decoupled/2p/test_mpfa2p.input b/test/decoupled/2p/test_mpfa2p.input
new file mode 100644
index 0000000000..2603d92ec5
--- /dev/null
+++ b/test/decoupled/2p/test_mpfa2p.input
@@ -0,0 +1,72 @@
+############################################################
+# Parameter file for test_mpfa2p.cc
+# Everything behind a '#' is a comment.
+# Type "./test_transport --help" for more information.
+############################################################
+
+############################################################
+# Mandatory arguments
+############################################################
+[TimeManager]
+TEnd = 5e4# [s]
+DtInitial = 0# [s]
+
+[Grid]
+#For nine-spot problem
+NumberOfCellsX = 10# [-] level 0 resolution in x-direction
+NumberOfCellsY = 10# [-] level 0 resolution in y-direction
+
+UpperRightX = 20# [m] dimension of the grid
+UpperRightY = 10# [m] dimension of the grid # 2D
+
+############################################################
+# Optional arguments
+############################################################
+[Problem]
+EnableGravity = true
+OutputInterval =0
+OutputTimeInterval = 1e4
+
+InletWidth = 2
+InjectionFlux = 0.1
+
+#OutputfileName =
+
+[SpatialParams]
+BackgroundPermeabilityXX = 1e-10
+BackgroundPermeabilityXY = 0
+BackgroundPermeabilityYX = 0
+BackgroundPermeabilityYY = 1e-10
+
+LensPermeabilityXX = 1e-14
+LensPermeabilityXY = 0
+LensPermeabilityYX = 0
+LensPermeabilityYY = 1e-14
+
+LensOneLowerLeftX = 7
+LensOneLowerLeftY = 6
+LensOneUpperRightX = 13
+LensOneUpperRightY = 7
+
+LensTwoLowerLeftX = 2
+LensTwoLowerLeftY = 4
+LensTwoUpperRightX = 8
+LensTwoUpperRightY = 5
+
+LensThreeLowerLeftX = 10
+LensThreeLowerLeftY = 2
+LensThreeUpperRightX = 18
+LensThreeUpperRightY = 3
+
+BackgroundEntryPressure = 500
+LenseEntryPressure = 5000
+
+[GridAdapt]
+MinLevel = 0# [-] minimum level of refinement
+MaxLevel = 1# [-] maximum level of refinement
+RefineTolerance = 0.2 # threshold for refinement criterion
+CoarsenTolerance = 0.01 # threshold for coarsening criterion
+EnableInitializationIndicator = true#
+RefineAtDirichletBC = false
+RefineAtFluxBC = true
+RefineAtSource = false
\ No newline at end of file
diff --git a/test/decoupled/2p/test_mpfa2pproblem.hh b/test/decoupled/2p/test_mpfa2pproblem.hh
new file mode 100644
index 0000000000..53fa2f3f8f
--- /dev/null
+++ b/test/decoupled/2p/test_mpfa2pproblem.hh
@@ -0,0 +1,351 @@
+/*****************************************************************************
+ * Copyright (C) 2012 by Markus Wolff *
+ * 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_TEST_MPFA2P_PROBLEM_HH
+#define DUMUX_TEST_MPFA2P_PROBLEM_HH
+
+#include <dune/grid/alugrid/2d/alugrid.hh>
+
+#include <dumux/common/cubegridcreator.hh>
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/simplednapl.hh>
+
+#include <dumux/decoupled/2p/diffusion/fv/fvpressureproperties2p.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvpressureproperties2padaptive.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfaopressureproperties2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfalpressureproperties2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfalpressureproperties2padaptive.hh>
+#include <dumux/decoupled/2p/transport/fv/fvtransportproperties2p.hh>
+#include <dumux/decoupled/2p/impes/impesproblem2p.hh>
+
+
+#include<dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh>
+#include<dumux/decoupled/2p/impes/gridadaptionindicator2plocal.hh>
+
+#include "test_mpfa2pspatialparams.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class MPFATwoPTestProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+
+NEW_TYPE_TAG(MPFATwoPTestProblem, INHERITS_FROM(Test2PSpatialParams));
+
+// Set the grid type
+SET_PROP(MPFATwoPTestProblem, Grid)
+{
+ typedef Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming> type;
+};
+
+// set the GridCreator property
+SET_TYPE_PROP(MPFATwoPTestProblem, GridCreator, CubeGridCreator<TypeTag>);
+
+
+// Set the problem property
+SET_PROP(MPFATwoPTestProblem, Problem)
+{
+public:
+ typedef Dumux::MPFATwoPTestProblem<TypeTag> type;
+};
+
+// Set the wetting phase
+SET_PROP(MPFATwoPTestProblem, 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(MPFATwoPTestProblem, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleDNAPL<Scalar> > type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(MPFATwoPTestProblem, ProblemEnableGravity, true);
+
+SET_TYPE_PROP(MPFATwoPTestProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>);
+SET_SCALAR_PROP(MPFATwoPTestProblem, ImpetCFLFactor, 1.0);
+SET_TYPE_PROP(MPFATwoPTestProblem, AdaptionIndicator, Dumux::GridAdaptionIndicator2PLocal<TypeTag>);
+
+NEW_TYPE_TAG(FVTwoPTestProblem, INHERITS_FROM(FVPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTestProblem));
+NEW_TYPE_TAG(FVAdaptiveTwoPTestProblem, INHERITS_FROM(FVPressureTwoPAdaptive, FVTransportTwoP, IMPESTwoPAdaptive, MPFATwoPTestProblem));
+NEW_TYPE_TAG(MPFAOTwoPTestProblem, INHERITS_FROM(FVMPFAOPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTestProblem));
+NEW_TYPE_TAG(MPFALTwoPTestProblem, INHERITS_FROM(FVMPFALPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTestProblem));
+NEW_TYPE_TAG(MPFALAdaptiveTwoPTestProblem, INHERITS_FROM(FVMPFALPressureTwoPAdaptive, FVTransportTwoP, IMPESTwoPAdaptive, MPFATwoPTestProblem));
+
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag>
+class MPFATwoPTestProblem: public IMPESProblem2P<TypeTag>
+{
+typedef IMPESProblem2P<TypeTag> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+typedef typename GET_PROP_TYPE(TypeTag, TwoPIndices) Indices;
+
+typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
+
+typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+typedef typename GET_PROP(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
+
+typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
+typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+
+enum
+{
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ pWIdx = Indices::pwIdx,
+ SwIdx = Indices::SwIdx,
+ eqIdxPress = Indices::pressEqIdx,
+ eqIdxSat = Indices::satEqIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+MPFATwoPTestProblem(TimeManager &timeManager,const GridView &gridView) :
+ParentType(timeManager, gridView)
+{
+ this->setGrid(GridCreator::grid());
+
+ Scalar inletWidth = 1.0;
+ if (ParameterTree::tree().hasKey("Problem.InletWidth"))
+ {
+ inletWidth = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InletWidth);
+ }
+ GlobalPosition inletCenter = this->bboxMax();
+ inletCenter[0] *= 0.5;
+
+ inletLeftCoord_ = inletCenter;
+ inletLeftCoord_[0] -=0.5*inletWidth;
+ inletRightCoord_ = inletCenter;
+ inletRightCoord_[0] +=0.5*inletWidth;
+
+ inFlux_ = 1e-4;
+ if (ParameterTree::tree().hasKey("Problem.InjectionFlux"))
+ {
+ inFlux_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, InjectionFlux);
+ }
+
+ int outputInterval = 0;
+ if (ParameterTree::tree().hasKey("Problem.OutputInterval"))
+ {
+ outputInterval = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval);
+ }
+ this->setOutputInterval(outputInterval);
+
+ Scalar outputTimeInterval = 1e6;
+ if (ParameterTree::tree().hasKey("Problem.OutputTimeInterval"))
+ {
+ outputTimeInterval = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, OutputTimeInterval);
+ }
+ this->setOutputTimeInterval(outputTimeInterval);
+}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+/*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+const char *name() const
+{
+ if (ParameterTree::tree().hasKey("Problem.OutputfileName"))
+ {
+ std::string fileName(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, OutputfileName));
+ return fileName.c_str();
+ }
+ else
+ {
+ return "testmpfa2p";
+ }
+}
+
+bool shouldWriteRestartFile() const
+{
+ return false;
+}
+
+
+/*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+Scalar temperatureAtPos(const GlobalPosition& globalPos) const
+{
+ return 273.15 + 10; // -> 10°C
+}
+
+// \}
+
+
+Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
+{
+ return 1e5; // -> 10°C
+}
+
+void source(PrimaryVariables &values,const Element& element) const
+{
+ values = 0;
+}
+
+void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
+{
+ if (isInlet(globalPos))
+ {
+ bcTypes.setNeumann(eqIdxPress);
+ bcTypes.setDirichlet(eqIdxSat);
+ }
+ else if (isBottom(globalPos) || isTop(globalPos))
+ {
+ bcTypes.setAllNeumann();
+ }
+ else
+ {
+ bcTypes.setDirichlet(eqIdxPress);
+ bcTypes.setOutflow(eqIdxSat);
+ }
+}
+
+void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
+{
+ Scalar pRef = referencePressureAtPos(globalPos);
+ Scalar temp = temperatureAtPos(globalPos);
+
+ values[pWIdx] = (1e5 - (this->bboxMax()- globalPos) * this->gravity() * WettingPhase::density(temp, pRef));
+ values[SwIdx] = 1.0;
+
+ if (isInlet(globalPos))
+ {
+ values[SwIdx] = 0.0;
+ }
+}
+
+void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
+{
+ values = 0;
+ if (isInlet(globalPos))
+ {
+ values[nPhaseIdx] = -inFlux_;
+ }
+
+}
+
+void initialAtPos(PrimaryVariables &values,
+ const GlobalPosition& globalPos) const
+{
+ values[pWIdx] = 0;
+ values[SwIdx] = 1.0;
+}
+
+private:
+
+bool isInlet(const GlobalPosition& globalPos) const
+{
+ if (!isTop(globalPos))
+ return false;
+
+ for (int i = 0; i < dim; i++)
+ {
+ if (globalPos[i] < inletLeftCoord_[i] - eps_)
+ return false;
+ if (globalPos[i] > inletRightCoord_[i] + eps_)
+ return false;
+ }
+ return true;
+}
+
+bool isTop(const GlobalPosition& globalPos) const
+{
+ if (dim == 2)
+ {
+ if (globalPos[1] > this->bboxMax()[1] - eps_)
+ return true;
+ }
+ if (dim == 3)
+ {
+ if (globalPos[2] > this->bboxMax()[2] - eps_)
+ return true;
+ }
+ return false;
+}
+
+bool isBottom(const GlobalPosition& globalPos) const
+{
+ if (dim == 2)
+ {
+ if (globalPos[1] < eps_)
+ return true;
+ }
+ if (dim == 3)
+ {
+ if (globalPos[2] < eps_)
+ return true;
+ }
+ return false;
+}
+
+Scalar inFlux_;
+GlobalPosition inletLeftCoord_;
+GlobalPosition inletRightCoord_;
+static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/decoupled/2p/test_mpfa2pspatialparams.hh b/test/decoupled/2p/test_mpfa2pspatialparams.hh
new file mode 100644
index 0000000000..c05699e21b
--- /dev/null
+++ b/test/decoupled/2p/test_mpfa2pspatialparams.hh
@@ -0,0 +1,283 @@
+/*****************************************************************************
+ * Copyright (C) 2012 by Markus Wolff *
+ * 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/>. *
+ *****************************************************************************/
+#ifndef TEST_2P_SPATIALPARAMETERS_HH
+#define TEST_2P_SPATIALPARAMETERS_HH
+
+#include <dumux/material/spatialparams/fvspatialparams.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class Test2PSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(Test2PSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(Test2PSpatialParams, SpatialParams, Dumux::Test2PSpatialParams<TypeTag>);
+
+// Set the material law
+SET_PROP(Test2PSpatialParams, MaterialLaw)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+public:
+ typedef EffToAbsLaw<RawMaterialLaw> type;
+};
+}
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class Test2PSpatialParams: public FVSpatialParams<TypeTag>
+{
+ typedef FVSpatialParams<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 = 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 typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
+
+public:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ const FieldMatrix& intrinsicPermeabilityAtPos(const GlobalPosition& globalPos) const
+ {
+ if (isLensOne(globalPos))
+ return permLenses_;
+ else if (isLensTwo(globalPos))
+ return permLenses_;
+ else if (isLensThree(globalPos))
+ return permLenses_;
+ else
+ return permBackground_;
+ }
+
+ Scalar porosity(const Element& element) const
+ {
+ return 0.4;
+ }
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ {
+ if (isLensOne(globalPos))
+ return materialLawParamsLenses_;
+ else if (isLensTwo(globalPos))
+ return materialLawParamsLenses_;
+ else if (isLensThree(globalPos))
+ return materialLawParamsLenses_;
+ else
+ return materialLawParamsBackground_;
+ }
+
+ Test2PSpatialParams(const GridView& gridView) :
+ ParentType(gridView), permBackground_(0), permLenses_(0),
+ lensOneLowerLeft_(0), lensOneUpperRight_(0), lensTwoLowerLeft_(0), lensTwoUpperRight_(0), lensThreeLowerLeft_(0), lensThreeUpperRight_(0)
+ {
+ // residual saturations
+ materialLawParamsBackground_.setSwr(0.);
+ materialLawParamsBackground_.setSnr(0.);
+
+ materialLawParamsLenses_.setSwr(0.);
+ materialLawParamsLenses_.setSnr(0.);
+
+ //parameters for Brooks-Corey law
+ // entry pressures function
+ materialLawParamsBackground_.setPe(0.);
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundEntryPressure"))
+ {
+ materialLawParamsBackground_.setPe(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundEntryPressure));
+ }
+
+ materialLawParamsLenses_.setPe(0.);
+ if (ParameterTree::tree().hasKey("SpatialParams.LenseEntryPressure"))
+ {
+ materialLawParamsLenses_.setPe(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LenseEntryPressure));
+ }
+
+ // Brooks-Corey shape parameters
+ materialLawParamsBackground_.setLambda(3);
+ materialLawParamsLenses_.setLambda(2);
+
+ permBackground_[0][0] = 1e-10;
+ permBackground_[1][1] = 1e-10;
+ permLenses_[0][0] = 1e-10;
+ permLenses_[1][1] = 1e-10;
+
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundPermeabilityXX"))
+ {
+ permBackground_[0][0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundPermeabilityXX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundPermeabilityXY"))
+ {
+ permBackground_[0][1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundPermeabilityXY);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundPermeabilityYX"))
+ {
+ permBackground_[1][0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundPermeabilityYX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundPermeabilityYY"))
+ {
+ permBackground_[1][1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundPermeabilityYY);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensPermeabilityXX"))
+ {
+ permLenses_[0][0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensPermeabilityXX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensPermeabilityXY"))
+ {
+ permLenses_[0][1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensPermeabilityXY);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensPermeabilityYX"))
+ {
+ permLenses_[1][0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensPermeabilityYX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensPermeabilityYY"))
+ {
+ permLenses_[1][1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensPermeabilityYY);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensOneLowerLeftX"))
+ {
+ lensOneLowerLeft_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensOneLowerLeftX);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensOneUpperRightX"))
+ {
+ lensOneUpperRight_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensOneUpperRightX);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensTwoLowerLeftX"))
+ {
+ lensTwoLowerLeft_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensTwoLowerLeftX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensTwoUpperRightX"))
+ {
+ lensTwoUpperRight_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensTwoUpperRightX);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensThreeLowerLeftX"))
+ {
+ lensThreeLowerLeft_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensThreeLowerLeftX);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensThreeUpperRightX"))
+ {
+ lensThreeUpperRight_[0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensThreeUpperRightX);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensOneLowerLeftY"))
+ {
+ lensOneLowerLeft_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensOneLowerLeftY);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensOneUpperRightY"))
+ {
+ lensOneUpperRight_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensOneUpperRightY);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensTwoLowerLeftY"))
+ {
+ lensTwoLowerLeft_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensTwoLowerLeftY);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensTwoUpperRightY"))
+ {
+ lensTwoUpperRight_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensTwoUpperRightY);
+ }
+
+ if (ParameterTree::tree().hasKey("SpatialParams.LensThreeLowerLeftY"))
+ {
+ lensThreeLowerLeft_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensThreeLowerLeftY);
+ }
+ if (ParameterTree::tree().hasKey("SpatialParams.LensThreeUpperRightY"))
+ {
+ lensThreeUpperRight_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensThreeUpperRightY);
+ }
+ }
+
+private:
+
+ bool isLensOne(const GlobalPosition& globalPos) const
+ {
+ for (int i = 0; i < dim; i++)
+ {
+ if (globalPos[i] < lensOneLowerLeft_[i] || globalPos[i] > lensOneUpperRight_[i])
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+ bool isLensTwo(const GlobalPosition& globalPos) const
+ {
+ for (int i = 0; i < dim; i++)
+ {
+ if (globalPos[i] < lensTwoLowerLeft_[i] || globalPos[i] > lensTwoUpperRight_[i])
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+ bool isLensThree(const GlobalPosition& globalPos) const
+ {
+ for (int i = 0; i < dim; i++)
+ {
+ if (globalPos[i] < lensThreeLowerLeft_[i] || globalPos[i] > lensThreeUpperRight_[i])
+ {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ MaterialLawParams materialLawParamsBackground_;
+ MaterialLawParams materialLawParamsLenses_;
+ FieldMatrix permBackground_;
+ FieldMatrix permLenses_;
+ GlobalPosition lensOneLowerLeft_;
+ GlobalPosition lensOneUpperRight_;
+ GlobalPosition lensTwoLowerLeft_;
+ GlobalPosition lensTwoUpperRight_;
+ GlobalPosition lensThreeLowerLeft_;
+ GlobalPosition lensThreeUpperRight_;
+};
+
+} // end namespace
+#endif
--
GitLab