From 912d07b8c4108777fd6ed3cccf869f0dc04a7ed6 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de> Date: Thu, 9 Sep 2021 15:55:03 +0200 Subject: [PATCH] [test][md] Add freeflow-porousmedium coupling test --- test/multidomain/boundary/CMakeLists.txt | 1 + .../freeflowporousmedium/1p_1p/CMakeLists.txt | 30 ++ .../freeflowporousmedium/1p_1p/main.cc | 213 ++++++++++++ .../1p_1p/params_horizontalflow.input | 33 ++ .../1p_1p/params_verticalflow.input | 37 +++ .../1p_1p/problem_darcy.hh | 212 ++++++++++++ .../1p_1p/problem_freeflow.hh | 302 ++++++++++++++++++ .../freeflowporousmedium/1p_1p/properties.hh | 129 ++++++++ .../1p_1p/spatialparams.hh | 91 ++++++ .../freeflowporousmedium/CMakeLists.txt | 1 + 10 files changed, 1049 insertions(+) create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh create mode 100644 test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh create mode 100644 test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt diff --git a/test/multidomain/boundary/CMakeLists.txt b/test/multidomain/boundary/CMakeLists.txt index addd460feb..5ffc5873b2 100644 --- a/test/multidomain/boundary/CMakeLists.txt +++ b/test/multidomain/boundary/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(darcydarcy) +add_subdirectory(freeflowporousmedium) add_subdirectory(stokesdarcy) diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt new file mode 100644 index 0000000000..da2b6d53d6 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/CMakeLists.txt @@ -0,0 +1,30 @@ +add_input_file_links() + +add_executable(test_md_boundary_ff1p_pm1p EXCLUDE_FROM_ALL main.cc) + +dumux_add_test(NAME test_md_boundary_ff1p_pm1p_horizontal + LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes + TARGET test_md_boundary_ff1p_pm1p + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_horizontal_freeflow-00000.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_horizontal_darcy-00000.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p params_horizontalflow.input + -Vtk.OutputName test_md_boundary_ff1p_pm1p_horizontal") + +dumux_add_test(NAME test_md_boundary_ff1p_pm1p_vertical + LABELS multidomain multidomain_boundary stokesdarcy 1p navierstokes + TARGET test_md_boundary_ff1p_pm1p + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_vertical_freeflow-00000.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p_vertical_darcy-00000.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_ff1p_pm1p params_verticalflow.input + -Vtk.OutputName test_md_boundary_ff1p_pm1p_vertical" + --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17}) diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc b/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc new file mode 100644 index 0000000000..17abbaac06 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/main.cc @@ -0,0 +1,213 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief A test problem for the coupled Stokes/Darcy problem (1p). + */ + +#include <config.h> + +#include <ctime> +#include <iostream> + +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/timer.hh> + +#include <dumux/common/properties.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/dumuxmessage.hh> +#include <dumux/linear/seqsolverbackend.hh> +#include <dumux/assembly/fvassembler.hh> +#include <dumux/assembly/diffmethod.hh> +#include <dumux/discretization/method.hh> +#include <dumux/io/vtkoutputmodule.hh> +#include <dumux/io/grid/gridmanager_yasp.hh> + +#include <dumux/multidomain/fvassembler.hh> +#include <dumux/multidomain/newtonsolver.hh> +#include <dumux/freeflow/navierstokes/velocityoutput.hh> + +#include "properties.hh" + +int main(int argc, char** argv) +{ + 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 FreeFlowMomentumTypeTag = Properties::TTag::FreeFlowOnePMomentum; + using FreeFlowMassTypeTag = Properties::TTag::FreeFlowOnePMass; + using DarcyTypeTag = Properties::TTag::DarcyOneP; + + // try to create a grid (from the given grid file or the input file) + // for both sub-domains + using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>; + DarcyGridManager darcyGridManager; + darcyGridManager.init("Darcy"); // pass parameter group + + using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowMomentumTypeTag, Properties::Grid>>; + FreeFlowGridManager freeFlowGridManager; + freeFlowGridManager.init("FreeFlow"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using FreeFlowMomentumGridGeometry = GetPropType<FreeFlowMomentumTypeTag, Properties::GridGeometry>; + auto freeFlowMomentumGridGeometry = std::make_shared<FreeFlowMomentumGridGeometry>(freeFlowGridView); + using FreeFlowMassGridGeometry = GetPropType<FreeFlowMassTypeTag, Properties::GridGeometry>; + auto freeFlowMassGridGeometry = std::make_shared<FreeFlowMassGridGeometry>(freeFlowGridView); + using DarcyGridGeometry = GetPropType<DarcyTypeTag, Properties::GridGeometry>; + auto darcyGridGeometry = std::make_shared<DarcyGridGeometry>(darcyGridView); + + using Traits = MultiDomainTraits<FreeFlowMomentumTypeTag, FreeFlowMassTypeTag, DarcyTypeTag>; + using CouplingManager = FreeFlowPorousMediumCouplingManager<Traits>; + auto couplingManager = std::make_shared<CouplingManager>(); + + // the indices + constexpr auto freeFlowMomentumIndex = CouplingManager::freeFlowMomentumIndex; + constexpr auto freeFlowMassIndex = CouplingManager::freeFlowMassIndex; + constexpr auto porousMediumIndex = CouplingManager::porousMediumIndex; + + // the problem (initial and boundary conditions) + using FreeFlowMomentumProblem = GetPropType<FreeFlowMomentumTypeTag, Properties::Problem>; + auto freeFlowMomentumProblem = std::make_shared<FreeFlowMomentumProblem>(freeFlowMomentumGridGeometry, couplingManager); + using FreeFlowMassProblem = GetPropType<FreeFlowMassTypeTag, Properties::Problem>; + auto freeFlowMassProblem = std::make_shared<FreeFlowMassProblem>(freeFlowMassGridGeometry, couplingManager); + using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>; + auto darcyProblem = std::make_shared<DarcyProblem>(darcyGridGeometry, couplingManager); + + // the solution vector + Traits::SolutionVector sol; + sol[freeFlowMomentumIndex].resize(freeFlowMomentumGridGeometry->numDofs()); + sol[freeFlowMassIndex].resize(freeFlowMassGridGeometry->numDofs()); + sol[porousMediumIndex].resize(darcyGridGeometry->numDofs()); + + // the grid variables + using FreeFlowMomentumGridVariables = GetPropType<FreeFlowMomentumTypeTag, Properties::GridVariables>; + auto freeFlowMomentumGridVariables = std::make_shared<FreeFlowMomentumGridVariables>(freeFlowMomentumProblem, freeFlowMomentumGridGeometry); + using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>; + using FreeFlowMassGridVariables = GetPropType<FreeFlowMassTypeTag, Properties::GridVariables>; + auto freeFlowMassGridVariables = std::make_shared<FreeFlowMassGridVariables>(freeFlowMassProblem, freeFlowMassGridGeometry); + using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>; + auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyGridGeometry); + + couplingManager->init(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem, + std::make_tuple(freeFlowMomentumGridVariables, freeFlowMassGridVariables, darcyGridVariables), + sol); + + freeFlowMomentumGridVariables->init(sol[freeFlowMomentumIndex]); + freeFlowMassGridVariables->init(sol[freeFlowMassIndex]); + darcyGridVariables->init(sol[porousMediumIndex]); + + // initialize the vtk output module + VtkOutputModule vtkWriterFF(*freeFlowMassGridVariables, sol[freeFlowMassIndex], freeFlowMassProblem->name()); + GetPropType<FreeFlowMassTypeTag, Properties::IOFields>::initOutputModule(vtkWriterFF); // Add model specific output fields + vtkWriterFF.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<FreeFlowMassGridVariables>>()); + + const bool addAnalyticSolutionForV = getParamFromGroup<bool>(freeFlowMassProblem->paramGroup(), "Problem.AddAnalyticSolutionForV", false); + + if (!freeFlowMassProblem->verticalFlow() && addAnalyticSolutionForV) + { + using Scalar = typename Traits::Scalar; + using std::sqrt; + std::vector<Scalar> analyticalVelocityX(freeFlowMassGridGeometry->gridView().size(0)); + const Scalar dPdX = -freeFlowMomentumProblem->pressureDifference() / (freeFlowMassGridGeometry->bBoxMax()[0] - freeFlowMassGridGeometry->bBoxMin()[0]); + static const Scalar mu = getParam<Scalar>("Component.LiquidKinematicViscosity") * getParam<Scalar>("Component.LiquidDensity"); + static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph"); + static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability"); + static const Scalar sqrtK = sqrt(K); + const Scalar sigma = (freeFlowMassGridGeometry->bBoxMax()[1] - freeFlowMassGridGeometry->bBoxMin()[1])/sqrtK; + + const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX; + + for (const auto& element : elements(freeFlowMassGridGeometry->gridView())) + { + const auto eIdx = freeFlowMassGridGeometry->gridView().indexSet().index(element); + const Scalar y = element.geometry().center()[1] - freeFlowMassGridGeometry->bBoxMin()[1]; + + const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX; + analyticalVelocityX[eIdx] = u; + } + + vtkWriterFF.addField(analyticalVelocityX, "analyticalV_x"); + } + + VtkOutputModule pmVtkWriter(*darcyGridVariables, sol[porousMediumIndex], darcyProblem->name()); + GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(pmVtkWriter); + pmVtkWriter.addVelocityOutput(std::make_shared<GetPropType<DarcyTypeTag, Properties::VelocityOutput>>(*darcyGridVariables)); + + // the assembler for a stationary problem + using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(std::make_tuple(freeFlowMomentumProblem, freeFlowMassProblem, darcyProblem), + std::make_tuple(freeFlowMomentumGridGeometry, + freeFlowMassGridGeometry, + darcyGridGeometry), + std::make_tuple(freeFlowMomentumGridVariables, + freeFlowMassGridVariables, + darcyGridVariables), + couplingManager); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>; + NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); + + const auto dofsVS = freeFlowMomentumGridGeometry->numDofs(); + const auto dofsPS = freeFlowMassGridGeometry->numDofs(); + const auto dofsPD = darcyGridGeometry->numDofs(); + std::cout << "Stokes velocity dofs: " << dofsVS << std::endl; + std::cout << "Stokes pressure dofs: " << dofsPS << std::endl; + std::cout << "Darcy pressure dofs: " << dofsPD << std::endl; + std::cout << "Total number of dofs: " << dofsVS + dofsPS + dofsPD << std::endl; + + // solve the non-linear system + nonLinearSolver.solve(sol); + + // write vtk output + vtkWriterFF.write(1.0); + pmVtkWriter.write(1.0); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input new file mode 100644 index 0000000000..3da8b5649d --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_horizontalflow.input @@ -0,0 +1,33 @@ +[Darcy.Grid] +UpperRight = 1 1 +Cells = 20 20 + +[FreeFlow.Grid] +LowerLeft = 0 1 +UpperRight = 1 2 +Cells = 20 20 + +[FreeFlow.Problem] +Name = freeflow +PressureDifference = 1e-9 + +[Darcy.Problem] +Name = darcy + +[Darcy.SpatialParams] +Permeability = 1e-6 # m^2 +AlphaBeaversJoseph = 1.0 + +[Component] +LiquidKinematicViscosity = 1e-6 +LiquidDensity = 1e3 + +[Vtk] +OutputName = test_md_boundary_freeflow1p_darcy1p_horizontal + +[Problem] +EnableGravity = false +EnableInertiaTerms = false + +[Vtk] +AddVelocity = 1 \ No newline at end of file diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input new file mode 100644 index 0000000000..54358b805b --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/params_verticalflow.input @@ -0,0 +1,37 @@ +[Darcy.Grid] +UpperRight = 2.0 2.0 +Cells = 20 20 + +[FreeFlow.Grid] +LowerLeft = 0.0 2.0 +UpperRight = 2.0 4.0 +Cells = 20 20 + +[FreeFlow.Problem] +Name = freeflow +Velocity = -1e-6 + +[Darcy.Problem] +Name = darcy +Pressure = 0.0 + +[Darcy.SpatialParams] +Permeability = 1e-10 # m^2 +AlphaBeaversJoseph = 1.0 + +[Problem] +EnableGravity = false +EnableInertiaTerms = false + +[Component] +LiquidKinematicViscosity = 1e-6 +LiquidDensity = 1e3 + +[Vtk] +OutputName = test_md_boundary_freeflow1p_darcy1p_vertical + +[Vtk] +AddVelocity = 1 + +[Assembly.NumericDifference] +BaseEpsilon = 1e-6 diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh new file mode 100644 index 0000000000..add1c15b12 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_darcy.hh @@ -0,0 +1,212 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief A simple Darcy test problem (cell-centered finite volume method). + */ + +#ifndef DUMUX_DARCY_SUBPROBLEM_HH +#define DUMUX_DARCY_SUBPROBLEM_HH + +#include <dumux/common/boundarytypes.hh> +#include <dumux/common/numeqvector.hh> +#include <dumux/common/properties.hh> +#include <dumux/common/parameters.hh> + +#include <dumux/porousmediumflow/problem.hh> + +namespace Dumux { + +template <class TypeTag> +class DarcySubProblem : public PorousMediumFlowProblem<TypeTag> +{ + using ParentType = PorousMediumFlowProblem<TypeTag>; + using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; + using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; + using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; + using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; + using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + + using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; + +public: + DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry, + std::shared_ptr<CouplingManager> couplingManager) + : ParentType(gridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) + { + problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name"); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; + } + + /*! + * \brief The problem name. + */ + const std::string& name() const + { return problemName_; } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the temperature within the domain in [K]. + * + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + // \} + + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary control volume. + * + * \param element The element + * \param scvf The boundary sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const + { + BoundaryTypes values; + values.setAllNeumann(); + + if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf)) + values.setAllCouplingNeumann(); + + if (verticalFlow_) + { + if (onLowerBoundary_(scvf.center())) + values.setAllDirichlet(); + } + + return values; + } + + /*! + * \brief Evaluates 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 + { return initial(element); } + + /*! + * \brief Evaluates the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeometry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFluxVarsCache Flux variables caches for all faces in stencil + * \param scvf The boundary sub control volume face + * + * For this method, the \a values variable stores primary variables. + */ + template<class ElementVolumeVariables, class ElementFluxVarsCache> + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFluxVarsCache& elemFluxVarsCache, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if (couplingManager().isCoupled(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, scvf)) + values[Indices::conti0EqIdx] = couplingManager().massCouplingCondition(CouplingManager::porousMediumIndex, CouplingManager::freeFlowMassIndex, + fvGeometry, scvf, elemVolVars); + + return values; + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + /*! + * \brief Evaluates the source term for all phases within a given + * sub control volume. + * + * \param element The element for which the source term is set + * \param fvGeometry The fvGeometry + * \param elemVolVars The element volume variables + * \param scv The sub control volume + */ + template<class ElementVolumeVariables> + NumEqVector source(const Element &element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const + { return NumEqVector(0.0); } + + // \} + + /*! + * \brief Evaluates 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 + { + return PrimaryVariables(0.0); + } + + // \} + + //! Get the coupling manager + const CouplingManager& couplingManager() const + { return *couplingManager_; } + +private: + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; } + + Scalar eps_; + std::shared_ptr<CouplingManager> couplingManager_; + std::string problemName_; + bool verticalFlow_; +}; + +} // end namespace Dumux + +#endif diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh new file mode 100644 index 0000000000..4feebadeed --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh @@ -0,0 +1,302 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup NavierStokesTests + * \brief Test for the staggered grid Navier-Stokes model with analytical solution (Kovasznay 1948, \cite Kovasznay1948) + */ + +#ifndef DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH +#define DUMUX_KOVASZNAY_TEST_PROBLEM_NEW_HH + + +#include <dumux/freeflow/navierstokes/problem.hh> +#include <dumux/common/properties.hh> +#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh> +#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh> + +namespace Dumux { + +/*! + * \ingroup NavierStokesTests + * \brief Test problem for the staggered grid (Kovasznay 1948, \cite Kovasznay1948) + * + * A two-dimensional Navier-Stokes flow with a periodicity in one direction + * is considered. The set-up represents a wake behind a two-dimensional grid + * and is chosen in a way such that an exact solution is available. + */ +template <class TypeTag> +class FreeFlowOnePTestProblem : public NavierStokesProblem<TypeTag> +{ + using ParentType = NavierStokesProblem<TypeTag>; + + using BoundaryTypes = typename ParentType::BoundaryTypes; + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; + using NumEqVector = typename ParentType::NumEqVector; + using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; + using PrimaryVariables = typename ParentType::PrimaryVariables; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; + + static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using VelocityVector = Dune::FieldVector<Scalar, dimWorld>; + + using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; + + // static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>(); + +public: + FreeFlowOnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager) + : ParentType(gridGeometry, couplingManager, "FreeFlow") + , couplingManager_(couplingManager) + { + problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name"); + deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference", 0.0); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief The problem name. + */ + const std::string& name() const + { + return problemName_; + } + + /*! + * \brief Returns the temperature within the domain in [K]. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperature() const + { return 298.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 values; + const auto& globalPos = scvf.center(); //avoid ambiguities at corners + + if (verticalFlow_) + { + if constexpr (ParentType::isMomentumProblem()) + { + if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)) + { + values.setCouplingNeumann(Indices::momentumYBalanceIdx); + values.setCouplingNeumann(Indices::momentumXBalanceIdx); + } + else + values.setAllDirichlet(); + } + else + { + if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf)) + values.setAllCouplingNeumann(); + else + values.setAllNeumann(); + } + } + else // horizontal flow + { + if constexpr (ParentType::isMomentumProblem()) + { + if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) + values.setAllNeumann(); + else if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)) + { + values.setCouplingNeumann(Indices::momentumYBalanceIdx); + values.setCouplingNeumann(Indices::momentumXBalanceIdx); + } + else + values.setAllDirichlet(); + } + else + { + if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf)) + values.setAllCouplingNeumann(); + else + values.setAllNeumann(); + } + } + + return values; + } + + /*! + * \brief Returns Dirichlet boundary values at a given position. + * + * \param globalPos The global position + */ + PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const + { + PrimaryVariables values(0.0); + if constexpr (ParentType::isMomentumProblem()) + { + if (verticalFlow_) + { + static const Scalar inletVelocity = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity"); + values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]); + } + } + return values; + } + + /*! + * \brief Evaluates the boundary conditions for a Neumann control volume. + * + * \param element The element for which the Neumann boundary condition is set + * \param fvGeometry The fvGeometry + * \param elemVolVars The element volume variables + * \param elemFaceVars The element face variables + * \param scvf The boundary sub control volume face + */ + template<class ElementVolumeVariables, class ElementFluxVariablesCache> + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFluxVariablesCache& elemFluxVarsCache, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + const auto& globalPos = scvf.ipGlobal(); + using FluxHelper = NavierStokesMomentumBoundaryFluxHelper; + + if constexpr (ParentType::isMomentumProblem()) + { + + if (couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf)) + { + values += couplingManager_->momentumCouplingCondition(CouplingManager::freeFlowMomentumIndex, + CouplingManager::porousMediumIndex, + fvGeometry, scvf, elemVolVars); + + values += FluxHelper::slipVelocityMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache); + } + else + { + const Scalar pressure = onLeftBoundary_(globalPos) ? deltaP_ : 0.0; + values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, pressure, true /*zeroNormalVelocityGradient*/); + } + } + else + { + if (couplingManager_->isCoupled(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, scvf)) + { + values = couplingManager_->massCouplingCondition(CouplingManager::freeFlowMassIndex, CouplingManager::porousMediumIndex, + fvGeometry, scvf, elemVolVars); + } + else + { + using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>; + values = FluxHelper::scalarOutflowFlux(*this, element, fvGeometry, scvf, elemVolVars); + } + } + + return values; + } + + /*! + * \brief Returns true if the scvf lies on a porous slip boundary + */ + bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const + { + assert(scvf.isFrontal()); + return scvf.boundary() && couplingManager_->isCoupled(CouplingManager::freeFlowMomentumIndex, CouplingManager::porousMediumIndex, scvf); + } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const + { return couplingManager_->darcyPermeability(fvGeometry, scvf); } + + /*! + * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition + */ + Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const + { return couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); } + + /*! + * \brief Returns the pressure difference between inlet and outlet for the horizontal flow scenario + */ + Scalar pressureDifference() const + { + assert(!verticalFlow_); + return deltaP_; + } + + /*! + * \brief Returns true if a vertical flow scenario is considered + */ + bool verticalFlow() const + { return verticalFlow_; } + + // \} + +private: + + bool onLeftBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; } + + bool onRightBoundary_(const GlobalPosition &globalPos) const + { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } + + bool onLowerBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; } + + bool onUpperBoundary_(const GlobalPosition &globalPos) const + { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; } + + std::string problemName_; + static constexpr Scalar eps_ = 1e-6; + Scalar deltaP_; + bool verticalFlow_; + + std::shared_ptr<CouplingManager> couplingManager_; +}; +} // end namespace Dumux + +#endif diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh new file mode 100644 index 0000000000..c60ecddd51 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/properties.hh @@ -0,0 +1,129 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \brief The properties for a simple Darcy test (cell-centered finite volume method) + */ +#ifndef DUMUX_DARCYSTOKES_PROPERTIES_HH +#define DUMUX_DARCYSTOKES_PROPERTIES_HH + +#include <dune/grid/yaspgrid.hh> + +#include <dumux/material/fluidsystems/1pliquid.hh> +#include <dumux/material/components/constant.hh> + +#include <dumux/porousmediumflow/1p/model.hh> +#include <dumux/freeflow/navierstokes/mass/1p/model.hh> +#include <dumux/freeflow/navierstokes/momentum/model.hh> +#include <dumux/discretization/cctpfa.hh> +#include <dumux/discretization/fcstaggered.hh> + +#include <dumux/multidomain/boundary/freeflowporousmedium/couplingmanager.hh> + +#include "spatialparams.hh" +#include "problem_darcy.hh" +#include "problem_freeflow.hh" + +namespace Dumux::Properties { + +// Create new type tags +namespace TTag { +struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; }; +} // end namespace TTag + +// Set the problem property +template<class TypeTag> +struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; }; + +// the fluid system +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::DarcyOneP> +{ + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ; +}; + +// Set the grid type +template<class TypeTag> +struct Grid<TypeTag, TTag::DarcyOneP> { using type = Dune::YaspGrid<2>; }; + +template<class TypeTag> +struct SpatialParams<TypeTag, TTag::DarcyOneP> +{ + using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using type = OnePSpatialParams<GridGeometry, Scalar>; +}; + +// Create new type tags +namespace TTag { +struct FreeFlowOneP {}; +struct FreeFlowOnePMass { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMassOneP, CCTpfaModel>; }; +struct FreeFlowOnePMomentum { using InheritsFrom = std::tuple<FreeFlowOneP, NavierStokesMomentum, FaceCenteredStaggeredModel>; }; +} // end namespace TTag + +// the fluid system +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::FreeFlowOneP> +{ + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ; +}; + +// Set the grid type +template<class TypeTag> +struct Grid<TypeTag, TTag::FreeFlowOneP> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; }; + +// Set the problem property +template<class TypeTag> +struct Problem<TypeTag, TTag::FreeFlowOneP> { using type = Dumux::FreeFlowOnePTestProblem<TypeTag> ; }; + +template<class TypeTag> +struct EnableGridGeometryCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowOneP> { static constexpr bool value = true; }; + +template<class TypeTag> +struct CouplingManager<TypeTag, TTag::DarcyOneP> +{ + using Traits = MultiDomainTraits<Properties::TTag::FreeFlowOnePMomentum, Properties::TTag::FreeFlowOnePMass, TTag::DarcyOneP>; + using type = Dumux::FreeFlowPorousMediumCouplingManager<Traits>; +}; + +template<class TypeTag> +struct CouplingManager<TypeTag, TTag::FreeFlowOneP> +{ + using Traits = MultiDomainTraits<Properties::TTag::FreeFlowOnePMomentum, Properties::TTag::FreeFlowOnePMass, TTag::DarcyOneP>; + using type = Dumux::FreeFlowPorousMediumCouplingManager<Traits>; +}; + +// template<class TypeTag> +// struct CouplingManager<TypeTag, TTag::StokesOneP> +// { +// using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>; +// using type = Dumux::StokesDarcyCouplingManager<Traits>; +// }; + + +} // end namespace Dumux::Properties + +#endif diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh new file mode 100644 index 0000000000..54339aae00 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/spatialparams.hh @@ -0,0 +1,91 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup BoundaryTests + * \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 <dumux/material/spatialparams/fv1p.hh> + +namespace Dumux { + +/*! + * \ingroup BoundaryTests + * \brief The spatial parameters class for the test problem using the + * 1p cc model. + */ +template<class GridGeometry, class Scalar> +class OnePSpatialParams +: public FVSpatialParamsOneP<GridGeometry, Scalar, + OnePSpatialParams<GridGeometry, Scalar>> +{ + using GridView = typename GridGeometry::GridView; + using ParentType = FVSpatialParamsOneP<GridGeometry, Scalar, + OnePSpatialParams<GridGeometry, Scalar>>; + + 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<const GridGeometry> gridGeometry) + : ParentType(gridGeometry) + { + permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability"); + alphaBJ_ = getParam<Scalar>("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 Defines the porosity in [-]. + * + * \param globalPos The global position + */ + Scalar porosityAtPos(const GlobalPosition& globalPos) const + { return 0.4; } + + /*! \brief Defines the Beavers-Joseph coefficient in [-]. + * + * \param globalPos The global position + */ + Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const + { return alphaBJ_; } + + +private: + Scalar permeability_; + Scalar alphaBJ_; +}; + +} // end namespace Dumux + +#endif diff --git a/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt b/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt new file mode 100644 index 0000000000..f62f49df14 --- /dev/null +++ b/test/multidomain/boundary/freeflowporousmedium/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(1p_1p) -- GitLab