Kernels_test_util.hpp

/*
 * Kernels_unit_tests.hpp
 *
 *  Created on: Feb 16, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_
#define OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_

#include "Kernels.hpp"
#include "Space/Ghost.hpp"
#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
#include "Decomposition/CartDecomposition.hpp"

struct PSEError
{
	double l2_error;
	double linf_error;
};


template<typename T> T f_xex2(T x)
{
	return x*exp(-x*x);
}

template<typename T> T f_xex2(Point<1,T> & x)
{
	return x.get(0)*exp(-x.get(0)*x.get(0));
}

template<typename T> T Lapf_xex2(Point<1,T> & x)
{
	return 2.0*x.get(0)*(2.0*x.get(0)*x.get(0) - 3.0)*exp(-x.get(0)*x.get(0));
}

/*
 * Given the Number of particles, it calculate the error produced by the standard
 * second order PSE kernel to approximate the laplacian of x*e^{-x^2} in the point
 * 0.5 (1D)
 *
 * The spacing h is calculated as
 *
 * \$ h = 1.0 / N_{part} \$
 *
 * epsilon of the 1D kernel is calculated as
 *
 * \$ \epsilon=overlap * h \$
 *
 * this mean that
 *
 * \$ overlap = \frac{\epsilon}{h} \$
 *
 * While the particles that are considered in the neighborhood of 0.5 are
 * calculated as
 *
 * \$ \N_contrib = 20*eps/spacing \$
 *
 * \param N_part Number of particles in the domain
 * \param overlap overlap parameter
 *
 *
 */
template<typename T, typename Kernel> void PSE_test(size_t Npart, size_t overlap,struct PSEError & error)
{
	// The domain
	Box<1,T> box({0.0},{4.0});

	// Calculated spacing
	T spacing = box.getHigh(0) / Npart;

	// Epsilon of the particle kernel
	const T eps = overlap*spacing;

	// Laplacian PSE kernel 1 dimension, on double, second order
	Kernel lker(eps);

	// For convergence we need less particles
	Npart = static_cast<size_t>(20*eps/spacing);

	// Middle particle
	long int mp = Npart / 2;

    size_t bc[1]={NON_PERIODIC};
	Ghost<1,T> g(20*eps);

	vector_dist<1,T, aggregate<T> > vd(Npart,box,bc,g);

	auto it2 = vd.getIterator();

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

		// set the position of the particles
		vd.getPos(key)[0] = 0.448000 - ((long int)key.getKey() - mp) * spacing;
		//set the property of the particles
		T d = vd.getPos(key)[0];

		vd.template getProp<0>(key) = f_xex2(d);

		++it2;
	}

	vect_dist_key_dx key;
	key.setKey(mp);

    // get and construct the Cell list
    auto cl = vd.getCellList(0.5);

    // Maximum infinity norm
    double linf = 0.0;

    // PSE integration accumulator
    T pse = 0;

    // Get the position of the particle
    Point<1,T> p = vd.getPos(key);

    // Get f(x) at the position of the particle
    T prp_x = vd.template getProp<0>(key);

    // Get the neighborhood of the particles
    auto NN = cl.template getNNIteratorBox(cl.getCell(p));
    while(NN.isNext())
    {
    	auto nnp = NN.get();

    	// Calculate contribution given by the kernel value at position p,
    	// given by the Near particle, exclude itself
    	if (nnp != key.getKey())
    	{
    		Point<1,T> xp = vd.getPos(nnp);

    		// W(x-y)
    		T ker = lker.value(p,xp);

    		// f(y)
    		T prp_y = vd.template getProp<0>(nnp);

    		// 1.0/(eps)^2 [f(y)-f(x)]*W(x,y)*V_q
    		T prp = 1.0/eps/eps * spacing * (prp_y - prp_x);
    		pse += prp * ker;
    	}

    	// Next particle
    	++NN;
    }

	// Here we calculate the L_infinity norm or the maximum difference
	// of the analytic solution from the PSE calculated

	T sol = Lapf_xex2(p);

	if (fabs(pse - sol) > linf)
		linf = static_cast<double>(fabs(pse - sol));


    error.linf_error = (double)linf;
}


#endif /* OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_ */