From 29a94c3b8641118d66fc56357ef212a1cce0c63d Mon Sep 17 00:00:00 2001 From: aabrown100-git <71424733+aabrown100-git@users.noreply.github.com> Date: Wed, 6 Mar 2024 21:32:06 -0500 Subject: [PATCH] struct/ustruct 0D coupling (#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 (https://github.com/SimVascular/svFSIplus/pull/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: https://github.com/SimVascular/svFSIplus/pull/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 | 3 +- .gitignore | 8 + Code/Source/svFSI/ComMod.h | 8 +- Code/Source/svFSI/all_fun.cpp | 137 +++- Code/Source/svFSI/all_fun.h | 6 +- Code/Source/svFSI/baf_ini.cpp | 10 +- Code/Source/svFSI/consts.h | 7 + Code/Source/svFSI/eq_assem.cpp | 90 ++- Code/Source/svFSI/eq_assem.h | 4 +- Code/Source/svFSI/fluid.cpp | 8 +- Code/Source/svFSI/lhsa.cpp | 7 +- Code/Source/svFSI/ls.cpp | 1 + Code/Source/svFSI/main.cpp | 5 +- Code/Source/svFSI/nn.cpp | 41 +- Code/Source/svFSI/nn.h | 3 +- Code/Source/svFSI/pic.cpp | 6 + Code/Source/svFSI/set_bc.cpp | 230 ++++-- Code/Source/svFSI/sv_struct.cpp | 16 +- Code/Source/svFSI/ustruct.cpp | 13 +- Code/Source/svFSI/vtk_xml.cpp | 1 + Code/Source/svFSILS/add_bc_mul.cpp | 28 +- Code/Source/svFSILS/bc.cpp | 57 +- Code/Source/svFSILS/fsils_api.hpp | 2 + .../struct/LV_Holzapfel_passive/svFSI.xml | 16 +- .../struct/LV_NeoHookean_passive/README.md | 2 + .../LV_NeoHookean_passive/endo_pressure.dat | 3 + .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../LV_NeoHookean_passive/result_005.vtu | 3 + .../struct/LV_NeoHookean_passive/svFSI.xml | 100 +++ .../LV_NeoHookean_passive_genBC/README.md | 29 + .../V3D_vs_V0D.png | 3 + .../dVdt3D_vs_dVdt0D.png | 3 + .../genBC/Makefile | 79 ++ .../genBC/Makefile.in | 165 +++++ .../genBC/include/initial_values_final.f | 4 + .../genBC/include/parameters_final.f | 15 + .../genBC/src/GenBC.f | 230 ++++++ .../genBC/src/Modules.f | 17 + .../genBC/src/USER.f | 143 ++++ .../genBC_svFSIplus/Makefile | 79 ++ .../genBC_svFSIplus/Makefile.in | 165 +++++ .../include/initial_values_final.f | 4 + .../include/parameters_final.f | 15 + .../genBC_svFSIplus/src/GenBC.f | 232 ++++++ .../genBC_svFSIplus/src/Modules.f | 17 + .../genBC_svFSIplus/src/USER.f | 143 ++++ .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../post_processing_functions.py | 700 ++++++++++++++++++ .../process_results.py | 154 ++++ .../LV_NeoHookean_passive_genBC/pv_plot.png | 3 + .../result_003.vtu | 3 + .../LV_NeoHookean_passive_genBC/svFSI.xml | 103 +++ .../LV_NeoHookean_passive_genBC/README.md | 29 + .../V3D_vs_V0D.png | 3 + .../dVdt3D_vs_dVdt0D.png | 3 + .../genBC/Makefile | 79 ++ .../genBC/Makefile.in | 165 +++++ .../genBC/include/initial_values_final.f | 4 + .../genBC/include/parameters_final.f | 15 + .../genBC/src/GenBC.f | 230 ++++++ .../genBC/src/Modules.f | 17 + .../genBC/src/USER.f | 143 ++++ .../genBC_svFSIplus/Makefile | 79 ++ .../genBC_svFSIplus/Makefile.in | 165 +++++ .../include/initial_values_final.f | 4 + .../include/parameters_final.f | 15 + .../genBC_svFSIplus/src/GenBC.f | 232 ++++++ .../genBC_svFSIplus/src/Modules.f | 17 + .../genBC_svFSIplus/src/USER.f | 143 ++++ .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/endo.vtp | 3 + .../mesh/mesh-surfaces/epi.vtp | 3 + .../mesh/mesh-surfaces/top.vtp | 3 + .../post_processing_functions.py | 700 ++++++++++++++++++ .../process_results.py | 154 ++++ .../LV_NeoHookean_passive_genBC/pv_plot.png | 3 + .../result_003.vtu | 3 + .../LV_NeoHookean_passive_genBC/svFSI.xml | 104 +++ tests/conftest.py | 7 + tests/test_struct.py | 26 +- tests/test_ustruct.py | 20 + 90 files changed, 5404 insertions(+), 122 deletions(-) create mode 100644 tests/cases/struct/LV_NeoHookean_passive/README.md create mode 100755 tests/cases/struct/LV_NeoHookean_passive/endo_pressure.dat create mode 100644 tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive/result_005.vtu create mode 100644 tests/cases/struct/LV_NeoHookean_passive/svFSI.xml create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/README.md create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile.in create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/USER.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/post_processing_functions.py create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/process_results.py create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/pv_plot.png create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/result_003.vtu create mode 100644 tests/cases/struct/LV_NeoHookean_passive_genBC/svFSI.xml create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/README.md create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile.in create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/USER.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/post_processing_functions.py create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/process_results.py create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/pv_plot.png create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/result_003.vtu create mode 100644 tests/cases/ustruct/LV_NeoHookean_passive_genBC/svFSI.xml diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 02726edb..1a92e608 100644 --- a/.github/workflows/test.yml +++ b/.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 diff --git a/.gitignore b/.gitignore index d90f540e..24ad342d 100644 --- a/.gitignore +++ b/.gitignore @@ -20,3 +20,11 @@ __pycache__ # Doxygen output Documentation/html/ + +# genBC files +**/genBC/obj/ +**/genBC_svFSIplus/obj/ +genBC.exe +AllData +GenBC.int +InitialData diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index ef4d5c3b..2123299a 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/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 Val; - /// @brief Position vector + /// @brief Position vector of mesh nodes (in ref config) Array x; /// @brief Old variables (velocity) diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index 94a02979..1285894c 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -34,12 +34,15 @@ #include "mat_fun.h" #include "nn.h" #include "utils.h" +#include "consts.h" #include #include 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& 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& s, bool pFlag) +double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector& 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 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& s) +double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, + const Array& 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,13 +759,20 @@ 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) { @@ -726,8 +780,10 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } } + // 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 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& s, const int l, std::optional uo, bool THflag) + const Array& s, const int l, std::optional 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 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 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"); } diff --git a/Code/Source/svFSI/all_fun.h b/Code/Source/svFSI/all_fun.h index 8f146ee4..89377920 100644 --- a/Code/Source/svFSI/all_fun.h +++ b/Code/Source/svFSI/all_fun.h @@ -58,12 +58,12 @@ namespace all_fun { bool pFlag=false); double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector& s, - bool pFlag=false); + bool pFlag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, - const int l, std::optional uo=std::nullopt, bool THflag=false); + const int l, std::optional uo=std::nullopt, bool THflag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s); + double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys); diff --git a/Code/Source/svFSI/baf_ini.cpp b/Code/Source/svFSI/baf_ini.cpp index b4fa389f..12e215a4 100644 --- a/Code/Source/svFSI/baf_ini.cpp +++ b/Code/Source/svFSI/baf_ini.cpp @@ -121,7 +121,8 @@ void baf_ini(Simulation* simulation) int iEq = 0; com_mod.cplBC.fa.resize(com_mod.cplBC.nFa); com_mod.cplBC.xn.resize(com_mod.cplBC.nX); - + + // Assign cplBC internal variables if (com_mod.cplBC.coupled) { auto& eq = com_mod.eq[iEq]; for (int iBc = 0; iBc < eq.nBc; iBc++) { @@ -139,6 +140,8 @@ void baf_ini(Simulation* simulation) } else if (utils::btest(bc.bType, iBC_Neu)) { com_mod.cplBC.fa[i].bGrp = CplBCType::cplBC_Neu; + // For implicit or semi-implicit (not explicit) Neumann 0D coupling scheme, + // set bType to resistance if (com_mod.cplBC.schm != CplBCType::cplBC_E) { bc.bType= utils::ibset(bc.bType, iBC_res); } @@ -483,6 +486,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) dmsg << "Flag: " << flag; #endif + // Compute integral of normal vector over surface element if (!flag) { Vector nV(nsd); for (int e = 0; e < lFa.nEl; e++) { @@ -702,6 +706,7 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp Array sV(nsd,tnNo); Vector gNodes(nNo); + // Copy mesh node id corresponding to face node id to gNodes for (int a= 0; a < nNo; a++) { gNodes(a) = lFa.gN(a); } @@ -733,6 +738,7 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp } } else if (btest(lBc.bType, iBC_Neu)) { + // Compute integral of normal vector over the face (needed for resistance BC/0D-coupling) if (btest(lBc.bType, iBC_res)) { sV = 0.0; for (int e = 0; e < lFa.nEl; e++) { @@ -762,6 +768,8 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp lsPtr = lsPtr + 1; lBc.lsPtr = lsPtr; + + // Fills lhs.face(i) variables, including val is sVl exists fsils_bc_create(com_mod.lhs, lsPtr, lFa.nNo, nsd, BcType::BC_TYPE_Neu, gNodes, sVl); } else { lBc.lsPtr = -1; diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h index 5ff9ffb5..1e978648 100644 --- a/Code/Source/svFSI/consts.h +++ b/Code/Source/svFSI/consts.h @@ -460,6 +460,13 @@ std::ostream& operator<<(typename std::enable_if::value, std::os return stream << static_cast::type>(e); } +//// Mechanical configurations +enum class MechanicalConfigurationType +{ + reference, // reference configuration + old_timestep, // old timestep (n) configuration + new_timestep // new timestep (n+1) configuration +}; }; #endif diff --git a/Code/Source/svFSI/eq_assem.cpp b/Code/Source/svFSI/eq_assem.cpp index 7da5096c..fe274e30 100644 --- a/Code/Source/svFSI/eq_assem.cpp +++ b/Code/Source/svFSI/eq_assem.cpp @@ -49,6 +49,8 @@ #include "sv_struct.h" #include "ustruct.h" +#include + #include #ifdef WITH_TRILINOS @@ -184,17 +186,22 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& } -/// @brief For struct/ustruct - construct follower pressure load. -/// +/// @brief For struct/ustruct - construct follower pressure load contribution +/// to the residual vector and stiffness matrix. /// We use Nanson's formula to take change in normal direction with /// deformation into account. Additional calculations based on mesh /// need to be performed. /// /// Reproduces 'SUBROUTINE BNEUFOLWP(lFa, hg, Dg)' -// -void b_neu_folw_p(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Dg) +/// @param com_mod +/// @param lBc +/// @param lFa +/// @param hg Pressure magnitude +/// @param Dg +void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg) { using namespace consts; + using namespace utils; #define n_debug_b_neu_folw_p #ifdef debug_b_neu_folw_p @@ -277,6 +284,7 @@ void b_neu_folw_p(ComMod& com_mod, const faceType& lFa, const Vector& hg nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, ksix); } + // Get surface normal vector Vector nV(nsd); auto Nx_g = lFa.Nx.slice(g); nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx_g, nV); @@ -284,6 +292,7 @@ void b_neu_folw_p(ComMod& com_mod, const faceType& lFa, const Vector& hg nV = nV / Jac; double w = lFa.w(g)*Jac; + // Compute residual and tangent contributions if (cPhys == EquationType::phys_ustruct) { if (nsd == 3) { ustruct::b_ustruct_3d(com_mod, eNoN, w, N, Nx, dl, hl, nV, lR, lK, lKd); @@ -316,8 +325,81 @@ void b_neu_folw_p(ComMod& com_mod, const faceType& lFa, const Vector& hg } #endif } + + // Now update surface integrals involved in coupled/resistance BC + // contribution to stiffness matrix to reflect deformed geometry. + if (btest(lBc.bType, iBC_res)) { + fsi_ls_upd(com_mod, lBc, lFa); + } } +/// @brief Update the surface integral involved in the coupled/resistance BC +/// contribution to the stiffness matrix to reflect deformed geometry, if using +/// a follower pressure load. +/// The value of this integral is stored in lhs%face%val. +/// This integral is sV = int_Gammat (Na * n_i) (See Brown et al. 2024, Eq. 56) +/// where Na is the shape function and n_i is the normal vector. +/// +/// This function updates the variable lhs%face%val with the new value, which +/// is eventually used in ADDBCMUL() in the linear solver to add the contribution +/// from the resistance BC to the matrix-vector product of the tangent matrix and +/// an arbitrary vector. +void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa) +{ + using namespace consts; + using namespace utils; + using namespace fsi_linear_solver; + + #define n_debug_fsi_ls_upd + #ifdef debug_fsi_ls_upd + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "lFa.name: " << lFa.name; + #endif + + auto& cm = com_mod.cm; + int nsd = com_mod.nsd; + int tnNo = com_mod.tnNo; + + int iM = lFa.iM; + int nNo = lFa.nNo; + + Array sVl(nsd,nNo); + Array sV(nsd,tnNo); + + // Updating the value of the surface integral of the normal vector + // using the deformed configuration ('n' = new = timestep n+1) + sV = 0.0; + for (int e = 0; e < lFa.nEl; e++) { + if (lFa.eType == ElementType::NRB) { + // CALL NRBNNXB(msh(iM),lFa,e) + } + for (int g = 0; g < lFa.nG; g++) { + Vector n(nsd); + auto Nx = lFa.Nx.slice(g); + + auto cfg = MechanicalConfigurationType::new_timestep; + + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, cfg); + // + for (int a = 0; a < lFa.eNoN; a++) { + int Ac = lFa.IEN(a,e); + for (int i = 0; i < nsd; i++) { + sV(i,Ac) = sV(i,Ac) + lFa.N(a,g)*lFa.w(g)*n(i); + } + } + } + } + + if (sVl.size() != 0) { + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN(a); + sVl.set_col(a, sV.col(Ac)); + } + } + // Update lhs.face(i).val with the new value of the surface integral + fsils_bc_update(com_mod.lhs, lBc.lsPtr, lFa.nNo, nsd, sVl); +}; /// @brief This routine assembles the equation on a given mesh. /// diff --git a/Code/Source/svFSI/eq_assem.h b/Code/Source/svFSI/eq_assem.h index 55029260..a1f4a20a 100644 --- a/Code/Source/svFSI/eq_assem.h +++ b/Code/Source/svFSI/eq_assem.h @@ -38,7 +38,9 @@ namespace eq_assem { void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Yg); -void b_neu_folw_p(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Dg); +void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg); + +void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa); void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Yg, const Array& Dg); diff --git a/Code/Source/svFSI/fluid.cpp b/Code/Source/svFSI/fluid.cpp index f0262ddc..bf640f57 100644 --- a/Code/Source/svFSI/fluid.cpp +++ b/Code/Source/svFSI/fluid.cpp @@ -93,7 +93,7 @@ void b_fluid(ComMod& com_mod, const int eNoN, const double w, const Vector& incL, const Vecto #ifdef WITH_TRILINOS + // Set up resistance/coupled-Neumann BC for Trilinos linear solver if (lEq.useTLS) { init_dir_and_coup_neu(com_mod, incL, res); } diff --git a/Code/Source/svFSI/main.cpp b/Code/Source/svFSI/main.cpp index c6d6969a..3a1fca6a 100644 --- a/Code/Source/svFSI/main.cpp +++ b/Code/Source/svFSI/main.cpp @@ -243,7 +243,7 @@ void iterate_solution(Simulation* simulation) set_bc::set_bc_dir(com_mod, An, Yn, Dn); - // Inner loop for iteration + // Inner loop for Newton iteration // int inner_count = 1; int reply; @@ -454,6 +454,7 @@ void iterate_solution(Simulation* simulation) for (int iBc = 0; iBc < eq.nBc; iBc++) { int i = eq.bc[iBc].lsPtr; if (i != -1) { + // Resistance term for coupled Neumann BC tangent contribution res(i) = eq.gam * dt * eq.bc[iBc].r; incL(i) = 1; } @@ -472,7 +473,7 @@ void iterate_solution(Simulation* simulation) com_mod.Val.write("Val_solve"+ istr); com_mod.R.write("R_solve"+ istr); - // Solution is obtained, now updating (Corrector) + // Solution is obtained, now updating (Corrector) and check for convergence // // Modifies: com_mod.An com_mod.Dn com_mod.Yn cep_mod.Xion com_mod.pS0 com_mod.pSa // com_mod.pSn com_mod.cEq eq.FSILS.RI.iNorm eq.pNorm diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp index 6310a685..2c864f27 100644 --- a/Code/Source/svFSI/nn.cpp +++ b/Code/Source/svFSI/nn.cpp @@ -527,14 +527,17 @@ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Arra } } -/// @brief This routine returns a vector at element "e" and Gauss point -/// 'g' of face 'lFa' that is the normal weigthed by Jac, i.e. -/// Jac = SQRT(NORM(n)). +/// @brief This routine returns a surface normal vector at element "e" and Gauss point +/// 'g' of face 'lFa' that is the normal weighted by Jac, i.e. +/// Jac = SQRT(NORM(n)), the Jacobian of the mapping from parent surface element to +/// reference/old/new configuration. +/// +/// cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference /// /// Reproduce Fortran 'GNNB'. // void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n) + const int eNoNb, const Array& Nx, Vector& n, MechanicalConfigurationType cfg) { auto& cm = com_mod.cm; @@ -547,6 +550,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, dmsg << "nsd: " << nsd; dmsg << "insd: " << insd; dmsg << "eNoNb: " << eNoNb; + dmsg << "cfg: " << cfg; #endif int iM = lFa.iM; @@ -606,19 +610,42 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, } } - // Correcting the position vector if mesh is moving + // Correcting the geometry if mesh is moving or if we want the + // area-weighted normal in a different configuration // for (int a = 0; a < eNoN; a++) { int Ac = msh.IEN(a,Ec); for (int i = 0; i < lX.nrows(); i++) { + // Get position vector lX(i,a) = com_mod.x(i,Ac); } - if (com_mod.mvMsh) { for (int i = 0; i < lX.nrows(); i++) { + // Add mesh displacement lX(i,a) = lX(i,a) + com_mod.Do(i+nsd+1,Ac); } } + else { + switch (cfg) { + case MechanicalConfigurationType::reference: + // Do nothing + break; + case MechanicalConfigurationType::old_timestep: + for (int i = 0; i < lX.nrows(); i++) { + // Add displacement at timestep n + lX(i,a) = lX(i,a) + com_mod.Do(i,Ac); + } + break; + case MechanicalConfigurationType::new_timestep: + for (int i = 0; i < lX.nrows(); i++) { + // Add displacement at timestep n+1 + lX(i,a) = lX(i,a) + com_mod.Dn(i,Ac); + } + break; + default: + throw std::runtime_error("gnnb: invalid MechanicalConfigurationType provided"); + } + } } // Calculating surface deflation @@ -699,7 +726,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, } // Changing the sign if neccessary. 'a' locates on the face and 'b' - // outside of the face, in the parent element + // in the interior of the element. v points outward along ba // a = ptr(0); int b = ptr(lFa.eNoN); diff --git a/Code/Source/svFSI/nn.h b/Code/Source/svFSI/nn.h index d48caa94..43172783 100644 --- a/Code/Source/svFSI/nn.h +++ b/Code/Source/svFSI/nn.h @@ -33,6 +33,7 @@ #include "Simulation.h" #include "ComMod.h" +#include "consts.h" namespace nn { @@ -64,7 +65,7 @@ namespace nn { double& Jac, Array& ks); void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n); + const int eNoNb, const Array& Nx, Vector& n, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); void gnns(const int nsd, const int eNoN, const Array& Nxi, Array& xl, Vector& nV, Array& gCov, Array& gCnv); diff --git a/Code/Source/svFSI/pic.cpp b/Code/Source/svFSI/pic.cpp index 69d07e51..de6f6356 100644 --- a/Code/Source/svFSI/pic.cpp +++ b/Code/Source/svFSI/pic.cpp @@ -29,6 +29,11 @@ */ // The code here replicates the Fortran code in PIC.f. +// +// See the publication below, section 4.4 for theory and derivation: +// Bazilevs, et al. "Isogeometric fluid-structure interaction: +// theory, algorithms, and computations.", Computational Mechanics, +// 43 (2008): 3-37. doi: 10.1007/s00466-008-0315-x #include "pic.h" @@ -221,6 +226,7 @@ void picc(Simulation* simulation) // IB treatment //if (ibFlag) CALL IB_PICC() + // Computes norms and check for convergence of Newton iterations double eps = std::numeric_limits::epsilon(); if (utils::is_zero(eq.FSILS.RI.iNorm)) { diff --git a/Code/Source/svFSI/set_bc.cpp b/Code/Source/svFSI/set_bc.cpp index f7872733..acf0f9d1 100644 --- a/Code/Source/svFSI/set_bc.cpp +++ b/Code/Source/svFSI/set_bc.cpp @@ -50,6 +50,12 @@ namespace set_bc { +/// @brief This function calculates updated cplBC pressures or flowrates from 0D, +/// as well as the resistance matrix M ~ dP/dQ from 0D using finite difference. +/// Updates the pressure or flowrates stored in cplBC.fa[i].y and the resistance +/// matrix M ~ dP/dQ stored in eq.bc[iBc].r. +/// @param com_mod +/// @param cm_mod void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) { using namespace consts; @@ -68,6 +74,8 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) auto& eq = com_mod.eq[iEq]; auto& cplBC = com_mod.cplBC; + // If coupling is all with Dirichlet faces, no derivative calculation is needed + // (see Moghadam et al. 2013 Section 2.2.2) if (std::count_if(cplBC.fa.begin(),cplBC.fa.end(),[](cplFaceType& fa){return fa.bGrp == CplBCType::cplBC_Dir;}) == cplBC.fa.size()) { #ifdef debug_calc_der_cpl_bc dmsg << "all cplBC_Dir " << std::endl; @@ -76,8 +84,17 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) } // if (ALL(cplBC.fa.bGrp .EQ. cplBC_Dir)) RETURN + // Determine current physics + auto& cDmn = com_mod.cDmn; + auto cPhys = eq.dmn[cDmn].phys; + + // Mechanical configuration in which to compute flowrate + auto cfg_o = MechanicalConfigurationType::reference; + auto cfg_n = MechanicalConfigurationType::reference; + bool RCRflag = false; + // Loop over BCs for (int iBc = 0; iBc < eq.nBc; iBc++) { #ifdef debug_calc_der_cpl_bc dmsg << "----- iBc " << iBc+1 << " -----" << std::endl; @@ -100,10 +117,29 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) if (ptr != -1) { auto& fa = com_mod.msh[iM].fa[iFa]; - + // Compute flowrates at 3D Neumann boundaries at timesteps n and n+1 if (utils::btest(bc.bType, iBC_Neu)) { - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 0, nsd-1); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, 0, nsd-1); + // If struct or ustruct, use old and new configurations to compute flowrate integral + if ((cPhys == EquationType::phys_struct) || (cPhys == EquationType::phys_ustruct)) { + + // Must use follower pressure load for 0D coupling with struct/ustruct + if (!bc.flwP) { + throw std::runtime_error("[calc_der_cpl_bc] Follower pressure load must be used for 0D coupling with struct/ustruct"); + } + cfg_o = MechanicalConfigurationType::old_timestep; + cfg_n = MechanicalConfigurationType::new_timestep; + } + // If fluid, FSI, or CMM, use reference configuration to compute flowrate integral + // Note that for FSI, mvMsh will modify geometry in gnnb() + else if ((cPhys == EquationType::phys_fluid) || (cPhys == EquationType::phys_FSI) || (cPhys == EquationType::phys_CMM)) { + cfg_o = MechanicalConfigurationType::reference; + cfg_n = MechanicalConfigurationType::reference; + } + else { + throw std::runtime_error("[calc_der_cpl_bc] Invalid physics type for 0D coupling"); + } + cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 0, nsd-1, false, cfg_o); + cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, 0, nsd-1, false, cfg_n); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -112,7 +148,9 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) dmsg << "cplBC.fa[ptr].Qn: " << cplBC.fa[ptr].Qn; #endif - } else if (utils::btest(bc.bType, iBC_Dir)) { + } + // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 + else if (utils::btest(bc.bType, iBC_Dir)) { double area = fa.area; cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd) / area; cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, nsd) / area; @@ -131,12 +169,15 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) dmsg << "RCRflag: " << RCRflag; #endif + // Call genBC or cplBC to get updated pressures or flowrates. if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "D"); } else { set_bc::cplBC_Integ_X(com_mod, cm_mod, RCRflag); } + // Compute the epsilon parameter (diff) for the finite difference calculation + // of the resistance matrix M ~ dP/dQ. Slightly different from Eq. 30 in Moghadam et al. 2013 int j = 0; double diff = 0.0; @@ -144,7 +185,7 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) auto& bc = eq.bc[iBc]; int i = bc.cplBCptr; if (i != -1 && utils::btest(bc.bType, iBC_Neu)) { - diff = diff + (cplBC.fa[i].Qo * cplBC.fa[i].Qo); + diff = diff + (cplBC.fa[i].Qn * cplBC.fa[i].Qn); j = j + 1; } } @@ -156,24 +197,39 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) diff = diff*relTol; } + // Store the original pressures and flowrates + std::vector orgY(cplBC.fa.size()); + std::vector orgQ(cplBC.fa.size()); + + for (size_t i = 0; i < cplBC.fa.size(); i++) { + orgY[i] = cplBC.fa[i].y; + orgQ[i] = cplBC.fa[i].Qn; + } + for (int iBc = 0; iBc < eq.nBc; iBc++) { auto& bc = eq.bc[iBc]; int i = bc.cplBCptr; if (i != -1 && utils::btest(bc.bType, iBC_Neu)) { - double orgY = cplBC.fa[i].y; - double orgQ = cplBC.fa[i].Qn; - cplBC.fa[i].Qn = cplBC.fa[i].Qn + diff; + // Finite difference perturbation in flowrate + cplBC.fa[i].Qn = orgQ[i] + diff; + + // 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); } - bc.r = (cplBC.fa[i].y - orgY) / diff; - cplBC.fa[i].y = orgY; - cplBC.fa[i].Qn = orgQ; + // Finite difference calculation of the resistance dP/dQ + bc.r = (cplBC.fa[i].y - orgY[i]) / diff; + + // Restore the original pressures and flowrates + for (size_t j = 0; j < cplBC.fa.size(); j++) { + cplBC.fa[j].y = orgY[j]; + cplBC.fa[j].Qn = orgQ[j]; +} } } } @@ -259,9 +315,12 @@ void cplBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const bool RCRflag) } } -/// @brief Interface to call 0D code (genBC/gcode) -/// -/// \todo [NOTE] not fully implemented. + +//--------------- +// genBC_Integ_X +//--------------- +// Interface to call 0D code (genBC/gcode) +// // void genBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const std::string& genFlag) { @@ -269,9 +328,17 @@ void genBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const std::string& genF int nDir = 0; int nNeu = 0; + double dt = com_mod.dt; auto& cplBC = com_mod.cplBC; auto& cm = com_mod.cm; + int len = cplBC.binPath.size() + cplBC.commuName.size() + 1; + char command[len]; + strcpy(command, cplBC.binPath.c_str()); + strcat(command, " "); + strcat(command, cplBC.commuName.c_str()); + + // If this process is the master process on the communicator if (cm.mas(cm_mod)) { for (int iFa = 0; iFa < cplBC.nFa; iFa++) { auto& fa = cplBC.fa[iFa]; @@ -283,65 +350,106 @@ void genBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const std::string& genF } } - int fid = 1; - //OPEN(fid, FILE=cplBC.commuName, FORM='UNFORMATTED') - //WRITE(fid) genFlag - //WRITE(fid) dt - //WRITE(fid) nDir - //WRITE(fid) nNeu + // Write coupling info (number of Dirichlet and Neumann surfaces, pressure, + // flow rate) from 3D to cplBC communication file (for GenBC, usually + // called GenBC.int) + int int_size = sizeof(int); + int double_size = sizeof(double); + int flag_size = genFlag.length(); + + std::ofstream genBC_writer; + genBC_writer.open(cplBC.commuName, std::ios::out|std::ios::binary); + if (!genBC_writer.is_open()) { + throw std::runtime_error("Failed to open the genBC initialization file '" + cplBC.commuName + "' to write."); + } + // Flag for how genBC behaves (I: Initializing, T: Iteration loop, L: Last iteration, D: Derivative) + genBC_writer.write( (char*)&flag_size, int_size); + genBC_writer.write(genFlag.c_str(), flag_size); + genBC_writer.write( (char*)&flag_size, int_size); + + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&dt, double_size); + genBC_writer.write( (char*)&double_size, int_size); + + genBC_writer.write( (char*)&int_size, int_size); + genBC_writer.write( (char*)&nDir, int_size); + genBC_writer.write( (char*)&int_size, int_size); + + genBC_writer.write( (char*)&int_size, int_size); + genBC_writer.write( (char*)&nNeu, int_size); + genBC_writer.write( (char*)&int_size, int_size); for (int iFa = 0; iFa < cplBC.nFa; iFa++) { if (cplBC.fa[iFa].bGrp == CplBCType::cplBC_Dir) { - //WRITE(fid) cplBC.fa(iFa).Po, cplBC.fa(iFa).Pn + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&cplBC.fa[iFa].Po, double_size); + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&cplBC.fa[iFa].Pn, double_size); + genBC_writer.write( (char*)&double_size, int_size); } } - + for (int iFa = 0; iFa < cplBC.nFa; iFa++) { if (cplBC.fa[iFa].bGrp == CplBCType::cplBC_Neu) { - //WRITE(fid) cplBC.fa(iFa).Qo, cplBC.fa(iFa).Qn + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&cplBC.fa[iFa].Qo, double_size); + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&double_size, int_size); + genBC_writer.write( (char*)&cplBC.fa[iFa].Qn, double_size); + genBC_writer.write( (char*)&double_size, int_size); } } - //CLOSE(fid) + genBC_writer.close(); + + // Call genBC executable which reads the communication file GenBC.int + system(command); - //CALL SYSTEM(TRIM(cplBC.binPath)//" "//TRIM(cplBC.commuName)) + // Read outputs from genBC, which are in the same GenBC.int + std::ifstream genBC_reader(cplBC.commuName, std::ios::out | std::ios::binary); + if (!genBC_reader.is_open()) { + throw std::runtime_error("Failed to open the genBC interface file '" + cplBC.commuName + "' to read."); + } - //OPEN(fid,FILE=cplBC.commuName,STATUS='OLD',FORM='UNFORMATTED') + int size_buffer; for (int iFa = 0; iFa < cplBC.nFa; iFa++) { if (cplBC.fa[iFa].bGrp == CplBCType::cplBC_Dir) { - //READ(fid) cplBC.fa(iFa).y + genBC_reader.read( (char*)&size_buffer, int_size ); + genBC_reader.read( (char*)&cplBC.fa[iFa].y, size_buffer ); + genBC_reader.read( (char*)&size_buffer, int_size ); } } for (int iFa = 0; iFa < cplBC.nFa; iFa++) { if (cplBC.fa[iFa].bGrp == CplBCType::cplBC_Neu) { - //READ(fid) cplBC.fa(iFa).y + genBC_reader.read( (char*)&size_buffer, int_size ); + genBC_reader.read( (char*)&cplBC.fa[iFa].y, size_buffer ); + genBC_reader.read( (char*)&size_buffer, int_size ); } } - //CLOSE(fid) + + genBC_reader.close(); } + // If there are multiple procs (not sequential), broadcast genBC outputs to + // follower procs if (!cm.seq()) { Vector y(cplBC.nFa); - //ALLOCATE(y(cplBC.nFa)) if (cm.mas(cm_mod)) { for (int i = 0; i < cplBC.nFa; i++) { y(i) = cplBC.fa[i].y; } - //y = cplBC.fa.y; } cm.bcast(cm_mod, y); - //CALL cm.bcast(y) if (cm.slv(cm_mod)) { for (int i = 0; i < cplBC.nFa; i++) { cplBC.fa[i].y = y(i); } - // cplBC.fa.y = y; } - //DEALLOCATE(y) } } @@ -558,7 +666,8 @@ void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, con } -/// @brief Reproduces the Fortran 'SETBCCPL()' subrotutine. +/// @brief Coupled BC quantities are computed here. +/// Reproduces the Fortran 'SETBCCPL()' subrotutine. // void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) { @@ -573,9 +682,21 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) const int iEq = 0; auto& eq = com_mod.eq[iEq]; + // Determine current physics + auto& cDmn = com_mod.cDmn; + auto cPhys = eq.dmn[cDmn].phys; + + // Configuration in which to compute flowrate + auto cfg_o = MechanicalConfigurationType::reference; + auto cfg_n = MechanicalConfigurationType::reference; + + // If coupling scheme is implicit, calculate updated pressure and flowrate + // from 0D, as well as resistance from 0D using finite difference. if (cplBC.schm == CplBCType::cplBC_I) { calc_der_cpl_bc(com_mod, cm_mod); + // If coupling scheme is semi-implicit or explicit, only calculated updated + // pressure and flowrate from 0D } else { bool RCRflag = false; @@ -593,12 +714,35 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) if (ptr != -1) { + // Compute flowrates at 3D Neumann boundaries at timesteps n and n+1 if (utils::btest(bc.bType,iBC_Neu)) { - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 0, nsd-1); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, nsd-1); + // If struct or ustruct, use old and new configurations to compute flowrate integral + if ((cPhys == EquationType::phys_struct) || (cPhys == EquationType::phys_ustruct)) { + + // Must use follower pressure load for 0D coupling with struct/ustruct + if (!bc.flwP) { + throw std::runtime_error("[set_bc_cpl] Follower pressure load must be used for 0D coupling with struct/ustruct"); + } + cfg_o = MechanicalConfigurationType::old_timestep; + cfg_n = MechanicalConfigurationType::new_timestep; + } + // If fluid, FSI, or CMM, use reference configuration to compute flowrate integral + // Note that for FSI, mvMsh will modify geometry in gnnb() + else if ((cPhys == EquationType::phys_fluid) || (cPhys == EquationType::phys_FSI) || (cPhys == EquationType::phys_CMM)) { + cfg_o = MechanicalConfigurationType::reference; + cfg_n = MechanicalConfigurationType::reference; + } + else { + throw std::runtime_error("[set_bc_cpl] Invalid physics type for 0D coupling"); + } + + cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 0, nsd-1, false, cfg_o); + cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, nsd-1, false, cfg_n); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; - } else if (utils::btest(bc.bType,iBC_Dir)) { + } + // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 + else if (utils::btest(bc.bType,iBC_Dir)) { double area = com_mod.msh[iM].fa[iFa].area; cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, nsd) / area; cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd) / area; @@ -608,6 +752,8 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) } } + // Call genBC or cplBC to get updated pressures or flowrates. + // Updates pressure or flowrates stored in cplBC.fa[i].y if (cplBC.useGenBC) { set_bc::genBC_Integ_X(com_mod, cm_mod, "D"); } else { @@ -1296,12 +1442,10 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } - // Add Neumann BCs contribution to the LHS/RHS - // - // if follower pressure load. + // Add Neumann BCs contribution to the residual (and tangent if flwP) // if (lBc.flwP) { - eq_assem::b_neu_folw_p(com_mod, lFa, hg, Dg); + eq_assem::b_neu_folw_p(com_mod, lBc, lFa, hg, Dg); } else { eq_assem::b_assem_neu_bc(com_mod, lFa, hg, Yg); @@ -1314,7 +1458,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } -/// @brief Set Robin BC +/// @brief Set Robin BC contribution to residual and tangent // void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const double cs, const bool isN, const Array& Yg, const Array& Dg) diff --git a/Code/Source/svFSI/sv_struct.cpp b/Code/Source/svFSI/sv_struct.cpp index 161c727c..57918d01 100644 --- a/Code/Source/svFSI/sv_struct.cpp +++ b/Code/Source/svFSI/sv_struct.cpp @@ -106,7 +106,17 @@ void b_struct_2d(const ComMod& com_mod, const int eNoN, const double w, const Ve } } - +/// @brief Add follower pressure load contributions to the local residual and stiffness matrix. +/// @param com_mod +/// @param eNoN +/// @param w Gauss point weight times reference configuration area +/// @param N Shape function values at the Gauss point +/// @param Nx Shape function derivatives at the Gauss point +/// @param dl Displacement vector +/// @param hl Magnitude of pressure +/// @param nV Normal vector (in reference configuration) +/// @param lR Local residual +/// @param lK Local stiffness matrix void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, Array& lR, Array3& lK) @@ -144,6 +154,7 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve double h = 0.0; + // Compute deformation gradient tensor F for (int a = 0; a < eNoN; a++) { h = h + N(a)*hl(a); F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); @@ -165,7 +176,7 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve NxFi(1,a) = Nx(0,a)*Fi(0,1) + Nx(1,a)*Fi(1,1) + Nx(2,a)*Fi(2,1); NxFi(2,a) = Nx(0,a)*Fi(0,2) + Nx(1,a)*Fi(1,2) + Nx(2,a)*Fi(2,2); } - + // Compute N.F^-1, used for Nanson' formula da.n = J*dA*N.F^-1 nFi(0) = nV(0)*Fi(0,0) + nV(1)*Fi(1,0) + nV(2)*Fi(2,0); nFi(1) = nV(0)*Fi(0,1) + nV(1)*Fi(1,1) + nV(2)*Fi(2,1); nFi(2) = nV(0)*Fi(0,2) + nV(1)*Fi(1,2) + nV(2)*Fi(2,2); @@ -179,6 +190,7 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve debug << "wl: " << wl; #endif + // Compute follower pressure load contributions to the local residual and stiffness matrix for (int a = 0; a < eNoN; a++) { lR(0,a) = lR(0,a) - wl*N(a)*nFi(0); lR(1,a) = lR(1,a) - wl*N(a)*nFi(1); diff --git a/Code/Source/svFSI/ustruct.cpp b/Code/Source/svFSI/ustruct.cpp index 0a57d42f..8e661ca1 100644 --- a/Code/Source/svFSI/ustruct.cpp +++ b/Code/Source/svFSI/ustruct.cpp @@ -98,7 +98,18 @@ void b_ustruct_2d(const ComMod& com_mod, const int eNoN, const double w, const V } } - +/// @brief Add follower pressure load contributions to the local residual and stiffness matrix. +/// @param com_mod +/// @param eNoN +/// @param w Gauss point weight times reference configuration area +/// @param N Shape function values at the Gauss point +/// @param Nx Shape function derivatives at the Gauss point +/// @param dl Displacement vector +/// @param hl Magnitude of pressure +/// @param nV Normal vector (in reference configuration) +/// @param lR Local residual +/// @param lK Local stiffness matrix +/// @param lKd Local stiffness matrix (displacement) void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& dl, const Vector& hl, const Vector& nV, Array& lR, Array3& lK, Array3& lKd) diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index 122f9a9a..a9ac6999 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -1041,6 +1041,7 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& X, Array& Y) { Vector coef(lhs.nFaces); @@ -64,18 +73,20 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const int nsd = std::min(face.dof, dof); if (face.coupledFlag) { + // If face is shared across procs if (face.sharedFlag) { v = 0.0; - + // Setting vector v = int{N_A n_i} dGamma for (int a = 0; a < face.nNo; a++) { int Ac = face.glob(a); for (int i = 0; i < nsd; i++) { v(i,Ac) = face.valM(i,a); } } - + // Computing S = coef * v^T * X double S = coef(faIn) * dot::fsils_dot_v(dof, lhs.mynNo, lhs.commu, v, X); + // Computing Y = Y + v * S for (int a = 0; a < face.nNo; a++) { int Ac = face.glob(a); for (int i = 0; i < nsd; i++) { @@ -83,7 +94,10 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const } } - } else { + } + // If face is not shared across procs + else { + // Computing S = coef * v^T * X double S = 0.0; for (int a = 0; a < face.nNo; a++) { int Ac = face.glob(a); @@ -91,9 +105,9 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const S = S + face.valM(i,a)*X(i,Ac); } } - S = coef(faIn) * S; - + + // Computing Y = Y + v * S for (int a = 0; a < face.nNo; a++) { int Ac = face.glob(a); for (int i = 0; i < nsd; i++) { diff --git a/Code/Source/svFSILS/bc.cpp b/Code/Source/svFSILS/bc.cpp index 119644e8..6924b5c4 100644 --- a/Code/Source/svFSILS/bc.cpp +++ b/Code/Source/svFSILS/bc.cpp @@ -33,7 +33,8 @@ namespace fsi_linear_solver { -/// @brief Modifies: +/// @brief Sets up 0D contribution to tangent matrix for 0D-coupled boundary condition. +/// Modifies /// lhs.face[faIn].nNo /// lhs.face[faIn].dof /// lhs.face[faIn].bGrp @@ -93,6 +94,7 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t lhs.face[faIn].val = 0.0; } + // Synchronize val on nodes on boundary between procs if (lhs.commu.nTasks > 1) { int a = 0; @@ -136,6 +138,9 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t } +/// @brief Frees FSILS boundary condition structure +/// @param lhs +/// @param faIn void fsils_bc_free(FSILS_lhsType& lhs, int faIn) { //IF (.NOT.lhs%face(faIn)%foC) THEN @@ -153,5 +158,53 @@ void fsils_bc_free(FSILS_lhsType& lhs, int faIn) } -}; +/// @brief Updates lhs.face[faIn].val with new values from 'Val' parameter. +/// Val should contain the integral of the normal vector in the current (n+1) configuration, +/// used in the tangent contribution of the 0D-coupled boundary condition. +/// Since fsils_bc_create() was already called at initialization, we don't +/// need to reallocate data structures, just update values. +/// @param lhs +/// @param faIn +/// @param nNo +/// @param dof +/// @param Val +void fsils_bc_update(FSILS_lhsType& lhs, int faIn, int nNo, int dof, const Array& Val) +{ + using namespace consts; + + // Set lhs.face[faIn].val with new values from 'Val' parameter + if (Val.size() != 0) { + for (int a = 0; a < nNo; a++) { + for (int i = 0; i < Val.nrows(); i++) { + lhs.face[faIn].val(i,a) = Val(i,a); + } + } + } else { + lhs.face[faIn].val = 0.0; + } + + // Communicate update among procs + if (lhs.face[faIn].sharedFlag){ + Array v(dof,lhs.nNo); + v = 0.0; + + for (int a = 0; a < nNo; a++) { + int Ac = lhs.face[faIn].glob(a); + for (int i = 0; i < dof; i++) { + v(i,Ac) = lhs.face[faIn].val(i,a); + } + } + + fsils_commuv(lhs, dof, v); + + for (int a = 0; a < nNo; a++) { + int Ac = lhs.face[faIn].glob(a); + for (int i = 0; i < dof; i++) { + lhs.face[faIn].val(i,a) = v(i,Ac); + } + } + } + +} +}; \ No newline at end of file diff --git a/Code/Source/svFSILS/fsils_api.hpp b/Code/Source/svFSILS/fsils_api.hpp index 3795dbfd..47bac5f0 100644 --- a/Code/Source/svFSILS/fsils_api.hpp +++ b/Code/Source/svFSILS/fsils_api.hpp @@ -45,6 +45,8 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t void fsils_bc_free(FSILS_lhsType& lhs, int faIn); +void fsils_bc_update(FSILS_lhsType& lhs, int faIn, int nNo, int dof, const Array& Val); + void fsils_commus(const FSILS_lhsType& lhs, Vector& R); void fsils_commuv(const FSILS_lhsType& lhs, const int dof, Array& R); diff --git a/tests/cases/struct/LV_Holzapfel_passive/svFSI.xml b/tests/cases/struct/LV_Holzapfel_passive/svFSI.xml index e854a060..e08d1e22 100644 --- a/tests/cases/struct/LV_Holzapfel_passive/svFSI.xml +++ b/tests/cases/struct/LV_Holzapfel_passive/svFSI.xml @@ -44,10 +44,10 @@ mesh/fibersLongCells.vtu mesh/fibersSheetCells.vtu - 100.0 + 100.0 - + true @@ -55,16 +55,16 @@ 4 1e-10 - 1.0 - 1.0e6 + 1.0 + 1.0e6 0.483333 0.0 - 590.0 - 8.023 + 590.0 + 8.023 184720.0 16.026 24810.0 @@ -96,8 +96,8 @@ Robin - 1.0e7 - 5.0e2 + 1.0e7 + 5.0e2 1 diff --git a/tests/cases/struct/LV_NeoHookean_passive/README.md b/tests/cases/struct/LV_NeoHookean_passive/README.md new file mode 100644 index 00000000..0516e92c --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/README.md @@ -0,0 +1,2 @@ +This test case simulates an idealized left ventricle (LV) with a NeoHookean material model +being inflated with a linearly ramped pressure on the endocardial surface. \ No newline at end of file diff --git a/tests/cases/struct/LV_NeoHookean_passive/endo_pressure.dat b/tests/cases/struct/LV_NeoHookean_passive/endo_pressure.dat new file mode 100755 index 00000000..da125712 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/endo_pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3b2d4b24ddefba1359ed39aa892ab35b510a78259fe1caefa1bdf6e77f997955 +size 22 diff --git a/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.exterior.vtp b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/endo.vtp b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/epi.vtp b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/top.vtp b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/struct/LV_NeoHookean_passive/result_005.vtu b/tests/cases/struct/LV_NeoHookean_passive/result_005.vtu new file mode 100644 index 00000000..e01aa6d8 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/result_005.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3ecec14df7ebec99b5267e4f6bc97fff8e86b4cd7d897430c12c61e79bceb86f +size 248941 diff --git a/tests/cases/struct/LV_NeoHookean_passive/svFSI.xml b/tests/cases/struct/LV_NeoHookean_passive/svFSI.xml new file mode 100644 index 00000000..f820ef98 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive/svFSI.xml @@ -0,0 +1,100 @@ + + + + + 0 + 3 + 5 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 10 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + 100.0 + + + + + + + true + 3 + 10 + 1e-9 + + 1.0 + 1.0e5 + 0.483333 + + + ST91 + + 0.0 + + + + + true + true + true + true + true + true + true + true + + + + FSILS + 1e-9 + 1000 + 50 + + + + Dirichlet + 0.0 + + + + Neu + Unsteady + endo_pressure.dat + true + true + + + + + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/README.md b/tests/cases/struct/LV_NeoHookean_passive_genBC/README.md new file mode 100644 index 00000000..dfcc0365 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/README.md @@ -0,0 +1,29 @@ +This test case simulates an idealized left ventricle (LV) with a NeoHookean material model +coupled to a lumped-parameter network (LPN), implemented in genBC. The LPN consists of a large pressure source and large resistor, which together produce an approximately constant flowrate into +the LV. This inflates the LV at an approximately constant rate of change of volume. + +Before running svFSIplus, genBC_svFSIplus must be compiled. Navigate to genBC_svFSIplus +and run `make clean` then `make`. + +The results can be post-processed by running `process_results.py`. If run for 100 +timesteps, the results should be the following. + +![V3D vs. V0D](V3D_vs_V0D.png) + +*Plot of volume vs. time for this simulation. V_3D is computed from results.vtu +by extracting the endocardial surface, warping it according to displacement at +each timestep, capping it, then measuring the enclosed volume. V_0D is read +from AllData (produced by genBC).* + + +![dVdt3D vs. dVdt0D](dVdt3D_vs_dVdt0D.png) + +*Plot of flowrate (dVdt) vs. time for this simulation. dVdt_3D is computed from results.vtu +by extracting the endocardial surface, warping it according to displacement at +each timestep, then computing the velocity flux integral. +V_0D is computed from AllData (produced by genBC). As seen in the plot, we +ramp the flowrate from 0 to 100 cm^3/s over 10 timesteps.* + +![Pressure vs. Volume](pv_plot.png) + +*Plot of pressure vs. volume as the LV is inflated.* \ No newline at end of file diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png b/tests/cases/struct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png new file mode 100644 index 00000000..19b86258 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:63c0dccfa18fd514053bd6ff2ef226ca2c2b69176131c1b751fe8703223d2a01 +size 23969 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png b/tests/cases/struct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png new file mode 100644 index 00000000..c2bfc23b --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:16bf0c37fe2087828f83da5226dbefc09b9975fef1148d9965eebb39180f3198 +size 22077 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile new file mode 100644 index 00000000..384f82a2 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile @@ -0,0 +1,79 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the Makefile to build GenBC. +# +#-------------------------------------------------------------------- + +# HERE COMES THE DEFINITION + +include Makefile.in + +INCLUDES = -I./include + +MYFUN = Modules.f \ + USER.f \ + GenBC.f + +GenBC_EXE = genBC.exe + +DSYM_DIR = +ifeq ($(debug),1) + ifeq ($(OS),Darwin) + DSYM_DIR = $(GenBC_EXE).dSYM + endif +endif +############################################################################# +# AND HERE ARE THE RULES + +SRC = $(patsubst %,$(SRC_DIR)/%,$(MYFUN)) +OBJ = $(patsubst %.f,$(OBJ_DIR)/%.o,$(MYFUN)) + +#.PHONY: $(TEST) +#$(TEST): $(TEST:.exe=.f) $(GenBC_EXE) +# $(FORTRAN) $< $(FFLAGS) $(INCLUDES) -o $@ + +.PHONY: $(GenBC_EXE) +$(GenBC_EXE): $(OBJ) + $(FORTRAN) $(OBJ) $(FFLAGS) -o $@ + +$(OBJ): | $(OBJ_DIR) + +$(OBJ_DIR): + mkdir -p $@ + +$(OBJ_DIR)/%.o : $(SRC_DIR)/%.f + $(FORTRAN) $(FFLAGS) $(INCLUDES) -c $< -o $@ + +clean: + rm -r -f $(OBJ_DIR) $(GenBC_EXE) $(TEST) $(DSYM_DIR) diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile.in b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile.in new file mode 100644 index 00000000..7a434ad0 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/Makefile.in @@ -0,0 +1,165 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the definitions for building process. +# +#-------------------------------------------------------------------- + + +.SUFFIXES: .f .o + +OS:= $(shell uname) + +AR = ar rv +RANLIB = ar -ts + +CFLAGS = -O3 -DAdd_ +CXXFLAGS = -O3 -std=c++11 +FFLAGS = -O3 -cpp +FCFLAGS = -lstdc++ -cpp +#FCFLAGS = -lc++ -cpp ## This may work if lstdc++ is not found +OBJ_DIR = obj +SRC_DIR = src + +MAKE_GUI = 0 + +debug = 0 +ifeq ($(debug),1) + CFLAGS = -O0 -DAdd_ + CXXFLAGS = -O0 -std=c++11 + FFLAGS = -O0 -cpp +endif + +seq = 1 +ifeq ($(seq),1) + CC = gcc + CXX = g++ + FORTRAN = gfortran + CFLAGS += -DSEQ +else + CC = mpicc + CXX = mpicxx + FORTRAN = mpif90 +endif + +# ---------------------------------------------------------------- +# Normally you would not need to change any line beyond this point +# ---------------------------------------------------------------- + +# Here I am finding the compiler group +ifeq ($(seq),1) + F_COMP = $(FORTRAN) +else + F_COMP = $(firstword $(shell $(FORTRAN) -show)) +endif +ifeq ($(F_COMP),gfortran) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),gcc) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),g77) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),f95) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),ifort) + COMP_GRP = intel +endif +ifeq ($(F_COMP),pgf90) + COMP_GRP = pgi +endif +ifeq ($(F_COMP),pgf77) + COMP_GRP = pgi +endif + +ifeq ($(OS),Darwin) + COMP_GRP = gnu +endif + +# If profiling is requested +prof = 0 +ifeq ($(prof),1) + FFLAGS += -pg -g +endif + +# If debuging is requested +ifeq ($(debug),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -g -Wall -Wconversion -Wline-truncation -pedantic -fimplicit-none -fbacktrace -fbounds-check -p -fcheck=all #-ffpe-trap=invalid,zero,overflow,underflow + CXXFLAGS += -g -Wall -pedantic -fbounds-check + CFLAGS += -g -Wall -pedantic -fbounds-check + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -g -traceback -check all -check bounds -check uninit -ftrapuv -debug all -fpe0 + CFLAGS += -g -traceback -check-uninit -fpe0 + CXXFLAGS += -g -traceback -check-uninit -fpe0 + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += # You need to add debuging flag for pgi compiler here + endif + OBJ_DIR = dbg +endif + +# If openMP parallelization is requested +mp = 0 +ifeq ($(mp),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -fopenmp + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -openmp + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += -mp + endif +endif + +# To make directories a bit cleaner with intel compiler +ifeq ($(COMP_GRP),intel) + FFLAGS += -module $(OBJ_DIR) +endif +ifeq ($(COMP_GRP),gnu) + FFLAGS += -J $(OBJ_DIR) + CFLAGS += -J $(OBJ_DIR) + CXXFLAGS += -J $(OBJ_DIR) +endif + +#LAPACK Library +ifeq ($(COMP_GRP),gnu) + LAPACK_INC = + LAPACK_LIB = + LAPACK_LFLAGS = -llapack +endif + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f new file mode 100644 index 00000000..991bbb50 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f @@ -0,0 +1,4 @@ +c Value of your unknown at time equal to zero (This is going to be used +c ONLY when you are initiating simulation) + + tZeroX = (/ 0.00D+00 /) diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f new file mode 100644 index 00000000..af40e47d --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f @@ -0,0 +1,15 @@ +! Specifing the constants + REAL(KIND=8), PARAMETER :: ! Order of surfaces: 3~8 + & Rmyo = 4.51000000D-02, + & Lv = 2.74000000D-05, + & Kv = 8.16000000D-04, + & Kao = 8.59000000D-05, + & pi = 3.14159265D+00, + & Afill = 6D+01, + & tfill = 5D+00 +! The parameters which might change through the code + REAL(KIND=8) Clung, Rlung + Clung= 2.56000000D-02 + Rlung= 7.29000000D-01 + + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f new file mode 100644 index 00000000..3be431da --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f @@ -0,0 +1,230 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c Here data are received from Phasta and it setup the data required for +c integration of ODE's inside FINDX + PROGRAM GenBC + USE COM + USE, INTRINSIC :: IEEE_ARITHMETIC + IMPLICIT NONE + + INTEGER i, n, nTimeStep + REAL(KIND=8) t, dt, tFinal, pMin, temp + CHARACTER(LEN=2048) string + CHARACTER(LEN=32) sTmp + CHARACTER flag + + REAL(KIND=8), ALLOCATABLE :: Qi(:), Qf(:), Pi(:), Pf(:), X(:), + 2 Xo(:), f(:,:) + +c flag = I : Initializing +c flag = T : Iteration loop +c flag = L : Last iteration +c flag = D : Derivative +c +c Here is the period of time that you should integrate for +c each of these flags: +c +c Flags I D&T&L +c ^ ^ +c 3D code time step: N.............................................N+1 +c 0D code time step: 1......................................nTimeStep+1 +c Flowrates: Qi............................ ................Qf +c Time, t: tInitial..............................tInitial+tFinal + + + CALL INITIALIZE(nTimeStep) + +c******************************************************************** +c Block for reading the data from phasta + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED'); + READ (1) flag + READ (1) tFinal + READ (1) i + IF (i .NE. nDirichletSrfs) THEN + PRINT *, 'Error: Number of Dirichlet Surfaces from Phasta is:', + 2 i + PRINT *, 'While nDirichletSrfs is equal to:', nDirichletSrfs + PRINT * + PRINT *, 'Number of Dirichlet surfaces defined in solver.inp', + 2 ' should match with nDirichletSrfs defined in USER.f' + STOP + END IF + READ (1) i + IF (i .NE. nNeumannSrfs) THEN + PRINT *, 'Error: Number of Neumann Surfaces from Phasta is:', + 2 i + PRINT *, 'While nNeumannSrfs is equal to:', nNeumannSrfs + PRINT * + PRINT *, 'Number of Neumann surfaces defined in solver.inp', + 2 ' should match with nNeumannSrfs defined in USER.f' + STOP + END IF + IF (nDirichletSrfs .GT. 0) THEN + IF (qCorr .OR. pCorr) THEN + PRINT *, 'You should only use P/Q correction when all', + 2 ' the surfaces are Neumann surfaces' + STOP + END IF + END IF + + ALLOCATE (Pi(nDirichletSrfs), Pf(nDirichletSrfs), + 2 PDirichlet(nDirichletSrfs,4)) + ALLOCATE (Qi(nNeumannSrfs), Qf(nNeumannSrfs), + 2 QNeumann(nNeumannSrfs,4)) + + DO i=1, nDirichletSrfs + READ(1) Pi(i), Pf(i) + END DO + Pi = Pi/pConv + Pf = Pf/pConv + + DO i=1,nNeumannSrfs + READ (1) Qi(i), Qf(i) +c Print to console and file pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Flowrates to genBC: ', Qi(i), Qf(i) + END DO + CLOSE(1) + Qi = Qi/qConv + Qf = Qf/qConv + +c******************************************************************** +c Block for initializing the unknowns + + IF (flag .EQ. 'I') THEN + nTimeStep = 0 + ELSE + dt = tFinal/REAL(nTimeStep,8) + END IF + ALLOCATE (X(nUnknowns), Xo(nUnknowns), f(nUnknowns,4), + 2 offset(nUnknowns), Xprint(nXprint)) + Xprint = 0D0 + + OPEN (1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED') + READ (1) t + DO i=1,nUnknowns + READ (1) Xo(i) + END DO + CLOSE (1) + +c******************************************************************** +c Setting up the system of equations + offset = 0D0 + DO n=1, nTimeStep + DO i=1, 4 + temp = (REAL(n-1,8) + REAL(i-1,8)/3D0)/REAL(nTimeStep,8) + + QNeumann(:,i) = Qi + (Qf-Qi)*temp + PDirichlet(:,i) = Pi + (Pf-Pi)*temp + + IF (qCorr) THEN + temp = SUM(QNeumann(:,i))/REAL(nNeumannSrfs,8) + QNeumann(:,i) = QNeumann(:,i) - temp + END IF + END DO + + X = Xo + CALL FINDF (t, X, f(:,1), QNeumann(:,1), + 2 PDirichlet(:,1)) + X = Xo + dt*f(:,1)/3D0 + + CALL FINDF (t+dt/3D0, X, f(:,2), QNeumann(:,2), + 2 PDirichlet(:,2)) + X = Xo - dt*f(:,1)/3D0 + dt*f(:,2) + + CALL FINDF (t+dt*2D0/3D0, X, f(:,3), QNeumann(:,3), + 2 PDirichlet(:,3)) + X = Xo + dt*f(:,1) - dt*f(:,2) + dt*f(:,3) + + CALL FINDF (t+dt, X, f(:,4), QNeumann(:,4), + 2 PDirichlet(:,4)) + + f(:,1) = (f(:,1) + 3D0*f(:,2) + 3D0*f(:,3) + f(:,4))/8D0 + Xo = Xo + dt*f(:,1) + t = t + dt + END DO + +c******************************************************************** +c Time to write the results + X = Xo + IF (pCorr .AND. flag.NE.'D') THEN + pMin = X(srfToXPtr(1)) + DO i=2, nNeumannSrfs + IF (X(srfToXPtr(i)) .LT. pMin) THEN + pMin = X(srfToXPtr(i)) + END IF + END DO + ELSE + pMin = 0D0 + END IF + +c Writing nDirichlet flowrates here + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED') + DO i=1, nDirichletSrfs + IF(IEEE_IS_NAN(X(srfToXdPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF + WRITE (1) X(srfToXdPtr(i))*qConv + END DO + +c Writing nNeumannSrfs pressures here + DO i=1, nNeumannSrfs + IF(IEEE_IS_NAN(X(srfToXPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF +c Print out pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Pressure from genBC: ', offset(srfToXPtr(i)) + +c Here, saving offset value to X itself. Necessary for keeping track of +c pressure from one time step to the next + X(srfToXPtr(i)) = offset(srfToXPtr(i)) + +c WRITE (1) (X(srfToXPtr(i)) - pMin +c 2 + offset(srfToXPtr(i)))*pConv + WRITE (1) (X(srfToXPtr(i)))*pConv + END DO + CLOSE(1) + + IF (flag .EQ. 'L') THEN + OPEN(1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED'); + WRITE (1) t + DO i=1,nUnknowns + WRITE (1) X(i) + END DO + CLOSE(1) + +c PRINT *,'Before AllData' + + OPEN(1, FILE='AllData', STATUS='UNKNOWN', ACCESS='APPEND'); + string = '' + DO i=1, nUnknowns + WRITE (sTmp,"(ES14.6E2)") X(i) + string = TRIM(string)//sTmp + END DO + DO i=1, nXprint + WRITE (sTmp,"(ES14.6E2)") Xprint(i) + string = TRIM(string)//sTmp + END DO + WRITE (1,"(A)") TRIM(string) + CLOSE(1) + END IF + + DEALLOCATE (Pi) + DEALLOCATE (Pf) + DEALLOCATE (PDirichlet) + DEALLOCATE (Qi) + DEALLOCATE (Qf) + DEALLOCATE (QNeumann) + DEALLOCATE (X) + DEALLOCATE (Xo) + DEALLOCATE (f) + DEALLOCATE (srfToXdPtr) + DEALLOCATE (srfToXPtr) + DEALLOCATE (offset) + DEALLOCATE (Xprint) + + END PROGRAM GenBC diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f new file mode 100644 index 00000000..6e472f3b --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f @@ -0,0 +1,17 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + + MODULE COM + + LOGICAL pCorr, qCorr + + INTEGER nUnknowns, nDirichletSrfs, nNeumannSrfs, nXprint + + REAL(KIND=8) pConv, qConv + + INTEGER, ALLOCATABLE :: srfToXPtr(:), srfToXdPtr(:) + + REAL(KIND=8), ALLOCATABLE :: QNeumann(:,:), PDirichlet(:,:), + 2 offset(:), Xprint(:) + + END MODULE COM diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/USER.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/USER.f new file mode 100644 index 00000000..e48f7c19 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC/src/USER.f @@ -0,0 +1,143 @@ +c MUSC 2 Immediate postop +c Ethan Kung keo@ucsd.edu + +c Created by Mahdi Esmaily Moghadam 12-01-2010 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c This subroutine initializes the parameters, you need to read the +c comments and specify them carefuly +c-------------------------------------------------------------------- +c This is an example for RCR boundary condition with parameters +c Rd, R, and C which are distal, and proximal resistance and capacitor. + + SUBROUTINE INITIALIZE(nTimeStep) + USE COM + IMPLICIT NONE + INTENT(OUT) nTimeStep + + LOGICAl ierr + INTEGER i, nTimeStep + REAL(KIND=8), ALLOCATABLE :: tZeroX(:) +c +c******************************************************************** +c For instance if pressure in 3D solver is in cgs and here mmHg +c pConv=1334=1.334D3, also if flowrate in 3D solver is mm^3/s and +c here is mL/s qConv=1000=1D3. In the case both solver are using same +c unites you can set these two conversion coefficients to 1D0 + pConv = 1D0 + qConv = 1D0 + +c Only when all the surfaces of you model are coupled with NeumannSrfs +c you may set this to .TRUE. + pCorr = .FALSE. + qCorr = .FALSE. + +c******************************************************************** +c Block for your inputs + +c These two value should match to that of defined in solver.inp + nDirichletSrfs = 0 + nNeumannSrfs = 1 +c Number of unknowns that you need inside your lumped parameter network + nUnknowns = 2 +c Number of time step between N and N+alpha, for higher values this code +c would be more accurate and more costly + nTimeStep = 10 + +c Number of parameters to be printed in AllData file (the first +c nUnknowns columns are the "X" values) + nXprint = 2 + +c-------------------------------------------------------------------- +c You don't need to change this part ! + ALLOCATE (tZeroX(nUnknowns), srfToXdPtr(nDirichletSrfs)) ! + ALLOCATE (srfToXPtr(nNeumannSrfs)) ! + tZeroX = 0D0 +c-------------------------------------------------------------------- + + INCLUDE "initial_values_final.f" + +c-------------------------------------------------------------------- +c You don't need to change this part ! + INQUIRE (FILE='InitialData', EXIST=ierr) ! + IF (.NOT.ierr) THEN ! +c PRINT *, 'Initializing unknowns in LPM' ! + OPEN (1, FILE='InitialData',STATUS='NEW',FORM='UNFORMATTED')! + WRITE (1) 0D0 ! + DO i=1, nUnknowns ! + WRITE (1) tZeroX(i) ! + END DO ! + CLOSE(1) ! + END IF ! +c-------------------------------------------------------------------- + +c Surface to X pointer: this defines which Unknown corresponds to which +c suface in "List of Neumann Surfaces" inside solver.inp +c For example, if you have "List of Neumann Surfaces= 2 8 4 ...." +c and you set "srfToXPtr = (/5,3,9,.../)" +C this means X(5) corresponds to surface 2, X(3) <-> surface 8, +c and X(9) <-> surface 4 +c Also, Q(1) corresponds to the first element in srfToXPtr, in this +c example, Q(1)<=>X(5), Q(2)<=>X(3) + srfToXPtr = (/1/) +c srfToXdPtr = (/1/) + + END SUBROUTINE INITIALIZE + +c#################################################################### +c Here you should find the f_i=dx_i/dt, based on the following parameters: +c current x_i: x(i) +c Current time: t +c Flowrates from 3D code: Q(i) +c Pressure from Dirichlet faces: P(i) + + SUBROUTINE FINDF(t, x, f, Q, P) + USE COM + IMPLICIT NONE + INTENT(IN) t, Q + INTENT(OUT) f + + REAL(KIND=8) t, x(nUnknowns), f(nUnknowns), Q(nNeumannSrfs), + 2 P(nDirichletSrfs) + +! These are the dumy variables + REAL(KIND=8) Pivc,Ppwc,Pcon,Qout, Phigh, Pmax, tramp, R + + INCLUDE "parameters_final.f" + +!******************************************************************** +! This is the only case that x is directly manipulated, ugly, but the +! only way!!! + +! Implement a resistance filling network so that the ventricle pressure +! Pv = Phigh - Q*R +! This fills the sphere at an approximately constant rate + +! x(1): Pressure applied to endo surface. Actually, we will assign the pressure +! by setting the offset instead +! x(2): Accumulated volume during inflation. Not necessary, but for insight + +! Define a pressure ramp for Phigh (1e7 dyne/cm^2 in 0.11 second, then constant) + Pmax = 1.0D+7 + tramp = 0.1 + IF (t .LT. tramp) THEN + Phigh = Pmax * t / tramp + ELSE + Phigh = Pmax + ENDIF + R = 1.0D+5 + +! The main body of equations + f(1) = 0 ! The endo surface pressure (we will set it with offset) + f(2) = -Q(1) ! Integrate the flow rate to compute accumulated volume + +! Define offset pressure. + offset(1) = Phigh + Q(1) * R + +c Assign the additional parameters to be printed + Xprint(1) = offset(1) ! This is the actual value of pressure for this sim + Xprint(2) = Q(1) ! The flowrate/dVdt + RETURN + END SUBROUTINE FINDF + + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile new file mode 100644 index 00000000..384f82a2 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile @@ -0,0 +1,79 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the Makefile to build GenBC. +# +#-------------------------------------------------------------------- + +# HERE COMES THE DEFINITION + +include Makefile.in + +INCLUDES = -I./include + +MYFUN = Modules.f \ + USER.f \ + GenBC.f + +GenBC_EXE = genBC.exe + +DSYM_DIR = +ifeq ($(debug),1) + ifeq ($(OS),Darwin) + DSYM_DIR = $(GenBC_EXE).dSYM + endif +endif +############################################################################# +# AND HERE ARE THE RULES + +SRC = $(patsubst %,$(SRC_DIR)/%,$(MYFUN)) +OBJ = $(patsubst %.f,$(OBJ_DIR)/%.o,$(MYFUN)) + +#.PHONY: $(TEST) +#$(TEST): $(TEST:.exe=.f) $(GenBC_EXE) +# $(FORTRAN) $< $(FFLAGS) $(INCLUDES) -o $@ + +.PHONY: $(GenBC_EXE) +$(GenBC_EXE): $(OBJ) + $(FORTRAN) $(OBJ) $(FFLAGS) -o $@ + +$(OBJ): | $(OBJ_DIR) + +$(OBJ_DIR): + mkdir -p $@ + +$(OBJ_DIR)/%.o : $(SRC_DIR)/%.f + $(FORTRAN) $(FFLAGS) $(INCLUDES) -c $< -o $@ + +clean: + rm -r -f $(OBJ_DIR) $(GenBC_EXE) $(TEST) $(DSYM_DIR) diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in new file mode 100644 index 00000000..aa00a611 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in @@ -0,0 +1,165 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the definitions for building process. +# +#-------------------------------------------------------------------- + + +.SUFFIXES: .f .o + +OS:= $(shell uname) + +AR = ar rv +RANLIB = ar -ts + +CFLAGS = -O3 -DAdd_ +CXXFLAGS = -O3 -std=c++11 +FFLAGS = -O3 -cpp +FCFLAGS = -lstdc++ -cpp +#FCFLAGS = -lc++ -cpp ## This may work if lstdc++ is not found +OBJ_DIR = obj +SRC_DIR = src + +MAKE_GUI = 0 + +debug = 0 +ifeq ($(debug),1) + CFLAGS = -O0 -DAdd_ + CXXFLAGS = -O0 -std=c++11 + FFLAGS = -O0 -cpp +endif + +seq = 0 +ifeq ($(seq),1) + CC = gcc + CXX = g++ + FORTRAN = gfortran + CFLAGS += -DSEQ +else + CC = mpicc + CXX = mpicxx + FORTRAN = mpif90 +endif + +# ---------------------------------------------------------------- +# Normally you would not need to change any line beyond this point +# ---------------------------------------------------------------- + +# Here I am finding the compiler group +ifeq ($(seq),1) + F_COMP = $(FORTRAN) +else + F_COMP = $(firstword $(shell $(FORTRAN) -show)) +endif +ifeq ($(F_COMP),gfortran) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),gcc) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),g77) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),f95) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),ifort) + COMP_GRP = intel +endif +ifeq ($(F_COMP),pgf90) + COMP_GRP = pgi +endif +ifeq ($(F_COMP),pgf77) + COMP_GRP = pgi +endif + +ifeq ($(OS),Darwin) + COMP_GRP = gnu +endif + +# If profiling is requested +prof = 0 +ifeq ($(prof),1) + FFLAGS += -pg -g +endif + +# If debuging is requested +ifeq ($(debug),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -g -Wall -Wconversion -Wline-truncation -pedantic -fimplicit-none -fbacktrace -fbounds-check -p -fcheck=all #-ffpe-trap=invalid,zero,overflow,underflow + CXXFLAGS += -g -Wall -pedantic -fbounds-check + CFLAGS += -g -Wall -pedantic -fbounds-check + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -g -traceback -check all -check bounds -check uninit -ftrapuv -debug all -fpe0 + CFLAGS += -g -traceback -check-uninit -fpe0 + CXXFLAGS += -g -traceback -check-uninit -fpe0 + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += # You need to add debuging flag for pgi compiler here + endif + OBJ_DIR = dbg +endif + +# If openMP parallelization is requested +mp = 0 +ifeq ($(mp),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -fopenmp + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -openmp + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += -mp + endif +endif + +# To make directories a bit cleaner with intel compiler +ifeq ($(COMP_GRP),intel) + FFLAGS += -module $(OBJ_DIR) +endif +ifeq ($(COMP_GRP),gnu) + FFLAGS += -J $(OBJ_DIR) + CFLAGS += -J $(OBJ_DIR) + CXXFLAGS += -J $(OBJ_DIR) +endif + +#LAPACK Library +ifeq ($(COMP_GRP),gnu) + LAPACK_INC = + LAPACK_LIB = + LAPACK_LFLAGS = -llapack +endif + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f new file mode 100644 index 00000000..991bbb50 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f @@ -0,0 +1,4 @@ +c Value of your unknown at time equal to zero (This is going to be used +c ONLY when you are initiating simulation) + + tZeroX = (/ 0.00D+00 /) diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f new file mode 100644 index 00000000..af40e47d --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f @@ -0,0 +1,15 @@ +! Specifing the constants + REAL(KIND=8), PARAMETER :: ! Order of surfaces: 3~8 + & Rmyo = 4.51000000D-02, + & Lv = 2.74000000D-05, + & Kv = 8.16000000D-04, + & Kao = 8.59000000D-05, + & pi = 3.14159265D+00, + & Afill = 6D+01, + & tfill = 5D+00 +! The parameters which might change through the code + REAL(KIND=8) Clung, Rlung + Clung= 2.56000000D-02 + Rlung= 7.29000000D-01 + + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f new file mode 100644 index 00000000..1d3d8908 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f @@ -0,0 +1,232 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c Here data are received from Phasta and it setup the data required for +c integration of ODE's inside FINDX + PROGRAM GenBC + USE COM + USE, INTRINSIC :: IEEE_ARITHMETIC + IMPLICIT NONE + + INTEGER i, n, nTimeStep + REAL(KIND=8) t, dt, tFinal, pMin, temp + CHARACTER(LEN=2048) string + CHARACTER(LEN=32) sTmp + CHARACTER flag + + REAL(KIND=8), ALLOCATABLE :: Qi(:), Qf(:), Pi(:), Pf(:), X(:), + 2 Xo(:), f(:,:) + +c flag = I : Initializing +c flag = T : Iteration loop +c flag = L : Last iteration +c flag = D : Derivative +c +c Here is the period of time that you should integrate for +c each of these flags: +c +c Flags I D&T&L +c ^ ^ +c 3D code time step: N.............................................N+1 +c 0D code time step: 1......................................nTimeStep+1 +c Flowrates: Qi............................ ................Qf +c Time, t: tInitial..............................tInitial+tFinal + + + CALL INITIALIZE(nTimeStep) + +c******************************************************************** +c Block for reading the data from phasta + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED'); + READ (1) flag + READ (1) tFinal + READ (1) i + IF (i .NE. nDirichletSrfs) THEN + PRINT *, 'Error: Number of Dirichlet Surfaces from Phasta is:', + 2 i + PRINT *, 'While nDirichletSrfs is equal to:', nDirichletSrfs + PRINT * + PRINT *, 'Number of Dirichlet surfaces defined in solver.inp', + 2 ' should match with nDirichletSrfs defined in USER.f' + STOP + END IF + READ (1) i + IF (i .NE. nNeumannSrfs) THEN + PRINT *, 'Error: Number of Neumann Surfaces from Phasta is:', + 2 i + PRINT *, 'While nNeumannSrfs is equal to:', nNeumannSrfs + PRINT * + PRINT *, 'Number of Neumann surfaces defined in solver.inp', + 2 ' should match with nNeumannSrfs defined in USER.f' + STOP + END IF + IF (nDirichletSrfs .GT. 0) THEN + IF (qCorr .OR. pCorr) THEN + PRINT *, 'You should only use P/Q correction when all', + 2 ' the surfaces are Neumann surfaces' + STOP + END IF + END IF + + ALLOCATE (Pi(nDirichletSrfs), Pf(nDirichletSrfs), + 2 PDirichlet(nDirichletSrfs,4)) + ALLOCATE (Qi(nNeumannSrfs), Qf(nNeumannSrfs), + 2 QNeumann(nNeumannSrfs,4)) + + DO i=1, nDirichletSrfs + READ(1) Pi(i) + READ(1) Pf(i) + END DO + Pi = Pi/pConv + Pf = Pf/pConv + + DO i=1,nNeumannSrfs + READ (1) Qi(i) + READ (1) Qf(i) +c Print to console and file pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Flowrates to genBC: ', Qi(i), Qf(i) + END DO + CLOSE(1) + Qi = Qi/qConv + Qf = Qf/qConv + +c******************************************************************** +c Block for initializing the unknowns + + IF (flag .EQ. 'I') THEN + nTimeStep = 0 + ELSE + dt = tFinal/REAL(nTimeStep,8) + END IF + ALLOCATE (X(nUnknowns), Xo(nUnknowns), f(nUnknowns,4), + 2 offset(nUnknowns), Xprint(nXprint)) + Xprint = 0D0 + + OPEN (1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED') + READ (1) t + DO i=1,nUnknowns + READ (1) Xo(i) + END DO + CLOSE (1) + +c******************************************************************** +c Setting up the system of equations + offset = 0D0 + DO n=1, nTimeStep + DO i=1, 4 + temp = (REAL(n-1,8) + REAL(i-1,8)/3D0)/REAL(nTimeStep,8) + + QNeumann(:,i) = Qi + (Qf-Qi)*temp + PDirichlet(:,i) = Pi + (Pf-Pi)*temp + + IF (qCorr) THEN + temp = SUM(QNeumann(:,i))/REAL(nNeumannSrfs,8) + QNeumann(:,i) = QNeumann(:,i) - temp + END IF + END DO + + X = Xo + CALL FINDF (t, X, f(:,1), QNeumann(:,1), + 2 PDirichlet(:,1)) + X = Xo + dt*f(:,1)/3D0 + + CALL FINDF (t+dt/3D0, X, f(:,2), QNeumann(:,2), + 2 PDirichlet(:,2)) + X = Xo - dt*f(:,1)/3D0 + dt*f(:,2) + + CALL FINDF (t+dt*2D0/3D0, X, f(:,3), QNeumann(:,3), + 2 PDirichlet(:,3)) + X = Xo + dt*f(:,1) - dt*f(:,2) + dt*f(:,3) + + CALL FINDF (t+dt, X, f(:,4), QNeumann(:,4), + 2 PDirichlet(:,4)) + + f(:,1) = (f(:,1) + 3D0*f(:,2) + 3D0*f(:,3) + f(:,4))/8D0 + Xo = Xo + dt*f(:,1) + t = t + dt + END DO + +c******************************************************************** +c Time to write the results + X = Xo + IF (pCorr .AND. flag.NE.'D') THEN + pMin = X(srfToXPtr(1)) + DO i=2, nNeumannSrfs + IF (X(srfToXPtr(i)) .LT. pMin) THEN + pMin = X(srfToXPtr(i)) + END IF + END DO + ELSE + pMin = 0D0 + END IF + +c Writing nDirichlet flowrates here + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED') + DO i=1, nDirichletSrfs + IF(IEEE_IS_NAN(X(srfToXdPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF + WRITE (1) X(srfToXdPtr(i))*qConv + END DO + +c Writing nNeumannSrfs pressures here + DO i=1, nNeumannSrfs + IF(IEEE_IS_NAN(X(srfToXPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF +c Print out pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Pressure from genBC: ', offset(srfToXPtr(i)) + +c Here, saving offset value to X itself. Necessary for keeping track of +c pressure from one time step to the next + X(srfToXPtr(i)) = offset(srfToXPtr(i)) + +c WRITE (1) (X(srfToXPtr(i)) - pMin +c 2 + offset(srfToXPtr(i)))*pConv + WRITE (1) (X(srfToXPtr(i)))*pConv + END DO + CLOSE(1) + + IF (flag .EQ. 'L') THEN + OPEN(1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED'); + WRITE (1) t + DO i=1,nUnknowns + WRITE (1) X(i) + END DO + CLOSE(1) + +c PRINT *,'Before AllData' + + OPEN(1, FILE='AllData', STATUS='UNKNOWN', ACCESS='APPEND'); + string = '' + DO i=1, nUnknowns + WRITE (sTmp,"(ES14.6E2)") X(i) + string = TRIM(string)//sTmp + END DO + DO i=1, nXprint + WRITE (sTmp,"(ES14.6E2)") Xprint(i) + string = TRIM(string)//sTmp + END DO + WRITE (1,"(A)") TRIM(string) + CLOSE(1) + END IF + + DEALLOCATE (Pi) + DEALLOCATE (Pf) + DEALLOCATE (PDirichlet) + DEALLOCATE (Qi) + DEALLOCATE (Qf) + DEALLOCATE (QNeumann) + DEALLOCATE (X) + DEALLOCATE (Xo) + DEALLOCATE (f) + DEALLOCATE (srfToXdPtr) + DEALLOCATE (srfToXPtr) + DEALLOCATE (offset) + DEALLOCATE (Xprint) + + END PROGRAM GenBC diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f new file mode 100644 index 00000000..6e472f3b --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f @@ -0,0 +1,17 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + + MODULE COM + + LOGICAL pCorr, qCorr + + INTEGER nUnknowns, nDirichletSrfs, nNeumannSrfs, nXprint + + REAL(KIND=8) pConv, qConv + + INTEGER, ALLOCATABLE :: srfToXPtr(:), srfToXdPtr(:) + + REAL(KIND=8), ALLOCATABLE :: QNeumann(:,:), PDirichlet(:,:), + 2 offset(:), Xprint(:) + + END MODULE COM diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f new file mode 100644 index 00000000..e48f7c19 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f @@ -0,0 +1,143 @@ +c MUSC 2 Immediate postop +c Ethan Kung keo@ucsd.edu + +c Created by Mahdi Esmaily Moghadam 12-01-2010 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c This subroutine initializes the parameters, you need to read the +c comments and specify them carefuly +c-------------------------------------------------------------------- +c This is an example for RCR boundary condition with parameters +c Rd, R, and C which are distal, and proximal resistance and capacitor. + + SUBROUTINE INITIALIZE(nTimeStep) + USE COM + IMPLICIT NONE + INTENT(OUT) nTimeStep + + LOGICAl ierr + INTEGER i, nTimeStep + REAL(KIND=8), ALLOCATABLE :: tZeroX(:) +c +c******************************************************************** +c For instance if pressure in 3D solver is in cgs and here mmHg +c pConv=1334=1.334D3, also if flowrate in 3D solver is mm^3/s and +c here is mL/s qConv=1000=1D3. In the case both solver are using same +c unites you can set these two conversion coefficients to 1D0 + pConv = 1D0 + qConv = 1D0 + +c Only when all the surfaces of you model are coupled with NeumannSrfs +c you may set this to .TRUE. + pCorr = .FALSE. + qCorr = .FALSE. + +c******************************************************************** +c Block for your inputs + +c These two value should match to that of defined in solver.inp + nDirichletSrfs = 0 + nNeumannSrfs = 1 +c Number of unknowns that you need inside your lumped parameter network + nUnknowns = 2 +c Number of time step between N and N+alpha, for higher values this code +c would be more accurate and more costly + nTimeStep = 10 + +c Number of parameters to be printed in AllData file (the first +c nUnknowns columns are the "X" values) + nXprint = 2 + +c-------------------------------------------------------------------- +c You don't need to change this part ! + ALLOCATE (tZeroX(nUnknowns), srfToXdPtr(nDirichletSrfs)) ! + ALLOCATE (srfToXPtr(nNeumannSrfs)) ! + tZeroX = 0D0 +c-------------------------------------------------------------------- + + INCLUDE "initial_values_final.f" + +c-------------------------------------------------------------------- +c You don't need to change this part ! + INQUIRE (FILE='InitialData', EXIST=ierr) ! + IF (.NOT.ierr) THEN ! +c PRINT *, 'Initializing unknowns in LPM' ! + OPEN (1, FILE='InitialData',STATUS='NEW',FORM='UNFORMATTED')! + WRITE (1) 0D0 ! + DO i=1, nUnknowns ! + WRITE (1) tZeroX(i) ! + END DO ! + CLOSE(1) ! + END IF ! +c-------------------------------------------------------------------- + +c Surface to X pointer: this defines which Unknown corresponds to which +c suface in "List of Neumann Surfaces" inside solver.inp +c For example, if you have "List of Neumann Surfaces= 2 8 4 ...." +c and you set "srfToXPtr = (/5,3,9,.../)" +C this means X(5) corresponds to surface 2, X(3) <-> surface 8, +c and X(9) <-> surface 4 +c Also, Q(1) corresponds to the first element in srfToXPtr, in this +c example, Q(1)<=>X(5), Q(2)<=>X(3) + srfToXPtr = (/1/) +c srfToXdPtr = (/1/) + + END SUBROUTINE INITIALIZE + +c#################################################################### +c Here you should find the f_i=dx_i/dt, based on the following parameters: +c current x_i: x(i) +c Current time: t +c Flowrates from 3D code: Q(i) +c Pressure from Dirichlet faces: P(i) + + SUBROUTINE FINDF(t, x, f, Q, P) + USE COM + IMPLICIT NONE + INTENT(IN) t, Q + INTENT(OUT) f + + REAL(KIND=8) t, x(nUnknowns), f(nUnknowns), Q(nNeumannSrfs), + 2 P(nDirichletSrfs) + +! These are the dumy variables + REAL(KIND=8) Pivc,Ppwc,Pcon,Qout, Phigh, Pmax, tramp, R + + INCLUDE "parameters_final.f" + +!******************************************************************** +! This is the only case that x is directly manipulated, ugly, but the +! only way!!! + +! Implement a resistance filling network so that the ventricle pressure +! Pv = Phigh - Q*R +! This fills the sphere at an approximately constant rate + +! x(1): Pressure applied to endo surface. Actually, we will assign the pressure +! by setting the offset instead +! x(2): Accumulated volume during inflation. Not necessary, but for insight + +! Define a pressure ramp for Phigh (1e7 dyne/cm^2 in 0.11 second, then constant) + Pmax = 1.0D+7 + tramp = 0.1 + IF (t .LT. tramp) THEN + Phigh = Pmax * t / tramp + ELSE + Phigh = Pmax + ENDIF + R = 1.0D+5 + +! The main body of equations + f(1) = 0 ! The endo surface pressure (we will set it with offset) + f(2) = -Q(1) ! Integrate the flow rate to compute accumulated volume + +! Define offset pressure. + offset(1) = Phigh + Q(1) * R + +c Assign the additional parameters to be printed + Xprint(1) = offset(1) ! This is the actual value of pressure for this sim + Xprint(2) = Q(1) ! The flowrate/dVdt + RETURN + END SUBROUTINE FINDF + + diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/post_processing_functions.py b/tests/cases/struct/LV_NeoHookean_passive_genBC/post_processing_functions.py new file mode 100644 index 00000000..503dbeb5 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/post_processing_functions.py @@ -0,0 +1,700 @@ +# Before running on Sherlock, need to have installed: +# +# pyvista: pip3 install pyvista +# ffmpy: pip3 install ffmpy +# +# Also, need to load the following required modules on Sherlock: +# +# ml system mesa ffmpeg +import glob +import os +import pyvista as pv # for VTK mesh manipulations +import numpy as np +#import ffmpy # to convert .avi movie to .mp4 movie + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start time, end time, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSI simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_time, end_time, step): A tuple of 3 integers, giving the first + time step of results to process, the last time step, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_time = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start time:', start_time) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_time = int("".join([i for i in end_file_name if i.isdigit()])) + print('End time:', end_time) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_time + print('Step:', step) + + return (start_time, end_time, step) + +def calc_volume_struct(start_time, end_time, step, results_folder, reference_surface): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI struct simulation, in which a model of the myocardium is inflated. + + Calculate the volume in the following steps + 1) Sample the result.vtu file onto the reference surface + 2) Warp the samples surface by the Displacement + 3) Flat fill any holes in the warped surface + 4) Calculate the volume of the warped and filled surface + + The units of volume are whatever units used in .vtu files, cubed. For example, + if units of length in the .vtu files are microns, then the volume calculated + here is cubic microns. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + reference_surface: The absolute file path of the .vtp file containing + the undeformed surface corresponding to the deformed surface of which + we want to compute the volume. + + output_folder: The absolute file path of the folder to which we output + intermediary results of the volume calculation process. + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_struct_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating volumes ##') + + # Load reference surface onto which we sample + ref_surface = pv.read(f"{reference_surface}") + + # Compute initial lumen volume + t = [] + vol = [] + k = 0 + # The initial volume is obtain by computing the lumen volume of the reference + # configuration, which is simply volume of the filled in reference surface. + ref_lumen = ref_surface.fill_holes(100) # 100 is the largest size of hole to fill + + # -------------------------------------------------------------------- + + # Recompute normals, incase the normals of the cap are opposite + ref_lumen.compute_normals(inplace=True) + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + ref_lumen.save(f'{output_folder}/resampled_warped_and_filled_{k:03d}.vtp') + + # Add time and volume to arrays + t.append(k) + vol.append(ref_lumen.volume) + + print(f"Iteration: {k}, Volume: {ref_lumen.volume}") + + + # Loop through results files at each time > 0 + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}_cpp.vtu")) + + # Sample result onto ref_surface + resampled = ref_surface.sample(result) + + # Warp resampled surface by displacement (needed for current configuration + # normals, as well volume calculation) + warped = resampled.warp_by_vector('Displacement') + + ## ------ Compute volume by warping resampled, then filling the hole ----- ## + + #warped.save(f'{output_folder}/resampled_warped_{k:03d}.vtp') + + # Fill the holes + # EITHER + # ----------------------------------------------------------------------- + # Fix the mesh (i.e. fill the holes where the inlet and outlet caps are) + #meshfix = mf.MeshFix(warped) + #meshfix.repair() + + # Convert back to pyvista polydata mesh + #lumen = meshfix.mesh + + ## PyMeshFix fills the hole, but leads to normals oriented inconsistently ## + ## The PyVista function .fill_holes() produces a bad cap mesh, but with ## + ## consistent normal vectors. The bad mesh is not important if we're only ## + ## computing volume ## + + # ---------------------------- OR ------------------------------------ + lumen = warped.fill_holes(100) # 100 is the largest size of hole to fill + + # -------------------------------------------------------------------- + + # Recompute normals, incase the normals of the cap are opposite + lumen.compute_normals(inplace=True) + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + lumen.save(f'{output_folder}/resampled_warped_and_filled_{k:03d}.vtp') + + # Add time and volume to arrays + t.append(k) + vol.append(lumen.volume) + + print(f"Iteration: {k}, Volume: {lumen.volume}") + + return (t, vol) +def calc_dVdt_struct(start_time, end_time, step, results_folder, reference_surface): + ''' + Computes the rate of change of volume of a closed (or partially closed) + surface. Intended to be used to calculate the rate of change of ventricular + volume from a struct simulation of an inflating ventricle. + + Compute the rate of change of volume by computing the velocity flux over the + reference surface, which should be the endo surface of the ventricle. + + dV/dt = int_{Gamma_t} (u . n) dA + + where Gamma_t is the reference surface in the current configuration, u is + the velocity on the reference surface, and n is the surface normal vector on + the reference surface + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + reference_surface: The absolute file path of the .vtp file containing + the undeformed surface corresponding to the deformed surface of which + we want to compute the volume. + + output_folder: The absolute file path of the folder to which we output + intermediary results of the volume calculation process. + + Returns: (t, dVdt), a tuple of lists of length number of time steps. t + contains the time step, and Q contains the rate of change of volume at that + time step. + ''' + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_struct_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + # Initialize volume array and time step lists + dVdt = [0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}_cpp.vtu")) + + # Load reference surface onto which we sample + ref_surface = pv.read(f"{reference_surface}") + + # Sample result onto ref_surface + ref= ref_surface.sample(result) + + # Warp ref surface by displacement (needed for current configuration + # normals, as well volume calculation) + warped = ref.warp_by_vector('Displacement') + + # Compute average velocity for each element + warped = warped.point_data_to_cell_data() + ref = ref.point_data_to_cell_data() + + # Recompute normals. Need this after point_data_to_cell_data(), because + # the average of point normals is not the same as a cell normal + warped.compute_normals(inplace=True) + ref.compute_normals(inplace=True) + + # Compute element areas + warped = warped.compute_cell_sizes() + ref = ref.compute_cell_sizes() + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + warped.save(f'{output_folder}/warped_{k:03d}.vtp') + ref.save(f'{output_folder}/ref_{k:03d}.vtp') + + # Compute dVdt from warped by computing velocity flux. Also, compute + # velocity flux using quantities from reference surface to compare. + # + # !! Looks like the velocity flux output in B_ST_Velocity_flux.txt is + # computed using reference surface normals and reference surface areas !! + dVdt_t = 0 + + cell_normals = warped.cell_data['Normals'] + cell_vels = warped.cell_data['Velocity'] + cell_areas = warped.cell_data['Area'] + + cell_normals_ref = ref.cell_data['Normals'] + cell_vels_ref = ref.cell_data['Velocity'] + cell_areas_ref = ref.cell_data['Area'] + + u_dot_n = np.sum(cell_vels * cell_normals, axis = 1) + dVdt_t = np.sum(u_dot_n * cell_areas) + + u_dot_n_ref = np.sum(cell_vels_ref * cell_normals_ref, axis = 1) + dVdt_t_ref = np.sum(u_dot_n_ref * cell_areas_ref) + + # Add time and volume to arrays + t.append(k) + dVdt.append(dVdt_t) + + print(f"Iteration: {k}, dVdt: {dVdt_t}, dVdt_ref: {dVdt_t_ref}") + + return (t, dVdt) + +def calc_pressure_struct(input_file, pressure_dat_file, t): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI struct simulation, in which a model of the myocardium is inflated. + + Calculates pressure on the endo surface by reading the input file and + pressure load file. The units of pressure are whatever units are used in the + pressure load file (usually Pa) + + Args: + input_file: The svFSI input file of the simulation whose results we are + now processing. + + pressure_dat_file: The file containing the pressure load information for + the simulation whose results we are now processing. + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + pressure: Pressure evaluated at list of time steps in t + """ + print('\n## Calculating pressure ##') + + # Compute pressure at each time step by reading input file and pressure file + with open(input_file) as f: + lines = f.readlines() + for line in lines: + if 'Number of time steps' in line: + num_time_steps = int(line.split(':')[-1]) + + with open(pressure_dat_file) as f: + f.readline() + line = f.readline() + vals = line.split() + t0 = float(vals[0]) + p0 = float(vals[1]) + line = f.readline() + vals = line.split() + t1 = float(vals[0]) + p1 = float(vals[1]) + + pressure = np.linspace(p0, p1, num_time_steps+1) + + # Sample pressure onto t + pressure = pressure[t] + + return pressure + +def calc_pressure_struct_genBC(genBC_out_file, col, t): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI struct simulation with genBC, in which a model of the myocardium is + inflated. + + Calculates pressure on the endo surface by reading the genBC output file, + usually called AllData. The units of pressure are whatever units are used in + the genBC input file. + + Args: + genBC_out_file: The file containing the pressure load information for + the simulation whose results we are now processing. Usually AllData + + col: The column in AllData that contains the pressure at each time step + col = 1 refers to the first column + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + pressure: Pressure evaluated at list of time steps in t + """ + print('\n## Calculating pressure ##') + + # Read pressure from AllData + pressure = [0] + with open(genBC_out_file) as f: + lines = f.readlines() + for line in lines: + vals = line.split() + pressure.append(float(vals[col-1])) + + pressure = np.array(pressure) + + # Sample pressure onto t + pressure = pressure[t] + + return pressure +def calc_volume_struct_genBC(genBC_out_file, col, t): + """ + Calculate the accumulated lumen volume at each time step from the results of + an svFSI struct simulation with genBC, in which a model of the myocardium is + inflated. + + Calculates accumulated volumeby reading the genBC output file, + usually called AllData. In genBC, you must define an additional unknown + x(i) such that f(i) = -Q(1). This will compute x(i) as the integral of the + flow rate, i.e. the accumulated volume + + Args: + genBC_out_file: The file containing the pressure load information for + the simulation whose results we are now processing. Usually AllData + + col: The column in AllData that contains the pressure at each time step + col = 1 refers to the first column + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + volume: Volume evaluated at list of time steps in t + """ + print('\n## Calculating volume from AllData ##') + + # Read pressure from AllData + volume = [0] + with open(genBC_out_file) as f: + lines = f.readlines() + for line in lines: + vals = line.split() + volume.append(float(vals[col-1])) + + volume = np.array(volume) + + # Sample pressure onto t + volume = volume[t] + + return volume + +def calc_volume_FSI(volume_file): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the volume by reading the output file V_FS_Volume_average.txt. This + is produced if the following is included in svFSI.inp + + Output: Volume_integral { + Volume: t + } + + This gives the volume of the fluid domain at each time step. + + The units of volume are whatever units used in svFSI.inp, cubed. For example, + if units of length in svFSI.inp are meters, then the volume calculated + here is cubic meters. Note that the mesh scale factor affects these units. If + the .vtu file has units of microns, but we apply a mesh scale factor of 1e-6 + to convert to meters in svFSI.inp, then the volume output here will be in meters^3 + + Args: + volume_file: Path to file V_FS_Volume_average.txt from simulation you + want to process + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + print('\n## Calculating volume ##') + + # Initialize time and volume arrays + t = [] + vol = [] + + # Read volume data from V_FS_Volume_average.txt + with open(volume_file) as f: + header = f.readline() # Read header + initial_volumes = f.readline() # Read domain initial volumes + f.readline() # Skip blank line + + # Read volume data + # (For some reason, the first volume is the initial volume, even though + # this corresponds to after the first time step) + lines = f.readlines() + k = 0 + for line in lines: + + # Add time iteration to t array + t.append(k) + k += 1 + + # Add fluid domain (Domain 0) volume to volume array + V_t = float(line.split()[0]) + vol.append(V_t) + + return (t, vol) + + +def calc_volume_FSI_2(start_time, end_time, step, results_folder): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the volume in the following steps + 1) Threshold to extract fluid domain + 2) Warp the fluid domain by the Displacement + 3) Calculate the volume of the warped fluid domain + + The units of volume are whatever units used in .vtu files, cubed. For example, + if units of length in the .vtu files are microns, then the volume calculated + here is cubic microns. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_FSI_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating volumes ##') + + # Initialize volume array and time step lists + vol = [0.0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Threshold to fluid domain + fluid = result.threshold(value = (1,1), scalars = 'Domain_ID') + + # Warp fluid domain by displacement (yields current fluid domain) + warped = fluid.warp_by_vector('Displacement') + + # Save warped fluid domain (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + #warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtp') + warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtu') + + + + # Add time and volume to arrays + t.append(k) + vol.append(warped.volume) + + print(f"Iteration: {k}, Volume: {warped.volume}") + + return (t, vol) + + +def calc_pressure_FSI(pressure_file): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the pressure by reading the output file V_FS_Pressure_average.txt. This + is produced if the following is included in svFSI.inp + + Output: Volume_integral { + Pressure: t + } + This gives the volume averaged fluid pressure at each time step. + + + The units of pressure are deduced by whatever units system you are using in + the input file. For SI (m, kg, s), the units of pressure are Pa. + + Args: + pressure_file: Path to file V_FS_Pressure_average.txt from simulation you + want to process + + Returns: (t, pressure), a tuple of lists of length number of time steps. t + contains the time step, and pressure contains the pressure at that time step. + """ + print('\n## Calculating pressure ##') + + # Initialize time and pressure arrays + t = [0] + pressure = [0] + + # Read pressure data from V_FS_Pressure_average.txt + with open(pressure_file) as f: + header = f.readline() # Read header + initial_volumes = f.readline() # Read domain initial volumes + f.readline() # Skip blank line + + # Read pressure data + lines = f.readlines() + for line in lines: + + # Add next time iteration to t array + t.append(t[-1] + 1) + + # Add fluid domain (Domain 0) pressure to pressure array + p_t = float(line.split()[0]) + pressure.append(p_t) + + return (t, pressure) + +def calc_pressure_FSI_2(start_time, end_time, step, results_folder): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the pressure in the following steps + 1) Threshold to extract fluid domain + 2) Warp the fluid domain by the Displacement + 3) Compute the volume of each cell + 4) Compute the pressure in each cell + 5) Compute the volume-averaged pressure + + The units of pressure are deduced by whatever units system you are using in + the input file. For SI (m, kg, s), the units of pressure are Pa. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + Returns: (t, pressure), a tuple of lists of length number of time steps. t + contains the time step, and pressure contains the pressure at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_pressure_FSI_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating pressures ##') + + # Initialize volume array and time step lists + pressure = [0.0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Threshold to fluid domain + fluid = result.threshold(value = (1,1), scalars = 'Domain_ID') + + # Warp fluid domain by displacement (yields current fluid domain) + warped = fluid.warp_by_vector('Displacement') + + # Compute volumes and areas + sized = warped.compute_cell_sizes() + + + # Compute pressure in each cell from pressure at each point + sized = sized.point_data_to_cell_data() + + + + # Save processed fluid domain (to check geometry, normals, and values) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + #warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtp') + sized.save(f'{output_folder}/fluid_processed_{k:03d}.vtu') + + # Grab volumes for all elements in mesh + elem_volumes = sized.cell_data['Volume'] + # Grab pressure for all elements in mesh + elem_pressures = sized.cell_data['Pressure'] + + # Compute volume averaged pressure by looping over fluid domain elements + p_avg = 0 + for i in range(len(elem_volumes)): + p_avg += elem_pressures[i] * elem_volumes[i] + p_avg = p_avg / sum(elem_volumes) + + + # Add time and pressure to arrays + t.append(k) + pressure.append(p_avg) + + print(f"Iteration: {k}, Pressure: {p_avg}") + + return (t, pressure) + +#def avi2mp4(input_movie_file, output_movie_file): +# ff = ffmpy.FFmpeg( +# inputs={input_movie_file: None}, +# outputs={output_movie_file: '-pix_fmt yuv420p -y'} +# # flag -pix_fmt yuv420p needed so that movie is compatible with QuickTime +# # flag -y needed to overwrite mp4 file if it already exists +# ) +# ff.run() \ No newline at end of file diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/process_results.py b/tests/cases/struct/LV_NeoHookean_passive_genBC/process_results.py new file mode 100644 index 00000000..ad94716d --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/process_results.py @@ -0,0 +1,154 @@ +# Python script to calculate lumen volume and fluid pressure on endo surface from +# results of svFSI struct +# +# Requires two meshes +# - svFSI result.vtu (with Displacement array) +# - reference surface file (polygonal mesh), representing the undeformed surface +# corresponding to the deformed surface of which we want to compute the volume +# +# We calculate the volume in the following steps +# 1) Sample the result.vtu file onto the reference surface +# 2) Warp the samples surface by the Displacement +# 3) Flat fill any holes in the warped surface +# 4) Calculate the volume of the warped and filled surface +# +# We calculate the fluid pressure on the endo surface by reading the input file +# and pressure load file. + +import os # for checking and creating directories, and loading modules +import matplotlib.pyplot as plt +import numpy as np +import sys + +# Import custom functions useful for post-processing +from post_processing_functions import * + +## -------- PARAMETERS TO CHANGE -------------------- ## + +# Simulation folder file path +sim_folder = os.path.dirname(os.path.realpath(__file__)) + +# Optionally read sim_folder as command line argument +if len(sys.argv) == 2: + print(sys.argv) + sim_folder = os.path.abspath(sys.argv[1]) + print(sim_folder) + + +# Undeformed surface, corresponding to deformed surface on result.vtu of which we want to compute the volume +reference_surface = os.path.join(sim_folder, 'mesh/mesh-surfaces/endo.vtp') + +# Undeformed vtu file (the original volume mesh) +reference_volume = os.path.join(sim_folder, 'mesh/mesh-complete.mesh.vtu') + +# svFSI results folder, containing results.vtu +#results_folder = os.path.join(sim_folder, 'results_svfsi/') +results_folder = os.path.join(sim_folder, '4-procs/') + +# File containing genBC output +#alldata_file = 'AllData_svfsi_vvedula22' +alldata_file = 'AllData' + +# Input file path +input_file = os.path.join(sim_folder, 'svFSI.inp') + +# Pressure file (unused) +pressure_dat_file = os.path.join(sim_folder, 'pressure.dat') + + +## -------------------- END PARAMETERS TO CHANGE ------------------------ ## + + +# Automatically determine the start time, end time, and step size based on all +# results file in results_folder +print(results_folder) +(start_time, end_time, step) = get_start_end_step(results_folder) +# Option to manually set start time, end time, and time step of results files to process +#start_time = 5 +#end_time = 85 +#step = 5 + +# Compute lumen volume from simulation results +(t_3D, vol_3D) = calc_volume_struct(start_time, end_time, step, results_folder, reference_surface) +vol_3D_cm3 = np.array(vol_3D) * (100)**3 # cm^3 + +# Compute lumen dVdt from simulation results +(t_dVdt_3D, dVdt_3D) = calc_dVdt_struct(start_time, end_time, step, results_folder, reference_surface) +# Convert volume to cm^3/s +dVdt_3D = np.array(dVdt_3D) # cm/s * m^2 +dVdt_3D_cm3 = dVdt_3D * (100)**2 # cm^3/s + +# Compute lumen volume from AllData file +vol_0D_cm3 = calc_volume_struct_genBC(os.path.join(sim_folder, alldata_file), 2, t_3D) +# Add on initial volume +vol_0D_cm3 += vol_3D_cm3[0] # cm^3 + +# Compute flow rate = dVdt from AllData file +dVdt_0D_cm3 = calc_volume_struct_genBC(os.path.join(sim_folder, alldata_file), 4, t_dVdt_3D) + + +# Compute lumen pressure at iterations in t +pressure = calc_pressure_struct_genBC(os.path.join(sim_folder, "AllData"), 1, t_3D) # dynes/cm^2 + +# Convert pressure from dynes/cm^2 to mmHg +#pressure_mmHg = pressure * 0.000750062 + + + +# Combine time step, pressure, and volume into one array +#PV = np.column_stack((t_3D, pressure, vol_3D_cm3)) + +#print('\n## Outputing pressure-volume data and plot ##') + +# Write pressure and volume to file +#np.savetxt(results_folder + '/../' + "pv.txt", PV, header = 'Timestep Pressure[mmHg] Volume[m^3]') + + +# Plot pressure vs. volume +fig, ax = plt.subplots() +ax.plot(vol_3D_cm3, pressure, linewidth=2.0, marker = 'o') +ax.set_xlabel('Volume [cm^3]') +ax.set_ylabel('Pressure [dyne/cm^2]') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'pv_plot')) + +# Plot dVdt vs. time +#fig, ax = plt.subplots() +#ax.plot(t_dVdt, dVdt_m3, linewidth=2.0, marker = 'o') +#ax.set_xlabel('Time') +#ax.set_ylabel('dVdt') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +#plt.savefig(os.path.join(sim_folder, 'dVdt_plot')) + +# Plot dVdt vs. Pressure +#fig, ax = plt.subplots() +#ax.plot(dVdt_m3, pressure_mmHg, linewidth=2.0, marker = 'o') +#ax.set_xlabel('dVdt') +#ax.set_ylabel('Pressure') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +#plt.savefig(os.path.join(sim_folder, 'dVdtvsP_plot')) + +# Plot 3D and 0D dVdt +fig, ax = plt.subplots() +ax.plot(dVdt_3D_cm3, label = 'dVdt_3D') +ax.plot(dVdt_0D_cm3, label = 'dVdt_0D', linestyle = '--') +ax.set_xlabel('Timestep') +ax.set_ylabel('dVdt (cm^3/s)') +ax.legend() +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'dVdt3D_vs_dVdt0D')) + +# Plot 3D and 0D volume +fig, ax = plt.subplots() +ax.plot(vol_3D_cm3, label = 'V_3D') +ax.plot(vol_0D_cm3, label = 'V_0D', linestyle = '--') +ax.set_xlabel('Timestep') +ax.set_ylabel('Volume (cm^3)') +ax.legend() +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'V3D_vs_V0D')) diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/pv_plot.png b/tests/cases/struct/LV_NeoHookean_passive_genBC/pv_plot.png new file mode 100644 index 00000000..f2aab43a --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/pv_plot.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d67aa33c8e9584754ab736140f84ab04ba831aeaf096dbdcd639d7faf1181841 +size 19560 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/result_003.vtu b/tests/cases/struct/LV_NeoHookean_passive_genBC/result_003.vtu new file mode 100644 index 00000000..96062588 --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/result_003.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3e7afe2fc92e79c64d89de3aaecdd55a2c2f05b052321e2276a059b4d9c50ba1 +size 248337 diff --git a/tests/cases/struct/LV_NeoHookean_passive_genBC/svFSI.xml b/tests/cases/struct/LV_NeoHookean_passive_genBC/svFSI.xml new file mode 100644 index 00000000..d8a0a9ee --- /dev/null +++ b/tests/cases/struct/LV_NeoHookean_passive_genBC/svFSI.xml @@ -0,0 +1,103 @@ + + + + + 0 + 3 + 3 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 10 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + 100.0 + + + + + + + true + 3 + 10 + 1e-9 + + 1.0 + 1.0e5 + 0.483333 + + + ST91 + + 0.0 + + + + + true + true + true + true + true + true + true + true + + + + FSILS + 1e-9 + 1000 + 50 + + + + + genBC_svFSIplus/genBC.exe + + + + Dirichlet + 0.0 + + + + Neu + Coupled + true + + + + + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/README.md b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/README.md new file mode 100644 index 00000000..dfcc0365 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/README.md @@ -0,0 +1,29 @@ +This test case simulates an idealized left ventricle (LV) with a NeoHookean material model +coupled to a lumped-parameter network (LPN), implemented in genBC. The LPN consists of a large pressure source and large resistor, which together produce an approximately constant flowrate into +the LV. This inflates the LV at an approximately constant rate of change of volume. + +Before running svFSIplus, genBC_svFSIplus must be compiled. Navigate to genBC_svFSIplus +and run `make clean` then `make`. + +The results can be post-processed by running `process_results.py`. If run for 100 +timesteps, the results should be the following. + +![V3D vs. V0D](V3D_vs_V0D.png) + +*Plot of volume vs. time for this simulation. V_3D is computed from results.vtu +by extracting the endocardial surface, warping it according to displacement at +each timestep, capping it, then measuring the enclosed volume. V_0D is read +from AllData (produced by genBC).* + + +![dVdt3D vs. dVdt0D](dVdt3D_vs_dVdt0D.png) + +*Plot of flowrate (dVdt) vs. time for this simulation. dVdt_3D is computed from results.vtu +by extracting the endocardial surface, warping it according to displacement at +each timestep, then computing the velocity flux integral. +V_0D is computed from AllData (produced by genBC). As seen in the plot, we +ramp the flowrate from 0 to 100 cm^3/s over 10 timesteps.* + +![Pressure vs. Volume](pv_plot.png) + +*Plot of pressure vs. volume as the LV is inflated.* \ No newline at end of file diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png new file mode 100644 index 00000000..68cf7cab --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/V3D_vs_V0D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd353855d1499573894e7bdb86dd2cfd99a9a80e79324a1433532ec15e2fa47f +size 23958 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png new file mode 100644 index 00000000..f8cba954 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/dVdt3D_vs_dVdt0D.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:27fa8d4244e55ece061e29f22afc9e207003190a66efbf8a86e86ede53688b58 +size 22082 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile new file mode 100644 index 00000000..384f82a2 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile @@ -0,0 +1,79 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the Makefile to build GenBC. +# +#-------------------------------------------------------------------- + +# HERE COMES THE DEFINITION + +include Makefile.in + +INCLUDES = -I./include + +MYFUN = Modules.f \ + USER.f \ + GenBC.f + +GenBC_EXE = genBC.exe + +DSYM_DIR = +ifeq ($(debug),1) + ifeq ($(OS),Darwin) + DSYM_DIR = $(GenBC_EXE).dSYM + endif +endif +############################################################################# +# AND HERE ARE THE RULES + +SRC = $(patsubst %,$(SRC_DIR)/%,$(MYFUN)) +OBJ = $(patsubst %.f,$(OBJ_DIR)/%.o,$(MYFUN)) + +#.PHONY: $(TEST) +#$(TEST): $(TEST:.exe=.f) $(GenBC_EXE) +# $(FORTRAN) $< $(FFLAGS) $(INCLUDES) -o $@ + +.PHONY: $(GenBC_EXE) +$(GenBC_EXE): $(OBJ) + $(FORTRAN) $(OBJ) $(FFLAGS) -o $@ + +$(OBJ): | $(OBJ_DIR) + +$(OBJ_DIR): + mkdir -p $@ + +$(OBJ_DIR)/%.o : $(SRC_DIR)/%.f + $(FORTRAN) $(FFLAGS) $(INCLUDES) -c $< -o $@ + +clean: + rm -r -f $(OBJ_DIR) $(GenBC_EXE) $(TEST) $(DSYM_DIR) diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile.in b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile.in new file mode 100644 index 00000000..7a434ad0 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/Makefile.in @@ -0,0 +1,165 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the definitions for building process. +# +#-------------------------------------------------------------------- + + +.SUFFIXES: .f .o + +OS:= $(shell uname) + +AR = ar rv +RANLIB = ar -ts + +CFLAGS = -O3 -DAdd_ +CXXFLAGS = -O3 -std=c++11 +FFLAGS = -O3 -cpp +FCFLAGS = -lstdc++ -cpp +#FCFLAGS = -lc++ -cpp ## This may work if lstdc++ is not found +OBJ_DIR = obj +SRC_DIR = src + +MAKE_GUI = 0 + +debug = 0 +ifeq ($(debug),1) + CFLAGS = -O0 -DAdd_ + CXXFLAGS = -O0 -std=c++11 + FFLAGS = -O0 -cpp +endif + +seq = 1 +ifeq ($(seq),1) + CC = gcc + CXX = g++ + FORTRAN = gfortran + CFLAGS += -DSEQ +else + CC = mpicc + CXX = mpicxx + FORTRAN = mpif90 +endif + +# ---------------------------------------------------------------- +# Normally you would not need to change any line beyond this point +# ---------------------------------------------------------------- + +# Here I am finding the compiler group +ifeq ($(seq),1) + F_COMP = $(FORTRAN) +else + F_COMP = $(firstword $(shell $(FORTRAN) -show)) +endif +ifeq ($(F_COMP),gfortran) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),gcc) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),g77) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),f95) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),ifort) + COMP_GRP = intel +endif +ifeq ($(F_COMP),pgf90) + COMP_GRP = pgi +endif +ifeq ($(F_COMP),pgf77) + COMP_GRP = pgi +endif + +ifeq ($(OS),Darwin) + COMP_GRP = gnu +endif + +# If profiling is requested +prof = 0 +ifeq ($(prof),1) + FFLAGS += -pg -g +endif + +# If debuging is requested +ifeq ($(debug),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -g -Wall -Wconversion -Wline-truncation -pedantic -fimplicit-none -fbacktrace -fbounds-check -p -fcheck=all #-ffpe-trap=invalid,zero,overflow,underflow + CXXFLAGS += -g -Wall -pedantic -fbounds-check + CFLAGS += -g -Wall -pedantic -fbounds-check + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -g -traceback -check all -check bounds -check uninit -ftrapuv -debug all -fpe0 + CFLAGS += -g -traceback -check-uninit -fpe0 + CXXFLAGS += -g -traceback -check-uninit -fpe0 + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += # You need to add debuging flag for pgi compiler here + endif + OBJ_DIR = dbg +endif + +# If openMP parallelization is requested +mp = 0 +ifeq ($(mp),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -fopenmp + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -openmp + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += -mp + endif +endif + +# To make directories a bit cleaner with intel compiler +ifeq ($(COMP_GRP),intel) + FFLAGS += -module $(OBJ_DIR) +endif +ifeq ($(COMP_GRP),gnu) + FFLAGS += -J $(OBJ_DIR) + CFLAGS += -J $(OBJ_DIR) + CXXFLAGS += -J $(OBJ_DIR) +endif + +#LAPACK Library +ifeq ($(COMP_GRP),gnu) + LAPACK_INC = + LAPACK_LIB = + LAPACK_LFLAGS = -llapack +endif + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f new file mode 100644 index 00000000..991bbb50 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/initial_values_final.f @@ -0,0 +1,4 @@ +c Value of your unknown at time equal to zero (This is going to be used +c ONLY when you are initiating simulation) + + tZeroX = (/ 0.00D+00 /) diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f new file mode 100644 index 00000000..af40e47d --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/include/parameters_final.f @@ -0,0 +1,15 @@ +! Specifing the constants + REAL(KIND=8), PARAMETER :: ! Order of surfaces: 3~8 + & Rmyo = 4.51000000D-02, + & Lv = 2.74000000D-05, + & Kv = 8.16000000D-04, + & Kao = 8.59000000D-05, + & pi = 3.14159265D+00, + & Afill = 6D+01, + & tfill = 5D+00 +! The parameters which might change through the code + REAL(KIND=8) Clung, Rlung + Clung= 2.56000000D-02 + Rlung= 7.29000000D-01 + + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f new file mode 100644 index 00000000..3be431da --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/GenBC.f @@ -0,0 +1,230 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c Here data are received from Phasta and it setup the data required for +c integration of ODE's inside FINDX + PROGRAM GenBC + USE COM + USE, INTRINSIC :: IEEE_ARITHMETIC + IMPLICIT NONE + + INTEGER i, n, nTimeStep + REAL(KIND=8) t, dt, tFinal, pMin, temp + CHARACTER(LEN=2048) string + CHARACTER(LEN=32) sTmp + CHARACTER flag + + REAL(KIND=8), ALLOCATABLE :: Qi(:), Qf(:), Pi(:), Pf(:), X(:), + 2 Xo(:), f(:,:) + +c flag = I : Initializing +c flag = T : Iteration loop +c flag = L : Last iteration +c flag = D : Derivative +c +c Here is the period of time that you should integrate for +c each of these flags: +c +c Flags I D&T&L +c ^ ^ +c 3D code time step: N.............................................N+1 +c 0D code time step: 1......................................nTimeStep+1 +c Flowrates: Qi............................ ................Qf +c Time, t: tInitial..............................tInitial+tFinal + + + CALL INITIALIZE(nTimeStep) + +c******************************************************************** +c Block for reading the data from phasta + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED'); + READ (1) flag + READ (1) tFinal + READ (1) i + IF (i .NE. nDirichletSrfs) THEN + PRINT *, 'Error: Number of Dirichlet Surfaces from Phasta is:', + 2 i + PRINT *, 'While nDirichletSrfs is equal to:', nDirichletSrfs + PRINT * + PRINT *, 'Number of Dirichlet surfaces defined in solver.inp', + 2 ' should match with nDirichletSrfs defined in USER.f' + STOP + END IF + READ (1) i + IF (i .NE. nNeumannSrfs) THEN + PRINT *, 'Error: Number of Neumann Surfaces from Phasta is:', + 2 i + PRINT *, 'While nNeumannSrfs is equal to:', nNeumannSrfs + PRINT * + PRINT *, 'Number of Neumann surfaces defined in solver.inp', + 2 ' should match with nNeumannSrfs defined in USER.f' + STOP + END IF + IF (nDirichletSrfs .GT. 0) THEN + IF (qCorr .OR. pCorr) THEN + PRINT *, 'You should only use P/Q correction when all', + 2 ' the surfaces are Neumann surfaces' + STOP + END IF + END IF + + ALLOCATE (Pi(nDirichletSrfs), Pf(nDirichletSrfs), + 2 PDirichlet(nDirichletSrfs,4)) + ALLOCATE (Qi(nNeumannSrfs), Qf(nNeumannSrfs), + 2 QNeumann(nNeumannSrfs,4)) + + DO i=1, nDirichletSrfs + READ(1) Pi(i), Pf(i) + END DO + Pi = Pi/pConv + Pf = Pf/pConv + + DO i=1,nNeumannSrfs + READ (1) Qi(i), Qf(i) +c Print to console and file pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Flowrates to genBC: ', Qi(i), Qf(i) + END DO + CLOSE(1) + Qi = Qi/qConv + Qf = Qf/qConv + +c******************************************************************** +c Block for initializing the unknowns + + IF (flag .EQ. 'I') THEN + nTimeStep = 0 + ELSE + dt = tFinal/REAL(nTimeStep,8) + END IF + ALLOCATE (X(nUnknowns), Xo(nUnknowns), f(nUnknowns,4), + 2 offset(nUnknowns), Xprint(nXprint)) + Xprint = 0D0 + + OPEN (1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED') + READ (1) t + DO i=1,nUnknowns + READ (1) Xo(i) + END DO + CLOSE (1) + +c******************************************************************** +c Setting up the system of equations + offset = 0D0 + DO n=1, nTimeStep + DO i=1, 4 + temp = (REAL(n-1,8) + REAL(i-1,8)/3D0)/REAL(nTimeStep,8) + + QNeumann(:,i) = Qi + (Qf-Qi)*temp + PDirichlet(:,i) = Pi + (Pf-Pi)*temp + + IF (qCorr) THEN + temp = SUM(QNeumann(:,i))/REAL(nNeumannSrfs,8) + QNeumann(:,i) = QNeumann(:,i) - temp + END IF + END DO + + X = Xo + CALL FINDF (t, X, f(:,1), QNeumann(:,1), + 2 PDirichlet(:,1)) + X = Xo + dt*f(:,1)/3D0 + + CALL FINDF (t+dt/3D0, X, f(:,2), QNeumann(:,2), + 2 PDirichlet(:,2)) + X = Xo - dt*f(:,1)/3D0 + dt*f(:,2) + + CALL FINDF (t+dt*2D0/3D0, X, f(:,3), QNeumann(:,3), + 2 PDirichlet(:,3)) + X = Xo + dt*f(:,1) - dt*f(:,2) + dt*f(:,3) + + CALL FINDF (t+dt, X, f(:,4), QNeumann(:,4), + 2 PDirichlet(:,4)) + + f(:,1) = (f(:,1) + 3D0*f(:,2) + 3D0*f(:,3) + f(:,4))/8D0 + Xo = Xo + dt*f(:,1) + t = t + dt + END DO + +c******************************************************************** +c Time to write the results + X = Xo + IF (pCorr .AND. flag.NE.'D') THEN + pMin = X(srfToXPtr(1)) + DO i=2, nNeumannSrfs + IF (X(srfToXPtr(i)) .LT. pMin) THEN + pMin = X(srfToXPtr(i)) + END IF + END DO + ELSE + pMin = 0D0 + END IF + +c Writing nDirichlet flowrates here + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED') + DO i=1, nDirichletSrfs + IF(IEEE_IS_NAN(X(srfToXdPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF + WRITE (1) X(srfToXdPtr(i))*qConv + END DO + +c Writing nNeumannSrfs pressures here + DO i=1, nNeumannSrfs + IF(IEEE_IS_NAN(X(srfToXPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF +c Print out pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Pressure from genBC: ', offset(srfToXPtr(i)) + +c Here, saving offset value to X itself. Necessary for keeping track of +c pressure from one time step to the next + X(srfToXPtr(i)) = offset(srfToXPtr(i)) + +c WRITE (1) (X(srfToXPtr(i)) - pMin +c 2 + offset(srfToXPtr(i)))*pConv + WRITE (1) (X(srfToXPtr(i)))*pConv + END DO + CLOSE(1) + + IF (flag .EQ. 'L') THEN + OPEN(1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED'); + WRITE (1) t + DO i=1,nUnknowns + WRITE (1) X(i) + END DO + CLOSE(1) + +c PRINT *,'Before AllData' + + OPEN(1, FILE='AllData', STATUS='UNKNOWN', ACCESS='APPEND'); + string = '' + DO i=1, nUnknowns + WRITE (sTmp,"(ES14.6E2)") X(i) + string = TRIM(string)//sTmp + END DO + DO i=1, nXprint + WRITE (sTmp,"(ES14.6E2)") Xprint(i) + string = TRIM(string)//sTmp + END DO + WRITE (1,"(A)") TRIM(string) + CLOSE(1) + END IF + + DEALLOCATE (Pi) + DEALLOCATE (Pf) + DEALLOCATE (PDirichlet) + DEALLOCATE (Qi) + DEALLOCATE (Qf) + DEALLOCATE (QNeumann) + DEALLOCATE (X) + DEALLOCATE (Xo) + DEALLOCATE (f) + DEALLOCATE (srfToXdPtr) + DEALLOCATE (srfToXPtr) + DEALLOCATE (offset) + DEALLOCATE (Xprint) + + END PROGRAM GenBC diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f new file mode 100644 index 00000000..6e472f3b --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/Modules.f @@ -0,0 +1,17 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + + MODULE COM + + LOGICAL pCorr, qCorr + + INTEGER nUnknowns, nDirichletSrfs, nNeumannSrfs, nXprint + + REAL(KIND=8) pConv, qConv + + INTEGER, ALLOCATABLE :: srfToXPtr(:), srfToXdPtr(:) + + REAL(KIND=8), ALLOCATABLE :: QNeumann(:,:), PDirichlet(:,:), + 2 offset(:), Xprint(:) + + END MODULE COM diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/USER.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/USER.f new file mode 100644 index 00000000..e48f7c19 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC/src/USER.f @@ -0,0 +1,143 @@ +c MUSC 2 Immediate postop +c Ethan Kung keo@ucsd.edu + +c Created by Mahdi Esmaily Moghadam 12-01-2010 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c This subroutine initializes the parameters, you need to read the +c comments and specify them carefuly +c-------------------------------------------------------------------- +c This is an example for RCR boundary condition with parameters +c Rd, R, and C which are distal, and proximal resistance and capacitor. + + SUBROUTINE INITIALIZE(nTimeStep) + USE COM + IMPLICIT NONE + INTENT(OUT) nTimeStep + + LOGICAl ierr + INTEGER i, nTimeStep + REAL(KIND=8), ALLOCATABLE :: tZeroX(:) +c +c******************************************************************** +c For instance if pressure in 3D solver is in cgs and here mmHg +c pConv=1334=1.334D3, also if flowrate in 3D solver is mm^3/s and +c here is mL/s qConv=1000=1D3. In the case both solver are using same +c unites you can set these two conversion coefficients to 1D0 + pConv = 1D0 + qConv = 1D0 + +c Only when all the surfaces of you model are coupled with NeumannSrfs +c you may set this to .TRUE. + pCorr = .FALSE. + qCorr = .FALSE. + +c******************************************************************** +c Block for your inputs + +c These two value should match to that of defined in solver.inp + nDirichletSrfs = 0 + nNeumannSrfs = 1 +c Number of unknowns that you need inside your lumped parameter network + nUnknowns = 2 +c Number of time step between N and N+alpha, for higher values this code +c would be more accurate and more costly + nTimeStep = 10 + +c Number of parameters to be printed in AllData file (the first +c nUnknowns columns are the "X" values) + nXprint = 2 + +c-------------------------------------------------------------------- +c You don't need to change this part ! + ALLOCATE (tZeroX(nUnknowns), srfToXdPtr(nDirichletSrfs)) ! + ALLOCATE (srfToXPtr(nNeumannSrfs)) ! + tZeroX = 0D0 +c-------------------------------------------------------------------- + + INCLUDE "initial_values_final.f" + +c-------------------------------------------------------------------- +c You don't need to change this part ! + INQUIRE (FILE='InitialData', EXIST=ierr) ! + IF (.NOT.ierr) THEN ! +c PRINT *, 'Initializing unknowns in LPM' ! + OPEN (1, FILE='InitialData',STATUS='NEW',FORM='UNFORMATTED')! + WRITE (1) 0D0 ! + DO i=1, nUnknowns ! + WRITE (1) tZeroX(i) ! + END DO ! + CLOSE(1) ! + END IF ! +c-------------------------------------------------------------------- + +c Surface to X pointer: this defines which Unknown corresponds to which +c suface in "List of Neumann Surfaces" inside solver.inp +c For example, if you have "List of Neumann Surfaces= 2 8 4 ...." +c and you set "srfToXPtr = (/5,3,9,.../)" +C this means X(5) corresponds to surface 2, X(3) <-> surface 8, +c and X(9) <-> surface 4 +c Also, Q(1) corresponds to the first element in srfToXPtr, in this +c example, Q(1)<=>X(5), Q(2)<=>X(3) + srfToXPtr = (/1/) +c srfToXdPtr = (/1/) + + END SUBROUTINE INITIALIZE + +c#################################################################### +c Here you should find the f_i=dx_i/dt, based on the following parameters: +c current x_i: x(i) +c Current time: t +c Flowrates from 3D code: Q(i) +c Pressure from Dirichlet faces: P(i) + + SUBROUTINE FINDF(t, x, f, Q, P) + USE COM + IMPLICIT NONE + INTENT(IN) t, Q + INTENT(OUT) f + + REAL(KIND=8) t, x(nUnknowns), f(nUnknowns), Q(nNeumannSrfs), + 2 P(nDirichletSrfs) + +! These are the dumy variables + REAL(KIND=8) Pivc,Ppwc,Pcon,Qout, Phigh, Pmax, tramp, R + + INCLUDE "parameters_final.f" + +!******************************************************************** +! This is the only case that x is directly manipulated, ugly, but the +! only way!!! + +! Implement a resistance filling network so that the ventricle pressure +! Pv = Phigh - Q*R +! This fills the sphere at an approximately constant rate + +! x(1): Pressure applied to endo surface. Actually, we will assign the pressure +! by setting the offset instead +! x(2): Accumulated volume during inflation. Not necessary, but for insight + +! Define a pressure ramp for Phigh (1e7 dyne/cm^2 in 0.11 second, then constant) + Pmax = 1.0D+7 + tramp = 0.1 + IF (t .LT. tramp) THEN + Phigh = Pmax * t / tramp + ELSE + Phigh = Pmax + ENDIF + R = 1.0D+5 + +! The main body of equations + f(1) = 0 ! The endo surface pressure (we will set it with offset) + f(2) = -Q(1) ! Integrate the flow rate to compute accumulated volume + +! Define offset pressure. + offset(1) = Phigh + Q(1) * R + +c Assign the additional parameters to be printed + Xprint(1) = offset(1) ! This is the actual value of pressure for this sim + Xprint(2) = Q(1) ! The flowrate/dVdt + RETURN + END SUBROUTINE FINDF + + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile new file mode 100644 index 00000000..384f82a2 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile @@ -0,0 +1,79 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the Makefile to build GenBC. +# +#-------------------------------------------------------------------- + +# HERE COMES THE DEFINITION + +include Makefile.in + +INCLUDES = -I./include + +MYFUN = Modules.f \ + USER.f \ + GenBC.f + +GenBC_EXE = genBC.exe + +DSYM_DIR = +ifeq ($(debug),1) + ifeq ($(OS),Darwin) + DSYM_DIR = $(GenBC_EXE).dSYM + endif +endif +############################################################################# +# AND HERE ARE THE RULES + +SRC = $(patsubst %,$(SRC_DIR)/%,$(MYFUN)) +OBJ = $(patsubst %.f,$(OBJ_DIR)/%.o,$(MYFUN)) + +#.PHONY: $(TEST) +#$(TEST): $(TEST:.exe=.f) $(GenBC_EXE) +# $(FORTRAN) $< $(FFLAGS) $(INCLUDES) -o $@ + +.PHONY: $(GenBC_EXE) +$(GenBC_EXE): $(OBJ) + $(FORTRAN) $(OBJ) $(FFLAGS) -o $@ + +$(OBJ): | $(OBJ_DIR) + +$(OBJ_DIR): + mkdir -p $@ + +$(OBJ_DIR)/%.o : $(SRC_DIR)/%.f + $(FORTRAN) $(FFLAGS) $(INCLUDES) -c $< -o $@ + +clean: + rm -r -f $(OBJ_DIR) $(GenBC_EXE) $(TEST) $(DSYM_DIR) diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in new file mode 100644 index 00000000..aa00a611 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/Makefile.in @@ -0,0 +1,165 @@ +# +# Copyright (c) Stanford University, The Regents of the University of +# California, and others. +# +# All Rights Reserved. +# +# See Copyright-SimVascular.txt for additional details. +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files (the +# "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to +# permit persons to whom the Software is furnished to do so, subject +# to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +#-------------------------------------------------------------------- +# +# This is the definitions for building process. +# +#-------------------------------------------------------------------- + + +.SUFFIXES: .f .o + +OS:= $(shell uname) + +AR = ar rv +RANLIB = ar -ts + +CFLAGS = -O3 -DAdd_ +CXXFLAGS = -O3 -std=c++11 +FFLAGS = -O3 -cpp +FCFLAGS = -lstdc++ -cpp +#FCFLAGS = -lc++ -cpp ## This may work if lstdc++ is not found +OBJ_DIR = obj +SRC_DIR = src + +MAKE_GUI = 0 + +debug = 0 +ifeq ($(debug),1) + CFLAGS = -O0 -DAdd_ + CXXFLAGS = -O0 -std=c++11 + FFLAGS = -O0 -cpp +endif + +seq = 0 +ifeq ($(seq),1) + CC = gcc + CXX = g++ + FORTRAN = gfortran + CFLAGS += -DSEQ +else + CC = mpicc + CXX = mpicxx + FORTRAN = mpif90 +endif + +# ---------------------------------------------------------------- +# Normally you would not need to change any line beyond this point +# ---------------------------------------------------------------- + +# Here I am finding the compiler group +ifeq ($(seq),1) + F_COMP = $(FORTRAN) +else + F_COMP = $(firstword $(shell $(FORTRAN) -show)) +endif +ifeq ($(F_COMP),gfortran) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),gcc) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),g77) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),f95) + COMP_GRP = gnu +endif +ifeq ($(F_COMP),ifort) + COMP_GRP = intel +endif +ifeq ($(F_COMP),pgf90) + COMP_GRP = pgi +endif +ifeq ($(F_COMP),pgf77) + COMP_GRP = pgi +endif + +ifeq ($(OS),Darwin) + COMP_GRP = gnu +endif + +# If profiling is requested +prof = 0 +ifeq ($(prof),1) + FFLAGS += -pg -g +endif + +# If debuging is requested +ifeq ($(debug),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -g -Wall -Wconversion -Wline-truncation -pedantic -fimplicit-none -fbacktrace -fbounds-check -p -fcheck=all #-ffpe-trap=invalid,zero,overflow,underflow + CXXFLAGS += -g -Wall -pedantic -fbounds-check + CFLAGS += -g -Wall -pedantic -fbounds-check + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -g -traceback -check all -check bounds -check uninit -ftrapuv -debug all -fpe0 + CFLAGS += -g -traceback -check-uninit -fpe0 + CXXFLAGS += -g -traceback -check-uninit -fpe0 + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += # You need to add debuging flag for pgi compiler here + endif + OBJ_DIR = dbg +endif + +# If openMP parallelization is requested +mp = 0 +ifeq ($(mp),1) + ifeq ($(COMP_GRP),gnu) + FFLAGS += -fopenmp + endif + ifeq ($(COMP_GRP),intel) + FFLAGS += -openmp + endif + ifeq ($(COMP_GRP),pgi) + FFLAGS += -mp + endif +endif + +# To make directories a bit cleaner with intel compiler +ifeq ($(COMP_GRP),intel) + FFLAGS += -module $(OBJ_DIR) +endif +ifeq ($(COMP_GRP),gnu) + FFLAGS += -J $(OBJ_DIR) + CFLAGS += -J $(OBJ_DIR) + CXXFLAGS += -J $(OBJ_DIR) +endif + +#LAPACK Library +ifeq ($(COMP_GRP),gnu) + LAPACK_INC = + LAPACK_LIB = + LAPACK_LFLAGS = -llapack +endif + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f new file mode 100644 index 00000000..991bbb50 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/initial_values_final.f @@ -0,0 +1,4 @@ +c Value of your unknown at time equal to zero (This is going to be used +c ONLY when you are initiating simulation) + + tZeroX = (/ 0.00D+00 /) diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f new file mode 100644 index 00000000..af40e47d --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/include/parameters_final.f @@ -0,0 +1,15 @@ +! Specifing the constants + REAL(KIND=8), PARAMETER :: ! Order of surfaces: 3~8 + & Rmyo = 4.51000000D-02, + & Lv = 2.74000000D-05, + & Kv = 8.16000000D-04, + & Kao = 8.59000000D-05, + & pi = 3.14159265D+00, + & Afill = 6D+01, + & tfill = 5D+00 +! The parameters which might change through the code + REAL(KIND=8) Clung, Rlung + Clung= 2.56000000D-02 + Rlung= 7.29000000D-01 + + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f new file mode 100644 index 00000000..1d3d8908 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/GenBC.f @@ -0,0 +1,232 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c Here data are received from Phasta and it setup the data required for +c integration of ODE's inside FINDX + PROGRAM GenBC + USE COM + USE, INTRINSIC :: IEEE_ARITHMETIC + IMPLICIT NONE + + INTEGER i, n, nTimeStep + REAL(KIND=8) t, dt, tFinal, pMin, temp + CHARACTER(LEN=2048) string + CHARACTER(LEN=32) sTmp + CHARACTER flag + + REAL(KIND=8), ALLOCATABLE :: Qi(:), Qf(:), Pi(:), Pf(:), X(:), + 2 Xo(:), f(:,:) + +c flag = I : Initializing +c flag = T : Iteration loop +c flag = L : Last iteration +c flag = D : Derivative +c +c Here is the period of time that you should integrate for +c each of these flags: +c +c Flags I D&T&L +c ^ ^ +c 3D code time step: N.............................................N+1 +c 0D code time step: 1......................................nTimeStep+1 +c Flowrates: Qi............................ ................Qf +c Time, t: tInitial..............................tInitial+tFinal + + + CALL INITIALIZE(nTimeStep) + +c******************************************************************** +c Block for reading the data from phasta + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED'); + READ (1) flag + READ (1) tFinal + READ (1) i + IF (i .NE. nDirichletSrfs) THEN + PRINT *, 'Error: Number of Dirichlet Surfaces from Phasta is:', + 2 i + PRINT *, 'While nDirichletSrfs is equal to:', nDirichletSrfs + PRINT * + PRINT *, 'Number of Dirichlet surfaces defined in solver.inp', + 2 ' should match with nDirichletSrfs defined in USER.f' + STOP + END IF + READ (1) i + IF (i .NE. nNeumannSrfs) THEN + PRINT *, 'Error: Number of Neumann Surfaces from Phasta is:', + 2 i + PRINT *, 'While nNeumannSrfs is equal to:', nNeumannSrfs + PRINT * + PRINT *, 'Number of Neumann surfaces defined in solver.inp', + 2 ' should match with nNeumannSrfs defined in USER.f' + STOP + END IF + IF (nDirichletSrfs .GT. 0) THEN + IF (qCorr .OR. pCorr) THEN + PRINT *, 'You should only use P/Q correction when all', + 2 ' the surfaces are Neumann surfaces' + STOP + END IF + END IF + + ALLOCATE (Pi(nDirichletSrfs), Pf(nDirichletSrfs), + 2 PDirichlet(nDirichletSrfs,4)) + ALLOCATE (Qi(nNeumannSrfs), Qf(nNeumannSrfs), + 2 QNeumann(nNeumannSrfs,4)) + + DO i=1, nDirichletSrfs + READ(1) Pi(i) + READ(1) Pf(i) + END DO + Pi = Pi/pConv + Pf = Pf/pConv + + DO i=1,nNeumannSrfs + READ (1) Qi(i) + READ (1) Qf(i) +c Print to console and file pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Flowrates to genBC: ', Qi(i), Qf(i) + END DO + CLOSE(1) + Qi = Qi/qConv + Qf = Qf/qConv + +c******************************************************************** +c Block for initializing the unknowns + + IF (flag .EQ. 'I') THEN + nTimeStep = 0 + ELSE + dt = tFinal/REAL(nTimeStep,8) + END IF + ALLOCATE (X(nUnknowns), Xo(nUnknowns), f(nUnknowns,4), + 2 offset(nUnknowns), Xprint(nXprint)) + Xprint = 0D0 + + OPEN (1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED') + READ (1) t + DO i=1,nUnknowns + READ (1) Xo(i) + END DO + CLOSE (1) + +c******************************************************************** +c Setting up the system of equations + offset = 0D0 + DO n=1, nTimeStep + DO i=1, 4 + temp = (REAL(n-1,8) + REAL(i-1,8)/3D0)/REAL(nTimeStep,8) + + QNeumann(:,i) = Qi + (Qf-Qi)*temp + PDirichlet(:,i) = Pi + (Pf-Pi)*temp + + IF (qCorr) THEN + temp = SUM(QNeumann(:,i))/REAL(nNeumannSrfs,8) + QNeumann(:,i) = QNeumann(:,i) - temp + END IF + END DO + + X = Xo + CALL FINDF (t, X, f(:,1), QNeumann(:,1), + 2 PDirichlet(:,1)) + X = Xo + dt*f(:,1)/3D0 + + CALL FINDF (t+dt/3D0, X, f(:,2), QNeumann(:,2), + 2 PDirichlet(:,2)) + X = Xo - dt*f(:,1)/3D0 + dt*f(:,2) + + CALL FINDF (t+dt*2D0/3D0, X, f(:,3), QNeumann(:,3), + 2 PDirichlet(:,3)) + X = Xo + dt*f(:,1) - dt*f(:,2) + dt*f(:,3) + + CALL FINDF (t+dt, X, f(:,4), QNeumann(:,4), + 2 PDirichlet(:,4)) + + f(:,1) = (f(:,1) + 3D0*f(:,2) + 3D0*f(:,3) + f(:,4))/8D0 + Xo = Xo + dt*f(:,1) + t = t + dt + END DO + +c******************************************************************** +c Time to write the results + X = Xo + IF (pCorr .AND. flag.NE.'D') THEN + pMin = X(srfToXPtr(1)) + DO i=2, nNeumannSrfs + IF (X(srfToXPtr(i)) .LT. pMin) THEN + pMin = X(srfToXPtr(i)) + END IF + END DO + ELSE + pMin = 0D0 + END IF + +c Writing nDirichlet flowrates here + OPEN (1, FILE='GenBC.int', STATUS='OLD', FORM='UNFORMATTED') + DO i=1, nDirichletSrfs + IF(IEEE_IS_NAN(X(srfToXdPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF + WRITE (1) X(srfToXdPtr(i))*qConv + END DO + +c Writing nNeumannSrfs pressures here + DO i=1, nNeumannSrfs + IF(IEEE_IS_NAN(X(srfToXPtr(i)))) THEN + PRINT*, 'Error! NAN encountered..' + STOP + END IF +c Print out pressure from genBC to svFSI for debugging + PRINT*, 'BC i = ', i, + 2 'Pressure from genBC: ', offset(srfToXPtr(i)) + +c Here, saving offset value to X itself. Necessary for keeping track of +c pressure from one time step to the next + X(srfToXPtr(i)) = offset(srfToXPtr(i)) + +c WRITE (1) (X(srfToXPtr(i)) - pMin +c 2 + offset(srfToXPtr(i)))*pConv + WRITE (1) (X(srfToXPtr(i)))*pConv + END DO + CLOSE(1) + + IF (flag .EQ. 'L') THEN + OPEN(1, FILE='InitialData', STATUS='OLD', FORM='UNFORMATTED'); + WRITE (1) t + DO i=1,nUnknowns + WRITE (1) X(i) + END DO + CLOSE(1) + +c PRINT *,'Before AllData' + + OPEN(1, FILE='AllData', STATUS='UNKNOWN', ACCESS='APPEND'); + string = '' + DO i=1, nUnknowns + WRITE (sTmp,"(ES14.6E2)") X(i) + string = TRIM(string)//sTmp + END DO + DO i=1, nXprint + WRITE (sTmp,"(ES14.6E2)") Xprint(i) + string = TRIM(string)//sTmp + END DO + WRITE (1,"(A)") TRIM(string) + CLOSE(1) + END IF + + DEALLOCATE (Pi) + DEALLOCATE (Pf) + DEALLOCATE (PDirichlet) + DEALLOCATE (Qi) + DEALLOCATE (Qf) + DEALLOCATE (QNeumann) + DEALLOCATE (X) + DEALLOCATE (Xo) + DEALLOCATE (f) + DEALLOCATE (srfToXdPtr) + DEALLOCATE (srfToXPtr) + DEALLOCATE (offset) + DEALLOCATE (Xprint) + + END PROGRAM GenBC diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f new file mode 100644 index 00000000..6e472f3b --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/Modules.f @@ -0,0 +1,17 @@ +c Created by Mahdi Esmaily Moghadam 05-25-2011 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + + MODULE COM + + LOGICAL pCorr, qCorr + + INTEGER nUnknowns, nDirichletSrfs, nNeumannSrfs, nXprint + + REAL(KIND=8) pConv, qConv + + INTEGER, ALLOCATABLE :: srfToXPtr(:), srfToXdPtr(:) + + REAL(KIND=8), ALLOCATABLE :: QNeumann(:,:), PDirichlet(:,:), + 2 offset(:), Xprint(:) + + END MODULE COM diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f new file mode 100644 index 00000000..e48f7c19 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/genBC_svFSIplus/src/USER.f @@ -0,0 +1,143 @@ +c MUSC 2 Immediate postop +c Ethan Kung keo@ucsd.edu + +c Created by Mahdi Esmaily Moghadam 12-01-2010 +c Please report any problem to mesmaily@ucsd.edu, memt63@yahoo.com + +c This subroutine initializes the parameters, you need to read the +c comments and specify them carefuly +c-------------------------------------------------------------------- +c This is an example for RCR boundary condition with parameters +c Rd, R, and C which are distal, and proximal resistance and capacitor. + + SUBROUTINE INITIALIZE(nTimeStep) + USE COM + IMPLICIT NONE + INTENT(OUT) nTimeStep + + LOGICAl ierr + INTEGER i, nTimeStep + REAL(KIND=8), ALLOCATABLE :: tZeroX(:) +c +c******************************************************************** +c For instance if pressure in 3D solver is in cgs and here mmHg +c pConv=1334=1.334D3, also if flowrate in 3D solver is mm^3/s and +c here is mL/s qConv=1000=1D3. In the case both solver are using same +c unites you can set these two conversion coefficients to 1D0 + pConv = 1D0 + qConv = 1D0 + +c Only when all the surfaces of you model are coupled with NeumannSrfs +c you may set this to .TRUE. + pCorr = .FALSE. + qCorr = .FALSE. + +c******************************************************************** +c Block for your inputs + +c These two value should match to that of defined in solver.inp + nDirichletSrfs = 0 + nNeumannSrfs = 1 +c Number of unknowns that you need inside your lumped parameter network + nUnknowns = 2 +c Number of time step between N and N+alpha, for higher values this code +c would be more accurate and more costly + nTimeStep = 10 + +c Number of parameters to be printed in AllData file (the first +c nUnknowns columns are the "X" values) + nXprint = 2 + +c-------------------------------------------------------------------- +c You don't need to change this part ! + ALLOCATE (tZeroX(nUnknowns), srfToXdPtr(nDirichletSrfs)) ! + ALLOCATE (srfToXPtr(nNeumannSrfs)) ! + tZeroX = 0D0 +c-------------------------------------------------------------------- + + INCLUDE "initial_values_final.f" + +c-------------------------------------------------------------------- +c You don't need to change this part ! + INQUIRE (FILE='InitialData', EXIST=ierr) ! + IF (.NOT.ierr) THEN ! +c PRINT *, 'Initializing unknowns in LPM' ! + OPEN (1, FILE='InitialData',STATUS='NEW',FORM='UNFORMATTED')! + WRITE (1) 0D0 ! + DO i=1, nUnknowns ! + WRITE (1) tZeroX(i) ! + END DO ! + CLOSE(1) ! + END IF ! +c-------------------------------------------------------------------- + +c Surface to X pointer: this defines which Unknown corresponds to which +c suface in "List of Neumann Surfaces" inside solver.inp +c For example, if you have "List of Neumann Surfaces= 2 8 4 ...." +c and you set "srfToXPtr = (/5,3,9,.../)" +C this means X(5) corresponds to surface 2, X(3) <-> surface 8, +c and X(9) <-> surface 4 +c Also, Q(1) corresponds to the first element in srfToXPtr, in this +c example, Q(1)<=>X(5), Q(2)<=>X(3) + srfToXPtr = (/1/) +c srfToXdPtr = (/1/) + + END SUBROUTINE INITIALIZE + +c#################################################################### +c Here you should find the f_i=dx_i/dt, based on the following parameters: +c current x_i: x(i) +c Current time: t +c Flowrates from 3D code: Q(i) +c Pressure from Dirichlet faces: P(i) + + SUBROUTINE FINDF(t, x, f, Q, P) + USE COM + IMPLICIT NONE + INTENT(IN) t, Q + INTENT(OUT) f + + REAL(KIND=8) t, x(nUnknowns), f(nUnknowns), Q(nNeumannSrfs), + 2 P(nDirichletSrfs) + +! These are the dumy variables + REAL(KIND=8) Pivc,Ppwc,Pcon,Qout, Phigh, Pmax, tramp, R + + INCLUDE "parameters_final.f" + +!******************************************************************** +! This is the only case that x is directly manipulated, ugly, but the +! only way!!! + +! Implement a resistance filling network so that the ventricle pressure +! Pv = Phigh - Q*R +! This fills the sphere at an approximately constant rate + +! x(1): Pressure applied to endo surface. Actually, we will assign the pressure +! by setting the offset instead +! x(2): Accumulated volume during inflation. Not necessary, but for insight + +! Define a pressure ramp for Phigh (1e7 dyne/cm^2 in 0.11 second, then constant) + Pmax = 1.0D+7 + tramp = 0.1 + IF (t .LT. tramp) THEN + Phigh = Pmax * t / tramp + ELSE + Phigh = Pmax + ENDIF + R = 1.0D+5 + +! The main body of equations + f(1) = 0 ! The endo surface pressure (we will set it with offset) + f(2) = -Q(1) ! Integrate the flow rate to compute accumulated volume + +! Define offset pressure. + offset(1) = Phigh + Q(1) * R + +c Assign the additional parameters to be printed + Xprint(1) = offset(1) ! This is the actual value of pressure for this sim + Xprint(2) = Q(1) ! The flowrate/dVdt + RETURN + END SUBROUTINE FINDF + + diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..37ab77b3 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:516ca5f67f51c432410c65a1f98d73baca653b0378664607a20eba3fd5e94db4 +size 19832 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..cf145a1d --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9b7c72acb6c365fae7c33091f8d932c95863d95ffee7fed934a2f832aeba976d +size 27192 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 00000000..e284a4c6 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4028808552b03a2854ce33bf5016243c6b994e62a61928625691b72938f1480d +size 10038 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 00000000..66272c04 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f125c7001fc4a5ad90313fcff8ef24749ccd054a6b045521dfc3d06169160a09 +size 11982 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp new file mode 100644 index 00000000..2d90c222 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2050d191642b719f651de6587571a7021eff67801eddd97c0fec397f3e9b62cd +size 3940 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/post_processing_functions.py b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/post_processing_functions.py new file mode 100644 index 00000000..503dbeb5 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/post_processing_functions.py @@ -0,0 +1,700 @@ +# Before running on Sherlock, need to have installed: +# +# pyvista: pip3 install pyvista +# ffmpy: pip3 install ffmpy +# +# Also, need to load the following required modules on Sherlock: +# +# ml system mesa ffmpeg +import glob +import os +import pyvista as pv # for VTK mesh manipulations +import numpy as np +#import ffmpy # to convert .avi movie to .mp4 movie + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start time, end time, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSI simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_time, end_time, step): A tuple of 3 integers, giving the first + time step of results to process, the last time step, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_time = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start time:', start_time) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_time = int("".join([i for i in end_file_name if i.isdigit()])) + print('End time:', end_time) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_time + print('Step:', step) + + return (start_time, end_time, step) + +def calc_volume_struct(start_time, end_time, step, results_folder, reference_surface): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI struct simulation, in which a model of the myocardium is inflated. + + Calculate the volume in the following steps + 1) Sample the result.vtu file onto the reference surface + 2) Warp the samples surface by the Displacement + 3) Flat fill any holes in the warped surface + 4) Calculate the volume of the warped and filled surface + + The units of volume are whatever units used in .vtu files, cubed. For example, + if units of length in the .vtu files are microns, then the volume calculated + here is cubic microns. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + reference_surface: The absolute file path of the .vtp file containing + the undeformed surface corresponding to the deformed surface of which + we want to compute the volume. + + output_folder: The absolute file path of the folder to which we output + intermediary results of the volume calculation process. + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_struct_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating volumes ##') + + # Load reference surface onto which we sample + ref_surface = pv.read(f"{reference_surface}") + + # Compute initial lumen volume + t = [] + vol = [] + k = 0 + # The initial volume is obtain by computing the lumen volume of the reference + # configuration, which is simply volume of the filled in reference surface. + ref_lumen = ref_surface.fill_holes(100) # 100 is the largest size of hole to fill + + # -------------------------------------------------------------------- + + # Recompute normals, incase the normals of the cap are opposite + ref_lumen.compute_normals(inplace=True) + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + ref_lumen.save(f'{output_folder}/resampled_warped_and_filled_{k:03d}.vtp') + + # Add time and volume to arrays + t.append(k) + vol.append(ref_lumen.volume) + + print(f"Iteration: {k}, Volume: {ref_lumen.volume}") + + + # Loop through results files at each time > 0 + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}_cpp.vtu")) + + # Sample result onto ref_surface + resampled = ref_surface.sample(result) + + # Warp resampled surface by displacement (needed for current configuration + # normals, as well volume calculation) + warped = resampled.warp_by_vector('Displacement') + + ## ------ Compute volume by warping resampled, then filling the hole ----- ## + + #warped.save(f'{output_folder}/resampled_warped_{k:03d}.vtp') + + # Fill the holes + # EITHER + # ----------------------------------------------------------------------- + # Fix the mesh (i.e. fill the holes where the inlet and outlet caps are) + #meshfix = mf.MeshFix(warped) + #meshfix.repair() + + # Convert back to pyvista polydata mesh + #lumen = meshfix.mesh + + ## PyMeshFix fills the hole, but leads to normals oriented inconsistently ## + ## The PyVista function .fill_holes() produces a bad cap mesh, but with ## + ## consistent normal vectors. The bad mesh is not important if we're only ## + ## computing volume ## + + # ---------------------------- OR ------------------------------------ + lumen = warped.fill_holes(100) # 100 is the largest size of hole to fill + + # -------------------------------------------------------------------- + + # Recompute normals, incase the normals of the cap are opposite + lumen.compute_normals(inplace=True) + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + lumen.save(f'{output_folder}/resampled_warped_and_filled_{k:03d}.vtp') + + # Add time and volume to arrays + t.append(k) + vol.append(lumen.volume) + + print(f"Iteration: {k}, Volume: {lumen.volume}") + + return (t, vol) +def calc_dVdt_struct(start_time, end_time, step, results_folder, reference_surface): + ''' + Computes the rate of change of volume of a closed (or partially closed) + surface. Intended to be used to calculate the rate of change of ventricular + volume from a struct simulation of an inflating ventricle. + + Compute the rate of change of volume by computing the velocity flux over the + reference surface, which should be the endo surface of the ventricle. + + dV/dt = int_{Gamma_t} (u . n) dA + + where Gamma_t is the reference surface in the current configuration, u is + the velocity on the reference surface, and n is the surface normal vector on + the reference surface + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + reference_surface: The absolute file path of the .vtp file containing + the undeformed surface corresponding to the deformed surface of which + we want to compute the volume. + + output_folder: The absolute file path of the folder to which we output + intermediary results of the volume calculation process. + + Returns: (t, dVdt), a tuple of lists of length number of time steps. t + contains the time step, and Q contains the rate of change of volume at that + time step. + ''' + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_struct_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + # Initialize volume array and time step lists + dVdt = [0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}_cpp.vtu")) + + # Load reference surface onto which we sample + ref_surface = pv.read(f"{reference_surface}") + + # Sample result onto ref_surface + ref= ref_surface.sample(result) + + # Warp ref surface by displacement (needed for current configuration + # normals, as well volume calculation) + warped = ref.warp_by_vector('Displacement') + + # Compute average velocity for each element + warped = warped.point_data_to_cell_data() + ref = ref.point_data_to_cell_data() + + # Recompute normals. Need this after point_data_to_cell_data(), because + # the average of point normals is not the same as a cell normal + warped.compute_normals(inplace=True) + ref.compute_normals(inplace=True) + + # Compute element areas + warped = warped.compute_cell_sizes() + ref = ref.compute_cell_sizes() + + # Save warped and filled lumen (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + warped.save(f'{output_folder}/warped_{k:03d}.vtp') + ref.save(f'{output_folder}/ref_{k:03d}.vtp') + + # Compute dVdt from warped by computing velocity flux. Also, compute + # velocity flux using quantities from reference surface to compare. + # + # !! Looks like the velocity flux output in B_ST_Velocity_flux.txt is + # computed using reference surface normals and reference surface areas !! + dVdt_t = 0 + + cell_normals = warped.cell_data['Normals'] + cell_vels = warped.cell_data['Velocity'] + cell_areas = warped.cell_data['Area'] + + cell_normals_ref = ref.cell_data['Normals'] + cell_vels_ref = ref.cell_data['Velocity'] + cell_areas_ref = ref.cell_data['Area'] + + u_dot_n = np.sum(cell_vels * cell_normals, axis = 1) + dVdt_t = np.sum(u_dot_n * cell_areas) + + u_dot_n_ref = np.sum(cell_vels_ref * cell_normals_ref, axis = 1) + dVdt_t_ref = np.sum(u_dot_n_ref * cell_areas_ref) + + # Add time and volume to arrays + t.append(k) + dVdt.append(dVdt_t) + + print(f"Iteration: {k}, dVdt: {dVdt_t}, dVdt_ref: {dVdt_t_ref}") + + return (t, dVdt) + +def calc_pressure_struct(input_file, pressure_dat_file, t): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI struct simulation, in which a model of the myocardium is inflated. + + Calculates pressure on the endo surface by reading the input file and + pressure load file. The units of pressure are whatever units are used in the + pressure load file (usually Pa) + + Args: + input_file: The svFSI input file of the simulation whose results we are + now processing. + + pressure_dat_file: The file containing the pressure load information for + the simulation whose results we are now processing. + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + pressure: Pressure evaluated at list of time steps in t + """ + print('\n## Calculating pressure ##') + + # Compute pressure at each time step by reading input file and pressure file + with open(input_file) as f: + lines = f.readlines() + for line in lines: + if 'Number of time steps' in line: + num_time_steps = int(line.split(':')[-1]) + + with open(pressure_dat_file) as f: + f.readline() + line = f.readline() + vals = line.split() + t0 = float(vals[0]) + p0 = float(vals[1]) + line = f.readline() + vals = line.split() + t1 = float(vals[0]) + p1 = float(vals[1]) + + pressure = np.linspace(p0, p1, num_time_steps+1) + + # Sample pressure onto t + pressure = pressure[t] + + return pressure + +def calc_pressure_struct_genBC(genBC_out_file, col, t): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI struct simulation with genBC, in which a model of the myocardium is + inflated. + + Calculates pressure on the endo surface by reading the genBC output file, + usually called AllData. The units of pressure are whatever units are used in + the genBC input file. + + Args: + genBC_out_file: The file containing the pressure load information for + the simulation whose results we are now processing. Usually AllData + + col: The column in AllData that contains the pressure at each time step + col = 1 refers to the first column + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + pressure: Pressure evaluated at list of time steps in t + """ + print('\n## Calculating pressure ##') + + # Read pressure from AllData + pressure = [0] + with open(genBC_out_file) as f: + lines = f.readlines() + for line in lines: + vals = line.split() + pressure.append(float(vals[col-1])) + + pressure = np.array(pressure) + + # Sample pressure onto t + pressure = pressure[t] + + return pressure +def calc_volume_struct_genBC(genBC_out_file, col, t): + """ + Calculate the accumulated lumen volume at each time step from the results of + an svFSI struct simulation with genBC, in which a model of the myocardium is + inflated. + + Calculates accumulated volumeby reading the genBC output file, + usually called AllData. In genBC, you must define an additional unknown + x(i) such that f(i) = -Q(1). This will compute x(i) as the integral of the + flow rate, i.e. the accumulated volume + + Args: + genBC_out_file: The file containing the pressure load information for + the simulation whose results we are now processing. Usually AllData + + col: The column in AllData that contains the pressure at each time step + col = 1 refers to the first column + + t: List of time steps (iterations) at which we want to output the + pressure. + Returns: + volume: Volume evaluated at list of time steps in t + """ + print('\n## Calculating volume from AllData ##') + + # Read pressure from AllData + volume = [0] + with open(genBC_out_file) as f: + lines = f.readlines() + for line in lines: + vals = line.split() + volume.append(float(vals[col-1])) + + volume = np.array(volume) + + # Sample pressure onto t + volume = volume[t] + + return volume + +def calc_volume_FSI(volume_file): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the volume by reading the output file V_FS_Volume_average.txt. This + is produced if the following is included in svFSI.inp + + Output: Volume_integral { + Volume: t + } + + This gives the volume of the fluid domain at each time step. + + The units of volume are whatever units used in svFSI.inp, cubed. For example, + if units of length in svFSI.inp are meters, then the volume calculated + here is cubic meters. Note that the mesh scale factor affects these units. If + the .vtu file has units of microns, but we apply a mesh scale factor of 1e-6 + to convert to meters in svFSI.inp, then the volume output here will be in meters^3 + + Args: + volume_file: Path to file V_FS_Volume_average.txt from simulation you + want to process + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + print('\n## Calculating volume ##') + + # Initialize time and volume arrays + t = [] + vol = [] + + # Read volume data from V_FS_Volume_average.txt + with open(volume_file) as f: + header = f.readline() # Read header + initial_volumes = f.readline() # Read domain initial volumes + f.readline() # Skip blank line + + # Read volume data + # (For some reason, the first volume is the initial volume, even though + # this corresponds to after the first time step) + lines = f.readlines() + k = 0 + for line in lines: + + # Add time iteration to t array + t.append(k) + k += 1 + + # Add fluid domain (Domain 0) volume to volume array + V_t = float(line.split()[0]) + vol.append(V_t) + + return (t, vol) + + +def calc_volume_FSI_2(start_time, end_time, step, results_folder): + """ + Calculate the ventricular lumen volume at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the volume in the following steps + 1) Threshold to extract fluid domain + 2) Warp the fluid domain by the Displacement + 3) Calculate the volume of the warped fluid domain + + The units of volume are whatever units used in .vtu files, cubed. For example, + if units of length in the .vtu files are microns, then the volume calculated + here is cubic microns. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + Returns: (t, vol), a tuple of lists of length number of time steps. t + contains the time step, and vol contains the volume at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_volume_FSI_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating volumes ##') + + # Initialize volume array and time step lists + vol = [0.0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Threshold to fluid domain + fluid = result.threshold(value = (1,1), scalars = 'Domain_ID') + + # Warp fluid domain by displacement (yields current fluid domain) + warped = fluid.warp_by_vector('Displacement') + + # Save warped fluid domain (to check geometry and normals) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + #warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtp') + warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtu') + + + + # Add time and volume to arrays + t.append(k) + vol.append(warped.volume) + + print(f"Iteration: {k}, Volume: {warped.volume}") + + return (t, vol) + + +def calc_pressure_FSI(pressure_file): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the pressure by reading the output file V_FS_Pressure_average.txt. This + is produced if the following is included in svFSI.inp + + Output: Volume_integral { + Pressure: t + } + This gives the volume averaged fluid pressure at each time step. + + + The units of pressure are deduced by whatever units system you are using in + the input file. For SI (m, kg, s), the units of pressure are Pa. + + Args: + pressure_file: Path to file V_FS_Pressure_average.txt from simulation you + want to process + + Returns: (t, pressure), a tuple of lists of length number of time steps. t + contains the time step, and pressure contains the pressure at that time step. + """ + print('\n## Calculating pressure ##') + + # Initialize time and pressure arrays + t = [0] + pressure = [0] + + # Read pressure data from V_FS_Pressure_average.txt + with open(pressure_file) as f: + header = f.readline() # Read header + initial_volumes = f.readline() # Read domain initial volumes + f.readline() # Skip blank line + + # Read pressure data + lines = f.readlines() + for line in lines: + + # Add next time iteration to t array + t.append(t[-1] + 1) + + # Add fluid domain (Domain 0) pressure to pressure array + p_t = float(line.split()[0]) + pressure.append(p_t) + + return (t, pressure) + +def calc_pressure_FSI_2(start_time, end_time, step, results_folder): + """ + Calculate the ventricular lumen pressure at each time step from the results of + an svFSI FSI simulation, in which a model of the myocardium is inflated. + + Calculate the pressure in the following steps + 1) Threshold to extract fluid domain + 2) Warp the fluid domain by the Displacement + 3) Compute the volume of each cell + 4) Compute the pressure in each cell + 5) Compute the volume-averaged pressure + + The units of pressure are deduced by whatever units system you are using in + the input file. For SI (m, kg, s), the units of pressure are Pa. + + Args: + start_time: The first svFSI result file to process + + end_time: The last svFSI result file to process + + step: The step in svFSI result files to process + + results_folder: The absolute file path of the svFSI results folder + (usually something/something/16-procs/) + + Returns: (t, pressure), a tuple of lists of length number of time steps. t + contains the time step, and pressure contains the pressure at that time step. + """ + + # Create folder to contain intermediary meshes (mostly for checking for errors) + output_folder = results_folder + '/../' + 'calc_pressure_FSI_output' + # checking if the directory exists + if not os.path.exists(output_folder): + # if the directory is not present then create it. + os.makedirs(output_folder) + + print('\n## Calculating pressures ##') + + # Initialize volume array and time step lists + pressure = [0.0] + t = [0] + + # Loop through results files at each time + for k in range(start_time, end_time+1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Threshold to fluid domain + fluid = result.threshold(value = (1,1), scalars = 'Domain_ID') + + # Warp fluid domain by displacement (yields current fluid domain) + warped = fluid.warp_by_vector('Displacement') + + # Compute volumes and areas + sized = warped.compute_cell_sizes() + + + # Compute pressure in each cell from pressure at each point + sized = sized.point_data_to_cell_data() + + + + # Save processed fluid domain (to check geometry, normals, and values) + # (Hopefully the normals on the filled cap will be consistent with the normals + # on the rest of the surface, but you should check to make sure.) + #warped.save(f'{output_folder}/fluid_warped_{k:03d}.vtp') + sized.save(f'{output_folder}/fluid_processed_{k:03d}.vtu') + + # Grab volumes for all elements in mesh + elem_volumes = sized.cell_data['Volume'] + # Grab pressure for all elements in mesh + elem_pressures = sized.cell_data['Pressure'] + + # Compute volume averaged pressure by looping over fluid domain elements + p_avg = 0 + for i in range(len(elem_volumes)): + p_avg += elem_pressures[i] * elem_volumes[i] + p_avg = p_avg / sum(elem_volumes) + + + # Add time and pressure to arrays + t.append(k) + pressure.append(p_avg) + + print(f"Iteration: {k}, Pressure: {p_avg}") + + return (t, pressure) + +#def avi2mp4(input_movie_file, output_movie_file): +# ff = ffmpy.FFmpeg( +# inputs={input_movie_file: None}, +# outputs={output_movie_file: '-pix_fmt yuv420p -y'} +# # flag -pix_fmt yuv420p needed so that movie is compatible with QuickTime +# # flag -y needed to overwrite mp4 file if it already exists +# ) +# ff.run() \ No newline at end of file diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/process_results.py b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/process_results.py new file mode 100644 index 00000000..ad94716d --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/process_results.py @@ -0,0 +1,154 @@ +# Python script to calculate lumen volume and fluid pressure on endo surface from +# results of svFSI struct +# +# Requires two meshes +# - svFSI result.vtu (with Displacement array) +# - reference surface file (polygonal mesh), representing the undeformed surface +# corresponding to the deformed surface of which we want to compute the volume +# +# We calculate the volume in the following steps +# 1) Sample the result.vtu file onto the reference surface +# 2) Warp the samples surface by the Displacement +# 3) Flat fill any holes in the warped surface +# 4) Calculate the volume of the warped and filled surface +# +# We calculate the fluid pressure on the endo surface by reading the input file +# and pressure load file. + +import os # for checking and creating directories, and loading modules +import matplotlib.pyplot as plt +import numpy as np +import sys + +# Import custom functions useful for post-processing +from post_processing_functions import * + +## -------- PARAMETERS TO CHANGE -------------------- ## + +# Simulation folder file path +sim_folder = os.path.dirname(os.path.realpath(__file__)) + +# Optionally read sim_folder as command line argument +if len(sys.argv) == 2: + print(sys.argv) + sim_folder = os.path.abspath(sys.argv[1]) + print(sim_folder) + + +# Undeformed surface, corresponding to deformed surface on result.vtu of which we want to compute the volume +reference_surface = os.path.join(sim_folder, 'mesh/mesh-surfaces/endo.vtp') + +# Undeformed vtu file (the original volume mesh) +reference_volume = os.path.join(sim_folder, 'mesh/mesh-complete.mesh.vtu') + +# svFSI results folder, containing results.vtu +#results_folder = os.path.join(sim_folder, 'results_svfsi/') +results_folder = os.path.join(sim_folder, '4-procs/') + +# File containing genBC output +#alldata_file = 'AllData_svfsi_vvedula22' +alldata_file = 'AllData' + +# Input file path +input_file = os.path.join(sim_folder, 'svFSI.inp') + +# Pressure file (unused) +pressure_dat_file = os.path.join(sim_folder, 'pressure.dat') + + +## -------------------- END PARAMETERS TO CHANGE ------------------------ ## + + +# Automatically determine the start time, end time, and step size based on all +# results file in results_folder +print(results_folder) +(start_time, end_time, step) = get_start_end_step(results_folder) +# Option to manually set start time, end time, and time step of results files to process +#start_time = 5 +#end_time = 85 +#step = 5 + +# Compute lumen volume from simulation results +(t_3D, vol_3D) = calc_volume_struct(start_time, end_time, step, results_folder, reference_surface) +vol_3D_cm3 = np.array(vol_3D) * (100)**3 # cm^3 + +# Compute lumen dVdt from simulation results +(t_dVdt_3D, dVdt_3D) = calc_dVdt_struct(start_time, end_time, step, results_folder, reference_surface) +# Convert volume to cm^3/s +dVdt_3D = np.array(dVdt_3D) # cm/s * m^2 +dVdt_3D_cm3 = dVdt_3D * (100)**2 # cm^3/s + +# Compute lumen volume from AllData file +vol_0D_cm3 = calc_volume_struct_genBC(os.path.join(sim_folder, alldata_file), 2, t_3D) +# Add on initial volume +vol_0D_cm3 += vol_3D_cm3[0] # cm^3 + +# Compute flow rate = dVdt from AllData file +dVdt_0D_cm3 = calc_volume_struct_genBC(os.path.join(sim_folder, alldata_file), 4, t_dVdt_3D) + + +# Compute lumen pressure at iterations in t +pressure = calc_pressure_struct_genBC(os.path.join(sim_folder, "AllData"), 1, t_3D) # dynes/cm^2 + +# Convert pressure from dynes/cm^2 to mmHg +#pressure_mmHg = pressure * 0.000750062 + + + +# Combine time step, pressure, and volume into one array +#PV = np.column_stack((t_3D, pressure, vol_3D_cm3)) + +#print('\n## Outputing pressure-volume data and plot ##') + +# Write pressure and volume to file +#np.savetxt(results_folder + '/../' + "pv.txt", PV, header = 'Timestep Pressure[mmHg] Volume[m^3]') + + +# Plot pressure vs. volume +fig, ax = plt.subplots() +ax.plot(vol_3D_cm3, pressure, linewidth=2.0, marker = 'o') +ax.set_xlabel('Volume [cm^3]') +ax.set_ylabel('Pressure [dyne/cm^2]') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'pv_plot')) + +# Plot dVdt vs. time +#fig, ax = plt.subplots() +#ax.plot(t_dVdt, dVdt_m3, linewidth=2.0, marker = 'o') +#ax.set_xlabel('Time') +#ax.set_ylabel('dVdt') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +#plt.savefig(os.path.join(sim_folder, 'dVdt_plot')) + +# Plot dVdt vs. Pressure +#fig, ax = plt.subplots() +#ax.plot(dVdt_m3, pressure_mmHg, linewidth=2.0, marker = 'o') +#ax.set_xlabel('dVdt') +#ax.set_ylabel('Pressure') +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +#plt.savefig(os.path.join(sim_folder, 'dVdtvsP_plot')) + +# Plot 3D and 0D dVdt +fig, ax = plt.subplots() +ax.plot(dVdt_3D_cm3, label = 'dVdt_3D') +ax.plot(dVdt_0D_cm3, label = 'dVdt_0D', linestyle = '--') +ax.set_xlabel('Timestep') +ax.set_ylabel('dVdt (cm^3/s)') +ax.legend() +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'dVdt3D_vs_dVdt0D')) + +# Plot 3D and 0D volume +fig, ax = plt.subplots() +ax.plot(vol_3D_cm3, label = 'V_3D') +ax.plot(vol_0D_cm3, label = 'V_0D', linestyle = '--') +ax.set_xlabel('Timestep') +ax.set_ylabel('Volume (cm^3)') +ax.legend() +#plt.xlim([0,0.4]) +#plt.ylim([-2, 14]) +plt.savefig(os.path.join(sim_folder, 'V3D_vs_V0D')) diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/pv_plot.png b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/pv_plot.png new file mode 100644 index 00000000..7953e628 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/pv_plot.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0cf62cb0d22137f3cb43473ee02c02fc64686b9f3fac79507b6aa1946fd34dd6 +size 21807 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/result_003.vtu b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/result_003.vtu new file mode 100644 index 00000000..23e80e4c --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/result_003.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8cfd6cb09faaa2b69537f9e3c08c3bc47cafa8cf2557f04c0518b1daf3067e9d +size 254881 diff --git a/tests/cases/ustruct/LV_NeoHookean_passive_genBC/svFSI.xml b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/svFSI.xml new file mode 100644 index 00000000..669fee21 --- /dev/null +++ b/tests/cases/ustruct/LV_NeoHookean_passive_genBC/svFSI.xml @@ -0,0 +1,104 @@ + + + + + 0 + 3 + 3 + 1e-2 + 0.50 + STOP_SIM + + 1 + result + + 1 + 1 + + 10 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + 100.0 + + + + + + + true + 3 + 10 + 1e-9 + + 1.0 + 1.0e5 + 0.483333 + + + ST91 + + 0.0 + + + + + true + true + true + true + true + true + true + true + true + + + + FSILS + 1e-9 + 1000 + 50 + + + + + genBC_svFSIplus/genBC.exe + + + + Dirichlet + 0.0 + + + + Neu + Coupled + true + + + + + diff --git a/tests/conftest.py b/tests/conftest.py index 9ea79161..8b6b3f81 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,6 +2,7 @@ import pytest import os +import shutil import subprocess import meshio @@ -49,6 +50,12 @@ def run_by_name(folder, name, t_max, n_proc=1): Returns: Simulation results """ + + # remove old results folders if they exist + dir_path = os.path.join(folder, str(n_proc) + "-procs") + if os.path.exists(dir_path): + shutil.rmtree(dir_path) + # run simulation cmd = " ".join( [ diff --git a/tests/test_struct.py b/tests/test_struct.py index 5376207d..b6658479 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -1,4 +1,6 @@ from .conftest import run_with_reference +import os +import subprocess # Common folder for all tests in this file base_folder = "struct" @@ -31,7 +33,29 @@ def test_block_compression(n_proc): test_folder = "block_compression" run_with_reference(base_folder, test_folder, fields, n_proc) - def test_robin(n_proc): test_folder = "robin" run_with_reference(base_folder, test_folder, fields, n_proc) + +def test_LV_NeoHookean_passive(n_proc): + test_folder = "LV_NeoHookean_passive" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=5) + +def test_LV_NeoHookean_passive_genBC(n_proc): + test_folder = "LV_NeoHookean_passive_genBC" + + # Remove old genBC output + os.chdir(os.path.join("cases", base_folder, test_folder)) + for name in ["AllData", "InitialData", "GenBC.int"]: + if os.path.isfile(name): + os.remove(name) + + # Compile genBC + os.chdir("genBC_svFSIplus") + subprocess.run(["make", "clean"], check=True) + subprocess.run(["make"], check=True) + + # Change back to original directory + os.chdir("../../../../") + + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=3) diff --git a/tests/test_ustruct.py b/tests/test_ustruct.py index b276edeb..f60b853d 100644 --- a/tests/test_ustruct.py +++ b/tests/test_ustruct.py @@ -1,5 +1,6 @@ import os import pytest +import subprocess from .conftest import run_with_reference @@ -32,3 +33,22 @@ def test_tensile_adventitia_HGO(n_proc): def test_LV_Guccione_active(n_proc): test_folder = "LV_Guccione_active" run_with_reference(base_folder, test_folder, fields, n_proc) + +def test_LV_NeoHookean_passive_genBC(n_proc): + test_folder = "LV_NeoHookean_passive_genBC" + + # Remove old genBC output + os.chdir(os.path.join("cases", base_folder, test_folder)) + for name in ["AllData", "InitialData", "GenBC.int"]: + if os.path.isfile(name): + os.remove(name) + + # Compile genBC + os.chdir("genBC_svFSIplus") + subprocess.run(["make", "clean"], check=True) + subprocess.run(["make"], check=True) + + # Change back to original directory + os.chdir("../../../../") + + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=3)