diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
index 5e9271041430048fcec8f4cce81b441584b9e6fe..d40b0679b2e784de9fb9d9b9c3463e05097494a7 100644
--- a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/CMakeLists.txt
@@ -1,29 +1,11 @@
# 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_mono_flux
-# SOURCES main_ff-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
-# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
-
-#dune_add_test(NAME test_pm_reversed_mono_flux
-# SOURCES main_pm-reversed.cc ../precice/preciceadapter.cc ../precice/dumuxpreciceindexwrapper.cc
-# COMPILE_DEFINITIONS ENABLEMONOLITHIC=0)
-
-add_executable(test_ff_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_ff-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
-add_executable(test_pm_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_pm-reversed.cc ../../precice-adapter/src/preciceadapter.cc ../../precice-adapter/src/dumuxpreciceindexwrapper.cc)
+add_executable(test_ff_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_ff-mono-flux.cc )
+add_executable(test_pm_reversed_mono_flux-no-precice EXCLUDE_FROM_ALL main_pm-mono-flux.cc )
target_compile_definitions(test_ff_reversed_mono_flux-no-precice PUBLIC "ENABLEMONOLITHIC=0")
target_compile_definitions(test_pm_reversed_mono_flux-no-precice PUBLIC "ENABLEMONOLITHIC=0")
-target_compile_definitions(test_ff_reversed_mono_flux-no-precice PRIVATE BOOST_ALL_DYN_LINK)
-target_link_libraries(test_ff_reversed_mono_flux-no-precice PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
-target_include_directories(test_ff_reversed_mono_flux-no-precice PRIVATE ${Boost_INCLUDE_DIRS})
-
-target_compile_definitions(test_pm_reversed_mono_flux-no-precice PRIVATE BOOST_ALL_DYN_LINK)
-target_link_libraries(test_pm_reversed_mono_flux-no-precice PRIVATE ${Boost_LIBRARIES} ${PRECICE_LIB})
-target_include_directories(test_pm_reversed_mono_flux-no-precice PRIVATE ${Boost_INCLUDE_DIRS})
# add a symlink for each input file
add_input_file_links()
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-mono-flux.hh
similarity index 91%
rename from appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh
rename to appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-mono-flux.hh
index 0a9ef2af63083faefd885bd4ebe6fee7ab0b73c9..b6a8f94d5d275aae4fad29b0636168a932c2119d 100644
--- a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/ffproblem-mono-flux.hh
@@ -27,6 +27,8 @@
#define ENABLEMONOLITHIC 0
#endif
+#include<algorithm>
+
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
@@ -36,7 +38,7 @@
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
-#include "../../precice-adapter/include/preciceadapter.hh"
+#include "monolithicdata.hh"
namespace Dumux
{
@@ -123,16 +125,10 @@ public:
#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)
+ eps_(1e-6)
#endif
{
deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
-// pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
-// velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
}
/*!
@@ -177,6 +173,12 @@ public:
const auto& globalPos = scvf.dofPosition();
const auto faceId = scvf.index();
+
+ const bool isCoupledEntity
+ = std::find( MonolithicSolution::velocityFaceIdx.begin(),
+ MonolithicSolution::velocityFaceIdx.end(), faceId )
+ != MonolithicSolution::velocityFaceIdx.end();
+
// left/right wall
if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
{
@@ -194,10 +196,8 @@ public:
}
#else
- else if ( couplingInterface_.isCoupledEntity(faceId) )
+ else if ( isCoupledEntity )
{
- // // 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);
@@ -229,10 +229,16 @@ public:
values = initialAtPos(scvf.center());
const auto faceId = scvf.index();
- if( couplingInterface_.isCoupledEntity( faceId ) )
+
+
+ const bool isCoupledEntity
+ = std::find( MonolithicSolution::velocityFaceIdx.begin(),
+ MonolithicSolution::velocityFaceIdx.end(), faceId )
+ != MonolithicSolution::velocityFaceIdx.end();
+ if( isCoupledEntity )
{
- values[Indices::velocityYIdx] =
- couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+ const int arrayIdx = (faceId-2) / 4;
+ values[Indices::velocityYIdx] = MonolithicSolution::velocity[arrayIdx];
}
@@ -258,6 +264,8 @@ public:
{
NumEqVector values(0.0);
+ throw std::runtime_error("You are not suppoed to call neumann boundary condition!");
+
#if ENABLEMONOLITHIC
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
@@ -265,15 +273,21 @@ public:
values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
}
#else
- assert( dataIdsWereSet_ );
const auto faceId = scvf.index();
- if( couplingInterface_.isCoupledEntity( faceId ) )
+
+ const bool isCoupledEntity
+ = std::find( MonolithicSolution::velocityFaceIdx.begin(),
+ MonolithicSolution::velocityFaceIdx.end(), faceId )
+ != MonolithicSolution::velocityFaceIdx.end();
+
+ if( isCoupledEntity )
{
+ const int arrayIdx = (faceId-2) / 4;
const Scalar density = 1000; // TODO how to handle compressible fluids?
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
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]) ;
+ values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( MonolithicSolution::pressure[arrayIdx] - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
}
#endif
@@ -371,16 +385,6 @@ public:
return analyticalVelocityX_;
}
-#if !ENABLEMONOLITHIC
- void updatePreciceDataIds()
- {
- pressureId_ = couplingInterface_.getIdFromName( "Pressure" );
- velocityId_ = couplingInterface_.getIdFromName( "Velocity" );
- dataIdsWereSet_ = true;
- }
-#endif
-
- // \}
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
@@ -401,10 +405,7 @@ private:
#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_;
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-mono-flux.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-mono-flux.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1ea5b151c6a56c472600a6f78f7aacf2814b5d70
--- /dev/null
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-mono-flux.cc
@@ -0,0 +1,183 @@
+// -*- 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-mono-flux.hh"
+#include "monolithicdata.hh"
+
+
+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());
+
+
+ // 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();
+ std::vector<double> coords; //( dim * vertexSize );
+ std::vector<int> coupledScvfIndices;
+
+ for (const auto& element : elements(freeFlowGridView))
+ {
+ auto fvGeometry = localView(*freeFlowFvGridGeometry);
+ fvGeometry.bindElement(element);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ static constexpr auto eps = 1e-7;
+ const auto& pos = scvf.center();
+ if (pos[1] < freeFlowFvGridGeometry->bBoxMin()[1] + eps)
+ {
+ if (pos[0] > xMin - eps && pos[0] < xMax + eps)
+ {
+ coupledScvfIndices.push_back(scvf.index());
+ std::cout << "Coupled index: " << scvf.index() << std::endl;
+ }
+ }
+ }
+ }
+
+
+ // the grid variables
+ using FreeFlowGridVariables = GetPropType<FreeFlowTypeTag, Properties::GridVariables>;
+ auto freeFlowGridVariables = std::make_shared<FreeFlowGridVariables>(freeFlowProblem, freeFlowFvGridGeometry);
+ freeFlowGridVariables->init(sol);
+
+ // intialize the vtk output module
+ StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(sol)> freeFlowVtkWriter(*freeFlowGridVariables, sol, freeFlowProblem->name());
+ GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter);
+ freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x");
+ freeFlowVtkWriter.write(0.0);
+
+ using FluxVariables = GetPropType<FreeFlowTypeTag, Properties::FluxVariables>;
+
+ // the assembler for a stationary problem
+ using Assembler = StaggeredFVAssembler<FreeFlowTypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(freeFlowProblem, freeFlowFvGridGeometry, freeFlowGridVariables);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ freeFlowVtkWriter.write(1.0);
+
+
+
+ // 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;
+}
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc
deleted file mode 100644
index 847c08d3d919dcef093b24a35c9e44a93cca7ff6..0000000000000000000000000000000000000000
--- a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_ff-reversed.cc
+++ /dev/null
@@ -1,584 +0,0 @@
-// -*- 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 "monolithicdata.hh"
-
-
-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 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" );
-
- size_t i = 0;
- 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);
- const auto p = MonolithicSolution::pressure[i];
- ++i;
- couplingInterface.writeScalarQuantityOnFace( pressureId, scvf.index(), p );
- }
- }
- }
-
-}
-
-
-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 );
- }
- }
- }
-
-}
-
-template<class Problem, class GridVariables, class SolutionVector>
-std::tuple<double,double,double> writeVelocitiesOnInterfaceToFile( const std::string& filename,
- 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());
-
- const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
-
- std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
- ofs << "x,y,";
- if ( couplingInterface.getDimensions() == 3 )
- ofs << "z,";
- ofs << "velocityY" << "\n";
-
-
- double min = std::numeric_limits<double>::max();
- double max = std::numeric_limits<double>::min();
- double sum = 0.;
- for (const auto& element : elements(fvGridGeometry.gridView()))
- {
- fvGeometry.bind(element);
- elemVolVars.bind(element, fvGeometry, sol);
- elemFaceVars.bindElement(element, fvGeometry, sol);
-
- for (const auto& scvf : scvfs(fvGeometry))
- {
-
- if ( couplingInterface.isCoupledEntity( scvf.index() ) )
- {
- const auto& pos = scvf.center();
- for (int i = 0; i < couplingInterface.getDimensions(); ++i )
- {
- ofs << pos[i] << ",";
- }
- const double v = problem.dirichlet( element, scvf )[1];
- max = std::max( v, max );
- min = std::min( v, min );
- sum += v;
- //velocityAtInterface(elemFaceVars, scvf)[scvf.directionIndex()];
- const int prec = ofs.precision();
- ofs << std::setprecision(std::numeric_limits<double>::digits10 + 1) << v << "\n";
- ofs.precision( prec );
- }
- }
- }
-
- ofs.close();
-
- return std::make_tuple(min, max, sum);
-}
-
-
-
-template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
-void writePressuresOnInterfaceToFile( const std::string& filename,
- 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());
-
- const auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
-
- std::ofstream ofs( filename+".csv", std::ofstream::out | std::ofstream::trunc);
- ofs << "x,y,";
- if ( couplingInterface.getDimensions() == 3 )
- ofs << "z,";
- ofs << "pressure" << "\n";
- 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() ) )
- {
- const auto& pos = scvf.center();
- for (int i = 0; i < couplingInterface.getDimensions(); ++i )
- {
- ofs << pos[i] << ",";
- }
- const double p = pressureAtInterface<FluxVariables>(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
- ofs << p << "\n";
- }
- }
- }
-
- ofs.close();
-}
-
-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)
-// preciceConfigFilename = argv[2];
- if (argc > 2)
- preciceConfigFilename = argv[argc-1];
-
- auto& couplingInterface =
- precice_adapter::PreciceAdapter::getInstance();
- couplingInterface.announceSolver( "FreeFlow", preciceConfigFilename,
- mpiHelper.rank(), mpiHelper.size() );
-
- const int dim = couplingInterface.getDimensions();
- std::cout << dim << " " << int(FreeFlowFVGridGeometry::GridView::dimension) << std::endl;
- 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();
- std::vector<double> coords; //( dim * vertexSize );
- std::vector<int> coupledScvfIndices;
-
- for (const auto& element : elements(freeFlowGridView))
- {
- auto fvGeometry = localView(*freeFlowFvGridGeometry);
- fvGeometry.bindElement(element);
-
- for (const auto& scvf : scvfs(fvGeometry))
- {
- static constexpr auto eps = 1e-7;
- const auto& pos = scvf.center();
- if (pos[1] < freeFlowFvGridGeometry->bBoxMin()[1] + eps)
- {
- if (pos[0] > xMin - eps && pos[0] < xMax + eps)
- {
- coupledScvfIndices.push_back(scvf.index());
- for (const auto p : pos)
- coords.push_back(p);
- }
- }
- }
- }
-
- const auto numberOfPoints = coords.size() / dim;
- const double preciceDt = couplingInterface.setMeshAndInitialize( "FreeFlowMesh",
- numberOfPoints,
- coords);
- couplingInterface.createIndexMapping( coupledScvfIndices );
-
- const auto velocityId = couplingInterface.announceQuantity( "Velocity" );
- const auto pressureId = couplingInterface.announceQuantity( "Pressure" );
-
- freeFlowProblem->updatePreciceDataIds();
-
- // apply initial solution for instationary problems
- freeFlowProblem->applyInitialSolution(sol);
-
- // the grid variables
- using FreeFlowGridVariables = GetPropType<FreeFlowTypeTag, Properties::GridVariables>;
- auto freeFlowGridVariables = std::make_shared<FreeFlowGridVariables>(freeFlowProblem, freeFlowFvGridGeometry);
- freeFlowGridVariables->init(sol);
-
- // intialize the vtk output module
- StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(sol)> freeFlowVtkWriter(*freeFlowGridVariables, sol, freeFlowProblem->name());
- GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter);
- freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x");
- freeFlowVtkWriter.write(0.0);
-
- using FluxVariables = GetPropType<FreeFlowTypeTag, Properties::FluxVariables>;
-
- if ( couplingInterface.hasToWriteInitialData() )
- {
- //TODO
-// couplingInterface.writeQuantityVector( pressureId );
-
- 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();
- }
- couplingInterface.initializeData();
-
- // the assembler for a stationary problem
- using Assembler = StaggeredFVAssembler<FreeFlowTypeTag, DiffMethod::numeric>;
- auto assembler = std::make_shared<Assembler>(freeFlowProblem, freeFlowFvGridGeometry, freeFlowGridVariables);
-
- // the linear solver
- using LinearSolver = UMFPackBackend;
- auto linearSolver = std::make_shared<LinearSolver>();
-
- // the non-linear solver
- using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver);
-
-
- auto dt = preciceDt;
- auto sol_checkpoint = sol;
-
- double vtkTime = 1.0;
- size_t iter = 0;
-
- while ( couplingInterface.isCouplingOngoing() )
- {
- if ( couplingInterface.hasToWriteIterationCheckpoint() )
- {
- //DO CHECKPOINTING
- sol_checkpoint = sol;
- couplingInterface.announceIterationCheckpointWritten();
- }
-
- // 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;
-// }
-
- // solve the non-linear system
- nonLinearSolver.solve(sol);
-
- // TODO
- 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 );
-
- //
- {
- double min = std::numeric_limits<double>::max();
- double max = std::numeric_limits<double>::min();
- double sum = 0.;
- const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-velocity-" + std::to_string(iter);
- std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile( filename,
- *freeFlowProblem,
- *freeFlowGridVariables,
- sol );
- const int prec = std::cout.precision();
- std::cout << "Velocity statistics:" << std::endl
- << std::setprecision(std::numeric_limits<double>::digits10 + 1)
- << " min: " << min << std::endl
- << " max: " << max << std::endl
- << " sum: " << sum << std::endl;
- std::cout.precision( prec );
- {
- const std::string filenameFlow="free-flow-statistics-" + std::to_string(iter);
- std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
- const auto prec = ofs.precision();
- ofs << "Velocity statistics (free flow):" << std::endl
- << std::setprecision(std::numeric_limits<double>::digits10 + 1)
- << " min: " << min << std::endl
- << " max: " << max << std::endl
- << " sum: " << sum << std::endl;
- ofs.precision( prec );
- ofs.close();
- }
- }
- {
- const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-pressure-" + std::to_string(iter);
- writePressuresOnInterfaceToFile<FluxVariables>( filename,
- *freeFlowProblem,
- *freeFlowGridVariables,
- sol );
- }
- ++iter;
-
- if ( couplingInterface.hasToReadIterationCheckpoint() )
- {
-// //Read checkpoint
-// freeFlowVtkWriter.write(vtkTime);
-// vtkTime += 1.;
- sol = sol_checkpoint;
- freeFlowGridVariables->update(sol);
- freeFlowGridVariables->advanceTimeStep();
- //freeFlowGridVariables->init(sol);
- couplingInterface.announceIterationCheckpointRead();
- }
- else // coupling successful
- {
- // write vtk output
- freeFlowVtkWriter.write(vtkTime);
- }
-
- }
- ////////////////////////////////////////////////////////////
- // finalize, print dumux message to say goodbye
- ////////////////////////////////////////////////////////////
-
- {
- double min = std::numeric_limits<double>::max();
- double max = std::numeric_limits<double>::min();
- double sum = 0.;
- const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-velocity";
- std::tie(min, max, sum) = writeVelocitiesOnInterfaceToFile( filename,
- *freeFlowProblem,
- *freeFlowGridVariables,
- sol );
- const int prec = std::cout.precision();
- std::cout << "Velocity statistics:" << std::endl
- << std::setprecision(std::numeric_limits<double>::digits10 + 1)
- << " min: " << min << std::endl
- << " max: " << max << std::endl
- << " sum: " << sum << std::endl;
- std::cout.precision( prec );
- {
- const std::string filenameFlow="free-flow-statistics";
- std::ofstream ofs( filenameFlow+".txt", std::ofstream::out | std::ofstream::trunc);
- const auto prec = ofs.precision();
- ofs << "Velocity statistics (free flow):" << std::endl
- << std::setprecision(std::numeric_limits<double>::digits10 + 1)
- << " min: " << min << std::endl
- << " max: " << max << std::endl
- << " sum: " << sum << std::endl;
- ofs.precision( prec );
- ofs.close();
- }
- }
- {
- const std::string filename = getParam<std::string>("Problem.Name") + "-" + freeFlowProblem->name() + "-interface-pressure";
- writePressuresOnInterfaceToFile<FluxVariables>( filename,
- *freeFlowProblem,
- *freeFlowGridVariables,
- sol );
- }
-
- couplingInterface.finalize();
-
- // 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;
-}
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-reversed.cc b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-mono-flux.cc
similarity index 100%
rename from appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-reversed.cc
rename to appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/main_pm-mono-flux.cc
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
index 5653a5951272867713e99856d9e353790958975e..65209ecbadc427f803f42605e9e2f1da965b75f9 100644
--- a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
+++ b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/monolithicdata.hh
@@ -29,6 +29,31 @@ static const std::array<double, 20> velocity =
2.142947374771422e-12
};
+
+static const std::array<int, 20> velocityFaceIdx =
+ {
+ 2,
+ 6,
+ 10,
+ 14,
+ 18,
+ 22,
+ 26,
+ 30,
+ 34,
+ 38,
+ 42,
+ 46,
+ 50,
+ 54,
+ 58,
+ 62,
+ 66,
+ 70,
+ 74,
+ 78
+};
+
static const std::array<double, 20> pressure = {
9.749674884723169e-10,
9.249954142532818e-10,
diff --git a/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-mono-flux.hh
similarity index 100%
rename from appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-reversed.hh
rename to appl/coupling-ff-pm/iterative-reversed-mono-flux-no-precice/pmproblem-mono-flux.hh