From ed9e6044d1e106cedfbbae37dfd3077121b12d2a Mon Sep 17 00:00:00 2001 From: Thomas Fetzer Date: Fri, 20 Jul 2018 16:05:23 +0200 Subject: [PATCH 01/10] TEMP [multidomain][ransdarcy] First draft runs, TODOs * other models * tests --- dumux/freeflow/rans/oneeq/problem.hh | 10 +- dumux/freeflow/rans/twoeq/kepsilon/problem.hh | 10 +- dumux/freeflow/rans/twoeq/komega/problem.hh | 8 +- .../rans/twoeq/lowrekepsilon/problem.hh | 11 +- test/multidomain/boundary/CMakeLists.txt | 1 + .../boundary/ransdarcy/CMakeLists.txt | 61 +++ .../boundary/ransdarcy/darcyproblem.hh | 447 +++++++++++++++++ .../boundary/ransdarcy/ransproblem.hh | 463 ++++++++++++++++++ .../boundary/ransdarcy/spatialparams.hh | 140 ++++++ .../ransdarcy/test_rans1p2cnidarcy2p2cni.cc | 303 ++++++++++++ .../test_rans1p2cnidarcy2p2cni.input | 64 +++ 11 files changed, 1508 insertions(+), 10 deletions(-) create mode 100644 test/multidomain/boundary/ransdarcy/CMakeLists.txt create mode 100644 test/multidomain/boundary/ransdarcy/darcyproblem.hh create mode 100644 test/multidomain/boundary/ransdarcy/ransproblem.hh create mode 100644 test/multidomain/boundary/ransdarcy/spatialparams.hh create mode 100644 test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc create mode 100644 test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input diff --git a/dumux/freeflow/rans/oneeq/problem.hh b/dumux/freeflow/rans/oneeq/problem.hh index 6788b7506b..730eacb9da 100644 --- a/dumux/freeflow/rans/oneeq/problem.hh +++ b/dumux/freeflow/rans/oneeq/problem.hh @@ -64,9 +64,13 @@ class OneEqProblem : public RANSProblem using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; public: - //! The constructor sets the gravity, if desired by the user. - OneEqProblem(std::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry) + /* + * \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) { useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", false); diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh index 8333634b66..13b685c451 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh @@ -82,9 +82,13 @@ class KEpsilonProblem : public RANSProblem public: static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); - //! The constructor sets the gravity, if desired by the user. - KEpsilonProblem(std::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry) + /* + * \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) { yPlusThreshold_ = getParamFromGroup(this->paramGroup(), "KEpsilon.YPlusThreshold", 30); useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", false); diff --git a/dumux/freeflow/rans/twoeq/komega/problem.hh b/dumux/freeflow/rans/twoeq/komega/problem.hh index 05c6d95d00..0814aa4b52 100644 --- a/dumux/freeflow/rans/twoeq/komega/problem.hh +++ b/dumux/freeflow/rans/twoeq/komega/problem.hh @@ -70,7 +70,13 @@ class KOmegaProblem : public RANSProblem using DimMatrix = Dune::FieldMatrix; public: - KOmegaProblem(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) + /* + * \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) { useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", false); } diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh index 742b6e790c..bdd67fd27e 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh @@ -61,9 +61,14 @@ class LowReKEpsilonProblem : public RANSProblem using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; public: - //! The constructor sets the gravity, if desired by the user. - LowReKEpsilonProblem(std::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry) + + /* + * \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) { useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", true); } diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt index addd460feb..7a82048d4d 100644 --- a/test/multidomain/boundary/CMakeLists.txt +++ b/test/multidomain/boundary/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(darcydarcy) +add_subdirectory(ransdarcy) add_subdirectory(stokesdarcy) diff --git a/test/multidomain/boundary/ransdarcy/CMakeLists.txt b/test/multidomain/boundary/ransdarcy/CMakeLists.txt new file mode 100644 index 0000000000..d2d569c9be --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/CMakeLists.txt @@ -0,0 +1,61 @@ +add_input_file_links() + +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") +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") +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") +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") +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") +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 0000000000..b6e3232021 --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/darcyproblem.hh @@ -0,0 +1,447 @@ +// -*- 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(); + 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'"); + if (plotStorage_) + gnuplotStorage_.plot("temp"); + } + + 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_.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 0000000000..110b5035cf --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/ransproblem.hh @@ -0,0 +1,463 @@ +// -*- 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, OneEqNCNI)); +#elif KOMEGA +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +#elif LOWREKEPSILON +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +#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 OneEqProblem +{ + using ParentType = OneEqProblem; +#elif KOMEGA +class FreeFlowSubProblem : public OneEqProblem +{ + using ParentType = OneEqProblem; +#elif LOWREKEPSILON +class FreeFlowSubProblem : public OneEqProblem +{ + using ParentType = OneEqProblem; +#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 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"); + inletMoleFrac_ = getParamFromGroup(this->paramGroup(), "Problem.InletMoleFrac"); + 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.setMoleFraction(0, 1, inletMoleFrac_); + fluidState.setMoleFraction(0, 0, 1.0 - inletMoleFrac_); + fluidState.setTemperature(inletTemperature_); + Scalar density = FluidSystem::density(fluidState, 0); + Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, 0) / density; + Scalar radius = 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_, radius*2.0, kinematicViscosity, true)); + turbulentKineticEnergy_ = getParam("Problem.InletTurbulentKineticEnergy", + turbulenceProperties.turbulentKineticEnergy(inletVelocity_, radius*2.0, kinematicViscosity, true)); +#if KOMEGA + dissipation_ = getParam("Problem.InletDissipationRate", + turbulenceProperties.dissipationRate(inletVelocity_, radius*2.0, kinematicViscosity, true)); +#else + dissipation_ = getParam("Problem.InletDissipation", + turbulenceProperties.dissipation(inletVelocity_, radius*2.0, kinematicViscosity, true)); +#endif + std::cout << std::endl; + } + + /*! + * \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 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 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(Indices::momentumYBalanceIdx); + bTypes.setBJS(Indices::momentumXBalanceIdx); + } + return bTypes; + } + + /*! + * \brief Evaluate the boundary conditions for a Dirichlet control volume. + * + * \param element The element + * \param scvf The subcontrolvolume face + */ + PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const + { + PrimaryVariables values(0.0); + values = initialAtPos(pos); + + 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 ONEEQ + values[Indices::viscosityTildeIdx] = viscosityTilde_; + if (isOnWall(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 0000000000..4bbd873bb8 --- /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(std::max(threshold,getParam("Darcy.SpatialParams.PcLowSw", 0.0))); + params_.setPcHighSw(std::min(1.0-threshold,getParam("Darcy.SpatialParams.PcHighSw", 1.0))); + } + + /*! + * \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/test_rans1p2cnidarcy2p2cni.cc b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc new file mode 100644 index 0000000000..ce23bd3fdb --- /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); + + // post time step treatment of Darcy problem + darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables); + + // advance grid variables to the next time step + ransGridVariables->advanceTimeStep(); + darcyGridVariables->advanceTimeStep(); + + // advance to the time loop to the next step + timeLoop->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 0000000000..e91d99650c --- /dev/null +++ b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input @@ -0,0 +1,64 @@ +[TimeLoop] +DtInitial = 1e-1 # [s] +MaxTimeStepSize = 43200 # [s] (12 hours) +TEnd = 864000 # [s] (6 days) + +[RANS.Grid] +Positions0 = -0.25 0.0 0.25 0.5 +Positions1 = 0.25 0.5 +Grading0 = 1.0 1.0 1.0 +Grading1 = 1.5 +Cells0 = 2 10 2 +Cells1 = 15 +Verbosity = true + +[Darcy.Grid] +Positions0 = 0.0 0.25 +Positions1 = 0.0 0.25 +Cells0 = 10 +Cells1 = 15 +Grading0 = 1.0 +Grading1 = -1.5 +Verbosity = true + +[RANS.Problem] +Name = rans +InletVelocity = 3.5 # [m/s] +InletPressure = 1e5 # [Pa] +InletMoleFrac = 0 # [-] +InletTemperature = 298.15 # [K] + +[Darcy.Problem] +Name = darcy +Pressure = 1e5 +Saturation = 0.95 # 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.371e-4 +VgN = 6.9 + +[Problem] +Name = test_rans1p2cnidarcy2p2cni +EnableGravity = true +InterfaceDiffusionCoefficientAvg = FreeFlowOnly +PlotFluxes = true +PlotStorage = true + +[Vtk] +AddVelocity = true +WriteFaceData = false + +[Newton] +MaxSteps = 12 +MaxRelativeShift = 1e-5 + +[Assembly] +NumericDifferenceMethod = 0 +NumericDifference.BaseEpsilon = 1e-8 -- GitLab From 1cd70da7c8590b172a48f3de0cc8db31289d1870 Mon Sep 17 00:00:00 2001 From: Thomas Fetzer Date: Fri, 20 Jul 2018 16:05:23 +0200 Subject: [PATCH 02/10] TEMP [multidomain][ransdarcy] First draft runs, TODOs * tests --- .../boundary/ransdarcy/CMakeLists.txt | 16 ++- .../boundary/ransdarcy/darcyproblem.hh | 10 +- .../boundary/ransdarcy/ransproblem.hh | 113 +++++++++++++----- .../boundary/ransdarcy/spatialparams.hh | 4 +- test/multidomain/boundary/ransdarcy/temp.sh | 5 + .../ransdarcy/test_rans1p2cnidarcy2p2cni.cc | 6 +- .../test_rans1p2cnidarcy2p2cni.input | 38 +++--- 7 files changed, 133 insertions(+), 59 deletions(-) create mode 100644 test/multidomain/boundary/ransdarcy/temp.sh diff --git a/test/multidomain/boundary/ransdarcy/CMakeLists.txt b/test/multidomain/boundary/ransdarcy/CMakeLists.txt index d2d569c9be..8b84477d3e 100644 --- a/test/multidomain/boundary/ransdarcy/CMakeLists.txt +++ b/test/multidomain/boundary/ransdarcy/CMakeLists.txt @@ -1,4 +1,5 @@ add_input_file_links() +dune_symlink_to_source_files(FILES temp.sh) dune_add_test(NAME test_kepsilon1p2cnidarcy2p2cni SOURCES test_rans1p2cnidarcy2p2cni.cc @@ -9,7 +10,8 @@ dune_add_test(NAME test_kepsilon1p2cnidarcy2p2cni ${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") + --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 @@ -21,7 +23,8 @@ dune_add_test(NAME test_komega1p2cnidarcy2p2cni ${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") + --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 @@ -33,7 +36,8 @@ dune_add_test(NAME test_lowrekepsilon1p2cnidarcy2p2cni ${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") + --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 @@ -45,7 +49,8 @@ dune_add_test(NAME test_oneeq1p2cnidarcy2p2cni ${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") + --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 @@ -57,5 +62,6 @@ dune_add_test(NAME test_zeroeq1p2cnidarcy2p2cni ${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") + --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 index b6e3232021..97b3777f3d 100644 --- a/test/multidomain/boundary/ransdarcy/darcyproblem.hh +++ b/test/multidomain/boundary/ransdarcy/darcyproblem.hh @@ -155,6 +155,8 @@ public: 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); @@ -167,8 +169,13 @@ public: 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("temp"); + gnuplotStorage_.plot("evapRate_"); } template @@ -258,6 +265,7 @@ public: } 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)]"); diff --git a/test/multidomain/boundary/ransdarcy/ransproblem.hh b/test/multidomain/boundary/ransdarcy/ransproblem.hh index 110b5035cf..f56fbb3d7f 100644 --- a/test/multidomain/boundary/ransdarcy/ransproblem.hh +++ b/test/multidomain/boundary/ransdarcy/ransproblem.hh @@ -34,14 +34,14 @@ #include #include #elif KEPSILON -#include -#include +#include +#include #elif KOMEGA -#include -#include +#include +#include #elif LOWREKEPSILON -#include -#include +#include +#include #else #include #include @@ -57,11 +57,11 @@ namespace Properties #if ONEEQ NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); #elif KEPSILON -NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KEpsilonNCNI)); #elif KOMEGA -NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmegaNCNI)); #elif LOWREKEPSILON -NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI)); +NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, LowReKEpsilonNCNI)); #else NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI)); #endif @@ -101,17 +101,17 @@ class FreeFlowSubProblem : public OneEqProblem { using ParentType = OneEqProblem; #elif KEPSILON -class FreeFlowSubProblem : public OneEqProblem +class FreeFlowSubProblem : public KEpsilonProblem { - using ParentType = OneEqProblem; + using ParentType = KEpsilonProblem; #elif KOMEGA -class FreeFlowSubProblem : public OneEqProblem +class FreeFlowSubProblem : public KOmegaProblem { - using ParentType = OneEqProblem; + using ParentType = KOmegaProblem; #elif LOWREKEPSILON -class FreeFlowSubProblem : public OneEqProblem +class FreeFlowSubProblem : public LowReKEpsilonProblem { - using ParentType = OneEqProblem; + using ParentType = LowReKEpsilonProblem; #else class FreeFlowSubProblem : public ZeroEqProblem { @@ -126,6 +126,7 @@ class FreeFlowSubProblem : public ZeroEqProblem 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; @@ -150,7 +151,7 @@ public: { inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.InletVelocity"); inletPressure_ = getParamFromGroup(this->paramGroup(), "Problem.InletPressure"); - inletMoleFrac_ = getParamFromGroup(this->paramGroup(), "Problem.InletMoleFrac"); + auto inletMassFrac = getParamFromGroup(this->paramGroup(), "Problem.InletMassFrac"); inletTemperature_ = getParamFromGroup(this->paramGroup(), "Problem.InletTemperature"); diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{}, @@ -159,23 +160,24 @@ public: Dumux::TurbulenceProperties turbulenceProperties; FluidState fluidState; fluidState.setPressure(0, inletPressure_); - fluidState.setMoleFraction(0, 1, inletMoleFrac_); - fluidState.setMoleFraction(0, 0, 1.0 - inletMoleFrac_); + 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 radius = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; + 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_, radius*2.0, kinematicViscosity, true)); + 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, charLength, kinematicViscosity, true)); turbulentKineticEnergy_ = getParam("Problem.InletTurbulentKineticEnergy", - turbulenceProperties.turbulentKineticEnergy(inletVelocity_, radius*2.0, kinematicViscosity, true)); + turbulenceProperties.turbulentKineticEnergy(inletVelocity_, charLength, kinematicViscosity, true)); #if KOMEGA dissipation_ = getParam("Problem.InletDissipationRate", - turbulenceProperties.dissipationRate(inletVelocity_, radius*2.0, kinematicViscosity, true)); + turbulenceProperties.dissipationRate(inletVelocity_, charLength, kinematicViscosity, true)); #else dissipation_ = getParam("Problem.InletDissipation", - turbulenceProperties.dissipation(inletVelocity_, radius*2.0, kinematicViscosity, true)); + turbulenceProperties.dissipation(inletVelocity_, charLength, kinematicViscosity, true)); #endif std::cout << std::endl; } @@ -219,8 +221,13 @@ public: const auto& globalPos = scvf.center(); +#if LOWREKEPSILON || KEPSILON || KOMEGA + bTypes.setDirichlet(Indices::turbulentKineticEnergyIdx); + bTypes.setDirichlet(Indices::dissipationIdx); +#endif + #if ONEEQ - bTypes.setDirichlet(Indices::viscosityTildeIdx); + bTypes.setDirichlet(Indices::viscosityTildeIdx); #endif if (onLeftBoundary_(globalPos)) @@ -250,6 +257,10 @@ public: 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 @@ -263,20 +274,46 @@ public: bTypes.setCouplingNeumann(Indices::momentumYBalanceIdx); bTypes.setBJS(Indices::momentumXBalanceIdx); } + +#if KOMEGA + // set a fixed dissipation (omega) in one cell + if (isOnWall(globalPos)) + bTypes.setDirichletCell(Indices::dissipationIdx); +#endif + return bTypes; } - /*! - * \brief Evaluate the boundary conditions for a Dirichlet control volume. - * - * \param element The element - * \param scvf The subcontrolvolume face - */ - PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const + /*! + * \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 { - PrimaryVariables values(0.0); - values = initialAtPos(pos); + 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; } @@ -353,6 +390,16 @@ public: if(onLowerBoundary_(globalPos)) values[Indices::velocityXIdx] = 0.0; +#if LOWREKEPSILON || KEPSILON || KOMEGA + values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; + values[Indices::dissipationIdx] = dissipation_; + if (isOnWall(globalPos)) + { + values[Indices::turbulentKineticEnergyIdx] = 0.0; + values[Indices::dissipationIdx] = 0.0; + } +#endif + #if ONEEQ values[Indices::viscosityTildeIdx] = viscosityTilde_; if (isOnWall(globalPos)) diff --git a/test/multidomain/boundary/ransdarcy/spatialparams.hh b/test/multidomain/boundary/ransdarcy/spatialparams.hh index 4bbd873bb8..dd70bbfbd3 100644 --- a/test/multidomain/boundary/ransdarcy/spatialparams.hh +++ b/test/multidomain/boundary/ransdarcy/spatialparams.hh @@ -75,8 +75,8 @@ public: params_.setVgAlpha(getParam("Darcy.SpatialParams.VgAlpha")); params_.setVgn(getParam("Darcy.SpatialParams.VgN")); Scalar threshold = 0.01 * (1.0 - params_.swr() - params_.snr()); - params_.setPcLowSw(std::max(threshold,getParam("Darcy.SpatialParams.PcLowSw", 0.0))); - params_.setPcHighSw(std::min(1.0-threshold,getParam("Darcy.SpatialParams.PcHighSw", 1.0))); + params_.setPcLowSw(getParam("Darcy.SpatialParams.PcLowSw", threshold)); + params_.setPcHighSw(getParam("Darcy.SpatialParams.PcHighSw", 1.0-threshold)); } /*! diff --git a/test/multidomain/boundary/ransdarcy/temp.sh b/test/multidomain/boundary/ransdarcy/temp.sh new file mode 100644 index 0000000000..fb35996564 --- /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 index ce23bd3fdb..5d174304a2 100644 --- a/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc +++ b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.cc @@ -239,6 +239,9 @@ int main(int argc, char** argv) try 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); @@ -246,9 +249,6 @@ int main(int argc, char** argv) try ransGridVariables->advanceTimeStep(); darcyGridVariables->advanceTimeStep(); - // advance to the time loop to the next step - timeLoop->advanceTimeStep(); - // write vtk output ransVtkWriter.write(timeLoop->time()); darcyVtkWriter.write(timeLoop->time()); diff --git a/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input index e91d99650c..e955974c5d 100644 --- a/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input +++ b/test/multidomain/boundary/ransdarcy/test_rans1p2cnidarcy2p2cni.input @@ -4,34 +4,34 @@ MaxTimeStepSize = 43200 # [s] (12 hours) TEnd = 864000 # [s] (6 days) [RANS.Grid] -Positions0 = -0.25 0.0 0.25 0.5 +Positions0 = -1.0 0.0 0.5 Positions1 = 0.25 0.5 -Grading0 = 1.0 1.0 1.0 -Grading1 = 1.5 -Cells0 = 2 10 2 -Cells1 = 15 +Grading0 = 1.0 1.0 +Grading1 = 1.57539 +Cells0 = 4 16 +Cells1 = 16 Verbosity = true [Darcy.Grid] -Positions0 = 0.0 0.25 +Positions0 = 0.0 0.5 Positions1 = 0.0 0.25 -Cells0 = 10 -Cells1 = 15 +Cells0 = 16 +Cells1 = 16 Grading0 = 1.0 -Grading1 = -1.5 +Grading1 = -1.57539 Verbosity = true [RANS.Problem] Name = rans InletVelocity = 3.5 # [m/s] InletPressure = 1e5 # [Pa] -InletMoleFrac = 0 # [-] +InletMassFrac = 0.008 # [-] InletTemperature = 298.15 # [K] [Darcy.Problem] Name = darcy Pressure = 1e5 -Saturation = 0.95 # initial Sw +Saturation = 0.8 # initial Sw Temperature = 298.15 # [K] InitPhasePresence = 3 # bothPhases @@ -41,21 +41,29 @@ Permeability = 2.65e-10 AlphaBeaversJoseph = 1.0 Swr = 0.005 Snr = 0.01 -VgAlpha = 6.371e-4 -VgN = 6.9 +VgAlpha = 6.37e-4 +VgN = 8.0 [Problem] Name = test_rans1p2cnidarcy2p2cni EnableGravity = true -InterfaceDiffusionCoefficientAvg = FreeFlowOnly -PlotFluxes = 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 -- GitLab From 28684a83a51e87210383749bfa069148b129aca8 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Fri, 12 Oct 2018 08:19:43 +0200 Subject: [PATCH 03/10] Update wall bools, rebase --- .../boundary/ransdarcy/ransproblem.hh | 32 +++++++++++++++---- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/test/multidomain/boundary/ransdarcy/ransproblem.hh b/test/multidomain/boundary/ransdarcy/ransproblem.hh index f56fbb3d7f..0ce62b22ba 100644 --- a/test/multidomain/boundary/ransdarcy/ransproblem.hh +++ b/test/multidomain/boundary/ransdarcy/ransproblem.hh @@ -89,7 +89,7 @@ SET_BOOL_PROP(RANSTypeTag, EnableFVGridGeometryCache, true); SET_BOOL_PROP(RANSTypeTag, EnableGridFluxVariablesCache, true); SET_BOOL_PROP(RANSTypeTag, EnableGridVolumeVariablesCache, true); -SET_BOOL_PROP(RANSTypeTag, EnableInertiaTerms, true); +// SET_BOOL_PROP(RANSTypeTag, EnableInertiaTerms, true); } /*! @@ -180,8 +180,25 @@ public: 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 */ @@ -271,13 +288,14 @@ public: bTypes.setCouplingNeumann(Indices::conti0EqIdx); bTypes.setCouplingNeumann(Indices::conti0EqIdx + 1); bTypes.setCouplingNeumann(Indices::energyBalanceIdx); - bTypes.setCouplingNeumann(Indices::momentumYBalanceIdx); - bTypes.setBJS(Indices::momentumXBalanceIdx); + + bTypes.setCouplingNeumann(scvf.directionIndex()); + bTypes.setBJS(1 - scvf.directionIndex()); } #if KOMEGA // set a fixed dissipation (omega) in one cell - if (isOnWall(globalPos)) + if (isOnWallAtPos(globalPos)) bTypes.setDirichletCell(Indices::dissipationIdx); #endif @@ -363,7 +381,6 @@ public: { return (onLowerBoundary_(globalPos)); } - /*! * \name Volume terms */ @@ -393,7 +410,8 @@ public: #if LOWREKEPSILON || KEPSILON || KOMEGA values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; values[Indices::dissipationIdx] = dissipation_; - if (isOnWall(globalPos)) + + if (isOnWallAtPos(globalPos)) { values[Indices::turbulentKineticEnergyIdx] = 0.0; values[Indices::dissipationIdx] = 0.0; @@ -402,7 +420,7 @@ public: #if ONEEQ values[Indices::viscosityTildeIdx] = viscosityTilde_; - if (isOnWall(globalPos)) + if (isOnWallAtPos(globalPos)) { values[Indices::viscosityTildeIdx] = 0.0; } -- GitLab From 3e96922bcffc5ec9606637040fb60837757288d1 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Fri, 12 Oct 2018 15:22:42 +0200 Subject: [PATCH 04/10] Develop a Komega-Forcheimer coupled model for Guang's paper --- test/multidomain/boundary/CMakeLists.txt | 1 + .../boundary/ransforchheimer/CMakeLists.txt | 7 + .../boundary/ransforchheimer/darcyproblem.hh | 289 ++++++++++++++ .../boundary/ransforchheimer/ransproblem.hh | 371 ++++++++++++++++++ .../boundary/ransforchheimer/spatialparams.hh | 95 +++++ .../test_rans1pdarcy1pforchheimer.cc | 299 ++++++++++++++ .../test_rans1pdarcy1pforchheimer.input | 61 +++ 7 files changed, 1123 insertions(+) create mode 100644 test/multidomain/boundary/ransforchheimer/CMakeLists.txt create mode 100644 test/multidomain/boundary/ransforchheimer/darcyproblem.hh create mode 100644 test/multidomain/boundary/ransforchheimer/ransproblem.hh create mode 100644 test/multidomain/boundary/ransforchheimer/spatialparams.hh create mode 100644 test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.cc create mode 100644 test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt index 7a82048d4d..4d35f107a2 100644 --- a/test/multidomain/boundary/CMakeLists.txt +++ b/test/multidomain/boundary/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(darcydarcy) add_subdirectory(ransdarcy) +add_subdirectory(ransforchheimer) add_subdirectory(stokesdarcy) diff --git a/test/multidomain/boundary/ransforchheimer/CMakeLists.txt b/test/multidomain/boundary/ransforchheimer/CMakeLists.txt new file mode 100644 index 0000000000..d204e1e993 --- /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 0000000000..542e382f53 --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh @@ -0,0 +1,289 @@ +// -*- 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) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + 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 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"); + } + + /*! + * \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(); + + 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); + + 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_; + + 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 0000000000..78583bc327 --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/ransproblem.hh @@ -0,0 +1,371 @@ +// -*- 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) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + 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"); + + 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 isOnWallAtPos(const GlobalPosition& globalPos) const + { + return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)); + } + // \} + + /*! + * \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 (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + bTypes.setCouplingNeumann(Indices::conti0EqIdx); + bTypes.setCouplingNeumann(scvf.directionIndex()); + bTypes.setBJS(1 - scvf.directionIndex()); + } + + // set a fixed dissipation (omega) in one cell + if (isOnWallAtPos(globalPos)) + 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_; } + + 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 + { + 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 (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_; + + 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 0000000000..db03f5787f --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/spatialparams.hh @@ -0,0 +1,95 @@ +// -*- 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"); + } + + /*! + * \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_; } + + +private: + Scalar permeability_; + Scalar porosity_; + Scalar alphaBJ_; +}; + +} // 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 0000000000..a188b01c66 --- /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 0000000000..816d35d053 --- /dev/null +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input @@ -0,0 +1,61 @@ +[TimeLoop] +DtInitial = 5e-4 # [s] +MaxTimeStepSize = 1 # [s] +TEnd = 2 # [s] + +[RANS.Grid] +Positions0 = -1.0 0.0 1.0 2.0 +Positions1 = 0.5 1.0 1.5 +Grading0 = 1.0 1.0 1.0 +Grading1 = 1.3 -1.3 +Cells0 = 10 20 5 +Cells1 = 20 20 +Verbosity = true + +[Darcy.Grid] +Positions0 = 0.0 1.0 +Positions1 = 0.0 0.5 +Cells0 = 20 +Cells1 = 20 +Grading0 = 1.0 +Grading1 = -1.3 +Verbosity = true + +[RANS.Problem] +Name = rans +InletReynoldsNumber = 1e6 # [m/s] +InletPressure = 1e5 # [Pa] +Temperature = 298.15 # [K] + +[Darcy.Problem] +Name = darcy +Pressure = 1e5 +Temperature = 298.15 # [K] + +[Darcy.SpatialParams] +Porosity = 0.41 +Permeability = 1e-6 +AlphaBeaversJoseph = 1.0 + +[SpatialParams] +ForchCoeff = 0.55 + +[Problem] +Name = test_rans1pdarcy1pforchheimer +EnableGravity = false + +[RANS.RANS] +UseStoredEddyViscosity = false + +[Vtk] +AddVelocity = true +WriteFaceData = false + +[Newton] +TargetSteps = 7 +MaxSteps = 12 +MaxRelativeShift = 1e-5 + +[Assembly] +NumericDifferenceMethod = 0 +NumericDifference.BaseEpsilon = 1e-8 -- GitLab From 4d2bcd12a5978b267807e1a861c7b0b2d9679ef7 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Thu, 18 Oct 2018 14:13:47 +0000 Subject: [PATCH 05/10] [stokesdarcy][couplingdata] Add pressure reconstruction for Forchheimer (cherry picked from commit 094ecf4ca8b1fab5037c4a4097ceb45e0c7b2c80) --- .../boundary/stokesdarcy/couplingdata.hh | 88 ++++++++++++++++--- 1 file changed, 76 insertions(+), 12 deletions(-) diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh index 9dcb8dce71..e835707b63 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_; -- GitLab From 84641077ec1d9ba1f2ff08541cc662a76e14562f Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Wed, 24 Oct 2018 15:24:39 +0200 Subject: [PATCH 06/10] Added slip zone to wall near the coupled interface --- .../boundary/ransforchheimer/darcyproblem.hh | 6 ++- .../boundary/ransforchheimer/ransproblem.hh | 41 ++++++++++++++---- .../test_rans1pdarcy1pforchheimer.input | 43 ++++++++++++------- 3 files changed, 63 insertions(+), 27 deletions(-) diff --git a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh index 542e382f53..ce30689f61 100644 --- a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh +++ b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh @@ -31,7 +31,7 @@ #include #include -#include +#include #if FORCHHEIMER #include @@ -54,8 +54,10 @@ SET_TYPE_PROP(DarcyOnePTypeTag, Problem, Dumux::DarcySubProblem); // the fluid system SET_PROP(DarcyOnePTypeTag, FluidSystem) { +private: using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using type = FluidSystems::OnePGas >; +public: + using type = FluidSystems::OnePGas >; }; // Set the grid type diff --git a/test/multidomain/boundary/ransforchheimer/ransproblem.hh b/test/multidomain/boundary/ransforchheimer/ransproblem.hh index 78583bc327..cdb42a6d54 100644 --- a/test/multidomain/boundary/ransforchheimer/ransproblem.hh +++ b/test/multidomain/boundary/ransforchheimer/ransproblem.hh @@ -26,7 +26,7 @@ #include #include -#include +#include #include #include @@ -46,8 +46,10 @@ NEW_TYPE_TAG(RANSTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmega)); // the fluid system SET_PROP(RANSTypeTag, FluidSystem) { +private: using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using type = FluidSystems::OnePGas >; +public: + using type = FluidSystems::OnePGas >; }; // Set the grid type @@ -100,6 +102,11 @@ public: 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_); @@ -126,12 +133,24 @@ public: 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 */ @@ -186,7 +205,7 @@ public: bTypes.setDirichlet(Indices::velocityYIdx); } - if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + if (isOnCouplingWall(scvf)) { bTypes.setCouplingNeumann(Indices::conti0EqIdx); bTypes.setCouplingNeumann(scvf.directionIndex()); @@ -194,7 +213,7 @@ public: } // set a fixed dissipation (omega) in one cell - if (isOnWallAtPos(globalPos)) + if (isOnWall(scvf)) bTypes.setDirichletCell(Indices::dissipationIdx); return bTypes; @@ -270,11 +289,6 @@ public: const CouplingManager& couplingManager() const { return *couplingManager_; } - bool isOnWall(const GlobalPosition& globalPos) const - { - return (onLowerBoundary_(globalPos)); - } - /*! * \name Volume terms */ @@ -298,6 +312,11 @@ public: values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; values[Indices::dissipationIdx] = dissipation_; + if (isSmoothingZoneAtPos(globalPos)) + { + values[Indices::velocityXIdx] = smoothingSlipVelocity_; + } + if (isOnWallAtPos(globalPos)) { values[Indices::turbulentKineticEnergyIdx] = 0.0; @@ -360,6 +379,10 @@ private: Scalar temperature_; Scalar turbulentKineticEnergy_; Scalar dissipation_; + Scalar darcyEnd_; + Scalar darcyStart_; + Scalar smoothingZoneDistance_; + Scalar smoothingSlipVelocity_; TimeLoopPtr timeLoop_; diff --git a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input index 816d35d053..afebae9183 100644 --- a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input @@ -1,31 +1,35 @@ [TimeLoop] -DtInitial = 5e-4 # [s] -MaxTimeStepSize = 1 # [s] -TEnd = 2 # [s] +DtInitial = 1e-4 # [s] +MaxTimeStepSize = 2 # [s] +TEnd = 10 # [s] [RANS.Grid] -Positions0 = -1.0 0.0 1.0 2.0 +Positions0 = -1.0 0.0 0.5 1.0 2.0 3.0 Positions1 = 0.5 1.0 1.5 -Grading0 = 1.0 1.0 1.0 +Grading0 = -1.40 1.25 -1.25 1.40 1.0 Grading1 = 1.3 -1.3 -Cells0 = 10 20 5 +Cells0 = 15 20 20 15 5 Cells1 = 20 20 Verbosity = true +DarcyEnd = 1.0 +DarcyStart = 0.0 [Darcy.Grid] -Positions0 = 0.0 1.0 +Positions0 = 0.0 0.5 1.0 Positions1 = 0.0 0.5 -Cells0 = 20 +Cells0 = 20 20 Cells1 = 20 -Grading0 = 1.0 +Grading0 = 1.25 -1.25 Grading1 = -1.3 Verbosity = true [RANS.Problem] Name = rans -InletReynoldsNumber = 1e6 # [m/s] +InletReynoldsNumber = 1e5 # [m/s] InletPressure = 1e5 # [Pa] Temperature = 298.15 # [K] +SmoothingZoneDistance = 0.02 +SmoothingSlipVelocity = 0.05 [Darcy.Problem] Name = darcy @@ -33,12 +37,19 @@ Pressure = 1e5 Temperature = 298.15 # [K] [Darcy.SpatialParams] -Porosity = 0.41 -Permeability = 1e-6 -AlphaBeaversJoseph = 1.0 +Porosity = 0.75 +Permeability = 3.1e-7 +AlphaBeaversJoseph = 2.0 +FittedAlpha = true +BaseClosed = true + [SpatialParams] -ForchCoeff = 0.55 +ForchCoeff = 2.65e-2 # + +[Component] +GasDensity = 1.2 +GasKinematicViscosity = 1.82e-5 [Problem] Name = test_rans1pdarcy1pforchheimer @@ -52,9 +63,9 @@ AddVelocity = true WriteFaceData = false [Newton] -TargetSteps = 7 +TargetSteps = 6 MaxSteps = 12 -MaxRelativeShift = 1e-5 +MaxRelativeShift = 1e-8 [Assembly] NumericDifferenceMethod = 0 -- GitLab From ecd596d9eadd90a5d0792b6995941033f10350c5 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Wed, 24 Oct 2018 15:24:54 +0200 Subject: [PATCH 07/10] added fitted function for Guang's Alphas --- .../boundary/ransforchheimer/spatialparams.hh | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/test/multidomain/boundary/ransforchheimer/spatialparams.hh b/test/multidomain/boundary/ransforchheimer/spatialparams.hh index db03f5787f..bc7b04c9a2 100644 --- a/test/multidomain/boundary/ransforchheimer/spatialparams.hh +++ b/test/multidomain/boundary/ransforchheimer/spatialparams.hh @@ -58,6 +58,8 @@ public: permeability_ = getParam("Darcy.SpatialParams.Permeability"); porosity_ = getParam("Darcy.SpatialParams.Porosity"); alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); + fittedAlpha_ = getParam("Darcy.SpatialParams.FittedAlpha"); + baseClosed_ = getParam("Darcy.SpatialParams.BaseClosed"); } /*! @@ -81,13 +83,23 @@ public: * \param globalPos The global position */ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const - { return alphaBJ_; } - + { + 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 + return alphaBJ_; + } private: Scalar permeability_; Scalar porosity_; Scalar alphaBJ_; + bool fittedAlpha_; + bool baseClosed_; }; } // end namespace -- GitLab From 0582ccae3a159c16fd9c986ecbeb7a6217761d68 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Thu, 8 Nov 2018 15:32:38 +0100 Subject: [PATCH 08/10] Add flow from the base --- .../boundary/ransforchheimer/darcyproblem.hh | 16 ++++++++++++++++ .../boundary/ransforchheimer/spatialparams.hh | 2 +- .../test_rans1pdarcy1pforchheimer.input | 6 +++--- 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh index ce30689f61..4b361cf589 100644 --- a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh +++ b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh @@ -114,6 +114,8 @@ public: { pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature"); + baseClosed_ = getParamFromGroup(this->paramGroup(),"Problem.BaseClosed"); + inletRelPressure_ = getParamFromGroup(this->paramGroup(),"Problem.InletRelPressure"); } /*! @@ -165,6 +167,11 @@ public: if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann(); + if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + { + values.setDirichlet(Indices::pressureIdx); + } + return values; } @@ -181,6 +188,11 @@ public: PrimaryVariables values(0.0); values = initial(element); + if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + { + values[pressureIdx] = pressure_ + inletRelPressure_; + } + return values; } @@ -203,7 +215,9 @@ public: NumEqVector values(0.0); if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + { values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf); + } return values; } @@ -282,6 +296,8 @@ private: Scalar pressure_; Scalar temperature_; Scalar eps_; + Scalar inletRelPressure_; + bool baseClosed_; TimeLoopPtr timeLoop_; std::shared_ptr couplingManager_; diff --git a/test/multidomain/boundary/ransforchheimer/spatialparams.hh b/test/multidomain/boundary/ransforchheimer/spatialparams.hh index bc7b04c9a2..b007c5706b 100644 --- a/test/multidomain/boundary/ransforchheimer/spatialparams.hh +++ b/test/multidomain/boundary/ransforchheimer/spatialparams.hh @@ -59,7 +59,7 @@ public: porosity_ = getParam("Darcy.SpatialParams.Porosity"); alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); fittedAlpha_ = getParam("Darcy.SpatialParams.FittedAlpha"); - baseClosed_ = getParam("Darcy.SpatialParams.BaseClosed"); + baseClosed_ = getParam("Darcy.Problem.BaseClosed"); } /*! diff --git a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input index afebae9183..b3f7f4f707 100644 --- a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input @@ -1,5 +1,5 @@ [TimeLoop] -DtInitial = 1e-4 # [s] +DtInitial = 1e-5 # [s] MaxTimeStepSize = 2 # [s] TEnd = 10 # [s] @@ -35,14 +35,14 @@ SmoothingSlipVelocity = 0.05 Name = darcy Pressure = 1e5 Temperature = 298.15 # [K] +BaseClosed = false +InletRelPressure = 10 [Darcy.SpatialParams] Porosity = 0.75 Permeability = 3.1e-7 AlphaBeaversJoseph = 2.0 FittedAlpha = true -BaseClosed = true - [SpatialParams] ForchCoeff = 2.65e-2 # -- GitLab From 4cfa62908e36aadd0bd3c8a6e68c5bc80cca5306 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Tue, 13 Nov 2018 09:39:22 +0100 Subject: [PATCH 09/10] [ransforchheimer] Adapt inflow from bottom --- .../boundary/ransforchheimer/darcyproblem.hh | 30 ++++++++++++------- .../test_rans1pdarcy1pforchheimer.input | 22 +++++++------- 2 files changed, 30 insertions(+), 22 deletions(-) diff --git a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh index 4b361cf589..dc549bba86 100644 --- a/test/multidomain/boundary/ransforchheimer/darcyproblem.hh +++ b/test/multidomain/boundary/ransforchheimer/darcyproblem.hh @@ -93,6 +93,8 @@ class DarcySubProblem : public PorousMediumFlowProblem 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; @@ -115,7 +117,12 @@ public: pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); temperature_ = getParamFromGroup(this->paramGroup(), "Problem.Temperature"); baseClosed_ = getParamFromGroup(this->paramGroup(),"Problem.BaseClosed"); - inletRelPressure_ = getParamFromGroup(this->paramGroup(),"Problem.InletRelPressure"); + bottomVelocity_ = getParamFromGroup(this->paramGroup(),"Problem.BottomVelocity"); + + FluidState fluidState; + fluidState.setPressure(0, pressure_); + fluidState.setTemperature(temperature_); + density_ = FluidSystem::density(fluidState, 0); } /*! @@ -167,10 +174,10 @@ public: if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann(); - if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) - { - values.setDirichlet(Indices::pressureIdx); - } + // if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + // { + // values.setDirichlet(Indices::pressureIdx); + // } return values; } @@ -188,11 +195,6 @@ public: PrimaryVariables values(0.0); values = initial(element); - if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) - { - values[pressureIdx] = pressure_ + inletRelPressure_; - } - return values; } @@ -219,6 +221,11 @@ public: values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf); } + if (!baseClosed_ && onLowerBoundary_(scvf.ipGlobal())) + { + values[Indices::conti0EqIdx] = -bottomVelocity_ * density_; + } + return values; } @@ -296,8 +303,9 @@ private: Scalar pressure_; Scalar temperature_; Scalar eps_; - Scalar inletRelPressure_; bool baseClosed_; + Scalar bottomVelocity_; + Scalar density_; TimeLoopPtr timeLoop_; std::shared_ptr couplingManager_; diff --git a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input index b3f7f4f707..972f69a574 100644 --- a/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input +++ b/test/multidomain/boundary/ransforchheimer/test_rans1pdarcy1pforchheimer.input @@ -4,12 +4,12 @@ MaxTimeStepSize = 2 # [s] TEnd = 10 # [s] [RANS.Grid] -Positions0 = -1.0 0.0 0.5 1.0 2.0 3.0 +Positions0 = -10.0 0.0 0.5 1.0 11.0 Positions1 = 0.5 1.0 1.5 -Grading0 = -1.40 1.25 -1.25 1.40 1.0 -Grading1 = 1.3 -1.3 -Cells0 = 15 20 20 15 5 -Cells1 = 20 20 +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 @@ -18,9 +18,9 @@ DarcyStart = 0.0 Positions0 = 0.0 0.5 1.0 Positions1 = 0.0 0.5 Cells0 = 20 20 -Cells1 = 20 +Cells1 = 40 Grading0 = 1.25 -1.25 -Grading1 = -1.3 +Grading1 = -1.15 Verbosity = true [RANS.Problem] @@ -35,8 +35,8 @@ SmoothingSlipVelocity = 0.05 Name = darcy Pressure = 1e5 Temperature = 298.15 # [K] -BaseClosed = false -InletRelPressure = 10 +BaseClosed = true +BottomVelocity = 0.606 # [m/s] [Darcy.SpatialParams] Porosity = 0.75 @@ -49,10 +49,10 @@ ForchCoeff = 2.65e-2 # [Component] GasDensity = 1.2 -GasKinematicViscosity = 1.82e-5 +GasKinematicViscosity = 1.5166e-5 [Problem] -Name = test_rans1pdarcy1pforchheimer +Name = test_rans1pdarcy1pforchheimer_closed EnableGravity = false [RANS.RANS] -- GitLab From 804b740b03de4d52c0f30885777e5c72d8949872 Mon Sep 17 00:00:00 2001 From: Ned Coltman Date: Wed, 14 Nov 2018 09:12:27 +0100 Subject: [PATCH 10/10] added averaged parameters to spatial params --- test/multidomain/boundary/ransforchheimer/spatialparams.hh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/multidomain/boundary/ransforchheimer/spatialparams.hh b/test/multidomain/boundary/ransforchheimer/spatialparams.hh index b007c5706b..af4ba63d06 100644 --- a/test/multidomain/boundary/ransforchheimer/spatialparams.hh +++ b/test/multidomain/boundary/ransforchheimer/spatialparams.hh @@ -90,6 +90,10 @@ public: 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_; } -- GitLab