diff --git a/dumux/freeflow/rans/oneeq/problem.hh b/dumux/freeflow/rans/oneeq/problem.hh index 9198797847e6fcb3e3b7ba4a5cd5825520b1a83d..730eacb9daa9dc559ed40b2433b458c91323be41 100644 --- a/dumux/freeflow/rans/oneeq/problem.hh +++ b/dumux/freeflow/rans/oneeq/problem.hh @@ -64,7 +64,11 @@ class OneEqProblem : public RANSProblem using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; public: - //! The constructor sets the gravity, if desired by the user. + /* + * \brief The constructor + * \param fvGridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ OneEqProblem(std::shared_ptr fvGridGeometry, const std::string& paramGroup = "") : ParentType(fvGridGeometry, paramGroup) { diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh index c1b2bde7a231a22945fe5923c5599bc1c3f52113..056a0b3a1df7065318c61c27513e6c91143476f1 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh @@ -80,8 +80,13 @@ class KEpsilonProblem : public RANSProblem static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes"); public: + static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); - //! The constructor sets the gravity, if desired by the user. + /* + * \brief The constructor + * \param fvGridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ KEpsilonProblem(std::shared_ptr fvGridGeometry, const std::string& paramGroup = "") : ParentType(fvGridGeometry, paramGroup) { diff --git a/dumux/freeflow/rans/twoeq/komega/problem.hh b/dumux/freeflow/rans/twoeq/komega/problem.hh index a4ee39561e6793931e99417861f72b8fede6260e..dd1aac3d0adec0bb9191cca706efc674d8456625 100644 --- a/dumux/freeflow/rans/twoeq/komega/problem.hh +++ b/dumux/freeflow/rans/twoeq/komega/problem.hh @@ -63,6 +63,11 @@ class KOmegaProblem : public RANSProblem using DimVector = typename Element::Geometry::GlobalCoordinate; public: + /* + * \brief The constructor + * \param fvGridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ KOmegaProblem(std::shared_ptr fvGridGeometry, const std::string& paramGroup = "") : ParentType(fvGridGeometry, paramGroup) { diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh index 34d30266e6ad82f903ba145fa269c0623d1eb041..5b0d6b02d59b947a9b1314b357abaeb76aaa9813 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh @@ -60,7 +60,11 @@ class LowReKEpsilonProblem : public RANSProblem using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; public: - //! The constructor sets the gravity, if desired by the user. + /* + * \brief The constructor + * \param fvGridGeometry The finite volume grid geometry + * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") + */ LowReKEpsilonProblem(std::shared_ptr fvGridGeometry, const std::string& paramGroup = "") : ParentType(fvGridGeometry, paramGroup) { diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh index 9dcb8dce71e63c67fe3973b5d89bac1fd0f60f84..e835707b63f9be4a28bd445413755999273c77bc 100644 --- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh +++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh @@ -181,6 +181,15 @@ struct IndexHelper { return FFFS::compIdx(coupledCompdIdx); } }; +//! forward declare +template +class DarcysLawImplementation; + +//! forward declare +template +class ForchheimersLawImplementation; + + template class StokesDarcyCouplingDataImplementation; @@ -218,9 +227,16 @@ class StokesDarcyCouplingDataImplementationBase template using FluidSystem = typename GET_PROP_TYPE(SubDomainTypeTag, FluidSystem); template using ModelTraits = typename GET_PROP_TYPE(SubDomainTypeTag, ModelTraits); + static constexpr auto stokesIdx = CouplingManager::stokesIdx; static constexpr auto darcyIdx = CouplingManager::darcyIdx; + using AdvectionType = typename GET_PROP_TYPE(SubDomainTypeTag, AdvectionType); + using DarcysLaw = DarcysLawImplementation, GET_PROP_TYPE(SubDomainTypeTag, FVGridGeometry)::discMethod>; + using ForchheimersLaw = ForchheimersLawImplementation, GET_PROP_TYPE(SubDomainTypeTag, FVGridGeometry)::discMethod>; + static constexpr bool darcyUsed = std::is_same::value; + static constexpr bool forchheimerUsed = std::is_same::value; + static constexpr bool adapterUsed = ModelTraits::numPhases() > 1; using IndexHelper = Dumux::IndexHelper, adapterUsed>; @@ -289,20 +305,10 @@ public: const Scalar darcyPressure = stokesContext.volVars.pressure(darcyPhaseIdx); if(numPhasesDarcy > 1) - { momentumFlux = darcyPressure; - } else // use pressure reconstruction for single phase models - { - // v = -K/mu * (gradP + rho*g) - const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf(); - const Scalar mu = stokesContext.volVars.viscosity(darcyPhaseIdx); - const Scalar rho = stokesContext.volVars.density(darcyPhaseIdx); - const Scalar distance = (stokesContext.element.geometry().center() - scvf.center()).two_norm(); - const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()]; - const Scalar interfacePressure = ((scvf.directionSign() * velocity * (mu/darcyPermeability(scvf))) + rho * g) * distance + darcyPressure; - momentumFlux = interfacePressure; - } + momentumFlux = pressureAtInterface_(scvf, stokesElemFaceVars, stokesContext); + // TODO: generalize for permeability tensors // normalize pressure if(GET_PROP_VALUE(SubDomainTypeTag, NormalizePressure)) @@ -424,6 +430,64 @@ protected: return volVars.effectiveThermalConductivity(); } + /*! + * \brief Returns the pressure at the interface using Darcy's law for reconstruction + */ + template = 0> + Scalar pressureAtInterface_(const SubControlVolumeFace& scvf, + const ElementFaceVariables& elemFaceVars, + const CouplingContext& context) const + { + const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); + const Scalar cellCenterPressure = context.volVars.pressure(darcyPhaseIdx); + + // v = -K/mu * (gradP + rho*g) + const Scalar velocity = elemFaceVars[scvf].velocitySelf(); + const Scalar mu = context.volVars.viscosity(darcyPhaseIdx); + const Scalar rho = context.volVars.density(darcyPhaseIdx); + const Scalar distance = (context.element.geometry().center() - scvf.center()).two_norm(); + const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()]; + const Scalar interfacePressure = ((scvf.directionSign() * velocity * (mu/darcyPermeability(scvf))) + rho * g) * distance + cellCenterPressure; + return interfacePressure; + } + + /*! + * \brief Returns the pressure at the interface using Forchheimers's law for reconstruction + */ + template = 0> + Scalar pressureAtInterface_(const SubControlVolumeFace& scvf, + const ElementFaceVariables& elemFaceVars, + const CouplingContext& context) const + { + const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx); + const Scalar cellCenterPressure = context.volVars.pressure(darcyPhaseIdx); + using std::abs; + + // v + cF * sqrt(K) * rho/mu * v * abs(v) + K/mu grad(p + rho z) + const Scalar velocity = elemFaceVars[scvf].velocitySelf(); + const Scalar mu = context.volVars.viscosity(darcyPhaseIdx); + const Scalar rho = context.volVars.density(darcyPhaseIdx); + const Scalar distance = (context.element.geometry().center() - scvf.center()).two_norm(); + const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()]; + + // get the Forchheimer coefficiencient + Scalar cF = 0.0; + for (const auto& darcyScvf : scvfs(context.fvGeometry)) + { + if (darcyScvf.index() == context.darcyScvfIdx) + cF = couplingManager_.problem(darcyIdx).spatialParams().forchCoeff(darcyScvf); + } + + const Scalar interfacePressure = ((scvf.directionSign() * velocity * (mu/darcyPermeability(scvf))) + + (scvf.directionSign() * velocity * abs(velocity) * rho * 1.0/sqrt(darcyPermeability(scvf)) * cF) + + rho * g) * distance + cellCenterPressure; + return interfacePressure; + } + + + private: const CouplingManager& couplingManager_; diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt index addd460feb4d2b502e6a92a0630ad445b11020d9..4d35f107a2aa48051852827bc87ad0e2d159fbe9 100644 --- a/test/multidomain/boundary/CMakeLists.txt +++ b/test/multidomain/boundary/CMakeLists.txt @@ -1,2 +1,4 @@ add_subdirectory(darcydarcy) +add_subdirectory(ransdarcy) +add_subdirectory(ransforchheimer) add_subdirectory(stokesdarcy) diff --git a/test/multidomain/boundary/ransdarcy/CMakeLists.txt b/test/multidomain/boundary/ransdarcy/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..8b84477d3e56133ec980c504db21bbb74df9596d --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/CMakeLists.txt @@ -0,0 +1,67 @@ +add_input_file_links() +dune_symlink_to_source_files(FILES temp.sh) + +dune_add_test(NAME test_kepsilon1p2cnidarcy2p2cni + SOURCES test_rans1p2cnidarcy2p2cni.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_kepsilon1p2cnidarcy2p2cni_kepsilon-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_kepsilon1p2cnidarcy2p2cni_kepsilon-00040.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_kepsilon1p2cnidarcy2p2cni_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_kepsilon1p2cnidarcy2p2cni_darcy-00040.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_kepsilon1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input + -Problem.Name test_kepsilon1p2cnidarcy2p2cni") +target_compile_definitions(test_kepsilon1p2cnidarcy2p2cni PUBLIC "KEPSILON=1") + +dune_add_test(NAME test_komega1p2cnidarcy2p2cni + SOURCES test_rans1p2cnidarcy2p2cni.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_komega1p2cnidarcy2p2cni_komega-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_komega1p2cnidarcy2p2cni_komega-00040.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_komega1p2cnidarcy2p2cni_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_komega1p2cnidarcy2p2cni_darcy-00040.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_komega1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input + -Problem.Name test_komega1p2cnidarcy2p2cni") +target_compile_definitions(test_komega1p2cnidarcy2p2cni PUBLIC "KOMEGA=1") + +dune_add_test(NAME test_lowrekepsilon1p2cnidarcy2p2cni + SOURCES test_rans1p2cnidarcy2p2cni.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_lowrekepsilon1p2cnidarcy2p2cni_lowrekepsilon-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_lowrekepsilon1p2cnidarcy2p2cni_lowrekepsilon-00040.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_lowrekepsilon1p2cnidarcy2p2cni_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_lowrekepsilon1p2cnidarcy2p2cni_darcy-00040.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_lowrekepsilon1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input + -Problem.Name test_lowrekepsilon1p2cnidarcy2p2cni") +target_compile_definitions(test_lowrekepsilon1p2cnidarcy2p2cni PUBLIC "LOWREKEPSILON=1") + +dune_add_test(NAME test_oneeq1p2cnidarcy2p2cni + SOURCES test_rans1p2cnidarcy2p2cni.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_oneeq1p2cnidarcy2p2cni_oneeq-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_oneeq1p2cnidarcy2p2cni_oneeq-00040.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_oneeq1p2cnidarcy2p2cni_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_oneeq1p2cnidarcy2p2cni_darcy-00040.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_oneeq1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input + -Problem.Name test_oneeq1p2cnidarcy2p2cni") +target_compile_definitions(test_oneeq1p2cnidarcy2p2cni PUBLIC "ONEEQ=1") + +dune_add_test(NAME test_zeroeq1p2cnidarcy2p2cni + SOURCES test_rans1p2cnidarcy2p2cni.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_zeroeq1p2cnidarcy2p2cni_zeroeq-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_zeroeq1p2cnidarcy2p2cni_zeroeq-00040.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_zeroeq1p2cnidarcy2p2cni_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_zeroeq1p2cnidarcy2p2cni_darcy-00040.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_zeroeq1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input + -Problem.Name test_zeroeq1p2cnidarcy2p2cni") +target_compile_definitions(test_zeroeq1p2cnidarcy2p2cni PUBLIC "ZEROEQ=1") diff --git a/test/multidomain/boundary/ransdarcy/darcyproblem.hh b/test/multidomain/boundary/ransdarcy/darcyproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..97b3777f3de2d862be28497a32f644d511ba55df --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/darcyproblem.hh @@ -0,0 +1,455 @@ +// -*- 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 + * + * \brief The porous medium sub problem + */ +#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH +#define DUMUX_DARCY2P2C_SUBPROBLEM_HH + +#include + +#include +#include +#include +#include + +#include + +#include "spatialparams.hh" + +namespace Dumux +{ +template +class DarcySubProblem; + +namespace Properties +{ +NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI)); + +// Set the problem property +SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem); + +// the fluid system +SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir); + +//! Set the default formulation to pw-Sn: This can be over written in the problem. +SET_PROP(DarcyTwoPTwoCTypeTag, Formulation) +{ static constexpr auto value = TwoPFormulation::p1s0; }; + +// The gas component balance (air) is replaced by the total mass balance +SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3); + +// Set the grid type +SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates >); + +SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true); + +SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams); +} + +/*! + * \brief The porous medium sub problem + */ +template +class DarcySubProblem : public PorousMediumFlowProblem +{ + using ParentType = PorousMediumFlowProblem; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; + + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + + // copy some indices for convenience + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + enum { + // primary variable indices + conti0EqIdx = Indices::conti0EqIdx, + contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx, + contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx, + pressureIdx = Indices::pressureIdx, + switchIdx = Indices::switchIdx + }; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager); + using TimeLoopPtr = std::shared_ptr>; + + using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType; + +public: + DarcySubProblem(std::shared_ptr fvGridGeometry, + std::shared_ptr couplingManager) + : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) + { + pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); + initialSw_ = getParamFromGroup(this->paramGroup(), "Problem.Saturation"); + temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature"); + initialPhasePresence_ = getParamFromGroup(this->paramGroup(), "Problem.InitPhasePresence"); + + diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{}, + getParamFromGroup(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg")); + + // initialize output file + plotFluxes_ = getParamFromGroup(this->paramGroup(), "Problem.PlotFluxes", false); + plotStorage_ = getParamFromGroup(this->paramGroup(), "Problem.PlotStorage", false); + storageFileName_ = "storage_" + getParam("Problem.Name") + "_" + this->name() + ".csv"; + storageFile_.open(storageFileName_); + storageFile_ << "#Time[s]" << ";" + << "WaterMass[kg]" << ";" + << "WaterMassLoss[kg]" << ";" + << "EvaporationRate[mm/d]" + << std::endl; + } + + /*! + * \name Simulation steering + */ + // \{ + + /*! + * \brief Initialize the problem. + */ + template + void init(const SolutionVector& curSol, + const GridVariables& gridVariables) + { + initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables); + lastWaterMass_ = initialWaterContent_; + } + + template + void postTimeStep(const SolutionVector& curSol, + const GridVariables& gridVariables) + + { + evaluateWaterMassStorageTerm(curSol, gridVariables); + evaluateInterfaceFluxes(curSol, gridVariables); + + gnuplotStorage_.resetPlot(); + if (timeLoop_->time() < getParam("TimeLoop.TEnd")) + gnuplotStorage_.setCreateImage(false); + gnuplotStorage_.setDatafileSeparator(';'); + gnuplotStorage_.setXlabel("time [d]"); + gnuplotStorage_.setXRange(0.0, getParam("TimeLoop.TEnd") / 86400); + gnuplotStorage_.setYlabel("evaporation rate [mm/d]"); + gnuplotStorage_.setOption("set yrange [0.0:15.0]"); + gnuplotStorage_.setOption("set y2label 'cumulative mass loss [kg]'"); + gnuplotStorage_.setOption("set y2range [0.0:15.0]"); + gnuplotStorage_.setOption("set ytics nomirror"); + gnuplotStorage_.setOption("set y2tics"); + + gnuplotStorage_.addFileToPlot(storageFileName_, "using ($1/86400):4 with lines title 'evaporation rate'"); + gnuplotStorage_.addFileToPlot(storageFileName_, "using ($1/86400):3 axes x1y2 with lines title 'cumulative mass loss'"); + gnuplotStorage_.addFileToPlot("storage_test_kepsilon1p2cnidarcy2p2cni_darcy.csv", "using ($1/86400):4 with lines title 'kepsilon'"); + gnuplotStorage_.addFileToPlot("storage_test_komega1p2cnidarcy2p2cni_darcy.csv", "using ($1/86400):4 with lines title 'komega'"); + gnuplotStorage_.addFileToPlot("storage_test_lowrekepsilon1p2cnidarcy2p2cni_darcy.csv", "using ($1/86400):4 with lines title 'lowrekepsilon'"); + gnuplotStorage_.addFileToPlot("storage_test_oneeq1p2cnidarcy2p2cni_darcy.csv", "using ($1/86400):4 with lines title 'oneeq'"); + gnuplotStorage_.addFileToPlot("storage_test_zeroeq1p2cnidarcy2p2cni_darcy.csv", "using ($1/86400):4 with lines title 'zeroeq'"); + if (plotStorage_) + gnuplotStorage_.plot("evapRate_"); + } + + template + Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol, + const GridVariables& gridVariables) + + { + // compute the mass in the entire domain + Scalar waterMass = 0.0; + + for (const auto& element : elements(this->fvGridGeometry().gridView())) + { + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(element); + + auto elemVolVars = localView(gridVariables.curGridVolVars()); + elemVolVars.bindElement(element, fvGeometry, curSol); + + for (auto&& scv : scvs(fvGeometry)) + { + const auto& volVars = elemVolVars[scv]; + for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) + { + // insert calculation of the water mass here + waterMass += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx) * volVars.density(phaseIdx) + * volVars.saturation(phaseIdx) * volVars.porosity() + * scv.volume() * volVars.extrusionFactor(); + } + } + } + + Scalar cumMassLoss = initialWaterContent_ - waterMass; + Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400 + / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]) + / timeLoop_->timeStepSize(); + lastWaterMass_ = waterMass; + + storageFile_ << timeLoop_->time() << ";" + << waterMass << ";" + << cumMassLoss << ";" + << evaporationRate + << std::endl; + + std::cout << "Mass of water is: " << waterMass << std::endl; + std::cout << "Evaporation rate is: " << evaporationRate << std::endl; + + return waterMass; + } + + template + void evaluateInterfaceFluxes(const SolutionVector& curSol, + const GridVariables& gridVariables) + + { + using std::max; + using std::min; + std::vector x; + std::vector y; + static Scalar maxFlux = -9e9; + static Scalar minFlux = 9e9; + + for (const auto& element : elements(this->fvGridGeometry().gridView())) + { + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(element); + + auto elemVolVars = localView(gridVariables.curGridVolVars()); + elemVolVars.bindElement(element, fvGeometry, curSol); + + for (auto&& scvf : scvfs(fvGeometry)) + { + if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + continue; + + // NOTE: binding the coupling context is necessary + couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element); + const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_); + NumEqVector flux(0.0); + for(int i = 0; i< massFlux.size(); ++i) + flux[i] = massFlux[i]; + + x.push_back(scvf.center()[0]); + y.push_back(flux[contiWEqIdx]); + maxFlux = max(maxFlux,flux[contiWEqIdx]); + minFlux = min(minFlux,flux[contiWEqIdx]); + } + } + + gnuplotInterfaceFluxes_.resetPlot(); + gnuplotInterfaceFluxes_.setCreateImage(false); + gnuplotInterfaceFluxes_.setXlabel("x-position [m]"); + gnuplotInterfaceFluxes_.setXRange(this->fvGridGeometry().bBoxMin()[0], this->fvGridGeometry().bBoxMax()[0]); + gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]"); + gnuplotInterfaceFluxes_.setYRange(minFlux, maxFlux); + gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 "); + std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) + + "_" + getParam("Problem.Name") + "_" + this->name() + ".csv"; + gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'"); + if (plotFluxes_) + gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex())); + } + + /*! + * \brief Return the temperature within the domain in [K]. + */ + Scalar temperature() const + { return temperature_; } + // \} + + /*! + * \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 scvf The boundary sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + values.setAllNeumann(); + + if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + values.setAllCouplingNeumann(); + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Dirichlet control volume. + * + * \param element The element for which the Dirichlet boundary condition is set + * \param scvf The boundary subcontrolvolumeface + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + PrimaryVariables values(0.0); + values = initialAtPos(scvf.center()); + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann + * control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param scvf The boundary sub control volume face + * + */ + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + { + const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_); + + for(int i = 0; i< massFlux.size(); ++i) + values[i] = massFlux[i]; + + values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_); + } + + 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 + * + * For this method, the \a values variable stores the rate mass + * of a component is generated or annihilate per volume + * unit. Positive values mean that mass is created, negative ones + * mean that it vanishes. + */ + NumEqVector source(const Element &element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const + { return NumEqVector(0.0); } + + // \} + + /*! + * \brief Evaluate the initial value for a control volume. + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + PrimaryVariables values(0.0); + values.setState(initialPhasePresence_); + + values[pressureIdx] = pressure_ + 1. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]); + values[switchIdx] = initialSw_; + values[Indices::temperatureIdx] = temperature_; + + return values; + } + + // \} + + /*! + * \brief Set the coupling manager + */ + void setCouplingManager(std::shared_ptr cm) + { couplingManager_ = cm; } + + /*! + * \brief Get the coupling manager + */ + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } + + Scalar pressure_; + Scalar initialSw_; + Scalar temperature_; + int initialPhasePresence_; + + TimeLoopPtr timeLoop_; + + Scalar eps_; + + std::shared_ptr couplingManager_; + DiffusionCoefficientAveragingType diffCoeffAvgType_; + + std::string storageFileName_; + std::ofstream storageFile_; + bool plotFluxes_; + bool plotStorage_; + Scalar initialWaterContent_ = 0.0; + Scalar lastWaterMass_ = 0.0; + Dumux::GnuplotInterface gnuplotInterfaceFluxes_; + Dumux::GnuplotInterface gnuplotStorage_; +}; +} //end namespace + +#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/ransdarcy/ransproblem.hh b/test/multidomain/boundary/ransdarcy/ransproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..0ce62b22bab704e6e943eed4d6d9ed262558bc90 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/ransproblem.hh @@ -0,0 +1,528 @@ +// -*- 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 + * \brief The free-flow sub problem + */ +#ifndef DUMUX_RANS1P2C_SUBPROBLEM_HH +#define DUMUX_RANS1P2C_SUBPROBLEM_HH + +#include + +#include +#include +#include +#include + +#if ONEEQ +#include +#include +#elif KEPSILON +#include +#include +#elif KOMEGA +#include +#include +#elif LOWREKEPSILON +#include +#include +#else +#include +#include +#endif + +namespace Dumux +{ +template +class FreeFlowSubProblem; + +namespace Properties +{ +#if ONEEQ +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +#elif KEPSILON +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KEpsilonNCNI)); +#elif KOMEGA +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmegaNCNI)); +#elif LOWREKEPSILON +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, LowReKEpsilonNCNI)); +#else +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI)); +#endif + +// Set the grid type +SET_TYPE_PROP(RANSTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates >); + +// The fluid system +SET_PROP(RANSTypeTag, FluidSystem) +{ + using H2OAir = FluidSystems::H2OAir; + static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase + using type = FluidSystems::OnePAdapter; +}; + +SET_INT_PROP(RANSTypeTag, ReplaceCompEqIdx, 3); + +// Use formulation based on mass fractions +SET_BOOL_PROP(RANSTypeTag, UseMoles, true); + +// Set the problem property +SET_TYPE_PROP(RANSTypeTag, Problem, Dumux::FreeFlowSubProblem ); + +SET_BOOL_PROP(RANSTypeTag, EnableFVGridGeometryCache, true); +SET_BOOL_PROP(RANSTypeTag, EnableGridFluxVariablesCache, true); +SET_BOOL_PROP(RANSTypeTag, EnableGridVolumeVariablesCache, true); + +// SET_BOOL_PROP(RANSTypeTag, EnableInertiaTerms, true); +} + +/*! + * \brief The free-flow sub problem + */ +template +#if ONEEQ +class FreeFlowSubProblem : public OneEqProblem +{ + using ParentType = OneEqProblem; +#elif KEPSILON +class FreeFlowSubProblem : public KEpsilonProblem +{ + using ParentType = KEpsilonProblem; +#elif KOMEGA +class FreeFlowSubProblem : public KOmegaProblem +{ + using ParentType = KOmegaProblem; +#elif LOWREKEPSILON +class FreeFlowSubProblem : public LowReKEpsilonProblem +{ + using ParentType = LowReKEpsilonProblem; +#else +class FreeFlowSubProblem : public ZeroEqProblem +{ + using ParentType = ZeroEqProblem; +#endif + + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; + using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView; + using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState); + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + + using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager); + using TimeLoopPtr = std::shared_ptr>; + + using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType; + + static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles(); + +public: + FreeFlowSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) + : ParentType(fvGridGeometry, "RANS"), eps_(1e-6), couplingManager_(couplingManager) + { + inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.InletVelocity"); + inletPressure_ = getParamFromGroup(this->paramGroup(), "Problem.InletPressure"); + auto inletMassFrac = getParamFromGroup(this->paramGroup(), "Problem.InletMassFrac"); + inletTemperature_ = getParamFromGroup(this->paramGroup(), "Problem.InletTemperature"); + + diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{}, + getParamFromGroup(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg")); + + Dumux::TurbulenceProperties turbulenceProperties; + FluidState fluidState; + fluidState.setPressure(0, inletPressure_); + fluidState.setMassFraction(0, 1, inletMassFrac); + fluidState.setMassFraction(0, 0, 1.0 - inletMassFrac); + fluidState.setTemperature(inletTemperature_); + inletMoleFrac_ = fluidState.moleFraction(0, 1); + Scalar density = FluidSystem::density(fluidState, 0); + Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density; + Scalar charLength = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; + // ideally the viscosityTilde parameter as inflow for the Spalart-Allmaras model should be zero + viscosityTilde_ = getParam("Problem.InletViscosityTilde", + 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, charLength, kinematicViscosity, true)); + turbulentKineticEnergy_ = getParam("Problem.InletTurbulentKineticEnergy", + turbulenceProperties.turbulentKineticEnergy(inletVelocity_, charLength, kinematicViscosity, true)); +#if KOMEGA + dissipation_ = getParam("Problem.InletDissipationRate", + turbulenceProperties.dissipationRate(inletVelocity_, charLength, kinematicViscosity, true)); +#else + dissipation_ = getParam("Problem.InletDissipation", + turbulenceProperties.dissipation(inletVelocity_, charLength, kinematicViscosity, true)); +#endif + std::cout << std::endl; + + } + /*! + * \name Boundary Locations + */ + // \{ + + bool isOnWall(const SubControlVolumeFace& scvf) const + { + GlobalPosition globalPos = scvf.ipGlobal(); + return isOnWallAtPos(globalPos); + } + + bool isOnWallAtPos(const GlobalPosition& globalPos) const + { + return onLowerBoundary_(globalPos); + } + // \} + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Return the temperature within the domain in [K]. + */ + Scalar temperature() const + { return inletTemperature_; } + + /*! + * \brief Return the sources within the domain. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition &globalPos) const + { return NumEqVector(0.0); } + + // \} + /*! + * \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 bTypes; + + const auto& globalPos = scvf.center(); + +#if LOWREKEPSILON || KEPSILON || KOMEGA + bTypes.setDirichlet(Indices::turbulentKineticEnergyIdx); + bTypes.setDirichlet(Indices::dissipationIdx); +#endif + +#if ONEEQ + bTypes.setDirichlet(Indices::viscosityTildeIdx); +#endif + + if (onLeftBoundary_(globalPos)) + { + bTypes.setDirichlet(Indices::velocityXIdx); + bTypes.setDirichlet(Indices::velocityYIdx); + bTypes.setDirichlet(Indices::conti0EqIdx + 1); + bTypes.setDirichlet(Indices::energyBalanceIdx); + } + + if (onLowerBoundary_(globalPos)) + { + bTypes.setDirichlet(Indices::velocityXIdx); + bTypes.setDirichlet(Indices::velocityYIdx); + bTypes.setNeumann(Indices::conti0EqIdx); + bTypes.setNeumann(Indices::conti0EqIdx + 1); + bTypes.setNeumann(Indices::energyBalanceIdx); + } + + if (onUpperBoundary_(globalPos)) + { + bTypes.setAllSymmetry(); + } + + if (onRightBoundary_(globalPos)) + { + bTypes.setDirichlet(Indices::pressureIdx); + bTypes.setOutflow(Indices::conti0EqIdx + 1); + bTypes.setOutflow(Indices::energyBalanceIdx); +#if LOWREKEPSILON || KEPSILON || KOMEGA + bTypes.setOutflow(Indices::turbulentKineticEnergyEqIdx); + bTypes.setOutflow(Indices::dissipationEqIdx); +#endif +#if ONEEQ + bTypes.setOutflow(Indices::viscosityTildeIdx); +#endif + } + + if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + bTypes.setCouplingNeumann(Indices::conti0EqIdx); + bTypes.setCouplingNeumann(Indices::conti0EqIdx + 1); + bTypes.setCouplingNeumann(Indices::energyBalanceIdx); + + bTypes.setCouplingNeumann(scvf.directionIndex()); + bTypes.setBJS(1 - scvf.directionIndex()); + } + +#if KOMEGA + // set a fixed dissipation (omega) in one cell + if (isOnWallAtPos(globalPos)) + bTypes.setDirichletCell(Indices::dissipationIdx); +#endif + + return bTypes; + } + + /*! + * \brief Evaluate the boundary conditions for a dirichlet values at the boundary. + * + * \param element The finite element + * \param scvf the sub control volume face + * \note used for cell-centered discretization schemes + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + return initialAtPos(scvf.ipGlobal()); + } + + /*! + * \brief Evaluate the boundary conditions for fixed values at cell centers + * + * \param element The finite element + * \param scv the sub control volume + * \note used for cell-centered discretization schemes + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const + { + const auto globalPos = scv.center(); + PrimaryVariables values(initialAtPos(globalPos)); +#if KOMEGA + using std::pow; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + const auto wallDistance = ParentType::wallDistance_[elementIdx]; + values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx] + / (ParentType::betaOmega() * pow(wallDistance, 2)); +#endif + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFaceVars The element face variables + * \param scvf The boundary sub control volume face + */ + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars, + const SubControlVolumeFace& scvf) const + { + PrimaryVariables values(0.0); + if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf); + + const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_); + values[Indices::conti0EqIdx] = massFlux[0]; + values[Indices::conti0EqIdx + 1] = massFlux[1]; + values[Indices::energyBalanceIdx] = couplingManager().couplingData().energyCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_); + } + return values; + } + + // \} + + /*! + * \brief Set the coupling manager + */ + void setCouplingManager(std::shared_ptr cm) + { couplingManager_ = cm; } + + /*! + * \brief Get the coupling manager + */ + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + bool isOnWall(const GlobalPosition& globalPos) const + { + return (onLowerBoundary_(globalPos)); + } + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + FluidState fluidState; + updateFluidStateForBC_(fluidState, inletTemperature(), inletPressure(), inletMoleFrac()); + + const Scalar density = FluidSystem::density(fluidState, 0); + + PrimaryVariables values(0.0); + values[Indices::pressureIdx] = inletPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]); + values[Indices::conti0EqIdx + 1] = inletMoleFrac(); + values[Indices::velocityXIdx] = inletVelocity(); + values[Indices::temperatureIdx] = inletTemperature(); + + if(onLowerBoundary_(globalPos)) + values[Indices::velocityXIdx] = 0.0; + +#if LOWREKEPSILON || KEPSILON || KOMEGA + values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; + values[Indices::dissipationIdx] = dissipation_; + + if (isOnWallAtPos(globalPos)) + { + values[Indices::turbulentKineticEnergyIdx] = 0.0; + values[Indices::dissipationIdx] = 0.0; + } +#endif + +#if ONEEQ + values[Indices::viscosityTildeIdx] = viscosityTilde_; + if (isOnWallAtPos(globalPos)) + { + values[Indices::viscosityTildeIdx] = 0.0; + } +#endif + + return values; + } + + //! \brief Returns the inlet velocity. + const Scalar inletVelocity() const + { return inletVelocity_ ;} + + //! \brief Returns the inlet pressure. + const Scalar inletPressure() const + { return inletPressure_; } + + //! \brief Returns the inlet mass fraction. + const Scalar inletMoleFrac() const + { return inletMoleFrac_; } + + //! \brief Returns the inlet temperature. + const Scalar inletTemperature() const + { return inletTemperature_; } + + + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar permeability(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center()); + } + + /*! + * \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::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); + } + + // \} + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } + + //! \brief updates the fluid state to obtain required quantities for IC/BC + void updateFluidStateForBC_(FluidState& fluidState, const Scalar temperature, + const Scalar pressure, const Scalar moleFraction) const + { + fluidState.setTemperature(temperature); + fluidState.setPressure(0, pressure); + fluidState.setSaturation(0, 1.0); + fluidState.setMoleFraction(0, 1, moleFraction); + fluidState.setMoleFraction(0, 0, 1.0 - moleFraction); + + typename FluidSystem::ParameterCache paramCache; + paramCache.updatePhase(fluidState, 0); + + const Scalar density = FluidSystem::density(fluidState, paramCache, 0); + fluidState.setDensity(0, density); + + const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0); + fluidState.setMolarDensity(0, molarDensity); + + const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0); + fluidState.setEnthalpy(0, enthalpy); + } + + // the height of the free-flow domain + const Scalar height_() const + { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; } + + Scalar eps_; + + Scalar inletVelocity_; + Scalar inletPressure_; + Scalar inletMoleFrac_; + Scalar inletTemperature_; + Scalar viscosityTilde_; + Scalar turbulentKineticEnergy_; + Scalar dissipation_; + + TimeLoopPtr timeLoop_; + + std::shared_ptr couplingManager_; + + DiffusionCoefficientAveragingType diffCoeffAvgType_; +}; +} //end namespace + +#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/ransdarcy/spatialparams.hh b/test/multidomain/boundary/ransdarcy/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..dd70bbfbd3ca5f6a1d7d899827aa8a8f229d7388 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/spatialparams.hh @@ -0,0 +1,140 @@ +// -*- 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 TwoPTests + * \brief The spatial parameters class for a test problem using the 2p cc model + */ +#ifndef DUMUX_SPATIAL_PARAMS_HH +#define DUMUX_SPATIAL_PARAMS_HH + +#include +#include +#include +#include + +namespace Dumux +{ + +/*! + * \ingroup TwoPModel + * \ingroup ImplicitTestProblems + * + * \brief The spatial parameters class for the test problem using the 2p cc model + */ +template +class TwoPSpatialParams +: public FVSpatialParams> +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using ParentType = FVSpatialParams>; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using EffectiveLaw = RegularizedVanGenuchten; + +public: + using MaterialLaw = EffToAbsLaw; + using MaterialLawParams = typename MaterialLaw::Params; + using PermeabilityType = Scalar; + + TwoPSpatialParams(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) + { + permeability_ = getParam("Darcy.SpatialParams.Permeability"); + porosity_ = getParam("Darcy.SpatialParams.Porosity"); + alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); + + // residual saturations + params_.setSwr(getParam("Darcy.SpatialParams.Swr")); + params_.setSnr(getParam("Darcy.SpatialParams.Snr")); + // parameters for the vanGenuchten law + params_.setVgAlpha(getParam("Darcy.SpatialParams.VgAlpha")); + params_.setVgn(getParam("Darcy.SpatialParams.VgN")); + Scalar threshold = 0.01 * (1.0 - params_.swr() - params_.snr()); + params_.setPcLowSw(getParam("Darcy.SpatialParams.PcLowSw", threshold)); + params_.setPcHighSw(getParam("Darcy.SpatialParams.PcHighSw", 1.0-threshold)); + } + + /*! + * \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 permeability_; } + + /*! \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_; } + + /*! + * \brief Returns the parameter object for the Brooks-Corey material law. + * In this test, we use element-wise distributed material parameters. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the material parameters object + */ + template + const MaterialLawParams& materialLawParams(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { return params_; } + + /*! + * \brief Function for defining which phase is to be considered as the wetting phase. + * + * \return the wetting phase index + * \param globalPos The global position + */ + template + int wettingPhaseAtPos(const GlobalPosition& globalPos) const + { return FluidSystem::phase0Idx; } + +private: + Scalar permeability_; + Scalar porosity_; + Scalar alphaBJ_; + MaterialLawParams params_; + static constexpr Scalar eps_ = 1.0e-7; +}; + +} // end namespace + +#endif diff --git a/test/multidomain/boundary/ransdarcy/temp.sh b/test/multidomain/boundary/ransdarcy/temp.sh new file mode 100644 index 0000000000000000000000000000000000000000..fb359965647347c00ebeabe28c3664a719afd963 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/temp.sh @@ -0,0 +1,5 @@ +./test_zeroeq1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input -Problem.Name test_zeroeq1p2cnidarcy2p2cni +./test_oneeq1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input -Problem.Name test_oneeq1p2cnidarcy2p2cni +./test_kepsilon1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input -Problem.Name test_kepsilon1p2cnidarcy2p2cni +./test_komega1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input -Problem.Name test_komega1p2cnidarcy2p2cni +./test_lowrekepsilon1p2cnidarcy2p2cni test_rans1p2cnidarcy2p2cni.input -Problem.Name test_lowrekepsilon1p2cnidarcy2p2cni diff --git a/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc new file mode 100644 index 0000000000000000000000000000000000000000..5d174304a2fe3e183780deb40c3c0577bd3f9b48 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc @@ -0,0 +1,303 @@ +// -*- 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 + * + * \brief A test problem for the isothermal coupled RANS/Darcy problem (1p2c/2p2c) + */ +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "darcyproblem.hh" +#include "ransproblem.hh" + +namespace Dumux { +namespace Properties { + +SET_PROP(RANSTypeTag, CouplingManager) +{ + using Traits = StaggeredMultiDomainTraits; + using type = Dumux::StokesDarcyCouplingManager; +}; + +SET_PROP(DarcyTwoPTwoCTypeTag, 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 RANSTypeTag = TTAG(RANSTypeTag); + using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag); + + // 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 RANSGridManager = Dumux::GridManager; + RANSGridManager ransGridManager; + ransGridManager.init("RANS"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + const auto& ransGridView = ransGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using RANSFVGridGeometry = typename GET_PROP_TYPE(RANSTypeTag, FVGridGeometry); + auto ransFvGridGeometry = std::make_shared(ransGridView); + ransFvGridGeometry->update(); + using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry); + auto darcyFvGridGeometry = std::make_shared(darcyGridView); + darcyFvGridGeometry->update(); + + using Traits = StaggeredMultiDomainTraits; + + // the coupling manager + using CouplingManager = StokesDarcyCouplingManager; + auto couplingManager = std::make_shared(ransFvGridGeometry, darcyFvGridGeometry); + + // the indices + constexpr auto ransCellCenterIdx = CouplingManager::stokesCellCenterIdx; + constexpr auto ransFaceIdx = CouplingManager::stokesFaceIdx; + constexpr auto darcyIdx = CouplingManager::darcyIdx; + + // initialize the fluidsystem (tabulation) + GET_PROP_TYPE(RANSTypeTag, FluidSystem)::init(); + + // the problem (initial and boundary conditions) + using RANSProblem = typename GET_PROP_TYPE(RANSTypeTag, Problem); + auto ransProblem = std::make_shared(ransFvGridGeometry, couplingManager); + using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem); + auto darcyProblem = std::make_shared(darcyFvGridGeometry, couplingManager); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(RANSTypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("TimeLoop.DtInitial"); + + // check if we are about to restart a previously interrupted simulation + Scalar restartTime = 0; + if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) + restartTime = getParam("TimeLoop.Restart"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // set timeloop for the subproblems, needed for boundary value variations + ransProblem->setTimeLoop(timeLoop); + darcyProblem->setTimeLoop(timeLoop); + + // the solution vector + Traits::SolutionVector sol; + sol[ransCellCenterIdx].resize(ransFvGridGeometry->numCellCenterDofs()); + sol[ransFaceIdx].resize(ransFvGridGeometry->numFaceDofs()); + sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); + + // apply initial solution for instationary problems + // auxiliary free flow solution vector + typename GET_PROP_TYPE(RANSTypeTag, SolutionVector) ransSol; + ransSol[ransCellCenterIdx].resize(sol[ransCellCenterIdx].size()); + ransSol[ransFaceIdx].resize(sol[ransFaceIdx].size()); + ransProblem->applyInitialSolution(ransSol); + ransProblem->updateStaticWallProperties(); + ransProblem->updateDynamicWallProperties(ransSol); + + auto solRANSOld = ransSol; + sol[ransCellCenterIdx] = ransSol[ransCellCenterIdx]; + sol[ransFaceIdx] = ransSol[ransFaceIdx]; + + darcyProblem->applyInitialSolution(sol[darcyIdx]); + auto solDarcyOld = sol[darcyIdx]; + + auto solOld = sol; + + // intialize the coupling manager + couplingManager->init(ransProblem, darcyProblem, sol); + + // intializethe grid variables and the subproblems + using RANSGridVariables = typename GET_PROP_TYPE(RANSTypeTag, GridVariables); + auto ransGridVariables = std::make_shared(ransProblem, ransFvGridGeometry); + ransGridVariables->init(ransSol, solRANSOld); + using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables); + auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry); + darcyGridVariables->init(sol[darcyIdx], solDarcyOld); + darcyProblem->init(sol[darcyIdx], *darcyGridVariables); + + // intialize the vtk output module + const auto ransName = getParam("Problem.Name") + "_" + ransProblem->name(); + const auto darcyName = getParam("Problem.Name") + "_" + darcyProblem->name(); + + StaggeredVtkOutputModule ransVtkWriter(*ransGridVariables, ransSol, ransName); + GET_PROP_TYPE(RANSTypeTag, VtkOutputFields)::init(ransVtkWriter); + ransVtkWriter.write(0.0); + + VtkOutputModule darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName); + using DarcyVelocityOutput = typename GET_PROP_TYPE(DarcyTypeTag, VelocityOutput); + darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); + GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter); + darcyVtkWriter.write(0.0); + + // the assembler with time loop for instationary problem + using Assembler = MultiDomainFVAssembler; + auto assembler = std::make_shared(std::make_tuple(ransProblem, ransProblem, darcyProblem), + std::make_tuple(ransFvGridGeometry->cellCenterFVGridGeometryPtr(), + ransFvGridGeometry->faceFVGridGeometryPtr(), + darcyFvGridGeometry), + std::make_tuple(ransGridVariables->cellCenterGridVariablesPtr(), + ransGridVariables->faceGridVariablesPtr(), + darcyGridVariables), + couplingManager, + timeLoop); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + + // the primary variable switches used by the sub models + using PriVarSwitchTuple = std::tuple; + + // the non-linear solver + using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(solOld); + + // solve the non-linear system with time step control + nonLinearSolver.solve(sol, *timeLoop); + + // make the new solution the old solution + solOld = sol; + + // update the auxiliary free flow solution vector + ransSol[ransCellCenterIdx] = sol[ransCellCenterIdx]; + ransSol[ransFaceIdx] = sol[ransFaceIdx]; + ransProblem->updateDynamicWallProperties(ransSol); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + // post time step treatment of Darcy problem + darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables); + + // advance grid variables to the next time step + ransGridVariables->advanceTimeStep(); + darcyGridVariables->advanceTimeStep(); + + // write vtk output + ransVtkWriter.write(timeLoop->time()); + darcyVtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton solver + timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); + + } while (!timeLoop->finished()); + + timeLoop->finalize(ransGridView.comm()); + timeLoop->finalize(darcyGridView.comm()); + + //////////////////////////////////////////////////////////// + // 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/ransdarcy/test_rans1p2cnidarcy2p2cni.input b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input new file mode 100644 index 0000000000000000000000000000000000000000..e955974c5d9ba1cdb962f12de316b636bda18267 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input @@ -0,0 +1,72 @@ +[TimeLoop] +DtInitial = 1e-1 # [s] +MaxTimeStepSize = 43200 # [s] (12 hours) +TEnd = 864000 # [s] (6 days) + +[RANS.Grid] +Positions0 = -1.0 0.0 0.5 +Positions1 = 0.25 0.5 +Grading0 = 1.0 1.0 +Grading1 = 1.57539 +Cells0 = 4 16 +Cells1 = 16 +Verbosity = true + +[Darcy.Grid] +Positions0 = 0.0 0.5 +Positions1 = 0.0 0.25 +Cells0 = 16 +Cells1 = 16 +Grading0 = 1.0 +Grading1 = -1.57539 +Verbosity = true + +[RANS.Problem] +Name = rans +InletVelocity = 3.5 # [m/s] +InletPressure = 1e5 # [Pa] +InletMassFrac = 0.008 # [-] +InletTemperature = 298.15 # [K] + +[Darcy.Problem] +Name = darcy +Pressure = 1e5 +Saturation = 0.8 # initial Sw +Temperature = 298.15 # [K] +InitPhasePresence = 3 # bothPhases + +[Darcy.SpatialParams] +Porosity = 0.41 +Permeability = 2.65e-10 +AlphaBeaversJoseph = 1.0 +Swr = 0.005 +Snr = 0.01 +VgAlpha = 6.37e-4 +VgN = 8.0 + +[Problem] +Name = test_rans1p2cnidarcy2p2cni +EnableGravity = true +InterfaceDiffusionCoefficientAvg = Harmonic # FreeFlowOnly Harmonic +PlotFluxes = false +PlotStorage = true + +[RANS.RANS] +EddyViscosityModel = "vanDriest" +UseStoredEddyViscosity = false + +[RANS.KEpsilon] +YPlusThreshold = 50. # should be large (30-60) for fine grids + +[Vtk] +AddVelocity = true +WriteFaceData = false + +[Newton] +TargetSteps = 5 +MaxSteps = 12 +MaxRelativeShift = 1e-5 + +[Assembly] +NumericDifferenceMethod = 0 +NumericDifference.BaseEpsilon = 1e-8 diff --git a/test/multidomain/boundary/ransforchheimer/CMakeLists.txt b/test/multidomain/boundary/ransforchheimer/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d204e1e99322dfa410bcba2cdc382c699801cdae --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/CMakeLists.txt @@ -0,0 +1,7 @@ +dune_add_test(NAME test_rans1pdarcy1pforchheimer + SOURCES test_rans1pdarcy1pforchheimer.cc + COMPILE_DEFINITIONS FORCHHEIMER=1 + CMAKE_GUARD HAVE_UMFPACK + CMD_ARGS test_rans1pdarcy1pforchheimer.input) + +dune_symlink_to_source_files(FILES "test_rans1pdarcy1pforchheimer.input") diff --git a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..dc549bba86063a8fe596b0781df8456b34601b9b --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh @@ -0,0 +1,315 @@ +// -*- 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 + * + * \brief The porous medium sub problem + */ +#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH +#define DUMUX_DARCY2P2C_SUBPROBLEM_HH + +#include + +#include +#include +#include + +#include +#include + +#if FORCHHEIMER +#include +#endif + +#include "spatialparams.hh" + +namespace Dumux +{ +template +class DarcySubProblem; + +namespace Properties +{ +NEW_TYPE_TAG(DarcyOnePTypeTag, INHERITS_FROM(CCTpfaModel, OneP)); + +// Set the problem property +SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem); + +// the fluid system +SET_PROP(DarcyOnePTypeTag, FluidSystem) +{ +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); +public: + using type = FluidSystems::OnePGas >; +}; + +// Set the grid type +SET_TYPE_PROP(DarcyOnePTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates >); + +SET_PROP(DarcyOnePTypeTag, SpatialParams) +{ + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using type = OnePSpatialParams; +}; + +#ifdef FORCHHEIMER +SET_TYPE_PROP(DarcyOnePTypeTag, AdvectionType, ForchheimersLaw); +#endif +} + +/*! + * \brief The porous medium sub problem + */ + +template +class DarcySubProblem : public PorousMediumFlowProblem +{ + using ParentType = PorousMediumFlowProblem; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState); + + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + + enum { + // index of the primary variable + pressureIdx = Indices::pressureIdx + }; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using TimeLoopPtr = std::shared_ptr>; + using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager); + +public: + DarcySubProblem(std::shared_ptr fvGridGeometry, + std::shared_ptr couplingManager) + : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) + { + pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); + temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature"); + baseClosed_ = getParamFromGroup(this->paramGroup(),"Problem.BaseClosed"); + bottomVelocity_ = getParamFromGroup(this->paramGroup(),"Problem.BottomVelocity"); + + FluidState fluidState; + fluidState.setPressure(0, pressure_); + fluidState.setTemperature(temperature_); + density_ = FluidSystem::density(fluidState, 0); + } + + /*! + * \name Simulation steering + */ + // \{ + + /*! + * \brief Returns true if a restart file should be written to + * disk. + */ + bool shouldWriteRestartFile() const + { return false; } + + // \} + + /*! + * \name Problem parameters + */ + // \{ + + bool shouldWriteOutput() const // define output + { return true; } + + /*! + * \brief Return the temperature within the domain in [K]. + */ + Scalar temperature() const + { return temperature_; } + + // \} + + /*! + * \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 scvf The boundary sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + values.setAllNeumann(); + + if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + values.setAllCouplingNeumann(); + + // if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + // { + // values.setDirichlet(Indices::pressureIdx); + // } + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Dirichlet control volume. + * + * \param element The element for which the Dirichlet boundary condition is set + * \param scvf The boundary subcontrolvolumeface + * + * For this method, the \a values parameter stores primary variables. + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + PrimaryVariables values(0.0); + values = initial(element); + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry The fvGeometry + * \param elemVolVars The element volume variables + * \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 SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + { + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf); + } + + if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + { + values[Indices::conti0EqIdx] = -bottomVelocity_ * density_; + } + + 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 + { return NumEqVector(0.0); } + + // \} + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param element The element + * + * For this method, the \a priVars parameter stores primary + * variables. + */ + PrimaryVariables initial(const Element &element) const + { + PrimaryVariables priVars(0.0); + priVars[pressureIdx] = pressure_; + return priVars; + } + + + // \} + + /*! + * \brief Set the coupling manager + */ + void setCouplingManager(std::shared_ptr cm) + { couplingManager_ = cm; } + + /*! + * \brief Get the coupling manager + */ + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } + + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } + + Scalar pressure_; + Scalar temperature_; + Scalar eps_; + bool baseClosed_; + Scalar bottomVelocity_; + Scalar density_; + + TimeLoopPtr timeLoop_; + std::shared_ptr couplingManager_; +}; +} //end namespace + +#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/ransforchheimer/ransproblem.hh b/test/multidomain/boundary/ransforchheimer/ransproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..cdb42a6d54620170319d0635247afdc85cbc9cf1 --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/ransproblem.hh @@ -0,0 +1,394 @@ +// -*- 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 + * \brief The free-flow sub problem + */ +#ifndef DUMUX_RANS1P_SUBPROBLEM_HH +#define DUMUX_RANS1P_SUBPROBLEM_HH + +#include + +#include +#include + +#include +#include + +#include +#include + +namespace Dumux { + +template +class FreeFlowSubProblem; + +namespace Properties { + +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmega)); + +// the fluid system +SET_PROP(RANSTypeTag, FluidSystem) +{ +private: + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); +public: + using type = FluidSystems::OnePGas >; +}; + +// Set the grid type +SET_TYPE_PROP(RANSTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates >); + +// Set the problem property +SET_TYPE_PROP(RANSTypeTag, Problem, Dumux::FreeFlowSubProblem ); + +SET_BOOL_PROP(RANSTypeTag, EnableFVGridGeometryCache, true); +SET_BOOL_PROP(RANSTypeTag, EnableGridFluxVariablesCache, true); +SET_BOOL_PROP(RANSTypeTag, EnableGridVolumeVariablesCache, true); +} + +/*! + * \brief The free-flow sub problem + */ +template +class FreeFlowSubProblem : public KOmegaProblem +{ + using ParentType = KOmegaProblem; + + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState); + + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + + using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager); + using TimeLoopPtr = std::shared_ptr>; + +public: + FreeFlowSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) + : ParentType(fvGridGeometry, "RANS"), eps_(1e-6), couplingManager_(couplingManager) + { + inletReynoldsNumber_ = getParamFromGroup(this->paramGroup(), "Problem.InletReynoldsNumber"); + inletPressure_ = getParamFromGroup(this->paramGroup(), "Problem.InletPressure"); + temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature"); + + darcyStart_ = getParamFromGroup(this->paramGroup(), "Grid.DarcyStart"); + darcyEnd_ = getParamFromGroup(this->paramGroup(), "Grid.DarcyEnd"); + smoothingZoneDistance_ = getParamFromGroup(this->paramGroup(), "Problem.SmoothingZoneDistance"); + smoothingSlipVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.SmoothingSlipVelocity"); + + Dumux::TurbulenceProperties turbulenceProperties; + FluidState fluidState; + fluidState.setPressure(0, inletPressure_); + fluidState.setTemperature(temperature_); + Scalar density = FluidSystem::density(fluidState, 0); + Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density; + Scalar charLength = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; + inletVelocity_ = inletReynoldsNumber_ * kinematicViscosity / charLength; + turbulentKineticEnergy_ = getParam("Problem.InletTurbulentKineticEnergy", + turbulenceProperties.turbulentKineticEnergy(inletVelocity_, charLength, kinematicViscosity, true)); + dissipation_ = getParam("Problem.InletDissipationRate", + turbulenceProperties.dissipationRate(inletVelocity_, charLength, kinematicViscosity, true)); + std::cout << std::endl; + } + + /*! + * \name Boundary Locations + */ + // \{ + + bool isOnWall(const SubControlVolumeFace& scvf) const + { + GlobalPosition globalPos = scvf.ipGlobal(); + return isOnWallAtPos(globalPos); + } + + bool isOnCouplingWall(const SubControlVolumeFace& scvf) const + { + return couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf); + } + + bool isOnWallAtPos(const GlobalPosition& globalPos) const + { + return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)); + } + // \} + + bool isSmoothingZoneAtPos(const GlobalPosition& globalPos) const + { + return (onLowerBoundary_(globalPos) && + (((globalPos[0] > ((darcyStart_ - smoothingZoneDistance_) - eps_)) && (globalPos[0] < (darcyStart_ + eps_))) || + ((globalPos[0] < ((darcyEnd_ + smoothingZoneDistance_) + eps_)) && (globalPos[0] > (darcyEnd_ - eps_))))); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Return the temperature within the domain in [K]. + */ + Scalar temperature() const + { return temperature_; } + + /*! + * \brief Return the sources within the domain. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition &globalPos) const + { return NumEqVector(0.0); } + + // \} + /*! + * \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 bTypes; + + const auto& globalPos = scvf.center(); + + bTypes.setDirichlet(Indices::turbulentKineticEnergyIdx); + bTypes.setDirichlet(Indices::dissipationIdx); + + if (onRightBoundary_(globalPos)) + { + bTypes.setDirichlet(Indices::pressureIdx); + bTypes.setOutflow(Indices::turbulentKineticEnergyEqIdx); + bTypes.setOutflow(Indices::dissipationEqIdx); + } + else + { + bTypes.setDirichlet(Indices::velocityXIdx); + bTypes.setDirichlet(Indices::velocityYIdx); + } + + if (isOnCouplingWall(scvf)) + { + bTypes.setCouplingNeumann(Indices::conti0EqIdx); + bTypes.setCouplingNeumann(scvf.directionIndex()); + bTypes.setBJS(1 - scvf.directionIndex()); + } + + // set a fixed dissipation (omega) in one cell + if (isOnWall(scvf)) + bTypes.setDirichletCell(Indices::dissipationIdx); + + return bTypes; + } + + /*! + * \brief Evaluate the boundary conditions for a dirichlet values at the boundary. + * + * \param element The finite element + * \param scvf the sub control volume face + * \note used for cell-centered discretization schemes + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + return initialAtPos(scvf.ipGlobal()); + } + + /*! + * \brief Evaluate the boundary conditions for fixed values at cell centers + * + * \param element The finite element + * \param scv the sub control volume + * \note used for cell-centered discretization schemes + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const + { + const auto globalPos = scv.center(); + PrimaryVariables values(initialAtPos(globalPos)); + using std::pow; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + const auto wallDistance = ParentType::wallDistance_[elementIdx]; + values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx] + / (ParentType::betaOmega() * pow(wallDistance, 2)); + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeomentry 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 + { + PrimaryVariables values(0.0); + if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf); + values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf); + } + return values; + } + + // \} + + /*! + * \brief Set the coupling manager + */ + void setCouplingManager(std::shared_ptr cm) + { couplingManager_ = cm; } + + /*! + * \brief Get the coupling manager + */ + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + PrimaryVariables values(0.0); + + values[Indices::pressureIdx] = inletPressure(); + values[Indices::velocityXIdx] = inletVelocity(); + + if(onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)) + values[Indices::velocityXIdx] = 0.0; + + values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; + values[Indices::dissipationIdx] = dissipation_; + + if (isSmoothingZoneAtPos(globalPos)) + { + values[Indices::velocityXIdx] = smoothingSlipVelocity_; + } + + if (isOnWallAtPos(globalPos)) + { + values[Indices::turbulentKineticEnergyIdx] = 0.0; + } + + return values; + } + + //! \brief Returns the inlet velocity. + const Scalar inletVelocity() const + { return inletVelocity_ ;} + + //! \brief Returns the inlet pressure. + const Scalar inletPressure() const + { return inletPressure_; } + + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar permeability(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center()); + } + + /*! + * \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::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); + } + + // \} + +private: + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } + + // the height of the free-flow domain + const Scalar height_() const + { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; } + + Scalar eps_; + + Scalar inletReynoldsNumber_; + Scalar inletVelocity_; + Scalar inletPressure_; + Scalar temperature_; + Scalar turbulentKineticEnergy_; + Scalar dissipation_; + Scalar darcyEnd_; + Scalar darcyStart_; + Scalar smoothingZoneDistance_; + Scalar smoothingSlipVelocity_; + + TimeLoopPtr timeLoop_; + + std::shared_ptr couplingManager_; + +}; +} //end namespace Dumux + +#endif // DUMUX_RANS1P_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/ransforchheimer/spatialparams.hh b/test/multidomain/boundary/ransforchheimer/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..af4ba63d06363727cab443b2eef06d5624852a9a --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/spatialparams.hh @@ -0,0 +1,111 @@ +// -*- 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 + * \ingroup ImplicitTestProblems + * + * \brief The spatial parameters class for the test problem using the + * 1p cc model + */ +template +class OnePSpatialParams +: public FVSpatialParamsOneP> +{ + using GridView = typename FVGridGeometry::GridView; + using ParentType = FVSpatialParamsOneP>; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + // export permeability type + using PermeabilityType = Scalar; + + OnePSpatialParams(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) + { + permeability_ = getParam("Darcy.SpatialParams.Permeability"); + porosity_ = getParam("Darcy.SpatialParams.Porosity"); + alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); + fittedAlpha_ = getParam("Darcy.SpatialParams.FittedAlpha"); + baseClosed_ = getParam("Darcy.Problem.BaseClosed"); + } + + /*! + * \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 permeability_; } + + /*! \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 + { + using std::pow; + Scalar x = globalPos[0]; + if(fittedAlpha_ && baseClosed_) + return 39.0*pow(x,4) - 79.0*pow(x,3) + 54.0*pow(x,2) - 15.0*x + 4.0; + else if (fittedAlpha_ && !baseClosed_) + return 44.4*pow(x,4) - 71.3*pow(x,3) + 20.5*pow(x,2) + 7.9*x + 3.0; + else if (!fittedAlpha_ && baseClosed_) + return 2.9; + else if (!fittedAlpha_ && !baseClosed_) + return 4.4; + else + return alphaBJ_; + } + +private: + Scalar permeability_; + Scalar porosity_; + Scalar alphaBJ_; + bool fittedAlpha_; + bool baseClosed_; +}; + +} // end namespace + +#endif diff --git a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.cc b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.cc new file mode 100644 index 0000000000000000000000000000000000000000..a188b01c66e959b5d01fe6d49186c6ba15c8e565 --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.cc @@ -0,0 +1,299 @@ +// -*- 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 + * + * \brief A test problem for the coupled RANS/Darcy problem (1p/1pForcheimer) + */ +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "darcyproblem.hh" +#include "ransproblem.hh" + +namespace Dumux { +namespace Properties { + +SET_PROP(RANSTypeTag, CouplingManager) +{ + using Traits = StaggeredMultiDomainTraits; + using type = Dumux::StokesDarcyCouplingManager; +}; + +SET_PROP(DarcyOnePTypeTag, 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 RANSTypeTag = TTAG(RANSTypeTag); + using DarcyTypeTag = TTAG(DarcyOnePTypeTag); + + // 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 RANSGridManager = Dumux::GridManager; + RANSGridManager ransGridManager; + ransGridManager.init("RANS"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + const auto& ransGridView = ransGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using RANSFVGridGeometry = typename GET_PROP_TYPE(RANSTypeTag, FVGridGeometry); + auto ransFvGridGeometry = std::make_shared(ransGridView); + ransFvGridGeometry->update(); + using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry); + auto darcyFvGridGeometry = std::make_shared(darcyGridView); + darcyFvGridGeometry->update(); + + using Traits = StaggeredMultiDomainTraits; + + // the coupling manager + using CouplingManager = StokesDarcyCouplingManager; + auto couplingManager = std::make_shared(ransFvGridGeometry, darcyFvGridGeometry); + + // the indices + constexpr auto ransCellCenterIdx = CouplingManager::stokesCellCenterIdx; + constexpr auto ransFaceIdx = CouplingManager::stokesFaceIdx; + constexpr auto darcyIdx = CouplingManager::darcyIdx; + + // initialize the fluidsystem (tabulation) + GET_PROP_TYPE(RANSTypeTag, FluidSystem)::init(); + + // the problem (initial and boundary conditions) + using RANSProblem = typename GET_PROP_TYPE(RANSTypeTag, Problem); + auto ransProblem = std::make_shared(ransFvGridGeometry, couplingManager); + using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem); + auto darcyProblem = std::make_shared(darcyFvGridGeometry, couplingManager); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(RANSTypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("TimeLoop.DtInitial"); + + // check if we are about to restart a previously interrupted simulation + Scalar restartTime = 0; + if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) + restartTime = getParam("TimeLoop.Restart"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // set timeloop for the subproblems, needed for boundary value variations + ransProblem->setTimeLoop(timeLoop); + darcyProblem->setTimeLoop(timeLoop); + + // the solution vector + Traits::SolutionVector sol; + sol[ransCellCenterIdx].resize(ransFvGridGeometry->numCellCenterDofs()); + sol[ransFaceIdx].resize(ransFvGridGeometry->numFaceDofs()); + sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); + + // apply initial solution for instationary problems + // auxiliary free flow solution vector + typename GET_PROP_TYPE(RANSTypeTag, SolutionVector) ransSol; + ransSol[ransCellCenterIdx].resize(sol[ransCellCenterIdx].size()); + ransSol[ransFaceIdx].resize(sol[ransFaceIdx].size()); + ransProblem->applyInitialSolution(ransSol); + ransProblem->updateStaticWallProperties(); + ransProblem->updateDynamicWallProperties(ransSol); + + auto solRANSOld = ransSol; + sol[ransCellCenterIdx] = ransSol[ransCellCenterIdx]; + sol[ransFaceIdx] = ransSol[ransFaceIdx]; + + darcyProblem->applyInitialSolution(sol[darcyIdx]); + auto solDarcyOld = sol[darcyIdx]; + + auto solOld = sol; + + // intialize the coupling manager + couplingManager->init(ransProblem, darcyProblem, sol); + + // intializethe grid variables and the subproblems + using RANSGridVariables = typename GET_PROP_TYPE(RANSTypeTag, GridVariables); + auto ransGridVariables = std::make_shared(ransProblem, ransFvGridGeometry); + ransGridVariables->init(ransSol, solRANSOld); + using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables); + auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry); + darcyGridVariables->init(sol[darcyIdx], solDarcyOld); + + // intialize the vtk output module + const auto ransName = getParam("Problem.Name") + "_" + ransProblem->name(); + const auto darcyName = getParam("Problem.Name") + "_" + darcyProblem->name(); + + StaggeredVtkOutputModule ransVtkWriter(*ransGridVariables, ransSol, ransName); + GET_PROP_TYPE(RANSTypeTag, VtkOutputFields)::init(ransVtkWriter); + ransVtkWriter.write(0.0); + + VtkOutputModule darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName); + using DarcyVelocityOutput = typename GET_PROP_TYPE(DarcyTypeTag, VelocityOutput); + darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); + GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter); + darcyVtkWriter.write(0.0); + + // the assembler with time loop for instationary problem + using Assembler = MultiDomainFVAssembler; + auto assembler = std::make_shared(std::make_tuple(ransProblem, ransProblem, darcyProblem), + std::make_tuple(ransFvGridGeometry->cellCenterFVGridGeometryPtr(), + ransFvGridGeometry->faceFVGridGeometryPtr(), + darcyFvGridGeometry), + std::make_tuple(ransGridVariables->cellCenterGridVariablesPtr(), + ransGridVariables->faceGridVariablesPtr(), + darcyGridVariables), + couplingManager, + timeLoop); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + + // the primary variable switches used by the sub models + using PriVarSwitchTuple = std::tuple; + + // the non-linear solver + using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(solOld); + + // solve the non-linear system with time step control + nonLinearSolver.solve(sol, *timeLoop); + + // make the new solution the old solution + solOld = sol; + + // update the auxiliary free flow solution vector + ransSol[ransCellCenterIdx] = sol[ransCellCenterIdx]; + ransSol[ransFaceIdx] = sol[ransFaceIdx]; + ransProblem->updateDynamicWallProperties(ransSol); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + // advance grid variables to the next time step + ransGridVariables->advanceTimeStep(); + darcyGridVariables->advanceTimeStep(); + + // write vtk output + ransVtkWriter.write(timeLoop->time()); + darcyVtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton solver + timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); + + } while (!timeLoop->finished()); + + timeLoop->finalize(ransGridView.comm()); + timeLoop->finalize(darcyGridView.comm()); + + //////////////////////////////////////////////////////////// + // 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/ransforchheimer/test_rans1pdarcy1pforchheimer.input b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input new file mode 100644 index 0000000000000000000000000000000000000000..972f69a57435246184b786c2aeb7009c9636734f --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input @@ -0,0 +1,72 @@ +[TimeLoop] +DtInitial = 1e-5 # [s] +MaxTimeStepSize = 2 # [s] +TEnd = 10 # [s] + +[RANS.Grid] +Positions0 = -10.0 0.0 0.5 1.0 11.0 +Positions1 = 0.5 1.0 1.5 +Grading0 = -1.15 1.25 -1.25 1.15 +Grading1 = 1.15 -1.3 +Cells0 = 30 20 20 30 +Cells1 = 40 20 +Verbosity = true +DarcyEnd = 1.0 +DarcyStart = 0.0 + +[Darcy.Grid] +Positions0 = 0.0 0.5 1.0 +Positions1 = 0.0 0.5 +Cells0 = 20 20 +Cells1 = 40 +Grading0 = 1.25 -1.25 +Grading1 = -1.15 +Verbosity = true + +[RANS.Problem] +Name = rans +InletReynoldsNumber = 1e5 # [m/s] +InletPressure = 1e5 # [Pa] +Temperature = 298.15 # [K] +SmoothingZoneDistance = 0.02 +SmoothingSlipVelocity = 0.05 + +[Darcy.Problem] +Name = darcy +Pressure = 1e5 +Temperature = 298.15 # [K] +BaseClosed = true +BottomVelocity = 0.606 # [m/s] + +[Darcy.SpatialParams] +Porosity = 0.75 +Permeability = 3.1e-7 +AlphaBeaversJoseph = 2.0 +FittedAlpha = true + +[SpatialParams] +ForchCoeff = 2.65e-2 # + +[Component] +GasDensity = 1.2 +GasKinematicViscosity = 1.5166e-5 + +[Problem] +Name = test_rans1pdarcy1pforchheimer_closed +EnableGravity = false + +[RANS.RANS] +UseStoredEddyViscosity = false + +[Vtk] +AddVelocity = true +WriteFaceData = false + +[Newton] +TargetSteps = 6 +MaxSteps = 12 +MaxRelativeShift = 1e-8 + +[Assembly] +NumericDifferenceMethod = 0 +NumericDifference.BaseEpsilon = 1e-8