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/adapter.hpp"
       6             : #include "issm-precice/confreader.hpp"
       7             : #include "issm-precice/coupler.hpp"
       8             : #include "issm-precice/issm.hpp"
       9             : #include "issm-precice/logging.hpp"
      10             : #include "issm-precice/ranges.hpp"
      11             : #include "issm-precice/timer.hpp"
      12             : #include "classes/classes.h"
      13             : #include "shared/shared.h"
      14             : #include <algorithm>
      15             : #include <iterator>
      16             : #include <ranges>
      17             : #include <span>
      18             : 
      19             : namespace stdv = std::views;
      20             : namespace stdr = std::ranges;
      21             : 
      22             : namespace ipc
      23             : {
      24           4 : void run_coupling(Adapter& adapter)
      25             : {
      26           8 :     IPC_LOG_INFO_0(adapter.get_mpi_comm(), "Starting coupling");
      27             : 
      28           4 :     ipc::ScopedTimer timer(ipc::TimingScope::Coupling);
      29             : 
      30             :     {
      31           4 :         ipc::ScopedTimer init_timer(ipc::TimingScope::Initialize);
      32           4 :         adapter.initialize();
      33           4 :         adapter.initialize_data();
      34           4 :     }
      35             : 
      36          11 :     while (adapter.is_coupling_ongoing())
      37             :     {
      38           7 :         auto dt = adapter.get_time_step();
      39           7 :         adapter.read_data(dt);
      40           7 :         adapter.solve(dt);
      41           7 :         adapter.write_data();
      42           7 :         adapter.advance(dt);
      43             :     }
      44             : 
      45           8 :     IPC_LOG_INFO_0(adapter.get_mpi_comm(), "Finished coupling.");
      46           4 : }
      47             : 
      48           4 : Adapter::Adapter(
      49             :     const Config& config, 
      50             :     std::unique_ptr<IIssm>&& issm, 
      51             :     std::unique_ptr<ICoupler>&& coupler, 
      52           4 :     MPI_Comm comm)
      53          32 :     : m_map_functions({
      54           4 :           {PressureEnum, nullptr},
      55           4 :           {TemperatureEnum, nullptr},
      56           4 :           {VelocityEnum, nullptr},
      57           4 :           {ThicknessEnum, nullptr},
      58           4 :           {MaskOceanLevelsetEnum, nullptr},
      59           8 :           {SealevelEnum, nullptr},
      60             :       })
      61           4 :     , m_vertices()
      62           4 :     , m_variables(config.variables)
      63           4 :     , m_issm(std::move(issm))
      64           4 :     , m_coupler(std::move(coupler))
      65           4 :     , m_comm(comm)
      66             : {
      67           8 :     IPC_LOG_INFO_0(m_comm, "Creating coupling adapter");
      68           4 : }
      69             : 
      70           4 : void Adapter::initialize()
      71             : {
      72           4 :     setup_mesh();
      73           4 :     m_coupler->initialize_mesh();
      74           4 : }
      75             : 
      76           4 : void Adapter::initialize_data()
      77             : {
      78           4 :     write_all_data(/*initialize = */ true);
      79           4 :     m_coupler->initialize_data();
      80           4 : }
      81             : 
      82          11 : bool Adapter::is_coupling_ongoing() const
      83             : {
      84          11 :     return m_coupler->is_coupling_ongoing();
      85             : }
      86             : 
      87           7 : void Adapter::read_data(double dt)
      88             : {
      89           7 :     ScopedTimer timer(TimingScope::Read);
      90           7 :     read_all_data(dt);
      91           7 : }
      92             : 
      93           7 : void Adapter::write_data()
      94             : {
      95           7 :     ScopedTimer timer(TimingScope::Write);
      96           7 :     write_all_data();
      97           7 : }
      98             : 
      99           7 : void Adapter::solve(double dt)
     100             : {
     101           7 :     ScopedTimer timer(TimingScope::Solve);
     102           7 :     m_issm->solve(dt);
     103           7 : }
     104             : 
     105           7 : void Adapter::advance(double dt)
     106             : {
     107           7 :     ScopedTimer timer(TimingScope::Advance);
     108           7 :     m_coupler->advance(dt);
     109           7 : }
     110             : 
     111           7 : void Adapter::read_all_data(double dt)
     112             : {
     113          14 :     IPC_LOG_INFO_0(m_comm, "Reading all data");
     114          25 :     auto read_variables = m_variables | stdv::filter([](auto& v) { return v.direction == DataDirection::Read; });
     115          16 :     for (VariableDefinition& var : read_variables)
     116             :     {
     117           9 :         auto map_func = m_map_functions.find(var.id);
     118           9 :         read_data(var, dt, map_func != m_map_functions.end() ? map_func->second : nullptr);
     119             :     }
     120          14 :     IPC_LOG_INFO_0(m_comm, "Finished reading all data");
     121           7 : }
     122             : 
     123          11 : void Adapter::write_all_data(bool initialize)
     124             : {
     125          22 :     IPC_LOG_INFO_0(m_comm, "Writing all data");
     126          39 :     auto write_variables = m_variables | stdv::filter([](auto& v) { return v.direction == DataDirection::Write; });
     127          25 :     for (VariableDefinition& var : write_variables)
     128             :     {
     129          14 :         if (initialize)
     130             :         {
     131           5 :             m_issm->initialize(var.id);
     132             :         }
     133             : 
     134          14 :         auto map_func = m_map_functions.find(var.id);
     135          14 :         write_data(var, map_func != m_map_functions.end() ? map_func->second : nullptr);
     136             :     }
     137          22 :     IPC_LOG_INFO_0(m_comm, "Finished writing all data");
     138          11 : }
     139             : 
     140             : // set up the mesh: calculate coordinates for preCICE by IDs
     141           4 : void Adapter::setup_mesh()
     142             : {
     143           8 :     IPC_LOG_INFO_0(m_comm, "Setting up Mesh");
     144           4 :     m_vertices = m_issm->get_vertices();
     145             : 
     146             :     // sort the vertices so vertices of elements are continuous (probably already sorted)
     147        1122 :     stdr::sort(m_vertices, {}, [](auto& v) { return v.element_id; });
     148             : 
     149             :     // coordinates of mesh vertices
     150           4 :     auto dim = m_coupler->get_dimensions();
     151             :     auto coordinates = to_vector(
     152           4 :         m_vertices                                                                    //
     153         100 :         | stdv::transform([dim](auto& v) { return v.coordinates | stdv::take(dim); }) //
     154           4 :         | stdv::join);
     155           4 :     m_coupler->setup_coordinates(coordinates);
     156             : 
     157             :     // mesh connectivity
     158             :     // vertices are sorted by element, so a continuous sequence defines the elements
     159           4 :     auto element_vertex_indices = to_vector(stdv::iota(size_t(0), m_vertices.size()));
     160           4 :     m_coupler->setup_simplex_elements(element_vertex_indices);
     161           4 : }
     162             : 
     163           9 : void Adapter::read_data(const VariableDefinition& var, double dt, std::function<void(double*)> map_func)
     164             : {
     165             :     // read
     166           9 :     IPC_LOG_INFO_0(m_comm, "Receiving {} from coupling", var.data_name);
     167           9 :     std::vector<double> values(m_vertices.size());
     168           9 :     m_coupler->read_data(var.data_name, dt, values);
     169             : 
     170             :     // apply map
     171           9 :     if (map_func != nullptr)
     172             :     {
     173           0 :         map_func(values.data());
     174             :     }
     175             : 
     176             :     // set in issm
     177           9 :     if (var.flags & VariableFlags::Constraint)
     178             :     {
     179           0 :         m_issm->set_constraints(var.id, m_vertices, values);
     180             :     }
     181             :     else
     182             :     {
     183           9 :         m_issm->set_input(var.id, m_vertices, values);
     184           9 :         if (var.flags & VariableFlags::Extrude)
     185             :         {
     186           1 :             m_issm->extrude(var.id, m_vertices);
     187             :         }
     188             :     }
     189           9 : }
     190             : 
     191          14 : void Adapter::write_data(const VariableDefinition& var, std::function<void(double*)> map_func)
     192             : {
     193             :     // get data values
     194          14 :     IPC_LOG_INFO_0(m_comm, "Get ISSM input values for {}", var.data_name);
     195          14 :     auto values = std::vector<double>(m_vertices.size());
     196             : 
     197          14 :     auto input = var.id;
     198          14 :     if (var.flags & VariableFlags::DepthAverage)
     199             :     {
     200           2 :         input = m_issm->depth_average(input, m_vertices);
     201             :     }
     202          14 :     m_issm->get_input(input, m_vertices, values);
     203             : 
     204             :     // apply map
     205          14 :     if (map_func != nullptr)
     206             :     {
     207           0 :         map_func(values.data());
     208             :     }
     209             : 
     210             :     // write
     211          14 :     IPC_LOG_INFO_0(m_comm, "Sending {} to coupling", var.data_name);
     212          14 :     m_coupler->write_data(var.data_name, values);
     213          14 : }
     214             : 
     215           7 : double Adapter::get_time_step() const
     216             : {
     217           7 :     return std::min(m_coupler->get_time_step(), m_issm->get_max_time_step());
     218             : }
     219             : 
     220             : } // namespace ipc