diff --git a/docs/use.md b/docs/use.md index 7bdaf08..1dcff93 100644 --- a/docs/use.md +++ b/docs/use.md @@ -7,6 +7,30 @@ summary: How to use the DuMuX adapter for building your own coupled solver. To understand how the adapter is used, see the tutorial [free-flow-over-porous-media](https://github.com/precice/tutorials/tree/develop/free-flow-over-porous-media). Additionally, the solver [macro-dumux](https://github.com/precice/tutorials/tree/develop/two-scale-heat-conduction/macro-dumux) uses the adapter to couple a heat conduction problem to micro-scale simulations. +### Configuration + The adapter is configured via the group called `[precice-adapter-config]` in the DuMux runtime parameter file with the default name `params.input`. The input files for the dummy simulation under `examples/dummysolver` provided examples of the parameters to be configured. This configuration follows the nomenclature of the [preCICE adapter configuration schema](https://github.com/precice/preeco-orga/tree/main/adapter-config-schema). It does not adhere to the schema completely because the configuration is done via a `.input` file instead of a JSON or YAML file. +### Use the API + +#### Set coupling mesh + +To inform preCICE the coupling mesh, vertices coordinates and `meshName` are to be set with `setSurfaceMesh()` or `setVolumeMesh()`. + +Both functions require: + +- `meshName`, which is the name of the mesh to be set +- `coupledDumuxIDs`, which can be indices of the finite-volume face in a surface coupling case or indices of elements in a volume coupling case, and +- `positions`,which are the spatial coordinates that the previous IDs represent. + +In addition, `setVolumeMesh()` requires `gridView`, which is used to filter the `overlap` and `ghost` cells in a distributed setting. This is meaningful for defining the correct number of elements that are coupled. + +In a distributed surface coupling case, there can also be `overlap` and `ghost` cells on the coupling interface. The filtering of these indices are not supported in the adapter yet, as this also requires information of the geometry. It's required, for user, to only add cell faces belonging to `interior` elements to the `coupledDumuxIDs` before calling `setSurfaceMesh()`. + +#### Exchange data + +In a parallel setting, the values on the `overlap` and `ghost` elements need to be updated manually via subdomain communications, e.g. `communicate()` provided by DUNE-Grid, after the values of the coupling variables have been set by data from another solver. This is needed since only `interior` elements or cell faces attached to `interior` cells are included in coupled mesh. + +### Additional build guideline + To use the adapter in a separate DUNE module, call `dune_enable_all_packages()` in the root `CMakeLists.txt` of the application module. If `libdumux-precice` is built as a static library, preCICE needs to be explicitly discovered with `find_package` as done in the root `CMakeLists.txt` of the adapter. To build the adapter library as a dynamic library, use the CMake option `-DBUILD_SHARED_LIBS=ON` to build `dumux-precice` and upstream modules. diff --git a/dumux-precice/CMakeLists.txt b/dumux-precice/CMakeLists.txt index 559b84f..4c71e46 100644 --- a/dumux-precice/CMakeLists.txt +++ b/dumux-precice/CMakeLists.txt @@ -3,4 +3,4 @@ install(FILES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux-precice) dune_library_add_sources(dumux-precice SOURCES couplingadapter.cc dumuxpreciceindexmapper.cc solverstate.cc) -target_link_libraries(dumux-precice PRIVATE precice::precice) +target_link_libraries(dumux-precice PUBLIC precice::precice) diff --git a/dumux-precice/couplingadapter.cc b/dumux-precice/couplingadapter.cc index 643a3bb..804eaae 100644 --- a/dumux-precice/couplingadapter.cc +++ b/dumux-precice/couplingadapter.cc @@ -1,7 +1,6 @@ #include "couplingadapter.hh" #include -#include #include #include #include @@ -34,6 +33,7 @@ void CouplingAdapter::announceConfig(const int rank, const int size) precice_ = std::make_unique( participantName_, preciceConfigName_, rank, size); wasCreated_ = true; + inParallel_ = size > 1; int interfaceIndex = 0; do { @@ -101,6 +101,7 @@ void CouplingAdapter::announceSolver(const std::string &name, precice_ = std::make_unique( name, configurationFileName, rank, size); wasCreated_ = true; + inParallel_ = size > 1; } int CouplingAdapter::getMeshDimensions(const std::string &meshName) const @@ -109,13 +110,16 @@ int CouplingAdapter::getMeshDimensions(const std::string &meshName) const return precice_->getMeshDimensions(meshName); } -void CouplingAdapter::setMesh(const std::string &meshName, - const std::vector &positions) +void CouplingAdapter::setSurfaceMesh(const std::string &meshName, + std::vector &coupledDumuxIDs, + std::vector &positions) { assert(wasCreated_); - vertexIDs_.resize(positions.size() / getMeshDimensions(meshName)); + + vertexIDs_.resize(coupledDumuxIDs.size()); precice_->setMeshVertices(meshName, positions, vertexIDs_); meshWasCreated_ = true; + createIndexMapping(coupledDumuxIDs); // compute size of data vectors for coupling data on this mesh auto dataToReadOnMesh = getReadDataNamesOnMesh(meshName); @@ -313,12 +317,6 @@ void CouplingAdapter::readQuantityFromOtherSolver(const std::string &meshName, double relativeReadTime) { auto key = std::make_pair(meshName, dataName); - for (std::map, - std::vector>::const_iterator it = dataRead_.begin(); - it != dataRead_.end(); ++it) { - std::cout << it->first.first << " " << it->first.second << " " - << it->second.size() << "\n"; - } std::vector &dataVector = dataRead_[key]; precice_->readData(meshName, dataName, vertexIDs_, relativeReadTime, diff --git a/dumux-precice/couplingadapter.hh b/dumux-precice/couplingadapter.hh index a9d930d..993fc5d 100644 --- a/dumux-precice/couplingadapter.hh +++ b/dumux-precice/couplingadapter.hh @@ -1,6 +1,8 @@ #ifndef PRECICEWRAPPER_HH #define PRECICEWRAPPER_HH +#include +#include #include #include #include @@ -55,6 +57,8 @@ private: std::string preciceConfigName_; //! Participant or solver name std::string participantName_; + //! If field needs to be synchronized after reading + bool inParallel_ = false; //! Constructor CouplingAdapter(); /*! @@ -100,6 +104,17 @@ public: const std::string &configurationFileName, const int rank, const int size); + /*! + * @brief Filter onon-interior element indices and positions out. + * + * @param[in] gv The grid view to find interior cells. + * @param[in] ids The dumuxIDs of the vertices. + * @param[in] positions The coordinates of the vertices. + */ + template + void filterInteriorEntities(const GridView &gv, + std::vector &ids, + std::vector &positions); /*! * @brief Get the number of spatial dimensions * @@ -176,13 +191,32 @@ public: * @return false No further action is needed. */ bool requiresToWriteInitialData(); - /*! - * @brief Adds mesh for coupling of solvers. With the mesh size, the data maps inside the adapter initialize the relevant data vector to size of meshSize*dataDimension. + * @brief Adds surface-coupled mesh for coupling of solvers. With the mesh size, the data maps inside the adapter initialize the relevant data vector to size of meshSize*dataDimension. * * @param[in] meshName The name of the mesh to add the vertices to. - * @param[in] positions A span to the coordinates of the vertices. + * @param[in] coupledDumuxIDs The dumux IDs of the elements to be coupled. + * @param[in] positions The coordinates of the vertices. + * @return dumux IDs of the actually coupled IDs + * \note The coordinates need to be stored consecutively + * according to their spatial coordinates as.\n + * Example 2D:\n + * [x_1, y_1, x_2, y_2,...x_numPoints, y_numPoints]\n + * Example 3D:\n + * [x_1, y_1, z_1, x_2, y_2, z_2,...x_numPoints, y_numPoints, z_numPoints] + */ + void setSurfaceMesh(const std::string &meshName, + std::vector &coupledDumuxIDs, + std::vector &positions); + + /*! + * @brief Adds volume-coupled mesh for coupling of solvers. With the mesh size, the data maps inside the adapter initialize the relevant data vector to size of meshSize*dataDimension. * + * @param[in] meshName The name of the mesh to add the vertices to. + * @param[in] coupledDumuxIDs The dumux IDs of the elements to be coupled. + * @param[in] positions The coordinates of the vertices. + * @param[in] gridView The gird view used to check repeating cells in parallel simulation + * @return dumux IDs of the actually coupled IDs after filtering out repeating vertices in a distributed case * \note The coordinates need to be stored consecutively * according to their spatial coordinates as.\n * Example 2D:\n @@ -190,8 +224,11 @@ public: * Example 3D:\n * [x_1, y_1, z_1, x_2, y_2, z_2,...x_numPoints, y_numPoints, z_numPoints] */ - void setMesh(const std::string &meshName, - const std::vector &positions); + template + std::vector setVolumeMesh(const std::string &meshName, + std::vector &coupledDumuxIDs, + std::vector &positions, + const GridView &gridView); /*! * @brief Initializes the coupling * @@ -335,5 +372,68 @@ void CouplingAdapter::initializeCheckpoint(SolutionVector &x) { states_.emplace_back(std::make_unique>(x)); } + +template +void CouplingAdapter::filterInteriorEntities(const GridView &gv, + std::vector &ids, + std::vector &positions) +{ + const auto &iset = gv.indexSet(); + const int dim = positions.size() / ids.size(); + + std::vector interior; + for (const auto &e : elements(gv)) { + if (e.partitionType() == Dune::InteriorEntity) + interior.push_back(static_cast(iset.index(e))); + } + + for (std::size_t i = 0; i < ids.size();) { + if (std::find(interior.begin(), interior.end(), ids[i]) == + interior.end()) { + ids.erase(ids.begin() + i); + positions.erase(positions.begin() + i * dim, + positions.begin() + (i + 1) * dim); + } else { + ++i; + } + } +} + +template +std::vector CouplingAdapter::setVolumeMesh( + const std::string &meshName, + std::vector &coupledDumuxIDs, + std::vector &positions, + const GridView &gridView) +{ + assert(wasCreated_); + + if (inParallel_) { + filterInteriorEntities(gridView, coupledDumuxIDs, positions); + } + + vertexIDs_.resize(coupledDumuxIDs.size()); + precice_->setMeshVertices(meshName, positions, vertexIDs_); + meshWasCreated_ = true; + createIndexMapping(coupledDumuxIDs); + + // compute size of data vectors for coupling data on this mesh + auto dataToReadOnMesh = getReadDataNamesOnMesh(meshName); + auto dataToWriteOnMesh = getWriteDataNamesOnMesh(meshName); + + for (auto dataName : dataToReadOnMesh) { + int dataDimension = precice_->getDataDimensions(meshName, dataName); + dataRead_[std::make_pair(meshName, dataName)].resize(vertexIDs_.size() * + dataDimension); + } + + for (auto dataName : dataToWriteOnMesh) { + int dataDimension = precice_->getDataDimensions(meshName, dataName); + dataWrite_[std::make_pair(meshName, dataName)].resize( + vertexIDs_.size() * dataDimension); + } + + return coupledDumuxIDs; +} } // namespace Dumux::Precice #endif diff --git a/examples/dummysolver/participantOne.cc b/examples/dummysolver/participantOne.cc index a7c53ac..70562f4 100644 --- a/examples/dummysolver/participantOne.cc +++ b/examples/dummysolver/participantOne.cc @@ -36,10 +36,10 @@ int main(int argc, char **argv) couplingParticipant.announceConfig(mpiHelper.rank(), mpiHelper.size()); const std::string meshName = couplingParticipant.getMeshNames()[0]; - const int dimensions = couplingParticipant.getMeshDimensions(meshName); - assert(dimensions == 3); + assert(couplingParticipant.getMeshDimensions(meshName) == 3); const int numberOfVertices = 3; + static constexpr int dimensions = 3; std::vector writeScalarData(numberOfVertices); std::vector readScalarData(numberOfVertices); @@ -60,11 +60,7 @@ int main(int argc, char **argv) std::cout << "DUMMY (" << mpiHelper.rank() << "): Initialize preCICE and set mesh\n"; - couplingParticipant.setMesh(meshName, vertices); - - // Create index mapping between DuMuX's index numbering and preCICE's numbering - std::cout << "DUMMY (" << mpiHelper.rank() << "): Create index mapping\n"; - couplingParticipant.createIndexMapping(dumuxVertexIDs); + couplingParticipant.setSurfaceMesh(meshName, dumuxVertexIDs, vertices); if (couplingParticipant.requiresToWriteInitialData()) { std::cout << "DUMMY (" << mpiHelper.rank() diff --git a/examples/dummysolver/participantTwo.cc b/examples/dummysolver/participantTwo.cc index b0f3e13..ce97db8 100644 --- a/examples/dummysolver/participantTwo.cc +++ b/examples/dummysolver/participantTwo.cc @@ -37,9 +37,6 @@ int main(int argc, char **argv) couplingParticipant.announceConfig(mpiHelper.rank(), mpiHelper.size()); const std::string meshName = couplingParticipant.getMeshNames()[0]; - const int dimensions = couplingParticipant.getMeshDimensions(meshName); - assert(dimensions == 3); - const std::string scalarDataWriteName = couplingParticipant.getWriteDataNamesOnMesh(meshName)[0]; const std::string scalarDataReadName = @@ -50,6 +47,7 @@ int main(int argc, char **argv) couplingParticipant.getReadDataNamesOnMesh(meshName)[1]; const int numberOfVertices = 3; + static constexpr int dimensions = 3; std::vector writeScalarData(numberOfVertices); std::vector readScalarData(numberOfVertices); @@ -70,11 +68,7 @@ int main(int argc, char **argv) std::cout << "DUMMY (" << mpiHelper.rank() << "): Initialize preCICE and set mesh\n"; - couplingParticipant.setMesh(meshName, vertices); - - // Create index mapping between DuMuX's index numbering and preCICE's numbering - std::cout << "DUMMY (" << mpiHelper.rank() << "): Create index mapping\n"; - couplingParticipant.createIndexMapping(dumuxVertexIDs); + couplingParticipant.setSurfaceMesh(meshName, dumuxVertexIDs, vertices); if (couplingParticipant.requiresToWriteInitialData()) { std::cout << "DUMMY (" << mpiHelper.rank() diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ec4a070..fc01f24 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(singleCheckpointTest) +add_subdirectory(filterInteriorTest) diff --git a/test/filterInteriorTest/CMakeLists.txt b/test/filterInteriorTest/CMakeLists.txt new file mode 100644 index 0000000..803d1cc --- /dev/null +++ b/test/filterInteriorTest/CMakeLists.txt @@ -0,0 +1,8 @@ +add_executable(dumuxprecice_filterInteriorTest EXCLUDE_FROM_ALL filterInteriorTest.cc) + +dumux_add_test( + NAME dumuxprecice_filterInteriorTest + TARGET dumuxprecice_filterInteriorTest + LABELS unit grid + MPI 2 +) diff --git a/test/filterInteriorTest/filterInteriorTest.cc b/test/filterInteriorTest/filterInteriorTest.cc new file mode 100644 index 0000000..8874a8e --- /dev/null +++ b/test/filterInteriorTest/filterInteriorTest.cc @@ -0,0 +1,87 @@ +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include "dumux-precice/couplingadapter.hh" + +int main(int argc, char **argv) +{ + auto &mpi = Dune::MPIHelper::instance(argc, argv); + Dune::TestSuite t; + + static constexpr int dim = 2; + + // Build a 8*8 grid; + Dune::FieldVector L(1.0); + std::array N{8, 8}; + std::bitset periodic(false); + int overlap = 1; + + Dune::YaspGrid grid(L, N, periodic, overlap); + auto gv = grid.leafGridView(); + const auto &iset = gv.indexSet(); + + std::vector ids; + std::vector pos; + for (const auto &e : elements(gv)) { + ids.push_back(static_cast(iset.index(e))); + const auto c = e.geometry().center(); + pos.push_back(c[0]); + pos.push_back(c[1]); + } + + // Reference: local interior element indices + std::vector refInterior; + for (const auto &e : elements(gv)) { + if (e.partitionType() == Dune::InteriorEntity) + refInterior.push_back(static_cast(iset.index(e))); + } + + std::cout << "\n=== Rank " << mpi.rank() << " ===\n"; + std::cout << "Total local elements: " << gv.size(0) << "\n"; + std::cout << "Initial ids.size(): " << ids.size() << "\n"; + std::cout << "Initial pos.size(): " << pos.size() << "\n"; + std::cout << "Reference interior count: " << refInterior.size() << "\n"; + + std::cout << "Initial local IDs on rank " << mpi.rank() << ":" << "\n"; + for (auto id : ids) + std::cout << id << " "; + std::cout << "\n"; + + auto &couplingParticipant = Dumux::Precice::CouplingAdapter::getInstance(); + couplingParticipant.filterInteriorEntities(gv, ids, pos); + + std::cout << "After filtering ids.size(): " << ids.size() << " on rank " + << mpi.rank() << "\n"; + std::cout << "After filtering pos.size(): " << pos.size() << " on rank " + << mpi.rank() << "\n"; + + std::cout << "Local IDs after filtering on rank " << mpi.rank() << ":" + << "\n"; + for (auto id : ids) + std::cout << id << " "; + std::cout << "\n"; + + t.check(pos.size() == ids.size() * dim) + << "pos size inconsistent after filtering"; + + t.check(ids.size() == refInterior.size()) << "filtered id count mismatch"; + + for (auto id : ids) { + t.check(std::find(refInterior.begin(), refInterior.end(), id) != + refInterior.end()) + << "filtered id not interior: " << id; + } + + std::cout << "Rank " << mpi.rank() << " test finished.\n"; + + return t.exit(); +} diff --git a/test/singleCheckpointTest/main.cc b/test/singleCheckpointTest/main.cc index 4f0c2a3..52a8264 100644 --- a/test/singleCheckpointTest/main.cc +++ b/test/singleCheckpointTest/main.cc @@ -13,6 +13,8 @@ #include #include +#include +#include #include "dumux-precice/couplingadapter.hh" @@ -72,10 +74,10 @@ int main(int argc, char **argv) const std::string meshName = (couplingParticipant.getSolverName() == "SolverOne") ? "MeshOne" : "MeshTwo"; - const int dimensions = couplingParticipant.getMeshDimensions(meshName); - assert(dimensions == 3); + assert(couplingParticipant.getMeshDimensions(meshName) == 3); const int numberOfVertices = 3; + static constexpr int dimensions = 3; std::vector writeScalarData(numberOfVertices); std::vector readScalarData(numberOfVertices); @@ -94,11 +96,7 @@ int main(int argc, char **argv) std::cout << "DUMMY (" << mpiHelper.rank() << "): Initialize preCICE and set mesh\n"; - couplingParticipant.setMesh(meshName, vertices); - - // Create index mapping between DuMuX's index numbering and preCICE's numbering - std::cout << "DUMMY (" << mpiHelper.rank() << "): Create index mapping\n"; - couplingParticipant.createIndexMapping(dumuxVertexIDs); + couplingParticipant.setSurfaceMesh(meshName, dumuxVertexIDs, vertices); if (couplingParticipant.requiresToWriteInitialData()) { std::cout << "DUMMY (" << mpiHelper.rank()