diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh index 3195589c599e6f61cd17e8b821dade9b4b3a914b..efc3cf647e76efe2eb86f42e55839efa6a6fe544 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) { diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh index 61948a2ac8b2fd3c558110734d18e7825906f856..65514ea3c3b785c13c68f051e6c48de53ca4a2c0 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 abs(x) > eps; } ) - vector.begin(); } /*! diff --git a/dumux/flux/porenetwork/fourierslaw.hh b/dumux/flux/porenetwork/fourierslaw.hh index 6db79c10cc4b2405b5f1880f1364ed5ae54ed418..2f091f6a20db667fce4a24f20ba5c3b608aed6cd 100644 --- a/dumux/flux/porenetwork/fourierslaw.hh +++ b/dumux/flux/porenetwork/fourierslaw.hh @@ -26,13 +26,22 @@ #define DUMUX_FLUX_PNM_FOURIERS_LAW_HH #include +#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 +77,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/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/localrulesforplatonicbody.hh index e13600e8b6969f5ef85c130d609107e4fb78c0c4..c97b28fcc132065e4704ec0b226b47a7aa272033 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_; } @@ -321,11 +317,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 +326,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 +338,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); } /*! @@ -516,6 +510,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 +547,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 +575,7 @@ private: const BaseLawParams* baseLawParamsPtr_; mutable std::optional> optionalPcSpline_; bool highSwSplineZeroSlope_; + Scalar highSwFixedSlope_; }; /*! diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh b/dumux/material/fluidmatrixinteractions/porenetwork/pore/2p/multishapelocalrules.hh index e248654efbd0c1dc622956f53e5c3b0118bf0240..5f3bab10a793a4fb85b8f5de947f634259a9d2f6 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 @@ -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 6f51c8ca857cf9dc9ea3c64baaeb1f73082249e9..032a9bb18da15dc63cffa720afd1744e2f168d5a 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. @@ -168,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")); } diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh index 0b94c8eb018ba6dce38b67f9016a6185123ae6f2..656a8ce65819576b7d775e0712d5437244884d5d 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) { diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 61ccca1595aa0d84bbfc9830d0587730ab40d888..5b11b9cc70b9ad6a2f98084eac6f1a250b3b7397 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; } diff --git a/dumux/porenetwork/2p/gridfluxvariablescache.hh b/dumux/porenetwork/2p/gridfluxvariablescache.hh index af5d4afa8805c910009d0b34556fd7302026a36b..a6a1969aa39ec9a7224baa0851f62ac3807bec44 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 e6fce845b847e3f6482120d1682d872156fba694..c27a8d307c247fc9b6f659487e1963d41c2c0f64 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,15 @@ class TwoPInvasionState { using Problem = P; + template + using GlobalCapillaryPressureDetector = decltype(std::declval().globalCapillaryPressure()); + + template + static constexpr bool hasGlobalCapillaryPressure() + { return Dune::Std::is_detected::value; } + + enum class EventType {invasion, snapOff, none}; + public: TwoPInvasionState(const Problem& problem) : problem_(problem) @@ -71,10 +67,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 constexpr (hasGlobalCapillaryPressure()) { if (restrictToGlobalCapillaryPressure_) std::cout << "\n *** Invasion behavior is restricted by a global capillary pressure defined in the problem! *** \n" << std::endl; @@ -115,11 +111,15 @@ 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; - 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]); + } } } } @@ -194,36 +194,46 @@ 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 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; + 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(); 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) { @@ -232,7 +242,7 @@ private: } invadedCurrentIteration_[eIdx] = false; - return false; + return Result{}; //nothing happened } if (*pcMax > pcEntry) @@ -242,68 +252,51 @@ 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 "); - } - return invadedBeforeSwitch != invadedAfterSwitch; - } + Result result; + result.localScvIdxWithCriticalPc = std::distance(pc.begin(), pcMax); + result.criticalPc = *pcMax; + result.event = !invadedBeforeSwitch && invadedAfterSwitch ? EventType::invasion : EventType::snapOff; - //! 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)); - } + 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; + } + } - //! 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) {} + return result; + } + } //! 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 8d676cf2f0e81b14cff8683c430b6961fd2bfcda..16591c003b4b0813a75e0586e7f5a3993b7dffeb 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; }; diff --git a/dumux/porenetwork/properties.hh b/dumux/porenetwork/properties.hh index 628182f78bb93359ae9ecf54b2919fde000102a3..dec06d8a0a3f017e099cf77919a7cbaa78d946cb 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 diff --git a/test/porousmediumflow/1p/incompressible/CMakeLists.txt b/test/porousmediumflow/1p/incompressible/CMakeLists.txt index 9eac64ac66d7312a9f2a8bca7ed142051abee048..02adbcc24585db8f2d13c44769a9fdaea54a2f93 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 3f8e191a6cd493c0997fc8d9b5c5ff82ea568d97..f0d7a5e666cc11a9814ac82f9a44567c8971127c 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 d333cff6fa25cd4a5a02c058b41585b11716afa6..7df2043848a99e25dc9b7a29e0414efef3b00ffd 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 {