Skip to content
Permalink
Browse files

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

…grid_refactor
  • Loading branch information...
KaiSzuttor committed Mar 15, 2019
2 parents 09b24dd + e0db0e1 commit 1c13b9a3e0eb154d371fdd49e10636ebde4910e3
@@ -38,7 +38,7 @@
lb_fluid = espressomd.lb.LBFluidGPU(
agrid=1.0, dens=1.0, visc=1.0, tau=0.01, ext_force_density=[0, 0, 0.15], kT=1.0, seed=32)
system.actors.add(lb_fluid)
system.thermostat.set_lb(LB_fluid=lb_fluid, friction=1.0, seed=23)
system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23)
fluid_obs = espressomd.observables.CylindricalLBVelocityProfile(
center=[5.0, 5.0, 0.0],
axis='z',
@@ -45,33 +45,27 @@
box_l = 50
system = espressomd.System(box_l=[box_l] * 3)
system.set_random_state_PRNG()
#system.seed = system.cell_system.get_state()['n_nodes'] * [1234]

system.time_step = 0.01
system.cell_system.skin = 0.1

system.part.add(id=0, pos=[box_l / 2.0, box_l /
2.0, box_l / 2.0], fix=[1, 1, 1])
system.part.add(pos=[box_l / 2.0] * 3, fix=[1, 1, 1])


lb_params = {'agrid': 1, 'dens': 1, 'visc': 1, 'tau': 0.01,
'ext_force_density': [0, 0, -1.0 / (box_l**3)]}
#lbf = espressomd.lb.LBFluidGPU(**lb_params)
lbf = espressomd.lb.LBFluid(**lb_params)
system.actors.add(lbf)
system.thermostat.set_lb(LB_fluid=lbf, friction=1.0)
print(system.actors)
system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0)
print(lbf.get_params())

f_list = []
f_list = np.zeros((10, 3))
for i in range(10):
f_list.append(system.part[0].f)
f_list[i] = system.part[0].f
system.integrator.run(steps=10)
print(i)

f_list = np.array(f_list)


fig1 = plt.figure()
ax = fig1.add_subplot(111)
ax.plot(f_list[:, 0], label=r"$F_x$")
@@ -45,7 +45,6 @@ set(EspressoCore_SRC
swimmer_reaction.cpp
SystemInterface.cpp
thermostat.cpp
topology.cpp
tuning.cpp
virtual_sites.cpp
accumulators/Correlator.cpp
@@ -384,21 +384,21 @@ void coldet_do_three_particle_bond(Particle &p, Particle &p1, Particle &p2) {

#ifdef VIRTUAL_SITES_RELATIVE
void place_vs_and_relate_to_particle(const int current_vs_pid,
const Vector3d &pos, const int relate_to,
const Vector3d &pos, int relate_to,
const Vector3d &initial_pos) {

// The virtual site is placed at initial_pos which will be in the local
// node's domain. It will then be moved to its final position.
// A resort occurs after vs-based collisions anyway, which will move the vs
// into the right cell.
added_particle(current_vs_pid);
local_place_particle(current_vs_pid, initial_pos.data(), 1);
local_particles[current_vs_pid]->r.p = pos;
local_vs_relate_to(current_vs_pid, relate_to);
auto p_vs = local_place_particle(current_vs_pid, initial_pos.data(), 1);
p_vs->r.p = pos;

(local_particles[max_seen_particle])->p.is_virtual = 1;
(local_particles[max_seen_particle])->p.type =
collision_params.vs_particle_type;
local_vs_relate_to(p_vs, &get_part(relate_to));

p_vs->p.is_virtual = 1;
p_vs->p.type = collision_params.vs_particle_type;
}

void bind_at_poc_create_bond_between_vs(const int current_vs_pid,
@@ -613,35 +613,28 @@ void handle_collisions() {
// Positions of the virtual sites
bind_at_point_of_collision_calc_vs_pos(p1, p2, pos1, pos2);

auto handle_particle = [&](Particle *p, Vector3d const &pos) {
if (not p->l.ghost) {
place_vs_and_relate_to_particle(current_vs_pid, pos,
p->identity(), initial_pos);
// Particle storage locations may have changed due to
// added particle
p1 = local_particles[c.pp1];
p2 = local_particles[c.pp2];
} else {
added_particle(current_vs_pid);
}
};

// place virtual sites on the node where the base particle is not a
// ghost
if (!p1->l.ghost) {
place_vs_and_relate_to_particle(current_vs_pid, pos1, c.pp1,
initial_pos);
// Particle storage locations may have changed due to
// added particle
p1 = local_particles[c.pp1];
p2 = local_particles[c.pp2];
} else // update the books
added_particle(current_vs_pid);

handle_particle(p1, pos1);
// Increment counter
current_vs_pid++;

// Same for particle 2
if (!p2->l.ghost) {
place_vs_and_relate_to_particle(current_vs_pid, pos2, c.pp2,
initial_pos);
// Particle storage locations may have changed due to
// added particle
p1 = local_particles[c.pp1];
p2 = local_particles[c.pp2];
} else // update the books
added_particle(current_vs_pid);

handle_particle(p2, pos2);
// Increment counter
current_vs_pid++;

// Create bonds between the vs.

bind_at_poc_create_bond_between_vs(current_vs_pid, c);
@@ -687,7 +680,7 @@ void handle_collisions() {
// Vs placement happens on the node that has p1
if (!attach_vs_to.l.ghost) {
place_vs_and_relate_to_particle(
current_vs_pid, pos, attach_vs_to.p.identity, initial_pos);
current_vs_pid, pos, attach_vs_to.identity(), initial_pos);
// Particle storage locations may have changed due to
// added particle
p1 = local_particles[c.pp1];
@@ -68,7 +68,6 @@
#include "statistics.hpp"
#include "statistics_chain.hpp"
#include "swimmer_reaction.hpp"
#include "topology.hpp"
#include "virtual_sites.hpp"

#include "serialization/IA_parameters.hpp"
@@ -127,8 +126,6 @@ int n_nodes = -1;
CB(mpi_rescale_particles_slave) \
CB(mpi_bcast_cell_structure_slave) \
CB(mpi_bcast_nptiso_geom_slave) \
CB(mpi_update_mol_ids_slave) \
CB(mpi_sync_topo_part_info_slave) \
CB(mpi_send_exclusion_slave) \
CB(mpi_bcast_lb_params_slave) \
CB(mpi_bcast_cuda_global_part_vars_slave) \
@@ -776,94 +773,6 @@ void mpi_bcast_nptiso_geom_slave(int, int) {
MPI_Bcast(&nptiso.non_const_dim, 1, MPI_INT, 0, comm_cart);
}

/***************REQ_UPDATE_MOL_IDS *********************/

void mpi_update_mol_ids() {
mpi_call(mpi_update_mol_ids_slave, -1, 0);
mpi_update_mol_ids_slave(-1, 0);
}

void mpi_update_mol_ids_slave(int, int) { update_mol_ids_setchains(); }

/******************* REQ_SYNC_TOPO ********************/
int mpi_sync_topo_part_info() {
int i;
int molsize = 0;
int moltype = 0;

mpi_call(mpi_sync_topo_part_info_slave, -1, 0);
int n_mols = topology.size();
MPI_Bcast(&n_mols, 1, MPI_INT, 0, comm_cart);

for (i = 0; i < n_mols; i++) {
molsize = topology[i].part.n;
moltype = topology[i].type;

#ifdef MOLFORCES
MPI_Bcast(&(topology[i].trap_flag), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(topology[i].trap_center, 3, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].trap_spring_constant), 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].drag_constant), 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].noforce_flag), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&(topology[i].isrelative), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&(topology[i].favcounter), 1, MPI_INT, 0, comm_cart);
if (topology[i].favcounter == -1)
MPI_Bcast(topology[i].fav, 3, MPI_DOUBLE, 0, comm_cart);
/* check if any molecules are trapped */
if ((topology[i].trap_flag != 32) && (topology[i].noforce_flag != 32)) {
IsTrapped = 1;
}
#endif

MPI_Bcast(&molsize, 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&moltype, 1, MPI_INT, 0, comm_cart);
MPI_Bcast(topology[i].part.e, topology[i].part.n, MPI_INT, 0, comm_cart);
MPI_Bcast(&topology[i].type, 1, MPI_INT, 0, comm_cart);
}

sync_topo_part_info();

return 1;
}

void mpi_sync_topo_part_info_slave(int, int) {
int i;
int molsize = 0;
int moltype = 0;
int n_mols = 0;

MPI_Bcast(&n_mols, 1, MPI_INT, 0, comm_cart);
realloc_topology(n_mols);
for (i = 0; i < n_mols; i++) {

#ifdef MOLFORCES
MPI_Bcast(&(topology[i].trap_flag), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(topology[i].trap_center, 3, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].trap_spring_constant), 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].drag_constant), 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&(topology[i].noforce_flag), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&(topology[i].isrelative), 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&(topology[i].favcounter), 1, MPI_INT, 0, comm_cart);
if (topology[i].favcounter == -1)
MPI_Bcast(topology[i].fav, 3, MPI_DOUBLE, 0, comm_cart);
/* check if any molecules are trapped */
if ((topology[i].trap_flag != 32) && (topology[i].noforce_flag != 32)) {
IsTrapped = 1;
}
#endif

MPI_Bcast(&molsize, 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&moltype, 1, MPI_INT, 0, comm_cart);
topology[i].type = moltype;
topology[i].part.resize(molsize);

MPI_Bcast(topology[i].part.e, topology[i].part.n, MPI_INT, 0, comm_cart);
MPI_Bcast(&topology[i].type, 1, MPI_INT, 0, comm_cart);
}

sync_topo_part_info();
}

/******************* REQ_BCAST_LBPAR ********************/

void mpi_bcast_lb_params(LBParam field, int value) {
@@ -280,11 +280,6 @@ void mpi_bcast_nptiso_geom(void);
* a single molecule */
void mpi_update_mol_ids(void);

/** Issue REQ_SYNC_TOPO: Update the molecules ids to that they correspond to
* the topology
*/
int mpi_sync_topo_part_info(void);

/** Issue REQ_BCAST_LBPAR: Broadcast a parameter for lattice Boltzmann.
* @param[in] field References the parameter field to be broadcasted.
* The references are defined in lb.hpp
@@ -36,14 +36,12 @@ class DensityProfile : public PidProfileObservable {
std::make_pair(min_z, max_z)}};
Utils::Histogram<double, 3> histogram(n_bins, 1, limits);
for (auto const &id : ids()) {
auto const ppos = ::Vector3d(folded_position(partCfg[id]));
histogram.update(ppos);
histogram.update(folded_position(partCfg[id]));
}
histogram.normalize();
return histogram.get_histogram();
}
};

} // Namespace Observables

#endif
@@ -37,8 +37,7 @@ class FluxDensityProfile : public PidProfileObservable {
Utils::Histogram<double, 3> histogram(n_bins, 3, limits);
for (auto const &id : ids()) {
auto const ppos = ::Vector3d(folded_position(partCfg[id]));
histogram.update(ppos, ::Vector3d{{partCfg[id].m.v[0], partCfg[id].m.v[1],
partCfg[id].m.v[2]}});
histogram.update(ppos, partCfg[id].m.v);
}
histogram.normalize();
return histogram.get_histogram();
@@ -38,8 +38,7 @@ class ForceDensityProfile : public PidProfileObservable {
Utils::Histogram<double, 3> histogram(n_bins, 3, limits);
for (auto const &id : ids()) {
auto const ppos = ::Vector3d(folded_position(partCfg[id]));
histogram.update(ppos, ::Vector3d{{partCfg[id].f.f[0], partCfg[id].f.f[1],
partCfg[id].f.f[2]}});
histogram.update(ppos, partCfg[id].f.f);
}
histogram.normalize();
return histogram.get_histogram();
@@ -1185,7 +1185,7 @@ void local_remove_particle(int part) {
Particle p_destroy = extract_indexed_particle(cell, n);
}

void local_place_particle(int part, const double p[3], int _new) {
Particle *local_place_particle(int part, const double p[3], int _new) {
Particle *pt;

Vector3i i{};
@@ -1222,6 +1222,10 @@ void local_place_particle(int part, const double p[3], int _new) {
#ifdef BOND_CONSTRAINT
pt->r.p_old = pp;
#endif

assert(local_particles[part] == pt);

return pt;
}

void local_remove_all_particles() {
@@ -299,14 +299,6 @@ struct ParticleLocal {
int ghost = 0;
};

#ifdef LB
/** Data related to the lattice Boltzmann hydrodynamic coupling */
struct ParticleLatticeCoupling {
/** fluctuating part of the coupling force */
Vector3d f_random;
};
#endif

struct ParticleParametersSwimming {
// ifdef inside because we need this type for some MPI prototypes
#ifdef ENGINE
@@ -364,9 +356,6 @@ struct Particle {
ret.m = m;
ret.f = f;
ret.l = l;
#ifdef LB
ret.lc = lc;
#endif
#ifdef ENGINE
ret.swim = swim;
#endif
@@ -387,10 +376,7 @@ struct Particle {
ParticleForce f;
///
ParticleLocal l;
///
#ifdef LB
ParticleLatticeCoupling lc;
#endif
///
/** Bonded interactions list
*
* The format is pretty simple: just the bond type, and then the particle
@@ -893,8 +879,10 @@ void remove_all_bonds_to(int part);
@param part the identity of the particle to move
@param p its new position
@param _new if true, the particle is allocated, else has to exists already
@return Pointer to the particle.
*/
void local_place_particle(int part, const double p[3], int _new);
Particle *local_place_particle(int part, const double p[3], int _new);

/** Used by \ref mpi_place_particle, should not be used elsewhere.
Called if on a different node a new particle was added.
Oops, something went wrong.

0 comments on commit 1c13b9a

Please sign in to comment.
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.