diff --git a/src/core/grid_based_algorithms/electrokinetics.hpp b/src/core/grid_based_algorithms/electrokinetics.hpp index 2dfea2379f8..c7799cbfcd9 100644 --- a/src/core/grid_based_algorithms/electrokinetics.hpp +++ b/src/core/grid_based_algorithms/electrokinetics.hpp @@ -147,31 +147,31 @@ int ek_print_vtk_lbforce_density(char *filename); int ek_lb_print_vtk_density(char *filename); int ek_lb_print_vtk_velocity(char *filename); 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); -int ek_set_electrostatics_coupling(bool electrostatics_coupling); +void ek_set_agrid(float agrid); +void ek_set_lb_density(float lb_density); +void ek_set_viscosity(float viscosity); +void ek_set_lb_ext_force_density(float lb_ext_force_dens_x, + float lb_ext_force_dens_y, + float lb_ext_force_dens_z); +void ek_set_friction(float friction); +void ek_set_T(float T); +void ek_set_prefactor(float prefactor); +void ek_set_electrostatics_coupling(bool electrostatics_coupling); void ek_calculate_electrostatic_coupling(); -int ek_set_bulk_viscosity(float bulk_viscosity); -int ek_set_gamma_odd(float gamma_odd); -int ek_set_gamma_even(float gamma_even); -int ek_set_density(int species, float density); -int ek_set_D(int species, float D); -int ek_set_valency(int species, float valency); -int ek_set_ext_force_density(int species, float ext_force_density_x, - float ext_force_density_y, - float ext_force_density_z); -int ek_set_stencil(int stencil); -int ek_set_advection(bool advection); -int ek_set_fluidcoupling(bool ideal_contribution); -int ek_set_fluctuations(bool fluctuations); -int ek_set_fluctuation_amplitude(float fluctuation_amplitude); +void ek_set_bulk_viscosity(float bulk_viscosity); +void ek_set_gamma_odd(float gamma_odd); +void ek_set_gamma_even(float gamma_even); +void ek_set_density(int species, float density); +void ek_set_D(int species, float D); +void ek_set_valency(int species, float valency); +void ek_set_ext_force_density(int species, float ext_force_density_x, + float ext_force_density_y, + float ext_force_density_z); +void ek_set_stencil(int stencil); +void ek_set_advection(bool advection); +void ek_set_fluidcoupling(bool ideal_contribution); +void ek_set_fluctuations(bool fluctuations); +void ek_set_fluctuation_amplitude(float fluctuation_amplitude); void ek_set_rng_state(uint64_t counter); int ek_node_get_density(int species, int x, int y, int z, double *density); int ek_node_get_flux(int species, int x, int y, int z, double *flux); @@ -179,8 +179,6 @@ int ek_node_get_potential(int x, int y, int z, double *potential); int ek_node_set_density(int species, int x, int y, int z, double density); float ek_calculate_net_charge(); int ek_neutralize_system(int species); -int ek_save_checkpoint(char *filename, char *lb_filename); -int ek_load_checkpoint(char *filename); #ifdef EK_BOUNDARIES void ek_gather_wallcharge_species_density(float *wallcharge_species_density, diff --git a/src/core/grid_based_algorithms/electrokinetics_cuda.cu b/src/core/grid_based_algorithms/electrokinetics_cuda.cu index 4a93e624b60..b7ad6e3ff1c 100644 --- a/src/core/grid_based_algorithms/electrokinetics_cuda.cu +++ b/src/core/grid_based_algorithms/electrokinetics_cuda.cu @@ -50,6 +50,7 @@ #include #include #include +#include #include #include @@ -3640,149 +3641,130 @@ void ek_print_lbpar() { printf("}\n"); } -int ek_set_agrid(float agrid) { +inline void ek_setter_throw_if_initialized() { + if (ek_initialized) + throw std::runtime_error( + "Electrokinetics parameters cannot be set after initialisation"); +}; +void ek_set_agrid(float agrid) { + ek_setter_throw_if_initialized(); ek_parameters.agrid = agrid; - return 0; } -int ek_set_lb_density(float lb_density) { - +void ek_set_lb_density(float lb_density) { + ek_setter_throw_if_initialized(); ek_parameters.lb_density = lb_density; - return 0; } -int ek_set_prefactor(float prefactor) { - +void ek_set_prefactor(float prefactor) { + ek_setter_throw_if_initialized(); ek_parameters.prefactor = prefactor; - return 0; } -int ek_set_electrostatics_coupling(bool electrostatics_coupling) { +void ek_set_electrostatics_coupling(bool electrostatics_coupling) { + ek_setter_throw_if_initialized(); ek_parameters.es_coupling = electrostatics_coupling; - return 0; } -int ek_set_viscosity(float viscosity) { - +void ek_set_viscosity(float viscosity) { + ek_setter_throw_if_initialized(); ek_parameters.viscosity = 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) { +void 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_setter_throw_if_initialized(); 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) { - +void ek_set_friction(float friction) { + ek_setter_throw_if_initialized(); ek_parameters.friction = friction; - return 0; } -int ek_set_bulk_viscosity(float bulk_viscosity) { - +void ek_set_bulk_viscosity(float bulk_viscosity) { + ek_setter_throw_if_initialized(); ek_parameters.bulk_viscosity = bulk_viscosity; - return 0; } -int ek_set_gamma_odd(float gamma_odd) { - +void ek_set_gamma_odd(float gamma_odd) { + ek_setter_throw_if_initialized(); ek_parameters.gamma_odd = gamma_odd; - return 0; } -int ek_set_gamma_even(float gamma_even) { +void ek_set_gamma_even(float gamma_even) { + ek_setter_throw_if_initialized(); ek_parameters.gamma_even = gamma_even; - return 0; } -int ek_set_stencil(int stencil) { +void ek_set_stencil(int stencil) { + ek_setter_throw_if_initialized(); if (!ek_parameters.fluidcoupling_ideal_contribution) - return 1; // combination not implemented - + throw std::runtime_error( + "Combination of stencil and fluid coupling not implmented."); ek_parameters.stencil = stencil; - return 0; } -int ek_set_advection(bool advection) { +void ek_set_advection(bool advection) { + ek_setter_throw_if_initialized(); ek_parameters.advection = advection; - return 0; } -int ek_set_fluctuations(bool fluctuations) { +void ek_set_fluctuations(bool fluctuations) { + ek_setter_throw_if_initialized(); ek_parameters.fluctuations = fluctuations; - return 0; } -int ek_set_fluctuation_amplitude(float fluctuation_amplitude) { +void ek_set_fluctuation_amplitude(float fluctuation_amplitude) { + ek_setter_throw_if_initialized(); ek_parameters.fluctuation_amplitude = fluctuation_amplitude; - return 0; } -int ek_set_fluidcoupling(bool ideal_contribution) { +void ek_set_fluidcoupling(bool ideal_contribution) { + ek_setter_throw_if_initialized(); if (ek_parameters.stencil != 0) - return 1; // combination not implemented - + throw std::runtime_error( + "Combination of stencil and fluid coupling not implemented."); ek_parameters.fluidcoupling_ideal_contribution = ideal_contribution; - return 0; } -int ek_set_density(int species, float density) { +void ek_set_T(float T) { + ek_setter_throw_if_initialized(); + ek_parameters.T = T; +} +void ek_set_density(int species, float density) { ek_init_species(species); - ek_parameters.density[ek_parameters.species_index[species]] = density; - - return 0; } -int ek_set_D(int species, float D) { - +void ek_set_D(int species, float D) { ek_init_species(species); - ek_parameters.D[ek_parameters.species_index[species]] = D; ek_parameters.d[ek_parameters.species_index[species]] = D / (1.0f + 2.0f * sqrt(2.0f)); - - return 0; } -int ek_set_T(float T) { - - ek_parameters.T = T; - - return 0; -} - -int ek_set_valency(int species, float valency) { - +void ek_set_valency(int species, float valency) { ek_init_species(species); - ek_parameters.valency[ek_parameters.species_index[species]] = valency; - - return 0; } -int ek_set_ext_force_density(int species, float ext_force_density_x, - float ext_force_density_y, - float ext_force_density_z) { - +void ek_set_ext_force_density(int species, float ext_force_density_x, + float ext_force_density_y, + float ext_force_density_z) { ek_init_species(species); - ek_parameters.ext_force_density[0][ek_parameters.species_index[species]] = ext_force_density_x; ek_parameters.ext_force_density[1][ek_parameters.species_index[species]] = ext_force_density_y; ek_parameters.ext_force_density[2][ek_parameters.species_index[species]] = ext_force_density_z; - - return 0; } struct ek_charge_of_particle { @@ -3866,57 +3848,6 @@ int ek_neutralize_system(int species) { return 0; } -int ek_save_checkpoint(char *filename, char *lb_filename) { - std::ofstream fout(filename, std::ofstream::binary); - std::vector densities(ek_parameters.number_of_nodes); - auto const nchars = - static_cast(densities.size() * sizeof(float)); - - for (unsigned i = 0; i < ek_parameters.number_of_species; i++) { - cuda_safe_mem(cudaMemcpy(densities.data(), ek_parameters.rho[i], - densities.size() * sizeof(float), - cudaMemcpyDeviceToHost)); - - if (!fout.write(reinterpret_cast(densities.data()), nchars)) { - fout.close(); - return 1; - } - } - - fout.close(); - - lb_lbfluid_save_checkpoint(lb_filename, true); - return 0; -} - -int ek_load_checkpoint(char *filename) { - std::string fname(filename); - std::ifstream fin((const char *)(fname + ".ek").c_str(), - std::ifstream::binary); - std::vector densities(ek_parameters.number_of_nodes); - auto const nchars = - static_cast(densities.size() * sizeof(float)); - - for (unsigned i = 0; i < ek_parameters.number_of_species; i++) { - if (!fin.read(reinterpret_cast(densities.data()), nchars)) { - fin.close(); - return 1; - } - - cuda_safe_mem(cudaMemcpy(ek_parameters.rho[i], densities.data(), - densities.size() * sizeof(float), - cudaMemcpyHostToDevice)); - } - - fin.close(); - - lb_lbfluid_load_checkpoint((char *)(fname + ".lb").c_str(), true); - - ek_integrate_electrostatics(); - - return 0; -} - void ek_set_rng_state(uint64_t counter) { if (ek_initialized) philox_counter = Utils::Counter(counter); diff --git a/src/python/espressomd/electrokinetics.pxd b/src/python/espressomd/electrokinetics.pxd index 40cad158358..0fe5ef57001 100644 --- a/src/python/espressomd/electrokinetics.pxd +++ b/src/python/espressomd/electrokinetics.pxd @@ -136,33 +136,31 @@ IF ELECTROKINETICS and CUDA: int ek_lb_print_vtk_density(char * filename) int ek_lb_print_vtk_velocity(char * filename) 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_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) - int ek_set_gamma_odd(float gamma_odd) - int ek_set_gamma_even(float gamma_even) - int ek_set_density(int species, float density) - int ek_set_D(int species, float D) - int ek_set_valency(int species, float valency) - int ek_set_ext_force_density(int species, float ext_force_density_x, float ext_force_density_y, float ext_force_density_z) - int ek_set_stencil(int stencil) - int ek_set_advection(bool advection) - int ek_set_fluctuations(bool fluctuations) - int ek_set_fluctuation_amplitude(float fluctuation_amplitude) - int ek_set_fluidcoupling(bool ideal_contribution) + void ek_set_agrid(float agrid) except + + void ek_set_lb_density(float lb_density) except + + void ek_set_viscosity(float viscosity) except + + void ek_set_friction(float friction) except + + void ek_set_lb_ext_force_density(float lb_ext_force_dens_x, float lb_ext_force_dens_y, float lb_ext_force_dens_z) except + + void ek_set_T(float T) except + + void ek_set_prefactor(float prefactor) except + + void ek_set_bulk_viscosity(float bulk_viscosity) except + + void ek_set_gamma_odd(float gamma_odd) except + + void ek_set_gamma_even(float gamma_even) except + + void ek_set_density(int species, float density) + void ek_set_D(int species, float D) + void ek_set_valency(int species, float valency) + void ek_set_ext_force_density(int species, float ext_force_density_x, float ext_force_density_y, float ext_force_density_z) + void ek_set_stencil(int stencil) except + + void ek_set_advection(bool advection) except + + void ek_set_fluctuations(bool fluctuations) except + + void ek_set_fluctuation_amplitude(float fluctuation_amplitude) except + + void ek_set_fluidcoupling(bool ideal_contribution) except + + void ek_set_electrostatics_coupling(bool electrostatics_coupling) except + int ek_node_get_density(int species, int x, int y, int z, double * density) int ek_node_get_flux(int species, int x, int y, int z, double * flux) int ek_node_get_potential(int x, int y, int z, double * potential) int ek_node_set_density(int species, int x, int y, int z, double density) float ek_calculate_net_charge() int ek_neutralize_system(int species) - int ek_save_checkpoint(char * filename, char * lb_filename) - int ek_load_checkpoint(char * filename) - int ek_set_electrostatics_coupling(bool electrostatics_coupling) int ek_print_vtk_particle_potential(char * filename) diff --git a/src/python/espressomd/electrokinetics.pyx b/src/python/espressomd/electrokinetics.pyx index 730a144ee56..9d4a6fe0f8c 100644 --- a/src/python/espressomd/electrokinetics.pyx +++ b/src/python/espressomd/electrokinetics.pyx @@ -343,25 +343,16 @@ IF ELECTROKINETICS: raise Exception("'es_coupling' is not active.") def save_checkpoint(self, path): - tmp_path = path + ".__tmp__" - tmpfile_ek = tempfile.NamedTemporaryFile(mode='w+b') - tmpfile_lb = tempfile.NamedTemporaryFile(mode='w+b') - ek_save_checkpoint(utils.to_char_pointer(tmpfile_ek.name), - utils.to_char_pointer(tmpfile_lb.name)) - ek_path = tmp_path + ".ek" - lb_path = tmp_path + ".lb" - shutil.copy(tmpfile_ek.name, path + ".ek") - shutil.copy(tmpfile_lb.name, path + ".lb") + raise RuntimeError("EK does not support checkpointing") def load_checkpoint(self, path): - self._activate_method() - ek_load_checkpoint(utils.to_char_pointer(path)) + raise RuntimeError("EK does not support checkpointing") def add_reaction(self, shape): - raise Exception("This method is not implemented yet.") + raise NotImplementedError("This method is not implemented yet.") def add_boundary(self, shape): - raise Exception("This method is not implemented yet.") + raise NotImplementedError("This method is not implemented yet.") cdef class ElectrokineticsRoutines(LBFluidRoutines): @@ -387,16 +378,10 @@ IF ELECTROKINETICS: # __getstate__ and __setstate__ define the pickle interaction def __getstate__(self): - odict = {} - odict["id"] = self.id - odict.update(self._params.copy()) - return odict + raise RuntimeError("EK does not support checkpointing") def __setstate__(self, params): - self.id = params["id"] - params.pop("id") - self._params = params - self._set_params_in_es_core() + raise RuntimeError("EK does not support checkpointing") def __str__(self): return f"{self.__class__.__name__}({self.get_params()})" diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 1cc3a321fef..2f0ad8eb4d0 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -57,7 +57,7 @@ endfunction(PYTHON_TEST) # Separate features with hyphens, use a period to add an optional flag. foreach( TEST_COMBINATION - lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.elc-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt;lb.off-int.sd;lb.off-dp3m.cpu-therm.langevin-int.nvt;lb.off-therm.dpd-int.nvt;lb.off-scafacos-therm.bd-int.bd;lb.off-therm.sdm-int.sdm;lb.cpu-p3m.cpu-lj-therm.lb-mpi.1.core + lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.elc-lj-therm.lb;lb.off-therm.npt-int.npt;lb.off-int.sd;lb.off-dp3m.cpu-therm.langevin-int.nvt;lb.off-therm.dpd-int.nvt;lb.off-scafacos-therm.bd-int.bd;lb.off-therm.sdm-int.sdm;lb.cpu-p3m.cpu-lj-therm.lb-mpi.1.core ) if(${TEST_COMBINATION} MATCHES "mpi\\.1\\.core") set(TEST_NPROCS 1) diff --git a/testsuite/python/ek_charged_plate.py b/testsuite/python/ek_charged_plate.py index 4a078336d1e..6b13be4381b 100644 --- a/testsuite/python/ek_charged_plate.py +++ b/testsuite/python/ek_charged_plate.py @@ -163,6 +163,12 @@ def test(self): positive_ions[i, j, 10].density = 0.0 negative_ions[i, j, 30].density = 0.0 + # Test error when trying to change ekin parameters after initialisation + ek._params.update({'agrid': 3, + 'T': 0.01}) + with self.assertRaises(RuntimeError): + ek._set_params_in_es_core() + if __name__ == "__main__": ut.main() diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 5855fe53c06..d44f94701bf 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -28,7 +28,6 @@ import espressomd.observables import espressomd.lb import espressomd.lbboundaries -import espressomd.electrokinetics import espressomd.shapes import espressomd.constraints @@ -78,26 +77,6 @@ system.lbboundaries.add( espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=(-1, 0, 0), dist=-(system.box_l[0] - 0.5)), velocity=(0, 0, 0))) -EK_implementation = None -if 'EK.GPU' in modes and espressomd.gpu_available( -) and espressomd.has_features('ELECTROKINETICS'): - EK_implementation = espressomd.electrokinetics - ek = EK_implementation.Electrokinetics( - agrid=0.5, - lb_density=26.15, - viscosity=1.7, - friction=0.0, - T=1.1, - prefactor=0.88, - stencil="linkcentered") - ek_species = EK_implementation.Species( - density=0.4, - D=0.02, - valency=0.3, - ext_force_density=[0.01, -0.08, 0.06]) - ek.add_species(ek_species) - system.actors.add(ek) - p1 = system.part.add(id=0, pos=[1.0] * 3) p2 = system.part.add(id=1, pos=[1.0, 1.0, 2.0]) @@ -312,23 +291,6 @@ lbf_cpt_path = checkpoint.checkpoint_dir + "/lb.cpt" lbf.save_checkpoint(lbf_cpt_path, cpt_mode) -if EK_implementation: - m = np.pi / 12 - nx = int(np.round(system.box_l[0] / ek.get_params()["agrid"])) - ny = int(np.round(system.box_l[1] / ek.get_params()["agrid"])) - nz = int(np.round(system.box_l[2] / ek.get_params()["agrid"])) - # Create a 3D grid with deterministic values to fill the LB fluid lattice - grid_3D = np.fromfunction( - lambda i, j, k: np.cos(i * m) * np.cos(j * m) * np.cos(k * m), - (nx, ny, nz), dtype=float) - for i in range(nx): - for j in range(ny): - for k in range(nz): - ek_species[i, j, k].density = grid_3D[i, j, k] - # save LB checkpoint file - ek_cpt_path = checkpoint.checkpoint_dir + "/ek" - ek.save_checkpoint(ek_cpt_path) - # save checkpoint file checkpoint.save(0) diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 8f215e179fc..0c8c031259d 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -35,9 +35,6 @@ LB = ('LB.CPU' in modes or 'LB.GPU' in modes and espressomd.gpu_available()) -EK = ('EK.GPU' in modes and espressomd.gpu_available() - and espressomd.has_features('ELECTROKINETICS')) - class CheckpointTest(ut.TestCase): @@ -104,45 +101,6 @@ def test_lb_fluid(self): self.assertIn(key, state) self.assertAlmostEqual(reference[key], state[key], delta=1E-7) - @utx.skipIfMissingFeatures('ELECTROKINETICS') - @ut.skipIf(not EK, "Skipping test due to missing mode.") - def test_EK(self): - ek = system.actors[0] - ek_species = ek.get_params()['species'][0] - cpt_path = self.checkpoint.checkpoint_dir + "/ek" - ek.load_checkpoint(cpt_path) - precision = 5 - m = np.pi / 12 - nx = int(np.round(system.box_l[0] / ek.get_params()["agrid"])) - ny = int(np.round(system.box_l[1] / ek.get_params()["agrid"])) - nz = int(np.round(system.box_l[2] / ek.get_params()["agrid"])) - grid_3D = np.fromfunction( - lambda i, j, k: np.cos(i * m) * np.cos(j * m) * np.cos(k * m), - (nx, ny, nz), dtype=float) - for i in range(nx): - for j in range(ny): - for k in range(nz): - np.testing.assert_almost_equal( - np.copy(ek_species[i, j, k].density), - grid_3D[i, j, k], - decimal=precision) - state = ek.get_params() - reference = {'agrid': 0.5, 'lb_density': 26.15, - 'viscosity': 1.7, 'friction': 0.0, - 'T': 1.1, 'prefactor': 0.88, 'stencil': "linkcentered"} - for key in reference: - self.assertIn(key, state) - self.assertAlmostEqual(reference[key], state[key], delta=1E-5) - state_species = ek_species.get_params() - reference_species = {'density': 0.4, 'D': 0.02, 'valency': 0.3, - 'ext_force_density': [0.01, -0.08, 0.06]} - for key in reference_species: - self.assertIn(key, state_species) - np.testing.assert_allclose( - reference_species[key], - state_species[key], - atol=1E-5) - def test_system_variables(self): cell_system_params = system.cell_system.get_state() self.assertTrue(cell_system_params['use_verlet_list']) @@ -484,8 +442,8 @@ def test_exclusions(self): self.assertEqual(list(system.part[1].exclusions), [2]) self.assertEqual(list(system.part[2].exclusions), [0, 1]) - @ut.skipIf(not LB or EK or not (espressomd.has_features("LB_BOUNDARIES") - or espressomd.has_features("LB_BOUNDARIES_GPU")), "Missing features") + @ut.skipIf(not LB or not (espressomd.has_features("LB_BOUNDARIES") + or espressomd.has_features("LB_BOUNDARIES_GPU")), "Missing features") @ut.skipIf(n_nodes > 1, "only runs for 1 MPI rank") def test_lb_boundaries(self): # check boundaries agree on all MPI nodes