Skip to content
Snippets Groups Projects
Commit e4f32ebd authored by yaskovet's avatar yaskovet
Browse files

DCPSE on GPU: add additional tests

parent 1371f3eb
No related branches found
No related tags found
No related merge requests found
......@@ -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
DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests.cu
DCPSE/DCPSE_op/tests/DCPSE_op_test_temporal.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)
......@@ -39,7 +41,7 @@ if ( CUDA_ON_BACKEND STREQUAL "HIP" AND HIP_FOUND )
OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
OdeIntegrators/tests/OdeIntegrator_grid_tests.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests.cpp
FiniteDifference/FD_Solver_test.cpp
FiniteDifference/FD_op_Tests.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
......@@ -80,7 +82,7 @@ else()
OdeIntegrators/tests/OdeIntegratores_base_tests.cpp
OdeIntegrators/tests/OdeIntegrator_grid_tests.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_subset_test.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests
DCPSE/DCPSE_op/tests/DCPSE_op_test_base_tests.cpp
FiniteDifference/FD_Solver_test.cpp
FiniteDifference/FD_op_Tests.cpp
DCPSE/DCPSE_op/tests/DCPSE_op_test3d.cpp
......
......@@ -24,127 +24,6 @@
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests_cu)
BOOST_AUTO_TEST_CASE(dcpse_op_vec3d_gpu) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
size_t edgeSemiSize = 257;
const size_t sz[3] = {edgeSemiSize, edgeSemiSize,edgeSemiSize};
Box<3, double> box({0, 0,0}, {1,1,1});
size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
double rCut = 3.1 * spacing;
Ghost<3, double> ghost(rCut);
BOOST_TEST_MESSAGE("Init vector_dist...");
double sigma2 = spacing * spacing/ (2 * 4);
vector_dist_gpu<3, double, aggregate<double, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>,double,double>> domain(
0, box, bc, ghost);
//Init_DCPSE(domain)
BOOST_TEST_MESSAGE("Init domain...");
auto it = domain.getGridIterator(sz);
size_t pointId = 0;
size_t counter = 0;
double minNormOne = 999;
while (it.isNext()) {
domain.add();
auto key = it.get();
mem_id k0 = key.get(0);
double x = k0 * spacing;
domain.getLastPos()[0] = x;//+ gaussian(rng);
mem_id k1 = key.get(1);
double y = k1 * spacing;
domain.getLastPos()[1] = y;//+gaussian(rng);
mem_id k2 = key.get(2);
double z = k2 * spacing;
domain.getLastPos()[2] = z;//+gaussian(rng);
// Here fill the function value
domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]) + sin(domain.getLastPos()[2]) ;
domain.template getLastProp<1>()[0] = cos(domain.getLastPos()[0]);
domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[1]) ;
domain.template getLastProp<1>()[2] = cos(domain.getLastPos()[2]);
// Here fill the validation value for Df/Dx
domain.template getLastProp<2>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
domain.template getLastProp<2>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
domain.template getLastProp<3>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
domain.template getLastProp<3>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
domain.template getLastProp<3>()[2] = 0;
domain.template getLastProp<4>()[0] = -cos(domain.getLastPos()[0]) * sin(domain.getLastPos()[0]);
domain.template getLastProp<4>()[1] = -cos(domain.getLastPos()[1]) * sin(domain.getLastPos()[1]);
domain.template getLastProp<4>()[2] = -cos(domain.getLastPos()[2]) * sin(domain.getLastPos()[2]);
/* domain.template getLastProp<4>()[0] = cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
domain.template getLastProp<4>()[2] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));*/
domain.template getLastProp<5>() = cos(domain.getLastPos()[0]) * cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]) * cos(domain.getLastPos()[1])+cos(domain.getLastPos()[2]) * cos(domain.getLastPos()[2]) ;
++counter;
++it;
}
BOOST_TEST_MESSAGE("Sync domain across processors...");
domain.map();
domain.ghost_get<0>();
Advection_gpu Adv(domain, 2, rCut, 1.9,support_options::RADIUS);
auto v = getV<1>(domain);
auto P = getV<0>(domain);
auto dv = getV<3>(domain);
auto dP = getV<6>(domain);
// typedef boost::mpl::int_<std::is_fundamental<point_expression_op<Point<2U, double>, point_expression<double>, Point<2U, double>, 3>>::value>::blabla blabla;
// std::is_fundamental<decltype(o1.value(key))>
domain.ghost_get<1>();
dv = Adv(v, v);
auto it2 = domain.getDomainIterator();
double worst1 = 0.0;
while (it2.isNext()) {
auto p = it2.get();
if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
}
++it2;
}
//std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
BOOST_REQUIRE(worst1 < 0.03);
//Adv.checkMomenta(domain);
//Adv.DrawKernel<2>(domain,0);
//domain.deleteGhost();
dP = Adv(v, P);//+Dy(P);
auto it3 = domain.getDomainIterator();
double worst2 = 0.0;
while (it3.isNext()) {
auto p = it3.get();
if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
}
++it3;
}
domain.deleteGhost();
BOOST_REQUIRE(worst2 < 0.03);
}
BOOST_AUTO_TEST_CASE(dcpse_op_solver) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
......
This diff is collapsed.
/*
* DCPSE_op_test_temporal.cpp
*
* Created on: Sep 15, 2020
* Author: i-bird
*/
#include "config.h"
#ifdef HAVE_EIGEN
#ifdef HAVE_PETSC
#ifndef DCPSE_OP_TEST_TEMPORAL_CPP_
#define DCPSE_OP_TEST_TEMPORAL_CPP_
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "Operators/Vector/vector_dist_operators.hpp"
#include "../DCPSE_op.hpp"
BOOST_AUTO_TEST_SUITE(temporal_test_suite_cu)
BOOST_AUTO_TEST_CASE(temporal_test)
{
size_t grd_sz = 12;
double boxsize = 10;
const size_t sz[3] = {grd_sz, grd_sz, grd_sz};
Box<3, double> box({0, 0, 0}, {boxsize, boxsize, boxsize});
size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
double Lx = box.getHigh(0);
double Ly = box.getHigh(1);
double Lz = box.getHigh(2);
double spacing = box.getHigh(0) / (sz[0] - 1);
double rCut = 3.9 * spacing;
double rCut2 = 3.9 * spacing;
int ord = 2;
int ord2 = 2;
double sampling_factor = 4.0;
double sampling_factor2 = 2.4;
Ghost<3, double> ghost(rCut);
auto &v_cl = create_vcluster();
vector_dist_gpu<3, double, aggregate<double,VectorS<3, double>,double[3][3],double,VectorS<3,double>,double[3][3]> > Particles(0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
while (it.isNext())
{
Particles.add();
auto key = it.get();
double x = key.get(0) * it.getSpacing(0);
Particles.getLastPos()[0] = x;
double y = key.get(1) * it.getSpacing(1);
Particles.getLastPos()[1] = y;
double z = key.get(2) * it.getSpacing(1);
Particles.getLastPos()[2] = z;
++it;
}
Particles.map();
Particles.ghost_get<0>();
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
constexpr int sScalar = 0;
constexpr int sVector = 1;
constexpr int sTensor = 2;
constexpr int dScalar = 3;
constexpr int dVector = 4;
constexpr int dTensor = 5;
auto Pos = getV<PROP_POS>(Particles);
auto sS = getV<sScalar>(Particles);
auto sV = getV<sVector>(Particles);
auto sT = getV<sTensor>(Particles);
auto dS = getV<dScalar>(Particles);
auto dV = getV<dVector>(Particles);
auto dT = getV<dTensor>(Particles);
//Particles_subset.write("Pars");
Derivative_x_gpu Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Bulk_Dx(Particles, ord,
rCut, sampling_factor,
support_options::RADIUS);
texp_v<double> TVx,TdxVx;
texp_v<VectorS<3, double>> TV;
texp_v<double[3][3]> TT;
auto it3 = Particles.getDomainIterator();
sS = 5;
sV[0] = 1;
sV[1] = 2;
sV[2] = 3;
sT[0][0] = 1;
sT[0][1] = 2;
sT[0][2] = 3;
sT[1][0] = 4;
sT[1][1] = 5;
sT[1][2] = 6;
sT[2][0] = 7;
sT[2][1] = 8;
sT[2][2] = 9;
TVx=sS;
dS = TVx;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
match &= Particles.template getProp<sScalar>(key) == Particles.template getProp<dScalar>(key);
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TVx=sS*sS;
dS = TVx;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
match &= Particles.template getProp<sScalar>(key)*Particles.template getProp<sScalar>(key) == Particles.template getProp<dScalar>(key);
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TVx=sV[0];
dS = TVx;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
match &= Particles.template getProp<sVector>(key)[0] == Particles.template getProp<dScalar>(key);
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TVx=sT[0][0];
dS = TVx;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
match &= Particles.template getProp<sVector>(key)[0] == Particles.template getProp<dScalar>(key);
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TV=sV;
dV=TV;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
match &= Particles.template getProp<sVector>(key)[0] == Particles.template getProp<dVector>(key)[0];
match &= Particles.template getProp<sVector>(key)[1] == Particles.template getProp<dVector>(key)[1];
match &= Particles.template getProp<sVector>(key)[2] == Particles.template getProp<dVector>(key)[2];
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TV=0.5*sV+sV;
dV=TV;
//Pol_bulk = dPol + (0.5 * dt) * k1;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
double x1=Particles.template getProp<sVector>(key)[0];
double y1=Particles.template getProp<dVector>(key)[0];
double x2=Particles.template getProp<sVector>(key)[1];
double y2=Particles.template getProp<dVector>(key)[1];
double x3=Particles.template getProp<sVector>(key)[2];
double y3=Particles.template getProp<dVector>(key)[2];
match &= 1.5*Particles.template getProp<sVector>(key)[0] == Particles.template getProp<dVector>(key)[0];
match &= 1.5*Particles.template getProp<sVector>(key)[1] == Particles.template getProp<dVector>(key)[1];
match &= 1.5*Particles.template getProp<sVector>(key)[2] == Particles.template getProp<dVector>(key)[2];
++it3;
}
BOOST_REQUIRE_EQUAL(match,true);
}
TV=sV;
dV=pmul(TV,TV);
//Pol_bulk = dPol + (0.5 * dt) * k1;
{
auto it3 = Particles.getDomainIterator();
bool match = true;
while (it3.isNext())
{
auto key = it3.get();
double x1=Particles.template getProp<sVector>(key)[0];
double y1=Particles.template getProp<dVector>(key)[0];
double x2=Particles.template getProp<sVector>(key)[1];
double y2=Particles.template getProp<dVector>(key)[1];
double x3=Particles.template getProp<sVector>(key)[2];
double y3=Particles.template getProp<dVector>(key)[2];
match &= Particles.template getProp<sVector>(key)[0]*Particles.template getProp<sVector>(key)[0] == Particles.template getProp<dVector>(key)[0];
match &= Particles.template getProp<sVector>(key)[1]*Particles.template getProp<sVector>(key)[1] == Particles.template getProp<dVector>(key)[1];
match &= Particles.template getProp<sVector>(key)[2]*Particles.template getProp<sVector>(key)[2] == Particles.template getProp<dVector>(key)[2];
++it3;
}
//THERE IS A BUG HERE IT IS SUMMING THE VECTORS.
BOOST_REQUIRE_EQUAL(match,true);
}
/*
TdxVx=Dx(sV[x]);
TT[0][0] = Dx(sV[x]);*/
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* DCPSE_OP_TEST_TEMPORAL_CPP_ */
#endif
#endif
\ No newline at end of file
......@@ -66,7 +66,7 @@ private:
openfpm::vector_custd<T> localEpsInvPow; // Each MPI rank has just access to the local ones
openfpm::vector_custd<T> calcKernels;
openfpm::vector<size_t> subsetKeyPid;
openfpm::vector_custd<size_t> subsetKeyPid;
vector_type & particles;
double rCut;
......@@ -432,7 +432,7 @@ public:
size_t localKey = subsetKeyPid.get(key.getKey());
double eps = localEps.get(localKey);
double epsInvPow = localEpsInvPow(localKey);
double epsInvPow = localEpsInvPow.get(localKey);
auto &particles = o1.getVector();
......@@ -482,6 +482,99 @@ public:
initializeStaticSize(particles, convergenceOrder, rCut, supportSizeFactor);
}
/*
*
Save computed DCPSE coefficients to file
*
*/
void save(const std::string &file){
auto & v_cl=create_vcluster();
size_t req = 0;
Packer<decltype(supportRefs),CudaMemory>::packRequest(supportRefs,req);
Packer<decltype(kerOffsets),CudaMemory>::packRequest(kerOffsets,req);
Packer<decltype(supportKeys1D),CudaMemory>::packRequest(supportKeys1D,req);
Packer<decltype(localEps),CudaMemory>::packRequest(localEps,req);
Packer<decltype(localEpsInvPow),CudaMemory>::packRequest(localEpsInvPow,req);
Packer<decltype(calcKernels),CudaMemory>::packRequest(calcKernels,req);
Packer<decltype(subsetKeyPid),CudaMemory>::packRequest(subsetKeyPid,req);
// allocate the memory
CudaMemory pmem;
ExtPreAlloc<CudaMemory> mem(req,pmem);
//Packing
Pack_stat sts;
Packer<decltype(supportRefs),CudaMemory>::pack(mem,supportRefs,sts);
Packer<decltype(kerOffsets),CudaMemory>::pack(mem,kerOffsets,sts);
Packer<decltype(supportKeys1D),CudaMemory>::pack(mem,supportKeys1D,sts);
Packer<decltype(localEps),CudaMemory>::pack(mem,localEps,sts);
Packer<decltype(localEpsInvPow),CudaMemory>::pack(mem,localEpsInvPow,sts);
Packer<decltype(calcKernels),CudaMemory>::pack(mem,calcKernels,sts);
Packer<decltype(subsetKeyPid),CudaMemory>::pack(mem,subsetKeyPid,sts);
// Save into a binary file
std::ofstream dump (file+"_"+std::to_string(v_cl.rank()), std::ios::out | std::ios::binary);
if (dump.is_open() == false)
{ std::cerr << __FILE__ << ":" << __LINE__ <<" Unable to write since dump is open at rank "<<v_cl.rank()<<std::endl;
return;
}
dump.write ((const char *)pmem.getPointer(), pmem.size());
return;
}
/*! \brief Load the DCPSE computations
*
*
*/
void load(const std::string & file)
{
auto & v_cl=create_vcluster();
std::ifstream fs (file+"_"+std::to_string(v_cl.rank()), std::ios::in | std::ios::binary | std::ios::ate );
if (fs.is_open() == false)
{
std::cerr << __FILE__ << ":" << __LINE__ << " error, opening file: " << file << std::endl;
return;
}
// take the size of the file
size_t sz = fs.tellg();
fs.close();
// reopen the file without ios::ate to read
std::ifstream input (file+"_"+std::to_string(v_cl.rank()), std::ios::in | std::ios::binary );
if (input.is_open() == false)
{//some message here maybe
return;}
// Create the CudaMemory memory
size_t req = 0;
req += sz;
CudaMemory pmem;
ExtPreAlloc<CudaMemory> mem(req,pmem);
mem.allocate(pmem.size());
// read
input.read((char *)pmem.getPointer(), sz);
//close the file
input.close();
//Unpacking
Unpack_stat ps;
Unpacker<decltype(supportRefs),CudaMemory>::unpack(mem,supportRefs,ps);
Unpacker<decltype(kerOffsets),CudaMemory>::unpack(mem,kerOffsets,ps);
Unpacker<decltype(supportKeys1D),CudaMemory>::unpack(mem,supportKeys1D,ps);
Unpacker<decltype(localEps),CudaMemory>::unpack(mem,localEps,ps);
Unpacker<decltype(localEpsInvPow),CudaMemory>::unpack(mem,localEpsInvPow,ps);
Unpacker<decltype(calcKernels),CudaMemory>::unpack(mem,calcKernels,ps);
Unpacker<decltype(subsetKeyPid),CudaMemory>::unpack(mem,subsetKeyPid,ps);
}
private:
template <typename U>
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment