diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index d274336b4955895eddf0205ffe4d09ee5e9fb64d..b6efa39179f1d4007fe4689b825b5ea66e79280f 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -118,8 +118,6 @@ public:
int numDofs = asImp_().numDofs();
uCur_.resize(numDofs);
- if (isBox)
- boxVolume_.resize(numDofs);
// apply initial solution
// for compositional models initial the phase presence herein
@@ -235,25 +233,6 @@ public:
storage = gridView_().comm().sum(storage);
}
- /*!
- * \brief Returns the volume \f$\mathrm{[m^3]}\f$ of a given control volume.
- *
- * \param vIdxGlobal The global index of the control volume's
- * associated vertex
- */
- Scalar boxVolume(const int vIdxGlobal) const
- {
- if (isBox)
- {
- return boxVolume_[vIdxGlobal][0];
- }
- else
- {
- DUNE_THROW(Dune::InvalidStateException,
- "requested box volume for cell-centered model");
- }
- }
-
void adaptVariableSize()
{
uCur_.resize(numDofs());
@@ -430,9 +409,6 @@ public:
// update boundary indices
updateBoundaryIndices_();
- if (isBox)
- boxVolume_.resize(asImp_().numDofs());
-
// update the fv geometry
globalFvGeometryPtr_->update(problem_());
@@ -727,62 +703,34 @@ protected:
{
// first set the whole domain to zero
uCur_ = Scalar(0.0);
- boxVolume_ = Scalar(0.0);
- // iterate through leaf grid and evaluate initial
- // condition at the center of each sub control volume
- for (const auto& element : elements(gridView_()))
+ // set the initial values
+ if(isBox)
{
- // deal with the current element, bind FVGeometry to it
- auto fvGeometry = localView(globalFvGeometry());
- fvGeometry.bindElement(element);
-
- // loop over sub control volumes
- for (auto&& scv : scvs(fvGeometry))
+ for (const auto& vertex : vertices(problem_().gridView()))
{
- // let the problem do the dirty work of nailing down
- // the initial solution.
- auto initPriVars = problem_().initial(scv);
-
- auto dofIdxGlobal = scv.dofIndex();
- if (isBox)
- {
- // add up the initial values of all sub-control
- // volumes. If the initial values disagree for
- // different sub control volumes, the initial value
- // will be the arithmetic mean.
- initPriVars *= scv.volume();
- boxVolume_[dofIdxGlobal] += scv.volume();
- }
-
- uCur_[dofIdxGlobal] += initPriVars;
+ const auto dofIdxGlobal = dofMapper().index(vertex);
+ uCur_[dofIdxGlobal] = problem_().initial(vertex);
+ }
+ }
+ else
+ {
+ for (const auto& element : elements(problem_().gridView()))
+ {
+ const auto dofIdxGlobal = dofMapper().index(element);
+ uCur_[dofIdxGlobal] = problem_().initial(element);
}
}
- // add up the primary variables and the volumes of the boxes
- // which cross process borders
- if (isBox && gridView_().comm().size() > 1) {
- VertexHandleSum<Dune::FieldVector<Scalar, 1>,
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> >,
- VertexMapper> sumVolumeHandle(boxVolume_, vertexMapper());
- gridView_().communicate(sumVolumeHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
-
+ // add up the primary variables which cross process borders
+ if (isBox && gridView_().comm().size() > 1)
+ {
VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
sumPVHandle(uCur_, vertexMapper());
gridView_().communicate(sumPVHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
}
-
- if (isBox)
- {
- // divide all primary variables by the volume of their boxes
- for (unsigned int i = 0; i < uCur_.size(); ++i) {
- uCur_[i] /= boxVolume(i);
- }
- }
}
/*!
@@ -860,9 +808,6 @@ protected:
// the finite volume element geometries
std::shared_ptr<GlobalFVGeometry> globalFvGeometryPtr_;
- // container to store the box volumes
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
-
private:
/*!
* \brief Returns whether messages should be printed
diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh
index 28ea314e1808957c8d4e0fd2bf97dfad5b597dc0..29b03f5ad8dcb4bb637a4aeff6aa701ffd7e3e14 100644
--- a/dumux/implicit/problem.hh
+++ b/dumux/implicit/problem.hh
@@ -422,32 +422,24 @@ public:
void pointSourceAtPos(PointSource& pointSource,
const GlobalPosition &globalPos) const {}
+
/*!
- * \brief Evaluate the initial value for a control volume.
+ * \brief Evaluate the initial value for an element (for cc models) or vertex (for box models)
*
- * \param values The initial values for the primary variables
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param scvIdx The local subcontrolvolume index
- *
- * For this method, the \a values parameter stores primary
- * variables.
+ * \param entity The entity (element or vertex)
*/
- PrimaryVariables initial(const SubControlVolume &scv) const
+ template<class Entity>
+ PrimaryVariables initial(const Entity& entity)
{
- // forward to generic interface
- return asImp_().initialAtPos(scv.dofPosition());
+ static_assert(int(Entity::codimension) == 0 || int(Entity::codimension) == dim, "Entity must be element or vertex");
+
+ return asImp_().initialAtPos(entity.geometry().center());
}
/*!
* \brief Evaluate the initial value for a control volume.
*
- * \param values The initial values for the primary variables
- * \param globalPos The position of the center of the finite volume
- * for which the initial values ought to be
- * set (in global coordinates)
- *
- * For this method, the \a values parameter stores primary variables.
+ * \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
@@ -458,33 +450,6 @@ public:
"an initial() or an initialAtPos() method.");
}
- /*!
- * \brief Evaluate the initial phase state inside a control volume.
- *
- * \param vertex The vertex
- * \param vIdxGlobal The global index of the vertex
- * \param globalPos The global position
- */
- int initialPhasePresence(const Vertex &vertex,
- int &vIdxGlobal,
- const GlobalPosition &globalPos) const
- {
- // forward to generic interface
- return asImp_().initialPhasePresenceAtPos(globalPos);
- }
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param globalPos The global position
- */
- int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
- {
- //! As a default, i.e. if the user's problem does not overload any initialPhasePresence method
- //! return 0 (the default phase state is depending on the model context)
- return 0;
- }
-
/*!
* \brief Return how much the domain is extruded at a given sub-control volume.
*
diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
index e9485dd050ee2cd87babf38611d7bbb3db7d823c..6a031ea3e8adcc7f604908abfd50679307b0e556 100644
--- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
@@ -54,6 +54,8 @@ class PrimaryVariableSwitch
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
void init(Problem& problem)
@@ -62,16 +64,20 @@ public:
phasePresence_.resize(numDofs);
wasSwitched_.resize(numDofs, false);
- for (const auto& element : elements(problem.gridView()))
+ if(isBox)
{
- // make sure FVElementGeometry is bound to the element
- auto fvGeometry = localView(problem.model().globalFvGeometry());
- fvGeometry.bindElement(element);
-
- for (auto&& scv : scvs(fvGeometry))
+ for (const auto& vertex : vertices(problem.gridView()))
{
- auto dofIdxGlobal = scv.dofIndex();
- phasePresence_[dofIdxGlobal] = problem.initialPhasePresence(scv);
+ const auto dofIdxGlobal = problem.model().dofMapper().index(vertex);
+ phasePresence_[dofIdxGlobal] = problem.initialPhasePresence(vertex);
+ }
+ }
+ else
+ {
+ for (const auto& element : elements(problem.gridView()))
+ {
+ const auto dofIdxGlobal = problem.model().dofMapper().index(element);
+ phasePresence_[dofIdxGlobal] = problem.initialPhasePresence(element);
}
}
diff --git a/dumux/porousmediumflow/implicit/problem.hh b/dumux/porousmediumflow/implicit/problem.hh
index 10ecfdd98ee0dc75bfe45ab56cb2d9d63c2fdff2..085ab189ce7d04a00334292646a22aec753a85b4 100644
--- a/dumux/porousmediumflow/implicit/problem.hh
+++ b/dumux/porousmediumflow/implicit/problem.hh
@@ -41,23 +41,25 @@ NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered i
template<class TypeTag>
class ImplicitPorousMediaProblem : public ImplicitProblem<TypeTag>
{
- typedef ImplicitProblem<TypeTag> ParentType;
+ using ParentType = ImplicitProblem<TypeTag>;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
- typedef typename GridView::ctype CoordScalar;
- typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
+
+ using CoordScalar = typename GridView::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
@@ -98,6 +100,32 @@ public:
*/
// \{
+ /*!
+ * \brief Evaluate the initial phase state at an element (for cc models) or vertex (for box models)
+ *
+ * \param entity The entity (element or vertex)
+ */
+ template<class Entity>
+ int initialPhasePresence(const Entity& entity)
+ {
+ static_assert(int(Entity::codimension) == 0 || int(Entity::codimension) == dim, "Entity must be element or vertex");
+
+ return asImp_().initialPhasePresenceAtPos(entity.geometry().center());
+ }
+
+ /*!
+ * \brief Evaluate the initial phase state at a given position
+ *
+ * \param globalPos The global position
+ */
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
+ {
+ DUNE_THROW(Dune::InvalidStateException,
+ "The problem does not provide a initialPhasePresence() "
+ "or initialPhasePresenceAtPos() method.");
+ return 0;
+ }
+
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at a given global position.
*
diff --git a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
index 2b94ffc26ff469142398931edd6cb02d0d644ebb..5168fb183673a0ed8609b552dc0b6958896680bd 100644
--- a/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/injectionproblem.hh
@@ -323,13 +323,11 @@ public:
}
/*!
- * \brief Returns the initial phase state for a control volume.
+ * \brief Evaluate the initial phase state at a given position
*
- * \param vertex The vertex
- * \param vIdxGlobal The global index of the vertex
* \param globalPos The global position
*/
- int initialPhasePresence(const SubControlVolume& scv) const
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
{ return Indices::wPhaseOnly; }
// \}
diff --git a/test/porousmediumflow/2p2c/implicit/waterairproblem.hh b/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
index 407ff58d52be72c81b2e73212cf09be1c434ac58..ade3380bab18227e7b0b6a92cd222e12c558bb61 100644
--- a/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
+++ b/test/porousmediumflow/2p2c/implicit/waterairproblem.hh
@@ -310,16 +310,12 @@ public:
}
/*!
- * \brief Return the initial phase state inside a control volume.
+ * \brief Evaluate the initial phase state at a given position
*
- * \param vertex The vertex
- * \param vIdxGlobal The global index of the vertex
* \param globalPos The global position
*/
- int initialPhasePresence(const SubControlVolume& scv) const
- {
- return wPhaseOnly;
- }
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
+ { return wPhaseOnly; }
private:
// internal method for the initial condition (reused for the
diff --git a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
index 9733d712667625cfc601521e6c6a600cbc60a381..1bf5cf015e20a73eb3b2ed1131253315d21fb447 100644
--- a/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
+++ b/test/porousmediumflow/2pnc/implicit/fuelcellproblem.hh
@@ -254,14 +254,14 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{ return initial_(globalPos); }
- /*!
- * \brief Return the initial phase state inside a sub control volume.
- *
- * \param scv The sub control volume
- */
- int initialPhasePresence(const SubControlVolume& scv) const
- { return Indices::bothPhases; }
+ /*!
+ * \brief Evaluate the initial phase state at a given position
+ *
+ * \param globalPos The global position
+ */
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
+ { return Indices::bothPhases; }
/*!
* \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
index 2f7313849d03d8e883d44f3bc6e5eef22fff72f4..0665d375f11e8b36d828f789c431faf51f2faa2e 100644
--- a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh
@@ -261,18 +261,13 @@ public:
}
/*!
- * \brief Return the initial phase state inside a control volume.
+ * \brief Evaluate the initial phase state at a given position
*
- * \param vertex The vertex
- * \param vIdxGlobal The global index of the vertex
* \param globalPos The global position
*/
- int initialPhasePresence(const Vertex &vertex,
- int &vIdxGlobal,
- const GlobalPosition &globalPos) const
- {
- return threePhases;
- }
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
+ { return threePhases; }
+
/*!
* \brief Append all quantities of interest which can be derived
diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
index 878a0b9f59fbb4c9370bbadd9d14b8b06aa13b60..a258c8d71fa193a7eb1a4d89e7086b64707af72a 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
@@ -318,14 +318,12 @@ public:
}
/*!
- * \brief Return the initial phase state inside a control volume.
+ * \brief Evaluate the initial phase state at a given position
*
* \param globalPos The global position
*/
- int initialPhasePresence(const SubControlVolume& scv) const
- {
- return wgPhaseOnly;
- }
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
+ { return wgPhaseOnly; }
private:
// internal method for the initial condition
diff --git a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
index 836c372ccb9c23469cfa2105630aba10f6cd14b1..82a3b8a190f464701d8db7b1f5fcc7202058bf0c 100644
--- a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh
@@ -268,15 +268,11 @@ public:
}
/*!
- * \brief Return the initial phase state inside a control volume.
+ * \brief Evaluate the initial phase state at a given position
*
- * \param vertex The vertex
- * \param vIdxGlobal The global index of the vertex
* \param globalPos The global position
*/
- int initialPhasePresence(const Vertex &vertex,
- const int &vIdxGlobal,
- const GlobalPosition &globalPos) const
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos)
{
if (isInContaminationZone(globalPos))
return threePhases;
@@ -284,7 +280,7 @@ public:
return wgPhaseOnly;
}
- /*!
+ /*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
* writer. Adjust this in case of anisotropic permeabilities.