...
 
Commits (7)
......@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
project(openfpm_numerics LANGUAGES C CXX)
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/CMakeFiles/)
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/../cmake_modules/)
set(BOOST_INCLUDE ${Boost_INCLUDE_DIR} CACHE PATH "Include directory for BOOST")
set(PETSC_ROOT CACHE PATH "If compiling with linear algebra indicate the PETSC root directory")
......@@ -29,10 +29,23 @@ set (CMAKE_CXX_STANDARD 11)
set (CMAKE_CUDA_STANDARD 11)
if(ENABLE_GPU)
enable_language(CUDA)
find_package(CUDA)
enable_language(CUDA)
find_package(CUDA)
if (CUDA_VERSION_MAJOR EQUAL 9 AND CUDA_VERSION_MINOR EQUAL 2)
message("CUDA is compatible")
set(WARNING_SUPPRESSION_AND_OPTION_NVCC -Xcudafe "--display_error_number --diag_suppress=611 --diag_suppress=2885 --diag_suppress=2886 --diag_suppress=2887 --diag_suppress=2888 --diag_suppress=186 --diag_suppress=111" --expt-extended-lambda)
FILE(WRITE cuda_options " -Xcudafe \"--display_error_number --diag_suppress=611 --diag_suppress=2885 --diag_suppress=2886 --diag_suppress=2887 --diag_suppress=2888 --diag_suppress=186 --diag_suppress=111\" --expt-extended-lambda ")
elseif ( CUDA_VERSION_MAJOR EQUAL 10 AND CUDA_VERSION_MINOR EQUAL 1 )
message("CUDA is compatible")
set(WARNING_SUPPRESSION_AND_OPTION_NVCC -Xcudafe "--display_error_number --diag_suppress=2915 --diag_suppress=2914 --diag_suppress=2912 --diag_suppress=2913 --diag_suppress=111 --diag_suppress=186 --diag_suppress=611 " --expt-extended-lambda )
FILE(WRITE cuda_options "-Xcudafe \"--display_error_number --diag_suppress=2915 --diag_suppress=2914 --diag_suppress=2912 --diag_suppress=2913 --diag_suppress=111 --diag_suppress=186 --diag_suppress=611 \" --expt-extended-lambda")
else()
message(FATAL_ERROR "CUDA is incompatible, version 9.2 and 10.1 is only supported")
endif()
endif()
find_package(Boost 1.66.0 REQUIRED COMPONENTS unit_test_framework iostreams program_options)
find_package(MPI REQUIRED)
find_package(PETSc)
......@@ -51,7 +64,7 @@ endif()
if(CUDA_FOUND)
set(OPENFPM_INIT_FILE "initialize/initialize_wrapper_cuda.cu")
else()
set(OPENFPM_INIT_FILE "initialize/initialize_wrapper_cuda.cpp")
set(OPENFPM_INIT_FILE "initialize/initialize_wrapper_cpu.cpp")
endif()
###### CONFIG.h FILE ######
......@@ -90,6 +103,10 @@ else()
message( FATAL_ERROR "MPI is required in order to install OpenFPM" )
endif()
if(PETSC_FOUND)
set(DEFINE_HAVE_PETSC "#define HAVE_PETSC")
endif()
if (Boost_FOUND)
set(DEFINE_HAVE_BOOST "#define HAVE_BOOST")
set(DEFINE_HAVE_BOOST_IOSTREAMS "#define HAVE_BOOST_IOSTREAMS")
......@@ -113,10 +130,20 @@ if(EIGEN_FOUND)
set(DEFINE_HAVE_EIGEN "#define HAVE_EIGEN")
endif()
if(EIGEN3_FOUND)
set(DEFINE_HAVE_EIGEN "#define HAVE_EIGEN")
endif()
if(LIBHILBERT_FOUND)
set(DEFINE_HAVE_LIBHILBERT "#define HAVE_LIBHILBERT 1")
set(DEFINE_HAVE_LIBHILBERT "#define HAVE_LIBHILBERT 1")
else()
message( FATAL_ERROR "LibHilbert is required in order to install OpenFPM")
file(WRITE error_code "210")
message( FATAL_ERROR "LibHilbert is required in order to install OpenFPM")
endif()
if(SUITESPARSE_FOUND AND SuiteSparse_UMFPACK_FOUND)
set(DEFINE_HAVE_SUITESPARSE "#define HAVE_SUITESPARSE")
endif()
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/src/config/config_cmake.h.in ${CMAKE_CURRENT_SOURCE_DIR}/src/config/config.h)
......
......@@ -2,13 +2,34 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
########################### Executables
add_executable(numerics main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp)
if (CUDA_FOUND)
set(CUDA_SOURCES Operators/Vector/vector_dist_operators_unit_tests.cu
Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cu)
endif()
add_executable(numerics ${OPENFPM_INIT_FILE} ${CUDA_SOURCES}
main.cpp
Matrix/SparseMatrix_unit_tests.cpp
interpolation/interpolation_unit_tests.cpp
Vector/Vector_unit_tests.cpp
Solvers/petsc_solver_unit_tests.cpp
FiniteDifference/FDScheme_unit_tests.cpp
FiniteDifference/eq_unit_test_3d.cpp
FiniteDifference/eq_unit_test.cpp
Operators/Vector/vector_dist_operators_unit_tests.cpp
Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cpp
../../openfpm_vcluster/src/VCluster/VCluster.cpp
../../openfpm_devices/src/memory/CudaMemory.cu
../../openfpm_devices/src/memory/HeapMemory.cpp
../../openfpm_devices/src/memory/PtrMemory.cpp
../../openfpm_devices/src/Memleak_check.cpp
../../src/lib/pdata.cpp)
###########################
if(CUDA_FOUND)
target_compile_options(numerics PUBLIC $<$<COMPILE_LANGUAGE:CUDA>:-Xcudafe "--display_error_number --diag_suppress=2885 --diag_suppress=2887 --diag_suppress=2888 --diag_suppress=186 --diag_suppress=111" --expt-extended-lambda>)
target_compile_options(numerics PUBLIC $<$<COMPILE_LANGUAGE:CUDA>: ${WARNING_SUPPRESSION_AND_OPTION_NVCC} >)
endif()
target_include_directories (numerics PUBLIC ${CUDA_INCLUDE_DIRS})
......@@ -28,8 +49,9 @@ if(EIGEN3_FOUND)
target_include_directories (numerics PUBLIC ${EIGEN3_INCLUDE_DIR})
endif()
link_directories(${PARMETIS_ROOT} ${METIS_ROOT})
target_link_libraries(numerics ${Boost_LIBRARIES})
target_link_libraries(numerics -L${PARMETIS_ROOT}/lib parmetis)
target_link_libraries(numerics ${PARMETIS_LIBRARIES})
target_link_libraries(numerics -L${HDF5_ROOT}/lib hdf5 hdf5_hl)
target_link_libraries(numerics -L${LIBHILBERT_LIBRARY_DIRS} ${LIBHILBERT_LIBRARIES})
if(PETSC_FOUND)
......
......@@ -249,7 +249,13 @@ template<typename solver_type,typename lid_nn> void lid_driven_cavity_2d()
#else
#if __GNUC__ == 7
#if __GNUC__ == 8
std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC7.vtk";
std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
#elif __GNUC__ == 7
std::string file1 = std::string("test/") + s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC7.vtk";
std::string file2 = s + "lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
......
......@@ -214,7 +214,12 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
#else
#if __GNUC__ == 7
#if __GNUC__ == 8
std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC8.vtk";
std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
#elif __GNUC__ == 7
std::string file1 = std::string("test/") + s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test_GCC7.vtk";
std::string file2 = s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()) + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk";
......
/*
* vector_dist_operators_cuda.cuh
*
* Created on: May 31, 2019
* Author: i-bird
*/
#ifndef VECTOR_DIST_OPERATORS_CUDA_CUH_
#define VECTOR_DIST_OPERATORS_CUDA_CUH_
constexpr unsigned int PROP_POS =(unsigned int)-1;
/*! \brief selector for position or properties left side expression
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector, unsigned int prp>
struct pos_or_propL
{
//! return the value (position or property) of the particle k in the vector v
static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(v.template getProp<prp>(k))
{
return v.template getProp<prp>(k);
}
};
/*! \brief selector for position or properties left side expression
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector, unsigned int prp>
struct pos_or_propL_ker
{
//! return the value (position or property) of the particle k in the vector v
__device__ static inline auto value(vector & v, const unsigned int & k) -> decltype(v.template getProp<prp>(k))
{
return v.template getProp<prp>(k);
}
};
/*! \brief selector for position or properties left side
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector>
struct pos_or_propL<vector,PROP_POS>
{
#ifdef SE_CLASS3
//! return the value (position or property) of the particle k in the vector v
static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(getExprL(v.getPos(k).getReference()))
{
return getExprL(v.getPos(k).getReference());
}
#else
//! return the value (position or property) of the particle k in the vector v
static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(ger<vector::dims,typename vector::stype>::getExprL(v.getPos(k)))
{
return ger<vector::dims,typename vector::stype>::getExprL(v.getPos(k));
}
#endif
};
/*! \brief selector for position or properties left side
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector>
struct pos_or_propL_ker<vector,PROP_POS>
{
#ifdef SE_CLASS3
//! return the value (position or property) of the particle k in the vector v
static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(getExprL(v.getPos(k).getReference()))
{
return getExprL(v.getPos(k).getReference());
}
#else
//! return the value (position or property) of the particle k in the vector v
__device__ static inline auto value(vector & v, const unsigned int & k) -> decltype(ger<vector::dims,typename vector::stype>::getExprL(v.getPos(k)))
{
return ger<vector::dims,typename vector::stype>::getExprL(v.getPos(k));
}
#endif
};
/*! \brief selector for position or properties right side position
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector, unsigned int prp>
struct pos_or_propR
{
//! return the value (position or property) of the particle k in the vector v
__device__ __host__ static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(v.template getProp<prp>(k))
{
return v.template getProp<prp>(k);
}
//! return the value (position or property) of the particle k in the vector v
__device__ __host__ static inline auto value(vector & v, const unsigned int & k) -> decltype(v.template getProp<prp>(k))
{
return v.template getProp<prp>(k);
}
};
/*! \brief selector for position or properties right side
*
* \tparam vector type of the original vector
*
* \tparam prp property id
*
*/
template <typename vector>
struct pos_or_propR<vector,PROP_POS>
{
//! return the value (position or property) of the particle k in the vector v
__device__ __host__ static inline auto value(vector & v, const vect_dist_key_dx & k) -> decltype(ger<vector::dims,typename vector::stype>::getExprR(v.getPos(k)))
{
return ger<vector::dims,typename vector::stype>::getExprR(v.getPos(k));
}
//! return the value (position or property) of the particle k in the vector v
__device__ __host__ static inline auto value(vector & v, const unsigned int & k) -> decltype(ger<vector::dims,typename vector::stype>::getExprR(v.getPos(k)))
{
return ger<vector::dims,typename vector::stype>::getExprR(v.getPos(k));
}
};
template<unsigned int prp ,bool is_sort, int impl>
struct vector_dist_op_compute_op
{
template<typename vector, typename expr>
static void compute_expr(vector & v,expr & v_exp)
{}
template<typename vector>
static void compute_const(vector & v,double d)
{}
};
template<unsigned int prp>
struct vector_dist_op_compute_op<prp,false,comp_host>
{
template<typename vector, typename expr>
static void compute_expr(vector & v,expr & v_exp)
{
v_exp.init();
auto it = v.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
pos_or_propL<vector,prp>::value(v,key) = v_exp.value(key);
++it;
}
}
template<typename vector>
static void compute_const(vector & v,double d)
{
auto it = v.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
pos_or_propL<vector,prp>::value(v,key) = d;
++it;
}
}
};
#ifdef __NVCC__
template<unsigned int prp, unsigned int dim ,typename vector, typename expr>
__global__ void compute_expr_ker_vv(vector vd, expr v_exp)
{
int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= vd.size()) {return;}
for (unsigned int i = 0 ; i < dim ; i++)
{
vd.template get<prp>(p)[i] = v_exp.value(p).get(i);
}
}
template<unsigned int prp, typename vector, typename expr>
__global__ void compute_expr_ker_v(vector vd, expr v_exp)
{
int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= vd.size()) {return;}
vd.template get<prp>(p) = v_exp.value(p);
}
template<unsigned int prp, typename vector, typename expr>
__global__ void compute_expr_ker(vector vd, expr v_exp)
{
int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= vd.size_local()) {return;}
pos_or_propL_ker<vector,prp>::value(vd,p) = v_exp.value(p);
}
template<unsigned int prp, typename vector>
__global__ void compute_double_ker(vector vd, double d)
{
int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= vd.size_local()) {return;}
pos_or_propL_ker<vector,prp>::value(vd,p) = d;
}
/////////// SORTED VERSION //
template<unsigned int prp, unsigned int dim ,typename vector, typename NN_type, typename expr>
__global__ void compute_expr_ker_sort_vv(vector vd, NN_type NN, expr v_exp)
{
int p;
GET_PARTICLE_SORT(p,NN);
for (unsigned int i = 0 ; i < dim ; i++)
{
vd.template get<prp>(p)[i] = v_exp.value(p).get(i);
}
}
template<unsigned int prp, typename vector, typename NN_type, typename expr>
__global__ void compute_expr_ker_sort_v(vector vd, NN_type NN, expr v_exp)
{
int p;
GET_PARTICLE_SORT(p,NN);
vd.template get<prp>(p) = v_exp.value(p);
}
template<unsigned int prp, typename vector, typename expr, typename NN_type>
__global__ void compute_expr_ker_sort(vector vd, NN_type NN, expr v_exp)
{
int p;
GET_PARTICLE_SORT(p,NN);
pos_or_propL_ker<vector,prp>::value(vd,p) = v_exp.value(p);
}
/////////////////////////////
template<unsigned int prp>
struct vector_dist_op_compute_op<prp,false,comp_dev>
{
template<typename vector, typename expr>
static void compute_expr(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getDomainIteratorGPU(256);
compute_expr_ker<prp><<<ite.wthr,ite.thr>>>(v,v_exp);
}
template<typename vector, typename expr>
static void compute_expr_v(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getGPUIterator(256);
compute_expr_ker_v<prp><<<ite.wthr,ite.thr>>>(v,v_exp);
}
template<unsigned int dim, typename vector, typename expr>
static void compute_expr_vv(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getGPUIterator(256);
compute_expr_ker_vv<prp,dim><<<ite.wthr,ite.thr>>>(v,v_exp);
}
template<typename vector>
static void compute_const(vector & v,double d)
{
auto ite = v.getDomainIteratorGPU(256);
compute_double_ker<prp><<<ite.wthr,ite.thr>>>(v,d);
}
};
template<unsigned int prp>
struct vector_dist_op_compute_op<prp,true,comp_dev>
{
template<typename vector, typename expr>
static void compute_expr(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getDomainIteratorGPU(256);
auto NN = v_exp.getNN();
compute_expr_ker_sort<prp><<<ite.wthr,ite.thr>>>(v,*NN,v_exp);
}
template<typename vector, typename expr>
static void compute_expr_v(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getGPUIterator(256);
auto NN = v_exp.getNN();
compute_expr_ker_sort_v<prp><<<ite.wthr,ite.thr>>>(v,*NN,v_exp);
}
template<unsigned int dim, typename vector, typename expr>
static void compute_expr_vv(vector & v,expr & v_exp)
{
v_exp.init();
auto ite = v.getGPUIterator(256);
auto NN = v_exp.getNN();
compute_expr_ker_sort_vv<prp,dim><<<ite.wthr,ite.thr>>>(v,*NN,v_exp);
}
};
#endif
#endif /* VECTOR_DIST_OPERATORS_CUDA_CUH_ */
/*
* vector_dist_operators_apply_kernel_unit_tests.hpp
*
* Created on: Jun 19, 2016
* Author: i-bird
*/
#include "config.h"
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "Vector/vector_dist.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Operators/Vector/tests/vector_dist_operators_tests_util.hpp"
BOOST_AUTO_TEST_SUITE( vector_dist_operators_apply_kernel_test_cpu )
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
{
if (create_vcluster().getProcessingUnits() > 3)
return;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
Ghost<3,float> ghost(0.05);
vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
check_all_apply_ker<comp_host>::check(vd);
}
BOOST_AUTO_TEST_SUITE_END()
/*
* vector_dist_operators_apply_kernel_unit_tests.hpp
*
* Created on: Jun 19, 2016
* Author: i-bird
*/
#include "config.h"
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "Vector/vector_dist.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Operators/Vector/tests/vector_dist_operators_tests_util.hpp"
BOOST_AUTO_TEST_SUITE( vector_dist_operators_apply_kernel_test_gpu )
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_host_gpu_test )
{
if (create_vcluster().getProcessingUnits() > 3)
return;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
Ghost<3,float> ghost(0.05);
vector_dist_gpu<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
check_all_apply_ker<comp_host>::check(vd);
}
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_gpu_test )
{
if (create_vcluster().getProcessingUnits() > 3)
return;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
Ghost<3,float> ghost(0.05);
vector_dist_gpu<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
check_all_apply_ker<comp_dev>::check(vd);
}
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_gpu_sort_test )
{
if (create_vcluster().getProcessingUnits() > 3)
return;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
Ghost<3,float> ghost(0.05);
vector_dist_gpu<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
check_all_apply_ker_sort::check(vd);
}
BOOST_AUTO_TEST_SUITE_END()
/*
* vector_dist_operators_apply_kernel_unit_tests.hpp
*
* Created on: Jun 19, 2016
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_UNIT_TESTS_HPP_
#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_UNIT_TESTS_HPP_
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel(vector & vd, Kernel & ker, NN_type & NN)
{
bool ret = true;
auto it = vd.getDomainIterator();
// ### WIKI 13 ###
//
// Check that apply kernel work
//
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
float base2 = 0.0;
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
float base1 = vd.template getProp<A>(p);
float prp_x = vd.template getProp<VC>(p) * vd.template getProp<VB>(p) + norm(vd.template getProp<VB>(p));
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
float prp_y = vd.template getProp<VC>(q) * vd.template getProp<VB>(q) + norm(vd.template getProp<VB>(q));
base2 += ker.value(xp,xq,prp_x,prp_y);
++Np;
}
base2 += vd.template getProp<C>(p);
ret &= base1 == base2;
++it2;
}
BOOST_REQUIRE_EQUAL(ret,true);
return ret;
}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel_reduce(vector & vd, Kernel & ker, NN_type & NN)
{
bool ret = true;
auto it = vd.getDomainIterator();
float base1 = 0.0;
float base2 = 0.0;
float base3 = 0.0;
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
float ker_accu = 0.0;
float prp_x = vd.template getProp<VC>(p) * vd.template getProp<VB>(p) + norm(vd.template getProp<VB>(p));
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
float prp_y = vd.template getProp<VC>(q) * vd.template getProp<VB>(q) + norm(vd.template getProp<VB>(q));
ker_accu += ker.value(xp,xq,prp_x,prp_y);
++Np;
}
base2 += ker_accu;
++it2;
}
auto it3 = vd.getDomainIterator();
while (it3.isNext())
{
auto p = it3.get();
base1 = vd.template getProp<A>(p);
base3 = vd.template getProp<C>(p) + base2;
ret &= base1 == base3;
++it3;
}
BOOST_REQUIRE_EQUAL(ret,true);
return ret;
}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel2(vector & vd, Kernel & ker, NN_type & NN)
{
bool ret = true;
auto it = vd.getDomainIterator();
// ### WIKI 13 ###
//
// Check that apply kernel work
//
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
Point<3,float> base2 = 0.0;
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
Point<3,float> base1 = vd.template getProp<VA>(p);
Point<3,float> prp_x = 2.0 * vd.template getProp<VC>(p) + vd.template getProp<VB>(p);
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
Point<3,float> prp_y = 2.0 * vd.template getProp<VC>(q) + vd.template getProp<VB>(q);
base2 += ker.value(xp,xq,prp_x,prp_y);
++Np;
}
base2 += vd.template getProp<VC>(p);
ret &= base1 == base2;
++it2;
}
BOOST_REQUIRE_EQUAL(ret,true);
return ret;
}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3(vector & vd, Kernel & ker, NN_type & NN)
{
bool ret = true;
auto it = vd.getDomainIterator();
// ### WIKI 13 ###
//
// Check that apply kernel work
//
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
Point<3,float> base2 = 0.0;
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
Point<3,float> base1 = vd.template getProp<VA>(p);
Point<3,float> prp_x = vd.template getProp<VC>(p);
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
Point<3,float> prp_y = vd.template getProp<VC>(q);
base2 += ker.value(xp,xq,prp_x,prp_y);
++Np;
}
base2 += vd.template getProp<VC>(p);
ret &= base1 == base2;
++it2;
}
BOOST_REQUIRE_EQUAL(ret,true);
return ret;
}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel2_reduce(vector & vd, Kernel & ker, NN_type & NN)
{
bool ret = true;
auto it = vd.getDomainIterator();
Point<3,float> base1 = 0.0;
Point<3,float> base2 = 0.0;
Point<3,float> base3 = 0.0;
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
Point<3,float> ker_accu = 0.0;
Point<3,float> prp_x = 2.0f*vd.template getProp<VC>(p) + vd.template getProp<VB>(p);
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
Point<3,float> prp_y = 2.0f*vd.template getProp<VC>(q) + vd.template getProp<VB>(q);
ker_accu += ker.value(xp,xq,prp_x,prp_y);
++Np;
}
base2 += ker_accu;
++it2;
}
auto it3 = vd.getDomainIterator();
while (it3.isNext())
{
auto p = it3.get();
base1 = vd.template getProp<VA>(p);
base3 = vd.template getProp<VC>(p) + base2;
ret &= base1 == base3;
++it3;
}
BOOST_REQUIRE_EQUAL(ret,true);
return ret;
}
template <typename vector,typename Kernel, typename NN_type> bool check_values_apply_kernel3_reduce(vector & vd, Kernel & ker, NN_type & NN, const Point<2,float> & p)
{
bool ret = true;
auto it = vd.getDomainIterator();
Point<2,float> base2 = 0.0;
auto it2 = vd.getDomainIterator();
while (it2.isNext())
{
auto p = it2.get();
Point<3,float> xp = vd.getPos(p);
Point<2,float> ker_accu = 0.0;
// For each neighborhood particle
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(xp));
while (Np.isNext())
{
// Neighborhood particle q
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
// position q
Point<3,float> xq = vd.getPos(q);
ker_accu += ker.value(xp,xq);
++Np;
}
base2 += ker_accu;
++it2;
}
BOOST_REQUIRE_EQUAL(p.get(0),base2.get(0));
BOOST_REQUIRE_EQUAL(p.get(1),base2.get(1));
return ret;
}
BOOST_AUTO_TEST_CASE( vector_dist_operators_apply_kernel_test )
{
if (create_vcluster().getProcessingUnits() > 3)
return;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// ghost
Ghost<3,float> ghost(0.05);
vector_dist<3,float,aggregate<float,float,float,VectorS<3,float>,VectorS<3,float>,VectorS<3,float>,float>> vd(512,box,bc,ghost);
auto vA = getV<A>(vd);
auto vC = getV<C>(vd);
auto vVA = getV<VA>(vd);
auto vVB = getV<VB>(vd);
auto vVC = getV<VC>(vd);
// fill vd with some value
fill_values(vd);
vd.map();
vd.ghost_get<0,1,2,3,4,5,6>();
// we apply an exponential kernel to calculate something
auto cl = vd.getCellList(0.05);
exp_kernel ker(0.2);
vA = applyKernel_in(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
check_values_apply_kernel(vd,ker,cl);
vVA = applyKernel_in(2.0*vVC + vVB ,vd,cl,ker) + vVC;
check_values_apply_kernel2(vd,ker,cl);
vA = rsum(applyKernel_in(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
check_values_apply_kernel_reduce(vd,ker,cl);
vVA = rsum(applyKernel_in(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
check_values_apply_kernel2_reduce(vd,ker,cl);
vA = applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker) + vC;
check_values_apply_kernel(vd,ker,cl);
vVA = applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker) + vVC;
check_values_apply_kernel2(vd,ker,cl);
vA = rsum(applyKernel_in_gen(vVC * vVB + norm(vVB),vd,cl,ker),vd) + vC;
check_values_apply_kernel_reduce(vd,ker,cl);
vVA = rsum(applyKernel_in_gen(2.0*vVC + vVB ,vd,cl,ker),vd) + vVC;
check_values_apply_kernel2_reduce(vd,ker,cl);
// Check it compile the code is the same
vVA = applyKernel_in_gen(vVC,vd,cl,ker) + vVC;
check_values_apply_kernel3(vd,ker,cl);
vVA = applyKernel_in (vVC,vd,cl,ker) + vVC;
check_values_apply_kernel3(vd,ker,cl);
// just check it compile
Point<2,float> p = rsum(applyKernel_in_sim(vd,cl,ker),vd).get();
check_values_apply_kernel3_reduce(vd,ker,cl,p);
}
#endif /* OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_APPLY_KERNEL_UNIT_TESTS_HPP_ */
......@@ -15,7 +15,8 @@
* \param v
*
*/
template <unsigned int dim, typename T> inline vector_dist_expression<16384,Point<dim,T> > getVExpr(Point<dim,T> & v)
template <unsigned int dim, typename T>
inline vector_dist_expression<16384,Point<dim,T> > getVExpr(Point<dim,T> & v)
{
vector_dist_expression<(unsigned int)16384,Point<dim,T>> exp_v(v);
......@@ -36,6 +37,14 @@ class vector_dist_expression<16384,point>
public:
typedef void vtype;
//! result for is sort
typedef boost::mpl::bool_<false> is_sort;
//! NN_type
typedef void NN_type;
//! vector expression from a constant point
vector_dist_expression(point p)
:p(p)
......
This diff is collapsed.
......@@ -186,6 +186,8 @@ public:
static Vector<double,EIGEN_BASE> try_solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error in order to use umfpack you must compile OpenFPM with linear algebra support" << "/n";
return Vector<double,EIGEN_BASE> ();
}
};
......
/*
* initialize_vcl.hpp
*
* Created on: Aug 21, 2018
* Author: i-bird
*/
#ifndef INITIALIZE_VCL_HPP_
#define INITIALIZE_VCL_HPP_
/*! \brief If openfpm has to work on GPU we have to be sure openfpm_init is called on a file compiled with NVCC
*
* There are two implementation initialize.cpp and initialize.cu. In configuration stage the second implementation is chosen
* if the test has to run on GPU
*
*/
void openfpm_init_wrapper(int * argc, char *** argv);
void openfpm_finalize_wrapper();
#endif /* INITIALIZE_VCL_HPP_ */
#include "initialize_wrapper.hpp"
#include "VCluster/VCluster.hpp"
void openfpm_init_wrapper(int * argc, char *** argv)
{
openfpm_init(argc,argv);
}
void openfpm_finalize_wrapper()
{
openfpm_finalize();
}
This diff is collapsed.
This diff is collapsed.