Add velocity limiter

Author: shkodm

models/multiphase/d3q27_pf_velocity/Dynamics.R

@@ -259,6 +259,7 @@ if (Options$thermo){
 	AddSetting(name="BuoyancyX", default=0.0, comment='applied (rho_h-rho)*BuoyancyX')
 	AddSetting(name="BuoyancyY", default=0.0, comment='applied (rho_h-rho)*BuoyancyY')
 	AddSetting(name="BuoyancyZ", default=0.0, comment='applied (rho_h-rho)*BuoyancyZ')
+    AddSetting(name="vlimiter", default=0.0, comment="Velocity limiter value")
 ##################################
 ########TRACKING VARIABLES########
 ##################################
@@ -321,6 +322,7 @@ if (Options$thermo){
 	AddGlobal(name="LiqTotalVelocityX", comment='use to determine avg velocity of droplets', unit="m/s")
 	AddGlobal(name="LiqTotalVelocityY", comment='use to determine avg velocity of droplets', unit="m/s")
 	AddGlobal(name="LiqTotalVelocityZ", comment='use to determine avg velocity of droplets', unit="m/s")
+	AddGlobal(name="LimitedCells", comment='Number of cells where we limited the velocity', unit="1")
     AddGlobal(name="NumFluidCells", comment='Number of fluid cells')
     AddGlobal(name="NumSpecialPoints", comment='Number of special points')
     AddGlobal(name="NumWallBoundaryPoints", comment='Number of boundary nodes')

models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt

@@ -687,6 +687,18 @@ CudaDeviceFunction void CollisionMRT()
 
 	<?R 	C( u, m0[2:4] + 0.5 * F_total * rho.inv) ?>
 	} while (j<force_fixed_iterator);
+
+   if (vlimiter != 0) {
+        real_t vl = U*U + V*V + W*W;
+        if (vl > vlimiter*vlimiter){
+            real_t reducer = vlimiter / sqrt(vl);
+            AddToLimitedCells(1);
+            U = U * reducer;
+            V = V * reducer;
+            W = W * reducer;
+        }
+    }
+
 	// PHASE FIELD POPULATION UPDATE:
 	tmp1 = (1.0 - 4.0*(C - 0.5)*(C - 0.5))/IntWidth;
 	<?R