diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index a9a1f2d4e08c1f8cc1bc9f897fe9b455136fdd75..c4689f2c284ca47e4a01ab86a06ee11dd47e3611 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -54,7 +54,7 @@ public: boundaryInfo_[eqIdx].visited = false; boundaryInfo_[eqIdx].isSymmetry = false; - boundaryInfo_[eqIdx].isBeaversJoseph = false; + boundaryInfo_[eqIdx].isSlipCondition = false; } /*! @@ -94,34 +94,67 @@ public: static_assert(AlwaysFalse::value, "Setting all boundary types to Neumann not permitted!"); } + // TODO: deprecate! /*! * \brief Set a boundary condition for a single equation to * Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.). */ + [[deprecated("This method will be removed after release (3.n). Use setSlipCondition instead!")]] void setBeaversJoseph(unsigned eqIdx) + { + setSlipCondition(eqIdx); + } + /*! + * \brief Set a boundary condition for a single equation to + * slip condition: bjs,bj or the new one by Elissa Eggenweiler + */ + void setSlipCondition(unsigned eqIdx) { resetEq(eqIdx); boundaryInfo_[eqIdx].visited = true; - boundaryInfo_[eqIdx].isBeaversJoseph = true; + boundaryInfo_[eqIdx].isSlipCondition = true; } + + // TODO: Deprecate/Remove? Where used? /*! * \brief Returns true if an equation is used to specify a * Beavers-Joseph(-Saffman) boundary condition. * * \param eqIdx The index of the equation */ + [[deprecated("This method will be removed after release (3.n). Use isSlipCondition instead!")]] bool isBeaversJoseph(unsigned eqIdx) const - { return boundaryInfo_[eqIdx].isBeaversJoseph; } + { return isSlipCondition(eqIdx); } + + /*! + * \brief Returns true if an equation is used to specify a + * slip condition + * + * \param eqIdx The index of the equation + */ + bool isSlipCondition(unsigned eqIdx) const + { return boundaryInfo_[eqIdx].isSlipCondition; } + // TODO: Deprecate/Remove? Where used? /*! * \brief Returns true if some equation is used to specify a * Beavers-Joseph(-Saffman) boundary condition. */ + [[deprecated("This method will be removed after release (3.n). Use hasSlipCondition instead!")]] bool hasBeaversJoseph() const + { + return hasSlipCondition(); + } + + /*! + * \brief Returns true if some equation is used to specify a + * slip condition + */ + bool hasSlipCondition() const { for (int i = 0; i < numEq; ++i) - if (boundaryInfo_[i].isBeaversJoseph) + if (boundaryInfo_[i].isSlipCondition) return true; return false; } @@ -131,7 +164,7 @@ protected: { bool visited; bool isSymmetry; - bool isBeaversJoseph; + bool isSlipCondition; }; std::array boundaryInfo_; diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index ab62ad440c08e765fadf82b31a0baeb9cadaef35..8dd33c0052412cd3d392b0e0824a8ec0848bdfc0 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -190,7 +190,7 @@ public: */ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the permeability must be returned in the acutal problem"); + DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the permeability must be returned in the actual problem"); } /*! @@ -200,7 +200,7 @@ public: */ Scalar alphaBJ(const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the alpha value must be returned in the acutal problem"); + DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the alpha value must be returned in the actual problem"); } /*! @@ -212,6 +212,23 @@ public: return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf)); } + /*! + * \brief Returns the scale separation value eps = (macroscopicLength/LengthOfPeriodicityCell). + */ + Scalar epsInterface(const SubControlVolumeFace& scvf) const + { + DUNE_THROW(Dune::NotImplemented, "When using the interface conditions for arbitrary flows to the interface, the epsInterface value must be returned in the actual problem"); + } + + /*! + * \brief Returns the boundary layer constant N^bl computed based on the exact pore geometry. + */ + Scalar factorNTangential(const SubControlVolumeFace& scvf) const + { + DUNE_THROW(Dune::NotImplemented, "When using the interface conditions for arbitrary flows to the interface, the factorNTangential value must be returned in the actual problem"); + } + + /*! * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann). */ @@ -220,26 +237,33 @@ public: return VelocityVector(0.0); } - //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used - const Scalar beaversJosephVelocity(const Element& element, + //! helper function to evaluate the slip velocity on the boundary when the slipCondition is used + const Scalar slipVelocity(const Element& element, const SubControlVolume& scv, const SubControlVolumeFace& ownScvf, const SubControlVolumeFace& faceOnPorousBoundary, - const Scalar velocitySelf, - const Scalar tangentialVelocityGradient) const + const Scalar velocitySelf, //vel perpendicular to tangential vel + const Scalar tangentialVelocityGradient) const //dv/dx (=0) { - // du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM) - // beta = alpha/sqrt(K) - const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary); + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + // du/dy + dv/dx = factor * (u_boundary-uPM) + Scalar factor; + if (newIc_) + { + factor = -1.0 / asImp_().epsInterface(faceOnPorousBoundary) / asImp_().factorNTangential(faceOnPorousBoundary); + } + else + { + factor = asImp_().betaBJ(element, faceOnPorousBoundary); //beta = alpha/sqrt(K) + } const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); // create a unit normal vector oriented in positive coordinate direction GlobalPosition orientation = ownScvf.unitOuterNormal(); orientation[ownScvf.directionIndex()] = 1.0; - return (tangentialVelocityGradient*distanceNormalToBoundary - + asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * betaBJ * distanceNormalToBoundary - + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0); + + asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * factor * distanceNormalToBoundary + + velocitySelf) / (factor*distanceNormalToBoundary + 1.0); } private: diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index bbd5e1981277bc3b39973cbc79608c72b1fa8246..42993bd320a5b5cdb77e3fe5d50e0bfb72c69746 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -310,9 +310,9 @@ public: lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralScvf)); } - // If the current scvf is on a bounary and if there is a Neumann or Beavers-Joseph-(Saffmann) BC for the stress in tangential direction, + // If the current scvf is on a bounary and if there is a Neumann or SlipCondition BC for the stress in tangential direction, // assign this value for the lateral momentum flux. No further calculations are required. We assume that all lateral faces - // have the same type of BC (Neumann or Beavers-Joseph-(Saffmann)), but we sample the value at their actual positions. + // have the same type of BC (Neumann or SlipCondition), but we sample the value at their actual positions. if (currentScvfBoundaryTypes) { // Handle Neumann BCs. @@ -331,8 +331,9 @@ public: // It is not clear how to evaluate the BJ condition here. // For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face. // TODO: We should clarify if this is the correct approach. - if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())) && lateralFaceBoundaryTypes && - lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex()))) + bool slipCondition = currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex())); + bool neumann = lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())); + if ( slipCondition && lateralFaceBoundaryTypes && neumann) { const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, @@ -375,7 +376,7 @@ public: std::bitset<3> admittableBcTypes; admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)); admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()))); - admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))); + admittableBcTypes.set(2, lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()))); if (admittableBcTypes.count() != 1) { DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf " @@ -502,9 +503,9 @@ private: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralFace.directionIndex())]; } - else if (bcTypes.isBeaversJoseph(Indices::velocity(lateralFace.directionIndex()))) + else if (bcTypes.isSlipCondition(Indices::velocity(lateralFace.directionIndex()))) { - return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, + return VelocityGradients::slipVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } else @@ -573,7 +574,8 @@ private: { if (!scvf.boundary() || currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) || - currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex()))) + currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralFace.directionIndex())) + ) { const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); // Account for the orientation of the staggered normal face's outer normal vector. diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh index bcae100b68b201298423b08b28f74d0c2b5c3cbb..30316395dd0acc2665ba19275a02486f61f2c603 100644 --- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh @@ -583,7 +583,7 @@ private: // Find out what boundary type is set on the lateral face const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)); - const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())); + const bool lateralFaceHasSlipCondition = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex())); const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())); const Scalar velocitySelf = faceVars.velocitySelf(); @@ -592,10 +592,9 @@ private: if (useZeroGradient) return velocitySelf; - if (lateralFaceHasBJS) - return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, + if (lateralFaceHasSlipCondition) + return VelocityGradients::slipVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - else if(lateralFaceHasDirichletVelocity) { // ________________ diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh index a9c3c452b31e9d5ed9ba383d29cf3a592bd87986..d1e347e5f276e01fd852368bfce628fb5d14d2ae 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -121,10 +121,9 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())]; } - else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))) - { - return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + else if (lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()))){ + return slipVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } else DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center()); @@ -200,10 +199,9 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())]; } - else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex()))) - { - return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + else if (currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex()))){ + return slipVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } else DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center()); @@ -218,10 +216,10 @@ public: } /*! - * \brief Returns the Beavers-Jospeh slip velocity for a scvf which lies on the boundary itself. + * \brief Returns the slip velocity for a scvf which lies on the boundary itself. * * \verbatim - * in.norm. B-J slip + * in.norm. slip * vel. vel. * ^ ^ * | | @@ -238,7 +236,7 @@ public: * */ template - static Scalar beaversJosephVelocityAtCurrentScvf(const Problem& problem, + static Scalar slipVelocityAtCurrentScvf(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, @@ -253,38 +251,48 @@ public: const auto tangentialVelocityGradient = [&]() { - // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for - // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face + // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a slip condition + // is set there, assume a tangential velocity gradient of zero along the lateral face // (towards the current scvf). - static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + static const bool unsymmetrizedGradientForBeaversJoseph = [&]() + { + const bool tmp = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + if (tmp) + { + std::cerr << std::endl << "Warning: You are using the deprecated parameter 'EnableUnsymmetrizedVelocityGradientForBeaversJoseph'. Use 'EnableUnsymmetrizedVelocityGradientForIC' instead." <(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); - if (unsymmetrizedGradientForBJ) + if (unsymmetrizedGradientForIC) return 0.0; if (lateralScvf.boundary()) { if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) || - lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))) + lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()))) return 0.0; } return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); }(); - return problem.beaversJosephVelocity(element, - fvGeometry.scv(scvf.insideScvIdx()), - lateralScvf, - scvf, /*on boundary*/ - innerLateralVelocity, - tangentialVelocityGradient); + return problem.slipVelocity(element, + fvGeometry.scv(scvf.insideScvIdx()), + lateralScvf, + scvf, /*on boundary*/ + innerLateralVelocity, + tangentialVelocityGradient); } /*! - * \brief Returns the Beavers-Jospeh slip velocity for a lateral scvf which lies on the boundary. + * \brief Returns the slip velocity for a lateral scvf which lies on the boundary. * * \verbatim - * B-J slip * boundary + * slip * boundary * ************** vel. ***** * scvf ---------##### ~~~~> :::: || and # staggered half-control-volume (own element) * | || | curr. :: @@ -298,7 +306,7 @@ public: * \endverbatim */ template - static Scalar beaversJosephVelocityAtLateralScvf(const Problem& problem, + static Scalar slipVelocityAtLateralScvf(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, @@ -313,31 +321,40 @@ public: const auto tangentialVelocityGradient = [&]() { - // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for - // the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face + // If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a slip condition + // is set there, assume a tangential velocity gradient of zero along the lateral face // (towards the current scvf). - static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - - if (unsymmetrizedGradientForBJ) + static const bool unsymmetrizedGradientForBeaversJoseph = [&]() + { + const bool tmp = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + if (tmp) + { + std::cout << std::endl << "Warning: You are using the deprecated parameter 'EnableUnsymmetrizedVelocityGradientForBeaversJoseph'. Use 'EnableUnsymmetrizedVelocityGradientForIC' instead." <(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); + if (unsymmetrizedGradientForIC) return 0.0; if (scvf.boundary()) { if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) || - currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex()))) + currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex()))) return 0.0; } return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); }(); - return problem.beaversJosephVelocity(element, - fvGeometry.scv(scvf.insideScvIdx()), - scvf, - lateralScvf, /*on boundary*/ - innerParallelVelocity, - tangentialVelocityGradient); + return problem.slipVelocity(element, + fvGeometry.scv(scvf.insideScvIdx()), + scvf, + lateralScvf, /*on boundary*/ + innerParallelVelocity, + tangentialVelocityGradient); } private: diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh index 1b5721d002b88025d55ba53f71f97a1be25684d4..c08f15d44e992060b6124550d1025c5d6a3daa63 100644 --- a/dumux/freeflow/rans/problem.hh +++ b/dumux/freeflow/rans/problem.hh @@ -324,10 +324,10 @@ public: } } - // Calculate the BJS-velocity by accounting for all sub faces. - std::vector bjsNumFaces(dim, 0); - std::vector bjsNeighbor(dim, 0); - DimVector bjsVelocityAverage(0.0); + // Calculate the slip-velocity by accounting for all sub faces. + std::vector slipNumFaces(dim, 0); + std::vector slipNeighbor(dim, 0); + DimVector slipVelocityAverage(0.0); DimVector normalNormCoordinate(0.0); unsigned int velIdx = Indices::velocity(scvfNormDim); const int numSubFaces = scvf.pairData().size(); @@ -335,33 +335,33 @@ public: { const auto& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localLateralFaceIdx); - // adapt calculations for Beavers-Joseph-Saffman condition + // adapt calculations for slip condition unsigned int normalNormDim = lateralFace.directionIndex(); - if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBeaversJoseph(Indices::velocity(velIdx)))) + if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isSlipCondition(Indices::velocity(velIdx)))) { unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0]; if (lateralFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim]) neighborIdx = neighborIdx_[elementIdx][normalNormDim][1]; const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx()); - bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0); - if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim]) + slipVelocityAverage[normalNormDim] += ParentType::slipVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0); + if (slipNumFaces[normalNormDim] > 0 && neighborIdx != slipNeighbor[normalNormDim]) DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur"); - bjsNeighbor[normalNormDim] = neighborIdx; + slipNeighbor[normalNormDim] = neighborIdx; normalNormCoordinate[normalNormDim] = lateralFace.center()[normalNormDim]; - bjsNumFaces[normalNormDim]++; + slipNumFaces[normalNormDim]++; } } for (unsigned dirIdx = 0; dirIdx < dim; ++dirIdx) { - if (bjsNumFaces[dirIdx] == 0) + if (slipNumFaces[dirIdx] == 0) continue; - unsigned int neighborIdx = bjsNeighbor[dirIdx]; - bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx]; + unsigned int neighborIdx = slipNeighbor[dirIdx]; + slipVelocityAverage[dirIdx] /= slipNumFaces[dirIdx]; velocityGradients_[elementIdx][velIdx][dirIdx] - = (velocity_[neighborIdx][velIdx] - bjsVelocityAverage[dirIdx]) + = (velocity_[neighborIdx][velIdx] - slipVelocityAverage[dirIdx]) / (cellCenter_[neighborIdx][dirIdx] - normalNormCoordinate[dirIdx]); } diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index b78243f07f1409f1e7491c851f67809c0b81d41f..fd09eaffaf0950aec86d56eed661ac214848e51d 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -28,7 +28,11 @@ #include #include -#include +#include //TODO: Why is this the right coupling manager, not the one in this folder? + +#include +#include +#include namespace Dumux { /*! @@ -58,6 +62,10 @@ class StokesDarcyCouplingDataBoxBase : public StokesDarcyCouplingDataImplementat static constexpr auto stokesIdx = CouplingManager::stokesIdx; static constexpr auto darcyIdx = CouplingManager::darcyIdx; + using VelocityVector = typename Element::Geometry::GlobalCoordinate; + template using BoundaryTypes = GetPropType, Properties::BoundaryTypes>; + using StokesVelocityGradients = StaggeredVelocityGradients, BoundaryTypes, Indices>; + using AdvectionType = GetPropType, Properties::AdvectionType>; using DarcysLaw = DarcysLawImplementation, GridGeometry::discMethod>; using ForchheimersLaw = ForchheimersLawImplementation, GridGeometry::discMethod>; @@ -72,8 +80,11 @@ public: /*! * \brief Returns the momentum flux across the coupling boundary. * + * Calculates the classical or new (normal)momentumCouplingCondition depending on the value of the parameter "Problem.NewIc" + * Defaults to classical momentumCouplingCondition. + * * For the normal momentum coupling, the porous medium side of the coupling condition - * is evaluated, i.e. -[p n]^pm. + * is evaluated, i.e. -[p n]^pm. For the new normal momentum coupling, the stokes term -viscosity*Nsbl*tau.T*grad(v_ff)*n is added. * */ template @@ -83,9 +94,12 @@ public: const ElementFaceVariables& stokesElemFaceVars, const SubControlVolumeFace& scvf) const { + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + Scalar momentumFlux(0.0); - const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); + //######## darcy contribution ################# + const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); // integrate darcy pressure over each coupling segment and average for (const auto& data : stokesContext) { @@ -123,10 +137,165 @@ public: if(getPropValue, Properties::NormalizePressure>()) momentumFlux -= this->couplingManager().problem(stokesIdx).initial(scvf)[Indices::pressureIdx]; - momentumFlux *= scvf.directionSign(); + if (newIc_) + { + //######## New stokes contribution ################# + static const bool unsymmetrizedGradientForBeaversJoseph = [&]() + { + const bool tmp = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + if (tmp) + { + std::cerr << "Warning: You are using the deprecated parameter 'EnableUnsymmetrizedVelocityGradientForBeaversJoseph'. Use 'EnableUnsymmetrizedVelocityGradientForIC' instead." << std::endl; + } + return tmp; + }(); + + // TODO: Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); + // NewIc not verified for symmetrized gradient + static bool once = [](){ + if(!unsymmetrizedGradientForIC) std::cerr << "Warning: The interface conditions for arbitrary flows to the interface are not verified for symmetrized stress tensors" <> currentScvfBoundaryTypes; + if (scvf.boundary()) + { + currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); + } + + std::optional> lateralFaceBoundaryTypes; + if (lateralScvf.boundary()) + { + lateralFaceBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, lateralScvf)); + } + + // Get velocity gradients + const Scalar velocityGrad_ji = StokesVelocityGradients::velocityGradJI( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + const Scalar velocityGrad_ij = unsymmetrizedGradientForIC ? 0.0 : StokesVelocityGradients::velocityGradIJ( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + + // Calculate stokes contribution to momentum flux: N_s^{bl} \tau T n + const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); + const Scalar viscosity = stokesElemVolVars[scvf.insideScvIdx()].viscosity(); + // Averaging the gradients over the subfaces to get evaluation at the center + momentumFlux -= 1.0/numSubFaces * viscosity * Nsbl * (velocityGrad_ji + velocityGrad_ij); + } + } + momentumFlux *= scvf.directionSign(); return momentumFlux; } + + /*! + * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law + * + * The tangential porous medium velocity needs to be evaluated for the classical and new tangential coupling (slipCondition) at the + * stokes-darcy interface, "Effective coupling conditions for arbitrary flows in Stokes-Darcy systems" by Elissa Eggenweiler. + * We use darcys law and perform an integral average over all coupling segments. + * + * Depending on the parameter "Problem.NewIc" the standard permeability tensor K + * or an altered permeability tensor M is used to evaluate the velocity. The method is using K per default. + * + */ + VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + + static constexpr int darcyDim = GridGeometry::GridView::dimension; + using JacobianType = Dune::FieldMatrix; + std::vector shapeDerivatives; + std::vector> shapeValues; + + VelocityVector velocity(0.0); // velocity darcy + VelocityVector gradP(0.0); // pressure gradient darcy + Scalar rho(0.0); // density darcy + Scalar intersectionLength = 0.0; // (total)intersection length, could differ from scfv length + + //Getting needed information from the darcy domain + const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); + static const bool enableGravity = getParamFromGroup(this->couplingManager().problem(darcyIdx).paramGroup(), "Problem.EnableGravity"); + + // Iteration over the different coupling segments + for (const auto& data : stokesContext) + { + //We are on (one of) the correct scvf(s) + if (scvf.index() == data.stokesScvfIdx) + { + const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); + const auto& elemVolVars = *(data.elementVolVars); + const auto& darcyFvGeometry = data.fvGeometry; + const auto& localBasis = darcyFvGeometry.feLocalBasis(); + + + // Darcy Permeability + const auto& K = data.volVars.permeability(); + + // INTEGRATION, second order as box provides linear functions + const auto& rule = Dune::QuadratureRules::rule(data.segmentGeometry.type(), 2); + //Loop over all quadrature points in the rule + for (const auto& qp : rule) + { + const auto& ipLocal = qp.position(); + const auto& ipGlobal = data.segmentGeometry.global(ipLocal); + const auto& ipElementLocal = data.element.geometry().local(ipGlobal); + + //reset pressure gradient and rho at this qp + gradP=0.0; + rho=0.0; + + //calculate the shape and derivative values at the qp + localBasis.evaluateFunction(ipElementLocal, shapeValues); + localBasis.evaluateJacobian(ipElementLocal, shapeDerivatives); + + //calc pressure gradient and rho at qp, every scv belongs to one node + for (const auto& scv : scvs(data.fvGeometry)){ + //gradP += p_i* (J^-T * L'_i) + data.element.geometry().jacobianInverseTransposed(ipElementLocal).usmv(elemVolVars[scv].pressure(darcyPhaseIdx), shapeDerivatives[scv.indexInElement()][0], gradP); + if (enableGravity) + { + rho += elemVolVars[scv].density(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; + } + } + //account for gravity + if (enableGravity) + { + gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); + } + + if(newIc_) + { + // darcy spatial dependent parameters + const auto& epsInterface = this->couplingManager().problem(darcyIdx).spatialParams().epsInterfaceAtPos(ipGlobal); + const auto& M = this->couplingManager().problem(darcyIdx).spatialParams().matrixNTangentialAtPos(ipGlobal); + //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*eps**2*M/mu*gradP + M.usmv(qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); + } + else + { + //Add the integrated segment velocity to the sum: v+= -weight_k * sqrt(det(A^T*A))*K/mu*gradP + K.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), gradP, velocity); + } + } + intersectionLength += data.segmentGeometry.volume(); + } + } + velocity /= intersectionLength; //averaging + return velocity; + } }; /*! diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh index 4bf72d67c66478338577e53cb23805c544967a7e..56181d1153703be82d09fb6ef828376dc55501c5 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh @@ -34,7 +34,8 @@ #include #include -#include "couplingdata.hh" +#include "../couplingmanager.hh" +#include "../couplingdata.hh" #include "couplingmapper.hh" namespace Dumux { @@ -45,18 +46,16 @@ namespace Dumux { */ template class StokesDarcyCouplingManagerImplementation -: public virtual StaggeredCouplingManager +: public virtual StaggeredCouplingManager, public FreeFlowPorousMediumCouplingManagerBase { using Scalar = typename MDTraits::Scalar; using ParentType = StaggeredCouplingManager; public: - static constexpr auto stokesCellCenterIdx = typename MDTraits::template SubDomain<0>::Index(); - static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<1>::Index(); - static constexpr auto cellCenterIdx = typename MDTraits::template SubDomain<0>::Index(); - static constexpr auto faceIdx = typename MDTraits::template SubDomain<1>::Index(); - static constexpr auto stokesIdx = stokesCellCenterIdx; - static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index(); + static constexpr auto freeFlowFaceIdx = FreeFlowPorousMediumCouplingManagerBase::freeFlowFaceIdx; + static constexpr auto freeFlowCellCenterIdx = FreeFlowPorousMediumCouplingManagerBase::freeFlowCellCenterIdx; + static constexpr auto freeFlowIdx = FreeFlowPorousMediumCouplingManagerBase::freeFlowIdx; + static constexpr auto porousMediumIdx = FreeFlowPorousMediumCouplingManagerBase::porousMediumIdx; private: @@ -92,32 +91,31 @@ private: using CellCenterSolutionVector = GetPropType; - using VelocityVector = typename Element::Geometry::GlobalCoordinate; + using VelocityVector = typename Element::Geometry::GlobalCoordinate; using CouplingMapper = StokesDarcyCouplingMapperBox; // Each context object contains the data related to one coupling segment struct StationaryStokesCouplingContext { - Element element; - FVElementGeometry fvGeometry; + Element element; + FVElementGeometry fvGeometry; std::size_t darcyScvfIdx; std::size_t stokesScvfIdx; - VolumeVariables volVars; - std::unique_ptr< ElementVolumeVariables > elementVolVars; - std::unique_ptr< ElementFluxVariablesCache > elementFluxVarsCache; + VolumeVariables volVars; + std::unique_ptr< ElementVolumeVariables > elementVolVars; + std::unique_ptr< ElementFluxVariablesCache > elementFluxVarsCache; typename CouplingMapper::CouplingSegment::Geometry segmentGeometry; }; - // Each context object contains the data related to one coupling segment struct StationaryDarcyCouplingContext { - Element element; - FVElementGeometry fvGeometry; + Element element; + FVElementGeometry fvGeometry; std::size_t stokesScvfIdx; std::size_t darcyScvfIdx; VelocityVector velocity; - VolumeVariables volVars; + VolumeVariables volVars; typename CouplingMapper::CouplingSegment::Geometry segmentGeometry; }; public: @@ -127,8 +125,8 @@ public: using CouplingData = StokesDarcyCouplingData>; //! Constructor - StokesDarcyCouplingManagerImplementation(std::shared_ptr> stokesFvGridGeometry, - std::shared_ptr> darcyFvGridGeometry) + StokesDarcyCouplingManagerImplementation(std::shared_ptr> stokesFvGridGeometry, + std::shared_ptr> darcyFvGridGeometry) { } /*! @@ -137,8 +135,8 @@ public: // \{ //! Initialize the coupling manager - void init(std::shared_ptr> stokesProblem, - std::shared_ptr> darcyProblem, + void init(std::shared_ptr> stokesProblem, + std::shared_ptr> darcyProblem, const SolutionVector& curSol) { if (Dune::FloatCmp::ne(stokesProblem->gravity(), darcyProblem->spatialParams().gravity({}))) @@ -192,12 +190,12 @@ public: /*! * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Darcy information) */ - template = 0> - void bindCouplingContext(Dune::index_constant domainI, const Element& element, const Assembler& assembler) const + template = 0> + void bindCouplingContext(Dune::index_constant domainI, const Element& element, const Assembler& assembler) const { stokesCouplingContext_.clear(); - const auto stokesElementIdx = this->problem(stokesIdx).gridGeometry().elementMapper().index(element); + const auto stokesElementIdx = this->problem(freeFlowIdx).gridGeometry().elementMapper().index(element); boundStokesElemIdx_ = stokesElementIdx; // do nothing if the element is not coupled to the other domain @@ -206,30 +204,30 @@ public: // prepare the coupling context const auto& couplingSegments = couplingMapper_.stokesElementToDarcyElementMap().at(stokesElementIdx); - auto darcyFvGeometry = localView(this->problem(darcyIdx).gridGeometry()); + auto darcyFvGeometry = localView(this->problem(porousMediumIdx).gridGeometry()); for(const auto& couplingSegment : couplingSegments) { - const auto& darcyElement = this->problem(darcyIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx); + const auto& darcyElement = this->problem(porousMediumIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx); darcyFvGeometry.bind(darcyElement); const auto& scvf = darcyFvGeometry.scvf(couplingSegment.scvfIdx); const auto& scv = darcyFvGeometry.scv(scvf.insideScvIdx()); - auto darcyElemVolVars = localView(assembler.gridVariables(darcyIdx).curGridVolVars()); - auto darcyElemFluxVarsCache = localView(assembler.gridVariables(darcyIdx).gridFluxVarsCache()); + auto darcyElemVolVars = localView(assembler.gridVariables(porousMediumIdx).curGridVolVars()); + auto darcyElemFluxVarsCache = localView(assembler.gridVariables(porousMediumIdx).gridFluxVarsCache()); - darcyElemVolVars.bind(darcyElement, darcyFvGeometry, this->curSol()[darcyIdx]); - darcyElemVolVars.bindElement(darcyElement, darcyFvGeometry, this->curSol()[darcyIdx]); + darcyElemVolVars.bind(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]); + darcyElemVolVars.bindElement(darcyElement, darcyFvGeometry, this->curSol()[porousMediumIdx]); darcyElemFluxVarsCache.bind(darcyElement, darcyFvGeometry, darcyElemVolVars); - const auto darcyElemSol = elementSolution(darcyElement, this->curSol()[darcyIdx], this->problem(darcyIdx).gridGeometry()); - VolumeVariables darcyVolVars; - darcyVolVars.update(darcyElemSol, this->problem(darcyIdx), darcyElement, scv); + const auto darcyElemSol = elementSolution(darcyElement, this->curSol()[porousMediumIdx], this->problem(porousMediumIdx).gridGeometry()); + VolumeVariables darcyVolVars; + darcyVolVars.update(darcyElemSol, this->problem(porousMediumIdx), darcyElement, scv); // add the context stokesCouplingContext_.push_back({darcyElement, darcyFvGeometry, couplingSegment.scvfIdx, couplingSegment.flipScvfIdx, darcyVolVars, - std::make_unique< ElementVolumeVariables >( std::move(darcyElemVolVars)), - std::make_unique< ElementFluxVariablesCache >( std::move(darcyElemFluxVarsCache)), + std::make_unique< ElementVolumeVariables >( std::move(darcyElemVolVars)), + std::make_unique< ElementFluxVariablesCache >( std::move(darcyElemFluxVarsCache)), couplingSegment.geometry}); } } @@ -238,11 +236,11 @@ public: * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) */ template - void bindCouplingContext(Dune::index_constant domainI, const Element& element, const Assembler& assembler) const + void bindCouplingContext(Dune::index_constant domainI, const Element& element, const Assembler& assembler) const { darcyCouplingContext_.clear(); - const auto darcyElementIdx = this->problem(darcyIdx).gridGeometry().elementMapper().index(element); + const auto darcyElementIdx = this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element); boundDarcyElemIdx_ = darcyElementIdx; // do nothing if the element is not coupled to the other domain @@ -251,11 +249,11 @@ public: // prepare the coupling context const auto& couplingSegments = couplingMapper_.darcyElementToStokesElementMap().at(darcyElementIdx); - auto stokesFvGeometry = localView(this->problem(stokesIdx).gridGeometry()); + auto stokesFvGeometry = localView(this->problem(freeFlowIdx).gridGeometry()); for(const auto& couplingSegment : couplingSegments) { - const auto& stokesElement = this->problem(stokesIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx); + const auto& stokesElement = this->problem(freeFlowIdx).gridGeometry().boundingBoxTree().entitySet().entity(couplingSegment.eIdx); stokesFvGeometry.bindElement(stokesElement); VelocityVector faceVelocity(0.0); @@ -263,16 +261,16 @@ public: for(const auto& scvf : scvfs(stokesFvGeometry)) { if(scvf.index() == couplingSegment.scvfIdx) - faceVelocity[scvf.directionIndex()] = this->curSol()[stokesFaceIdx][scvf.dofIndex()]; + faceVelocity[scvf.directionIndex()] = this->curSol()[freeFlowFaceIdx][scvf.dofIndex()]; } - using PriVarsType = typename VolumeVariables::PrimaryVariables; - const auto& cellCenterPriVars = this->curSol()[stokesCellCenterIdx][couplingSegment.eIdx]; + using PriVarsType = typename VolumeVariables::PrimaryVariables; + const auto& cellCenterPriVars = this->curSol()[freeFlowCellCenterIdx][couplingSegment.eIdx]; const auto elemSol = makeElementSolutionFromCellCenterPrivars(cellCenterPriVars); - VolumeVariables stokesVolVars; + VolumeVariables stokesVolVars; for(const auto& scv : scvs(stokesFvGeometry)) - stokesVolVars.update(elemSol, this->problem(stokesIdx), stokesElement, scv); + stokesVolVars.update(elemSol, this->problem(freeFlowIdx), stokesElement, scv); // add the context darcyCouplingContext_.push_back({stokesElement, stokesFvGeometry, @@ -285,11 +283,11 @@ public: * \brief Update the coupling context for the Darcy residual w.r.t. Darcy DOFs */ template - void updateCouplingContext(Dune::index_constant domainI, + void updateCouplingContext(Dune::index_constant domainI, const LocalAssemblerI& localAssemblerI, - Dune::index_constant domainJ, + Dune::index_constant domainJ, std::size_t dofIdxGlobalJ, - const PrimaryVariables& priVarsJ, + const PrimaryVariables& priVarsJ, int pvIdxJ) { this->curSol()[domainJ][dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; @@ -299,27 +297,27 @@ public: * \brief Update the coupling context for the Darcy residual w.r.t. the Stokes cell-center DOFs (DarcyToCC) */ template - void updateCouplingContext(Dune::index_constant domainI, + void updateCouplingContext(Dune::index_constant domainI, const LocalAssemblerI& localAssemblerI, - Dune::index_constant domainJ, + Dune::index_constant domainJ, const std::size_t dofIdxGlobalJ, - const PrimaryVariables& priVars, + const PrimaryVariables& priVars, int pvIdxJ) { this->curSol()[domainJ][dofIdxGlobalJ] = priVars; for (auto& data : darcyCouplingContext_) { - const auto stokesElemIdx = this->problem(stokesIdx).gridGeometry().elementMapper().index(data.element); + const auto stokesElemIdx = this->problem(freeFlowIdx).gridGeometry().elementMapper().index(data.element); if(stokesElemIdx != dofIdxGlobalJ) continue; - using PriVarsType = typename VolumeVariables::PrimaryVariables; + using PriVarsType = typename VolumeVariables::PrimaryVariables; const auto elemSol = makeElementSolutionFromCellCenterPrivars(priVars); for(const auto& scv : scvs(data.fvGeometry)) - data.volVars.update(elemSol, this->problem(stokesIdx), data.element, scv); + data.volVars.update(elemSol, this->problem(freeFlowIdx), data.element, scv); } } @@ -327,9 +325,9 @@ public: * \brief Update the coupling context for the Darcy residual w.r.t. the Stokes face DOFs (DarcyToFace) */ template - void updateCouplingContext(Dune::index_constant domainI, + void updateCouplingContext(Dune::index_constant domainI, const LocalAssemblerI& localAssemblerI, - Dune::index_constant domainJ, + Dune::index_constant domainJ, const std::size_t dofIdxGlobalJ, const PrimaryVariables<1>& priVars, int pvIdxJ) @@ -349,12 +347,12 @@ public: /*! * \brief Update the coupling context for the Stokes cc residual w.r.t. the Darcy DOFs (FaceToDarcy) */ - template = 0> + template = 0> void updateCouplingContext(Dune::index_constant domainI, const LocalAssemblerI& localAssemblerI, - Dune::index_constant domainJ, + Dune::index_constant domainJ, const std::size_t dofIdxGlobalJ, - const PrimaryVariables& priVars, + const PrimaryVariables& priVars, int pvIdxJ) { this->curSol()[domainJ][dofIdxGlobalJ] = priVars; @@ -362,19 +360,19 @@ public: for (auto& data : stokesCouplingContext_) { //Derivatives are assumed to be always calculated with respect to unknowns associated with its own element - const auto darcyElemSol = elementSolution(data.element, this->curSol()[darcyIdx], this->problem(darcyIdx).gridGeometry()); + const auto darcyElemSol = elementSolution(data.element, this->curSol()[porousMediumIdx], this->problem(porousMediumIdx).gridGeometry()); const auto& scvf = data.fvGeometry.scvf(data.darcyScvfIdx); const auto& scv = data.fvGeometry.scv(scvf.insideScvIdx()); if(scv.dofIndex() == dofIdxGlobalJ) - data.volVars.update(darcyElemSol, this->problem(darcyIdx), data.element, scv); + data.volVars.update(darcyElemSol, this->problem(porousMediumIdx), data.element, scv); for (auto&& scv : scvs(data.fvGeometry)) { if(scv.dofIndex() == dofIdxGlobalJ) { auto& volVars = (*data.elementVolVars)[scv]; - volVars.update(darcyElemSol, this->problem(darcyIdx), data.element, scv); + volVars.update(darcyElemSol, this->problem(porousMediumIdx), data.element, scv); } } } @@ -393,7 +391,7 @@ public: /*! * \brief Access the coupling context needed for the Stokes domain */ - const auto& stokesCouplingContext(const Element& element, const SubControlVolumeFace& scvf) const + const auto& stokesCouplingContext(const Element& element, const SubControlVolumeFace& scvf) const { if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx()) DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); @@ -410,9 +408,9 @@ public: /*! * \brief Access the coupling context needed for the Darcy domain */ - const auto& darcyCouplingContext(const Element& element, const SubControlVolumeFace& scvf) const + const auto& darcyCouplingContext(const Element& element, const SubControlVolumeFace& scvf) const { - if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(darcyIdx).gridGeometry().elementMapper().index(element)) + if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element)) DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); for(const auto& context : darcyCouplingContext_) @@ -427,9 +425,9 @@ public: /*! * \brief Access the coupling context needed for the Darcy domain */ - const auto& darcyCouplingContextVector(const Element& element, const SubControlVolumeFace& scvf) const + const auto& darcyCouplingContextVector(const Element& element, const SubControlVolumeFace& scvf) const { - if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(darcyIdx).gridGeometry().elementMapper().index(element)) + if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element)) DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); return darcyCouplingContext_; @@ -438,7 +436,7 @@ public: /*! * \brief Access the coupling context needed for the Stokes domain */ - const auto& stokesCouplingContextVector(const Element& element, const SubControlVolumeFace& scvf) const + const auto& stokesCouplingContextVector(const Element& element, const SubControlVolumeFace& scvf) const { if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx()) DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); @@ -453,9 +451,9 @@ public: /*! * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs */ - const CouplingStencil& couplingStencil(Dune::index_constant domainI, - const Element& element, - Dune::index_constant domainJ) const + const CouplingStencil& couplingStencil(Dune::index_constant domainI, + const Element& element, + Dune::index_constant domainJ) const { const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); if(stokesCellCenterCouplingStencils_.count(eIdx)) @@ -478,9 +476,9 @@ public: * \brief The coupling stencil of domain I, i.e. which domain J dofs * the given domain I element's residual depends on. */ - const CouplingStencil& couplingStencil(Dune::index_constant domainI, - const Element& element, - Dune::index_constant domainJ) const + const CouplingStencil& couplingStencil(Dune::index_constant domainI, + const Element& element, + Dune::index_constant domainJ) const { const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); if(darcyToStokesCellCenterCouplingStencils_.count(eIdx)) @@ -493,9 +491,9 @@ public: * \brief The coupling stencil of domain I, i.e. which domain J dofs * the given domain I element's residual depends on. */ - const CouplingStencil& couplingStencil(Dune::index_constant domainI, - const Element& element, - Dune::index_constant domainJ) const + const CouplingStencil& couplingStencil(Dune::index_constant domainI, + const Element& element, + Dune::index_constant domainJ) const { const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); if (darcyToStokesFaceCouplingStencils_.count(eIdx)) @@ -510,16 +508,16 @@ public: */ template const CouplingStencil& couplingStencil(Dune::index_constant domainI, - const SubControlVolumeFace& scvf, + const SubControlVolumeFace& scvf, Dune::index_constant domainJ) const { return emptyStencil_; } /*! * \brief The coupling stencil of a Stokes face w.r.t. Darcy DOFs */ - const CouplingStencil& couplingStencil(Dune::index_constant domainI, - const SubControlVolumeFace& scvf, - Dune::index_constant domainJ) const + const CouplingStencil& couplingStencil(Dune::index_constant domainI, + const SubControlVolumeFace& scvf, + Dune::index_constant domainJ) const { const auto faceDofIdx = scvf.dofIndex(); if(stokesFaceCouplingStencils_.count(faceDofIdx)) @@ -547,7 +545,7 @@ public: /*! * \brief Returns whether a given free flow scvf is coupled to the other domain */ - bool isCoupledEntity(Dune::index_constant, const SubControlVolumeFace& scvf) const + bool isCoupledEntity(Dune::index_constant, const SubControlVolumeFace& scvf) const { return stokesFaceCouplingStencils_.count(scvf.dofIndex()); } @@ -555,9 +553,9 @@ public: /*! * \brief Returns whether a given free flow scvf is coupled to the other domain */ - bool isCoupledEntity(Dune::index_constant, const Element& element, const SubControlVolumeFace& scvf) const + bool isCoupledEntity(Dune::index_constant, const Element& element, const SubControlVolumeFace& scvf) const { - return couplingMapper_.isCoupledDarcyScvf(this->problem(darcyIdx).gridGeometry().elementMapper().index(element), + return couplingMapper_.isCoupledDarcyScvf(this->problem(porousMediumIdx).gridGeometry().elementMapper().index(element), scvf.index()); } diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh index 68af4896cd7930826398ec4efed1dcddb8ff900e..3cf3e9b51cf07fdfc5d24b35396ff2cbdb380f22 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh @@ -197,7 +197,7 @@ public: typename IntersectionAlgorithm::Intersection rawIs; if(IntersectionAlgorithm::intersection(darcyScvfGeometry, stokesScvf.geometry(), rawIs)) { - const auto is = typename CouplingSegment::Geometry(Dune::GeometryTypes::simplex(dim), rawIs); + const auto is = typename CouplingSegment::Geometry(Dune::GeometryTypes::simplex(dim-1), rawIs); isCoupledDarcyScvf_[darcyEIdx][darcyScvf.index()] = true; darcyElementToStokesElementMap_[darcyEIdx].push_back({stokesEIdx, stokesScvf.index(), darcyScvf.index(), is}); stokesElementToDarcyElementMap_[stokesEIdx].push_back({darcyEIdx, darcyScvf.index(), stokesScvf.index(), is}); diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt index 45acb89ecd3ed8e3f2b58c3260ba52309dd9ee7f..84140c7d0e7c7e60fc3a16afe85d7165a6e95ca7 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(convergencetest) +add_subdirectory(bj_newic) add_input_file_links() diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/1pspatialparams.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/1pspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..a861e5f18e4b467b32f0e268a640df9ba773afbe --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/1pspatialparams.hh @@ -0,0 +1,127 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup OnePTests + * \brief The spatial parameters class for the test problem using the 1p cc model + */ +#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH +#define DUMUX_1P_TEST_SPATIALPARAMS_HH + +#include + +namespace Dumux { + +/*! + * \ingroup OnePModel + * + * \brief The spatial parameters class for the test problem using the + * 1p cc model + */ +template +class OnePSpatialParams +: public FVSpatialParamsOneP> +{ + using GridView = typename GridGeometry::GridView; + using ParentType = FVSpatialParamsOneP>; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + // export permeability type + static constexpr int dim = GridView::dimension; + using PermeabilityType = Dune::FieldMatrix; + + OnePSpatialParams(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) + { + K_ = PermeabilityType(0.0); + M_ = PermeabilityType(0.0); + + K_[0][0] = getParam("Darcy.SpatialParams.Permeability_xx"); + K_[0][1] = K_[1][0] = getParam("Darcy.SpatialParams.Permeability_xy"); + K_[1][1] = getParam("Darcy.SpatialParams.Permeability_yy"); + + porosity_ = getParam("Darcy.SpatialParams.Porosity"); + alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); + + //for new ic + epsInterface_ = getParam("Darcy.InterfaceParams.epsInterface"); + N_s_bl = getParam("Darcy.InterfaceParams.N_s_bl"); + N_1_bl = getParam("Darcy.InterfaceParams.N_1_bl"); + + M_[0][0]= getParam("Darcy.InterfaceParams.M_1_bl"); + M_[1][1]= getParam("Darcy.InterfaceParams.M_2_bl"); + } + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. + * + * \param globalPos The global position + * \return the intrinsic permeability + */ + PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const + { return K_; } + + /*! \brief Define the porosity in [-]. + * + * \param globalPos The global position + */ + Scalar porosityAtPos(const GlobalPosition& globalPos) const + { return porosity_; } + + /*! \brief Define the Beavers-Joseph coefficient in [-]. + * + * \param globalPos The global position + */ + Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const + { return alphaBJ_; } + + Scalar epsInterfaceAtPos(const GlobalPosition& globalPos) const + { return epsInterface_;} + + Scalar factorNMomentumAtPos(const GlobalPosition& globalPos) const + { return N_s_bl;} + + Scalar factorNTangentialAtPos(const GlobalPosition& globalPos) const + { return N_1_bl;} + + Dune::FieldMatrix matrixNTangentialAtPos(const GlobalPosition& globalPos) const + { return M_; } + + + +private: + PermeabilityType K_; + + Scalar porosity_; + Scalar alphaBJ_; + + //for new ic + Scalar epsInterface_; + Scalar N_s_bl; + Scalar N_1_bl; + Dune::FieldMatrix M_; + +}; + +} // end namespace + +#endif diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4c86072741f069566f46b4f8a3f933d462b056ea --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt @@ -0,0 +1,21 @@ +add_input_file_links() +dune_symlink_to_source_files(FILES "./../convergencetest/convergencetest.py") + +dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_bj + SOURCES main.cc + LABELS multidomain multidomain_boundary stokesdarcy + TIMEOUT 1000 + CMAKE_GUARD HAVE_UMFPACK + COMMAND ./convergencetest.py + CMD_ARGS test_md_boundary_darcy1p_stokes1p_bj_newic params_bj.input) + +dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_newic + SOURCES main.cc + LABELS multidomain multidomain_boundary stokesdarcy + TIMEOUT 1000 + CMAKE_GUARD HAVE_UMFPACK + COMMAND ./convergencetest.py + CMD_ARGS test_md_boundary_darcy1p_stokes1p_bj_newic params_newic.input) + + +add_executable(test_md_boundary_darcy1p_stokes1p_bj_newic EXCLUDE_FROM_ALL main.cc) \ No newline at end of file diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc new file mode 100644 index 0000000000000000000000000000000000000000..d74a5d46998edbd193795795391f7fb06221bbc8 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc @@ -0,0 +1,399 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief A test problem for the coupled Stokes/Darcy problem (1p). + */ + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include "problem_darcy.hh" +#include "problem_stokes.hh" + +namespace Dumux { +namespace Properties { + +template +struct CouplingManager +{ + using Traits = StaggeredMultiDomainTraits; + using type = Dumux::StokesDarcyCouplingManager; +}; + +template +struct CouplingManager +{ + using Traits = StaggeredMultiDomainTraits; + using type = Dumux::StokesDarcyCouplingManager; +}; + +} // end namespace Properties +} // end namespace Dumux + +/*! +* \brief Creates analytical solution. +* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces +* \param problem the problem for which to evaluate the analytical solution +*/ +template +auto createStokesAnalyticalSolution(const Problem& problem) +{ + const auto& gridGeometry = problem.gridGeometry(); + using GridView = typename std::decay_t::GridView; + + static constexpr auto dim = GridView::dimension; + static constexpr auto dimWorld = GridView::dimensionworld; + + using VelocityVector = Dune::FieldVector; + + std::vector analyticalPressure; + std::vector analyticalVelocity; + std::vector analyticalVelocityOnFace; + + analyticalPressure.resize(gridGeometry.numCellCenterDofs()); + analyticalVelocity.resize(gridGeometry.numCellCenterDofs()); + analyticalVelocityOnFace.resize(gridGeometry.numFaceDofs()); + + using Indices = typename Problem::Indices; + for (const auto& element : elements(gridGeometry.gridView())) + { + auto fvGeometry = localView(gridGeometry); + fvGeometry.bindElement(element); + for (auto&& scv : scvs(fvGeometry)) + { + auto ccDofIdx = scv.dofIndex(); + auto ccDofPosition = scv.dofPosition(); + auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition); + + // velocities on faces + for (auto&& scvf : scvfs(fvGeometry)) + { + const auto faceDofIdx = scvf.dofIndex(); + const auto faceDofPosition = scvf.center(); + const auto dirIdx = scvf.directionIndex(); + const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition); + analyticalVelocityOnFace[faceDofIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)]; + } + + analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx]; + + for (int dirIdx = 0; dirIdx < dim; ++dirIdx) + analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)]; + } + } + + return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace); +} + +/*! +* \brief Creates analytical solution. +* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces +* \param problem the problem for which to evaluate the analytical solution +*/ +template +auto createDarcyAnalyticalSolution(const Problem& problem) +{ + const auto& gridGeometry = problem.gridGeometry(); + using GridView = typename std::decay_t::GridView; + + static constexpr auto dim = GridView::dimension; + static constexpr auto dimWorld = GridView::dimensionworld; + + using VelocityVector = Dune::FieldVector; + + std::vector analyticalPressure; + std::vector analyticalVelocity; + + analyticalPressure.resize(gridGeometry.numDofs()); + analyticalVelocity.resize(gridGeometry.numDofs()); + + for (const auto& element : elements(gridGeometry.gridView())) + { + auto fvGeometry = localView(gridGeometry); + fvGeometry.bindElement(element); + for (auto&& scv : scvs(fvGeometry)) + { + const auto ccDofIdx = scv.dofIndex(); + const auto ccDofPosition = scv.dofPosition(); + const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition); + analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim]; + for (int dirIdx = 0; dirIdx < dim; ++dirIdx) + analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx]; + } + } + + return std::make_tuple(analyticalPressure, analyticalVelocity); +} + +template +void printStokesL2Error(const Problem& problem, const SolutionVector& x) +{ + using namespace Dumux; + using Scalar = double; + using TypeTag = Properties::TTag::StokesOneP; + using Indices = typename GetPropType::Indices; + using ModelTraits = GetPropType; + using PrimaryVariables = GetPropType; + + using L2Error = NavierStokesTestL2Error; + const auto l2error = L2Error::calculateL2Error(problem, x); + const int numCellCenterDofs = problem.gridGeometry().numCellCenterDofs(); + const int numFaceDofs = problem.gridGeometry().numFaceDofs(); + std::ostream tmpOutputObject(std::cout.rdbuf()); // create temporary output with fixed formatting without affecting std::cout + tmpOutputObject << std::setprecision(8) << "** L2 error (abs/rel) for " + << std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): " + << std::scientific + << "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx] + << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] + << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx] + << std::endl; + + // write the norm into a log file + std::ofstream logFile; + logFile.open(problem.name() + ".log", std::ios::app); + logFile << "[ConvergenceTest] L2(p) = " << l2error.first[Indices::pressureIdx] << " L2(vx) = " << l2error.first[Indices::velocityXIdx] << " L2(vy) = " << l2error.first[Indices::velocityYIdx] << std::endl; + logFile.close(); +} + +template +void printDarcyL2Error(const Problem& problem, const SolutionVector& x) +{ + using namespace Dumux; + using Scalar = double; + + Scalar l2error = 0.0; + + for (const auto& element : elements(problem.gridGeometry().gridView())) + { + auto fvGeometry = localView(problem.gridGeometry()); + fvGeometry.bindElement(element); + + for (auto&& scv : scvs(fvGeometry)) + { + const auto dofIdx = scv.dofIndex(); + const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.dofPosition())[2/*pressureIdx*/]; + l2error += scv.volume()*(delta*delta); + } + } + using std::sqrt; + l2error = sqrt(l2error); + + const auto numDofs = problem.gridGeometry().numDofs(); + std::ostream tmp(std::cout.rdbuf()); + tmp << std::setprecision(8) << "** L2 error (abs) for " + << std::setw(6) << numDofs << " cc dofs " + << std::scientific + << "L2 error = " << l2error + << std::endl; + + // write the norm into a log file + std::ofstream logFile; + logFile.open(problem.name() + ".log", std::ios::app); + logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl; + logFile.close(); +} + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv); + + // Define the sub problem type tags + using StokesTypeTag = Properties::TTag::StokesOneP; + using DarcyTypeTag = Properties::TTag::DarcyOneP; + + // try to create a grid (from the given grid file or the input file) + // for both sub-domains + using DarcyGridManager = Dumux::GridManager>; + DarcyGridManager darcyGridManager; + darcyGridManager.init("Darcy"); // pass parameter group + + using StokesGridManager = Dumux::GridManager>; + StokesGridManager stokesGridManager; + stokesGridManager.init("Stokes"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + const auto& stokesGridView = stokesGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using StokesGridGeometry = GetPropType; + auto stokesGridGeometry = std::make_shared(stokesGridView); + stokesGridGeometry->update(); + using DarcyGridGeometry = GetPropType; + auto darcyGridGeometry = std::make_shared(darcyGridView); + darcyGridGeometry->update(); + + using Traits = StaggeredMultiDomainTraits; + + // the coupling manager + using CouplingManager = StokesDarcyCouplingManager; + auto couplingManager = std::make_shared(stokesGridGeometry, darcyGridGeometry); + + // the indices + constexpr auto freeFlowCellCenterIdx = CouplingManager::freeFlowCellCenterIdx; + constexpr auto freeFlowFaceIdx = CouplingManager::freeFlowFaceIdx; + constexpr auto porousMediumIdx = CouplingManager::porousMediumIdx; + + // the problem (initial and boundary conditions) + using StokesProblem = GetPropType; + auto stokesProblem = std::make_shared(stokesGridGeometry, couplingManager); + using DarcyProblem = GetPropType; + auto darcyProblem = std::make_shared(darcyGridGeometry, couplingManager); + + // the solution vector + Traits::SolutionVector sol; + sol[freeFlowCellCenterIdx].resize(stokesGridGeometry->numCellCenterDofs()); + sol[freeFlowFaceIdx].resize(stokesGridGeometry->numFaceDofs()); + sol[porousMediumIdx].resize(darcyGridGeometry->numDofs()); + + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, freeFlowFaceIdx, freeFlowCellCenterIdx); + + couplingManager->init(stokesProblem, darcyProblem, sol); + + // the grid variables + using StokesGridVariables = GetPropType; + auto stokesGridVariables = std::make_shared(stokesProblem, stokesGridGeometry); + stokesGridVariables->init(stokesSol); + using DarcyGridVariables = GetPropType; + auto darcyGridVariables = std::make_shared(darcyProblem, darcyGridGeometry); + darcyGridVariables->init(sol[porousMediumIdx]); + + // intialize the vtk output module + using Scalar = typename Traits::Scalar; + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + GetPropType::initOutputModule(stokesVtkWriter); + const auto stokesAnalyticalSolution = createStokesAnalyticalSolution(*stokesProblem); + stokesVtkWriter.addField(std::get<0>(stokesAnalyticalSolution), "pressureExact"); + stokesVtkWriter.addField(std::get<1>(stokesAnalyticalSolution), "velocityExact"); + stokesVtkWriter.addFaceField(std::get<2>(stokesAnalyticalSolution), "faceVelocityExact"); + stokesVtkWriter.write(0.0); + + VtkOutputModule> darcyVtkWriter(*darcyGridVariables, sol[porousMediumIdx], darcyProblem->name()); + using DarcyVelocityOutput = GetPropType; + darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); + GetPropType::initOutputModule(darcyVtkWriter); + const auto darcyAnalyticalSolution = createDarcyAnalyticalSolution(*darcyProblem); + darcyVtkWriter.addField(std::get<0>(darcyAnalyticalSolution), "pressureExact"); + darcyVtkWriter.addField(std::get<1>(darcyAnalyticalSolution), "velocityExact"); + darcyVtkWriter.write(0.0); + + // the assembler for a stationary problem + using Assembler = MultiDomainFVAssembler; + auto assembler = std::make_shared(std::make_tuple(stokesProblem, stokesProblem, darcyProblem), + std::make_tuple(stokesGridGeometry->faceFVGridGeometryPtr(), + stokesGridGeometry->cellCenterFVGridGeometryPtr(), + darcyGridGeometry), + std::make_tuple(stokesGridVariables->faceGridVariablesPtr(), + stokesGridVariables->cellCenterGridVariablesPtr(), + darcyGridVariables), + couplingManager); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + + // the non-linear solver + using NewtonSolver = MultiDomainNewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); + + // solve the non-linear system + nonLinearSolver.solve(sol); + + // write vtk output + stokesVtkWriter.write(1.0); + darcyVtkWriter.write(1.0); + + printStokesL2Error(*stokesProblem, stokesSol); + printDarcyL2Error(*darcyProblem, sol[porousMediumIdx]); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_bj.input b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_bj.input new file mode 100644 index 0000000000000000000000000000000000000000..431620690ca4d001acd53a22c15ed8fccf73b34a --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_bj.input @@ -0,0 +1,47 @@ +[Darcy.Grid] +UpperRight = 1 1 +Cells = 20 20 + +[Stokes.Grid] +LowerLeft = 0 1 +UpperRight = 1 2 +Cells = 20 20 + +[Stokes.Problem] +Name = stokes + +[Darcy.Problem] +Name = darcy + +[Darcy.SpatialParams] +Porosity = 1.0 +Permeability_xx = 1.0 # m^2 +Permeability_xy = 2.0 # m^2 +Permeability_yy = 1.0 # m^2 +AlphaBeaversJoseph = 1.0 + +#remove later +[Darcy.InterfaceParams] +epsInterface = 1.0 +N_s_bl = 0.5 +N_1_bl = -3.0 +M_1_bl = 1.0 +M_2_bl = 0.0 + +[Vtk] +OutputName = test_md_boundary_stokes1p_darcy1p_BJ + +[Problem] +EnableGravity = false +EnableInertiaTerms = false +TestCase = BJSymmetrized + +[Vtk] +AddVelocity = 1 + +[Component] +LiquidDensity = 1.0 +LiquidKinematicViscosity = 1.0 + +[Assembly] +NumericDifference.BaseEpsilon = 1e-4 diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input new file mode 100644 index 0000000000000000000000000000000000000000..2e0495a12c2bdde25876af1b0de124ed4bacf9c0 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input @@ -0,0 +1,49 @@ +[Darcy.Grid] +UpperRight = 1 1 +Cells = 20 20 + +[Stokes.Grid] +LowerLeft = 0 1 +UpperRight = 1 2 +Cells = 20 20 + +[Stokes.Problem] +Name = stokes + +[Darcy.Problem] +Name = darcy + +[Darcy.SpatialParams] +Porosity = 1.0 +Permeability_xx = 1.0 # m^2 +Permeability_xy = 0.0 # m^2 +Permeability_yy = 1.0 # m^2 +# remove later! +AlphaBeaversJoseph = 1.0 + +[Darcy.InterfaceParams] +epsInterface = 1.0 +N_s_bl = 0.5 +N_1_bl = -3.0 +M_1_bl = 1.0 +M_2_bl = 1.0 + +[Vtk] +OutputName = test_md_boundary_stokes1p_darcy1p_newic + +[Problem] +EnableGravity = false +EnableInertiaTerms = false +NewIc = true +TestCase = NewICNonSymmetrized + +[Stokes.FreeFlow] +EnableUnsymmetrizedVelocityGradientForIC = true +EnableUnsymmetrizedVelocityGradient = true + +[Vtk] +AddVelocity = 1 + +[Component] +LiquidDensity = 1.0 +LiquidKinematicViscosity = 2.0 \ No newline at end of file diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh new file mode 100644 index 0000000000000000000000000000000000000000..23a536e857bd295586e79b4b97886d7bade4587d --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh @@ -0,0 +1,365 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief The Darcy sub-problem of coupled Stokes-Darcy BJ/NewIC test + */ + +#ifndef DUMUX_DARCY_SUBPROBLEM_HH +#define DUMUX_DARCY_SUBPROBLEM_HH + +#include +#include + +#include + +#include +#include + +#include "1pspatialparams.hh" + +#include +#include + +namespace Dumux { +template +class DarcySubProblem; + +namespace Properties { +// Create new type tags +namespace TTag { +struct DarcyOneP { using InheritsFrom = std::tuple; }; +} // end namespace TTag + +// Set the problem property +template +struct Problem { using type = Dumux::DarcySubProblem; }; + +// the fluid system +template +struct FluidSystem +{ + using Scalar = GetPropType; + using type = FluidSystems::OnePLiquid > ; +}; + +// Set the grid type +template +struct Grid { using type = Dune::YaspGrid<2>; }; + +template +struct SpatialParams +{ + using GridGeometry = GetPropType; + using Scalar = GetPropType; + using type = OnePSpatialParams; +}; +} // end namespace Properties + +/*! + * \ingroup BoundaryTests + * \brief The Darcy sub-problem of coupled Stokes-Darcy BJ/NewIC test + */ +template +class DarcySubProblem : public PorousMediumFlowProblem +{ + using ParentType = PorousMediumFlowProblem; + using GridView = typename GetPropType::GridView; + using Scalar = GetPropType; + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using BoundaryTypes = GetPropType; + using VolumeVariables = GetPropType; + using FVElementGeometry = typename GetPropType::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using GridGeometry = GetPropType; + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using CouplingManager = GetPropType; + + static constexpr auto velocityXIdx = 0; + static constexpr auto velocityYIdx = 1; + static constexpr auto pressureIdx = 2; + + enum class TestCase + { + BJSymmetrized, NewICNonSymmetrized + }; + +public: + //! export the Indices + using Indices = typename GetPropType::Indices; + + DarcySubProblem(std::shared_ptr gridGeometry, + std::shared_ptr couplingManager) + : ParentType(gridGeometry, "Darcy"), couplingManager_(couplingManager) + { + problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + const auto testCaseInput = getParamFromGroup(this->paramGroup(), "Problem.TestCase", "BJSymmetrized"); + if (testCaseInput == "BJSymmetrized") + testCase_ = TestCase::BJSymmetrized; + else if (testCaseInput == "NewICNonSymmetrized") + testCase_ = TestCase::NewICNonSymmetrized; + else + DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case"); + } + + /*! + * \brief The problem name. + */ + const std::string& name() const + { + return problemName_; + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the temperature within the domain in [K]. + * + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary control volume. + * + * \param element The element + * \param scv The boundary sub control volume + */ + BoundaryTypes boundaryTypes(const Element &element, const SubControlVolume &scv) const + { + BoundaryTypes values; + values.setAllDirichlet(); + + //corners will not be handled as coupled + if(onLeftBoundary_(scv.dofPosition()) || onRightBoundary_(scv.dofPosition())) + { + values.setAllDirichlet(); + return values; + } + + + auto fvGeometry = localView(this->gridGeometry()); + fvGeometry.bindElement(element); + for (auto&& scvf : scvfs(fvGeometry)) + { + if (couplingManager().isCoupledEntity(CouplingManager::porousMediumIdx, element, scvf)) + { + values.setAllCouplingNeumann(); + } + } + return values; + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet control volume. + * + * \param globalPos The position for which the Dirichlet value is set + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const + { + const auto p = analyticalSolution(globalPos)[pressureIdx]; + return PrimaryVariables(p); + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeometry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFluxVarsCache Flux variables caches for all faces in stencil + * \param scvf The boundary sub control volume face + * + * For this method, the \a values variable stores primary variables. + */ + template + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFluxVarsCache& elemFluxVarsCache, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if (couplingManager().isCoupledEntity(CouplingManager::porousMediumIdx, element, scvf)) + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); + + return values; + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + /*! + * \brief Evaluate the source term for all phases within a given + * sub-control-volume. + * + * \param element The element for which the source term is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param scv The subcontrolvolume + */ + template + NumEqVector source(const Element &element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const + { + GlobalPosition globalPos = scv.geometry().center(); + switch (testCase_) + { + case TestCase::BJSymmetrized: + return rhsBJSymmetrized_(globalPos); + case TestCase::NewICNonSymmetrized: + return rhsNewICNonSymmetrized_(globalPos); + default: + DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); + } + } + + // \} + + /*! + * \brief Evaluates the initial value for a control volume. + * + * \param globalPos The global position + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + return PrimaryVariables(0.0); + } + + /*! + * \brief Returns the analytical solution of the problem at a given position. + * + * \param globalPos The global position + */ + auto analyticalSolution(const GlobalPosition& globalPos) const + { + switch (testCase_) + { + case TestCase::BJSymmetrized: + return analyticalSolutionBJSymmetrized_(globalPos); + case TestCase::NewICNonSymmetrized: + return analyticalSolutionNewICNonSymmetrized_(globalPos); + default: + DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); + } + } + + // \} + + //! Set the coupling manager + void setCouplingManager(std::shared_ptr cm) + { couplingManager_ = cm; } + + //! Get the coupling manager + const CouplingManager& couplingManager() const + { return *couplingManager_; } + +private: + + // see exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler) + Dune::FieldVector analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const + { + Dune::FieldVector sol(0.0); + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + + using std::exp; using std::sin; using std::cos; + sol[velocityXIdx] = 2.0*x*x + 2.0*x*y - 6.0*x - 2.0*y*y + 3.0*y - 5.0; + sol[velocityYIdx] = 1.0*x*x + 4.0*x*y - 9.0*x - 1.0*y*y - 1.0; + sol[pressureIdx] = -x*(x-1.0)*(y-1.0) + 1.0/3.0*(y-1.0)*(y-1.0)*(y-1.0) + 4.0 + 2.0*y + x*x; + return sol; + } + + // see exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler) + NumEqVector rhsBJSymmetrized_(const GlobalPosition& globalPos) const + { + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + using std::sin; + return NumEqVector(8.0*x - 6.0); + } + + + // see exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler) + Dune::FieldVector analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const + { + Dune::FieldVector sol(0.0); + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + + using std::exp; using std::sin; using std::cos; + sol[velocityXIdx] = x*(y-1.0) + (x-1.0)*(y-1.0) - 3.0; + sol[velocityYIdx] = x*(x-1.0) - (y-1.0)*(y-1.0) - 2.0; + sol[pressureIdx] = x*(1.0-x)*(y-1.0) + 1.0/3.0*(y-1.0)*(y-1.0)*(y-1.0) + 3.0*x + 2.0*y + 1.0; + return sol; + } + + // see exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler) + NumEqVector rhsNewICNonSymmetrized_(const GlobalPosition& globalPos) const + { + // const Scalar x = globalPos[0]; + // const Scalar y = globalPos[1]; + // using std::sin; + return NumEqVector(0.0); + } + + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; } + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; } + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; } + + static constexpr Scalar eps_ = 1e-7; + std::shared_ptr couplingManager_; + std::string problemName_; + TestCase testCase_; +}; +} // end namespace Dumux + +#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh new file mode 100644 index 0000000000000000000000000000000000000000..219c1d3a0048cd0ba7f86392cf815909e46ec453 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh @@ -0,0 +1,400 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief The Stokes sub-problem of coupled Stokes-Darcy BJ/NewIC test + */ + +#ifndef DUMUX_STOKES_SUBPROBLEM_HH +#define DUMUX_STOKES_SUBPROBLEM_HH + +#include +#include + +#include +#include + +#include +#include +#include + +namespace Dumux { +template +class StokesSubProblem; + +namespace Properties { +// Create new type tags +namespace TTag { +struct StokesOneP { using InheritsFrom = std::tuple; }; +} // end namespace TTag + +// the fluid system +template +struct FluidSystem +{ + using Scalar = GetPropType; + using type = FluidSystems::OnePLiquid > ; +}; + +// Set the grid type +template +struct Grid { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates, 2> >; }; + +// Set the problem property +template +struct Problem { using type = Dumux::StokesSubProblem ; }; + +template +struct EnableGridGeometryCache { static constexpr bool value = true; }; +template +struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; +template +struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; +} // end namespace Properties + +/*! + * \ingroup BoundaryTests + * \brief The Stokes sub-problem of coupled Stokes-Darcy convergence test + */ +template +class StokesSubProblem : public NavierStokesProblem +{ + using ParentType = NavierStokesProblem; + using Scalar = GetPropType; + using BoundaryTypes = GetPropType; + using GridGeometry = GetPropType; + using GridView = typename GridGeometry::GridView; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using ModelTraits = GetPropType; + + using CouplingManager = GetPropType; + + using VelocityVector = typename Element::Geometry::GlobalCoordinate; + + + enum class TestCase + { + BJSymmetrized, NewICNonSymmetrized + }; + +public: + //! export the Indices + using Indices = typename ModelTraits::Indices; + + StokesSubProblem(std::shared_ptr gridGeometry, std::shared_ptr couplingManager) + : ParentType(gridGeometry, "Stokes"), couplingManager_(couplingManager) + { + problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + const auto testCaseInput = getParamFromGroup(this->paramGroup(), "Problem.TestCase", "BJSymmetrized"); + if (testCaseInput == "BJSymmetrized") + testCase_ = TestCase::BJSymmetrized; + else if (testCaseInput == "NewICNonSymmetrized") + testCase_ = TestCase::NewICNonSymmetrized; + else + DUNE_THROW(Dune::InvalidStateException, testCaseInput + " is not a valid test case"); + } + + /*! + * \brief The problem name. + */ + const std::string& name() const + { + return problemName_; + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the temperature within the domain in [K]. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + + /*! + * \brief Returns the sources within the domain. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition &globalPos) const + { + switch (testCase_) + { + case TestCase::BJSymmetrized: + return rhsBJSymmetrized_(globalPos); + case TestCase::NewICNonSymmetrized: + return rhsNewICNonSymmetrized_(globalPos); + default: + DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); + } + } + + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment. + * + * \param element The finite element + * \param scvf The sub control volume face + */ + BoundaryTypes boundaryTypes(const Element& element, + const SubControlVolumeFace& scvf) const + { + BoundaryTypes values; + + const auto& globalPos = scvf.dofPosition(); + + if(onUpperBoundary_(globalPos)) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + + // left/right wall + if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos))) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + + if (couplingManager().isCoupledEntity(CouplingManager::freeFlowIdx, scvf)) + { + values.setCouplingNeumann(Indices::conti0EqIdx); + values.setCouplingNeumann(Indices::momentumYBalanceIdx); + values.setSlipCondition(Indices::momentumXBalanceIdx); + } + + return values; + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet control volume. + */ + PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const + { + return analyticalSolution(globalPos); + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeometry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFaceVars The element face variables + * \param scvf The boundary sub control volume face + */ + template + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if(couplingManager().isCoupledEntity(CouplingManager::freeFlowIdx, scvf)) + { + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); + values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); + } + return values; + } + + // \} + + //! Get the coupling manager + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + /*! + * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann). + */ + VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + return couplingManager().couplingData().porousMediumVelocity(element, scvf); + } + + + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial value for a control volume. + * + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const + { + return PrimaryVariables(0.0); + } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter + for the Beavers-Joseph(-Saffman) boundary condition + */ + Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const + { + return couplingManager().couplingData().darcyPermeability(element, scvf); + } + + /*! + * \brief Returns the alpha value required as input parameter for the + Beavers-Joseph-Saffman boundary condition. + */ + Scalar alphaBJ(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); + } + + + //for new IC, input parameters that stay in the tangential coupling condition + /*! + * \brief Returns the scale separation parameter epsilon required as input parameter for the + new coupling conditions + */ + Scalar epsInterface(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().epsInterfaceAtPos(scvf.center()); + } + + /*! + * \brief Returns the boundary layer constant N_1_bl required as input parameter for the + new coupling condition for the tangential component + */ + Scalar factorNTangential(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::porousMediumIdx).spatialParams().factorNTangentialAtPos(scvf.center()); + } + + + /*! + * \brief Returns the analytical solution of the problem at a given position. + * + * \param globalPos The global position + */ + PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const + { + switch (testCase_) + { + case TestCase::BJSymmetrized: + return analyticalSolutionBJSymmetrized_(globalPos); + case TestCase::NewICNonSymmetrized: + return analyticalSolutionNewICNonSymmetrized_(globalPos); + default: + DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); + } + } + + // \} + +private: + + // see exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler) + PrimaryVariables analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const + { + PrimaryVariables sol(0.0); + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + + using std::sin; using std::cos; + sol[Indices::velocityXIdx] = (y-1.0)*(y-1.0) + x*(y-1.0) - x - 7.0 + 2.0*x*x + 2.0*y*y; + sol[Indices::velocityYIdx] = (x-1.0)*x - 0.5*(y-1.0)*(y-1.0) - 3.0*y - 3.0 -4.0*y*(x-1.0); + sol[Indices::pressureIdx] = -8.0*x + 1.0*y + 7.0 + x*x; + return sol; + } + + // see exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler) + NumEqVector rhsBJSymmetrized_(const GlobalPosition& globalPos) const + { + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + NumEqVector source(0.0); + using std::sin; using std::cos; + source[Indices::momentumXBalanceIdx] = 2.0*x - 18.0; + source[Indices::momentumYBalanceIdx] = 0.0; + return source; + } + + // see exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler) + PrimaryVariables analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const + { + PrimaryVariables sol(0.0); + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + + using std::exp; using std::sin; using std::cos; + sol[Indices::velocityXIdx] = (y-1.0)*(y-1.0) + x*(y-1.0) + 3.0*x + 1.5; + sol[Indices::velocityYIdx] = 0.5*x*(x-1.0) - 0.5*(y-1.0)*(y-1.0) - 3.0*y + 2.0; + sol[Indices::pressureIdx] = 2.0*x + y - 4.0; + return sol; + } + + // see exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler) + NumEqVector rhsNewICNonSymmetrized_(const GlobalPosition& globalPos) const + { + const Scalar x = globalPos[0]; + const Scalar y = globalPos[1]; + using std::exp; using std::sin; using std::cos; + NumEqVector source(0.0); + source[Indices::momentumXBalanceIdx] = -2.0; + source[Indices::momentumYBalanceIdx] = 1.0; + return source; + } + + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; } + + static constexpr Scalar eps_ = 1e-7; + std::string problemName_; + std::shared_ptr couplingManager_; + TestCase testCase_; +}; +} // end namespace Dumux + +#endif // DUMUX_STOKES_SUBPROBLEM_HH