Correct second order term for forces in LB (#3885)

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.

Description of changes:
- Replace `gamma_bulk` by `gamma_shear` in the corresponding terms in `lb_apply_forces`

src/core/grid_based_algorithms/lb.cpp

@@ -894,13 +894,13 @@ std::array<T, 19> 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] =

src/core/grid_based_algorithms/lbgpu_cuda.cu

@@ -1012,21 +1012,21 @@ __device__ void apply_forces(unsigned int index, Utils::Array<float, 19> &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] +