From 26eb9f2bd8db4f43ffd810e882b277b3e817e880 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Thu, 1 Sep 2022 13:42:56 +0200 Subject: [PATCH 1/7] [pnm][pnm-regularization] Use the upwinded pore to determine conducitivty reg interval when reg with satruation --- dumux/porenetwork/2p/fluxvariablescache.hh | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index 56c52d6b08..23b3f95d98 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -70,10 +70,18 @@ public: pc_ = std::max(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); + auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, outsideScv, elemSol); + if ( (!invaded && pcInsidePore > pcOutsidePore) || (invaded && pcInsidePore < pcOutsidePore) ) + 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; -- GitLab From d58d52fa74d45c00c3021ad6a2cd1429abdaeea4 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Fri, 2 Sep 2022 11:15:57 +0200 Subject: [PATCH 2/7] [pnm] Use minimum of two adjacent pore combined with histroy invasion state to decidede if snapoff happens --- dumux/porenetwork/2p/fluxvariablescache.hh | 5 +++++ dumux/porenetwork/2p/invasionstate.hh | 6 +++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index 23b3f95d98..f2daaa7882 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -82,6 +82,11 @@ public: 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()); + // 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; diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh index 6322574eeb..f8c1967e07 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 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) -- GitLab From 99cb1491fa1e39d6a06f371f2644c20bc6451987 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Mon, 12 Sep 2022 20:00:15 +0200 Subject: [PATCH 3/7] [pnm] Use averaged pc to calculate the avg radius/wetting layer and conductivity --- dumux/porenetwork/2p/fluxvariablescache.hh | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index f2daaa7882..0092657430 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -68,7 +68,8 @@ 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()]; @@ -251,7 +252,7 @@ public: * \brief Returns the curvature radius within the throat. */ Scalar curvatureRadius() const - { return surfaceTension_ / pc_;} + { return surfaceTension_ / pcAvg_;} /*! * \brief Returns the cross-sectional area of a wetting layer within @@ -379,6 +380,8 @@ private: Scalar throatInscribedRadius_; Scalar pcEntry_; Scalar pcSnapoff_; + Scalar pcMax_; + Scalar pcAvg_; Scalar pc_; Scalar surfaceTension_; bool invaded_; -- GitLab From ff8755a0342ad7af44594507fec4fb82031f4605 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 29 Nov 2022 09:24:11 +0100 Subject: [PATCH 4/7] [pnm] reduce output and use pc instead of pcAvg to calculate k --- dumux/porenetwork/2p/fluxvariablescache.hh | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/dumux/porenetwork/2p/fluxvariablescache.hh b/dumux/porenetwork/2p/fluxvariablescache.hh index 0092657430..b46588f930 100644 --- a/dumux/porenetwork/2p/fluxvariablescache.hh +++ b/dumux/porenetwork/2p/fluxvariablescache.hh @@ -88,10 +88,6 @@ public: else pc_ = std::min(elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure()); - // 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; - regInvasionInterval_[0] = pcEntry_; regSnapoffInterval_[1] = pcSnapoff_; @@ -252,7 +248,7 @@ public: * \brief Returns the curvature radius within the throat. */ Scalar curvatureRadius() const - { return surfaceTension_ / pcAvg_;} + { return surfaceTension_ / pc_;} /*! * \brief Returns the cross-sectional area of a wetting layer within -- GitLab From 752ffb99cf8beb173ad74993526cbf680c862064 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 29 Nov 2022 09:26:59 +0100 Subject: [PATCH 5/7] [boxassembler][pnm] Update FluxVarsCache --- dumux/assembly/boxlocalassembler.hh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh index 6cc4d0f712..dc03466bde 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(); }; -- GitLab From cfd09d558fd77c549db0e53625a4d0481021a576 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 29 Nov 2022 09:30:17 +0100 Subject: [PATCH 6/7] [pnm][debug] Update invasion state only after converged --- dumux/nonlinear/newtonsolver.hh | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 9ef27335cb..b4717a7010 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", -- GitLab From 79ae66bc105c78523b17f383a76a77760a7f1cf2 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 29 Nov 2022 09:34:11 +0100 Subject: [PATCH 7/7] [pnm] Better regularize k with Sw --- test/porenetwork/2p/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/porenetwork/2p/CMakeLists.txt b/test/porenetwork/2p/CMakeLists.txt index c02a388b1c..b0be575a0a 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 -- GitLab