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