From 057166129705420f11df38879a8e6130714921ff Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Tue, 2 Mar 2021 18:07:11 +0100 Subject: [PATCH 01/15] [hack] Change free flow block structure to old version: cc face-coupling face cc-coupling --- dumux/nonlinear/newtonsolver.hh | 262 ++++++++++++++++++++++++++++++-- 1 file changed, 251 insertions(+), 11 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 61ccca1595..5b11b9cc70 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -42,6 +42,7 @@ #include #include +#include #include #include #include @@ -158,6 +159,117 @@ constexpr std::size_t blockSize() { return 1; } template>::value, int> = 0> constexpr std::size_t blockSize() { return std::decay_t::size(); } + +//! This is a partial copy of Dune::MultiTypeBlockMatrix. TODO: It can be deleted once Dune::MultiTypeBlockMatrix +//! exposes std::tuple's constructor. +template +class MultiTypeBlockMatrix : std::tuple +{ + using ParentType = std::tuple; + public: + + using ParentType::ParentType; + + /** + * own class' type + */ + using type = MultiTypeBlockMatrix; + + /** \brief Type used for sizes */ + using size_type = std::size_t; + + using field_type = typename FirstRow::field_type; + + /** \brief Return the number of matrix rows */ + static constexpr size_type N() + { + return 1+sizeof...(Args); + } + + /** \brief Return the number of matrix columns */ + static constexpr size_type M() + { + return FirstRow::size(); + } + + template< size_type index > + auto + operator[] ( const std::integral_constant< size_type, index > indexVariable ) -> decltype(std::get(*this)) + { + DUNE_UNUSED_PARAMETER(indexVariable); + return std::get(*this); + } + + /** \brief Const random-access operator + * + * This is the const version of the random-access operator. See the non-const version for a full + * explanation of how to use it. + */ + template< size_type index > + auto + operator[] ( const std::integral_constant< size_type, index > indexVariable ) const -> decltype(std::get(*this)) + { + DUNE_UNUSED_PARAMETER(indexVariable); + return std::get(*this); + } + + /** \brief y = A x + */ + template + void mv (const X& x, Y& y) const { + static_assert(X::size() == M(), "length of x does not match row length"); + static_assert(Y::size() == N(), "length of y does not match row count"); + y = 0; //reset y (for mv uses umv) + umv(x,y); + } + + /** \brief y += A x + */ + template + void umv (const X& x, Y& y) const { + static_assert(X::size() == M(), "length of x does not match row length"); + static_assert(Y::size() == N(), "length of y does not match row count"); + using namespace Dune::Hybrid; + forEach(integralRange(Dune::Hybrid::size(y)), [&](auto&& i) { + using namespace Dune::Hybrid; // needed for icc, see issue #31 + forEach(integralRange(Dune::Hybrid::size(x)), [&](auto&& j) { + (*this)[i][j].umv(x[j], y[i]); + }); + }); + } + + /** \brief y += alpha A x + */ + template + void usmv (const AlphaType& alpha, const X& x, Y& y) const { + static_assert(X::size() == M(), "length of x does not match row length"); + static_assert(Y::size() == N(), "length of y does not match row count"); + using namespace Dune::Hybrid; + forEach(integralRange(Dune::Hybrid::size(y)), [&](auto&& i) { + using namespace Dune::Hybrid; // needed for icc, see issue #31 + forEach(integralRange(Dune::Hybrid::size(x)), [&](auto&& j) { + (*this)[i][j].usmv(alpha, x[j], y[i]); + }); + }); + } + +}; + +/*! + * \brief a function to get a MultiTypeBlockVector with const references to some entries of another MultiTypeBlockVector + * \param v a MultiTypeBlockVector + * \param indices the indices of the entries that should be referenced + * TODO can be removed if it gets implemented in dumux-master + */ +template +auto partial(const Dune::MultiTypeBlockVector& v, Dune::index_constant... indices) +{ + return Dune::MultiTypeBlockVector>>>...>(v[indices]...); +} + + + + } // end namespace Detail /*! @@ -966,6 +1078,8 @@ private: solveLinearSystem(deltaU); solveTimer.stop(); + std::cout << "Solve took " << solveTimer.elapsed() << std::endl; + /////////////// // update /////////////// @@ -1183,16 +1297,142 @@ private: { assert(this->checkSizesOfSubMatrices(A) && "Sub-blocks of MultiTypeBlockMatrix have wrong sizes!"); - // create the bcrs matrix the IterativeSolver backend can handle - const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); - - // get the new matrix sizes - const std::size_t numRows = M.N(); - assert(numRows == M.M()); + using namespace Dune::Indices; + + const auto& A00 = A[_0][_0]; + const auto& A11 = A[_1][_1]; + const auto& A22 = A[_2][_2]; + + const auto& A01 = A[_0][_1]; + const auto& A02 = A[_0][_2]; + + const auto& A10 = A[_1][_0]; + const auto& A12 = A[_1][_2]; + + const auto& A20 = A[_2][_0]; + const auto& A21 = A[_2][_1]; + + + using FirstRow = Dune::MultiTypeBlockVector&, const std::decay_t&, const std::decay_t&>; + using SecondRow = Dune::MultiTypeBlockVector&, const std::decay_t&, const std::decay_t&>; + using ThirdRow = Dune::MultiTypeBlockVector&, const std::decay_t&, const std::decay_t&>; + + using NewA = Detail::MultiTypeBlockMatrix; + + FirstRow firstRow(A11, A10, A12); + SecondRow secondRow(A01, A00, A02); + ThirdRow thirdRow(A21, A20, A22); + + const NewA newA = NewA(firstRow, secondRow, thirdRow); + + const auto newM = MatrixConverter::multiTypeToBCRSMatrix(newA); + + using namespace Dune::Indices; + + auto newB = partial(b, _1, _0, _2); + const auto newBtmp = VectorConverter::multiTypeToBlockVector(newB); + + auto newX = partial(x, _1, _0, _2); + + + // const auto& firstRow = partial(A[Dune::index_constant<0>{}], Dune::index_constant<1>{}, Dune::index_constant<0>{}, Dune::index_constant<2>{}); + // const auto& secondRow = partial(A[Dune::index_constant<1>{}], Dune::index_constant<1>{}, Dune::index_constant<0>{}, Dune::index_constant<2>{}); + // const auto& thirdRow = partial(A[Dune::index_constant<2>{}], Dune::index_constant<1>{}, Dune::index_constant<0>{}, Dune::index_constant<2>{}); + // + // + // using NewA = Detail::MultiTypeBlockMatrix; + // + // const auto newA = NewA(firstRow, secondRow, thirdRow); + + + +// // create the bcrs matrix the IterativeSolver backend can handle +// const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); +// +// // get the new matrix sizes + const std::size_t numRows = newM.N(); +// assert(numRows == M.M()); +// +// // create the vector the IterativeSolver backend can handle +// const auto bTmp = VectorConverter::multiTypeToBlockVector(b); +// assert(bTmp.size() == numRows); +// +// static const bool printmatrix = getParam("Problem.PrintMatrix", false); +// +// if (printmatrix) +// { +// static int counter = 0; +// +// std::ofstream logFile; +// const auto rank = Dune::MPIHelper::getCollectiveCommunication().rank(); +// // logFile.open("solver_log_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log"); +// // +// // const auto& A00 = A[Dune::index_constant<0>{}][Dune::index_constant<0>{}]; +// // const auto& A11 = A[Dune::index_constant<1>{}][Dune::index_constant<1>{}]; +// // const auto& A22 = A[Dune::index_constant<2>{}][Dune::index_constant<2>{}]; +// // +// // const auto& A01 = A[Dune::index_constant<0>{}][Dune::index_constant<1>{}]; +// // const auto& A02 = A[Dune::index_constant<0>{}][Dune::index_constant<2>{}]; +// // +// // +// // const auto& A10 = A[Dune::index_constant<1>{}][Dune::index_constant<0>{}]; +// // const auto& A12 = A[Dune::index_constant<1>{}][Dune::index_constant<2>{}]; +// // +// // +// // const auto& A20 = A[Dune::index_constant<2>{}][Dune::index_constant<0>{}]; +// // const auto& A21 = A[Dune::index_constant<2>{}][Dune::index_constant<1>{}]; +// +// +// +// +// // if (A02.frobenius_norm() > 1e-18) +// // DUNE_THROW(Dune::InvalidStateException, "A02 > 0"); +// // +// // if (A12.frobenius_norm() > 1e-18) +// // DUNE_THROW(Dune::InvalidStateException, "A12 > 0"); +// // +// // if (A20.frobenius_norm() > 1e-18) +// // DUNE_THROW(Dune::InvalidStateException, "A20 > 0"); +// // +// // if (A21.frobenius_norm() > 1e-18) +// // DUNE_THROW(Dune::InvalidStateException, "A21 > 0"); +// +// +// +// Dune::writeMatrixToMatlab(A00, "matrix00_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// Dune::writeMatrixToMatlab(A11, "matrix11_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// Dune::writeMatrixToMatlab(A22, "matrix22_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// +// +// Dune::writeMatrixToMatlab(A01, "matrix01_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// Dune::writeMatrixToMatlab(A02, "matrix02_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// +// Dune::writeMatrixToMatlab(A10, "matrix10_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// Dune::writeMatrixToMatlab(A12, "matrix12_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// +// Dune::writeMatrixToMatlab(A20, "matrix20_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// Dune::writeMatrixToMatlab(A21, "matrix21_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".mat"); +// +// // +// logFile.open("solver_residual_log0_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log"); +// Dune::printvector(logFile, b[Dune::index_constant<1>{}], "", ""); +// logFile.close(); +// +// logFile.open("solver_residual_log1_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log"); +// Dune::printvector(logFile, b[Dune::index_constant<0>{}], "", ""); +// logFile.close(); +// +// logFile.open("solver_residual_log2_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log"); +// Dune::printvector(logFile, b[Dune::index_constant<2>{}], "", ""); +// logFile.close(); +// +// // logFile.open("solver_residual_log_" + std::to_string(rank) + "_iter_" + std::to_string(counter) + ".log"); +// // Dune::printvector(logFile, b[Dune::index_constant<0>{}], "", ""); +// +// ++counter; +// +// } - // create the vector the IterativeSolver backend can handle - const auto bTmp = VectorConverter::multiTypeToBlockVector(b); - assert(bTmp.size() == numRows); // create a blockvector to which the linear solver writes the solution using VectorBlock = typename Dune::FieldVector; @@ -1200,11 +1440,11 @@ private: BlockVector y(numRows); // solve - const bool converged = ls.solve(M, y, bTmp); + const bool converged = ls.solve(newM, y, newBtmp); // copy back the result y into x if(converged) - VectorConverter::retrieveValues(x, y); + VectorConverter::retrieveValues(newX, y); return converged; } -- GitLab From 30756a360530073550b289cf4444fe7b0dae61db Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sun, 7 Mar 2021 15:58:27 +0100 Subject: [PATCH 02/15] [pnm2p] Clean up invasionstate --- .../porenetwork/2p/gridfluxvariablescache.hh | 11 ++- dumux/porenetwork/2p/invasionstate.hh | 75 ++++++------------- dumux/porenetwork/2p/model.hh | 4 +- 3 files changed, 29 insertions(+), 61 deletions(-) diff --git a/dumux/porenetwork/2p/gridfluxvariablescache.hh b/dumux/porenetwork/2p/gridfluxvariablescache.hh index af5d4afa88..a6a1969aa3 100644 --- a/dumux/porenetwork/2p/gridfluxvariablescache.hh +++ b/dumux/porenetwork/2p/gridfluxvariablescache.hh @@ -36,13 +36,12 @@ namespace Dumux::PoreNetwork { * \ingroup PNMTwoPModel * \brief Flux variable caches traits */ -template +template> struct PNMTwoPDefaultGridFVCTraits { using Problem = P; using FluxVariablesCache = FVC; - using Indices = I; - using Labels = L; + using InvasionState = IS; template using LocalView = PNMTwoPElementFluxVariablesCache; @@ -69,9 +68,9 @@ class PNMTwoPGridFluxVariablesCache { using Problem = typename Traits::Problem; using ThisType = PNMTwoPGridFluxVariablesCache; - using InvasionState = TwoPInvasionState; + using InvasionState = typename Traits::InvasionState; - public: +public: //! export the flux variable cache type using FluxVariablesCache = typename Traits::FluxVariablesCache; @@ -140,7 +139,7 @@ class PNMTwoPGridFluxVariablesCache { using Problem = typename Traits::Problem; using ThisType = PNMTwoPGridFluxVariablesCache; - using InvasionState = TwoPInvasionState; + using InvasionState = typename Traits::InvasionState; public: //! export the flux variable cache type diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh index e6fce845b8..ca9e9d065b 100644 --- a/dumux/porenetwork/2p/invasionstate.hh +++ b/dumux/porenetwork/2p/invasionstate.hh @@ -26,25 +26,12 @@ #include #include +#include #include -#include #include namespace Dumux::PoreNetwork { -#ifndef DOXYGEN -namespace Detail { -// helper struct detecting if the user-defined problem params class has a globalCapillaryPressure function -struct hasGlobalCapillaryPressure -{ - template - auto operator()(const Problem& p) - -> decltype(p.globalCapillaryPressure()) - {} -}; -} // end namespace Detail -#endif - /*! * \ingroup PNMTwoPModel * \brief The grid flux variables cache for the two-phase PNM hodling the invasion state of the throats @@ -55,6 +42,13 @@ class TwoPInvasionState { using Problem = P; + template + using GlobalCapillaryPressureDetector = decltype(std::declval().globalCapillaryPressure()); + + template + static constexpr bool hasGlobalCapillaryPressure() + { return Dune::Std::is_detected::value; } + public: TwoPInvasionState(const Problem& problem) : problem_(problem) @@ -71,10 +65,10 @@ public: } numThroatsInvaded_ = std::count(invadedCurrentIteration_.begin(), invadedCurrentIteration_.end(), true); - verbose_ = getParamFromGroup(problem.paramGroup(), "InvasionState.InvasionStateVerbosity", true); + verbose_ = getParamFromGroup(problem.paramGroup(), "InvasionState.Verbosity", true); restrictToGlobalCapillaryPressure_ = getParamFromGroup(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false); - if (decltype(isValid(Detail::hasGlobalCapillaryPressure())(problem)){}) + if (hasGlobalCapillaryPressure()) { if (restrictToGlobalCapillaryPressure_) std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl; @@ -118,8 +112,11 @@ public: if (invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf])) { hasChangedInCurrentIteration_ = true; - static constexpr auto cachingEnabled = std::integral_constant{}; - updateChangedFluxVarCache_(element, fvGeometry, gridFluxVarsCache, elemVolVars, scvf, cachingEnabled); + if constexpr (GridFluxVariablesCache::cachingEnabled) + { + const auto eIdx = problem_.gridGeometry().elementMapper().index(element); + gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invadedCurrentIteration_[eIdx]); + } } } } @@ -208,8 +205,8 @@ private: const auto wPhaseIdx = spatialParams.template wettingPhase(element, elemVolVars); // Block non-wetting phase flux out of the outlet - static const bool blockOutlet = getParamFromGroup(problem_.paramGroup(), "InvasionState.BlockOutletForNonwettingPhase", true); - if (blockOutlet && gridGeometry.throatLabel(eIdx) == Labels::outlet) + static const auto blockNonwettingPhase = getParamFromGroup>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector{Labels::outlet}); + if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end()) { invadedCurrentIteration_[eIdx] = false; return false; @@ -223,7 +220,7 @@ private: const Scalar pcSnapoff = fluxVarsCache.pcSnapoff(); // check if there is a user-specified global capillary pressure which needs to be obeyed - if (maybeRestrictToGlobalCapillaryPressure_(pcEntry, decltype(isValid(Detail::hasGlobalCapillaryPressure())(problem_)){})) + if (maybeRestrictToGlobalCapillaryPressure_(pcEntry)) { if (*pcMax > pcEntry) { @@ -267,43 +264,17 @@ private: return invadedBeforeSwitch != invadedAfterSwitch; } - //! Update the fluxVarsCache after an invasion/snap-off occured (only for caching enabled) - template - void updateChangedFluxVarCache_(const Element& element, - const FVElementGeometry& fvGeometry, - GridFluxVariablesCache& gridFluxVarsCache, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf, - std::true_type) - { - const auto eIdx = problem_.gridGeometry().elementMapper().index(element); - gridFluxVarsCache.cache(eIdx, scvf.index()).update(problem_, element, fvGeometry, elemVolVars, scvf, invaded(element)); - } - - //! Do nothing if the fluxVarCaches are not cached globally - template - void updateChangedFluxVarCache_(const Element& element, - const FVElementGeometry& fvGeometry, - GridFluxVariablesCache& gridFluxVarsCache, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf, - std::false_type) {} - //! If the user has specified a global capillary pressure, check if it is lower than the given entry capillary pressure. //! This may be needed to exactly reproduce pc-S curves given by static network models. template - bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry, std::true_type) const + bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pcEntry) const { - return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure()); + if constexpr (hasGlobalCapillaryPressure()) + return restrictToGlobalCapillaryPressure_ && (pcEntry > problem_.globalCapillaryPressure()); + else + return false; } - //! Do nothing here. - template - bool maybeRestrictToGlobalCapillaryPressure_(const Scalar pc, std::false_type) const - { return false; } - std::vector invadedCurrentIteration_; std::vector invadedPreviousTimeStep_; bool hasChangedInCurrentIteration_ = false; diff --git a/dumux/porenetwork/2p/model.hh b/dumux/porenetwork/2p/model.hh index 8d676cf2f0..16591c003b 100644 --- a/dumux/porenetwork/2p/model.hh +++ b/dumux/porenetwork/2p/model.hh @@ -130,10 +130,8 @@ struct GridFluxVariablesCache private: static constexpr bool enableCache = getPropValue(); using Problem = GetPropType; - using Indices = typename GetPropType::Indices; - using Labels = GetPropType; using FluxVariablesCache = GetPropType; - using Traits = Dumux::PoreNetwork::PNMTwoPDefaultGridFVCTraits; + using Traits = Dumux::PoreNetwork::PNMTwoPDefaultGridFVCTraits; public: using type = Dumux::PoreNetwork::PNMTwoPGridFluxVariablesCache; }; -- GitLab From 141be322a57b45d8b71cd8a7643ce00bfd02e78c Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Mon, 8 Mar 2021 08:32:43 +0100 Subject: [PATCH 03/15] [pnm2p] Further cleanup of invasion state --- dumux/porenetwork/2p/invasionstate.hh | 78 +++++++++++++++++---------- 1 file changed, 50 insertions(+), 28 deletions(-) diff --git a/dumux/porenetwork/2p/invasionstate.hh b/dumux/porenetwork/2p/invasionstate.hh index ca9e9d065b..c27a8d307c 100644 --- a/dumux/porenetwork/2p/invasionstate.hh +++ b/dumux/porenetwork/2p/invasionstate.hh @@ -49,6 +49,8 @@ class TwoPInvasionState static constexpr bool hasGlobalCapillaryPressure() { return Dune::Std::is_detected::value; } + enum class EventType {invasion, snapOff, none}; + public: TwoPInvasionState(const Problem& problem) : problem_(problem) @@ -68,7 +70,7 @@ public: verbose_ = getParamFromGroup(problem.paramGroup(), "InvasionState.Verbosity", true); restrictToGlobalCapillaryPressure_ = getParamFromGroup(problem.paramGroup(), "InvasionState.RestrictInvasionToGlobalCapillaryPressure", false); - if (hasGlobalCapillaryPressure()) + if constexpr (hasGlobalCapillaryPressure()) { if (restrictToGlobalCapillaryPressure_) std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl; @@ -109,8 +111,9 @@ public: for (auto&& scvf : scvfs(fvGeometry)) { // checks if invasion or snap-off occured after Newton iteration step - if (invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf])) + if (const auto invasionResult = invasionSwitch_(element, elemVolVars, elemFluxVarsCache[scvf]); invasionResult) { + const auto localScvIdxWithCriticalPc = invasionResult.localScvIdxWithCriticalPc; hasChangedInCurrentIteration_ = true; if constexpr (GridFluxVariablesCache::cachingEnabled) { @@ -191,29 +194,39 @@ private: //! The switch for determining the invasion state of a pore throat. Called at the end of each Newton step. template - bool invasionSwitch_(const Element& element, + auto invasionSwitch_(const Element& element, const ElementVolumeVariables& elemVolVars, const FluxVariablesCache& fluxVarsCache) { + using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type; const auto& gridGeometry = problem_.gridGeometry(); const auto& spatialParams = problem_.spatialParams(); const auto eIdx = gridGeometry.elementMapper().index(element); bool invadedBeforeSwitch = invadedCurrentIteration_[eIdx]; bool invadedAfterSwitch = invadedBeforeSwitch; - const auto wPhaseIdx = spatialParams.template wettingPhase(element, elemVolVars); + // Result type, containing the local scv index of the pore from which the invasion/snap-off occurred + // Evaluates to 'false' if no invasion/snap-off occurred + struct Result + { + std::uint8_t localScvIdxWithCriticalPc; + Scalar criticalPc; + EventType event = EventType::none; + + operator bool() const + { return event != EventType::none; } + }; // Block non-wetting phase flux out of the outlet static const auto blockNonwettingPhase = getParamFromGroup>(problem_.paramGroup(), "InvasionState.BlockNonwettingPhaseAtThroatLabel", std::vector{Labels::outlet}); if (!blockNonwettingPhase.empty() && std::find(blockNonwettingPhase.begin(), blockNonwettingPhase.end(), gridGeometry.throatLabel(eIdx)) != blockNonwettingPhase.end()) { invadedCurrentIteration_[eIdx] = false; - return false; + return Result{}; // nothing happened } //Determine whether throat gets invaded or snap-off occurs - using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type; const std::array pc = { elemVolVars[0].capillaryPressure(), elemVolVars[1].capillaryPressure() }; const auto pcMax = std::max_element(pc.begin(), pc.end()); const Scalar pcEntry = fluxVarsCache.pcEntry(); @@ -229,7 +242,7 @@ private: } invadedCurrentIteration_[eIdx] = false; - return false; + return Result{}; //nothing happened } if (*pcMax > pcEntry) @@ -239,29 +252,38 @@ private: invadedCurrentIteration_[eIdx] = invadedAfterSwitch; - if (invadedAfterSwitch != invadedBeforeSwitch && verbose_) + if (invadedBeforeSwitch == invadedAfterSwitch) + return Result{}; // nothing happened + else { - const std::array sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) }; - const auto scvIdx = pcMax - pc.begin(); - const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, scvIdx, 1); - if (!invadedBeforeSwitch && invadedAfterSwitch) - { - std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :"; - std::cout << " pc: " << *pcMax; - std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars); - std::cout << ", sw: " << sw[scvIdx] << std::endl; - } - else if (invadedBeforeSwitch && !invadedAfterSwitch) - { - std::cout << "Snap-off occured at: " << eIdx << " from pore " << vIdx << " :"; - std::cout << " pc: " << *pcMax; - std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars); - std::cout << ", sw: " << sw[scvIdx] << std::endl; - } - else - DUNE_THROW(Dune::InvalidStateException, "Invalid Process "); + Result result; + result.localScvIdxWithCriticalPc = std::distance(pc.begin(), pcMax); + result.criticalPc = *pcMax; + result.event = !invadedBeforeSwitch && invadedAfterSwitch ? EventType::invasion : EventType::snapOff; + + if (verbose_) + { + const auto wPhaseIdx = spatialParams.template wettingPhase(element, elemVolVars); + const std::array sw = { elemVolVars[0].saturation(wPhaseIdx), elemVolVars[1].saturation(wPhaseIdx) }; + const auto vIdx = gridGeometry.gridView().indexSet().subIndex(element, result.localScvIdxWithCriticalPc, 1); + if (result.event == EventType::invasion) + { + std::cout << "Throat " << eIdx << " was invaded from pore " << vIdx << " :"; + std::cout << " pc: " << *pcMax; + std::cout << ", pcEntry: " << spatialParams.pcEntry(element, elemVolVars); + std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl; + } + else + { + std::cout << "Snap-off occured at throat " << eIdx << " from pore " << vIdx << " :"; + std::cout << " pc: " << *pcMax; + std::cout << ", pcSnapoff: " << spatialParams.pcSnapoff(element, elemVolVars); + std::cout << ", sw: " << sw[result.localScvIdxWithCriticalPc] << std::endl; + } + } + + return result; } - return invadedBeforeSwitch != invadedAfterSwitch; } //! If the user has specified a global capillary pressure, check if it is lower than the given entry capillary pressure. -- GitLab From b51ecac9a35160c1bab0bd62ba2d71fd8032f68b Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Tue, 23 Mar 2021 17:57:53 +0100 Subject: [PATCH 04/15] [localrulesforplatonicbody] Add option to set fixed high sw slope --- .../porenetwork/pore/2p/localrulesforplatonicbody.hh | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh index e13600e8b6..a2f24899e0 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh @@ -516,6 +516,9 @@ private: baseLawParamsPtr_ = &m->basicParams(); + static const auto highSwFixedSlopeInput = getParamFromGroup(paramGroup, "RegularizationHighSwFixedSlope", std::numeric_limits::quiet_NaN()); + highSwFixedSlope_ = highSwFixedSlopeInput; + // Note: we do not pre-calculate all end-point values (e.g., pcLowSwPcValue_ and pcDerivativeLowSw_) // as done in, e.g., VanGenuchten. This is because this law will generally be instantiated only as // a temporary object since all pore bodies generally have different parameters. @@ -550,7 +553,11 @@ private: // dpc_dsw at Sw = 1.0 Scalar pcDerivativeHighSwEnd_() const { - return (0.0 - pcHighSwPcValue_()) / (1.0 - pcHighSw_); + using std::isnan; + if (!isnan(highSwFixedSlope_)) + return highSwFixedSlope_; + else + return (0.0 - pcHighSwPcValue_()) / (1.0 - pcHighSw_); } const Spline& pcSpline_() const @@ -574,6 +581,7 @@ private: const BaseLawParams* baseLawParamsPtr_; mutable std::optional> optionalPcSpline_; bool highSwSplineZeroSlope_; + Scalar highSwFixedSlope_; }; /*! -- GitLab From 74c6fc240a6f09931e9df5b773f2f7ceb01c306d Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sun, 28 Mar 2021 16:21:02 +0200 Subject: [PATCH 05/15] [localrules] Fix linear regularisation --- .../porenetwork/pore/2p/localrulesforplatonicbody.hh | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh index a2f24899e0..c5d9455c24 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh @@ -321,11 +321,6 @@ private: if (sw < pcLowSw_) return pcLowSwPcValue_() + pcDerivativeLowSw_() * (sw - pcLowSw_); - auto linearCurveForHighSw = [&]() - { - return pcDerivativeHighSwEnd_()*(sw - 1.0); - }; - if (sw <= pcHighSw_) return {}; // standard else if (sw < 1.0) // regularized part below sw = 1.0 @@ -335,7 +330,10 @@ private: return pcHighSwPcValue_() * pow(((1.0-sw)/(1.0-pcHighSw_)), 1.0/3.0); else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::linear) - return linearCurveForHighSw(); + { + const Scalar slope = -pcHighSwPcValue_() / (1.0 - pcHighSw_); + return pcHighSwPcValue_() + (sw - pcHighSw_) * slope; + } else if (highSwRegularizationMethod_ == HighSwRegularizationMethod::spline) return pcSpline_().eval(sw); @@ -344,7 +342,7 @@ private: DUNE_THROW(Dune::NotImplemented, "Regularization not method not implemented"); } else // regularized part above sw = 1.0 - return linearCurveForHighSw(); + return pcDerivativeHighSwEnd_()*(sw - 1.0); } /*! -- GitLab From 3fd1673bc197e725c6479446049a44a97626b0f4 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sun, 4 Apr 2021 17:19:19 +0200 Subject: [PATCH 06/15] [PNM] Add diffusive enthalpy transport --- dumux/flux/porenetwork/fourierslaw.hh | 42 +++++++++++++++++++++++++++ dumux/porenetwork/properties.hh | 2 +- 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/dumux/flux/porenetwork/fourierslaw.hh b/dumux/flux/porenetwork/fourierslaw.hh index 6db79c10cc..87d66d7fc1 100644 --- a/dumux/flux/porenetwork/fourierslaw.hh +++ b/dumux/flux/porenetwork/fourierslaw.hh @@ -26,13 +26,21 @@ #define DUMUX_FLUX_PNM_FOURIERS_LAW_HH #include +#include namespace Dumux::PoreNetwork { +namespace Detail { + +struct NoDiffusionType {}; + +} // end namespace Detail + /*! * \ingroup PoreNetworkModels * \brief Specialization of Fourier's Law for the pore-network model. */ +template struct PNMFouriersLaw { @@ -68,7 +76,41 @@ struct PNMFouriersLaw const Scalar gradT = deltaT/fluxVarsCache.throatLength(); heatflux += thermalConductivity*gradT*area; + + if constexpr (!std::is_same_v) + heatflux += componentEnthalpyFlux_(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache, phaseIdx); + } + + return heatflux; + } + +private: + template + static auto componentEnthalpyFlux_(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const typename FVElementGeometry::SubControlVolumeFace& scvf, + const ElementFluxVariablesCache& elemFluxVarsCache, + const int phaseIdx) + { + using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type; + Scalar heatflux = 0.0; + using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem; + const auto diffusiveFlux = MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache); + for (int compIdx = 0; compIdx < ElementVolumeVariables::VolumeVariables::numFluidComponents(); ++compIdx) + { + const bool insideIsUpstream = diffusiveFlux[compIdx] > 0.0; + const auto& upstreamVolVars = insideIsUpstream ? elemVolVars[scvf.insideScvIdx()] : elemVolVars[scvf.outsideScvIdx()]; + const Scalar componentEnthalpy = FluidSystem::componentEnthalpy(upstreamVolVars.fluidState(), phaseIdx, compIdx); + + if (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) + heatflux += diffusiveFlux[compIdx] * componentEnthalpy; + else + heatflux += diffusiveFlux[compIdx] * FluidSystem::molarMass(compIdx) * componentEnthalpy; } + return heatflux; } }; diff --git a/dumux/porenetwork/properties.hh b/dumux/porenetwork/properties.hh index 628182f78b..dec06d8a0a 100644 --- a/dumux/porenetwork/properties.hh +++ b/dumux/porenetwork/properties.hh @@ -68,7 +68,7 @@ public: }; template -struct HeatConductionType { using type = Dumux::PoreNetwork::PNMFouriersLaw; }; +struct HeatConductionType { using type = Dumux::PoreNetwork::PNMFouriersLaw<>; }; //! The labels template -- GitLab From e3820aad3e2e1bfb5547f218dbe5baaf1c543b39 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Sun, 4 Apr 2021 17:22:01 +0200 Subject: [PATCH 07/15] [PNM] Add surface tension to spatialParams --- .../pore/2p/localrulesforplatonicbody.hh | 10 +++------- .../spatialparams/porenetwork/porenetwork2p.hh | 18 ++++++++++++++++-- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh index c5d9455c24..c97b28fcc1 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh @@ -57,10 +57,8 @@ struct PlatonicBodyParams const ElemSol& elemSol) : shape_(spatialParams.gridGeometry().poreGeometry()[scv.dofIndex()]) , radius_(spatialParams.poreInscribedRadius(element, scv, elemSol)) - { - static const Scalar surfaceTension = getParam("SpatialParams.SurfaceTension", 0.0725); // TODO - surfaceTension_ = surfaceTension; - } + , surfaceTension_(spatialParams.surfaceTension(element, scv, elemSol)) + {} template void update(const SpatialParams& spatialParams, @@ -71,9 +69,7 @@ struct PlatonicBodyParams const auto& gridGeometry = spatialParams.gridGeometry(); shape_ = gridGeometry.poreGeometry()[scv.dofIndex()]; radius_ = spatialParams.poreInscribedRadius(element, scv, elemSol); - - static const Scalar surfaceTension = getParam("SpatialParams.SurfaceTension", 0.0725); // TODO - surfaceTension_ = surfaceTension; + surfaceTension_ = spatialParams.surfaceTension(element, scv, elemSol); } Pore::Shape poreShape() const { return shape_; } diff --git a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh index 6f51c8ca85..feab65a281 100644 --- a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh +++ b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh @@ -154,9 +154,23 @@ public: } /*! - * \brief Returns the parameter object for the Brooks-Corey material law. + * \brief Returns the surface tension \f$ N/m\f$ * - * In this test, we use element-wise distributed material parameters. + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + */ + template + Scalar surfaceTension(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + static const Scalar gamma = getParam("SpatialParams.SurfaceTension", 0.0725); // default to surface tension of water/air + return gamma; + } + + /*! + * \brief Returns the parameter object for the pore-local pc-Sw law * * \param element The current element * \param scv The sub-control volume inside the element. -- GitLab From 36a312e9c831f6922b8750d1eb60343dc05e5ea4 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 7 Apr 2021 08:11:47 +0200 Subject: [PATCH 08/15] [staggeredgeohelper] Using std::abs --- .../staggered/freeflow/staggeredgeometryhelper.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh index 61948a2ac8..dd7ba55e5f 100644 --- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh +++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh @@ -133,7 +133,8 @@ template inline static unsigned int directionIndex(Vector&& vector) { const auto eps = 1e-8; - return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin(); + using std::abs; + return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return x > eps; } ) - vector.begin(); } /*! -- GitLab From 9ab60fd3f2f6701656e1c29a99b98e5acacf6bab Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 7 Apr 2021 08:13:26 +0200 Subject: [PATCH 09/15] [multishapelocalules] Pass ScalarT --- .../porenetwork/pore/2p/multishapelocalrules.hh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh index e248654efb..68ea4c5a67 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh @@ -34,11 +34,11 @@ namespace Dumux::PoreNetwork::FluidMatrix { template struct LocalRulesTraits { - using Tetrahedron = TwoPLocalRulesPlatonicBodyDefault; - using Cube = TwoPLocalRulesPlatonicBodyDefault; - using Octahedron = TwoPLocalRulesPlatonicBodyDefault; - using Icosahedron = TwoPLocalRulesPlatonicBodyDefault; - using Dodecahedron = TwoPLocalRulesPlatonicBodyDefault; + using Tetrahedron = TwoPLocalRulesPlatonicBodyDefault; + using Cube = TwoPLocalRulesPlatonicBodyDefault; + using Octahedron = TwoPLocalRulesPlatonicBodyDefault; + using Icosahedron = TwoPLocalRulesPlatonicBodyDefault; + using Dodecahedron = TwoPLocalRulesPlatonicBodyDefault; }; template -- GitLab From 1eaf20d4358f7f4d2adacf8d9bf3525b0ca243d2 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 7 Apr 2021 09:35:42 +0200 Subject: [PATCH 10/15] fixup geohelper --- .../staggered/freeflow/staggeredgeometryhelper.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh index dd7ba55e5f..65514ea3c3 100644 --- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh +++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh @@ -134,7 +134,7 @@ inline static unsigned int directionIndex(Vector&& vector) { const auto eps = 1e-8; using std::abs; - return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return x > eps; } ) - vector.begin(); + return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return abs(x) > eps; } ) - vector.begin(); } /*! -- GitLab From 3b2b42fac9a49ea56c8eaae61600b98a491051a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 7 Apr 2021 14:36:40 +0200 Subject: [PATCH 11/15] [fvassembler] add missing query if parallel --- dumux/assembly/fvassembler.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh index 3195589c59..efc3cf647e 100644 --- a/dumux/assembly/fvassembler.hh +++ b/dumux/assembly/fvassembler.hh @@ -170,7 +170,7 @@ public: // issue a warning if the caluclation is used in parallel with overlap static bool warningIssued = false; - if (gridView().overlapSize(0) == 0) + if (gridView().comm().size() > 1 && gridView().overlapSize(0) == 0) { if constexpr (isBox) { -- GitLab From 21aaa86d089d780a260db6eb5e0b7fe5c2145a7b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 7 Apr 2021 14:36:56 +0200 Subject: [PATCH 12/15] [md][fvassembler] add missing query if parallel --- dumux/multidomain/fvassembler.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh index 0b94c8eb01..656a8ce658 100644 --- a/dumux/multidomain/fvassembler.hh +++ b/dumux/multidomain/fvassembler.hh @@ -239,7 +239,7 @@ public: const auto gridGeometry = std::get(gridGeometryTuple_); const auto& gridView = gridGeometry->gridView(); - if (gridView.overlapSize(0) == 0) + if (gridView.comm().size() && gridView.overlapSize(0) == 0) { if constexpr (GridGeometry::discMethod == DiscretizationMethod::box) { -- GitLab From 8ffca9ee0fd42e36fa927e00ff17ddeda9e2cf6d Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 7 Apr 2021 14:14:44 +0200 Subject: [PATCH 13/15] [test] Add test with nonoverlapping grid manager that cannot communicate --- .../1p/incompressible/CMakeLists.txt | 8 +++ .../1p/incompressible/main.cc | 3 + .../1p/incompressible/properties.hh | 61 +++++++++++++++++++ 3 files changed, 72 insertions(+) diff --git a/test/porousmediumflow/1p/incompressible/CMakeLists.txt b/test/porousmediumflow/1p/incompressible/CMakeLists.txt index 9eac64ac66..02adbcc245 100644 --- a/test/porousmediumflow/1p/incompressible/CMakeLists.txt +++ b/test/porousmediumflow/1p/incompressible/CMakeLists.txt @@ -145,3 +145,11 @@ dumux_add_test(NAME test_1p_incompressible_mpfa_extrude_distorted -Problem.CheckIsConstantVelocity true -Problem.EnableGravity false -Grid.File ./grids/randomlydistorted.dgf) + +# check grids without communicate method (using box and numeric differentiation) +dumux_add_test(NAME test_1p_incompressible_box_numdiff_no_communicate + LABELS porousmediumflow 1p + SOURCES main.cc + COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox NUMDIFFMETHOD=DiffMethod::numeric GRIDTYPE=Dumux::NoCommunicateGrid<2> + COMMAND ./test_1p_incompressible_box_numdiff_no_communicate + CMD_ARGS -Grid.Overlap 0) diff --git a/test/porousmediumflow/1p/incompressible/main.cc b/test/porousmediumflow/1p/incompressible/main.cc index 3f8e191a6c..f0d7a5e666 100644 --- a/test/porousmediumflow/1p/incompressible/main.cc +++ b/test/porousmediumflow/1p/incompressible/main.cc @@ -149,6 +149,9 @@ int main(int argc, char** argv) // output result to vtk vtkWriter.write(1.0); + // output residual norm (test assembler interface) + std::cout << "Residual norm: " << assembler->residualNorm(x) << std::endl; + timer.stop(); const bool checkIsConstantVelocity = getParam("Problem.CheckIsConstantVelocity", false); diff --git a/test/porousmediumflow/1p/incompressible/properties.hh b/test/porousmediumflow/1p/incompressible/properties.hh index d333cff6fa..7df2043848 100644 --- a/test/porousmediumflow/1p/incompressible/properties.hh +++ b/test/porousmediumflow/1p/incompressible/properties.hh @@ -50,6 +50,67 @@ #define GRIDTYPE Dune::YaspGrid<2> #endif +//////////////////////////////////////////////////// +// A fake grid that cannot communicate. // +// Can be used to make sure that compilation and // +// sequential execution also work for such grids. // +//////////////////////////////////////////////////// +namespace Dumux { + +template +class NoCommunicateGrid; + +template +class NoCommunicateGridLeafGridView +: public Dune::YaspGrid::LeafGridView +{ + using ParentType = typename Dune::YaspGrid::LeafGridView; +public: + using ParentType::ParentType; + + struct Traits : public ParentType::Traits + { using Grid = NoCommunicateGrid; }; +}; + +template +class NoCommunicateGrid : public Dune::YaspGrid +{ + using ParentType = Dune::YaspGrid; +public: + using ParentType::ParentType; + struct Traits : public ParentType::Traits + { using LeafGridView = NoCommunicateGridLeafGridView; }; + + using LeafGridView = NoCommunicateGridLeafGridView; + + typename Traits::LeafGridView leafGridView() const + { return NoCommunicateGridLeafGridView(*this); } +private: + using ParentType::communicate; + using ParentType::communicateCodim; +}; + +template +class GridManager> : public GridManager> +{ + using ParentType = GridManager>; +public: + using ParentType::ParentType; + using Grid = NoCommunicateGrid; + Grid& grid() { return static_cast&>(ParentType::grid()); } +}; + +} // end namespace Dumux + +namespace Dune::Capabilities { + +template +struct canCommunicate, codim> +{ static constexpr bool v = false; }; + +} // end namespace Dune::Capabilities +///////////////////////////////////////// + namespace Dumux::Properties { // Create new type tags namespace TTag { -- GitLab From 776dfaddb163ef889f622475ad495ebfa4409bfb Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 14 Apr 2021 15:16:00 +0200 Subject: [PATCH 14/15] fixup fouriers law --- dumux/flux/porenetwork/fourierslaw.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/dumux/flux/porenetwork/fourierslaw.hh b/dumux/flux/porenetwork/fourierslaw.hh index 87d66d7fc1..2f091f6a20 100644 --- a/dumux/flux/porenetwork/fourierslaw.hh +++ b/dumux/flux/porenetwork/fourierslaw.hh @@ -26,6 +26,7 @@ #define DUMUX_FLUX_PNM_FOURIERS_LAW_HH #include +#include #include namespace Dumux::PoreNetwork { -- GitLab From dddf016eb30dab76fd089cba3ea6513e159943e8 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 14 Apr 2021 15:17:47 +0200 Subject: [PATCH 15/15] [pnm][multishapelocalrules] Deprecate makeParams() * function doesn't add any value * use ctor of BasicParams directly --- .../porenetwork/pore/2p/multishapelocalrules.hh | 1 + dumux/material/spatialparams/porenetwork/porenetwork2p.hh | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh index 68ea4c5a67..5f3bab10a7 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh @@ -102,6 +102,7 @@ public: { return 2; } template + [[deprecated("Use constructor of BasicParams directly. Will be removed before release 3.4")]] static BasicParams makeParams(const SpatialParams& spatialParams, const Element& element, const SubControlVolume& scv, diff --git a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh index feab65a281..032a9bb18d 100644 --- a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh +++ b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh @@ -182,7 +182,7 @@ public: const SubControlVolume& scv, const ElementSolution& elemSol) const { - const auto params = LocalRules::makeParams(*this, element, scv, elemSol); + const auto params = typename LocalRules::BasicParams(*this, element, scv, elemSol); return makeFluidMatrixInteraction(LocalRules(params, "SpatialParams")); } -- GitLab