Fixing slow nonlinear convergence for URIS valve

Code/Source/solver/ComMod.h

@@ -1209,6 +1209,9 @@ class eqType
     /// @brief URIS: Outputs
     std::vector<outputType> outURIS;
 
+    /// @brief Explicit geometry coupling
+    bool expl_geom_cpl = false;
+
     /// @brief Body force associated with this equation
     std::vector<bfType> bf;
 };

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);
 
 }
 

Code/Source/solver/Parameters.h

@@ -1378,6 +1378,11 @@ class EquationParameters : public ParameterLists
     Parameter<std::string> type;
     Parameter<bool> use_taylor_hood_type_basis;
 
+    // Explicit geometric coupling for FSI simulations: the fluid-structure equations
+    // are solved to convergence using the mesh displacement from the previous time step,
+    // and only then is the mesh equation solved.
+    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;
 

Code/Source/solver/distribute.cpp

@@ -1470,6 +1470,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;

Code/Source/solver/fluid.cpp

@@ -2239,11 +2239,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]+up[0]);
+                        + Res*DDir*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]+up[1]);
+                        + Res*DDir*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]+up[2]);
+                        + Res*DDir*w*Nw(a)*u[2];
   }
 
 }

Code/Source/solver/pic.cpp

@@ -143,7 +143,12 @@ void picc(Simulation* simulation)
     pic_eth(simulation);
   }
 
-  if (eq.phys == Equation_FSI) {
+  // If explicit geometric coupling is not used, update the acceleration, 
+  // velocity, and displacement for the FSI equation as usual
+  if (eq.phys == Equation_FSI && !eq.expl_geom_cpl) {
+    if (com_mod.eq[1].phys != Equation_mesh) {
+      throw std::runtime_error("PICC: FSI simulation requires a mesh motion equation as the second equation.");
+    }
     int s = com_mod.eq[1].s;
     int e = com_mod.eq[1].e;
     #ifdef debug_picc
@@ -263,30 +268,62 @@ void picc(Simulation* simulation)
   //if (ALL(eq.ok)) RETURN
 
   if (eq.coupled) {
-    cEq = cEq + 1;
-    #ifdef debug_picc
-    dmsg << "eq " << " coupled ";
-    dmsg << "1st update cEq: " << cEq;
-    #endif
+    if (!eq.expl_geom_cpl) {
+      // For coupled equations, if explicit geometric coupling is not used, 
+      // increment the equation counter after each linear solve
+      cEq = cEq + 1;
+      #ifdef debug_picc
+      dmsg << "eq " << " coupled ";
+      dmsg << "1st update cEq: " << cEq;
+      #endif
+
+      auto& eqs = com_mod.eq;
+      if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return !eq.coupled || eq.ok;}) == eqs.size()) {
+        while (cEq < com_mod.nEq) {
+          if (!eqs[cEq].coupled) {
+            break;
+          }
+          cEq = cEq + 1;
+        }
 
-    auto& eqs = com_mod.eq;
-    if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return !eq.coupled || eq.ok;}) == eqs.size()) {
-      while (cEq < com_mod.nEq) {
-        if (!eqs[cEq].coupled) {
-          break;
+      } else {
+        if (cEq >= com_mod.nEq) {
+          cEq = 0;
         }
-        cEq = cEq + 1;
-      }
 
-    } else {
-      if (cEq >= com_mod.nEq) {
-        cEq = 0;
+        while (!eqs[cEq].coupled) {
+          cEq = cEq + 1;
+          if (cEq >= com_mod.nEq) {
+            cEq = 0;
+          }
+        }
       }
-
-      while (!eqs[cEq].coupled) {
+    } else {
+      // If explicit geometric coupling is used for coupled equations, 
+      // only update the equation counter if the current equation is converged
+      if (eq.ok) {
         cEq = cEq + 1;
         if (cEq >= com_mod.nEq) {
-          cEq = 0;
+            cEq = 0;
+        }
+      }
+      // Update the acceleration, velocity, and displacement when the FSI equation is ok
+      if (eq.ok && eq.phys == Equation_FSI) {
+        if (com_mod.eq[1].phys != Equation_mesh) {
+          throw std::runtime_error("PICC: FSI simulation requires a mesh motion equation as the second equation.");
+        }
+        int s = com_mod.eq[1].s;
+        int e = com_mod.eq[1].e;
+        for (int Ac = 0; Ac < tnNo; Ac++) {
+          if (all_fun::is_domain(com_mod, eq, Ac, Equation_struct) || 
+              all_fun::is_domain(com_mod, eq, Ac, Equation_ustruct) || 
+              all_fun::is_domain(com_mod, eq, Ac, Equation_lElas)) {
+            for (int i = 0; i < e-s+1; i++) {
+              An(i+s,Ac) = An(i,Ac);
+              Yn(i+s,Ac) = Yn(i,Ac);
+              Dn(i+s,Ac) = Dn(i,Ac);
+            }
+          }
         }
       }
     }

Code/Source/solver/read_files.cpp

@@ -1381,6 +1381,7 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq)
   lEq.minItr = eq_params->min_iterations.value();
   lEq.maxItr = eq_params->max_iterations.value();
   lEq.tol = eq_params->tolerance.value();
+  lEq.expl_geom_cpl = eq_params->explicit_geometric_coupling.value();
 
   // Initialize coupled BC.
   //
@@ -1459,6 +1460,13 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq)
 
   set_equation_properties(simulation, eq_params, lEq, propL, outPuts, nDOP);
 
+  // Check if explicit geometric coupling is allowed for the equation
+  if (lEq.expl_geom_cpl) {
+    if (lEq.phys != EquationType::phys_FSI) {
+      throw std::runtime_error("Explicit geometric coupling is only allowed for FSI equation.");
+    }
+  }
+
   // Read VTK files or boundaries. [TODO:DaveP] this is not a correct comment.
   read_outputs(simulation, eq_params, lEq, nDOP, outPuts);
 
@@ -1796,7 +1804,9 @@ void read_files(Simulation* simulation, const std::string& file_name)
     if (eq.phys == EquationType::phys_mesh) {   
       if (!com_mod.mvMsh) {
         throw std::runtime_error("mesh equation can only be specified after FSI equation");
-      }     
+      }
+      // Use the explicit geometry coupling flag of the FSI equation.
+      eq.expl_geom_cpl = com_mod.eq[0].expl_geom_cpl; 
     }     
   }
   #ifdef debug_read_files