From fb9d0b46110d9365b918679108f17a76cec1961b Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 21 Mar 2011 14:55:19 +0000
Subject: [PATCH] box models: add a few additional valgrind checks, allow
problems to cap time step sizes more easily
to cap the time step size a problem can overload the maxTimeStepSize() method.
Conflicts:
dumux/boxmodels/common/boxproblem.hh
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5439 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/common/boxlocaljacobian.hh | 3 +-
dumux/boxmodels/common/boxlocalresidual.hh | 49 +++++++-
dumux/boxmodels/common/boxproblem.hh | 119 ++++++++++++++++++-
dumux/boxmodels/common/boxvolumevariables.hh | 24 +++-
dumux/common/valgrind.hh | 7 +-
dumux/nonlinear/newtoncontroller.hh | 3 +-
6 files changed, 190 insertions(+), 15 deletions(-)
diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index a90ecb235e..c3b557a0b2 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -353,6 +353,7 @@ protected:
int pvIdx)
{
int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
+
PrimaryVariables priVars(model_().curSol()[globalIdx]);
VolumeVariables origVolVars(curVolVars_[scvIdx]);
@@ -445,7 +446,7 @@ protected:
Scalar numericEpsilon_(int scvIdx,
int pvIdx) const
{
- Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
+ Scalar pv = this->curVolVars_[scvIdx].primaryVar(pvIdx);
return 1e-9*(std::abs(pv) + 1);
}
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index 495670bc47..ffa9008bf2 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -269,10 +269,10 @@ public:
Valgrind::CheckDefined(prevVolVars);
Valgrind::CheckDefined(curVolVars);
-#if 0 // HAVE_VALGRIND
+#if !defined NDEBUG && HAVE_VALGRIND
for (int i=0; i < fvGeom.numVertices; i++) {
- Valgrind::CheckDefined(prevVolVars[i]);
- Valgrind::CheckDefined(curVolVars[i]);
+ prevVolVars[i].checkDefined();
+ curVolVars[i].checkDefined();
}
#endif // HAVE_VALGRIND
@@ -287,13 +287,29 @@ public:
residual_.resize(numVerts);
residual_ = 0;
+#if !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
+
asImp_().evalFluxes_();
+
+#if !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
+
asImp_().evalVolumeTerms_();
+#if !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
+
// evaluate the boundary
asImp_().evalBoundary_();
-#if HAVE_VALGRIND
+#if !defined NDEBUG && HAVE_VALGRIND
for (int i=0; i < fvElemGeom_().numVertices; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
@@ -340,8 +356,18 @@ protected:
{
if (bcTypes_().hasNeumann())
asImp_().evalNeumann_();
+#if !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
+
if (bcTypes_().hasDirichlet())
asImp_().evalDirichlet_();
+
+#if !defined NDEBUG && HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
}
/*!
@@ -365,8 +391,12 @@ protected:
if (!bcTypes.isDirichlet(eqIdx))
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ assert(0 <= pvIdx && pvIdx < numEq);
+ Valgrind::CheckDefined(pvIdx);
+ Valgrind::CheckDefined(curPrimaryVar_(i, pvIdx));
+ Valgrind::CheckDefined(tmp[pvIdx]);
this->residual_[i][eqIdx] =
- curPrimaryVars_(i)[pvIdx] - tmp[pvIdx];
+ curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
};
};
}
@@ -583,6 +613,15 @@ protected:
return curVolVars_(i).primaryVars();
}
+ /*!
+ * \brief Returns the j'th primary of the i'th sub-control volume
+ * of the current element.
+ */
+ Scalar curPrimaryVar_(int i, int j) const
+ {
+ return curVolVars_(i).primaryVar(j);
+ }
+
/*!
* \brief Returns a reference to the current volume variables of
* all sub-control volumes of the current element.
diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index 6040bda0f4..070a30f518 100644
--- a/dumux/boxmodels/common/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -65,11 +65,18 @@ private:
dimWorld = GridView::dimensionworld
};
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename GridView::Intersection Intersection;
+
typedef typename GridView::Grid::ctype CoordScalar;
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
// copying a problem is not a good idea
BoxProblem(const BoxProblem &);
@@ -128,6 +135,14 @@ public:
model().init(asImp_());
}
+ /*!
+ * \brief Returns the maximum allowed time step size [s]
+ *
+ * By default this the time step size is unrestricted.
+ */
+ Scalar maxTimeStepSize() const
+ { return std::numeric_limits<Scalar>::infinity(); }
+
/*!
* \brief If model coupling is used, this updates the parameters
* required to calculate the coupling fluxes between the
@@ -352,6 +367,108 @@ public:
{ return model_; }
// \}
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * 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
+ */
+ void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ {
+ // Throw an exception
+ DUNE_THROW(Dune::NotImplemented,
+ "The problem specifies does not provide "
+ "a boundaryTypes() method.");
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param values The dirichlet values for the primary variables
+ * \param pos The position of the center of the finite volume
+ * for which the dirichlet condition ought to be
+ * set in global coordinates
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVariables &values,
+ const GlobalPosition &pos) const
+ {
+ // Throw an exception (there is no reasonable default value
+ // for Dirichlet conditions)
+ DUNE_THROW(Dune::NotImplemented,
+ "The problem specifies that some boundary "
+ "segments are dirichlet, but does not provide "
+ "a dirichlet() method.");
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * \param values The neumann values for the conservation equations [kg / (m^2 *s )]
+ * \param element The finite element
+ * \param fvElemGeom The finite-volume geometry in the box scheme
+ * \param is The intersection between element and boundary
+ * \param scvIdx The local vertex index
+ * \param boundaryFaceIdx The index of the boundary face
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ // Throw an exception
+ DUNE_THROW(Dune::NotImplemented,
+ "The problem specifies that some boundary "
+ "segments are dirichlet, but does not provide "
+ "a neumann() method.");
+ };
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ // Throw an exception
+ DUNE_THROW(Dune::NotImplemented,
+ "The problem specifies does not provide "
+ "a initial() method.");
+ }
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ // Throw an exception
+ DUNE_THROW(Dune::NotImplemented,
+ "The problem specifies does not provide "
+ "a source() method.");
+ }
+
/*!
* \name Restart mechanism
*/
diff --git a/dumux/boxmodels/common/boxvolumevariables.hh b/dumux/boxmodels/common/boxvolumevariables.hh
index 724a6b903d..cd6dd02439 100644
--- a/dumux/boxmodels/common/boxvolumevariables.hh
+++ b/dumux/boxmodels/common/boxvolumevariables.hh
@@ -115,8 +115,8 @@ public:
int scvIdx,
bool isOldSol)
{
+ Valgrind::CheckDefined(priVars);
primaryVars_ = priVars;
- Valgrind::CheckDefined(primaryVars_);
}
/*!
@@ -131,7 +131,23 @@ public:
* \param pvIdx The index of the primary variable of interest
*/
Scalar primaryVar(int pvIdx) const
- { return primaryVars_[pvIdx]; }
+ {
+ return primaryVars_[pvIdx];
+ }
+
+ /*!
+ * \brief If running in valgrind this makes sure that all
+ * quantities in the volume variables are defined.
+ */
+ void checkDefined() const
+ {
+#if !defined NDEBUG && HAVE_VALGRIND
+ Valgrind::CheckDefined(primaryVars_);
+ Valgrind::CheckDefined(evalPoint_);
+ if (evalPoint_ && evalPoint_ != this)
+ evalPoint_->checkDefined();
+#endif
+ };
protected:
const Implementation &asImp_() const
@@ -139,10 +155,10 @@ protected:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
- PrimaryVariables primaryVars_;
-
// the evaluation point of the local jacobian
const Implementation *evalPoint_;
+
+ PrimaryVariables primaryVars_;
};
} // end namepace
diff --git a/dumux/common/valgrind.hh b/dumux/common/valgrind.hh
index e4073a8ce7..1bc26efda3 100644
--- a/dumux/common/valgrind.hh
+++ b/dumux/common/valgrind.hh
@@ -44,7 +44,7 @@
#define SetNoAccess(a) foo();
namespace Valgrind
{
-inline void foo() {}; // dummy function
+inline bool foo() { return true; }; // dummy function
}
#else
@@ -67,10 +67,11 @@ namespace Valgrind
* \param value the object which valgrind should check
*/
template <class T>
-inline void CheckDefined(const T &value)
+inline bool CheckDefined(const T &value)
{
#if HAVE_VALGRIND
- VALGRIND_CHECK_MEM_IS_DEFINED(&value, sizeof(T));
+ unsigned int tmp = VALGRIND_CHECK_MEM_IS_DEFINED(&value, sizeof(T));
+ return tmp == 0;
#endif
}
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 8b7dbec540..428a67d4a0 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -622,7 +622,8 @@ public:
return oldTimeStep*(1 + percent);
*/
Scalar percent = (targetSteps_ - n)/targetSteps_;
- return oldTimeStep*(1.0 + percent/1.2);
+ return std::min(oldTimeStep*(1.0 + percent/1.2),
+ this->problem_().maxTimeStepSize());
}
}
--
GitLab