except the standing wave

src/shared/particle_dynamics/fluid_dynamics/fluid_integration.hpp

@@ -88,7 +88,7 @@ void Integration1stHalf<Inner<>, RiemannSolverType, KernelCorrectionType>::inter
         rho_dissipation += riemann_solver_.DissipativeUJump(p_[index_i] - p_[index_j]) * dW_ijV_j;
     }
     force_[index_i] += force * Vol_[index_i];
-    drho_dt_[index_i] = rho_dissipation * mass_[index_i] / Vol_[index_i];
+    drho_dt_[index_i] = rho_dissipation * rho_[index_i];
 }
 //=================================================================================================//
 template <class RiemannSolverType, class KernelCorrectionType>
@@ -121,7 +121,7 @@ void Integration1stHalf<Contact<Wall>, RiemannSolverType, KernelCorrectionType>:
         }
     }
     force_[index_i] += force * Vol_[index_i];
-    drho_dt_[index_i] += rho_dissipation * mass_[index_i] / Vol_[index_i];
+    drho_dt_[index_i] += rho_dissipation * rho_[index_i];
 }
 //=================================================================================================//
 template <class RiemannSolverType, class KernelCorrectionType>
@@ -165,7 +165,7 @@ void Integration1stHalf<Contact<>, RiemannSolverType, KernelCorrectionType>::
         }
     }
     this->force_[index_i] += force * this->Vol_[index_i];
-    this->drho_dt_[index_i] += rho_dissipation * this->mass_[index_i] / this->Vol_[index_i];
+    this->drho_dt_[index_i] += rho_dissipation * this->rho_[index_i];
 }
 //=================================================================================================//
 template <class RiemannSolverType>
@@ -203,7 +203,7 @@ void Integration2ndHalf<Inner<>, RiemannSolverType>::interaction(size_t index_i,
         density_change_rate += u_jump * dW_ijV_j;
         p_dissipation += riemann_solver_.DissipativePJump(u_jump) * dW_ijV_j * e_ij;
     }
-    drho_dt_[index_i] += density_change_rate * mass_[index_i] / Vol_[index_i];
+    drho_dt_[index_i] += density_change_rate * rho_[index_i];
     force_[index_i] = p_dissipation * Vol_[index_i];
 };
 //=================================================================================================//
@@ -236,7 +236,7 @@ void Integration2ndHalf<Contact<Wall>, RiemannSolverType>::interaction(size_t in
             p_dissipation += riemann_solver_.DissipativePJump(u_jump) * dW_ijV_j * n_k[index_j];
         }
     }
-    drho_dt_[index_i] += density_change_rate * this->mass_[index_i] / this->Vol_[index_i];
+    drho_dt_[index_i] += density_change_rate * this->rho_[index_i];
     force_[index_i] += p_dissipation * this->Vol_[index_i];
 }
 //=================================================================================================//
@@ -277,7 +277,7 @@ void Integration2ndHalf<Contact<>, RiemannSolverType>::interaction(size_t index_
             p_dissipation += riemann_solver_k.DissipativePJump(u_jump) * dW_ijV_j * e_ij;
         }
     }
-    this->drho_dt_[index_i] += density_change_rate * this->mass_[index_i] / this->Vol_[index_i];
+    this->drho_dt_[index_i] += density_change_rate * this->rho_[index_i];
     this->force_[index_i] += p_dissipation * this->Vol_[index_i];
 }
 //=================================================================================================//

src/shared/particle_dynamics/general_dynamics/kernel_correction.cpp

@@ -20,7 +20,7 @@ void LinearGradientCorrectionMatrix<Inner<>>::interaction(size_t index_i, Real d
 //=================================================================================================//
 void LinearGradientCorrectionMatrix<Inner<>>::update(size_t index_i, Real dt)
 {
-    Real det_sqr = SMAX(alpha_ - B_[index_i].determinant(), 0.0);
+    Real det_sqr = SMAX(alpha_ - B_[index_i].determinant(), Real(0));
     Matd inverse = B_[index_i].inverse();
     Real weight1_ = B_[index_i].determinant() / (B_[index_i].determinant() + det_sqr);
     Real weight2_ = det_sqr / (B_[index_i].determinant() + det_sqr);

tests/2d_examples/test_2d_impact_patch/impact_patch.cpp

@@ -100,7 +100,6 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     BoundingBox system_domain_bounds(Vec2d(-LH, -LH), Vec2d(LL + BW, LH + BW));
     SPHSystem sph_system(system_domain_bounds, particle_spacing_ref);
-    sph_system.setGenerateRegressionData(true);
     sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment();
     //----------------------------------------------------------------------
     //	Creating bodies with corresponding materials and particles.

tests/2d_examples/test_2d_rotation_patch/rotation_patch.cpp

@@ -114,7 +114,6 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     BoundingBox system_domain_bounds(Vec2d(-LL, -LL), Vec2d(LL + 10 * BW, LH + 10 * BW));
     SPHSystem sph_system(system_domain_bounds, particle_spacing_ref);
-    //sph_system.setGenerateRegressionData(true);
     sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment();
     //----------------------------------------------------------------------
     //	Creating bodies with corresponding materials and particles.

tests/2d_examples/test_2d_standing_wave/standing_wave.cpp

@@ -112,7 +112,6 @@ int main(int ac, char *av[])
     sph_system.setRunParticleRelaxation(false);
     /** Tag for computation start with relaxed body fitted particles distribution. */
     sph_system.setReloadParticles(false);
-    //sph_system.setGenerateRegressionData(true);
     sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment();
     //----------------------------------------------------------------------
     //	Creating bodies with corresponding materials and particles.
@@ -271,7 +270,7 @@ int main(int ac, char *av[])
         {
             /** outer loop for dual-time criteria time-stepping. */
             time_instance = TickCount::now();
-            Real advection_dt = fluid_advection_time_step.exec();
+            Real advection_dt = 0.3 * fluid_advection_time_step.exec();
             fluid_density_by_summation.exec();
             corrected_configuration_fluid.exec();
             interval_computing_time_step += TickCount::now() - time_instance;