Commit 75891010 authored by Martin Schneider's avatar Martin Schneider
Browse files

[tests][IMPES] Implementation of first working Impes test

parent 5bb904ea
add_subdirectory("impes")
add_input_file_links()
dune_symlink_to_source_files(FILES grids)
......
dune_symlink_to_source_files(FILES "test_impes.input")
dune_add_test(NAME test_impes_fv
SOURCES test_impes_fv.cc
COMPILE_DEFINITIONS TYPETAG=TwoPTransport
CMD_ARGS --command "${CMAKE_CURRENT_BINARY_DIR}/test_impes_fv test_impes.input")
#install sources
install(FILES
problem.hh
spatialparams.hh
test_transport_fv.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/porousmediumflow/2p/sequential/impes)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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/>. *
*****************************************************************************/
/*!
* \ingroup TwoPTests
* \brief The properties for the incompressible 2p test
*/
#ifndef DUMUX_IMPES_PRESSURE_TEST_PROBLEM_HH
#define DUMUX_IMPES_PRESSURE_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/box/properties.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/cellcentered/mpfa/properties.hh>
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/2p/sequential/pressure/model.hh>
#include "pressurespatialparams.hh"
namespace Dumux {
// forward declarations
template<class TypeTag> class TwoPTestProblem;
namespace Properties {
NEW_TYPE_TAG(TwoPImpes, INHERITS_FROM(CCTpfaModel,Pressure));
// Set the grid type
SET_TYPE_PROP(TwoPImpes, Grid, Dune::YaspGrid<2>);
// Set the problem type
SET_TYPE_PROP(TwoPImpes, Problem, TwoPTestProblem<TypeTag>);
// Set the fluid system
SET_PROP(TwoPImpes, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
// Set the spatial parameters
SET_PROP(TwoPImpes, SpatialParams)
{
private:
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
};
// Enable caching
SET_BOOL_PROP(TwoPImpes, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(TwoPImpes, EnableGridFluxVariablesCache, false);
SET_BOOL_PROP(TwoPImpes, EnableFVGridGeometryCache, false);
} // end namespace Properties
/*!
* \ingroup TwoPTests
* \brief The incompressible 2p test problem.
*/
template<class TypeTag>
class TwoPTestProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
enum {
pressureEqIdx = Indices::pressureEqIdx,
pressureIdx = Indices::pressureIdx
};
public:
TwoPTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry) {}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment
*
* \param values Stores the value of the boundary type
* \param globalPos The global position
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param values Stores the Dirichlet values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} ] \f$
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(1.0e5);
if (onLeftBoundary_(globalPos))
values[pressureIdx] = 2.0e5;
return values;
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values Stores the Neumann values for the conservation equations in
* \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param globalPos The position of the integration point of the boundary segment.
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
return values;
}
/*!
* \brief Evaluates the initial values for a control volume
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[pressureIdx] = 1e5;
return values;
}
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
*
* This is not specific to the discretization. By default it just
* throws an exception so it must be overloaded by the problem if
* no energy equation is used.
*/
Scalar temperature() const
{
return 293.15; // 10°C
}
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
}
bool onRightBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
}
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_;
}
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
}
static constexpr Scalar eps_ = 1e-6;
};
} // end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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/>. *
*****************************************************************************/
/*!
* \ingroup TwoPTests
* \brief The spatial params for the incompressible 2p test
*/
#ifndef DUMUX_IMPES_PRESSURE_TEST_SPATIAL_PARAMS_HH
#define DUMUX_IMPES_PRESSURE_TEST_SPATIAL_PARAMS_HH
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/porousmediumflow/2p/boxmaterialinterfaceparams.hh>
namespace Dumux {
/*!
* \ingroup TwoPTests
* \brief The spatial params for the incompressible 2p test
*/
template<class FVGridGeometry, class Scalar>
class TwoPTestSpatialParams
: public FVSpatialParams<FVGridGeometry, Scalar, TwoPTestSpatialParams<FVGridGeometry, Scalar>>
{
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ThisType = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using EffectiveLaw = RegularizedBrooksCorey<Scalar>;
public:
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
TwoPTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
// residual saturations
lensMaterialParams_.setSwr(0.0);
lensMaterialParams_.setSnr(0.0);
outerMaterialParams_.setSwr(0.0);
outerMaterialParams_.setSnr(0.0);
// parameters for the Van Genuchten law
// alpha and n
lensMaterialParams_.setLambda(2.0);
lensMaterialParams_.setPe(0);
outerMaterialParams_.setLambda(2.0);
outerMaterialParams_.setPe(0);
lensK_ = 9.05e-12;
outerK_ = 4.6e-10;
wettingSaturation_.resize(fvGridGeometry->numDofs(),0.0);
}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
* In this test, we use element-wise distributed permeabilities.
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return permeability
*/
template<class ElementSolution>
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// do not use a less permeable lens in the test with inverted wettability
if (isInLens_(element.geometry().center()))
return lensK_;
return outerK_;
}
/*!
* \brief Returns the porosity \f$[-]\f$
*
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
/*!
* \brief Returns the parameter object for the Brooks-Corey material law.
* In this test, we use element-wise distributed material parameters.
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return the material parameters object
*/
template<class ElementSolution>
const MaterialLawParams& materialLawParams(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
// do not use different parameters in the test with inverted wettability
if (isInLens_(element.geometry().center()))
return lensMaterialParams_;
return outerMaterialParams_;
}
/*!
* \brief Function for defining which phase is to be considered as the wetting phase.
*
* \return the wetting phase index
* \param globalPos The global position
*/
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{
return FluidSystem::phase0Idx;
}
Scalar wettingSaturation(const Element& element,
const SubControlVolume& scv) const
{
return wettingSaturation_[scv.dofIndex()];
}
template<class Scvf>
Scalar wettingSaturation(const Element& element,
const SubControlVolume& scv,
const Scvf& scvf) const
{
return wettingSaturation_[scvf.outsideScvIdx()];
}
void setWettingSaturation(const std::vector<Scalar>& sat)
{ wettingSaturation_ = sat;}
private:
bool isInLens_(const GlobalPosition &globalPos) const
{
for (int i = 0; i < dimWorld; ++i) {
if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
return false;
}
return true;
}
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
Scalar lensK_;
Scalar outerK_;
MaterialLawParams lensMaterialParams_;
MaterialLawParams outerMaterialParams_;
std::vector<Scalar> wettingSaturation_;
static constexpr Scalar eps_ = 1.5e-7;
};
} // end namespace Dumux
#endif
[TimeLoop]
DtInitial = 1# [s]
MaxTimeStepSize = 10
TEnd = 400 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 100 10
Cells = 100 1
[SpatialParams]
LensLowerLeft = 20.0 0.0 # [m] coordinates of the lower left lens corner
LensUpperRight = 40.0 1.0 # [m] coordinates of the upper right lens corner
[Problem]
Name = test_impes
EnableGravity = false
[Component]
LiquidDensity = 1000
LiquidKinematicViscosity = 1e-5
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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 porousmedium flow model
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include "pressureproblem.hh"
#include "transportproblem.hh"
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
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-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.LowerLeft Lower left corner coordinates\n"
"\t-Grid.UpperRight Upper right corner coordinates\n"
"\t-Grid.Cells Number of cells in respective coordinate directions\n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
"\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
"\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
std::cout << errorMessageOut
<< "\n";
}
}
int main(int argc, char** argv) try
{
using namespace Dumux;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
using TwoPImpesTT = TTAG(TwoPImpes);
using TwoPTransportTT = TTAG(TwoPTransport);
// we simply extract the grid creator from one of the type tags
using GridManager = Dumux::GridManager<typename GET_PROP_TYPE(TwoPImpesTT, Grid)>;
GridManager gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometries
using TwoPImpesFVGridGeometry = typename GET_PROP_TYPE(TwoPImpesTT, FVGridGeometry);
using TwoPTransportFVGridGeometry = typename GET_PROP_TYPE(TwoPTransportTT, FVGridGeometry);
auto twoPImpesFvGridGeometry = std::make_shared<TwoPImpesFVGridGeometry>(leafGridView);
auto twoPTransportFvGridGeometry = std::make_shared<TwoPTransportFVGridGeometry>(leafGridView);
twoPImpesFvGridGeometry->update();
twoPTransportFvGridGeometry->update();
// the problems (boundary conditions)
using TwoPImpes = typename GET_PROP_TYPE(TwoPImpesTT, Problem);
using TwoPTransportProblem = typename GET_PROP_TYPE(TwoPTransportTT, Problem);
auto twoPImpesProblem = std::make_shared<TwoPImpes>(twoPImpesFvGridGeometry);
auto twoPTransportProblem = std::make_shared<TwoPTransportProblem>(twoPTransportFvGridGeometry);
// the solution vector
using SVTwoPImpes = typename GET_PROP_TYPE(TwoPImpesTT, SolutionVector);
SVTwoPImpes x_p(twoPImpesFvGridGeometry->numDofs());
twoPImpesProblem->applyInitialSolution(x_p);
auto x_p_old = x_p;
using SVTwoPTransport = typename GET_PROP_TYPE(TwoPTransportTT, SolutionVector);
SVTwoPTransport x_s(twoPTransportFvGridGeometry->numDofs());
twoPTransportProblem->applyInitialSolution(x_s);
auto x_s_old = x_s;
// the grid variables
using GVTwoPImpes = typename GET_PROP_TYPE(TwoPImpesTT, GridVariables);
auto gridVarTwoPImpes = std::make_shared<GVTwoPImpes>(twoPImpesProblem, twoPImpesFvGridGeometry);
gridVarTwoPImpes->init(x_p, x_p_old);