Skip to content

Commit

Permalink
Merge branch 'v2.x' into complex_bvp
Browse files Browse the repository at this point in the history
  • Loading branch information
@RaulPPelaez
RaulPPelaez committed Feb 4, 2023
2 parents 5ad238d + 93dc669 commit 73335f7
Show file tree
Hide file tree
Showing 10 changed files with 319 additions and 78 deletions.
172 changes: 109 additions & 63 deletions src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include <stdexcept>
#include "misc/LanczosAlgorithm.cuh"
#include"StokesSlab/spreadInterp.cuh"
#include<thrust/functional.h>

namespace uammd{
namespace DPStokesSlab_ns{
Expand Down Expand Up @@ -140,21 +141,23 @@ namespace uammd{

};

struct SaruTransform{
struct SaruTransform{
uint s1, s2;
SaruTransform(uint s1, uint s2):s1(s1), s2(s2){}
real std;
SaruTransform(uint s1, uint s2, real std = 1.0):
s1(s1), s2(s2), std(std){}

__device__ real3 operator()(uint id){
Saru rng(s1, s2, id);
return make_real3(rng.gf(0.0f,1.0f), rng.gf(0.0f,1.0f).x);
return make_real3(rng.gf(0.0f,std), rng.gf(0.0f,std).x);

}
};

auto fillRandomVectorReal3(int numberParticles, uint s1, uint s2){
auto fillRandomVectorReal3(int numberParticles, uint s1, uint s2, real std = 1.0){
cached_vector<real3> noise(numberParticles);
auto cit = thrust::make_counting_iterator<uint>(0);
auto tr = thrust::make_transform_iterator(cit, SaruTransform(s1, s2));
auto tr = thrust::make_transform_iterator(cit, SaruTransform(s1, s2, std));
thrust::copy(thrust::cuda::par,
tr, tr + numberParticles,
noise.begin());
Expand All @@ -175,32 +178,28 @@ namespace uammd{
struct BDIntegrate{
real4* pos;
real3* mf;
real3* sqrtmdwn;
real3* sqrtmdwnm1;
real3* mpw;
real3* mmw;
real rfdPrefactor;
real3* noise;
real3* noisePrev;
real3* rfd;
real dt;
real noisePrefactor;
BDIntegrate(real4* pos,
real3* mf,
real3* sqrtmdwn, real3* sqrtmdwnm1,
real3* mpw, real3* mmw, real deltaRFD,
real3* noise, real3* noisePrev,
real3* rfd,
real dt, real temperature):
pos(pos), mf(mf), sqrtmdwn(sqrtmdwn),sqrtmdwnm1(sqrtmdwnm1),
mpw(mpw), mmw(mmw), dt(dt){
this->noisePrefactor = sqrt(2*temperature*dt);
this->rfdPrefactor = dt*temperature/deltaRFD;
pos(pos), mf(mf), noise(noise),noisePrev(noisePrev),
rfd(rfd), dt(dt){
}
__device__ void operator()(int id){
real3 fluct;
if(sqrtmdwnm1)
fluct = noisePrefactor*real(0.5)*(sqrtmdwn[id] + sqrtmdwnm1[id]);
else
fluct = noisePrefactor*sqrtmdwn[id];
real3 displacement = mf[id]*dt +
fluct
+rfdPrefactor * (mpw[id] - mmw[id]);
real3 displacement = mf[id]*dt;
if(noise){
real3 fluct;
if(noisePrev)
fluct = dt*real(0.5)*(noise[id] + noisePrev[id]);
else
fluct = dt*noise[id];
displacement += fluct + rfd[id];
}
pos[id] += make_real4(displacement);
}
};
Expand All @@ -209,7 +208,7 @@ namespace uammd{

class DPStokesIntegrator: public Integrator{
int steps = 0;
uint seed;
uint seed, seedRFD;
std::shared_ptr<DPStokes> dpstokes;
std::shared_ptr<lanczos::Solver> lanczos;
thrust::device_vector<real3> previousNoise;
Expand All @@ -228,9 +227,68 @@ namespace uammd{
System::log<System::MESSAGE>("[DPStokes] dt %g", par.dt);
System::log<System::MESSAGE>("[DPStokes] temperature: %g", par.temperature);
this->seed = pd->getSystem()->rng().next32();
this->seedRFD = pd->getSystem()->rng().next32();
this->deltaRFD = 1e-6*par.hydrodynamicRadius;
}

//Returns the thermal drift term: temperature*dt*(\partial_q \cdot M)
auto computeThermalDrift(){
auto pos = pd->getPos(access::gpu, access::read);
const int numberParticles = pd->getNumParticles();
if(par.temperature){
auto noise2 = detail::fillRandomVectorReal3(numberParticles, seedRFD, steps);
auto cit = thrust::make_counting_iterator(0);
auto posp = thrust::make_transform_iterator(cit,
detail::SumPosAndNoise(pos.raw(),
noise2.data().get(),
deltaRFD*0.5));
auto mpw = dpstokes->Mdot(posp, noise2.data().get(), numberParticles, 0);
auto posm = thrust::make_transform_iterator(cit,
detail::SumPosAndNoise(pos.raw(),
noise2.data().get(),
-deltaRFD*0.5));
auto mmw = dpstokes->Mdot(posm, noise2.data().get(), numberParticles, 0);
using namespace thrust::placeholders;
thrust::transform(mpw.begin(), mpw.end(),
mmw.begin(),
mpw.begin(),
make_real3(par.dt*par.temperature/deltaRFD)*(_1 - _2));
return mpw;
}
else{
return cached_vector<real3>();
}
}

//Returns sqrt(2*M*temperature/dt)dW
auto computeFluctuations(){
auto pos = pd->getPos(access::gpu, access::read);
const int numberParticles = pd->getNumParticles();
cached_vector<real3> bdw(numberParticles);
thrust::fill(bdw.begin(), bdw.end(), real3());
if(par.temperature){
detail::LanczosAdaptor dot(dpstokes, pos.raw(), numberParticles);
auto noise = detail::fillRandomVectorReal3(numberParticles, seed, steps,
sqrt(2*par.temperature/par.dt));
lanczos->run(dot,
(real*)bdw.data().get(), (real*)noise.data().get(),
par.tolerance, 3*numberParticles);
}
return bdw;
}

//Returns the product of the forces and the mobility matrix, M F
auto computeDeterministicDisplacements(){
auto pos = pd->getPos(access::gpu, access::read);
auto force = pd->getForce(access::gpu, access::read);
if(pd->isTorqueAllocated()){
System::log<System::EXCEPTION>("[DPStokes] Torques are not yet implemented");
throw std::runtime_error("Operation not implemented");
}
const int numberParticles = pd->getNumParticles();
return dpstokes->Mdot(pos.raw(), force.raw(), numberParticles, 0);
}

void forwardTime(){
System::log<System::DEBUG2>("[DPStokes] Running step %d", steps);
if(steps==0) setUpInteractors();
Expand All @@ -239,46 +297,34 @@ namespace uammd{
i->updateSimulationTime(par.dt*steps);
i->sum({.force=true});
}
auto pos = pd->getPos(access::gpu, access::readwrite);
auto force = pd->getForce(access::gpu, access::read);
auto torque = pd->getTorqueIfAllocated(access::gpu, access::read);
const int numberParticles = pd->getNumParticles();
auto mf = dpstokes->Mdot(pos.raw(), force.raw(), numberParticles, 0);
if(torque.raw()){
System::log<System::EXCEPTION>("[DPStokes] Torques are not yet implemented");
throw std::runtime_error("Operation not implemented");
auto mf = computeDeterministicDisplacements();
if(par.temperature){
auto bdw = computeFluctuations();
if(par.useLeimkuhler and previousNoise.size() != numberParticles){
previousNoise.resize(numberParticles);
thrust::copy(bdw.begin(), bdw.end(), previousNoise.begin());
}
auto rfd = computeThermalDrift();
real3* d_prevNoise = par.useLeimkuhler?previousNoise.data().get():nullptr;
auto pos = pd->getPos(access::gpu, access::readwrite);
detail::BDIntegrate bd(pos.raw(), mf.data().get(),
bdw.data().get(), d_prevNoise,
rfd.data().get(),
par.dt, par.temperature);
auto cit = thrust::make_counting_iterator(0);
thrust::for_each_n(thrust::cuda::par, cit, numberParticles, bd);
if(par.useLeimkuhler)
thrust::copy(bdw.begin(), bdw.end(), previousNoise.begin());
}
cached_vector<real3> bdw(numberParticles);
thrust::fill(bdw.begin(), bdw.end(), real3());
detail::LanczosAdaptor dot(dpstokes, pos.raw(), numberParticles);
auto noise = detail::fillRandomVectorReal3(numberParticles, seed, 2*steps);
lanczos->run(dot,
(real*)bdw.data().get(), (real*)noise.data().get(),
par.tolerance, 3*numberParticles);
if(par.useLeimkuhler and previousNoise.size() != numberParticles){
previousNoise.resize(numberParticles);
thrust::copy(bdw.begin(), bdw.end(), previousNoise.begin());
else{
auto pos = pd->getPos(access::gpu, access::readwrite);
detail::BDIntegrate bd(pos.raw(), mf.data().get(),
nullptr, nullptr, nullptr,
par.dt, 0);
auto cit = thrust::make_counting_iterator(0);
thrust::for_each_n(thrust::cuda::par, cit, numberParticles, bd);
}
auto noise2 = detail::fillRandomVectorReal3(numberParticles, seed, 2*steps+1);
auto cit = thrust::make_counting_iterator(0);
auto posp = thrust::make_transform_iterator(cit,
detail::SumPosAndNoise(pos.raw(),
noise2.data().get(),
deltaRFD*0.5));
auto mpw = dpstokes->Mdot(posp, noise2.data().get(), numberParticles, 0);
auto posm = thrust::make_transform_iterator(cit,
detail::SumPosAndNoise(pos.raw(),
noise2.data().get(),
-deltaRFD*0.5));
auto mmw = dpstokes->Mdot(posm, noise2.data().get(), numberParticles, 0);
real3* d_prevNoise = par.useLeimkuhler?previousNoise.data().get():nullptr;
detail::BDIntegrate bd(pos.raw(), mf.data().get(),
bdw.data().get(), d_prevNoise,
mpw.data().get(), mmw.data().get(),
deltaRFD, par.dt, par.temperature);
thrust::for_each_n(thrust::cuda::par, cit, numberParticles, bd);
if(par.useLeimkuhler)
thrust::copy(bdw.begin(), bdw.end(), previousNoise.begin());
steps++;
}

Expand Down
10 changes: 5 additions & 5 deletions src/Integrator/BDHI/DoublyPeriodic/StokesSlab/utils.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,16 +202,16 @@ namespace uammd{
m_z.resize(newSize);
}

void fillWithZero() const{
void fillWithZero(){
thrust::fill(m_x.begin(), m_x.end(), T());
thrust::fill(m_y.begin(), m_y.end(), T());
thrust::fill(m_z.begin(), m_z.end(), T());
}

using Iterator = T*;
Iterator x()const{return thrust::raw_pointer_cast(m_x.data());}
Iterator y()const{return thrust::raw_pointer_cast(m_y.data());}
Iterator z()const{return thrust::raw_pointer_cast(m_z.data());}
Iterator x()const{return (Iterator)m_x.data().get();}
Iterator y()const{return (Iterator)m_y.data().get();}
Iterator z()const{return (Iterator)m_z.data().get();}

auto xyz() const{
auto zip = thrust::make_zip_iterator(thrust::make_tuple(x(), y(), z()));
Expand Down Expand Up @@ -292,7 +292,7 @@ namespace uammd{
FluidPointers(){}
template<class Container>
FluidPointers(const Container &pressure, const DataXYZ<T> &vel):
pressure(thrust::raw_pointer_cast(pressure.data())),
pressure((T*)(pressure.data().get())),
velocityX(vel.x()), velocityY(vel.y()), velocityZ(vel.z()){}
T* pressure;
typename DataXYZ<T>::Iterator velocityX, velocityY, velocityZ;
Expand Down
10 changes: 4 additions & 6 deletions src/Integrator/BDHI/FCM/FCM_impl.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,10 +81,6 @@ namespace uammd{
System::log<System::EXCEPTION>("FCM_impl requires instances of the spreading kernels");
throw std::runtime_error("Invalid arguments");
}
if(par.cells.x <= 0){
System::log<System::EXCEPTION>("FCM_impl requires valid grid dimensions");
throw std::runtime_error("Invalid arguments");
}
printMessages(par);
initCuFFT();
CudaSafeCall(cudaDeviceSynchronize());
Expand Down Expand Up @@ -149,8 +145,10 @@ namespace uammd{
auto rh = this->getHydrodynamicRadius();
auto M0 = this->getSelfMobility();
System::log<System::MESSAGE>("[BDHI::FCM] Using kernel: %s", type_name<Kernel>().c_str());
System::log<System::MESSAGE>("[BDHI::FCM] Closest possible hydrodynamic radius: %g", rh);
System::log<System::MESSAGE>("[BDHI::FCM] Self mobility: %g", (double)M0);
if(rh>0){
System::log<System::MESSAGE>("[BDHI::FCM] Hydrodynamic radius: %g", rh);
System::log<System::MESSAGE>("[BDHI::FCM] Self mobility: %g", (double)M0);
}
if(box.boxSize.x != box.boxSize.y || box.boxSize.y != box.boxSize.z || box.boxSize.x != box.boxSize.z){
System::log<System::WARNING>("[BDHI::FCM] Self mobility will be different for non cubic boxes!");
}
Expand Down
2 changes: 1 addition & 1 deletion test/BDHI/DPStokes/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ set(CMAKE_CUDA_STANDARD_REQUIRED ON)
set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF)
set(CUDA_ARCHITECTURES OFF)

set(UAMMD_INCLUDE ../../../../uammd/src ../../../../uammd/src/third_party)
set(UAMMD_INCLUDE ../../../src ../../../src/third_party)


include(FetchContent)
Expand Down
64 changes: 64 additions & 0 deletions test/BDHI/FCM/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
cmake_minimum_required(VERSION 3.14)
project(fcm_test)
enable_language(CUDA)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/")
set(CMAKE_BUILD_TYPE Release)
#set( CMAKE_VERBOSE_MAKEFILE on )
#add_compile_definitions(PUBLIC MAXLOGLEVEL=15)
#GoogleTest requires at least C++14
set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CUDA_STANDARD 14)
set(CMAKE_CUDA_STANDARD_REQUIRED ON)
set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF)
set(CUDA_ARCHITECTURES OFF)

set(UAMMD_INCLUDE ../../../src ../../../src/third_party)


include(FetchContent)
FetchContent_Declare(
googletest
GIT_REPOSITORY https://github.com/google/googletest.git
GIT_TAG release-1.12.1
)
# For Windows: Prevent overriding the parent project's compiler/linker settings
set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
FetchContent_MakeAvailable(googletest)
enable_testing()



#set(BLA_VENDOR Intel10_64lp)
find_package(BLAS)
if(BLAS_FOUND)
message("mkl environment detected")
add_compile_definitions(PUBLIC USE_MKL)
else()
unset(BLA_VENDOR)
find_package(BLAS REQUIRED)
find_package(LAPACK REQUIRED)
find_package(LAPACKE REQUIRED)
endif()

add_compile_definitions(PUBLIC DOUBLE_PRECISION)
include_directories(${UAMMD_INCLUDE})
link_libraries(${LAPACK_LIBRARIES} ${LAPACKE_LIBRARIES} ${BLAS_LIBRARIES} ${CUDA_LIBRARY} cufft cublas cusolver)

add_executable(
fcm_test
fcm_test.cu
)

target_link_libraries(
fcm_test
GTest::gtest_main
GTest::gmock_main
)

include(GoogleTest)
gtest_discover_tests(fcm_test)

IF (CMAKE_BUILD_TYPE MATCHES "Debug")
set(CMAKE_CUDA_FLAGS "-g -G -lineinfo -src-in-ptx")
set(CMAKE_CXX_FLAGS "-g -Wall -Wextra -pedantic")
ENDIF()
Loading

0 comments on commit 73335f7

Please sign in to comment.