Updates for URIS valve with nonzero valve velocity and improved nonlinear convergence

Code/Source/solver/ComMod.h

@@ -1209,6 +1209,9 @@ class eqType
     /// @brief URIS: Outputs
     std::vector<outputType> outURIS;
 
+    /// @brief Explicit geometry coupling for FSI problem
+    bool expl_geom_cpl = false;
+
     /// @brief Body force associated with this equation
     std::vector<bfType> bf;
 };
@@ -1459,17 +1462,38 @@ class urisType
     // Position coordinates (2D array: rows x columns).
     Array<double> x;
 
+    // Position coordinates of last time step.
+    Array<double> x_prev;
+
+    // Time derivative of the position coordinates (2D array).
+    Array<double> v;
+
     // Displacement (new) (2D array).
     Array<double> Yd;
 
     // Default signed distance value away from the valve.
     double sdf_default = 1000.0;
 
-    // Default distance value of the valve boundary (valve thickness).
-    double sdf_deps = 0.25;
+    // Default distance value of the valve boundary (half valve thickness).
+    double sdf_deps = 0.2;
+
+    // Default half valve thickness when valve is closed.
+    double sdf_deps_close = 0.2;
 
-    // Default distance value of the valve boundary when the valve is closed.
-    double sdf_deps_close = 0.25;
+    // Default half scaffold thickness.
+    double sdf_deps_scaffold = 0.1;
+
+    /// @brief URIS resistance
+    double resistance;
+
+    /// @brief URIS resistance when the valve is closed
+    double resistance_close;
+
+    // Whether to use the valve velocity.
+    bool use_valve_velocity = false;
+
+    // Reverse the surface normal vector
+    bool reverse_normal = false;
 
     // Displacements of the valve when it opens (3D array).
     Array3<double> DxOpen;
@@ -1486,9 +1510,18 @@ class urisType
     // Iteration count.
     int cnt = 1000000;
 
-    // URIS: signed distance function of each node to the uris (1D array).
+    // pressurization time
+    double pressurization_time = 0.0;
+
+    // URIS: signed distance function of each node to the uris valves (1D array).
+    Vector<double> dirac_delta_func;
+
+    // URIS: signed distance function of each node to the uris valves (1D array).
     Vector<double> sdf;
 
+    // URIS: time derivative of the signed distance field       
+    Array<double> sdf_t;
+
     // Mesh scale factor.
     double scF = 1.0;
 
@@ -1501,16 +1534,22 @@ class urisType
     // Relaxation factor to compute weighted averages of pressure values.
     double relax_factor = 0.5;
 
-    // Array to store the fluid mesh elements that the uris node is in (2D array).
-    Array<int> elemId;
-
     // Array to count how many times a uris node is found in the fluid mesh of a processor (1D array).
     Vector<int> elemCounter;
 
+    // Array to store the fluid mesh elements that the uris node is in (2D array).
+    Array<int> elemId;
+
     // Derived type variables
     // IB meshes
     std::vector<mshType> msh;
 
+    // Scaffold mesh
+    bool scaffold_flag = false;
+    mshType scaffold_mesh;
+
+    Vector<double> sdf_scaffold;
+
 };
 
 /// @brief The ComMod class duplicates the data structures in the Fortran COMMOD module
@@ -1598,12 +1637,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 Whether to use precomputed state-variable solutions
     bool usePrecomp = false;
     //----- int members -----//

Code/Source/solver/Parameters.cpp

@@ -2154,6 +2154,7 @@ EquationParameters::EquationParameters()
 
   set_parameter("Tolerance", 0.5, !required, tolerance);
   set_parameter("Use_taylor_hood_type_basis", false, !required, use_taylor_hood_type_basis);
+  set_parameter("Explicit_geometric_coupling", false, !required, explicit_geometric_coupling);
 
 }
 
@@ -2828,12 +2829,16 @@ 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("Thickness", 0.25,  !required, thickness);
   set_parameter("Closed_thickness", 0.25,  !required, close_thickness);
   set_parameter("Resistance", 1.0e5,  !required, resistance);
   set_parameter("Closed_resistance", 1.0e5,  !required, resistance_close);
+  set_parameter("Use_valve_velocity", false,  !required, use_valve_velocity);
   set_parameter("Valve_starts_as_closed", true,  !required, valve_starts_as_closed);
   set_parameter("Positive_flow_normal_file_path", "",  !required, positive_flow_normal_file_path);
+  set_parameter("Scaffold_file_path", "",  !required, scaffold_file_path);
+  set_parameter("Pressurization_time", 0.0,  !required, pressurization_time);
+  set_parameter("Reverse_surface_normal", false,  !required, reverse_surface_normal);
 }
 
 void URISMeshParameters::print_parameters()
@@ -2948,7 +2953,6 @@ void URISFaceParameters::set_values(tinyxml2::XMLElement* face_elem)
   }
 }
 
-
 //////////////////////////////////////////////////////////
 //                 LinearAlgebraParameters              //
 //////////////////////////////////////////////////////////

Code/Source/solver/Parameters.h

@@ -1377,6 +1377,7 @@ class EquationParameters : public ParameterLists
 
     Parameter<std::string> type;
     Parameter<bool> use_taylor_hood_type_basis;
+    Parameter<bool> explicit_geometric_coupling;
 
     // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman
     Parameter<double> inverse_darcy_permeability;
@@ -1648,9 +1649,12 @@ class URISMeshParameters : public ParameterLists
     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> use_valve_velocity; // Whether to use the valve velocity
     Parameter<bool> valve_starts_as_closed; // Whether the valve starts as closed
     Parameter<std::string> positive_flow_normal_file_path; // File path for the positive flow normal
-
+    Parameter<std::string> scaffold_file_path; // File path for the scaffold
+    Parameter<double> pressurization_time; // Time for pressurization
+    Parameter<bool> reverse_surface_normal; // Whether to reverse the surface normal vector
 };
 
 

Code/Source/solver/distribute.cpp

@@ -343,8 +343,7 @@ 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;
@@ -1169,10 +1168,42 @@ 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].sdf_deps_scaffold);
+    cm.bcast(cm_mod, &uris[iUris].resistance);
+    cm.bcast(cm_mod, &uris[iUris].resistance_close);
+    cm.bcast(cm_mod, &uris[iUris].pressurization_time);
     cm.bcast(cm_mod, &uris[iUris].clsFlg);
     cm.bcast(cm_mod, &uris[iUris].cnt);
     cm.bcast(cm_mod, &uris[iUris].scF);
     cm.bcast(cm_mod, uris[iUris].nrm);
+
+    cm.bcast(cm_mod, &uris[iUris].scaffold_flag);
+    cm.bcast(cm_mod, &uris[iUris].reverse_normal);
+    cm.bcast(cm_mod, &uris[iUris].use_valve_velocity);
+    if (uris[iUris].scaffold_flag) {
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.lShl);
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.nEl);
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.gnEl);
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.eNoN);
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.nNo);
+      cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.gnNo);
+    }
+  }
+
+  if (cm.slv(cm_mod)) {
+    for (int iUris = 0; iUris < com_mod.nUris; iUris++) {
+      if (uris[iUris].scaffold_flag) {
+        uris[iUris].scaffold_mesh.x.resize(com_mod.nsd, uris[iUris].scaffold_mesh.gnNo);
+        uris[iUris].scaffold_mesh.IEN.resize(uris[iUris].scaffold_mesh.eNoN, uris[iUris].scaffold_mesh.gnEl);
+      }
+    }
+  }
+
+  for (int iUris = 0; iUris < com_mod.nUris; iUris++) {
+    if (uris[iUris].scaffold_flag) {
+      cm.bcast(cm_mod, uris[iUris].scaffold_mesh.x);
+      cm.bcast(cm_mod, uris[iUris].scaffold_mesh.IEN);
+    }
   }
 
   std::vector<Vector<int>> lM_gN_flat(com_mod.nUris);
@@ -1328,6 +1359,8 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) {
       uris[iUris].DxOpen.resize(OpenN(iUris), com_mod.nsd, uris[iUris].tnNo);
       uris[iUris].DxClose.resize(CloseN(iUris), com_mod.nsd, uris[iUris].tnNo);
       uris[iUris].x.resize(com_mod.nsd, uris[iUris].tnNo);
+      uris[iUris].x_prev.resize(com_mod.nsd, uris[iUris].tnNo);
+      uris[iUris].v.resize(com_mod.nsd, uris[iUris].tnNo);
       uris[iUris].Yd.resize(com_mod.nsd, uris[iUris].tnNo);
       DxOpenFlat[iUris].resize(DxOpenFlatSize(iUris));
       DxCloseFlat[iUris].resize(DxCloseFlatSize(iUris));
@@ -1336,6 +1369,8 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) {
 
   for (int iUris = 0; iUris < com_mod.nUris; iUris++) {
     cm.bcast(cm_mod, uris[iUris].x);
+    cm.bcast(cm_mod, uris[iUris].x_prev);
+    cm.bcast(cm_mod, uris[iUris].v);
     cm.bcast(cm_mod, uris[iUris].Yd);
     cm.bcast(cm_mod, DxOpenFlat[iUris]);
     cm.bcast(cm_mod, DxCloseFlat[iUris]);
@@ -1470,6 +1505,7 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std::
   cm.bcast(cm_mod, &lEq.tol);
   cm.bcast(cm_mod, &lEq.useTLS);
   cm.bcast(cm_mod, &lEq.assmTLS);
+  cm.bcast(cm_mod, &lEq.expl_geom_cpl);
 
   #ifdef dist_eq
   dmsg << "lEq.nOutput: " << lEq.nOutput;