Commit 03d95e3c authored by Markus Wolff's avatar Markus Wolff
Browse files

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
parent f16ca59c
......@@ -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();
}
/*!
......
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
/*****************************************************************************
* 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;
}
}
############################################################
# 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
/*****************************************************************************
* 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
/*****************************************************************************
* Copyright (C) 2012 by Markus Wolff *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *