From 23bc012a8055bb47f4394497507bdb6215f65097 Mon Sep 17 00:00:00 2001
From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>
Date: Wed, 28 Jan 2015 12:34:26 +0000
Subject: [PATCH] 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
---
 dumux/common/math.hh                          |  2 +-
 dumux/common/splinecommon_.hh                 |  3 +-
 dumux/common/timemanager.hh                   |  3 +-
 .../decoupled/2p/diffusion/fv/fvpressure2p.hh | 26 ++++++------
 .../2p/diffusion/fv/fvpressure2padaptive.hh   | 14 ++++---
 .../decoupled/2p/diffusion/fv/fvvelocity2p.hh | 25 +++++------
 .../2p/diffusion/fv/fvvelocity2padaptive.hh   | 41 ++++++++++---------
 .../lmethod/fvmpfal2dpressurevelocity2p.hh    |  6 ++-
 .../fvmpfal2dpressurevelocity2padaptive.hh    |  6 ++-
 .../lmethod/fvmpfal3dpressurevelocity2p.hh    |  5 ++-
 .../fvmpfal3dpressurevelocity2padaptive.hh    |  5 ++-
 .../omethod/fvmpfao2dpressurevelocity2p.hh    |  5 ++-
 .../2p/transport/fv/evalcflfluxcoats.hh       |  5 ++-
 dumux/decoupled/2p2c/celldata2p2cadaptive.hh  |  4 +-
 .../2p2c/fv2dtransport2p2cadaptive.hh         |  5 ++-
 .../2p2c/fv3dtransport2p2cadaptive.hh         |  5 ++-
 dumux/decoupled/2p2c/fvpressure2p2c.hh        | 17 ++++----
 .../2p2c/fvpressure2p2cmultiphysics.hh        | 17 ++++----
 .../decoupled/2p2c/fvpressurecompositional.hh |  3 +-
 dumux/decoupled/2p2c/fvtransport2p2c.hh       |  6 ++-
 dumux/decoupled/common/impetproblem.hh        |  4 +-
 dumux/freeflow/stokes/stokeslocalresidual.hh  | 13 +++---
 dumux/implicit/2pdfm/2pdfmfluxvariables.hh    |  4 +-
 .../common/implicitdarcyfluxvariables.hh      |  3 +-
 dumux/implicit/mpnc/diffusion/diffusion.hh    |  4 +-
 dumux/linear/amgparallelhelpers.hh            |  4 +-
 dumux/material/components/co2tablereader.hh   |  5 ++-
 .../2cnistokes2p2cniproblem.hh                |  3 +-
 .../2cnistokes2p2cni/2p2cnisubproblem.hh      |  6 ++-
 .../2cstokes2p2c/2cstokes2p2cproblem.hh       |  3 +-
 .../2cstokes2p2c/2p2csubproblem.hh            |  6 ++-
 31 files changed, 150 insertions(+), 108 deletions(-)

diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index 56981485ed..45f703c9f2 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -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.
  *
diff --git a/dumux/common/splinecommon_.hh b/dumux/common/splinecommon_.hh
index dafa1bd0d2..4dab68c951 100644
--- a/dumux/common/splinecommon_.hh
+++ b/dumux/common/splinecommon_.hh
@@ -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)
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index 1927c442fb..3667e2bd4c 100644
--- a/dumux/common/timemanager.hh
+++ b/dumux/common/timemanager.hh
@@ -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());
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index 1646fed9dd..fef5b51e0f 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
index aabcab61a0..db4aa45a3f 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
@@ -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
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index 4da38ca88e..95caa1c972 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
index 1a4d1cdf8a..594c1af858 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
@@ -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];
         }
 
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
index c89587200c..2552f39ebe 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2p.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
index a82e278b65..a108c812da 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal2dpressurevelocity2padaptive.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
index ab3ff1fb2b..e5d455a620 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2p.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
index 60029247dc..d51618404f 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/lmethod/fvmpfal3dpressurevelocity2padaptive.hh
@@ -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();
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
index 7642da0776..45aa6303c9 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/omethod/fvmpfao2dpressurevelocity2p.hh
@@ -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();
diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
index ed670d7e92..cf0e3733b2 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
@@ -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;
diff --git a/dumux/decoupled/2p2c/celldata2p2cadaptive.hh b/dumux/decoupled/2p2c/celldata2p2cadaptive.hh
index 048650b773..a1232ee81e 100644
--- a/dumux/decoupled/2p2c/celldata2p2cadaptive.hh
+++ b/dumux/decoupled/2p2c/celldata2p2cadaptive.hh
@@ -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;
diff --git a/dumux/decoupled/2p2c/fv2dtransport2p2cadaptive.hh b/dumux/decoupled/2p2c/fv2dtransport2p2cadaptive.hh
index 703170e81d..2b6bb69961 100644
--- a/dumux/decoupled/2p2c/fv2dtransport2p2cadaptive.hh
+++ b/dumux/decoupled/2p2c/fv2dtransport2p2cadaptive.hh
@@ -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;
diff --git a/dumux/decoupled/2p2c/fv3dtransport2p2cadaptive.hh b/dumux/decoupled/2p2c/fv3dtransport2p2cadaptive.hh
index cfd652aa84..a9645ef941 100644
--- a/dumux/decoupled/2p2c/fv3dtransport2p2cadaptive.hh
+++ b/dumux/decoupled/2p2c/fv3dtransport2p2cadaptive.hh
@@ -20,6 +20,7 @@
 #define DUMUX_FV3DTRANSPORT2P2C_ADAPTIVE_HH
 
 #include <dune/grid/common/gridenums.hh>
+#include <dune/common/float_cmp.hh>
 
 #include <dumux/common/math.hh>
 #include <dumux/linear/vectorexchange.hh>
@@ -543,9 +544,9 @@ void FV3dTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2
             else    // i.e. restrictFluxInTransport == 1
             {
                //check if harmonic weithing is necessary
-                if (potential[phaseIdx] > 0. && cellDataJ.mobility(phaseIdx) != 0.) // check if outflow induce neglected (i.e. mob=0) phase flux
+                if (potential[phaseIdx] > 0. && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataJ.mobility(phaseIdx), 0.0, 1.0e-30)) // check if outflow induce neglected (i.e. mob=0) phase flux
                     lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
-                else if (potential[phaseIdx] < 0. && cellDataI.mobility(phaseIdx) != 0.) // check if inflow induce neglected phase flux
+                else if (potential[phaseIdx] < 0. && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataI.mobility(phaseIdx), 0.0, 1.0e-30)) // check if inflow induce neglected phase flux
                     lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
                 else
                     doUpwinding[phaseIdx] = false;
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 124eb18f22..4b9610372e 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -24,6 +24,7 @@
 #include <dune/istl/operators.hh>
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
+#include <dune/common/float_cmp.hh>
 
 // dumux environment
 #include <dumux/decoupled/2p2c/fvpressurecompositional.hh>
@@ -268,7 +269,7 @@ void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEn
     Scalar timestep_ = problem().timeManager().timeStepSize();
 
     // compressibility term
-    if (!first && timestep_ != 0)
+    if (!first && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(timestep_, 0.0, 1.0e-30))
     {
         Scalar compress_term = cellDataI.dv_dp() / timestep_;
 
@@ -783,11 +784,11 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
 
             if (potentialW >= 0.)
             {
-                densityW = (potentialW == 0) ? rhoMeanW : cellDataI.density(wPhaseIdx);
+                densityW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialW, 0.0, 1.0e-30)) ? rhoMeanW : cellDataI.density(wPhaseIdx);
                 dV_w = (dv_dC1 * cellDataI.massFraction(wPhaseIdx, wCompIdx)
                                    + dv_dC2 * cellDataI.massFraction(wPhaseIdx, nCompIdx));
                 dV_w *= densityW;
-                lambdaW = (potentialW == 0) ? 0.5 * (cellDataI.mobility(wPhaseIdx) + lambdaWBound)
+                lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialW, 0.0, 1.0e-30)) ? 0.5 * (cellDataI.mobility(wPhaseIdx) + lambdaWBound)
                                             : cellDataI.mobility(wPhaseIdx);
             }
             else
@@ -800,11 +801,11 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
             }
             if (potentialNW >= 0.)
             {
-                densityNW = (potentialNW == 0) ? rhoMeanNW : cellDataI.density(nPhaseIdx);
+                densityNW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNW, 0.0, 1.0e-30)) ? rhoMeanNW : cellDataI.density(nPhaseIdx);
                 dV_n = (dv_dC1 * cellDataI.massFraction(nPhaseIdx, wCompIdx)
                         + dv_dC2 * cellDataI.massFraction(nPhaseIdx, nCompIdx));
                 dV_n *= densityNW;
-                lambdaNW = (potentialNW == 0) ? 0.5 * (cellDataI.mobility(nPhaseIdx) + lambdaNWBound)
+                lambdaNW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialNW, 0.0, 1.0e-30)) ? 0.5 * (cellDataI.mobility(nPhaseIdx) + lambdaNWBound)
                                               : cellDataI.mobility(nPhaseIdx);
             }
             else
@@ -910,7 +911,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
                     + cellData.massConcentration(nCompIdx));
 
     // make shure only physical quantities enter flash calculation
-    if(Z1<0. || Z1 > 1.)
+    if(Z1 < 0. || Z1 > 1.)
     {
         Dune::dgrave << "Feed mass fraction unphysical: Z1 = " << Z1
                << " at global Idx " << eIdxGlobal
@@ -918,7 +919,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
                << cellData.totalConcentration(wCompIdx)
                << " and totalConcentration(nCompIdx) = "
                << cellData.totalConcentration(nCompIdx)<< std::endl;
-        if(Z1<0.)
+        if(Z1 < 0.)
             {
             Z1 = 0.;
             cellData.setTotalConcentration(wCompIdx, 0.);
@@ -1026,7 +1027,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
     if ((cellData.density(wPhaseIdx)*cellData.density(nPhaseIdx)) == 0)
         DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
     Scalar vol = massw / cellData.density(wPhaseIdx) + massn / cellData.density(nPhaseIdx);
-    if (problem().timeManager().timeStepSize() != 0)
+    if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30))
     {
         cellData.volumeError()=(vol - problem().spatialParams().porosity(elementI));
 
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index 5b3793b274..2fa19766d2 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -19,6 +19,8 @@
 #ifndef DUMUX_FVPRESSURE2P2C_MULTIPHYSICS_HH
 #define DUMUX_FVPRESSURE2P2C_MULTIPHYSICS_HH
 
+#include <dune/common/float_cmp.hh>
+
 // dumux environment
 #include <dumux/decoupled/2p2c/fvpressure2p2c.hh>
 #include <dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh>
@@ -724,8 +726,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFluxOnBoundary(Dune::FieldVector<
 
                         if (potential >= 0.)
                         {
-                            density = (potential == 0) ? rhoMean : cellDataI.density(phaseIdx);
-                            lambda = (potential == 0) ? 0.5 * (lambdaI + lambdaBound) : lambdaI;
+                            density = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potential, 0.0, 1.0e-30)) ? rhoMean : cellDataI.density(phaseIdx);
+                            lambda = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potential, 0.0, 1.0e-30)) ? 0.5 * (lambdaI + lambdaBound) : lambdaI;
                         }
                         else
                         {
@@ -784,7 +786,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
     Scalar maxError = 0.;
 
     // next subdomain map
-    if (problem().timeManager().time() == 0.)
+    if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().time(), 0.0, 1.0e-30))
         nextSubdomain = 2;  // start with complicated sub in initialization
     else
         nextSubdomain = -1;  // reduce complexity after first TS
@@ -808,7 +810,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
             PrimaryVariables source(NAN);
             problem().source(source, *eIt);
 
-            if (cellData.saturation(wPhaseIdx) != (1. || 0.) or source.one_norm()!= 0.) // cell still 2p
+            if ((cellData.saturation(wPhaseIdx) > 0.0 && cellData.saturation(wPhaseIdx) < 1.0)
+            		|| Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(source.one_norm(), 0.0, 1.0e-30)) // cell still 2p
             {
                 // mark this element
                 nextSubdomain[eIdxGlobal] = 2;
@@ -827,9 +830,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
             }
             else if(nextSubdomain[eIdxGlobal] != 2)// update next subdomain if possible
             {
-                if(cellData.saturation(wPhaseIdx) != 0.)
+                if(Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(wPhaseIdx), 0.0, 1.0e-30))
                     nextSubdomain[eIdxGlobal] = wPhaseIdx;
-                else if (cellData.saturation(nPhaseIdx) != 0.)
+                else if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(nPhaseIdx), 0.0, 1.0e-30))
                     nextSubdomain[eIdxGlobal] = nPhaseIdx;
             }
             timer_.stop();
@@ -974,7 +977,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const El
     Scalar vol(0.);
     vol = sumConc / pseudoFluidState.density(presentPhaseIdx);
 
-    if (problem().timeManager().timeStepSize() != 0)
+    if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30))
         cellData.volumeError() = (vol - problem().spatialParams().porosity(elementI));
     return;
 }
diff --git a/dumux/decoupled/2p2c/fvpressurecompositional.hh b/dumux/decoupled/2p2c/fvpressurecompositional.hh
index 9105dee35c..ace74fc50d 100644
--- a/dumux/decoupled/2p2c/fvpressurecompositional.hh
+++ b/dumux/decoupled/2p2c/fvpressurecompositional.hh
@@ -20,6 +20,7 @@
 #define DUMUX_FVPRESSURECOMPOSITIONAL_HH
 
 #include <cmath>
+#include <dune/common/float_cmp.hh>
 
 // dumux environment
 #include <dumux/common/math.hh>
@@ -137,7 +138,7 @@ public:
 
         // if we just started a new episode, the TS size of the update Estimate is a better
         // estimate then the size of the last time step
-        if(problem_.timeManager().time() == problem_.timeManager().episodeStartTime()
+        if(Dune::FloatCmp::eq<Scalar>(problem_.timeManager().time(), problem_.timeManager().episodeStartTime())
                 && problem_.timeManager().episodeIndex() > 1)
             problem_.timeManager().setTimeStepSize(dt_estimate*GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, CFLFactor));
 
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index 343442da83..24b8e1b3e7 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -23,6 +23,8 @@
 #include <unordered_map>
 
 #include <dune/grid/common/gridenums.hh>
+#include <dune/common/float_cmp.hh>
+
 #include <dumux/decoupled/2p2c/2p2cproperties.hh>
 #include <dumux/material/constraintsolvers/compositionalflash.hh>
 #include <dumux/common/math.hh>
@@ -723,10 +725,10 @@ void FVTransport2P2C<TypeTag>::getFlux(ComponentVector& fluxEntries,
             else    // i.e. restrictFluxInTransport == 1
             {
                //check if harmonic weighting is necessary
-                if (potential[phaseIdx] > 0. && (cellDataJ.mobility(phaseIdx) != 0.   // check if outflow induce neglected (i.e. mob=0) phase flux
+                if (potential[phaseIdx] > 0. && (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataJ.mobility(phaseIdx), 0.0, 1.0e-30)   // check if outflow induce neglected (i.e. mob=0) phase flux
                        or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
                     lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
-                else if (potential[phaseIdx] < 0. && (cellDataI.mobility(phaseIdx) != 0. // check if inflow induce neglected phase flux
+                else if (potential[phaseIdx] < 0. && (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellDataI.mobility(phaseIdx), 0.0, 1.0e-30) // check if inflow induce neglected phase flux
                         or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
                     lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
                 else
diff --git a/dumux/decoupled/common/impetproblem.hh b/dumux/decoupled/common/impetproblem.hh
index 09487d2a70..8c617c1728 100644
--- a/dumux/decoupled/common/impetproblem.hh
+++ b/dumux/decoupled/common/impetproblem.hh
@@ -19,6 +19,7 @@
 #ifndef DUMUX_IMPETPROBLEM_HH
 #define DUMUX_IMPETPROBLEM_HH
 
+#include <dune/common/float_cmp.hh>
 #include "impetproperties.hh"
 #include <dumux/io/vtkmultiwriter.hh>
 #include <dumux/io/restart.hh>
@@ -371,7 +372,8 @@ public:
         dt = std::min(dt, timeManager().episodeMaxTimeStepSize());
 
         // check if we are in first TS and an initialDt was assigned
-        if (t==0. && timeManager().timeStepSize()!=0.)
+        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(t, 0.0, 1.0e-30) 
+        	&& Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(timeManager().timeStepSize(), 0.0, 1.0e-30))
         {
             if (this->gridView().comm().size() > 1)
                 dt = this->gridView().comm().min(dt);
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index 55e10312c2..4be5661c97 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -26,6 +26,7 @@
 #ifndef DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
 #define DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
 
+#include <dune/common/float_cmp.hh>
 #include <dune/grid/common/grid.hh>
 
 #include <dumux/implicit/common/implicitmodel.hh>
@@ -190,7 +191,7 @@ protected:
         if(!useMoles)
         {
             // mass balance with upwinded density
-            if (massUpwindWeight_ == 1.0) // fully upwind
+            if (Dune::FloatCmp::eq<Scalar>(massUpwindWeight_, 1.0)) // fully upwind
                 massBalanceResidual *= up.density();
             else
                 massBalanceResidual *= massUpwindWeight_ * up.density() +
@@ -199,7 +200,7 @@ protected:
         else
         {
             // mass balance with upwinded molarDensity
-            if (massUpwindWeight_ == 1.0) // fully upwind
+            if (Dune::FloatCmp::eq<Scalar>(massUpwindWeight_, 1.0)) // fully upwind
                 massBalanceResidual *= up.molarDensity();
             else
                 massBalanceResidual *= massUpwindWeight_ * up.molarDensity() +
@@ -493,7 +494,7 @@ protected:
                 // calculate  mu grad v t t
                 // center in the face of the reference element
                 DimVector tangent;
-                if (stabilizationBeta_ != 0)
+                if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(stabilizationBeta_, 0.0, 1.0e-30))
                 {
                     const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
                     tangent[0] = elementUnitNormal[1];  //TODO: 3D
@@ -564,7 +565,7 @@ protected:
 
             // beta-stabilization at the boundary in case of outflow condition
             // for the momentum balance
-            if(momentumBalanceOutflow_(bcTypes) && stabilizationBeta_ != 0)
+            if(momentumBalanceOutflow_(bcTypes) && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(stabilizationBeta_, 0.0, 1.0e-30))
             {
                 // calculate  mu grad v t t for beta-stabilization
                 // center in the face of the reference element
@@ -603,7 +604,7 @@ protected:
      */
     void removeStabilizationAtBoundary_(const int scvIdx)
     {
-        if (stabilizationAlpha_ != 0)
+        if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(stabilizationAlpha_, 0.0, 1.0e-30))
         {
             // loop over the edges of the element
             for (int fIdx = 0; fIdx < this->fvGeometry_().numScvf; fIdx++)
@@ -679,7 +680,7 @@ protected:
                                      DimVector& averagedNormal,
                                      const int scvIdx)
     {
-        assert(averagedNormal.two_norm() != 0.0);
+        assert( (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(averagedNormal.two_norm(), 0.0, 1.0e-30)) );
 
         // divide averagedNormal by its length
         averagedNormal /= averagedNormal.two_norm();
diff --git a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
index 9bd7cbcf48..fa38e3b148 100644
--- a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
@@ -23,6 +23,8 @@
 #ifndef DUMUX_MODELS_2PDFM_FLUX_VARIABLES_HH
 #define DUMUX_MODELS_2PDFM_FLUX_VARIABLES_HH
 
+#include <dune/common/float_cmp.hh>
+
 #include <dumux/common/math.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/valgrind.hh>
@@ -229,7 +231,7 @@ private:
                         Scalar rhoJ = elemVolVars[this->face().j].fluidState().density(phaseIdx);
                         Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
                         Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
-                        if (fI + fJ == 0)
+                        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
                             // doesn't matter because no wetting phase is present in
                             // both cells!
                             fI = fJ = 0.5;
diff --git a/dumux/implicit/common/implicitdarcyfluxvariables.hh b/dumux/implicit/common/implicitdarcyfluxvariables.hh
index 31196ec5a7..f31687db6f 100644
--- a/dumux/implicit/common/implicitdarcyfluxvariables.hh
+++ b/dumux/implicit/common/implicitdarcyfluxvariables.hh
@@ -27,6 +27,7 @@
 #define DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
 
 #include "implicitproperties.hh"
+#include <dune/common/float_cmp.hh>
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/math.hh>
@@ -210,7 +211,7 @@ protected:
                 Scalar rhoJ = elemVolVars[face().j].fluidState().density(phaseIdx);
                 Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
                 Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
-                if (fI + fJ == 0)
+                if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
                     // doesn't matter because no wetting phase is present in
                     // both cells!
                     fI = fJ = 0.5;
diff --git a/dumux/implicit/mpnc/diffusion/diffusion.hh b/dumux/implicit/mpnc/diffusion/diffusion.hh
index 808d12e654..413f2c45f2 100644
--- a/dumux/implicit/mpnc/diffusion/diffusion.hh
+++ b/dumux/implicit/mpnc/diffusion/diffusion.hh
@@ -25,9 +25,11 @@
 #ifndef DUMUX_MPNC_DIFFUSION_HH
 #define DUMUX_MPNC_DIFFUSION_HH
 
+#include <dune/common/float_cmp.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
+
 #include <dumux/implicit/mpnc/mpncproperties.hh>
 
 namespace Dumux {
@@ -130,7 +132,7 @@ protected:
             for (int compIIdx = 0; compIIdx < numComponents - 1; ++compIIdx) {
                 for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
                     Scalar Dij = fluxVars.porousDiffCoeffG(compIIdx, compJIdx);
-                    if (Dij) {
+                    if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(Dij, 0.0, 1.0e-30)) {
                         M[compIIdx][compJIdx] += fluxVars.moleFraction(nPhaseIdx, compIIdx) / Dij;
                         M[compIIdx][compIIdx] -= fluxVars.moleFraction(nPhaseIdx, compJIdx) / Dij;
                     }
diff --git a/dumux/linear/amgparallelhelpers.hh b/dumux/linear/amgparallelhelpers.hh
index 6c41c10910..cfed547a8c 100644
--- a/dumux/linear/amgparallelhelpers.hh
+++ b/dumux/linear/amgparallelhelpers.hh
@@ -932,7 +932,7 @@ void ParallelISTLHelper<TypeTag>::createIndexSetAndProjectForAMG(M& m, C& c)
     auto owned=owner_.begin();
     
     for(auto v=sharedDofs.begin(), vend=sharedDofs.end(); v != vend; ++v, ++owned)
-        if(*v && *owned==1.0)
+        if(*v && *owned==1)
             ++count;
 
     Dune::dverb<<gridview.comm().rank()<<": shared count is "<< count.touint()
@@ -956,7 +956,7 @@ void ParallelISTLHelper<TypeTag>::createIndexSetAndProjectForAMG(M& m, C& c)
     auto index=scalarIndices.begin();
 
     for(auto i=owner_.begin(), iend=owner_.end(); i!=iend; ++i, ++shared, ++index)
-        if(*i==1.0 && *shared){
+        if(*i==1 && *shared){
             *index=start;
             ++start;
         }
diff --git a/dumux/material/components/co2tablereader.hh b/dumux/material/components/co2tablereader.hh
index 30a7693c95..e19ace9398 100644
--- a/dumux/material/components/co2tablereader.hh
+++ b/dumux/material/components/co2tablereader.hh
@@ -25,6 +25,7 @@
 #ifndef DUMUX_TABULATED_CO2_HH
 #define DUMUX_TABULATED_CO2_HH
 
+#include <dune/common/float_cmp.hh>
 #include <dumux/common/exceptions.hh>
 
 #include <assert.h>
@@ -119,7 +120,7 @@ public:
 protected:
     int findTempIdx_(Scalar temperature) const
     {
-        if (temperature == maxTemp())
+        if (Dune::FloatCmp::eq<Scalar>(temperature, maxTemp()))
             return numTempSteps - 2;
         const int result = static_cast<int>((temperature - minTemp())/(maxTemp() - minTemp())*(numTempSteps - 1));
         return std::max(0, std::min(result, numTempSteps - 2));
@@ -127,7 +128,7 @@ protected:
 
     int findPressIdx_(Scalar pressure) const
     {
-        if (pressure == maxPress())
+        if (Dune::FloatCmp::eq<Scalar>(pressure, maxPress()))
             return numPressSteps - 2;
         const int result = static_cast<int>((pressure - minPress())/(maxPress() - minPress())*(numPressSteps - 1));
         return std::max(0, std::min(result, numPressSteps - 2));
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
index 33a126d79b..4c6bdfae1f 100644
--- a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
@@ -32,6 +32,7 @@
 #ifndef DUMUX_2CNISTOKES2P2CNIPROBLEM_HH
 #define DUMUX_2CNISTOKES2P2CNIPROBLEM_HH
 
+#include <dune/common/float_cmp.hh>
 #include <dune/grid/common/gridinfo.hh>
 #include <dune/grid/multidomaingrid.hh>
 #include <dune/grid/io/file/dgfparser.hh>
@@ -262,7 +263,7 @@ public:
         ParentType::init();
 
         std::cout << "Writing flux data at interface\n";
-        if (this->timeManager().time() == 0)
+        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(this->timeManager().time(), 0.0, 1.0e-30))
         {
             fluxFile_.open("fluxes.out");
             fluxFile_ << "Time;"
diff --git a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
index 5f03262e77..5556e8bd6d 100644
--- a/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2p2cnisubproblem.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_2P2CNISUB_PROBLEM_HH
 #define DUMUX_2P2CNISUB_PROBLEM_HH
 
+#include <dune/common/float_cmp.hh>
+
 #include <dumux/implicit/common/implicitporousmediaproblem.hh>
 #include <dumux/implicit/2p2c/2p2cmodel.hh>
 #include <dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh>
@@ -367,7 +369,7 @@ public:
                 PrimaryVariables storageChange(0.);
                 storageChange = storageLastTimestep_ - storage;
 
-                assert(time - lastMassOutputTime_ != 0);
+                assert( (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(time - lastMassOutputTime_, 0.0, 1.0e-30)) );
                 storageChange /= (time - lastMassOutputTime_);
                 // 2d: interface length has to be accounted for
                 // in order to obtain kg/m²s
@@ -379,7 +381,7 @@ public:
                           << " IntEnergy: " << storage[energyEqIdx]
                           << " WaterMassChange: " << storageChange[contiWEqIdx]
                           << std::endl;
-                if (this->timeManager().time() != 0.)
+                if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(this->timeManager().time(), 0.0, 1.0e-30))
                     outfile << time << ";"
                             << storageChange[contiTotalMassIdx] << ";"
                             << storageChange[contiWEqIdx] << ";"
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
index 8730ad1ba1..6039597b56 100644
--- a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
@@ -32,6 +32,7 @@
 #ifndef DUMUX_2CSTOKES2P2CPROBLEM_HH
 #define DUMUX_2CSTOKES2P2CPROBLEM_HH
 
+#include <dune/common/float_cmp.hh>
 #include <dune/grid/common/gridinfo.hh>
 #include <dune/grid/multidomaingrid.hh>
 #include <dune/grid/io/file/dgfparser.hh>
@@ -255,7 +256,7 @@ public:
         ParentType::init();
 
         std::cout << "Writing flux data at interface\n";
-        if (this->timeManager().time() == 0)
+        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(this->timeManager().time(), 0.0, 1.0e-30))
         {
             fluxFile_.open("fluxes.out");
             fluxFile_ << "Time;"
diff --git a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
index a920e71fc1..ca8b208150 100644
--- a/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2p2csubproblem.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_2P2C_SUBPROBLEM_HH
 #define DUMUX_2P2C_SUBPROBLEM_HH
 
+#include <dune/common/float_cmp.hh>
+
 #include <dumux/implicit/2p2c/2p2cindices.hh>
 #include <dumux/implicit/common/implicitporousmediaproblem.hh>
 #include <dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh>
@@ -347,7 +349,7 @@ public:
                 PrimaryVariables storageChange(0.);
                 storageChange = storageLastTimestep_ - storage;
 
-                assert(time - lastMassOutputTime_ != 0);
+                assert( (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(time - lastMassOutputTime_, 0.0, 1.0e-30)) );
                 storageChange /= (time - lastMassOutputTime_);
                 // 2d: interface length has to be accounted for
                 // in order to obtain kg/m²s
@@ -358,7 +360,7 @@ public:
                           << " WaterMass: " << storage[contiWEqIdx]
                           << " WaterMassChange: " << storageChange[contiWEqIdx]
                           << std::endl;
-                if (this->timeManager().time() != 0.)
+                if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(this->timeManager().time(), 0.0, 1.0e-30))
                     outfile << time << ";"
                             << storageChange[contiTotalMassIdx] << ";"
                             << storageChange[contiWEqIdx] << ";"
-- 
GitLab