Skip to content

Commit

Permalink
Fix initialization bug in PSE (#40)
Browse files Browse the repository at this point in the history 
* Fix initialization bug in PSE

* Add stochastic tests to PSE

* Enable full warnings when compiling examples in Debug mode
  • Loading branch information
@RaulPPelaez
RaulPPelaez authored Nov 14, 2024
1 parent 0ca7ce7 commit e1b43bc
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 47 deletions.
3 changes: 1 addition & 2 deletions examples/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@ set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CUDA_STANDARD 14)
set(CMAKE_CUDA_STANDARD_REQUIRED ON)
set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF)

#add_compile_definitions(PUBLIC DOUBLE_PRECISION)
set(UAMMD_INCLUDE ../src ../src/third_party)
set(BLA_STATIC OFF)
Expand Down Expand Up @@ -70,7 +69,7 @@ add_subdirectory(uammd_as_a_library)

IF (CMAKE_BUILD_TYPE MATCHES "Debug")
if(CMAKE_CUDA_COMPILER_ID STREQUAL "NVIDIA")
set(CMAKE_CUDA_FLAGS "-g -G -lineinfo -src-in-ptx")
set(CMAKE_CUDA_FLAGS "-g -G -lineinfo -src-in-ptx -Wall -Wextra -pedantic")
else()
set(CMAKE_CUDA_FLAGS "-g")
endif()
Expand Down
2 changes: 1 addition & 1 deletion src/Integrator/BDHI/PSE/initialization.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ namespace uammd{
}

PSE::PSE(shared_ptr<ParticleGroup> pg, Parameters par):
pg(pg), hydrodynamicRadius(par.hydrodynamicRadius){
pg(pg), hydrodynamicRadius(par.hydrodynamicRadius), temperature(par.temperature), dt(par.dt){
System::log<System::MESSAGE>("[BDHI::PSE] Initialized");
this->M0 = pse_ns::computeSelfMobility(par);
System::log<System::MESSAGE>("[BDHI::PSE] Self mobility: %f", M0);
Expand Down
186 changes: 142 additions & 44 deletions test/BDHI/PSE/pse_test.cu
View file
Original file line number Diff line number Diff line change
@@ -1,31 +1,31 @@
/*Raul P. Pelaez 2022. Tests for the FCM algorithm and integrator
*/
#include <fstream>
#include <gtest/gtest.h>
#include "gmock/gmock.h"
#include "Integrator/BDHI/BDHI_PSE.cuh"
#include "utils/container.h"
#include "gmock/gmock.h"
#include <fstream>
#include <gtest/gtest.h>
#include <memory>
#include<random>
#include<thrust/iterator/constant_iterator.h>
#include<thrust/host_vector.h>
#include <random>
#include <thrust/host_vector.h>
#include <thrust/iterator/constant_iterator.h>

using namespace uammd;
using BDHI::PSE;
//Different kernels can achieve different maximum accuracies. For
//instance, the Gaussian kernel can safely achieve 8 digits of
//accuracy in the self mobility, but the Peskin 3pt kernel will only
//give about 3. Place the maximum expected accuracy in this variable
//if you want to check a new kernel.
// Different kernels can achieve different maximum accuracies. For
// instance, the Gaussian kernel can safely achieve 8 digits of
// accuracy in the self mobility, but the Peskin 3pt kernel will only
// give about 3. Place the maximum expected accuracy in this variable
// if you want to check a new kernel.
constexpr real expectedAccuracy = 1e-8;
//Although the Gaussian kernel should be able to achieve even more
//accuracy, beyond 8 digits the approximate solution for the periodic
//correction of the self mobility I am using starts to fail.
// Although the Gaussian kernel should be able to achieve even more
// accuracy, beyond 8 digits the approximate solution for the periodic
// correction of the self mobility I am using starts to fail.

template<class T> using cached_vector = uninitialized_cached_vector<T>;
template <class T> using cached_vector = uninitialized_cached_vector<T>;

TEST(PSE, CanBeCreated){
TEST(PSE, CanBeCreated) {
PSE::Parameters par;
real hydrodynamicRadius = 1;
real3 L = {128, 128, 128};
Expand All @@ -35,33 +35,39 @@ TEST(PSE, CanBeCreated){
par.dt = 1;
par.box = Box(L);
par.psi = 1.0;
par.hydrodynamicRadius = hydrodynamicRadius;
auto pd = std::make_shared<ParticleData>(1);
auto pse = std::make_shared<PSE>(pd, par);
}

real selfMobility(real hydrodynamicRadius, real viscosity, real L){
//O(a^8) accuracy. See Hashimoto 1959.
//With a Gaussian this expression has a minimum deviation from measuraments of 7e-7*rh at L=64*rh.
//The translational invariance of the hydrodynamic radius however decreases arbitrarily with the tolerance.
//Seems that this deviation decreases with L, so probably is due to the correction below missing something.
real selfMobility(real hydrodynamicRadius, real viscosity, real L) {
// O(a^8) accuracy. See Hashimoto 1959.
// With a Gaussian this expression has a minimum deviation from measuraments
// of 7e-7*rh at L=64*rh. The translational invariance of the hydrodynamic
// radius however decreases arbitrarily with the tolerance. Seems that this
// deviation decreases with L, so probably is due to the correction below
// missing something.
long double rh = hydrodynamicRadius;
long double a = rh/L;
long double a2= a*a; long double a3 = a2*a;
long double a = rh / L;
long double a2 = a * a;
long double a3 = a2 * a;
long double c = 2.83729747948061947666591710460773907l;
long double b = 0.19457l;
long double a6pref = 16.0l*M_PIl*M_PIl/45.0l + 630.0L*b*b;
return 1.0l/(6.0l*M_PIl*viscosity*rh)*(1.0l-c*a+(4.0l/3.0l)*M_PIl*a3-a6pref*a3*a3);
long double a6pref = 16.0l * M_PIl * M_PIl / 45.0l + 630.0L * b * b;
return 1.0l / (6.0l * M_PIl * viscosity * rh) *
(1.0l - c * a + (4.0l / 3.0l) * M_PIl * a3 - a6pref * a3 * a3);
}

//Check that the self mobility stays below tolerance at a series of random points inside the domain
// in every direction when pulling a particle
TEST(PSE, SelfMobilityIsCorrectUpToTolerance){
// Check that the self mobility stays below tolerance at a series of random
// points inside the domain
// in every direction when pulling a particle
TEST(PSE, SelfMobilityIsCorrectUpToTolerance) {
PSE::Parameters par;
real hydrodynamicRadius = 1.012312;
par.viscosity = 1.12321;
par.tolerance = expectedAccuracy;
par.dt = 1;
real3 L = make_real3(128*hydrodynamicRadius);
real3 L = make_real3(128 * hydrodynamicRadius);
par.box = Box(L);
par.hydrodynamicRadius = hydrodynamicRadius;
par.psi = 1.0;
Expand All @@ -75,33 +81,125 @@ TEST(PSE, SelfMobilityIsCorrectUpToTolerance){
real3 m0;
int ntest = 20;
Saru rng(1234);
for(int j = 0; j<ntest; j++){
for (int j = 0; j < ntest; j++) {
{
auto pos = pd->getPos(access::cpu, access::write);
real3 randomPos = make_real3(rng.f(-0.5, 0.5), rng.f(-0.5, 0.5), rng.f(-0.5, 0.5))*L;
real3 randomPos =
make_real3(rng.f(-0.5, 0.5), rng.f(-0.5, 0.5), rng.f(-0.5, 0.5)) * L;
pos[0] = make_real4(randomPos);
}
for(int i = 0; i<3; i++){
switch(i){
for (int i = 0; i < 3; i++) {
switch (i) {
case 0:
force[0] = make_real4(1,0,0,0);
m0 = make_real3(selfMobility(hydrodynamicRadius, par.viscosity, L.x),0,0);
break;
force[0] = make_real4(1, 0, 0, 0);
m0 = make_real3(selfMobility(hydrodynamicRadius, par.viscosity, L.x), 0,
0);
break;
case 1:
force[0] = make_real4(0,1,0,0);
m0 = make_real3(0, selfMobility(hydrodynamicRadius, par.viscosity, L.x),0);
break;
force[0] = make_real4(0, 1, 0, 0);
m0 = make_real3(0, selfMobility(hydrodynamicRadius, par.viscosity, L.x),
0);
break;
case 2:
force[0] = make_real4(0,0,1,0);
m0 = make_real3(0, 0, selfMobility(hydrodynamicRadius, par.viscosity, L.x));
break;
force[0] = make_real4(0, 0, 1, 0);
m0 = make_real3(0, 0,
selfMobility(hydrodynamicRadius, par.viscosity, L.x));
break;
}
pse->computeHydrodynamicDisplacements(force.data().get(),
MF.data().get(), temperature, prefactor, 0);
pse->computeHydrodynamicDisplacements(force.data().get(), MF.data().get(),
temperature, prefactor, 0);
real3 dx = MF[0];
ASSERT_THAT(dx.x, ::testing::DoubleNear(m0.x, par.tolerance));
ASSERT_THAT(dx.y, ::testing::DoubleNear(m0.y, par.tolerance));
ASSERT_THAT(dx.z, ::testing::DoubleNear(m0.z, par.tolerance));
}
}
}

// Self diffusion is D=kT*M0
// <dx^2> = 2*D*dt
TEST(PSE, SelfDiffusionIsCorrectUpToToleranceHydroDisp) {
PSE::Parameters par;
real hydrodynamicRadius = 1.012312;
par.viscosity = 1.12321;
par.tolerance = 1e-4;
par.dt = 121312312231; // Should not be used
real3 L = make_real3(32 * hydrodynamicRadius);
par.box = Box(L);
par.hydrodynamicRadius = hydrodynamicRadius;
par.psi = 1.0;
par.temperature = 12312312; // Should not be used
int numberParticles = 1;
auto pd = std::make_shared<ParticleData>(numberParticles);
auto pse = std::make_shared<PSE>(pd, par);
thrust::device_vector<real3> BdW(numberParticles);
thrust::fill(thrust::cuda::par, BdW.begin(), BdW.end(), make_real3(0));
real temperature = 1;
real prefactor = 1;
real m0 = selfMobility(hydrodynamicRadius, par.viscosity, L.x);
int ntest = 1000;
real3 dx2 = make_real3(0);
Saru rng(1234);
for (int j = 0; j < ntest; j++) {
{
auto pos = pd->getPos(access::cpu, access::write);
real3 randomPos =
make_real3(rng.f(-0.5, 0.5), rng.f(-0.5, 0.5), rng.f(-0.5, 0.5)) * L;
pos[0] = make_real4(randomPos);
}
pse->computeHydrodynamicDisplacements(nullptr, BdW.data().get(),
temperature, prefactor, 0);
real3 dx = BdW[0];
dx2 += dx * dx;
}
auto diffusion = dx2 / (ntest);
ASSERT_THAT(diffusion.x, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
ASSERT_THAT(diffusion.y, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
ASSERT_THAT(diffusion.z, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
}

TEST(PSE, SelfDiffusionIsCorrectUpToToleranceComputeMF) {
PSE::Parameters par;
real hydrodynamicRadius = 1.012312;
real temperature = 1.12312;
par.viscosity = 1.12321;
par.tolerance = 1e-4;
par.dt = 1;
real3 L = make_real3(32 * hydrodynamicRadius);
par.box = Box(L);
par.hydrodynamicRadius = hydrodynamicRadius;
par.psi = 1.0;
par.temperature = temperature;
int numberParticles = 1;
auto pd = std::make_shared<ParticleData>(numberParticles);
auto pse = std::make_shared<PSE>(pd, par);
thrust::device_vector<real3> MF(numberParticles);
thrust::device_vector<real3> BdW(numberParticles);
thrust::fill(thrust::cuda::par, BdW.begin(), BdW.end(), make_real3(0));
thrust::fill(thrust::cuda::par, MF.begin(), MF.end(), make_real3(0));
real m0 = selfMobility(hydrodynamicRadius, par.viscosity, L.x);
int ntest = 1000;
real3 dx2 = make_real3(0);
Saru rng(1234);
for (int j = 0; j < ntest; j++) {
{
auto pos = pd->getPos(access::cpu, access::write);
auto forces = pd->getForce(access::cpu, access::write);
real3 randomPos =
make_real3(rng.f(-0.5, 0.5), rng.f(-0.5, 0.5), rng.f(-0.5, 0.5)) * L;
pos[0] = make_real4(randomPos);
forces[0] = make_real4(0, 0, 0, 0);
}
pse->computeMF(MF.data().get(), 0); // Far field noise is computed here
pse->computeBdW(BdW.data().get(), 0);
// Sum MF to BdW
thrust::transform(thrust::cuda::par, BdW.begin(), BdW.end(), MF.begin(),
BdW.begin(), thrust::plus<real3>());
real3 dx = BdW[0];
dx2 += dx * dx;
}
auto diffusion = dx2 / (ntest*par.dt);
ASSERT_THAT(diffusion.x, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
ASSERT_THAT(diffusion.y, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
ASSERT_THAT(diffusion.z, ::testing::DoubleNear(2.0*temperature*m0, 1e-2));
}

0 comments on commit e1b43bc

Please sign in to comment.