Permalink
Browse files

Merge pull request #976 from sschoell/ek_remove_total_momentum

Ek remove total momentum
  • Loading branch information...
2 parents b36e54d + 1c012ec commit 803902ef3865c82104f3522032768c8c340903f2 @rempferg rempferg committed on GitHub Jan 31, 2017
@@ -126,7 +126,11 @@ 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;
+ 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;
#endif
#ifdef ROTATION
@@ -237,6 +241,10 @@ 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,8 +88,12 @@ typedef struct {
float q;
#endif
+#ifdef MASS
+ float mass;
+#endif
+
unsigned int fixed;
-
+
#ifdef IMMERSED_BOUNDARY
bool isVirtual;
#endif
@@ -1,3 +1,38 @@
#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,6 +188,12 @@ 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,6 +3880,92 @@ 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 < 2)
+ if(argc < 1)
{
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,6 +108,7 @@ 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"))
@@ -615,6 +616,10 @@ 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 803902e

Please sign in to comment.