Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Daniel Abele <daniel.abele@dlr.de>
       2             : //
       3             : // SPDX-License-Identifier: BSD-3-Clause
       4             : 
       5             : #include "issm-precice/coupler.hpp"
       6             : #include "issm-precice/compiler.hpp"
       7             : #include "issm-precice/logging.hpp"
       8             : #include "issm-precice/ranges.hpp"
       9             : #include <cassert>
      10             : #include <map>
      11             : #include <ranges>
      12             : #include <string>
      13             : #include <vector>
      14             : 
      15             : namespace stdr = std::ranges;
      16             : namespace stdv = std::views;
      17             : 
      18             : namespace ipc
      19             : {
      20           8 : ICoupler::~ICoupler()
      21             : {
      22           8 : }
      23             : 
      24           4 : PreciceCoupler::PreciceCoupler(
      25             :     const std::vector<VariableDefinition>& variables,
      26             :     const std::string& config_file,
      27             :     const std::string& solver_name,
      28             :     const std::string& mesh_name,
      29           4 :     MPI_Comm comm)
      30           4 :     : m_initialized_mesh(false)
      31           4 :     , m_initialized_data(false)
      32           4 :     , m_finalized(false)
      33           4 :     , m_vertex_ids()
      34           4 :     , m_precice(
      35             :           solver_name,
      36             :           config_file,
      37           4 :           [comm]
      38             :           {
      39             :               int r;
      40           4 :               MPI_Comm_rank(comm, &r);
      41           4 :               return r;
      42             :           }(),
      43           4 :           [comm]
      44             :           {
      45             :               int sz;
      46           4 :               MPI_Comm_size(comm, &sz);
      47           4 :               return sz;
      48             :           }(),
      49             :           &comm)
      50           4 :     , m_mesh_name(mesh_name)
      51           4 :     , m_comm(comm)
      52             : {
      53             :     // check valid dimensions
      54           4 :     maybe_unused(variables);
      55           4 :     auto dimensions = m_precice.getMeshDimensions(m_mesh_name);
      56           4 :     if (dimensions < 2 || dimensions > 3)
      57             :     {
      58           0 :         throw std::domain_error("Unsupported number of dimensions.");
      59             :     }
      60           4 : }
      61             : 
      62           4 : PreciceCoupler::~PreciceCoupler()
      63             : {
      64           4 :     if (!m_finalized)
      65             :     {
      66           4 :         finalize();
      67             :     }
      68           4 : }
      69             : 
      70           4 : void PreciceCoupler::finalize()
      71             : {
      72           8 :     IPC_LOG_INFO_0(m_comm, "Finalizing preCICE coupler.");
      73           4 :     assert(!m_finalized);
      74           4 :     m_precice.finalize();
      75           4 :     m_finalized = true;
      76           4 : }
      77             : 
      78          12 : int PreciceCoupler::get_dimensions() const
      79             : {
      80          12 :     return m_precice.getMeshDimensions(m_mesh_name);
      81             : }
      82             : 
      83           4 : void PreciceCoupler::initialize_mesh()
      84             : {
      85             :     //nothing to do since precice 3
      86           4 :     m_initialized_mesh = true;
      87           4 : }
      88             : 
      89           4 : void PreciceCoupler::initialize_data()
      90             : {
      91           4 :     m_initialized_data = true;
      92           4 :     if (!m_precice.requiresInitialData())
      93             :     {
      94             :         // it's very unlikely that no data initialization is required, so we warn here.
      95             :         // also stops the exception thrown by precice when data initialization is enabled but requiresInitialData is
      96             :         // never called.
      97           0 :         IPC_LOG_WARN_0(m_comm, "No data initialization required. This is usually an error, did you forget to enable "
      98             :             "data initialization in your precice config? (<exchange ... initialize=\"1\">");
      99             :     }
     100           4 :     m_precice.initialize();
     101           4 : }
     102             : 
     103           4 : void PreciceCoupler::setup_coordinates(std::span<const double> coordinates)
     104             : {
     105           4 :     assert(coordinates.size() % get_dimensions() == 0);
     106             : 
     107           4 :     auto num_vertices = coordinates.size() / get_dimensions();
     108             : 
     109           4 :     m_vertex_ids.resize(num_vertices);
     110           4 :     m_precice.setMeshVertices(m_mesh_name, coordinates, m_vertex_ids);
     111           4 : }
     112             : 
     113           4 : void PreciceCoupler::setup_simplex_elements(std::span<const size_t> vertex_idxs)
     114             : {
     115           4 :     if (get_dimensions() == 2)
     116             :     {
     117             :         //triangles
     118           4 :         assert(vertex_idxs.size() % 3 == 0);
     119          16 :         for (auto i = size_t(0); i < vertex_idxs.size(); i += 3)
     120             :         {
     121          12 :             m_precice.setMeshTriangle(
     122          12 :                 m_mesh_name,
     123          12 :                 m_vertex_ids[vertex_idxs[i]],
     124          12 :                 m_vertex_ids[vertex_idxs[i + 1]],
     125          12 :                 m_vertex_ids[vertex_idxs[i + 2]]);
     126             :         }
     127             :     }
     128           0 :     else if (get_dimensions() == 3)
     129             :     {
     130             :         //tetrahedrons
     131           0 :         assert(vertex_idxs.size() % 4 == 0);
     132           0 :         for (auto i = size_t(0); i < vertex_idxs.size(); i += 4)
     133             :         {
     134           0 :             m_precice.setMeshTetrahedron(
     135           0 :                 m_mesh_name,
     136           0 :                 m_vertex_ids[vertex_idxs[i]],
     137           0 :                 m_vertex_ids[vertex_idxs[i + 1]],
     138           0 :                 m_vertex_ids[vertex_idxs[i + 2]],
     139           0 :                 m_vertex_ids[vertex_idxs[i + 3]]);
     140             :         }
     141             :     }
     142           4 : }
     143             : 
     144           0 : bool PreciceCoupler::is_coupling_ongoing() const
     145             : {
     146           0 :     return m_precice.isCouplingOngoing();
     147             : }
     148             : 
     149           4 : void PreciceCoupler::read_data(const std::string& data_name, double dt, std::span<double> data) const
     150             : {
     151           4 :     m_precice.readData(m_mesh_name, data_name, m_vertex_ids, dt, data);
     152           4 : }
     153             : 
     154           4 : void PreciceCoupler::write_data(const std::string& data_name, std::span<const double> data)
     155             : {
     156           4 :     m_precice.writeData(m_mesh_name, data_name, m_vertex_ids, data);
     157           4 : }
     158             : 
     159           6 : void PreciceCoupler::advance(double dt)
     160             : {
     161           6 :     IPC_LOG_INFO_0(m_comm, "Advance preCICE by {}", dt);
     162           6 :     m_precice.advance(dt);
     163           6 : }
     164             : 
     165           4 : double PreciceCoupler::get_time_step() const
     166             : {
     167           4 :     return m_precice.getMaxTimeStepSize();
     168             : }
     169             : 
     170             : } // namespace ipc