diff --git a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp index cc65c8f7e5..7e8f79cc82 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/viscous_dynamics.hpp @@ -37,11 +37,12 @@ void ViscousForce, ViscosityType, KernelCorrectionType>::interaction(siz for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; + const Vecd &e_ij = inner_neighborhood.e_ij_[n]; // viscous force Vecd vel_derivative = (vel_[index_i] - vel_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); - force += (kernel_correction_(index_i) + kernel_correction_(index_j)) * + force += e_ij.dot((kernel_correction_(index_i) + kernel_correction_(index_j)) * e_ij) * mu_(index_i, index_j) * vel_derivative * inner_neighborhood.dW_ij_[n] * Vol_[index_j]; } @@ -60,11 +61,11 @@ void ViscousForce, ViscosityType, KernelCorrectionTyp interaction(size_t index_i, Real dt) { Vecd force = Vecd::Zero(); - Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - Vecd &e_ij = inner_neighborhood.e_ij_[n]; + const Vecd &e_ij = inner_neighborhood.e_ij_[n]; Real r_ij = inner_neighborhood.r_ij_[n]; /** The following viscous force is given in Monaghan 2005 (Rep. Prog. Phys.), it seems that @@ -72,7 +73,7 @@ void ViscousForce, ViscosityType, KernelCorrectionTyp Real v_r_ij = (vel_[index_i] - vel_[index_j]).dot(e_ij); Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_j) * v_r_ij / (r_ij + 0.01 * smoothing_length_); - force += (kernel_correction_(index_i) + kernel_correction_(index_j)) * + force += e_ij.dot((kernel_correction_(index_i) + kernel_correction_(index_j)) * e_ij) * eta_ij * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * e_ij; } @@ -94,15 +95,16 @@ void ViscousForce, ViscosityType, KernelCorrectionType>:: { StdLargeVec &vel_ave_k = *(wall_vel_ave_[k]); StdLargeVec &wall_Vol_k = *(wall_Vol_[k]); - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { size_t index_j = contact_neighborhood.j_[n]; Real r_ij = contact_neighborhood.r_ij_[n]; + const Vecd &e_ij = contact_neighborhood.e_ij_[n]; Vecd vel_derivative = 2.0 * (vel_[index_i] - vel_ave_k[index_j]) / (r_ij + 0.01 * smoothing_length_); - force += 2.0 * kernel_correction_(index_i) * mu_(index_i, index_i) * + force += 2.0 * e_ij.dot(kernel_correction_(index_i) * e_ij) * mu_(index_i, index_i) * vel_derivative * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j]; } } @@ -125,18 +127,18 @@ void ViscousForce, ViscosityType, KernelCorre { StdLargeVec &vel_ave_k = *(wall_vel_ave_[k]); StdLargeVec &wall_Vol_k = *(wall_Vol_[k]); - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { size_t index_j = contact_neighborhood.j_[n]; - Vecd &e_ij = contact_neighborhood.e_ij_[n]; + const Vecd &e_ij = contact_neighborhood.e_ij_[n]; Real r_ij = contact_neighborhood.r_ij_[n]; Real v_r_ij = 2.0 * (vel_[index_i] - vel_ave_k[index_j]).dot(e_ij); Real eta_ij = Real(Dimensions + 2) * mu_(index_i, index_i) * v_r_ij / (r_ij + 0.01 * smoothing_length_); - force += 2.0 * kernel_correction_(index_i) * eta_ij * - contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij; + force += 2.0 * e_ij.dot(kernel_correction_(index_i) * e_ij) * mu_(index_i, index_i) * + eta_ij * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j] * e_ij; } } @@ -169,13 +171,15 @@ void ViscousForce, ViscosityType, KernelCorrectionType>:: StdLargeVec &vel_k = *(contact_vel_[k]); StdLargeVec &wall_Vol_k = *(wall_Vol_[k]); KernelCorrectionType &kernel_correction_k = contact_kernel_corrections_[k]; - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + const Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { size_t index_j = contact_neighborhood.j_[n]; + const Vecd &e_ij = contact_neighborhood.e_ij_[n]; + Vecd vel_derivative = (vel_[index_i] - vel_k[index_j]) / (contact_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); - force += (kernel_correction_(index_i) + kernel_correction_k(index_j)) * + force += e_ij.dot((kernel_correction_(index_i) + kernel_correction_k(index_j)) * e_ij) * contact_mu_k(index_i, index_j) * vel_derivative * contact_neighborhood.dW_ij_[n] * wall_Vol_k[index_j]; } diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml index 21f59ae38f..6ac4c2a0b0 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_0_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml index 0f3d27c3b9..ba8360a519 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_3_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml new file mode 100644 index 0000000000..834c8af597 --- /dev/null +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml deleted file mode 100644 index 19ffb9658b..0000000000 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_Run_6_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml index 36b6bf2534..8018a3a600 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat index 34de7daedb..de07de18db 100644 --- a/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat +++ b/tests/2d_examples/test_2d_poiseuille_flow/regression_test_tool/WaterBody_TotalKineticEnergy_runtimes.dat @@ -1,3 +1,3 @@ true -7 +6 4 \ No newline at end of file