Fixes to particle coupling on arbitrary grid

models/flow/d3q27_cumulant_part/Dynamics.c.Rt

@@ -20,34 +20,78 @@
 #	EQ = MRT_eq(U, rho, J);
 ?>
 
+#define PI 3.141592653589793
+
 CudaDeviceFunction void CalcF(){
 	for (auto p : SyncParticleIterator(X,Y,Z)) {
-	if ((NodeType & NODE_BOUNDARY) == 0) {
-	real_t d = getRho();
+        real_t d = getRho();
         vector_t u;
-	vector_t force;
-<?R C(PV(c("u.x","u.y", "u.z")), f %*% U) ?>
-	if (p.dist < p.rad) {
-		real_t cvel_len = sqrt(p.cvel.x*p.cvel.x + p.cvel.y*p.cvel.y + p.cvel.z*p.cvel.z);
-		if (cvel_len > ParticleVelocityLimit) {
-			p.cvel.x = p.cvel.x / cvel_len * ParticleVelocityLimit;
-			p.cvel.y = p.cvel.y / cvel_len * ParticleVelocityLimit;
-			p.cvel.z = p.cvel.z / cvel_len * ParticleVelocityLimit;
-		}
-		force.x = u.x - d*p.cvel.x;
-		force.y = u.y - d*p.cvel.y;
-		force.z = u.z - d*p.cvel.z;
-		p.applyForce(force);
+        <?R C(PV(c("u.x","u.y", "u.z")), f %*% U) ?>
+
+        // if particle can be shared etween two cells
+        // can't apply this approach
+        bool is_subgrid_particle =  p.rad <= 0.5;
+        // apply subgrid forcing model
+        if (is_subgrid_particle) {
+            bool found_on_boundary = false;
+            if ((NodeType & NODE_BOUNDARY) == NODE_Wall) {
+                // might not work in case of many particles
+                bool owned_by_cell = floor(p.pos.x) == floor(X) && floor(p.pos.y) == floor(Y) && floor(p.pos.z) == floor(Z);
+                if (owned_by_cell) {
+                    real_t wrong_value = nanf("");
+                    vector_t nan_vec = {wrong_value, wrong_value, wrong_value};
+                    p.applyForce(nan_vec);
+                }
+            }
+            // no need to apply any other forces
+            if (found_on_boundary) continue;
+            real_t dist = sqrt(p.diff.x*p.diff.x + p.diff.y*p.diff.y + p.diff.z*p.diff.z);
+            bool owned_by_cell = floor(p.pos.x) == floor(X) && floor(p.pos.y) == floor(Y) && floor(p.pos.z) == floor(Z);
+
+            if (owned_by_cell) {
+                vector_t u_particle = {p.cvel.x, p.cvel.y, p.cvel.z};
+
+                vector_t u_diff;
+                u_diff.x = u_particle.x - u.x;
+                u_diff.y = u_particle.y - u.y;
+                u_diff.z = u_particle.z - u.z;
+
+                vector_t force;
+                real_t u_diff_mag = sqrt(u_diff.x*u_diff.x + u_diff.y*u_diff.y + u_diff.z*u_diff.z);
+
+                real_t Re = u_diff_mag * (2*p.rad) / nu;
+                real_t drag_coeff = 1 + 0.15*pow(Re, 0.687);
+
+                force.x = - 3*PI * (p.rad*2)*nu* d*u_diff.x * drag_coeff;
+                force.y = - 3*PI * (p.rad*2)*nu* d*u_diff.y * drag_coeff;
+                force.z = - 3*PI * (p.rad*2)*nu* d*u_diff.z * drag_coeff;
+                p.applyForce(force);
+            }
+        } else {
+            if ((NodeType & NODE_BOUNDARY) == 0) {
+                vector_t force;
+                if (p.dist < p.rad) {
+                    real_t cvel_len = sqrt(p.cvel.x*p.cvel.x + p.cvel.y*p.cvel.y + p.cvel.z*p.cvel.z);
+                    if (cvel_len > ParticleVelocityLimit) {
+                        p.cvel.x = p.cvel.x / cvel_len * ParticleVelocityLimit;
+                        p.cvel.y = p.cvel.y / cvel_len * ParticleVelocityLimit;
+                        p.cvel.z = p.cvel.z / cvel_len * ParticleVelocityLimit;
+                    }
+                    force.x = u.x - d*p.cvel.x;
+                    force.y = u.y - d*p.cvel.y;
+                    force.z = u.z - d*p.cvel.z;
+                    p.applyForce(force);
 //		p.force.x += +force.x;
-		fx         = -force.x;
+                    fx         = -force.x;
 //		p.force.y += +force.y;
-		fy         = -force.y;
+                    fy         = -force.y;
 //		p.force.z += +force.z;
-		fz         = -force.z;
-		sol += 1;
-	}
-	}
-	}
+                    fz         = -force.z;
+                    sol += 1;
+                }
+            }
+        }
+    }
 }
 
 

src/ArbLattice.cpp

@@ -476,7 +476,7 @@ void ArbLattice::initContainer() {
 }
 
 int ArbLattice::fullLatticePos(double pos) const {
-    const auto retval = std::lround(pos / connect.grid_size - .5);
+    const auto retval = std::lround(pos - .5);
     assert(retval <= std::numeric_limits<int>::max() && retval >= std::numeric_limits<int>::min());
     return static_cast<int>(retval);
 }
@@ -895,3 +895,8 @@ void ArbLattice::resetAverage(){
         }
     }
 }
+
+ArbLattice::~ArbLattice()
+{
+	RFI.Close();
+}

src/ArbLattice.hpp

@@ -79,30 +79,31 @@ class ArbLattice : public LatticeBase {
     ArbLattice(ArbLattice&&) = delete;
     ArbLattice& operator=(const ArbLattice&) = delete;
     ArbLattice& operator=(ArbLattice&&) = delete;
-    virtual ~ArbLattice() = default;
+    virtual ~ArbLattice();
 
     int reinitialize(size_t num_snaps_, const std::map<std::string, int>& setting_zones, pugi::xml_node arb_node);  /// Init if passed args differ from those passed at construction or the last call to reinitialize (avoid duplicating work)
     size_t getLocalSize() const final { return connect.chunk_end - connect.chunk_begin; }
     size_t getGlobalSize() const final { return connect.num_nodes_global; }
 
 
-    virtual std::vector<int> shape() const { return {static_cast<int>(getLocalSize())}; };
-    virtual std::vector<real_t> getQuantity(const Model::Quantity& q, real_t scale = 1) ;
-    virtual std::vector<big_flag_t> getFlags() const;
-    virtual std::vector<real_t> getField(const Model::Field& f);
-    virtual std::vector<real_t> getFieldAdj(const Model::Field& f);
-    virtual std::vector<real_t> getCoord(const Model::Coord& q, real_t scale = 1);
+    virtual std::vector<int> shape() const override { return {static_cast<int>(getLocalSize())}; };
+    virtual std::vector<real_t> getQuantity(const Model::Quantity& q, real_t scale = 1) override ;
+    virtual std::vector<big_flag_t> getFlags() const override;
+    virtual std::vector<real_t> getField(const Model::Field& f) override;
+    virtual std::vector<real_t> getFieldAdj(const Model::Field& f) override;
+    virtual std::vector<real_t> getCoord(const Model::Coord& q, real_t scale = 1) override;
 
-    virtual void setFlags(const std::vector<big_flag_t>& x);
-    virtual void setField(const Model::Field& f, const std::vector<real_t>& x);
-    virtual void setFieldAdjZero(const Model::Field& f);
+    virtual void setFlags(const std::vector<big_flag_t>& x) override;
+    virtual void setField(const Model::Field& f, const std::vector<real_t>& x) override;
+    virtual void setFieldAdjZero(const Model::Field& f) override;
 
     const ArbVTUGeom& getVTUGeom() const { return vtu_geom; }
     Span<const flag_t> getNodeTypes() const { return {node_types_host.data(), node_types_host.size()}; }  /// Get host view of node types (permuted)
     const ArbLatticeConnectivity& getConnectivity() const { return connect; }
     const std::vector<unsigned>& getLocalPermutation() const { return local_permutation; }
 
     void resetAverage();
+    lbRegion getLocalBoundingBox() const override;                                                                                          /// Compute local bounding box, assuming the arbitrary lattice is a subset of a Cartesian lattice 
 
    protected:
     ArbLatticeLauncher launcher;  /// Launcher responsible for running CUDA kernels on the lattice
@@ -152,7 +153,6 @@ class ArbLattice : public LatticeBase {
     void initCommManager();                                                                                                        /// Compute which fields need to be sent to/received from which neighbors
     void initContainer();                                                                                                          /// Initialize the data residing in launcher.container
     int fullLatticePos(double pos) const;                                                                                          /// Compute the position (in terms of lattice offsets) of a node, assuming the arbitrary lattice is a subset of a Cartesian lattice
-    lbRegion getLocalBoundingBox() const;                                                                                          /// Compute local bounding box, assuming the arbitrary lattice is a subset of a Cartesian lattice
     ArbVTUGeom makeVTUGeom() const;                                                                                                /// Compute VTU geometry
     void communicateBorder();                                                                                                      /// Send and receive border values in snap (overlapped with interior computation)
     unsigned lookupLocalGhostIndex(ArbLatticeConnectivity::Index gid) const;                                                       /// For a given ghost gid, look up its local id

src/CartLattice.h.Rt

@@ -85,6 +85,13 @@ public:
 
   const lbRegion& getLocalRegion() const { return connectivity.getLocalRegion(); }
   const lbRegion& getGlobalRegion() const { return connectivity.global_region; }
+
+  lbRegion getLocalBoundingBox() const override {
+    lbRegion reg = getLocalRegion();
+    return lbRegion(px + reg.dx, py + reg.dy, pz + reg.dz,
+                    reg.nx, reg.ny, reg.nz);
+  }
+
   size_t getLocalSize() const override { return getLocalRegion().sizeL(); }
   size_t getGlobalSize() const override { return getGlobalRegion().sizeL(); }
 

src/Handlers/acRemoteForceInterface.cpp

@@ -77,10 +77,9 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
           }
         }
 
-        const auto lattice = solver->getCartLattice();
         if (stats) {
           output("Asking for stats on RFI ( %s every %d it)\n", stats_prefix.c_str(), stats_iter);
-          lattice->RFI.enableStats(stats_prefix.c_str(), stats_iter);
+          solver->lattice->RFI.enableStats(stats_prefix.c_str(), stats_iter);
         }
 
         inter = MPMD[integrator_];
@@ -91,21 +90,18 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
         integrator = integrator_;
 
         if (use_box) {
-          lbRegion reg = lattice->getLocalRegion();
-          double px = lattice->px;
-          double py = lattice->py;
-          double pz = lattice->pz;
-          lattice->RFI.DeclareSimpleBox(
-            px + reg.dx - PART_MAR_BOX,
-            px + reg.dx + reg.nx + PART_MAR_BOX,
-            py + reg.dy - PART_MAR_BOX,
-            py + reg.dy + reg.ny + PART_MAR_BOX,
-            pz + reg.dz - PART_MAR_BOX,
-            pz + reg.dz + reg.nz + PART_MAR_BOX);
+          lbRegion reg = solver->lattice->getLocalBoundingBox();
+          solver->lattice->RFI.DeclareSimpleBox(
+            reg.dx - PART_MAR_BOX,
+            reg.dx + reg.nx + PART_MAR_BOX,
+            reg.dy - PART_MAR_BOX,
+            reg.dy + reg.ny + PART_MAR_BOX,
+            reg.dz - PART_MAR_BOX,
+            reg.dz + reg.nz + PART_MAR_BOX);
         }
 
         MPI_Barrier(MPMD.local);
-        lattice->RFI.Connect(MPMD.work,inter.work);
+        solver->lattice->RFI.Connect(MPMD.work,inter.work);
 
         return 0;
 }

src/LatticeBase.hpp

@@ -10,6 +10,7 @@
 #include "SolidContainer.h"
 #include "SyntheticTurbulence.h"
 #include "ZoneSettings.h"
+#include "Region.h"
 #include "cross.h"
 #include "unit.h"
 
@@ -93,6 +94,9 @@ class LatticeBase {
    public:
     virtual size_t getLocalSize() const = 0;
     virtual size_t getGlobalSize() const = 0;
+
+    virtual lbRegion getLocalBoundingBox() const = 0;
+
     void initLattice();  /// Called by handlers
 
     template <class F>

src/Particle.hpp

@@ -8,6 +8,7 @@
 
 struct Particle {
     real_t* particle_data;
+    vector_t pos; 
 	vector_t cvel, diff;
 	real_t rad; real_t dist;
 	CudaDeviceFunction bool in() {
@@ -19,7 +20,7 @@ struct ParticleI : Particle {
 	size_t i;
 	CudaDeviceFunction ParticleI(real_t* particle_data_, const size_t& i_, const real_t node[3]) : i(i_) {
         particle_data = particle_data_;
-		vector_t pos, vel, angvel;
+		vector_t vel, angvel;
 		rad = particle_data[i*RFI_DATA_SIZE+RFI_DATA_R];
 		pos.x = particle_data[i*RFI_DATA_SIZE+RFI_DATA_POS+0];
 		pos.y = particle_data[i*RFI_DATA_SIZE+RFI_DATA_POS+1];

src/toArb.cpp

@@ -90,9 +90,9 @@ static int writeArbLatticeNodes(const Geometry& geo,
                 if (lin_to_arb_index_map.find(current_lin_pos) == lin_to_arb_index_map.end()) continue;  // void node
 
                 // Coordinates
-                const double x_coord = (static_cast<double>(x) + .5) * spacing;
-                const double y_coord = (static_cast<double>(y) + .5) * spacing;
-                const double z_coord = (static_cast<double>(z) + .5) * spacing;
+                const double x_coord = (static_cast<double>(x) + .5) ;
+                const double y_coord = (static_cast<double>(y) + .5) ;
+                const double z_coord = (static_cast<double>(z) + .5) ;
                 file << x_coord << ' ' << y_coord << ' ' << z_coord << ' ';
 
                 // Neighbors