Commit 741cce5b authored by Alexander Jaust's avatar Alexander Jaust
Browse files

files for reverse coupling of ff-pm coupling

Change coupling to:
pm receives pressure and communicates velocity

WIP
parent 59a673a4
# executables for the iterative case
string(REPLACE ":" ";" LIBRARY_DIRS $ENV{LD_LIBRARY_PATH})
find_library(PRECICE_LIB NAMES precice HINTS ${LIBRARY_DIRS})
find_package(Boost 1.65.1 REQUIRED COMPONENTS log system) #Require same version as preCICE
#dune_add_test(NAME test_ff_reversed
# SOURCES main_ff-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
#dune_add_test(NAME test_pm_reversed
# SOURCES main_pm-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
add_executable(test_ff_reversed EXCLUDE_FROM_ALL main_ff-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
add_executable(test_pm_reversed EXCLUDE_FROM_ALL main_pm-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
target_compile_definitions(test_ff_reversed PUBLIC "ENABLEMONOLITHIC=0")
target_compile_definitions(test_pm_reversed PUBLIC "ENABLEMONOLITHIC=0")
target_compile_definitions(test_ff_reversed PRIVATE BOOST_ALL_DYN_LINK)
target_link_libraries(test_ff_reversed PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
target_include_directories(test_ff_reversed PRIVATE ${Boost_INCLUDE_DIRS})
target_compile_definitions(test_pm_reversed PRIVATE BOOST_ALL_DYN_LINK)
target_link_libraries(test_pm_reversed PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
target_include_directories(test_pm_reversed PRIVATE ${Boost_INCLUDE_DIRS})
# add a symlink for each input file
add_input_file_links()
// -*- 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 The free flow sub problem
*/
#ifndef DUMUX_STOKES_SUBPROBLEM_HH
#define DUMUX_STOKES_SUBPROBLEM_HH
#ifndef ENABLEMONOLITHIC
#define ENABLEMONOLITHIC 0
#endif
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
#include "../../precice-adapter/include/preciceadapter.hh"
namespace Dumux
{
template <class TypeTag>
class StokesSubProblem;
namespace Properties
{
// Create new type tags
namespace TTag {
struct FreeFlowModel { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
} // end namespace TTag
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::FreeFlowModel>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
};
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::FreeFlowModel>
{
static constexpr auto dim = 2;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
//****** comment out for the last exercise *****//
using type = TensorGrid;
//****** uncomment for the last exercise *****//
// using HostGrid = TensorGrid;
// using type = Dune::SubGrid<dim, HostGrid>;
};
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::FreeFlowModel> { using type = Dumux::StokesSubProblem<TypeTag> ; };
template<class TypeTag>
struct EnableFVGridGeometryCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; };
}
/*!
* \brief The free flow sub problem
*/
template <class TypeTag>
class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
#if ENABLEMONOLITHIC
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
#endif
public:
#if ENABLEMONOLITHIC
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
#else
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry, "FreeFlow"),
eps_(1e-6),
couplingInterface_(precice_adapter::PreciceAdapter::getInstance() ),
pressureId_(0),
velocityId_(0),
dataIdsWereSet_(false)
#endif
{
deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
// pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
// velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \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();
// left/right wall
if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
{
values.setDirichlet(Indices::pressureIdx);
}
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
// coupling interface
#if ENABLEMONOLITHIC
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values.setCouplingNeumann(Indices::conti0EqIdx);
values.setCouplingNeumann(Indices::momentumYBalanceIdx);
values.setBJS(Indices::momentumXBalanceIdx);
}
#else
// // TODO do preCICE stuff in analogy to heat transfer
assert( dataIdsWereSet_ );
const auto faceId = scvf.index();
if ( couplingInterface_.isCoupledEntity(faceId) )
{
//TODO What do I want to do here?
values.setNeumann(Indices::conti0EqIdx);
values.setNeumann(Indices::momentumYBalanceIdx);
//values.setCouplingNeumann(Indices::conti0EqIdx);
//values.setCouplingNeumann(Indices::momentumYBalanceIdx);
//values.setBJS(Indices::momentumXBalanceIdx);
}
#endif
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Dirichlet control volume.
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
values = initialAtPos(globalPos);
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFaceVars The element face variables
* \param scvf The boundary sub control volume face
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
#if ENABLEMONOLITHIC
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);
}
#else
assert( dataIdsWereSet_ );
const auto faceId = scvf.index();
if( couplingInterface_.isCoupledEntity( faceId ) )
{
//const Scalar density = 1000; // TODO how to handle compressible fluids?
const auto velocity = couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
// const Scalar ccPressure = volVars.pressure();
// const Scalar mobility = volVars.mobility();
const Scalar density = volVars.density();
// const auto K = volVars.permeability();
// // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
// const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, this->gravity() );
// auto distanceVector = scvf.center() - element.geometry().center();
// distanceVector /= distanceVector.two_norm2();
// const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
// const Scalar pressure = (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
// + ccPressure;
values[Indices::conti0EqIdx] = velocity * density;
// values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( pressure - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
}
#endif
return values;
}
// \}
#if ENABLEMONOLITHIC
//! Get the coupling manager
const CouplingManager& couplingManager() const
{ return *couplingManager_; }
#endif
/*!
* \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] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
if(onLeftBoundary_(globalPos))
values[Indices::pressureIdx] = deltaP_;
if(onRightBoundary_(globalPos))
values[Indices::pressureIdx] = 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
{
#if ENABLEMONOLITHIC
return couplingManager().couplingData().darcyPermeability(element, scvf);
#else
return 1e-10; // TODO transfer information or just use constant value
#endif
}
/*!
* \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
*/
Scalar alphaBJ(const SubControlVolumeFace& scvf) const
{
#if ENABLEMONOLITHIC
return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
#else
return 1.0; // TODO transfer information or just use constant value
#endif
}
/*!
* \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
*/
void calculateAnalyticalVelocityX() const
{
analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
using std::sqrt;
const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability");
static const Scalar sqrtK = sqrt(K);
const Scalar sigma = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])/sqrtK;
const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX;
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().bBoxMin()[1];
const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX;
analyticalVelocityX_[eIdx] = u;
}
}
/*!
* \brief Get the analytical velocity in x direction
*/
const std::vector<Scalar>& getAnalyticalVelocityX() const
{
if(analyticalVelocityX_.empty())
calculateAnalyticalVelocityX();
return analyticalVelocityX_;
}
#if !ENABLEMONOLITHIC
void updatePreciceDataIds()
{
pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
dataIdsWereSet_ = true;
}
#endif
// \}
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 deltaP_;
#if ENABLEMONOLITHIC
std::shared_ptr<CouplingManager> couplingManager_;
#else
precice_adapter::PreciceAdapter& couplingInterface_;
size_t pressureId_;
size_t velocityId_;
bool dataIdsWereSet_;
#endif
mutable std::vector<Scalar> analyticalVelocityX_;
};
} //end namespace
#endif // DUMUX_STOKES_SUBPROBLEM_HH
// -*- 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/partial.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include "ffproblem-reversed.hh"
#include "../../precice-adapter/include/preciceadapter.hh"
//TODO
// Helper function to put pressure on interface
template<class ElementFaceVariables, class SubControlVolumeFace>
auto velocityAtInterface(const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf)
{
const double scalarVelocity = elemFaceVars[scvf].velocitySelf();
auto velocity = scvf.unitOuterNormal();
velocity[scvf.directionIndex()] = scalarVelocity;
return velocity;
}
template<class Problem, class GridVariables, class SolutionVector>
void setInterfaceVelocities(const Problem& problem,
const GridVariables& gridVars,
const SolutionVector& sol)
{
const auto& fvGridGeometry = problem.fvGridGeometry();
auto fvGeometry = localView(fvGridGeometry);
auto elemVolVars = localView(gridVars.curGridVolVars());
auto elemFaceVars = localView(gridVars.curGridFaceVars());
auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
const auto velocityId = couplingInterface.getIdFromName( "Velocity" );
for (const auto& element : elements(fvGridGeometry.gridView()))
{
fvGeometry.bindElement(element);
elemVolVars.bindElement(element, fvGeometry, sol);
elemFaceVars.bindElement(element, fvGeometry, sol);
for (const auto& scvf : scvfs(fvGeometry))
{
if ( couplingInterface.isCoupledEntity( scvf.index() ) )
{
//TODO: What to do here?
const auto v = velocityAtInterface(elemFaceVars, scvf)[scvf.directionIndex()];
couplingInterface.writeScalarQuantityOnFace( velocityId, scvf.index(), v );
}
}
}
}
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 FreeFlowTypeTag = Properties::TTag::FreeFlowModel;
// try to create a grid (from the given grid file or the input file)
using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowTypeTag, Properties::Grid>>;
FreeFlowGridManager freeFlowGridManager;
freeFlowGridManager.init("FreeFlow"); // pass parameter group
// we compute on the leaf grid view
const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView();
// create the finite volume grid geometry
using FreeFlowFVGridGeometry = GetPropType<FreeFlowTypeTag, Properties::FVGridGeometry>;
auto freeFlowFvGridGeometry = std::make_shared<FreeFlowFVGridGeometry>(freeFlowGridView);
freeFlowFvGridGeometry->update();
// the problem (initial and boundary conditions)
using FreeFlowProblem = GetPropType<FreeFlowTypeTag, Properties::Problem>;
auto freeFlowProblem = std::make_shared<FreeFlowProblem>(freeFlowFvGridGeometry);
// the solution vector
GetPropType<FreeFlowTypeTag, Properties::SolutionVector> sol;
sol[FreeFlowFVGridGeometry::cellCenterIdx()].resize(freeFlowFvGridGeometry->numCellCenterDofs());
sol[FreeFlowFVGridGeometry::faceIdx()].resize(freeFlowFvGridGeometry->numFaceDofs());
// Initialize preCICE.Tell preCICE about:
// - Name of solver
// - What rank of how many ranks this instance is
// Configure preCICE. For now the config file is hardcoded.
//couplingInterface.createInstance( "FreeFlow", mpiHelper.rank(), mpiHelper.size() );
std::string preciceConfigFilename = "precice-config.xml";
if (argc == 3)