diff --git a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
index 97390739a2de1c352e551fae7f35c70dca01f6cf..668305d1fc8ac325a80f2bad43e17d986ad99842 100644
--- a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
@@ -202,10 +202,12 @@ public:
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.setCouplingNeumann(Indices::conti0EqIdx);
+ // values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+
+// values.setNeumann(Indices::conti0EqIdx);
+// values.setNeumann(Indices::momentumYBalanceIdx);
values.setBJS(Indices::momentumXBalanceIdx);
}
#endif
@@ -218,10 +220,21 @@ public:
*
* \param globalPos The global position
*/
- PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ using ParentType::dirichlet;
+ PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
PrimaryVariables values(0.0);
- values = initialAtPos(globalPos);
+ values = initialAtPos(scvf.center());
+
+ const auto faceId = scvf.index();
+ if( couplingInterface_.isCoupledEntity( faceId ) )
+ {
+ values[Indices::velocityYIdx] =
+ couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+ }
+
+
+
return values;
}
@@ -254,28 +267,11 @@ public:
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 Scalar density = 1000; // TODO how to handle compressible fluids?
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]) ;
-// values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( pressure - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
+ const Scalar density_ = volVars.density();
+ values[Indices::conti0EqIdx] = density * elemFaceVars[scvf].velocitySelf() * scvf.directionSign();
+ values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
}
#endif
diff --git a/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc b/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
index d0d3c6b5562943a2320e77e54c52f1886e6e2cca..3f59d6b97f7b4b8500b808a6537c8d16cf976a4f 100644
--- a/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
+++ b/appl/coupling-ff-pm/iterative-reversed/main_ff-reversed.cc
@@ -62,6 +62,62 @@ auto velocityAtInterface(const ElementFaceVariables& elemFaceVars,
return velocity;
}
+template<class FluxVariables,
+ class Problem,
+ class Element,
+ class SubControlVolumeFace,
+ class FVElementGeometry,
+ class ElementVolumeVariables,
+ class ElementFaceVariables,
+ class ElementFluxVariablesCache>
+auto pressureAtInterface(const Problem& problem,
+ const Element& element,
+ const SubControlVolumeFace& scvf,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
+{
+ FluxVariables fluxVars;
+ return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache()) /scvf.area();
+}
+
+template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+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 elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto pressureId = couplingInterface.getIdFromName( "Pressure" );
+
+ for (const auto& element : elements(fvGridGeometry.gridView()))
+ {
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, sol);
+ elemFaceVars.bind(element, fvGeometry, sol);
+ elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+
+ if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+ const auto p = pressureAtInterface<FluxVariables>(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+ couplingInterface.writeScalarQuantityOnFace( pressureId, scvf.index(), p );
+ }
+ }
+ }
+
+}
+
+
template<class Problem, class GridVariables, class SolutionVector>
void setInterfaceVelocities(const Problem& problem,
const GridVariables& gridVars,
@@ -154,6 +210,7 @@ int main(int argc, char** argv) try
if (dim != int(FreeFlowFVGridGeometry::GridView::dimension))
DUNE_THROW(Dune::InvalidStateException, "Dimensions do not match");
+
// GET mesh corodinates
const double xMin = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0")[0];
const double xMax = getParamFromGroup<std::vector<double>>("Darcy", "Grid.Positions0").back();
@@ -206,11 +263,22 @@ int main(int argc, char** argv) try
freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x");
freeFlowVtkWriter.write(0.0);
+ using FluxVariables = GetPropType<FreeFlowTypeTag, Properties::FluxVariables>;
+
if ( couplingInterface.hasToWriteInitialData() )
{
//TODO
// couplingInterface.writeQuantityVector( pressureId );
- setInterfacePressures( *freeFlowProblem, *freeFlowGridVariables, sol );
+
+ setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+ //For testing
+ {
+ std::cout << "Pressures to be sent to pm" << std::endl;
+ const auto p = couplingInterface.getQuantityVector( pressureId );
+ for (size_t i = 0; i < p.size(); ++i) {
+ std::cout << "p[" << i << "]=" <<p[i] << std::endl;
+ }
+ }
couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
couplingInterface.announceInitialDataWritten();
}
@@ -232,6 +300,8 @@ int main(int argc, char** argv) try
auto dt = preciceDt;
auto sol_checkpoint = sol;
+ double vtkTime = 1.0;
+
while ( couplingInterface.isCouplingOngoing() )
{
if ( couplingInterface.hasToWriteIterationCheckpoint() )
@@ -244,25 +314,41 @@ int main(int argc, char** argv) try
// TODO
couplingInterface.readScalarQuantityFromOtherSolver( velocityId );
// For testing
-// {
-// const auto v = couplingInterface.getQuantityVector( velocityId );
-// const double sum = std::accumulate( v.begin(), v.end(), 0. );
-// std::cout << "Sum of velocities over boundary to pm: \n" << sum << std::endl;
-// }
+ {
+ const auto v = couplingInterface.getQuantityVector( velocityId );
+ const double sum = std::accumulate( v.begin(), v.end(), 0. );
+ std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
+ }
// solve the non-linear system
nonLinearSolver.solve(sol);
// TODO
- setInterfacePressures( *freeFlowProblem, *freeFlowGridVariables, sol );
+ setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+ // For testing
+ {
+ const auto p = couplingInterface.getQuantityVector( pressureId );
+ const double sum = std::accumulate( p.begin(), p.end(), 0. );
+ std::cout << "Pressures to be sent to pm" << std::endl;
+ for (size_t i = 0; i < p.size(); ++i) {
+ std::cout << "p[" << i << "]=" << p[i] << std::endl;
+ }
+ std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
+ }
couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+
+ //Read checkpoint
+ freeFlowVtkWriter.write(vtkTime);
+ vtkTime += 1.;
const double preciceDt = couplingInterface.advance( dt );
dt = std::min( preciceDt, dt );
if ( couplingInterface.hasToReadIterationCheckpoint() )
{
- //Read checkpoint
+// //Read checkpoint
+// freeFlowVtkWriter.write(vtkTime);
+// vtkTime += 1.;
sol = sol_checkpoint;
freeFlowGridVariables->update(sol);
freeFlowGridVariables->advanceTimeStep();
@@ -272,7 +358,7 @@ int main(int argc, char** argv) try
else // coupling successful
{
// write vtk output
- freeFlowVtkWriter.write(1.0);
+ freeFlowVtkWriter.write(vtkTime);
}
}
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 a49cfbe3ceb5fc83ce4785e6049253e3bae854e5..b860427f148b5c1c49fe56256779b1f57eda89fc 100644
--- a/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
+++ b/appl/coupling-ff-pm/iterative-reversed/main_pm-reversed.cc
@@ -69,14 +69,15 @@
const Scalar boundaryFlux = problem.neumann(element, fvGeometry, elemVolVars, scvf);
+ const auto K = volVars.permeability();
+ const Scalar ccPressure = volVars.pressure();
+ const Scalar mobility = volVars.mobility();
+ const Scalar density = volVars.density();
+
// v = -K/mu * (gradP + rho*g)
auto velocity = scvf.unitOuterNormal();
velocity *= boundaryFlux; // TODO check sign
-
- const Scalar ccPressure = volVars.pressure();
- const Scalar mobility = volVars.mobility();
- const Scalar density = volVars.density();
- const auto K = volVars.permeability();
+ velocity /= density;
// v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
@@ -120,6 +121,60 @@
}
+ /*!
+ * \brief Returns the pressure 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 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 elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+ auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+ const auto velocityId = couplingInterface.getIdFromName( "Velocity" );
+
+ 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 ( couplingInterface.isCoupledEntity( scvf.index() ) )
+ {
+ //TODO: What to do here?
+ const double v = velocityAtInterface<FluxVariables>(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ couplingInterface.writeScalarQuantityOnFace( velocityId, scvf.index(), v );
+ }
+ }
+ }
+
+ }
+
int main(int argc, char** argv) try
{
using namespace Dumux;
@@ -228,12 +283,20 @@ int main(int argc, char** argv) try
GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
darcyVtkWriter.write(0.0);
-
+ using FluxVariables = GetPropType<DarcyTypeTag, Properties::FluxVariables>;
if ( couplingInterface.hasToWriteInitialData() )
{
//TODO
//couplingInterface.writeQuantityVector(velocityId);
- setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+ setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+ // For testing
+ {
+ const auto v = couplingInterface.getQuantityVector( velocityId );
+ std::cout << "velocities to be sent to ff" << std::endl;
+ for (size_t i = 0; i < v.size(); ++i) {
+ std::cout << "v[" << i << "]=" << v[i] << std::endl;
+ }
+ }
couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
couplingInterface.announceInitialDataWritten();
}
@@ -254,6 +317,8 @@ int main(int argc, char** argv) try
auto dt = preciceDt;
auto sol_checkpoint = sol;
+ double vtkTime = 1.0;
+
while ( couplingInterface.isCouplingOngoing() )
{
if ( couplingInterface.hasToWriteIterationCheckpoint() )
@@ -264,18 +329,32 @@ int main(int argc, char** argv) try
}
// TODO
- couplingInterface.readScalarQuantityFromOtherSolver( velocityId );
+ couplingInterface.readScalarQuantityFromOtherSolver( pressureId );
// For testing
{
- const auto v = couplingInterface.getQuantityVector( velocityId );
- const double sum = std::accumulate( v.begin(), v.end(), 0. );
- std::cout << "Sum of fluxes over boundary to pm: \n" << sum << std::endl;
+ const auto p = couplingInterface.getQuantityVector( pressureId );
+ const double sum = std::accumulate( p.begin(), p.end(), 0. );
+ std::cout << "Sum of pressures over boundary to ff: \n" << sum << std::endl;
+ std::cout << "Pressure received from ff" << std::endl;
+ for (size_t i = 0; i < p.size(); ++i) {
+ std::cout << "p[" << i << "]=" << p[i] << std::endl;
+ }
}
// solve the non-linear system
nonLinearSolver.solve(sol);
- setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
- couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+ setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+ // For testing
+ {
+ const auto v = couplingInterface.getQuantityVector( velocityId );
+ const double sum = std::accumulate( v.begin(), v.end(), 0. );
+ std::cout << "Velocities to be sent to ff" << std::endl;
+ for (size_t i = 0; i < v.size(); ++i) {
+ std::cout << "v[" << i << "]=" << v[i] << std::endl;
+ }
+ std::cout << "Sum of velocities over boundary to ff: \n" << sum << std::endl;
+ }
+ couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
const double preciceDt = couplingInterface.advance( dt );
dt = std::min( preciceDt, dt );
@@ -283,6 +362,8 @@ int main(int argc, char** argv) try
if ( couplingInterface.hasToReadIterationCheckpoint() )
{
//Read checkpoint
+ darcyVtkWriter.write(vtkTime);
+ vtkTime += 1.;
sol = sol_checkpoint;
darcyGridVariables->update(sol);
darcyGridVariables->advanceTimeStep();
@@ -292,7 +373,7 @@ int main(int argc, char** argv) try
else // coupling successful
{
// write vtk output
- darcyVtkWriter.write(1.0);
+ darcyVtkWriter.write(vtkTime);
}
}
diff --git a/appl/coupling-ff-pm/iterative-reversed/params.input b/appl/coupling-ff-pm/iterative-reversed/params.input
index c0fa10ab8a8eadbe9163164568179b5e75dabf91..2b8584a8de67878eb562f9c93f9656f850b1640f 100644
--- a/appl/coupling-ff-pm/iterative-reversed/params.input
+++ b/appl/coupling-ff-pm/iterative-reversed/params.input
@@ -4,7 +4,7 @@
Verbosity = true
Positions0 = 0.0 1.0
Positions1 = 1.0 2.0
-Cells0 = 20
+Cells0 = 57
Cells1 = 100
Grading1 = 1
diff --git a/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
index dfca6d5cf12703aeb9fcd6de21d010a8d1fbb65d..40e7d7f153f3f4f50fa84ee6cb9a67c98e073694 100644
--- a/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
@@ -171,16 +171,9 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
#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.setCouplingNeumann(Indices::conti0EqIdx);
- values.setCouplingNeumann(Indices::momentumYBalanceIdx);
- values.setBJS(Indices::momentumXBalanceIdx);
- }
+ values.setAllDirichlet();
#endif
return values;
}
@@ -199,6 +192,10 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
PrimaryVariables values(0.0);
values = initial(element);
+ const auto faceId = scvf.index();
+ if ( couplingInterface_.isCoupledEntity(faceId) )
+ values = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+
return values;
}
@@ -226,34 +223,13 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
#else
- zassert( dataIdsWereSet_ );
+ assert( dataIdsWereSet_ );
const auto faceId = scvf.index();
- if( couplingInterface_.isCoupledEntity( faceId ) )
+ if ( couplingInterface_.isCoupledEntity(faceId) )
{
- const Scalar densityOnFace = 1000; // TODO how to handle compressible fluids?
-
- const Scalar pressure = couplingInterface_.getScalarQuantityOnFace( pressureId_, 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;
-// mobility * ( (pressure - ccPressure) * ti - + density * alpha ) / scvf.unitOuterNormal()= velocity ;
- const Scalar velocity = mobility * ( (pressure - ccPressure) * ti - + density * alpha ) / scvf.unitOuterNormal();
-
- values[Indices::conti0EqIdx] = velocity * scvf.directionSign() * densityOnFace;
- values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
+ const Scalar density = 1000.;
+ values[Indices::conti0EqIdx] = density * couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+ std::cout << "pm: values[Indices::conti0EqIdx] = " << values << std::endl;
}
#endif
return values;
diff --git a/appl/coupling-ff-pm/monolithic/params.input b/appl/coupling-ff-pm/monolithic/params.input
index 455dd8fbc09300fdc031e51ed0e4a2cb35e34f96..93eb547d328d24313a50df35e2e092ee72f91a51 100644
--- a/appl/coupling-ff-pm/monolithic/params.input
+++ b/appl/coupling-ff-pm/monolithic/params.input
@@ -27,8 +27,8 @@ Grading1 = 1
[Stokes.Problem]
Name = stokes
-PressureDifference = 1e2
-EnableInertiaTerms = true
+PressureDifference = 1e-9
+EnableInertiaTerms = false
[Darcy.Problem]
Name = darcy