Commit 2b71199c authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[ex3] implement executable solution

parent 0468f5a2
add_subdirectory(ex1)
add_subdirectory(ex2)
add_subdirectory(ex3)
// -*- 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 Tutorial problem for a fully coupled two phase-two component box model.
*/
#ifndef DUMUX_EXERCISE_THREE_B_PROBLEM_HH
#define DUMUX_EXERCISE_THREE_B_PROBLEM_HH
// The numerical model
#include <dumux/porousmediumflow/2p2c/model.hh>
// The box discretization
#include <dumux/discretization/box/properties.hh>
// The base porous media box problem
#include <dumux/porousmediumflow/problem.hh>
// Spatially dependent parameters
#include "spatialparams.hh"
// The fluid system that is created in this exercise
#include "fluidsystems/h2omycompressiblecomponent.hh"
namespace Dumux{
// Forward declaration of the problem class
template <class TypeTag> class ExerciseThreeProblemTwoPTwoC;
namespace Properties {
// Create a new type tag for the problem
NEW_TYPE_TAG(ExerciseThreeTwoPTwoCTypeTag, INHERITS_FROM(TwoPTwoC));
NEW_TYPE_TAG(ExerciseThreeBoxTwoPTwoCTypeTag, INHERITS_FROM(BoxModel, ExerciseThreeTwoPTwoCTypeTag));
// Set the "Problem" property
SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Problem, ExerciseThreeProblemTwoPTwoC<TypeTag>);
// Set the spatial parameters
SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, SpatialParams,
ExerciseThreeSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Set grid and the grid creator to be used
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
#elif HAVE_UG
SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2>);
#endif // HAVE_DUNE_ALUGRID
// The fluid system property
SET_PROP(ExerciseThreeTwoPTwoCTypeTag, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
};
}
/*!
* \ingroup TwoPBoxModel
*
* \brief Tutorial problem for a fully coupled two phase-two component box model.
*/
template <class TypeTag>
class ExerciseThreeProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
// Grid dimension
enum { dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// Dumux specific types
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
public:
ExerciseThreeProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, eps_(3e-6)
{
#if !(HAVE_DUNE_ALUGRID || HAVE_UG)
std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl;
#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG)
// initialize the fluid system
FluidSystem::init();
// set the depth of the bottom of the reservoir
depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the temperature \f$ K \f$
*/
Scalar temperature() const
{ return 283.15; }
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
if (globalPos[0] < eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
bcTypes.setAllDirichlet();
else
bcTypes.setAllNeumann();
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
// use initial as Dirichlet conditions
return initialAtPos(globalPos);
}
/*!
* \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.
*/
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
PrimaryVariables values(0.0);
Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
// extraction of oil (30 g/m/s) on a segment of the upper boundary
if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
{
// we solve for the mole balance, so we have to divide by the molar mass
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = -3e-2/FluidSystem::MyCompressibleComponent::molarMass();
}
else
{
// no-flow on the remaining Neumann-boundaries.
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = 0;
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
// tell the primary variables the phase state, i.e. which phase/phases
// is/are present, because this changes the meaning of the primary variable
// value at the index Indices::switchIdx
values.setState(Indices::firstPhaseOnly);
// use hydrostatic pressure distribution with 2 bar at the top and zero saturation
values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar
values[Indices::switchIdx] = 0.0;
return values;
}
// \}
/*!
* \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
*/
PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
{
// we do not define any sources
PrimaryVariables values(0.0);
return values;
}
private:
Scalar eps_; //! small epsilon value
Scalar depthBOR_; //! depth at the bottom of the reservoir
};
} // 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/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Tutorial problem for a fully coupled twophase box model.
*/
#ifndef DUMUX_EXERCISE_THREE_A_PROBLEM_HH
#define DUMUX_EXERCISE_THREE_A_PROBLEM_HH
// The numerical model
#include <dumux/porousmediumflow/2p/model.hh>
// The box discretization
#include <dumux/discretization/box/properties.hh>
// The grid managers
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#elif HAVE_UG
#include <dune/grid/uggrid.hh>
#else
#include <dune/grid/yaspgrid.hh>
#endif
// The porous media base problem
#include <dumux/porousmediumflow/problem.hh>
// Spatially dependent parameters
#include "spatialparams.hh"
// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>
// The components that will be created in this exercise
#include "components/myincompressiblecomponent.hh"
#include "components/mycompressiblecomponent.hh"
// We will only have liquid phases here
#include <dumux/material/fluidsystems/1pliquid.hh>
// The two-phase immiscible fluid system
#include <dumux/material/fluidsystems/2pimmiscible.hh>
namespace Dumux{
// Forward declaration of the problem class
template <class TypeTag> class ExerciseThreeProblemTwoP;
namespace Properties {
// Create a new type tag for the problem
NEW_TYPE_TAG(ExerciseThreeTwoPTypeTag, INHERITS_FROM(TwoP));
NEW_TYPE_TAG(ExerciseThreeBoxTwoPTypeTag, INHERITS_FROM(BoxModel, ExerciseThreeTwoPTypeTag));
// Set the "Problem" property
SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Problem, ExerciseThreeProblemTwoP<TypeTag>);
// Set the spatial parameters
SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, SpatialParams,
ExerciseThreeSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Set grid to be used
#if HAVE_DUNE_ALUGRID
SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
#elif HAVE_UG
SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::YaspGrid<2>);
#endif // HAVE_DUNE_ALUGRID
// we use the immiscible fluid system here
SET_PROP(ExerciseThreeTwoPTypeTag, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using TabulatedH2O = Components::TabulatedComponent<Components::H2O<Scalar>>;
using WettingPhase = typename FluidSystems::OnePLiquid<Scalar, TabulatedH2O>;
/*!
* Uncomment first line and comment second line for using the incompressible component
* Uncomment second line and comment first line for using the compressible component
*/
using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyIncompressibleComponent<Scalar> >;
// using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyCompressibleComponent<Scalar> >;
public:
using type = typename FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonWettingPhase>;
};
}
/*!
* \ingroup TwoPBoxModel
* \brief Tutorial problem for a fully coupled twophase box model.
*/
template <class TypeTag>
class ExerciseThreeProblemTwoP : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
// Grid dimension
enum { dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// Dumux specific types
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
enum {
waterPressureIdx = Indices::pressureIdx,
naplSaturationIdx = Indices::saturationIdx,
contiWEqIdx = Indices::conti0EqIdx + FluidSystem::comp0Idx, // water transport equation index
contiNEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx // napl transport equation index
};
public:
ExerciseThreeProblemTwoP(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, eps_(3e-6)
{
#if !(HAVE_DUNE_ALUGRID || HAVE_UG)
std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl;
#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG)
// initialize the tables for the water properties
std::cout << "Initializing the tables for the water properties" << std::endl;
Components::TabulatedComponent<Components::H2O<Scalar>>::init(/*tempMin=*/273.15,
/*tempMax=*/623.15,
/*numTemp=*/100,
/*pMin=*/0.0,
/*pMax=*/20e6,
/*numP=*/200);
// set the depth of the bottom of the reservoir
depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the temperature \f$ K \f$
*/
Scalar temperature() const
{ return 283.15; }
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
if (globalPos[0] < eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
bcTypes.setAllDirichlet();
else
bcTypes.setAllNeumann();
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
// use the initial values as Dirichlet values
return initialAtPos(globalPos);
}
/*!
* \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.
*/
PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
PrimaryVariables values(0.0);
Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
// influx of oil (30 g/m/s) over a segment of the top boundary
if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
{
values[contiWEqIdx] = 0;
values[contiNEqIdx] = -3e-2;
}
else
{
// no-flow on the remaining Neumann-boundaries.
values[contiWEqIdx] = 0;
values[contiNEqIdx] = 0;
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
// use hydrostatic pressure distribution with 2 bar at the top and zero saturation
values[waterPressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]);
values[naplSaturationIdx] = 0.0;
return values;
}
// \}
/*!
* \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
*/
PrimaryVariables sourceAtPos(const GlobalPosition& globalPos) const
{
// we define no source term here
PrimaryVariables values(0.0);
return values;
}
private:
Scalar eps_; //! small epsilon value
Scalar depthBOR_; //! depth at the bottom of the reservoir
};
}
#endif
# executables for exercise part a & b
#part a: 2pproblem
dune_add_test(NAME exercise3_a_solution
SOURCES exercise3.cc
COMPILE_DEFINITIONS TYPETAG=ExerciseThreeBoxTwoPTypeTag
CMD_ARGS exercise3_a.input)
#part b: 2p2cproblem
dune_add_test(NAME exercise3_b_solution
SOURCES exercise3.cc
COMPILE_DEFINITIONS TYPETAG=ExerciseThreeBoxTwoPTwoCTypeTag
CMD_ARGS exercise3_b.input)
# add tutorial to the common target
add_dependencies(test_exercises exercise3_a exercise3_b)
# add symlinks for the input files
add_input_file_links()
// -*- 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 Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a.
*/
#ifndef DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH
#define DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH
namespace Dumux
{
namespace BinaryCoeff
{
/*!
* \brief Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a
* The implementation of the missing methods in this file is part of exercise 3b.
*/