diff --git a/Makefile b/Makefile
index c0e74cf883cb7b7572865dbf9ff54f282e08e031..3f07f7168d85a17f57344cb146370db50a146310 100644
--- a/Makefile
+++ b/Makefile
@@ -1,27 +1,27 @@
-.PHONY: all build clean lj_ns_cpu lj_ns_gpu
+.PHONY: all build clean
-all: build lj_ns_cpu lj_ns_gpu
+all: build lj_cpu lj_gpu
@echo "Everything was done!"
build:
@echo "Building pairs package..."
python3 setup.py build && python3 setup.py install --user
-lj_ns_cpu:
+lj_ns.cpp:
@echo "Generating and compiling Lennard-Jones example for CPU..."
python3 examples/lj_func.py cpu
-lj_ns_gpu:
+lj_ns.cu:
@echo "Generating and compiling Lennard-Jones example for GPU..."
python3 examples/lj_func.py gpu
# Targets
-cpu: build lj_ns_cpu
- g++ -o lj_ns lj_ns.cpp
+lj_cpu: lj_ns.cpp
+ g++ -o lj_cpu lj_ns.cpp
-gpu: build lj_ns_gpu
- nvcc -o lj_ns lj_ns.cu
+lj_gpu: lj_ns.cu
+ nvcc -o lj_gpu lj_ns.cu
clean:
@echo "Cleaning..."
- rm -rf build lj_ns lj_ns.cpp lj_ns.cu dist pairs.egg-info functions functions.pdf
+ rm -rf build lj_cpu lj_gpu lj_ns.cpp lj_ns.cu dist pairs.egg-info functions functions.pdf
diff --git a/examples/lj_func.py b/examples/lj_func.py
index b75cd272dde75b166bf979fb05da6726ba48cdcd..9dc2fdf2613fb451409ca52fe03b649bdd223d93 100644
--- a/examples/lj_func.py
+++ b/examples/lj_func.py
@@ -32,13 +32,10 @@ psim.add_position('position')
psim.add_vector_property('velocity')
psim.add_vector_property('force', vol=True)
psim.from_file("data/minimd_setup_4x4x4.input", ['mass', 'position', 'velocity'])
-psim.create_cell_lists(2.8, 2.8)
-psim.create_neighbor_lists()
-psim.periodic(2.8)
+psim.build_neighbor_lists(cutoff_radius + skin)
psim.vtk_output(f"output/test_{target}")
psim.compute(lj, cutoff_radius, {'sigma6': sigma6, 'epsilon': epsilon})
psim.compute(euler, symbols={'dt': dt})
-psim.target(pairs.target_cpu())
if target == 'gpu':
psim.target(pairs.target_gpu())
diff --git a/runtime/domain/domain_partitioning.hpp b/runtime/domain/domain_partitioning.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d5dbdfc0e28e5d28c2d22407aa7270ab434eca71
--- /dev/null
+++ b/runtime/domain/domain_partitioning.hpp
@@ -0,0 +1,60 @@
+//---
+#include "pairs.hpp"
+
+#pragma once
+
+namespace pairs {
+
+template<int ndims>
+class DomainPartitioning {
+protected:
+ int world_size, rank;
+ real_t grid_min[ndims];
+ real_t grid_max[ndims];
+
+public:
+ DomainPartitioning(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
+ static_assert(ndims == 3, "Constructor called from DomainPartitioning class mismatches number of dimensions!");
+ PAIRS_ASSERT(xmax > xmin);
+ PAIRS_ASSERT(ymax > ymin);
+ PAIRS_ASSERT(zmax > zmin);
+
+ grid_min[0] = xmin;
+ grid_max[0] = xmax;
+ grid_min[1] = ymin;
+ grid_max[1] = ymax;
+ grid_min[2] = zmin;
+ grid_max[2] = zmax;
+ }
+
+ virtual void initialize();
+ virtual void finalize();
+};
+
+template<int ndims>
+class DimensionRanges : DomainPartitioning<ndims> {
+protected:
+ int nranks[ndims];
+ int prev[ndims];
+ int next[ndims];
+ int pbc_prev[ndims];
+ int pbc_next[ndims];
+ real_t subdom_min[ndims];
+ real_t subdom_max[ndims];
+
+public:
+ void fillArrays(int neighbor_ranks[], int pbc[], real_t subdom[]) {
+ for(int d = 0; d < ndims; d++) {
+ neighbor_ranks[d * 2 + 0] = this->prev[d];
+ neighbor_ranks[d * 2 + 1] = this->next[d];
+ pbc[d * 2 + 0] = this->pbc_prev[d];
+ pbc[d * 2 + 1] = this->pbc_next[d];
+ subdom[d * 2 + 0] = this->subdom_min[d];
+ subdom[d * 2 + 1] = this->subdom_max[d];
+ }
+ }
+};
+
+class ListOfBoxes : DomainPartitioning {};
+
+}
diff --git a/runtime/domain/no_partition.hpp b/runtime/domain/no_partition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..db9561f3ccf8e672816c3979f5a31aa09a9714d6
--- /dev/null
+++ b/runtime/domain/no_partition.hpp
@@ -0,0 +1,30 @@
+//---
+#include "pairs.hpp"
+#include "domain_partitioning.hpp"
+
+#pragma once
+
+namespace pairs {
+
+template <int ndims>
+class NoPartitioning : DimensionRanges<ndims> {
+public:
+ void initialize(int *argc, const char **argv) {
+ this->world_size = 1;
+ this->rank = 0;
+
+ for(int d = 0; d < ndims; d++) {
+ this->nranks[d] = 1;
+ this->prev[d] = 0;
+ this->next[d] = 0;
+ this->pbc_prev[d] = 1;
+ this->pbc_next[d] = -1;
+ this->subdom_min[d] = this->grid_min[d];
+ this->subdom_max[d] = this->grid_max[d];
+ }
+ }
+
+ void finalize() {}
+};
+
+}
diff --git a/runtime/domain/regular_6d_stencil.hpp b/runtime/domain/regular_6d_stencil.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a852dcbfee069f12932e4c7af98cf9801637eb21
--- /dev/null
+++ b/runtime/domain/regular_6d_stencil.hpp
@@ -0,0 +1,81 @@
+//---
+#include "pairs.hpp"
+#include "domain_partitioning.hpp"
+
+#pragma once
+
+namespace pairs {
+
+template <int ndims>
+class Regular6DStencil : DimensionRanges<ndims> {
+public:
+ void setConfig() {
+ static_assert(ndims == 3, "setConfig() only implemented for three dimensions!");
+ real_t area[ndims];
+ real_t best_surf = 0.0;
+
+ for(int d = 0; d < ndims; d++) {
+ this->nranks[d] = 1;
+ area[d] = (this->grid_max[d] - this->grid_min[d]) * (this->grid_max[d] - this->grid_min[d]);
+ best_surf += 2.0 * area[d];
+ }
+
+ for(int i = 1; i < world_size; i++) {
+ if(world_size % i == 0) {
+ const int rem_yz = world_size / i;
+
+ for(int j = 0; j < rem_yz; j++) {
+ if(rem_yz % j == 0) {
+ const int k = rem_yz / j;
+ const real_t surf = (area[0] / i / j) + (area[1] / i / k) + (area[2] / j / k);
+ if(surf < best_surf) {
+ this->nranks[0] = i;
+ this->nranks[1] = j;
+ this->nranks[2] = k;
+ bestsurf = surf;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void setBoundingBox() {
+ MPI_Comm cartesian;
+ int myloc[ndims];
+ int periods[ndims];
+ int rank_length[ndims];
+ int reorder = 0;
+
+ for(int d = 0; d < ndims; d++) {
+ periods[d] = 1;
+ rank_length[d] = (this->grid_max[d] - this->grid_min[d]) / this->nranks[d];
+ }
+
+ MPI_Cart_create(MPI_COMM_WORLD, ndims, this->nranks, periods, reorder, &cartesian);
+ MPI_Cart_get(cartesian, ndims, this->nranks, periods, myloc);
+ for(int d = 0; d < ndims; d++) {
+ MPI_Cart_shift(cartesian, d, 1, &(this->prev[d]), &(this->next[d]));
+ this->pbc_prev[d] = (myloc[d] == 0) ? 1 : 0;
+ this->pbc_next[d] = (myloc[d] == this->nranks[d]) ? -1 : 0;
+ this->subdom_min[d] = this->grid_min[d] + rank_length[d] * myloc[d];
+ this->subdom_max[d] = this->subdom_min[d] + rank_length[d];
+ }
+
+ MPI_Comm_free(cartesian);
+ }
+
+ void initialize(int *argc, const char **argv) {
+ MPI_Init(argc, argv);
+ MPI_Comm_size(MPI_COMM_WORLD, &(this->world_size));
+ MPI_Comm_rank(MPI_COMM_WORLD, &(this->rank));
+ this->setConfig();
+ this->setBoundingBox();
+ }
+
+ void finalize() {
+ MPI_Finalize();
+ }
+};
+
+}
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 8c188b73ab439c6b3b8b11b53862d889467e5f13..a79fea8167dedb88342ea83017a282aa01d73e1e 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -243,19 +243,19 @@ public:
}
};
-class PairsSim {
+class PairsSimulation {
private:
std::vector<Property> properties;
std::vector<Array> arrays;
DeviceFlags *prop_flags, *array_flags;
int nprops, narrays;
public:
- PairsSim(int nprops_, int narrays_) {
+ PairsSimulation(int nprops_, int narrays_) {
prop_flags = new DeviceFlags(nprops_);
array_flags = new DeviceFlags(narrays_);
}
- ~PairsSim() {
+ ~PairsSimulation() {
delete prop_flags;
delete array_flags;
}
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index d28d2ef6e46436f9585bf3d3e0e21e24eb0428f9..ddbc65e517cac5e13bbc56548cd8dbc10771d63f 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -9,12 +9,12 @@
namespace pairs {
-size_t read_particle_data(PairsSim *ps, const char *filename, double *grid_buffer, const property_t properties[], size_t nprops) {
+size_t read_particle_data(PairsSimulation *ps, const char *filename, double *grid_buffer, const property_t properties[], size_t nprops) {
std::ifstream in_file(filename, std::ifstream::in);
std::string line;
size_t n = 0;
int read_grid_data = 0;
- // TODO: store this from PairsSim class
+ // TODO: store this from PairsSimulation class
const int num_dims = 3;
if(in_file.is_open()) {
diff --git a/runtime/vtk.hpp b/runtime/vtk.hpp
index d7d84f11a9b40958230177f7f95bf40614aac502..33d598af1bb63400851d5dff4f0ef3ddd1c76958 100644
--- a/runtime/vtk.hpp
+++ b/runtime/vtk.hpp
@@ -8,7 +8,7 @@
namespace pairs {
-void vtk_write_data(PairsSim *ps, const char *filename, int start, int end, int timestep) {
+void vtk_write_data(PairsSimulation *ps, const char *filename, int start, int end, int timestep) {
std::string output_filename(filename);
std::ostringstream filename_oss;
filename_oss << filename << "_" << timestep << ".vtk";
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index ad50e3e7235c6184ebda49ffc68c776f755ee4f5..dfec0dbf1248447669b8ef6488cb5c036e6609f7 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -90,7 +90,7 @@ class CGen:
nprops = module.sim.properties.nprops()
narrays = module.sim.arrays.narrays()
self.print("int main() {")
- self.print(f" PairsSim *ps = new PairsSim({nprops}, {narrays});")
+ self.print(f" PairsSimulation *pairs = new PairsSimulation({nprops}, {narrays});")
self.generate_statement(module.block)
self.print(" return 0;")
self.print("}")
@@ -263,54 +263,54 @@ class CGen:
array_name = ast_node.array.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->copyArrayToDevice({array_id}); // {array_name}")
+ self.print(f"pairs->copyArrayToDevice({array_id}); // {array_name}")
else:
- self.print(f"ps->copyArrayToHost({array_id}); // {array_name}")
+ self.print(f"pairs->copyArrayToHost({array_id}); // {array_name}")
if isinstance(ast_node, CopyProperty):
prop_id = ast_node.prop.id()
prop_name = ast_node.prop.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->copyPropertyToDevice({prop_id}); // {prop_name}")
+ self.print(f"pairs->copyPropertyToDevice({prop_id}); // {prop_name}")
else:
- self.print(f"ps->copyPropertyToHost({prop_id}); // {prop_name}")
+ self.print(f"pairs->copyPropertyToHost({prop_id}); // {prop_name}")
if isinstance(ast_node, ClearArrayFlag):
array_id = ast_node.array.id()
array_name = ast_node.array.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->clearArrayDeviceFlag({array_id}); // {array_name}")
+ self.print(f"pairs->clearArrayDeviceFlag({array_id}); // {array_name}")
else:
- self.print(f"ps->clearArrayHostFlag({array_id}); // {array_name}")
+ self.print(f"pairs->clearArrayHostFlag({array_id}); // {array_name}")
if isinstance(ast_node, ClearPropertyFlag):
prop_id = ast_node.prop.id()
prop_name = ast_node.prop.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->clearPropertyDeviceFlag({prop_id}); // {prop_name}")
+ self.print(f"pairs->clearPropertyDeviceFlag({prop_id}); // {prop_name}")
else:
- self.print(f"ps->clearPropertyHostFlag({prop_id}); // {prop_name}")
+ self.print(f"pairs->clearPropertyHostFlag({prop_id}); // {prop_name}")
if isinstance(ast_node, SetArrayFlag):
array_id = ast_node.array.id()
array_name = ast_node.array.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->setArrayDeviceFlag({array_id}); // {array_name}")
+ self.print(f"pairs->setArrayDeviceFlag({array_id}); // {array_name}")
else:
- self.print(f"ps->setArrayHostFlag({array_id}); // {array_name}")
+ self.print(f"pairs->setArrayHostFlag({array_id}); // {array_name}")
if isinstance(ast_node, SetPropertyFlag):
prop_id = ast_node.prop.id()
prop_name = ast_node.prop.name()
if ast_node.context() == Contexts.Device:
- self.print(f"ps->setPropertyDeviceFlag({prop_id}); // {prop_name}")
+ self.print(f"pairs->setPropertyDeviceFlag({prop_id}); // {prop_name}")
else:
- self.print(f"ps->setPropertyHostFlag({prop_id}); // {prop_name}")
+ self.print(f"pairs->setPropertyHostFlag({prop_id}); // {prop_name}")
if isinstance(ast_node, For):
iterator = self.generate_expression(ast_node.iterator)
@@ -414,7 +414,7 @@ class CGen:
size = self.generate_expression(ast_node.size())
if a.is_static():
- self.print(f"ps->addStaticArray({a.id()}, \"{a.name()}\", {ptr}, {d_ptr}, {size});")
+ self.print(f"pairs->addStaticArray({a.id()}, \"{a.name()}\", {ptr}, {d_ptr}, {size});")
else:
if self.target.is_gpu() and a.device_flag:
@@ -423,7 +423,7 @@ class CGen:
else:
self.print(f"{tkw} *{ptr};")
- self.print(f"ps->addArray({a.id()}, \"{a.name()}\", &{ptr}, {d_ptr}, {size});")
+ self.print(f"pairs->addArray({a.id()}, \"{a.name()}\", &{ptr}, {d_ptr}, {size});")
if isinstance(ast_node, RegisterProperty):
p = ast_node.property()
@@ -449,7 +449,7 @@ class CGen:
else:
self.print(f"{tkw} *{ptr};")
- self.print(f"ps->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
+ self.print(f"pairs->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
if isinstance(ast_node, Timestep):
self.generate_statement(ast_node.block)
@@ -459,16 +459,16 @@ class CGen:
ptr = p.name()
d_ptr_addr = f"&d_{ptr}" if self.target.is_gpu() and p.device_flag else "nullptr"
sizes = ", ".join([str(self.generate_expression(BinOp.inline(size))) for size in ast_node.sizes()])
- self.print(f"ps->reallocProperty({p.id()}, &{ptr}, {d_ptr_addr}, {sizes});")
- #self.print(f"ps->reallocProperty({p.id()}, (void **) &{ptr}, (void **) &d_{ptr}, {sizes});")
+ self.print(f"pairs->reallocProperty({p.id()}, &{ptr}, {d_ptr_addr}, {sizes});")
+ #self.print(f"pairs->reallocProperty({p.id()}, (void **) &{ptr}, (void **) &d_{ptr}, {sizes});")
if isinstance(ast_node, ReallocArray):
a = ast_node.array()
size = self.generate_expression(ast_node.size())
ptr = a.name()
d_ptr_addr = f"&d_{ptr}" if self.target.is_gpu() and a.device_flag else "nullptr"
- self.print(f"ps->reallocArray({a.id()}, &{ptr}, {d_ptr_addr}, {size});")
- #self.print(f"ps->reallocArray({a.id()}, (void **) &{ptr}, (void **) &d_{ptr}, {size});")
+ self.print(f"pairs->reallocArray({a.id()}, &{ptr}, {d_ptr_addr}, {size});")
+ #self.print(f"pairs->reallocArray({a.id()}, (void **) &{ptr}, (void **) &d_{ptr}, {size});")
if isinstance(ast_node, VarDecl):
tkw = Types.c_keyword(ast_node.var.type())
@@ -518,7 +518,7 @@ class CGen:
return f"e{ast_node.id()}"
if isinstance(ast_node, Call):
- params = ", ".join(["ps"] + [str(self.generate_expression(p)) for p in ast_node.parameters()])
+ params = ", ".join(["pairs"] + [str(self.generate_expression(p)) for p in ast_node.parameters()])
return f"{ast_node.name()}({params})"
if isinstance(ast_node, Cast):
diff --git a/src/pairs/ir/loops.py b/src/pairs/ir/loops.py
index aa86b7788e0bbb0c842a5d267efddcf8d1645632..46abc6188a129bdcde80e4df5ea8e0f53cc52bf9 100644
--- a/src/pairs/ir/loops.py
+++ b/src/pairs/ir/loops.py
@@ -69,14 +69,14 @@ class For(ASTNode):
class ParticleFor(For):
def __init__(self, sim, block=None, local_only=True):
- super().__init__(sim, 0, sim.nlocal if local_only else sim.nlocal + sim.pbc.npbc, block)
+ super().__init__(sim, 0, sim.nlocal if local_only else sim.nlocal + sim.comm.nghost, block)
self.local_only = local_only
def __str__(self):
return f"ParticleFor<self.iterator>"
def children(self):
- return [self.block, self.sim.nlocal] + ([] if self.local_only else [self.sim.pbc.npbc])
+ return [self.block, self.sim.nlocal] + ([] if self.local_only else [self.sim.comm.nghost])
class While(ASTNode):
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
new file mode 100644
index 0000000000000000000000000000000000000000..240305dece1e0394d09c558fe56315fab64724bd
--- /dev/null
+++ b/src/pairs/sim/comm.py
@@ -0,0 +1,98 @@
+from pairs.ir.block import pairs_device_block, pairs_host_block
+from pairs.ir.branches import Branch, Filter
+from pairs.ir.loops import For, ParticleFor
+from pairs.ir.utils import Print
+from pairs.ir.select import Select
+from pairs.ir.types import Types
+from pairs.sim.lowerable import Lowerable
+
+
+class Comm:
+ def __init__(self, sim):
+ self.sim = sim
+ self.nghost = sim.add_var('nghost', Types.Int32)
+ self.ghost_capacity = sim.add_var('ghost_capacity', Types.Int32, 100)
+ self.ghost_map = sim.add_array('ghost_map', [self.ghost_capacity], Types.Int32)
+ self.ghost_mult = sim.add_array('ghost_mult', [self.ghost_capacity, sim.ndims()], Types.Int32)
+
+
+class DetermineGhostParticles(Lowerable):
+ def __init__(self, sim, comm):
+ super().__init__(sim)
+ self.comm = comm
+
+ @pairs_host_block
+ def lower(self):
+ sim = self.sim
+ ndims = sim.ndims()
+ grid = self.sim.grid
+ cutneigh = self.sim.cell_spacing()
+ nghost = self.comm.nghost
+ ghost_capacity = self.comm.ghost_capacity
+ ghost_map = self.comm.ghost_map
+ ghost_mult = self.comm.ghost_mult
+ positions = self.sim.position()
+ nlocal = self.sim.nlocal
+ sim.module_name("setup_comm")
+ sim.check_resize(ghost_capacity, nghost)
+
+ nghost.set(0)
+ for d in range(0, ndims):
+ for i in For(sim, 0, nlocal + nghost):
+ pos = positions[i]
+ grid_length = grid.length(d)
+
+ # TODO: VecFilter?
+ for _ in Filter(sim, pos[d] < grid.min(d) + cutneigh):
+ last_pos = positions[nlocal + nghost]
+ ghost_map[nghost].set(i)
+ ghost_mult[nghost][d].set(1)
+ last_pos[d].set(pos[d] + grid_length)
+
+ for d_ in [x for x in range(0, ndims) if x != d]:
+ ghost_mult[nghost][d_].set(0)
+ last_pos[d_].set(pos[d_])
+
+ nghost.add(1)
+
+ for _ in Filter(sim, pos[d] > grid.max(d) - cutneigh):
+ last_pos = positions[nlocal + nghost]
+ ghost_map[nghost].set(i)
+ ghost_mult[nghost][d].set(-1)
+ last_pos[d].set(pos[d] - grid_length)
+
+ for d_ in [x for x in range(0, ndims) if x != d]:
+ ghost_mult[nghost][d_].set(0)
+ last_pos[d_].set(pos[d_])
+
+ nghost.add(1)
+
+
+class UpdateGhostParticles(Lowerable):
+ def __init__(self, sim, comm):
+ super().__init__(sim)
+ self.comm = comm
+
+ @pairs_device_block
+ def lower(self):
+ sim = self.sim
+ ndims = sim.ndims()
+ grid = self.sim.grid
+ nghost = self.comm.nghost
+ ghost_map = self.comm.ghost_map
+ ghost_mult = self.comm.ghost_mult
+ positions = self.sim.position()
+ nlocal = self.sim.nlocal
+ sim.module_name("update_comm")
+
+ for i in For(sim, 0, nghost):
+ # TODO: allow syntax:
+ # positions[nlocal + i].set(positions[ghost_map[i]] + ghost_mult[i] * grid.length)
+ for d in range(0, ndims):
+ positions[nlocal + i][d].set(positions[ghost_map[i]][d] + ghost_mult[i][d] * grid.length(d))
+
+
+class ExchangeParticles(Lowerable):
+ def __init__(self, sim, comm):
+ super().__init__(sim)
+ self.comm = comm
diff --git a/src/pairs/sim/pbc.py b/src/pairs/sim/pbc.py
index 3c812b395df635ea1cfbeb77c0695ad86b8c509a..47bb6676623c5d27d44ef27ba4ad0f56c740bc5c 100644
--- a/src/pairs/sim/pbc.py
+++ b/src/pairs/sim/pbc.py
@@ -7,52 +7,15 @@ from pairs.ir.types import Types
from pairs.sim.lowerable import Lowerable
-class PBC:
- def __init__(self, sim, grid, cutneigh, pbc_flags=[1, 1, 1]):
- self.sim = sim
- self.grid = grid
- self.cutneigh = cutneigh
- self.pbc_flags = pbc_flags
- self.npbc = sim.add_var('npbc', Types.Int32)
- self.pbc_capacity = sim.add_var('pbc_capacity', Types.Int32, 100)
- self.pbc_map = sim.add_array('pbc_map', [self.pbc_capacity], Types.Int32)
- self.pbc_mult = sim.add_array('pbc_mult', [self.pbc_capacity, sim.ndims()], Types.Int32)
-
-
-class UpdatePBC(Lowerable):
- def __init__(self, sim, pbc):
- super().__init__(sim)
- self.pbc = pbc
-
- @pairs_device_block
- def lower(self):
- sim = self.sim
- ndims = sim.ndims()
- grid = self.pbc.grid
- npbc = self.pbc.npbc
- pbc_map = self.pbc.pbc_map
- pbc_mult = self.pbc.pbc_mult
- positions = self.pbc.sim.position()
- nlocal = self.pbc.sim.nlocal
- sim.module_name("update_pbc")
-
- for i in For(sim, 0, npbc):
- # TODO: allow syntax:
- # positions[nlocal + i].set(positions[pbc_map[i]] + pbc_mult[i] * grid.length)
- for d in range(0, ndims):
- positions[nlocal + i][d].set(positions[pbc_map[i]][d] + pbc_mult[i][d] * grid.length(d))
-
-
class EnforcePBC(Lowerable):
- def __init__(self, sim, pbc):
+ def __init__(self, sim):
super().__init__(sim)
- self.pbc = pbc
@pairs_device_block
def lower(self):
sim = self.sim
+ grid = sim.grid
ndims = sim.ndims()
- grid = self.pbc.grid
positions = sim.position()
sim.module_name("enforce_pbc")
@@ -64,55 +27,3 @@ class EnforcePBC(Lowerable):
for _ in Filter(sim, positions[i][d] > grid.max(d)):
positions[i][d].sub(grid.length(d))
-
-
-class SetupPBC(Lowerable):
- def __init__(self, sim, pbc):
- super().__init__(sim)
- self.pbc = pbc
-
- @pairs_host_block
- def lower(self):
- sim = self.sim
- ndims = sim.ndims()
- grid = self.pbc.grid
- cutneigh = self.pbc.cutneigh
- npbc = self.pbc.npbc
- pbc_capacity = self.pbc.pbc_capacity
- pbc_map = self.pbc.pbc_map
- pbc_mult = self.pbc.pbc_mult
- positions = self.pbc.sim.position()
- nlocal = self.pbc.sim.nlocal
- sim.module_name("setup_pbc")
- sim.check_resize(pbc_capacity, npbc)
-
- npbc.set(0)
- for d in range(0, ndims):
- for i in For(sim, 0, nlocal + npbc):
- pos = positions[i]
- grid_length = grid.length(d)
-
- # TODO: VecFilter?
- for _ in Filter(sim, pos[d] < grid.min(d) + cutneigh):
- last_pos = positions[nlocal + npbc]
- pbc_map[npbc].set(i)
- pbc_mult[npbc][d].set(1)
- last_pos[d].set(pos[d] + grid_length)
-
- for d_ in [x for x in range(0, ndims) if x != d]:
- pbc_mult[npbc][d_].set(0)
- last_pos[d_].set(pos[d_])
-
- npbc.add(1)
-
- for _ in Filter(sim, pos[d] > grid.max(d) - cutneigh):
- last_pos = positions[nlocal + npbc]
- pbc_map[npbc].set(i)
- pbc_mult[npbc][d].set(-1)
- last_pos[d].set(pos[d] - grid_length)
-
- for d_ in [x for x in range(0, ndims) if x != d]:
- pbc_mult[npbc][d_].set(0)
- last_pos[d_].set(pos[d_])
-
- npbc.add(1)
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 22f96f8c60ceae9a4d40034136a4a1b5d8ad5173..971eda45ccf3835be693d66242842c17ccd94626 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -12,10 +12,11 @@ from pairs.graph.graphviz import ASTGraph
from pairs.mapping.funcs import compute
from pairs.sim.arrays import ArraysDecl
from pairs.sim.cell_lists import CellLists, CellListsBuild, CellListsStencilBuild
+from pairs.sim.comm import Comm, DetermineGhostParticles, ExchangeParticles, UpdateGhostParticles
from pairs.sim.grid import Grid2D, Grid3D
from pairs.sim.lattice import ParticleLattice
from pairs.sim.neighbor_lists import NeighborLists, NeighborListsBuild
-from pairs.sim.pbc import PBC, UpdatePBC, EnforcePBC, SetupPBC
+from pairs.sim.pbc import EnforcePBC
from pairs.sim.properties import PropertiesAlloc, PropertiesResetVolatile
from pairs.sim.read_from_file import ReadFromFile
from pairs.sim.timestep import Timestep
@@ -34,12 +35,10 @@ class Simulation:
self.arrays = Arrays(self)
self.particle_capacity = self.add_var('particle_capacity', Types.Int32, particle_capacity)
self.nlocal = self.add_var('nlocal', Types.Int32)
- self.nghost = self.add_var('nghost', Types.Int32)
self.resizes = self.add_array('resizes', 3, Types.Int32, arr_sync=False)
self.grid = None
self.cell_lists = None
self.neighbor_lists = None
- self.pbc = None
self.scope = []
self.nested_count = 0
self.nest = False
@@ -58,7 +57,8 @@ class Simulation:
self.iter_id = 0
self.vtk_file = None
self._target = None
- self.nparticles = self.nlocal + self.nghost
+ self.comm = Comm(self)
+ self.nparticles = self.nlocal + self.comm.nghost
self.properties.add_capacity(self.particle_capacity)
def add_module(self, module):
@@ -145,19 +145,15 @@ class Simulation:
self.setups.add_statement(read_object)
self.grid = read_object.grid
- def create_cell_lists(self, spacing, cutoff_radius):
- self.cell_lists = CellLists(self, self.grid, spacing, cutoff_radius)
+ def build_cell_lists(self, spacing):
+ self.cell_lists = CellLists(self, self.grid, spacing, spacing)
return self.cell_lists
- def create_neighbor_lists(self):
+ def build_neighbor_lists(self, spacing):
+ self.cell_lists = CellLists(self, self.grid, spacing, spacing)
self.neighbor_lists = NeighborLists(self.cell_lists)
return self.neighbor_lists
- def periodic(self, cutneigh, flags=[1, 1, 1]):
- self.pbc = PBC(self, self.grid, cutneigh, flags)
- self.properties.add_capacity(self.pbc.pbc_capacity)
- return self.pbc
-
def compute(self, func, cutoff_radius=None, symbols={}):
return compute(self, func, cutoff_radius, symbols)
@@ -220,12 +216,15 @@ class Simulation:
self._target = target
self.code_gen.assign_target(target)
+ def cell_spacing(self):
+ return self.cell_lists.cutoff_radius
+
def generate(self):
assert self._target is not None, "Target not specified!"
timestep = Timestep(self, self.ntimesteps, [
- (EnforcePBC(self, self.pbc), 20),
- (SetupPBC(self, self.pbc), UpdatePBC(self, self.pbc), 20),
+ (EnforcePBC(self), 20),
+ (DetermineGhostParticles(self, self.comm), UpdateGhostParticles(self, self.comm), 20),
(CellListsBuild(self, self.cell_lists), 20),
(NeighborListsBuild(self, self.neighbor_lists), 20),
PropertiesResetVolatile(self),
diff --git a/src/pairs/sim/vtk.py b/src/pairs/sim/vtk.py
index dc19a8ed88576c9d2df799322fb8c6ef13647484..0faf8ec739273bf58e9c96d0d5f3c13a854ffcaa 100644
--- a/src/pairs/sim/vtk.py
+++ b/src/pairs/sim/vtk.py
@@ -14,7 +14,7 @@ class VTKWrite(Lowerable):
@pairs_inline
def lower(self):
nlocal = self.sim.nlocal
- npbc = self.sim.pbc.npbc
- nall = nlocal + npbc
+ nghost = self.sim.comm.nghost
+ nall = nlocal + nghost
Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_local", 0, nlocal, self.timestep])
- Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_pbc", nlocal, nall, self.timestep])
+ Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_ghost", nlocal, nall, self.timestep])