From 081ee98c6827f4834277a19b8f232a7158f38e27 Mon Sep 17 00:00:00 2001
From: Alexander Jaust <alexander.jaust@ipvs.uni-stuttgart.de>
Date: Tue, 28 May 2019 14:53:57 +0200
Subject: [PATCH] [coupling-ff-pm] Fix coupling with Kilian - Write routines to
extrapolate pressure. - Solvers use data from preCICE now - Some hacks to
make it work for now
---
appl/coupling-ff-pm/iterative/main_ff.cc | 54 +++++++++++++++++++--
appl/coupling-ff-pm/iterative/main_pm.cc | 44 +++++++++--------
appl/coupling-ff-pm/iterative/params.input | 18 ++-----
appl/coupling-ff-pm/monolithic/ffproblem.hh | 18 +++----
appl/coupling-ff-pm/monolithic/pmproblem.hh | 12 ++---
5 files changed, 88 insertions(+), 58 deletions(-)
diff --git a/appl/coupling-ff-pm/iterative/main_ff.cc b/appl/coupling-ff-pm/iterative/main_ff.cc
index c7f3942..31e5f0b 100644
--- a/appl/coupling-ff-pm/iterative/main_ff.cc
+++ b/appl/coupling-ff-pm/iterative/main_ff.cc
@@ -52,6 +52,50 @@
//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.writeQuantityOnFace( velocityId, scvf.index(), v );
+ }
+ }
+ }
+
+}
+
+
int main(int argc, char** argv) try
{
using namespace Dumux;
@@ -166,7 +210,8 @@ int main(int argc, char** argv) try
{
//TODO
// couplingInterface.writeQuantityVector( pressureId );
- couplingInterface.writeQuantityToOtherSolver( pressureId );
+ setInterfaceVelocities( *freeFlowProblem, *freeFlowGridVariables, sol );
+ couplingInterface.writeQuantityToOtherSolver( velocityId );
couplingInterface.announceInitialDataWritten();
}
couplingInterface.initializeData();
@@ -197,15 +242,14 @@ int main(int argc, char** argv) try
}
// TODO
- couplingInterface.readQuantityFromOtherSolver( velocityId );
+ couplingInterface.readQuantityFromOtherSolver( pressureId );
// solve the non-linear system
nonLinearSolver.solve(sol);
// TODO
- // FILL DATA vector
- // couplingInterface.writeQuantityVector( pressureId );
- couplingInterface.writeQuantityToOtherSolver( pressureId );
+ setInterfaceVelocities( *freeFlowProblem, *freeFlowGridVariables, sol );
+ couplingInterface.writeQuantityToOtherSolver( velocityId );
const double preciceDt = couplingInterface.advance( dt );
dt = std::min( preciceDt, dt );
diff --git a/appl/coupling-ff-pm/iterative/main_pm.cc b/appl/coupling-ff-pm/iterative/main_pm.cc
index 0eb1935..225ed24 100644
--- a/appl/coupling-ff-pm/iterative/main_pm.cc
+++ b/appl/coupling-ff-pm/iterative/main_pm.cc
@@ -32,6 +32,7 @@
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
+#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
@@ -56,14 +57,17 @@
/*!
* \brief Returns the pressure at the interface using Darcy's law for reconstruction
*/
- template<class Problem, class Element, class SubControlVolumeFace, class VolumeVariables, class Scalar>
- Scalar pressureAtInterface(const Problem& problem,
+ template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class SubControlVolumeFace>
+ auto pressureAtInterface(const Problem& problem,
const Element& element,
- const SubControlVolumeFace& scvf,
- const VolumeVariables& volVars,
- const Scalar boundaryFlux)
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf)
{
- const Scalar cellCenterPressure = volVars.pressure();
+ using Scalar = typename ElementVolumeVariables::VolumeVariables::PrimaryVariables::value_type;
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+
+ const Scalar boundaryFlux = problem.neumann(element, fvGeometry, elemVolVars, scvf);
// v = -K/mu * (gradP + rho*g)
auto velocity = scvf.unitOuterNormal();
@@ -75,34 +79,32 @@
const auto K = volVars.permeability();
// v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
-
- const auto alpha = vtmv(scvf.unitOuterNormal(), K, problem().spatialParams.gravity(scvf.center()));
+ const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
auto distanceVector = scvf.center() - element.geometry().center();
distanceVector /= distanceVector.two_norm2();
- const Scalar ti = vtmv(distanceVector, K, scvf.unitOuterNormal());
+ const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
return (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
+ ccPressure;
}
template<class Problem, class GridVariables, class SolutionVector>
- void setBoundaryHeatFluxes(const Problem& problem,
+ void setInterfacePressures(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 pressureId = couplingInterface.getIdFromName( "Pressure" );
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))
{
@@ -110,7 +112,8 @@
if ( couplingInterface.isCoupledEntity( scvf.index() ) )
{
//TODO: What to do here?
-// couplingInterface.pressureAtInterface( problem, element, scvf, elemVolVars, ??? );
+ const double p = pressureAtInterface(problem, element, fvGridGeometry, elemVolVars, scvf);
+ couplingInterface.writeQuantityOnFace( pressureId, scvf.index(), p );
}
}
}
@@ -186,7 +189,7 @@ int main(int argc, char** argv) try
{
static constexpr auto eps = 1e-7;
const auto& pos = scvf.center();
- if (pos[1] < darcyFvGridGeometry->bBoxMin()[1] + eps)
+ if (pos[1] > darcyFvGridGeometry->bBoxMax()[1] - eps)
{
if (pos[0] > xMin - eps && pos[0] < xMax + eps)
{
@@ -219,7 +222,7 @@ int main(int argc, char** argv) try
// intialize the vtk output module
const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
- VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol, darcyName);
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol, darcyProblem->name());
using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
@@ -230,7 +233,8 @@ int main(int argc, char** argv) try
{
//TODO
//couplingInterface.writeQuantityVector(velocityId);
- couplingInterface.writeQuantityToOtherSolver( velocityId );
+ setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+ couplingInterface.writeQuantityToOtherSolver( pressureId );
couplingInterface.announceInitialDataWritten();
}
couplingInterface.initializeData();
@@ -260,14 +264,12 @@ int main(int argc, char** argv) try
}
// TODO
- couplingInterface.readQuantityFromOtherSolver( pressureId );
+ couplingInterface.readQuantityFromOtherSolver( velocityId );
// solve the non-linear system
nonLinearSolver.solve(sol);
- // TODO
- // FILL DATA vector
- // couplingInterface.writeQuantityVector( pressureId );
- couplingInterface.writeQuantityToOtherSolver( velocityId );
+ setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+ couplingInterface.writeQuantityToOtherSolver( pressureId );
const double preciceDt = couplingInterface.advance( dt );
dt = std::min( preciceDt, dt );
diff --git a/appl/coupling-ff-pm/iterative/params.input b/appl/coupling-ff-pm/iterative/params.input
index 4117b7f..ab27df7 100644
--- a/appl/coupling-ff-pm/iterative/params.input
+++ b/appl/coupling-ff-pm/iterative/params.input
@@ -1,15 +1,6 @@
- # for dune-subgrid
- [Grid]
-Positions0 = 0 1
-Positions1 = 0 0.2 0.3 0.65
-Cells0 = 100
-Cells1 = 10 50 18
-Baseline = 0.25 # [m]
-Amplitude = 0.04 # [m]
-Offset = 0.5 # [m]
-Scaling = 0.2 #[m]
-
-[Stokes.Grid]
+
+
+[FreeFlow.Grid]
Verbosity = true
Positions0 = 0.0 1.0
Positions1 = 1.0 2.0
@@ -25,13 +16,14 @@ Cells0 = 20
Cells1 = 20
Grading1 = 1
-[Stokes.Problem]
+[FreeFlow.Problem]
Name = stokes
PressureDifference = 1e-9
EnableInertiaTerms = false
[Darcy.Problem]
Name = darcy
+InitialP = 0.5e-9
[Darcy.SpatialParams]
Permeability = 1e-6 # m^2
diff --git a/appl/coupling-ff-pm/monolithic/ffproblem.hh b/appl/coupling-ff-pm/monolithic/ffproblem.hh
index d336848..130bcc1 100644
--- a/appl/coupling-ff-pm/monolithic/ffproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/ffproblem.hh
@@ -24,7 +24,7 @@
#define DUMUX_STOKES_SUBPROBLEM_HH
#ifndef ENABLEMONOLITHIC
-#define ENABLEMONOLITHIC 1
+#define ENABLEMONOLITHIC 0
#endif
#include <dune/grid/yaspgrid.hh>
@@ -122,7 +122,7 @@ public:
: ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
#else
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
- : ParentType(fvGridGeometry, "Stokes"),
+ : ParentType(fvGridGeometry, "FreeFlow"),
eps_(1e-6),
couplingInterface_(precice_adapter::PreciceAdapter::getInstance() ),
pressureId_(0),
@@ -202,9 +202,9 @@ public:
if ( couplingInterface_.isCoupledEntity(faceId) )
{
//TODO What do I want to do here?
- // values.setCouplingNeumann(Indices::conti0EqIdx);
- // values.setCouplingNeumann(Indices::momentumYBalanceIdx);
- // values.setBJS(Indices::momentumXBalanceIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
}
#endif
@@ -254,14 +254,8 @@ public:
{
const Scalar density = 1000; // TODO how to handle compressible fluids?
values[Indices::conti0EqIdx] = elemFaceVars[scvf].velocitySelf() * scvf.directionSign() * density;
- values[Indices::momentumYBalanceIdx] = couplingInterface_.getQuantityOnFace( pressureId_, faceId );
+ values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
}
- // if(/*preCICE*/)
- // {
- // const Scalar density = 1000; // TODO how to handle compressible fluids?
- // values[Indices::conti0EqIdx] = elemFaceVars[scvf].velocitySelf() * scvf.directionSign() * density;
- // values[Indices::momentumYBalanceIdx] = /*pressure from Darcy*/
- // }
#endif
return values;
diff --git a/appl/coupling-ff-pm/monolithic/pmproblem.hh b/appl/coupling-ff-pm/monolithic/pmproblem.hh
index bbec14c..7ed57e3 100644
--- a/appl/coupling-ff-pm/monolithic/pmproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/pmproblem.hh
@@ -25,7 +25,7 @@
#define DUMUX_DARCY_SUBPROBLEM_HH
#ifndef ENABLEMONOLITHIC
-#define ENABLEMONOLITHIC 1
+#define ENABLEMONOLITHIC 0
#endif
#include <dune/grid/yaspgrid.hh>
@@ -219,12 +219,9 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
const auto faceId = scvf.index();
if ( couplingInterface_.isCoupledEntity(faceId) )
{
- values[Indices::conti0EqIdx] = - couplingInterface_.getQuantityOnFace( velocityId_, faceId );
+ const Scalar density = 1000.;
+ values[Indices::conti0EqIdx] = density * couplingInterface_.getQuantityOnFace( velocityId_, faceId );
}
- // if (/*preCICE*/)
- // {
- // values[Indices::conti0EqIdx] = /*mass flux from stokes*/;
- // }
#endif
return values;
}
@@ -263,7 +260,8 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
*/
PrimaryVariables initial(const Element &element) const
{
- return PrimaryVariables(0.0);
+ static const Scalar p = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialP");
+ return PrimaryVariables(p);
}
// \}
--
GitLab