Commit 1b0b176a authored by Andreas Lauser's avatar Andreas Lauser
Browse files

boxmodels: fix a few issues found while porting the benchmark3 problem

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@3967 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent bf7ffe18
......@@ -90,7 +90,6 @@ class TwoPModel : public BoxModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSecondaryVars)) ElementSecondaryVars;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
......@@ -108,13 +107,15 @@ public:
* and get minimum and maximum values of primary variables
*
*/
void calculateMass(const SolutionVector &sol, Dune::FieldVector<Scalar, 2> &mass)
void calculateMass(Dune::FieldVector<Scalar, 2> &mass)
{
const SolutionVector &sol = this->curSol();
mass = 0;
ElementIterator elementIt =
this->model().gridView().template begin<0> ();
ElementIterator endit = this->model().gridView().template end<0> ();
this->gridView_().template begin<0> ();
ElementIterator endit = this->gridView_().template end<0> ();
ElementSecondaryVars elemDat;
......@@ -129,15 +130,12 @@ public:
FVElementGeometry fvElemGeom;
SecondaryVars secVars;
ElementBoundaryTypes elemBcTypes;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt)
{
int idx = this->problem_.model().elementMapper().map(*elemIt);
fvElemGeom.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
int numLocalVerts = elementIt->template count<dim> ();
for (int i = 0; i < numLocalVerts; ++i)
......@@ -174,16 +172,16 @@ public:
}
}
mass = this->gridView_.comm().sum(mass);
mass = this->gridView_().comm().sum(mass);
Scalar minS = this->gridView_.comm().min(minSat);
Scalar maxS = this->gridView_.comm().max(maxSat);
Scalar minPr = this->gridView_.comm().min(minP);
Scalar maxPr = this->gridView_.comm().max(maxP);
Scalar minT = this->gridView_.comm().min(minTe);
Scalar maxT = this->gridView_.comm().max(maxTe);
Scalar minS = this->gridView_().comm().min(minSat);
Scalar maxS = this->gridView_().comm().max(maxSat);
Scalar minPr = this->gridView_().comm().min(minP);
Scalar maxPr = this->gridView_().comm().max(maxP);
Scalar minT = this->gridView_().comm().min(minTe);
Scalar maxT = this->gridView_().comm().max(maxTe);
if (this->problem_.gridView().comm().rank() == 0) // IF PARALLEL: only print by processor with rank() == 0
if (this->gridView_().comm().rank() == 0) // IF PARALLEL: only print by processor with rank() == 0
{
// print minimum and maximum values
std::cout << "nonwetting phase saturation: min = " << minS
......@@ -227,7 +225,6 @@ public:
FVElementGeometry fvElemGeom;
SecondaryVars secVars;
ElementBoundaryTypes elemBcTypes;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
......@@ -237,7 +234,6 @@ public:
(*rank)[idx] = this->gridView_().comm().rank();
fvElemGeom.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
int numVerts = elemIt->template count<dim> ();
for (int i = 0; i < numVerts; ++i)
......@@ -275,24 +271,6 @@ public:
writer.addVertexData(Te, "temperature");
writer.addCellData(rank, "process rank");
}
/*!
* \brief Calculate the flux of the nonwetting phase across a given
* layer for the current timestep
*/
void calculateFluxAcrossLayer(Dune::FieldVector<Scalar, 2> &flux, int coord, Scalar coordVal)
{
this->localJacobian().calculateFluxAcrossLayer(this->curSol(), flux, coord, coordVal);
}
/*!
* \brief Calculate the phase masses in the system for the current
* timestep.
*/
void calculateMass(Dune::FieldVector<Scalar, 2> &mass)
{
this->localJacobian().calculateMass(this->curSol(), mass);
}
};
}
......
......@@ -226,19 +226,19 @@ public:
* and get minimum and maximum values of primary variables
*
*/
void calculateMass(const SolutionVector &sol, Dune::FieldVector<
Scalar, 2> &massGas, Dune::FieldVector<Scalar, 2> &massLiquid)
void calculateMass(Dune::FieldVector<Scalar, 2> &massGas,
Dune::FieldVector<Scalar, 2> &massLiquid)
{
const SolutionVector &sol = this->curSol();
massGas = 0;
massLiquid = 0;
ElementIterator elemIt =
this->gridView().template begin<0> ();
ElementIterator endit = this->gridView().template end<0> ();
ElementIterator elemIt = this->gridView_().template begin<0> ();
ElementIterator endit = this->gridView_().template end<0> ();
FVElementGeometry fvElemGeom;
SecondaryVars secVars;
ElementBoundaryTypes elemBcTypes;
Scalar minSat = 1e100;
Scalar maxSat = -1e100;
......@@ -255,9 +255,7 @@ public:
if (elemIt->partitionType() != Dune::InteriorEntity)
continue;
fvElemGeom.update(this->gridView_()(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
fvElemGeom.update(this->gridView_(), *elemIt);
// Loop over element vertices
for (int i = 0; i < fvElemGeom.numVertices; ++i)
{
......@@ -270,7 +268,7 @@ public:
false);
const FluidState &fs = secVars.fluidState();
Scalar vol = this->fvElemGeom_().subContVol[i].volume;
Scalar vol = fvElemGeom.subContVol[i].volume;
Scalar satN = fs.saturation(gPhaseIdx);
Scalar xAW = fs.massFrac(lPhaseIdx, gCompIdx);
......@@ -396,7 +394,6 @@ public:
FVElementGeometry fvElemGeom;
SecondaryVars secVars;
ElementBoundaryTypes elemBcTypes;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
......@@ -404,9 +401,7 @@ public:
{
int idx = this->problem_().elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
fvElemGeom.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
int numVerts = elemIt->template count<dim> ();
for (int i = 0; i < numVerts; ++i)
......
......@@ -79,7 +79,6 @@ namespace Dumux {
* specified by setting the <tt>Formulation</tt> property to either <tt>TwoPNIIndices::pWsN</tt> or <tt>TwoPIndices::pNsW</tt>. By
* default, the model uses \f$p_w\f$, \f$S_n\f$ and \f$T\f$.
*/
template<class TypeTag>
class TwoPNIModel: public TwoPModel<TypeTag>
{
......
......@@ -75,8 +75,6 @@ public:
timeManager_(&timeManager),
resultWriter_(asImp_().name())
{
wasRestarted_ = false;
// calculate the bounding box of the grid view
VertexIterator vIt = gridView.template begin<dim>();
const VertexIterator vEndIt = gridView.template end<dim>();
......@@ -181,8 +179,7 @@ public:
*/
bool doSerialize() const
{
return !restarted() &&
timeManager().timeStepIndex() > 0 &&
return timeManager().timeStepIndex() > 0 &&
(timeManager().timeStepIndex() % 10 == 0);
}
......@@ -203,6 +200,15 @@ public:
void postProcess()
{ }
/*!
* \brief Called when the end of an simulation episode is reached.
*/
void episodeEnd()
{
std::cerr << "The end of an episode is reached, but the problem "
<< "does not override the episodeEnd() method. "
<< "Doing nothing!\n";
};
// \}
/*!
......@@ -291,13 +297,6 @@ public:
*/
// \{
/*!
* \brief Returns true, if the current state of the problem was
* loaded from a restart file.
*/
bool restarted() const
{ return wasRestarted_; }
/*!
* \brief This method writes the complete state of the simulation
* to the harddisk.
......@@ -363,8 +362,6 @@ public:
{
resultWriter_.deserialize(res);
model().deserialize(res);
wasRestarted_ = true;
};
// \}
......@@ -417,8 +414,6 @@ private:
NewtonController newtonCtl_;
VtkMultiWriter resultWriter_;
bool wasRestarted_;
};
// definition of the static class member simname_,
// which is necessary because it is of type string.
......
......@@ -20,6 +20,8 @@
#ifndef DUMUX_TIME_MANAGER_HH
#define DUMUX_TIME_MANAGER_HH
#include <dumux/common/propertysystem.hh>
#include <boost/format.hpp>
#include <dune/common/timer.hh>
......@@ -29,6 +31,11 @@
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Problem);
}
/*!
* \addtogroup SimControl Simulation Supervision
*/
......@@ -241,13 +248,21 @@ public:
* simulated time.
*/
Scalar episodeLength() const
{ return std::min(episodeLength_, endTime_ - episodeStartTime_); }
{ return episodeLength_; }
/*!
* \brief Returns true if the current episode is over.
* \brief Returns true if the current episode is finished at the
* current time.
*/
bool episodeIsOver() const
{ return time() >= episodeStartTime_ + (1 - 1e-14)*episodeLength(); }
{ return time() >= episodeStartTime_ + episodeLength(); }
/*!
* \brief Returns true if the current episode will be finished
* after the current time step.
*/
bool episodeWillBeOver() const
{ return time() + timeStepSize() >= episodeStartTime_ + episodeLength(); }
/*!
......@@ -308,6 +323,9 @@ public:
if (problem_->doSerialize())
problem_->serialize();
if (episodeIsOver())
problem_->episodeEnd();
// notify the problem that the timestep is done and ask it
// for a suggestion for the next timestep size
......
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