isolation.cu

#include <iostream>
#include <thread>

size_t debug_tot_call = 0;

#define PRINT_STACKTRACE
#define CHECKFOR_POSNAN
#define CHECKFOR_POSINF
#define CHECKFOR_PROPNAN
#define CHECKFOR_PROPINF

#define NO_WARNING
#include "Graph/CartesianGraphFactory.hpp"

void timeout_cycle()
{
	// 6 seconds
	std::this_thread::sleep_for (std::chrono::seconds(900));

	std::cout << "Time Out" << std::endl;
	std::exit(1);
}


#define BOOST_DISABLE_ASSERTS


#include "config.h"
#undef VERSION

#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "VCluster/VCluster.hpp"
#include <Vector/vector_dist.hpp>
#include "Vector/tests/vector_dist_util_unit_tests.hpp"

// initialization function:
bool init_unit_test()
{
//  std::thread to (timeout_cycle);
//  to.detach();
  return true;
}

// entry point
int main(int argc, char* argv[])
{
	openfpm_init(&argc,&argv);

  return boost::unit_test::unit_test_main( &init_unit_test, argc, argv );
}



BOOST_AUTO_TEST_CASE( vector_dist_ghost_put_gpu )
{
	

	Vcluster<> & v_cl = create_vcluster();

	long int k = 25*25*25*create_vcluster().getProcessingUnits();
	k = std::pow(k, 1/3.);

	if (v_cl.getProcessingUnits() > 48)
		return;

	BOOST_TEST_CHECKPOINT( "Testing 3D periodic ghost put k=" << k );

	long int big_step = k / 30;
	big_step = (big_step == 0)?1:big_step;
	long int small_step = 21;

	// 3D test
	for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
	{
		float r_cut = 1.3 / k;
		float r_g = 1.5 / k;

		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});

		// Boundary conditions
		size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

		// ghost
		Ghost<3,float> ghost(r_g);

		typedef  aggregate<float,float,float> part_prop;

		// Distributed vector
		vector_dist_gpu<3,float, part_prop > vd(0,box,bc,ghost);

		auto it = vd.getGridIterator({(size_t)k,(size_t)k,(size_t)k});

		while (it.isNext())
		{
			auto key = it.get();

			vd.add();

			vd.getLastPosWrite()[0] = key.get(0)*it.getSpacing(0);
			vd.getLastPosWrite()[1] = key.get(1)*it.getSpacing(1);
			vd.getLastPosWrite()[2] = key.get(2)*it.getSpacing(2);

			// Fill some properties randomly

			vd.getLastPropWrite<0>() = 0.0;

			vd.getLastPropWrite<2>() = 0.0;

			++it;
		}

		vd.map();

		vd.hostToDevicePos();
		vd.template hostToDeviceProp<0,2>();
		// sync the ghost
		vd.ghost_get<0,2>(RUN_ON_DEVICE);
		vd.template deviceToHostProp<0,2>();
		vd.deviceToHostPos();

		{
			auto NN = vd.getCellList(r_cut);
			float a = 1.0f*k*k;

			// run trough all the particles + ghost

			auto it2 = vd.getDomainIterator();

			while (it2.isNext())
			{
				// particle p
				auto p = it2.get();
				Point<3,float> xp = vd.getPos(p);

				// Get an iterator over the neighborhood particles of p
				auto Np = NN.getNNIterator<NO_CHECK>(NN.getCell(xp));

				// For each neighborhood particle ...
				while (Np.isNext())
				{
					auto q = Np.get();
					Point<3,float> xq = vd.getPosRead(q);

					float dist = xp.distance(xq);

					if (dist < r_cut)
					{
						vd.getPropWrite<0>(q) += a*(-dist*dist+r_cut*r_cut);
						vd.getPropWrite<2>(q) += a*(-dist*dist+r_cut*r_cut);
					}

					++Np;
				}

				++it2;
			}

			vd.hostToDevicePos();
			vd.template hostToDeviceProp<0,2>();
			vd.template ghost_put<add_atomic_,0,2>(RUN_ON_DEVICE);
			vd.template deviceToHostProp<0,2>();
			vd.deviceToHostPos();

			bool ret = true;
			auto it3 = vd.getDomainIterator();

			float constant = vd.getProp<0>(it3.get());
			float eps = 0.001;

			while (it3.isNext())
			{
				float constant2 = vd.getProp<0>(it3.get());
				float constant3 = vd.getProp<2>(it3.get());
				if (fabs(constant - constant2)/constant > eps || fabs(constant - constant3)/constant > eps)
				{
					Point<3,float> p = vd.getPosRead(it3.get());

					std::cout << p.toString() << "    " <<  constant2 << "/" << constant << "/" << constant3 << "    " << v_cl.getProcessUnitID() << std::endl;
					ret = false;
					break;
				}

				++it3;
			}
			BOOST_REQUIRE_EQUAL(ret,true);
		}

		auto itp = vd.getDomainAndGhostIterator();
		while (itp.isNext())
		{
			auto key = itp.get();

			vd.getPropWrite<0>(key) = 0.0;
			vd.getPropWrite<2>(key) = 0.0;

			++itp;
		}

		{
			auto NN = vd.getCellList(r_cut);
			float a = 1.0f*k*k;

			// run trough all the particles + ghost

			auto it2 = vd.getDomainIterator();

			while (it2.isNext())
			{
				// particle p
				auto p = it2.get();
				Point<3,float> xp = vd.getPosRead(p);

				// Get an iterator over the neighborhood particles of p
				auto Np = NN.getNNIterator<NO_CHECK>(NN.getCell(xp));

				// For each neighborhood particle ...
				while (Np.isNext())
				{
					auto q = Np.get();
					Point<3,float> xq = vd.getPosRead(q);

					float dist = xp.distance(xq);

					if (dist < r_cut)
					{
						vd.getPropWrite<0>(q) += a*(-dist*dist+r_cut*r_cut);
						vd.getPropWrite<2>(q) += a*(-dist*dist+r_cut*r_cut);
					}

					++Np;
				}

				++it2;
			}

			vd.hostToDevicePos();
			vd.template hostToDeviceProp<0,2>();
			vd.template ghost_put<add_atomic_,0>(RUN_ON_DEVICE);
			vd.template ghost_put<add_atomic_,2>(RUN_ON_DEVICE);
			vd.template deviceToHostProp<0,2>();
			vd.deviceToHostPos();

			bool ret = true;
			auto it3 = vd.getDomainIterator();

			float constant = vd.getPropRead<0>(it3.get());
			float eps = 0.001;

			while (it3.isNext())
			{
				float constant2 = vd.getPropRead<0>(it3.get());
				float constant3 = vd.getPropRead<0>(it3.get());
				if (fabs(constant - constant2)/constant > eps || fabs(constant - constant3)/constant > eps)
				{
					Point<3,float> p = vd.getPosRead(it3.get());

					std::cout << p.toString() << "    " <<  constant2 << "/" << constant << "/" << constant3 << "  " << it3.get().getKey() << "    " << v_cl.getProcessUnitID() << std::endl;
					ret = false;
					break;
				}

				++it3;
			}
			BOOST_REQUIRE_EQUAL(ret,true);
		}
	}

	openfpm_finalize();
}