diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index 6cc4d0f7129c8546be49cb8efaa542d94ad74536..dc03466bdef197f5a33ed534184b94058d9c5add 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -368,6 +368,7 @@ public:
         const auto& fvGeometry = this->fvGeometry();
         const auto& curSol = this->curSol();
         auto&& curElemVolVars = this->curElemVolVars();
+        auto&& elemFluxVarsCache = this->elemFluxVarsCache();
 
         // get the vector of the actual element residuals
         const auto origResiduals = this->evalLocalResidual();
@@ -417,6 +418,7 @@ public:
                     // update the volume variables and compute element residual
                     elemSol[scv.localDofIndex()][pvIdx] = priVar;
                     deflectionHelper.deflect(elemSol, scv, this->problem());
+                    elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
                     return this->evalLocalResidual();
                 };
 
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 9ef27335cb4278331e5ec32483822054333488f3..b4717a70105a7c804839f4500fc126186e6d5866 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -659,6 +659,7 @@ public:
 
             if (enableShiftCriterion_)
                 std::cout << Fmt::format(", maximum relative shift = {:.4e}", shift_);
+
             if (enableResidualCriterion_ && enableAbsoluteResidualCriterion_)
                 std::cout << Fmt::format(", residual = {:.4e}", residualNorm_);
             else if (enableResidualCriterion_)
@@ -1030,9 +1031,6 @@ private:
                 converged = newtonConverged();
             }
 
-            // tell solver we are done
-            newtonEnd(vars, uLastIter);
-
             // reset state if Newton failed
             if (!newtonConverged())
             {
@@ -1047,6 +1045,9 @@ private:
             // tell solver we converged successfully
             newtonSucceed();
 
+            // tell solver we are done
+            newtonEnd(vars, uLastIter);
+
             if (verbosity_ >= 1) {
                 const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
                 std::cout << Fmt::format("Assemble/solve/update time: {:.2g}({:.2f}%)/{:.2g}({:.2f}%)/{:.2g}({:.2f}%)\n",
diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh
index 56c52d6b08c15063355286e06a2e51c924f087f5..b46588f930aedf58fecc38a30a599a5d76907b52 100644
--- a/dumux/porenetwork/2p/fluxvariablescache.hh
+++ b/dumux/porenetwork/2p/fluxvariablescache.hh
@@ -68,15 +68,25 @@ public:
         // take the average surface tension of both adjacent pores TODO: is this correct?
         surfaceTension_ = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension());
 
-        pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure());
+        pcMax_ = pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure());
+        pcAvg_ = 0.5 * (pc_ + std::min(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()));
 
+        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+        auto pcOutsidePore = outsideVolVars.capillaryPressure();
+        auto pcInsidePore = insideVolVars.capillaryPressure();
         const auto& elemSol = elementSolution(element, elemVolVars, fvGeometry);
         const auto& spatialParams = problem.spatialParams();
-        auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, insideScv, elemSol);
-        // open a file to write the analytical pcEntry and Sw for comparison purposes   
-        std::ofstream file("analytical_SwEntry");
-        file << std::setprecision(15) << "pcEntry: " << pcEntry_ << " Sw: " << fluidMatrixInteraction.sw(pcEntry_) << std::endl;
+        auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, outsideScv, elemSol);
+        if (  (!invaded &&  pcInsidePore > pcOutsidePore) || (invaded && pcInsidePore < pcOutsidePore) )
+            fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, insideScv, elemSol);
+
+        if (!invaded)
+            pc_ = std::max(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure());
+        else
+            pc_ = std::min(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure());
 
         regInvasionInterval_[0] = pcEntry_;
         regSnapoffInterval_[1] = pcSnapoff_;
@@ -366,6 +376,8 @@ private:
     Scalar throatInscribedRadius_;
     Scalar pcEntry_;
     Scalar pcSnapoff_;
+    Scalar pcMax_;
+    Scalar pcAvg_;
     Scalar pc_;
     Scalar surfaceTension_;
     bool invaded_;
diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh
index 6322574eebde91a86f7e08e646193a745bce94f2..f8c1967e07c23943b3cc631d661e2670dffafbaf 100644
--- a/dumux/porenetwork/2p/invasionstate.hh
+++ b/dumux/porenetwork/2p/invasionstate.hh
@@ -217,6 +217,7 @@ private:
         //Determine whether throat gets invaded or snap-off occurs
         const std::array<Scalar, 2> pc = { elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure() };
         const auto pcMax = std::max_element(pc.begin(), pc.end());
+        const auto pcMin = std::min_element(pc.begin(), pc.end());
         const Scalar pcEntry = fluxVarsCache.pcEntry();
         const Scalar pcSnapoff = fluxVarsCache.pcSnapoff();
 
@@ -233,11 +234,10 @@ private:
             return Result{}; //nothing happened
         }
 
-        if (*pcMax > pcEntry)
+        if (!invadedBeforeSwitch && *pcMax > pcEntry)
            invadedAfterSwitch = true;
-        else if (*pcMax <= pcSnapoff)
+        else if (invadedBeforeSwitch && *pcMin <= pcSnapoff)
            invadedAfterSwitch = false;
-
         invadedCurrentTimeStep_[eIdx] = invadedAfterSwitch;
 
         if (invadedBeforeSwitch == invadedAfterSwitch)
diff --git a/test/porenetwork/2p/CMakeLists.txt b/test/porenetwork/2p/CMakeLists.txt
index c02a388b1ccbae9da847d88f71afad92ef065172..b0be575a0a54541f478d50a11a48e4a99cf5b2b9 100644
--- a/test/porenetwork/2p/CMakeLists.txt
+++ b/test/porenetwork/2p/CMakeLists.txt
@@ -5,7 +5,7 @@ add_input_file_links()
 dumux_add_test(NAME test_pnm_2p
                SOURCES main.cc
                LABELS porenetwork
-               COMPILE_DEFINITIONS ISOTHERMAL=1 REGULARIZATIONWITHPRESSURE=1
+               COMPILE_DEFINITIONS ISOTHERMAL=1 REGULARIZATIONWITHPRESSURE=0
                CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
                TIMEOUT "1000"
                COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py