Commit 50e9c45f authored by Lisa Hug's avatar Lisa Hug
Browse files

update box-interface test cases

parent 328c50bc
......@@ -146,7 +146,7 @@ public:
satNFracture_ = priVars[saturationIdx];
satWFracture_ = 1.0 - satNFracture_;
}
bool useInterfaceCondition_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, UseInterfaceCondition);
//update fracture saturation using interface condition
const VertexMapper &vertexMapper = problem.vertexMapper();
int globalIdx = vertexMapper.subIndex(element, scvIdx, dim);
......@@ -170,40 +170,41 @@ public:
pe = MaterialLaw::pc(materialParamsMatrix, 1.0-materialParamsMatrix.snr());
}
FVElementGeometry minPcElemFvGeometry;
Element minPcElem = vertIdxToMinPcMapper.vertexElement(globalIdx);
minPcElemFvGeometry.update(problem.gridView(), minPcElem);
//choose materialLawParamsFracture as fracture has lower pe
// we set an scvIdx of 0, because we don't know the actual scv idx and assume element wise parameters anyway
MaterialLawParams minPcElemMaterialParams = problem.spatialParams().materialLawParamsFracture(minPcElem, minPcElemFvGeometry, /*dummy*/0);
//calculate capillary pressure based on the pc-sw relation of the minPc element
Scalar pcmin;
if (hasFractureFaces)
pcmin = MaterialLaw::pc(minPcElemMaterialParams, satWFracture_);
else
pcmin = MaterialLaw::pc(minPcElemMaterialParams, satWMatrix_);
// update fracture saturation for scvs containing fracture scvs, otherwise update matrix saturations
if (hasFractureFaces)
{
if (std::abs(pc-pcmin) < 1e-6) {}
else if (pcmin < pe)
satNFracture_ = std::min(materialParamsFracture.snr(), 0.0 /* SnInitial*/);
if (useInterfaceCondition_){
FVElementGeometry minPcElemFvGeometry;
Element minPcElem = vertIdxToMinPcMapper.vertexElement(globalIdx);
minPcElemFvGeometry.update(problem.gridView(), minPcElem);
//choose materialLawParamsFracture as fracture has lower pe
// we set an scvIdx of 0, because we don't know the actual scv idx and assume element wise parameters anyway
MaterialLawParams minPcElemMaterialParams = problem.spatialParams().materialLawParamsFracture(minPcElem, minPcElemFvGeometry, /*dummy*/0);
//calculate capillary pressure based on the pc-sw relation of the minPc element
Scalar pcmin;
if (hasFractureFaces)
pcmin = MaterialLaw::pc(minPcElemMaterialParams, satWFracture_);
else
satNFracture_ = 1.0 - MaterialLaw::sw(materialParamsFracture, pcmin);
satWFracture_ = 1.0 - satNFracture_;
}
else
{
if (std::abs(pc-pcmin) < 1e-6) {}
else if (pcmin < pe)
satNMatrix_ = std::min(materialParamsMatrix.snr(), 0.0 /* SnInitial*/);
pcmin = MaterialLaw::pc(minPcElemMaterialParams, satWMatrix_);
// update fracture saturation for scvs containing fracture scvs, otherwise update matrix saturations
if (hasFractureFaces)
{
if (std::abs(pc-pcmin) < 1e-6) {}
else if (pcmin < pe)
satNFracture_ = std::min(materialParamsFracture.snr(), 0.0 /* SnInitial*/);
else
satNFracture_ = 1.0 - MaterialLaw::sw(materialParamsFracture, pcmin);
satWFracture_ = 1.0 - satNFracture_;
}
else
satNMatrix_ = 1.0 - MaterialLaw::sw(materialParamsMatrix, pcmin);
satWMatrix_ = 1.0 - satNMatrix_;
{
if (std::abs(pc-pcmin) < 1e-6) {}
else if (pcmin < pe)
satNMatrix_ = std::min(materialParamsMatrix.snr(), 0.0 /* SnInitial*/);
else
satNMatrix_ = 1.0 - MaterialLaw::sw(materialParamsMatrix, pcmin);
satWMatrix_ = 1.0 - satNMatrix_;
}
}
// update fracture mobilities if we are on an scv containing fracture scvs
if (hasFractureFaces)
{
......@@ -226,8 +227,7 @@ public:
MaterialLaw::krn(materialParamsFracture, fluidStateFracture_.saturation(wPhaseIdx))
/ fluidStateFracture_.viscosity(nPhaseIdx);
// use interface condition - extended capillary pressure inteface condition
if (problem.useInterfaceCondition())
// use interface condition - extended capillary pressure inteface condition
// updated solution in the fracture is used for the interface condition
interfaceCondition(priVarsFracture, materialParamsMatrix, materialParamsFracture);
}
......
add_subdirectory(pointsources)
add_subdirectory(test1d_saturationcondition)
add_subdirectory(lensbox)
add_input_file_links()
# isothermal tests
......
add_input_file_links()
# isothermal tests
add_dumux_test(lensbox lensbox lensbox.cc
python ${CMAKE_SOURCE_DIR}/bin/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/satcondition-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/satcondition-00009.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test1d_satcondition")
#install sources
install(FILES
lensboxproblem.hh
lensboxparams.hh
lensbox.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/2p/lensbox)
// -*- 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 box model
*/
#include <config.h>
#include "lensboxproblem.hh"
#include <dumux/common/start.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 options for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
"\t-Problem.Name String for naming of the output files \n"
"\n";
std::cout << errorMessageOut << std::endl;
}
}
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv)
{
typedef TTAG(SatConditionBoxProblem) TypeTag;
return Dumux::start<TypeTag>(argc, argv, usage);
}
[TimeManager]
DtInitial = 100 # [s]
TEnd = 6000 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 6 4
Cells = 30 20
[SpatialParams]
LensLowerLeftX = 1.0 # [m] x-coordinate of the lower left lens corner
LensLowerLeftY = 2.0 # [m] y-coordinate of the lower left lens corner
LensUpperRightX = 4.0 # [m] x-coordinate of the upper right lens corner
LensUpperRightY = 3.0 # [m] y-coordinate of the upper right lens corner
[Problem]
Name = lensbox_classicbox # name passed to the output routines
SnInitial = 0
[Implicit]
EnablePartialReassemble = 1 # enable partial reassembly of the jacobian matrix?
EnableJacobianRecycling = 1 # Enable reuse of jacobian matrices?
UseInterfaceCondition = 1
// -*- 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_SatConditionProblem_HH
#define DUMUX_SatConditionProblem_HH
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/dnapl.hh>
#include <dumux/porousmediumflow/2p/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/implicit/cellcentered/propertydefaults.hh>
#include <dumux/porousmediumflow/2p/implicit/gridadaptindicator.hh>
#include <dumux/implicit/adaptive/gridadaptinitializationindicator.hh>
#include <dumux/implicit/box/vertidxtoscvneighbormapper.hh>
#include <dumux/porousmediumflow/2p/implicit/vertextominpcelemmapper.hh>
#include "lensboxspatialparams.hh"
namespace Dumux
{
template <class TypeTag>
class SatConditionProblem;
//////////
// Specify the properties for the satcondition problem
//////////
namespace Properties
{
NEW_TYPE_TAG(SatConditionProblem, INHERITS_FROM(TwoP, SatConditionSpatialParams));
NEW_TYPE_TAG(SatConditionBoxProblem, INHERITS_FROM(BoxModel, SatConditionProblem));
NEW_TYPE_TAG(SatConditionBoxAdaptiveProblem, INHERITS_FROM(BoxModel, SatConditionProblem));
NEW_TYPE_TAG(SatConditionCCProblem, INHERITS_FROM(CCModel, SatConditionProblem));
NEW_TYPE_TAG(SatConditionCCAdaptiveProblem, INHERITS_FROM(CCModel, SatConditionProblem));
#if HAVE_UG
SET_TYPE_PROP(SatConditionCCProblem, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(SatConditionBoxProblem, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(SatConditionCCProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(SatConditionBoxProblem, Grid, Dune::YaspGrid<2>);
#endif
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(SatConditionBoxAdaptiveProblem, Grid, Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>);
SET_TYPE_PROP(SatConditionCCAdaptiveProblem, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
#endif
// Set the problem property
SET_TYPE_PROP(SatConditionProblem, Problem, Dumux::SatConditionProblem<TypeTag>);
// Set the wetting phase
SET_PROP(SatConditionProblem, WettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Dumux::FluidSystems::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
};
// Set the non-wetting phase
SET_PROP(SatConditionProblem, NonwettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Dumux::FluidSystems::LiquidPhase<Scalar, Dumux::DNAPL<Scalar> > type;
};
// Linear solver settings
SET_TYPE_PROP(SatConditionCCProblem, LinearSolver, Dumux::ILU0BiCGSTABBackend<TypeTag> );
SET_TYPE_PROP(SatConditionBoxProblem, LinearSolver, Dumux::ILU0BiCGSTABBackend<TypeTag> );
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(SatConditionCCAdaptiveProblem, LinearSolver, Dumux::ILU0BiCGSTABBackend<TypeTag> );
SET_TYPE_PROP(SatConditionBoxAdaptiveProblem, LinearSolver, Dumux::ILU0BiCGSTABBackend<TypeTag> );
SET_BOOL_PROP(SatConditionCCAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(SatConditionCCAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
SET_TYPE_PROP(SatConditionCCAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
SET_BOOL_PROP(SatConditionBoxAdaptiveProblem, AdaptiveGrid, true);
SET_TYPE_PROP(SatConditionBoxAdaptiveProblem, AdaptionIndicator, TwoPImplicitGridAdaptIndicator<TypeTag>);
SET_TYPE_PROP(SatConditionBoxAdaptiveProblem, AdaptionInitializationIndicator, ImplicitGridAdaptInitializationIndicator<TypeTag>);
#endif
NEW_PROP_TAG(BaseProblem);
SET_TYPE_PROP(SatConditionBoxProblem, BaseProblem, ImplicitPorousMediaProblem<TypeTag>);
SET_TYPE_PROP(SatConditionCCProblem, BaseProblem, ImplicitPorousMediaProblem<TypeTag>);
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(SatConditionCCAdaptiveProblem, BaseProblem, ImplicitPorousMediaProblem<TypeTag>);
SET_TYPE_PROP(SatConditionBoxAdaptiveProblem, BaseProblem, ImplicitPorousMediaProblem<TypeTag>);
#endif
}
/*!
* \ingroup TwoPModel
* \ingroup ImplicitTestProblems
* \brief Soil contamination problem where DNAPL infiltrates a fully
* water saturated medium.
*
* The domain is a vertical column sized 0.1m times 0.5m and features a
* low permeablility zone which spans from a height of 0.155 m to 0.355m
* and is bounded above and below by a medium with higher permability.
* The problem is simulated quasi 1D so there is only one cell in
* x-direction.
*
* On the top, left and right of the domain neumann boundary
* conditions are used, while dirichlet conditions apply on the bottom.
*
* DNAPL is injected at the top boundary at a rate of
* 0.05 kg/(s m^2), the remaining neumann boundaries are no-flow
* boundaries.
*
*/
template <class TypeTag >
class SatConditionProblem : 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;
typedef typename GridView::Grid Grid;
enum {dim=Grid::dimension};
typedef typename Grid::template Codim<0>::Entity Element;
typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;
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
dimWorld = GridView::dimensionworld
};
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 Dumux::VertIdxToScvNeighborMapper<GridView> VertIdxToScvNeighborMapper;
typedef typename Dumux::VertIdxToMinPcMapper<TypeTag> VertIdxToMinPcMapper;
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
SatConditionProblem(TimeManager &timeManager,
const GridView &gridView)
: ParentType(timeManager, gridView),
vertIdxToScvNeighborMapper_(gridView),
vertIdxToMinPcMapper_(gridView, this->spatialParams())
{
eps_ = 3e-6;
temperature_ = 273.15 + 20; // -> 20°C
name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
std::string,
Problem,
Name);
vertIdxToMinPcMapper_.update();
vertIdxToScvNeighborMapper_.update();
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the vertIdxToScvNeighborMapper
*
* This is used as a prefix for files generated by the simulation.
*/
const VertIdxToScvNeighborMapper &vertIdxToScvNeighborMapper() const
{
return vertIdxToScvNeighborMapper_;
}
const VertIdxToMinPcMapper &vertIdxToMinPcMapper() const
{
return vertIdxToMinPcMapper_;
}
/*!
* \brief Returns the problem name
*
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
{
return name_.c_str();
}
/*!
* \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 &values,
const GlobalPosition &globalPos) const
{
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) {
values.setAllDirichlet();
}
else {
values.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
{
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1.0e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1.0e5);
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
Scalar height = this->bBoxMax()[1] - this->bBoxMin()[1];
Scalar depth = this->bBoxMax()[1] - globalPos[1];
Scalar alpha = 1 + 1.5/height;
Scalar width = this->bBoxMax()[0] - this->bBoxMin()[0];
Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
// hydrostatic pressure scaled by alpha
values[pwIdx] = 1.0e5 - factor*densityW*this->gravity()[1]*depth;
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;
if (onInlet_(globalPos)) {
values[contiNEqIdx] = -0.04; // kg / (m * s)
}
}
// \}
/*!
* \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
*/
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(FluidSystem::wPhaseIdx, /*pressure=*/1.0e5);
fluidState.setPressure(FluidSystem::nPhaseIdx, /*pressure=*/1.0e5);
Scalar densityW = FluidSystem::density(fluidState, FluidSystem::wPhaseIdx);
Scalar depth = this->bBoxMax()[1] - globalPos[1];
// hydrostatic pressure
values[pwIdx] = 1.0e5 - densityW*this->gravity()[1]*depth;
values[snIdx] = 0.0;
}
// \}
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{