Fixed URIS contribution for more than one valve in fluid equation

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<bool> urisFluidNodeMask;

Code/Source/solver/Parameters.cpp

@@ -2984,10 +2984,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);

Code/Source/solver/Parameters.h

@@ -1783,7 +1783,6 @@ class URISMeshParameters : public ParameterLists
     Parameter<double> thickness; // Thickness of the valve
     Parameter<double> close_thickness; // Thickness of the valve when it is closed
     Parameter<double> resistance; // Resistance of the valve
-    Parameter<double> resistance_close; // Resistance of the valve when it is closed
     Parameter<bool> valve_starts_as_closed; // Whether the valve starts as closed
     Parameter<bool> invert_normal; // Whether to invert the valve surface normal vector
     Parameter<std::string> positive_flow_normal_file_path; // File path for the positive flow normal

Code/Source/solver/distribute.cpp

@@ -345,8 +345,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;
@@ -1178,6 +1176,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);

Code/Source/solver/fluid.cpp

@@ -10,6 +10,7 @@
 #include "nn.h"
 #include "utils.h"
 #include "ris.h"
+#include "uris.h"
 
 #include <array>
 #include <iomanip>
@@ -536,9 +537,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s
 
   // local tangent matrix (for a single element)
   Array3<double> 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<double> ksix(nsd,nsd);
+    // Total resistance factor value of the RIS valves for the current element 
+    // at different quadrature points
+    Vector<double> 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<double> 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<int>(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<double>& Kxi, const Vector<double>& Nw, const Vector<double>& Nq, const Array<double>& Nwx, 
     const Array<double>& Nqx, const Array<double>& Nwxx, const Array<double>& al, const Array<double>& yl, 
-    const Array<double>& bfl, Array<double>& lR, Array3<double>& lK, double K_inverse_darcy_permeability, double DDir)
+    const Array<double>& bfl, Array<double>& lR, Array3<double>& 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<double>& Kxi, const Vector<double>& Nw, const Vector<double>& Nq, const Array<double>& Nwx,
     const Array<double>& Nqx, const Array<double>& Nwxx, const Array<double>& al, const Array<double>& yl,
-    const Array<double>& bfl, Array<double>& lR, Array3<double>& lK, double K_inverse_darcy_permeability, double DDir)
+    const Array<double>& bfl, Array<double>& lR, Array3<double>& 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<double>::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];
   }
 
 }

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<double>& Kxi, 
     const Vector<double>& Nw, const Vector<double>& Nq, const Array<double>& Nwx, const Array<double>& Nqx, 
     const Array<double>& Nwxx, const Array<double>& al, const Array<double>& yl, const Array<double>& bfl, 
-    Array<double>& lR, Array3<double>& lK, double K_inverse_darcy_permeability, double DDir=0.0);
+    Array<double>& lR, Array3<double>& 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<double>& Kxi, 
     const Vector<double>& Nw, const Vector<double>& Nq, const Array<double>& Nwx, const Array<double>& Nqx, 
     const Array<double>& Nwxx, const Array<double>& al, const Array<double>& yl, const Array<double>& bfl, 
-    Array<double>& lR, Array3<double>& lK, double K_inverse_darcy_permeability, double DDir=0.0);
+    Array<double>& lR, Array3<double>& 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);