Skip to content

Commit

Permalink
struct/ustruct 0D coupling (SimVascular#188)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@aabrown100-git
aabrown100-git authored Mar 7, 2024
1 parent 308791b commit 29a94c3
Show file tree
Hide file tree
Showing 90 changed files with 5,404 additions and 122 deletions.
3 changes: 2 additions & 1 deletion .github/workflows/test.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 8 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,3 +20,11 @@ __pycache__

# Doxygen output
Documentation/html/

# genBC files
**/genBC/obj/
**/genBC_svFSIplus/obj/
genBC.exe
AllData
GenBC.int
InitialData
8 changes: 5 additions & 3 deletions Code/Source/svFSI/ComMod.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
137 changes: 107 additions & 30 deletions Code/Source/svFSI/all_fun.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
//--------------
Expand Down Expand Up @@ -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)
Expand All @@ -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);
}
}

Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand All @@ -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);
}
}

Expand Down Expand Up @@ -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) {
Expand All @@ -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);
Expand All @@ -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;
}
Expand All @@ -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;

Expand All @@ -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();
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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;
}
Expand All @@ -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;

Expand All @@ -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 {
Expand All @@ -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");
}
Expand Down
Loading

0 comments on commit 29a94c3

Please sign in to comment.