From 49a261a8f61738bb50a5d0bbae819514d2848bae Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 5 Dec 2012 13:47:51 +0000
Subject: [PATCH] implicit branch: add cell-centered treatment to the 3p3c
model, unify the test.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9783 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/2p2c/2p2cmodel.hh | 4 +-
dumux/implicit/2p2c/2p2cvolumevariables.hh | 8 +-
dumux/implicit/3p3c/3p3cmodel.hh | 316 +++++++++++-------
dumux/implicit/3p3c/3p3cvolumevariables.hh | 11 +-
test/implicit/3p3c/Makefile.am | 5 +-
test/implicit/3p3c/infiltrationproblem.hh | 48 +--
test/implicit/3p3c/test_box3p3c.cc | 59 ++++
test/implicit/3p3c/test_box3p3c.input | 31 ++
.../3p3c/{test_3p3c.cc => test_cc3p3c.cc} | 2 +-
.../{test_3p3c.input => test_cc3p3c.input} | 3 +
10 files changed, 331 insertions(+), 156 deletions(-)
create mode 100644 test/implicit/3p3c/test_box3p3c.cc
create mode 100644 test/implicit/3p3c/test_box3p3c.input
rename test/implicit/3p3c/{test_3p3c.cc => test_cc3p3c.cc} (98%)
rename test/implicit/3p3c/{test_3p3c.input => test_cc3p3c.input} (96%)
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index 82a91f9e16..19e9475954 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -266,9 +266,9 @@ public:
}
/*!
- * \brief Returns the phase presence of the current or the old solution of a vertex.
+ * \brief Returns the phase presence of the current or the old solution of a degree of freedom.
*
- * \param globalIdx The global vertex index
+ * \param globalIdx The global index of the degree of freedom
* \param oldSol Evaluate function with solution of current or previous time step
*/
int phasePresence(int globalIdx, bool oldSol) const
diff --git a/dumux/implicit/2p2c/2p2cvolumevariables.hh b/dumux/implicit/2p2c/2p2cvolumevariables.hh
index ac9719c58b..9591cae3f7 100644
--- a/dumux/implicit/2p2c/2p2cvolumevariables.hh
+++ b/dumux/implicit/2p2c/2p2cvolumevariables.hh
@@ -101,6 +101,8 @@ class TwoPTwoCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem> MiscibleMultiPhaseComposition;
typedef Dumux::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
//! The type of the object returned by the fluidState() method
typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
@@ -186,10 +188,10 @@ public:
unsigned numVertices = problem.gridView().size(dim);
int globalIdx;
- if (numDofs != numVertices) // element data
- globalIdx = problem.model().dofMapper().map(element);
- else
+ if (isBox) // vertex data
globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ else
+ globalIdx = problem.model().dofMapper().map(element);
int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
diff --git a/dumux/implicit/3p3c/3p3cmodel.hh b/dumux/implicit/3p3c/3p3cmodel.hh
index c12cadbb93..bc8da5efe6 100644
--- a/dumux/implicit/3p3c/3p3cmodel.hh
+++ b/dumux/implicit/3p3c/3p3cmodel.hh
@@ -139,6 +139,8 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief Initialize the static data with the initial solution.
@@ -149,25 +151,47 @@ public:
{
ParentType::init(problem);
- staticVertexDat_.resize(this->gridView_().size(dim));
+ staticDat_.resize(this->numDofs());
setSwitched_(false);
- VertexIterator it = this->gridView_().template begin<dim> ();
- const VertexIterator &endit = this->gridView_().template end<dim> ();
- for (; it != endit; ++it)
+ if (isBox)
{
- int globalIdx = this->dofMapper().map(*it);
- const GlobalPosition &globalPos = it->geometry().corner(0);
+ VertexIterator it = this->gridView_().template begin<dim> ();
+ const VertexIterator &endit = this->gridView_().template end<dim> ();
+ for (; it != endit; ++it)
+ {
+ int globalIdx = this->dofMapper().map(*it);
+ const GlobalPosition &globalPos = it->geometry().corner(0);
- // initialize phase presence
- staticVertexDat_[globalIdx].phasePresence
- = this->problem_().initialPhasePresence(*it, globalIdx,
+ // initialize phase presence
+ staticDat_[globalIdx].phasePresence
+ = this->problem_().initialPhasePresence(*it, globalIdx,
globalPos);
- staticVertexDat_[globalIdx].wasSwitched = false;
+ staticDat_[globalIdx].wasSwitched = false;
- staticVertexDat_[globalIdx].oldPhasePresence
- = staticVertexDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
+ }
+ }
+ else
+ {
+ ElementIterator it = this->gridView_().template begin<0> ();
+ const ElementIterator &endit = this->gridView_().template end<0> ();
+ for (; it != endit; ++it)
+ {
+ int globalIdx = this->dofMapper().map(*it);
+ const GlobalPosition &globalPos = it->geometry().center();
+
+ // initialize phase presence
+ staticDat_[globalIdx].phasePresence
+ = this->problem_().initialPhasePresence(*this->gridView_().template begin<dim> (),
+ globalIdx, globalPos);
+ staticDat_[globalIdx].wasSwitched = false;
+
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
+ }
}
}
@@ -175,26 +199,27 @@ public:
* \brief Compute the total storage inside one phase of all
* conservation quantities.
*
- * \param dest Contains the storage of each component for one phase
+ * \param storage Contains the storage of each component for one phase
* \param phaseIdx The phase index
*/
- void globalPhaseStorage(PrimaryVariables &dest, int phaseIdx)
+ void globalPhaseStorage(PrimaryVariables &storage, const int phaseIdx)
{
- dest = 0;
+ storage = 0;
- ElementIterator elemIt = this->gridView_().template begin<0>();
- const ElementIterator elemEndIt = this->gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- this->localResidual().evalPhaseStorage(*elemIt, phaseIdx);
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ const ElementIterator eEndIt = this->gridView_().template end<0>();
+ for (; eIt != eEndIt; ++eIt) {
+ this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
- for (int i = 0; i < elemIt->template count<dim>(); ++i)
- dest += this->localResidual().residual(i);
+ for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+ storage += this->localResidual().storageTerm()[i];
}
if (this->gridView_().comm().size() > 1)
- dest = this->gridView_().comm().sum(dest);
+ storage = this->gridView_().comm().sum(storage);
}
+
/*!
* \brief Called by the update() method if applying the newton
* method was unsuccessful.
@@ -233,17 +258,17 @@ public:
}
/*!
- * \brief Returns the phase presence of the current or the old solution of a vertex.
+ * \brief Returns the phase presence of the current or the old solution of a degree of freedom.
*
- * \param globalVertexIdx The global vertex index
+ * \param globalIdx The global index of the degree of freedom
* \param oldSol Evaluate function with solution of current or previous time step
*/
- int phasePresence(int globalVertexIdx, bool oldSol) const
+ int phasePresence(int globalIdx, bool oldSol) const
{
return
oldSol
- ? staticVertexDat_[globalVertexIdx].oldPhasePresence
- : staticVertexDat_[globalVertexIdx].phasePresence;
+ ? staticDat_[globalIdx].oldPhasePresence
+ : staticDat_[globalIdx].phasePresence;
}
/*!
@@ -260,27 +285,28 @@ public:
{
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
- // create the required scalar fields
- unsigned numVertices = this->problem_().gridView().size(dim);
+ // get the number of degrees of freedom
+ unsigned numDofs = this->numDofs();
+ // create the required scalar fields
ScalarField *saturation[numPhases];
ScalarField *pressure[numPhases];
ScalarField *density[numPhases];
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
- saturation[phaseIdx] = writer.allocateManagedBuffer(numVertices);
- pressure[phaseIdx] = writer.allocateManagedBuffer(numVertices);
- density[phaseIdx] = writer.allocateManagedBuffer(numVertices);
+ saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
}
- ScalarField *phasePresence = writer.allocateManagedBuffer (numVertices);
+ ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
ScalarField *moleFraction[numPhases][numComponents];
for (int i = 0; i < numPhases; ++i)
for (int j = 0; j < numComponents; ++j)
- moleFraction[i][j] = writer.allocateManagedBuffer (numVertices);
- ScalarField *temperature = writer.allocateManagedBuffer (numVertices);
- ScalarField *poro = writer.allocateManagedBuffer(numVertices);
- ScalarField *perm = writer.allocateManagedBuffer(numVertices);
+ moleFraction[i][j] = writer.allocateManagedBuffer (numDofs);
+ ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer(numDofs);
+ ScalarField *perm = writer.allocateManagedBuffer(numDofs);
unsigned numElements = this->gridView_().size(0);
ScalarField *rank =
@@ -297,15 +323,19 @@ public:
(*rank)[idx] = this->gridView_().comm().rank();
fvGeometry.update(this->gridView_(), *elemIt);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx;
+ if (numDofs == numElements) // element data
+ globalIdx = idx;
+ else
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
fvGeometry,
- i,
+ scvIdx,
false);
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
@@ -327,79 +357,113 @@ public:
(*poro)[globalIdx] = volVars.porosity();
(*perm)[globalIdx] = volVars.permeability();
(*temperature)[globalIdx] = volVars.temperature();
- (*phasePresence)[globalIdx] = staticVertexDat_[globalIdx].phasePresence;
+ (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
}
}
- writer.attachVertexData(*saturation[wPhaseIdx], "Sw");
- writer.attachVertexData(*saturation[nPhaseIdx], "Sn");
- writer.attachVertexData(*saturation[gPhaseIdx], "Sg");
- writer.attachVertexData(*pressure[wPhaseIdx], "pw");
- writer.attachVertexData(*pressure[nPhaseIdx], "pn");
- writer.attachVertexData(*pressure[gPhaseIdx], "pg");
- writer.attachVertexData(*density[wPhaseIdx], "rhow");
- writer.attachVertexData(*density[nPhaseIdx], "rhon");
- writer.attachVertexData(*density[gPhaseIdx], "rhog");
-
- for (int i = 0; i < numPhases; ++i)
+ if (isBox) // vertex data
{
- for (int j = 0; j < numComponents; ++j)
- {
- std::ostringstream oss;
- oss << "x^"
- << FluidSystem::phaseName(i)
- << "_"
- << FluidSystem::componentName(j);
- writer.attachVertexData(*moleFraction[i][j], oss.str().c_str());
+ writer.attachVertexData(*saturation[wPhaseIdx], "Sw");
+ writer.attachVertexData(*saturation[nPhaseIdx], "Sn");
+ writer.attachVertexData(*saturation[gPhaseIdx], "Sg");
+ writer.attachVertexData(*pressure[wPhaseIdx], "pw");
+ writer.attachVertexData(*pressure[nPhaseIdx], "pn");
+ writer.attachVertexData(*pressure[gPhaseIdx], "pg");
+ writer.attachVertexData(*density[wPhaseIdx], "rhow");
+ writer.attachVertexData(*density[nPhaseIdx], "rhon");
+ writer.attachVertexData(*density[gPhaseIdx], "rhog");
+
+ for (int i = 0; i < numPhases; ++i)
+ {
+ for (int j = 0; j < numComponents; ++j)
+ {
+ std::ostringstream oss;
+ oss << "x^"
+ << FluidSystem::phaseName(i)
+ << "_"
+ << FluidSystem::componentName(j);
+ writer.attachVertexData(*moleFraction[i][j], oss.str().c_str());
+ }
}
+ writer.attachVertexData(*poro, "porosity");
+ writer.attachVertexData(*perm, "permeability");
+ writer.attachVertexData(*temperature, "temperature");
+ writer.attachVertexData(*phasePresence, "phase presence");
}
- writer.attachVertexData(*poro, "porosity");
- writer.attachVertexData(*perm, "permeability");
- writer.attachVertexData(*temperature, "temperature");
- writer.attachVertexData(*phasePresence, "phase presence");
+ else // cell data
+ {
+ writer.attachCellData(*saturation[wPhaseIdx], "Sw");
+ writer.attachCellData(*saturation[nPhaseIdx], "Sn");
+ writer.attachCellData(*saturation[gPhaseIdx], "Sg");
+ writer.attachCellData(*pressure[wPhaseIdx], "pw");
+ writer.attachCellData(*pressure[nPhaseIdx], "pn");
+ writer.attachCellData(*pressure[gPhaseIdx], "pg");
+ writer.attachCellData(*density[wPhaseIdx], "rhow");
+ writer.attachCellData(*density[nPhaseIdx], "rhon");
+ writer.attachCellData(*density[gPhaseIdx], "rhog");
+
+ for (int i = 0; i < numPhases; ++i)
+ {
+ for (int j = 0; j < numComponents; ++j)
+ {
+ std::ostringstream oss;
+ oss << "x^"
+ << FluidSystem::phaseName(i)
+ << "_"
+ << FluidSystem::componentName(j);
+ writer.attachCellData(*moleFraction[i][j], oss.str().c_str());
+ }
+ }
+ writer.attachCellData(*poro, "porosity");
+ writer.attachCellData(*perm, "permeability");
+ writer.attachCellData(*temperature, "temperature");
+ writer.attachCellData(*phasePresence, "phase presence");
+ }
+
writer.attachCellData(*rank, "process rank");
}
/*!
* \brief Write the current solution to a restart file.
*
- * \param outStream The output stream of one vertex for the restart file
- * \param vert The vertex
+ * \param outStream The output stream of one entity for the restart file
+ * \param entity The entity, either a vertex or an element
*/
- void serializeEntity(std::ostream &outStream, const Vertex &vert)
+ template<class Entity>
+ void serializeEntity(std::ostream &outStream, const Entity &entity)
{
// write primary variables
- ParentType::serializeEntity(outStream, vert);
+ ParentType::serializeEntity(outStream, entity);
- int vertIdx = this->dofMapper().map(vert);
+ int globalIdx = this->dofMapper().map(entity);
if (!outStream.good())
- DUNE_THROW(Dune::IOError, "Could not serialize vertex " << vertIdx);
+ DUNE_THROW(Dune::IOError, "Could not serialize entity " << globalIdx);
- outStream << staticVertexDat_[vertIdx].phasePresence << " ";
+ outStream << staticDat_[globalIdx].phasePresence << " ";
}
/*!
- * \brief Reads the current solution for a vertex from a restart
- * file.
+ * \brief Reads the current solution from a restart file.
*
- * \param inStream The input stream of one vertex from the restart file
- * \param vert The vertex
+ * \param inStream The input stream of one entity from the restart file
+ * \param entity The entity, either a vertex or an element
*/
- void deserializeEntity(std::istream &inStream, const Vertex &vert)
+ template<class Entity>
+ void deserializeEntity(std::istream &inStream, const Entity &entity)
{
// read primary variables
- ParentType::deserializeEntity(inStream, vert);
+ ParentType::deserializeEntity(inStream, entity);
// read phase presence
- int vertIdx = this->dofMapper().map(vert);
+ int globalIdx = this->dofMapper().map(entity);
if (!inStream.good())
DUNE_THROW(Dune::IOError,
- "Could not deserialize vertex " << vertIdx);
+ "Could not deserialize entity " << globalIdx);
- inStream >> staticVertexDat_[vertIdx].phasePresence;
- staticVertexDat_[vertIdx].oldPhasePresence
- = staticVertexDat_[vertIdx].phasePresence;
+ inStream >> staticDat_[globalIdx].phasePresence;
+ staticDat_[globalIdx].oldPhasePresence
+ = staticDat_[globalIdx].phasePresence;
}
@@ -414,36 +478,44 @@ public:
{
bool wasSwitched = false;
- for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
- staticVertexDat_[i].visited = false;
+ for (unsigned i = 0; i < staticDat_.size(); ++i)
+ staticDat_[i].visited = false;
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
- ElementIterator it = this->gridView_().template begin<0> ();
- const ElementIterator &endit = this->gridView_().template end<0> ();
- for (; it != endit; ++it)
+ ElementIterator eIt = this->gridView_().template begin<0> ();
+ const ElementIterator &eEndIt = this->gridView_().template end<0> ();
+ for (; eIt != eEndIt; ++eIt)
{
- fvGeometry.update(this->gridView_(), *it);
- for (int i = 0; i < fvGeometry.numVertices; ++i)
+ fvGeometry.update(this->gridView_(), *eIt);
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*it, i, dim);
+ int globalIdx;
+
+ if (isBox)
+ globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*eIt);
- if (staticVertexDat_[globalIdx].visited)
+ if (staticDat_[globalIdx].visited)
continue;
- staticVertexDat_[globalIdx].visited = true;
+ staticDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
- *it,
+ *eIt,
fvGeometry,
- i,
+ scvIdx,
false);
- const GlobalPosition &global = it->geometry().corner(i);
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
- global))
+ globalPos))
+ {
+ this->jacobianAssembler().markDofRed(globalIdx);
wasSwitched = true;
+ }
}
}
@@ -477,12 +549,11 @@ protected:
*/
void resetPhasePresence_()
{
- int numVertices = this->gridView_().size(dim);
- for (int i = 0; i < numVertices; ++i)
+ for (int i = 0; i < this->numDofs(); ++i)
{
- staticVertexDat_[i].phasePresence
- = staticVertexDat_[i].oldPhasePresence;
- staticVertexDat_[i].wasSwitched = false;
+ staticDat_[i].phasePresence
+ = staticDat_[i].oldPhasePresence;
+ staticDat_[i].wasSwitched = false;
}
}
@@ -491,12 +562,11 @@ protected:
*/
void updateOldPhasePresence_()
{
- int numVertices = this->gridView_().size(dim);
- for (int i = 0; i < numVertices; ++i)
+ for (int i = 0; i < this->numDofs(); ++i)
{
- staticVertexDat_[i].oldPhasePresence
- = staticVertexDat_[i].phasePresence;
- staticVertexDat_[i].wasSwitched = false;
+ staticDat_[i].oldPhasePresence
+ = staticDat_[i].phasePresence;
+ staticDat_[i].wasSwitched = false;
}
}
@@ -518,14 +588,14 @@ protected:
{
// evaluate primary variable switch
bool wouldSwitch = false;
- int phasePresence = staticVertexDat_[globalIdx].phasePresence;
+ int phasePresence = staticDat_[globalIdx].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
if (phasePresence == threePhases)
{
Scalar Smin = 0;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(gPhaseIdx) <= Smin)
@@ -583,7 +653,7 @@ protected:
Scalar xgMax = 1.0;
if (xwg + xgg + xng > xgMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xgMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -609,7 +679,7 @@ protected:
Scalar xnMax = 1.0;
if (xnn > xnMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xnMax *= 1.02;
// if the sum of the hypothetical mole fractions would be larger than
@@ -649,7 +719,7 @@ protected:
bool wettingFlag = 0;
Scalar Smin = 0.0;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
@@ -672,7 +742,7 @@ protected:
Scalar xwMax = 1.0;
if (xww > xwMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xwMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -714,7 +784,7 @@ protected:
bool gasFlag = 0;
Scalar Smin = 0.0;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
@@ -740,7 +810,7 @@ protected:
Scalar xgMax = 1.0;
if (xwg + xgg + xng > xgMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xgMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -791,7 +861,7 @@ protected:
Scalar xnMax = 1.0;
if (xnn > xnMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xnMax *= 1.02;
// if the sum of the hypothetical mole fraction would be larger than
@@ -813,7 +883,7 @@ protected:
Scalar xwMax = 1.0;
if (xww > xwMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xwMax *= 1.02;
// if the sum of the mole fractions would be larger than
@@ -862,7 +932,7 @@ protected:
Scalar xnMax = 1.0;
if (xnn > xnMax)
wouldSwitch = true;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
xnMax *= 1.02;
// if the sum of the hypothetical mole fraction would be larger than
@@ -877,7 +947,7 @@ protected:
}
Scalar Smin = -1.e-6;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(gPhaseIdx) <= Smin)
@@ -891,7 +961,7 @@ protected:
}
Smin = 0.0;
- if (staticVertexDat_[globalIdx].wasSwitched)
+ if (staticDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(wPhaseIdx) <= Smin)
@@ -935,13 +1005,13 @@ protected:
}
}
- staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
- staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
+ staticDat_[globalIdx].phasePresence = newPhasePresence;
+ staticDat_[globalIdx].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
// parameters given in constructor
- std::vector<StaticVars> staticVertexDat_;
+ std::vector<StaticVars> staticDat_;
bool switchFlag_;
};
diff --git a/dumux/implicit/3p3c/3p3cvolumevariables.hh b/dumux/implicit/3p3c/3p3cvolumevariables.hh
index 0049fe5290..edbbf9ebc8 100644
--- a/dumux/implicit/3p3c/3p3cvolumevariables.hh
+++ b/dumux/implicit/3p3c/3p3cvolumevariables.hh
@@ -100,6 +100,8 @@ class ThreePThreeCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
static const Scalar R; // universial gas constant
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
//! The type of the object returned by the fluidState() method
typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
@@ -128,8 +130,13 @@ public:
const MaterialLawParams &materialParams =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
- int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
- int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
+ int globalIdx;
+ if (isBox) // vertex data
+ globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ else
+ globalIdx = problem.model().dofMapper().map(element);
+
+ int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
Scalar temp = Implementation::temperature_(priVars, problem, element, fvGeometry, scvIdx);
fluidState_.setTemperature(temp);
diff --git a/test/implicit/3p3c/Makefile.am b/test/implicit/3p3c/Makefile.am
index 8d4fe42cbc..ca3ec4c029 100644
--- a/test/implicit/3p3c/Makefile.am
+++ b/test/implicit/3p3c/Makefile.am
@@ -1,8 +1,9 @@
-check_PROGRAMS = test_3p3c
+check_PROGRAMS = test_box3p3c test_cc3p3c
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_3p3c_SOURCES = test_3p3c.cc
+test_box3p3c_SOURCES = test_box3p3c.cc
+test_cc3p3c_SOURCES = test_cc3p3c.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/3p3c/infiltrationproblem.hh b/test/implicit/3p3c/infiltrationproblem.hh
index 3c5736ec11..a168f9d629 100644
--- a/test/implicit/3p3c/infiltrationproblem.hh
+++ b/test/implicit/3p3c/infiltrationproblem.hh
@@ -43,7 +43,9 @@ class InfiltrationProblem;
namespace Properties
{
-NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(BoxThreePThreeC, InfiltrationSpatialParams));
+NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(ThreePThreeC, InfiltrationSpatialParams));
+NEW_TYPE_TAG(InfiltrationBoxProblem, INHERITS_FROM(BoxModel, InfiltrationProblem));
+NEW_TYPE_TAG(InfiltrationCCProblem, INHERITS_FROM(CCModel, InfiltrationProblem));
// Set the grid type
SET_TYPE_PROP(InfiltrationProblem, Grid, Dune::YaspGrid<2>);
@@ -143,6 +145,8 @@ class InfiltrationProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief The constructor
@@ -161,6 +165,8 @@ public:
/*pressMin=*/0.8*1e5,
/*pressMax=*/3*1e5,
/*nPress=*/200);
+
+ name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name);
}
/*!
@@ -173,8 +179,8 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const char *name() const
- { return "infiltration"; }
+ const std::string name() const
+ { return name_; }
/*!
* \brief Returns the temperature within the domain.
@@ -210,12 +216,11 @@ public:
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
- * \param vertex The vertex for which the boundary type is set
+ * \param globalPos The position for which the bc type should be evaluated
*/
- void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
if(globalPos[0] > 500. - eps_)
values.setAllDirichlet();
else if(globalPos[0] < eps_)
@@ -229,14 +234,12 @@ public:
* boundary segment.
*
* \param values The dirichlet values for the primary variables
- * \param vertex The vertex for which the boundary type is set
+ * \param globalPos The position for which the bc type should be evaluated
*
* For this method, the \a values parameter stores primary variables.
*/
- void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
Scalar y = globalPos[1];
Scalar x = globalPos[0];
Scalar Sw, Swr=0.12, Sgr=0.03;
@@ -286,8 +289,13 @@ public:
int scvIdx,
const int boundaryFaceIdx) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
values = 0;
+
+ GlobalPosition globalPos;
+ if (isBox)
+ globalPos = element.geometry().corner(scvIdx);
+ else
+ globalPos = is.geometry().center();
// negative values for injection
if ((globalPos[0] <= 75.+eps_) && (globalPos[0] >= 50.+eps_) && (globalPos[1] >= 10.-eps_))
@@ -309,21 +317,14 @@ public:
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry in the box scheme
- * \param scvIdx The local vertex index
+ * \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
* variables.
*/
- void initial(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
- initial_(values, globalPos, element);
+ initial_(values, globalPos);
}
/*!
@@ -345,7 +346,7 @@ private:
// internal method for the initial condition (reused for the
// dirichlet conditions!)
void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos, const Element &element) const
+ const GlobalPosition &globalPos) const
{
Scalar y = globalPos[1];
Scalar x = globalPos[0];
@@ -400,6 +401,7 @@ private:
Scalar temperature_;
Scalar eps_;
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/3p3c/test_box3p3c.cc b/test/implicit/3p3c/test_box3p3c.cc
new file mode 100644
index 0000000000..49e2db460a
--- /dev/null
+++ b/test/implicit/3p3c/test_box3p3c.cc
@@ -0,0 +1,59 @@
+// -*- 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 3p3c box model
+ */
+#include "config.h"
+#include "infiltrationproblem.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";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(InfiltrationBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
+
diff --git a/test/implicit/3p3c/test_box3p3c.input b/test/implicit/3p3c/test_box3p3c.input
new file mode 100644
index 0000000000..954b95f8f8
--- /dev/null
+++ b/test/implicit/3p3c/test_box3p3c.input
@@ -0,0 +1,31 @@
+###############################################################
+# Parameter file for test_3p3c.
+# Everything behind a '#' is a comment.
+# Type "./test_3p3c --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 60 # [s]
+TEnd = 864000 # [s]
+
+[Grid]
+File = ./grids/test_3p3c.dgf
+
+[Problem]
+Name = infiltrationbox
+
+###############################################################
+# Simulation restart
+#
+# DuMux simulations can be restarted from *.drs files
+# Set Restart to the value of a specific file,
+# e.g.: 'Restart = 27184.1' for the restart file
+# name_time=27184.1_rank=0.drs
+# Please comment in the two lines below, if restart is desired.
+###############################################################
+# [TimeManager]
+# Restart = ...
diff --git a/test/implicit/3p3c/test_3p3c.cc b/test/implicit/3p3c/test_cc3p3c.cc
similarity index 98%
rename from test/implicit/3p3c/test_3p3c.cc
rename to test/implicit/3p3c/test_cc3p3c.cc
index eaa2a7b5e5..d60dcee9b2 100644
--- a/test/implicit/3p3c/test_3p3c.cc
+++ b/test/implicit/3p3c/test_cc3p3c.cc
@@ -53,7 +53,7 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(InfiltrationProblem) ProblemTypeTag;
+ typedef TTAG(InfiltrationCCProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/3p3c/test_3p3c.input b/test/implicit/3p3c/test_cc3p3c.input
similarity index 96%
rename from test/implicit/3p3c/test_3p3c.input
rename to test/implicit/3p3c/test_cc3p3c.input
index 2d36d05c82..b6ddf71376 100644
--- a/test/implicit/3p3c/test_3p3c.input
+++ b/test/implicit/3p3c/test_cc3p3c.input
@@ -15,6 +15,9 @@ TEnd = 864000 # [s]
[Grid]
File = ./grids/test_3p3c.dgf
+[Problem]
+Name = infiltrationcc
+
###############################################################
# Simulation restart
#
--
GitLab