Revert "Ek remove total momentum"

src/core/cuda_interface.cpp

@@ -126,11 +126,7 @@ void cuda_mpi_get_particles(CUDA_particle_data *particle_data_host)
 #endif
 
 #ifdef ELECTROSTATICS
-                particle_data_host[i+g].q = (float)part[i].p.q;
-#endif
-
-#ifdef MASS
-                particle_data_host[i+g].mass = (float)part[i].p.mass;
+		particle_data_host[i+g].q = (float)part[i].p.q;
 #endif
 
 #ifdef ROTATION
@@ -241,10 +237,6 @@ static void cuda_mpi_get_particles_slave(){
             particle_data_host_sl[i+g].q = (float)part[i].p.q;
   #endif
 
-  #ifdef MASS
-      particle_data_host_sl[i+g].mass = (float)part[i].p.mass;
-  #endif
-
 #ifdef ROTATION
           particle_data_host_sl[i+g].quatu[0] = (float)part[i].r.quatu[0];
           particle_data_host_sl[i+g].quatu[1] = (float)part[i].r.quatu[1];

src/core/cuda_interface.hpp

@@ -66,7 +66,7 @@ typedef struct {
 #ifdef ENGINE
   CUDA_ParticleParametersSwimming swim;
 #endif
-
+  
   /** particle position given from md part*/
   float p[3];
   /** particle momentum struct velocity p.m->v*/
@@ -88,12 +88,8 @@ typedef struct {
   float q;
 #endif
 
-#ifdef MASS
-  float mass;
-#endif
-
   unsigned int fixed;
-
+  
 #ifdef IMMERSED_BOUNDARY
   bool isVirtual;
 #endif

src/core/electrokinetics.cpp

@@ -1,38 +1,3 @@
 #include "config.hpp"
 #include "electrokinetics.hpp"
 
-#include "particle_data.hpp"
-#include "statistics.hpp"
-
-#include "cuda_interface.hpp"
-#include "lbgpu.hpp"
-
-extern CUDA_particle_data *particle_data_host;
-extern LB_parameters_gpu lbpar_gpu;
-
-void ek_particles_add_momentum ( ekfloat velocity[3] ) {
-  for(int i=0; i<lbpar_gpu.number_of_particles; i++){
-    double mass;
-#ifdef MASS
-    mass = particle_data_host[i].mass;
-#else
-    mass = 1.;
-#endif
-
-    double new_velocity[3] = {
-      particle_data_host[i].v[0] + velocity[0],
-      particle_data_host[i].v[1] + velocity[1],
-      particle_data_host[i].v[2] + velocity[2]
-    };
-    set_particle_v( i, new_velocity );
-  }
-}
-
-void ek_get_total_momentum ( ekfloat momentum[3] ) {
-  std::vector<double> total_momentum;
-  total_momentum = calc_linear_momentum(1,1);
-  momentum[0] = total_momentum[0];
-  momentum[1] = total_momentum[1];
-  momentum[2] = total_momentum[2];
-  return;
-}

src/core/electrokinetics.hpp

@@ -188,12 +188,6 @@ int ek_node_print_flux(int species, int x, int y, int z, double* flux);
 int ek_node_set_density(int species, int x, int y, int z, double density);
 ekfloat ek_calculate_net_charge(); 
 int ek_neutralize_system(int species); 
-int ek_remove_total_momentum();
-void ek_fluid_add_momentum(ekfloat momentum[3]);
-void ek_get_total_momentum(ekfloat momentum[3]);
-void ek_particles_add_momentum(ekfloat velocity[3]);
-
-
 int ek_save_checkpoint(char* filename);
 int ek_load_checkpoint(char* filename);
 

src/core/electrokinetics_cuda.cu

@@ -3880,92 +3880,6 @@ int ek_neutralize_system(int species) {
   return 0;
 }
 
-struct ek_mass_of_particle
-{
-  __device__ float operator()(CUDA_particle_data particle)
-  {
-#ifdef MASS
-    return particle.mass;
-#else
-    return 1.;
-#endif
-  };
-};
-
-int ek_remove_total_momentum() {
-  ekfloat total_momentum[3];
-  ekfloat velocity[3];
-  ekfloat fluid_mass;
-  ekfloat particles_mass;
-
-  // calculate momentum of fluid and particles
-  ek_get_total_momentum(total_momentum);
-
-  particle_data_gpu = gpu_get_particle_pointer();
-  thrust::device_ptr<CUDA_particle_data> ptr(particle_data_gpu);
-  particles_mass = thrust::transform_reduce(
-      ptr,
-      ptr + lbpar_gpu.number_of_particles,
-      ek_mass_of_particle(),
-      0.0f,
-      thrust::plus<ekfloat>());
-
-  lb_calc_fluid_mass_GPU( &fluid_mass );
-
-  velocity[0] = -total_momentum[0]/(fluid_mass + particles_mass);
-  velocity[1] = -total_momentum[1]/(fluid_mass + particles_mass);
-  velocity[2] = -total_momentum[2]/(fluid_mass + particles_mass);
-
-  ekfloat momentum_fluid[3] = {
-    velocity[0]*fluid_mass,
-    velocity[1]*fluid_mass,
-    velocity[2]*fluid_mass
-  };
-
-  ek_particles_add_momentum( velocity );
-  ek_fluid_add_momentum( momentum_fluid );
-
-  return 0;
-}
-
-__global__ void ek_fluid_add_momentum_kernel(
-                                            ekfloat momentum[3],
-                                            LB_nodes_gpu n_a,
-                                            LB_parameters_gpu *ek_lbparameters,
-                                            LB_node_force_gpu node_f,
-                                            LB_rho_v_gpu *d_v
-                                          ){
-  unsigned int index = ek_getThreadIndex();
-  unsigned int number_of_nodes = ek_parameters_gpu.number_of_nodes;
-#ifdef EK_BOUNDARIES
-  number_of_nodes -= ek_parameters_gpu.number_of_boundary_nodes;
-#endif
-  if( index < ek_parameters_gpu.number_of_nodes ) {
-    if( n_a.boundary[index] == 0 ){
-      float force_factor=powf(ek_lbparameters->agrid,2)*ek_lbparameters->tau*ek_lbparameters->tau;
-      for(int ii=0 ; ii < LB_COMPONENTS ; ii++ ) { 
-        // add force density onto each node (momentum / time_step / Volume)
-        node_f.force[(0+ii*3)*ek_parameters_gpu.number_of_nodes + index] += momentum[0] / ek_lbparameters->tau / (number_of_nodes * powf(ek_lbparameters->agrid,3)) * force_factor;
-        node_f.force[(1+ii*3)*ek_parameters_gpu.number_of_nodes + index] += momentum[1] / ek_lbparameters->tau / (number_of_nodes * powf(ek_lbparameters->agrid,3)) * force_factor;
-        node_f.force[(2+ii*3)*ek_parameters_gpu.number_of_nodes + index] += momentum[2] / ek_lbparameters->tau / (number_of_nodes * powf(ek_lbparameters->agrid,3)) * force_factor;
-      }
-    }
-  }
-}
-
-void ek_fluid_add_momentum( ekfloat momentum_host[3] ){
-  ekfloat* momentum_device;
-  cuda_safe_mem(cudaMalloc((void**)&momentum_device,3*sizeof(ekfloat)));
-  cuda_safe_mem(cudaMemcpy(momentum_device, momentum_host, 3*sizeof(ekfloat), cudaMemcpyHostToDevice));
-
-  int threads_per_block = 64;
-  int blocks_per_grid_y = 4;
-  int blocks_per_grid_x = (lbpar_gpu.number_of_nodes + threads_per_block * blocks_per_grid_y - 1)/(threads_per_block * blocks_per_grid_y);
-  dim3 dim_grid = make_uint3(blocks_per_grid_x, blocks_per_grid_y, 1);
-
-  KERNELCALL( ek_fluid_add_momentum_kernel, dim_grid, threads_per_block, (momentum_device, *current_nodes, ek_lbparameters_gpu, node_f, ek_lb_device_values));
-}
-
 int ek_save_checkpoint(char* filename) {
   std::string fname(filename);
   std::ofstream fout((const char *) (fname + ".ek").c_str(), std::ofstream::binary);

src/tcl/electrokinetics_tcl.cpp

@@ -61,7 +61,7 @@ int tclcommand_electrokinetics(ClientData data, Tcl_Interp *interp, int argc, ch
          mass_product1;
 #endif
 
-  if(argc < 1) 
+  if(argc < 2) 
   {
     Tcl_AppendResult(interp, "Usage of \"electrokinetics\":\n\n", (char *)NULL);
     Tcl_AppendResult(interp, "electrokinetics [agrid #float] [lb_density #float] [viscosity #float] [friction #float]\n", (char *)NULL);
@@ -108,7 +108,6 @@ int tclcommand_electrokinetics(ClientData data, Tcl_Interp *interp, int argc, ch
     Tcl_AppendResult(interp, "electrokinetics #int node #int #int #int print <density|flux>\n", (char *)NULL);
     Tcl_AppendResult(interp, "electrokinetics pdb-parse #string #string #double\n", (char *)NULL);
     Tcl_AppendResult(interp, "electrokinetics checkpoint <save|load> #string\n", (char *)NULL);
-    Tcl_AppendResult(interp, "electrokinetics remove_total_momentum\n", (char *)NULL);
     return TCL_ERROR;
   }
   else if(ARG0_IS_S("checkpoint"))
@@ -616,10 +615,6 @@ int tclcommand_electrokinetics(ClientData data, Tcl_Interp *interp, int argc, ch
       }
     }
   }
-  else if(ARG0_IS_S("remove_total_momentum"))
-  {
-    err = ek_remove_total_momentum();
-  }
   else 
   {
     Tcl_ResetResult(interp);