...
 
Commits (62)
......@@ -27,14 +27,12 @@
#include <dumux/material/components/co2tablereader.hh>
namespace Dumux
{
namespace CO2TablesBenchmarkThree
{
namespace Dumux {
namespace CO2TablesBenchmarkThree {
// the real work is done by some external program which provides
// ready-to-use tables.
#include "co2values.inc"
}
}
} // end namespace CO2TablesBenchmarkThree
} // end namespace Dumux
#endif
#endif // DUMUX_HETEROGENEOUS_CO2TABLES_HH
......@@ -24,8 +24,7 @@
#ifndef DUMUX_VIPLAB_OUTPUT_HH
#define DUMUX_VIPLAB_OUTPUT_HH
namespace Dumux
{
namespace Dumux {
/**
* \brief Writes output to use with ViPLab.
*
......@@ -35,8 +34,8 @@ namespace Dumux
template<class TypeTag>
class ViplabOutput
{
using Problem = GET_PROP_TYPE(TypeTag, Problem);
using Scalar = GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum{ wPhaseIdx = Indices::wPhaseIdx };
......@@ -92,9 +91,9 @@ public:
{
setColor(color);
const auto cellNumber = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::vector<int>, Grid, Cells);
const auto cellNumber = getParam<std::vector<int>>("Grid.Cells");
assert(cellNumber[0] == solution.size());
upperRight_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::vector<double>, Grid, UpperRight);
upperRight_ = getParam<std::vector<double>>("Grid.UpperRight");
const Scalar discretizationLength = upperRight_[0] / cellNumber[0];
dataFile_.open(outputName_, std::fstream::app);
......@@ -108,8 +107,9 @@ public:
dataFile_.close();
}
void writeViplabOutput2d(const std::vector<double> &solution){
void writeViplabOutput2d(const std::vector<double> &solution)
{
// TODO: check if this is used anymore:
// double ymin= *min_element(solution.begin(),solution.end());
// Scalar zmax, zmin;
// zmax=problem_.variables().cellData(0).pressure()/(Fluid::density(0,0)*9.81);
......@@ -254,6 +254,7 @@ private:
std::string outputName_;
std::vector<double> upperRight_;
};
}//end namespace Dumux
#endif
#endif // DUMUX_VIPLAB_OUTPUT_HH
......@@ -5,21 +5,21 @@ add_dumux_test(lens2pexercise3 lens2pexercise3 lens2pexercise3.cc
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/lens-2p-exercise3-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lens-2p-00011.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2pexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input")
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2pexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input -Problem.Name lens-2p")
add_dumux_test(lens2p2cexercise3 lens2p2cexercise3 lens2p2cexercise3.cc
python ${dumux_INCLUDE_DIRS}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/lens-2p2c-exercise3-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lens-2p2c-00011.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2p2cexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input")
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2p2cexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input -Problem.Name lens-2p2c")
add_dumux_test(lens1p2cexercise3 lens1p2cexercise3 lens1p2cexercise3.cc
python ${dumux_INCLUDE_DIRS}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/lens-1p2c-exercise3-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lens-1p2c-00011.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/lens1p2cexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input")
--command "${CMAKE_CURRENT_BINARY_DIR}/lens1p2cexercise3 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise3.input -Problem.Name lens-1p2c")
# headers for installation and headercheck
install(FILES
......
......@@ -40,7 +40,7 @@
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
/*!
......@@ -52,9 +52,6 @@
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
......@@ -86,12 +83,10 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
int main(int argc, char** argv) try
{
// typedef TTAG(LensOnePTwoCProblem) TypeTag;
// return Dumux::start<TypeTag>(argc, argv, usage);
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(LensOnePTwoCProblem);
using TypeTag = TTAG(LensOnePTwoCProblemTypeTag);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -103,17 +98,19 @@ int main(int argc, char** argv) try
// 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 GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
/////////////////////////////////////////////////////////////////////
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -148,7 +145,7 @@ int main(int argc, char** argv) try
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
......@@ -208,8 +205,7 @@ int main(int argc, char** argv) try
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
catch (Dumux::ParameterException &e)
......@@ -217,6 +213,7 @@ catch (Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
......@@ -226,11 +223,13 @@ catch (Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
catch (Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
......
......@@ -22,21 +22,19 @@
* \brief Exercise to show the diffusive spreading of contaminants.
*/
#ifndef DUMUX_LENS_1P2C_PROBLEM_HH
#define DUMUX_LENS_1P2C_PROBLEM_HH
#include <dumux/porousmediumflow/1pnc/model.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/box/properties.hh>
// #include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include <dumux/material/fluidsystems/1padapter.hh>
#include "lens1p2cspatialparams.hh"
namespace Dumux
{
namespace Dumux {
template <class TypeTag>
class LensOnePTwoCProblem;
......@@ -44,29 +42,28 @@ class LensOnePTwoCProblem;
//////////
// Specify the properties for the lens problem
//////////
namespace Properties
{
NEW_TYPE_TAG(LensOnePTwoCProblem, INHERITS_FROM(BoxModel, OnePNC, Lens1p2cSpatialParams));
namespace Properties {
NEW_TYPE_TAG(LensOnePTwoCProblemTypeTag, INHERITS_FROM(BoxModel, OnePNC, Lens1p2cSpatialParamsTypeTag));
// Set the grid type
SET_TYPE_PROP(LensOnePTwoCProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(LensOnePTwoCProblemTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(LensOnePTwoCProblem, Problem, LensOnePTwoCProblem<TypeTag>);
SET_TYPE_PROP(LensOnePTwoCProblemTypeTag, Problem, LensOnePTwoCProblem<TypeTag>);
// Set fluid configuration
SET_PROP(LensOnePTwoCProblem, FluidSystem)
{private:
SET_PROP(LensOnePTwoCProblemTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
static const bool useComplexRelations = false;
public:
using type = FluidSystems::H2ON2<Scalar, useComplexRelations>;
using H2ON2 = FluidSystems::H2ON2<Scalar, FluidSystems::H2ON2DefaultPolicy</*simplified=*/true>>;
using type = FluidSystems::OnePAdapter<H2ON2, H2ON2::liquidPhaseIdx>;
};
// Define whether mole(true) or mass(false) fractions are used
SET_BOOL_PROP(LensOnePTwoCProblem, UseMoles, true);
SET_BOOL_PROP(LensOnePTwoCProblemTypeTag, UseMoles, true);
}
} // end namespace Properties
/*!
* \ingroup TwoPBoxProblems
......@@ -97,43 +94,33 @@ template <class TypeTag >
class LensOnePTwoCProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
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 Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
enum {
H2OIdx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::H2OIdx),
N2Idx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::N2Idx)
};
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
LensOnePTwoCProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry )
: ParentType(fvGridGeometry)
: ParentType(fvGridGeometry)
{
// the boundary condition data
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure");
......@@ -145,16 +132,6 @@ public:
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
std::string name() const
{
return "lens-1p2c";
}
/*!
* \brief Returns the temperature within the domain.
* This problem assumes a temperature of 10 degrees Celsius.
......@@ -162,7 +139,7 @@ public:
Scalar temperature() const
{
return 273.15 + 10; // -> 10 degrees Celsius
};
}
// \}
......@@ -184,6 +161,7 @@ public:
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -207,10 +185,12 @@ public:
{
values[Indices::pressureIdx] = upperPressure_;
}
else if (onLowerBoundary_(globalPos))
{
values[Indices::pressureIdx] = lowerPressure_;
}
return values;
}
......@@ -268,7 +248,7 @@ public:
const Scalar height = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
values[Indices::conti0EqIdx] = upperPressure_ - depth/height*(upperPressure_ - lowerPressure_);
values[FluidSystem::N2Idx] = (isInitial_(globalPos)) ? 9.12825137e-6 : 0.0;
values[N2Idx] = (isInitial_(globalPos)) ? 9.12825137e-6 : 0.0;
return values;
}
......@@ -312,7 +292,9 @@ private:
Scalar infiltrationEndTime_;
bool paraviewOutput_;
std::vector<int> numCells_;
};
} //end namespace
#endif
} //end namespace Dumux
#endif // DUMUX_LENS_1P2C_PROBLEM_HH
......@@ -28,20 +28,18 @@
* @author Bernd Flemisch, Klaus Mosthaf, Markus Wolff
*/
namespace Dumux
{
namespace Dumux {
//forward declaration
template<class FVGridGeometry, class Scalar>
class Lens1p2cSpatialParams;
namespace Properties
{
namespace Properties {
// The spatial parameters TypeTag
NEW_TYPE_TAG(Lens1p2cSpatialParams);
NEW_TYPE_TAG(Lens1p2cSpatialParamsTypeTag);
// Set the spatial parameters
SET_PROP(Lens1p2cSpatialParams, SpatialParams)
SET_PROP(Lens1p2cSpatialParamsTypeTag, SpatialParams)
{
private:
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -49,7 +47,8 @@ private:
public:
using type = Lens1p2cSpatialParams<FVGridGeometry, Scalar>;
};
}
} // end namespace Properties
/** \todo Please doc me! */
......@@ -58,12 +57,10 @@ class Lens1p2cSpatialParams : public FVSpatialParamsOneP<FVGridGeometry, Scalar,
{
using ThisType = Lens1p2cSpatialParams<FVGridGeometry, Scalar>;
using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, ThisType>;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
......@@ -109,6 +106,7 @@ public:
{
if (isInLens_(globalPos))
return lensK_;
return outerK_;
}
......@@ -116,6 +114,7 @@ public:
{
if (isInLens_(globalPos))
return lensPorosity_;
return outerPorosity_;
}
/*!
......@@ -125,8 +124,7 @@ public:
*/
//! Set the bounding box of the fine-sand lens
void setLensCoords(const GlobalPosition& lensLowerLeft,
const GlobalPosition& lensUpperRight)
void setLensCoords(const GlobalPosition& lensLowerLeft, const GlobalPosition& lensUpperRight)
{
lensLowerLeft_ = lensLowerLeft;
lensUpperRight_ = lensUpperRight;
......@@ -140,8 +138,8 @@ public:
* \param scvIdx The local index of the sub-control volume where
*/
Dune::FieldVector<Scalar,2> dispersivity(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
Dune::FieldVector<Scalar,2> values (0);
values[0] = longitudinalDispersivity_; // alpha_l
......@@ -155,14 +153,11 @@ public:
* \param vertexI DOC ME!
* \param vertexJ DOC ME!
*/
bool useTwoPointGradient(const Element &elem,
int vertexI,
int vertexJ) const
bool useTwoPointGradient(const Element &elem, int vertexI, int vertexJ) const
{
return false;
}
private:
/*!
* \brief DOC ME!
......@@ -170,10 +165,12 @@ private:
*/
bool isInLens_(const GlobalPosition &pos) const
{
for (int i = 0; i < dimWorld; ++i) {
for (int i = 0; i < dimWorld; ++i)
{
if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
return false;
}
return true;
}
......@@ -188,5 +185,6 @@ private:
Scalar transverseDispersivity_;
};
} // end namespace
#endif
} // end namespace Dumux
#endif // DUMUX_1P2C_SPATIALPARAMS_HH
......@@ -49,10 +49,9 @@
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
// #include "config.h"
#include "lens2p2cproblem.hh"
/*!
* \brief Provides an interface for customizing error messages associated with
......@@ -62,7 +61,6 @@
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
// TODO: do we need this void-function anymore?
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
......@@ -106,12 +104,10 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
int main(int argc, char** argv) try
{
// typedef TTAG(LensOnePTwoCProblem) TypeTag;
// return Dumux::start<TypeTag>(argc, argv, usage);
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(LensTwoPTwoCProblem);
using TypeTag = TTAG(LensTwoPTwoCProblemTypeTag);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -123,17 +119,19 @@ int main(int argc, char** argv) try
// 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 GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
/////////////////////////////////////////////////////////////////////
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -168,7 +166,7 @@ int main(int argc, char** argv) try
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
......@@ -228,6 +226,8 @@ int main(int argc, char** argv) try
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
catch (Dumux::ParameterException &e)
......@@ -235,6 +235,7 @@ catch (Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
......@@ -244,11 +245,13 @@ catch (Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
catch (Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
......
......@@ -28,7 +28,7 @@
#define DUMUX_LENS_2P2C_PROBLEM_HH
#include <dumux/porousmediumflow/2p2c/model.hh>
// #include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include "lens2pspatialparams.hh"
......@@ -36,11 +36,10 @@
#include <dumux/porousmediumflow/1pnc/model.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/box/properties.hh>
// #include <dumux/porousmediumflow/implicit/problem.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
namespace Dumux
{
namespace Dumux {
template<class TypeTag>
class LensTwoPTwoCProblem;
......@@ -48,25 +47,22 @@ class LensTwoPTwoCProblem;
//////////
// Specify the properties for the lens problem
//////////
namespace Properties
{
NEW_TYPE_TAG(LensTwoPTwoCProblem, INHERITS_FROM(TwoPTwoC, BoxModel, Lens2pSpatialParams));
namespace Properties {
NEW_TYPE_TAG(LensTwoPTwoCProblemTypeTag, INHERITS_FROM(TwoPTwoC, BoxModel, Lens2pSpatialParamsTypeTag));
// Set the grid type
SET_TYPE_PROP(LensTwoPTwoCProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(LensTwoPTwoCProblemTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(LensTwoPTwoCProblem, Problem, LensTwoPTwoCProblem<TypeTag>);
SET_TYPE_PROP(LensTwoPTwoCProblemTypeTag, Problem, LensTwoPTwoCProblem<TypeTag>);
// Set fluid configuration
SET_PROP(LensTwoPTwoCProblem, FluidSystem)
{private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
static const bool useComplexRelations = false;
public:
using type = FluidSystems::H2ON2<Scalar, useComplexRelations>;
};
}
SET_TYPE_PROP(LensTwoPTwoCProblemTypeTag,
FluidSystem,
FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/true>>);
} // end namespace Properties
/*!
* \ingroup TwoPTwoCModel
......@@ -87,35 +83,28 @@ class LensTwoPTwoCProblem: public PorousMediumFlowProblem<TypeTag>
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
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 SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
LensTwoPTwoCProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry )
: ParentType(fvGridGeometry)
{
......@@ -139,15 +128,6 @@ public:
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
std::string name() const
{
return "lens-2p2c";
}
/*!
* \brief Returns the temperature within the domain.
......@@ -181,6 +161,7 @@ public:
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -205,6 +186,7 @@ public:
values[Indices::pressureIdx] = upperPressure_;
values[Indices::switchIdx] = 0.0;
}
else if (onLowerBoundary_(globalPos))
{
values[Indices::pressureIdx] = lowerPressure_;
......@@ -275,6 +257,7 @@ public:
if (isInitial_(globalPos))
values.setState(Indices::bothPhases);
else
values.setState(Indices::firstPhaseOnly);
......@@ -289,9 +272,7 @@ public:
}
// \}
private:
Scalar eps_;
Scalar episodeLength_; // [s]
Scalar temperature_; // [K]
......@@ -322,10 +303,10 @@ private:
*/
bool isInitial_(const GlobalPosition &globalPos) const
{
return (globalPos[1] > 1.0 && globalPos[1] < 1.25 && globalPos[0] > 1.7
&& globalPos[0] < 2.3);
return (globalPos[1] > 1.0 && globalPos[1] < 1.25 && globalPos[0] > 1.7 && globalPos[0] < 2.3);
}
};
} //end namespace
#endif
} //end namespace Dumux
#endif // DUMUX_LENS_2P2C_PROBLEM_HH
......@@ -46,7 +46,7 @@
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
/*!
......@@ -57,7 +57,7 @@
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
// TODO is this void-function still in use?
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
......@@ -101,12 +101,10 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
int main(int argc, char** argv) try
{
// typedef TTAG(LensOnePTwoCProblem) TypeTag;
// return Dumux::start<TypeTag>(argc, argv, usage);
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(LensTwoPProblem);
using TypeTag = TTAG(LensTwoPProblemTypeTag);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -118,17 +116,19 @@ int main(int argc, char** argv) try
// 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 GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
/////////////////////////////////////////////////////////////////////
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -163,7 +163,7 @@ int main(int argc, char** argv) try
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
......@@ -223,8 +223,7 @@ int main(int argc, char** argv) try
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
catch (Dumux::ParameterException &e)
......@@ -232,6 +231,7 @@ catch (Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
......@@ -241,11 +241,13 @@ catch (Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
catch (Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
......
......@@ -38,8 +38,7 @@
#include "lens2pspatialparams.hh"
namespace Dumux
{
namespace Dumux {
template <class TypeTag>
class LensTwoPProblem;
......@@ -47,25 +46,26 @@ class LensTwoPProblem;
//////////
// Specify the properties for the lens problem
//////////
namespace Properties
{
NEW_TYPE_TAG(LensTwoPProblem, INHERITS_FROM(TwoP, BoxModel, Lens2pSpatialParams));
namespace Properties {
NEW_TYPE_TAG(LensTwoPProblemTypeTag, INHERITS_FROM(TwoP, BoxModel, Lens2pSpatialParamsTypeTag));
// Set the grid type
SET_TYPE_PROP(LensTwoPProblem, Grid, Dune::YaspGrid<2>);
SET_TYPE_PROP(LensTwoPProblemTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(LensTwoPProblem, Problem, LensTwoPProblem<TypeTag>);
SET_TYPE_PROP(LensTwoPProblemTypeTag, Problem, LensTwoPProblem<TypeTag>);
// Set the fluid system
SET_PROP(LensTwoPProblem, FluidSystem)
SET_PROP(LensTwoPProblemTypeTag, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::H2O<Scalar> >;
using NonwettingPhase = FluidSystems::OnePGas<Scalar, Components::N2<Scalar> >;
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
}
} // end namespace Properties
/*!
* \ingroup TwoPBoxProblems
......@@ -108,14 +108,10 @@ class LensTwoPProblem : public PorousMediumFlowProblem<TypeTag>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
......@@ -123,18 +119,13 @@ class LensTwoPProblem : public PorousMediumFlowProblem<TypeTag>
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
LensTwoPProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
// this is placed down to the declaration of the private varialbes as it is done in all other cases i have seen eps_ = 3e-6;
// the boundary condition data
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, LowerPressure);
upperPressure_ = getParam<Scalar>("Boundary.UpperPressure"); // GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Boundary, UpperPressure);
lowerPressure_ = getParam<Scalar>("Boundary.LowerPressure");
upperPressure_ = getParam<Scalar>("Boundary.UpperPressure");
}
/*!
......@@ -142,16 +133,6 @@ public:
*/
// \{
/*!
* \brief The problem name.
*
* This is used as a prefix for files generated by the simulation.
*/
std::string name() const
{
return "lens-2p";
}
/*!
* \brief Returns the temperature within the domain.
* \param temperature DOC ME!
......@@ -181,6 +162,7 @@ public:
BoundaryTypes values;
if (onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -209,6 +191,7 @@ public:
values[Indices::pressureIdx] = upperPressure_;
values[Indices::saturationIdx] = 0.0;
}
else if (onLowerBoundary_(globalPos))
{
values[Indices::pressureIdx] = lowerPressure_;
......@@ -286,7 +269,6 @@ public:
// \}
private:
static constexpr Scalar eps_ = 3e-6;
Scalar episodeLength_;
......@@ -320,6 +302,7 @@ private:
&& globalPos[0] < 2.3 - eps_);
}
};
} //end namespace
#endif
} //end namespace Dumux
#endif // DUMUX_LENS2P_PROBLEM_HH
......@@ -26,29 +26,25 @@
#ifndef DUMUX_LENS2P_SPATIALPARAMS_HH
#define DUMUX_LENS2P_SPATIALPARAMS_HH
// #include <dumux/material/spatialparams/implicit.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/spatialparams/fv.hh>
namespace Dumux
{
namespace Dumux {
//forward declaration
template<class FVGridGeometry, class Scalar>
class Lens2pSpatialParams;
namespace Properties
{
namespace Properties {
// The spatial parameters TypeTag
NEW_TYPE_TAG(Lens2pSpatialParams);
NEW_TYPE_TAG(Lens2pSpatialParamsTypeTag);
// Set the spatial parameters
SET_PROP(Lens2pSpatialParams, SpatialParams)
SET_PROP(Lens2pSpatialParamsTypeTag, SpatialParams)
{
private:
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
......@@ -56,25 +52,23 @@ private:
public:
using type = Lens2pSpatialParams<FVGridGeometry, Scalar>;
};
}
} // end namespace Properties
/** \todo Please doc me! */
template<class FVGridGeometry, class Scalar>
class Lens2pSpatialParams
: public FVSpatialParams<FVGridGeometry, Scalar, Lens2pSpatialParams<FVGridGeometry, Scalar>>
: public FVSpatialParams<FVGridGeometry, Scalar, Lens2pSpatialParams<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 ThisType = Lens2pSpatialParams<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:
......@@ -131,6 +125,7 @@ public:
{
if (isInLens_(element.geometry().center()))
return lensK_;
else
return outerK_;
}
......@@ -152,6 +147,7 @@ public:
{
if (isInLens_(element.geometry().center()))
return lensPorosity_;
else
return outerPorosity_;
}
......@@ -172,6 +168,7 @@ public:
// const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
if (isInLense_(element.geometry().center() ))
return std::cbrt(lensPorosity_);
else
return std::cbrt(outerPorosity_);
}
......@@ -192,6 +189,7 @@ public:
{
if (isInLens_(element.geometry().center()))
return lensMaterialParams_;
else
return outerMaterialParams_;
}
......@@ -225,6 +223,7 @@ private:
{
// points on the lens boundary are considered inside
for (int i = 0; i < dimWorld; ++i)
if (globalPos[i] < lensLowerLeft_[i] - eps_ || globalPos[i] > lensUpperRight_[i] + eps_)
return false;
......@@ -243,5 +242,6 @@ private:
Scalar eps_ = 1e-6;
};
} // end namespace
#endif
} // end namespace Dumux
#endif // DUMUX_LENS2P_SPATIALPARAMS_HH
......@@ -5,14 +5,14 @@ add_dumux_test(lens1p2cexercise1 lens1p2cexercise1 lens1p2cexercise1.cc
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/lens-1p2c-exercise1-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lens-1p2c-00101.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/lens1p2cexercise1 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise1.input")
--command "${CMAKE_CURRENT_BINARY_DIR}/lens1p2cexercise1 ${CMAKE_CURRENT_SOURCE_DIR}/exercise1.input -Problem.Name lens-1p2c")
add_dumux_test(lens2pexercise1 lens2pexercise1 lens2pexercise1.cc
python ${dumux_INCLUDE_DIRS}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/lens-2p-exercise1-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/lens-2p-00101.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2pexercise1 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/exercise1.input")
--command "${CMAKE_CURRENT_BINARY_DIR}/lens2pexercise1 ${CMAKE_CURRENT_SOURCE_DIR}/exercise1.input -Problem.Name lens-2p")
# headers for installation and headercheck
install(FILES
......
[TimeManager]
[TimeLoop]
MaxTimeStepSize = 5.0e1 # maximal time step size [s]
TEnd = 5.0e3 # end time of the simulation [s]
DtInitial = 1e1 # initial time step for the simulation [s]
......
......@@ -16,10 +16,35 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#include "config.h"
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.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>
#include "lens1p2cproblem.hh"
#include <dumux/common/start.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
......@@ -60,8 +85,155 @@ void usage(const char *progName, const std::string &errorMsg)
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv)
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(LensOnePTwoCProblemTypeTag);
// 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)
/////////////////////////////////////////////////////////////////////
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> 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 geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(fvGridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
problem->setTimeLoop(timeLoop);
// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
using LinearSolver = AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the Newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
catch (Dumux::ParameterException &e)