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

[2p][tests] Add SPE10 test case

parent d5db5983
add_subdirectory(pointsources)
add_subdirectory(spe10)
add_input_file_links()
# isothermal tests
......
add_input_file_links()
add_dumux_test(test_spe10 test_spe10 test_spe10.cc
${CMAKE_CURRENT_BINARY_DIR}/test_spe10)
\ No newline at end of file
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* 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/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Implementation of quadratic relative-permeability-saturation relations.
* @author Markus Wolff
*/
#ifndef DUMUX_QUADRATICLAW_HH
#define DUMUX_QUADRATICLAW_HH
#include <algorithm>
namespace Dumux
{
/*!
* \ingroup fluidmatrixinteractionslaws
*
* \brief Implementation of quadratic relative-permeability-saturation relations.
*
* This class provides quadratic relative-permeability-saturation relations:
*
* \f[
* k_{r_w} = \overline{S}_w^2
* \f]
*
* \f[
* k_{r_n} = \overline{S}_n^2 = (1 - \overline{S}_n)^2
* \f]
*
* A class providing functions for capillary-pressure-saturation relations has to be given as template argument.
*
* @tparam ScalarT The type of a scalar
* @tparam PCLawT The class providing the capillary-pressure-saturation relations
*
*/
template <class ScalarT, class PCLawT>
class QuadraticLaw
{
public:
typedef PCLawT PCLaw;
typedef typename PCLaw::Params Params;
typedef ScalarT Scalar;
/*!
* \brief The capillary pressure-saturation curve.
*
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Capillary pressure.
*/
static Scalar pc(const Params &params, Scalar Swe)
{
return PCLaw::pc(params, Swe);
}
DUNE_DEPRECATED_MSG("use pc() (uncapitalized 'c') instead")
static Scalar pC(const Params &params, Scalar Swe)
{
return PCLaw::pc(params, Swe);
}
/*!
* \brief The saturation-capillary pressure curve.
*
*
* \param pC Capillary pressure \f$p_C\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Effective wetting phase saturation.
*/
static Scalar sw(const Params &params, Scalar pC)
{
return PCLaw::sw(params, pC);
}
DUNE_DEPRECATED_MSG("use sw() (uncapitalized 's') instead")
static Scalar Sw(const Params &params, Scalar pC)
{
return PCLaw::sw(params, pC);
}
/*!
* \brief The partial derivative of the capillary
* pressure w.r.t. the effective saturation.
*
* \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of \f$p_c\f$ w.r.t. effective saturation.
*/
static Scalar dpc_dsw(const Params &params, Scalar Swe)
{
return PCLaw::dpc_dsw(params, Swe);
}
DUNE_DEPRECATED_MSG("use dpc_dsw() (uncapitalized 'c', 's') instead")
static Scalar dpC_dSw(const Params &params, Scalar Swe)
{
return PCLaw::dpc_dsw(params, Swe);
}
/*!
* \brief The partial derivative of the effective
* saturation w.r.t. the capillary pressure.
*
* \param pC Capillary pressure \f$p_C\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of effective saturation w.r.t. \f$p_c\f$.
*/
static Scalar dsw_dpc(const Params &params, Scalar pC)
{
return PCLaw::dsw_dpc(params, pC);
}
DUNE_DEPRECATED_MSG("use dsw_dpc() (uncapitalized 's', 'c') instead")
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return PCLaw::dsw_dpc(params, pC);
}
/*!
* \brief The relative permeability for the wetting phase of
* the medium.
*
* \param Swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the wetting phase.
*/
static Scalar krw(const Params &params, Scalar Swe)
{
if (Swe <= 0.0)
return 0.0;
else if (Swe >= 1.0)
return 1.0;
return Swe*Swe;
};
/*!
* \brief The derivative of the relative permeability for the
* wetting phase with regard to the wetting saturation of the
* medium.
*
* \param Swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the wetting phase w.r.t. effective wetting phase saturation calculated as implied by Brooks & Corey.
*/
static Scalar dkrw_dsw(const Params &params, Scalar Swe)
{
assert(0 <= Swe && Swe <= 1);
return 2*Swe;
};
DUNE_DEPRECATED_MSG("use dkrw_dsw() (uncapitalized 's') instead")
static Scalar dkrw_dSw(const Params &params, Scalar Swe)
{
assert(0 <= Swe && Swe <= 1);
return 2*Swe;
};
/*!
* \brief The relative permeability for the non-wetting phase of
* the medium.
*
* \param Swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the non-wetting phase.
*/
static Scalar krn(const Params &params, Scalar Swe)
{
if (Swe >= 1.0)
return 0.0;
else if (Swe <= 0.0)
return 1.0;
return (1-Swe)*(1-Swe);
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium.
*
* \param Swe The mobile saturation of the wetting phase.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the non-wetting phase w.r.t. effective wetting phase saturation calculated as implied by Brooks & Corey.
*/
static Scalar dkrn_dsw(const Params &params, Scalar Swe)
{
assert(0 <= Swe && Swe <= 1);
return -2*(1-Swe);
};
DUNE_DEPRECATED_MSG("use dkrn_dsw() (uncapitalized 's') instead")
static Scalar dkrn_dSw(const Params &params, Scalar Swe)
{
assert(0 <= Swe && Swe <= 1);
return -2*(1-Swe);
};
};
}
#endif // QUADRATICLAW_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*\ingroup Components
* \brief Rough estimate of some oil for testing purposes.
*/
#ifndef DUMUX_OIL_HH
#define DUMUX_OIL_HH
#include <dumux/material/components/component.hh>
namespace Dumux
{
/*!
* \ingroup Components
*
* \brief Rough estimate for testing purposes of some oil.
*
* \tparam Scalar The type used for scalar values
*/
template <class Scalar>
class Oil : public Component<Scalar, Oil<Scalar> >
{
public:
/*!
* \brief A human readable name for the water.
*/
static const char *name()
{ return "Oil"; }
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
static bool liquidIsCompressible()
{ return false; }
/*!
* \brief Rough estimate of the density of oil \f$\mathrm{[kg/m^3]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
return 890;
}
/*!
* \brief Rough estimate of the viscosity of oil in \f$\mathrm{[Pa*s]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{
return 5e-3;
};
};
} // end namespace
#endif
This diff is collapsed.
This diff is collapsed.
// -*- 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/>. *
*****************************************************************************/
#ifndef DUMUX_SPE10_PROBLEM_HH
#define DUMUX_SPE10_PROBLEM_HH
#include <dumux/material/components/simpleh2o.hh>
#include "spe10oil.hh"
#include <dumux/porousmediumflow/2p/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/implicit/cellcentered/propertydefaults.hh>
#include "spe10spatialparams.hh"
namespace Dumux
{
template <class TypeTag>
class CCTwoPSpe10Problem;
namespace Properties
{
NEW_TYPE_TAG(TwoPSpe10Problem, INHERITS_FROM(TwoP, Spe10SpatialParams));
NEW_TYPE_TAG(CCTwoPSpe10Problem, INHERITS_FROM(CCModel, TwoPSpe10Problem));
// Set the grid type
SET_TYPE_PROP(TwoPSpe10Problem, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(TwoPSpe10Problem, Problem, Dumux::CCTwoPSpe10Problem<TypeTag>);
SET_TYPE_PROP(TwoPSpe10Problem, LinearSolver, Dumux::SuperLUBackend<TypeTag> );
// Set the wetting phase
SET_PROP(TwoPSpe10Problem, WettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef FluidSystems::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
};
// Set the non-wetting phase
SET_PROP(TwoPSpe10Problem, NonwettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef FluidSystems::LiquidPhase<Scalar, Dumux::Oil<Scalar> > type;
};
// Enable gravity
SET_BOOL_PROP(TwoPSpe10Problem, ProblemEnableGravity, false);
}
template <class TypeTag >
class CCTwoPSpe10Problem : public ImplicitPorousMediaProblem<TypeTag>
{
typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum
{
// primary variable indices
pwIdx = Indices::pwIdx,
snIdx = Indices::snIdx,
// equation indices
contiNEqIdx = Indices::contiNEqIdx,
contiWEqIdx = Indices::contiWEqIdx,
// phase indices
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
// world dimension
dimWorld = GridView::dimensionworld,
dim = GridView::dimension
};
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dimWorld,dimWorld> DimWorldMatrix;
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
CCTwoPSpe10Problem(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, gridView)
{
eps_ = 3e-6;
temperature_ = 273.15 + 20; // -> 20°C
name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, EpisodeLength);
this->timeManager().startNextEpisode(episodeLength_);
pIn_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, PressureIn);
pOut_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, PressureOut);
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the problem name
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{
return name_;
}
/*!
* \brief User defined output after the time integration
*
* Will be called diretly after the time integration.
*/
void postTimeStep()
{
// Calculate storage terms
PrimaryVariables storage;
this->model().globalStorage(storage);
// Write mass balance information for rank 0
if (this->gridView().comm().rank() == 0) {
std::cout<<"Storage: " << storage << std::endl;
}
}
/*!
* \brief Returns the temperature \f$ K \f$
*
* This problem assumes a uniform temperature of 20 degrees Celsius.
*/
Scalar temperature() const
{ return temperature_; }
/*!
* \brief Returns the source term
*
* \param values Stores the source values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
* \param globalPos The global position
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0.0;
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \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
*/
void boundaryTypesAtPos(BoundaryTypes &bcTypes,
const GlobalPosition &globalPos) const
{
if (globalPos[1] < eps_ || globalPos[1] > this->bBoxMax()[1] - eps_)
{
bcTypes.setAllDirichlet();
}
else
{
bcTypes.setAllNeumann();
}
}
/*!
* \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
*/
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0;
if (globalPos[1] > this->bBoxMax()[1] - eps_)
{
values[pwIdx] = pOut_;
values[snIdx] = 1.0;
}
else
{
values[pwIdx] = pIn_;
values[snIdx] = 0.0;
}
}
/*!
* \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.
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
values = 0.0;