Skip to content

GitLab

Commit 8c865a7d authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[2p] make incompressible tests the new isothermal 2p tests

parent 7a0a17b8
Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/2p/implicit/properties.hh
View file @ 8c865a7d
......@@ -154,9 +154,6 @@ public:
typedef ThermalConductivitySomerton<Scalar, Indices> type;
};
}
}
......
test/porousmediumflow/2p/implicit/CMakeLists.txt
View file @ 8c865a7d
......@@ -3,20 +3,6 @@ add_subdirectory(incompressible)
add_input_file_links()
# isothermal tests
add_dumux_test(test_box2p test_box2p test_box2p.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/lensbox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lensbox-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_box2p")
add_dumux_test(test_cc2p test_cc2p test_cc2p.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/lenscc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lenscc-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_cc2p")
add_dumux_test(test_ccadaptive2p test_ccadaptive2p test_ccadaptive2p.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
......@@ -110,10 +96,8 @@ generalizeddirichletspatialparams.hh
injectionproblem2pni.hh
lensproblem.hh
lensspatialparams.hh
test_box2p.cc
test_box2pni.cc
test_boxadaptive2p.cc
test_cc2p.cc
test_cc2pcornerpoint.cc
test_cc2pni.cc
test_ccadaptive2p.cc
......
test/porousmediumflow/2p/implicit/incompressible/CMakeLists.txt
View file @ 8c865a7d
......@@ -4,17 +4,29 @@ dune_symlink_to_source_files(FILES "test_2p.input")
dune_add_test(NAME test_2p_incompressible_tpfa
SOURCES test_2p_cc.cc
COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleTpfa
CMD_ARGS test_2p.input -Problem.Name 2p_tpfa)
COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/lenscc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/2p_tpfa-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_tpfa test_2p.input -Problem.Name 2p_tpfa")
# using box
dune_add_test(NAME test_2p_incompressible_box
SOURCES test_2p_box.cc
CMD_ARGS test_2p.input -Problem.Name 2p_box)
COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/lensbox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/2p_box-00007.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_box test_2p.input -Problem.Name 2p_box")
# using mpfa
dune_add_test(NAME test_2p_incompressible_mpfa
SOURCES test_2p_cc.cc
COMPILE_DEFINITIONS TYPETAG=TwoPIncompressibleMpfa
CMD_ARGS test_2p.input -Problem.Name 2p_mpfa)
COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/lenscc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/2p_mpfa-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_incompressible_mpfa test_2p.input -Problem.Name 2p_mpfa")
set(CMAKE_BUILD_TYPE Release)
test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh
View file @ 8c865a7d
......@@ -100,35 +100,23 @@ public:
}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
* \brief Returns the scalar intrinsic permeability \f$[m^2]\f$
*
* \param element The element
* \param scv The sub control volume
* \param elemSol The element solution vector
* \return the intrinsic permeability
* \param globalPos The global position
*/
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
{
if (isInLens_(scv.dofPosition()))
if (isInLens_(globalPos))
return lensK_;
else
return outerK_;
return outerK_;
}
/*!
* \brief Function for defining the porosity.
* That is possibly solution dependent.
* \brief Returns the porosity \f$[-]\f$
*
* \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 porosity
* \param globalPos The global position
*/
Scalar porosity(const Element &element,
const SubControlVolume &scv,
const ElementSolutionVector &elemSol) const
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
/*!
......
test/porousmediumflow/2p/implicit/incompressible/test_2p_cc.cc
View file @ 8c865a7d
......@@ -122,9 +122,8 @@ int main(int argc, char** argv) try
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(leafGridView.size(GridView::dimension));
SolutionVector x(leafGridView.size(0));
problem->applyInitialSolution(x);
auto xOld = x;
......
test/porousmediumflow/2p/implicit/lensproblem.hh deleted 100644 → 0
View file @ 7a0a17b8
// -*- 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 Soil contamination problem where DNAPL infiltrates a fully
* water saturated medium.
*/
#ifndef DUMUX_LENSPROBLEM_HH
#define DUMUX_LENSPROBLEM_HH
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/dnapl.hh>
#include <dumux/discretization/cellcentered/box/properties.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/2p/implicit/model.hh>
#include <dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh>
#include <dumux/porousmediumflow/2p/implicit/adaptionhelper.hh>
#include <dumux/implicit/adaptive/gridadaptinitializationindicator.hh>
#include "lensspatialparams.hh"
namespace Dumux
{
template <class TypeTag>
class LensProblem;
//////////
// Specify the properties for the lens problem
//////////
namespace Properties
{
NEW_PROP_TAG(AdaptiveGrid);
NEW_TYPE_TAG(LensProblem, INHERITS_FROM(TwoP, LensSpatialParams));
NEW_TYPE_TAG(LensBoxProblem, INHERITS_FROM(BoxModel, LensProblem));
NEW_TYPE_TAG(LensBoxAdaptiveProblem, INHERITS_FROM(BoxModel, LensProblem));
NEW_TYPE_TAG(LensCCProblem, INHERITS_FROM(CCTpfaModel, LensProblem));
NEW_TYPE_TAG(LensCCAdaptiveProblem, INHERITS_FROM(CCTpfaModel, LensProblem));
#if HAVE_UG
SET_TYPE_PROP(LensCCProblem, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(LensBoxProblem, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(LensBoxAdaptiveProblem, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(LensCCProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(LensBoxProblem, Grid, Dune::YaspGrid<2>);
#endif
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(LensCCAdaptiveProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
#endif
// Set the problem property
SET_TYPE_PROP(LensProblem, Problem, LensProblem<TypeTag>);
// Set the wetting phase
SET_PROP(LensProblem, WettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar> > type;
};
// Set the non-wetting phase
SET_PROP(LensProblem, NonwettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef FluidSystems::LiquidPhase<Scalar, DNAPL<Scalar> > type;
};
#if HAVE_DUNE_ALUGRID
SET_BOOL_PROP(LensCCAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
SET_TYPE_PROP(LensCCAdaptiveProblem, AdaptionHelper, TwoPAdaptionHelper<TypeTag>);
#endif
#if HAVE_UG
SET_BOOL_PROP(LensBoxAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
SET_TYPE_PROP(LensBoxAdaptiveProblem, AdaptionHelper, TwoPAdaptionHelper<TypeTag>);
#endif
NEW_PROP_TAG(BaseProblem);
SET_TYPE_PROP(LensBoxProblem, BaseProblem, PorousMediaProblem<TypeTag>);
SET_TYPE_PROP(LensCCProblem, BaseProblem, PorousMediaProblem<TypeTag>);
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(LensCCAdaptiveProblem, BaseProblem, PorousMediaProblem<TypeTag>);
#endif
#if HAVE_UG
SET_TYPE_PROP(LensBoxAdaptiveProblem, BaseProblem, PorousMediaProblem<TypeTag>);
#endif
}
/*!
* \ingroup TwoPModel
* \ingroup ImplicitTestProblems
* \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 TwoPModel.
*
* 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_box2p -parameterFile test_box2p.input</tt> or
* <tt>./test_cc2p -parameterFile test_cc2p.input</tt>
*/
template <class TypeTag >
class LensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
{
typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
enum {
// primary variable indices
pwIdx = Indices::pwIdx,
snIdx = Indices::snIdx,
// equation indices
contiNEqIdx = Indices::contiNEqIdx,
// phase indices
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
// world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
static const bool adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid);
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
LensProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
temperature_ = 273.15 + 20; // -> 20°C
name_ = getParam<std::string>("Problem.Name");
}
/*!
* \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 before the time integration
*/
void preTimeStep()
{
if(adaptiveGrid)
{
PrimaryVariables totalMass = calculateTotalMass();
std::cout << "Total mass before grid adaption: " << std::endl;
std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << std::endl;
ParentType::preTimeStep();
totalMass = calculateTotalMass();
std::cout << "Total mass after grid adaption: " << std::endl;
std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << std::endl;
}
else
{
ParentType::preTimeStep();
PrimaryVariables totalMass = calculateTotalMass();
std::cout << "Total mass: " << std::endl;
std::cout << "wPhaseIdx: " << totalMass[wPhaseIdx] << " " << "nPhaseIdx: " << totalMass[nPhaseIdx] << 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
*/
Sources sourceAtPos(const GlobalPosition &globalPos) const
{
return PrimaryVariables(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
*/
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;
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
Scalar height = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
Scalar alpha = 1 + 1.5/height;
Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
// hydrostatic pressure scaled by alpha
values[pwIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
values[snIdx] = 0.0;
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.
*/
NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
{
NeumannFluxes values(0.0);
if (onInlet_(globalPos)) {
values[contiNEqIdx] = -0.04; // kg / (m * s)
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \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;
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1e5);
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
// hydrostatic pressure
values[pwIdx] = 1e5 - densityW*this->gravity()[1]*depth;
values[snIdx] = 0.0;
return values;
}
// \}
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_;
}
bool onInlet_(const GlobalPosition &globalPos) const
{
Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
Scalar lambda = (this->fvGridGeometry().bBoxMax()[0] - globalPos[0])/width;
return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
}
PrimaryVariables calculateTotalMass()
{
PrimaryVariables totalMass(0);
for (const auto& element : elements(this->gridView(), Dune::Partitions::interior))
{
auto fvGeometry = localView(this->model().fvGridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(this->model().curGlobalVolVars());
elemVolVars.bindElement(element, fvGeometry, this->model().curSol());
for (auto&& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
totalMass[nPhaseIdx] += scv.volume()*volVars.density(nPhaseIdx)
*volVars.porosity()*volVars.saturation(nPhaseIdx);
totalMass[wPhaseIdx] += scv.volume()*volVars.density(wPhaseIdx)
*volVars.porosity()*volVars.saturation(wPhaseIdx);
}
}
// communicate global sum if we are running mpi parallel
if (this->gridView().comm().size() > 1)
totalMass = this->gridView().comm().sum(totalMass);
return totalMass;
}
Scalar temperature_;
static constexpr Scalar eps_ = 3e-6;
std::string name_;
};
} //end namespace
#endif
test/porousmediumflow/2p/implicit/lensspatialparams.hh deleted 100644 → 0
View file @ 7a0a17b8
// -*- 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 *