From 2604e2bcae47214be7c3970585d28c08196a7432 Mon Sep 17 00:00:00 2001
From: Behzad Safaei <iwia103h@a0325.nhr.fau.de>
Date: Mon, 10 Mar 2025 22:16:38 +0100
Subject: [PATCH] Fix initDomain and PBCs
---
examples/modular/force_reduction.cpp | 4 ++--
examples/modular/sd_1.cpp | 2 +-
examples/modular/sd_2.cpp | 2 +-
examples/modular/sd_3_GPU.cu | 2 +-
examples/modular/sd_4.cpp | 3 +--
runtime/domain/block_forest.cpp | 18 +++++++-------
runtime/domain/block_forest.hpp | 3 +--
runtime/pairs.cpp | 3 +--
runtime/pairs.hpp | 2 +-
src/pairs/code_gen/cgen.py | 12 ----------
src/pairs/code_gen/interface.py | 17 ++++++--------
src/pairs/sim/domain.py | 8 -------
src/pairs/sim/domain_partitioning.py | 21 +----------------
src/pairs/sim/simulation.py | 35 ++++------------------------
14 files changed, 29 insertions(+), 103 deletions(-)
diff --git a/examples/modular/force_reduction.cpp b/examples/modular/force_reduction.cpp
index e884b3d..7e770a5 100644
--- a/examples/modular/force_reduction.cpp
+++ b/examples/modular/force_reduction.cpp
@@ -15,8 +15,8 @@ int main(int argc, char **argv) {
// Create bodies
pairs::id_t pUid = pairs::create_sphere(pairs_runtime, 0.0499, 0.0499, 0.07, 0.5, 0.5, 0 , 1000, 0.0045, 0, 0);
- // setup_sim after creating all bodies
- pairs_sim->setup_sim();
+ // setup_cells after creating all bodies
+ pairs_sim->setup_cells();
pairs_sim->update_mass_and_inertia();
// Track particle
diff --git a/examples/modular/sd_1.cpp b/examples/modular/sd_1.cpp
index e56cfc5..d7e0dca 100644
--- a/examples/modular/sd_1.cpp
+++ b/examples/modular/sd_1.cpp
@@ -21,7 +21,7 @@ int main(int argc, char **argv) {
pairs::create_sphere(pairs_runtime, 0.6, 0.6, 0.7, -2, -2, 0, 1000, 0.05, 0, 0);
pairs::create_sphere(pairs_runtime, 0.4, 0.4, 0.68, 2, 2, 0, 1000, 0.05, 0, 0);
- pairs_sim->setup_sim(0.1, 0.1, 0.1, 0.1);
+ pairs_sim->setup_cells(0.1, 0.1, 0.1, 0.1);
pairs_sim->update_mass_and_inertia();
int num_timesteps = 2000;
diff --git a/examples/modular/sd_2.cpp b/examples/modular/sd_2.cpp
index 4bdb4c8..3ae1269 100644
--- a/examples/modular/sd_2.cpp
+++ b/examples/modular/sd_2.cpp
@@ -42,7 +42,7 @@ int main(int argc, char **argv) {
pairs::create_sphere(pairs_runtime, 0.6, 0.6, 0.7, -2, -2, 0, 1000, 0.05, 0, 0);
pairs::create_sphere(pairs_runtime, 0.4, 0.4, 0.68, 2, 2, 0, 1000, 0.05, 0, 0);
- pairs_sim->setup_sim(0.1, 0.1, 0.1, 0.1);
+ pairs_sim->setup_cells(0.1, 0.1, 0.1, 0.1);
pairs_sim->update_mass_and_inertia();
int num_timesteps = 2000;
diff --git a/examples/modular/sd_3_GPU.cu b/examples/modular/sd_3_GPU.cu
index b44af84..bf95b20 100644
--- a/examples/modular/sd_3_GPU.cu
+++ b/examples/modular/sd_3_GPU.cu
@@ -65,7 +65,7 @@ int main(int argc, char **argv) {
auto pIsLocalInMyRank = [&](pairs::id_t uid){return ac->uidToIdxLocal(uid) != ac->getInvalidIdx();};
- pairs_sim->setup_sim(0.1, 0.1, 0.1, 0.1);
+ pairs_sim->setup_cells(0.1, 0.1, 0.1, 0.1);
pairs_sim->update_mass_and_inertia();
pairs_sim->communicate(0);
diff --git a/examples/modular/sd_4.cpp b/examples/modular/sd_4.cpp
index 80ec403..e83e66d 100644
--- a/examples/modular/sd_4.cpp
+++ b/examples/modular/sd_4.cpp
@@ -41,7 +41,6 @@ int main(int argc, char **argv) {
pairs_runtime->initDomain(&argc, &argv,
0, 0, 0, 20, 20, 20, // Domain bounds
- false, false, false, // PBCs --------------> TODO: runtime pbc
true // Enable dynamic load balancing (does initial refinement on a <1,1,1> blockforest)
);
@@ -61,7 +60,7 @@ int main(int argc, char **argv) {
double lcw = diameter_max * 1.01; // Linked-cell width
double interaction_radius = diameter_max;
- pairs_sim->setup_sim(lcw, lcw, lcw, interaction_radius);
+ pairs_sim->setup_cells(lcw, lcw, lcw, interaction_radius);
pairs_sim->update_mass_and_inertia();
diff --git a/runtime/domain/block_forest.cpp b/runtime/domain/block_forest.cpp
index 0851f2c..d19fcc6 100644
--- a/runtime/domain/block_forest.cpp
+++ b/runtime/domain/block_forest.cpp
@@ -24,8 +24,8 @@ namespace pairs {
BlockForest::BlockForest(
PairsRuntime *ps_,
- real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax, bool pbcx, bool pbcy, bool pbcz, bool balance_workload_) :
- DomainPartitioner(xmin, xmax, ymin, ymax, zmin, zmax), ps(ps_), globalPBC{pbcx, pbcy, pbcz}, balance_workload(balance_workload_) {
+ real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax, bool balance_workload_) :
+ DomainPartitioner(xmin, xmax, ymin, ymax, zmin, zmax), ps(ps_), balance_workload(balance_workload_) {
subdom = new real_t[ndims * 2];
}
@@ -35,8 +35,7 @@ BlockForest::BlockForest(PairsRuntime *ps_, const std::shared_ptr<walberla::bloc
DomainPartitioner(bf->getDomain().xMin(), bf->getDomain().xMax(),
bf->getDomain().yMin(), bf->getDomain().yMax(),
bf->getDomain().zMin(), bf->getDomain().zMax()),
- ps(ps_),
- globalPBC{bf->isXPeriodic(), bf->isYPeriodic(), bf->isZPeriodic()} {
+ ps(ps_) {
subdom = new real_t[ndims * 2];
mpiManager = walberla::mpi::MPIManager::instance();
world_size = mpiManager->numProcesses();
@@ -60,10 +59,8 @@ void BlockForest::updateNeighborhood() {
for(uint neigh = 0; neigh < block->getNeighborhoodSize(); ++neigh) {
auto neighbor_rank = walberla::int_c(block->getNeighborProcess(neigh));
- // Neighbor blocks that belong to the same rank should be added to
- // neighboorhood only if there's PBC along any dim, otherwise they should be skipped.
// TODO: Make PBCs work with runtime load balancing
- if((neighbor_rank != me) || globalPBC[0] || globalPBC[1] || globalPBC[2]) {
+ // if(neighbor_rank != me) {
const walberla::BlockID& neighbor_id = block->getNeighborId(neigh);
walberla::math::AABB neighbor_aabb = block->getNeighborAABB(neigh);
auto begin = blocks_pushed[neighbor_rank].begin();
@@ -73,7 +70,7 @@ void BlockForest::updateNeighborhood() {
neighborhood[neighbor_rank].push_back(neighbor_aabb);
blocks_pushed[neighbor_rank].push_back(neighbor_id);
}
- }
+ // }
}
}
@@ -264,7 +261,8 @@ void BlockForest::initialize(int *argc, char ***argv) {
auto ref_level = balance_workload ? getInitialRefinementLevel(procs) : 0;
- walberla::Vector3<bool> pbc(globalPBC[0], globalPBC[1], globalPBC[2]);
+ // PBC's are forced to true here and sperately handled when determining ghosts
+ walberla::Vector3<bool> pbc(true, true, true);
forest = walberla::blockforest::createBlockForest(domain, block_config, pbc, procs, ref_level);
@@ -298,7 +296,7 @@ void BlockForest::update() {
// PAIRS_DEBUG("Rebalance\n");
if (rank==0) std::cout << "Rebalance" << std::endl;
forest->refresh();
-}
+ }
this->updateNeighborhood();
this->setBoundingBox();
diff --git a/runtime/domain/block_forest.hpp b/runtime/domain/block_forest.hpp
index d814d02..bfda71b 100644
--- a/runtime/domain/block_forest.hpp
+++ b/runtime/domain/block_forest.hpp
@@ -40,14 +40,13 @@ private:
std::vector<double> aabbs;
PairsRuntime *ps;
real_t *subdom;
- const bool globalPBC[3];
int world_size, rank, nranks, total_aabbs;
bool balance_workload = false;
public:
BlockForest(
PairsRuntime *ps_,
- real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax, bool pbcx, bool pbcy, bool pbcz, bool balance_workload_);
+ real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax, bool balance_workload_);
BlockForest(PairsRuntime *ps_, const std::shared_ptr<walberla::blockforest::BlockForest> &bf);
diff --git a/runtime/pairs.cpp b/runtime/pairs.cpp
index d1c2d63..351aa36 100644
--- a/runtime/pairs.cpp
+++ b/runtime/pairs.cpp
@@ -15,7 +15,6 @@ namespace pairs {
void PairsRuntime::initDomain(
int *argc, char ***argv,
real_t xmin, real_t ymin, real_t zmin, real_t xmax, real_t ymax, real_t zmax,
- bool pbcx, bool pbcy, bool pbcz,
bool balance_workload) {
int mpi_initialized=0;
@@ -40,7 +39,7 @@ void PairsRuntime::initDomain(
#ifdef USE_WALBERLA
else if(dom_part_type == BlockForestPartitioning) {
- dom_part = new BlockForest(this, xmin, xmax, ymin, ymax, zmin, zmax, pbcx, pbcy, pbcz, balance_workload);
+ dom_part = new BlockForest(this, xmin, xmax, ymin, ymax, zmin, zmax, balance_workload);
}
#endif
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 14ffc5a..d582d27 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -321,7 +321,7 @@ public:
void initDomain(
int *argc, char ***argv,
real_t xmin, real_t ymin, real_t zmin, real_t xmax, real_t ymax, real_t zmax,
- bool pbcx = 0, bool pbcy = 0, bool pbcz = 0, bool balance_workload = 0);
+ bool balance_workload = false);
template<typename Domain_T>
void useDomain(const std::shared_ptr<Domain_T> &domain_ptr);
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 028c090..c8078a8 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -186,12 +186,6 @@ class CGen:
if not module.interface:
module_params += ["PairsRuntime *pairs_runtime", "struct PairsObjects *pobj"]
- if module.name=="initialize" and self.sim.create_domain_at_initialization:
- module_params += ["int argc", "char **argv"]
-
- if module.name=="set_domain":
- module_params += ["int argc", "char **argv"]
-
module_params += [f"{Types.c_keyword(self.sim, param.type())} {param.name()}" for param in module.parameters()]
print_params = ", ".join(module_params)
@@ -923,9 +917,6 @@ class CGen:
if isinstance(ast_node, ModuleCall):
module_params = ["pairs_runtime", "pobj"]
- if ast_node.module.name=="init_domain":
- module_params += ["argc", "argv"]
-
module_params += [f"{param.name()}" for param in ast_node.module.parameters()]
print_params = ", ".join(module_params)
@@ -1106,9 +1097,6 @@ class CGen:
if ast_node.name().startswith("pairs::"):
extra_params += ["pairs_runtime"]
- if ast_node.name() == "pairs_runtime->initDomain":
- extra_params += ["&argc", "&argv"]
-
params = ", ".join(extra_params + [str(self.generate_expression(p)) for p in ast_node.parameters()])
return f"{ast_node.name()}({params})"
diff --git a/src/pairs/code_gen/interface.py b/src/pairs/code_gen/interface.py
index cf59ca2..78d5a85 100644
--- a/src/pairs/code_gen/interface.py
+++ b/src/pairs/code_gen/interface.py
@@ -28,7 +28,7 @@ class InterfaceModules:
def create_all(self):
self.initialize()
- self.setup_sim()
+ self.setup_cells()
self.update_domain()
self.update_cells(self.sim.reneighbor_frequency)
self.communicate(self.sim.reneighbor_frequency)
@@ -60,6 +60,9 @@ class InterfaceModules:
PrintCode(self.sim, f"pairs_runtime = new PairsRuntime({nprops}, {ncontactprops}, {narrays}, {part});")
PrintCode(self.sim, f"pobj = new PairsObjects();")
+ if self.sim.grid is None:
+ self.sim.grid = MutableGrid(self.sim, self.sim.dims)
+
inits = Block.from_list(self.sim, [
DeclareVariables(self.sim),
DeclareArrays(self.sim),
@@ -70,16 +73,11 @@ class InterfaceModules:
RegisterMarkers(self.sim)
])
- if self.sim.create_domain_at_initialization:
- self.sim.add_statement(Block.merge_blocks(inits, self.sim.create_domain))
- else:
- assert self.sim.grid is None, "A grid already exists"
- self.sim.grid = MutableGrid(self.sim, self.sim.dims)
- self.sim.add_statement(inits)
+ self.sim.add_statement(inits)
@pairs_interface_block
- def setup_sim(self):
- self.sim.module_name('setup_sim')
+ def setup_cells(self):
+ self.sim.module_name('setup_cells')
if self.sim.cell_lists.runtime_spacing:
for d in range(self.sim.dims):
@@ -88,7 +86,6 @@ class InterfaceModules:
if self.sim.cell_lists.runtime_cutoff_radius:
Assign(self.sim, self.sim.cell_lists.cutoff_radius, Parameter(self.sim, 'cutoff_radius', Types.Real))
- self.sim.add_statement(self.sim.setup_particles)
# This update assumes all particles have been created exactly in the rank that contains them
self.sim.add_statement(UpdateDomain(self.sim))
self.sim.add_statement(BuildCellListsStencil(self.sim, self.sim.cell_lists))
diff --git a/src/pairs/sim/domain.py b/src/pairs/sim/domain.py
index 2560f1a..dfa73de 100644
--- a/src/pairs/sim/domain.py
+++ b/src/pairs/sim/domain.py
@@ -1,14 +1,6 @@
from pairs.ir.block import pairs_inline
from pairs.sim.lowerable import Lowerable
-class InitializeDomain(Lowerable):
- def __init__(self, sim):
- super().__init__(sim)
-
- @pairs_inline
- def lower(self):
- self.sim.domain_partitioning().initialize()
-
class UpdateDomain(Lowerable):
def __init__(self, sim):
super().__init__(sim)
diff --git a/src/pairs/sim/domain_partitioning.py b/src/pairs/sim/domain_partitioning.py
index e72e445..10fac21 100644
--- a/src/pairs/sim/domain_partitioning.py
+++ b/src/pairs/sim/domain_partitioning.py
@@ -45,10 +45,6 @@ class DimensionRanges:
def reduce_sum_step_indexes(self, step, array):
return sum([array[i] for i in self.step_indexes(step)])
- def initialize(self):
- grid_array = [self.sim.grid.min(d) for d in range(self.sim.ndims())] + [self.sim.grid.max(d) for d in range(self.sim.ndims())]
- Call_Void(self.sim, "pairs_runtime->initDomain", grid_array)
-
def update(self):
Call_Void(self.sim, "pairs_runtime->updateDomain", [])
Assign(self.sim, self.rank, Call_Int(self.sim, "pairs_runtime->getDomainPartitioner()->getRank", []))
@@ -140,25 +136,10 @@ class BlockForest:
return self.reduce_step
- def initialize(self):
- grid_array = [self.sim.grid.min(d) for d in range(self.sim.ndims())] + [self.sim.grid.max(d) for d in range(self.sim.ndims())]
-
- Call_Void(self.sim, "pairs_runtime->initDomain",
- grid_array + self.sim._pbc + ([True] if self.load_balancer is not None else []))
-
- if self.load_balancer is not None:
- PrintCode(self.sim, "pairs_runtime->getDomainPartitioner()->initWorkloadBalancer"
- f"({LoadBalancingAlgorithms.c_keyword(self.load_balancer)}, {self.regrid_min}, {self.regrid_max});")
-
- # Call_Void(self.sim, "pairs_runtime->getDomainPartitioner()->initWorkloadBalancer",
- # [self.load_balancer, self.regrid_min, self.regrid_max])
-
def update(self):
Call_Void(self.sim, "pairs_runtime->updateDomain", [])
Assign(self.sim, self.rank, Call_Int(self.sim, "pairs_runtime->getDomainPartitioner()->getRank", []))
Assign(self.sim, self.nranks, Call_Int(self.sim, "pairs_runtime->getNumberOfNeighborRanks", []))
-
- # for _ in Filter(self.sim, ScalarOp.neq(self.nranks, 0)): # TODO: Test different block configs with PBC
Assign(self.sim, self.ntotal_aabbs, Call_Int(self.sim, "pairs_runtime->getNumberOfNeighborAABBs", []))
for _ in Filter(self.sim, self.nranks_capacity < self.nranks):
@@ -229,7 +210,7 @@ class BlockForest:
for _ in Filter(self.sim, full_cond):
yield i, r, self.ranks[r], pbc_shifts
- for _ in Filter(self.sim, ScalarOp.cmp(self.ranks[r] , self.rank)): # if my neighbor is me (cuz I'm the only rank in a dimension that has pbc)
+ for _ in Filter(self.sim, ScalarOp.cmp(self.ranks[r] , self.rank)): # if my neighbor is me
pbc_shifts = []
isghost = Lit(self.sim, 0)
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 432d7e5..9ac9977 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -17,7 +17,7 @@ from pairs.sim.comm import Comm, Synchronize, Borders, Exchange, ReverseComm
from pairs.sim.contact_history import ContactHistory, BuildContactHistory, ClearUnusedContactHistory, ResetContactHistoryUsageStatus
from pairs.sim.copper_fcc_lattice import CopperFCCLattice
from pairs.sim.dem_sc_grid import DEMSCGrid
-from pairs.sim.domain import InitializeDomain, UpdateDomain
+from pairs.sim.domain import UpdateDomain
from pairs.sim.domain_partitioners import DomainPartitioners
from pairs.sim.domain_partitioning import BlockForest, DimensionRanges
from pairs.sim.load_balancing_algorithms import LoadBalancingAlgorithms
@@ -93,11 +93,6 @@ class Simulation:
self._capture_statements = True
self._block = Block(self, [])
- # Different segments of particle code/functions
- self.create_domain = Block(self, [])
- self.create_domain_at_initialization = False
-
- self.setup_particles = Block(self, [])
self.module_list = []
self.kernel_list = []
@@ -311,35 +306,12 @@ class Simulation:
return self.vars.find(var_name)
def set_domain(self, grid):
- """Set domain bounds.
- If the domain is set through this function, the 'set_domain' module won't be generated in the modular version.
- Use this function only if you do not need to set domain at runtime.
- This function is required only for whole-program generation."""
- self.create_domain_at_initialization = True
+ """Set domain bounds if they are known at P4IRS compile time"""
self.grid = Grid3D(self, grid[0], grid[1], grid[2], grid[3], grid[4], grid[5])
- self.create_domain.add_statement(InitializeDomain(self))
def reneighbor_every(self, frequency):
self.reneighbor_frequency = frequency
- def create_particle_lattice(self, grid, spacing, props={}):
- self.setup_particles.add_statement(ParticleLattice(self, grid, spacing, props, self.position()))
-
- def read_particle_data(self, filename, prop_names, shape_id):
- """Generate statement to read particle data from file"""
- props = [self.property(prop_name) for prop_name in prop_names]
- self.setup_particles.add_statement(ReadParticleData(self, filename, props, shape_id))
-
- def copper_fcc_lattice(self, nx, ny, nz, rho, temperature, ntypes):
- """Specific initialization for MD Copper FCC lattice case"""
- self.setup_particles.add_statement(CopperFCCLattice(self, nx, ny, nz, rho, temperature, ntypes))
-
- def dem_sc_grid(self, xmax, ymax, zmax, spacing, diameter, min_diameter, max_diameter, initial_velocity, particle_density, ntypes):
- """Specific initialization for DEM grid"""
- self.setup_particles.add_statement(
- DEMSCGrid(self, xmax, ymax, zmax, spacing, diameter, min_diameter, max_diameter,
- initial_velocity, particle_density, ntypes))
-
def build_cell_lists(self, spacing=None, store_neighbors_per_cell=False):
"""Add routines to build the linked-cells acceleration structure.
Leave spacing as None so it can be set at runtime."""
@@ -546,7 +518,7 @@ class Simulation:
# Generate getters for the runtime functions
self.code_gen.generate_interfaces()
-
+"""
def generate_program(self):
assert self.grid, "No domain is created. Set domain bounds with 'set_domain'."
@@ -652,3 +624,4 @@ class Simulation:
# Generate getters for the runtime functions
self.code_gen.generate_interfaces()
+"""
\ No newline at end of file
--
GitLab