From 8229905f3713ab4aa9729a1b5f5412a9616e7223 Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Wed, 25 Feb 2026 15:42:09 +1000 Subject: [PATCH 1/6] Adjust for particles computation --- src/ArbLattice.hpp | 2 +- src/CartLattice.h.Rt | 7 +++++++ src/Handlers/acRemoteForceInterface.cpp | 24 ++++++++++-------------- src/LatticeBase.hpp | 4 ++++ 4 files changed, 22 insertions(+), 15 deletions(-) diff --git a/src/ArbLattice.hpp b/src/ArbLattice.hpp index 4883ff4b6..62658cb9c 100644 --- a/src/ArbLattice.hpp +++ b/src/ArbLattice.hpp @@ -101,6 +101,7 @@ class ArbLattice : public LatticeBase { const std::vector& getLocalPermutation() const { return local_permutation; } void resetAverage(); + lbRegion getLocalBoundingBox() const; /// 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 @@ -149,7 +150,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 diff --git a/src/CartLattice.h.Rt b/src/CartLattice.h.Rt index bebc6f05f..783fc6c05 100644 --- a/src/CartLattice.h.Rt +++ b/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 { + 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(); } diff --git a/src/Handlers/acRemoteForceInterface.cpp b/src/Handlers/acRemoteForceInterface.cpp index ca0a0ce14..88d94b814 100644 --- a/src/Handlers/acRemoteForceInterface.cpp +++ b/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; } diff --git a/src/LatticeBase.hpp b/src/LatticeBase.hpp index 275d3ec70..292952ef0 100644 --- a/src/LatticeBase.hpp +++ b/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 From b6dc1f63712b4b5c6f72684519e5dd15faccefa6 Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Fri, 27 Feb 2026 10:29:08 +1000 Subject: [PATCH 2/6] Fix coordinate units between lattices On arbitrary lattice getCoord returns in real units, whereas in cartesian lattice it is lattice units. Then X, Y,Z macros on node return in real units which is not expected, and maybe it impossible to find particles (because particles are using lattice units). Need to also tweak vtk export --- src/ArbLattice.cpp | 2 +- src/Particle.hpp | 3 ++- src/toArb.cpp | 6 +++--- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/ArbLattice.cpp b/src/ArbLattice.cpp index 23129f962..bf9575899 100644 --- a/src/ArbLattice.cpp +++ b/src/ArbLattice.cpp @@ -425,7 +425,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::max() && retval >= std::numeric_limits::min()); return static_cast(retval); } diff --git a/src/Particle.hpp b/src/Particle.hpp index 91b70841c..0d2913862 100644 --- a/src/Particle.hpp +++ b/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]; diff --git a/src/toArb.cpp b/src/toArb.cpp index 7244fd721..0885b0395 100644 --- a/src/toArb.cpp +++ b/src/toArb.cpp @@ -88,9 +88,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(x) + .5) * spacing; - const double y_coord = (static_cast(y) + .5) * spacing; - const double z_coord = (static_cast(z) + .5) * spacing; + const double x_coord = (static_cast(x) + .5) ; + const double y_coord = (static_cast(y) + .5) ; + const double z_coord = (static_cast(z) + .5) ; file << x_coord << ' ' << y_coord << ' ' << z_coord << ' '; // Neighbors From 113e9adacdb07d441012a22b0e5d29760feae9b8 Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Wed, 4 Mar 2026 16:53:31 +1000 Subject: [PATCH 3/6] Add subgrid particles to the cumulant model --- models/flow/d3q27_cumulant_part/Dynamics.c.Rt | 88 ++++++++++++++----- 1 file changed, 66 insertions(+), 22 deletions(-) diff --git a/models/flow/d3q27_cumulant_part/Dynamics.c.Rt b/models/flow/d3q27_cumulant_part/Dynamics.c.Rt index c77d11e03..7dd3d01f1 100644 --- a/models/flow/d3q27_cumulant_part/Dynamics.c.Rt +++ b/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; - - 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); + + + // 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; + } + } + } + } } From cf47d890703d198a29e2a9bb7ed532d103d0bc6a Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Wed, 4 Mar 2026 17:17:36 +1000 Subject: [PATCH 4/6] Add destructor to arblattice --- src/ArbLattice.cpp | 7 ++++++- src/ArbLattice.hpp | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/ArbLattice.cpp b/src/ArbLattice.cpp index bf9575899..cd5473805 100644 --- a/src/ArbLattice.cpp +++ b/src/ArbLattice.cpp @@ -843,4 +843,9 @@ void ArbLattice::resetAverage(){ CudaMemset(&getSnapPtr(Snap)[f.id*sizes.snaps_pitch], 0, sizes.snaps_pitch*sizeof(real_t)); } } -} \ No newline at end of file +} + +ArbLattice::~ArbLattice() +{ + RFI.Close(); +} diff --git a/src/ArbLattice.hpp b/src/ArbLattice.hpp index 62658cb9c..58ec0a5b5 100644 --- a/src/ArbLattice.hpp +++ b/src/ArbLattice.hpp @@ -77,7 +77,7 @@ 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& 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; } From c2be9f85f11bd251470497d77685b238a1430fb2 Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Wed, 11 Mar 2026 16:13:20 +1000 Subject: [PATCH 5/6] Mark override --- src/ArbLattice.hpp | 2 +- src/CartLattice.h.Rt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ArbLattice.hpp b/src/ArbLattice.hpp index 58ec0a5b5..e02fbce00 100644 --- a/src/ArbLattice.hpp +++ b/src/ArbLattice.hpp @@ -101,7 +101,7 @@ class ArbLattice : public LatticeBase { const std::vector& getLocalPermutation() const { return local_permutation; } void resetAverage(); - lbRegion getLocalBoundingBox() const; /// Compute local bounding box, assuming the arbitrary lattice is a subset of a Cartesian lattice + 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 diff --git a/src/CartLattice.h.Rt b/src/CartLattice.h.Rt index 783fc6c05..5572a17e3 100644 --- a/src/CartLattice.h.Rt +++ b/src/CartLattice.h.Rt @@ -86,7 +86,7 @@ public: const lbRegion& getLocalRegion() const { return connectivity.getLocalRegion(); } const lbRegion& getGlobalRegion() const { return connectivity.global_region; } - lbRegion getLocalBoundingBox() const { + lbRegion getLocalBoundingBox() const override { lbRegion reg = getLocalRegion(); return lbRegion(px + reg.dx, py + reg.dy, pz + reg.dz, reg.nx, reg.ny, reg.nz); From 88597408d2b04d0edae71bc7f40f2834bce64ec4 Mon Sep 17 00:00:00 2001 From: Dmytro Sashko Date: Thu, 12 Mar 2026 16:03:03 +1000 Subject: [PATCH 6/6] Fixup overrides of virtual functions Were not triggered previously because destructor was virtual. Now more consistent with CartLattice --- src/ArbLattice.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ArbLattice.hpp b/src/ArbLattice.hpp index a2fea87be..5349db6ca 100644 --- a/src/ArbLattice.hpp +++ b/src/ArbLattice.hpp @@ -86,16 +86,16 @@ class ArbLattice : public LatticeBase { size_t getGlobalSize() const final { return connect.num_nodes_global; } - virtual std::vector shape() const { return {static_cast(getLocalSize())}; }; - virtual std::vector getQuantity(const Model::Quantity& q, real_t scale = 1) ; - virtual std::vector getFlags() const; - virtual std::vector getField(const Model::Field& f); - virtual std::vector getFieldAdj(const Model::Field& f); - virtual std::vector getCoord(const Model::Coord& q, real_t scale = 1); - - virtual void setFlags(const std::vector& x); - virtual void setField(const Model::Field& f, const std::vector& x); - virtual void setFieldAdjZero(const Model::Field& f); + virtual std::vector shape() const override { return {static_cast(getLocalSize())}; }; + virtual std::vector getQuantity(const Model::Quantity& q, real_t scale = 1) override ; + virtual std::vector getFlags() const override; + virtual std::vector getField(const Model::Field& f) override; + virtual std::vector getFieldAdj(const Model::Field& f) override; + virtual std::vector getCoord(const Model::Coord& q, real_t scale = 1) override; + + virtual void setFlags(const std::vector& x) override; + virtual void setField(const Model::Field& f, const std::vector& x) override; + virtual void setFieldAdjZero(const Model::Field& f) override; const ArbVTUGeom& getVTUGeom() const { return vtu_geom; } Span getNodeTypes() const { return {node_types_host.data(), node_types_host.size()}; } /// Get host view of node types (permuted)