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           0 : bool PreciceCoupler::is_mesh_connectivity_required() const
     104             : {
     105           0 :     return m_precice.requiresMeshConnectivityFor(m_mesh_name);
     106             : }
     107             : 
     108           4 : void PreciceCoupler::setup_coordinates(std::span<const double> coordinates)
     109             : {
     110           4 :     assert(coordinates.size() % get_dimensions() == 0);
     111             : 
     112           4 :     auto num_vertices = coordinates.size() / get_dimensions();
     113             : 
     114           4 :     m_vertex_ids.resize(num_vertices);
     115           4 :     m_precice.setMeshVertices(m_mesh_name, coordinates, m_vertex_ids);
     116           4 : }
     117             : 
     118           4 : void PreciceCoupler::setup_simplex_elements(std::span<const size_t> vertex_idxs)
     119             : {
     120           4 :     if (get_dimensions() == 2)
     121             :     {
     122             :         //triangles
     123           4 :         assert(vertex_idxs.size() % 3 == 0);
     124          16 :         for (auto i = size_t(0); i < vertex_idxs.size(); i += 3)
     125             :         {
     126          12 :             m_precice.setMeshTriangle(
     127          12 :                 m_mesh_name,
     128          12 :                 m_vertex_ids[vertex_idxs[i]],
     129          12 :                 m_vertex_ids[vertex_idxs[i + 1]],
     130          12 :                 m_vertex_ids[vertex_idxs[i + 2]]);
     131             :         }
     132             :     }
     133           0 :     else if (get_dimensions() == 3)
     134             :     {
     135             :         //tetrahedrons
     136           0 :         assert(vertex_idxs.size() % 4 == 0);
     137           0 :         for (auto i = size_t(0); i < vertex_idxs.size(); i += 4)
     138             :         {
     139           0 :             m_precice.setMeshTetrahedron(
     140           0 :                 m_mesh_name,
     141           0 :                 m_vertex_ids[vertex_idxs[i]],
     142           0 :                 m_vertex_ids[vertex_idxs[i + 1]],
     143           0 :                 m_vertex_ids[vertex_idxs[i + 2]],
     144           0 :                 m_vertex_ids[vertex_idxs[i + 3]]);
     145             :         }
     146             :     }
     147           4 : }
     148             : 
     149           0 : bool PreciceCoupler::is_coupling_ongoing() const
     150             : {
     151           0 :     return m_precice.isCouplingOngoing();
     152             : }
     153             : 
     154           4 : void PreciceCoupler::read_data(const std::string& data_name, double dt, std::span<double> data) const
     155             : {
     156           4 :     m_precice.readData(m_mesh_name, data_name, m_vertex_ids, dt, data);
     157           4 : }
     158             : 
     159           4 : void PreciceCoupler::write_data(const std::string& data_name, std::span<const double> data)
     160             : {
     161           4 :     m_precice.writeData(m_mesh_name, data_name, m_vertex_ids, data);
     162           4 : }
     163             : 
     164           6 : void PreciceCoupler::advance(double dt)
     165             : {
     166           6 :     IPC_LOG_INFO_0(m_comm, "Advance preCICE by {}", dt);
     167           6 :     m_precice.advance(dt);
     168           6 : }
     169             : 
     170           4 : double PreciceCoupler::get_time_step() const
     171             : {
     172           4 :     return m_precice.getMaxTimeStepSize();
     173             : }
     174             : 
     175             : } // namespace ipc