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