From 3e1b7a6a7d53aef13d23b281f89b0ae640953ed9 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 20 May 2026 16:56:36 -0700 Subject: [PATCH 1/4] Updated URIS valve test input files and results --- tests/cases/uris/pipe_uris_cfd/result_005.vtu | 4 ++-- tests/cases/uris/pipe_uris_cfd/solver.xml | 2 -- tests/cases/uris/pipe_uris_fsi/result_005.vtu | 4 ++-- tests/cases/uris/pipe_uris_fsi/solver.xml | 2 -- 4 files changed, 4 insertions(+), 8 deletions(-) diff --git a/tests/cases/uris/pipe_uris_cfd/result_005.vtu b/tests/cases/uris/pipe_uris_cfd/result_005.vtu index f8d9f2a0e..feb20dc3a 100644 --- a/tests/cases/uris/pipe_uris_cfd/result_005.vtu +++ b/tests/cases/uris/pipe_uris_cfd/result_005.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d50afef3020ef7541fe422fbe75b5ddaf0df8d06d1815a035ed836836341dda4 -size 1460447 +oid sha256:1908c879949ea025d53905e92bec5807d39a26012e716d1180ba33f61b0bea31 +size 1460399 diff --git a/tests/cases/uris/pipe_uris_cfd/solver.xml b/tests/cases/uris/pipe_uris_cfd/solver.xml index 18f238a91..f57184175 100644 --- a/tests/cases/uris/pipe_uris_cfd/solver.xml +++ b/tests/cases/uris/pipe_uris_cfd/solver.xml @@ -59,7 +59,6 @@ 0.2 0.2 1.0e5 - 1.0e5 false meshes/normal.dat @@ -84,7 +83,6 @@ 0.2 0.2 1.0e5 - 1.0e5 false meshes/normal.dat diff --git a/tests/cases/uris/pipe_uris_fsi/result_005.vtu b/tests/cases/uris/pipe_uris_fsi/result_005.vtu index e49bdc462..95e5b971a 100644 --- a/tests/cases/uris/pipe_uris_fsi/result_005.vtu +++ b/tests/cases/uris/pipe_uris_fsi/result_005.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:c626bca307ffeae844a71a4109d46da37b4e7776739528ffefb071e021a8ae37 -size 411402 +oid sha256:f7cc3986cd696000d86321a643809429510e8e0adcb1f3a6d2d3ec15fc51b23b +size 392425 diff --git a/tests/cases/uris/pipe_uris_fsi/solver.xml b/tests/cases/uris/pipe_uris_fsi/solver.xml index bf9106540..ff0b9bcd6 100644 --- a/tests/cases/uris/pipe_uris_fsi/solver.xml +++ b/tests/cases/uris/pipe_uris_fsi/solver.xml @@ -81,7 +81,6 @@ 0.2 0.2 1.0e5 - 1.0e5 false meshes/normal.dat @@ -106,7 +105,6 @@ 0.2 0.2 1.0e5 - 1.0e5 false meshes/normal.dat From a134f47dd8fe075fb72983f060d68b842daa5397 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 20 May 2026 16:58:28 -0700 Subject: [PATCH 2/4] Removed closed resistance input parameter for URIS valve --- Code/Source/solver/ComMod.h | 9 +++------ Code/Source/solver/Parameters.cpp | 5 ++--- Code/Source/solver/Parameters.h | 1 - Code/Source/solver/distribute.cpp | 3 +-- 4 files changed, 6 insertions(+), 12 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 26410a084..8b41e5c6a 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1482,6 +1482,9 @@ class urisType // if the fully closed position alone is not able to prevent backflow. double sdf_deps_close; + // Resistance value of the valve. + double resistance; + // Whether to invert the valve surface normal vector. Default is false. // // Valve normal vectors are assumed to point downstream, so that the @@ -1622,12 +1625,6 @@ class ComMod { /// @brief Number of URIS surfaces (uninitialized, to be set later) int nUris; - /// @brief URIS resistance - double urisRes; - - /// @brief URIS resistance when the valve is closed - double urisResClose; - /// @brief Fluid-related node mask for URIS SDF. Built once when /// consistent with tnNo; rebuilt automatically if tnNo changes. std::vector urisFluidNodeMask; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 996f155e6..da63f6785 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2998,10 +2998,9 @@ URISMeshParameters::URISMeshParameters() // Parameters under Add_mesh element. // set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); - set_parameter("Thickness", 0.04, !required, thickness); - set_parameter("Closed_thickness", 0.25, !required, close_thickness); + set_parameter("Thickness", 0.2, !required, thickness); + set_parameter("Closed_thickness", 0.2, !required, close_thickness); set_parameter("Resistance", 1.0e5, !required, resistance); - set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Invert_normal", false, !required, invert_normal); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 422cfce55..63b39b369 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1728,7 +1728,6 @@ class URISMeshParameters : public ParameterLists Parameter thickness; // Thickness of the valve Parameter close_thickness; // Thickness of the valve when it is closed Parameter resistance; // Resistance of the valve - Parameter resistance_close; // Resistance of the valve when it is closed Parameter valve_starts_as_closed; // Whether the valve starts as closed Parameter invert_normal; // Whether to invert the valve surface normal vector Parameter positive_flow_normal_file_path; // File path for the positive flow normal diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 5ddc95369..775b3e8a0 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -344,8 +344,6 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.ris0DFlag); cm.bcast(cm_mod, &com_mod.urisFlag); cm.bcast(cm_mod, &com_mod.urisActFlag); - cm.bcast(cm_mod, &com_mod.urisRes); - cm.bcast(cm_mod, &com_mod.urisResClose); cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; @@ -1177,6 +1175,7 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].sdf_default); cm.bcast(cm_mod, &uris[iUris].sdf_deps); cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); + cm.bcast(cm_mod, &uris[iUris].resistance); cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].invert_normal); cm.bcast(cm_mod, &uris[iUris].cnt); From d16987bf3e82face1ea39cc6dcf74709afb10713 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 20 May 2026 17:00:07 -0700 Subject: [PATCH 3/4] Updated URIS valve contribution to fluid equation for more than one valve --- Code/Source/solver/fluid.cpp | 121 +++++++++++++++-------------------- Code/Source/solver/fluid.h | 4 +- Code/Source/solver/fsi.cpp | 69 ++++++++------------ Code/Source/solver/uris.cpp | 73 +++++++++++++++++++-- Code/Source/solver/uris.h | 3 + 5 files changed, 149 insertions(+), 121 deletions(-) diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 2bf56dfb0..a2638315f 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -10,6 +10,7 @@ #include "nn.h" #include "utils.h" #include "ris.h" +#include "uris.h" #include #include @@ -536,9 +537,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s // local tangent matrix (for a single element) Array3 lK(dof*dof,eNoN,eNoN); - - double DDir = 0.0; - + // Loop over all elements of mesh // int num_c = lM.nEl / 10; @@ -620,6 +619,12 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s double Jac{0.0}; Array ksix(nsd,nsd); + // Total resistance factor value of the RIS valves for the current element + // at different quadrature points + Vector risFactorTotalEl; + if (com_mod.urisFlag) { + uris::uris_compute_ris_factor(com_mod, lM, fs[0], e, risFactorTotalEl); + } for (int g = 0; g < fs[0].nG; g++) { #ifdef debug_construct_fluid @@ -650,43 +655,18 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s dmsg << "w: " << w; #endif - // Plot the coordinates of the quad point in the current configuration - if (com_mod.urisFlag) { - Vector distSrf(com_mod.nUris); - distSrf = 0.0; - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a,e); - for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - distSrf(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); - } - } - - DDir = 0.0; - double DDirTmp = 0.0; - double sdf_deps_temp = 0.0; - for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - if (com_mod.uris[iUris].clsFlg) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - } else { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - } - if (distSrf(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} - } - } - - if (!com_mod.urisActFlag) {DDir = 0.0;} - } - // Compute momentum residual and tangent matrix. // if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); + double risFactorTotal = 0.0; + if (com_mod.urisFlag) { + risFactorTotal = risFactorTotalEl(g); + } fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, - Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, + risFactorTotal); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -711,6 +691,14 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s dmsg << "fs[2].lShpF: " << fs[1].lShpF; #endif + // If the number of quadrature points is different for the continuity and + // momentum function spaces, recompute the RIS factor + if (com_mod.urisFlag) { + if (static_cast(risFactorTotalEl.size()) != fs[1].nG) { + uris::uris_compute_ris_factor(com_mod, lM, fs[1], e, risFactorTotalEl); + } + } + for (int g = 0; g < fs[1].nG; g++) { if (g == 0 || !fs[0].lShpF) { auto Nx = fs[0].Nx.rslice(g); @@ -736,12 +724,19 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); - fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); + double risFactorTotal = 0.0; + if (com_mod.urisFlag) { + risFactorTotal = risFactorTotalEl(g); + } + fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, + risFactorTotal); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); - fluid_2d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability); + fluid_2d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability); } } // g: loop @@ -1443,7 +1438,8 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + double risFactorTotal) { #define n_debug_fluid3d_c #ifdef debug_fluid3d_c @@ -1469,14 +1465,6 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e const double ctM = 1.0; const double ctC = 36.0; - double Res; - if (!com_mod.urisFlag) { - Res = 0.0; - } else { - Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} - } - double rho = dmn.prop[PhysicalProperyType::fluid_density]; double f[3]; f[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -1666,7 +1654,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param - kT = kT + pow(Res*DDir, 2.0); + kT = kT + pow(risFactorTotal, 2.0); double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) @@ -1694,11 +1682,11 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2]); up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0] - + (Res*DDir)*u[0]); + + risFactorTotal*u[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1] - + (Res*DDir)*u[1]); + + risFactorTotal*u[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2] - + (Res*DDir)*u[2]); + + risFactorTotal*u[2]); for (int a = 0; a < eNoNw; a++) { double uNx = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); @@ -1707,7 +1695,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a) - - (Res*DDir)*Nw(a); + - risFactorTotal*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -1775,7 +1763,8 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + double risFactorTotal) { #define n_debug_fluid_3d_m #ifdef debug_fluid_3d_m @@ -1799,14 +1788,6 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double Res; - if (!com_mod.urisFlag) { - Res = 0.0; - } else { - Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} - } - double ctM = 1.0; double ctC = 36.0; @@ -2019,7 +2000,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param - kT = kT + pow(Res*DDir, 2.0); + kT = kT + pow(risFactorTotal, 2.0); double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) @@ -2054,11 +2035,11 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2]); up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0] - + (Res*DDir)*u[0]); + + risFactorTotal * u[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1] - + (Res*DDir)*u[1]); + + risFactorTotal * u[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2] - + (Res*DDir)*u[2]); + + risFactorTotal * u[2]); double tauC, tauB, pa; double eps = std::numeric_limits::epsilon(); @@ -2140,7 +2121,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a) - - (Res*DDir)*Nw(a); + - risFactorTotal*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -2177,7 +2158,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(0,a,b) = lK(0,a,b) + wl*(T2 + T1); // lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + risFactorTotal*wl*Nw(b)*Nw(a); // dRm_a1/du_b2 T2 = mu*rM[1][0] + tauC*rM[0][1] + esNx[0][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][0][b]; @@ -2196,7 +2177,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(5,a,b) = lK(5,a,b) + wl*(T2 + T1); // lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + risFactorTotal*wl*Nw(b)*Nw(a); // dRm_a2/du_b3 T2 = mu*rM[2][1] + tauC*rM[1][2] + esNx[1][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][1][b]; @@ -2215,7 +2196,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(10,a,b) = lK(10,a,b) + wl*(T2 + T1); // lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + risFactorTotal*wl*Nw(b)*Nw(a); //dmsg << "lK(10,a,b): " << lK(10,a,b); } } @@ -2241,11 +2222,11 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Local residue for (int a = 0; a < eNoNw; a++) { lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) - + Res*DDir*w*Nw(a)*u[0]; + + risFactorTotal*w*Nw(a)*u[0]; lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) - + Res*DDir*w*Nw(a)*u[1]; + + risFactorTotal*w*Nw(a)*u[1]; lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) - + Res*DDir*w*Nw(a)*u[2]; + + risFactorTotal*w*Nw(a)*u[2]; } } diff --git a/Code/Source/solver/fluid.h b/Code/Source/solver/fluid.h index f4bf26abe..02fab7de3 100644 --- a/Code/Source/solver/fluid.h +++ b/Code/Source/solver/fluid.h @@ -36,12 +36,12 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double risFactorTotal); void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double risFactorTotal); void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, double& mu, double& mu_s, double& mu_x); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 913827a26..c74156410 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -13,6 +13,7 @@ #include "ustruct.h" #include "utils.h" #include "ris.h" +#include "uris.h" #include #include @@ -80,7 +81,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So // double struct_3d_time = 0.0; double fluid_3d_time = 0.0; - double DDir = 0.0; for (int e = 0; e < lM.nEl; e++) { // setting globals @@ -165,6 +165,12 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So // double Jac{0.0}; Array ksix(nsd,nsd); + // Total resistance factor value of the RIS valves for the current element + // at different quadrature points + Vector risFactorTotalEl; + if (com_mod.urisFlag) { + uris::uris_compute_ris_factor(com_mod, lM, fs_1[0], e, risFactorTotalEl); + } for (int g = 0; g < fs_1[0].nG; g++) { if (g == 0 || !fs_1[1].lShpF) { @@ -188,54 +194,18 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So double w = fs_1[0].w(g) * Jac; - // Plot the coordinates of the quad point in the current configuration - if (com_mod.urisFlag) { - Vector distSrf(com_mod.nUris); - distSrf = 0.0; - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a,e); - for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - distSrf(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); - } - } - - DDir = 0.0; - double sdf_deps_temp = 0; - double DDirTmp = 0.0; - for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - // if (distSrf(iUris) <= com_mod.uris[iUris].sdf_deps) { - // DDirTmp = (1 + cos(pi*distSrf(iUris)/com_mod.uris[iUris].sdf_deps))/ - // (2*com_mod.uris[iUris].sdf_deps*com_mod.uris[iUris].sdf_deps); - // if (DDirTmp > DDir) {DDir = DDirTmp;} - // } - - if (com_mod.uris[iUris].clsFlg) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - } else { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - } - if (distSrf(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} - } - } - - if (!com_mod.urisActFlag) {DDir = 0.0;} - - // std::cout << "===== DDir: " << DDir << std::endl; - } - - if (nsd == 3) { switch (cPhys) { case Equation_fluid: { auto N0 = fs_1[0].N.col(g); auto N1 = fs_1[1].N.col(g); + double risFactorTotal = 0.0; + if (com_mod.urisFlag) { + risFactorTotal = risFactorTotalEl(g); + } // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman - // fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); + fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, risFactorTotal); } break; @@ -283,6 +253,14 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } } // g: loop + // If the number of quadrature points is different for the continuity and + // momentum function spaces, recompute the RIS factor + if (com_mod.urisFlag) { + if (static_cast(risFactorTotalEl.size()) != fs_2[1].nG) { + uris::uris_compute_ris_factor(com_mod, lM, fs_2[1], e, risFactorTotalEl); + } + } + // Gauss integration 2 // for (int g = 0; g < fs_2[1].nG; g++) { @@ -310,10 +288,13 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_fluid: { auto N0 = fs_2[0].N.col(g); auto N1 = fs_2[1].N.col(g); + double risFactorTotal = 0.0; + if (com_mod.urisFlag) { + risFactorTotal = risFactorTotalEl(g); + } // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman - //fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); + fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, risFactorTotal); } break; case Equation_ustruct: diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 714c03511..35e738894 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -4,6 +4,7 @@ #include "uris.h" #include "all_fun.h" +#include "consts.h" #include "lhsa.h" #include "mat_fun.h" #include "nn.h" @@ -557,11 +558,6 @@ void uris_read_msh(Simulation* simulation) { com_mod.urisActFlag = true; auto param = simulation->parameters.URIS_mesh_parameters[0]; - com_mod.urisRes = param->resistance(); - com_mod.urisResClose = param->resistance_close(); - - std::cout << "URIS resistance: " << com_mod.urisRes << std::endl; - std::cout << "URIS resistance when the valve is closed: " << com_mod.urisResClose << std::endl; int nUris = simulation->parameters.URIS_mesh_parameters.size(); com_mod.nUris = nUris; @@ -603,6 +599,7 @@ void uris_read_msh(Simulation* simulation) { // Use large default value for the signed distance function to indicate that // the fluid node is far away from the valve. uris_obj.sdf_default = param->close_thickness() * 1e6; + uris_obj.resistance = param->resistance(); uris_obj.clsFlg = param->valve_starts_as_closed(); uris_obj.invert_normal = param->invert_normal(); @@ -1286,4 +1283,70 @@ double uris_compute_sdf_sign(const urisType& uris_obj, const Vector& xp, return (dotP < 0.0) ? -1.0 : 1.0; } +/// @brief Compute the RIS factor for the current element at different quadrature points +void uris_compute_ris_factor(const ComMod& com_mod, const mshType& lM, const fsType& fs, + const int e, Vector& ris_factor_total_el) { + #define n_dbg_uris_compute_ris_factor + #ifdef dbg_uris_compute_ris_factor + DebugMsg dmsg(__func__, 0); + dmsg.banner(); + dmsg << "computing RIS factor"; + #endif + + const int eNoN = lM.eNoN; + const int nUris = com_mod.nUris; + ris_factor_total_el.resize(fs.nG); + ris_factor_total_el = 0.0; + + for (int g = 0; g < fs.nG; g++) { + Vector dist_srf(nUris); + dist_srf = 0.0; + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + for (int iUris = 0; iUris < nUris; iUris++) { + dist_srf(iUris) += fs.N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + } + } + + double sdf_deps; + double delta_eps; + for (int iUris = 0; iUris < nUris; iUris++) { + sdf_deps = 0.0; + delta_eps = 0.0; + double start_deps, end_deps; + int n_steps; + if (com_mod.uris[iUris].clsFlg) { + start_deps = com_mod.uris[iUris].sdf_deps; + end_deps = com_mod.uris[iUris].sdf_deps_close; + n_steps = com_mod.uris[iUris].DxClose.nrows(); + } else { + start_deps = com_mod.uris[iUris].sdf_deps_close; + end_deps = com_mod.uris[iUris].sdf_deps; + n_steps = com_mod.uris[iUris].DxOpen.nrows(); + } + + if (n_steps <= 0) { + sdf_deps = end_deps; + } else if (com_mod.uris[iUris].cnt >= n_steps) { + sdf_deps = end_deps; + } else if (com_mod.uris[iUris].cnt <= 0) { + sdf_deps = start_deps; + } else { + // Linear ramping: start_deps -> end_deps over n_steps. + // This alleviates the sudden change in the resistance factor when the + // valve status changes from open to close and reduces oscillations + double progress = static_cast(com_mod.uris[iUris].cnt) / static_cast(n_steps); + sdf_deps = start_deps + progress * (end_deps - start_deps); + } + + if (dist_srf(iUris) < sdf_deps && sdf_deps > 0.0) { + delta_eps = (1 + cos(consts::pi*dist_srf(iUris)/sdf_deps))/(2*sdf_deps*sdf_deps); + } + ris_factor_total_el(g) += com_mod.uris[iUris].resistance * delta_eps; + } // iUris: loop + + if (!com_mod.urisActFlag) {ris_factor_total_el(g) = 0.0;} + } +} + } \ No newline at end of file diff --git a/Code/Source/solver/uris.h b/Code/Source/solver/uris.h index d44428bc2..b001d44d6 100644 --- a/Code/Source/solver/uris.h +++ b/Code/Source/solver/uris.h @@ -46,6 +46,9 @@ double uris_compute_sdf_sign(const urisType& uris_obj, const Vector& xp, void uris_build_fluid_node_mask(ComMod& com_mod); +void uris_compute_ris_factor(const ComMod& com_mod, const mshType& lM, const fsType& fs, const int e, + Vector& ris_factor_total_el); + } #endif From c5f1454190121fe366a94b3f95abcbbe603cdd89 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Thu, 21 May 2026 13:53:21 -0700 Subject: [PATCH 4/4] Updated RIS factor computation function with safe number of shape functions --- Code/Source/solver/uris.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 35e738894..109060fe1 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -557,8 +557,6 @@ void uris_read_msh(Simulation* simulation) { com_mod.urisFlag = true; com_mod.urisActFlag = true; - auto param = simulation->parameters.URIS_mesh_parameters[0]; - int nUris = simulation->parameters.URIS_mesh_parameters.size(); com_mod.nUris = nUris; @@ -1293,15 +1291,14 @@ void uris_compute_ris_factor(const ComMod& com_mod, const mshType& lM, const fsT dmsg << "computing RIS factor"; #endif - const int eNoN = lM.eNoN; const int nUris = com_mod.nUris; ris_factor_total_el.resize(fs.nG); ris_factor_total_el = 0.0; + Vector dist_srf(nUris); for (int g = 0; g < fs.nG; g++) { - Vector dist_srf(nUris); dist_srf = 0.0; - for (int a = 0; a < eNoN; a++) { + for (int a = 0; a < fs.eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < nUris; iUris++) { dist_srf(iUris) += fs.N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac));