Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize the classes to use complex or real numbers #21

Open
wants to merge 12 commits into
base: v2.x
Choose a base branch
from
4 changes: 2 additions & 2 deletions src/Integrator/BDHI/BDHI_Cholesky.cu
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ namespace uammd{
//Curand fill with gaussian numbers with mean 0 and var 1
/*This shit is obscure, curand will only work with an even number of elements*/
auto d_noise = thrust::raw_pointer_cast(noise.data());
curandGenerateNormal(curng, d_noise, 3*numberParticles + ((3*numberParticles)%2), real(0.0), real(1.0));
curandgeneratenormal(curng, d_noise, 3*numberParticles + ((3*numberParticles)%2), real(0.0), real(1.0));
isMup2date = false;
}

Expand Down Expand Up @@ -219,7 +219,7 @@ namespace uammd{
3*numberParticles, d_M, 3*numberParticles, d_work, h_work_size, d_info));
curandSetStream(curng, st);
/*Gen new noise in BdW*/
curandGenerateNormal(curng,
curandgeneratenormal(curng,
(real*) BdW,
3*numberParticles + ((3*numberParticles)%2),
real(0.0), real(1.0));
Expand Down
6 changes: 3 additions & 3 deletions src/Integrator/BDHI/BDHI_Lanczos.cu
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,8 @@ namespace uammd{
thrust::device_vector<real> noise(30000);
auto noise_ptr = thrust::raw_pointer_cast(noise.data());
//Warm cuRNG
curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0);
curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0);
curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0);
curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0);
}

namespace Lanczos_ns{
Expand Down Expand Up @@ -172,7 +172,7 @@ namespace uammd{
Lanczos_ns::Dotctor<real3> Mdot(rpy, this->hydrodynamicRadius, radius_ptr, nbody, st);
//Filling V instead of an external array (for v in sqrt(M)·v) is faster
uninitialized_cached_vector<real3> noise(numberParticles);
curandGenerateNormal(curng, (real*)noise.data().get(),
curandgeneratenormal(curng, (real*)noise.data().get(),
3*numberParticles + (3*numberParticles)%2,
real(0.0), real(1.0));
//lanczosAlgorithm->solve(Mdot, (real*) BdW, noise, numberParticles, st);
Expand Down
2 changes: 1 addition & 1 deletion src/Integrator/BDHI/DoublyPeriodic/DPStokesSlab.cuh
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ else throw std::runtime_error("[DPStokesSlab] Can only average in direction X (0
real gw;
real tolerance;
WallMode mode;
shared_ptr<BVP::BatchedBVPHandler> bvpSolver;
shared_ptr<BVP::BatchedBVPHandler<real>> bvpSolver;

};

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -506,7 +506,7 @@ namespace uammd{
}

auto computeZeroModeBoundaryConditions(int nz, real H, WallMode mode){
BVP::SchurBoundaryCondition bcs(nz, H);
BVP::SchurBoundaryCondition<real> bcs(nz, H);
if(mode == WallMode::bottom){
correction_ns::TopBoundaryConditions top(0, H);
correction_ns::BottomBoundaryConditions bot(0, H);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ namespace uammd{
auto botBC = thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), botdispatch);
int numberSystems = (nk.x/2+1)*nk.y;
int nz = grid.cellDim.z;
this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler>(klist, topBC, botBC, numberSystems, halfH, nz);
this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler<real>>(klist, topBC, botBC, numberSystems, halfH, nz);
CudaCheckError();
}

Expand Down
10 changes: 5 additions & 5 deletions src/Integrator/BDHI/FIB/FIB.cu
Original file line number Diff line number Diff line change
Expand Up @@ -129,8 +129,8 @@ namespace uammd{
thrust::device_vector<real> noise(30000);
auto noise_ptr = thrust::raw_pointer_cast(noise.data());
//Warm cuRNG
CurandSafeCall(curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
}
catch(thrust::system_error &e){
sys->log<System::CRITICAL>("[BDHI::FIB] Thrust could not allocate necessary arrays at initialization with error: %s", e.what());
Expand Down Expand Up @@ -871,7 +871,7 @@ namespace uammd{
sys->log<System::DEBUG2>("[BDHI::FIB] Random advection");
if(temperature!=real(0.0)){
sys->log<System::DEBUG2>("[BDHI::FIB] Generate random numbers");
CurandSafeCall(curandGenerateNormal(curng,
CurandSafeCall(curandgeneratenormal(curng,
thrust::raw_pointer_cast(random.data()),
random.size(),
0.0, 1.0));
Expand Down Expand Up @@ -906,7 +906,7 @@ namespace uammd{
sys->log<System::DEBUG2>("[BDHI::FIB] Random advection");
if(temperature!=real(0.0)){
sys->log<System::DEBUG2>("[BDHI::FIB] Generate random numbers");
CurandSafeCall(curandGenerateNormal(curng,
CurandSafeCall(curandgeneratenormal(curng,
thrust::raw_pointer_cast(random.data()),
random.size(),
0.0, 1.0));
Expand All @@ -932,7 +932,7 @@ namespace uammd{
double dV = grid.cellSize.x*grid.cellSize.y*grid.cellSize.z;
real noisePrefactor = sqrt(viscosity*temperature/(dt*dV));
if(temperature!=real(0.0))
CurandSafeCall(curandGenerateNormal(curng,
CurandSafeCall(curandgeneratenormal(curng,
thrust::raw_pointer_cast(random.data()),
random.size(),
0.0, 1.0));
Expand Down
6 changes: 3 additions & 3 deletions src/Integrator/Hydro/ICM.cu
Original file line number Diff line number Diff line change
Expand Up @@ -935,8 +935,8 @@ namespace uammd{
thrust::device_vector<real> noise(30000);
auto noise_ptr = thrust::raw_pointer_cast(noise.data());
//Warm cuRNG
CurandSafeCall(curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandGenerateNormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
CurandSafeCall(curandgeneratenormal(curng, noise_ptr, noise.size(), 0.0, 1.0));
}
//Sum S·F term using the current particle positions
void ICM::spreadParticleForces(){
Expand Down Expand Up @@ -993,7 +993,7 @@ namespace uammd{
if(temperature!=real(0.0)){
CurandSafeCall(curandSetStream(curng, st));
sys->log<System::DEBUG2>("[Hydro::ICM] Generate random numbers");
CurandSafeCall(curandGenerateNormal(curng, thrust::raw_pointer_cast(random.data()), random.size(), 0.0, 1.0));
CurandSafeCall(curandgeneratenormal(curng, thrust::raw_pointer_cast(random.data()), random.size(), 0.0, 1.0));
}
real3* d_gridVels = (real3*)thrust::raw_pointer_cast(gridVels.data());
real3* d_cellAdvection = (real3*)thrust::raw_pointer_cast(cellAdvection.data());
Expand Down
5 changes: 0 additions & 5 deletions src/Integrator/VerletNVE.cu
Original file line number Diff line number Diff line change
Expand Up @@ -28,11 +28,6 @@

#include"VerletNVE.cuh"


#ifndef SINGLE_PRECISION
#define curandGenerateNormal curandGenerateNormalDouble
#endif

namespace uammd{

VerletNVE::VerletNVE(shared_ptr<ParticleGroup> pg, VerletNVE::Parameters par):
Expand Down
4 changes: 2 additions & 2 deletions src/Interactor/DoublyPeriodic/PoissonSlab/BVPPoisson.cuh
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ namespace uammd{
}

class BVPPoissonSlab{
std::shared_ptr<BVP::BatchedBVPHandler> bvpSolver;
std::shared_ptr<BVP::BatchedBVPHandler<real>> bvpSolver;
real2 Lxy;
real H;
int3 cellDim;
Expand Down Expand Up @@ -150,7 +150,7 @@ namespace uammd{
BoundaryConditionsDispatch<BottomBoundaryConditions, decltype(klist)>(klist, H));

int numberSystems = (nk.x/2+1)*nk.y;
this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler>(klist, topBC, bottomBC, numberSystems, H, cellDim.z);
this->bvpSolver = std::make_shared<BVP::BatchedBVPHandler<real>>(klist, topBC, bottomBC, numberSystems, H, cellDim.z);
CudaCheckError();
}
};
Expand Down
25 changes: 14 additions & 11 deletions src/misc/BoundaryValueProblem/BVPSchurComplementMatrices.cuh
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ namespace uammd{
}
};

template <typename T>
class SchurBoundaryCondition{
int nz;
real H;
Expand All @@ -125,8 +126,8 @@ namespace uammd{
SchurBoundaryCondition(int nz, real H):nz(nz), H(H), bcs(nz){}

template<class TopBC, class BottomBC>
std::vector<real> computeBoundaryConditionMatrix(const TopBC &top, const BottomBC &bottom){
std::vector<real> CandD(2*nz+4, 0);
std::vector<T> computeBoundaryConditionMatrix(const TopBC &top, const BottomBC &bottom){
std::vector<T> CandD(2*nz+4, 0);
auto topRow = computeTopRow(top, bottom);
auto bottomRow = computeBottomRow(top, bottom);
std::copy(topRow.begin(), topRow.end()-2, CandD.begin());
Expand All @@ -141,8 +142,8 @@ namespace uammd{
private:

template<class TopBC, class BottomBC>
std::vector<real> computeTopRow(const TopBC &top, const BottomBC &bottom){
std::vector<real> topRow(nz+2, 0);
std::vector<T> computeTopRow(const TopBC &top, const BottomBC &bottom){
std::vector<T> topRow(nz+2, 0);
auto tfi = bcs.topFirstIntegral();
auto tsi = bcs.topSecondIntegral();
auto tfiFactor = top.getFirstIntegralFactor();
Expand All @@ -156,8 +157,8 @@ namespace uammd{
}

template<class TopBC, class BottomBC>
std::vector<real> computeBottomRow(const TopBC &top, const BottomBC &bottom){
std::vector<real> bottomRow(nz+2, 0);
std::vector<T> computeBottomRow(const TopBC &top, const BottomBC &bottom){
std::vector<T> bottomRow(nz+2, 0);
auto bfi = bcs.bottomFirstIntegral();
auto bsi = bcs.bottomSecondIntegral();
auto bfiFactor = bottom.getFirstIntegralFactor();
Expand All @@ -172,10 +173,11 @@ namespace uammd{

};

std::vector<real> computeSecondIntegralMatrix(real k, real H, int nz){
std::vector<real> A(nz*nz, 0);
template<class T>
std::vector<T> computeSecondIntegralMatrix(T k, real H, int nz){
std::vector<T> A(nz*nz, 0);
SecondIntegralMatrix sim(nz);
real kH2 = k*k*H*H;
T kH2 = k*k*H*H;
fori(0, nz){
forj(0,nz){
A[i+nz*j] = (i==j) - kH2*sim.getElement(i, j);
Expand All @@ -184,9 +186,10 @@ namespace uammd{
return std::move(A);
}

std::vector<real> computeInverseSecondIntegralMatrix(real k, real H, int nz){
template<class T>
std::vector<T> computeInverseSecondIntegralMatrix(T k, real H, int nz){
if(k==0){
std::vector<real> invA(nz*nz, 0);
std::vector<T> invA(nz*nz, 0);
fori(0, nz){
invA[i+nz*i] = 1;
}
Expand Down