Add filter for interior elements in volume coupling

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.

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)

dumux-precice/couplingadapter.cc

@@ -1,7 +1,6 @@
 #include "couplingadapter.hh"
 #include <dumux/common/parameters.hh>
 
-#include <algorithm>
 #include <cassert>
 #include <exception>
 #include <limits>
@@ -34,6 +33,7 @@ void CouplingAdapter::announceConfig(const int rank, const int size)
     precice_ = std::make_unique<precice::Participant>(
         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<precice::Participant>(
         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<double> &positions)
+void CouplingAdapter::setSurfaceMesh(const std::string &meshName,
+                                     std::vector<int> &coupledDumuxIDs,
+                                     std::vector<double> &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::pair<std::string, std::string>,
-                  std::vector<double>>::const_iterator it = dataRead_.begin();
-         it != dataRead_.end(); ++it) {
-        std::cout << it->first.first << " " << it->first.second << " "
-                  << it->second.size() << "\n";
-    }
     std::vector<double> &dataVector = dataRead_[key];
 
     precice_->readData(meshName, dataName, vertexIDs_, relativeReadTime,

dumux-precice/couplingadapter.hh

@@ -1,6 +1,8 @@
 #ifndef PRECICEWRAPPER_HH
 #define PRECICEWRAPPER_HH
 
+#include <algorithm>
+#include <dune/grid/common/partitionset.hh>
 #include <map>
 #include <ostream>
 #include <precice/precice.hpp>
@@ -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<class GridView>
+    void filterInteriorEntities(const GridView &gv,
+                                std::vector<int> &ids,
+                                std::vector<double> &positions);
     /*!
      * @brief Get the number of spatial dimensions
      *
@@ -176,22 +191,44 @@ 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<int> &coupledDumuxIDs,
+                        std::vector<double> &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
      *       [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 setMesh(const std::string &meshName,
-                 const std::vector<double> &positions);
+    template<class GridView>
+    std::vector<int> setVolumeMesh(const std::string &meshName,
+                                   std::vector<int> &coupledDumuxIDs,
+                                   std::vector<double> &positions,
+                                   const GridView &gridView);
     /*!
      * @brief Initializes the coupling
      *
@@ -335,5 +372,68 @@ void CouplingAdapter::initializeCheckpoint(SolutionVector &x)
 {
     states_.emplace_back(std::make_unique<SolverStateOnly<SolutionVector>>(x));
 }
+
+template<class GridView>
+void CouplingAdapter::filterInteriorEntities(const GridView &gv,
+                                             std::vector<int> &ids,
+                                             std::vector<double> &positions)
+{
+    const auto &iset = gv.indexSet();
+    const int dim = positions.size() / ids.size();
+
+    std::vector<int> interior;
+    for (const auto &e : elements(gv)) {
+        if (e.partitionType() == Dune::InteriorEntity)
+            interior.push_back(static_cast<int>(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<class GridView>
+std::vector<int> CouplingAdapter::setVolumeMesh(
+    const std::string &meshName,
+    std::vector<int> &coupledDumuxIDs,
+    std::vector<double> &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

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<double> writeScalarData(numberOfVertices);
     std::vector<double> 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()

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<double> writeScalarData(numberOfVertices);
     std::vector<double> 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()

test/CMakeLists.txt

@@ -1 +1,2 @@
 add_subdirectory(singleCheckpointTest)
+add_subdirectory(filterInteriorTest)

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
+)