From 98f5d422c04d2f42b3203433c4ac505ef3159f05 Mon Sep 17 00:00:00 2001
From: Alexander Jaust <alexander.jaust@ipvs.uni-stuttgart.de>
Date: Tue, 22 Oct 2019 11:17:39 +0200
Subject: [PATCH] WIP: precice-dumux coupling
Updated code to current version on my laptop. Need to cleanup up code. It has a lot of debugging output.
Good news: Partitioned results no looke very similar to monolithic result after Kilian has changed the boundary conditions in the monolithic case
---
.../iterative-reversed/ffproblem-reversed.hh | 22 ++--
.../iterative-reversed/main_pm-reversed.cc | 5 +-
.../iterative-reversed/params.input | 2 +-
appl/coupling-ff-pm/monolithic/ffproblem.hh | 4 +-
appl/coupling-ff-pm/monolithic/main.cc | 120 ++++++++++++++++++
appl/coupling-ff-pm/monolithic/params.input | 2 +-
appl/coupling-ff-pm/monolithic/pmproblem.hh | 4 +-
7 files changed, 144 insertions(+), 15 deletions(-)
diff --git a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
index 668305d..0a9ef2a 100644
--- a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
@@ -175,32 +175,29 @@ public:
BoundaryTypes values;
const auto& globalPos = scvf.dofPosition();
+ const auto faceId = scvf.index();
// 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))
+ else 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) )
+
+ else if ( couplingInterface_.isCoupledEntity(faceId) )
{
+ // // TODO do preCICE stuff in analogy to heat transfer
+ assert( dataIdsWereSet_ );
//TODO What do I want to do here?
// values.setCouplingNeumann(Indices::conti0EqIdx);
// values.setCouplingNeumann(Indices::momentumYBalanceIdx);
@@ -211,6 +208,11 @@ public:
values.setBJS(Indices::momentumXBalanceIdx);
}
#endif
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
return values;
}
diff --git a/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc b/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
index 1e1881a..2178fbd 100644
--- a/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
+++ b/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
@@ -28,6 +28,8 @@
#include <iostream>
#include <string>
+bool printstuff = false;
+
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
@@ -109,6 +111,7 @@
fvGeometry.bindElement(element);
elemVolVars.bindElement(element, fvGeometry, sol);
+ //sstd::cout << "Pressure by reconstruction" << std::endl;
for (const auto& scvf : scvfs(fvGeometry))
{
@@ -393,7 +396,7 @@ int main(int argc, char** argv) try
std::cout << "| v[" << i << "]=" << v[i] << std::endl;
}
}
- couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+ couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
couplingInterface.announceInitialDataWritten();
}
couplingInterface.initializeData();
diff --git a/appl/coupling-ff-pm/iterative-reversed/params.input b/appl/coupling-ff-pm/iterative-reversed/params.input
index 00515fb..31a7880 100644
--- a/appl/coupling-ff-pm/iterative-reversed/params.input
+++ b/appl/coupling-ff-pm/iterative-reversed/params.input
@@ -5,7 +5,7 @@ Verbosity = true
Positions0 = 0.0 1.0
Positions1 = 1.0 2.0
Cells0 = 20
-Cells1 = 100
+Cells1 = 20
Grading1 = 1
[Darcy.Grid]
diff --git a/appl/coupling-ff-pm/monolithic/ffproblem.hh b/appl/coupling-ff-pm/monolithic/ffproblem.hh
index 668305d..c5b642c 100644
--- a/appl/coupling-ff-pm/monolithic/ffproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/ffproblem.hh
@@ -226,13 +226,15 @@ public:
PrimaryVariables values(0.0);
values = initialAtPos(scvf.center());
+#if ENABLEMONOLITHIC
+#else
const auto faceId = scvf.index();
if( couplingInterface_.isCoupledEntity( faceId ) )
{
values[Indices::velocityYIdx] =
couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
}
-
+#endif
return values;
diff --git a/appl/coupling-ff-pm/monolithic/main.cc b/appl/coupling-ff-pm/monolithic/main.cc
index 6d228a9..01af755 100644
--- a/appl/coupling-ff-pm/monolithic/main.cc
+++ b/appl/coupling-ff-pm/monolithic/main.cc
@@ -69,6 +69,106 @@ struct CouplingManager<TypeTag, TTag::DarcyOneP>
} // end namespace Properties
} // end namespace Dumux
+
+/*!
+ * \brief Returns the velocity at the interface using Darcy's law for reconstruction
+ */
+template<class FluxVariables, class Problem, class Element, class FVElementGeometry,
+ class ElementVolumeVariables, class SubControlVolumeFace,
+ class ElementFluxVariablesCache>
+auto velocityAtInterface(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+{
+ const int phaseIdx = 0;
+ FluxVariables fluxVars;
+ fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); };
+ const auto scalarVelocity = fluxVars.advectiveFlux(phaseIdx, upwindTerm)/scvf.area();
+ return scalarVelocity;
+}
+
+template<class FluxVariables, class CouplingManager, class Problem, class GridVariables, class SolutionVector>
+void writeVelocitiesOnInterfaceToFile( const std::string& filename,
+ const CouplingManager& couplingManager,
+ const Problem& problem,
+ const GridVariables& gridVars,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ auto elemVolVars = localView(gridVars.curGridVolVars());
+ auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ ofs << "velocityY" << "\n";
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if (couplingManager.isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < 2; ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double v = couplingManager.couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
+ ofs << v / 1e3 << "\n";
+ }
+ }
+ }
+
+ ofs.close();
+}
+
+template<class CouplingManager, class Problem, class SolutionVector>
+void writeStokesVelocitiesOnInterfaceToFile( const std::string& filename,
+ const CouplingManager& couplingManager,
+ const Problem& problem,
+ const SolutionVector& sol)
+{
+ const auto& fvGridGeometry = problem.fvGridGeometry();
+ auto fvGeometry = localView(fvGridGeometry);
+ using FVGridGeometry = std::decay_t<decltype (fvGridGeometry)>;
+
+ std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
+ ofs << "x,y,";
+ ofs << "velocityY" << "\n";
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if (couplingManager.isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ const auto& pos = scvf.center();
+ for (int i = 0; i < 2; ++i )
+ {
+ ofs << pos[i] << ",";
+ }
+ const double v = sol[scvf.dofIndex()];
+ ofs << v << "\n";
+ }
+ }
+ }
+
+ ofs.close();
+}
+
int main(int argc, char** argv) try
{
using namespace Dumux;
@@ -241,6 +341,26 @@ int main(int argc, char** argv) try
stokesVtkWriter.write(1.0);
darcyVtkWriter.write(1.0);
+ {
+ using FluxVariables = GetPropType<DarcyTypeTag, Properties::FluxVariables>;
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + darcyProblem->name() + "-interface-velocity";
+ writeVelocitiesOnInterfaceToFile<FluxVariables>( filename,
+ *couplingManager,
+ *darcyProblem,
+ *darcyGridVariables,
+ sol[darcyIdx] );
+ }
+
+ //TODO make freeflow
+ {
+ const std::string filename = getParam<std::string>("Problem.Name") + "-" + stokesProblem->name() + "-interface-velocity";
+ writeStokesVelocitiesOnInterfaceToFile( filename,
+ *couplingManager,
+ *stokesProblem,
+ sol[stokesFaceIdx] );
+ }
+
+
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
diff --git a/appl/coupling-ff-pm/monolithic/params.input b/appl/coupling-ff-pm/monolithic/params.input
index 843e4fb..053c1da 100644
--- a/appl/coupling-ff-pm/monolithic/params.input
+++ b/appl/coupling-ff-pm/monolithic/params.input
@@ -14,7 +14,7 @@ Verbosity = true
Positions0 = 0.0 1.0
Positions1 = 1.0 2.0
Cells0 = 20
-Cells1 = 100
+Cells1 = 20
Grading1 = 1
[Darcy.Grid]
diff --git a/appl/coupling-ff-pm/monolithic/pmproblem.hh b/appl/coupling-ff-pm/monolithic/pmproblem.hh
index 40e7d7f..2506477 100644
--- a/appl/coupling-ff-pm/monolithic/pmproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/pmproblem.hh
@@ -192,10 +192,12 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
PrimaryVariables values(0.0);
values = initial(element);
+#if ENABLEMONOLITHIC
+#else
const auto faceId = scvf.index();
if ( couplingInterface_.isCoupledEntity(faceId) )
values = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
-
+#endif
return values;
}
--
GitLab