From 3f16ad727616a461a709a8480ca25b9f970f7ac2 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 19 May 2020 15:48:48 +0200 Subject: [PATCH 01/40] Added nMomentum-IC, porousMediumVel --- .../boundary/stokesdarcy/box/couplingdata.hh | 267 +++++++++++++++++- 1 file changed, 266 insertions(+), 1 deletion(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index b78243f07f..b7b32ac5c5 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -28,7 +28,12 @@ #include #include -#include +#include //TODO: Why is this the right coupling manager, not the one in this folder? + +//Needed for nMomentumCouplingCondition +#include +#include +#include namespace Dumux { /*! @@ -58,6 +63,12 @@ class StokesDarcyCouplingDataBoxBase : public StokesDarcyCouplingDataImplementat static constexpr auto stokesIdx = CouplingManager::stokesIdx; static constexpr auto darcyIdx = CouplingManager::darcyIdx; + // Needed for velocityPorousMedium + using VelocityVector = typename Element::Geometry::GlobalCoordinate; + // Needed for nMomentum + 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>; @@ -127,6 +138,260 @@ public: return momentumFlux; } + + //TODO: Review! + /*! + * \brief Returns the new interface condition momentum flux across the coupling boundary. + * + * For the new momentum coupling, the porous medium side and a difference to the usual momentum flux is calculated + * + */ + template + Scalar nMomentumCouplingCondition(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& stokesElemVolVars, + const ElementFaceVariables& stokesElemFaceVars, + const SubControlVolumeFace& scvf) const + { + Scalar momentumFlux(0.0); + const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); + //######## darcy contribution ################# + // integrate darcy pressure over each coupling segment and average + for (const auto& data : stokesContext) + { + 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(); + + // do second order integration as box provides linear functions + static constexpr int darcyDim = GridGeometry::GridView::dimension; + const auto& rule = Dune::QuadratureRules::rule(data.segmentGeometry.type(), 2); + 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); + + std::vector> shapeValues; + localBasis.evaluateFunction(ipElementLocal, shapeValues); + + Scalar pressure = 0.0; + for (const auto& scv : scvs(data.fvGeometry)) + pressure += elemVolVars[scv].pressure(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; + + momentumFlux += pressure*data.segmentGeometry.integrationElement(qp.position())*qp.weight(); + } + } + } + momentumFlux /= scvf.area(); + + // normalize pressure + if(getPropValue, Properties::NormalizePressure>()) + momentumFlux -= this->couplingManager().problem(stokesIdx).initial(scvf)[Indices::pressureIdx]; + + //######## New (n) stokes contribution ################# + const std::size_t numSubFaces = scvf.pairData().size(); //numer of adjacent sub faces? + + // Account for all sub faces. + for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) + { + const auto eIdx = scvf.insideScvIdx(); + const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); + + // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction. + // Create a boundaryTypes object (will be empty if not at a boundary). + Dune::Std::optional currentScvfBoundaryTypes; + if (scvf.boundary()) + currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); + + // Getting boundary type for lateral face + Dune::Std::optional> lateralFaceBoundaryTypes; + if (lateralScvf.boundary()) + { + lateralFaceBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, lateralScvf)); + } + + + // Getting velocity gradients + const Scalar velocityGrad_ji = StokesVelocityGradients::velocityGradJI( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + const Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + + // Calculating additional term for momentum flux + //TODO: inverted sign... + const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); + //TODO: viscosity? + //TODO: ij should be 0 for unsymm, is this fullfilled? yes, but just if nTangential/bj/pressure bc is used + // Averaging the gradients to get evaluation at the center + momentumFlux += 1.0/numSubFaces*Nsbl * (velocityGrad_ji + velocityGrad_ij); + } + momentumFlux *= scvf.directionSign(); + + return momentumFlux; + } + + // TODO: Review! + /*! + * \brief Returns the velocity vector at the interface of the porous medium according to darcys law + * + * For the tangential (bj(s) and nTangential) coupling, the tangential porous medium velocity needs to + * be evaluated. We use darcys law and perform an integral average over all coupling segments + * + */ + VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + VelocityVector velocity(0.0); // velocity darcy + VelocityVector gradP(0.0); // pressure gradient darcy + Scalar rho(0.0); // density darcy + + //Getting needed information from the darcy domain + const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); + + //TODO: //Gravity changes darcy's law for calculating the velocity: ...-rho*g, + or -? + static const bool enableGravity = getParamFromGroup(this->couplingManager().problem(darcyIdx).paramGroup(), "Problem.EnableGravity"); + + // Iteraton 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 + static constexpr int darcyDim = GridGeometry::GridView::dimension; + 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; + //initialize the shape values + //TODO: move definitions outside the loop? + std::vector> shapeValues; + localBasis.evaluateFunction(ipElementLocal, shapeValues); + //and derivate values + using JacobianType = Dune::FieldMatrix; + std::vector shapeDerivates; + localBasis.evaluateJacobian(ipElementLocal, shapeDerivates); + + //calc pressure gradient and rho at qp, every scv belongs to one node + for (const auto& scv : scvs(data.fvGeometry)){ + gradP.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivates[scv.indexInElement()][0]); + if (enableGravity){ + rho += elemVolVars[scv].density(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; + } + } + //account gravity + if (enableGravity){ + gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); + } + + //Add the integrated segment velocity to the sum + velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()/data.volVars.viscosity(data.darcyScvfIdx), mv(K,gradP)); + } + } + } + //The integration is performed to get the average of the darcy velocity over one stokes face + velocity /= scvf.area(); + return velocity; + } + + + // TODO: Review! + /*! + * \brief Returns the velocity vector with a different permeability tensor + * + * For the new interface condition by Elissa Eggenweiler, a pm-velocity with altered permeability tensor needs to be evaluated + * + */ + VelocityVector newPorousMediumInterfaceVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + VelocityVector velocity(0.0); // velocity darcy + VelocityVector gradP(0.0); // pressure gradient darcy + Scalar rho(0.0); // density darcy + + //Getting needed information from the darcy domain + const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); + + //TODO: //Gravity changes darcy's law for calculating the velocity: ...-rho*g, + or -? + static const bool enableGravity = getParamFromGroup(this->couplingManager().problem(darcyIdx).paramGroup(), "Problem.EnableGravity"); + + // Iteraton 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(); + + // INTEGRATION, second order as box provides linear functions + static constexpr int darcyDim = GridGeometry::GridView::dimension; + 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; + //Darcy parameters + const auto& epsInterface = this->couplingManager().problem(darcyIdx).spatialParams().epsInterfaceAtPos(ipGlobal); + const auto& M = this->couplingManager().problem(darcyIdx).spatialParams().matrixNTangentialAtPos(ipGlobal); + + //initialize the shape values + //TODO: move definitions outside the loop? + std::vector> shapeValues; + localBasis.evaluateFunction(ipElementLocal, shapeValues); + //and derivate values + using JacobianType = Dune::FieldMatrix; + std::vector shapeDerivates; + localBasis.evaluateJacobian(ipElementLocal, shapeDerivates); + + //calc pressure gradient and rho at qp, every scv belongs to one node + for (const auto& scv : scvs(data.fvGeometry)){ + gradP.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivates[scv.indexInElement()][0]); + if (enableGravity){ + rho += elemVolVars[scv].density(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; + } + } + //account gravity + if (enableGravity){ + gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); + } + + //Add the integrated segment velocity to the sum + velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()*epsInterface*epsInterface/data.volVars.viscosity(data.darcyScvfIdx), mv(M,gradP)); + } + } + } + //The integration is performed to get the average of the darcy velocity over one stokes face + velocity /= scvf.area(); + return velocity; + } }; /*! -- GitLab From 54b3fc2875c9d465c7a136365793823215dd6c37 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 19 May 2020 15:51:27 +0200 Subject: [PATCH 02/40] Added nTangential(-IC) boundary type --- .../staggered/freeflow/boundarytypes.hh | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index a9a1f2d4e0..c5ff3031e5 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -55,6 +55,7 @@ public: boundaryInfo_[eqIdx].visited = false; boundaryInfo_[eqIdx].isSymmetry = false; boundaryInfo_[eqIdx].isBeaversJoseph = false; + boundaryInfo_[eqIdx].isNTangential = false; } /*! @@ -104,6 +105,16 @@ public: boundaryInfo_[eqIdx].visited = true; boundaryInfo_[eqIdx].isBeaversJoseph = true; } + /*! + * \brief Set a boundary condition for a single equation to + * new slip condition by Elissa Eggenweiler + */ + void setNTangential(unsigned eqIdx) + { + resetEq(eqIdx); + boundaryInfo_[eqIdx].visited = true; + boundaryInfo_[eqIdx].isNTangential = true; + } /*! * \brief Returns true if an equation is used to specify a @@ -114,6 +125,15 @@ public: bool isBeaversJoseph(unsigned eqIdx) const { return boundaryInfo_[eqIdx].isBeaversJoseph; } + /*! + * \brief Returns true if an equation is used to specify a + * nTangential boundary condition. + * + * \param eqIdx The index of the equation + */ + bool isNTangential(unsigned eqIdx) const + { return boundaryInfo_[eqIdx].isNTangential; } + /*! * \brief Returns true if some equation is used to specify a * Beavers-Joseph(-Saffman) boundary condition. @@ -126,12 +146,25 @@ public: return false; } + /*! + * \brief Returns true if some equation is used to specify a + * nTangential boundary condition. + */ + bool hasNTangential() const + { + for (int i = 0; i < numEq; ++i) + if (boundaryInfo_[i].isNTangential) + return true; + return false; + } + protected: struct StaggeredFreeFlowBoundaryInfo { bool visited; bool isSymmetry; bool isBeaversJoseph; + bool isNTangential; }; std::array boundaryInfo_; -- GitLab From 8213a6d0b3148f878384093d17993bf54a41d0aa Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 19 May 2020 15:53:27 +0200 Subject: [PATCH 03/40] Duplicated code from bj for nTangential --- dumux/freeflow/navierstokes/problem.hh | 20 ++++ .../navierstokes/staggered/fluxvariables.hh | 17 +++- .../staggered/staggeredupwindfluxvariables.hh | 4 + .../staggered/velocitygradients.hh | 94 +++++++++++++++++++ 4 files changed, 131 insertions(+), 4 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index ab62ad440c..6684f32536 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -242,6 +242,26 @@ public: + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0); } + //TODO: I assume an inverted sign for factorNTangential... + //! helper function to evaluate the slip velocity on the boundary when the new tangential condition is used + const Scalar nTangentialVelocity(const Element& element, + const SubControlVolume& scv, + const SubControlVolumeFace& ownScvf, //stehendes + const SubControlVolumeFace& faceOnPorousBoundary, //liegendes + const Scalar velocitySelf, //vel auf stehendem + const Scalar tangentialVelocityGradient) const //dv/dx=0 + { + const Scalar factor = 1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); + 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_().newPorousMediumInterfaceVelocity(element, faceOnPorousBoundary) * orientation * factor * distanceNormalToBoundary + + velocitySelf) / (factor*distanceNormalToBoundary + 1.0); + } + private: //! Returns the implementation of the problem (i.e. static polymorphism) Implementation &asImp_() diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index bbd5e19812..121d9104d2 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -331,8 +331,10 @@ 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 bj = currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())); + bool nTangential = currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex())); + bool neumann = lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())); + if ( (bj || nTangential) && lateralFaceBoundaryTypes && neumann) { const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, @@ -372,10 +374,11 @@ public: // Check the consistency of the boundary conditions, exactly one of the following must be set if (lateralFaceBoundaryTypes) { - std::bitset<3> admittableBcTypes; + std::bitset<4> 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(3, lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()))); if (admittableBcTypes.count() != 1) { DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf " @@ -507,6 +510,10 @@ private: return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } + else if (bcTypes.isNTangential(Indices::velocity(lateralFace.directionIndex()))){ + return VelocityGradients::nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + } else return faceVars.velocityLateralInside(localSubFaceIdx); } @@ -573,7 +580,9 @@ private: { if (!scvf.boundary() || currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) || - currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex()))) + currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())) || + currentScvfBoundaryTypes->isNTangential(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 bcae100b68..3241255d26 100644 --- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh @@ -584,6 +584,7 @@ private: const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)); const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())); + const bool lateralFaceHasNTangential = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex())); const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())); const Scalar velocitySelf = faceVars.velocitySelf(); @@ -596,6 +597,9 @@ private: return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + else if(lateralFaceHasNTangential) + return VelocityGradients::nTangentialVelocityAtLateralScvf(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 a9c3c452b3..51dfabd4a6 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -121,6 +121,10 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())]; } + else if (lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()))){ + return nTangentialVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + } else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))) { return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, @@ -200,6 +204,10 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())]; } + else if (currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex()))){ + return nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + } else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex()))) { return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, @@ -280,6 +288,49 @@ public: tangentialVelocityGradient); } + template + static Scalar nTangentialVelocityAtCurrentScvf(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf, + const FaceVariables& faceVars, + const Dune::Std::optional& currentScvfBoundaryTypes, + const Dune::Std::optional& lateralFaceBoundaryTypes, + const std::size_t localSubFaceIdx) + { + const auto eIdx = scvf.insideScvIdx(); + const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); + const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx); + + 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 + // (towards the current scvf). + static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + + if (unsymmetrizedGradientForBJ) + return 0.0; + + if (lateralScvf.boundary()) + { + if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) || + lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()))) + return 0.0; + } + + return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + }(); + + return problem.nTangentialVelocity(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. * @@ -340,6 +391,49 @@ public: tangentialVelocityGradient); } + template + static Scalar nTangentialVelocityAtLateralScvf(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const SubControlVolumeFace& scvf, + const FaceVariables& faceVars, + const Dune::Std::optional& currentScvfBoundaryTypes, + const Dune::Std::optional& lateralFaceBoundaryTypes, + const std::size_t localSubFaceIdx) + { + const auto eIdx = scvf.insideScvIdx(); + const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); + const Scalar innerParallelVelocity = faceVars.velocitySelf(); + + 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 + // (towards the current scvf). + static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + + if (unsymmetrizedGradientForBJ) + return 0.0; + + if (scvf.boundary()) + { + if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) || + currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex()))) + return 0.0; + } + + return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + }(); + + return problem.nTangentialVelocity(element, + fvGeometry.scv(scvf.insideScvIdx()), + scvf, + lateralScvf, /*on boundary*/ + innerParallelVelocity, + tangentialVelocityGradient); + } + private: /*! -- GitLab From 9fbdd106617104c4f264d14695e28cfeb59b1aae Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 19 May 2020 15:54:51 +0200 Subject: [PATCH 04/40] Added functions "linking" to spatialParameters for nTangential and nMomentum --- dumux/freeflow/navierstokes/problem.hh | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 6684f32536..5cfd08eba3 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -212,6 +212,18 @@ public: return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf)); } + Scalar epsInterface(const SubControlVolumeFace& scvf) const + { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem");} + + Scalar factorNMomentum(const SubControlVolumeFace& scvf) const + { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNMomentum value must be returned in the acutal problem");} + + Scalar factorNTangential(const SubControlVolumeFace& scvf) const + { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the acutal problem");} + + Dune::FieldMatrix matrixNTangential(const SubControlVolumeFace& scvf) const + { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the matrixNTangential value must be returned in the acutal problem");} + /*! * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann). */ -- GitLab From a3e0ddef0f8f26d3f701c5cd799c2b2cb2a13b59 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 19 May 2020 15:59:55 +0200 Subject: [PATCH 05/40] Call actual implementation for porousMediumVelocity --- dumux/freeflow/navierstokes/problem.hh | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 5cfd08eba3..5503be390c 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -229,7 +229,15 @@ public: */ VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - return VelocityVector(0.0); + return asImp_().porousMediumVelocity(element, scvf); + } + + /*! + * \brief Returns the darcy velocity vector with a different permeability tensor + */ + VelocityVector newPorousMediumInterfaceVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + return asImp_().newPorousMediumInterfaceVelocity(element, scvf); } //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used -- GitLab From 95a8de4f468e92b96d59472119e7652ab11bee4f Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Wed, 20 May 2020 12:00:20 +0200 Subject: [PATCH 06/40] replaced Dune::Std::optional with std::optonal --- .../freeflow/navierstokes/staggered/velocitygradients.hh | 8 ++++---- .../multidomain/boundary/stokesdarcy/box/couplingdata.hh | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh index 51dfabd4a6..3fa41ae2f3 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -294,8 +294,8 @@ public: const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const FaceVariables& faceVars, - const Dune::Std::optional& currentScvfBoundaryTypes, - const Dune::Std::optional& lateralFaceBoundaryTypes, + const std::optional& currentScvfBoundaryTypes, + const std::optional& lateralFaceBoundaryTypes, const std::size_t localSubFaceIdx) { const auto eIdx = scvf.insideScvIdx(); @@ -397,8 +397,8 @@ public: const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const FaceVariables& faceVars, - const Dune::Std::optional& currentScvfBoundaryTypes, - const Dune::Std::optional& lateralFaceBoundaryTypes, + const std::optional& currentScvfBoundaryTypes, + const std::optional& lateralFaceBoundaryTypes, const std::size_t localSubFaceIdx) { const auto eIdx = scvf.insideScvIdx(); diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index b7b32ac5c5..9b2f287222 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -33,7 +33,7 @@ //Needed for nMomentumCouplingCondition #include #include -#include +#include namespace Dumux { /*! @@ -203,12 +203,12 @@ public: // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction. // Create a boundaryTypes object (will be empty if not at a boundary). - Dune::Std::optional currentScvfBoundaryTypes; + std::optional> currentScvfBoundaryTypes; if (scvf.boundary()) currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); // Getting boundary type for lateral face - Dune::Std::optional> lateralFaceBoundaryTypes; + std::optional> lateralFaceBoundaryTypes; if (lateralScvf.boundary()) { lateralFaceBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, lateralScvf)); -- GitLab From 38ce7af77fdb658edb2e5b4aaac2cab37d992cdb Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Wed, 20 May 2020 17:51:01 +0200 Subject: [PATCH 07/40] removed asImp_().... --- dumux/freeflow/navierstokes/problem.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 5503be390c..6f64fdca1e 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -229,7 +229,7 @@ public: */ VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - return asImp_().porousMediumVelocity(element, scvf); + return VelocityVector(0.0);//TODO: -> Dont force implementation? } /*! @@ -237,8 +237,8 @@ public: */ VelocityVector newPorousMediumInterfaceVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - return asImp_().newPorousMediumInterfaceVelocity(element, scvf); - } + DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the newPorousMediumInterfaceVelocity must be returned in the acutal problem"); + } //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used const Scalar beaversJosephVelocity(const Element& element, -- GitLab From e2353dc1208f2f54a29344248cef91acf7649fe2 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Wed, 20 May 2020 13:13:47 +0200 Subject: [PATCH 08/40] [md][ffpm][box] Update indices and correct include --- .../stokesdarcy/box/couplingmanager.hh | 166 +++++++++--------- 1 file changed, 82 insertions(+), 84 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingmanager.hh index 4bf72d67c6..56181d1153 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()); } -- GitLab From e259d213c688fd6bd96b4aa80dbb98e330347cdc Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Sat, 23 May 2020 22:58:03 +0200 Subject: [PATCH 09/40] "TODO" -> "TODO by Lars" ==> searchable --- .../boundary/stokesdarcy/box/couplingdata.hh | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 9b2f287222..decad5cfb2 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -28,7 +28,7 @@ #include #include -#include //TODO: Why is this the right coupling manager, not the one in this folder? +#include //TODO by Lars: Why is this the right coupling manager, not the one in this folder? //Needed for nMomentumCouplingCondition #include @@ -139,7 +139,8 @@ public: return momentumFlux; } - //TODO: Review! + //TODO by Lars: Review! + //TODO by Lars: Use momentumCouplingCondition(...)/introduce new method to remove dulplicated code? /*! * \brief Returns the new interface condition momentum flux across the coupling boundary. * @@ -224,10 +225,10 @@ public: currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); // Calculating additional term for momentum flux - //TODO: inverted sign... + //TODO by Lars: inverted sign... const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); - //TODO: viscosity? - //TODO: ij should be 0 for unsymm, is this fullfilled? yes, but just if nTangential/bj/pressure bc is used + //TODO by Lars: viscosity? + //TODO by Lars: ij should be 0 for unsymm, is this fullfilled? yes, but just if nTangential/bj/pressure bc is used // Averaging the gradients to get evaluation at the center momentumFlux += 1.0/numSubFaces*Nsbl * (velocityGrad_ji + velocityGrad_ij); } @@ -236,7 +237,7 @@ public: return momentumFlux; } - // TODO: Review! + // TODO by Lars: Review! /*! * \brief Returns the velocity vector at the interface of the porous medium according to darcys law * @@ -252,8 +253,6 @@ public: //Getting needed information from the darcy domain const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); - - //TODO: //Gravity changes darcy's law for calculating the velocity: ...-rho*g, + or -? static const bool enableGravity = getParamFromGroup(this->couplingManager().problem(darcyIdx).paramGroup(), "Problem.EnableGravity"); // Iteraton over the different coupling segments @@ -300,6 +299,7 @@ public: } } //account gravity + //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g correct? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } @@ -315,7 +315,7 @@ public: } - // TODO: Review! + // TODO by Lars: Review! /*! * \brief Returns the velocity vector with a different permeability tensor * @@ -331,7 +331,6 @@ public: //Getting needed information from the darcy domain const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); - //TODO: //Gravity changes darcy's law for calculating the velocity: ...-rho*g, + or -? static const bool enableGravity = getParamFromGroup(this->couplingManager().problem(darcyIdx).paramGroup(), "Problem.EnableGravity"); // Iteraton over the different coupling segments @@ -363,7 +362,7 @@ public: const auto& M = this->couplingManager().problem(darcyIdx).spatialParams().matrixNTangentialAtPos(ipGlobal); //initialize the shape values - //TODO: move definitions outside the loop? + //TODO by Lars: move definitions outside the loop? std::vector> shapeValues; localBasis.evaluateFunction(ipElementLocal, shapeValues); //and derivate values @@ -379,6 +378,7 @@ public: } } //account gravity + //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g correct? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } -- GitLab From 419cfef4351ad44231c1f0bcbcc1a9fb9ccfadde Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Sun, 5 Jul 2020 22:12:16 +0200 Subject: [PATCH 10/40] fixing call to viscosity --- dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index decad5cfb2..8a07091851 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -305,7 +305,7 @@ public: } //Add the integrated segment velocity to the sum - velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()/data.volVars.viscosity(data.darcyScvfIdx), mv(K,gradP)); + velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()/data.volVars.viscosity(darcyPhaseIdx), mv(K,gradP)); } } } @@ -384,7 +384,7 @@ public: } //Add the integrated segment velocity to the sum - velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()*epsInterface*epsInterface/data.volVars.viscosity(data.darcyScvfIdx), mv(M,gradP)); + velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()*epsInterface*epsInterface/data.volVars.viscosity(darcyPhaseIdx), mv(M,gradP)); } } } -- GitLab From c2b875f4e25efc9df6673a044aa3a081203d6996 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Thu, 9 Jul 2020 08:53:02 +0200 Subject: [PATCH 11/40] fixing Dimension of CouplingSegment::Geometry --- dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingmapper.hh index 68af4896cd..3cf3e9b51c 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}); -- GitLab From 2ba648c20dea87a76e46dc59c020db5ff4f4de67 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Thu, 9 Jul 2020 08:56:19 +0200 Subject: [PATCH 12/40] Finished function porousMediumVelocity --- .../boundary/stokesdarcy/box/couplingdata.hh | 56 ++++++++++--------- 1 file changed, 31 insertions(+), 25 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 8a07091851..e9affc009c 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -237,25 +237,30 @@ public: return momentumFlux; } - // TODO by Lars: Review! /*! - * \brief Returns the velocity vector at the interface of the porous medium according to darcys law + * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law * - * For the tangential (bj(s) and nTangential) coupling, the tangential porous medium velocity needs to - * be evaluated. We use darcys law and perform an integral average over all coupling segments + * The tangential porous medium velocity needs to be evaluated at the stokes domain for the tangential (bj and nTangential) + * coupling at the stokes-darcy interface. We use darcys law and perform an integral average over all coupling segments. * */ VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { + 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"); - // Iteraton over the different coupling segments + // Iteration over the different coupling segments for (const auto& data : stokesContext) { //We are on (one of) the correct scvf(s) @@ -266,11 +271,11 @@ public: 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 - static constexpr int darcyDim = GridGeometry::GridView::dimension; const auto& rule = Dune::QuadratureRules::rule(data.segmentGeometry.type(), 2); //Loop over all quadrature points in the rule for (const auto& qp : rule) @@ -282,35 +287,36 @@ public: //reset pressure gradient and rho at this qp gradP=0.0; rho=0.0; - //initialize the shape values - //TODO: move definitions outside the loop? - std::vector> shapeValues; + //TODO: Is this needed? + shapeValues.clear(); + shapeDerivatives.clear(); + + //calculate the shape and derivative values at the qp localBasis.evaluateFunction(ipElementLocal, shapeValues); - //and derivate values - using JacobianType = Dune::FieldMatrix; - std::vector shapeDerivates; - localBasis.evaluateJacobian(ipElementLocal, shapeDerivates); + 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.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivates[scv.indexInElement()][0]); + //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 gravity - //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g correct? + //account for gravity + //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g the correct sign? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } - - //Add the integrated segment velocity to the sum - velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()/data.volVars.viscosity(darcyPhaseIdx), mv(K,gradP)); + //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*K/mu*gradP + //TODO: which fits dumux style better? + K.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), gradP, velocity); + //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), mv(K,gradP)); } + intersectionLength += data.segmentGeometry.volume(); } } - //The integration is performed to get the average of the darcy velocity over one stokes face - velocity /= scvf.area(); + velocity /= intersectionLength; //averaging return velocity; } @@ -365,14 +371,14 @@ public: //TODO by Lars: move definitions outside the loop? std::vector> shapeValues; localBasis.evaluateFunction(ipElementLocal, shapeValues); - //and derivate values + //and derivative values using JacobianType = Dune::FieldMatrix; - std::vector shapeDerivates; - localBasis.evaluateJacobian(ipElementLocal, shapeDerivates); + std::vector shapeDerivatives; + 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.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivates[scv.indexInElement()][0]); + gradP.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivatives[scv.indexInElement()][0]); if (enableGravity){ rho += elemVolVars[scv].density(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; } -- GitLab From 4b0ba300356f5c6cb01076dfed4d227367573b1a Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 21 Jul 2020 10:38:35 +0200 Subject: [PATCH 13/40] newPorousMediumInterfaceVelocity finished --- .../boundary/stokesdarcy/box/couplingdata.hh | 54 +++++++++++-------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index e9affc009c..e9be87bee3 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -230,10 +230,11 @@ public: //TODO by Lars: viscosity? //TODO by Lars: ij should be 0 for unsymm, is this fullfilled? yes, but just if nTangential/bj/pressure bc is used // Averaging the gradients to get evaluation at the center + std::cout<<"velGrad ij, ji: "<< velocityGrad_ij << "; " << velocityGrad_ji <& element, const SubControlVolumeFace& scvf) const { + 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"); - // Iteraton over the different coupling segments + // Iteration over the different coupling segments for (const auto& data : stokesContext) { //We are on (one of) the correct scvf(s) @@ -350,8 +356,11 @@ public: 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 - static constexpr int darcyDim = GridGeometry::GridView::dimension; const auto& rule = Dune::QuadratureRules::rule(data.segmentGeometry.type(), 2); //Loop over all quadrature points in the rule for (const auto& qp : rule) @@ -363,39 +372,40 @@ public: //reset pressure gradient and rho at this qp gradP=0.0; rho=0.0; - //Darcy parameters + //TODO: Is this needed? + shapeValues.clear(); + shapeDerivatives.clear(); + + // darcy spatial dependent parameters const auto& epsInterface = this->couplingManager().problem(darcyIdx).spatialParams().epsInterfaceAtPos(ipGlobal); const auto& M = this->couplingManager().problem(darcyIdx).spatialParams().matrixNTangentialAtPos(ipGlobal); - //initialize the shape values - //TODO by Lars: move definitions outside the loop? - std::vector> shapeValues; + //calculate the shape and derivative values at the qp localBasis.evaluateFunction(ipElementLocal, shapeValues); - //and derivative values - using JacobianType = Dune::FieldMatrix; - std::vector shapeDerivatives; 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.axpy(elemVolVars[scv].pressure(darcyPhaseIdx),shapeDerivatives[scv.indexInElement()][0]); + //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 gravity - //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g correct? + //account for gravity + //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g the correct sign? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } - - //Add the integrated segment velocity to the sum - velocity.axpy(-data.segmentGeometry.integrationElement(qp.position())*qp.weight()*epsInterface*epsInterface/data.volVars.viscosity(darcyPhaseIdx), mv(M,gradP)); + //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*eps**2*M/mu*gradP + //TODO: which fits dumux style better? + M.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); + //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, mv(M,gradP)); } + intersectionLength += data.segmentGeometry.volume(); } } - //The integration is performed to get the average of the darcy velocity over one stokes face - velocity /= scvf.area(); + velocity /= intersectionLength; //averaging return velocity; } }; -- GitLab From 0252539ad9975a8ca5b9994a1eee12a5ba2f9a44 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 21 Jul 2020 10:38:54 +0200 Subject: [PATCH 14/40] Debugging activated --- cmake.opts | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake.opts b/cmake.opts index 8c5f184a1b..5dedd1d416 100644 --- a/cmake.opts +++ b/cmake.opts @@ -54,6 +54,6 @@ CMAKE_FLAGS="$SPECIFIC_COMPILER $SPECIFIC_GENERATOR $OPM_FLAGS -DCMAKE_CXX_FLAGS_RELEASE='$GXX_RELEASE_OPTS $GXX_RELEASE_WARNING_OPTS $DUMUX_DISABLE_PROP_MACROS' -DCMAKE_CXX_FLAGS_DEBUG='-O0 -g -ggdb -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare $DUMUX_DISABLE_PROP_MACROS' -DCMAKE_CXX_FLAGS_RELWITHDEBINFO='$GXX_RELEASE_OPTS $GXX_RELEASE_WARNING_OPTS -g -ggdb -Wall $DUMUX_DISABLE_PROP_MACROS' --DCMAKE_BUILD_TYPE=Release +-DCMAKE_BUILD_TYPE=Debug -DENABLE_HEADERCHECK=$DUMUX_ENABLE_HEADERCHECK " -- GitLab From cebe7d8c632d8b52ac613315fc22c09eb0f42a49 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 21 Jul 2020 10:43:18 +0200 Subject: [PATCH 15/40] removing newPorousMediumInterfaceVelocity in problem.hh --- dumux/freeflow/navierstokes/problem.hh | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 6f64fdca1e..a3d13a3747 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -232,14 +232,6 @@ public: return VelocityVector(0.0);//TODO: -> Dont force implementation? } - /*! - * \brief Returns the darcy velocity vector with a different permeability tensor - */ - VelocityVector newPorousMediumInterfaceVelocity(const Element& element, const SubControlVolumeFace& scvf) const - { - DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the newPorousMediumInterfaceVelocity must be returned in the acutal problem"); - } - //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used const Scalar beaversJosephVelocity(const Element& element, const SubControlVolume& scv, @@ -278,7 +270,7 @@ public: GlobalPosition orientation = ownScvf.unitOuterNormal(); orientation[ownScvf.directionIndex()] = 1.0; return (tangentialVelocityGradient*distanceNormalToBoundary - + asImp_().newPorousMediumInterfaceVelocity(element, faceOnPorousBoundary) * orientation * factor * distanceNormalToBoundary + + asImp_().porousMediumVelocity(element, faceOnPorousBoundary) * orientation * factor * distanceNormalToBoundary + velocitySelf) / (factor*distanceNormalToBoundary + 1.0); } -- GitLab From 8000d2c68f573a69deff1b029ba5bd948825cc08 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Wed, 22 Jul 2020 12:33:37 +0200 Subject: [PATCH 16/40] changed signs --- dumux/freeflow/navierstokes/problem.hh | 5 +++-- .../boundary/stokesdarcy/box/couplingdata.hh | 16 ++++++++++------ 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index a3d13a3747..6435fdbe0b 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -254,7 +254,7 @@ public: + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0); } - //TODO: I assume an inverted sign for factorNTangential... + //TODO: viscosity? //! helper function to evaluate the slip velocity on the boundary when the new tangential condition is used const Scalar nTangentialVelocity(const Element& element, const SubControlVolume& scv, @@ -263,7 +263,8 @@ public: const Scalar velocitySelf, //vel auf stehendem const Scalar tangentialVelocityGradient) const //dv/dx=0 { - const Scalar factor = 1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); + //const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].effectiveViscosity(); + const Scalar factor = -1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); // create a unit normal vector oriented in positive coordinate direction diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index e9be87bee3..1d8f67405f 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -220,18 +220,22 @@ public: const Scalar velocityGrad_ji = StokesVelocityGradients::velocityGradJI( this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - const Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( + Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + const bool unsymmetrizedGradientForBJ = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + //TODO: Remove calculation above in this case + if (unsymmetrizedGradientForBJ){ + velocityGrad_ij = 0.0; + } + // Calculating additional term for momentum flux - //TODO by Lars: inverted sign... const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); - //TODO by Lars: viscosity? - //TODO by Lars: ij should be 0 for unsymm, is this fullfilled? yes, but just if nTangential/bj/pressure bc is used // Averaging the gradients to get evaluation at the center - std::cout<<"velGrad ij, ji: "<< velocityGrad_ij << "; " << velocityGrad_ji < Date: Wed, 22 Jul 2020 12:34:01 +0200 Subject: [PATCH 17/40] changed contribution of viscosity --- dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 1d8f67405f..57af2a1b00 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -403,7 +403,7 @@ public: } //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*eps**2*M/mu*gradP //TODO: which fits dumux style better? - M.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); + M.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)*data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, mv(M,gradP)); } intersectionLength += data.segmentGeometry.volume(); -- GitLab From bdccbf53f2d8ab99b0db937281ac9a78d27e1432 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 24 Jul 2020 11:17:11 +0200 Subject: [PATCH 18/40] replaced effectiveViscosity by (dynamic)Viscosity --- dumux/freeflow/navierstokes/problem.hh | 2 +- dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 6435fdbe0b..44303a1e87 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -263,7 +263,7 @@ public: const Scalar velocitySelf, //vel auf stehendem const Scalar tangentialVelocityGradient) const //dv/dx=0 { - //const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].effectiveViscosity(); + //const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].Viscosity(); const Scalar factor = -1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 57af2a1b00..1f2307ec72 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -234,7 +234,7 @@ public: // Calculating additional term for momentum flux const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); // Averaging the gradients to get evaluation at the center - const Scalar viscosity = stokesElemVolVars[scvf.insideScvIdx()].effectiveViscosity(); + const Scalar viscosity = stokesElemVolVars[scvf.insideScvIdx()].viscosity(); momentumFlux -= 1.0/numSubFaces * viscosity * Nsbl * (velocityGrad_ji + velocityGrad_ij); } momentumFlux *= scvf.directionSign(); -- GitLab From 7486751e3d98d08621002c28219d20b8fa989d3e Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 24 Jul 2020 14:17:31 +0200 Subject: [PATCH 19/40] Removed unnecessary functions --- dumux/freeflow/navierstokes/problem.hh | 5 ----- 1 file changed, 5 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 44303a1e87..342b5a3292 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -215,14 +215,9 @@ public: Scalar epsInterface(const SubControlVolumeFace& scvf) const { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem");} - Scalar factorNMomentum(const SubControlVolumeFace& scvf) const - { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNMomentum value must be returned in the acutal problem");} - Scalar factorNTangential(const SubControlVolumeFace& scvf) const { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the acutal problem");} - Dune::FieldMatrix matrixNTangential(const SubControlVolumeFace& scvf) const - { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the matrixNTangential value must be returned in the acutal problem");} /*! * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann). -- GitLab From cb8ddf0eb87808c23b5c1b6e603dd4b4fbb3f7ec Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 24 Jul 2020 14:18:13 +0200 Subject: [PATCH 20/40] Comments changed --- dumux/freeflow/navierstokes/problem.hh | 6 +++- .../boundary/stokesdarcy/box/couplingdata.hh | 32 +++++++++---------- 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 342b5a3292..c118209903 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -212,6 +212,10 @@ public: return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf)); } + //TODO: Elissa, maybe you can write a small description, see above for these functions + /*! + * \brief Returns the ... value, + */ Scalar epsInterface(const SubControlVolumeFace& scvf) const { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem");} @@ -224,7 +228,7 @@ public: */ VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - return VelocityVector(0.0);//TODO: -> Dont force implementation? + return VelocityVector(0.0); } //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 1f2307ec72..878c1640ff 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -142,9 +142,10 @@ public: //TODO by Lars: Review! //TODO by Lars: Use momentumCouplingCondition(...)/introduce new method to remove dulplicated code? /*! - * \brief Returns the new interface condition momentum flux across the coupling boundary. + * \brief Returns the momentum flux across the coupling boundary which is calculated according to the new interface condition * - * For the new momentum coupling, the porous medium side and a difference to the usual momentum flux is calculated + * For the new normal momentum coupling, the porous medium side and also the stokes side is evaluated. + * [p]^pm + N_s^{bl} \tau T n * */ template @@ -193,22 +194,20 @@ public: if(getPropValue, Properties::NormalizePressure>()) momentumFlux -= this->couplingManager().problem(stokesIdx).initial(scvf)[Indices::pressureIdx]; - //######## New (n) stokes contribution ################# - const std::size_t numSubFaces = scvf.pairData().size(); //numer of adjacent sub faces? + //######## New stokes contribution ################# + const std::size_t numSubFaces = scvf.pairData().size(); - // Account for all sub faces. + // Account for all sub faces for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) { const auto eIdx = scvf.insideScvIdx(); const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction. - // Create a boundaryTypes object (will be empty if not at a boundary). + // Create a boundaryTypes object (will be empty if not at a boundary) std::optional> currentScvfBoundaryTypes; if (scvf.boundary()) currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); - // Getting boundary type for lateral face std::optional> lateralFaceBoundaryTypes; if (lateralScvf.boundary()) { @@ -216,7 +215,7 @@ public: } - // Getting velocity gradients + // Get velocity gradients const Scalar velocityGrad_ji = StokesVelocityGradients::velocityGradJI( this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); @@ -231,22 +230,21 @@ public: velocityGrad_ij = 0.0; } - // Calculating additional term for momentum flux + // Calculate stokes contribution to momentum flux: N_s^{bl} \tau T n const Scalar Nsbl = this->couplingManager().problem(darcyIdx).spatialParams().factorNMomentumAtPos(scvf.center()); - // Averaging the gradients to get evaluation at the 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(); - std::cout<<"momFlux" << momentumFlux<& element, const SubControlVolumeFace& scvf) const @@ -259,7 +257,7 @@ public: 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 + 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); @@ -313,7 +311,7 @@ public: if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } - //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*K/mu*gradP + //Add the integrated segment velocity to the sum: v+= -weight_k * sqrt(det(A^T*A))*K/mu*gradP //TODO: which fits dumux style better? K.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), gradP, velocity); //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), mv(K,gradP)); @@ -327,7 +325,7 @@ public: /*! - * \brief Returns the averaged velocity vector at the interface of the porous medium with a different permeability tensor according to darcys law + * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law with a different permeability tensor * * For the new tangential interface condition by Elissa Eggenweiler, a porous medium velocity with altered permeability tensor needs to be evaluated. * We use darcys law and perform an integral average over all coupling segments. -- GitLab From 5da23b771bff01537ead458abdde7029ec9044e1 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 24 Jul 2020 14:53:48 +0200 Subject: [PATCH 21/40] changed v_pm_int sign to fit paper --- dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 878c1640ff..f40bf96f36 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -401,7 +401,7 @@ public: } //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*eps**2*M/mu*gradP //TODO: which fits dumux style better? - M.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)*data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); + M.usmv(qp.weight()*data.segmentGeometry.integrationElement(ipLocal)*data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, mv(M,gradP)); } intersectionLength += data.segmentGeometry.volume(); -- GitLab From 487a88da8dd8d01045a62d5b4e4a355403cf3d0b Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 24 Jul 2020 16:26:42 +0200 Subject: [PATCH 22/40] Combining bj and nTangential to slip --- .../staggered/freeflow/boundarytypes.hh | 39 ++--- dumux/freeflow/navierstokes/problem.hh | 16 ++ .../navierstokes/staggered/fluxvariables.hh | 25 +-- .../staggered/staggeredupwindfluxvariables.hh | 11 +- .../staggered/velocitygradients.hh | 160 ++++-------------- dumux/freeflow/rans/problem.hh | 28 +-- .../boundary/stokesdarcy/box/couplingdata.hh | 49 +++++- 7 files changed, 136 insertions(+), 192 deletions(-) diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index c5ff3031e5..b55640dc5c 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -54,8 +54,7 @@ public: boundaryInfo_[eqIdx].visited = false; boundaryInfo_[eqIdx].isSymmetry = false; - boundaryInfo_[eqIdx].isBeaversJoseph = false; - boundaryInfo_[eqIdx].isNTangential = false; + boundaryInfo_[eqIdx].isSlipCondition = false; } /*! @@ -95,27 +94,28 @@ 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.). */ void setBeaversJoseph(unsigned eqIdx) { - resetEq(eqIdx); - boundaryInfo_[eqIdx].visited = true; - boundaryInfo_[eqIdx].isBeaversJoseph = true; + setSlipCondition(eqIdx); } /*! * \brief Set a boundary condition for a single equation to - * new slip condition by Elissa Eggenweiler + * slip condition: bjs,bj or the new one by Elissa Eggenweiler */ - void setNTangential(unsigned eqIdx) + void setSlipCondition(unsigned eqIdx) { resetEq(eqIdx); boundaryInfo_[eqIdx].visited = true; - boundaryInfo_[eqIdx].isNTangential = 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. @@ -123,37 +123,35 @@ public: * \param eqIdx The index of the equation */ bool isBeaversJoseph(unsigned eqIdx) const - { return boundaryInfo_[eqIdx].isBeaversJoseph; } + { return isSlipCondition(eqIdx); } /*! * \brief Returns true if an equation is used to specify a - * nTangential boundary condition. + * slip condition * * \param eqIdx The index of the equation */ - bool isNTangential(unsigned eqIdx) const - { return boundaryInfo_[eqIdx].isNTangential; } + 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. */ bool hasBeaversJoseph() const { - for (int i = 0; i < numEq; ++i) - if (boundaryInfo_[i].isBeaversJoseph) - return true; - return false; + return hasSlipCondition(); } /*! * \brief Returns true if some equation is used to specify a - * nTangential boundary condition. + * slip condition */ - bool hasNTangential() const + bool hasSlipCondition() const { for (int i = 0; i < numEq; ++i) - if (boundaryInfo_[i].isNTangential) + if (boundaryInfo_[i].isSlipCondition) return true; return false; } @@ -163,8 +161,7 @@ protected: { bool visited; bool isSymmetry; - bool isBeaversJoseph; - bool isNTangential; + bool isSlipCondition; }; std::array boundaryInfo_; diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index c118209903..7424b9dfb0 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -231,6 +231,22 @@ public: return VelocityVector(0.0); } + //! helper function to evaluate the slip velocity on the boundary + const Scalar slipVelocity(const Element& element, + const SubControlVolume& scv, + const SubControlVolumeFace& ownScvf, + const SubControlVolumeFace& faceOnPorousBoundary, + const Scalar velocitySelf, + const Scalar tangentialVelocityGradient) const + { + if (getParamFromGroup("Problem", "NewIc", false)){ + return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + } + else{ + return beaversJosephVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + } + } + //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used const Scalar beaversJosephVelocity(const Element& element, const SubControlVolume& scv, diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 121d9104d2..42993bd320 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,10 +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. - bool bj = currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())); - bool nTangential = currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex())); + bool slipCondition = currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex())); bool neumann = lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())); - if ( (bj || nTangential) && lateralFaceBoundaryTypes && neumann) + if ( slipCondition && lateralFaceBoundaryTypes && neumann) { const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, @@ -374,11 +373,10 @@ public: // Check the consistency of the boundary conditions, exactly one of the following must be set if (lateralFaceBoundaryTypes) { - std::bitset<4> admittableBcTypes; + 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(3, lateralFaceBoundaryTypes->isNTangential(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 " @@ -505,15 +503,11 @@ 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 if (bcTypes.isNTangential(Indices::velocity(lateralFace.directionIndex()))){ - return VelocityGradients::nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - } else return faceVars.velocityLateralInside(localSubFaceIdx); } @@ -580,8 +574,7 @@ private: { if (!scvf.boundary() || currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) || - currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())) || - currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralFace.directionIndex())) + currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralFace.directionIndex())) ) { const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh index 3241255d26..30316395dd 100644 --- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh @@ -583,8 +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 lateralFaceHasNTangential = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isNTangential(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(); @@ -593,12 +592,8 @@ private: if (useZeroGradient) return velocitySelf; - if (lateralFaceHasBJS) - return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - - else if(lateralFaceHasNTangential) - return VelocityGradients::nTangentialVelocityAtLateralScvf(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 3fa41ae2f3..a7344666b0 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -121,15 +121,10 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())]; } - else if (lateralFaceBoundaryTypes->isNTangential(Indices::velocity(scvf.directionIndex()))){ - return nTangentialVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, + else if (lateralFaceBoundaryTypes->isSlipCondition(Indices::velocity(scvf.directionIndex()))){ + return slipVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } - else if (lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))) - { - return beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - } else DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center()); }(); @@ -204,15 +199,10 @@ public: const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())]; } - else if (currentScvfBoundaryTypes->isNTangential(Indices::velocity(lateralScvf.directionIndex()))){ - return nTangentialVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, + else if (currentScvfBoundaryTypes->isSlipCondition(Indices::velocity(lateralScvf.directionIndex()))){ + return slipVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); } - else if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex()))) - { - return beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - } else DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center()); }(); @@ -226,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. * ^ ^ * | | @@ -246,50 +236,7 @@ public: * */ template - static Scalar beaversJosephVelocityAtCurrentScvf(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf, - const FaceVariables& faceVars, - const std::optional& currentScvfBoundaryTypes, - const std::optional& lateralFaceBoundaryTypes, - const std::size_t localSubFaceIdx) - { - const auto eIdx = scvf.insideScvIdx(); - const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx); - - 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 - // (towards the current scvf). - static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - - if (unsymmetrizedGradientForBJ) - return 0.0; - - if (lateralScvf.boundary()) - { - if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) || - lateralFaceBoundaryTypes->isBeaversJoseph(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); - } - - template - static Scalar nTangentialVelocityAtCurrentScvf(const Problem& problem, + static Scalar slipVelocityAtCurrentScvf(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, @@ -304,38 +251,38 @@ 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(), + static const bool unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - if (unsymmetrizedGradientForBJ) + if (unsymmetrizedGradientForIC) return 0.0; if (lateralScvf.boundary()) { if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) || - lateralFaceBoundaryTypes->isNTangential(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.nTangentialVelocity(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. :: @@ -349,50 +296,7 @@ public: * \endverbatim */ template - static Scalar beaversJosephVelocityAtLateralScvf(const Problem& problem, - const Element& element, - const FVElementGeometry& fvGeometry, - const SubControlVolumeFace& scvf, - const FaceVariables& faceVars, - const std::optional& currentScvfBoundaryTypes, - const std::optional& lateralFaceBoundaryTypes, - const std::size_t localSubFaceIdx) - { - const auto eIdx = scvf.insideScvIdx(); - const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - const Scalar innerParallelVelocity = faceVars.velocitySelf(); - - 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 - // (towards the current scvf). - static const bool unsymmetrizedGradientForBJ = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - - if (unsymmetrizedGradientForBJ) - return 0.0; - - if (scvf.boundary()) - { - if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) || - currentScvfBoundaryTypes->isBeaversJoseph(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); - } - - template - static Scalar nTangentialVelocityAtLateralScvf(const Problem& problem, + static Scalar slipVelocityAtLateralScvf(const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, @@ -407,31 +311,31 @@ 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(), + static const bool unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - if (unsymmetrizedGradientForBJ) + if (unsymmetrizedGradientForIC) return 0.0; if (scvf.boundary()) { if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) || - currentScvfBoundaryTypes->isNTangential(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.nTangentialVelocity(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 1b5721d002..c08f15d44e 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 f40bf96f36..f3176dfe5b 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -80,6 +80,28 @@ public: using ParentType::couplingPhaseIdx; + /*! + * \brief Returns the momentum flux across the coupling boundary. + * + * Calls (old/new)MomentumCouplingCondition depending on the value of the parameter "Problem.NewIc" + * Defaults to oldMomentumCouplingCondition + * + */ + template + Scalar momentumCouplingCondition(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& stokesElemVolVars, + const ElementFaceVariables& stokesElemFaceVars, + const SubControlVolumeFace& scvf) const + { + if (getParamFromGroup("Problem", "NewIc", false)){ + return newMomentumCouplingCondition(element, fvGeometry, stokesElemVolVars, stokesElemFaceVars, scvf); + } + else{ + return oldMomentumCouplingCondition(element, fvGeometry, stokesElemVolVars, stokesElemFaceVars, scvf); + } + } + /*! * \brief Returns the momentum flux across the coupling boundary. * @@ -88,7 +110,7 @@ public: * */ template - Scalar momentumCouplingCondition(const Element& element, + Scalar oldMomentumCouplingCondition(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& stokesElemVolVars, const ElementFaceVariables& stokesElemFaceVars, @@ -149,7 +171,7 @@ public: * */ template - Scalar nMomentumCouplingCondition(const Element& element, + Scalar newMomentumCouplingCondition(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& stokesElemVolVars, const ElementFaceVariables& stokesElemFaceVars, @@ -223,10 +245,10 @@ public: this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - const bool unsymmetrizedGradientForBJ = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); //TODO: Remove calculation above in this case - if (unsymmetrizedGradientForBJ){ + if (unsymmetrizedGradientForIC){ velocityGrad_ij = 0.0; } @@ -240,6 +262,23 @@ public: return momentumFlux; } + /*! + * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law + * + * Calls standardPorousMediumVelocity/newPorousMediumInterfaceVelocity depending on the value of the parameter "Problem.NewIc" + * Defaults to standardPorousMediumVelocity + * + */ + VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + { + if (getParamFromGroup("Problem", "NewIc", false)){ + return newPorousMediumInterfaceVelocity(element, scvf); + } + else{ + return standardPorousMediumVelocity(element, scvf); + } + } + /*! * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law * @@ -247,7 +286,7 @@ public: * stokes-darcy interface. We use darcys law and perform an integral average over all coupling segments. * */ - VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + VelocityVector standardPorousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { static constexpr int darcyDim = GridGeometry::GridView::dimension; using JacobianType = Dune::FieldMatrix; -- GitLab From dc8771bb4141fce25eea49dd6a525874af0fa19a Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 28 Jul 2020 15:54:21 +0200 Subject: [PATCH 23/40] Added description --- dumux/freeflow/navierstokes/problem.hh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 7424b9dfb0..05fa8c64da 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -212,13 +212,15 @@ public: return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf)); } - //TODO: Elissa, maybe you can write a small description, see above for these functions /*! - * \brief Returns the ... value, + * \brief Returns the scale separation value eps = (macroscopicLength/LengthOfPeriodicityCell). */ Scalar epsInterface(const SubControlVolumeFace& scvf) const { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal 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 new interface conditions, the factorNTangential value must be returned in the acutal problem");} -- GitLab From 7dbc148de29b1a42a4956d737f5bf4d152a2f775 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 10:33:14 +0200 Subject: [PATCH 24/40] uses static variables for getParam, comments --- dumux/freeflow/navierstokes/problem.hh | 5 ++--- .../boundary/stokesdarcy/box/couplingdata.hh | 17 ++++++----------- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 05fa8c64da..4d30d36d9b 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -241,7 +241,8 @@ public: const Scalar velocitySelf, const Scalar tangentialVelocityGradient) const { - if (getParamFromGroup("Problem", "NewIc", false)){ + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + if (newIc_){ return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); } else{ @@ -271,7 +272,6 @@ public: + velocitySelf) / (betaBJ*distanceNormalToBoundary + 1.0); } - //TODO: viscosity? //! helper function to evaluate the slip velocity on the boundary when the new tangential condition is used const Scalar nTangentialVelocity(const Element& element, const SubControlVolume& scv, @@ -280,7 +280,6 @@ public: const Scalar velocitySelf, //vel auf stehendem const Scalar tangentialVelocityGradient) const //dv/dx=0 { - //const Scalar viscosity = elemVolVars[scvf.insideScvIdx()].Viscosity(); const Scalar factor = -1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index f3176dfe5b..ea9b949212 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -94,7 +94,8 @@ public: const ElementFaceVariables& stokesElemFaceVars, const SubControlVolumeFace& scvf) const { - if (getParamFromGroup("Problem", "NewIc", false)){ + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + if (newIc_){ return newMomentumCouplingCondition(element, fvGeometry, stokesElemVolVars, stokesElemFaceVars, scvf); } else{ @@ -245,7 +246,7 @@ public: this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); //TODO: Remove calculation above in this case if (unsymmetrizedGradientForIC){ @@ -271,7 +272,8 @@ public: */ VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - if (getParamFromGroup("Problem", "NewIc", false)){ + static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + if (newIc_){ return newPorousMediumInterfaceVelocity(element, scvf); } else{ @@ -346,14 +348,11 @@ public: } } //account for gravity - //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g the correct sign? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } //Add the integrated segment velocity to the sum: v+= -weight_k * sqrt(det(A^T*A))*K/mu*gradP - //TODO: which fits dumux style better? K.usmv(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), gradP, velocity); - //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx), mv(K,gradP)); } intersectionLength += data.segmentGeometry.volume(); } @@ -434,15 +433,11 @@ public: } } //account for gravity - //TODO by Lars: Gravity changes darcy's law for calculating the velocity: Why is -rho*g the correct sign? if (enableGravity){ gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } //Add the integrated segment velocity to the sum: v+= -w_k * sqrt(det(A^T*A))*eps**2*M/mu*gradP - //TODO: which fits dumux style better? - M.usmv(qp.weight()*data.segmentGeometry.integrationElement(ipLocal)*data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); - //alternativ: velocity.axpy(-qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, mv(M,gradP)); - } + M.usmv(qp.weight()*data.segmentGeometry.integrationElement(ipLocal)/data.volVars.viscosity(darcyPhaseIdx)*epsInterface*epsInterface, gradP, velocity); } intersectionLength += data.segmentGeometry.volume(); } } -- GitLab From d3c0e0bcb0d0fc2f9d2f6a14f8ff1abe22aba902 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 16:03:33 +0200 Subject: [PATCH 25/40] Trying to deprecate BeaversJoseph methods --- dumux/discretization/staggered/freeflow/boundarytypes.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index b55640dc5c..f5cde7ccb5 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -99,6 +99,7 @@ public: * \brief Set a boundary condition for a single equation to * Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.). */ + [[deprecated("Use setSlipCondition")]] void setBeaversJoseph(unsigned eqIdx) { setSlipCondition(eqIdx); @@ -122,6 +123,7 @@ public: * * \param eqIdx The index of the equation */ + [[deprecated("Use isSlipCondition")]] bool isBeaversJoseph(unsigned eqIdx) const { return isSlipCondition(eqIdx); } @@ -139,6 +141,7 @@ public: * \brief Returns true if some equation is used to specify a * Beavers-Joseph(-Saffman) boundary condition. */ + [[deprecated("Use hasSlipCondition")]] bool hasBeaversJoseph() const { return hasSlipCondition(); -- GitLab From a2fc653df39912720322776d8748aec7e7958237 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 16:13:43 +0200 Subject: [PATCH 26/40] Combining BeaversJoseph and nTangential --- dumux/freeflow/navierstokes/problem.hh | 48 ++++++-------------------- 1 file changed, 10 insertions(+), 38 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 4d30d36d9b..ff99d954d6 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -233,54 +233,26 @@ public: return VelocityVector(0.0); } - //! helper function to evaluate the slip velocity on the boundary + //! 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) { static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); + // du/dy + dv/dx = factor * (u_boundary-uPM) + Scalar factor; if (newIc_){ - return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + //return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + factor = -1.0 / asImp_().epsInterface(faceOnPorousBoundary) / asImp_().factorNTangential(faceOnPorousBoundary); + } else{ - return beaversJosephVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + //return beaversJosephVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + factor = asImp_().betaBJ(element, faceOnPorousBoundary); //beta = alpha/sqrt(K) } - } - - //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used - const Scalar beaversJosephVelocity(const Element& element, - const SubControlVolume& scv, - const SubControlVolumeFace& ownScvf, - const SubControlVolumeFace& faceOnPorousBoundary, - const Scalar velocitySelf, - const Scalar tangentialVelocityGradient) const - { - // du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM) - // beta = alpha/sqrt(K) - const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary); - 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); - } - - //! helper function to evaluate the slip velocity on the boundary when the new tangential condition is used - const Scalar nTangentialVelocity(const Element& element, - const SubControlVolume& scv, - const SubControlVolumeFace& ownScvf, //stehendes - const SubControlVolumeFace& faceOnPorousBoundary, //liegendes - const Scalar velocitySelf, //vel auf stehendem - const Scalar tangentialVelocityGradient) const //dv/dx=0 - { - const Scalar factor = -1.0/(asImp_().epsInterface(faceOnPorousBoundary) * asImp_().factorNTangential(faceOnPorousBoundary)); const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); // create a unit normal vector oriented in positive coordinate direction -- GitLab From 1cb7b8be1bf0448d4ac6d6d378585345d39de107 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 16:23:32 +0200 Subject: [PATCH 27/40] Accepting both parameters: unsymmGradForBJ/IC --- .../staggered/velocitygradients.hh | 11 +++++-- .../boundary/stokesdarcy/box/couplingdata.hh | 32 ++++--------------- 2 files changed, 14 insertions(+), 29 deletions(-) diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh index a7344666b0..e8650c4b90 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -254,8 +254,11 @@ public: // 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 unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), + static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + //TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + static const bool unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); if (unsymmetrizedGradientForIC) return 0.0; @@ -314,9 +317,11 @@ public: // 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 unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), + static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - + // TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + static const bool unsymmetrizedGradientForIC = getParamFromGroup(problem.paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); if (unsymmetrizedGradientForIC) return 0.0; diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index ea9b949212..8ad85f561d 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -191,33 +191,13 @@ public: const auto& darcyFvGeometry = data.fvGeometry; const auto& localBasis = darcyFvGeometry.feLocalBasis(); - // do second order integration as box provides linear functions - static constexpr int darcyDim = GridGeometry::GridView::dimension; - const auto& rule = Dune::QuadratureRules::rule(data.segmentGeometry.type(), 2); - 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); - - std::vector> shapeValues; - localBasis.evaluateFunction(ipElementLocal, shapeValues); - - Scalar pressure = 0.0; - for (const auto& scv : scvs(data.fvGeometry)) - pressure += elemVolVars[scv].pressure(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; - - momentumFlux += pressure*data.segmentGeometry.integrationElement(qp.position())*qp.weight(); - } - } - } - momentumFlux /= scvf.area(); - - // normalize pressure - if(getPropValue, Properties::NormalizePressure>()) - momentumFlux -= this->couplingManager().problem(stokesIdx).initial(scvf)[Indices::pressureIdx]; - //######## New stokes contribution ################# + static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); + // TODO: how to deprecate unsymmBeaverJoseph? + // Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); const std::size_t numSubFaces = scvf.pairData().size(); // Account for all sub faces -- GitLab From 6a40bd50da05276c5ed0ca9f74cb9b0699dcf332 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 16:27:19 +0200 Subject: [PATCH 28/40] Simplified newMomentum with oldMomentum --- .../boundary/stokesdarcy/box/couplingdata.hh | 21 ++++++------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 8ad85f561d..2651f0de72 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -162,8 +162,6 @@ public: return momentumFlux; } - //TODO by Lars: Review! - //TODO by Lars: Use momentumCouplingCondition(...)/introduce new method to remove dulplicated code? /*! * \brief Returns the momentum flux across the coupling boundary which is calculated according to the new interface condition * @@ -179,17 +177,9 @@ public: const SubControlVolumeFace& scvf) const { Scalar momentumFlux(0.0); - const auto& stokesContext = this->couplingManager().stokesCouplingContextVector(element, scvf); //######## darcy contribution ################# - // integrate darcy pressure over each coupling segment and average - for (const auto& data : stokesContext) - { - 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(); + momentumFlux = oldMomentumCouplingCondition(element,fvGeometry,stokesElemVolVars,stokesElemFaceVars,scvf); + momentumFlux*= scvf.directionSign(); //Revert sign change, should be applied to total flux //######## New stokes contribution ################# static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), @@ -209,7 +199,9 @@ public: // Create a boundaryTypes object (will be empty if not at a boundary) std::optional> currentScvfBoundaryTypes; if (scvf.boundary()) + { currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); + } std::optional> lateralFaceBoundaryTypes; if (lateralScvf.boundary()) @@ -226,10 +218,9 @@ public: this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); //TODO: Remove calculation above in this case - if (unsymmetrizedGradientForIC){ + if (unsymmetrizedGradientForIC) + { velocityGrad_ij = 0.0; } -- GitLab From 0ebb6ba9cc9f71a5c3b1088943b6edd6633447fd Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 7 Aug 2020 16:28:45 +0200 Subject: [PATCH 29/40] Combining velocity functions --- .../boundary/stokesdarcy/box/couplingdata.hh | 321 +++++++----------- 1 file changed, 119 insertions(+), 202 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 2651f0de72..f554241d7d 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -176,242 +176,159 @@ public: const ElementFaceVariables& stokesElemFaceVars, const SubControlVolumeFace& scvf) const { - Scalar momentumFlux(0.0); - //######## darcy contribution ################# + Scalar momentumFlux(0.0); + //######## darcy contribution ################# momentumFlux = oldMomentumCouplingCondition(element,fvGeometry,stokesElemVolVars,stokesElemFaceVars,scvf); momentumFlux*= scvf.directionSign(); //Revert sign change, should be applied to total flux - //######## New stokes contribution ################# + //######## New stokes contribution ################# static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); // TODO: how to deprecate unsymmBeaverJoseph? // Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); - const std::size_t numSubFaces = scvf.pairData().size(); + const std::size_t numSubFaces = scvf.pairData().size(); - // Account for all sub faces - for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) - { - const auto eIdx = scvf.insideScvIdx(); - const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); + // Account for all sub faces + for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) + { + const auto eIdx = scvf.insideScvIdx(); + const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - // Create a boundaryTypes object (will be empty if not at a boundary) - std::optional> currentScvfBoundaryTypes; - if (scvf.boundary()) + // Create a boundaryTypes object (will be empty if not at a boundary) + std::optional> currentScvfBoundaryTypes; + if (scvf.boundary()) { - currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); + currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); } - std::optional> lateralFaceBoundaryTypes; - if (lateralScvf.boundary()) - { - lateralFaceBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, lateralScvf)); - } + 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); - Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( - this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + // Get velocity gradients + const Scalar velocityGrad_ji = StokesVelocityGradients::velocityGradJI( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); + Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( + this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], + currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - //TODO: Remove calculation above in this case + //TODO: Remove calculation above in this case if (unsymmetrizedGradientForIC) { - velocityGrad_ij = 0.0; - } - - // 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; - } + velocityGrad_ij = 0.0; + } - /*! - * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law - * - * Calls standardPorousMediumVelocity/newPorousMediumInterfaceVelocity depending on the value of the parameter "Problem.NewIc" - * Defaults to standardPorousMediumVelocity - * - */ - VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const - { - static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); - if (newIc_){ - return newPorousMediumInterfaceVelocity(element, scvf); - } - else{ - return standardPorousMediumVelocity(element, scvf); + // 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 tangential coupling at the - * stokes-darcy interface. We use darcys law and perform an integral average over all coupling segments. + * 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 standardPorousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const + VelocityVector porousMediumVelocity(const Element& element, const SubControlVolumeFace& scvf) const { - 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; - //TODO: Is this needed? - shapeValues.clear(); - shapeDerivatives.clear(); - - //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)); - } - //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; - } + 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; - /*! - * \brief Returns the averaged velocity vector at the interface of the porous medium according to darcys law with a different permeability tensor - * - * For the new tangential interface condition by Elissa Eggenweiler, a porous medium velocity with altered permeability tensor needs to be evaluated. - * We use darcys law and perform an integral average over all coupling segments. - * - */ - VelocityVector newPorousMediumInterfaceVelocity(const Element& element, const SubControlVolumeFace& scvf) const - { - 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) + 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) { - 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; - //TODO: Is this needed? - shapeValues.clear(); - shapeDerivatives.clear(); - - // darcy spatial dependent parameters - const auto& epsInterface = this->couplingManager().problem(darcyIdx).spatialParams().epsInterfaceAtPos(ipGlobal); - const auto& M = this->couplingManager().problem(darcyIdx).spatialParams().matrixNTangentialAtPos(ipGlobal); - - //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)); + //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; + //TODO: Is this needed? -> Numerisch ausprobieren, sollte überschrieben werden, siehe doxygen dune + shapeValues.clear(); + shapeDerivatives.clear(); + + //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(); } - //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); } - intersectionLength += data.segmentGeometry.volume(); } - } velocity /= intersectionLength; //averaging return velocity; } -- GitLab From 957a45d1c29013981dea764d13cc9c21f8452d70 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Sat, 8 Aug 2020 12:35:13 +0200 Subject: [PATCH 30/40] [Deprecated] Improved message --- dumux/discretization/staggered/freeflow/boundarytypes.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index f5cde7ccb5..c4689f2c28 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -99,7 +99,7 @@ public: * \brief Set a boundary condition for a single equation to * Beavers-Joseph(-Saffmann) (special case of Dirichlet b.c.). */ - [[deprecated("Use setSlipCondition")]] + [[deprecated("This method will be removed after release (3.n). Use setSlipCondition instead!")]] void setBeaversJoseph(unsigned eqIdx) { setSlipCondition(eqIdx); @@ -123,7 +123,7 @@ public: * * \param eqIdx The index of the equation */ - [[deprecated("Use isSlipCondition")]] + [[deprecated("This method will be removed after release (3.n). Use isSlipCondition instead!")]] bool isBeaversJoseph(unsigned eqIdx) const { return isSlipCondition(eqIdx); } @@ -141,7 +141,7 @@ public: * \brief Returns true if some equation is used to specify a * Beavers-Joseph(-Saffman) boundary condition. */ - [[deprecated("Use hasSlipCondition")]] + [[deprecated("This method will be removed after release (3.n). Use hasSlipCondition instead!")]] bool hasBeaversJoseph() const { return hasSlipCondition(); -- GitLab From de8144a8f25911e622642de62634c3f4b7cf5098 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Sat, 8 Aug 2020 12:37:39 +0200 Subject: [PATCH 31/40] [styleguide] braces and code in comments --- dumux/freeflow/navierstokes/problem.hh | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index ff99d954d6..3baf917571 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -216,13 +216,17 @@ public: * \brief Returns the scale separation value eps = (macroscopicLength/LengthOfPeriodicityCell). */ Scalar epsInterface(const SubControlVolumeFace& scvf) const - { DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem");} + { + DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal 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 new interface conditions, the factorNTangential value must be returned in the acutal problem");} + { + DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the acutal problem"); + } /*! @@ -244,13 +248,12 @@ public: static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); // du/dy + dv/dx = factor * (u_boundary-uPM) Scalar factor; - if (newIc_){ - //return nTangentialVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + if (newIc_) + { factor = -1.0 / asImp_().epsInterface(faceOnPorousBoundary) / asImp_().factorNTangential(faceOnPorousBoundary); - } - else{ - //return beaversJosephVelocity(element, scv,ownScvf, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient); + else + { factor = asImp_().betaBJ(element, faceOnPorousBoundary); //beta = alpha/sqrt(K) } const Scalar distanceNormalToBoundary = (faceOnPorousBoundary.center() - scv.center()).two_norm(); -- GitLab From 44e827f5c71633dca756cd3091c64c4b975a2ec3 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Sat, 8 Aug 2020 12:41:18 +0200 Subject: [PATCH 32/40] [combine] Combined porousMediumVelocity & cleanup --- .../boundary/stokesdarcy/box/couplingdata.hh | 149 ++++++------------ 1 file changed, 50 insertions(+), 99 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index f554241d7d..1556d58515 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -28,9 +28,8 @@ #include #include -#include //TODO by Lars: Why is this the right coupling manager, not the one in this folder? +#include //TODO: Why is this the right coupling manager, not the one in this folder? -//Needed for nMomentumCouplingCondition #include #include #include @@ -63,9 +62,7 @@ class StokesDarcyCouplingDataBoxBase : public StokesDarcyCouplingDataImplementat static constexpr auto stokesIdx = CouplingManager::stokesIdx; static constexpr auto darcyIdx = CouplingManager::darcyIdx; - // Needed for velocityPorousMedium using VelocityVector = typename Element::Geometry::GlobalCoordinate; - // Needed for nMomentum template using BoundaryTypes = GetPropType, Properties::BoundaryTypes>; using StokesVelocityGradients = StaggeredVelocityGradients, BoundaryTypes, Indices>; @@ -83,8 +80,11 @@ public: /*! * \brief Returns the momentum flux across the coupling boundary. * - * Calls (old/new)MomentumCouplingCondition depending on the value of the parameter "Problem.NewIc" - * Defaults to oldMomentumCouplingCondition + * 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. For the new normal momentum coupling, the stokes term -viscosity*Nsbl*tau.T*grad(v_ff)*n is added. * */ template @@ -95,31 +95,11 @@ public: const SubControlVolumeFace& scvf) const { static const bool newIc_ = getParamFromGroup("Problem", "NewIc", false); - if (newIc_){ - return newMomentumCouplingCondition(element, fvGeometry, stokesElemVolVars, stokesElemFaceVars, scvf); - } - else{ - return oldMomentumCouplingCondition(element, fvGeometry, stokesElemVolVars, stokesElemFaceVars, scvf); - } - } - /*! - * \brief Returns the momentum flux across the coupling boundary. - * - * For the normal momentum coupling, the porous medium side of the coupling condition - * is evaluated, i.e. -[p n]^pm. - * - */ - template - Scalar oldMomentumCouplingCondition(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& stokesElemVolVars, - const ElementFaceVariables& stokesElemFaceVars, - const SubControlVolumeFace& scvf) const - { 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) { @@ -157,78 +137,51 @@ public: if(getPropValue, Properties::NormalizePressure>()) momentumFlux -= this->couplingManager().problem(stokesIdx).initial(scvf)[Indices::pressureIdx]; - momentumFlux *= scvf.directionSign(); - - return momentumFlux; - } - - /*! - * \brief Returns the momentum flux across the coupling boundary which is calculated according to the new interface condition - * - * For the new normal momentum coupling, the porous medium side and also the stokes side is evaluated. - * [p]^pm + N_s^{bl} \tau T n - * - */ - template - Scalar newMomentumCouplingCondition(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& stokesElemVolVars, - const ElementFaceVariables& stokesElemFaceVars, - const SubControlVolumeFace& scvf) const - { - Scalar momentumFlux(0.0); - //######## darcy contribution ################# - momentumFlux = oldMomentumCouplingCondition(element,fvGeometry,stokesElemVolVars,stokesElemFaceVars,scvf); - momentumFlux*= scvf.directionSign(); //Revert sign change, should be applied to total flux - - //######## New stokes contribution ################# - static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + if (newIc_) + { + //######## New stokes contribution ################# + static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - // TODO: how to deprecate unsymmBeaverJoseph? - // Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired - static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), + // TODO: how to deprecate unsymmBeaverJoseph? + // Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + static const bool unsymmetrizedGradientForIC = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); - const std::size_t numSubFaces = scvf.pairData().size(); + const std::size_t numSubFaces = scvf.pairData().size(); - // Account for all sub faces - for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) - { - const auto eIdx = scvf.insideScvIdx(); - const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - - // Create a boundaryTypes object (will be empty if not at a boundary) - std::optional> currentScvfBoundaryTypes; - if (scvf.boundary()) + // Account for all sub faces + for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) { - currentScvfBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, scvf)); - } + const auto eIdx = scvf.insideScvIdx(); + const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); - std::optional> lateralFaceBoundaryTypes; - if (lateralScvf.boundary()) - { - lateralFaceBoundaryTypes.emplace(this->couplingManager().problem(stokesIdx).boundaryTypes(element, lateralScvf)); - } + // Create a boundaryTypes object (will be empty if not at a boundary) + std::optional> 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); - Scalar velocityGrad_ij = StokesVelocityGradients::velocityGradIJ( - this->couplingManager().problem(stokesIdx), element, fvGeometry, scvf , stokesElemFaceVars[scvf], - currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); - //TODO: Remove calculation above in this case - if (unsymmetrizedGradientForIC) - { - velocityGrad_ij = 0.0; - } + // 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); + // 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; @@ -290,30 +243,28 @@ public: //reset pressure gradient and rho at this qp gradP=0.0; rho=0.0; - //TODO: Is this needed? -> Numerisch ausprobieren, sollte überschrieben werden, siehe doxygen dune - shapeValues.clear(); - shapeDerivatives.clear(); //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){ + if (enableGravity) + { rho += elemVolVars[scv].density(darcyPhaseIdx)*shapeValues[scv.indexInElement()][0]; } } //account for gravity - if (enableGravity){ + if (enableGravity) + { gradP.axpy(-rho, this->couplingManager().problem(darcyIdx).spatialParams().gravity(ipGlobal)); } - if(newIc_){ + 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); -- GitLab From 389334618b207ad0f49f749707099e37be8913d7 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Fri, 14 Aug 2020 00:27:27 +0200 Subject: [PATCH 33/40] [Test] Add test for BJ and NewIC, wip Started to implement test Vtk-Data for comparison missing! --- .../1p_1p/BJandNewIC/1pspatialparams.hh | 127 ++++++ .../1p_1p/BJandNewIC/CMakeLists.txt | 30 ++ .../stokesdarcy/1p_1p/BJandNewIC/main.cc | 224 +++++++++++ .../1p_1p/BJandNewIC/params_BJ.input | 46 +++ .../1p_1p/BJandNewIC/params_newIC.input | 46 +++ .../1p_1p/BJandNewIC/problem_darcy.hh | 365 +++++++++++++++++ .../1p_1p/BJandNewIC/problem_stokes.hh | 377 ++++++++++++++++++ .../boundary/stokesdarcy/1p_1p/CMakeLists.txt | 1 + 8 files changed, 1216 insertions(+) create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/1pspatialparams.hh create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh create mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/1pspatialparams.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/1pspatialparams.hh new file mode 100644 index 0000000000..a861e5f18e --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/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/BJandNewIC/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt new file mode 100644 index 0000000000..fa10a840c6 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt @@ -0,0 +1,30 @@ +add_input_file_links() + +add_executable(test_md_boundary_darcy1p_stokes1p_New EXCLUDE_FROM_ALL main.cc) + +dumux_add_test(NAME test_md_boundary_darcy1p_stokes1p_BJ + LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes + TARGET test_md_boundary_darcy1p_stokes1p_New + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_BJ_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_BJ_stokes-00001.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_BJ_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_BJ_darcy-00001.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_New params_BJ.input + -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_BJ") + +dumux_add_test(NAME test_md_boundary_darcy1p_stokes1p_NewIC + LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes + TARGET test_md_boundary_darcy1p_stokes1p_New + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_NewIC_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_NewIC_stokes-00001.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_NewIC_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_NewIC_darcy-00001.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_New params_newIC.input + -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_NewIC" + --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17}) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc new file mode 100644 index 0000000000..09b96b3a24 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc @@ -0,0 +1,224 @@ +// -*- 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 "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 + +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 + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + GetPropType::initOutputModule(stokesVtkWriter); + //stokesVtkWriter.write(0.0); + + VtkOutputModule> darcyVtkWriter(*darcyGridVariables, sol[porousMediumIdx], darcyProblem->name()); + using DarcyVelocityOutput = GetPropType; + darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); + GetPropType::initOutputModule(darcyVtkWriter); + //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); + + //////////////////////////////////////////////////////////// + // 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/BJandNewIC/params_BJ.input b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input new file mode 100644 index 0000000000..6304ebfa85 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input @@ -0,0 +1,46 @@ +[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 + +[Vtk] +AddVelocity = 1 + +[Component] +LiquidDensity = 1.0 +LiquidKinematicViscosity = 1.0 + +[Assembly] +NumericDifference.BaseEpsilon = 1e-4 diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input new file mode 100644 index 0000000000..7d2e010528 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input @@ -0,0 +1,46 @@ +[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 = 0.0 + +[Vtk] +OutputName = test_md_boundary_stokes1p_darcy1p_newIC + +[Problem] +EnableGravity = false +EnableInertiaTerms = false + +[Vtk] +AddVelocity = 1 + +[Component] +LiquidDensity = 1.0 +LiquidKinematicViscosity = 2.0 + +[Assembly] +NumericDifference.BaseEpsilon = 1e-4 diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh new file mode 100644 index 0000000000..6989c3592a --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/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)[Indices::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) + PrimaryVariables analyticalSolutionBJSymmetrized_(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] = 2.0*x*x + 2.0*x*y - 6.0*x - 2.0*y*y + 3.0*y - 5.0; + //sol[Indices::velocityYIdx] = 1.0*x*x + 4.0*x*y - 9.0*x - 1.0*y*y - 1.0; + sol[Indices::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) + 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] = x*(y-1.0) + (x-1.0)*(y-1.0) - 3.0; + //sol[Indices::velocityYIdx] = x*(x-1.0) - (y-1.0)*(y-1.0) - 2.0; + sol[Indices::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/BJandNewIC/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh new file mode 100644 index 0000000000..4acdc3ae70 --- /dev/null +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh @@ -0,0 +1,377 @@ +// -*- 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; + + 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 rhsBJSymmetrized_(globalPos); + case TestCase::NewICNonSymmetrized: + return rhsNewICNonSymmetrized_(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; + } + + + static constexpr Scalar eps_ = 1e-7; + std::string problemName_; + std::shared_ptr couplingManager_; + TestCase testCase_; +}; +} // end namespace Dumux + +#endif // DUMUX_STOKES_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt index 45acb89ecd..84dfb3bf77 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(BJandNewIC) add_input_file_links() -- GitLab From 804698c2d08fb1e7d706a750b569b1903811c585 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Mon, 21 Sep 2020 01:15:02 +0200 Subject: [PATCH 34/40] [CleanUp] Fix typo --- dumux/freeflow/navierstokes/problem.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 3baf917571..2e3aeb3fa2 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"); } /*! @@ -217,7 +217,7 @@ public: */ Scalar epsInterface(const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the acutal problem"); + DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the actual problem"); } /*! @@ -225,7 +225,7 @@ public: */ Scalar factorNTangential(const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the acutal problem"); + DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the actual problem"); } -- GitLab From 4638307a17c75762376b2eb6646171011989965d Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Mon, 21 Sep 2020 01:26:49 +0200 Subject: [PATCH 35/40] [Warnings] deprecation and parameter combination Added deprecation warning for old parameters. Added warning if using newIc and symmetricGradient --- .../staggered/velocitygradients.hh | 26 ++++++++++++++----- .../boundary/stokesdarcy/box/couplingdata.hh | 21 ++++++++++++--- 2 files changed, 37 insertions(+), 10 deletions(-) diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh index e8650c4b90..d1e347e5f2 100644 --- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh +++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh @@ -254,9 +254,16 @@ public: // 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 unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - //TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + 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); @@ -317,9 +324,16 @@ public: // 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 unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(problem.paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - // TODO: Deprecate and replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + 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) diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index 1556d58515..c2f9c5a55e 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -140,12 +140,25 @@ public: if (newIc_) { //######## New stokes contribution ################# - static const bool unsymmetrizedGradientForBeaversJoseph = getParamFromGroup(this->couplingManager().problem(stokesIdx).paramGroup(), - "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false); - // TODO: how to deprecate unsymmBeaverJoseph? - // Replace unsymmetrizedGradientForBeaversJoseph below by false, when deprecation period expired + 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 new interface conditions are not verified for symmetrized stress tensors" < Date: Mon, 21 Sep 2020 01:29:44 +0200 Subject: [PATCH 36/40] [test] Add parameters and fix coefficient --- .../boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input | 1 + .../stokesdarcy/1p_1p/BJandNewIC/params_newIC.input | 8 +++++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input index 6304ebfa85..431620690c 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input @@ -34,6 +34,7 @@ OutputName = test_md_boundary_stokes1p_darcy1p_BJ [Problem] EnableGravity = false EnableInertiaTerms = false +TestCase = BJSymmetrized [Vtk] AddVelocity = 1 diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input index 7d2e010528..f389442fb7 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input @@ -26,7 +26,7 @@ epsInterface = 1.0 N_s_bl = 0.5 N_1_bl = -3.0 M_1_bl = 1.0 -M_2_bl = 0.0 +M_2_bl = 1.0 [Vtk] OutputName = test_md_boundary_stokes1p_darcy1p_newIC @@ -34,6 +34,12 @@ 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 -- GitLab From 36ca959098c0cf0cd9498a2ec77a040484e98872 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Mon, 21 Sep 2020 01:30:41 +0200 Subject: [PATCH 37/40] [test] Switching to convergence test --- .../1p_1p/BJandNewIC/CMakeLists.txt | 41 ++++++++----------- 1 file changed, 16 insertions(+), 25 deletions(-) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt index fa10a840c6..7c567ee59b 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt @@ -1,30 +1,21 @@ add_input_file_links() +dune_symlink_to_source_files(FILES "./../convergencetest/convergencetest.py") -add_executable(test_md_boundary_darcy1p_stokes1p_New EXCLUDE_FROM_ALL main.cc) - -dumux_add_test(NAME test_md_boundary_darcy1p_stokes1p_BJ - LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes - TARGET test_md_boundary_darcy1p_stokes1p_New +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 ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMD_ARGS --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_BJ_stokes-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_BJ_stokes-00001.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_BJ_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_BJ_darcy-00001.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_New params_BJ.input - -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_BJ") + COMMAND ./convergencetest.py + CMD_ARGS test_md_boundary_darcy1p_stokes1p_BJNewIC params_BJ.input) -dumux_add_test(NAME test_md_boundary_darcy1p_stokes1p_NewIC - LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes - TARGET test_md_boundary_darcy1p_stokes1p_New +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 ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMD_ARGS --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_NewIC_stokes-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_NewIC_stokes-00001.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_NewIC_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_NewIC_darcy-00001.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_New params_newIC.input - -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_NewIC" - --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17}) + COMMAND ./convergencetest.py + CMD_ARGS test_md_boundary_darcy1p_stokes1p_BJNewIC params_newIC.input) + + +add_executable(test_md_boundary_darcy1p_stokes1p_BJNewIC EXCLUDE_FROM_ALL main.cc) \ No newline at end of file -- GitLab From 148ad11249ad82ed4dfb1989eb54b343162a7bbd Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Mon, 21 Sep 2020 10:32:44 +0200 Subject: [PATCH 38/40] Revert "Debugging activated" This reverts commit 0252539ad9975a8ca5b9994a1eee12a5ba2f9a44. --- cmake.opts | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake.opts b/cmake.opts index 5dedd1d416..8c5f184a1b 100644 --- a/cmake.opts +++ b/cmake.opts @@ -54,6 +54,6 @@ CMAKE_FLAGS="$SPECIFIC_COMPILER $SPECIFIC_GENERATOR $OPM_FLAGS -DCMAKE_CXX_FLAGS_RELEASE='$GXX_RELEASE_OPTS $GXX_RELEASE_WARNING_OPTS $DUMUX_DISABLE_PROP_MACROS' -DCMAKE_CXX_FLAGS_DEBUG='-O0 -g -ggdb -Wall -Wextra -Wno-unused-parameter -Wno-sign-compare $DUMUX_DISABLE_PROP_MACROS' -DCMAKE_CXX_FLAGS_RELWITHDEBINFO='$GXX_RELEASE_OPTS $GXX_RELEASE_WARNING_OPTS -g -ggdb -Wall $DUMUX_DISABLE_PROP_MACROS' --DCMAKE_BUILD_TYPE=Debug +-DCMAKE_BUILD_TYPE=Release -DENABLE_HEADERCHECK=$DUMUX_ENABLE_HEADERCHECK " -- GitLab From a3f5fcb820bd79bd85b0bbf617eef9c18bb8392d Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Mon, 21 Sep 2020 23:59:04 +0200 Subject: [PATCH 39/40] [test] lowercase filenames, fix/complete test Switched to lowercase foldername and filenames. Added the l2error calculation and fixed some bugs. --- .../boundary/stokesdarcy/1p_1p/CMakeLists.txt | 2 +- .../1pspatialparams.hh | 0 .../{BJandNewIC => bj_newic}/CMakeLists.txt | 10 +- .../1p_1p/{BJandNewIC => bj_newic}/main.cc | 179 +++++++++++++++++- .../params_bj.input} | 0 .../params_newic.input} | 7 +- .../{BJandNewIC => bj_newic}/problem_darcy.hh | 30 +-- .../problem_stokes.hh | 31 ++- 8 files changed, 227 insertions(+), 32 deletions(-) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC => bj_newic}/1pspatialparams.hh (100%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC => bj_newic}/CMakeLists.txt (58%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC => bj_newic}/main.cc (55%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC/params_BJ.input => bj_newic/params_bj.input} (100%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC/params_newIC.input => bj_newic/params_newic.input} (84%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC => bj_newic}/problem_darcy.hh (92%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{BJandNewIC => bj_newic}/problem_stokes.hh (92%) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt index 84dfb3bf77..84140c7d0e 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt @@ -1,5 +1,5 @@ add_subdirectory(convergencetest) -add_subdirectory(BJandNewIC) +add_subdirectory(bj_newic) add_input_file_links() diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/1pspatialparams.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/1pspatialparams.hh similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/1pspatialparams.hh rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/1pspatialparams.hh diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt similarity index 58% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt index 7c567ee59b..4c86072741 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/CMakeLists.txt @@ -1,21 +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 +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_BJNewIC params_BJ.input) + CMD_ARGS test_md_boundary_darcy1p_stokes1p_bj_newic params_bj.input) -dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_NewIC +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_BJNewIC params_newIC.input) + CMD_ARGS test_md_boundary_darcy1p_stokes1p_bj_newic params_newic.input) -add_executable(test_md_boundary_darcy1p_stokes1p_BJNewIC EXCLUDE_FROM_ALL main.cc) \ No newline at end of file +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/BJandNewIC/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc similarity index 55% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc index 09b96b3a24..d74a5d4699 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/main.cc @@ -46,6 +46,7 @@ #include #include #include +#include #include @@ -72,6 +73,169 @@ struct CouplingManager } // 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; @@ -149,15 +313,23 @@ int main(int argc, char** argv) try darcyGridVariables->init(sol[porousMediumIdx]); // intialize the vtk output module + using Scalar = typename Traits::Scalar; StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); GetPropType::initOutputModule(stokesVtkWriter); - //stokesVtkWriter.write(0.0); + 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); - //darcyVtkWriter.write(0.0); + 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; @@ -185,6 +357,9 @@ int main(int argc, char** argv) try stokesVtkWriter.write(1.0); darcyVtkWriter.write(1.0); + printStokesL2Error(*stokesProblem, stokesSol); + printDarcyL2Error(*darcyProblem, sol[porousMediumIdx]); + //////////////////////////////////////////////////////////// // finalize, print dumux message to say goodbye //////////////////////////////////////////////////////////// diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_bj.input similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_BJ.input rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_bj.input diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input similarity index 84% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input index f389442fb7..2e0495a12c 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/params_newIC.input +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/params_newic.input @@ -29,7 +29,7 @@ M_1_bl = 1.0 M_2_bl = 1.0 [Vtk] -OutputName = test_md_boundary_stokes1p_darcy1p_newIC +OutputName = test_md_boundary_stokes1p_darcy1p_newic [Problem] EnableGravity = false @@ -46,7 +46,4 @@ AddVelocity = 1 [Component] LiquidDensity = 1.0 -LiquidKinematicViscosity = 2.0 - -[Assembly] -NumericDifference.BaseEpsilon = 1e-4 +LiquidKinematicViscosity = 2.0 \ No newline at end of file diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh similarity index 92% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh index 6989c3592a..23a536e857 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_darcy.hh +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_darcy.hh @@ -35,7 +35,7 @@ #include "1pspatialparams.hh" -#include +#include #include namespace Dumux { @@ -96,9 +96,9 @@ class DarcySubProblem : public PorousMediumFlowProblem using CouplingManager = GetPropType; - // static constexpr auto velocityXIdx = 0; - // static constexpr auto velocityYIdx = 1; - // static constexpr auto pressureIdx = 2; + static constexpr auto velocityXIdx = 0; + static constexpr auto velocityYIdx = 1; + static constexpr auto pressureIdx = 2; enum class TestCase { @@ -191,7 +191,7 @@ public: */ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const { - const auto p = analyticalSolution(globalPos)[Indices::pressureIdx]; + const auto p = analyticalSolution(globalPos)[pressureIdx]; return PrimaryVariables(p); } @@ -300,16 +300,16 @@ public: private: // see exact solution for BJ-IC with symmetrized stress tensor (by Elissa Eggenweiler) - PrimaryVariables analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const + Dune::FieldVector analyticalSolutionBJSymmetrized_(const GlobalPosition& globalPos) const { - PrimaryVariables sol(0.0); + 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[Indices::velocityXIdx] = 2.0*x*x + 2.0*x*y - 6.0*x - 2.0*y*y + 3.0*y - 5.0; - //sol[Indices::velocityYIdx] = 1.0*x*x + 4.0*x*y - 9.0*x - 1.0*y*y - 1.0; - sol[Indices::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; + 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; } @@ -324,16 +324,16 @@ private: // see exact solution for new IC with non-symmetrized stress tensor (by Elissa Eggenweiler) - PrimaryVariables analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const + Dune::FieldVector analyticalSolutionNewICNonSymmetrized_(const GlobalPosition& globalPos) const { - PrimaryVariables sol(0.0); + 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[Indices::velocityXIdx] = x*(y-1.0) + (x-1.0)*(y-1.0) - 3.0; - //sol[Indices::velocityYIdx] = x*(x-1.0) - (y-1.0)*(y-1.0) - 2.0; - sol[Indices::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; + 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; } diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh similarity index 92% rename from test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh rename to test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh index 4acdc3ae70..219c1d3a00 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/BJandNewIC/problem_stokes.hh +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/bj_newic/problem_stokes.hh @@ -175,8 +175,20 @@ public: { BoundaryTypes values; - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); + 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)) { @@ -302,9 +314,9 @@ public: switch (testCase_) { case TestCase::BJSymmetrized: - return rhsBJSymmetrized_(globalPos); + return analyticalSolutionBJSymmetrized_(globalPos); case TestCase::NewICNonSymmetrized: - return rhsNewICNonSymmetrized_(globalPos); + return analyticalSolutionNewICNonSymmetrized_(globalPos); default: DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); } @@ -366,6 +378,17 @@ private: 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_; -- GitLab From 5fd5eadc10434828c5163874be660ba3e9d11a19 Mon Sep 17 00:00:00 2001 From: Lars Kaiser Date: Tue, 22 Sep 2020 00:04:59 +0200 Subject: [PATCH 40/40] [warnings] introduce better term for new ic This commits makes a suggestion of what term can be used when printing information about the new interface conditions. --- dumux/freeflow/navierstokes/problem.hh | 4 ++-- dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 2e3aeb3fa2..8dd33c0052 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -217,7 +217,7 @@ public: */ Scalar epsInterface(const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the epsInterface value must be returned in the actual problem"); + 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"); } /*! @@ -225,7 +225,7 @@ public: */ Scalar factorNTangential(const SubControlVolumeFace& scvf) const { - DUNE_THROW(Dune::NotImplemented, "When using the new interface conditions, the factorNTangential value must be returned in the actual problem"); + 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"); } diff --git a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh index c2f9c5a55e..fd09eaffaf 100644 --- a/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/box/couplingdata.hh @@ -155,7 +155,7 @@ public: "FreeFlow.EnableUnsymmetrizedVelocityGradientForIC", unsymmetrizedGradientForBeaversJoseph); // NewIc not verified for symmetrized gradient static bool once = [](){ - if(!unsymmetrizedGradientForIC) std::cerr << "Warning: The new interface conditions are not verified for symmetrized stress tensors" <