svZeroD Solver Coupling Interface (SimVascular#217)

* Transfer from kharold23/svFSIplus_OLD * sv0D tests * More test case files * Add copyright headers * add cases to python test scripts * fix typo in python test scripts * Add test cases to macos workflow * Added new files to CMakeLists.txt * Update test_fluid.py * Update test_fluid.py * More test updates * More test case debugging * Fix versioning issue and sv0D test cases * Update struct and ustruct test case fields * Fix svZeroD library finding * Add sv0D build to ubuntu yaml * Try setting up conda for sv0D build * Trying to fix sv0D build in the docker * Update test_ubuntu.yml * Update test_ubuntu.yml * Update test_ubuntu.yml * Update test_ubuntu.yml * Update test_ubuntu.yml * Update test_ubuntu.yml * Add "Linear Algebra" element to coupling test cases * Decrease tolerances in ustruct sv0D test case
mrp089 · Jun 19, 2024 · e675676 · e675676
1 parent 64789b6
commit e675676
Show file tree

Hide file tree

Showing 71 changed files with 3,838 additions and 8 deletions.
diff --git a/.github/workflows/test_macos12.yml b/.github/workflows/test_macos12.yml
@@ -12,6 +12,14 @@ jobs:
  brew install lcov
  sudo ln -s /usr/local/opt/qt5/mkspecs /usr/local/mkspecs
  sudo ln -s /usr/local/opt/qt5/plugins /usr/local/plugins
+ - name: Build svZeroDSolver
+ run: |
+ git clone https://github.com/SimVascular/svZeroDSolver.git
+ cd svZeroDSolver
+ mkdir build
+ cd build
+ cmake ..
+ make -j2
  - name: Build svFSIplus
  run: |
  mkdir build

diff --git a/.github/workflows/test_ubuntu.yml b/.github/workflows/test_ubuntu.yml
@@ -16,6 +16,18 @@ jobs:
  container: dcodoni/lib:${{ matrix.image }}
  steps:
  - uses: actions/checkout@v3
+ - name: Build svZeroDSolver
+ shell: bash
+ run: |
+ conda init
+ source /conda/etc/profile.d/conda.sh
+ conda activate svfsiplus
+ git clone https://github.com/SimVascular/svZeroDSolver.git
+ cd svZeroDSolver
+ mkdir build
+ cd build
+ cmake ..
+ make -j2
  - name: Build svFSIplus
  run: |
  if [ -d "build" ]; then

diff --git a/.gitignore b/.gitignore
@@ -28,3 +28,10 @@ genBC.exe
 AllData
 GenBC.int
 InitialData
+
+# svZeroD files
+svZeroD_data
+Q_svZeroD
+P_svZeroD
+# local svZeroD build
+svZeroDSolver
diff --git a/Code/Source/svFSI/CMakeLists.txt b/Code/Source/svFSI/CMakeLists.txt
@@ -201,6 +201,7 @@ set(CSRCS
  shells.h shells.cpp
  stokes.h stokes.cpp
  sv_struct.h sv_struct.cpp
+ svZeroD_subroutines.h svZeroD_subroutines.cpp
  txt.h txt.cpp
  utils.h utils.cpp
  ustruct.h ustruct.cpp
@@ -218,6 +219,8 @@ set(CSRCS
  Timer.h
 
  SPLIT.c
+
+ svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp
 )
 
  # Set PETSc interace code.

diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h
@@ -741,7 +741,10 @@ class cplBCType
  /// @brief Whether to use genBC
  bool useGenBC = false;
 
- /// @brief Whether to initialize RCR from flow data
+ // Whether to use svZeroD
+ bool useSvZeroD = false;
+
+ // Whether to initialize RCR from flow data
  bool initRCR = false;
 
  /// @brief Number of coupled faces

diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp
@@ -892,6 +892,48 @@ void CoupleGenBCParameters::set_values(tinyxml2::XMLElement* xml_elem)
  value_set = true;
 }
 
+
+//////////////////////////////////////////////////////////
+// CoupleSvZeroDParameters //
+//////////////////////////////////////////////////////////
+
+// Coupling to svZeroD.
+
+// Define the XML element name for equation Couple_to_svZeroD parameters.
+const std::string CoupleSvZeroDParameters::xml_element_name_ = "Couple_to_svZeroD";
+
+CoupleSvZeroDParameters::CoupleSvZeroDParameters()
+{
+ // A parameter that must be defined.
+ bool required = true;
+
+ type = Parameter<std::string>("type", "", required);
+};
+
+void CoupleSvZeroDParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+ std::string error_msg = "Unknown Couple_to_svZeroD type=TYPE XML element '";
+
+ // Get the 'type' from the <Couple_to_genBC type=TYPE> element.
+ const char* stype;
+ auto result = xml_elem->QueryStringAttribute("type", &stype);
+ if (stype == nullptr) {
+ throw std::runtime_error("No TYPE given in the XML <Couple_to_svZeroD type=TYPE> element.");
+ }
+ type.set(std::string(stype));
+ auto item = xml_elem->FirstChildElement();
+
+ using std::placeholders::_1;
+ using std::placeholders::_2;
+ std::function<void(const std::string&, const std::string&)> ftpr = 
+ std::bind( &CoupleSvZeroDParameters::set_parameter_value, *this, _1, _2);
+
+ xml_util_set_parameters(ftpr, xml_elem, error_msg);
+
+ value_set = true;
+}
+
+
 //////////////////////////////////////////////////////////
 // OutputParameters //
 //////////////////////////////////////////////////////////
@@ -1701,6 +1743,9 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem)
  } else if (name == CoupleGenBCParameters::xml_element_name_) {
  couple_to_genBC.set_values(item);
 
+ } else if (name == CoupleSvZeroDParameters::xml_element_name_) {
+ couple_to_svZeroD.set_values(item);
+
  } else if (name == DomainParameters::xml_element_name_) {
  auto domain_params = new DomainParameters();
  domain_params->set_values(item);

diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h
@@ -583,6 +583,26 @@ class CoupleGenBCParameters : public ParameterLists
  bool value_set = false;
 };
 
+//-----------------------
+// CoupleSvZeroDParameters
+//-----------------------
+// Coupling to svZeroD.
+//
+class CoupleSvZeroDParameters : public ParameterLists
+{
+ public:
+ CoupleSvZeroDParameters();
+
+ static const std::string xml_element_name_;
+
+ bool defined() const { return value_set; };
+ void set_values(tinyxml2::XMLElement* xml_elem);
+
+ // attributes.
+ Parameter<std::string> type;
+
+ bool value_set = false;
+};
 /// @brief Body force over a mesh using the "Add_BF" command.
 ///
 /// \code {.xml}
@@ -1165,6 +1185,7 @@ class EquationParameters : public ParameterLists
 
  CoupleCplBCParameters couple_to_cplBC;
  CoupleGenBCParameters couple_to_genBC;
+ CoupleSvZeroDParameters couple_to_svZeroD;
 
  DomainParameters* default_domain = nullptr;
 

diff --git a/Code/Source/svFSI/baf_ini.cpp b/Code/Source/svFSI/baf_ini.cpp
@@ -38,6 +38,7 @@
 #include "nn.h"
 #include "set_bc.h"
 #include "utils.h"
+#include "svZeroD_subroutines.h"
 
 #include "fsils_api.hpp"
 #include "fils_struct.hpp"
@@ -166,6 +167,10 @@ void baf_ini(Simulation* simulation)
  set_bc::genBC_Integ_X(com_mod, cm_mod, "I");
  }
 
+ if (com_mod.cplBC.useSvZeroD) {
+ svZeroD::init_svZeroD(com_mod, cm_mod);
+ }
+
  if (com_mod.cplBC.schm != CplBCType::cplBC_E) {
  set_bc::calc_der_cpl_bc(com_mod, cm_mod);
  }

diff --git a/Code/Source/svFSI/distribute.cpp b/Code/Source/svFSI/distribute.cpp
@@ -526,13 +526,20 @@ void distribute(Simulation* simulation)
  cm.bcast(cm_mod, &cplBC.nFa);
  cm.bcast_enum(cm_mod, &cplBC.schm);
  cm.bcast(cm_mod, &cplBC.useGenBC);
+ cm.bcast(cm_mod, &cplBC.useSvZeroD);
 
  if (cplBC.useGenBC) { 
  if (cm.slv(cm_mod)) { 
  cplBC.nX = 0;
  cplBC.xo.resize(cplBC.nX);
  }
 
+ } else if (cplBC.useSvZeroD) { 
+ if (cm.slv(cm_mod)) { 
+ cplBC.nX = 0;
+ cplBC.xo.resize(cplBC.nX);
+ }
+
  } else { 
  cm.bcast(cm_mod, &cplBC.nX);
  if (cplBC.xo.size() == 0) {

diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp
@@ -1366,11 +1366,14 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq)
  if (eq_params->couple_to_genBC.defined()) {
  cplBC.useGenBC = true;
  cplbc_type_str = eq_params->couple_to_genBC.type.value();
+ } else if (eq_params->couple_to_svZeroD.defined()) {
+ cplBC.useSvZeroD = true;
+ cplbc_type_str = eq_params->couple_to_svZeroD.type.value();
  } else if (eq_params->couple_to_cplBC.defined()) {
  cplbc_type_str = eq_params->couple_to_cplBC.type.value();
  }
 
- if (eq_params->couple_to_genBC.defined() || eq_params->couple_to_cplBC.defined()) {
+ if (eq_params->couple_to_genBC.defined() || eq_params->couple_to_cplBC.defined() || eq_params->couple_to_svZeroD.defined()) {
  try {
  cplBC.schm = consts::cplbc_name_to_type.at(cplbc_type_str);
  } catch (const std::out_of_range& exception) {
@@ -1384,6 +1387,9 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq)
  cplBC.binPath = eq_params->couple_to_genBC.zerod_code_file_path.value();
  cplBC.commuName = "GenBC.int";
  cplBC.nX = 0;
+ } else if (cplBC.useSvZeroD) {
+ cplBC.commuName = "svZeroD_interface.dat";
+ cplBC.nX = 0;
  } else {
  auto& cplBC_params = eq_params->couple_to_cplBC;
  cplBC.nX = cplBC_params.number_of_unknowns.value();

diff --git a/Code/Source/svFSI/read_files.h b/Code/Source/svFSI/read_files.h
@@ -46,7 +46,7 @@ namespace read_files_ns {
  using EquationNdop = std::array<int, 4>;
  using EquationOutputs = std::array<consts::OutputType, maxOutput>;
  using EquationPhys = std::vector<consts::EquationType>;
- using EquationProps = std::array<std::array<consts::PhysicalProperyType, 10>, consts::maxNProp>;
+ using EquationProps = std::array<std::array<consts::PhysicalProperyType, consts::maxNProp>, 20>;
 
  void face_match(ComMod& com_mod, faceType& lFa, faceType& gFa, Vector<int>& ptr);
 

diff --git a/Code/Source/svFSI/set_bc.cpp b/Code/Source/svFSI/set_bc.cpp
@@ -43,6 +43,7 @@
 #include "ustruct.h"
 #include "utils.h"
 #include <math.h>
+#include "svZeroD_subroutines.h"
 
 namespace set_bc {
 
@@ -63,7 +64,9 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod)
  #endif
 
  const int iEq = 0;
- const double absTol = 1.0e-8;
+ // NOTE: For coupling with svZeroDPlus, absTol needs to be > 1e-8 to be compatible with the default convergence tolerance of svZeroDPlus (1e-8)
+ // If this is not true, the finite difference computation of bc.r below results in zero because the perturbation is below the svZeroDPlus tolerance 
+ const double absTol = 1.0e-7;
  const double relTol = 1.0e-5;
 
  int nsd = com_mod.nsd;
@@ -168,6 +171,8 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod)
  // Call genBC or cplBC to get updated pressures or flowrates.
  if (cplBC.useGenBC) {
  set_bc::genBC_Integ_X(com_mod, cm_mod, "D");
+ } else if (cplBC.useSvZeroD) {
+ svZeroD::calc_svZeroD(com_mod, cm_mod, 'D');
  } else {
  set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag);
  }
@@ -213,9 +218,11 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod)
 
  // Call genBC or cplBC again with perturbed flowrate
  if (cplBC.useGenBC) {
- set_bc::genBC_Integ_X(com_mod, cm_mod, "D");
- } else {
- set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag);
+ set_bc::genBC_Integ_X(com_mod, cm_mod, "D");
+ } else if (cplBC.useSvZeroD) {
+ svZeroD::calc_svZeroD(com_mod, cm_mod, 'D');
+ } else {
+ set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag);
  }
 
  // Finite difference calculation of the resistance dP/dQ
@@ -752,6 +759,8 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod)
  // Updates pressure or flowrates stored in cplBC.fa[i].y
  if (cplBC.useGenBC) {
  set_bc::genBC_Integ_X(com_mod, cm_mod, "D");
+ } else if (cplBC.useSvZeroD){
+ svZeroD::calc_svZeroD(com_mod, cm_mod, 'D');
  } else {
  set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag);
  }