Commit cc6d8292 authored by Jonas Landsgesell's avatar Jonas Landsgesell

Merge branch 'python' of https://github.com/espressomd/espresso into improve_test_coverage_WLRE

parents 56fe0156 87fd43c8
......@@ -205,6 +205,16 @@ empty:
- linux
- cuda
ubuntu:wo-dependencies:
stage: build
image: gitlab.icp.uni-stuttgart.de:4567/espressomd/docker/$CI_JOB_NAME
script:
- export myconfig=maxset make_check=false
- bash maintainer/CI/build_cmake.sh
tags:
- docker
- linux
### Builds with ROCm
rocm-maxset:
......
......@@ -71,7 +71,10 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#define SOFT_SPHERE
#define INTER_RF
#define OVERLAPPED
#ifdef P3M
#define THOLE
#endif
#define BOND_ANGLE
......
......@@ -306,33 +306,6 @@ void fold_and_reset(Particle &p) {
p.l.p_old = p.r.p;
}
/**
* @brief Extract an indexed particle from a list.
*
* Removes a particle from a particle list and
* from the particle index.
*/
Particle extract_indexed_particle(ParticleList *sl, int i) {
Particle *src = &sl->part[i];
Particle *end = &sl->part[sl->n - 1];
Particle p = std::move(*src);
assert(p.p.identity <= max_seen_particle);
local_particles[p.p.identity] = nullptr;
if (src != end) {
new (src) Particle(std::move(*end));
}
if (realloc_particlelist(sl, --sl->n)) {
update_local_particles(sl);
} else if (src != end) {
local_particles[src->p.identity] = src;
}
return p;
}
} // namespace
/**
......@@ -491,13 +464,11 @@ void cells_update_ghosts() {
}
Cell *find_current_cell(const Particle &p) {
auto c = cell_structure.position_to_cell(p.r.p);
if (c) {
return c;
} else if (!p.l.ghost) {
// Old pos must lie within the cell system
return cell_structure.position_to_cell(p.l.p_old);
} else {
assert(not resort_particles);
if (p.l.ghost) {
return nullptr;
}
return cell_structure.position_to_cell(p.l.p_old);
}
......@@ -298,13 +298,12 @@ void check_resort_particles();
/*@}*/
/* @brief Finds the cell in which a particle is stored
Uses position_to_cell on p.r.p. If this is not on the node's domain,
uses position at last Verlet list rebuild (p.l.p_old).
@return pointer to the cell or nullptr if the particle is not on the node
*/
/**
* @brief Finds the cell in which a particle is stored
*
*
* @return pointer to the cell or nullptr if the particle is not on the node
*/
Cell *find_current_cell(const Particle &p);
#endif
......@@ -2045,7 +2045,7 @@ void mpi_send_fluid_populations_slave(int node, int index) {
/****************************************************/
void mpi_bcast_max_mu() {
#ifdef DIPOLES
#if defined(DIPOLES) and defined(DP3M)
mpi_call(mpi_bcast_max_mu_slave, -1, 0);
calc_mu_max();
......@@ -2054,7 +2054,7 @@ void mpi_bcast_max_mu() {
}
void mpi_bcast_max_mu_slave(int node, int dummy) {
#ifdef DIPOLES
#if defined(DIPOLES) and defined(DP3M)
calc_mu_max();
......
......@@ -44,7 +44,7 @@
#include "particle_data.hpp"
#include "utils.hpp"
#ifdef DIPOLES
#if defined(DIPOLES) && defined(DP3M)
DLC_struct dlc_params = {1e100, 0, 0, 0, 0};
......
......@@ -42,7 +42,7 @@
#include "config.hpp"
#ifdef DIPOLES
#if defined(DIPOLES) && defined(DP3M)
/** parameters for the MDLC method */
typedef struct {
......
......@@ -245,10 +245,12 @@ void calc_long_range_energies() {
case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA:
energy.dipolar[1] = dawaanr_calculations(0, 1);
break;
#ifdef DP3M
case DIPOLAR_MDLC_DS:
energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(0, 1);
energy.dipolar[2] = add_mdlc_energy_corrections();
break;
#endif
case DIPOLAR_DS:
energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(0, 1);
break;
......
......@@ -292,9 +292,11 @@ void calc_long_range_forces() {
case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA:
dawaanr_calculations(1, 0);
break;
#ifdef DP3M
case DIPOLAR_MDLC_DS:
add_mdlc_force_corrections();
// fall through
// fall through
#endif
case DIPOLAR_DS:
magnetic_dipolar_direct_sum_calculations(1, 0);
break;
......
......@@ -214,19 +214,19 @@ Vector3d get_mi_vector(T const &a, U const &b) {
template <typename T1, typename T2, typename T3>
void fold_coordinate(T1 &pos, T2 &vel, T3 &image_box, int dir) {
if (PERIODIC(dir)) {
int img_count = (int)floor(pos[dir] * box_l_i[dir]);
image_box[dir] += img_count;
pos[dir] = pos[dir] - img_count * box_l[dir];
if (pos[dir] * box_l_i[dir] < -ROUND_ERROR_PREC ||
pos[dir] * box_l_i[dir] >= 1 + ROUND_ERROR_PREC) {
runtimeErrorMsg() << "particle coordinate out of range, pos = "
<< pos[dir] << ", image box = " << image_box[dir];
image_box[dir] = 0;
pos[dir] = 0;
return;
while ((pos[dir] < 0) && (image_box[dir] > INT_MIN)) {
pos[dir] += box_l[dir];
image_box[dir] -= 1;
}
while ((pos[dir] >= box_l[dir]) && (image_box[dir] < INT_MAX)) {
pos[dir] -= box_l[dir];
image_box[dir] += 1;
}
if ((image_box[dir] == INT_MIN) || (image_box[dir] == INT_MAX)) {
throw std::runtime_error(
"Overflow in the image box count while folding a particle coordinate "
"into the primary simulation box. Maybe a particle experienced a "
"huge force.");
}
}
}
......@@ -277,14 +277,9 @@ template <typename T1, typename T2> void fold_position(T1 &pos, T2 &image_box) {
*/
inline Vector3d folded_position(Particle const &p) {
Vector3d pos{p.r.p};
int tmp[3] = {0, 0, 0};
for (int dir = 0; dir < 3; dir++) {
if (PERIODIC(dir)) {
const int img_count =
static_cast<int>(std::floor(pos[dir] * box_l_i[dir]));
pos[dir] -= img_count * box_l[dir];
}
fold_coordinate(pos, tmp, dir);
}
return pos;
......
......@@ -472,9 +472,8 @@ int interactions_sanity_checks() {
}
#endif /* ifdef ELECTROSTATICS */
#ifdef DIPOLES
#if defined(DIPOLES) and defined(DP3M)
switch (coulomb.Dmethod) {
#ifdef DP3M
case DIPOLAR_MDLC_P3M:
if (mdlc_sanity_checks())
state = 0; // fall through
......@@ -482,7 +481,6 @@ int interactions_sanity_checks() {
if (dp3m_sanity_checks())
state = 0;
break;
#endif
case DIPOLAR_MDLC_DS:
if (mdlc_sanity_checks())
state = 0; // fall through
......
......@@ -76,6 +76,7 @@ void add_id_to_type_map(int part_id, int type);
int max_seen_particle = -1;
int n_part = 0;
bool swimming_particles_exist = false;
/**
* @brief id -> rank
*/
......@@ -220,6 +221,7 @@ void init_particlelist(ParticleList *pList) {
}
int realloc_particlelist(ParticleList *l, int size) {
assert(size >= 0);
int old_max = l->max;
Particle *old_start = l->part;
......@@ -280,6 +282,9 @@ Particle *append_indexed_particle(ParticleList *l, Particle &&part) {
}
Particle *move_unindexed_particle(ParticleList *dl, ParticleList *sl, int i) {
assert(sl->n > 0);
assert(i < sl->n);
realloc_particlelist(dl, ++dl->n);
auto dst = &dl->part[dl->n - 1];
auto src = &sl->part[i];
......@@ -290,12 +295,13 @@ Particle *move_unindexed_particle(ParticleList *dl, ParticleList *sl, int i) {
new (src) Particle(std::move(*end));
}
sl->n -= 1;
realloc_particlelist(sl, sl->n);
realloc_particlelist(sl, --sl->n);
return dst;
}
Particle *move_indexed_particle(ParticleList *dl, ParticleList *sl, int i) {
assert(sl->n > 0);
assert(i < sl->n);
int re = realloc_particlelist(dl, ++dl->n);
Particle *dst = &dl->part[dl->n - 1];
Particle *src = &sl->part[i];
......@@ -322,6 +328,41 @@ Particle *move_indexed_particle(ParticleList *dl, ParticleList *sl, int i) {
return dst;
}
/**
* @brief Extract an indexed particle from a list.
*
* Removes a particle from a particle list and
* from the particle index.
*
* @param i Index of particle to remove,
* needs to be valid.
* @param sl List to remove the particle from,
* needs to be non-empty.
* @return The extracted particle.
*/
Particle extract_indexed_particle(ParticleList *sl, int i) {
assert(sl->n > 0);
assert(i < sl->n);
Particle *src = &sl->part[i];
Particle *end = &sl->part[sl->n - 1];
Particle p = std::move(*src);
assert(p.p.identity <= max_seen_particle);
local_particles[p.p.identity] = nullptr;
if (src != end) {
new (src) Particle(std::move(*end));
}
if (realloc_particlelist(sl, --sl->n)) {
update_local_particles(sl);
} else if (src != end) {
local_particles[src->p.identity] = src;
}
return p;
}
namespace {
/* Limit cache to 100 MiB */
std::size_t const max_cache_size = (100ul * 1048576ul) / sizeof(Particle);
......@@ -803,42 +844,50 @@ int remove_particle(int p_id) {
return ES_OK;
}
void local_remove_particle(int part) {
int ind, c;
Particle *p = local_particles[part];
ParticleList *pl = nullptr, *tmp;
/* the tricky - say ugly - part: determine
the cell the particle is located in by checking
whether the particle address is inside the array */
for (c = 0; c < local_cells.n; c++) {
tmp = local_cells.cell[c];
ind = p - tmp->part;
if (ind >= 0 && ind < tmp->n) {
pl = tmp;
break;
namespace {
std::pair<Cell *, size_t> find_particle(Particle *p, Cell *c) {
for (int i = 0; i < c->n; ++i) {
if ((c->part + i) == p) {
return {c, i};
}
}
if (!pl) {
fprintf(stderr,
"%d: INTERNAL ERROR: could not find cell of particle %d, exiting\n",
this_node, part);
errexit();
}
return {nullptr, 0};
}
free_particle(p);
std::pair<Cell *, size_t> find_particle(Particle *p, CellPList cells) {
for (auto &c : cells) {
auto res = find_particle(p, c);
if (res.first) {
return res;
}
}
/* remove local_particles entry */
local_particles[p->p.identity] = nullptr;
return {nullptr, 0};
}
} // namespace
if (&pl->part[pl->n - 1] != p) {
/* move last particle to free position */
*p = pl->part[pl->n - 1];
void local_remove_particle(int part) {
Particle *p = local_particles[part];
assert(p);
assert(not p->l.ghost);
/* If the particles are sorted we can use the
* cell system to find the cell containing the
* particle. Otherwise we do a brute force search
* of the cells. */
Cell *cell = nullptr;
size_t n = 0;
if (Cells::RESORT_NONE == get_resort_particles()) {
std::tie(cell, n) = find_particle(p, find_current_cell(*p));
}
/* update the local_particles array for the moved particle */
local_particles[p->p.identity] = p;
if (not cell) {
std::tie(cell, n) = find_particle(p, local_cells);
}
pl->n--;
assert(cell && cell->part && (n < cell->n) && ((cell->part + n) == p));
Particle p_destroy = extract_indexed_particle(cell, n);
}
void local_place_particle(int part, const double p[3], int _new) {
......
......@@ -479,6 +479,8 @@ struct ParticleList {
extern int max_seen_particle;
/** total number of particles on all nodes. */
extern int n_part;
/** flag that active swimming particles exist */
extern bool swimming_particles_exist;
/** id->particle mapping on all nodes. This is used to find partners
of bonded interactions. */
......@@ -562,6 +564,8 @@ Particle *move_unindexed_particle(ParticleList *destList,
Particle *move_indexed_particle(ParticleList *destList,
ParticleList *sourceList, int ind);
Particle extract_indexed_particle(ParticleList *sl, int i);
/* Other Functions */
/************************************************/
......
......@@ -294,7 +294,7 @@ void convert_torques_propagate_omega() {
fprintf(stderr, "%d: convert_torques_propagate_omega:\n", this_node));
#if defined(LB_GPU) && defined(ENGINE)
if (lattice_switch & LATTICE_LB_GPU) {
if ((lattice_switch & LATTICE_LB_GPU) && swimming_particles_exist) {
copy_v_cs_from_GPU(local_cells.particles());
}
#endif
......
......@@ -73,6 +73,7 @@ cdef extern from "domain_decomposition.hpp":
cdef extern from "particle_data.hpp":
extern int n_part
extern bool swimming_particles_exist
cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
double dpd_gamma
......
......@@ -23,7 +23,7 @@ include "myconfig.pxi"
from espressomd.system cimport *
from espressomd.utils cimport *
IF DIPOLES == 1:
IF DIPOLES and DP3M:
cdef extern from "electrostatics_magnetostatics/mdlc_correction.hpp":
ctypedef struct dlc_struct "DLC_struct":
......
......@@ -22,7 +22,7 @@ from . cimport utils
include "myconfig.pxi"
from .actors import Actor
IF DIPOLES == 1:
IF DIPOLES and DP3M:
class MagnetostaticExtension(Actor):
pass
......
......@@ -28,7 +28,7 @@ from .interactions import BondedInteraction
from .interactions import BondedInteractions
from .interactions cimport bonded_ia_params
from copy import copy
from globals cimport max_seen_particle, time_step, box_l, n_part, n_rigidbonds, max_seen_particle_type
from globals cimport max_seen_particle, time_step, box_l, n_part, n_rigidbonds, max_seen_particle_type, swimming_particles_exist
import collections
import functools
import types
......@@ -1140,22 +1140,16 @@ cdef class ParticleHandle(object):
for e in self.exclusions:
self.delete_exclusion(e)
# Empty list? ony delete
if _partners:
nlvl = nesting_level(_partners)
if nlvl == 0: # Single item
self.add_exclusion(_partners)
elif nlvl == 1: # List of items
for partner in _partners:
self.add_exclusion(partner)
else:
raise ValueError(
"Exclusions have to specified as a lists of partners or a single partner.")
nlvl = nesting_level(_partners)
# Set new exclusion list
# self.add_exclusion(_partners)
if nlvl == 0: # Single item
self.add_exclusion(_partners)
elif nlvl == 1: # List of items
for partner in _partners:
self.add_exclusion(partner)
else:
raise ValueError(
"Exclusions have to be specified as a lists of partners or a single item.")
def __get__(self):
self.update_particle_data()
......@@ -1273,6 +1267,7 @@ cdef class ParticleHandle(object):
"""
def __set__(self, _params):
global swimming_particles_exist
cdef particle_parameters_swimming swim
swim.swimming = True
swim.v_swim = 0.0
......@@ -1325,6 +1320,9 @@ cdef class ParticleHandle(object):
swim.rotational_friction = _params[
'rotational_friction']
if swim.f_swim != 0 or swim.v_swim != 0:
swimming_particles_exist = True
if set_particle_swimming(self._id, swim) == 1:
raise Exception("Set particle position first.")
......@@ -1809,9 +1807,19 @@ cdef class ParticleList(object):
return odict
def __setstate__(self, params):
exclusions = collections.OrderedDict()
for particle_number in params.keys():
params[particle_number]["id"] = particle_number
IF EXCLUSIONS:
exclusions[
particle_number] = params[
particle_number][
"exclusions"]
del params[particle_number]["exclusions"]
self._place_new_particle(params[particle_number])
IF EXCLUSIONS:
for pid in exclusions:
self[pid].exclusions = exclusions[pid]
def __len__(self):
return n_part
......
......@@ -310,7 +310,7 @@ def nesting_level(obj):
"""
if not isinstance(obj, (list, tuple)):
if not isinstance(obj, (list, tuple, np.ndarray)):
return 0
obj = list(obj)
......
......@@ -40,6 +40,8 @@ if espressomd.has_features('LB'):
system.part.add(pos=[1.0] * 3)
system.part.add(pos=[1.0, 1.0, 2.0])
if espressomd.has_features('EXCLUSIONS'):
system.part.add(pos=[2.0] * 3, exclusions=[0, 1])
if espressomd.has_features('ELECTROSTATICS'):
system.part[0].q = 1
system.part[1].q = -1
......
......@@ -21,6 +21,7 @@ import numpy as np
import espressomd
import espressomd.checkpointing
import espressomd.virtual_sites
import tests_common
class CheckpointTest(ut.TestCase):
......@@ -120,6 +121,15 @@ class CheckpointTest(ut.TestCase):
self.assertAlmostEqual(coldet.distance, 0.11, delta=1E-9)
self.assertTrue(coldet.bond_centers, system.bonded_inter[0])
@ut.skipIf(not espressomd.has_features("EXCLUSIONS"), "Skipped because feature EXCLUSIONS missing.")
def test_exclusions(self):
self.assertTrue(
tests_common.lists_contain_same_elements(system.part[0].exclusions, [2]))
self.assertTrue(
tests_common.lists_contain_same_elements(system.part[1].exclusions, [2]))
self.assertTrue(
tests_common.lists_contain_same_elements(system.part[2].exclusions, [0, 1]))
if __name__ == '__main__':
ut.main()
......@@ -139,8 +139,8 @@ def lj_force_vector(v_d, d, lj_params):
def verify_lj_forces(system, tolerance, ids_to_skip=[]):
"""Goes over all pairs of paritcles in system and compares the forces on them
to what would be expected based on the systems lj parametes.
"""Goes over all pairs of particles in system and compares the forces on them
to what would be expected based on the systems LJ parametes.
Particle ids listed in ids_to_skip are not checked
Do not run this with a thermostat enabled."""
......@@ -594,3 +594,7 @@ def single_component_maxwell(x1, x2, kT):
"""Integrate the probability density from x1 to x2 using the trapezoidal rule"""
x = np.linspace(x1, x2, 1000)
return np.trapz(np.exp(-x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT)
def lists_contain_same_elements(list1, list2):
return len(list1) == len(list2) and sorted(list1) == sorted(list2)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment