From 998739777993ad7329835eaafcfb45c7d121c884 Mon Sep 17 00:00:00 2001 From: Ulf Schiller Date: Fri, 4 Sep 2020 17:39:46 -0400 Subject: [PATCH] Correct second order term for forces in LB The prefactor for the traceless part of the second order term is (1+gamma_shear)/2. Terms with this prefactor must cancel in the trace. Cf. Eq. (4.61) and (4.67) in my thesis. --- src/core/grid_based_algorithms/lb.cpp | 6 +++--- src/core/grid_based_algorithms/lbgpu_cuda.cu | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/grid_based_algorithms/lb.cpp b/src/core/grid_based_algorithms/lb.cpp index f465eebb450..1827d88553e 100644 --- a/src/core/grid_based_algorithms/lb.cpp +++ b/src/core/grid_based_algorithms/lb.cpp @@ -894,13 +894,13 @@ std::array lb_apply_forces(Lattice::index_t index, density; double C[6]; - C[0] = (1. + lb_parameters.gamma_bulk) * u[0] * f[0] + + C[0] = (1. + lb_parameters.gamma_shear) * u[0] * f[0] + 1. / 3. * (lb_parameters.gamma_bulk - lb_parameters.gamma_shear) * (u * f); - C[2] = (1. + lb_parameters.gamma_bulk) * u[1] * f[1] + + C[2] = (1. + lb_parameters.gamma_shear) * u[1] * f[1] + 1. / 3. * (lb_parameters.gamma_bulk - lb_parameters.gamma_shear) * (u * f); - C[5] = (1. + lb_parameters.gamma_bulk) * u[2] * f[2] + + C[5] = (1. + lb_parameters.gamma_shear) * u[2] * f[2] + 1. / 3. * (lb_parameters.gamma_bulk - lb_parameters.gamma_shear) * (u * f); C[1] = diff --git a/src/core/grid_based_algorithms/lbgpu_cuda.cu b/src/core/grid_based_algorithms/lbgpu_cuda.cu index c133bf97b3b..7f13c804148 100644 --- a/src/core/grid_based_algorithms/lbgpu_cuda.cu +++ b/src/core/grid_based_algorithms/lbgpu_cuda.cu @@ -1012,21 +1012,21 @@ __device__ void apply_forces(unsigned int index, Utils::Array &mode, u[1] = d_v[index].v[1]; u[2] = d_v[index].v[2]; - C[0] += (1.0f + para->gamma_bulk) * u[0] * + C[0] += (1.0f + para->gamma_shear) * u[0] * node_f.force_density[0 * para->number_of_nodes + index] + 1.0f / 3.0f * (para->gamma_bulk - para->gamma_shear) * (u[0] * node_f.force_density[0 * para->number_of_nodes + index] + u[1] * node_f.force_density[1 * para->number_of_nodes + index] + u[2] * node_f.force_density[2 * para->number_of_nodes + index]); - C[2] += (1.0f + para->gamma_bulk) * u[1] * + C[2] += (1.0f + para->gamma_shear) * u[1] * node_f.force_density[1 * para->number_of_nodes + index] + 1.0f / 3.0f * (para->gamma_bulk - para->gamma_shear) * (u[0] * node_f.force_density[0 * para->number_of_nodes + index] + u[1] * node_f.force_density[1 * para->number_of_nodes + index] + u[2] * node_f.force_density[2 * para->number_of_nodes + index]); - C[5] += (1.0f + para->gamma_bulk) * u[2] * + C[5] += (1.0f + para->gamma_shear) * u[2] * node_f.force_density[2 * para->number_of_nodes + index] + 1.0f / 3.0f * (para->gamma_bulk - para->gamma_shear) * (u[0] * node_f.force_density[0 * para->number_of_nodes + index] +