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/confreader.hpp"
7 : #include "issm-precice/logging.hpp"
8 : #include "issm-precice/math.hpp"
9 : #include "issm-precice/timer.hpp"
10 : #include "shared/shared.h"
11 : #include "classes/classes.h"
12 : #include "cores/cores.h"
13 : #include <algorithm>
14 : #include <cmath>
15 : #include <iterator>
16 : #include <limits>
17 : #include <ranges>
18 : #include <span>
19 : #include <type_traits>
20 :
21 : namespace stdr = std::ranges;
22 : namespace stdv = std::views;
23 :
24 : namespace ipc
25 : {
26 45923 : const char* issm_input_name(definitions input_id)
27 : {
28 45923 : return EnumToStringx(input_id);
29 : }
30 :
31 0 : definitions issm_input_id(const std::string& input_name)
32 : {
33 0 : return definitions(StringToEnumx(input_name.c_str()));
34 : }
35 :
36 : namespace
37 : {
38 :
39 : /**
40 : * Turn a ISSM DataSet into a C++ range.
41 : * Cast the void ptrs stored by the DataSet into the specified type.
42 : * @tparam Type of objects stored by the DataSet.
43 : * @param ds An ISSM DataSet representing a collection of abstract objects.
44 : * @returns Random access view of the range of objects in the DataSet as T*.
45 : */
46 : template<class T, class DS, class = std::enable_if_t<std::is_base_of_v<DataSet, std::decay_t<DS>>>>
47 108977 : auto as_range(DS&& ds)
48 : {
49 2353513 : return ds.objects | std::views::transform([](auto& o) { return static_cast<T*>(o); });
50 : }
51 :
52 : /**
53 : * Get an Issm parameter.
54 : */
55 : template<class T>
56 343 : T get_parameter(const ::FemModel& model, ::definitions par)
57 : {
58 : T value;
59 343 : const_cast<::FemModel&>(model).parameters->FindParam(&value, par);
60 343 : return value;
61 : }
62 :
63 : /**
64 : * Get an Issm enum parameter.
65 : */
66 : template<>
67 318 : ::definitions get_parameter<::definitions>(const ::FemModel& model, ::definitions par)
68 : {
69 : int value;
70 318 : const_cast<::FemModel&>(model).parameters->FindParam(&value, par);
71 318 : return (::definitions)value;
72 : }
73 :
74 : /**
75 : * Set an Issm parameter.
76 : */
77 : template<class T>
78 490 : void set_parameter(::FemModel& model, ::definitions par, T value)
79 : {
80 490 : model.parameters->SetParam(value, par);
81 490 : }
82 :
83 : } // namespace
84 :
85 118 : FemModelPtr make_fem_model(const std::filesystem::path& root_path, const std::string& model_name, MPI_Comm comm)
86 : {
87 118 : IPC_LOG_INFO_0(comm, "Creating ISSM Model using {} / {}", root_path.string(), model_name);
88 :
89 : //ISSM requires non-const char*, but doesn't modify, so cast is safe.
90 118 : char* args[] = {
91 : const_cast<char*>(""),
92 : const_cast<char*>("TransientSolution"),
93 118 : const_cast<char*>(root_path.c_str()),
94 118 : const_cast<char*>(model_name.c_str())};
95 118 : FemModel* model = new FemModel{4, args, comm};
96 236 : return FemModelPtr(model);
97 : }
98 :
99 : /**
100 : * Deleter that deletes an Element if it was spawned temporarily.
101 : * Elements are spawned e.g. as a the base triangle element of a prism element.
102 : * Otherwise the Element is still required by the Model and not deleted.
103 : */
104 : struct SpawnedElementDeleter
105 : {
106 65600 : void operator()(::Element* element)
107 : {
108 65600 : if (element->IsSpawnedElement())
109 : {
110 42400 : element->DeleteMaterials();
111 42400 : delete element;
112 : }
113 65600 : }
114 : };
115 :
116 133 : std::pair<Mesh, Mesh> make_mesh(const Mesh2dConfig& mesh2d, const ::FemModel& model)
117 : {
118 133 : auto layer = mesh2d.layer;
119 133 : if (layer != 0 && layer != -1)
120 : {
121 0 : throw std::runtime_error("Only coupling with base (0) and surface (-1) layers implemented.");
122 : }
123 :
124 192800 : auto is_on_layer = [layer](auto&& el)
125 : {
126 192800 : return (layer == 0 && el->IsOnBase()) || (layer == -1 && el->IsOnSurface());
127 133 : };
128 133 : auto layer_elements = as_range<::Element>(*model.elements) //
129 266 : | stdv::filter(is_on_layer);
130 :
131 133 : Mesh connected_mesh;
132 133 : Mesh unique_mesh;
133 133 : const size_t INVALID_IDX = std::numeric_limits<size_t>::max();
134 133 : auto vertex_index_map = std::vector<size_t>(model.vertices->Size(), INVALID_IDX);
135 65733 : for (auto&& element : layer_elements)
136 : {
137 : int element_idx;
138 65600 : model.elements->GetObjectById(&element_idx, element->Id());
139 :
140 65600 : connected_mesh.elements.push_back(Element{.id = (size_t)element_idx});
141 :
142 : auto layer_element = std::unique_ptr<::Element, SpawnedElementDeleter>(
143 65600 : layer == 0 ? element->SpawnBasalElement() : element->SpawnTopElement());
144 65600 : assert(layer_element->GetNumberOfVertices() == 3 && "Non triangle element.");
145 262400 : for (int i = 0; i < layer_element->GetNumberOfVertices(); ++i)
146 : {
147 196800 : auto issm_vertex = layer_element->vertices[i];
148 : int vertex_idx;
149 196800 : model.vertices->GetObjectById(&vertex_idx, issm_vertex->Id());
150 196800 : if (vertex_index_map[vertex_idx] == INVALID_IDX)
151 : {
152 : auto vertex = Vertex{
153 37318 : .id = size_t(vertex_idx),
154 37318 : .coordinates = {issm_vertex->x, issm_vertex->y, issm_vertex->z},
155 37318 : .element_id = (size_t)element_idx
156 37318 : };
157 37318 : connected_mesh.vertices.push_back(vertex);
158 37318 : if (!issm_vertex->clone)
159 : {
160 36162 : unique_mesh.vertices.push_back(vertex);
161 : }
162 37318 : IPC_LOG_TRACE("Vertex {} [{}, {}]", issm_vertex->id, issm_vertex->x, issm_vertex->y);
163 37318 : vertex_index_map[vertex_idx] = connected_mesh.vertices.size() - 1;
164 : }
165 196800 : connected_mesh.elements.back().vertices[i] = vertex_index_map[vertex_idx];
166 : }
167 65600 : }
168 :
169 133 : IPC_LOG_TRACE("Use {} vertices of {} elements.", connected_mesh.vertices.size(), connected_mesh.elements.size());
170 :
171 266 : return {connected_mesh, unique_mesh};
172 133 : }
173 :
174 0 : std::pair<Mesh, Mesh> make_mesh(const Mesh3dConfig& /* config */, const ::FemModel& /* model */)
175 : {
176 0 : throw std::runtime_error("3D Mesh not implemented.");
177 : }
178 :
179 118 : Issm::Issm(const std::filesystem::path& root_path, const std::string& model_name, const MeshConfig& mesh, MPI_Comm comm)
180 118 : : m_model(make_fem_model(root_path, model_name, comm))
181 118 : , m_coupling_mesh_config(mesh)
182 236 : , m_initialized_stress_balance(false)
183 : {
184 : //check for compatible mesh type
185 118 : auto domain_type = get_parameter<::definitions>(*m_model, DomainTypeEnum);
186 118 : if (domain_type == Domain2DverticalEnum)
187 : {
188 0 : throw std::runtime_error("Coupling for vertical domains not implemented.");
189 : }
190 :
191 : //time steps are controlled by precice
192 : //store the start time as the current time so it can be used in advance().
193 : //ISSM only sets the current time when simulation begins.
194 118 : auto start_time = get_parameter<double>(*m_model, TimesteppingStartTimeEnum);
195 118 : set_parameter(*m_model, TimeEnum, start_time);
196 :
197 : #ifdef ISSM_PRECICE_FEATURE_CONFIGURE_OUTPUT
198 : // don't save results at "final" time step
199 : // when coupled, there is a "final" time step at the end of every coupling window
200 118 : set_parameter(*m_model, SaveFinalResultsEnum, false);
201 : #endif
202 :
203 118 : if (get_parameter<::definitions>(*m_model, TimesteppingTypeEnum) == FixedTimesteppingEnum)
204 : {
205 : // save the time step set in the model setup;
206 : // we need it during solve and it may be overridden during simulation
207 117 : m_initial_dt = get_parameter<double>(*m_model, TimesteppingTimeStepEnum);
208 : }
209 : else
210 : {
211 : // time step not available if adaptive time stepping
212 1 : m_initial_dt = std::numeric_limits<double>::signaling_NaN();
213 : }
214 :
215 : //initialize meshes
216 118 : std::tie(m_connected_mesh, m_unique_mesh) =
217 354 : std::visit([&](auto&& mesh_type) { return make_mesh(mesh_type, *m_model); }, m_coupling_mesh_config.type);
218 118 : }
219 :
220 122 : IIssm::~IIssm()
221 : {
222 122 : }
223 :
224 102 : Mesh Issm::get_vertices(bool with_connectivity) const
225 : {
226 102 : return with_connectivity ? m_connected_mesh : m_unique_mesh;
227 : }
228 :
229 15 : Mesh Issm::get_vertices(const MeshConfig& config, bool with_connectivity) const
230 : {
231 15 : auto [connected_mesh, unique_mesh] =
232 30 : std::visit([this](auto& mesh_type) { return make_mesh(mesh_type, *this->m_model); }, config.type);
233 30 : return with_connectivity ? connected_mesh : unique_mesh;
234 15 : }
235 :
236 2878 : auto get_analysis(definitions input)
237 : {
238 2878 : switch (input)
239 : {
240 2878 : case VxEnum:
241 : case VyEnum:
242 2878 : return StressbalanceAnalysisEnum;
243 0 : default:
244 0 : throw std::runtime_error(fmt::format("Constraint coupling not implemented for {}.", EnumToStringx(input)));
245 : }
246 : }
247 :
248 4299 : auto get_constraint_dof_idx(definitions input)
249 : {
250 4299 : switch (input)
251 : {
252 1433 : case VxEnum:
253 1433 : return 0;
254 2866 : case VyEnum:
255 2866 : return 1;
256 0 : default:
257 0 : return 0;
258 : }
259 : }
260 :
261 36 : std::vector<double> Issm::synchronize_ghosts(const std::vector<Vertex>& vertices, std::span<const double> values) const
262 : {
263 : static_assert(std::numeric_limits<double>::has_quiet_NaN);
264 : auto all_values =
265 36 : std::vector<double>(m_model->vertices->NumberOfVerticesLocalAll(), std::numeric_limits<double>::quiet_NaN());
266 36 : auto issm_vertices = as_range<::Vertex>(*m_model->vertices);
267 36 : assert(all_values.size() == issm_vertices.size());
268 5464 : for (size_t vtx_idx = 0; vtx_idx < vertices.size(); ++vtx_idx)
269 : {
270 5428 : auto issm_vertex_idx = vertices[vtx_idx].id;
271 5428 : all_values[issm_vertices[issm_vertex_idx]->Lid()] = values[vtx_idx];
272 : }
273 36 : m_model->SyncLocalVectorWithClonesVertices(all_values.data());
274 72 : return all_values;
275 0 : }
276 :
277 12 : void Issm::set_constraints(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values)
278 : {
279 12 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Setting ISSM constraints for {}", issm_input_name(input));
280 :
281 12 : m_model->SetCurrentConfiguration(get_analysis(input));
282 :
283 : // TODO: this condition might be relaxed if we use coupling meshes based on nodes instead of vertices.
284 : static int allowed_fe_types[] = {P1Enum, P1P1Enum};
285 24 : if (stdr::find(allowed_fe_types, as_range<::Element>(*m_model->elements)[0]->element_type) == std::end(allowed_fe_types))
286 : {
287 0 : throw std::runtime_error("Constraint coupling only implemented for P1 elements (nodes == vertices).");
288 : }
289 :
290 : // set constraints
291 : // !! loop over ALL vertices of the model, vertices in the coupling mesh do not include ghosts/clones !!
292 12 : auto all_vertex_values = synchronize_ghosts(vertices, values);
293 14342 : for (size_t vtx_idx = 0; vtx_idx < all_vertex_values.size(); ++vtx_idx)
294 : {
295 14330 : auto vertex = as_range<::Vertex>(*m_model->vertices)[vtx_idx];
296 14330 : if (std::isnan(all_vertex_values[vertex->Lid()]))
297 : {
298 11464 : continue; // constraint not set
299 : }
300 :
301 : // find constraint on this vertex
302 : // This can probably be optimised since the constraint should have a similar order as the vertices
303 2866 : auto static_constraints = as_range<::Constraint>(*m_model->constraints) |
304 1880498 : stdv::transform([](auto cons) { return dynamic_cast<SpcStatic*>(cons); }) |
305 937383 : stdv::filter([](auto cons) { return cons != nullptr; });
306 2866 : auto constraint = stdr::find_if(
307 : static_constraints,
308 931651 : [&](auto static_constraint)
309 : {
310 935950 : return static_constraint->GetNodeId() == vertex->Sid() + 1 &&
311 935950 : static_constraint->GetDof() == get_constraint_dof_idx(input);
312 : });
313 2866 : if (constraint == end(static_constraints))
314 : {
315 0 : throw std::runtime_error(fmt::format(
316 : "Constraint not found for vertex {}, Sid {}, Value {}.",
317 : vtx_idx,
318 0 : vertex->Sid(),
319 0 : all_vertex_values[vertex->Lid()]));
320 : }
321 : // set constraint
322 8598 : **constraint = SpcStatic(
323 2866 : (**constraint).Id(),
324 : (**constraint).GetNodeId(),
325 : (**constraint).GetDof(),
326 2866 : all_vertex_values[vertex->Lid()],
327 5732 : get_analysis(input));
328 : }
329 12 : }
330 :
331 62 : void Issm::set_input(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values)
332 : {
333 62 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Setting ISSM input values for {}", issm_input_name(input));
334 :
335 62 : if (vertices.size() < m_connected_mesh.vertices.size())
336 : {
337 : // input array does not include ghosts
338 24 : auto all_vertex_values = synchronize_ghosts(vertices, values);
339 :
340 : // set value for each vertex.
341 : // Could probably be made more efficient by setting all values at the same time but
342 : // making this work for many different types of meshes is probably not worth it.
343 3078 : for (size_t i = 0; i < m_connected_mesh.vertices.size(); ++i)
344 : {
345 3054 : auto issm_vertex = as_range<::Vertex>(*m_model->vertices)[m_connected_mesh.vertices[i].id];
346 3054 : auto issm_element = as_range<::Element>(*m_model->elements)[m_connected_mesh.vertices[i].element_id];
347 3054 : IPC_LOG_TRACE(
348 : "Set value for {} at vertex local {} global {} ghost {} [{} {}] = {}",
349 : issm_input_name(input),
350 : issm_vertex->lid,
351 : issm_vertex->sid,
352 : issm_vertex->clone ? 1 : 0,
353 : issm_vertex->x,
354 : issm_vertex->y,
355 : values[i]);
356 3054 : issm_element->SetElementInput(
357 3054 : m_model->inputs,
358 : 1,
359 : &issm_vertex->lid,
360 3054 : const_cast<double*>(&all_vertex_values[issm_vertex->Lid()]),
361 : input);
362 : }
363 24 : }
364 : else
365 : {
366 : // input array include ghosts, set directly
367 : // set value for each vertex.
368 : // Could probably be made more efficient by setting all values at the same time but
369 : // making this work for many different types of meshes is probably not worth it.
370 11776 : for (size_t i = 0; i < vertices.size(); ++i)
371 : {
372 11738 : auto issm_vertex = as_range<::Vertex>(*m_model->vertices)[vertices[i].id];
373 11738 : auto issm_element = as_range<::Element>(*m_model->elements)[vertices[i].element_id];
374 11738 : IPC_LOG_TRACE(
375 : "Set value for {} at vertex local {} global {} ghost {} [{} {}] = {}",
376 : issm_input_name(input),
377 : issm_vertex->lid,
378 : issm_vertex->sid,
379 : issm_vertex->clone ? 1 : 0,
380 : issm_vertex->x,
381 : issm_vertex->y,
382 : values[i]);
383 : issm_element
384 11738 : ->SetElementInput(m_model->inputs, 1, &issm_vertex->lid, const_cast<double*>(&values[i]), input);
385 : }
386 : }
387 62 : }
388 :
389 125 : void Issm::get_input(definitions input, const std::vector<Vertex>& vertices, std::span<double> values) const
390 : {
391 125 : assert(values.size() == vertices.size());
392 :
393 125 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Get ISSM input values for {}", issm_input_name(input));
394 :
395 : //get value for each vertex.
396 : //Could probably be made more efficient by setting all values at the same time but
397 : //making this work for many different types of meshes is probably not worth it.
398 31008 : auto get_value = [&](auto& v)
399 : {
400 31008 : auto vertex = as_range<::Vertex>(*m_model->vertices)[v.id];
401 31008 : auto element = as_range<::Element>(*m_model->elements)[v.element_id];
402 31008 : IPC_LOG_TRACE("Get value of at vertex {} [{} {}]", vertex->lid, vertex->x, vertex->y);
403 : double value;
404 31008 : if (element->GetInput(input))
405 : {
406 31008 : element->GetInputValue(&value, vertex, input);
407 : }
408 : else
409 : {
410 0 : if (input == VelEnum)
411 : {
412 : //Vel doesn't exist yet during initialization
413 : double vx, vy, vz;
414 0 : element->GetInputValue(&vx, vertex, VxEnum);
415 0 : element->GetInputValue(&vy, vertex, VyEnum);
416 0 : element->GetInputValue(&vz, vertex, VzEnum);
417 0 : value = sqrt(vx * vx + vy * vy + vz * vz);
418 : }
419 : else
420 : {
421 0 : IPC_LOG_ERROR("Input {} not found.", issm_input_name(input));
422 0 : throw std::runtime_error("Input not found.");
423 : }
424 : }
425 31008 : 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);
426 :
427 31008 : return value;
428 125 : };
429 125 : std::ranges::transform(vertices, values.begin(), get_value);
430 125 : }
431 :
432 3 : void Issm::extrude(definitions input, const std::vector<Vertex>& /*vertices*/)
433 : {
434 3 : if (auto p_mesh2d = get_if<Mesh2dConfig>(&m_coupling_mesh_config.type); p_mesh2d)
435 : {
436 3 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Extruding ISSM input {}", issm_input_name(input));
437 3 : set_parameter(*m_model, InputToExtrudeEnum, input);
438 3 : if (p_mesh2d->layer == 0)
439 : {
440 2 : extrudefrombase_core(m_model.get());
441 : }
442 1 : else if (p_mesh2d->layer == -1)
443 : {
444 1 : extrudefromtop_core(m_model.get());
445 : }
446 : else
447 : {
448 0 : throw std::runtime_error("Extruding from layer other than base or surface not implemented.");
449 : }
450 : }
451 3 : }
452 :
453 5 : definitions depth_average_input(definitions input)
454 : {
455 5 : switch (input)
456 : {
457 2 : case VxEnum: return VxAverageEnum;
458 2 : case VyEnum: return VyAverageEnum;
459 0 : case WaterfractionDrainageEnum: return WaterfractionDrainageIntegratedEnum;
460 0 : case MaterialsRheologyBEnum: return MaterialsRheologyBbarEnum;
461 0 : case DamageDEnum: return DamageDbarEnum;
462 0 : case DamageDOldEnum: return DamageDbarOldEnum;
463 0 : case MaterialsRheologyEEnum: return MaterialsRheologyEbarEnum;
464 0 : case MaterialsRheologyEsEnum: return MaterialsRheologyEsbarEnum;
465 0 : case MaterialsRheologyEcEnum: return MaterialsRheologyEcbarEnum;
466 1 : default: throw std::domain_error("Depth averaging not implemented for this input.");
467 : }
468 : }
469 :
470 3 : definitions Issm::depth_average(definitions input, const std::vector<Vertex>& /*vertices*/)
471 : {
472 3 : auto average_input = depth_average_input(input);
473 2 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Depth averaging ISSM input {} into {}", issm_input_name(input), issm_input_name(average_input));
474 2 : if (auto p_mesh2d = get_if<Mesh2dConfig>(&m_coupling_mesh_config.type); !p_mesh2d || p_mesh2d->layer != 0)
475 : {
476 0 : IPC_LOG_WARN_0(IssmComm::GetComm(), "Values of depth averaged input {} only valid at base layer.", issm_input_name(average_input));
477 : }
478 :
479 2 : set_parameter(*m_model, InputToDepthaverageInEnum, input);
480 2 : set_parameter(*m_model, InputToDepthaverageOutEnum, average_input);
481 2 : depthaverage_core(m_model.get());
482 2 : return average_input;
483 : }
484 :
485 3 : void Issm::initialize(definitions input)
486 : {
487 :
488 3 : switch (input)
489 : {
490 0 : case MaskIceLevelsetEnum:
491 : case MaskOceanLevelsetEnum:
492 : case ThicknessEnum:
493 : case BaseEnum:
494 : case SurfaceEnum:
495 : case BedEnum:
496 : //need no initialization, must be initialized in a valid setup
497 0 : break;
498 3 : case VxEnum:
499 : case VyEnum:
500 : case VzEnum:
501 : case VelEnum:
502 : case PressureEnum:
503 3 : if (!m_initialized_stress_balance)
504 : {
505 3 : stressbalance_core(m_model.get());
506 3 : m_model->results->clear(); //don't store results
507 3 : m_initialized_stress_balance = true;
508 : }
509 3 : break;
510 0 : default:
511 0 : IPC_LOG_WARN_0(
512 : IssmComm::GetComm(),
513 : "Initializing coupling data of {}. "
514 : "Make sure this field has a valid initial value as the adapter can not compute it. "
515 : "A procedure to compute the initial value can be added to ipc::Issm::initialize(). ",
516 : issm_input_name(input));
517 0 : break;
518 : }
519 3 : }
520 :
521 7 : double Issm::get_max_time_step() const
522 : {
523 : #ifdef ISSM_PRECICE_FEATURE_SUBCYCLING
524 : int ts_type;
525 7 : m_model->parameters->FindParam(&ts_type, TimesteppingTypeEnum);
526 :
527 7 : if (ts_type == AdaptiveTimesteppingEnum)
528 : {
529 : // determine adaptive step size
530 : double adapt_dt;
531 1 : m_model->TimeAdaptx(&adapt_dt);
532 1 : return adapt_dt;
533 : }
534 6 : else if (ts_type == FixedTimesteppingEnum)
535 : {
536 : // use either the saved initial dt or override dt
537 : // don't use the time step that is currently set in ISSM, since it may have been overridden,
538 : // e.g., to match the exact final time of the last coupling window.
539 6 : return m_dt_override > 0 ? m_dt_override : m_initial_dt;
540 : }
541 : else
542 : {
543 0 : throw std::runtime_error("Unknown time stepping scheme.");
544 : }
545 : #else
546 : return std::numeric_limits<double>::infinity();
547 : #endif
548 : }
549 :
550 82 : void Issm::solve(double dt)
551 : {
552 82 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Advance ISSM by {}", dt);
553 :
554 : double t;
555 82 : m_model->parameters->FindParam(&t, TimeEnum);
556 :
557 : // set the time step if not using adaptive time steps.
558 : // must be set before every solve because it may have been overridden internally,
559 : // e.g., to match the final time exactly.
560 82 : if (get_parameter<::definitions>(*m_model, TimesteppingTypeEnum) == FixedTimesteppingEnum)
561 : {
562 81 : set_parameter(*m_model, TimesteppingTimeStepEnum, m_dt_override > 0.0 ? m_dt_override : m_initial_dt);
563 : }
564 :
565 : // set or disable output frequency
566 82 : if (m_output_frequency_override > 0)
567 : {
568 1 : set_parameter(*m_model, SettingsOutputFrequencyEnum, m_output_frequency_override);
569 : }
570 81 : else if (m_output_frequency_override == 0)
571 : {
572 : // disable output
573 1 : set_parameter(*m_model, SettingsOutputFrequencyEnum, std::numeric_limits<int>::max());
574 : #ifndef ISSM_PRECICE_FEATURE_CONFIGURE_OUTPUT
575 : IPC_LOG_WARN_0(IssmComm::GetComm(), "Disabling output entirely requires ISSM version after 2026.2.");
576 : #endif
577 : }
578 :
579 82 : set_parameter(*m_model, TimesteppingStartTimeEnum, t);
580 82 : auto final_t = t + dt;
581 82 : set_parameter(*m_model, TimesteppingFinalTimeEnum, final_t);
582 :
583 82 : IPC_LOG_INFO_0(IssmComm::GetComm(), "Run ISSM from t = {} to {}", t, t + dt);
584 :
585 82 : m_model->Solve();
586 :
587 : // check that ISSM has hit the exact final time.
588 : // there is a bug in issm that makes it overshoot.
589 : // workaround: set the time step so the coupling window
590 : // is an integer multiple of the time step.
591 : // bug should be fixed in ISSM version 2026.1
592 82 : auto new_t = get_parameter<double>(*m_model, TimeEnum);
593 82 : if (!fp_equal(new_t, final_t, 0, 1e-10))
594 : {
595 0 : IPC_LOG_WARN_0(
596 : IssmComm::GetComm(),
597 : "ISSM overshot the expected final time by {} s. Check time step settings.",
598 : new_t - final_t);
599 : }
600 82 : }
601 :
602 6 : void Issm::set_time_step_override(double dt)
603 : {
604 6 : m_dt_override = dt;
605 6 : }
606 :
607 2 : void Issm::set_output_frequency_override(int f)
608 : {
609 2 : m_output_frequency_override = f;
610 2 : }
611 :
612 13 : double Issm::get_time() const
613 : {
614 13 : return get_parameter<double>(*m_model, TimeEnum);
615 : }
616 :
617 13 : int Issm::get_step() const
618 : {
619 13 : return get_parameter<int>(*m_model, StepEnum);
620 : }
621 :
622 : } // namespace ipc
|
