Fix interfacing with svSolver/svFSI (#102)

Co-authored-by: Karthik Menon <kmenon13@sh03-ln06.stanford.edu>
Co-authored-by: Karthik Menon <kmenon13@sh02-ln01.stanford.edu>

src/model/Parameter.h

@@ -36,6 +36,7 @@
 
 #include <math.h>
 
+#include <iostream>
 #include <numeric>
 #include <vector>
 

src/solve/SimulationParameters.cpp

@@ -354,8 +354,10 @@ void create_external_coupling(
     std::string coupling_loc = coupling_config["location"];
     bool periodic = coupling_config.value("periodic", true);
     const auto& coupling_values = coupling_config["values"];
+    const bool internal = false;
 
-    generate_block(model, coupling_values, coupling_type, coupling_name);
+    generate_block(model, coupling_values, coupling_type, coupling_name,
+                   internal, periodic);
 
     // Determine the type of connected block
     std::string connected_block = coupling_config["connected_block"];

tests/test_interface/CMakeLists.txt

@@ -5,3 +5,4 @@ project(svZeroDSolverInterfaceTests)
 
 add_subdirectory("test_01/")
 add_subdirectory("test_02/")
+add_subdirectory("test_03/")

tests/test_interface/test_01/main.cpp

@@ -1,4 +1,7 @@
 // Test interfacing to svZeroSolver. 
+// This test mimics an external 3D solver (svSolver/svFSI) interfacing with svZeroDSolver
+// The model consists of a closed-loop heart model with coronary BCs
+// It is run for one time step of the external solver
 
 #include "../LPNSolverInterface/LPNSolverInterface.h" 
 #include <iostream>

tests/test_interface/test_02/main.cpp

@@ -1,4 +1,6 @@
-// Test interfacing to svZeroDSolver. 
+// Test interfacing to svZeroDSolver.
+// This test mimics the coupling of svZeroDSolver with an external parameter estimation code (eg. Tulip)
+// A coronary model is used and parameters of the BCs are updated and then compared to a reference solution
 
 #include "../LPNSolverInterface/LPNSolverInterface.h" 
 #include <iostream>

tests/test_interface/test_03/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_executable(svZeroD_interface_test03 ../LPNSolverInterface/LPNSolverInterface.cpp  main.cpp)
+target_link_libraries(svZeroD_interface_test03 ${CMAKE_DL_LIBS})

tests/test_interface/test_03/main.cpp

@@ -0,0 +1,151 @@
+// Test interfacing to svZeroSolver.
+// This test mimics an external 3D solver (svSolver/svFSI) interfacing with svZeroDSolver
+// The model consists of an RCR BC which acts as a Neumann BC for an external solver
+// It mimics two consecutive time steps of an external solver
+
+#include "../LPNSolverInterface/LPNSolverInterface.h" 
+#include <iostream>
+#include <map>
+#include <fstream>
+
+//------
+// main
+//------
+//
+int main(int argc, char** argv)
+{
+  LPNSolverInterface interface;
+
+  if (argc != 3) {
+    std::runtime_error("Usage: svZeroD_interface_test03 <path_to_svzeroDSolver_build_folder> <path_to_json_file>");
+  }
+
+  // Load shared library and get interface functions.
+  // File extension of the shared library depends on the system
+  std::string svzerod_build_path = std::string(argv[1]);
+  std::string interface_lib_path = svzerod_build_path + "/src/interface/libsvzero_interface";
+  std::string interface_lib_so = interface_lib_path + ".so";
+  std::string interface_lib_dylib = interface_lib_path + ".dylib";
+  std::ifstream lib_so_exists(interface_lib_so);
+  std::ifstream lib_dylib_exists(interface_lib_dylib);
+  if (lib_so_exists) {
+    interface.load_library(interface_lib_so);
+  } else if (lib_dylib_exists) {
+    interface.load_library(interface_lib_dylib);
+  } else {
+    throw std::runtime_error("Could not find shared libraries " + interface_lib_so + " or " + interface_lib_dylib);
+  } 
+
+  // Set up the svZeroD model
+  std::string file_name = std::string(argv[2]);
+  interface.initialize(file_name);
+
+  // Check number of variables and blocks
+  if (interface.system_size_ != 3) {
+    throw std::runtime_error("interface.system_size_ != 3");
+  }
+  if (interface.block_names_.size() != 2) {
+    throw std::runtime_error("interface.block_names_.size() != 2");
+  }
+
+  // Set external time step size (flow solver step size)
+  double external_step_size = 0.005;
+  interface.set_external_step_size(external_step_size);
+
+  // Save the initial condition
+  std::vector<double> init_state_y = {-6.2506662304695681e+01, -3.8067539421845140e+04, -3.0504233282976966e+04};
+  std::vector<double> init_state_ydot = {-3.0873806830951793e+01, -2.5267653962355386e+05, -2.4894080899699836e+05};
+
+  // Get variable IDs for inlet to RCR block
+  std::vector<int> IDs;
+  std::string block_name = "RCR"; 
+  interface.get_block_node_IDs(block_name, IDs);
+  int num_inlet_nodes = IDs[0];
+  int num_outlet_nodes = IDs[1+num_inlet_nodes*2];
+  if ((num_inlet_nodes != 1) || (num_outlet_nodes != 0)) {
+    throw std::runtime_error("Wrong number of inlets/outlets for RCR");
+  }
+  int rcr_inlet_flow_id = IDs[1];
+  int rcr_inlet_pressure_id = IDs[2];
+
+  // Update block parameters with current flow from 3D solver
+  std::vector<double> new_params(5);
+  std::vector<double> params = {-6.2506662041472836e+01, -6.2599344518688739e+01};
+  std::vector<double> interface_times = {1.9899999999999796e+00, 1.9949999999999795e+00};
+  // Format of new_params for flow/pressure blocks: 
+  // [N, time_1, time_2, ..., time_N, value1, value2, ..., value_N]
+  // where N is number of time points and value* is flow/pressure
+  new_params[0] = 2.0;
+  for (int i = 0; i < 2; i++) {
+    new_params[1+i] = interface_times[i];
+    new_params[3+i] = params[i];
+  }
+  interface.update_block_params("RCR_coupling", new_params); 
+
+  // Set up vectors to run svZeroD simulation
+  std::vector<double> solutions(interface.system_size_*interface.num_output_steps_);
+  std::vector<double> times(interface.num_output_steps_);
+  int error_code = 0;
+
+  // Run svZeroD simulation
+  interface.update_state(init_state_y, init_state_ydot); 
+  interface.run_simulation(interface_times[0], times, solutions, error_code);
+
+  // Parse output and calculate mean flow/pressure in aorta and coronary
+  int sol_idx = 0;
+  double mean_pressure = 0.0;
+  for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) {
+    for (int state = 0; state < interface.system_size_; state++) {
+      sol_idx = interface.system_size_*tstep + state;
+      if (state == rcr_inlet_pressure_id) {
+        mean_pressure += solutions[sol_idx];
+      }
+    }
+  }
+  mean_pressure /= (double)interface.num_output_steps_;
+  std::cout <<"Simulation output: " << std::endl;
+  std::cout <<"Mean pressure = " << mean_pressure << std::endl;
+
+  // Compare mean pressure with pre-computed ("correct") values
+  double error_limit = 0.05;
+  if (abs(-mean_pressure / 38690.2 - 1.0) > error_limit) {
+    throw std::runtime_error("Error in mean pressure at RCR inlet.");
+  }
+
+  // Get state vector to prepare for next 3D time step
+  interface.return_y(init_state_y);
+  interface.return_ydot(init_state_ydot);
+
+  // Update parameters for next 3D time step
+  params = {-6.2599344283486374e+01, -6.2630248751732964e+01};
+  interface_times = {1.9949999999999795e+00, 1.9999999999999793e+00};
+  for (int i = 0; i < 2; i++) {
+    new_params[1+i] = interface_times[i];
+    new_params[3+i] = params[i];
+  }
+  interface.update_block_params("RCR_coupling", new_params); 
+
+  // Run svZeroD simulation
+  interface.update_state(init_state_y, init_state_ydot); 
+  interface.run_simulation(interface_times[0], times, solutions, error_code);
+
+  // Parse output and calculate mean flow/pressure in aorta and coronary
+  sol_idx = 0;
+  mean_pressure = 0.0;
+  for (int tstep = 0; tstep < interface.num_output_steps_; tstep++) {
+    for (int state = 0; state < interface.system_size_; state++) {
+      sol_idx = interface.system_size_*tstep + state;
+      if (state == rcr_inlet_pressure_id) {
+        mean_pressure += solutions[sol_idx];
+      }
+    }
+  }
+  mean_pressure /= (double)interface.num_output_steps_;
+  std::cout <<"Simulation output: " << std::endl;
+  std::cout <<"Mean pressure = " << mean_pressure << std::endl;
+
+  // Compare mean pressure with pre-computed ("correct") values
+  if (abs(-mean_pressure / 39911.3 - 1.0) > error_limit) {
+    throw std::runtime_error("Error in mean pressure at RCR inlet.");
+  }
+}

tests/test_interface/test_03/svzerod_3Dcoupling.json

@@ -0,0 +1,35 @@
+{
+    "simulation_parameters": {
+        "coupled_simulation": true,
+        "number_of_time_pts": 50,
+        "output_all_cycles": true,
+        "steady_initial": false
+    },
+    "boundary_conditions": [
+        {
+            "bc_name": "RCR",
+            "bc_type": "RCR",
+            "bc_values": {
+                "Rp": 121.0,
+                "Rd": 1212.0,
+                "C": 1.5e-4,
+                "Pd": 0.0
+            }
+        }
+    ],
+    "external_solver_coupling_blocks": [
+        {
+            "name": "RCR_coupling",
+            "type": "FLOW",
+            "location": "inlet",
+            "connected_block": "RCR",
+            "periodic": false,
+            "values": {
+                "t": [0.0, 1.0],
+                "Q": [1.0, 1.0]
+            }
+        }
+    ],
+    "junctions": [],
+    "vessels": []
+}