diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index 82a91f9e163adf24780ba9419cebdc846fc84550..19e9475954ab17c3fcfda0fcce3a7217ba7b6c02 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 ac9719c58b76bc0c2f582a87ead85bc54d624d93..9591cae3f738fc3650f0394aa321e0e5278667c0 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 c12cadbb93c8ed1289cdab906c0185f883aaefbd..9d8d48cd64ae30b8ec8b88f2b1c8e9ea29fada30 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 (isBox) // vertex data
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = idx;
+
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 0049fe529009a255bbe97bd25515a25b578ce8ff..edbbf9ebc81228fc59e27e11112458c6fb53a570 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 8d4fe42cbc9af8902bcdc2ccea42beb45f1e9946..ca3ec4c029acf6c1ed751230ab9043db79f046b1 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 3c5736ec1199a8ba8602ea7645a0616451553103..a168f9d6295caa6478b7c9e9f55f219de3481a6f 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 0000000000000000000000000000000000000000..49e2db460a3bc7f2ddd596c8bd32591510357e07
--- /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 0000000000000000000000000000000000000000..954b95f8f805bf1b548315f716cc7f61d47f3ad0
--- /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 eaa2a7b5e57d79d94051907f04a3737eee540830..d60dcee9b20f6260ded10208e11cc11af9fcdcdc 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 2d36d05c820c2010b63209e6d0ec7d59a051e5d9..b6ddf71376303335aa4e4baeaa6510373673ae0c 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
#
diff --git a/test/implicit/3p3cni/Makefile.am b/test/implicit/3p3cni/Makefile.am
index f1eaf98125f5526a15be290645e4a2680f500c12..5eaa31dae3c76a6b2ee378380f7c6dad6bccfde2 100644
--- a/test/implicit/3p3cni/Makefile.am
+++ b/test/implicit/3p3cni/Makefile.am
@@ -1,8 +1,9 @@
-check_PROGRAMS = test_3p3cni
+check_PROGRAMS = test_box3p3cni test_cc3p3cni
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_3p3cni_SOURCES = test_3p3cni.cc
+test_box3p3cni_SOURCES = test_box3p3cni.cc
+test_cc3p3cni_SOURCES = test_cc3p3cni.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/3p3cni/columnxylolproblem.hh b/test/implicit/3p3cni/columnxylolproblem.hh
index 0bc6e1f9820ffc80c25ff7234a52fea6217f79a2..c9193430a62c324246dd4d2008886f6c57c52a03 100644
--- a/test/implicit/3p3cni/columnxylolproblem.hh
+++ b/test/implicit/3p3cni/columnxylolproblem.hh
@@ -46,8 +46,10 @@ class ColumnProblem;
namespace Properties
{
-NEW_TYPE_TAG(ColumnProblem, INHERITS_FROM(BoxThreePThreeCNI, ColumnSpatialParams));
-
+NEW_TYPE_TAG(ColumnProblem, INHERITS_FROM(ThreePThreeCNI, ColumnSpatialParams));
+NEW_TYPE_TAG(ColumnBoxProblem, INHERITS_FROM(BoxModel, ColumnProblem));
+NEW_TYPE_TAG(ColumnCCProblem, INHERITS_FROM(CCModel, ColumnProblem));
+
// Set the grid type
SET_PROP(ColumnProblem, Grid)
{
@@ -124,6 +126,8 @@ class ColumnProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief The constructor
@@ -135,6 +139,8 @@ public:
: ParentType(timeManager, gridView), eps_(1e-6)
{
FluidSystem::init();
+
+ name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name);
}
/*!
@@ -147,8 +153,8 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const char *name() const
- { return "columnxylol"; }
+ const std::string name() const
+ { return name_; }
void sourceAtPos(PrimaryVariables &values,
@@ -169,12 +175,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[1] < eps_)
values.setAllDirichlet();
else
@@ -186,14 +191,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();
-
initial_(values, globalPos);
}
@@ -218,9 +221,14 @@ public:
const int scvIdx,
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[1] > 1.2 - eps_)
{
@@ -242,22 +250,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 fvGeomtry 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 &fvGeomtry,
- const int scvIdx) const
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
initial_(values, globalPos);
-
}
/*!
@@ -274,48 +274,6 @@ public:
return threePhases;
}
- /*
- bool shouldWriteOutput() const
- {
- return
- this->timeManager().timeStepIndex() == 0 ||
- this->timeManager().timeStepIndex() == 10 ||
- this->timeManager().timeStepIndex() == 20 ||
- this->timeManager().timeStepIndex() == 50 ||
- this->timeManager().timeStepIndex() == 100 ||
- this->timeManager().timeStepIndex() == 200 ||
- this->timeManager().timeStepIndex() == 300 ||
- this->timeManager().timeStepIndex() == 400 ||
- this->timeManager().timeStepIndex() == 500 ||
- this->timeManager().timeStepIndex() == 600 ||
- this->timeManager().timeStepIndex() == 700 ||
- this->timeManager().timeStepIndex() == 800 ||
- this->timeManager().timeStepIndex() == 900 ||
- this->timeManager().timeStepIndex() == 1000 ||
- this->timeManager().timeStepIndex() == 1100 ||
- this->timeManager().timeStepIndex() == 1200 ||
- this->timeManager().timeStepIndex() == 1300 ||
- this->timeManager().timeStepIndex() == 1400 ||
- this->timeManager().timeStepIndex() == 1500 ||
- this->timeManager().timeStepIndex() == 1600 ||
- this->timeManager().timeStepIndex() == 1700 ||
- this->timeManager().timeStepIndex() == 1800 ||
- this->timeManager().timeStepIndex() == 1900 ||
- this->timeManager().timeStepIndex() == 2000 ||
- this->timeManager().timeStepIndex() == 2100 ||
- this->timeManager().timeStepIndex() == 2200 ||
- this->timeManager().timeStepIndex() == 2300 ||
- this->timeManager().timeStepIndex() == 2400 ||
- this->timeManager().timeStepIndex() == 2500 ||
- this->timeManager().timeStepIndex() == 2600 ||
- this->timeManager().timeStepIndex() == 2700 ||
- this->timeManager().timeStepIndex() == 2800 ||
- this->timeManager().timeStepIndex() == 2900 ||
- this->timeManager().willBeFinished();
- }
- */
-
-
private:
// internal method for the initial condition (reused for the
// dirichlet conditions!)
@@ -375,6 +333,7 @@ private:
}
const Scalar eps_;
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/3p3cni/kuevetteproblem.hh b/test/implicit/3p3cni/kuevetteproblem.hh
index f6ead4d8f7f68abc6faa7e83574254076c8ebde3..b193e5a610f894662eb0664888c9a4916feeaa35 100644
--- a/test/implicit/3p3cni/kuevetteproblem.hh
+++ b/test/implicit/3p3cni/kuevetteproblem.hh
@@ -46,8 +46,10 @@ class KuevetteProblem;
namespace Properties
{
-NEW_TYPE_TAG(KuevetteProblem, INHERITS_FROM(BoxThreePThreeCNI, KuevetteSpatialParams));
-
+NEW_TYPE_TAG(KuevetteProblem, INHERITS_FROM(ThreePThreeCNI, KuevetteSpatialParams));
+NEW_TYPE_TAG(KuevetteBoxProblem, INHERITS_FROM(BoxModel, KuevetteProblem));
+NEW_TYPE_TAG(KuevetteCCProblem, INHERITS_FROM(CCModel, KuevetteProblem));
+
// Set the grid type
SET_PROP(KuevetteProblem, Grid)
{
@@ -158,6 +160,8 @@ class KuevetteProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief The constructor
@@ -169,6 +173,8 @@ public:
: ParentType(timeManager, gridView), eps_(1e-6)
{
FluidSystem::init();
+
+ name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name);
}
/*!
@@ -181,9 +187,9 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const char *name() const
- { return "kuevette3p3cni"; }
-
+ const std::string name() const
+ { return name_; }
+
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
@@ -203,12 +209,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] > 1.5 - eps_)
values.setAllDirichlet();
else
@@ -221,14 +226,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();
-
initial_(values, globalPos);
}
@@ -253,9 +256,14 @@ public:
const int scvIdx,
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] < eps_)
{
@@ -277,22 +285,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 fvGeomtry 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 &fvGeomtry,
- const int scvIdx) const
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
initial_(values, globalPos);
-
}
/*!
@@ -332,6 +332,7 @@ private:
}
const Scalar eps_;
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/3p3cni/kuevettespatialparams.hh b/test/implicit/3p3cni/kuevettespatialparams.hh
index 40f502991f7a915e7778daa5a0d45066f77fcf61..d3487b383c2559e49fb9df29b002fbdbad146af6 100644
--- a/test/implicit/3p3cni/kuevettespatialparams.hh
+++ b/test/implicit/3p3cni/kuevettespatialparams.hh
@@ -159,14 +159,14 @@ public:
* potential gradient.
*
* \param element The current finite element
- * \param fvGeomtry The current finite volume geometry of the element
+ * \param fvGeometry The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
const Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvGeomtry,
+ const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
- const GlobalPosition &pos = fvGeomtry.subContVol[scvIdx].global;
+ const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return fineK_;
return coarseK_;
@@ -176,15 +176,15 @@ public:
* \brief Define the porosity \f$[-]\f$ of the spatial parameters
*
* \param element The finite element
- * \param fvGeomtry The finite volume geometry
+ * \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* the porosity needs to be defined
*/
double porosity(const Element &element,
- const FVElementGeometry &fvGeomtry,
+ const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
- const GlobalPosition &pos = fvGeomtry.subContVol[scvIdx].global;
+ const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return finePorosity_;
else
@@ -196,14 +196,14 @@ public:
* \brief return the parameter object for the Brooks-Corey material law which depends on the position
*
* \param element The current finite element
- * \param fvGeomtry The current finite volume geometry of the element
+ * \param fvGeometry The current finite volume geometry of the element
* \param scvIdx The index of the sub-control volume
*/
const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvGeomtry,
+ const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
- const GlobalPosition &pos = fvGeomtry.subContVol[scvIdx].global;
+ const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
if (isFineMaterial_(pos))
return fineMaterialParams_;
else
@@ -216,18 +216,18 @@ public:
* This is only required for non-isothermal models.
*
* \param element The finite element
- * \param fvGeomtry The finite volume geometry
+ * \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
* the heat capacity needs to be defined
*/
double heatCapacity(const Element &element,
- const FVElementGeometry &fvGeomtry,
+ const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return
850 // specific heat capacity [J / (kg K)]
* 2650 // density of sand [kg/m^3]
- * (1 - porosity(element, fvGeomtry, scvIdx));
+ * (1 - porosity(element, fvGeometry, scvIdx));
}
/*!
@@ -241,7 +241,7 @@ public:
* \param elemVolVars The volume variables
* \param tempGrad The temperature gradient
* \param element The current finite element
- * \param fvGeomtry The finite volume geometry of the current element
+ * \param fvGeometry The finite volume geometry of the current element
* \param faceIdx The local index of the sub-control volume face where
* the matrix heat flux should be calculated
*/
@@ -250,7 +250,7 @@ public:
const ElementVolumeVariables &elemVolVars,
const DimVector &tempGrad,
const Element &element,
- const FVElementGeometry &fvGeomtry,
+ const FVElementGeometry &fvGeometry,
const int faceIdx) const
{
static const Scalar ldry = 0.35;
@@ -258,8 +258,8 @@ public:
static const Scalar lSn1 = 0.65;
// arithmetic mean of the liquid saturation and the porosity
- const int i = fvGeomtry.subContVolFace[faceIdx].i;
- const int j = fvGeomtry.subContVolFace[faceIdx].j;
+ const int i = fvGeometry.subContVolFace[faceIdx].i;
+ const int j = fvGeometry.subContVolFace[faceIdx].j;
Scalar Sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
elemVolVars[j].saturation(wPhaseIdx)) / 2);
Scalar Sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
diff --git a/test/implicit/3p3cni/test_3p3cni.cc b/test/implicit/3p3cni/test_box3p3cni.cc
similarity index 95%
rename from test/implicit/3p3cni/test_3p3cni.cc
rename to test/implicit/3p3cni/test_box3p3cni.cc
index 3b5710201bb5d0fbdfde1652021930f1eab51103..6eb5ff6b879662d778649e22768718745949c4e8 100644
--- a/test/implicit/3p3cni/test_3p3cni.cc
+++ b/test/implicit/3p3cni/test_box3p3cni.cc
@@ -24,7 +24,6 @@
#include "config.h"
//#include "kuevetteproblem.hh"
#include "columnxylolproblem.hh"
-//#include "permdruckproblem.hh"
#include <dumux/common/start.hh>
@@ -56,7 +55,7 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
-// typedef TTAG(KuevetteProblem) ProblemTypeTag;
- typedef TTAG(ColumnProblem) ProblemTypeTag;
+// typedef TTAG(KuevetteBoxProblem) ProblemTypeTag;
+ typedef TTAG(ColumnBoxProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/3p3cni/test_box3p3cni.input b/test/implicit/3p3cni/test_box3p3cni.input
new file mode 100644
index 0000000000000000000000000000000000000000..c1e0432b0500fc891a5f70313e52827bd6a08b75
--- /dev/null
+++ b/test/implicit/3p3cni/test_box3p3cni.input
@@ -0,0 +1,31 @@
+###############################################################
+# Parameter file for test_3p3cni.
+# Everything behind a '#' is a comment.
+# Type "./test_3p3cni --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 200 # [s]
+
+[Grid]
+File = ./grids/column.dgf
+
+[Problem]
+Name = columnxylolbox # name passed to the output routines
+
+###############################################################
+# 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/3p3cni/test_cc3p3cni.cc b/test/implicit/3p3cni/test_cc3p3cni.cc
new file mode 100644
index 0000000000000000000000000000000000000000..69a2abe259781bbd8e76fc655faf5b73917a28e3
--- /dev/null
+++ b/test/implicit/3p3cni/test_cc3p3cni.cc
@@ -0,0 +1,61 @@
+// -*- 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 3p3cni box model
+ */
+#include "config.h"
+//#include "kuevetteproblem.hh"
+#include "columnxylolproblem.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(KuevetteCCProblem) ProblemTypeTag;
+ typedef TTAG(ColumnCCProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/3p3cni/test_3p3cni.input b/test/implicit/3p3cni/test_cc3p3cni.input
similarity index 93%
rename from test/implicit/3p3cni/test_3p3cni.input
rename to test/implicit/3p3cni/test_cc3p3cni.input
index 7be6df79dd5a8cc0ffbb55d069a25ea3f9f65fc4..d41669fddecc208d4b80db50b4c5099f86594ba4 100644
--- a/test/implicit/3p3cni/test_3p3cni.input
+++ b/test/implicit/3p3cni/test_cc3p3cni.input
@@ -15,6 +15,9 @@ TEnd = 200 # [s]
[Grid]
File = ./grids/column.dgf
+[Problem]
+Name = columnxylolcc # name passed to the output routines
+
###############################################################
# Simulation restart
#