From 3c0baf10269a9cedc2b291bcffc04e369718682b Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 21 Nov 2018 17:52:53 +0100 Subject: [PATCH 1/4] [staggered][problem] Make SolutionVector a template paramter --- dumux/common/staggeredfvproblem.hh | 4 +++- dumux/freeflow/navierstokes/problem.hh | 3 +-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh index 2597a72be5..cd7d68b0c0 100644 --- a/dumux/common/staggeredfvproblem.hh +++ b/dumux/common/staggeredfvproblem.hh @@ -57,7 +57,6 @@ class StaggeredFVProblem : public FVProblem using PrimaryVariables = GetPropType; using FVGridGeometry = GetPropType; - using SolutionVector = GetPropType; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; @@ -162,6 +161,7 @@ public: * \brief Applies the initial solution for all degrees of freedom of the grid. * */ + template void applyInitialSolution(SolutionVector& sol) const { sol[cellCenterIdx].resize(this->fvGridGeometry().numCellCenterDofs()); @@ -192,6 +192,7 @@ public: //! Applys the initial cell center solution + template void applyInitialCellCenterSolution(SolutionVector& sol, const SubControlVolume& scv, const PrimaryVariables& initSol) const @@ -205,6 +206,7 @@ public: } //! Applys the initial face solution + template void applyInitialFaceSolution(SolutionVector& sol, const SubControlVolumeFace& scvf, const PrimaryVariables& initSol) const diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 486e398a1d..3c37603c73 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -74,7 +74,6 @@ class NavierStokesProblem : public NavierStokesParentProblem using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using SolutionVector = GetPropType; using PrimaryVariables = GetPropType; using Indices = typename GetPropType::Indices; @@ -137,7 +136,7 @@ public: { return enableInertiaTerms_; } //! Applys the initial face solution (velocities on the faces). Specialization for staggered grid discretization. - template + template typename std::enable_if::type applyInitialFaceSolution(SolutionVector& sol, const SubControlVolumeFace& scvf, -- GitLab From 07b043d14a0a2fb65ec0a4a86de25783f178b557 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 28 Nov 2018 18:13:49 +0100 Subject: [PATCH 2/4] [test][stokesDarcy] Remove obsolete helper vector, pass sol directly --- .../1p2c_1p2c/diffusionlawcomparison/main.cc | 14 ++++---------- .../boundary/stokesdarcy/1p2c_2p2c/main.cc | 17 ++++++----------- .../boundary/stokesdarcy/1p3c_1p3c/main.cc | 14 ++++---------- .../boundary/stokesdarcy/1p_2p/main.cc | 19 +++++-------------- 4 files changed, 19 insertions(+), 45 deletions(-) diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc index 3d6e2c8de9..e73f2ea4e9 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/diffusionlawcomparison/main.cc @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -162,19 +163,12 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; stokesProblem->applyInitialSolution(stokesSol); - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; auto solOld = sol; @@ -192,7 +186,7 @@ int main(int argc, char** argv) try const auto stokesName = getParam("Problem.Name") + "_" + stokesProblem->name(); const auto darcyName = getParam("Problem.Name") + "_" + darcyProblem->name(); - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc index 5024433d0c..a6208789c2 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/main.cc @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -32,7 +33,9 @@ #include #include +#include #include +#include #include #include #include @@ -148,20 +151,12 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; stokesProblem->applyInitialSolution(stokesSol); - auto solStokesOld = stokesSol; - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; auto solOld = sol; @@ -176,7 +171,7 @@ int main(int argc, char** argv) try darcyGridVariables->init(sol[darcyIdx]); // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); diff --git a/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc b/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc index 6cfa8f9178..bb09aad03a 100644 --- a/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p3c_1p3c/main.cc @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -150,19 +151,12 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; stokesProblem->applyInitialSolution(stokesSol); - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; auto solOld = sol; @@ -180,7 +174,7 @@ int main(int argc, char** argv) try const auto stokesName = getParam("Problem.Name") + "_" + stokesProblem->name(); const auto darcyName = getParam("Problem.Name") + "_" + darcyProblem->name(); - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); diff --git a/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc b/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc index 776b53804c..b0c2d89cd8 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p_2p/main.cc @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -33,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -175,20 +177,9 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; - // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; - stokesProblem->applyInitialSolution(stokesSol); - auto solStokesOld = stokesSol; - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - + stokesProblem->applyInitialSolution(sol); darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; auto solOld = sol; @@ -197,7 +188,7 @@ int main(int argc, char** argv) try // the grid variables using StokesGridVariables = GetPropType; auto stokesGridVariables = std::make_shared(stokesProblem, stokesFvGridGeometry); - stokesGridVariables->init(stokesSol); + stokesGridVariables->init(sol); using DarcyGridVariables = GetPropType; auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry); darcyGridVariables->init(sol[darcyIdx]); @@ -209,7 +200,7 @@ int main(int argc, char** argv) try auto dt = getParam("TimeLoop.DtInitial"); // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, sol, stokesProblem->name()); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); -- GitLab From a1e232e121c54920a913bc471c0bcd972258338c Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 28 Nov 2018 22:08:24 +0100 Subject: [PATCH 3/4] [test][md][stokesdarcy] Clean-up 1p_1p * use only one executable --- .../boundary/stokesdarcy/1p_1p/CMakeLists.txt | 32 +- .../1p_1p/horizontalflow/CMakeLists.txt | 15 - .../1p_1p/{horizontalflow => }/main.cc | 15 +- ...rams.input => params_horizontalflow.input} | 0 ...params.input => params_verticalflow.input} | 0 .../{horizontalflow => }/problem_darcy.hh | 41 +-- .../{horizontalflow => }/problem_stokes.hh | 60 +++- .../1p_1p/verticalflow/CMakeLists.txt | 16 - .../stokesdarcy/1p_1p/verticalflow/main.cc | 233 ------------- .../1p_1p/verticalflow/problem_darcy.hh | 286 ---------------- .../1p_1p/verticalflow/problem_stokes.hh | 307 ------------------ 11 files changed, 89 insertions(+), 916 deletions(-) delete mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/CMakeLists.txt rename test/multidomain/boundary/stokesdarcy/1p_1p/{horizontalflow => }/main.cc (93%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{horizontalflow/params.input => params_horizontalflow.input} (100%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{verticalflow/params.input => params_verticalflow.input} (100%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{horizontalflow => }/problem_darcy.hh (90%) rename test/multidomain/boundary/stokesdarcy/1p_1p/{horizontalflow => }/problem_stokes.hh (84%) delete mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/CMakeLists.txt delete mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/main.cc delete mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_darcy.hh delete mode 100644 test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_stokes.hh diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt index 44146a2819..021f50dad7 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/CMakeLists.txt @@ -1,2 +1,30 @@ -add_subdirectory("verticalflow") -add_subdirectory("horizontalflow") +add_input_file_links() + +add_executable(test_md_boundary_darcy1p_stokes1p EXCLUDE_FROM_ALL main.cc) + +dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_horizontal + LABELS multidomain freeflow + TARGET test_md_boundary_darcy1p_stokes1p + 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_darcy1p_stokes1p_horizontal_stokes-00001.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-00001.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p params_horizontalflow.input + -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_horizontal") + +dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_vertical + LABELS multidomain freeflow + TARGET test_md_boundary_darcy1p_stokes1p + 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_darcy1p_stokes1p_vertical_stokes-00001.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_vertical_darcy-00001.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p params_verticalflow.input + -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_vertical" + --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17}) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/CMakeLists.txt deleted file mode 100644 index 48bfbeca35..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -add_input_file_links() - -dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_horizontal - LABELS multidomain freeflow - SOURCES main.cc - 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_darcy1p_stokes1p_horizontal_stokes-00001.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_horizontal_darcy-00001.vtu - - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_horizontal params.input - -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_horizontal") diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc similarity index 93% rename from test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/main.cc rename to test/multidomain/boundary/stokesdarcy/1p_1p/main.cc index f85c72ca83..8d4b638819 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/main.cc @@ -32,6 +32,7 @@ #include #include +#include #include #include #include @@ -133,16 +134,8 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; - - // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; - stokesProblem->applyInitialSolution(stokesSol); - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); couplingManager->init(stokesProblem, darcyProblem, sol); @@ -155,7 +148,7 @@ int main(int argc, char** argv) try darcyGridVariables->init(sol[darcyIdx]); // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/params.input b/test/multidomain/boundary/stokesdarcy/1p_1p/params_horizontalflow.input similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/params.input rename to test/multidomain/boundary/stokesdarcy/1p_1p/params_horizontalflow.input diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/params.input b/test/multidomain/boundary/stokesdarcy/1p_1p/params_verticalflow.input similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/params.input rename to test/multidomain/boundary/stokesdarcy/1p_1p/params_verticalflow.input diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/problem_darcy.hh similarity index 90% rename from test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_darcy.hh rename to test/multidomain/boundary/stokesdarcy/1p_1p/problem_darcy.hh index 81d52233f4..feec52ceda 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_darcy.hh +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/problem_darcy.hh @@ -31,7 +31,7 @@ #include #include -#include "./../spatialparams.hh" +#include "spatialparams.hh" #include #include @@ -101,6 +101,9 @@ public: : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) { problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; } /*! @@ -111,27 +114,11 @@ public: return problemName_; } - - /*! - * \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]. * @@ -160,6 +147,12 @@ public: if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann(); + if (verticalFlow_) + { + if (onLowerBoundary_(scvf.center())) + values.setAllDirichlet(); + } + return values; } @@ -173,10 +166,7 @@ public: */ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const { - PrimaryVariables values(0.0); - values = initial(element); - - return values; + return initial(element); } /*! @@ -247,21 +237,14 @@ public: { return *couplingManager_; } 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::shared_ptr couplingManager_; std::string problemName_; + bool verticalFlow_; }; } //end namespace diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/problem_stokes.hh similarity index 84% rename from test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_stokes.hh rename to test/multidomain/boundary/stokesdarcy/1p_1p/problem_stokes.hh index 3184e9ba98..f34fa2cd9f 100644 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/horizontalflow/problem_stokes.hh +++ b/test/multidomain/boundary/stokesdarcy/1p_1p/problem_stokes.hh @@ -103,8 +103,10 @@ public: StokesSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager) { - deltaP_ = getParamFromGroup(this->paramGroup(), "Problem.PressureDifference"); problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; } /*! @@ -120,9 +122,6 @@ public: */ // \{ - bool shouldWriteRestartFile() const - { return false; } - /*! * \brief Return the temperature within the domain in [K]. * @@ -159,12 +158,33 @@ public: const auto& globalPos = scvf.dofPosition(); - if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) - values.setDirichlet(Indices::pressureIdx); - else + 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 { - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); + // 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)) @@ -185,10 +205,7 @@ public: */ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const { - PrimaryVariables values(0.0); - values = initialAtPos(globalPos); - - return values; + return initialAtPos(globalPos); } /*! @@ -233,12 +250,21 @@ public: * * \param globalPos The global position */ - PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const { PrimaryVariables values(0.0); - if(onLeftBoundary_(globalPos)) - values[Indices::pressureIdx] = deltaP_; + if (verticalFlow_) + { + static const Scalar inletVelocity = getParamFromGroup(this->paramGroup(), "Problem.Velocity"); + values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]); + } + else // horizontal flow + { + static const Scalar deltaP = getParamFromGroup(this->paramGroup(), "Problem.PressureDifference"); + if(onLeftBoundary_(globalPos)) + values[Indices::pressureIdx] = deltaP; + } return values; } @@ -275,8 +301,8 @@ private: { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } Scalar eps_; - Scalar deltaP_; std::string problemName_; + bool verticalFlow_; std::shared_ptr couplingManager_; }; diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/CMakeLists.txt deleted file mode 100644 index 3ad02d9bb5..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/CMakeLists.txt +++ /dev/null @@ -1,16 +0,0 @@ -add_input_file_links() - -dune_add_test(NAME test_md_boundary_darcy1p_stokes1p_vertical - LABELS multidomain freeflow - SOURCES main.cc - 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_darcy1p_stokes1p_vertical_stokes-00001.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p_stokes1p_vertical_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_vertical_darcy-00001.vtu - - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p_stokes1p_vertical params.input - -Vtk.OutputName test_md_boundary_darcy1p_stokes1p_vertical" - --zeroThreshold {"velocity_liq \(m/s\)_0":6e-17}) diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/main.cc deleted file mode 100644 index f1bc479ced..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/main.cc +++ /dev/null @@ -1,233 +0,0 @@ -// -*- 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 Stokes/Darcy problem (1p) - */ -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -#include "problem_darcy.hh" -#include "problem_stokes.hh" - -namespace Dumux { -namespace Properties { - -template -struct CouplingManager -{ - using Traits = StaggeredMultiDomainTraits; - using type = Dumux::StokesDarcyCouplingManager; -}; - -template -struct CouplingManager -{ - using Traits = StaggeredMultiDomainTraits; - using type = Dumux::StokesDarcyCouplingManager; -}; - -} // end namespace Properties -} // end namespace Dumux - -int main(int argc, char** argv) try -{ - using namespace Dumux; - - // initialize MPI, finalize is done automatically on exit - const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); - - // print dumux start message - if (mpiHelper.rank() == 0) - DumuxMessage::print(/*firstCall=*/true); - - // parse command line arguments and input file - Parameters::init(argc, argv); - - // Define the sub problem type tags - using StokesTypeTag = Properties::TTag::StokesOneP; - using DarcyTypeTag = Properties::TTag::DarcyOneP; - - // try to create a grid (from the given grid file or the input file) - // for both sub-domains - using DarcyGridManager = Dumux::GridManager>; - DarcyGridManager darcyGridManager; - darcyGridManager.init("Darcy"); // pass parameter group - - using StokesGridManager = Dumux::GridManager>; - StokesGridManager stokesGridManager; - stokesGridManager.init("Stokes"); // pass parameter group - - // we compute on the leaf grid view - const auto& darcyGridView = darcyGridManager.grid().leafGridView(); - const auto& stokesGridView = stokesGridManager.grid().leafGridView(); - - // create the finite volume grid geometry - using StokesFVGridGeometry = GetPropType; - auto stokesFvGridGeometry = std::make_shared(stokesGridView); - stokesFvGridGeometry->update(); - using DarcyFVGridGeometry = GetPropType; - auto darcyFvGridGeometry = std::make_shared(darcyGridView); - darcyFvGridGeometry->update(); - - using Traits = StaggeredMultiDomainTraits; - - // the coupling manager - using CouplingManager = StokesDarcyCouplingManager; - auto couplingManager = std::make_shared(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; - auto stokesProblem = std::make_shared(stokesFvGridGeometry, couplingManager); - using DarcyProblem = GetPropType; - auto darcyProblem = std::make_shared(darcyFvGridGeometry, couplingManager); - - // the solution vector - Traits::SolutionVector sol; - sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs()); - sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); - sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; - - // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; - stokesProblem->applyInitialSolution(stokesSol); - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - - darcyProblem->applyInitialSolution(sol[darcyIdx]); - - - couplingManager->init(stokesProblem, darcyProblem, sol); - - // the grid variables - using StokesGridVariables = GetPropType; - auto stokesGridVariables = std::make_shared(stokesProblem, stokesFvGridGeometry); - stokesGridVariables->init(stokesSol); - using DarcyGridVariables = GetPropType; - auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry); - darcyGridVariables->init(sol[darcyIdx]); - - // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); - GetPropType::initOutputModule(stokesVtkWriter); - stokesVtkWriter.write(0.0); - - VtkOutputModule> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name()); - using DarcyVelocityOutput = GetPropType; - darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); - GetPropType::initOutputModule(darcyVtkWriter); - darcyVtkWriter.write(0.0); - - // the assembler for a stationary problem - using Assembler = MultiDomainFVAssembler; - auto assembler = std::make_shared(std::make_tuple(stokesProblem, stokesProblem, darcyProblem), - std::make_tuple(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(); - - // the non-linear solver - using NewtonSolver = MultiDomainNewtonSolver; - NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); - - // solve the non-linear system - nonLinearSolver.solve(sol); - - // write vtk output - stokesVtkWriter.write(1.0); - darcyVtkWriter.write(1.0); - - //////////////////////////////////////////////////////////// - // finalize, print dumux message to say goodbye - //////////////////////////////////////////////////////////// - - // print dumux end message - if (mpiHelper.rank() == 0) - { - Parameters::print(); - DumuxMessage::print(/*firstCall=*/false); - } - - return 0; -} // end main -catch (Dumux::ParameterException &e) -{ - std::cerr << std::endl << e << " ---> Abort!" << std::endl; - return 1; -} -catch (Dune::DGFException & e) -{ - std::cerr << "DGF exception thrown (" << e << - "). Most likely, the DGF file name is wrong " - "or the DGF file is corrupted, " - "e.g. missing hash at end of file or wrong number (dimensions) of entries." - << " ---> Abort!" << std::endl; - return 2; -} -catch (Dune::Exception &e) -{ - std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; - return 3; -} -catch (...) -{ - std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; - return 4; -} diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_darcy.hh deleted file mode 100644 index 8a1d4b1cff..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_darcy.hh +++ /dev/null @@ -1,286 +0,0 @@ -// -*- 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 simple Darcy test problem (cell-centered finite volume method). - */ -#ifndef DUMUX_DARCY_SUBPROBLEM_HH -#define DUMUX_DARCY_SUBPROBLEM_HH - -#include - -#include - -#include -#include - -#include "./../spatialparams.hh" - -#include -#include - -namespace Dumux -{ -template -class DarcySubProblem; - -namespace Properties -{ -// Create new type tags -namespace TTag { -struct DarcyOneP { using InheritsFrom = std::tuple; }; -} // end namespace TTag - -// Set the problem property -template -struct Problem { using type = Dumux::DarcySubProblem; }; - -// the fluid system -template -struct FluidSystem -{ - using Scalar = GetPropType; - using type = FluidSystems::OnePLiquid > ; -}; - -// Set the grid type -template -struct Grid { using type = Dune::YaspGrid<2>; }; - -template -struct SpatialParams -{ - using FVGridGeometry = GetPropType; - using Scalar = GetPropType; - using type = OnePSpatialParams; -}; -} - -template -class DarcySubProblem : public PorousMediumFlowProblem -{ - using ParentType = PorousMediumFlowProblem; - using GridView = GetPropType; - using Scalar = GetPropType; - using PrimaryVariables = GetPropType; - using NumEqVector = GetPropType; - using BoundaryTypes = GetPropType; - using VolumeVariables = GetPropType; - using FVElementGeometry = typename GetPropType::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; - using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using FVGridGeometry = GetPropType; - using ElementVolumeVariables = typename GetPropType::LocalView; - - // copy some indices for convenience - using Indices = typename GetPropType::Indices; - enum { - // grid and world dimension - dim = GridView::dimension, - dimworld = GridView::dimensionworld, - - // primary variable indices - conti0EqIdx = Indices::conti0EqIdx, - pressureIdx = Indices::pressureIdx, - }; - - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - using CouplingManager = GetPropType; - -public: - DarcySubProblem(std::shared_ptr fvGridGeometry, - std::shared_ptr couplingManager) - : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) - { - pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure", 0.0); - problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); - } - - /*! - * \brief The problem name. - */ - const std::string& name() const - { - return problemName_; - } - - /*! - * \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 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(); // right wall, top - - if (onLowerBoundary_(scvf.center())) - values.setAllDirichlet(); - - 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. - */ - 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[pressureIdx] = couplingManager().couplingData().massCouplingCondition(element, 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 - */ - 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 values(0.0); - values[pressureIdx] = pressure_; - - return values; - } - - // \} - - //! Get the coupling manager - const CouplingManager& couplingManager() const - { return *couplingManager_; } - -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_; - Scalar pressure_; - std::string problemName_; - std::shared_ptr couplingManager_; -}; -} //end namespace - -#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_stokes.hh deleted file mode 100644 index 835c33d927..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p_1p/verticalflow/problem_stokes.hh +++ /dev/null @@ -1,307 +0,0 @@ -// -*- 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 NavierStokesTests - * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model. - */ -#ifndef DUMUX_STOKES_SUBPROBLEM_HH -#define DUMUX_STOKES_SUBPROBLEM_HH - -#include - -#include -#include - -#include -#include -#include - -namespace Dumux -{ -template -class StokesSubProblem; - -namespace Properties -{ -// Create new type tags -namespace TTag { -struct StokesOneP { using InheritsFrom = std::tuple; }; -} // end namespace TTag - -// the fluid system -template -struct FluidSystem -{ - using Scalar = GetPropType; - using type = FluidSystems::OnePLiquid > ; -}; - -// Set the grid type -template -struct Grid { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates, 2> >; }; - -// Set the problem property -template -struct Problem { using type = Dumux::StokesSubProblem ; }; - -template -struct EnableFVGridGeometryCache { static constexpr bool value = true; }; -template -struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; -template -struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; - -} - -/*! - * \ingroup NavierStokesTests - * \brief Test problem for the one-phase (Navier-) Stokes problem. - * - * Vertical flow from top to bottom with a parabolic velocity profile. - */ -template -class StokesSubProblem : public NavierStokesProblem -{ - using ParentType = NavierStokesProblem; - - using GridView = GetPropType; - using Scalar = GetPropType; - - using Indices = typename GetPropType::Indices; - - using BoundaryTypes = GetPropType; - - using FVGridGeometry = GetPropType; - using FVElementGeometry = typename FVGridGeometry::LocalView; - using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using Element = typename GridView::template Codim<0>::Entity; - using ElementVolumeVariables = typename GetPropType::LocalView; - using ElementFaceVariables = typename GetPropType::LocalView; - - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - using PrimaryVariables = GetPropType; - using NumEqVector = GetPropType; - - using CouplingManager = GetPropType; - -public: - StokesSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) - : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager) - { - inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.Velocity"); - problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); - } - - /*! - * \brief The problem name. - */ - const std::string& name() const - { - return problemName_; - } - - /*! - * \name Problem parameters - */ - // \{ - - - bool shouldWriteRestartFile() const - { return false; } - - /*! - * \brief Return 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 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 values; - - const auto& globalPos = scvf.dofPosition(); - - // 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); - } - - // outflow/coupling - if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) - { - values.setCouplingNeumann(Indices::conti0EqIdx); - values.setCouplingNeumann(Indices::momentumYBalanceIdx); - values.setBJS(Indices::momentumXBalanceIdx); - } - - return values; - } - - /*! - * \brief Evaluate the boundary conditions for a Dirichlet control volume. - * - * \param element The element - * \param scvf The sub control volume 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 - * - * For this method, the \a values variable stores primary variables. - */ - 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)) - { - values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); - values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); - } - return values; - } - - // \} - - //! Get the coupling manager - const CouplingManager& couplingManager() const - { return *couplingManager_; } - - /*! - * \brief Check if on coupling interface - * - * \param globalPos The global position - * - * Returns true if globalPos is on coupling interface - * (here: lower boundary of Stokes domain) - */ - bool onCouplingInterface(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::velocityYIdx] = inletVelocity_ * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]); - - 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_; - Scalar inletVelocity_; - std::string problemName_; - std::shared_ptr couplingManager_; -}; -} //end namespace - -#endif // DUMUX_STOKES_SUBPROBLEM_HH -- GitLab From b391efedb261e9ec8fc62d277608f5345d51e816 Mon Sep 17 00:00:00 2001 From: Kilian Weishaupt Date: Wed, 28 Nov 2018 18:13:59 +0100 Subject: [PATCH 4/4] [test][md][stokesdarcy] Clean-up 1p2c_1p2c * use only one executable --- .../stokesdarcy/1p2c_1p2c/CMakeLists.txt | 47 ++- .../1p2c_1p2c/horizontalflow/CMakeLists.txt | 14 - .../1p2c_1p2c/{horizontalflow => }/main.cc | 37 +- ...rams.input => params_horizontalflow.input} | 3 - ...params.input => params_verticalflow.input} | 0 ...ut => params_verticalflow_diffusion.input} | 0 .../{horizontalflow => }/problem_darcy.hh | 107 +++--- .../{horizontalflow => }/problem_stokes.hh | 145 ++++--- .../1p2c_1p2c/verticalflow/CMakeLists.txt | 31 -- .../1p2c_1p2c/verticalflow/main.cc | 293 --------------- .../1p2c_1p2c/verticalflow/problem_darcy.hh | 315 ---------------- .../1p2c_1p2c/verticalflow/problem_stokes.hh | 355 ------------------ 12 files changed, 220 insertions(+), 1127 deletions(-) delete mode 100644 test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/CMakeLists.txt rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{horizontalflow => }/main.cc (89%) rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{horizontalflow/params.input => params_horizontalflow.input} (90%) rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{verticalflow/params.input => params_verticalflow.input} (100%) rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{verticalflow/params_diffusion.input => params_verticalflow_diffusion.input} (100%) rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{horizontalflow => }/problem_darcy.hh (82%) rename test/multidomain/boundary/stokesdarcy/1p2c_1p2c/{horizontalflow => }/problem_stokes.hh (70%) delete mode 100644 test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/CMakeLists.txt delete mode 100644 test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/main.cc delete mode 100644 test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_darcy.hh delete mode 100644 test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_stokes.hh diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/CMakeLists.txt index 05a41dcbc9..a393a4b349 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/CMakeLists.txt +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/CMakeLists.txt @@ -1,3 +1,46 @@ -add_subdirectory("verticalflow") -add_subdirectory("horizontalflow") add_subdirectory("diffusionlawcomparison") + +add_input_file_links() + +add_executable(test_md_boundary_darcy1p2c_stokes1p2c EXCLUDE_FROM_ALL main.cc) + +dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_horizontal + LABELS multidomain freeflow 1pnc + TARGET test_md_boundary_darcy1p2c_stokes1p2c + 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_darcy1p2c_stokes1p2c_horizontal_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_stokes-00020.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_darcy-00020.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c params_horizontalflow.input + -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_horizontal") + +dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion + LABELS multidomain freeflow 1pnc + TARGET test_md_boundary_darcy1p2c_stokes1p2c + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --zeroThreshold {"velocity_liq \(m/s\)":1e-20} + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_stokes-00003.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_darcy-00003.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c params_verticalflow_diffusion.input + -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion") + +dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection + LABELS multidomain freeflow 1pnc + TARGET test_md_boundary_darcy1p2c_stokes1p2c + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --zeroThreshold {"velocity_liq \(m/s\)":1e-15} + --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection_stokes-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_stokes-00030.vtu + ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection_darcy-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_darcy-00030.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c params_verticalflow.input + -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_vertical") diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/CMakeLists.txt deleted file mode 100644 index b5aba6a76e..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/CMakeLists.txt +++ /dev/null @@ -1,14 +0,0 @@ -add_input_file_links() - -dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_horizontal - LABELS multidomain freeflow 1pnc - SOURCES main.cc - 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_darcy1p2c_stokes1p2c_horizontal_stokes-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_stokes-00020.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_horizontal_darcy-00020.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_horizontal params.input - -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_horizontal") diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc similarity index 89% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/main.cc rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc index 156c1b88c1..af1c67dc73 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/main.cc +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/main.cc @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -32,7 +33,9 @@ #include #include +#include #include +#include #include #include #include @@ -137,17 +140,11 @@ int main(int argc, char** argv) try auto dt = getParam("TimeLoop.DtInitial"); // instantiate time loop - auto timeLoop = std::make_shared>(0, dt, tEnd); + auto timeLoop = std::make_shared>(0, dt, tEnd); timeLoop->setMaxTimeStepSize(maxDt); - // control the injection period - const Scalar injectionBegin = getParam("Stokes.Problem.InjectionBegin"); - const Scalar injectionEnd = getParam("Stokes.Problem.InjectionEnd"); - - if(injectionBegin > 0.0) - timeLoop->setCheckPoint({injectionBegin, injectionEnd}); - else - timeLoop->setCheckPoint({injectionEnd}); + stokesProblem->setTimeLoop(timeLoop); + darcyProblem->setTimeLoop(timeLoop); // the solution vector Traits::SolutionVector sol; @@ -155,20 +152,12 @@ int main(int argc, char** argv) try sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; + // get a solution vector storing references to the two Stokes solution vectors + auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx); // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; stokesProblem->applyInitialSolution(stokesSol); - auto solStokesOld = stokesSol; - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; auto solOld = sol; @@ -183,7 +172,7 @@ int main(int argc, char** argv) try darcyGridVariables->init(sol[darcyIdx]); // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); + StaggeredVtkOutputModule stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); GetPropType::initOutputModule(stokesVtkWriter); stokesVtkWriter.write(0.0); @@ -213,20 +202,12 @@ int main(int argc, char** argv) try using NewtonSolver = MultiDomainNewtonSolver; NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); - constexpr auto eps = 1e-6; - // time loop timeLoop->start(); do { // set previous solution for storage evaluations assembler->setPreviousSolution(solOld); - - if(timeLoop->time() > injectionBegin - eps && timeLoop->time() < injectionEnd + eps) - stokesProblem->setInjectionState(true); - else - stokesProblem->setInjectionState(false); - // solve the non-linear system with time step control nonLinearSolver.solve(sol, *timeLoop); diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/params.input b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_horizontalflow.input similarity index 90% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/params.input rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_horizontalflow.input index 36c9624800..302eb4565a 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/params.input +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_horizontalflow.input @@ -14,15 +14,12 @@ Cells = 20 20 [Stokes.Problem] Name = stokes Velocity = 1e-6 -Pressure = 1.0e5 InletMoleFraction = 1e-3 InjectionBegin = 0 InjectionEnd = 3e6 [Darcy.Problem] Name = darcy -Pressure = 1.0e5 -InitialMoleFraction = 0.0 [SpatialParams] AlphaBeaversJoseph = 1.0 diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/params.input b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_verticalflow.input similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/params.input rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_verticalflow.input diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/params_diffusion.input b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_verticalflow_diffusion.input similarity index 100% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/params_diffusion.input rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/params_verticalflow_diffusion.input diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_darcy.hh similarity index 82% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_darcy.hh rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_darcy.hh index 6a39e93357..0c895eb348 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_darcy.hh +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_darcy.hh @@ -31,7 +31,7 @@ #include #include -#include "./../spatialparams.hh" +#include "spatialparams.hh" #include #include @@ -93,42 +93,32 @@ template class DarcySubProblem : public PorousMediumFlowProblem { using ParentType = PorousMediumFlowProblem; - using GridView = GetPropType; + using FVGridGeometry = GetPropType; + using GridView = typename FVGridGeometry::GridView; using Scalar = GetPropType; using PrimaryVariables = GetPropType; - using FluidSystem = GetPropType; using NumEqVector = GetPropType; using BoundaryTypes = GetPropType; using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using FVGridGeometry = GetPropType; - - // copy some indices for convenience + using Element = typename GridView::template Codim<0>::Entity; using Indices = typename GetPropType::Indices; - enum { - // grid and world dimension - dim = GridView::dimension, - dimworld = GridView::dimensionworld, - - // primary variable indices - conti0EqIdx = Indices::conti0EqIdx, - pressureIdx = Indices::pressureIdx, - }; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = Dune::FieldVector; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using CouplingManager = GetPropType; + using TimeLoopPtr = std::shared_ptr>; public: DarcySubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) { - pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); - initialMoleFraction_ = getParamFromGroup(this->paramGroup(), "Problem.InitialMoleFraction"); problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; } /*! @@ -139,26 +129,6 @@ public: return problemName_; } - /*! - * \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]. * @@ -187,6 +157,46 @@ public: if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann(); + if (verticalFlow_) + { + if (onLowerBoundary_(scvf.center())) + values.setAllDirichlet(); + } + + 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 + * + */ + PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const + { + PrimaryVariables values(0.0); + values = initial(element); + + if (verticalFlow_) + { + // Check if this a pure diffusion problem. + static const bool isDiffusionProblem = problemName_.find("diffusion") != std::string::npos; + + Scalar bottomMoleFraction = 0.0; + + if (isDiffusionProblem) + { + // For the diffusion problem, change the top mole fraction after some time + // in order to revert the concentration gradient. + if (time() >= 1e10) + bottomMoleFraction = 1e-3; + } + + if(onLowerBoundary_(scvf.center())) + values[Indices::conti0EqIdx + 1] = bottomMoleFraction; + } + return values; } @@ -249,8 +259,7 @@ public: PrimaryVariables initial(const Element &element) const { PrimaryVariables values(0.0); - values[pressureIdx] = pressure_; - values[conti0EqIdx + 1] = initialMoleFraction_; + values[Indices::pressureIdx] = 1e5; return values; } @@ -261,6 +270,18 @@ public: const CouplingManager& couplingManager() const { return *couplingManager_; } + /*! + * \brief Sets the time loop pointer + */ + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } + + /*! + * \brief Returns the time + */ + Scalar time() const + { return timeLoop_->time(); } + private: bool onLeftBoundary_(const GlobalPosition &globalPos) const { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; } @@ -275,10 +296,10 @@ private: { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } Scalar eps_; - Scalar pressure_; - Scalar initialMoleFraction_; std::string problemName_; + bool verticalFlow_; std::shared_ptr couplingManager_; + TimeLoopPtr timeLoop_; }; } //end namespace diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_stokes.hh similarity index 70% rename from test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_stokes.hh rename to test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_stokes.hh index 41c6ef3fc4..2fbef2e814 100644 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/problem_stokes.hh +++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/problem_stokes.hh @@ -91,7 +91,6 @@ class StokesSubProblem : public NavierStokesProblem using GridView = GetPropType; using Scalar = GetPropType; using Indices = typename GetPropType::Indices; - using FluidSystem = GetPropType; using BoundaryTypes = GetPropType; using FVGridGeometry = GetPropType; using FVElementGeometry = typename FVGridGeometry::LocalView; @@ -103,15 +102,16 @@ class StokesSubProblem : public NavierStokesProblem using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using CouplingManager = GetPropType; + using TimeLoopPtr = std::shared_ptr>; public: StokesSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) - : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), injectionState_(false), couplingManager_(couplingManager) + : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager) { - inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.Velocity"); - pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure"); - inletMoleFraction_ = getParamFromGroup(this->paramGroup(), "Problem.InletMoleFraction"); problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); + + // determine whether to simulate a vertical or horizontal flow configuration + verticalFlow_ = problemName_.find("vertical") != std::string::npos; } /*! @@ -127,9 +127,6 @@ public: */ // \{ - bool shouldWriteRestartFile() const - { return false; } - /*! * \brief Return the temperature within the domain in [K]. * @@ -166,23 +163,46 @@ public: const auto& globalPos = scvf.dofPosition(); - if(onLeftBoundary_(globalPos)) - { - values.setDirichlet(Indices::conti0EqIdx + 1); - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); - } - else if(onRightBoundary_(globalPos)) + if (verticalFlow_) { - values.setDirichlet(Indices::pressureIdx); - values.setOutflow(Indices::conti0EqIdx + 1); + // inflow + if(onUpperBoundary_(globalPos)) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + values.setDirichlet(Indices::conti0EqIdx + 1); + } + + // left/right wall + if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos))) + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + values.setNeumann(Indices::conti0EqIdx); + values.setNeumann(Indices::conti0EqIdx + 1); + } + } - else + else // horizontal flow { - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); - values.setNeumann(Indices::conti0EqIdx); - values.setNeumann(Indices::conti0EqIdx + 1); + if (onLeftBoundary_(globalPos)) + { + values.setDirichlet(Indices::conti0EqIdx + 1); + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + } + else if (onRightBoundary_(globalPos)) + { + values.setDirichlet(Indices::pressureIdx); + values.setOutflow(Indices::conti0EqIdx + 1); + } + else + { + values.setDirichlet(Indices::velocityXIdx); + values.setDirichlet(Indices::velocityYIdx); + values.setNeumann(Indices::conti0EqIdx); + values.setNeumann(Indices::conti0EqIdx + 1); + } } if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) @@ -207,9 +227,37 @@ public: PrimaryVariables values(0.0); values = initialAtPos(globalPos); - // start injecting after the velocity field had enough time to initialize - if(globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_ && isInjectionPeriod()) - values[Indices::conti0EqIdx + 1] = inletMoleFraction_; + if (verticalFlow_) + { + // Check if this a pure diffusion problem. + static const bool isDiffusionProblem = problemName_.find("diffusion") != std::string::npos; + + Scalar topMoleFraction = 1e-3; + + if (isDiffusionProblem) + { + // For the diffusion problem, change the top mole fraction after some time + // in order to revert the concentration gradient. + if (time() >= 1e10) + topMoleFraction = 0.0; + } + else // advection problem + { + // reverse the flow direction after some time for the advection problem + if (time() >= 3e5) + values[Indices::velocityYIdx] *= -1.0; + } + + if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_) + values[Indices::conti0EqIdx + 1] = topMoleFraction; + } + else // horizontal flow + { + static const Scalar inletMoleFraction = getParamFromGroup(this->paramGroup(), "Problem.InletMoleFraction"); + if(globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_) + values[Indices::conti0EqIdx + 1] = inletMoleFraction; + } + return values; } @@ -259,14 +307,25 @@ public: * * \param globalPos The global position */ - PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const { PrimaryVariables values(0.0); - values[Indices::pressureIdx] = pressure_; - values[Indices::velocityXIdx] = inletVelocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1]) - * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1]) - / (0.25 * (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]) - * (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])); + values[Indices::pressureIdx] = 1e5; + + static const Scalar vMax = getParamFromGroup(this->paramGroup(), "Problem.Velocity", 0.0); + + auto parabolicProfile = [&](const GlobalPosition& globalPos, int coord) + { + return vMax * (globalPos[coord] - this->fvGridGeometry().bBoxMin()[coord]) + * (this->fvGridGeometry().bBoxMax()[coord] - globalPos[coord]) + / (0.25 * (this->fvGridGeometry().bBoxMax()[coord] - this->fvGridGeometry().bBoxMin()[coord]) + * (this->fvGridGeometry().bBoxMax()[coord] - this->fvGridGeometry().bBoxMin()[coord])); + }; + + if (verticalFlow_) + values[Indices::velocityYIdx] = parabolicProfile(globalPos, 0); + else // horizontal flow + values[Indices::velocityXIdx] = parabolicProfile(globalPos, 1); return values; } @@ -287,15 +346,17 @@ public: return 1.0; } - void setInjectionState(const bool yesNo) - { - injectionState_ = yesNo; - } + /*! + * \brief Sets the time loop pointer + */ + void setTimeLoop(TimeLoopPtr timeLoop) + { timeLoop_ = timeLoop; } - bool isInjectionPeriod() const - { - return injectionState_; - } + /*! + * \brief Returns the time + */ + Scalar time() const + { return timeLoop_->time(); } // \} @@ -313,12 +374,10 @@ private: { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } Scalar eps_; - Scalar inletVelocity_; - Scalar pressure_; - Scalar inletMoleFraction_; - bool injectionState_; + bool verticalFlow_; std::string problemName_; std::shared_ptr couplingManager_; + TimeLoopPtr timeLoop_; }; } //end namespace diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/CMakeLists.txt deleted file mode 100644 index 1574bd3746..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/CMakeLists.txt +++ /dev/null @@ -1,31 +0,0 @@ -add_input_file_links() - -add_executable(test_md_boundary_darcy1p2c_stokes1p2c_vertical EXCLUDE_FROM_ALL main.cc) - -dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion - LABELS multidomain freeflow 1pnc - TARGET test_md_boundary_darcy1p2c_stokes1p2c_vertical - CMAKE_GUARD HAVE_UMFPACK - COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMD_ARGS --script fuzzy - --zeroThreshold {"velocity_liq \(m/s\)":1e-20} - --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_stokes-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_stokes-00003.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion_darcy-00003.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical params_diffusion.input - -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_vertical_diffusion") - -dune_add_test(NAME test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection - LABELS multidomain freeflow 1pnc - TARGET test_md_boundary_darcy1p2c_stokes1p2c_vertical - CMAKE_GUARD HAVE_UMFPACK - COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - CMD_ARGS --script fuzzy - --zeroThreshold {"velocity_liq \(m/s\)":1e-15} - --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection_stokes-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_stokes-00030.vtu - ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy1p2c_stokes1p2c_vertical_advection_darcy-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical_darcy-00030.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy1p2c_stokes1p2c_vertical params.input - -Vtk.OutputName test_md_boundary_darcy1p2c_stokes1p2c_vertical") diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/main.cc b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/main.cc deleted file mode 100644 index 3b913a3ebb..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/main.cc +++ /dev/null @@ -1,293 +0,0 @@ -// -*- 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 Stokes/Darcy problem (1p) - */ -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -#include "problem_darcy.hh" -#include "problem_stokes.hh" - -namespace Dumux { -namespace Properties { - -template -struct CouplingManager -{ - using Traits = StaggeredMultiDomainTraits; - using type = Dumux::StokesDarcyCouplingManager; -}; - -template -struct CouplingManager -{ - using Traits = StaggeredMultiDomainTraits; - using type = Dumux::StokesDarcyCouplingManager; -}; - -} // end namespace Properties -} // end namespace Dumux - -int main(int argc, char** argv) try -{ - using namespace Dumux; - - // initialize MPI, finalize is done automatically on exit - const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); - - // print dumux start message - if (mpiHelper.rank() == 0) - DumuxMessage::print(/*firstCall=*/true); - - // parse command line arguments and input file - Parameters::init(argc, argv); - - // Define the sub problem type tags - using StokesTypeTag = Properties::TTag::StokesOnePTwoC; - using DarcyTypeTag = Properties::TTag::DarcyOnePTwoC; - - // try to create a grid (from the given grid file or the input file) - // for both sub-domains - using DarcyGridManager = Dumux::GridManager>; - DarcyGridManager darcyGridManager; - darcyGridManager.init("Darcy"); // pass parameter group - - using StokesGridManager = Dumux::GridManager>; - StokesGridManager stokesGridManager; - stokesGridManager.init("Stokes"); // pass parameter group - - // we compute on the leaf grid view - const auto& darcyGridView = darcyGridManager.grid().leafGridView(); - const auto& stokesGridView = stokesGridManager.grid().leafGridView(); - - // create the finite volume grid geometry - using StokesFVGridGeometry = GetPropType; - auto stokesFvGridGeometry = std::make_shared(stokesGridView); - stokesFvGridGeometry->update(); - using DarcyFVGridGeometry = GetPropType; - auto darcyFvGridGeometry = std::make_shared(darcyGridView); - darcyFvGridGeometry->update(); - - using Traits = StaggeredMultiDomainTraits; - - // the coupling manager - using CouplingManager = StokesDarcyCouplingManager; - auto couplingManager = std::make_shared(stokesFvGridGeometry, darcyFvGridGeometry); - - // the indices - constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx; - constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx; - constexpr auto darcyIdx = CouplingManager::darcyIdx; - - // get some time loop parameters - using Scalar = GetPropType; - const auto tEnd = getParam("TimeLoop.TEnd"); - const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); - auto dt = getParam("TimeLoop.DtInitial"); - - const bool isDiffusionProblem = getParam("Problem.OnlyDiffusion", false); - - // the problem (initial and boundary conditions) - using StokesProblem = GetPropType; - auto stokesProblem = std::make_shared(stokesFvGridGeometry, couplingManager); - using DarcyProblem = GetPropType; - auto darcyProblem = std::make_shared(darcyFvGridGeometry, couplingManager); - - // initialize the fluidsystem (tabulation) - GetPropType::init(); - - // instantiate time loop - auto timeLoop = std::make_shared>(0, dt, tEnd); - timeLoop->setMaxTimeStepSize(maxDt); - - if(!isDiffusionProblem) - timeLoop->setPeriodicCheckPoint(tEnd/30.0); - // else - // timeLoop->setCheckPoint({1e10, 1.1e10, 1.2e10}); - // TODO This does not work due to an issue with Dune's eq() - - stokesProblem->setTimeLoop(timeLoop); - - // the solution vector - Traits::SolutionVector sol; - sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs()); - sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs()); - sol[darcyIdx].resize(darcyFvGridGeometry->numDofs()); - - const auto& cellCenterSol = sol[stokesCellCenterIdx]; - const auto& faceSol = sol[stokesFaceIdx]; - - // apply initial solution for instationary problems - GetPropType stokesSol; - std::get<0>(stokesSol) = cellCenterSol; - std::get<1>(stokesSol) = faceSol; - stokesProblem->applyInitialSolution(stokesSol); - auto solStokesOld = stokesSol; - sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx]; - sol[stokesFaceIdx] = stokesSol[stokesFaceIdx]; - - darcyProblem->applyInitialSolution(sol[darcyIdx]); - auto solDarcyOld = sol[darcyIdx]; - - auto solOld = sol; - - couplingManager->init(stokesProblem, darcyProblem, sol); - - // the grid variables - using StokesGridVariables = GetPropType; - auto stokesGridVariables = std::make_shared(stokesProblem, stokesFvGridGeometry); - stokesGridVariables->init(stokesSol); - using DarcyGridVariables = GetPropType; - auto darcyGridVariables = std::make_shared(darcyProblem, darcyFvGridGeometry); - darcyGridVariables->init(sol[darcyIdx]); - - // intialize the vtk output module - StaggeredVtkOutputModule> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name()); - GetPropType::initOutputModule(stokesVtkWriter); - stokesVtkWriter.write(0.0); - - VtkOutputModule> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name()); - using DarcyVelocityOutput = GetPropType; - darcyVtkWriter.addVelocityOutput(std::make_shared(*darcyGridVariables)); - GetPropType::initOutputModule(darcyVtkWriter); - darcyVtkWriter.write(0.0); - - // the assembler with time loop for instationary problem - using Assembler = MultiDomainFVAssembler; - auto assembler = std::make_shared(std::make_tuple(stokesProblem, stokesProblem, darcyProblem), - std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(), - stokesFvGridGeometry->faceFVGridGeometryPtr(), - darcyFvGridGeometry), - std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(), - stokesGridVariables->faceGridVariablesPtr(), - darcyGridVariables), - couplingManager, - timeLoop); - - // the linear solver - using LinearSolver = UMFPackBackend; - auto linearSolver = std::make_shared(); - - // the non-linear solver - using NewtonSolver = MultiDomainNewtonSolver; - NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager); - - // time loop - timeLoop->start(); do - { - // set previous solution for storage evaluations - assembler->setPreviousSolution(solOld); - - // revert the concentration gradient after some time - if(isDiffusionProblem && timeLoop->time() > 1e10) - { - darcyProblem->setBottomMoleFraction(1e-3); - stokesProblem->setTopMoleFraction(0.0); - } - - // solve the non-linear system with time step control - nonLinearSolver.solve(sol, *timeLoop); - - // make the new solution the old solution - solOld = sol; - stokesGridVariables->advanceTimeStep(); - darcyGridVariables->advanceTimeStep(); - - // advance to the time loop to the next step - timeLoop->advanceTimeStep(); - - // write vtk output -// vtkWriter.write(timeLoop->time()); - stokesVtkWriter.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(stokesGridView.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/stokesdarcy/1p2c_1p2c/verticalflow/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_darcy.hh deleted file mode 100644 index 438b364dd4..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_darcy.hh +++ /dev/null @@ -1,315 +0,0 @@ -// -*- 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 TODO docme - */ -#ifndef DUMUX_DARCY_SUBPROBLEM_HH -#define DUMUX_DARCY_SUBPROBLEM_HH - -#include - -#include - -#include -#include - -#include "./../spatialparams.hh" - -#include -#include -#include - -namespace Dumux -{ -template -class DarcySubProblem; - -namespace Properties -{ -// Create new type tags -namespace TTag { -struct DarcyOnePTwoC { using InheritsFrom = std::tuple; }; -} // end namespace TTag - -// Set the problem property -template -struct Problem { using type = Dumux::DarcySubProblem; }; - -// The fluid system -template -struct FluidSystem -{ - using H2OAir = FluidSystems::H2OAir>; - static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase - using type = FluidSystems::OnePAdapter; -}; - -// Use moles -template -struct UseMoles { static constexpr bool value = true; }; - -// Do not replace one equation with a total mass balance -template -struct ReplaceCompEqIdx { static constexpr int value = 3; }; - -//! Use a model with constant tortuosity for the effective diffusivity -template -struct EffectiveDiffusivityModel -{ using type = DiffusivityConstantTortuosity>; }; - -// Set the grid type -template -struct Grid { using type = Dune::YaspGrid<2>; }; - -// Set the spatial paramaters type -template -struct SpatialParams -{ - using FVGridGeometry = GetPropType; - using Scalar = GetPropType; - using type = OnePSpatialParams; -}; -} - -template -class DarcySubProblem : public PorousMediumFlowProblem -{ - using ParentType = PorousMediumFlowProblem; - using GridView = GetPropType; - using Scalar = GetPropType; - using PrimaryVariables = GetPropType; - using FluidSystem = GetPropType; - using NumEqVector = GetPropType; - using BoundaryTypes = GetPropType; - using VolumeVariables = GetPropType; - using FVElementGeometry = typename GetPropType::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; - using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using FVGridGeometry = GetPropType; - using ElementVolumeVariables = typename GetPropType::LocalView; - - // copy some indices for convenience - using Indices = typename GetPropType::Indices; - enum { - // grid and world dimension - dim = GridView::dimension, - dimworld = GridView::dimensionworld, - - // primary variable indices - conti0EqIdx = Indices::conti0EqIdx, - pressureIdx = Indices::pressureIdx, - }; - - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - using CouplingManager = GetPropType; - -public: - DarcySubProblem(std::shared_ptr fvGridGeometry, - std::shared_ptr couplingManager) - : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager), xBottom_(0.0) - { - pressure_ = getParamFromGroup(this->paramGroup(), "Problem.Pressure", 0.0); - problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); - } - - /*! - * \brief The problem name. - */ - const std::string& name() const - { - return problemName_; - } - - /*! - * \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 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(); // left/right wall, top - - if(onLowerBoundary_(scvf.center())) - values.setAllDirichlet(); - - 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); - - if(onLowerBoundary_(scvf.center())) - values[Indices::conti0EqIdx + 1] = xBottom_; - - 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. - */ - 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 = couplingManager().couplingData().massCouplingCondition(element, 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 - */ - 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 values(0.0); - values[pressureIdx] = pressure_; - - return values; - } - - // \} - - /*! - * \brief Get the coupling manager - */ - const CouplingManager& couplingManager() const - { return *couplingManager_; } - - /*! - * \brief Set the fixed value of the mole fraction at the bottom of the domain - */ - void setBottomMoleFraction(const Scalar x) - { xBottom_ = x; } - -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::shared_ptr couplingManager_; - Scalar xBottom_; - Scalar pressure_; - std::string problemName_; -}; -} //end namespace - -#endif //DUMUX_DARCY_SUBPROBLEM_HH diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_stokes.hh deleted file mode 100644 index 30859be729..0000000000 --- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/problem_stokes.hh +++ /dev/null @@ -1,355 +0,0 @@ -// -*- 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 NavierStokesTests - * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model. - */ -#ifndef DUMUX_STOKES_1P2C_SUBPROBLEM_HH -#define DUMUX_STOKES_1P2C_SUBPROBLEM_HH - -#include - -#include -#include - -#include -#include -#include - -namespace Dumux -{ -template -class StokesSubProblem; - -namespace Properties -{ -// Create new type tags -namespace TTag { -struct StokesOnePTwoC { using InheritsFrom = std::tuple; }; -} // end namespace TTag - -// The fluid system -template -struct FluidSystem -{ - using H2OAir = FluidSystems::H2OAir>; - static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase - using type = FluidSystems::OnePAdapter; -}; - -// Set the grid type -template -struct Grid { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates, 2> >; }; - -// Set the problem property -template -struct Problem { using type = Dumux::StokesSubProblem ; }; - -// Enable all caches -template -struct EnableFVGridGeometryCache { static constexpr bool value = true; }; -template -struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; -template -struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; - -// Use moles -template -struct UseMoles { static constexpr bool value = true; }; - -// Do not replace one equation with a total mass balance -template -struct ReplaceCompEqIdx { static constexpr int value = 3; }; -} - -/*! - * \ingroup NavierStokesTests - * \brief Test problem for the one-phase (Navier-) Stokes problem. - * - * Vertical flow from top to bottom with a parabolic velocity profile. - */ -template -class StokesSubProblem : public NavierStokesProblem -{ - using ParentType = NavierStokesProblem; - - using GridView = GetPropType; - using Scalar = GetPropType; - using Indices = typename GetPropType::Indices; - using BoundaryTypes = GetPropType; - using FluidSystem = GetPropType; - using FVGridGeometry = GetPropType; - using FVElementGeometry = typename FVGridGeometry::LocalView; - using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; - using Element = typename GridView::template Codim<0>::Entity; - using ElementVolumeVariables = typename GetPropType::LocalView; - using ElementFaceVariables = typename GetPropType::LocalView; - - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - - using PrimaryVariables = GetPropType; - using NumEqVector = GetPropType; - - using CouplingManager = GetPropType; - using TimeLoopPtr = std::shared_ptr>; - -public: - StokesSubProblem(std::shared_ptr fvGridGeometry, std::shared_ptr couplingManager) - : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager), xTop_(1e-3) - { - inletVelocity_ = getParamFromGroup(this->paramGroup(), "Problem.Velocity"); - problemName_ = getParam("Vtk.OutputName") + "_" + getParamFromGroup(this->paramGroup(), "Problem.Name"); - } - - /*! - * \brief The problem name. - */ - const std::string& name() const - { - return problemName_; - } - - /*! - * \name Problem parameters - */ - // \{ - - - bool shouldWriteRestartFile() const - { return false; } - - /*! - * \brief Return 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 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 values; - - const auto& globalPos = scvf.dofPosition(); - - // inflow - if(onUpperBoundary_(globalPos)) - { - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); - values.setDirichlet(Indices::conti0EqIdx + 1); - } - - // left/right wall - if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos))) - { - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); - values.setNeumann(Indices::conti0EqIdx); - values.setNeumann(Indices::conti0EqIdx + 1); - } - - // outflow/coupling - if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) - { - values.setCouplingNeumann(Indices::conti0EqIdx); - values.setCouplingNeumann(Indices::conti0EqIdx + 1); - values.setCouplingNeumann(Indices::momentumYBalanceIdx); - values.setBJS(Indices::momentumXBalanceIdx); - } - - return values; - } - - /*! - * \brief Evaluate the boundary conditions for a Dirichlet control volume. - * - * \param element The element - * \param scvf The sub control volume face - */ - PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const - { - PrimaryVariables values(0.0); - values = initialAtPos(pos); - - if(pos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_) - values[Indices::conti0EqIdx + 1] = xTop_; - - // reverse the flow direction after some time - if(time() >= 3e5) - values[Indices::velocityYIdx] *= -1.0; - - 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 - * - * For this method, the \a values variable stores primary variables. - */ - 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)) - { - values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); - - const auto tmp = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); - values[Indices::conti0EqIdx] = tmp[0]; - values[Indices::conti0EqIdx + 1] = tmp[1]; - } - return values; - } - - // \} - - /*! - * \brief Get the coupling manager - */ - const CouplingManager& couplingManager() const - { return *couplingManager_; } - - /*! - * \brief Check if on coupling interface - * - * \param globalPos The global position - * - * Returns true if globalPos is on coupling interface - * (here: lower boundary of Stokes domain) - */ - bool onCouplingInterface(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] = 1e5; - values[Indices::velocityYIdx] = inletVelocity_ * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]); - values[Indices::conti0EqIdx + 1] = 0.0; - - 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()); - } - - /*! - * \brief Sets the time loop pointer - */ - void setTimeLoop(TimeLoopPtr timeLoop) - { timeLoop_ = timeLoop; } - - /*! - * \brief Returns the time - */ - Scalar time() const - { return timeLoop_->time(); } - - /*! - * \brief Set the fixed value of the mole fraction at the top of the domain - */ - void setTopMoleFraction(const Scalar x) - { xTop_ = x; } - - // \} - -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_; - Scalar inletVelocity_; - std::string problemName_; - - std::shared_ptr couplingManager_; - Scalar xTop_; - - TimeLoopPtr timeLoop_; -}; -} //end namespace - -#endif // DUMUX_STOKES_SUBPROBLEM_HH -- GitLab