From 4e5378942624afd9fad5da0520cd10278b3a1273 Mon Sep 17 00:00:00 2001 From: Pietro Incardona <i-bird@private-incardon-3.mpi-cbg.de> Date: Thu, 10 Mar 2016 06:49:22 -0500 Subject: [PATCH] Adding missing files --- src/PSE/Kernels_test_util.hpp | 161 +++++++++++++++++++++++++++++++++ src/PSE/Kernels_unit_tests.hpp | 95 +++++++++++++++++++ src/test_data/PSE_convergence | Bin 0 -> 2504 bytes src/unit_test_init_cleanup.hpp | 24 +++++ 4 files changed, 280 insertions(+) create mode 100644 src/PSE/Kernels_test_util.hpp create mode 100644 src/PSE/Kernels_unit_tests.hpp create mode 100644 src/test_data/PSE_convergence create mode 100644 src/unit_test_init_cleanup.hpp diff --git a/src/PSE/Kernels_test_util.hpp b/src/PSE/Kernels_test_util.hpp new file mode 100644 index 00000000..81101b9c --- /dev/null +++ b/src/PSE/Kernels_test_util.hpp @@ -0,0 +1,161 @@ +/* + * 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>, CartDecomposition<1,T> > vd(Npart,box,bc,g); + + auto it2 = vd.getIterator(); + + while (it2.isNext()) + { + auto key = it2.get(); + + // set the position of the particles + vd.template getPos<0>(key)[0] = 0.448000 - ((long int)key.getKey() - mp) * spacing; + //set the property of the particles + vd.template getProp<0>(key) = f_xex2(vd.template getPos<0>(key)[0]); + + ++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.template getPos<0>(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 getNNIterator<NO_CHECK>(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()) + { + // W(x-y) + T ker = lker.value(p,vd.template getPos<0>(nnp)); + + // 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_ */ diff --git a/src/PSE/Kernels_unit_tests.hpp b/src/PSE/Kernels_unit_tests.hpp new file mode 100644 index 00000000..80e9b105 --- /dev/null +++ b/src/PSE/Kernels_unit_tests.hpp @@ -0,0 +1,95 @@ +/* + * Kernels_unit_tests.hpp + * + * Created on: Feb 17, 2016 + * Author: i-bird + */ + +#ifndef OPENFPM_NUMERICS_SRC_PSE_KERNELS_UNIT_TESTS_HPP_ +#define OPENFPM_NUMERICS_SRC_PSE_KERNELS_UNIT_TESTS_HPP_ + +#include "PSE/Kernels_test_util.hpp" +#include <boost/multiprecision/float128.hpp> + +BOOST_AUTO_TEST_SUITE( pse_kernels_unit_tests ) + +BOOST_AUTO_TEST_CASE( pse_ker ) +{ + Vcluster & v_cl = *global_v_cluster; + + // This test is not made to run in parallel + if (v_cl.getProcessingUnits() > 1) + return; + + openfpm::vector<openfpm::vector<double>> y; + openfpm::vector<openfpm::vector<double>> y_res; + + // Load the result of the test + y_res.load("test_data/PSE_convergence"); + + // Every time increase the number of particles by 2 + for (size_t i = 250 ; i <= 2097152000 ; i*=2) + { + y.add(); + + PSEError err; + + /////// Order 2 ////////////// + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,2>>(i,2,err); + y.last().add(err.linf_error); + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,2>>(i,4,err); + y.last().add(err.linf_error); + + PSE_test<double,Lap_PSE<1,double,2>>(i,2,err); + y.last().add(err.linf_error); + + PSE_test<double,Lap_PSE<1,double,2>>(i,4,err); + y.last().add(err.linf_error); + + PSE_test<float,Lap_PSE<1,float,2>>(i,2,err); + y.last().add(err.linf_error); + + PSE_test<float,Lap_PSE<1,float,2>>(i,4,err); + y.last().add(err.linf_error); + + //////// Order 4 ///////////// + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,4>>(i,2,err); + y.last().add(err.linf_error); + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,4>>(i,4,err); + y.last().add(err.linf_error); + + + //////// Order 6 ///////////// + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,6>>(i,2,err); + y.last().add(err.linf_error); + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,6>>(i,4,err); + y.last().add(err.linf_error); + + //////// Order 8 ///////////// + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,8>>(i,8,err); + y.last().add(err.linf_error); + + PSE_test<boost::multiprecision::float128,Lap_PSE<1,boost::multiprecision::float128,8>>(i,16,err); + y.last().add(err.linf_error); + } + + // Check the result + for (size_t i = 0 ; i < y.size(); i++) + { + for (size_t j = 0 ; j < y.get(i).size(); j++) + { + BOOST_REQUIRE_CLOSE(y.get(i).get(j),y_res.get(i).get(j),0.01); + } + } +} + +BOOST_AUTO_TEST_SUITE_END() + +#endif /* OPENFPM_NUMERICS_SRC_PSE_KERNELS_UNIT_TESTS_HPP_ */ diff --git a/src/test_data/PSE_convergence b/src/test_data/PSE_convergence new file mode 100644 index 0000000000000000000000000000000000000000..9ae0642f606ceb86fe67f80d65a87f9867f1185e GIT binary patch literal 2504 zcmX}tc~leE9tZFnsv-&q8pKisLa0DJrC2lpW{nD{PzfsNQL2f8qJ6R`_W3{{K`@{| zK%pXwhD}?nqP1>4A_*W-3yOszWqAe>Ns+P$&jnfDy%Wg%^F6=GoZlzEd+wcC1az12 zRGLJcjw`ShQ-3!X9oT9vrNI7EObx21YV_PDF@?fqzLZj&ui6oRPerDd57dZDSQLBE zjkO`3LnZ~~zmf+^$zrlK!6=ST?Yr-@z~JCbYTDWD)zdh{5l|wL&I7_XObv3UJOs>N z!#FhbUX*n$HM=h*@r2Ti$TK?Z8L%$AieRz%%TLE}Nc+QM2G4^4*=3%&ue{S$?|kQZ z)Jko~p}ar5()Cw7{1nwd(XOkJ<Lt&b;O|w`UU5==Wvzfr%c~BZ2faCwU;NMJT>_E2 zr>?OpQAL=y<(&LAMM6w_)AOXIxrdM0v8dGfF^mU*$re-t&ErLk2N_AFzs_lPqh5EU zb6eFcGADlV=ym9sb!AdK&sajZrghX!CMG`HR`3jD1HCuxq!(YTG2(WT#nmeq4}jjc zhOh?UydUF1hi-_{IT9Z#e5R{#R-=l@-v71366o0<!xv`R#uC?uQZ_%9v51KKY1{ui zBf+hvylLm6C2nQHRgqD#S&Rq3Z26d-3OfBEvl$P7+i`h}BXT}My{K3691{}=ZSN<C z_-Yjq-QT}0!6ue)ygqViSXPM7evy0nL=M2sYIWYU^Q7*`pb58-tdk6UiaY=`7kups zVY(!LU_1cwi|#OvFfWq|8sa)PWq(7&&t=Ub)l#B1?jP?R&^g-H#u}Fia3EP#e=A3V zwan3Z)6Rc#LYz#wQgY~^*8uY1GF@57lfsnL`|4`I`K^O->`(nh36sWajGEnv3(u8h zlv+SsiY=@OfR63bvzsSnofrV3i~~6U0}+*V$eVVm+$=p!xnffOb-q3Fpu<-Bi;z>b z&y1-x*D4HXc%P2vy`btdA4OhK@`;i5!3DIMMO@BnD%}Jfu$9w0DU(76*d8&7_W(dB z3x?jblOE6^HRTG(i$z1rkq4hs7H}<4HrburiIOhZ`vJ|Pv`BqTZTQwDWb1i8aq$T5 zOsiFROT3dc;{e|nqa+S;!2WCBKnIADmc)5XJGteDL+5ch<h2Q}mB@p?JB@coqRjZ+ z!waRx)y$QKDG9n^LZ{9;+rO%SPmHfjUO}s+_}mt9Y$KWn0E-U)8>qjVhp!wZaCZpR zJIBm#P0r`C$m;=%d8W{Ve|uisnS9h==6E_9r6M@Xfrcq*)LPLsio2N{=qMU)E~nK3 zJd0N|JmoaG6=$LT)?ePch*D^mlA(9f&Wvp5a65@Zm;9l_Le4qwcm?X;rM#JjQfXPS zi-xIryByF1MVZky&}sPfKrpFh;p39hA)TX??6}Z?`uAog)uPm92Vd1Y3oqZU+s0K9 zkHAjobfWnseXVzpx9ojwDRO|gAcIEtkEU;--&a;_MUEsz_;P2h3gdMxo^l%4$07~U zej8H(Ep#OABDY6+XC&^foVG6{u5R819Tjxee$Z$eQUBR>err%l$}KC9qvYXdXuKjd z5II2GQPz?z#Y}(xIH;Sa)@J@wv=1CV%wikCw4EV8q18LAeGl`VIY^0@cDtcd2%Xjv zzxHm}KY8`*ot?})xz!!W(a?z%^p)}godcqK*Tzy<i1_`K<I;Nkry)2;l&uAC5p)J8 z8|o+Z&L1glzq|6p#NO}3(2+u?{q|x*_`U||2NRR~kfT|%WCA(Lw0}D2AnN#LW*#t3 zKVa>^6!z6R$JWp=MzII>Hw=`T2w~1WYC*W@&re&Q(}|LB0TKO9039)O1i6m~ec*Q? z_vxX)Wz4zxXs|>89mPt!33}teD<8&jbiMH;M~ZFUz?z!p_=3LClNivnw5i{RL&xf9 zUZha(bWao~?+vgc8opnN&JO|bW}7PaZz=Gjd)~9sb;zMDpLifg)XiQZD1Ya!o9D;X zfe%XM*uPG&b^G7}@0as?Ph%h>;&WL|A_i==cQgb0^v+ySieFUUOk&Y+4tyVQ2!LHy z-XD#iz}g2hfA$JNj>vx19^^>8LcIiY%<JXMJaojJUz-1oy{0z3W%fZCX)W_F!a!3) zPf1oL2Fzl{zKE6So&MJjo<8<(@IOQ@ohDosK_}T<*{Uq3l+U{h+LoST*3I27206BG z26qXXWn5=#wBJtlPY?^PVK4o^YG9m^z@)WdN(_kR7~Si=fC2lz?AHHvUhnwGfB0%e zSPzb`$k=JX?Zj_O&o#T|D=EoSqn<PTOU$|{)`v2VS%0*E?fQcuaum*KRmsJ@n4l?4 ziyTGBI?J$EBN%Y1Bn~JVF(7_e?VWH_@04qKLwx>A{F7zcyOkUjKH|a7AO@dPEt=il zS8u7Ac@A4-(rEq8LGOuiaL6L)gg#PiZ%&$ri^TcskrUK(T>kkHb6nW$8ZECf!h!MZ z<7HNTdWW7d{C(E0E}Zl3u}Rd-LfmXPaIVn-6Z+I=*)H~rwaBq8j-FuV%9)WtduNVq vDbD8;8sGFXLluU1t)ATB4ISgGe9eFFufg|}eC!!L>x{QM`lf^#v+@4}KRI2^ literal 0 HcmV?d00001 diff --git a/src/unit_test_init_cleanup.hpp b/src/unit_test_init_cleanup.hpp new file mode 100644 index 00000000..12b98f01 --- /dev/null +++ b/src/unit_test_init_cleanup.hpp @@ -0,0 +1,24 @@ +/* + * unit_test_init_cleanup.hpp + * + * Created on: Apr 17, 2015 + * Author: Pietro Incardona + */ + +#ifndef UNIT_TEST_INIT_CLEANUP_HPP_ +#define UNIT_TEST_INIT_CLEANUP_HPP_ + +#include "VCluster.hpp" + +struct ut_start { + ut_start() { BOOST_TEST_MESSAGE("Initialize global VCluster"); init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv); } + ~ut_start() { BOOST_TEST_MESSAGE("Delete global VClster");delete_global_v_cluster(); } +}; + +//____________________________________________________________________________// + +BOOST_GLOBAL_FIXTURE( ut_start ); + + + +#endif /* UNIT_TEST_INIT_CLEANUP_HPP_ */ -- GitLab