diff --git a/examples/BinaryLennardJones/SPARTIAN.cpp b/examples/BinaryLennardJones/SPARTIAN.cpp index 10a2f54c..ecaae34b 100644 --- a/examples/BinaryLennardJones/SPARTIAN.cpp +++ b/examples/BinaryLennardJones/SPARTIAN.cpp @@ -209,6 +209,7 @@ void spartian(YAML::Node& config, thermodynamicForce.getDensityProfile().min, thermodynamicForce.getDensityProfile().max, thermodynamicForce.getDensityProfile().numBins, + thermodynamicForce.getDensityProfile().binSize, COORD_X); densityProfile.scale( 1_r / (densityProfile.binSize * subdomain.diameter[1] * subdomain.diameter[2])); diff --git a/examples/DensityFluctuations/DensityFluctuations.cpp b/examples/DensityFluctuations/DensityFluctuations.cpp index e108d0c2..f9be26b1 100644 --- a/examples/DensityFluctuations/DensityFluctuations.cpp +++ b/examples/DensityFluctuations/DensityFluctuations.cpp @@ -247,6 +247,7 @@ void LJ(Config& config) subdomain.minCorner[0], subdomain.maxCorner[0], 10, + config.densityBinWidth, COORD_X); densityProfile.scale(1_r / (densityProfile.binSize * config.Ly * config.Lz * 4)); auto fluctuation = analysis::getFluctuation(densityProfile, rho, 0); diff --git a/mrmd/action/ThermodynamicForce.cpp b/mrmd/action/ThermodynamicForce.cpp index 0e23a563..cb099367 100644 --- a/mrmd/action/ThermodynamicForce.cpp +++ b/mrmd/action/ThermodynamicForce.cpp @@ -21,55 +21,121 @@ namespace mrmd { namespace action { -ThermodynamicForce::ThermodynamicForce(const std::vector& targetDensity, + +ThermodynamicForce::ThermodynamicForce(const std::vector& targetDensities, const data::Subdomain& subdomain, + const real_t& requestedDensityGridSpacing, const real_t& requestedDensityBinWidth, - const std::vector& thermodynamicForceModulation, + const std::vector& thermodynamicForceModulations, const bool enforceSymmetry, const bool usePeriodicity) : force_("thermodynamic-force", subdomain.minCorner[0], subdomain.maxCorner[0], - idx_c(std::ceil(subdomain.diameter[0] / requestedDensityBinWidth)), - idx_c(targetDensity.size())), - densityProfile_("density-profile", force_), + idx_c(std::ceil(subdomain.diameter[0] / requestedDensityGridSpacing)), + idx_c(targetDensities.size())), + densityProfile_("density-profile", + subdomain.minCorner[0], + subdomain.maxCorner[0], + idx_c(std::ceil(subdomain.diameter[0] / requestedDensityGridSpacing)), + idx_c(targetDensities.size())), binVolume_(subdomain.diameter[1] * subdomain.diameter[2] * densityProfile_.binSize), - targetDensity_(targetDensity), - thermodynamicForceModulation_(thermodynamicForceModulation), - forceFactor_("force-factor", targetDensity.size()), + targetDensities_(targetDensities), + densityBinWidth_(requestedDensityBinWidth), + thermodynamicForceModulations_(thermodynamicForceModulations), + forceFactor_("force-factor", targetDensities.size()), enforceSymmetry_(enforceSymmetry), usePeriodicity_(usePeriodicity) { - MRMD_HOST_CHECK_LESSEQUAL( - force_.binSize, requestedDensityBinWidth, "requested bin size is not achieved"); + MRMD_HOST_CHECK_LESSEQUAL(force_.binSize, + requestedDensityGridSpacing, + "requested thermodynamic force bin size is not achieved"); + MRMD_HOST_CHECK_LESSEQUAL(densityProfile_.binSize, + requestedDensityGridSpacing, + "requested density profile bin size is not achieved"); - MRMD_HOST_CHECK_EQUAL(targetDensity.size(), thermodynamicForceModulation.size()); - numTypes_ = idx_c(targetDensity.size()); + MRMD_HOST_CHECK_EQUAL(targetDensities.size(), thermodynamicForceModulations.size()); + numTypes_ = idx_c(targetDensities.size()); MRMD_HOST_CHECK_GREATER(numTypes_, 0); auto hForceFactor = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), forceFactor_); for (auto i = 0; i < numTypes_; ++i) { - hForceFactor(i) = thermodynamicForceModulation_[i] / targetDensity_[i]; + hForceFactor(i) = thermodynamicForceModulations_[i] / targetDensities_[i]; } Kokkos::deep_copy(forceFactor_, hForceFactor); + } -ThermodynamicForce::ThermodynamicForce(const real_t targetDensity, +ThermodynamicForce::ThermodynamicForce(const std::vector& targetDensities, const data::Subdomain& subdomain, - const real_t& requestedDensityBinWidth, - const real_t thermodynamicForceModulation, + const real_t& requestedDensityGridSpacing, + const std::vector& thermodynamicForceModulations, + const bool enforceSymmetry, + const bool usePeriodicity) + : ThermodynamicForce(targetDensities, + subdomain, + requestedDensityGridSpacing, + requestedDensityGridSpacing, + thermodynamicForceModulations, + enforceSymmetry, + usePeriodicity) +{ +} + +ThermodynamicForce::ThermodynamicForce(const real_t& targetDensity, + const data::Subdomain& subdomain, + const real_t& requestedDensityGridSpacing, + const real_t& thermodynamicForceModulation, const bool enforceSymmetry, const bool usePeriodicity) - : ThermodynamicForce(std::vector{targetDensity}, + : ThermodynamicForce({targetDensity}, subdomain, - requestedDensityBinWidth, + requestedDensityGridSpacing, + requestedDensityGridSpacing, {thermodynamicForceModulation}, enforceSymmetry, usePeriodicity) { } +ThermodynamicForce::ThermodynamicForce(const std::vector& targetDensities, + const data::Subdomain& subdomain, + const idx_t& requestedDensityBinNumber, + const real_t& requestedDensityBinWidth, + const std::vector& thermodynamicForceModulations, + const bool enforceSymmetry, + const bool usePeriodicity) + : force_("thermodynamic-force", + subdomain.minCorner[0], + subdomain.maxCorner[0], + requestedDensityBinNumber, + idx_c(targetDensities.size())), + densityProfile_("density-profile", + subdomain.minCorner[0], + subdomain.maxCorner[0], + requestedDensityBinNumber, + idx_c(targetDensities.size())), + binVolume_(subdomain.diameter[1] * subdomain.diameter[2] * densityProfile_.binSize), + targetDensities_(targetDensities), + densityBinWidth_(requestedDensityBinWidth), + thermodynamicForceModulations_(thermodynamicForceModulations), + forceFactor_("force-factor", targetDensities.size()), + enforceSymmetry_(enforceSymmetry), + usePeriodicity_(usePeriodicity) +{ + MRMD_HOST_CHECK_EQUAL(targetDensities.size(), thermodynamicForceModulations.size()); + numTypes_ = idx_c(targetDensities.size()); + MRMD_HOST_CHECK_GREATER(numTypes_, 0); + + auto hForceFactor = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), forceFactor_); + for (auto i = 0; i < numTypes_; ++i) + { + hForceFactor(i) = thermodynamicForceModulations_[i] / targetDensities_[i]; + } + Kokkos::deep_copy(forceFactor_, hForceFactor); +} + void ThermodynamicForce::sample(data::Atoms& atoms) { densityProfile_ += analysis::getAxialDensityProfile(atoms.numLocalAtoms, @@ -79,6 +145,7 @@ void ThermodynamicForce::sample(data::Atoms& atoms) densityProfile_.min, densityProfile_.max, densityProfile_.numBins, + densityBinWidth_, COORD_X); ++densityProfileSamples_; @@ -101,7 +168,36 @@ void ThermodynamicForce::update(const real_t& smoothingSigma, const real_t& smoo auto smoothedDensityGradient = data::gradient(smoothedDensityProfile, usePeriodicity_); smoothedDensityGradient.scale(forceFactor_); - force_ -= smoothedDensityGradient; + force_ -= util::interpolate(smoothedDensityGradient, force_.createGrid_d()); + + // reset sampling data + Kokkos::deep_copy(densityProfile_.data, 0_r); + densityProfileSamples_ = 0; +} + +void ThermodynamicForce::update(const real_t& smoothingSigma, + const real_t& smoothingIntensity, + const util::ApplicationRegion& applicationRegion) +{ + MRMD_HOST_CHECK_GREATER(densityProfileSamples_, 0); + + if (enforceSymmetry_) + { + densityProfile_.makeSymmetric(); + } + + auto normalizationFactor = 1_r / (binVolume_ * real_c(densityProfileSamples_)); + densityProfile_.scale(normalizationFactor); + + auto smoothedDensityProfile = + data::smoothen(densityProfile_, smoothingSigma, smoothingIntensity, usePeriodicity_); + auto smoothedDensityGradient = data::gradient(smoothedDensityProfile, usePeriodicity_); + smoothedDensityGradient.scale(forceFactor_); + + auto constrainedDensityGradient = + util::constrainToApplicationRegion(smoothedDensityGradient, applicationRegion); + + force_ -= util::interpolate(constrainedDensityGradient, force_.createGrid_d()); // reset sampling data Kokkos::deep_copy(densityProfile_.data, 0_r); diff --git a/mrmd/action/ThermodynamicForce.hpp b/mrmd/action/ThermodynamicForce.hpp index 5740d102..7dee8640 100644 --- a/mrmd/action/ThermodynamicForce.hpp +++ b/mrmd/action/ThermodynamicForce.hpp @@ -20,6 +20,7 @@ #include "data/Subdomain.hpp" #include "datatypes.hpp" #include "util/ApplicationRegion.hpp" +#include "util/interpolation.hpp" #include "weighting_function/Slab.hpp" namespace mrmd @@ -33,8 +34,9 @@ class ThermodynamicForce data::MultiHistogram densityProfile_; idx_t densityProfileSamples_ = 0; real_t binVolume_; - const std::vector targetDensity_; - const std::vector thermodynamicForceModulation_; + const std::vector targetDensities_; + const real_t densityBinWidth_; + const std::vector thermodynamicForceModulations_; idx_t numTypes_; ScalarView forceFactor_; ///< precalculated prefactor for force calculation @@ -61,23 +63,42 @@ class ThermodynamicForce void sample(data::Atoms& atoms); void update(const real_t& smoothingSigma, const real_t& smoothingIntensity); + void update(const real_t& smoothingSigma, + const real_t& smoothingIntensity, + const util::ApplicationRegion& applicationRegion); void apply(const data::Atoms& atoms) const; void apply(const data::Atoms& atoms, const weighting_function::Slab& slab) const; void apply(const data::Atoms& atoms, const util::ApplicationRegion& applicationRegion) const; std::vector getMuLeft() const; std::vector getMuRight() const; - ThermodynamicForce(const std::vector& targetDensity, + ThermodynamicForce(const std::vector& targetDensities, const data::Subdomain& subdomain, + const real_t& requestedDensityGridSpacing, const real_t& requestedDensityBinWidth, - const std::vector& thermodynamicForceModulation, + const std::vector& thermodynamicForceModulations, const bool enforceSymmetry = false, const bool usePeriodicity = false); - ThermodynamicForce(const real_t targetDensity, + ThermodynamicForce(const std::vector& targetDensities, const data::Subdomain& subdomain, + const real_t& requestedDensityGridSpacing, + const std::vector& thermodynamicForceModulations, + const bool enforceSymmetry = false, + const bool usePeriodicity = false); + + ThermodynamicForce(const real_t& targetDensity, + const data::Subdomain& subdomain, + const real_t& requestedDensityGridSpacing, + const real_t& thermodynamicForceModulation, + const bool enforceSymmetry = false, + const bool usePeriodicity = false); + + ThermodynamicForce(const std::vector& targetDensities, + const data::Subdomain& subdomain, + const idx_t& requestedDensityBinNumber, const real_t& requestedDensityBinWidth, - const real_t thermodynamicForceModulation, + const std::vector& thermodynamicForceModulations, const bool enforceSymmetry = false, const bool usePeriodicity = false); }; diff --git a/mrmd/analysis/AxialDensityProfile.cpp b/mrmd/analysis/AxialDensityProfile.cpp index 8064b060..2adab549 100644 --- a/mrmd/analysis/AxialDensityProfile.cpp +++ b/mrmd/analysis/AxialDensityProfile.cpp @@ -20,6 +20,7 @@ namespace mrmd { namespace analysis { + data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms, const data::Atoms::pos_t& positions, const data::Atoms::type_t& types, @@ -27,6 +28,7 @@ data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms, const real_t min, const real_t max, const int64_t numBins, + const real_t binWidth, const int64_t axis) { MRMD_HOST_CHECK_GREATEREQUAL(max, min); @@ -37,15 +39,18 @@ data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms, data::MultiHistogram histogram("density-profile", min, max, numBins, numTypes); MultiScatterView scatter(histogram.data); - auto policy = Kokkos::RangePolicy<>(0, numAtoms); - auto kernel = KOKKOS_LAMBDA(const idx_t idx) + auto policy = Kokkos::MDRangePolicy>({0, 0}, {numAtoms, numBins}); + auto kernel = KOKKOS_LAMBDA(const idx_t atomIdx, const idx_t binIdx) { - MRMD_DEVICE_ASSERT_GREATEREQUAL(types(idx), 0); - MRMD_DEVICE_ASSERT_LESS(types(idx), numTypes); - auto bin = histogram.getBin(positions(idx, axis)); - if (bin == -1) return; + MRMD_DEVICE_ASSERT_GREATEREQUAL(types(atomIdx), 0); + MRMD_DEVICE_ASSERT_LESS(types(atomIdx), numTypes); + auto binPos = histogram.getBinPosition(binIdx); + auto atomPos = positions(atomIdx, axis); auto access = scatter.access(); - access(bin, types(idx)) += 1_r; + if (atomPos >= binPos - binWidth / 2 && atomPos <= binPos + binWidth / 2) + { + access(binIdx, types(atomIdx)) += 1_r; + } }; Kokkos::parallel_for(policy, kernel); Kokkos::Experimental::contribute(histogram.data, scatter); diff --git a/mrmd/analysis/AxialDensityProfile.hpp b/mrmd/analysis/AxialDensityProfile.hpp index d6c7a639..c5c9c28f 100644 --- a/mrmd/analysis/AxialDensityProfile.hpp +++ b/mrmd/analysis/AxialDensityProfile.hpp @@ -35,6 +35,7 @@ data::MultiHistogram getAxialDensityProfile(const idx_t numAtoms, const real_t min, const real_t max, const int64_t numBins, + const real_t binWidth, const int64_t axis); } // namespace analysis diff --git a/mrmd/analysis/AxialDensityProfile.test.cpp b/mrmd/analysis/AxialDensityProfile.test.cpp index cfbe9e68..c9f3b0f1 100644 --- a/mrmd/analysis/AxialDensityProfile.test.cpp +++ b/mrmd/analysis/AxialDensityProfile.test.cpp @@ -60,7 +60,7 @@ TEST(AxialDensityProfile, histogram) auto atoms = initAtoms(); auto histogram = getAxialDensityProfile( - atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, COORD_X); + atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, 1_r, COORD_X); auto h_data = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogram.data); for (auto i = 0; i < 10; ++i) @@ -69,5 +69,34 @@ TEST(AxialDensityProfile, histogram) EXPECT_FLOAT_EQ(h_data(i, 1), real_c(11 - (i + 1))); } } + +TEST(AxialDensityProfile, overlappingBins) +{ + auto atoms = initAtoms(); + + auto histogram = getAxialDensityProfile( + atoms.numLocalAtoms, atoms.getPos(), atoms.getType(), 2, 0_r, 10_r, 10, 3_r, COORD_X); + auto h_data = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogram.data); + + for (auto i = 0; i < 10; ++i) + { + if (i == 0) + { + EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1) + real_c(i + 2))); + EXPECT_FLOAT_EQ(h_data(i, 1), (real_c(11 - (i + 1)) + real_c(11 - (i + 2)))); + continue; + } + if (i == 9) + { + EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1))); + EXPECT_FLOAT_EQ(h_data(i, 1), (real_c(11 - i) + real_c(11 - (i + 1)))); + continue; + } + EXPECT_FLOAT_EQ(h_data(i, 0), (real_c(i) + real_c(i + 1) + real_c(i + 2))); + EXPECT_FLOAT_EQ(h_data(i, 1), + (real_c(11 - i) + real_c(11 - (i + 1)) + real_c(11 - (i + 2)))); + } +} + } // namespace analysis } // namespace mrmd \ No newline at end of file diff --git a/mrmd/data/MultiHistogram.cpp b/mrmd/data/MultiHistogram.cpp index 9b184746..d92dc459 100644 --- a/mrmd/data/MultiHistogram.cpp +++ b/mrmd/data/MultiHistogram.cpp @@ -33,6 +33,24 @@ ScalarView::HostMirror MultiHistogram::createGrid() const return grid; } +ScalarView MultiHistogram::createGrid_d() const +{ + MRMD_HOST_CHECK_GREATEREQUAL(numBins, 0); + + ScalarView grid("grid", numBins); + const real_t min_ = min; + const real_t binSize_ = binSize; + auto policy = Kokkos::RangePolicy<>(0, numBins); + auto kernel = KOKKOS_LAMBDA(const idx_t idx) + { + grid[idx] = min_ + (real_c(idx) + 0.5_r) * binSize_; + }; + Kokkos::parallel_for("MultiHistogram::scale", policy, kernel); + Kokkos::fence(); + + return grid; +} + MultiHistogram& MultiHistogram::operator+=(const MultiHistogram& rhs) { transform(*this, rhs, *this, bin_op::Add()); diff --git a/mrmd/data/MultiHistogram.hpp b/mrmd/data/MultiHistogram.hpp index 469ec337..fc9ea6d9 100644 --- a/mrmd/data/MultiHistogram.hpp +++ b/mrmd/data/MultiHistogram.hpp @@ -31,12 +31,13 @@ struct MultiHistogram const real_t minArg, const real_t maxArg, idx_t numBinsArg, + real_t binSizeArg, idx_t numHistogramsArg) : min(minArg), max(maxArg), numBins(numBinsArg), numHistograms(numHistogramsArg), - binSize((max - min) / real_c(numBins)), + binSize(binSizeArg), inverseBinSize(1_r / binSize), data(label, numBins, numHistograms) { @@ -45,6 +46,20 @@ struct MultiHistogram MRMD_HOST_CHECK_GREATEREQUAL(numHistogramsArg, 0); } + MultiHistogram(const std::string& label, + const real_t minArg, + const real_t maxArg, + idx_t numBinsArg, + idx_t numHistogramsArg) + : MultiHistogram(label, + minArg, + maxArg, + numBinsArg, + (maxArg - minArg) / real_c(numBinsArg), + numHistogramsArg) + { + } + MultiHistogram(const std::string& label, const MultiHistogram& histogram) : MultiHistogram( label, histogram.min, histogram.max, histogram.numBins, histogram.numHistograms) @@ -73,6 +88,7 @@ struct MultiHistogram } ScalarView::HostMirror createGrid() const; + ScalarView createGrid_d() const; const real_t min; const real_t max; diff --git a/mrmd/data/MultiHistogram.test.cpp b/mrmd/data/MultiHistogram.test.cpp index 240c7b39..a788a369 100644 --- a/mrmd/data/MultiHistogram.test.cpp +++ b/mrmd/data/MultiHistogram.test.cpp @@ -37,6 +37,16 @@ TEST(MultiHistogram, getBinPosition) EXPECT_FLOAT_EQ(histogram.getBinPosition(5), 5.5_r); } +TEST(MultiHistogram, consistencyBinToPosition) +{ + MultiHistogram histogram("histogram", 0_r, 10_r, 10, 2); + + for (idx_t idx = 0; idx < histogram.numBins; ++idx) + { + EXPECT_FLOAT_EQ(histogram.getBin(histogram.getBinPosition(idx)), idx); + } +} + TEST(MultiHistogram, scale) { MultiHistogram histogram("histogram", 0_r, 10_r, 11, 2); diff --git a/mrmd/io/RestoreThermoForce.cpp b/mrmd/io/RestoreThermoForce.cpp index 4cda5f3d..0128d533 100644 --- a/mrmd/io/RestoreThermoForce.cpp +++ b/mrmd/io/RestoreThermoForce.cpp @@ -27,43 +27,33 @@ action::ThermodynamicForce restoreThermoForce( const std::vector& thermodynamicForceModulations, const bool enforceSymmetry, const bool usePeriodicity, - const idx_t maxNumForces) + const idx_t maxNumForces, + const real_t requestedDensityBinWidth) { std::string line; std::string word; - int binNum = 0; + int binNumForce = 0; int histNum = 0; - real_t grid0 = 0; - real_t grid1 = 0; std::ifstream fileThermoForce(filename); std::getline(fileThermoForce, line); std::stringstream gridLineStream(line); while (gridLineStream >> word) { - if (binNum == 0) - { - grid0 = std::stod(word); - } - if (binNum == 1) - { - grid1 = std::stod(word); - } - binNum++; + binNumForce++; } - MRMD_HOST_ASSERT_GREATER(binNum, 1); - real_t binWidth = grid1 - grid0; + MRMD_HOST_ASSERT_GREATER(binNumForce, 1); - MultiView::HostMirror h_forcesRead("h_forcesRead", binNum, maxNumForces); + MultiView::HostMirror h_forcesRead("h_forcesRead", binNumForce, maxNumForces); while (std::getline(fileThermoForce, line)) { - binNum = 0; + binNumForce = 0; std::stringstream forceLineStream(line); while (forceLineStream >> word) { - h_forcesRead(binNum, histNum) = std::stod(word); - binNum++; + h_forcesRead(binNumForce, histNum) = std::stod(word); + binNumForce++; } histNum++; @@ -71,14 +61,23 @@ action::ThermodynamicForce restoreThermoForce( } fileThermoForce.close(); - auto h_forces = - Kokkos::subview(h_forcesRead, Kokkos::make_pair(0, binNum), Kokkos::make_pair(0, histNum)); - MultiView d_forces("d_forces", binNum, histNum); + auto h_forces = Kokkos::subview( + h_forcesRead, Kokkos::make_pair(0, binNumForce), Kokkos::make_pair(0, histNum)); + MultiView d_forces("d_forces", binNumForce, histNum); Kokkos::deep_copy(d_forces, h_forces); + auto forceBinNumber = idx_c(binNumForce); + real_t densityBinWidth = requestedDensityBinWidth; + + if (requestedDensityBinWidth < 0_r) + { + densityBinWidth = (subdomain.diameter[0] / real_c(forceBinNumber)); + } + action::ThermodynamicForce thermodynamicForce(targetDensities, subdomain, - binWidth, + forceBinNumber, + densityBinWidth, thermodynamicForceModulations, enforceSymmetry, usePeriodicity); diff --git a/mrmd/io/RestoreThermoForce.hpp b/mrmd/io/RestoreThermoForce.hpp index 4d40babb..979f2640 100644 --- a/mrmd/io/RestoreThermoForce.hpp +++ b/mrmd/io/RestoreThermoForce.hpp @@ -27,6 +27,7 @@ action::ThermodynamicForce restoreThermoForce( const std::vector& thermodynamicForceModulations = {1_r}, const bool enforceSymmetry = false, const bool usePeriodicity = false, - const idx_t maxNumForces = 10); + const idx_t maxNumForces = 10, + const real_t requestedDensityBinWidth = -1_r); } // namespace io } // namespace mrmd \ No newline at end of file diff --git a/mrmd/io/ThermoForce.test.cpp b/mrmd/io/ThermoForce.test.cpp index 0c4fd62d..31bb682f 100644 --- a/mrmd/io/ThermoForce.test.cpp +++ b/mrmd/io/ThermoForce.test.cpp @@ -28,9 +28,9 @@ action::ThermodynamicForce createThermoForce(const idx_t& numBins, const std::vector& targetDensities, const std::vector& forceModulations) { - real_t binWidth = (subdomain.maxCorner[0] - subdomain.minCorner[0]) / real_c(numBins); + real_t gridSpacing = (subdomain.maxCorner[0] - subdomain.minCorner[0]) / real_c(numBins); action::ThermodynamicForce thermodynamicForce( - targetDensities, subdomain, binWidth, forceModulations); + targetDensities, subdomain, gridSpacing, forceModulations); MultiView forces("thermoForce", numBins, numForces); diff --git a/mrmd/util/ApplicationRegion.hpp b/mrmd/util/ApplicationRegion.hpp index 9f46642b..22eef592 100644 --- a/mrmd/util/ApplicationRegion.hpp +++ b/mrmd/util/ApplicationRegion.hpp @@ -34,7 +34,7 @@ class ApplicationRegion auto dx = x - center_[0]; auto absDx = std::abs(dx); - return (absDx >= applicationRegionMin_ && absDx <= applicationRegionMax_); + return (absDx > applicationRegionMin_ && absDx < applicationRegionMax_); } ApplicationRegion(const Point3D& center, @@ -45,6 +45,12 @@ class ApplicationRegion applicationRegionMax_(applicationRegionMax) { } + + Point3D getCenter() const { return center_; } + + real_t getApplicationRegionMin() const { return applicationRegionMin_; } + + real_t getApplicationRegionMax() const { return applicationRegionMax_; } }; } // namespace util } // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/ApplicationRegion.test.cpp b/mrmd/util/ApplicationRegion.test.cpp index 27decbbe..3165d903 100644 --- a/mrmd/util/ApplicationRegion.test.cpp +++ b/mrmd/util/ApplicationRegion.test.cpp @@ -27,15 +27,17 @@ TEST(ApplicationRegion, testRegion) const auto applicationRegionMax = 4_r; auto applicationRegion = ApplicationRegion(center, applicationRegionMin, applicationRegionMax); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(4_r, 10_r, 10_r)); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(5.7_r, 10_r, 10_r)); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(6_r, 10_r, 10_r)); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(-2_r, 10_r, 10_r)); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(-1.2_r, 10_r, 10_r)); - EXPECT_TRUE(applicationRegion.isInApplicationRegion(0_r, 10_r, 10_r)); - + EXPECT_FALSE(applicationRegion.isInApplicationRegion(4_r, 10_r, 10_r)); + EXPECT_FALSE(applicationRegion.isInApplicationRegion(6_r, 10_r, 10_r)); + EXPECT_FALSE(applicationRegion.isInApplicationRegion(0_r, 10_r, 10_r)); + EXPECT_FALSE(applicationRegion.isInApplicationRegion(-2_r, 10_r, 10_r)); EXPECT_FALSE(applicationRegion.isInApplicationRegion(3.1_r, 10_r, 10_r)); EXPECT_FALSE(applicationRegion.isInApplicationRegion(-2.1_r, 10_r, 10_r)); + + EXPECT_TRUE(applicationRegion.isInApplicationRegion(5.7_r, 10_r, 10_r)); + EXPECT_TRUE(applicationRegion.isInApplicationRegion(-1.2_r, 10_r, 10_r)); + EXPECT_TRUE(applicationRegion.isInApplicationRegion(4.1_r, 10_r, 10_r)); + EXPECT_TRUE(applicationRegion.isInApplicationRegion(-0.1_r, 10_r, 10_r)); } } // namespace util diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp new file mode 100644 index 00000000..e7102e75 --- /dev/null +++ b/mrmd/util/interpolation.cpp @@ -0,0 +1,81 @@ +// Copyright 2024 Sebastian Eibl +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "util/interpolation.hpp" + +namespace mrmd +{ +namespace util +{ +data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid) +{ + data::MultiHistogram output( + "interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms); + const ScalarView& inputGrid = input.createGrid_d(); + + auto policy = Kokkos::MDRangePolicy>( + {idx_t(0), idx_t(0)}, {idx_c(grid.extent(0)), input.numHistograms}); + auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx) + { + // find the two enclosing bins in the input histogram + auto rightBin = findRightBin(inputGrid, grid(binIdx)); + auto leftBin = rightBin - 1; + + // handle boundaries + if (leftBin < 0 || rightBin >= input.numBins) + { + output.data(binIdx, histogramIdx) = 0.0_r; // out of bounds, set to zero + return; + } + + auto inputDataLeft = input.data(leftBin, histogramIdx); + auto inputDataRight = input.data(rightBin, histogramIdx); + + output.data(binIdx, histogramIdx) = + lerp(inputDataLeft, + inputDataRight, + (grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize); + }; + Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); + Kokkos::fence(); + + return output; +} + +data::MultiHistogram constrainToApplicationRegion(const data::MultiHistogram& input, + const util::ApplicationRegion& applicationRegion) +{ + data::MultiHistogram constrainedProfile( + "constrained-profile", input.min, input.max, input.numBins, input.numHistograms); + + auto policy = Kokkos::MDRangePolicy>({idx_t(0), idx_t(0)}, + {input.numBins, input.numHistograms}); + auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx) + { + if (!applicationRegion.isInApplicationRegion(input.getBinPosition(binIdx), 0_r, 0_r)) + { + constrainedProfile.data(binIdx, histogramIdx) = 0_r; + } + else + { + constrainedProfile.data(binIdx, histogramIdx) = input.data(binIdx, histogramIdx); + } + }; + Kokkos::parallel_for("MultiHistogram::constrainToApplicationRegion", policy, kernel); + Kokkos::fence(); + + return constrainedProfile; +} +} // namespace util +} // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp new file mode 100644 index 00000000..e9db338e --- /dev/null +++ b/mrmd/util/interpolation.hpp @@ -0,0 +1,67 @@ +// Copyright 2024 Sebastian Eibl +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "data/MultiHistogram.hpp" +#include "util/ApplicationRegion.hpp" + +namespace mrmd +{ +namespace util +{ +/** + * Linear interpolation between two values. + * @param left left value + * @param right right value + * @param factor interpolation factor in [0, 1] + * @return interpolated value + */ +KOKKOS_INLINE_FUNCTION +real_t lerp(const real_t& left, const real_t& right, const real_t& factor) +{ + return left + (right - left) * factor; +} + +/** + * Find the index of the first bin in the grid that is greater than or equal to the value. + * @param grid grid to search in + * @param value value to find + * @return index of the right bin + */ +KOKKOS_INLINE_FUNCTION +idx_t findRightBin(const ScalarView& grid, const real_t& value) +{ + idx_t idx = 0; + for (; idx < idx_c(grid.extent(0)) && grid(idx) < value; ++idx) + { + // empty loop body + } + + return idx; +} + +/** + * Interpolate data values of MultiHistogram to a new grid. + * @param input input MultiHistogram + * @param grid grid to interpolate to + * @return MultiHistogram with interpolated values + */ +data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid); + +data::MultiHistogram constrainToApplicationRegion(const data::MultiHistogram& input, + const util::ApplicationRegion& applicationRegion); + +} // namespace util +} // namespace mrmd \ No newline at end of file diff --git a/pyMRMD/action.cpp b/pyMRMD/action.cpp index 0a485f0b..72188bf6 100644 --- a/pyMRMD/action.cpp +++ b/pyMRMD/action.cpp @@ -85,7 +85,8 @@ void init_action(py::module_& m) static_cast( &action::ThermodynamicForce::getDensityProfile)) .def("sample", &action::ThermodynamicForce::sample) - .def("update", &action::ThermodynamicForce::update) + .def("update", + static_cast(&action::ThermodynamicForce::update)) .def("apply", static_cast( &action::ThermodynamicForce::apply)) diff --git a/pyMRMD/analysis.cpp b/pyMRMD/analysis.cpp index 33722b75..48970b48 100644 --- a/pyMRMD/analysis.cpp +++ b/pyMRMD/analysis.cpp @@ -37,10 +37,11 @@ void init_analysis(py::module_ &m) const real_t min, const real_t max, const int64_t numBins, + const real_t binWidth, const int64_t axis) { return analysis::getAxialDensityProfile( - numAtoms, atoms.getPos(), atoms.getType(), numTypes, min, max, numBins, axis); + numAtoms, atoms.getPos(), atoms.getType(), numTypes, min, max, numBins, binWidth, axis); }); py::class_(m, "MeanSquareDisplacement")