Commit 177c9fce authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

finished transport test

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4062 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent b44b0f10
......@@ -322,8 +322,7 @@ public:
if (problem_->shouldWriteOutput())
problem_->writeOutput();
// advance the simulated time by the current time step
// size
// advance the simulated time by the current time step size
Scalar dt = timeStepSize();
time_ += timeStepSize_;
++timeStepIdx_;
......
......@@ -31,6 +31,18 @@ class DiffusivePart;
template<class TypeTag>
class ConvectivePart;
template<class TypeTag>
class TwoPCommonIndices;
template<class TypeTag>
class FluidSystem2P;
template<class TypeTag>
class TwoPFluidState;
template<class TypeTag>
class VariableClass2P;
namespace Properties
{
/*!
......@@ -51,10 +63,31 @@ NEW_TYPE_TAG(Transport);
NEW_PROP_TAG( DiffusivePart ); //!< The type of the diffusive part in a transport equation
NEW_PROP_TAG( ConvectivePart ); //!< The type of a convective part in a transport equation
NEW_PROP_TAG( Variables );
NEW_PROP_TAG( NumPhases );
NEW_PROP_TAG( NumComponents );
NEW_PROP_TAG( TwoPIndices );
NEW_PROP_TAG( FluidSystem );
NEW_PROP_TAG( FluidState );
NEW_PROP_TAG( EnableCompressibility );
NEW_PROP_TAG( PressureFormulation );
NEW_PROP_TAG( SaturationFormulation );
NEW_PROP_TAG( VelocityFormulation );
NEW_PROP_TAG( CFLFactor );
SET_TYPE_PROP(Transport, DiffusivePart, DiffusivePart<TypeTag>);
SET_TYPE_PROP(Transport, ConvectivePart, ConvectivePart<TypeTag>);
SET_TYPE_PROP(Transport, Variables, VariableClass2P<TypeTag>);
SET_INT_PROP(Transport, NumPhases, 2);
SET_INT_PROP(Transport, NumComponents, 1);
SET_TYPE_PROP(Transport, TwoPIndices, TwoPCommonIndices<TypeTag>);
SET_TYPE_PROP(Transport, FluidSystem, FluidSystem2P<TypeTag>);
SET_TYPE_PROP(Transport, FluidState, TwoPFluidState<TypeTag>);
SET_BOOL_PROP(Transport, EnableCompressibility, false);
SET_INT_PROP(Transport, PressureFormulation, TwoPCommonIndices<TypeTag>::pressureW);
SET_INT_PROP(Transport, SaturationFormulation, TwoPCommonIndices<TypeTag>::saturationW);
SET_INT_PROP(Transport, VelocityFormulation, TwoPCommonIndices<TypeTag>::velocityTotal);
SET_SCALAR_PROP(Transport, CFLFactor, 1.0);
}
}
......
......@@ -82,7 +82,7 @@ public:
typedef typename SolutionTypes::PhasePropertyElemFace PhasePropertyElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of scalars
typedef typename SolutionTypes::DimVecElemFace DimVecElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of vector (of size dimension) of scalars
private:
public:
const int codim_;
ScalarSolutionType saturation_;
......
......@@ -33,6 +33,11 @@
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(CFLFactor);
}
/*! \ingroup diffusion
* @brief base class that defines the parameters of loosely coupled diffusion and transport equations
*
......@@ -67,6 +72,9 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
typedef typename SolutionTypes::ScalarSolution Solution;
enum
{
dim = GridView::dimension,
......@@ -108,6 +116,7 @@ public:
// bboxMax_[i] = std::max(bboxMax_[i], vIt->geometry().corner(0)[i]);
// }
// }
cFLFactor_ = GET_PROP_VALUE(TypeTag, PTAG(CFLFactor));
model_ = new Model(asImp_()) ;
}
......@@ -126,7 +135,7 @@ public:
void init()
{
// set the initial condition of the model
model().initial();
model().initialize();
}
/*!
......@@ -150,7 +159,6 @@ public:
void timeIntegration()
{
// allocate temporary vectors for the updates
typedef typename Model::SolutionType Solution;
Solution k1 = asImp_().variables().saturation();
dt_ = 1e100;
......@@ -160,7 +168,7 @@ public:
model().update(t, dt_, k1);
//make sure t_old + dt is not larger than tend
dt_ = std::min(dt_, timeManager().episodeMaxTimeStepSize());
dt_ = std::min(dt_*cFLFactor_, timeManager().episodeMaxTimeStepSize());
timeManager().setTimeStepSize(dt_);
// explicit Euler: Sat <- Sat + dt*N(Sat)
......@@ -226,7 +234,7 @@ public:
if (gridView().comm().rank() == 0)
std::cout << "Writing result file for current time step\n";
resultWriter_.beginTimestep(timeManager_.time(), gridView());
resultWriter_.beginTimestep(timeManager_.time() + timeManager_.timeStepSize(), gridView());
model().addOutputVtkFields(resultWriter_);
resultWriter_.endTimestep();
}
......@@ -338,7 +346,7 @@ public:
typedef Dumux::Restart Restarter;
Restarter res;
res.serializeBegin(*asImp_());
res.serializeBegin(asImp_());
std::cerr << "Serialize to file " << res.fileName() << "\n";
timeManager_.serialize(res);
......@@ -405,6 +413,8 @@ private:
Scalar dt_;
Scalar cFLFactor_;
Model* model_;
VtkMultiWriter resultWriter_;
......
......@@ -141,19 +141,6 @@ public:
instream >> pressure_[globalIdx];
}
private:
void initializeGlobalVariables(Dune::FieldVector<Scalar, dim>& initialVel)
{
//resize to grid size
pressure_.resize(gridSize_);
velocity_.resize(gridSize_);//depends on pressure
potential_.resize(gridSize_);//depends on pressure
//initialise variables
pressure_ = 0;
velocity_ = initialVel;
}
void initializePotentials(Dune::FieldVector<Scalar, dim>& initialVel)
{
if (initialVel.two_norm())
......@@ -185,6 +172,50 @@ private:
return;
}
private:
void initializeGlobalVariables(Dune::FieldVector<Scalar, dim>& initialVel)
{
//resize to grid size
pressure_.resize(gridSize_);
velocity_.resize(gridSize_);//depends on pressure
potential_.resize(gridSize_);//depends on pressure
//initialise variables
pressure_ = 0;
velocity_ = initialVel;
}
// void initializePotentials(Dune::FieldVector<Scalar, dim>& initialVel)
// {
// if (initialVel.two_norm())
// {
// // compute update vector
// ElementIterator eItEnd = gridView_.template end<0> ();
// for (ElementIterator eIt = gridView_.template begin<0> (); eIt != eItEnd; ++eIt)
// {
// // cell index
// int globalIdxI = elementMapper_.map(*eIt);
//
// // run through all intersections with neighbors and boundary
// IntersectionIterator isItEnd = gridView_.iend(*eIt);
// for (IntersectionIterator isIt = gridView_.ibegin(*eIt); isIt != isItEnd; ++isIt)
// {
// // local number of facet
// int indexInInside = isIt->indexInInside();
//
// Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->centerUnitOuterNormal();
//
// for (int i = 0; i < numPhase; i++) {potential_[globalIdxI][indexInInside][i] = initialVel * unitOuterNormal;}
// }
// }
// }
// else
// {
// potential_ = Dune::FieldVector<Scalar, numPhase> (0);
// }
// return;
// }
//Write saturation and pressure into file
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
......
DGF
Interval
0 0 % first corner
1 1 % second corner
10 10 % cells in x and y direction
#
BOUNDARYDOMAIN
default 1 % all boundaries have id 1
#BOUNDARYDOMAIN
# unitcube.dgf
......@@ -33,7 +33,7 @@
////////////////////////
void usage(const char *progname)
{
std::cout << boost::format("usage: %s [--restart restartTime] tEnd\n")%progname;
std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd\n")%progname;
exit(1);
}
......@@ -54,7 +54,7 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////
// parse the command line arguments
////////////////////////////////////////////////////////////
if (argc < 2)
if (argc < 3)
usage(argv[0]);
// deal with the restart stuff
......@@ -68,28 +68,27 @@ int main(int argc, char** argv)
std::istringstream(argv[argPos++]) >> restartTime;
}
if (argc - argPos != 1) {
if (argc - argPos != 2) {
usage(argv[0]);
}
////////////////////////////////////////////////////////////
// create the grid
////////////////////////////////////////////////////////////
const char *dgfFileName = argv[argPos++];
Dune::GridPtr<Grid> gridPtr(dgfFileName);
// read the initial time step and the end time
double tEnd, dt;
std::istringstream(argv[argPos++]) >> tEnd;
dt = tEnd;
////////////////////////////////////////////////////////////
// create the grid
// instantiate and run the concrete problem
////////////////////////////////////////////////////////////
Dune::FieldVector<int,dim> N(4); N[0] = 4;
Dune::FieldVector<double ,dim> L(0);
Dune::FieldVector<double,dim> H(1);
Grid grid(N,L,H);
////////////////////////////////////////////////////////////
// instantiate and run the concrete problem
////////////////////////////////////////////////////////////
Problem problem(grid.leafView(), L, H);
Problem problem(gridPtr->leafView(), L, H);
// load restart file if necessarry
if (restart)
......
......@@ -22,8 +22,7 @@
#include <dune/grid/uggrid.hh>
#endif
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/sgrid.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include <dumux/material/components/unit.hh>
......@@ -51,8 +50,7 @@ NEW_TYPE_TAG(TransportTestProblem, INHERITS_FROM(DecoupledModel, Transport));
// Set the grid type
SET_PROP(TransportTestProblem, Grid)
{
// typedef Dune::YaspGrid<2> type;
typedef Dune::SGrid<2, 2> type;
typedef Dune::UGGrid<2> type;
};
// Set the problem property
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment