Commit a1e232e1 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[test][md][stokesdarcy] Clean-up 1p_1p

* use only one executable
parent 07b043d1
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})
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")
......@@ -32,6 +32,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/partial.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
......@@ -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<StokesTypeTag, Properties::SolutionVector> 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<StokesGridVariables, GetPropType<StokesTypeTag, Properties::SolutionVector>> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
stokesVtkWriter.write(0.0);
......
......@@ -31,7 +31,7 @@
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "./../spatialparams.hh"
#include "spatialparams.hh"
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
......@@ -101,6 +101,9 @@ public:
: ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
{
problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
// determine whether to simulate a vertical or horizontal flow configuration
verticalFlow_ = problemName_.find("vertical") != std::string::npos;
}
/*!
......@@ -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> couplingManager_;
std::string problemName_;
bool verticalFlow_;
};
} //end namespace
......
......@@ -103,8 +103,10 @@ public:
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
{
deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
problemName_ = getParam<std::string>("Vtk.OutputName") + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
// determine whether to simulate a vertical or horizontal flow configuration
verticalFlow_ = problemName_.find("vertical") != std::string::npos;
}
/*!
......@@ -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<Scalar>(this->paramGroup(), "Problem.Velocity");
values[Indices::velocityYIdx] = inletVelocity * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
}
else // horizontal flow
{
static const Scalar deltaP = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
if(onLeftBoundary_(globalPos))
values[Indices::pressureIdx] = deltaP;
}
return values;
}
......@@ -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> couplingManager_;
};
......
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})
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief A test problem for the coupled Stokes/Darcy problem (1p)
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "problem_darcy.hh"
#include "problem_stokes.hh"
namespace Dumux {
namespace Properties {
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::StokesOneP>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::DarcyOneP>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
} // end namespace Properties
} // end namespace Dumux
int main(int argc, char** argv) try
{
using namespace Dumux;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv);
// Define the sub problem type tags
using StokesTypeTag = Properties::TTag::StokesOneP;
using DarcyTypeTag = Properties::TTag::DarcyOneP;
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
// we compute on the leaf grid view
const auto& darcyGridView = darcyGridManager.grid().leafGridView();
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
stokesFvGridGeometry->update();
using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
darcyFvGridGeometry->update();
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// the solution vector
Traits::SolutionVector sol;
sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
const auto& cellCenterSol = sol[stokesCellCenterIdx];
const auto& faceSol = sol[stokesFaceIdx];
// apply initial solution for instationary problems
GetPropType<StokesTypeTag, Properties::SolutionVector> 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<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
stokesGridVariables->init(stokesSol);
using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
StaggeredVtkOutputModule<StokesGridVariables, GetPropType<StokesTypeTag, Properties::SolutionVector>> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesProblem->name());
GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
stokesVtkWriter.write(0.0);
VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
darcyVtkWriter.write(0.0);
// the assembler for a stationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
stokesFvGridGeometry->faceFVGridGeometryPtr(),
darcyFvGridGeometry),
std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
stokesGridVariables->faceGridVariablesPtr(),
darcyGridVariables),
couplingManager);
// the linear solver
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// solve the non-linear system
nonLinearSolver.solve(sol);
// write vtk output
stokesVtkWriter.write(1.0);
darcyVtkWriter.write(1.0);
////////////////////////////////////////////////////////////
// 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;
}
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief A simple Darcy test problem (cell-centered finite volume method).
*/
#ifndef DUMUX_DARCY_SUBPROBLEM_HH
#define DUMUX_DARCY_SUBPROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "./../spatialparams.hh"
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
namespace Dumux
{
template <class TypeTag>