Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • mosaic/software/parallel-computing/openfpm/openfpm_numerics
  • argupta/openfpm_numerics
2 results
Show changes
Commits on Source (52)
Hide whitespace changes
Inline Side-by-side
Showing
with 2520 additions and 681 deletions
run.sh
View file @ 1371f3eb
......@@ -5,7 +5,7 @@
#cd openfpm_numerics
workspace=$1
hostname=$2
hostname=$(cat hostname)
nproc=$3
ntask_per_node=$5
nodes=$4
......@@ -15,6 +15,10 @@ echo "Directory: $workspace"
echo "Machine: $hostname"
echo "Branch name: $branch"
if [ x"$hostname" == x"cifarm-centos-node" ]; then
mpi_options="--allow-run-as-root"
fi
# Get the branch of pdata
echo "RUN numerics test"
......@@ -24,7 +28,7 @@ source $HOME/openfpm_vars_$branch
cd openfpm_numerics
mpirun -np $3 ../build/openfpm_numerics/src/numerics
mpirun $mpi_options -np $3 ../build/openfpm_numerics/src/numerics
if [ $? -ne 0 ]; then
curl -X POST --data "payload={\"icon_emoji\": \":jenkins:\", \"username\": \"jenkins\" , \"attachments\":[{ \"title\":\"Error:\", \"color\": \"#FF0000\", \"text\":\"$2 failed to complete the openfpm_numerics test \" }] }" https://hooks.slack.com/services/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
exit 1 ;
......
src/CMakeLists.txt
View file @ 1371f3eb
......@@ -18,6 +18,8 @@ endif()
if (NOT CUDA_ON_BACKEND STREQUAL "None")
set(CUDA_SOURCES Operators/Vector/vector_dist_operators_unit_tests.cu
DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cu
OdeIntegrators/tests/Odeintegrators_test_gpu.cu
#DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cu
Operators/Vector/vector_dist_operators_apply_kernel_unit_tests.cu)
endif()
......@@ -50,7 +52,7 @@ if ( CUDA_ON_BACKEND STREQUAL "HIP" AND HIP_FOUND )
DCPSE/tests/Vandermonde_unit_tests.cpp
main.cpp
Matrix/SparseMatrix_unit_tests.cpp
interpolation/interpolation_unit_tests.cpp
#interpolation/interpolation_unit_tests.cpp
Vector/Vector_unit_tests.cpp
Solvers/petsc_solver_unit_tests.cpp
FiniteDifference/FDScheme_unit_tests.cpp
......@@ -91,7 +93,7 @@ else()
DCPSE/tests/Vandermonde_unit_tests.cpp
main.cpp
Matrix/SparseMatrix_unit_tests.cpp
interpolation/interpolation_unit_tests.cpp
#interpolation/interpolation_unit_tests.cpp
Vector/Vector_unit_tests.cpp
Solvers/petsc_solver_unit_tests.cpp
FiniteDifference/FDScheme_unit_tests.cpp
......@@ -250,7 +252,7 @@ install(FILES FiniteDifference/Average.hpp
FiniteDifference/FDScheme.hpp
FiniteDifference/Laplacian.hpp
FiniteDifference/mul.hpp
FiniteDifference/sum.hpp
FiniteDifference/sum.hpp
FiniteDifference/Upwind_gradient.hpp
FiniteDifference/Eno_Weno.hpp
FiniteDifference/FD_op.hpp
......@@ -278,7 +280,7 @@ install(FILES Operators/Vector/vector_dist_operators_extensions.hpp
COMPONENT OpenFPM)
install(FILES Operators/Vector/cuda/vector_dist_operators_cuda.cuh
DESTINATION openfpm_numerics/include/Operators/Vector/cuda
DESTINATION openfpm_numerics/include/Operators/Vector/cuda
COMPONENT OpenFPM)
install(FILES DCPSE/Dcpse.hpp
......@@ -313,7 +315,8 @@ install(FILES DCPSE/DCPSE_op/DCPSE_op.hpp
COMPONENT OpenFPM)
install(FILES OdeIntegrators/OdeIntegrators.hpp
OdeIntegrators/boost_vector_algebra_ofp.hpp
OdeIntegrators/vector_algebra_ofp.hpp
OdeIntegrators/vector_algebra_ofp_gpu.hpp
DESTINATION openfpm_numerics/include/OdeIntegrators
COMPONENT OpenFPM)
......@@ -322,7 +325,7 @@ install(FILES Draw/DrawParticles.hpp
Draw/PointIteratorSkin.hpp
Draw/DrawDisk.hpp
Draw/DrawSphere.hpp
DESTINATION openfpm_numerics/include/Draw
DESTINATION openfpm_numerics/include/Draw
COMPONENT OpenFPM)
install(FILES interpolation/interpolation.hpp
......
src/DCPSE/DCPSE_op/DCPSE_Solver.hpp
View file @ 1371f3eb
......@@ -125,11 +125,11 @@ class DCPSE_scheme {
} else if (opt == options_solver::LAGRANGE_MULTIPLIER) {
if (v_cl.rank() == v_cl.size() - 1) {
b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz + 1);
x_ig.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz + 1);
b.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * (sz + 1));
x_ig.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * (sz + 1));
} else {
b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz);
x_ig.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz);
b.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * sz);
x_ig.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * sz);
}
}
//Use Custom number of constraints using opt as an integer
......@@ -252,8 +252,9 @@ class DCPSE_scheme {
// Indicate all the non zero rows
openfpm::vector<unsigned char> nz_rows;
if (v_cl.rank() == v_cl.size()-1 && opt == options_solver::LAGRANGE_MULTIPLIER) {
nz_rows.resize(row_b+1);
nz_rows.resize(row_b+Sys_eqs::nvar);
}
else{
nz_rows.resize(row_b);
......@@ -273,7 +274,7 @@ class DCPSE_scheme {
if (v_cl.getProcessingUnits() == 1) {
openfpm::vector<unsigned> nz_cols;
if (v_cl.rank() == v_cl.size()-1 && opt == options_solver::LAGRANGE_MULTIPLIER) {
nz_cols.resize(row_b+1);
nz_cols.resize(row_b+Sys_eqs::nvar);
}
else{
nz_cols.resize(row_b);
......@@ -438,6 +439,54 @@ public:
copy_nested(x, comp, exps ...);
}
/*! \brief Successive Solve an equation
*
* \warning exp must be a scalar type
*
* \param Solver Manually created Solver instead from the Equation structure
* \param exp where to store the result
*
*/
template<typename SolverType, typename ... expr_type>
void solve_with_solver_successive(SolverType &solver,expr_type ... exps) {
#ifdef SE_CLASS1
if (sizeof...(exps) != Sys_eqs::nvar) {
std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
dimensionality, I am expecting " << Sys_eqs::nvar <<
" properties " << std::endl;
};
#endif
auto x = solver.solve_successive(getB(opt));
unsigned int comp = 0;
copy_nested(x, comp, exps ...);
}
/*! \brief Successive Solve an equation with inital guess
*
* \warning exp must be a scalar type
*
* \param Solver Manually created Solver instead from the Equation structure
* \param exp where to store the result
*
*/
template<typename SolverType, typename ... expr_type>
void solve_with_solver_ig_successive(SolverType &solver,expr_type ... exps) {
#ifdef SE_CLASS1
if (sizeof...(exps) != Sys_eqs::nvar) {
std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
dimensionality, I am expecting " << Sys_eqs::nvar <<
" properties " << std::endl;
};
#endif
auto x = solver.solve_successive(get_x_ig(opt),getB(opt));
unsigned int comp = 0;
copy_nested(x, comp, exps ...);
}
/*! \brief Solve an equation with a given Nullspace
*
* \warning exp must be a scalar type
......@@ -480,7 +529,7 @@ public:
" properties " << std::endl;
};
#endif
auto x = solver.with_constant_nullspace_solve(getA(opt), getB(opt));
auto x = solver.with_nullspace_solve(getA(opt), getB(opt));
unsigned int comp = 0;
copy_nested(x, comp, exps ...);
......@@ -758,6 +807,7 @@ public:
*/
template<typename options>
typename Sys_eqs::SparseMatrix_type &getA(options opt) {
if (A.isMatrixFilled()) return A;
if (opt == options_solver::STANDARD) {
A.resize(tot * Sys_eqs::nvar, tot * Sys_eqs::nvar,
p_map.size_local() * Sys_eqs::nvar,
......@@ -768,60 +818,46 @@ public:
openfpm::vector<triplet> &trpl = A.getMatrixTriplets();
if (v_cl.rank() == v_cl.size() - 1) {
A.resize(tot * Sys_eqs::nvar + 1, tot * Sys_eqs::nvar + 1,
p_map.size_local() * Sys_eqs::nvar + 1,
p_map.size_local() * Sys_eqs::nvar + 1);
for (int i = 0; i < tot * Sys_eqs::nvar; i++) {
triplet t1;
t1.row() = tot * Sys_eqs::nvar;
t1.col() = i;
t1.value() = 1;
trpl.add(t1);
}
for (int i = 0; i < p_map.size_local() * Sys_eqs::nvar; i++) {
triplet t2;
t2.row() = i + s_pnt * Sys_eqs::nvar;
t2.col() = tot * Sys_eqs::nvar;
t2.value() = 1;
trpl.add(t2);
A.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * (tot + 1),
Sys_eqs::nvar * (p_map.size_local() + 1),
Sys_eqs::nvar * (p_map.size_local() + 1));
for (int j = 0; j < Sys_eqs::nvar; j++) {
for (int i = 0; i < tot; i++) {
triplet t1;
t1.row() = tot * Sys_eqs::nvar + j;
t1.col() = i * Sys_eqs::nvar + j;
t1.value() = 1;
trpl.add(t1);
}
for (int i = 0; i < p_map.size_local(); i++) {
triplet t2;
t2.row() = s_pnt * Sys_eqs::nvar + i * Sys_eqs::nvar + j;
t2.col() = tot * Sys_eqs::nvar + j;
t2.value() = 1;
trpl.add(t2);
}
triplet t3;
t3.col() = tot * Sys_eqs::nvar + j;
t3.row() = tot * Sys_eqs::nvar + j;
t3.value() = 0;
trpl.add(t3);
}
triplet t3;
t3.col() = tot * Sys_eqs::nvar;
t3.row() = tot * Sys_eqs::nvar;
t3.value() = 0;
trpl.add(t3);
//row_b++;
//row++;
}
else {
A.resize(tot * Sys_eqs::nvar + 1, tot * Sys_eqs::nvar + 1,
} else {
A.resize(Sys_eqs::nvar * (tot + 1), Sys_eqs::nvar * (tot + 1),
p_map.size_local() * Sys_eqs::nvar,
p_map.size_local() * Sys_eqs::nvar);
for (int i = 0; i < p_map.size_local() * Sys_eqs::nvar; i++) {
triplet t2;
t2.row() = i + s_pnt * Sys_eqs::nvar;
t2.col() = tot * Sys_eqs::nvar;
t2.value() = 1;
trpl.add(t2);
for (int j = 0; j < Sys_eqs::nvar; j++) {
for (int i = 0; i < p_map.size_local(); i++) {
triplet t2;
t2.row() = s_pnt * Sys_eqs::nvar + i * Sys_eqs::nvar + j;
t2.col() = tot * Sys_eqs::nvar + j;
t2.value() = 1;
trpl.add(t2);
}
}
}
}
else{
auto &v_cl = create_vcluster();
if (v_cl.rank() == v_cl.size() - 1) {
......@@ -833,8 +869,8 @@ public:
A.resize(tot * Sys_eqs::nvar - offset, tot * Sys_eqs::nvar - offset,
p_map.size_local() * Sys_eqs::nvar,
p_map.size_local() * Sys_eqs::nvar);
}
}
}
#ifdef SE_CLASS1
consistency(opt);
#endif
......@@ -854,8 +890,8 @@ public:
if (opt == options_solver::LAGRANGE_MULTIPLIER) {
auto &v_cl = create_vcluster();
if (v_cl.rank() == v_cl.size() - 1) {
b(tot * Sys_eqs::nvar) = 0;
for(int j=0;j<Sys_eqs::nvar;j++)
{b(tot * Sys_eqs::nvar+j) = 0;}
}
}
return b;
......@@ -873,8 +909,8 @@ public:
if (opt == options_solver::LAGRANGE_MULTIPLIER) {
auto &v_cl = create_vcluster();
if (v_cl.rank() == v_cl.size() - 1) {
x_ig(tot * Sys_eqs::nvar) = 0;
for(int j=0;j<Sys_eqs::nvar;j++)
{x_ig(tot * Sys_eqs::nvar+j) = 0;}
}
}
return x_ig;
......@@ -959,14 +995,6 @@ public:
// get the particle
auto key = it.get();
/*
if (key == 298 && create_vcluster().rank() == 1)
{
int debug = 0;
debug++;
}
*/
// Calculate the non-zero colums
typename Sys_eqs::stype coeff = 1.0;
op.template value_nz<Sys_eqs>(p_map, key, cols, coeff, 0);
......@@ -975,11 +1003,11 @@ public:
bool is_diag = false;
// create the triplet
for (auto it = cols.begin(); it != cols.end(); ++it) {
for (auto it2 = cols.begin(); it2 != cols.end(); ++it2) {
trpl.add();
trpl.last().row() = p_map.template getProp<0>(key) * Sys_eqs::nvar + id;
trpl.last().col() = it->first;
trpl.last().value() = it->second;
trpl.last().col() = it2->first;
trpl.last().value() = it2->second;
if (trpl.last().row() == trpl.last().col())
{is_diag = true;}
}
......
src/DCPSE/DCPSE_op/DCPSE_op.hpp
View file @ 1371f3eb
......@@ -7,6 +7,7 @@
#ifndef DCPSE_OP_HPP_
#define DCPSE_OP_HPP_
#ifdef HAVE_EIGEN
#include "Decomposition/CartDecomposition.hpp"
#include "DCPSE/Dcpse.hpp"
......@@ -2664,5 +2665,5 @@ typedef Derivative_yyy_T<Dcpse_gpu> Derivative_yyy_gpu;
typedef Derivative_G_T<Dcpse_gpu> Derivative_G_gpu;
#endif
#endif /*EIGEN */
#endif /* DCPSE_OP_HPP_ */
src/DCPSE/DCPSE_op/DCPSE_surface_op.hpp
View file @ 1371f3eb
......@@ -4,6 +4,7 @@
#ifndef OPENFPM_PDATA_DCPSE_SURFACE_OP_HPP
#define OPENFPM_PDATA_DCPSE_SURFACE_OP_HPP
#ifdef HAVE_EIGEN
#include "DCPSE/DCPSE_op/DCPSE_op.hpp"
......@@ -91,7 +92,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
}
......@@ -181,7 +182,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -271,7 +272,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -363,8 +364,8 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
}
......@@ -453,8 +454,8 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
}
......@@ -544,7 +545,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -726,7 +727,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -818,7 +819,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -910,7 +911,7 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
......@@ -997,11 +998,11 @@ public:
template<typename particles_type>
void update(particles_type &particles) {
auto dcpse_temp = (Dcpse<particles_type::dims, particles_type> *) dcpse;
dcpse_temp->createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->template createNormalParticles<NORMAL_ID>(particles);
dcpse_temp->initializeUpdate(particles);
dcpse_temp->accumulateAndDeleteNormalParticles(particles);
}
};
#endif //Eigen
#endif //OPENFPM_PDATA_DCPSE_SURFACE_OP_HPP
src/DCPSE/DCPSE_op/EqnsStruct.hpp
View file @ 1371f3eb
......@@ -1228,6 +1228,7 @@ struct equations3d3EPz_gpu {
typedef umfpack_solver<double> solver_type;
};
struct equations3d1Pxy_gpu {
//! dimensionaly of the equation ( 3D problem ...)
static const unsigned int dims = 3;
......
src/DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cpp
View file @ 1371f3eb
......@@ -1011,7 +1011,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.solve_with_solver(solver,sol);
Solver.reset_b();
Solver.impose_b(bulk, prop_id<1>());
Solver.impose_b(up_p, prop_id<1>());
......@@ -1040,6 +1039,170 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
//domain.write("Neumann");
}
BOOST_AUTO_TEST_CASE(dcpse_poisson_Neumann2d) {
const size_t sz[2] = {31,31};
Box<2, double> box({0, 0}, {1.0, 1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
double rCut = 3.1 * spacing;
Ghost<2, double> ghost(spacing * 3.1);
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist<2, double, aggregate<double[2],double[2],double[2],double[2],double[2]>> domain(0, box, bc, ghost);
//Init_DCPSE(domain)
BOOST_TEST_MESSAGE("Init domain...");
auto it = domain.getGridIterator(sz);
while (it.isNext()) {
domain.add();
auto key = it.get();
double x = key.get(0) * it.getSpacing(0);
domain.getLastPos()[0] = x;
double y = key.get(1) * it.getSpacing(1);
domain.getLastPos()[1] = y;
++it;
}
BOOST_TEST_MESSAGE("Sync domain across processors...");
domain.map();
domain.ghost_get<0>();
Derivative_x Dx(domain, 2, rCut,1.9,support_options::N_PARTICLES);
Derivative_y Dy(domain, 2, rCut,1.9,support_options::N_PARTICLES);
Laplacian Lap(domain, 2, rCut, 1.9,support_options::N_PARTICLES);
petsc_solver<double> solver;
solver.setRestart(500);
solver.setSolver(KSPGMRES);
solver.setPreconditioner(PCSVD);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
openfpm::vector<aggregate<int>> dw_p;
openfpm::vector<aggregate<int>> l_p;
openfpm::vector<aggregate<int>> r_p;
auto v = getV<0>(domain);
auto RHS=getV<1>(domain);
auto sol = getV<2>(domain);
auto anasol = getV<3>(domain);
auto err = getV<4>(domain);
// Here fill me
Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
openfpm::vector<Box<2, double>> boxes;
boxes.add(up);
boxes.add(down);
boxes.add(left);
boxes.add(right);
// Create a writer and write
VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
vtk_box.add(boxes);
//vtk_box.write("vtk_box.vtk");
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = domain.getPos(p);
//domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
domain.getProp<1>(p)[0] = sin(5*xp.get(0));
domain.getProp<1>(p)[1] = sin(5*xp.get(0));
} else if (down.isInside(xp) == true) {
dw_p.add();
dw_p.last().get<0>() = p.getKey();
domain.getProp<1>(p)[0] = sin(5*xp.get(0));
domain.getProp<1>(p)[1] = sin(5*xp.get(0));
} else if (left.isInside(xp) == true) {
l_p.add();
l_p.last().get<0>() = p.getKey();
domain.getProp<1>(p)[0] = sin(5*xp.get(0));
domain.getProp<1>(p)[1] = sin(5*xp.get(0));
} else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
domain.getProp<1>(p)[0] = sin(5*xp.get(0));
domain.getProp<1>(p)[1] = sin(5*xp.get(0));
} else {
bulk.add();
bulk.last().get<0>() = p.getKey();
domain.getProp<1>(p)[0] = -10*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
domain.getProp<1>(p)[1] = -10*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
}
++it2;
}
DCPSE_scheme<equations2d2,decltype(domain)> Solver(domain,options_solver::LAGRANGE_MULTIPLIER);
eq_id vx,vy;
vx.setId(0);
vy.setId(1);
auto Poisson0 = -Lap(v[0]);
auto D_x0 = Dx(v[0]);
auto D_y0 = Dy(v[0]);
auto Poisson1 = -Lap(v[1]);
auto D_x1 = Dx(v[1]);
auto D_y1 = Dy(v[1]);
Solver.impose(Poisson0, bulk, RHS[0],vx);
Solver.impose(Poisson1, bulk, RHS[1],vy);
Solver.impose(D_y0, up_p, RHS[0],vx);
Solver.impose(-D_y0, dw_p, RHS[0],vx);
Solver.impose(-D_x0, l_p, RHS[0],vx);
Solver.impose(D_x0, r_p, RHS[0],vx);
Solver.impose(D_y1, up_p, RHS[1],vy);
Solver.impose(-D_y1, dw_p, RHS[1],vy);
Solver.impose(-D_x1, l_p, RHS[1],vy);
Solver.impose(D_x1, r_p, RHS[1],vy);
Solver.solve_with_solver(solver,sol[0],sol[1]);
// Solver.solve(sol);
domain.ghost_get<2>();
anasol[0]=-Lap(sol[0]);
anasol[1]=-Lap(sol[1]);
double worst1 = 0.0,worst2 = 0.0;
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
if (fabs(domain.getProp<3>(p)[0]- domain.getProp<1>(p)[0]) >= worst1) {
worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<1>(p)[0]);
}
if (fabs(domain.getProp<3>(p)[1]- domain.getProp<1>(p)[1]) >= worst2) {
worst2 = fabs(domain.getProp<3>(p)[1] - domain.getProp<1>(p)[1]);
}
domain.getProp<4>(p)[0] = fabs(domain.getProp<1>(p)[0] - domain.getProp<3>(p)[0]);
domain.getProp<4>(p)[1] = fabs(domain.getProp<1>(p)[1] - domain.getProp<3>(p)[1]);
}
//Auto Error
BOOST_REQUIRE(worst1 < 1.0);
BOOST_REQUIRE(worst2 < 1.0);
domain.write("Neumann2d");
}
BOOST_AUTO_TEST_CASE(dcpse_slice_solver) {
// int rank;
......
src/DCPSE/DCPSE_op/tests/DCPSE_op_Solver_test.cu
View file @ 1371f3eb
......@@ -17,7 +17,6 @@
#include <iostream>
#include "../DCPSE_op.hpp"
#include "../DCPSE_Solver.hpp"
#include "../DCPSE_Solver.cuh"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
#include "../EqnsStruct.hpp"
......@@ -154,7 +153,7 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
Ghost<2, double> ghost(spacing * 3);
double rCut = 3.1 * spacing;
double rCut = 2.0 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist_gpu<2, double, aggregate<double,double,double,double>> domain(0, box, bc, ghost);
......@@ -180,7 +179,7 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
domain.map();
domain.ghost_get<0>();
Laplacian_gpu Lap(domain, 2, rCut);
Laplacian_gpu Lap(domain, 2, rCut, 2,support_options::N_PARTICLES);
DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver( domain);
......@@ -415,10 +414,9 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
domain.map();
domain.ghost_get<0>();
Derivative_x_gpu Dx(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
Derivative_y_gpu Dy(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
Laplacian_gpu Lap(domain, 2, rCut / 3.0 ,1.9/*,support_options::RADIUS*/);
Derivative_x_gpu Dx(domain, 2, rCut/3.0 ,1.9,support_options::N_PARTICLES);
Derivative_y_gpu Dy(domain, 2, rCut/3.0,1.9,support_options::N_PARTICLES);
Laplacian_gpu Lap(domain, 2, rCut/3.0 ,1.9,support_options::N_PARTICLES);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
......@@ -552,7 +550,7 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
BOOST_AUTO_TEST_CASE(dcpse_poisson_Dirichlet_anal) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
const size_t sz[2] = {200,200};
const size_t sz[2] = {81,81};
Box<2, double> box({0, 0}, {1, 1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
......@@ -860,7 +858,7 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
Ghost<2, double> ghost(spacing * 3);
double rCut = 3.1 * spacing;
double rCut = 2.0 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist_gpu<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> domain(0, box, bc, ghost);
......@@ -886,8 +884,8 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
domain.map();
domain.ghost_get<0>();
Derivative_y_gpu Dy(domain, 2, rCut);
Laplacian_gpu Lap(domain, 2, rCut);
Derivative_y_gpu Dy(domain, 2, rCut,2,support_options::N_PARTICLES);
Laplacian_gpu Lap(domain, 2, rCut, 3,support_options::N_PARTICLES);
DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver(domain);
......@@ -1033,9 +1031,10 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
domain.map();
domain.ghost_get<0>();
Derivative_x_gpu Dx(domain, 2, rCut);
Derivative_y_gpu Dy(domain, 2, rCut);
Laplacian_gpu Lap(domain, 2, rCut);
Derivative_x_gpu Dx(domain, 2, rCut,1.9,support_options::N_PARTICLES);
Derivative_y_gpu Dy(domain, 2, rCut,1.9,support_options::N_PARTICLES);
Laplacian_gpu Lap(domain, 2, rCut, 1.9,support_options::N_PARTICLES);
petsc_solver<double> solver;
solver.setRestart(500);
solver.setSolver(KSPGMRES);
......@@ -1175,10 +1174,9 @@ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
domain.map();
domain.ghost_get<0>();
Derivative_x_gpu Dx(domain, 2, rCut,1.9,support_options::RADIUS);
Derivative_y_gpu Dy(domain, 2, rCut,1.9,support_options::RADIUS);
Laplacian_gpu Lap(domain, 2, rCut,1.9,support_options::RADIUS);
Derivative_x Dx(domain, 2, rCut,1.9,support_options::RADIUS);
Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
Laplacian Lap(domain, 2, rCut,1.9,support_options::RADIUS);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
......
src/DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cu 0 → 100644
View file @ 1371f3eb
/*
* DCPSE_op_test.cpp
*
* Created on: May 15, 2020
* Author: Abhinav Singh
*
*/
#include "config.h"
#ifdef HAVE_EIGEN
#ifdef HAVE_PETSC
#define BOOST_TEST_DYN_LINK
#include "util/util_debug.hpp"
#include <boost/test/unit_test.hpp>
#include <iostream>
#include "../DCPSE_op.hpp"
#include "../DCPSE_Solver.hpp"
#include "../DCPSE_Solver.cuh"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
#include "../EqnsStruct.hpp"
//template<typename T>
//struct Debug;
#if 0
BOOST_AUTO_TEST_SUITE(dcpse_op_subset_suite_tests_cu)
BOOST_AUTO_TEST_CASE(dcpse_op_subset_tests) {
size_t edgeSemiSize = 40;
const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
Box<2, double> box({0, 0}, {1.0,1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
double rCut = 3.9 * spacing[0];
int ord = 2;
double sampling_factor = 4.0;
Ghost<2, double> ghost(rCut);
BOOST_TEST_MESSAGE("Init vector_dist...");
double sigma2 = spacing[0] * spacing[1] / (2 * 4);
vector_dist_ws_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles(0, box,
bc,
ghost);
//Init_DCPSE(Particles)
BOOST_TEST_MESSAGE("Init Particles...");
std::mt19937 rng{6666666};
std::normal_distribution<> gaussian{0, sigma2};
// openfpm::vector<aggregate<int>> bulk;
// openfpm::vector<aggregate<int>> boundary;
auto it = Particles.getGridIterator(sz);
size_t pointId = 0;
size_t counter = 0;
double minNormOne = 999;
while (it.isNext())
{
Particles.add();
auto key = it.get();
mem_id k0 = key.get(0);
double x = k0 * spacing[0];
Particles.getLastPos()[0] = x;//+ gaussian(rng);
mem_id k1 = key.get(1);
double y = k1 * spacing[1];
Particles.getLastPos()[1] = y;//+gaussian(rng);
// Here fill the function value
Particles.template getLastProp<0>() = sin(Particles.getLastPos()[0]) + sin(Particles.getLastPos()[1]);
if (k0 != 0 && k1 != 0 && k0 != sz[0] -1 && k1 != sz[1] - 1)
{
// bulk.add();
// bulk.template get<0>(bulk.size()-1) = Particles.size_local() - 1;
Particles.getLastSubset(0);
}
else
{
// boundary.add();
// boundary.template get<0>(boundary.size()-1) = Particles.size_local() - 1;
Particles.getLastSubset(1);
}
++counter;
++it;
}
BOOST_TEST_MESSAGE("Sync Particles across processors...");
Particles.map();
Particles.ghost_get<0>();
auto git = Particles.getGhostIterator();
while (git.isNext())
{
auto p = git.get();
Particles.template getProp<0>(p) = std::numeric_limits<double>::quiet_NaN();
++git;
}
vector_dist_subset_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles_bulk(Particles,0);
vector_dist_subset_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles_boundary(Particles,1);
auto & boundary = Particles_boundary.getIds();
// move particles
auto P = getV<0>(Particles);
auto Out = getV<1>(Particles);
auto Pb = getV<2>(Particles);
auto Out_V = getV<3>(Particles);
auto P_bulk = getV<2>(Particles_bulk);
auto Out_bulk = getV<1>(Particles_bulk);
auto Out_V_bulk = getV<3>(Particles_bulk);
Out=10;
P_bulk = 5;
P_bulk=Pb+Out;
// Particles.write("Test_output_subset");
// Create the subset
/* Derivative_x Dx(Particles, 2, rCut);
Derivative_y Dy(Particles, 2, rCut);
Derivative_x Dx_bulk(Particles_bulk, 2, rCut);
*/
Derivative_x_gpu Dx_bulk(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_y_gpu Dy_bulk(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
Out_bulk = Dx_bulk(P);
Out_V_bulk[0] = P + Dx_bulk(P);
Out_V_bulk[1] = Out_V[0] +Dy_bulk(P);
// Check
bool is_nan = false;
auto & v_cl = create_vcluster();
if (v_cl.size() > 1)
{
auto it2 = Particles_bulk.getDomainIterator();
while (it2.isNext())
{
auto p = it2.get();
/* BOOST_REQUIRE_EQUAL(Particles_bulk.getProp<2>(p),15.0);
BOOST_REQUIRE(fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.005 );
BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[0] - Particles_bulk.getProp<0>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.001 );
BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[1] - Particles_bulk.getProp<3>(p)[0] - cos(Particles_bulk.getPos(p)[1])) < 0.001 );*/
is_nan |= std::isnan(Particles_bulk.template getProp<1>(p));
// Particles_bulk.template getProp<0>(p) = fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0]));
++it2;
}
BOOST_REQUIRE_EQUAL(is_nan,true);
}
// P_bulk = Dx_bulk(P_bulk); <------------ Incorrect produce error message
// P = Dx_bulk(P); <------- Incorrect produce overflow
Particles.ghost_get<0>();
for (int i = 0 ; i < boundary.size() ; i++)
{
Particles.template getProp<0>(boundary.template get<0>(i)) = std::numeric_limits<double>::quiet_NaN();
}
Particles.ghost_get<0>();
Out_bulk = Dx_bulk(P);
Out_V_bulk[0] = P + Dx_bulk(P);
Out_V_bulk[1] = Out_V[0] +Dy_bulk(P);
auto it2 = Particles_bulk.getDomainIterator();
while (it2.isNext())
{
auto p = it2.get();
BOOST_REQUIRE_EQUAL(Particles_bulk.getProp<2>(p),15.0);
BOOST_REQUIRE(fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.005 );
BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[0] - Particles_bulk.getProp<0>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.001 );
BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[1] - Particles_bulk.getProp<3>(p)[0] - cos(Particles_bulk.getPos(p)[1])) < 0.001 );
++it2;
}
}
BOOST_AUTO_TEST_CASE(dcpse_op_subset_PC_lid) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
constexpr int x = 0;
constexpr int y = 1;
size_t edgeSemiSize = 20;
const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
Box<2, double> box({0, 0}, {1.0,1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
double rCut = 3.9 * spacing[0];
int ord = 2;
double sampling_factor = 4.0;
Ghost<2, double> ghost(rCut);
BOOST_TEST_MESSAGE("Init vector_dist...");
double sigma2 = spacing[0] * spacing[1] / (2 * 4);
auto &v_cl = create_vcluster();
typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double> particle_type;
vector_dist_ws_gpu<2, double, particle_type> Particles(0, box,bc,ghost);
//Init_DCPSE(Particles)
BOOST_TEST_MESSAGE("Init Particles...");
// openfpm::vector<aggregate<int>> bulk;
// openfpm::vector<aggregate<int>> boundary;
auto it = Particles.getGridIterator(sz);
while (it.isNext())
{
Particles.add();
auto key = it.get();
mem_id k0 = key.get(0);
double xp0 = k0 * spacing[0];
Particles.getLastPos()[0] = xp0;
mem_id k1 = key.get(1);
double yp0 = k1 * spacing[1];
Particles.getLastPos()[1] = yp0;
++it;
}
BOOST_TEST_MESSAGE("Sync Particles across processors...");
Particles.map();
Particles.ghost_get<0>();
auto it2 = Particles.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
if (xp[0] != 0 && xp[1] != 0 && xp[0] != 1.0 && xp[1] != 1.0) {
// bulk.add();
// bulk.last().get<0>() = p.getKey();
Particles.setSubset(p,0);
Particles.getProp<3>(p)[x] = 3.0;
Particles.getProp<3>(p)[y] = 3.0;
} else {
// boundary.add();
// boundary.last().get<0>() = p.getKey();
Particles.setSubset(p,1);
Particles.getProp<3>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
Particles.getProp<3>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
}
Particles.getProp<6>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
Particles.getProp<6>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
Particles.getProp<7>(p) = xp[0]+xp[1]-1.0;
++it2;
}
vector_dist_subset_gpu<2, double, particle_type> Particles_bulk(Particles,0);
vector_dist_subset_gpu<2, double, particle_type> Particles_boundary(Particles,1);
auto & bulk = Particles_bulk.getIds();
auto & boundary = Particles_boundary.getIds();
auto P = getV<0>(Particles);
auto V = getV<1>(Particles);
auto RHS = getV<2>(Particles);
auto dV = getV<3>(Particles);
auto div = getV<4>(Particles);
auto V_star = getV<5>(Particles);
auto P_bulk = getV<0>(Particles_bulk);
auto RHS_bulk =getV<2>(Particles_bulk);
P_bulk = 0;
Derivative_x_gpu Dx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_xx_gpu Dxx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_yy_gpu Dyy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_y_gpu Dy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_x_gpu Bulk_Dx(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_y_gpu Bulk_Dy(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
int n = 0, nmax = 5, ctr = 0, errctr=1, Vreset = 0;
double V_err=1;
if (Vreset == 1) {
P_bulk = 0;
P = 0;
Vreset = 0;
}
P=0;
eq_id vx,vy;
vx.setId(0);
vy.setId(1);
double sum, sum1, sum_k,V_err_eps=1e-3,V_err_old;
auto Stokes1=Dxx(V[x])+Dyy(V[x]);
auto Stokes2=Dxx(V[y])+Dyy(V[y]);
petsc_solver<double> solverPetsc;
//solverPetsc.setSolver(KSPGMRES);
//solverPetsc.setRestart(250);
//solverPetsc.setPreconditioner(PCJACOBI);
V_star=0;
RHS[x] = dV[x];
RHS[y] = dV[y];
while (V_err >= V_err_eps && n <= nmax) {
Particles.ghost_get<0>(SKIP_LABELLING);
RHS_bulk[x] = dV[x] + Bulk_Dx(P);
RHS_bulk[y] = dV[y] + Bulk_Dy(P);
DCPSE_scheme_gpu<equations2d2_gpu, decltype(Particles)> Solver(Particles);
Solver.impose(Stokes1, bulk, RHS[0], vx);
Solver.impose(Stokes2, bulk, RHS[1], vy);
Solver.impose(V[x], boundary, RHS[0], vx);
Solver.impose(V[y], boundary, RHS[1], vy);
/*auto A=Solver.getA(options_solver::STANDARD);
//A.getMatrixTriplets().save("Tripletes");
A.write("Mat_lid");*/
Solver.solve_with_solver(solverPetsc, V[x], V[y]);
Particles.ghost_get<1>(SKIP_LABELLING);
div = -(Dx(V[x]) + Dy(V[y]));
P_bulk = P + div;
sum = 0;
sum1 = 0;
for (int j = 0; j < bulk.size(); j++) {
auto p = bulk.get<0>(j);
sum += (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) *
(Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) +
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]) *
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]);
sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1];
}
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_star = V;
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
V_err_old = V_err;
V_err = sum / sum1;
if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
errctr++;
//alpha_P -= 0.1;
} else {
errctr = 0;
}
if (n > 3) {
if (errctr > 3) {
std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
Vreset = 1;
break;
} else {
Vreset = 0;
}
}
n++;
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << std::endl;
}
}
double worst1 = 0.0;
double worst2 = 0.0;
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
if (fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]) >= worst1) {
worst1 = fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]);
}
}
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
if (fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]) >= worst2) {
worst2 = fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]);
}
}
//Particles.deleteGhost();
//Particles.write("PC_subset_lid");
std::cout << "Maximum Analytic Error in Vx: " << worst1 << std::endl;
std::cout << "Maximum Analytic Error in Vy: " << worst2 << std::endl;
BOOST_REQUIRE(worst1 < 0.03);
BOOST_REQUIRE(worst2 < 0.03);
}
BOOST_AUTO_TEST_CASE(dcpse_op_subset_PC_lid2) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
constexpr int x = 0;
constexpr int y = 1;
size_t edgeSemiSize = 20;
const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
Box<2, double> box({0, 0}, {1.0,1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
double rCut = 3.9 * spacing[0];
int ord = 2;
double sampling_factor = 4.0;
Ghost<2, double> ghost(rCut);
BOOST_TEST_MESSAGE("Init vector_dist...");
auto &v_cl = create_vcluster();
vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double>> Particles(0, box,
bc,
ghost);
vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double>> Particles_subset(Particles.getDecomposition(), 0);
//Init_DCPSE(Particles)
BOOST_TEST_MESSAGE("Init Particles...");
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> boundary;
auto it = Particles.getGridIterator(sz);
size_t pointId = 0;
double minNormOne = 999;
while (it.isNext())
{
Particles.add();
auto key = it.get();
mem_id k0 = key.get(0);
double xp0 = k0 * spacing[0];
Particles.getLastPos()[0] = xp0;
mem_id k1 = key.get(1);
double yp0 = k1 * spacing[1];
Particles.getLastPos()[1] = yp0;
++it;
}
BOOST_TEST_MESSAGE("Sync Particles across processors...");
Particles.map();
Particles.ghost_get<0>();
auto it2 = Particles.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
if (xp[0] != 0 && xp[1] != 0 && xp[0] != 1.0 && xp[1] != 1.0) {
bulk.add();
bulk.last().get<0>() = p.getKey();
Particles.getProp<3>(p)[x] = 3.0;
Particles.getProp<3>(p)[y] = 3.0;
} else {
boundary.add();
boundary.last().get<0>() = p.getKey();
Particles.getProp<3>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
Particles.getProp<3>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
}
Particles.getProp<6>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
Particles.getProp<6>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
Particles.getProp<7>(p) = xp[0]+xp[1]-1.0;
++it2;
}
for (int i = 0; i < bulk.size(); i++) {
Particles_subset.add();
Particles_subset.getLastPos()[0] = Particles.getPos(bulk.template get<0>(i))[0];
Particles_subset.getLastPos()[1] = Particles.getPos(bulk.template get<0>(i))[1];
}
Particles_subset.map();
Particles_subset.ghost_get<0>();
auto P = getV<0>(Particles);
auto V = getV<1>(Particles);
auto RHS = getV<2>(Particles);
auto dV = getV<3>(Particles);
auto div = getV<4>(Particles);
auto V_star = getV<5>(Particles);
auto P_bulk = getV<0>(Particles_subset);
auto Grad_bulk= getV<2>(Particles_subset);
P_bulk = 0;
Derivative_x Dx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_x Bulk_Dx(Particles_subset, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_xx Dxx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_yy Dyy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_y Dy(Particles, 2, rCut,sampling_factor, support_options::RADIUS),Bulk_Dy(Particles_subset, 2, rCut,sampling_factor, support_options::RADIUS);;
int n = 0, nmax = 5, ctr = 0, errctr=0, Vreset = 0;
double V_err=1;
if (Vreset == 1) {
P_bulk = 0;
P = 0;
Vreset = 0;
}
P=0;
eq_id vx,vy;
vx.setId(0);
vy.setId(1);
double sum, sum1, sum_k,V_err_eps=1e-3,V_err_old;
auto Stokes1=Dxx(V[x])+Dyy(V[x]);
auto Stokes2=Dxx(V[y])+Dyy(V[y]);
petsc_solver<double> solverPetsc;
//solverPetsc.setSolver(KSPGMRES);
//solverPetsc.setRestart(250);
//solverPetsc.setPreconditioner(PCJACOBI);
V_star=0;
while (V_err >= V_err_eps && n <= nmax) {
RHS[x] = dV[x];
RHS[y] = dV[y];
Particles_subset.ghost_get<0>(SKIP_LABELLING);
Grad_bulk[x] = Bulk_Dx(P_bulk);
Grad_bulk[y] = Bulk_Dy(P_bulk);
for (int i = 0; i < bulk.size(); i++) {
Particles.template getProp<2>(bulk.template get<0>(i))[x] += Particles_subset.getProp<2>(i)[x];
Particles.template getProp<2>(bulk.template get<0>(i))[y] += Particles_subset.getProp<2>(i)[y];
}
DCPSE_scheme<equations2d2_gpu, decltype(Particles)> Solver(Particles);
Solver.impose(Stokes1, bulk, RHS[0], vx);
Solver.impose(Stokes2, bulk, RHS[1], vy);
Solver.impose(V[x], boundary, RHS[0], vx);
Solver.impose(V[y], boundary, RHS[1], vy);
Solver.solve_with_solver(solverPetsc, V[x], V[y]);
Particles.ghost_get<1>(SKIP_LABELLING);
div = -(Dx(V[x]) + Dy(V[y]));
P = P + div;
for (int i = 0; i < bulk.size(); i++) {
Particles_subset.getProp<0>(i) = Particles.template getProp<0>(bulk.template get<0>(i));
}
sum = 0;
sum1 = 0;
for (int j = 0; j < bulk.size(); j++) {
auto p = bulk.get<0>(j);
sum += (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) *
(Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) +
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]) *
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]);
sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1];
}
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_star=V;
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
V_err_old = V_err;
V_err = sum / sum1;
if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
errctr++;
//alpha_P -= 0.1;
} else {
errctr = 0;
}
if (n > 3) {
if (errctr > 3) {
std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
Vreset = 1;
break;
} else {
Vreset = 0;
}
}
n++;
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << std::endl;
}
}
double worst1 = 0.0;
double worst2 = 0.0;
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
if (fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]) >= worst1) {
worst1 = fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]);
}
}
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
if (fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]) >= worst2) {
worst2 = fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]);
}
}
std::cout << "Maximum Analytic Error in slice x: " << worst1 << std::endl;
std::cout << "Maximum Analytic Error in slice y: " << worst2 << std::endl;
BOOST_REQUIRE(worst1 < 0.03);
BOOST_REQUIRE(worst2 < 0.03);
//Particles.write("PC_subset_lid2");
}
BOOST_AUTO_TEST_SUITE_END()
#endif
#endif
#endif //if 0
\ No newline at end of file
src/DCPSE/Dcpse.cuh
View file @ 1371f3eb
......@@ -484,9 +484,10 @@ public:
private:
template <typename U>
void initializeAdaptive(vector_type &particles,
unsigned int convergenceOrder,
double rCut) {
U rCut) {
// Still need to be tested
#ifdef SE_CLASS1
this->update_ctr=particles.getMapCtr();
......@@ -505,77 +506,7 @@ private:
if (!isSharedSupport) {
SupportBuilder<vector_type,vector_type>
supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
unsigned int requiredSupportSize = monomialBasis.size();
// need to resize supportKeys1D to yet unknown supportKeysTotalN
// add() takes too long
openfpm::vector<openfpm::vector<size_t>> tempSupportKeys(supportRefs.size());
auto it = particles.getDomainIterator();
while (it.isNext()) {
auto key_o = particles.getOriginKey(it.get());
subsetKeyPid.get(key_o.getKey()) = it.get().getKey();
Support support = supportBuilder.getSupport(it, requiredSupportSize, opt);
supportRefs.get(key_o.getKey()) = key_o.getKey();
tempSupportKeys.get(key_o.getKey()) = support.getKeys();
kerOffsets.get(key_o.getKey()) = supportKeysTotalN;
if (maxSupportSize < support.size())
maxSupportSize = support.size();
supportKeysTotalN += support.size();
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
// Vandermonde matrix computation
Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
vandermonde(support, monomialBasis, particles, particles,HOverEpsilon);
vandermonde.getMatrix(V);
T condV = conditionNumber(V, condVTOL);
T eps = vandermonde.getEps();
if (condV > condVTOL) {
requiredSupportSize *= 2;
std::cout << "INFO: Increasing, requiredSupportSize = " << requiredSupportSize << std::endl; // debug
continue;
} else requiredSupportSize = monomialBasis.size();
++it;
}
kerOffsets.get(supportRefs.size()) = supportKeysTotalN;
supportKeys1D.resize(supportKeysTotalN);
size_t offset = 0;
for (size_t i = 0; i < tempSupportKeys.size(); ++i)
for (size_t j = 0; j < tempSupportKeys.get(i).size(); ++j, ++offset)
supportKeys1D.get(offset) = tempSupportKeys.get(i).get(j);
}
kerOffsets.hostToDevice(); supportKeys1D.hostToDevice();
assembleLocalMatrices(cublasDgetrfBatched, cublasDtrsmBatched);
}
void initializeAdaptive(vector_type &particles,
unsigned int convergenceOrder,
float rCut) {
// Still need to be tested
#ifdef SE_CLASS1
this->update_ctr=particles.getMapCtr();
#endif
if (!isSharedSupport) {
subsetKeyPid.resize(particles.size_local_orig());
supportRefs.resize(particles.size_local());
}
localEps.resize(particles.size_local());
localEpsInvPow.resize(particles.size_local());
kerOffsets.resize(particles.size_local()+1);
const T condVTOL = 1e2;
if (!isSharedSupport) {
SupportBuilder<vector_type,vector_type>
supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
unsigned int requiredSupportSize = monomialBasis.size();
// need to resize supportKeys1D to yet unknown supportKeysTotalN
// add() takes too long
......@@ -598,13 +529,11 @@ private:
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
// Vandermonde matrix computation
Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
vandermonde(support, monomialBasis, particles, particles, HOverEpsilon);
vandermonde(support, monomialBasis, particles, particles,HOverEpsilon);
vandermonde.getMatrix(V);
T condV = conditionNumber(V, condVTOL);
T eps = vandermonde.getEps();
if (condV > condVTOL) {
requiredSupportSize *= 2;
std::cout << "INFO: Increasing, requiredSupportSize = " << requiredSupportSize << std::endl; // debug
......@@ -624,13 +553,14 @@ private:
}
kerOffsets.hostToDevice(); supportKeys1D.hostToDevice();
assembleLocalMatrices(cublasSgetrfBatched, cublasStrsmBatched);
assembleLocalMatrices_t(rCut);
}
template <typename U>
void initializeStaticSize(vector_type &particles,
unsigned int convergenceOrder,
double rCut,
double supportSizeFactor) {
U rCut,
U supportSizeFactor) {
#ifdef SE_CLASS1
this->update_ctr=particles.getMapCtr();
#endif
......@@ -697,82 +627,19 @@ std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution
std::chrono::duration<double> time_span2 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << "Support building took " << time_span2.count() * 1000. << " milliseconds." << std::endl;
assembleLocalMatrices(cublasDgetrfBatched, cublasDtrsmBatched);
assembleLocalMatrices_t(rCut);
}
// ad hoc solution to template specialization for float/double
void initializeStaticSize(vector_type &particles,
unsigned int convergenceOrder,
float rCut,
float supportSizeFactor) {
#ifdef SE_CLASS1
this->update_ctr=particles.getMapCtr();
#endif
this->rCut=rCut;
this->supportSizeFactor=supportSizeFactor;
this->convergenceOrder=convergenceOrder;
// Cublas subroutine selector: float or double
// rCut is needed to overcome limitation of nested template class specialization
template <typename U> void assembleLocalMatrices_t(U rCut) {
throw std::invalid_argument("DCPSE operator error: CUBLAS supports only float or double"); }
if (!isSharedSupport) {
subsetKeyPid.resize(particles.size_local_orig());
supportRefs.resize(particles.size_local());
}
localEps.resize(particles.size_local());
localEpsInvPow.resize(particles.size_local());
void assembleLocalMatrices_t(float rCut) {
assembleLocalMatrices(cublasSgetrfBatched, cublasStrsmBatched); }
std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
auto it = particles.getDomainIterator();
if (opt==support_options::RADIUS) {
if (!isSharedSupport) {
while (it.isNext()) {
auto key_o = it.get(); subsetKeyPid.get(particles.getOriginKey(key_o).getKey()) = key_o.getKey();
supportRefs.get(key_o.getKey()) = key_o.getKey();
++it;
}
SupportBuilderGPU<vector_type> supportBuilder(particles, rCut);
supportBuilder.getSupport(supportRefs.size(), kerOffsets, supportKeys1D, maxSupportSize, supportKeysTotalN);
}
} else {
if (!isSharedSupport){
openfpm::vector<openfpm::vector<size_t>> tempSupportKeys(supportRefs.size());
size_t requiredSupportSize = monomialBasis.size() * supportSizeFactor;
// need to resize supportKeys1D to yet unknown supportKeysTotalN
// add() takes too long
SupportBuilder<vector_type,vector_type> supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
kerOffsets.resize(supportRefs.size()+1);
while (it.isNext()) {
auto key_o = it.get(); subsetKeyPid.get(particles.getOriginKey(key_o).getKey()) = key_o.getKey();
Support support = supportBuilder.getSupport(it, requiredSupportSize, opt);
supportRefs.get(key_o.getKey()) = key_o.getKey();
tempSupportKeys.get(key_o.getKey()) = support.getKeys();
kerOffsets.get(key_o.getKey()) = supportKeysTotalN;
if (maxSupportSize < support.size()) maxSupportSize = support.size();
supportKeysTotalN += support.size();
++it;
}
kerOffsets.get(supportRefs.size()) = supportKeysTotalN;
supportKeys1D.resize(supportKeysTotalN);
size_t offset = 0;
for (size_t i = 0; i < tempSupportKeys.size(); ++i)
for (size_t j = 0; j < tempSupportKeys.get(i).size(); ++j, ++offset)
supportKeys1D.get(offset) = tempSupportKeys.get(i).get(j);
}
kerOffsets.hostToDevice(); supportKeys1D.hostToDevice();
}
std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> time_span2 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << "Support building took " << time_span2.count() * 1000. << " milliseconds." << std::endl;
assembleLocalMatrices(cublasSgetrfBatched, cublasStrsmBatched);
}
void assembleLocalMatrices_t(double rCut) {
assembleLocalMatrices(cublasDgetrfBatched, cublasDtrsmBatched); }
template<typename cublasLUDec_type, typename cublasTriangSolve_type>
void assembleLocalMatrices(cublasLUDec_type cublasLUDecFunc, cublasTriangSolve_type cublasTriangSolveFunc) {
......
src/DCPSE/Dcpse.hpp
View file @ 1371f3eb
......@@ -878,7 +878,7 @@ private:
localEpsInvPow.resize(particlesTo.size_local_orig());
kerOffsets.resize(particlesTo.size_local_orig());
kerOffsets.fill(-1);
T avgSpacingGlobal=0,maxSpacingGlobal=0,minSpacingGlobal=std::numeric_limits<T>::max();
T avgSpacingGlobal=0,avgSpacingGlobal2=0,maxSpacingGlobal=0,minSpacingGlobal=std::numeric_limits<T>::max();
size_t Counter=0;
auto it = particlesTo.getDomainIterator();
while (it.isNext()) {
......@@ -900,6 +900,7 @@ private:
T eps = vandermonde.getEps();
avgSpacingGlobal+=eps;
T tSpacing = vandermonde.getMinSpacing();
avgSpacingGlobal2+=tSpacing;
if(tSpacing>maxSpacingGlobal)
{
maxSpacingGlobal=tSpacing;
......@@ -948,12 +949,13 @@ private:
}
v_cl.sum(avgSpacingGlobal);
v_cl.sum(avgSpacingGlobal2);
v_cl.max(maxSpacingGlobal);
v_cl.min(minSpacingGlobal);
v_cl.sum(Counter);
v_cl.execute();
if(v_cl.rank()==0)
{std::cout<<"DCPSE Operator Construction Complete. The global avg spacing in the support <h> is: "<<HOverEpsilon*avgSpacingGlobal/(T(Counter))<<" (c="<<HOverEpsilon<<"). Min Spacing Range=["<<minSpacingGlobal<<","<<maxSpacingGlobal<<"]."<<std::endl;}
{std::cout<<"DCPSE Operator Construction Complete. The global avg spacing in the support <h> is: "<<HOverEpsilon*avgSpacingGlobal/(T(Counter))<<" (c="<<HOverEpsilon<<"). Avg:"<<avgSpacingGlobal2/(T(Counter))<<" Range:["<<minSpacingGlobal<<","<<maxSpacingGlobal<<"]."<<std::endl;}
}
T computeKernel(Point<dim, T> x, EMatrix<T, Eigen::Dynamic, 1> & a) const {
......
src/DCPSE/SupportBuilder.hpp
View file @ 1371f3eb
......@@ -209,6 +209,9 @@ private:
MinSpacing = rp.get(i).dist;
}
}
#ifdef SE_CLASS1
assert(MinSpacing !=0 && "You have multiple particles on the same position.");
#endif
for (int i = 0; i < rp.size(); i++) {
if (rp.get(i).dist < AdapFac * MinSpacing) {
points.push_back(rp.get(i).offset);
......@@ -221,6 +224,7 @@ private:
points.push_back(rp.get(i).offset);
}
}
//MinSpacing=MinSpacing/requiredSupportSize
return points;
}
......
src/Matrix/SparseMatrix_petsc.hpp
View file @ 1371f3eb
......@@ -372,6 +372,17 @@ public:
return 0;
}
/*! \brief Return the state of matrix
*
* Returns a bool flag that indicated whether the matrix
* has already been filled via MatSetValues
*
*/
bool isMatrixFilled()
{
return m_created;
}
/* Write matrix on vtk
*
* \param out file to write into
......
src/OdeIntegrators/OdeIntegrators.hpp
View file @ 1371f3eb
......@@ -30,7 +30,58 @@ namespace boost{
#include <boost/numeric/odeint.hpp>
#include "Operators/Vector/vector_dist_operators.hpp"
#include "FiniteDifference/FD_expressions.hpp"
#include "OdeIntegrators/boost_vector_algebra_ofp.hpp"
#include "OdeIntegrators/vector_algebra_ofp.hpp"
#ifdef __NVCC__
#include "OdeIntegrators/vector_algebra_ofp_gpu.hpp"
/*! \brief A 1d Odeint and Openfpm compatible structure.
*
* Use the method this.data.get<d>() to refer to property of all the particles in the dimension d.
*
* d starts with 0.
*
*/
struct state_type_1d_ofp_ker{
state_type_1d_ofp_ker(){
}
typedef decltype(std::declval<texp_v_gpu<double>>().getVector().toKernel()) state_kernel;
typedef size_t size_type;
typedef int is_state_vector;
aggregate<state_kernel> data;
__host__ __device__ size_t size() const
{ return data.get<0>().size(); }
};
/*! \brief A 1d Odeint and Openfpm compatible structure.
*
* Use the method this.data.get<d>() to refer to property of all the particles in the dimension d.
*
* d starts with 0.
*
*/
struct state_type_1d_ofp_gpu{
state_type_1d_ofp_gpu(){
}
typedef size_t size_type;
typedef int is_state_vector;
aggregate<texp_v_gpu<double>> data;
size_t size() const
{ return data.get<0>().size(); }
void resize(size_t n)
{
data.get<0>().resize(n);
}
state_type_1d_ofp_ker toKernel() const
{
state_type_1d_ofp_ker s1_ker;
s1_ker.data.get<0>()=data.get<0>().getVector().toKernel();
return s1_ker;
}
};
#endif
namespace boost { namespace numeric { namespace odeint {
......@@ -184,7 +235,7 @@ struct state_type_ofpm_add_elements
template<typename state_type, typename ... list>
struct state_type_ofpm_add_elements<0,state_type,list ...>
{
typedef aggregate<list ...> type;
typedef aggregate<list ...> type;
};
template<int n_state, typename state_type>
......@@ -199,8 +250,8 @@ struct state_type_ofpm_impl
type_data data;
FD::gdb_ext_plus_g_info<state_type::dims> size() const
{
return data.template get<0>().size();
{
return data.template get<0>().size();
}
......@@ -212,19 +263,6 @@ struct state_type_ofpm_impl
namespace boost {
// template<typename state_type>
// struct range_size<const state_type_ofpm_impl<3,state_type>>
// {
// typedef FD::gdb_ext_plus_g_info<state_type::dims> type;
// };
// template<typename state_type>
// struct range_size<const state_type_ofpm_impl<1,state_type>>
// {
// typedef FD::gdb_ext_plus_g_info<state_type::dims> type;
// };
namespace numeric {
namespace odeint {
......@@ -235,7 +273,13 @@ namespace boost {
typedef boost::true_type type;
static const bool value = type::value;
};
#ifdef __NVCC__
template<>
struct is_resizeable<state_type_1d_ofp_gpu> {
typedef boost::true_type type;
static const bool value = type::value;
};
#endif
template<>
struct is_resizeable<state_type_2d_ofp> {
typedef boost::true_type type;
......
src/OdeIntegrators/tests/OdeIntegrator_grid_tests.cpp
View file @ 1371f3eb
......@@ -11,7 +11,6 @@
#include "config.h"
#include "Grid/grid_dist_id.hpp"
#include "OdeIntegrators/OdeIntegrators.hpp"
#include "OdeIntegrators/boost_vector_algebra_ofp.hpp"
#include "FiniteDifference/FD_op.hpp"
#include "util/util_debug.hpp"
#include "util/common.hpp"
......
src/OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
View file @ 1371f3eb
......@@ -15,8 +15,10 @@
#include "Vector/vector_dist_subset.hpp"
#include "Decomposition/Distribution/SpaceDistribution.hpp"
#include "OdeIntegrators/OdeIntegrators.hpp"
#ifdef HAVE_EIGEN
#include "DCPSE/DCPSE_op/DCPSE_op.hpp"
#include "OdeIntegrators/boost_vector_algebra_ofp.hpp"
#endif
#include "OdeIntegrators/vector_algebra_ofp.hpp"
typedef texp_v<double> state_type;
const double a = 2.8e-4;
......@@ -465,7 +467,6 @@ BOOST_AUTO_TEST_CASE(odeint_base_test3)
}
#ifdef HAVE_EIGEN
BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
size_t edgeSemiSize = 5;
const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
......@@ -567,5 +568,4 @@ BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
}
}
#endif
BOOST_AUTO_TEST_SUITE_END()
src/OdeIntegrators/tests/Odeintegrators_test_gpu.cu 0 → 100644
View file @ 1371f3eb
//
// Created by abhinav on 2/28/23.
//
#include "config.h"
#include <type_traits>
#include <cstring>
#include "util/common.hpp"
#define BOOST_TEST_DYN_LINK
#include "util/util_debug.hpp"
#include <boost/test/unit_test.hpp>
#include <iostream>
#include "Operators/Vector/vector_dist_operators.hpp"
#include "OdeIntegrators/OdeIntegrators.hpp"
//#include "DCPSE/DCPSE_op/DCPSE_op.hpp"
#ifdef __NVCC__
typedef state_type_1d_ofp_gpu state_type;
//const double a = 2.8e-4;
//const double b = 5e-3;
//const double tau = .1;
//const double k = .005;
void ExponentialGPU( const state_type &x , state_type &dxdt , const double t )
{
dxdt.data.get<0>() = x.data.get<0>();
//x.data.get<0>().getVector().deviceToHost<0>();
//dxdt.data.get<0>().getVector().deviceToHost<0>();
}
BOOST_AUTO_TEST_SUITE(odeInt_BASE_tests)
BOOST_AUTO_TEST_CASE(odeint_base_test_gpu)
{
size_t edgeSemiSize = 512;
const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
double spacing[2];
spacing[0] = 1.0 / (sz[0] - 1);
spacing[1] = 1.0 / (sz[1] - 1);
double rCut = 3.9 * spacing[0];
Ghost<2, double> ghost(rCut);
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist_gpu<2, double, aggregate<double, double,double>> Particles(0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
while (it.isNext())
{
Particles.add();
auto key = it.get();
mem_id k0 = key.get(0);
double xp0 = k0 * spacing[0];
Particles.getLastPos()[0] = xp0;
mem_id k1 = key.get(1);
double yp0 = k1 * spacing[1];
Particles.getLastPos()[1] = yp0;
Particles.getLastProp<0>() = xp0*yp0*exp(-5);
Particles.getLastProp<1>() = xp0*yp0*exp(5);
++it;
}
Particles.map();
Particles.ghost_get<0>();
Particles.hostToDeviceProp<0,1,2>();
auto Init = getV<0,comp_dev>(Particles);
auto Sol = getV<1,comp_dev>(Particles);
auto OdeSol = getV<2,comp_dev>(Particles);
state_type x0;
x0.data.get<0>()=Init;
x0.data.get<0>().getVector().deviceToHost<0>();
// The rhs of x' = f(x)
double t0=-5,tf=5;
const double dt=0.01;
//This doesnt work Why?
//size_t steps=boost::numeric::odeint::integrate(Exponential,x0,0.0,tf,dt);
timer tt;
tt.start();
size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type, double, state_type, double, boost::numeric::odeint::vector_space_algebra_ofp_gpu,boost::numeric::odeint::ofp_operations>(),ExponentialGPU,x0,t0,tf,dt);
tt.stop();
OdeSol=x0.data.get<0>();
Particles.deviceToHostProp<0,1,2>();
auto it2 = Particles.getDomainIterator();
double worst = 0.0;
while (it2.isNext()) {
auto p = it2.get();
if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
}
++it2;
}
std::cout<<"WCT:"<<tt.getwct()<<std::endl;
std::cout<<"CPU:"<<tt.getcputime()<<std::endl;
std::cout<<worst<<std::endl;
BOOST_REQUIRE(worst < 1e-6);
}
BOOST_AUTO_TEST_SUITE_END()
#endif
\ No newline at end of file
src/OdeIntegrators/boost_vector_algebra_ofp.hpp → src/OdeIntegrators/vector_algebra_ofp.hpp
View file @ 1371f3eb
......@@ -2,126 +2,40 @@
// Created by Abhinav Singh on 18.02.21.
//
#ifndef OPENFPM_PDATA_BOOST_VECTOR_ALGEBRA_OFP_HPP
#define OPENFPM_PDATA_BOOST_VECTOR_ALGEBRA_OFP_HPP
#ifndef OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
#define OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
namespace boost {
namespace numeric {
namespace odeint {
/*
* This class template has to be overload in order to call vector_space_algebra::norm_inf
*/
// template< class State, class Enabler = void > struct vector_space_norm_inf;
/*
* Example: instantiation for sole doubles and complex
*/
/* template<>
struct vector_space_norm_inf< double >
/* It copy one element of the chunk for each property
*
*/
template<typename vector_type,typename index_type,typename op_type>
struct for_each_prop1
{
typedef double result_type;
double operator()( double x ) const
{
using std::abs;
return abs(x);
}
};
template<>
struct vector_space_norm_inf< float >
{
typedef float result_type;
result_type operator()( float x ) const
vector_type &v;
index_type &p;
op_type &op;
/*! \brief constructor
*
*
* \param src source encapsulated object
* \param dst destination encapsulated object
*
*/
__device__ __host__ inline for_each_prop1(vector_type &v,index_type &p,op_type &op)
:v(v),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
__device__ __host__ inline void operator()(T& t) const
{
using std::abs;
return abs(x);
}
};
template< typename T >
struct vector_space_norm_inf< std::complex<T> >
{
typedef T result_type;
result_type operator()( std::complex<T> x ) const
{
using std::abs;
return abs( x );
}
};*/
template<typename S1,typename S2>
struct for_each_prop_resize{
S1 &v1;
S2 &v2;
/*! \brief constructor
*
*
* \param src source encapsulated object
* \param dst destination encapsulated object
*
*/
inline for_each_prop_resize(S1 &v1,S2 &v2)
:v1(v1),v2(v2)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
{
v1.data.template get<T::value>().getVector().resize(v2.data.template get<T::value>().getVector().size());
}
};
/* It copy one element of the chunk for each property
*
*/
template<typename vector_type,typename index_type,typename op_type>
struct for_each_prop1
{
vector_type &v;
index_type &p;
op_type &op;
/*! \brief constructor
*
*
* \param src source encapsulated object
* \param dst destination encapsulated object
*
*/
inline for_each_prop1(vector_type &v,index_type &p,op_type &op)
:v(v),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
{
op(v.data.template get<T::value>().getVector().template get<0>(p));
}
};
struct vector_space_algebra_ofp
{
template< class S1 , class Op >
static void for_each1( S1 &s1 , Op op )
{
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop1<S1,typename S1::index_type,Op> cp(s1,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
op(v.data.template get<T::value>().getVector().template get<0>(p));
}
}
};
template<typename S1,typename S2,typename index_type,typename op_type>
struct for_each_prop2
{
......@@ -137,36 +51,16 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop2(S1 &v1,S2 &v2,index_type &p,op_type &op)
:v1(v1),v2(v2),p(p),op(op)
__device__ __host__ inline for_each_prop2(S1 &v1,S2 &v2,index_type &p,op_type &op)
:v1(v1),v2(v2),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class Op >
static void for_each2( S1 &s1 , S2 &s2 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
//s1.data.template get<0>().getVector().resize(s2.data.template get<0>().getVector().size());
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop2<S1,S2,typename S1::index_type,Op> cp(s1,s2,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename index_type,typename op_type>
struct for_each_prop3
......@@ -184,41 +78,21 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop3(S1 &v1,S2 &v2,S3 &v3,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),p(p),op(op)
__device__ __host__ inline for_each_prop3(S1 &v1,S2 &v2,S3 &v3,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
//std::cout<<v1.data.template get<T::value>().getVector().size()<<":"<<v2.data.template get<T::value>().getVector().size()<<":"<<v3.data.template get<T::value>().getVector().size()<<std::endl;
//printf("v2:%f,v3:%f \n",v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p));
//printf("2\n");
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class Op >
static void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
{
//
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
//printf("v1:%f, v2:%f,v3:%f \n",v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p));
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop3<S1,S2,S3,typename S1::index_type,Op> cp(s1,s2,s3,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
};
template<typename S1,typename S2,typename S3,typename S4,typename index_type,typename op_type>
struct for_each_prop4
{
......@@ -237,36 +111,16 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop4(S1 &v1,S2 &v2,S3 &v3,S4 &v4,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),p(p),op(op)
__device__ __host__ inline for_each_prop4(S1 &v1,S2 &v2,S3 &v3,S4 &v4,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4 , class Op >
static void for_each4( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop4<S1,S2,S3,S4,typename S1::index_type,Op> cp(s1,s2,s3,s4,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename index_type,typename op_type>
struct for_each_prop5
{
......@@ -286,36 +140,17 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop5(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop5(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5 , class Op >
static void for_each5( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop5<S1,S2,S3,S4,S5,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename index_type,typename op_type>
struct for_each_prop6
{
......@@ -337,34 +172,16 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop6(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop6(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 , class Op >
static void for_each6( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop6<S1,S2,S3,S4,S5,S6,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename index_type,typename op_type>
......@@ -389,36 +206,17 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop7(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop7(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7, class Op >
static void for_each7( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop7<S1,S2,S3,S4,S5,S6,S7,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename index_type,typename op_type>
struct for_each_prop8
{
......@@ -442,36 +240,17 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop8(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop8(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class Op >
static void for_each8( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop8<S1,S2,S3,S4,S5,S6,S7,S8,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9,typename index_type,typename op_type>
struct for_each_prop9
{
......@@ -496,38 +275,19 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop9(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop9(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class Op >
static void for_each9( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop9<S1,S2,S3,S4,S5,S6,S7,S8,S9,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10,typename index_type,typename op_type>
struct for_each_prop10
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10,typename index_type,typename op_type>
struct for_each_prop10
{
S1 &v1;
......@@ -551,36 +311,17 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop10(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop10(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class Op >
static void for_each10( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop10<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11,typename index_type,typename op_type>
struct for_each_prop11
{
......@@ -607,36 +348,18 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop11(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop11(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class Op >
static void for_each11( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop11<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12,typename index_type,typename op_type>
struct for_each_prop12
{
......@@ -664,36 +387,18 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop12(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop12(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class Op >
static void for_each12( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop12<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13,typename index_type,typename op_type>
struct for_each_prop13
{
......@@ -722,36 +427,18 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop13(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop13(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class Op >
static void for_each13( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12,S13 &s13, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop13<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13, typename S14,typename index_type,typename op_type>
struct for_each_prop14
{
......@@ -781,36 +468,18 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop14(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop14(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),v14(v14),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p),v14.data.template get<T::value>().getVector().template get<0>(p));
}
};
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class Op >
static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop14<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename S9, typename S10, typename S11, typename S12, typename S13, typename S14, typename S15,typename index_type,typename op_type>
struct for_each_prop15
{
......@@ -841,17 +510,349 @@ namespace boost {
* \param dst destination encapsulated object
*
*/
inline for_each_prop15(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,S15 &v15,index_type &p,op_type &op)
__device__ __host__ inline for_each_prop15(S1 &v1,S2 &v2,S3 &v3,S4 &v4,S5 &v5,S6 &v6,S7 &v7,S8 &v8,S9 &v9,S10 &v10, S11 &v11,S12 &v12,S13 &v13,S14 &v14,S15 &v15,index_type &p,op_type &op)
:v1(v1),v2(v2),v3(v3),v4(v4),v5(v5),v6(v6),v7(v7),v8(v8),v9(v9),v10(v10),v11(v11), v12(v12),v13(v13),v14(v14),v15(v15),p(p),op(op)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
__device__ __host__ inline void operator()(T& t) const
{
op(v1.data.template get<T::value>().getVector().template get<0>(p),v2.data.template get<T::value>().getVector().template get<0>(p),v3.data.template get<T::value>().getVector().template get<0>(p),v4.data.template get<T::value>().getVector().template get<0>(p),v5.data.template get<T::value>().getVector().template get<0>(p),v6.data.template get<T::value>().getVector().template get<0>(p),v7.data.template get<T::value>().getVector().template get<0>(p),v8.data.template get<T::value>().getVector().template get<0>(p),v9.data.template get<T::value>().getVector().template get<0>(p),v10.data.template get<T::value>().getVector().template get<0>(p),v11.data.template get<T::value>().getVector().template get<0>(p),v12.data.template get<T::value>().getVector().template get<0>(p),v13.data.template get<T::value>().getVector().template get<0>(p),v14.data.template get<T::value>().getVector().template get<0>(p),v15.data.template get<T::value>().getVector().template get<0>(p));
}
};
/*
* This class template has to be overload in order to call vector_space_algebra::norm_inf
*/
// template< class State, class Enabler = void > struct vector_space_norm_inf;
/*
* Example: instantiation for sole doubles and complex
*/
/* template<>
struct vector_space_norm_inf< double >
{
typedef double result_type;
double operator()( double x ) const
{
using std::abs;
return abs(x);
}
};
template<>
struct vector_space_norm_inf< float >
{
typedef float result_type;
result_type operator()( float x ) const
{
using std::abs;
return abs(x);
}
};
template< typename T >
struct vector_space_norm_inf< std::complex<T> >
{
typedef T result_type;
result_type operator()( std::complex<T> x ) const
{
using std::abs;
return abs( x );
}
};*/
template<typename S1,typename S2>
struct for_each_prop_resize{
S1 &v1;
S2 &v2;
/*! \brief constructor
*
*
* \param src source encapsulated object
* \param dst destination encapsulated object
*
*/
inline for_each_prop_resize(S1 &v1,S2 &v2)
:v1(v1),v2(v2)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
{
v1.data.template get<T::value>().getVector().resize(v2.data.template get<T::value>().getVector().size());
}
};
struct vector_space_algebra_ofp
{
template< class S1 , class Op >
static void for_each1( S1 &s1 , Op op )
{
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop1<S1,typename S1::index_type,Op> cp(s1,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class Op >
static void for_each2( S1 &s1 , S2 &s2 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
//s1.data.template get<0>().getVector().resize(s2.data.template get<0>().getVector().size());
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop2<S1,S2,typename S1::index_type,Op> cp(s1,s2,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class Op >
static void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
{
//
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop3<S1,S2,S3,typename S1::index_type,Op> cp(s1,s2,s3,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4 , class Op >
static void for_each4( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop4<S1,S2,S3,S4,typename S1::index_type,Op> cp(s1,s2,s3,s4,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5 , class Op >
static void for_each5( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop5<S1,S2,S3,S4,S5,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 , class Op >
static void for_each6( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop6<S1,S2,S3,S4,S5,S6,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7, class Op >
static void for_each7( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop7<S1,S2,S3,S4,S5,S6,S7,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class Op >
static void for_each8( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop8<S1,S2,S3,S4,S5,S6,S7,S8,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class Op >
static void for_each9( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop9<S1,S2,S3,S4,S5,S6,S7,S8,S9,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class Op >
static void for_each10( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop10<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class Op >
static void for_each11( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop11<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class Op >
static void for_each12( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop12<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class Op >
static void for_each13( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12,S13 &s13, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop13<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class Op >
static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_prop14<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,typename S1::index_type,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
++it;
}
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class S15, class Op >
static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14,S15 &s15, Op op )
......@@ -940,4 +941,4 @@ namespace boost {
#endif //OPENFPM_PDATA_BOOST_VECTOR_ALGEBRA_OFP_HPP
#endif //OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
src/OdeIntegrators/vector_algebra_ofp_gpu.hpp 0 → 100644
View file @ 1371f3eb
//
// Created by Abhinav Singh on 1.03.23.
//
#ifndef OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_GPU_HPP
#define OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_GPU_HPP
namespace boost {
namespace numeric {
namespace odeint {
template<typename S1, typename Op>
__global__ void for_each1_ker(S1 s1, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.size()) {return;}
for_each_prop1<S1,size_t,Op> cp(s1,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2, typename Op>
__global__ void for_each2_ker(S1 s1,S2 s2, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//printf("%f \n",s2.data.template get<0>().getVector().template get<0>(p));
//converting to boost vector ids.
for_each_prop2<S1,S2,unsigned int,Op> cp(s1,s2,p,op);
//s1.data.template get<0>().getVector().template get<0>(p)=1.0*s2.data.template get<0>().getVector().template get<0>(p)+0.05*s3.data.template get<0>().getVector().template get<0>(p);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3, typename Op>
__global__ void for_each3_ker(S1 s1,S2 s2,S3 s3, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//printf("%f \n",s2.data.template get<0>().getVector().template get<0>(p));
//converting to boost vector ids.
for_each_prop3<S1,S2,S3,unsigned int,Op> cp(s1,s2,s3,p,op);
//s1.data.template get<0>().getVector().template get<0>(p)=1.0*s2.data.template get<0>().getVector().template get<0>(p)+0.05*s3.data.template get<0>().getVector().template get<0>(p);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4, typename Op>
__global__ void for_each4_ker(S1 s1,S2 s2,S3 s3,S4 s4, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//printf("%f \n",s2.data.template get<0>().getVector().template get<0>(p));
//converting to boost vector ids.
for_each_prop4<S1,S2,S3,S4,unsigned int,Op> cp(s1,s2,s3,s4,p,op);
//s1.data.template get<0>().getVector().template get<0>(p)=1.0*s2.data.template get<0>().getVector().template get<0>(p)+0.05*s3.data.template get<0>().getVector().template get<0>(p);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5, typename Op>
__global__ void for_each5_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//printf("%f \n",s2.data.template get<0>().getVector().template get<0>(p));
//converting to boost vector ids.
for_each_prop5<S1,S2,S3,S4,S5,unsigned int,Op> cp(s1,s2,s3,s4,s5,p,op);
//s1.data.template get<0>().getVector().template get<0>(p)=1.0*s2.data.template get<0>().getVector().template get<0>(p)+0.05*s3.data.template get<0>().getVector().template get<0>(p);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6, typename Op>
__global__ void for_each6_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop6<S1,S2,S3,S4,S5,S6,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7, typename Op>
__global__ void for_each7_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop7<S1,S2,S3,S4,S5,S6,S7,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8, typename Op>
__global__ void for_each8_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop8<S1,S2,S3,S4,S5,S6,S7,S8,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9, typename Op>
__global__ void for_each9_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop9<S1,S2,S3,S4,S5,S6,S7,S8,S9,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10, typename Op>
__global__ void for_each10_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop10<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10,typename S11, typename Op>
__global__ void for_each11_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10,S11 s11, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop11<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10,typename S11,typename S12, typename Op>
__global__ void for_each12_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10,S11 s11,S12 s12, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop12<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10,typename S11,typename S12,typename S13, typename Op>
__global__ void for_each13_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10,S11 s11,S12 s12,S13 s13, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop13<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10,typename S11,typename S12,typename S13,typename S14, typename Op>
__global__ void for_each14_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10,S11 s11,S12 s12,S13 s13,S14 s14, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop14<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
template<typename S1,typename S2,typename S3,typename S4,typename S5,typename S6,typename S7,typename S8,typename S9,typename S10,typename S11,typename S12,typename S13,typename S14,typename S15, typename Op>
__global__ void for_each15_ker(S1 s1,S2 s2,S3 s3,S4 s4,S5 s5,S6 s6,S7 s7,S8 s8,S9 s9,S10 s10,S11 s11,S12 s12,S13 s13,S14 s14,S15 s15, Op op)
{
unsigned int p = threadIdx.x + blockIdx.x * blockDim.x;
if (p >= s1.data.template get<0>().size()) {return;}
//converting to boost vector ids.
for_each_prop15<S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15,unsigned int,Op> cp(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,p,op);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(cp);
}
struct vector_space_algebra_ofp_gpu
{
template< class S1 , class Op >
static void for_each1( S1 &s1 , Op op )
{
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getDomainIteratorGPU();
CUDA_LAUNCH((for_each1_ker),it,s1,op);
}
template< class S1 , class S2 , class Op >
static void for_each2( S1 &s1 , S2 &s2 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
//s1.data.template get<0>().getVector().resize(s2.data.template get<0>().getVector().size());
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each2_ker),it,s1.toKernel(),s2.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class Op >
static void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype( s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each3_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4 , class Op >
static void for_each4( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each4_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5 , class Op >
static void for_each5( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each5_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 , class Op >
static void for_each6( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each6_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7, class Op >
static void for_each7( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each7_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class Op >
static void for_each8( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each8_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class Op >
static void for_each9( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each9_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class Op >
static void for_each10( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each10_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),s10.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class Op >
static void for_each11( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each11_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class Op >
static void for_each12( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each12_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),s10.toKernel(),s11.toKernel(),s12.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class Op >
static void for_each13( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9 , S10 &s10,S11 &s11,S12 &s12,S13 &s13, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each13_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),s10.toKernel(),s11.toKernel(),s12.toKernel(),s13.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class Op >
static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
// ToDo : build checks, that the +-*/ operators are well defined
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each14_ker),it,it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),s10.toKernel(),s11.toKernel(),s12.toKernel(),s13.toKernel(),s14.toKernel(),op);
}
template< class S1 , class S2 , class S3 , class S4,class S5,class S6 ,class S7,class S8, class S9, class S10, class S11, class S12, class S13, class S14, class S15, class Op >
static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4,S5 &s5,S6 &s6,S7 &s7,S8 &s8, S9 &s9, S10 &s10,S11 &s11,S12 &s12,S13 &s13,S14 &s14,S15 &s15, Op op )
{
for_each_prop_resize<S1,S2> the_resize(s1,s2);
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s1.data)::max_prop>>(the_resize);
auto it=s1.data.template get<0>().getVector().getGPUIterator();
CUDA_LAUNCH((for_each15_ker),it,s1.toKernel(),s2.toKernel(),s3.toKernel(),s4.toKernel(),s5.toKernel(),s6.toKernel(),s7.toKernel(),s8.toKernel(),s9.toKernel(),s10.toKernel(),s11.toKernel(),s12.toKernel(),s13.toKernel(),s14.toKernel(),s15.toKernel(),op);
}
template<typename vector_type,typename index_type,typename norm_result_type>
struct for_each_norm
{
const vector_type &v;
index_type &p;
norm_result_type &n;
/*! \brief constructor
*
*
* \param src source encapsulated object
* \param dst destination encapsulated object
*
*/
inline for_each_norm(const vector_type &v,index_type &p,norm_result_type &n)
:v(v),p(p),n(n)
{};
//! It call the copy function for each property
template<typename T>
inline void operator()(T& t) const
{
if(fabs(v.data.template get<T::value>().getVector().template get<0>(p)) > n)
{
n=fabs(v.data.template get<T::value>().getVector().template get<0>(p));
}
}
};
template< class S >
static typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type norm_inf( const S &s )
{
typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type n=0;
auto it=s.data.template get<0>().getVector().getIterator();
while(it.isNext()){
auto p=it.get();
//converting to boost vector ids.
for_each_norm<S,size_t,typename boost::numeric::odeint::vector_space_norm_inf< S >::result_type> cp(s,p,n);
//creating an iterator on v_ids[0] [1] [2]
boost::mpl::for_each_ref<boost::mpl::range_c<int,0,decltype(s.data)::max_prop>>(cp);
++it;
}
auto &v_cl = create_vcluster();
v_cl.max(n);
v_cl.execute();
//std::max();
//std::cout<<n<<std::endl;
return n;
}
};
#include <algorithm>
#include <boost/config.hpp>
#include <boost/array.hpp>
#include <boost/numeric/odeint/util/unit_helper.hpp>
/*
* Notes:
*
* * the results structs are needed in order to work with fusion_algebra
*/
struct ofp_operations
{
template< class Fac1 = double >
struct scale
{
const Fac1 m_alpha1;
scale( Fac1 alpha1 ) : m_alpha1( alpha1 ) { }
template< class T1 >
__device__ __host__ void operator()( T1 &t1 ) const
{
t1 *= m_alpha1;
}
typedef void result_type;
};
template< class Fac1 = double >
struct scale_sum1
{
const Fac1 m_alpha1;
scale_sum1( Fac1 alpha1 ) : m_alpha1( alpha1 ) { }
template< class T1 , class T2 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 ) const
{
t1 = m_alpha1 * t2;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 >
struct scale_sum2
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
scale_sum2( Fac1 alpha1 , Fac2 alpha2 ) : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
template< class T1 , class T2 , class T3 >
BOOST_FUSION_GPU_ENABLED
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 >
struct scale_sum3
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
scale_sum3( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { }
template< class T1 , class T2 , class T3 , class T4 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 >
struct scale_sum4
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
scale_sum4( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 >
struct scale_sum5
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
scale_sum5( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 , Fac5 alpha5 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 >
struct scale_sum6
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
scale_sum6( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 , Fac5 alpha5 , Fac6 alpha6 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ){ }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 ,const T7 &t7) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 >
struct scale_sum7
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
scale_sum7( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 >
struct scale_sum8
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
scale_sum8( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 >
struct scale_sum9
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
scale_sum9( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 , class Fac10 = Fac9 >
struct scale_sum10
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
const Fac10 m_alpha10;
scale_sum10( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 , Fac10 alpha10 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) , m_alpha10( alpha10 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 , class T11 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 , const T11 &t11 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 , class Fac10 = Fac9 , class Fac11 = Fac10 >
struct scale_sum11
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
const Fac10 m_alpha10;
const Fac11 m_alpha11;
scale_sum11( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 ,
Fac10 alpha10 , Fac11 alpha11 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) , m_alpha10( alpha10 ) , m_alpha11( alpha11 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 , class T11 , class T12 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 , const T11 &t11 , const T12 &t12 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 + m_alpha11 * t12;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 , class Fac10 = Fac9 , class Fac11 = Fac10 , class Fac12 = Fac11 >
struct scale_sum12
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
const Fac10 m_alpha10;
const Fac11 m_alpha11;
const Fac12 m_alpha12;
scale_sum12( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 ,
Fac10 alpha10 , Fac11 alpha11 , Fac12 alpha12 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) , m_alpha10( alpha10 ) , m_alpha11( alpha11 ) , m_alpha12( alpha12 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 , class T11 , class T12 , class T13 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 , const T11 &t11 , const T12 &t12 , const T13 &t13 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 + m_alpha11 * t12 + m_alpha12 * t13;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 , class Fac10 = Fac9 , class Fac11 = Fac10 , class Fac12 = Fac11 , class Fac13 = Fac12 >
struct scale_sum13
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
const Fac10 m_alpha10;
const Fac11 m_alpha11;
const Fac12 m_alpha12;
const Fac13 m_alpha13;
scale_sum13( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 ,
Fac10 alpha10 , Fac11 alpha11 , Fac12 alpha12 , Fac13 alpha13 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) , m_alpha10( alpha10 ) , m_alpha11( alpha11 ) , m_alpha12( alpha12 ) , m_alpha13( alpha13 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 , class T11 , class T12 , class T13 , class T14 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 , const T11 &t11 , const T12 &t12 , const T13 &t13 , const T14 &t14 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 + m_alpha11 * t12 + m_alpha12 * t13 + m_alpha13 * t14;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 , class Fac3 = Fac2 , class Fac4 = Fac3 , class Fac5 = Fac4 , class Fac6 = Fac5 , class Fac7 = Fac6 , class Fac8 = Fac7 , class Fac9 = Fac8 , class Fac10 = Fac9 , class Fac11 = Fac10 , class Fac12 = Fac11 , class Fac13 = Fac12 , class Fac14 = Fac13 >
struct scale_sum14
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
const Fac3 m_alpha3;
const Fac4 m_alpha4;
const Fac5 m_alpha5;
const Fac6 m_alpha6;
const Fac7 m_alpha7;
const Fac8 m_alpha8;
const Fac9 m_alpha9;
const Fac10 m_alpha10;
const Fac11 m_alpha11;
const Fac12 m_alpha12;
const Fac13 m_alpha13;
const Fac14 m_alpha14;
scale_sum14( Fac1 alpha1 , Fac2 alpha2 , Fac3 alpha3 , Fac4 alpha4 ,
Fac5 alpha5 , Fac6 alpha6 , Fac7 alpha7 , Fac8 alpha8 , Fac9 alpha9 ,
Fac10 alpha10 , Fac11 alpha11 , Fac12 alpha12 , Fac13 alpha13 , Fac14 alpha14 )
: m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) , m_alpha6( alpha6 ) , m_alpha7( alpha7 ) , m_alpha8( alpha8 ) , m_alpha9( alpha9 ) , m_alpha10( alpha10 ) , m_alpha11( alpha11 ) , m_alpha12( alpha12 ) , m_alpha13( alpha13 ) , m_alpha14( alpha14 ) { }
template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 , class T7 , class T8 , class T9 , class T10 , class T11 , class T12 , class T13 , class T14 , class T15 >
__device__ __host__ void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 , const T7 &t7 , const T8 &t8 , const T9 &t9 , const T10 &t10 , const T11 &t11 , const T12 &t12 , const T13 &t13 , const T14 &t14 , const T15 &t15 ) const
{
t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6 + m_alpha6 * t7 + m_alpha7 * t8 + m_alpha8 * t9 + m_alpha9 * t10 + m_alpha10 * t11 + m_alpha11 * t12 + m_alpha12 * t13 + m_alpha13 * t14 + m_alpha14 * t15;
}
typedef void result_type;
};
template< class Fac1 = double , class Fac2 = Fac1 >
struct scale_sum_swap2
{
const Fac1 m_alpha1;
const Fac2 m_alpha2;
scale_sum_swap2( Fac1 alpha1 , Fac2 alpha2 ) : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
template< class T1 , class T2 , class T3 >
__device__ __host__ void operator()( T1 &t1 , T2 &t2 , const T3 &t3) const
{
const T1 tmp( t1 );
t1 = m_alpha1 * t2 + m_alpha2 * t3;
t2 = tmp;
}
typedef void result_type;
};
/*
* for usage in for_each2
*
* Works with boost::units by eliminating the unit
*/
template< class Fac1 = double >
struct rel_error
{
const Fac1 m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt;
rel_error( Fac1 eps_abs , Fac1 eps_rel , Fac1 a_x , Fac1 a_dxdt )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt ) { }
template< class T1 , class T2 , class T3 >
__device__ __host__ void operator()( T3 &t3 , const T1 &t1 , const T2 &t2 ) const
{
using std::abs;
set_unit_value( t3 , abs( get_unit_value( t3 ) ) / ( m_eps_abs + m_eps_rel * ( m_a_x * abs( get_unit_value( t1 ) ) + m_a_dxdt * abs( get_unit_value( t2 ) ) ) ) );
}
typedef void result_type;
};
/*
* for usage in for_each3
*
* used in the controller for the rosenbrock4 method
*
* Works with boost::units by eliminating the unit
*/
template< class Fac1 = double >
struct default_rel_error
{
const Fac1 m_eps_abs , m_eps_rel ;
default_rel_error( Fac1 eps_abs , Fac1 eps_rel )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) { }
/*
* xerr = xerr / ( eps_abs + eps_rel * max( x , x_old ) )
*/
template< class T1 , class T2 , class T3 >
__device__ __host__ void operator()( T3 &t3 , const T1 &t1 , const T2 &t2 ) const
{
BOOST_USING_STD_MAX();
using std::abs;
Fac1 x1 = abs( get_unit_value( t1 ) ) , x2 = abs( get_unit_value( t2 ) );
set_unit_value( t3 , abs( get_unit_value( t3 ) ) / ( m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION ( x1 , x2 ) ) );
}
typedef void result_type;
};
/*
* for usage in reduce
*/
template< class Value >
struct maximum
{
template< class Fac1 , class Fac2 >
__device__ __host__ Value operator()( Fac1 t1 , const Fac2 t2 ) const
{
using std::abs;
Value a1 = abs( get_unit_value( t1 ) ) , a2 = abs( get_unit_value( t2 ) );
return ( a1 < a2 ) ? a2 : a1 ;
}
typedef Value result_type;
};
template< class Fac1 = double >
struct rel_error_max
{
const Fac1 m_eps_abs , m_eps_rel;
rel_error_max( Fac1 eps_abs , Fac1 eps_rel )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel )
{ }
template< class Res , class T1 , class T2 , class T3 >
__device__ __host__ Res operator()( Res r , const T1 &x_old , const T2 &x , const T3 &x_err )
{
BOOST_USING_STD_MAX();
using std::abs;
Res tmp = abs( get_unit_value( x_err ) ) / ( m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( x_old ) , abs( x ) ) );
return max BOOST_PREVENT_MACRO_SUBSTITUTION ( r , tmp );
}
};
template< class Fac1 = double >
struct rel_error_max2
{
const Fac1 m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt;
rel_error_max2( Fac1 eps_abs , Fac1 eps_rel , Fac1 a_x , Fac1 a_dxdt )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
template< class Res , class T1 , class T2 , class T3 , class T4 >
__device__ __host__ Res operator()( Res r , const T1 &x_old , const T2 &/*x*/ , const T3 &dxdt_old , const T4 &x_err )
{
BOOST_USING_STD_MAX();
using std::abs;
Res tmp = abs( get_unit_value( x_err ) ) /
( m_eps_abs + m_eps_rel * ( m_a_x * abs( get_unit_value( x_old ) ) + m_a_dxdt * abs( get_unit_value( dxdt_old ) ) ) );
return max BOOST_PREVENT_MACRO_SUBSTITUTION ( r , tmp );
}
};
template< class Fac1 = double >
struct rel_error_l2
{
const Fac1 m_eps_abs , m_eps_rel;
rel_error_l2( Fac1 eps_abs , Fac1 eps_rel )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel )
{ }
template< class Res , class T1 , class T2 , class T3 >
__device__ __host__ Res operator()( Res r , const T1 &x_old , const T2 &x , const T3 &x_err )
{
BOOST_USING_STD_MAX();
using std::abs;
Res tmp = abs( get_unit_value( x_err ) ) / ( m_eps_abs + m_eps_rel * max BOOST_PREVENT_MACRO_SUBSTITUTION ( abs( x_old ) , abs( x ) ) );
return r + tmp * tmp;
}
};
template< class Fac1 = double >
struct rel_error_l2_2
{
const Fac1 m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt;
rel_error_l2_2( Fac1 eps_abs , Fac1 eps_rel , Fac1 a_x , Fac1 a_dxdt )
: m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
{ }
template< class Res , class T1 , class T2 , class T3 , class T4 >
__device__ __host__ Res operator()( Res r , const T1 &x_old , const T2 &/*x*/ , const T3 &dxdt_old , const T4 &x_err )
{
using std::abs;
Res tmp = abs( get_unit_value( x_err ) ) /
( m_eps_abs + m_eps_rel * ( m_a_x * abs( get_unit_value( x_old ) ) + m_a_dxdt * abs( get_unit_value( dxdt_old ) ) ) );
return r + tmp * tmp;
}
};
};
} // odeint
} // numeric
} // boost
#endif //OPENFPM_PDATA_VECTOR_ALGEBRA_OFP_HPP
src/Operators/Vector/cuda/vector_dist_operators_cuda.cuh
View file @ 1371f3eb
......@@ -80,7 +80,6 @@ struct pos_or_propL_ker
};
/*! \brief selector for position or properties left side
*
* \tparam vector type of the original vector
......@@ -428,7 +427,7 @@ struct vector_dist_op_compute_op<prp,false,comp_host>
}
};
#define NVCC
#ifdef __NVCC__
template<unsigned int prp, unsigned int dim ,typename vector, typename expr>
......