Permalink
Browse files

Merge pull request #1006 from espressomd/revert-976-ek_remove_total_m…

…omentum

Revert "Ek remove total momentum"
  • Loading branch information...
2 parents 803902e + 29fdfe7 commit 4440cfb1966178ac8ab08622ed5f6f6b8bceb7d3 @rempferg rempferg committed on GitHub Jan 31, 2017
@@ -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];
@@ -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
@@ -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;
-}
@@ -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);
@@ -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);
@@ -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);

0 comments on commit 4440cfb

Please sign in to comment.