From 141e3655988cd5d0f70b42ee5ce1eafc33e96844 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Wed, 21 Feb 2024 01:41:54 +0100
Subject: [PATCH] Fix bugs in blockforest code

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
 runtime/boundary_weights.hpp                  | 10 +++--
 ...taHanding.hpp => ParticleDataHandling.hpp} | 10 ++---
 runtime/domain/block_forest.cpp               | 16 +++----
 runtime/domain/block_forest.hpp               |  2 +-
 runtime/domain/domain_partitioning.hpp        | 11 ++++-
 runtime/domain/regular_6d_stencil.hpp         |  1 -
 runtime/interfaces/last_generated.hpp         | 27 +++++++++++
 runtime/pairs.cpp                             |  7 +--
 runtime/pairs.hpp                             | 45 ++++++++++++++-----
 runtime/tracked_variable.hpp                  | 18 ++++++++
 src/pairs/code_gen/cgen.py                    |  8 +++-
 src/pairs/ir/variables.py                     | 10 +++--
 src/pairs/sim/simulation.py                   |  6 +--
 13 files changed, 130 insertions(+), 41 deletions(-)
 rename runtime/domain/{ParticleDataHanding.hpp => ParticleDataHandling.hpp} (96%)
 create mode 100644 runtime/interfaces/last_generated.hpp
 create mode 100644 runtime/tracked_variable.hpp

diff --git a/runtime/boundary_weights.hpp b/runtime/boundary_weights.hpp
index f784d04..2c91aa1 100644
--- a/runtime/boundary_weights.hpp
+++ b/runtime/boundary_weights.hpp
@@ -6,12 +6,16 @@
 #include "pairs.hpp"
 #include "pairs_common.hpp"
 
+/*
 #define INTERFACE_DIR "interfaces/"
 #define INTERFACE_EXT ".hpp"
 #define INTERFACE_FILE(a, b, c) a ## b ## c
 #define INCLUDE_FILE(filename) #filename
-
 #include INCLUDE_FILE(INTERFACE_FILE(INTERFACE_DIR, APPLICATION_REFERENCE, INTERFACE_EXT))
+*/
+
+// Always include last generated interfaces
+#include "interfaces/last_generated.hpp"
 
 #pragma once
 
@@ -29,8 +33,8 @@ void compute_boundary_weights(
     int *comp_weight, int *comm_weight) {
 
     const int particle_capacity = ps->getParticleCapacity();
-    const int nlocal = ps->getNumberOfLocalParticles();
-    const int nghost = ps->getNumberOfGhostParticles();
+    const int nlocal = ps->getTrackedVariableAsInteger("nlocal");
+    const int nghost = ps->getTrackedVariableAsInteger("nghost");
     auto position_prop = ps->getPropertyByName("position");
 
     #ifndef PAIRS_TARGET_CUDA
diff --git a/runtime/domain/ParticleDataHanding.hpp b/runtime/domain/ParticleDataHandling.hpp
similarity index 96%
rename from runtime/domain/ParticleDataHanding.hpp
rename to runtime/domain/ParticleDataHandling.hpp
index 6ee2710..35d1ad3 100644
--- a/runtime/domain/ParticleDataHanding.hpp
+++ b/runtime/domain/ParticleDataHandling.hpp
@@ -101,7 +101,7 @@ private:
         }
 
         auto position = ps->getPropertyByName("position");
-        int nlocal = ps->getNumberOfLocalParticles();
+        int nlocal = ps->getTrackedVariableAsInteger("nlocal");
         int i = 0;
         int nserialized = 0;
 
@@ -163,12 +163,12 @@ private:
             }
         }
 
-        ps->setNumberOfLocalParticles(nlocal);
+        ps->setTrackedVariableAsInteger("nlocal", nlocal);
         *ptr = (uint_t) nserialized;
     }
 
     void deserializeImpl(IBlock *const, const BlockDataID&, mpi::RecvBuffer& buffer) {
-        int nlocal = ps->getNumberOfLocalParticles();
+        int nlocal = ps->getTrackedVariableAsInteger("nlocal");
         real_t real_tmp;
         int int_tmp;
         uint_t nrecv;
@@ -222,8 +222,8 @@ private:
             }
         }
 
-        ps->setNumberOfLocalParticles(nlocal + nrecv);
+        ps->setTrackedVariableAsInteger("nlocal", nlocal + nrecv);
     }
 };
 
-} // namespace walberla
\ No newline at end of file
+} // namespace walberla
diff --git a/runtime/domain/block_forest.cpp b/runtime/domain/block_forest.cpp
index 8996ed6..e6538a8 100644
--- a/runtime/domain/block_forest.cpp
+++ b/runtime/domain/block_forest.cpp
@@ -16,7 +16,7 @@
 #include "../pairs_common.hpp"
 #include "../devices/device.hpp"
 #include "regular_6d_stencil.hpp"
-#include "ParticleDataHandling.h"
+#include "ParticleDataHandling.hpp"
 
 namespace pairs {
 
@@ -78,13 +78,13 @@ void BlockForest::updateNeighborhood() {
 }
 
 void BlockForest::copyRuntimeArray(const std::string& name, void *dest, const int size) {
-    void *src = name.compare('ranks') ? ranks.data() :
-                name.compare('naabbs') ? vec_naabbs.data() :
-                name.compare('aabb_offsets') ? aabb_offsets.data() :
-                name.compare('aabbs') ? vec_aabbs.data() :
-                name.compare('subdom') ? subdom;
+    void *src = name.compare("ranks") ? ranks.data() :
+                name.compare("naabbs") ? vec_naabbs.data() :
+                name.compare("aabb_offsets") ? aabb_offsets.data() :
+                name.compare("aabbs") ? vec_aabbs.data() :
+                name.compare("subdom") ? subdom;
 
-    bool is_real = name.compare('aabbs') || name.compare('subdom');
+    bool is_real = name.compare("aabbs") || name.compare("subdom");
     int tsize = is_real ? sizeof(real_t) : sizeof(int);
     std::memcpy(dest, src, size * tsize);
 }
@@ -315,7 +315,7 @@ void BlockForest::finalize() {
     MPI_Finalize();
 }
 
-bool BlockForest::isWithinSubdomain(real_t x, real_t y, real_t z) {
+int BlockForest::isWithinSubdomain(real_t x, real_t y, real_t z) {
     for(auto& iblock: *forest) {
         auto block = static_cast<blockforest::Block *>(&iblock);
 
diff --git a/runtime/domain/block_forest.hpp b/runtime/domain/block_forest.hpp
index 7a11aa4..a0b5043 100644
--- a/runtime/domain/block_forest.hpp
+++ b/runtime/domain/block_forest.hpp
@@ -22,7 +22,7 @@ namespace pairs {
 
 class BlockForest : public DomainPartitioner {
 private:
-    std::shared_ptr<BlockForest> forest;
+    std::shared_ptr<walberla::BlockForest> forest;
     walberla::blockforest::InfoCollection info;
     std::map<int, std::vector<walberla::math::AABB>> neighborhood;
     std::map<int, std::vector<walberla::BlockID>> blocks_pushed;
diff --git a/runtime/domain/domain_partitioning.hpp b/runtime/domain/domain_partitioning.hpp
index e08e5ee..7384133 100644
--- a/runtime/domain/domain_partitioning.hpp
+++ b/runtime/domain/domain_partitioning.hpp
@@ -37,12 +37,21 @@ public:
     }
 
     virtual void initialize(int *argc, char ***argv) = 0;
-    virtual void fillArrays(int *neighbor_ranks, int *pbc, real_t *subdom) = 0;
+    virtual int getWorldSize() const = 0;
+    virtual int getRank() const = 0;
+    virtual int getNumberOfNeighborAABBs() = 0;
+    virtual int getNumberOfNeighborRanks() = 0;
+    virtual int isWithinSubdomain(real_t x, real_t y, real_t z) = 0;
+    virtual void copyRuntimeArray(const std::string& name, void *dest, const int size) = 0;
     virtual void communicateSizes(int dim, const int *nsend, int *nrecv) = 0;
     virtual void communicateData(
         int dim, int elem_size,
         const real_t *send_buf, const int *send_offsets, const int *nsend,
         real_t *recv_buf, const int *recv_offsets, const int *nrecv) = 0;
+    virtual void communicateAllData(
+        int ndims, int elem_size,
+        const real_t *send_buf, const int *send_offsets, const int *nsend,
+        real_t *recv_buf, const int *recv_offsets, const int *nrecv) = 0;
     virtual void finalize() = 0;
 };
 
diff --git a/runtime/domain/regular_6d_stencil.hpp b/runtime/domain/regular_6d_stencil.hpp
index 84c01a7..b507cd0 100644
--- a/runtime/domain/regular_6d_stencil.hpp
+++ b/runtime/domain/regular_6d_stencil.hpp
@@ -60,7 +60,6 @@ public:
 
     int isWithinSubdomain(real_t x, real_t y, real_t z);
     void copyRuntimeArray(const std::string& name, void *dest, const int size);
-    void fillArrays(int *neighbor_ranks, int *pbc, real_t *subdom);
     void communicateSizes(int dim, const int *send_sizes, int *recv_sizes);
     void communicateData(
         int dim, int elem_size,
diff --git a/runtime/interfaces/last_generated.hpp b/runtime/interfaces/last_generated.hpp
new file mode 100644
index 0000000..c8bb6b2
--- /dev/null
+++ b/runtime/interfaces/last_generated.hpp
@@ -0,0 +1,27 @@
+#include "../pairs.hpp"
+
+namespace pairs_host_interface {
+
+int get_uid(int *uid, int i) { return uid[i]; }
+int get_shape(int *shape, int i) { return shape[i]; }
+int get_flags(int *flags, int i) { return flags[i]; }
+double get_position(double *position, int i, int j, int capacity) { return position[i * 3 + j]; }
+double get_mass(double *mass, int i) { return mass[i]; }
+double get_linear_velocity(double *linear_velocity, int i, int j, int capacity) { return linear_velocity[i * 3 + j]; }
+double get_force(double *force, int i, int j, int capacity) { return force[i * 3 + j]; }
+int get_type(int *type, int i) { return type[i]; }
+
+}
+
+namespace pairs_cuda_interface {
+
+__inline__ __device__ int get_uid(int *uid, int i) { return uid[i]; }
+__inline__ __device__ int get_shape(int *shape, int i) { return shape[i]; }
+__inline__ __device__ int get_flags(int *flags, int i) { return flags[i]; }
+__inline__ __device__ double get_position(double *position, int i, int j, int capacity) { return position[i * 3 + j]; }
+__inline__ __device__ double get_mass(double *mass, int i) { return mass[i]; }
+__inline__ __device__ double get_linear_velocity(double *linear_velocity, int i, int j, int capacity) { return linear_velocity[i * 3 + j]; }
+__inline__ __device__ double get_force(double *force, int i, int j, int capacity) { return force[i * 3 + j]; }
+__inline__ __device__ int get_type(int *type, int i) { return type[i]; }
+
+}
diff --git a/runtime/pairs.cpp b/runtime/pairs.cpp
index e78e494..ab0ec47 100644
--- a/runtime/pairs.cpp
+++ b/runtime/pairs.cpp
@@ -7,6 +7,7 @@
 #include "pairs.hpp"
 #include "pairs_common.hpp"
 #include "devices/device.hpp"
+#include "domain/block_forest.hpp"
 #include "domain/regular_6d_stencil.hpp"
 
 namespace pairs {
@@ -15,13 +16,13 @@ void PairsSimulation::initDomain(
     int *argc, char ***argv,
     real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
 
-    if(dom_part_type == Regular) {
+    if(dom_part_type == RegularPartitioning) {
         const int flags[] = {1, 1, 1};
         dom_part = new Regular6DStencil(xmin, xmax, ymin, ymax, zmin, zmax, flags);
-    } else if(dom_part_type == RegularXY) {
+    } else if(dom_part_type == RegularXYPartitioning) {
         const int flags[] = {1, 1, 0};
         dom_part = new Regular6DStencil(xmin, xmax, ymin, ymax, zmin, zmax, flags);
-    } else if(dom_part_type == BlockForest) {
+    } else if(dom_part_type == BlockForestPartitioning) {
         dom_part = new BlockForest(xmin, xmax, ymin, ymax, zmin, zmax);
     } else {
         PAIRS_EXCEPTION("Domain partitioning type not implemented!\n");
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index f47e714..8edad7f 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -12,6 +12,7 @@
 #include "property.hpp"
 #include "runtime_var.hpp"
 #include "timers.hpp"
+#include "tracked_variable.hpp"
 #include "devices/device.hpp"
 #include "domain/block_forest.hpp"
 #include "domain/regular_6d_stencil.hpp"
@@ -27,13 +28,14 @@ namespace pairs {
 
 class PairsSimulation {
 private:
-    Regular6DStencil *dom_part;
-    //DomainPartitioner *dom_part;
+    //Regular6DStencil *dom_part;
+    DomainPartitioner *dom_part;
     DomainPartitioners dom_part_type;
     std::vector<Property> properties;
     std::vector<ContactProperty> contact_properties;
     std::vector<FeatureProperty> feature_properties;
     std::vector<Array> arrays;
+    std::vector<TrackedVariable> tracked_variables;
     DeviceFlags *prop_flags, *contact_prop_flags, *array_flags;
     Timers<double> *timers;
     int *nlocal, *nghost;
@@ -43,16 +45,12 @@ public:
         int nprops_,
         int ncontactprops_,
         int narrays_,
-        int *nlocal_,
-        int *nghost_,
         DomainPartitioners dom_part_type_) {
 
         dom_part_type = dom_part_type_;
         prop_flags = new DeviceFlags(nprops_);
         contact_prop_flags = new DeviceFlags(ncontactprops_);
         array_flags = new DeviceFlags(narrays_);
-        nlocal = nlocal_;
-        nghost = nghost_;
         timers = new Timers<double>(1e-6);
     }
 
@@ -70,10 +68,35 @@ public:
        return RuntimeVar<T>(h_ptr); 
     }
 
-    void setNumberOfLocalParticles(int n) { *nlocal = n; }
-    const int getNumberOfLocalParticles() { return *nlocal; }
-    void setNumberOfGhostParticles(int n) { *nghost = n; }
-    const int getNumberOfGhostParticles() { return *nghost; }
+    void trackVariable(std::string variable_name, void *ptr) {
+        auto v = std::find_if(
+            tracked_variables.begin(),
+            tracked_variables.end(),
+            [variable_name](TrackedVariable _v) { return _v.getName() == variable_name; });
+
+        PAIRS_ASSERT(v == std::end(tracked_variables));
+        tracked_variables.push_back(TrackedVariable(variable_name, ptr)); 
+    }
+
+    TrackedVariable &getTrackedVariable(std::string variable_name) {
+        auto v = std::find_if(
+            tracked_variables.begin(),
+            tracked_variables.end(),
+            [variable_name](TrackedVariable _v) { return _v.getName() == variable_name; });
+
+        PAIRS_ASSERT(v == std::end(tracked_variables));
+        return *v;
+    }
+
+    void setTrackedVariableAsInteger(std::string variable_name, int value) {
+        auto& tv = getTrackedVariable(variable_name);
+        *(static_cast<int *>(tv.getPointer())) = value;
+    }
+
+    const int getTrackedVariableAsInteger(std::string variable_name) {
+        auto& tv = getTrackedVariable(variable_name);
+        return *(static_cast<int *>(tv.getPointer()));
+    }
 
     // Arrays
     Array &getArray(array_t id);
@@ -280,7 +303,7 @@ public:
         int *argc, char ***argv,
         real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax);
 
-    Regular6DStencil *getDomainPartitioner() { return dom_part; }
+    DomainPartitioner *getDomainPartitioner() { return dom_part; }
     void communicateSizes(int dim, const int *send_sizes, int *recv_sizes);
 
     void communicateData(
diff --git a/runtime/tracked_variable.hpp b/runtime/tracked_variable.hpp
new file mode 100644
index 0000000..d985512
--- /dev/null
+++ b/runtime/tracked_variable.hpp
@@ -0,0 +1,18 @@
+#include "pairs_common.hpp"
+
+#pragma once
+
+namespace pairs {
+
+class TrackedVariable {
+protected:
+    std::string name;
+    void *ptr;
+
+public:
+    TrackedVariable(std::string name_, void *ptr_) : name(name_), ptr(ptr_) {}
+    std::string getName() { return name; }
+    void *getPointer() { return ptr; }
+};
+
+}
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 7973f44..6e1b310 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -55,7 +55,8 @@ class CGen:
         return Types.c_keyword(self.sim, Types.Real)
 
     def generate_interfaces(self):
-        self.print = Printer(f"runtime/interfaces/{self.ref}.hpp")
+        #self.print = Printer(f"runtime/interfaces/{self.ref}.hpp")
+        self.print = Printer("runtime/interfaces/last_generated.hpp")
         self.print.start()
         self.print("#include \"../pairs.hpp\"")
         self.generate_interface_namespace('pairs_host_interface')
@@ -97,7 +98,7 @@ class CGen:
         ext = ".cu" if self.target.is_gpu() else ".cpp"
         self.print = Printer(self.ref + ext)
         self.print.start()
-        self.print("#define APPLICATION_REFERENCE \"{self.ref}\"")
+        self.print(f"#define APPLICATION_REFERENCE \"{self.ref}\"")
 
         if self.target.is_gpu():
             self.print("#define PAIRS_TARGET_CUDA")
@@ -772,6 +773,9 @@ class CGen:
                 init = self.generate_expression(ast_node.var.init_value())
                 self.print(f"{tkw} {var} = {init};")
 
+                if ast_node.var.runtime_track():
+                    self.print(f"pairs->trackVariable(\"{var}\", &{var});")
+
             else:
                 for i in range(Types.number_of_elements(self.sim, ast_node.var.type())):
                     var = self.generate_expression(ast_node.var, index=i)
diff --git a/src/pairs/ir/variables.py b/src/pairs/ir/variables.py
index 3d96065..9d36cf6 100644
--- a/src/pairs/ir/variables.py
+++ b/src/pairs/ir/variables.py
@@ -18,8 +18,8 @@ class Variables:
         self.vars = []
         self.nvars = 0
 
-    def add(self, v_name, v_type, v_value=0):
-        var = Var(self.sim, v_name, v_type, v_value)
+    def add(self, v_name, v_type, v_value=0, v_runtime_track=False):
+        var = Var(self.sim, v_name, v_type, v_value, v_runtime_track)
         self.vars.append(var)
         return var
 
@@ -39,11 +39,12 @@ class Variables:
 
 
 class Var(ASTTerm):
-    def __init__(self, sim, var_name, var_type, init_value=0, temp=False):
+    def __init__(self, sim, var_name, var_type, init_value=0, runtime_track=False, temp=False):
         super().__init__(sim, OperatorClass.from_type(var_type))
         self.var_name = var_name
         self.var_type = var_type
         self.var_init_value = Lit.cvt(sim, init_value)
+        self.var_runtime_track = runtime_track
         self.var_temporary = temp
         self.mutable = True
         self.var_bonded_arrays = []
@@ -74,6 +75,9 @@ class Var(ASTTerm):
     def init_value(self):
         return self.var_init_value
 
+    def runtime_track(self):
+        return self.var_runtime_track
+
     def add_bonded_array(self, array):
         self.var_bonded_arrays.append(array)
 
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 048f309..1fe6bfc 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -57,8 +57,8 @@ class Simulation:
         self.contact_properties = ContactProperties(self)
         self.particle_capacity = self.add_var('particle_capacity', Types.Int32, particle_capacity)
         self.neighbor_capacity = self.add_var('neighbor_capacity', Types.Int32, neighbor_capacity)
-        self.nlocal = self.add_var('nlocal', Types.Int32)
-        self.nghost = self.add_var('nghost', Types.Int32)
+        self.nlocal = self.add_var('nlocal', Types.Int32, runtime=True)
+        self.nghost = self.add_var('nghost', Types.Int32, runtime=True)
         self.resizes = self.add_array('resizes', 3, Types.Int32, arr_sync=False)
         self.particle_uid = self.add_property('uid', Types.Int32, 0)
         self.particle_shape = self.add_property('shape', Types.Int32, 0)
@@ -215,7 +215,7 @@ class Simulation:
     def array(self, arr_name):
         return self.arrays.find(arr_name)
 
-    def add_var(self, var_name, var_type, init_value=0):
+    def add_var(self, var_name, var_type, init_value=0, runtime=False):
         assert self.var(var_name) is None, f"Variable already defined: {var_name}"
         return self.vars.add(var_name, var_type, init_value)
 
-- 
GitLab