diff --git a/appl/coupling-ff-pm/CMakeLists.txt b/appl/coupling-ff-pm/CMakeLists.txt index d3e96b83cfd100be34c5e8da7806d8089e6de69c..943da4e7824478d66bc07cf1140a887e2f6cfc13 100644 --- a/appl/coupling-ff-pm/CMakeLists.txt +++ b/appl/coupling-ff-pm/CMakeLists.txt @@ -3,5 +3,6 @@ add_custom_target(test_exercises) add_subdirectory(monolithic) +add_subdirectory(fvca-monolithic-reversed) add_subdirectory(iterative) add_subdirectory(iterative-reversed) diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/CMakeLists.txt b/appl/coupling-ff-pm/fvca-monolithic-reversed/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..685542c73737916e8022a5971f5c24867c20c81a --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/CMakeLists.txt @@ -0,0 +1,4 @@ + +add_input_file_links() + +add_executable(fvca-monolithic-reversed EXCLUDE_FROM_ALL main.cc) diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/main.cc b/appl/coupling-ff-pm/fvca-monolithic-reversed/main.cc new file mode 100644 index 0000000000000000000000000000000000000000..8cfe6fc893e773c195b1b60ce80a3a449b5da426 --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/main.cc @@ -0,0 +1,274 @@ +// -*- 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 <dune/istl/io.hh> + +#include <dumux/common/properties.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/partial.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/staggeredvtkoutputmodule.hh> +#include <dumux/io/grid/gridmanager.hh> + +#include <dumux/multidomain/staggeredtraits.hh> +#include <dumux/multidomain/fvassembler.hh> +#include <dumux/multidomain/newtonsolver.hh> + +#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh> + +#include "problem_darcy.hh" +#include "problem_stokes.hh" + +namespace Dumux { +namespace Properties { + +template<class TypeTag> +struct CouplingManager<TypeTag, TTag::StokesOneP> +{ + using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>; + using type = Dumux::StokesDarcyCouplingManager<Traits>; +}; + +template<class TypeTag> +struct CouplingManager<TypeTag, TTag::DarcyOneP> +{ + using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>; + using type = Dumux::StokesDarcyCouplingManager<Traits>; +}; + +} // end namespace Properties +} // end namespace Dumux + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv); + + // Define the sub problem type tags + using StokesTypeTag = Properties::TTag::StokesOneP; + using DarcyTypeTag = Properties::TTag::DarcyOneP; + + // try to create a grid (from the given grid file or the input file) + // for both sub-domains + using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>; + DarcyGridManager darcyGridManager; + darcyGridManager.init("Darcy"); // pass parameter group + + using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>; + StokesGridManager stokesGridManager; + stokesGridManager.init("Stokes"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + const auto& stokesGridView = stokesGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>; + auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView); + stokesFvGridGeometry->update(); + using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>; + auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView); + darcyFvGridGeometry->update(); + + using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>; + + // the coupling manager + using CouplingManager = StokesDarcyCouplingManager<Traits>; + auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry); + + // the indices + constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx; + constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx; + constexpr auto darcyIdx = CouplingManager::darcyIdx; + + // the problem (initial and boundary conditions) + using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>; + auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager); + using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>; + auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager); + + // the solution vector + Traits::SolutionVector sol; + sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs()); + sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); + sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); + + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); + + couplingManager->init(stokesProblem, darcyProblem, sol); + + // the grid variables + using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>; + auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry); + stokesGridVariables->init(stokesSol); + using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>; + auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry); + darcyGridVariables->init(sol[darcyIdx]); + + couplingManager->setGridVariables(std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(), + stokesGridVariables->faceGridVariablesPtr(), + darcyGridVariables)); + + // intialize the vtk output module + StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter); + stokesVtkWriter.write(0.0); + + VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name()); + using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>; + darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables)); + GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter); + darcyVtkWriter.write(0.0); + + // the assembler for a stationary problem + using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem), + std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(), + stokesFvGridGeometry->faceFVGridGeometryPtr(), + darcyFvGridGeometry), + std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(), + stokesGridVariables->faceGridVariablesPtr(), + 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); + + // solve the non-linear system + nonLinearSolver.solve(sol); + + // write vtk output + stokesVtkWriter.write(1.0); + darcyVtkWriter.write(1.0); + + for (const auto& element : elements(darcyGridView)) + { + auto fvGeometry = localView(*darcyFvGridGeometry); + fvGeometry.bind(element); + + for (const auto& scvf : scvfs(fvGeometry)) + { + if (couplingManager->isCoupledEntity(CouplingManager::darcyIdx, scvf)) + { + std::ostream tmp(std::cout.rdbuf()); + + tmp << std::scientific << "Interface pressure ff at " << scvf.center() << " is " << couplingManager->couplingData().freeFlowInterfacePressure(element, scvf) << std::endl; + } + } + } + + for (const auto& element : elements(stokesGridView)) + { + auto fvGeometry = localView(*stokesFvGridGeometry); + auto elemVolVars = localView(stokesGridVariables->cellCenterGridVariablesPtr()->curGridVolVars()); + auto elemFaceVars = localView(stokesGridVariables->faceGridVariablesPtr()->curGridFaceVars()); + // auto elemFluxVarsCache = localView(gridVars_(stokesIdx).gridFluxVarsCache()); + + fvGeometry.bind(element); + elemVolVars.bind(element, fvGeometry, sol[stokesCellCenterIdx]); + elemFaceVars.bind(element, fvGeometry, sol[stokesFaceIdx]); + + + for (const auto& scvf : scvfs(fvGeometry)) + { + if (couplingManager->isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + std::ostream tmp(std::cout.rdbuf()); + tmp << std::scientific << "Interface pressure pm at " << scvf.center() << " is " << stokesProblem->neumann(element, + fvGeometry, + elemVolVars, + elemFaceVars, + scvf)[scvf.directionIndex()] << ", vself " << elemFaceVars[scvf].velocitySelf() << ", vDarcy " << couplingManager->couplingData().darcyInterfaceVelocity(element, + fvGeometry, + elemVolVars, + elemFaceVars, + scvf) << std::endl; + } + } + } + + //////////////////////////////////////////////////////////// + // 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/appl/coupling-ff-pm/fvca-monolithic-reversed/params_horizontalflow.input b/appl/coupling-ff-pm/fvca-monolithic-reversed/params_horizontalflow.input new file mode 100644 index 0000000000000000000000000000000000000000..efc0963b125cf5c7583fa12eaa8fc58190f45d6f --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/params_horizontalflow.input @@ -0,0 +1,29 @@ +[Darcy.Grid] +UpperRight = 1 1 +Cells = 20 20 + +[Stokes.Grid] +LowerLeft = 0 1 +UpperRight = 1 2 +Cells = 20 20 + +[Stokes.Problem] +Name = stokes +PressureDifference = 1e-9 + +[Darcy.Problem] +Name = darcy + +[Darcy.SpatialParams] +Permeability = 1e-6 # m^2 +AlphaBeaversJoseph = 1.0 + +[Vtk] +OutputName = test_md_boundary_stokes1p_darcy1p_horizontal + +[Problem] +EnableGravity = false +EnableInertiaTerms = false + +[Vtk] +AddVelocity = 1 diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/params_verticalflow.input b/appl/coupling-ff-pm/fvca-monolithic-reversed/params_verticalflow.input new file mode 100644 index 0000000000000000000000000000000000000000..b5f39b858a7d2e31a64d74636487edf39f53b74f --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/params_verticalflow.input @@ -0,0 +1,33 @@ +[Darcy.Grid] +UpperRight = 2.0 2.0 +Cells = 20 20 + +[Stokes.Grid] +LowerLeft = 0.0 2.0 +UpperRight = 2.0 4.0 +Cells = 20 20 + +[Stokes.Problem] +Name = stokes +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 + +[Vtk] +OutputName = test_md_boundary_darcy1p_stokes1p_vertical + +[Vtk] +AddVelocity = 1 + +[Assembly.NumericDifference] +BaseEpsilon = 1e-6 diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_darcy.hh b/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_darcy.hh new file mode 100644 index 0000000000000000000000000000000000000000..b0c51e2969429cfc62b025d8ad155aabe21835d3 --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_darcy.hh @@ -0,0 +1,257 @@ +// -*- 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 <dune/grid/yaspgrid.hh> + +#include <dumux/discretization/cctpfa.hh> + +#include <dumux/porousmediumflow/1p/model.hh> +#include <dumux/porousmediumflow/problem.hh> + +#include "spatialparams.hh" + +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/fluidsystems/1pliquid.hh> + +namespace Dumux { +template <class TypeTag> +class DarcySubProblem; + +namespace 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::SimpleH2O<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 FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using type = OnePSpatialParams<FVGridGeometry, Scalar>; +}; +} // end namespace Properties + +template <class TypeTag> +class DarcySubProblem : public PorousMediumFlowProblem<TypeTag> +{ + using ParentType = PorousMediumFlowProblem<TypeTag>; + using GridView = GetPropType<TypeTag, Properties::GridView>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; + using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>; + using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>; + using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; + using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; + + 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 FVGridGeometry> fvGridGeometry, + std::shared_ptr<CouplingManager> couplingManager) + : ParentType(fvGridGeometry, "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().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + values.setAllCouplingDirichlet(); + // 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 + { + PrimaryVariables values = initial(element); + if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + values = couplingManager().couplingData().freeFlowInterfacePressure(element, scvf); + + 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 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().isCoupledEntity(CouplingManager::darcyIdx, scvf)) + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf); + + 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->fvGridGeometry().bBoxMin()[1] + eps_; } + + Scalar eps_; + std::shared_ptr<CouplingManager> couplingManager_; + std::string problemName_; + bool verticalFlow_; +}; +} // end namespace Dumux + +#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_stokes.hh b/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_stokes.hh new file mode 100644 index 0000000000000000000000000000000000000000..d88892f6697893433cf112248cfa51d9a4e4929b --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/problem_stokes.hh @@ -0,0 +1,353 @@ +// -*- 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 Stokes test problem for the staggered grid (Navier-)Stokes model. + */ + +#ifndef DUMUX_STOKES_SUBPROBLEM_HH +#define DUMUX_STOKES_SUBPROBLEM_HH + +#include <dune/grid/yaspgrid.hh> + +#include <dumux/material/fluidsystems/1pliquid.hh> +#include <dumux/material/components/simpleh2o.hh> + +#include <dumux/freeflow/navierstokes/problem.hh> +#include <dumux/discretization/staggered/freeflow/properties.hh> +#include <dumux/freeflow/navierstokes/model.hh> + +namespace Dumux { +template <class TypeTag> +class StokesSubProblem; + +namespace Properties { +// Create new type tags +namespace TTag { +struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; }; +} // end namespace TTag + +// the fluid system +template<class TypeTag> +struct FluidSystem<TypeTag, TTag::StokesOneP> +{ + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ; +}; + +// Set the grid type +template<class TypeTag> +struct Grid<TypeTag, TTag::StokesOneP> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; }; + +// Set the problem property +template<class TypeTag> +struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; }; + +template<class TypeTag> +struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +template<class TypeTag> +struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +} // end namespace Properties + +/*! + * \ingroup BoundaryTests + * \brief Test problem for the one-phase (Navier-) Stokes problem. + * + * Horizontal flow from left to right with a parabolic velocity profile. + */ +template <class TypeTag> +class StokesSubProblem : public NavierStokesProblem<TypeTag> +{ + using ParentType = NavierStokesProblem<TypeTag>; + + using GridView = GetPropType<TypeTag, Properties::GridView>; + using Scalar = GetPropType<TypeTag, Properties::Scalar>; + + using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; + + using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>; + + using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using Element = typename GridView::template Codim<0>::Entity; + + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; + using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>; + + using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; + +public: + StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager) + : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), 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]. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperature() const + { return 273.15 + 10; } // 10°C + + /*! + * \brief Returns the sources within the domain. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition &globalPos) const + { 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 values; + + const auto& globalPos = scvf.dofPosition(); + + if (verticalFlow_) + { + // inflow + if(onUpperBoundary_(globalPos)) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + + // left/right wall + if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos))) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + } + else // horizontal flow + { + // inlet / outlet + if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) + values.setDirichlet(Indices::pressureIdx); + else // wall + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + + } + + if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + values.setCouplingDirichlet(Indices::velocityYIdx); + // values.setCouplingNeumann(Indices::conti0EqIdx); + // values.setCouplingNeumann(Indices::momentumYBalanceIdx); + values.setBJS(Indices::momentumXBalanceIdx); + } + + return values; + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet control volume. + */ + PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const + { + return initialAtPos(globalPos); + } + + // /*! + // * \brief Evaluates the boundary conditions for a Neumann control volume. + // * + // * \param element The element for which the Neumann boundary condition is set + // * \param fvGeometry The fvGeometry + // * \param elemVolVars The element volume variables + // * \param elemFaceVars The element face variables + // * \param scvf The boundary sub control volume face + // */ + // template<class ElementVolumeVariables, class ElementFaceVariables> + // NumEqVector dirichlet(const Element& element, + // const FVElementGeometry& fvGeometry, + // const ElementVolumeVariables& elemVolVars, + // const ElementFaceVariables& elemFaceVars, + // const SubControlVolumeFace& scvf) const + // { + // NumEqVector values(0.0); + // + // if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + // { + // + // } + // } + + /*! + * \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 ElementFaceVariables> + NumEqVector neumann(const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const ElementFaceVariables& elemFaceVars, + const SubControlVolumeFace& scvf) const + { + NumEqVector values(0.0); + + if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) + { + // std::cout << "vel is " << elemFaceVars[scvf].velocitySelf() << std::endl; + values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); + values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); + + // std::cout << "v self is " << elemFaceVars[scvf].velocitySelf() << ", v darcy " << couplingManager().couplingData().darcyInterfaceVelocity(element, fvGeometry, elemVolVars, elemFaceVars, scvf) << std::endl; + + // + // using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; + // FluxVariables fluxVars; + // + // const auto& gridFluxVarsCache = couplingManager().gridVars_(CouplingManager::stokesIdx).gridFluxVarsCache(); + // + // std::cout << "p is " << fluxVars.computeMomentumFlux(*this, + // element, + // scvf, + // fvGeometry, + // elemVolVars, + // elemFaceVars, + // gridFluxVarsCache) << std::endl; + + + } + return values; + } + + // \} + + //! Get the coupling manager + const CouplingManager& couplingManager() const + { return *couplingManager_; } + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluates the initial value for a control volume. + * + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const + { + PrimaryVariables values(0.0); + + if (verticalFlow_) + { + static const Scalar inletVelocity = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity"); + values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]); + } + else // horizontal flow + { + static const Scalar deltaP = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference"); + if(onLeftBoundary_(globalPos)) + values[Indices::pressureIdx] = deltaP; + } + + return values; + } + + /*! + * \brief Returns the intrinsic permeability of required as input parameter + for the Beavers-Joseph-Saffman boundary condition + */ + Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const + { + return couplingManager().couplingData().darcyPermeability(element, scvf); + } + + /*! + * \brief Returns the alpha value required as input parameter for the + Beavers-Joseph-Saffman boundary condition. + */ + Scalar alphaBJ(const SubControlVolumeFace& scvf) const + { + return couplingManager().problem(CouplingManager::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_; } + + Scalar eps_; + std::string problemName_; + bool verticalFlow_; + + std::shared_ptr<CouplingManager> couplingManager_; +}; +} // end namespace Dumux + +#endif // DUMUX_STOKES_SUBPROBLEM_HH diff --git a/appl/coupling-ff-pm/fvca-monolithic-reversed/spatialparams.hh b/appl/coupling-ff-pm/fvca-monolithic-reversed/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..fb0adabe6240371a406469339a2099ca259eff03 --- /dev/null +++ b/appl/coupling-ff-pm/fvca-monolithic-reversed/spatialparams.hh @@ -0,0 +1,92 @@ +// -*- 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 FVGridGeometry, class Scalar> +class OnePSpatialParams +: public FVSpatialParamsOneP<FVGridGeometry, Scalar, + OnePSpatialParams<FVGridGeometry, Scalar>> +{ + using GridView = typename FVGridGeometry::GridView; + using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, + OnePSpatialParams<FVGridGeometry, 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 FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry) + { + 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/appl/coupling-ff-pm/iterative-reversed/precice-config-parallel-implicit-reversed.xml b/appl/coupling-ff-pm/iterative-reversed/precice-config-parallel-implicit-reversed.xml new file mode 100644 index 0000000000000000000000000000000000000000..b0ff6c0d2b15e16dc691c1c9b79816cfa9997aee --- /dev/null +++ b/appl/coupling-ff-pm/iterative-reversed/precice-config-parallel-implicit-reversed.xml @@ -0,0 +1,79 @@ +<?xml version="1.0"?> + +<precice-configuration> + <log> + <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" /> + <sink type="stream" output="stdout" enabled="false" /> + </log> + + <solver-interface dimensions="2"> + + <data:scalar name="Velocity"/> + <data:scalar name="Pressure"/> + + <mesh name="FreeFlowMesh"> + <use-data name="Velocity" /> + <use-data name="Pressure" /> + </mesh> + + <mesh name="DarcyMesh"> + <use-data name="Velocity" /> + <use-data name="Pressure" /> + </mesh> + + <participant name="FreeFlow"> + <use-mesh name="FreeFlowMesh" provide="yes"/> + <use-mesh name="DarcyMesh" from="Darcy"/> + + <read-data name="Velocity" mesh="FreeFlowMesh"/> + <write-data name="Pressure" mesh="FreeFlowMesh"/> + + <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/> + <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/> + + </participant> + + <participant name="Darcy"> + <use-mesh name="DarcyMesh" provide="yes"/> + + <read-data name="Pressure" mesh="DarcyMesh"/> + <write-data name="Velocity" mesh="DarcyMesh"/> + </participant> + + <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" /> + + + <coupling-scheme:parallel-implicit> + <max-time value="1"/> + <timestep-length value="1" /> + <max-iterations value="100"/> + + + <participants first="FreeFlow" second="Darcy"/> + <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="true" /> + <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="true" /> + + <relative-convergence-measure limit="1.0e-5" data="Pressure" mesh="DarcyMesh"/> + <relative-convergence-measure limit="1.0e-5" data="Velocity" mesh="DarcyMesh"/> + + + <!-- + <relative-convergence-measure limit="1.0e-2" data="Velocity" mesh="FreeFlowMesh"/> + --> + + <extrapolation-order value="0"/> + + <post-processing:IQN-ILS> + <data mesh="DarcyMesh" name="Velocity" /> + <data mesh="DarcyMesh" name="Pressure" /> + <preconditioner type="residual-sum" /> + <initial-relaxation value="0.1" /> + <max-used-iterations value="9" /> + <timesteps-reused value="10" /> + <filter type="QR2" limit="1e-3" /> + </post-processing:IQN-ILS> + + </coupling-scheme:parallel-implicit> + </solver-interface> +</precice-configuration> + diff --git a/appl/coupling-ff-pm/iterative-reversed/precice-config-serial-explicit-reversed.xml b/appl/coupling-ff-pm/iterative-reversed/precice-config-serial-explicit-reversed.xml new file mode 100644 index 0000000000000000000000000000000000000000..d78a37837f8b6f765a5b691c8ca0caf1dbff79b7 --- /dev/null +++ b/appl/coupling-ff-pm/iterative-reversed/precice-config-serial-explicit-reversed.xml @@ -0,0 +1,57 @@ +<?xml version="1.0"?> + +<precice-configuration> + <log> + <sink type="stream" output="stdout" filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="false" /> + <sink type="stream" output="stdout" enabled="false" /> + </log> + + <solver-interface dimensions="2"> + + <data:scalar name="Velocity"/> + <data:scalar name="Pressure"/> + + <mesh name="FreeFlowMesh"> + <use-data name="Velocity" /> + <use-data name="Pressure" /> + </mesh> + + <mesh name="DarcyMesh"> + <use-data name="Velocity" /> + <use-data name="Pressure" /> + </mesh> + + <participant name="FreeFlow"> + <use-mesh name="FreeFlowMesh" provide="yes"/> + <use-mesh name="DarcyMesh" from="Darcy"/> + + <read-data name="Velocity" mesh="FreeFlowMesh"/> + <write-data name="Pressure" mesh="FreeFlowMesh"/> + + <mapping:nearest-neighbor direction="write" from="FreeFlowMesh" to="DarcyMesh" constraint="consistent"/> + <mapping:nearest-neighbor direction="read" from="DarcyMesh" to="FreeFlowMesh" constraint="consistent"/> + + </participant> + + <participant name="Darcy"> + <use-mesh name="DarcyMesh" provide="yes"/> + + <read-data name="Pressure" mesh="DarcyMesh"/> + <write-data name="Velocity" mesh="DarcyMesh"/> + </participant> + + <m2n:sockets from="FreeFlow" to="Darcy" distribution-type="gather-scatter" network="lo" /> + + + <coupling-scheme:serial-explicit> + <max-time value="1"/> + <timestep-length value="1" /> + + <participants first="FreeFlow" second="Darcy"/> + <exchange data="Pressure" mesh="DarcyMesh" from="FreeFlow" to="Darcy" initialize="false" /> + <exchange data="Velocity" mesh="DarcyMesh" from="Darcy" to="FreeFlow" initialize="true" /> + + </coupling-scheme:serial-explicit> + </solver-interface> +</precice-configuration> + diff --git a/appl/coupling-ff-pm/monolithic/ffproblem.hh b/appl/coupling-ff-pm/monolithic/ffproblem.hh index c5b642c8b98d0e527858502b3264d6dd3ba97d69..a52c685050cc4a53c639f03ae6c09d5520b68cfe 100644 --- a/appl/coupling-ff-pm/monolithic/ffproblem.hh +++ b/appl/coupling-ff-pm/monolithic/ffproblem.hh @@ -36,8 +36,10 @@ #include <dumux/discretization/staggered/freeflow/properties.hh> #include <dumux/freeflow/navierstokes/model.hh> +#ifdef ENABLEMONOLITHIC +#else #include "../../precice-adapter/include/preciceadapter.hh" - +#endif namespace Dumux { template <class TypeTag> diff --git a/appl/coupling-ff-pm/monolithic/pmproblem.hh b/appl/coupling-ff-pm/monolithic/pmproblem.hh index 25064770ca158ecd78ba35f5dc7ecfbb82a82c58..ff7d69364d7e7de1fd593c9ec3906914455a0cfe 100644 --- a/appl/coupling-ff-pm/monolithic/pmproblem.hh +++ b/appl/coupling-ff-pm/monolithic/pmproblem.hh @@ -43,7 +43,10 @@ #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/fluidsystems/1pliquid.hh> +#ifdef ENABLEMONOLITHIC +#else #include "../../precice-adapter/include/preciceadapter.hh" +#endif namespace Dumux {