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/issm.hpp"
6 : #include "issm-precice/logging.hpp"
7 : #include "issm-precice/timer.hpp"
8 : #include "shared/shared.h"
9 : #include "classes/classes.h"
10 : #include "cores/cores.h"
11 : #include <algorithm>
12 : #include <cmath>
13 : #include <iterator>
14 : #include <limits>
15 : #include <ranges>
16 : #include <span>
17 :
18 : namespace stdr = std::ranges;
19 : namespace stdv = std::views;
20 :
21 : namespace ipc
22 : {
23 26091 : const char* issm_input_name(definitions input_id)
24 : {
25 26091 : return EnumToStringx(input_id);
26 : }
27 :
28 0 : definitions issm_input_id(const std::string& input_name)
29 : {
30 0 : return definitions(StringToEnumx(input_name.c_str()));
31 : }
32 :
33 : /**
34 : * Turn a ISSM DataSet into a C++ range.
35 : * Cast the void ptrs stored by the DataSet into the specified type.
36 : * @tparam Type of objects stored by the DataSet.
37 : * @param ds An ISSM DataSet representing a collection of abstract objects.
38 : * @returns Random access view of the range of objects in the DataSet as T*.
39 : */
40 : template<class T, class DS, class = std::enable_if_t<std::is_base_of_v<DataSet, std::decay_t<DS>>>>
41 62671 : auto as_range(DS&& ds)
42 : {
43 1832365 : return ds.objects | std::views::transform([](auto& o) { return static_cast<T*>(o); });
44 : }
45 :
46 34 : FemModelPtr make_fem_model(const std::filesystem::path& root_path, const std::string& model_name, MPI_Comm comm)
47 : {
48 34 : IPC_LOG_INFO_0(comm, "Creating ISSM Model using {} / {}", root_path.string(), model_name);
49 :
50 : //ISSM requires non-const char*, but doesn't modify, so cast is safe.
51 34 : char* args[] = {
52 : const_cast<char*>(""),
53 : const_cast<char*>("TransientSolution"),
54 34 : const_cast<char*>(root_path.c_str()),
55 34 : const_cast<char*>(model_name.c_str())};
56 34 : FemModel* model = new FemModel{4, args, comm};
57 68 : return FemModelPtr(model);
58 : }
59 :
60 34 : Issm::Issm(const std::filesystem::path& root_path, const std::string& model_name, const MeshConfig& mesh, MPI_Comm comm)
61 34 : : m_model(make_fem_model(root_path, model_name, comm))
62 34 : , m_coupling_mesh_config(mesh)
63 34 : , m_initialized_stress_balance(false)
64 : {
65 : //check for compatible mesh type
66 : int domain_type;
67 34 : m_model->parameters->FindParam(&domain_type, DomainTypeEnum);
68 :
69 34 : if (domain_type == Domain2DverticalEnum)
70 : {
71 0 : throw std::runtime_error("Coupling for vertical domains not implemented.");
72 : }
73 :
74 : //time steps are controlled by precice
75 : //store the start time as the current time so it can be used in advance().
76 : //ISSM only sets the current time when simulation begins.
77 : double start_time;
78 34 : m_model->parameters->FindParam(&start_time, TimesteppingStartTimeEnum);
79 34 : m_model->parameters->SetParam(start_time, TimeEnum);
80 :
81 : // ISSM can make any step required by the coupler by making multiple time steps.
82 : // We currently do not support subcycling (i.e. we don't call advance() for every ISSM time step)
83 : // because ISSM writes one ouput step per call to FemModel::Solve, and that could be data to handle.
84 34 : m_max_dt = std::numeric_limits<double>::max();
85 34 : }
86 :
87 : /**
88 : * Deleter that deletes an Element if it was spawned temporarily.
89 : * Elements are spawned e.g. as a the base triangle element of a prism element.
90 : * Otherwise the Element is still required by the Model.
91 : */
92 : struct SpawnedElementDeleter
93 : {
94 32800 : void operator()(::Element* element)
95 : {
96 32800 : if (element->IsSpawnedElement())
97 : {
98 25600 : element->DeleteMaterials();
99 25600 : delete element;
100 : }
101 32800 : }
102 : };
103 :
104 38 : IIssm::~IIssm()
105 : {
106 38 : }
107 :
108 30 : Mesh Issm::get_vertices() const
109 : {
110 30 : return get_vertices(m_coupling_mesh_config);
111 : }
112 :
113 41 : Mesh Issm::get_vertices(const MeshConfig& mesh) const
114 : {
115 82 : return std::visit([this](auto& mesh_type) { return this->get_vertices(mesh_type); }, mesh.type);
116 : }
117 :
118 41 : Mesh Issm::get_vertices(const Mesh2dConfig& mesh2d) const
119 : {
120 41 : auto layer = mesh2d.layer;
121 41 : if (layer != 0 && layer != -1)
122 : {
123 0 : throw std::runtime_error("Only coupling with base (0) and surface (-1) layers implemented.");
124 : }
125 :
126 196000 : auto is_on_layer = [layer](auto&& el)
127 : {
128 109600 : return (layer == 0 && el->IsOnBase()) || (layer == -1 && el->IsOnSurface());
129 41 : };
130 41 : auto layer_elements = as_range<::Element>(*m_model->elements) //
131 82 : | std::views::filter(is_on_layer);
132 :
133 41 : Mesh mesh;
134 41 : const size_t INVALID_IDX = std::numeric_limits<size_t>::max();
135 41 : auto vertex_index_map = std::vector<size_t>(m_model->vertices->Size(), INVALID_IDX);
136 32841 : for (auto&& element : layer_elements)
137 : {
138 : int element_idx;
139 32800 : m_model->elements->GetObjectById(&element_idx, element->Id());
140 :
141 32800 : mesh.elements.push_back(Element{.id = (size_t)element_idx});
142 :
143 : auto layer_element = std::unique_ptr<::Element, SpawnedElementDeleter>(
144 32800 : layer == 0 ? element->SpawnBasalElement() : element->SpawnTopElement());
145 32800 : assert(layer_element->GetNumberOfVertices() == 3 && "Non triangle element.");
146 131200 : for (int i = 0; i < layer_element->GetNumberOfVertices(); ++i)
147 : {
148 98400 : auto vertex = layer_element->vertices[i];
149 : int vertex_idx;
150 98400 : m_model->vertices->GetObjectById(&vertex_idx, vertex->Id());
151 98400 : if (vertex_index_map[vertex_idx] == INVALID_IDX)
152 : {
153 0 : mesh.vertices.push_back(Vertex{
154 18081 : .id = size_t(vertex_idx),
155 18081 : .coordinates = {vertex->x, vertex->y, vertex->z},
156 18081 : .element_id = (size_t)element_idx,
157 : });
158 18081 : IPC_LOG_TRACE("Vertex {} [{}, {}]", vertex->id, vertex->x, vertex->y);
159 18081 : vertex_index_map[vertex_idx] = mesh.vertices.size() - 1;
160 : }
161 98400 : mesh.elements.back().vertices[i] = vertex_index_map[vertex_idx];
162 : }
163 32800 : }
164 :
165 41 : IPC_LOG_TRACE("Use {} vertices of {} elements.", mesh.vertices.size(), mesh.elements.size());
166 :
167 82 : return mesh;
168 41 : }
169 :
170 0 : Mesh Issm::get_vertices(const Mesh3dConfig&) const
171 : {
172 0 : throw std::runtime_error("3D coupling interface not implemented.");
173 : }
174 :
175 1768 : auto get_analysis(definitions input)
176 : {
177 1768 : switch (input)
178 : {
179 1768 : case VxEnum:
180 : case VyEnum:
181 1768 : return StressbalanceAnalysisEnum;
182 0 : default:
183 0 : throw std::runtime_error(fmt::format("Constraint coupling not implemented for {}.", EnumToStringx(input)));
184 : }
185 : }
186 :
187 2646 : auto get_constraint_dof_idx(definitions input)
188 : {
189 2646 : switch (input)
190 : {
191 882 : case VxEnum:
192 882 : return 0;
193 1764 : case VyEnum:
194 1764 : return 1;
195 0 : default:
196 0 : return 0;
197 : }
198 : }
199 :
200 4 : void Issm::set_constraints(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values)
201 : {
202 4 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Setting ISSM constraints for {}", issm_input_name(input));
203 :
204 4 : m_model->SetCurrentConfiguration(get_analysis(input));
205 :
206 : // TODO: this condition might be relaxed if we use coupling meshes based on nodes instead of vertices.
207 : static int allowed_fe_types[] = {P1Enum, P1P1Enum};
208 4 : if (stdr::find(allowed_fe_types, as_range<::Element>(*m_model->elements)[0]->element_type) == std::end(allowed_fe_types))
209 : {
210 0 : throw std::runtime_error("Constraint coupling only implemented for P1 elements (nodes == vertices).");
211 : }
212 :
213 : // get values of constraints of ghost/clone vertices. NaN means no constraint
214 : static_assert(std::numeric_limits<double>::has_quiet_NaN);
215 : auto all_vertex_values =
216 4 : std::vector<double>(m_model->vertices->NumberOfVerticesLocalAll(), std::numeric_limits<double>::quiet_NaN());
217 4 : auto vs = as_range<::Vertex>(*m_model->vertices);
218 4 : assert(all_vertex_values.size() == vs.size());
219 1768 : for (size_t vtx_idx = 0; vtx_idx < vertices.size(); ++vtx_idx)
220 : {
221 1764 : all_vertex_values[vs[vertices[vtx_idx].id]->Lid()] = values[vtx_idx];
222 : }
223 4 : m_model->SyncLocalVectorWithClonesVertices(all_vertex_values.data());
224 :
225 : // set constraints
226 : // !! loop over ALL vertices of the model, vertices in the coupling mesh do not include ghosts/clones !!
227 8824 : for (size_t vtx_idx = 0; vtx_idx < all_vertex_values.size(); ++vtx_idx)
228 : {
229 8820 : auto vertex = as_range<::Vertex>(*m_model->vertices)[vtx_idx];
230 8820 : if (std::isnan(all_vertex_values[vertex->Lid()]))
231 : {
232 7056 : continue; // constraint not set
233 : }
234 :
235 : // find constraint on this vertex
236 : // This can probably be optimised since the constraint should have a similar order as the vertices
237 1764 : auto static_constraints = as_range<::Constraint>(*m_model->constraints) |
238 1568196 : stdv::transform([](auto cons) { return dynamic_cast<SpcStatic*>(cons); }) |
239 782334 : stdv::filter([](auto cons) { return cons != nullptr; });
240 1764 : auto constraint = stdr::find_if(
241 : static_constraints,
242 778806 : [&](auto static_constraint)
243 : {
244 781452 : return static_constraint->GetNodeId() == vertex->Sid() + 1 &&
245 781452 : static_constraint->GetDof() == get_constraint_dof_idx(input);
246 : });
247 1764 : if (constraint == end(static_constraints))
248 : {
249 0 : throw std::runtime_error(fmt::format(
250 : "Constraint not found for vertex {}, Sid {}, Value {}.",
251 : vtx_idx,
252 0 : vertex->Sid(),
253 0 : all_vertex_values[vertex->Lid()]));
254 : }
255 : // set constraint
256 5292 : **constraint = SpcStatic(
257 1764 : (**constraint).Id(),
258 : (**constraint).GetNodeId(),
259 : (**constraint).GetDof(),
260 1764 : all_vertex_values[vertex->Lid()],
261 3528 : get_analysis(input));
262 : }
263 4 : }
264 :
265 16 : void Issm::set_input(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values)
266 : {
267 16 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Setting ISSM input values for {}", issm_input_name(input));
268 :
269 : //set value for each vertex.
270 : //Could probably be made more efficient by setting all values at the same time but
271 : //making this work for many different types of meshes is probably not worth it.
272 7072 : for (size_t i = 0; i < vertices.size(); ++i)
273 : {
274 7056 : auto vertex = as_range<::Vertex>(*m_model->vertices)[vertices[i].id];
275 7056 : auto element = as_range<::Element>(*m_model->elements)[vertices[i].element_id];
276 7056 : IPC_LOG_TRACE("Set value for {} at vertex local {} global {} ghost {} [{} {}] = {}", issm_input_name(input), vertex->lid, vertex->sid, vertex->clone ? 1 : 0, vertex->x, vertex->y, values[i]);
277 7056 : element->SetElementInput(m_model->inputs, 1, &vertex->lid, const_cast<double*>(&values[i]), input);
278 : }
279 16 : }
280 :
281 43 : void Issm::get_input(definitions input, const std::vector<Vertex>& vertices, std::span<double> values) const
282 : {
283 43 : assert(values.size() == vertices.size());
284 :
285 43 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Get ISSM input values for {}", issm_input_name(input));
286 :
287 : //get value for each vertex.
288 : //Could probably be made more efficient by setting all values at the same time but
289 : //making this work for many different types of meshes is probably not worth it.
290 18963 : auto get_value = [&](auto& v)
291 : {
292 18963 : auto vertex = as_range<::Vertex>(*m_model->vertices)[v.id];
293 18963 : auto element = as_range<::Element>(*m_model->elements)[v.element_id];
294 18963 : IPC_LOG_TRACE("Get value of at vertex {} [{} {}]", vertex->lid, vertex->x, vertex->y);
295 : double value;
296 18963 : if (element->GetInput(input))
297 : {
298 18963 : element->GetInputValue(&value, vertex, input);
299 : }
300 : else
301 : {
302 0 : if (input == VelEnum)
303 : {
304 : double vx, vy, vz;
305 0 : element->GetInputValue(&vx, vertex, VxEnum);
306 0 : element->GetInputValue(&vy, vertex, VyEnum);
307 0 : element->GetInputValue(&vz, vertex, VzEnum);
308 0 : value = sqrt(vx * vx + vy * vy + vz * vz);
309 : }
310 : else
311 : {
312 0 : IPC_LOG_ERROR("Input {} not found.", issm_input_name(input));
313 0 : throw std::runtime_error("Input not found.");
314 : }
315 : }
316 18963 : IPC_LOG_TRACE("Got value for {} at vertex local {} global {} ghost {} [{} {}] = {}", issm_input_name(input), vertex->lid, vertex->sid, vertex->clone ? 1 : 0, vertex->x, vertex->y, value);
317 :
318 18963 : return value;
319 43 : };
320 43 : std::ranges::transform(vertices, values.begin(), get_value);
321 43 : }
322 :
323 3 : void Issm::extrude(definitions input, const std::vector<Vertex>& /*vertices*/)
324 : {
325 3 : if (auto p_mesh2d = get_if<Mesh2dConfig>(&m_coupling_mesh_config.type); p_mesh2d)
326 : {
327 3 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Extruding ISSM input {}", issm_input_name(input));
328 3 : m_model->parameters->SetParam(input, InputToExtrudeEnum);
329 3 : if (p_mesh2d->layer == 0)
330 : {
331 2 : extrudefrombase_core(m_model.get());
332 : }
333 1 : else if (p_mesh2d->layer == -1)
334 : {
335 1 : extrudefromtop_core(m_model.get());
336 : }
337 : else
338 : {
339 0 : throw std::runtime_error("Extruding from layer other than base or surface not implemented.");
340 : }
341 : }
342 3 : }
343 :
344 5 : definitions depth_average_input(definitions input)
345 : {
346 5 : switch (input)
347 : {
348 2 : case VxEnum: return VxAverageEnum;
349 2 : case VyEnum: return VyAverageEnum;
350 0 : case WaterfractionDrainageEnum: return WaterfractionDrainageIntegratedEnum;
351 0 : case MaterialsRheologyBEnum: return MaterialsRheologyBbarEnum;
352 0 : case DamageDEnum: return DamageDbarEnum;
353 0 : case DamageDOldEnum: return DamageDbarOldEnum;
354 0 : case MaterialsRheologyEEnum: return MaterialsRheologyEbarEnum;
355 0 : case MaterialsRheologyEsEnum: return MaterialsRheologyEsbarEnum;
356 0 : case MaterialsRheologyEcEnum: return MaterialsRheologyEcbarEnum;
357 1 : default: throw std::domain_error("Depth averaging not implemented for this input.");
358 : }
359 : }
360 :
361 3 : definitions Issm::depth_average(definitions input, const std::vector<Vertex>& /*vertices*/)
362 : {
363 3 : auto average_input = depth_average_input(input);
364 2 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Depth averaging ISSM input {} into {}", issm_input_name(input), issm_input_name(average_input));
365 2 : if (auto p_mesh2d = get_if<Mesh2dConfig>(&m_coupling_mesh_config.type); !p_mesh2d || p_mesh2d->layer != 0)
366 : {
367 0 : IPC_LOG_WARN_0(IssmComm::GetComm(), "Values of depth averaged input {} only valid at base layer.", issm_input_name(average_input));
368 : }
369 :
370 2 : m_model->parameters->SetParam(input, InputToDepthaverageInEnum);
371 2 : m_model->parameters->SetParam(average_input, InputToDepthaverageOutEnum);
372 2 : depthaverage_core(m_model.get());
373 2 : return average_input;
374 : }
375 :
376 3 : void Issm::initialize(definitions input)
377 : {
378 :
379 3 : switch (input)
380 : {
381 0 : case MaskIceLevelsetEnum:
382 : case MaskOceanLevelsetEnum:
383 : case ThicknessEnum:
384 : case BaseEnum:
385 : case SurfaceEnum:
386 : case BedEnum:
387 : //need no initialization, must be initialized in a valid setup
388 0 : break;
389 3 : case VxEnum:
390 : case VyEnum:
391 : case VzEnum:
392 : case VelEnum:
393 : case PressureEnum:
394 3 : if (!m_initialized_stress_balance)
395 : {
396 3 : stressbalance_core(m_model.get());
397 3 : m_model->results->clear(); //don't store results
398 3 : m_initialized_stress_balance = true;
399 : }
400 3 : break;
401 0 : default:
402 0 : IPC_LOG_WARN_0(
403 : IssmComm::GetComm(),
404 : "Initializing coupling data of {}. "
405 : "Make sure this field has a valid initial value as the adapter can not compute it. "
406 : "A procedure to compute the initial value can be added to ipc::Issm::initialize(). ",
407 : issm_input_name(input));
408 0 : break;
409 : }
410 3 : }
411 :
412 19 : double Issm::get_max_time_step() const
413 : {
414 19 : return m_max_dt;
415 : }
416 :
417 19 : void Issm::solve(double dt)
418 : {
419 19 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Advance ISSM by {}", dt);
420 :
421 19 : assert(dt < get_max_time_step() && "Time step too big.");
422 :
423 : double t;
424 19 : m_model->parameters->FindParam(&t, TimeEnum);
425 :
426 19 : if (m_dt_override > 0)
427 : {
428 0 : m_model->parameters->SetParam(m_dt_override, TimesteppingTimeStepEnum);
429 : }
430 19 : if (m_output_frequency_override > 0)
431 : {
432 0 : m_model->parameters->SetParam(m_output_frequency_override, SettingsOutputFrequencyEnum);
433 : //TODO: find a way to disable all output (including at the end of "solve") with output_frequency_override == 0!
434 : }
435 :
436 19 : m_model->parameters->SetParam(t, TimesteppingStartTimeEnum);
437 19 : m_model->parameters->SetParam(t + dt, TimesteppingFinalTimeEnum);
438 :
439 19 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Run ISSM from t = {} to {}", t, t + dt);
440 :
441 19 : m_model->Solve();
442 19 : }
443 :
444 0 : void Issm::set_time_step_override(double dt)
445 : {
446 0 : m_dt_override = dt;
447 0 : }
448 :
449 0 : void Issm::set_output_frequency_override(int f)
450 : {
451 0 : m_output_frequency_override = f;
452 0 : }
453 :
454 : } // namespace ipc
|
