Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Daniel Abele <daniel.abele@dlr.de>
       2             : //
       3             : // SPDX-License-Identifier: BSD-3-Clause
       4             : 
       5             : #pragma once
       6             : 
       7             : #include "issm-precice/confreader.hpp"
       8             : #include "mpi.h"
       9             : #include <filesystem>
      10             : #include <memory>
      11             : #include <ranges>
      12             : #include <string>
      13             : #include <vector>
      14             : #include <span>
      15             : 
      16             : class FemModel;
      17             : 
      18             : namespace ipc
      19             : {
      20             : 
      21             : struct NoopDeleter
      22             : {
      23          34 :     void operator()(FemModel*) { }
      24             : };
      25             : 
      26             : /**
      27             :  * Unique ptr to own an ISSM FemModel object.
      28             :  * Deleter does not do anything because FemModel causes a segfault on deletion.
      29             :  */
      30             : using FemModelPtr = std::unique_ptr<FemModel, NoopDeleter>;
      31             : 
      32             : /**
      33             :  * Create an ISSM FemModel by loading the specified model.
      34             :  * @param root_path Directory where the model files are stored.
      35             :  * @param model_name Name of the model to load.
      36             :  * @param comm MPI communicator to use by ISSM.
      37             :  */
      38             : FemModelPtr make_fem_model(const std::filesystem::path& root_path, const std::string& model_name, MPI_Comm comm);
      39             : 
      40             : /**
      41             :  * Get the ISSM input id of the given name.
      42             :  */
      43             : definitions issm_input_id(const std::string& input_name);
      44             : 
      45             : /**
      46             :  * Get the name of the ISSM input.
      47             :  */
      48             : const char* issm_input_name(definitions input_id);
      49             : 
      50             : /**
      51             :  * Get the ISSM input that stores depth averaged values for the specified input.
      52             :  */
      53             : definitions depth_average_input(definitions input_id);
      54             : 
      55             : /**
      56             :  * Vertex of the ISSM mesh.
      57             :  */
      58             : struct Vertex
      59             : {
      60             :     /**
      61             :      * Locally unique id of the vertex.
      62             :      */
      63       92160 :     size_t id = 0;
      64             : 
      65             :     /**
      66             :      * Coordinates of the vertex.
      67             :      */
      68       92160 :     std::array<double, 3> coordinates;
      69             : 
      70             :     /**
      71             :      * Locally unique id of the element (2d triangle or 3d tetraeder) this vertex belongs to.
      72             :      */    
      73        4480 :     size_t element_id = 0;
      74             : 
      75       92160 :     bool operator==(const Vertex&) const = default;
      76             : };
      77             : 
      78             : /**
      79             :  * Abstract interface for ISSM used for coupling.
      80             :  */
      81             : class IIssm
      82             : {
      83             : public:
      84             :     /**
      85             :      * Virtual destructor.
      86             :      */
      87             :     virtual ~IIssm();
      88             : 
      89             :     /**
      90             :      * Get the vertices of the coupled interface mesh.
      91             :      * @returns Vertices of the ISSM mesh and the elements they belong to.
      92             :      */
      93             :     virtual std::vector<Vertex> get_vertices() const = 0;
      94             : 
      95             :     /**
      96             :      * Maximum length of next time step allowed by the model (e.g. for numeric stability).
      97             :      */
      98             :     virtual double get_max_time_step() const = 0;
      99             : 
     100             :     /**
     101             :      * Set the specified input for the vertices.
     102             :      * @param input Input to set.
     103             :      * @param vertices Vertices and corresponding Element to set the values for.
     104             :      * @param values Values to set.
     105             :      */
     106             :     virtual void
     107             :     set_input(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values) = 0;
     108             : 
     109             :     /**
     110             :      * Set the specified constraints for the vertices.
     111             :      * @param input Input to set.
     112             :      * @param vertices Vertices and corresponding Element to set the values for.
     113             :      * @param values Values to set.
     114             :      */
     115             :     virtual void
     116             :     set_constraints(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values) = 0;
     117             : 
     118             :     /**
     119             :      * Extrude the specified input.
     120             :      * Values of the input are copied from base (or surface, depending on the used coupling interface)
     121             :      * to the other layers of a 3D ISSM model. No effect on 2D ISSM models.
     122             :      * @param input Input to extrude.
     123             :      * @param vertices Vertices to extrude.
     124             :      */
     125             :     virtual void extrude(definitions input, const std::vector<Vertex>& vertices) = 0;
     126             : 
     127             :     /**
     128             :      * Compute valid initial values for the input.
     129             :      * For most inputs, this does nothing since the input has valid values in
     130             :      * the ISSM model setup. Some inputs are not required to be have correct initial values as they are 
     131             :      * computed by ISSM without initial state. For these inputs, this function does the necessary computation that
     132             :      * would normally happen only during the simulation.
     133             :      * @param input Input to initialize.
     134             :      */
     135             :     virtual void initialize(definitions input) = 0;
     136             : 
     137             :     /**
     138             :      * Get the specified input for the vertices.
     139             :      * @param input Input to get.
     140             :      * @param vertices Vertices and corresponding Element to get the values for.
     141             :      * @param[out] values Vector to store the values in. Will be overwritten by the function.
     142             :      */
     143             :     virtual void
     144             :     get_input(definitions input, const std::vector<Vertex>& vertices, std::span<double> values) const = 0;
     145             : 
     146             :     /**
     147             :      * Run depth averaging for the specified input.
     148             :      * Computes the average of an input over all the layers of a 3D ISSM model and stores it
     149             :      * in a specific input that depends on the argument input. For every supported input there is one
     150             :      * corresponding average input.
     151             :      * @param input Input to depth average.
     152             :      * @param vertices vertices to depth average.
     153             :      * @returns Input where the depth averaged values are stored.
     154             :      */
     155             :     virtual definitions depth_average(definitions input, const std::vector<Vertex>& vertices) = 0;
     156             : 
     157             :     /**
     158             :      * Advance the simulation one time step.
     159             :      * @param dt Length of the time step.
     160             :      */
     161             :     virtual void solve(double dt) = 0;
     162             : };
     163             : 
     164             : /**
     165             :  * Concretization of the IIssm interface.
     166             :  * Gives access to ISSM model data required for coupling.
     167             :  */
     168             : class Issm : public IIssm
     169             : {
     170             : public:
     171             :     /**
     172             :      * Initialize ISSM by loading a model.
     173             :      * @param root_path Path to directory that stores the model.
     174             :      * @param model_name Name of the model.
     175             :      * @param mesh Definition of the coupling interface mesh.
     176             :      * @param comm MPI communicator to used by ISSM.
     177             :      */
     178             :     Issm(
     179             :         const std::filesystem::path& root_path,
     180             :         const std::string& model_name,
     181             :         const MeshDefinition& mesh,
     182             :         MPI_Comm comm);
     183             : 
     184             :     /**
     185             :      * @copydoc IIssm::get_vertices()
     186             :      */
     187             :     std::vector<Vertex> get_vertices() const override;
     188             : 
     189             :     /**
     190             :      * @copydoc IIssm::get_vertices(const MeshDefinition&)
     191             :      */
     192             :     std::vector<Vertex> get_vertices(const MeshDefinition& mesh) const;
     193             :     
     194             :     /**
     195             :      * @copydoc IIssm::get_time_step()
     196             :      */
     197             :     double get_max_time_step() const override;
     198             : 
     199             :     /**
     200             :      * @copydoc IIssm::initialize()
     201             :      */
     202             :     void initialize(definitions input) override;
     203             : 
     204             :     /**
     205             :      * @copydoc IIssm:set_input()
     206             :      */
     207             :     void set_input(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values) override;
     208             : 
     209             :     /**
     210             :      * @copydoc IIssm:set_constraints()
     211             :      */
     212             :     virtual void
     213             :     set_constraints(definitions input, const std::vector<Vertex>& vertices, std::span<const double> values) override;
     214             : 
     215             :     /**
     216             :      * @copydoc IIssm:extrude()
     217             :      */
     218             :     void extrude(definitions input, const std::vector<Vertex>& vertices) override;
     219             : 
     220             :     /**
     221             :      * @copydoc IIssm:get_input()
     222             :      */
     223             :     void get_input(definitions input, const std::vector<Vertex>& vertices, std::span<double> values) const override;
     224             : 
     225             :     /**
     226             :      * @copydoc IIssm:depth_average()
     227             :      */
     228             :     definitions depth_average(definitions input, const std::vector<Vertex>& vertices) override;
     229             : 
     230             :     /**
     231             :      * @copydoc IIssm:solve()
     232             :      */
     233             :     void solve(double dt) override;
     234             : 
     235             : private:
     236             :     bool try_set_constraint_value(const Vertex&, definitions input, double value) const;
     237             :     std::vector<Vertex> get_vertices(const Mesh2dDefinition& mesh_2d) const;
     238             :     std::vector<Vertex> get_vertices(const Mesh3dDefinition& mesh_3d) const;
     239             : 
     240             : private:
     241             :     FemModelPtr m_model;
     242             :     MeshDefinition m_mesh_definition;
     243             :     double m_max_dt;
     244             :     bool m_initialized_stress_balance;
     245             : };
     246             : } // namespace ipc