Merge pull request #347 from Exawind/restart_from_file

Added restart capability for CFD Interface.

src/interfaces/cfd/interface.cpp

@@ -7,8 +7,10 @@
 #include "node_data.hpp"
 #include "state/clone_state.hpp"
 #include "state/copy_state_data.hpp"
+#include "state/read_state_from_file.hpp"
 #include "state/set_node_external_loads.hpp"
 #include "state/state.hpp"
+#include "state/write_state_to_file.hpp"
 #include "step/step.hpp"
 #include "turbine.hpp"
 #include "turbine_input.hpp"
@@ -200,6 +202,26 @@ Interface::Interface(const InterfaceInput& input)
     );
 }
 
+void Interface::WriteRestart(const std::filesystem::path& filename) const {
+    auto output = std::ofstream(filename);
+    WriteStateToFile(output, state);
+}
+
+void Interface::ReadRestart(const std::filesystem::path& filename) {
+    auto input = std::ifstream(filename);
+    ReadStateFromFile(input, state);
+
+    Kokkos::deep_copy(this->host_state_x, this->state.x);
+    Kokkos::deep_copy(this->host_state_q, this->state.q);
+    Kokkos::deep_copy(this->host_state_v, this->state.v);
+    Kokkos::deep_copy(this->host_state_vd, this->state.vd);
+
+    GetTurbineMotion(
+        this->turbine, this->host_state_x, this->host_state_q, this->host_state_v,
+        this->host_state_vd
+    );
+}
+
 bool Interface::Step() {
     // Transfer loads to solver
     SetTurbineLoads(this->turbine, this->state);

src/interfaces/cfd/interface.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <filesystem>
+
 #include "interfaces/cfd/interface_input.hpp"
 #include "interfaces/cfd/turbine.hpp"
 #include "model/model.hpp"
@@ -19,6 +21,12 @@ class Interface {
     /// @brief Restore state for correction step
     void RestoreState();
 
+    /// @brief Write restart file
+    void WriteRestart(const std::filesystem::path& filename) const;
+
+    /// @brief Read restart file
+    void ReadRestart(const std::filesystem::path& filename);
+
     /// @brief  OpenTurbine class used for model construction
     Model model;
 

src/state/read_state_from_file.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <istream>
+
+#include <Kokkos_Core.hpp>
+
+#include "state.hpp"
+
+namespace openturbine {
+
+inline void ReadStateFromFile(std::istream& input, State& state) {
+    auto num_system_nodes = size_t{};
+    // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+    input.read(reinterpret_cast<char*>(&num_system_nodes), sizeof(size_t));
+
+    if (num_system_nodes != state.num_system_nodes) {
+        throw std::length_error("Number of system nodes in file is not the same as in model");
+    }
+
+    auto mirror_7 = Kokkos::View<double* [7]>::HostMirror("mirror_7", num_system_nodes);
+    auto out_7 = Kokkos::View<double* [7], Kokkos::HostSpace>("out_7", num_system_nodes);
+
+    auto mirror_6 = Kokkos::View<double* [6]>::HostMirror("mirror_6", num_system_nodes);
+    auto out_6 = Kokkos::View<double* [6], Kokkos::HostSpace>("out_6", num_system_nodes);
+
+    auto read_7 = [&](const Kokkos::View<double* [7]>& data) {
+        const auto stream_size = static_cast<long>(7U * num_system_nodes * sizeof(double));
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        input.read(reinterpret_cast<char*>(out_7.data()), stream_size);
+
+        Kokkos::deep_copy(mirror_7, out_7);
+        Kokkos::deep_copy(data, mirror_7);
+    };
+
+    auto read_6 = [&](const Kokkos::View<double* [6]>& data) {
+        const auto stream_size = static_cast<long>(6U * num_system_nodes * sizeof(double));
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        input.read(reinterpret_cast<char*>(out_6.data()), stream_size);
+
+        Kokkos::deep_copy(mirror_6, out_6);
+        Kokkos::deep_copy(data, mirror_6);
+    };
+
+    read_7(state.x0);
+    read_7(state.x);
+    read_6(state.q_delta);
+    read_7(state.q_prev);
+    read_7(state.q);
+    read_6(state.v);
+    read_6(state.vd);
+    read_6(state.a);
+    read_6(state.f);
+}
+
+}  // namespace openturbine

src/state/write_state_to_file.hpp

@@ -0,0 +1,51 @@
+#pragma once
+
+#include <ostream>
+
+#include <Kokkos_Core.hpp>
+
+#include "state.hpp"
+
+namespace openturbine {
+
+inline void WriteStateToFile(std::ostream& output, const State& state) {
+    auto num_system_nodes = state.num_system_nodes;
+    // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+    output.write(reinterpret_cast<char*>(&num_system_nodes), sizeof(size_t));
+
+    auto mirror_7 = Kokkos::View<double* [7]>::HostMirror("mirror_7", num_system_nodes);
+    auto out_7 = Kokkos::View<double* [7], Kokkos::HostSpace>("out_7", num_system_nodes);
+
+    auto mirror_6 = Kokkos::View<double* [6]>::HostMirror("mirror_6", num_system_nodes);
+    auto out_6 = Kokkos::View<double* [6], Kokkos::HostSpace>("out_6", num_system_nodes);
+
+    auto write_7 = [&](const Kokkos::View<double* [7]>& data) {
+        Kokkos::deep_copy(mirror_7, data);
+        Kokkos::deep_copy(out_7, mirror_7);
+
+        const auto stream_size = static_cast<long>(7U * num_system_nodes * sizeof(double));
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        output.write(reinterpret_cast<char*>(out_7.data()), stream_size);
+    };
+
+    auto write_6 = [&](const Kokkos::View<double* [6]>& data) {
+        Kokkos::deep_copy(mirror_6, data);
+        Kokkos::deep_copy(out_6, mirror_6);
+
+        const auto stream_size = static_cast<long>(6U * num_system_nodes * sizeof(double));
+        // NOLINTNEXTLINE(cppcoreguidelines-pro-type-reinterpret-cast)
+        output.write(reinterpret_cast<char*>(out_6.data()), stream_size);
+    };
+
+    write_7(state.x0);
+    write_7(state.x);
+    write_6(state.q_delta);
+    write_7(state.q_prev);
+    write_7(state.q);
+    write_6(state.v);
+    write_6(state.vd);
+    write_6(state.a);
+    write_6(state.f);
+}
+
+}  // namespace openturbine

tests/regression_tests/interfaces/test_cfd_interface.cpp

@@ -275,4 +275,124 @@ TEST(CFDInterfaceTest, FloatingPlatform) {
     }
 }
 
+TEST(CFDInterfaceTest, Restart) {
+    // Solution parameters
+    constexpr auto time_step = 0.01;
+    constexpr auto rho_inf = 0.0;
+    constexpr auto max_iter = 5;
+
+    // Construct platform mass matrix
+    constexpr auto platform_mass{1.419625E+7};                           // kg
+    constexpr Array_3 gravity{0., 0., -9.8124};                          // m/s/s
+    constexpr Array_3 platform_moi{1.2898E+10, 1.2851E+10, 1.4189E+10};  // kg*m*m
+    constexpr Array_3 platform_cm_position{0., 0., -7.53};               // m
+    constexpr Array_6x6 platform_mass_matrix{{
+        {platform_mass, 0., 0., 0., 0., 0.},    // Row 1
+        {0., platform_mass, 0., 0., 0., 0.},    // Row 2
+        {0., 0., platform_mass, 0., 0., 0.},    // Row 3
+        {0., 0., 0., platform_moi[0], 0., 0.},  // Row 4
+        {0., 0., 0., 0., platform_moi[1], 0.},  // Row 5
+        {0., 0., 0., 0., 0., platform_moi[2]},  // Row 6
+    }};
+
+    // Mooring line properties
+    constexpr auto mooring_line_stiffness{48.9e3};       // N
+    constexpr auto mooring_line_initial_length{55.432};  // m
+
+    auto interface_input = InterfaceInput{
+        gravity,
+        time_step,  // time step
+        rho_inf,    // rho infinity (numerical damping)
+        max_iter,   // max convergence iterations
+        TurbineInput{
+            FloatingPlatformInput{
+                true,  // enable
+                {
+                    platform_cm_position[0],
+                    platform_cm_position[1],
+                    platform_cm_position[2],
+                    1.,
+                    0.,
+                    0.,
+                    0.,
+                },                         // position
+                {0., 0., 0., 0., 0., 0.},  // velocity
+                {0., 0., 0., 0., 0., 0.},  // acceleration
+                platform_mass_matrix,
+                {
+                    {
+                        mooring_line_stiffness,
+                        mooring_line_initial_length,
+                        {-40.87, 0.0, -14.},    // Fairlead node coordinates
+                        {0., 0., 0.},           // Fairlead node velocity
+                        {0., 0., 0.},           // Fairlead node acceleration
+                        {-105.47, 0.0, -58.4},  // Anchor node coordinates
+                        {0., 0., 0.},           // Anchor node velocity
+                        {0., 0., 0.},           // Anchor node acceleration
+                    },
+                    {
+                        mooring_line_stiffness,
+                        mooring_line_initial_length,
+                        {20.43, -35.39, -14.},   // Fairlead node coordinates
+                        {0., 0., 0.},            // Fairlead node velocity
+                        {0., 0., 0.},            // Fairlead node acceleration
+                        {52.73, -91.34, -58.4},  // Anchor node coordinates
+                        {0., 0., 0.},            // Anchor node velocity
+                        {0., 0., 0.},            // Anchor node acceleration
+                    },
+                    {
+                        mooring_line_stiffness,
+                        mooring_line_initial_length,
+                        {20.43, 35.39, -14.},   // Fairlead node coordinates
+                        {0., 0., 0.},           // Fairlead node velocity
+                        {0., 0., 0.},           // Fairlead node acceleration
+                        {52.73, 91.34, -58.4},  // Anchor node coordinates
+                        {0., 0., 0.},           // Anchor node velocity
+                        {0., 0., 0.},           // Anchor node acceleration
+                    },
+                },
+            },
+        },
+    };
+
+    auto interface1 = Interface(interface_input);
+
+    // Take 10 initial steps
+    for (auto i = 0U; i < 100U; ++i) {
+        const auto t = static_cast<double>(i) * time_step;
+        interface1.turbine.floating_platform.node.loads[1] = 1e6 * sin(2. * M_PI / 20. * t);
+        auto converged = interface1.Step();
+        EXPECT_TRUE(converged);
+    }
+
+    interface1.WriteRestart("test_restart.dat");
+
+    // Take 10 more steps using original system
+    for (auto i = 0U; i < 100U; ++i) {
+        const auto t = static_cast<double>(i) * time_step;
+        interface1.turbine.floating_platform.node.loads[1] = 1e6 * sin(2. * M_PI / 20. * t);
+        auto converged = interface1.Step();
+        EXPECT_TRUE(converged);
+    }
+
+    auto interface2 = Interface(interface_input);
+    interface2.ReadRestart("test_restart.dat");
+
+    // Take 10 steps using restarted system
+    for (auto i = 0U; i < 100U; ++i) {
+        const auto t = static_cast<double>(i) * time_step;
+        interface2.turbine.floating_platform.node.loads[1] = 1e6 * sin(2. * M_PI / 20. * t);
+        auto converged = interface2.Step();
+        EXPECT_TRUE(converged);
+    }
+
+    const auto& platform_u_1 = interface1.turbine.floating_platform.node.displacement;
+    const auto& platform_u_2 = interface2.turbine.floating_platform.node.displacement;
+    // Ensure platform location is the same for original and restarted system
+    for (auto i = 0U; i < 7U; ++i) {
+        EXPECT_EQ(platform_u_1[i], platform_u_2[i]);
+    }
+
+    std::filesystem::remove("test_restart.dat");
+}
 }  // namespace openturbine::tests