Commit 23bc012a authored by Martin Schneider's avatar Martin Schneider
Browse files

Changed floating point comparisons. For comparisons with zero an absolute...

Changed floating point comparisons. For comparisons with zero an absolute criteria is used, otherwise a relative one.

Reviewed by Timo 

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14145 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent e1b42df9
......@@ -123,7 +123,7 @@ int invertLinearPolynomial(SolContainer &sol,
* The polynomial is defined as
* \f[ p(x) = a\; x^2 + + b\;x + c \f]
*
* This method teturns the number of solutions which are in the real
* This method returns the number of solutions which are in the real
* numbers. The "sol" argument contains the real roots of the parabola
* in order with the smallest root first.
*
......
......@@ -28,6 +28,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/tuples.hh>
#include <dune/common/float_cmp.hh>
#include "valgrind.hh"
#include "math.hh"
......@@ -236,7 +237,7 @@ public:
{
assert(applies(x0));
assert(applies(x1));
assert(x0 != x1);
assert(Dune::FloatCmp::ne<Scalar>(x0, x1));
// make sure that x0 is smaller than x1
if (x0 > x1)
......
......@@ -23,6 +23,7 @@
#ifndef DUMUX_TIME_MANAGER_HH
#define DUMUX_TIME_MANAGER_HH
#include <dune/common/float_cmp.hh>
#include <dune/common/timer.hh>
#include <dune/common/parallel/mpihelper.hh>
......@@ -410,7 +411,7 @@ public:
problem_->episodeEnd();
//check if a time step size was explicitly defined in problem->episodeEnd()
if (dt == timeStepSize())
if (Dune::FloatCmp::eq<Scalar>(dt, timeStepSize()))
{
// set the initial time step size of a an episode to the last real time step size before the episode
Scalar nextDt = std::max(previousTimeStepSize_, timeStepSize());
......
......@@ -19,6 +19,8 @@
#ifndef DUMUX_FVPRESSURE2P_HH
#define DUMUX_FVPRESSURE2P_HH
#include <dune/common/float_cmp.hh>
// dumux environment
#include <dumux/decoupled/common/fv/fvpressure.hh>
#include <dumux/decoupled/2p/diffusion/diffusionproperties2p.hh>
......@@ -748,8 +750,8 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] = (potentialDiffW == 0.) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw == 0.) ? rhoMeanNw : density_[nPhaseIdx];
density_[wPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? rhoMeanNw : density_[nPhaseIdx];
potentialDiffW = cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx);
potentialDiffNw = cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx);
......@@ -761,17 +763,17 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0) ? lambdaNwI : lambdaNwJ;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
if (compressibility_)
{
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] = (potentialDiffW == 0) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw == 0) ? rhoMeanNw : density_[nPhaseIdx];
density_[wPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? rhoMeanNw : density_[nPhaseIdx];
}
Scalar scalarPerm = permeability.two_norm();
......@@ -957,8 +959,8 @@ const Intersection& intersection, const CellData& cellData, const bool first)
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
density_[wPhaseIdx] = (potentialDiffW == 0.) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw == 0.) ? rhoMeanNw : density_[nPhaseIdx];
density_[wPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? rhoMeanNw : density_[nPhaseIdx];
}
//calculate potential gradient
......@@ -990,16 +992,16 @@ const Intersection& intersection, const CellData& cellData, const bool first)
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
if (compressibility_)
{
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
density_[wPhaseIdx] = (potentialDiffW == 0) ? rhoMeanW : density_[wPhaseIdx];
density_[wPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? rhoMeanW : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
density_[nPhaseIdx] = (potentialDiffNw == 0) ? rhoMeanNw : density_[nPhaseIdx];
density_[nPhaseIdx] = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? rhoMeanNw : density_[nPhaseIdx];
}
Scalar scalarPerm = permeability.two_norm();
......
......@@ -19,6 +19,8 @@
#ifndef DUMUX_FVPRESSURE2P_ADAPTIVE_HH
#define DUMUX_FVPRESSURE2P_ADAPTIVE_HH
#include <dune/common/float_cmp.hh>
// dumux environment
#include <dumux/decoupled/2p/2pproperties.hh>
#include <dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh>
......@@ -453,20 +455,20 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
lambdaWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
Scalar lambdaNwIJ = (potentialNwIJ > 0) ? lambdaNwI : lambdaNwJ;
lambdaNwIJ = (potentialNwIJ == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
lambdaNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
if (compressibility_)
{
densityWIJ = (potentialWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
densityNwIJ = (potentialNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (potentialNwIJ == 0) ? rhoMeanNwIJ : densityNwIJ;
densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
densityWIK = (potentialWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
densityNwIK = (potentialNwIK > 0.) ? cellData.density(nPhaseIdx) : densityNwIK;
densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
densityNwIK = (potentialNwIK == 0) ? rhoMeanNwIK : densityNwIK;
densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
}
// update diagonal entry and right hand side
......
......@@ -24,6 +24,7 @@
* @brief Velocity Field from a finite volume solution of a pressure equation.
*/
#include <dune/common/float_cmp.hh>
#include <dune/grid/common/gridenums.hh>
#include <dumux/decoupled/2p/diffusion/diffusionproperties2p.hh>
......@@ -384,10 +385,10 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
density_[wPhaseIdx];
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
density_[nPhaseIdx];
potentialDiffW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx));
......@@ -406,9 +407,9 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwJ;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
if (compressibility_)
{
......@@ -416,10 +417,10 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
density_[wPhaseIdx];
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
density_[nPhaseIdx];
}
......@@ -616,10 +617,10 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
{
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
}
//calculate potential gradient
......@@ -643,18 +644,18 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
if (compressibility_)
{
density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
}
Scalar scalarPerm = permeability.two_norm();
......
......@@ -24,6 +24,7 @@
* @brief Velocity Field from a finite volume solution of a pressure equation.
*/
#include <dune/common/float_cmp.hh>
#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
namespace Dumux
......@@ -315,10 +316,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
densityWIJ = (potentialDiffWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (potentialDiffNwIJ == 0.) ? rhoMeanNwIJ : densityNwIJ;
densityWIK = (potentialDiffWIK == 0.) ? rhoMeanWIK : densityWIK;
densityNwIK = (potentialDiffNwIK == 0.) ? rhoMeanNwIK : densityNwIK;
densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
}
Scalar fractionalWIJ = (potentialDiffWIJ > 0.) ? fractionalWI : fractionalWJ;
......@@ -326,10 +327,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
Scalar fractionalWIK = (potentialDiffWIK > 0.) ? fractionalWI : fractionalWK;
Scalar fractionalNwIK = (potentialDiffNwIK > 0.) ? fractionalNwI : fractionalNwK;
fractionalWIJ = (potentialDiffWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
fractionalNwIJ = (potentialDiffNwIJ == 0.) ? fMeanNwIJ : fractionalNwIJ;
fractionalWIK = (potentialDiffWIK == 0.) ? fMeanWIK : fractionalWIK;
fractionalNwIK = (potentialDiffNwIK == 0.) ? fMeanNwIK : fractionalNwIK;
fractionalWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? fMeanWIJ : fractionalWIJ;
fractionalNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? fMeanNwIJ : fractionalNwIJ;
fractionalWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? fMeanWIK : fractionalWIK;
fractionalNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? fMeanNwIK : fractionalNwIK;
switch (pressureType_)
{
......@@ -375,20 +376,20 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaWIJ = (potentialDiffWIJ > 0.) ? lambdaWI : lambdaWJ;
lambdaWIJ = (potentialDiffWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
lambdaWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
Scalar lambdaNwIJ = (potentialDiffNwIJ > 0.) ? lambdaNwI : lambdaNwJ;
lambdaNwIJ = (potentialDiffNwIJ == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
lambdaNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
if (compressibility_)
{
densityWIJ = (potentialDiffWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
densityNwIJ = (potentialDiffNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
densityWIJ = (potentialDiffWIJ == 0) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (potentialDiffNwIJ == 0) ? rhoMeanNwIJ : densityNwIJ;
densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
densityWIK = (potentialDiffWIK == 0) ? rhoMeanWIK : densityWIK;
densityNwIK = (potentialDiffNwIK == 0) ? rhoMeanNwIK : densityNwIK;
densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
}
//calculate velocities and the gravity term
......@@ -503,10 +504,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
density_[wPhaseIdx];
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
density_[nPhaseIdx];
potentialDiffW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx));
......@@ -526,9 +527,9 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30 )) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwJ;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
if (compressibility_)
{
......@@ -536,10 +537,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters
density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
density_[wPhaseIdx] =
(potentialDiffW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
density_[wPhaseIdx];
density_[nPhaseIdx] =
(potentialDiffNw == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
(Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
density_[nPhaseIdx];
}
......
......@@ -18,6 +18,8 @@
#ifndef DUMUX_MPFAL2DPRESSUREVELOCITY2P_HH
#define DUMUX_MPFAL2DPRESSUREVELOCITY2P_HH
#include <dune/common/float_cmp.hh>
#include "fvmpfal2dpressure2p.hh"
#include "fvmpfal2dvelocity2p.hh"
......@@ -498,9 +500,9 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
Scalar scalarPerm = permeability.two_norm();
......
......@@ -18,6 +18,8 @@
#ifndef DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH
#define DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH
#include <dune/common/float_cmp.hh>
#include "fvmpfal2dpressure2padaptive.hh"
#include "fvmpfal2dvelocity2padaptive.hh"
......@@ -635,9 +637,9 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
Scalar scalarPerm = permeability.two_norm();
......
......@@ -19,6 +19,7 @@
#ifndef DUMUX_FVMPFAL2PFABOUND3DVELOCITIES2P_HH
#define DUMUX_FVMPFAL2PFABOUND3DVELOCITIES2P_HH
#include <dune/common/float_cmp.hh>
#include "fvmpfal3dpressure2p.hh"
#include "fvmpfal3dvelocity2p.hh"
......@@ -496,9 +497,9 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
Scalar scalarPerm = permeability.two_norm();
......
......@@ -18,6 +18,7 @@
#ifndef DUMUX_FVMPFALPRESSUREVELOCITY2P_ADAPTIVE_HH
#define DUMUX_FVMPFALPRESSUREVELOCITY2P_ADAPTIVE_HH
#include <dune/common/float_cmp.hh>
#include "fvmpfal3dpressure2padaptive.hh"
#include "fvmpfal3dvelocity2padaptive.hh"
......@@ -598,9 +599,9 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
Scalar scalarPerm = permeability.two_norm();
......
......@@ -18,6 +18,7 @@
#ifndef DUMUX_MPFAO2DPRESSUREVELOCITIES2P_HH
#define DUMUX_MPFAO2DPRESSUREVELOCITIES2P_HH
#include <dune/common/float_cmp.hh>
#include "fvmpfao2dpressure2p.hh"
#include "fvmpfao2dvelocity2p.hh"
......@@ -499,9 +500,9 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
//do the upwinding of the mobility depending on the phase potentials
Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
lambdaW = (potentialDiffW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
lambdaNw = (potentialDiffNw == 0) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
Scalar scalarPerm = permeability.two_norm();
......
......@@ -24,6 +24,7 @@
* @brief Cfl-flux-function to evaluate a Cfl-Condition after Coats 2003
*/
#include <dune/common/float_cmp.hh>
#include <dumux/decoupled/common/impetproperties.hh>
#include "evalcflflux.hh"
......@@ -575,13 +576,13 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
bcValues[nPhaseIdx] *= faceArea;
bool hasPotWBound = false;
if (lambdaW != 0 && bcValues[wPhaseIdx] != 0)
if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(lambdaW, 0.0, 1.0e-30) && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(bcValues[wPhaseIdx], 0.0, 1.0e-30))
{
potWBound -= bcValues[wPhaseIdx] / (transmissibility * lambdaW);
hasPotWBound = true;
}
bool hasPotNwBound = false;
if (lambdaNw != 0 && bcValues[nPhaseIdx] != 0)
if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(lambdaNw, 0.0, 1.0e-30) && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(bcValues[nPhaseIdx], 0.0, 1.0e-30))
{
potNwBound -= bcValues[nPhaseIdx] / (transmissibility * lambdaNw);
hasPotNwBound = true;
......
......@@ -19,6 +19,8 @@
#ifndef DUMUX_ELEMENTDATA2P2C_ADAPTIVE_HH
#define DUMUX_ELEMENTDATA2P2C_ADAPTIVE_HH
#include <dune/common/float_cmp.hh>
#include <dumux/decoupled/2p2c/celldata2p2c.hh>
#include <dumux/decoupled/2p2c/celldata2p2cmultiphysics.hh>
/**
......@@ -154,7 +156,7 @@ public:
+= adaptedValues.cellVolume* adaptedValues.totalConcentration_[nCompIdx];
// if all cells are summed up, re-convert mass into total concentrations
Scalar fatherVolume = fatherElement.geometry().volume();
if(adaptedValuesFather.count == (pow(2,dim)))
if(adaptedValuesFather.count == 1 << dim)
{
adaptedValuesFather.totalConcentration_[wCompIdx] /= fatherVolume;
adaptedValuesFather.totalConcentration_[nCompIdx] /= fatherVolume;
......
......@@ -20,6 +20,7 @@
#define DUMUX_FV2DTRANSPORT2P2C_ADAPTIVE_HH
#include <dune/grid/common/gridenums.hh>
#include <dune/common/float_cmp.hh>
#include <dumux/common/math.hh>
#include <dumux/linear/vectorexchange.hh>
......@@ -537,10 +538,10 @@ void FV2dTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2
{
//check if harmonic weighting is necessary
// check if outflow induce neglected (i.e. mob=0) phase flux
if (potential[phaseIdx] > 0. && cellDataJ.mobility(phaseIdx) != 0.)
if (potential[phaseIdx] > 0. && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataJ.mobility(phaseIdx), 0.0, 1.0e-30))
lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
// check if inflow induce neglected phase flux
else if (potential[phaseIdx] < 0. && cellDataI.mobility(phaseIdx) != 0.)
else if (potential[phaseIdx] < 0. && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataI.mobility(phaseIdx), 0.0, 1.0e-30))
lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
else
doUpwinding[phaseIdx] = false;
......