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; + } + } + } + } } diff --git a/src/ArbLattice.cpp b/src/ArbLattice.cpp index a24f1dda0..e2521888f 100644 --- a/src/ArbLattice.cpp +++ b/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::max() && retval >= std::numeric_limits::min()); return static_cast(retval); } @@ -895,3 +895,8 @@ void ArbLattice::resetAverage(){ } } } + +ArbLattice::~ArbLattice() +{ + RFI.Close(); +} diff --git a/src/ArbLattice.hpp b/src/ArbLattice.hpp index 9857dc23b..5349db6ca 100644 --- a/src/ArbLattice.hpp +++ b/src/ArbLattice.hpp @@ -79,23 +79,23 @@ 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; } 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 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); - virtual void setField(const Model::Field& f, const std::vector& x); - virtual void setFieldAdjZero(const Model::Field& f); + 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) @@ -103,6 +103,7 @@ class ArbLattice : public LatticeBase { const std::vector& 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 diff --git a/src/CartLattice.h.Rt b/src/CartLattice.h.Rt index bebc6f05f..5572a17e3 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 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(); } 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 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 12384cff4..7621b4b4d 100644 --- a/src/toArb.cpp +++ b/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(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