Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Correct sign for efield...
  • Loading branch information
fweik committed Mar 25, 2015
1 parent 20263fe commit 126d7b5
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 33 deletions.
46 changes: 23 additions & 23 deletions src/core/electrokinetics_cuda.cu
Expand Up @@ -2066,60 +2066,60 @@ __global__ void ek_spread_particle_force( CUDA_particle_data * particle_data,
lowernode[1] = (lowernode[1] + ek_lbparameters_gpu->dim_y) % ek_lbparameters_gpu->dim_y;
lowernode[2] = (lowernode[2] + ek_lbparameters_gpu->dim_z) % ek_lbparameters_gpu->dim_z;

float force[3] = { 0., 0., 0. };
float efield[3] = { 0., 0., 0. };
#pragma unroll 3
for(unsigned int dim = 0; dim < 3; ++dim) {
// 0 0 0
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
lowernode[1],
lowernode[2]) + dim]
* particle_data[ index ].q * ( 1 - cellpos[0] ) * ( 1 - cellpos[1] ) * ( 1 - cellpos[2] );
* ( 1 - cellpos[0] ) * ( 1 - cellpos[1] ) * ( 1 - cellpos[2] );

// 0 0 1
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
lowernode[1],
(lowernode[2] + 1 ) % ek_lbparameters_gpu->dim_z ) + dim]
* particle_data[ index ].q * ( 1 - cellpos[0] ) * ( 1 - cellpos[1] ) * cellpos[2];
* ( 1 - cellpos[0] ) * ( 1 - cellpos[1] ) * cellpos[2];

// 0 1 0
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
(lowernode[1] + 1) % ek_lbparameters_gpu->dim_y,
lowernode[2] ) + dim]
* particle_data[ index ].q * ( 1 - cellpos[0] ) * cellpos[1] * ( 1 - cellpos[2] );
* ( 1 - cellpos[0] ) * cellpos[1] * ( 1 - cellpos[2] );

// 0 1 1
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( lowernode[0],
(lowernode[1] + 1) % ek_lbparameters_gpu->dim_y,
(lowernode[2] + 1 ) % ek_lbparameters_gpu->dim_z ) + dim]
* particle_data[ index ].q * ( 1 - cellpos[0] ) * cellpos[1] * cellpos[2];
* ( 1 - cellpos[0] ) * cellpos[1] * cellpos[2];

// 1 0 0
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
lowernode[1],
lowernode[2] ) + dim]
* particle_data[ index ].q * cellpos[0] * ( 1 - cellpos[1] ) * ( 1 - cellpos[2] );
* cellpos[0] * ( 1 - cellpos[1] ) * ( 1 - cellpos[2] );

// 1 0 1
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
lowernode[1],
(lowernode[2] + 1 ) % ek_lbparameters_gpu->dim_z ) + dim]
* particle_data[ index ].q * cellpos[0] * ( 1 - cellpos[1] ) * cellpos[2];
* cellpos[0] * ( 1 - cellpos[1] ) * cellpos[2];

// 1 1 0
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
(lowernode[1] + 1) % ek_lbparameters_gpu->dim_y,
lowernode[2] ) + dim]
* particle_data[ index ].q * cellpos[0] * cellpos[1] * ( 1 - cellpos[2] );
* cellpos[0] * cellpos[1] * ( 1 - cellpos[2] );

// 1 1 1
force[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
efield[dim] += ek_parameters_gpu.electric_field[3*rhoindex_cartesian2linear( (lowernode[0] + 1) % ek_lbparameters_gpu->dim_x,
(lowernode[1] + 1) % ek_lbparameters_gpu->dim_y,
(lowernode[2] + 1 ) % ek_lbparameters_gpu->dim_z ) + dim]
* particle_data[ index ].q * cellpos[0] * cellpos[1] * cellpos[2];
* cellpos[0] * cellpos[1] * cellpos[2];
}
particle_forces[index].f[0] += force[0];
particle_forces[index].f[1] += force[1];
particle_forces[index].f[2] += force[2];
particle_forces[index].f[0] += particle_data[ index ].q * efield[0];
particle_forces[index].f[1] += particle_data[ index ].q * efield[1];
particle_forces[index].f[2] += particle_data[ index ].q * efield[2];
}
}

Expand All @@ -2132,17 +2132,17 @@ __global__ void ek_calc_electric_field(const float *potential) {
rhoindex_linear2cartesian(index, coord);
const float agrid_inv = 1.0f / ek_parameters_gpu.agrid;

ek_parameters_gpu.electric_field[3*index + 0] = 0.5f * agrid_inv *
ek_parameters_gpu.electric_field[3*index + 0] = -0.5f * agrid_inv *
(
potential[rhoindex_cartesian2linear((coord[0] + 1) % ek_parameters_gpu.dim_x, coord[1], coord[2])]
- potential[rhoindex_cartesian2linear((coord[0] - 1 + ek_parameters_gpu.dim_x) % ek_parameters_gpu.dim_x, coord[1], coord[2])]
);
ek_parameters_gpu.electric_field[3*index + 1] = 0.5f * agrid_inv *
ek_parameters_gpu.electric_field[3*index + 1] = -0.5f * agrid_inv *
(
potential[rhoindex_cartesian2linear(coord[0], (coord[1] + 1) % ek_parameters_gpu.dim_y, coord[2])]
- potential[rhoindex_cartesian2linear(coord[0], (coord[1] - 1 + ek_parameters_gpu.dim_y) % ek_parameters_gpu.dim_y, coord[2])]
);
ek_parameters_gpu.electric_field[3*index + 2] = 0.5f * agrid_inv *
ek_parameters_gpu.electric_field[3*index + 2] = -0.5f * agrid_inv *
(
potential[rhoindex_cartesian2linear(coord[0], coord[1], (coord[2] + 1) % ek_parameters_gpu.dim_z)]
- potential[rhoindex_cartesian2linear(coord[0], coord[1], (coord[2] - 1 + ek_parameters_gpu.dim_z) % ek_parameters_gpu.dim_z)]
Expand Down
23 changes: 13 additions & 10 deletions testsuite/ek_electrostatics_coupling.tcl
Expand Up @@ -36,18 +36,17 @@ puts "###############################################################\n"

set box_x 6
set box_y 6
set width 50
set width 40

set padding 6
set padding 10
set box_z [expr $width+2*$padding]

setmd box_l $box_x $box_y $box_z

# Set the electrokinetic parameters

set agrid 1.0
set dt [expr 1.0/7.0]
set force 0.13
set dt [expr 1e-15]
set sigma -0.05
set viscosity_kinematic 2.3
set friction 4.3
Expand All @@ -72,19 +71,23 @@ set valency 1.0

# Set up the (LB) electrokinetics fluid

electrokinetics agrid $agrid lb_density $density_water viscosity $viscosity_kinematic friction $friction T $temperature bjerrum_length $bjerrum_length stencil linkcentered
electrokinetics agrid $agrid lb_density $density_water viscosity $viscosity_kinematic friction $friction T $temperature bjerrum_length $bjerrum_length stencil linkcentered electrostatics_coupling

electrokinetics 1 density 0.0 D 0.0 valency $valency 0 0
electrokinetics 1 density 0.0 D 1.0 valency $valency

# Set up the charged boundaries

electrokinetics boundary charge_density [expr $sigma/$agrid] wall normal 0 0 1 d $padding 0 0 direction outside
electrokinetics boundary charge_density [expr $sigma/$agrid] wall normal 0 0 -1 d -[expr $padding+$width] 0 0 direction outside
electrokinetics boundary charge_density [expr -$sigma/$agrid] wall normal 0 0 -1 d -[expr $padding+$width] 0 0 direction outside

# Integrate the system

for { set i 0 } { i < 100 } { incr i } {
part 0 pos [expr $padding + 0.01*$i * $width] 0 0 q 1.0
set chen [open "forces.dat" "w"]

for { set i 10 } { $i < 90 } { incr i } {
part 0 pos 3. 3. [expr $padding + 0.01*$i * $width] q 1.0
integrate 0
puts [part 0 pr pos f]
puts $chen [part 0 pr pos f]
}

close $chen

0 comments on commit 126d7b5

Please sign in to comment.