diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 1bbf4b3d56c..8cfedeaf2e7 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -1283,7 +1283,7 @@ Initialization system.time_step = 0.0 system.cell_system.skin = 0.4 ek = espressomd.electrokinetics.Electrokinetics(agrid=1.0, lb_density=1.0, - viscosity=1.0, friction=1.0, T=1.0, prefactor=1.0, + viscosity=1.0, ext_force_density = [1,0,0], friction=1.0, T=1.0, prefactor=1.0, stencil='linkcentered', advection=True, fluid_coupling='friction') system.actors.add(ek) @@ -1303,7 +1303,8 @@ that your computer contains a CUDA capable GPU which is sufficiently modern. To set up a proper LB fluid using this command one has to specify at -least the following options: ``agrid``, ``lb_density``, ``viscosity``, ``friction``, ``T``, and ``prefactor``. The other options can be +least the following options: ``agrid``, ``lb_density``, ``viscosity``, +``friction``, ``T``, and ``prefactor``. The other options can be used to modify the behavior of the LB fluid. Note that the command does not allow the user to set the time step parameter as is the case for the lattice-Boltzmann command, this parameter is instead taken directly from the value set for @@ -1382,7 +1383,7 @@ Boundaries ^^^^^^^^^^ :: - ek_boundary = espressomd.electrokinetics.EKBoundary(charge_density=1.0, shape=my_shape) + ek_boundary = espressomd.ekboundaries.EKBoundary(charge_density=1.0, shape=my_shape) system.ekboundaries.add(ek_boundary) .. note:: Feature ``EK_BOUNDARIES`` required diff --git a/src/core/grid_based_algorithms/electrokinetics.hpp b/src/core/grid_based_algorithms/electrokinetics.hpp index b0d2413ee06..0e5316211aa 100644 --- a/src/core/grid_based_algorithms/electrokinetics.hpp +++ b/src/core/grid_based_algorithms/electrokinetics.hpp @@ -48,7 +48,7 @@ typedef struct { float friction; float T; float prefactor; - float lb_force_density[3]; + float lb_ext_force_density[3]; unsigned int number_of_species; int reaction_species[3]; float rho_reactant_reservoir; @@ -150,6 +150,9 @@ int ek_init(); int ek_set_agrid(float agrid); int ek_set_lb_density(float lb_density); int ek_set_viscosity(float viscosity); +int ek_set_lb_ext_force_density(float lb_ext_force_dens_x, + float lb_ext_force_dens_y, + float lb_ext_force_dens_z); int ek_set_friction(float friction); int ek_set_T(float T); int ek_set_prefactor(float prefactor); diff --git a/src/core/grid_based_algorithms/electrokinetics_cuda.cu b/src/core/grid_based_algorithms/electrokinetics_cuda.cu index 069ea3579c5..1eeba89b6b5 100644 --- a/src/core/grid_based_algorithms/electrokinetics_cuda.cu +++ b/src/core/grid_based_algorithms/electrokinetics_cuda.cu @@ -103,7 +103,7 @@ EK_parameters ek_parameters = { -1.0, // prefactor -1.0, - // lb_force_density + // lb_ext_force_density {0.0, 0.0, 0.0}, // number_of_species 0, @@ -2241,49 +2241,33 @@ int ek_init() { lbpar_gpu.rho = (ek_parameters.lb_density < 0.0) ? 1.0f - : ek_parameters.lb_density * Utils::int_pow<3>(ek_parameters.agrid); + : ek_parameters.lb_density * Utils::int_pow<3>(lbpar_gpu.agrid); lbpar_gpu.is_TRT = true; lb_reinit_parameters_gpu(); lbpar_gpu.viscosity = ek_parameters.viscosity * lbpar_gpu.time_step / - Utils::sqr(ek_parameters.agrid); + Utils::sqr(lbpar_gpu.agrid); lbpar_gpu.bulk_viscosity = ek_parameters.bulk_viscosity * lbpar_gpu.time_step / - Utils::sqr(ek_parameters.agrid); - lb_reinit_parameters_gpu(); + Utils::sqr(lbpar_gpu.agrid); - lb_init_gpu(); + lbpar_gpu.external_force_density = + ek_parameters.lb_ext_force_density[0] != 0 || + ek_parameters.lb_ext_force_density[1] != 0 || + ek_parameters.lb_ext_force_density[2] != 0; + lbpar_gpu.ext_force_density = + Utils::Vector3f(ek_parameters.lb_ext_force_density) * + Utils::sqr(lbpar_gpu.agrid * lbpar_gpu.time_step); - if (ek_parameters.lb_force_density[0] != 0 || - ek_parameters.lb_force_density[1] != 0 || - ek_parameters.lb_force_density[2] != 0) { - lbpar_gpu.external_force_density = 1; - lbpar_gpu.ext_force_density[0] = - ek_parameters.lb_force_density[0] * ek_parameters.agrid * - ek_parameters.agrid * ek_parameters.time_step * - ek_parameters.time_step; - lbpar_gpu.ext_force_density[1] = - ek_parameters.lb_force_density[1] * ek_parameters.agrid * - ek_parameters.agrid * ek_parameters.time_step * - ek_parameters.time_step; - lbpar_gpu.ext_force_density[2] = - ek_parameters.lb_force_density[2] * ek_parameters.agrid * - ek_parameters.agrid * ek_parameters.time_step * - ek_parameters.time_step; - lb_reinit_extern_nodeforce_GPU(&lbpar_gpu); - } else { - lbpar_gpu.external_force_density = 0; - lbpar_gpu.ext_force_density[0] = 0; - lbpar_gpu.ext_force_density[1] = 0; - lbpar_gpu.ext_force_density[2] = 0; - } + lb_reinit_parameters_gpu(); + lb_init_gpu(); + ek_parameters.time_step = lbpar_gpu.time_step; ek_parameters.dim_x = lbpar_gpu.dim[0]; ek_parameters.dim_x_padded = (ek_parameters.dim_x / 2 + 1) * 2; ek_parameters.dim_y = lbpar_gpu.dim[1]; ek_parameters.dim_z = lbpar_gpu.dim[2]; - ek_parameters.time_step = lbpar_gpu.time_step; ek_parameters.number_of_nodes = ek_parameters.dim_x * ek_parameters.dim_y * ek_parameters.dim_z; @@ -3535,9 +3519,10 @@ void ek_print_parameters() { printf(" float friction = %f;\n", ek_parameters.friction); printf(" float T = %f;\n", ek_parameters.T); printf(" float prefactor = %f;\n", ek_parameters.prefactor); - printf(" float lb_force_density[] = {%f, %f, %f};\n", - ek_parameters.lb_force_density[0], ek_parameters.lb_force_density[1], - ek_parameters.lb_force_density[2]); + printf(" float lb_ext_force_density[] = {%f, %f, %f};\n", + ek_parameters.lb_ext_force_density[0], + ek_parameters.lb_ext_force_density[1], + ek_parameters.lb_ext_force_density[2]); printf(" unsigned int number_of_species = %d;\n", ek_parameters.number_of_species); printf(" int reaction_species[] = {%d, %d, %d};\n", @@ -3697,6 +3682,15 @@ int ek_set_viscosity(float viscosity) { return 0; } +int ek_set_lb_ext_force_density(float lb_ext_force_dens_x, + float lb_ext_force_dens_y, + float lb_ext_force_dens_z) { + ek_parameters.lb_ext_force_density[0] = lb_ext_force_dens_x; + ek_parameters.lb_ext_force_density[1] = lb_ext_force_dens_y; + ek_parameters.lb_ext_force_density[2] = lb_ext_force_dens_z; + return 0; +} + int ek_set_friction(float friction) { ek_parameters.friction = friction; diff --git a/src/python/espressomd/electrokinetics.pxd b/src/python/espressomd/electrokinetics.pxd index 4c6fa8bad7e..f42961044ae 100644 --- a/src/python/espressomd/electrokinetics.pxd +++ b/src/python/espressomd/electrokinetics.pxd @@ -39,7 +39,7 @@ IF ELECTROKINETICS and CUDA: float friction float T float prefactor - float lb_force_density[3] + float lb_ext_force_density[3] unsigned int number_of_species int reaction_species[3] float rho_reactant_reservoir @@ -88,7 +88,7 @@ IF ELECTROKINETICS and CUDA: float friction float T float prefactor - float lb_force_density[3] + float lb_ext_force_density[3] unsigned int number_of_species int reaction_species[3] float rho_reactant_reservoir @@ -140,6 +140,7 @@ IF ELECTROKINETICS and CUDA: int ek_set_lb_density(float lb_density) int ek_set_viscosity(float viscosity) int ek_set_friction(float friction) + int ek_set_lb_ext_force_density(float lb_ext_force_dens_x, float lb_ext_force_dens_y, float lb_ext_force_dens_z) int ek_set_T(float T) int ek_set_prefactor(float prefactor) int ek_set_bulk_viscosity(float bulk_viscosity) diff --git a/src/python/espressomd/electrokinetics.pyx b/src/python/espressomd/electrokinetics.pyx index 73edb490835..63c5ec5f910 100644 --- a/src/python/espressomd/electrokinetics.pyx +++ b/src/python/espressomd/electrokinetics.pyx @@ -66,7 +66,7 @@ IF ELECTROKINETICS: """ return ["agrid", "lb_density", "viscosity", "friction", - "bulk_viscosity", "gamma_even", "gamma_odd", "T", + "bulk_viscosity", "gamma_even", "gamma_odd", "T", "ext_force_density", "prefactor", "stencil", "advection", "fluid_coupling", "fluctuations", "fluctuation_amplitude", "es_coupling", "species"] @@ -90,6 +90,7 @@ IF ELECTROKINETICS: "bulk_viscosity": -1, "gamma_odd": 0.0, "gamma_even": 0.0, + "ext_force_density": [0., 0., 0.], "friction": 0.0, "T": -1, "prefactor": -1, @@ -120,6 +121,7 @@ IF ELECTROKINETICS: "bulk_viscosity": ek_parameters.bulk_viscosity, "gamma_odd": ek_parameters.gamma_odd, "gamma_even": ek_parameters.gamma_even, + "ext_force_density": ek_parameters.lb_ext_force_density, "friction": ek_parameters.friction, "T": ek_parameters.T, "prefactor": ek_parameters.prefactor, @@ -146,6 +148,9 @@ IF ELECTROKINETICS: ek_set_lb_density(self._params["lb_density"]) ek_set_viscosity(self._params["viscosity"]) ek_set_friction(self._params["friction"]) + ek_set_lb_ext_force_density(self._params["ext_force_density"][0], + self._params["ext_force_density"][1], + self._params["ext_force_density"][2]) ek_set_T(self._params["T"]) ek_set_prefactor(self._params["prefactor"]) ek_set_bulk_viscosity(self._params["bulk_viscosity"]) diff --git a/testsuite/python/lb_poiseuille.py b/testsuite/python/lb_poiseuille.py index 8bd38bd9511..39e514375d2 100644 --- a/testsuite/python/lb_poiseuille.py +++ b/testsuite/python/lb_poiseuille.py @@ -20,6 +20,7 @@ import espressomd.lb import espressomd.lbboundaries +import espressomd.electrokinetics import espressomd.shapes @@ -41,6 +42,14 @@ 'tau': TIME_STEP, 'ext_force_density': [0.0, 0.0, EXT_FORCE]} +EK_PARAMS = {'agrid': AGRID, + 'lb_density': DENS, + 'viscosity': VISC, + 'ext_force_density': [0.0, 0.0, EXT_FORCE], + 'friction': 0., + 'T': 1, + 'prefactor': 0.} + def poiseuille_flow(z, H, ext_force_density, dyn_visc): """ @@ -86,9 +95,7 @@ def prepare(self): self.system.lbboundaries.add(wall1) self.system.lbboundaries.add(wall2) - mid_indices = [int((self.system.box_l[0] / AGRID) / 2), - int((self.system.box_l[1] / AGRID) / 2), - int((self.system.box_l[2] / AGRID) / 2)] + mid_indices = (self.system.box_l / AGRID / 2).astype(int) diff = float("inf") old_val = self.lbf[mid_indices].velocity[2] while diff > 0.005: @@ -114,9 +121,7 @@ def test_profile(self): velocities[x, 0] = (x + 0.5) * AGRID v_measured = velocities[1:-1, 1] - v_expected = poiseuille_flow(velocities[1:-1, - 0] - 0.5 * self.system.box_l[ - 0], + v_expected = poiseuille_flow(velocities[1:-1, 0] - 0.5 * self.system.box_l[0], self.system.box_l[0] - 2.0 * AGRID, EXT_FORCE, VISC * DENS) @@ -143,6 +148,20 @@ def setUp(self): self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS) +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures( + ['LB_BOUNDARIES_GPU', "ELECTROKINETICS", "EXTERNAL_FORCES"]) +class LBEkinPoiseuille(ut.TestCase, LBPoiseuilleCommon): + + """Test the LB part of electrokinetics. """ + + def setUp(self): + self.lbf = espressomd.electrokinetics.Electrokinetics(**EK_PARAMS) + species = espressomd.electrokinetics.Species( + density=0., D=1., valency=0.) + self.lbf.add_species(species) + + @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures(['LB_BOUNDARIES_GPU', 'EXTERNAL_FORCES']) class LBGPUPoiseuilleInterpolation(ut.TestCase, LBPoiseuilleCommon):