struct/ustruct 0D coupling (SimVascular#188)

* struct 0D coupling changes and test, but still some bugs

* Bug fix, cleanup, and tests

Bug fix: In set_bc.cpp, when computing flowrates with integ, I was calling the wrong overloaded integ function, so the integration was not being performed in the 'o' or 'n' configurations. Instead, the flowrates were being computed with the 'r'eference configuration. The tangent contribution, however, correctly accounted for the 'n' configuration, so it was inconsistent with the flowrate. This lead to poor convergence after significant deformation of the LV model in the test case.

Cleanup: Removing some unnecessary lines and function arguments. Also, adding catch for invalid configuration character in gnnb(), and initializing v to 0 in fsils_bc_update()

Tests: 2 tests, LV_NeoHookean_passive_genBC and LV_NeoHookean_passive. LV_NeoHookean_passive_genBC is the test for struct - 0D coupling feature, but I figured I'd add LV_NeoHookean_passive (no genBC coupling) anyway. With this commit, there is very good (not exact) agreement of results between vvedula22/svFSI and aabrown100-git/svFSIplus for LV_NeoHookean_passive_genBC. The convergence histories, AllData files, and results.vtu are all nearly identical. There is also good agreement for LV_NeoHookean_passive.

* Extremely minor edit to LV_NeoHookean_passive input file

* Add folwP check, finish test cases.

* Bug for ustruct 0D coupling, ustruct test

Fixed bug in set_bc.cpp for ustruct 0D coupling. When using ustruct, the Yo or Yn has nsd+1 components (nsd velocities and 1 pressure). When computing flowrate, using Yo or Yn with nsd+1 components caused an error in integ(). Error is avoided by explicitly providing the start (0) and stop (nsd-1) indices to integ(). This tells integ to only use the velocity components in Yo or Yn.

Also, adding LV_NeoHookean_passive_genBC test. This is identical to the test in struct/, except uses ustruct.

* Remove previous results folders when testing

* Add CMM to valid physics for 0D coupling

* Update genBC Makefile to use mpif90, remove system calls in tests

* Add gcc to MacOS dependencies for testing

Adding gcc in hopes that gfortran will be available to compile genBC on GitHub Mac machine (SimVascular#188)

* Add debug statements for Github Mac (will revert)

* Modify debug statements for Github Mac

* Remove debug statements for Github Mac

Keeping brew reinstall -v gcc in test.yml, which seemingly properly installs gfortran and allows genBC to be compiled in my test cases.

* Addressing style comments

Addressing style comments from Matteo and Dave in PR: SimVascular#188

* Removing LV_Holzapfel_passive svFSI.inp

* Add all output fields to test cases and update ref sols

* Change enum MechConfigType names

.github/workflows/test.yml

@@ -20,7 +20,8 @@ jobs:
  - name: Install MacOS dependencies
  if: startsWith(matrix.os, 'macos')
  run: | 
- brew install cmake vtk openblas lapack mesa open-mpi qt
+ brew reinstall -v gcc
+ brew install -v cmake vtk openblas lapack mesa open-mpi qt
  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

.gitignore

@@ -20,3 +20,11 @@ __pycache__
 
 # Doxygen output
 Documentation/html/
+
+# genBC files
+**/genBC/obj/
+**/genBC_svFSIplus/obj/
+genBC.exe
+AllData
+GenBC.int
+InitialData

Code/Source/svFSI/ComMod.h

@@ -1392,7 +1392,8 @@ class ComMod {
  /// @brief Current equation degrees of freedom
  int dof = 0;
 
- /// @brief Global total number of nodes
+ /// @brief Global total number of nodes, across all meshes (total) and all 
+ /// procs (global)
  int gtnNo = 0;
 
  /// @brief Number of equations
@@ -1431,7 +1432,8 @@ class ComMod {
  /// @brief Total number of degrees of freedom per node
  int tDof = 0;
 
- /// @brief Total number of nodes
+ /// @brief Total number of nodes (number of nodes on current proc across
+ /// all meshes)
  int tnNo = 0;
 
  /// @brief Restart Time Step
@@ -1509,7 +1511,7 @@ class ComMod {
  /// @brief LHS matrix
  Array<double> Val;
 
- /// @brief Position vector
+ /// @brief Position vector of mesh nodes (in ref config)
  Array<double> x;
 
  /// @brief Old variables (velocity)

Code/Source/svFSI/all_fun.cpp

@@ -34,12 +34,15 @@
 #include "mat_fun.h"
 #include "nn.h"
 #include "utils.h"
+#include "consts.h"
 
 #include <bitset>
 #include <math.h>
 
 namespace all_fun {
 
+ using namespace consts;
+
 //--------------
 // aspect_ratio
 //--------------
@@ -292,10 +295,14 @@ global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Arra
  return result;
 }
 
-/// @brief This routine integrate an equation over a particular domain
+/// @brief This routine integrated a scalar field over a particular domain.
 ///
 /// Note that 'l' and 'u' should be 0-based and are used to index into 's'.
-//
+/// @param dId domain id
+/// @param s an array containing a scalar value for each node in the mesh
+/// @param l lower index of s
+/// @param u upper index of s (must be equal to l)
+/// @param pFlag flag for using Taylor-Hood function space for pressure
 /// Replicates 'FUNCTION vInteg(dId, s, l, u, pFlag)' defined in ALLFUN.f.
 //
 double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array<double>& s, int l, int u, bool pFlag)
@@ -319,10 +326,16 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array<do
  if (nNo != tnNo) {
  if (ibFlag) {
  if (nNo != com_mod.ib.tnNo) {
- throw std::runtime_error("Incompatible vector size in vInteg");
+ std::string msg = "Incompatible vector size in vInteg in domain: ";
+ msg += std::to_string(dId);
+ msg += "\nNumber of nodes in s must be equal to total number of nodes (immersed boundary).\n";
+ throw std::runtime_error(msg);
  }
  } else { 
- throw std::runtime_error("Incompatible vector size in vInteg");
+ std::string msg = "Incompatible vector size in vInteg in domain: ";
+ msg += std::to_string(dId);
+ msg += "\nNumber of nodes in s must be equal to total number of nodes.\n";
+ throw std::runtime_error(msg);
  } 
  }
 
@@ -535,11 +548,16 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array<do
  return result;
 }
 
-/// @brief This routine integrate s over the surface faId.
+/// @brief This routine integrate a scalar field s over the face lFa.
 ///
 /// Reproduces 'FUNCTION IntegS(lFa, s, pflag)'.
+///
+/// @param lFa face type, representing a face on the computational mesh
+/// @param s an array containing a scalar value for each node in the mesh
+/// @param pFlag flag for using Taylor-Hood function space for pressure
+/// @param cfg denotes which mechanical configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference.
 //
-double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector<double>& s, bool pFlag)
+double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector<double>& s, bool pFlag, MechanicalConfigurationType cfg)
 {
  using namespace consts;
  #define n_debug_integ_s
@@ -564,7 +582,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  insd = 0;
  }
 
- int nNo = s.size();
+ int nNo = s.size(); // Total number of nodes on a processor
  #ifdef debug_integ_s
  dmsg << "nNo: " << nNo;
  dmsg << "insd: " << insd;
@@ -574,10 +592,16 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  if (nNo != com_mod.tnNo) {
  if (com_mod.ibFlag) {
  if (nNo != com_mod.ib.tnNo) {
- throw std::runtime_error("Incompatible vector size in Integ");
+ std::string msg = "Incompatible vector size in integS on face: ";
+ msg += lFa.name;
+ msg += "\nNumber of nodes in s must be equal to total number of nodes (immersed boundary).\n";
+ throw std::runtime_error(msg);
  }
  } else {
- throw std::runtime_error("Incompatible vector size in vInteg");
+ std::string msg = "Incompatible vector size in integS on face: ";
+ msg += lFa.name;
+ msg += "\nNumber of nodes in s must be equal to total number of nodes.\n";
+ throw std::runtime_error(msg);
  }
  }
 
@@ -634,8 +658,11 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  dmsg << "fs.eType: " << fs.eType;
  dmsg << "fs.w: " << fs.w;
  #endif
+
+ // Initialize integral to 0
  double result = 0.0;
 
+ // Loop over elements on face
  for (int e = 0; e < lFa.nEl; e++) {
  // [TODO:DaveP] not implemented.
  if (lFa.eType == ElementType::NRB) {
@@ -649,16 +676,19 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  fs.Nx = lFa.Nx;
  }
 
+ // Loop over the Gauss points
  for (int g = 0; g < fs.nG; g++) {
  Vector<double> n(nsd);
  if (!isIB) {
+ // Get normal vector in cfg configuration
  auto Nx = fs.Nx.slice(g);
- nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n);
+ nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, cfg);
  }
 
+ // Calculating the Jacobian (encodes area of face element)
  double Jac = sqrt(utils::norm(n));
 
- // Calculating the function value
+ // Calculating the function value at Gauss point
  double sHat = 0.0;
  for (int a = 0; a < fs.eNoN; a++) {
  int Ac = lFa.IEN(a,e);
@@ -670,6 +700,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  }
  }
 
+ // If using multiple processors, add result from all processors
  if (com_mod.cm.seq() || isIB) {
  return result; 
  }
@@ -678,11 +709,19 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  return result; 
 }
 
-/// @brief This routine integrate s over the surface faId. 
+/// @brief This routine integrates vector field s dotted with the face normal n
+/// over the face lFa. For example, if s contains the velocity at each node on 
+/// the face, this function computed the velocity flux through the face.
 ///
 /// Reproduces 'FUNCTION IntegV(lFa, s)'
+///
+/// @param lFa face type, representing a face on the computational mesh
+/// @param s an array containing a vector value for each node in the mesh
+/// @param pFlag flag for using Taylor-Hood function space for pressure
+/// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference.
 //
-double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array<double>& s)
+double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, 
+ const Array<double>& s, MechanicalConfigurationType cfg)
 {
  using namespace consts;
 
@@ -699,8 +738,16 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  int insd = nsd - 1;
  int tnNo = com_mod.tnNo;
 
+ #ifdef debug_integ_V
+ dmsg << "s.nrows(): " << s.nrows();
+ dmsg << "nsd: " << nsd;
+ #endif
+
  if (s.nrows() != nsd) {
- throw std::runtime_error("Incompatible vector size in integ");
+ std::string msg = "Incompatible vector size in integV on face: "; 
+ msg += lFa.name;
+ msg += "\nNumber of rows in s must be equal to number of spatial dimensions.\n";
+ throw std::runtime_error(msg);
  }
 
  int nNo = s.ncols();
@@ -712,22 +759,31 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  if (nNo != tnNo) {
  if (com_mod.ibFlag) {
  if (nNo != com_mod.ib.tnNo) {
- throw std::runtime_error("Incompatible vector size in integ");
+ std::string msg = "Incompatible vector size in integV on face: ";
+ msg += lFa.name;
+ msg += "\nNumber of nodes in s must be equal to total number of nodes (immersed boundary).\n";
+ throw std::runtime_error(msg);
  }
  } else {
- throw std::runtime_error("Incompatible vector size in integ");
+ std::string msg = "Incompatible vector size in integV on face: ";
+ msg += lFa.name; 
+ msg += "\nNumber of nodes in s must be equal to total number of nodes.\n";
+ throw std::runtime_error(msg);
  }
  }
 
+ // If using Immersed Boundary Method
  bool isIB = false;
  if (com_mod.ibFlag) {
  if (nNo == com_mod.ib.tnNo) {
  isIB = true;
  }
  }
 
+ // Initialize integral to 0
  double result = 0.0;
 
+ // Loop over elements on face
  for (int e = 0; e < lFa.nEl; e++) {
  //dmsg << "----- e " << e+1 << " -----";
  // Updating the shape functions, if this is a NURB
@@ -739,23 +795,25 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  }
  }
 
+ // Loop over the Gauss points
  for (int g = 0; g < lFa.nG; g++) {
  //dmsg << ">>> g: " << g+1;
  Vector<double> n(nsd);
  if (!isIB) {
+ // Get normal vector in cfg configuration
  auto Nx = lFa.Nx.slice(g);
- nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n);
+ nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, cfg);
  //CALL GNNB(lFa, e, g, nsd-1, lFa.eNoN, lFa.Nx(:,:,g), n)
  } else {
  //CALL GNNIB(lFa, e, g, n)
  }
 
- // Calculating the function value
- //
+ // Calculating the function value (s dot n)dA at this Gauss point
  double sHat = 0.0;
 
  for (int a = 0; a < lFa.eNoN; a++) {
  int Ac = lFa.IEN(a,e);
+ // Compute s dot n
  for (int i = 0; i < nsd; i++) {
  sHat = sHat + lFa.N(a,g) * s(i,Ac) * n(i);
  //dmsg << "s(i,Ac): " << s(i,Ac);
@@ -769,6 +827,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  }
  }
 
+ // If using multiple processors, add result from all processors
  if (cm.seq() || isIB) {
  return result; 
  }
@@ -778,14 +837,25 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
  return result; 
 }
 
-/// @brief This routine integrate s(l:u,:) over the surface faId.
+/// @brief This routine integrate s(l:u,:) over the surface faId, where s is an
+/// array of scalars or an array of nsd-vectors. This routine calls overloaded
+/// functions to integrate scalars, if s is scalar (i.e. l=u), or vectors if s 
+/// is vector (i.e. l<u).
 ///
 /// Note that 'l' and 'u' should be 0-based and are used to index into 's'.
 ///
 /// Reproduces 'FUNCTION IntegG(lFa, s, l, uo, THflag)'.
+///
+/// @param lFa face type, representing a face on the computational mesh.
+/// @param s an array containing a vector value for each node in the mesh.
+/// @param l lower index of s
+/// @param uo optional: upper index of s. Default u = l.
+/// @param THlag flag for using Taylor-Hood function space for pressure.
+/// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference.
+///
 //
 double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, 
- const Array<double>& s, const int l, std::optional<int> uo, bool THflag)
+ const Array<double>& s, const int l, std::optional<int> uo, bool THflag, MechanicalConfigurationType cfg)
 {
  using namespace consts;
 
@@ -803,9 +873,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa,
  int insd = nsd - 1;
  int tnNo = com_mod.tnNo;
 
- // Set u if uo is given.
+ // Set u if uo is given. Else, set u = l.
  int u{0};
-
  if (uo.has_value()) { 
  u = uo.value();
  } else {
@@ -818,30 +887,38 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa,
  if (nNo != tnNo) {
  if (com_mod.ibFlag) {
  if (nNo != com_mod.tnNo) { 
- throw std::runtime_error("Incompatible vector size in integ");
+ std::string msg = "Incompatible vector size in integG on face: ";
+ msg += lFa.name; 
+ msg += "\nNumber of nodes in s must be equal to total number of nodes (immersed boundary).\n";
+ throw std::runtime_error(msg);
  }
  } else {
- throw std::runtime_error("Incompatible vector size in integ");
+ std::string msg = "Incompatible vector size in integG on face: ";
+ msg += lFa.name;
+ msg += "\nNumber of nodes in s must be equal to total number of nodes.\n";
+ throw std::runtime_error(msg);
  }
  }
 
+ // Initialize integral to 0 (not strictly necessary to initialize to 0)
  double result = 0.0; 
-
- if (u-l+1 == nsd) {
+
+ // If s vector, integrate as vector (dot with surface normal)
+ if (u-l+1 == nsd) { 
  Array<double> vec(nsd,nNo);
  for (int a = 0; a < nNo; a++) {
  for (int i = l, n = 0; i <= u; i++, n++) {
  vec(n,a) = s(i,a); 
  }
  }
- result = integ(com_mod, cm_mod, lFa, vec);
-
+ result = integ(com_mod, cm_mod, lFa, vec, cfg);
+ // If s scalar, integrate as scalar
  } else if (l == u) {
  Vector<double> sclr(nNo);
  for (int a = 0; a < nNo; a++) {
  sclr(a) = s(l,a);
  }
- result = integ(com_mod, cm_mod, lFa, sclr, flag);
+ result = integ(com_mod, cm_mod, lFa, sclr, flag, cfg);
  } else {
  throw std::runtime_error("Unexpected dof in integ");
  }