diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt index 787f6816c9..222e686217 100644 --- a/source/module_base/test/CMakeLists.txt +++ b/source/module_base/test/CMakeLists.txt @@ -1,4 +1,5 @@ remove_definitions(-D__MPI) +install(DIRECTORY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( TARGET base_blas_connector LIBS ${math_libs} diff --git a/source/module_base/test/data/bjo.bin b/source/module_base/test/data/bjo.bin new file mode 100644 index 0000000000..39c137c618 Binary files /dev/null and b/source/module_base/test/data/bjo.bin differ diff --git a/source/module_base/test/data/bjxo.bin b/source/module_base/test/data/bjxo.bin new file mode 100644 index 0000000000..0a93b2997f Binary files /dev/null and b/source/module_base/test/data/bjxo.bin differ diff --git a/source/module_base/test/math_sphbes_test.cpp b/source/module_base/test/math_sphbes_test.cpp index a56892f3d6..fd89950414 100644 --- a/source/module_base/test/math_sphbes_test.cpp +++ b/source/module_base/test/math_sphbes_test.cpp @@ -1,5 +1,6 @@ #include"../math_sphbes.h" #include +#include #ifdef __MPI #include"mpi.h" @@ -8,7 +9,7 @@ #include"gtest/gtest.h" #define doublethreshold 1e-7 - +#define MAX(x, y) ((x) > (y) ? (x) : (y)) /************************************************ * unit test of class Integral @@ -197,48 +198,139 @@ TEST_F(Sphbes,SphericalBesselRoots) } -TEST_F(Sphbes, SeriesAndRecurrence) +// TEST_F(Sphbes, SeriesAndRecurrence) +// { +// // This test checks whether Spherical_Bessel and sphbesj agree with each other +// // on a coarse grid for a range of l and q values. +// // +// // NOTE: this test should be removed once Spherical_Bessel is removed from the code. +// int lmax = 8; +// int nr = 5000; +// double rcut = 50; +// double dr = rcut / (nr - 1); +// double* r = new double[nr]; +// for (int i = 0; i < nr; ++i) +// { +// r[i] = i * dr; +// } + +// double* jl_old = new double[nr]; +// double* jl_new = new double[nr]; +// for (int l = 0; l <= lmax; ++l) +// { +// for (double q = 0.1; q < 10; q += 0.1) +// { +// ModuleBase::Sphbes::Spherical_Bessel(nr, r, q, l, jl_old); +// ModuleBase::Sphbes::sphbesj(nr, r, q, l, jl_new); +// for (int i = 0; i < nr; ++i) +// { +// EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); +// } + +// // derivative +// ModuleBase::Sphbes::dSpherical_Bessel_dx(nr, r, q, l, jl_old); +// ModuleBase::Sphbes::dsphbesj(nr, r, q, l, jl_new); +// for (int i = 0; i < nr; ++i) +// { +// EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); +// } +// } +// } + +// delete[] r; +// delete[] jl_old; +// delete[] jl_new; +// } + +TEST_F(Sphbes, SphericalBesselPrecisionGrid) { - // This test checks whether Spherical_Bessel and sphbesj agree with each other + // This test checks whether sphbesj agrees with the Octave implementation // on a coarse grid for a range of l and q values. - // - // NOTE: this test should be removed once Spherical_Bessel is removed from the code. - int lmax = 8; - int nr = 5000; - double rcut = 50; - double dr = rcut / (nr - 1); - double* r = new double[nr]; - for (int i = 0; i < nr; ++i) + const double q = 0.1; + const double rcut = 50; + const int nr = 5000; + const int l_lo = 0, l_hi = 11; + const double dr = rcut / nr; + double* r = new double[nr + 10]; + double* Y = new double[nr * (l_hi - l_lo + 1) + 10]; + + // case 0: x = 0, l = 0 + EXPECT_NEAR(ModuleBase::Sphbes::sphbesj(0,0), 1.0, 1e-12); + // case 1: x = 0, l = 1, ... , 11 + for (int l = 1; l <= l_hi; ++l) { - r[i] = i * dr; + EXPECT_NEAR(ModuleBase::Sphbes::sphbesj(l,0), 0.0, 1e-12); + } + // case 2: wide range of x and l + double y; + // read reference data + std::ifstream fin("data/bjo.bin", std::ios::binary); + int i = 0; + while (fin.read(reinterpret_cast(&y), sizeof(double))){ + Y[i] = y; ++i; } + fin.close(); - double* jl_old = new double[nr]; - double* jl_new = new double[nr]; - for (int l = 0; l <= lmax; ++l) + for(int i = 0; i < nr; ++i){ + r[i] = (i + 1) * dr; + } + + // test for new sphbesj + for (int l = l_lo; l <= l_hi; ++l) { - for (double q = 0.1; q < 10; q += 0.1) - { - ModuleBase::Sphbes::Spherical_Bessel(nr, r, q, 0, jl_old); - ModuleBase::Sphbes::sphbesj(nr, r, q, 0, jl_new); - for (int i = 0; i < nr; ++i) - { - EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); - } - - // derivative - ModuleBase::Sphbes::dSpherical_Bessel_dx(nr, r, q, 0, jl_old); - ModuleBase::Sphbes::dsphbesj(nr, r, q, 0, jl_new); - for (int i = 0; i < nr; ++i) - { - EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); - } + for(int i = 0; i< nr; ++i){ + EXPECT_NEAR(ModuleBase::Sphbes::sphbesj(l, r[i] * q), Y[l * nr + i], 1e-12); } + } + // test for old Bessel + // most of l cases precision failed to achieve 1e-12 + double* jl_old = new double[nr + 10]; + for (int l = l_lo; l <= l_hi; ++l) + { + ModuleBase::Sphbes::Spherical_Bessel(nr, r, q, l, jl_old); + double errs = 0.0; + for(int i = 0; i< nr; ++i){ + errs += MAX(jl_old[i] - Y[l * nr + i], Y[l * nr + i] - jl_old[i]); + } + // std::cout << " l = " << l << ", errors: " << std::scientific << errs / nr << std::endl; } - delete[] r; + delete[] Y; delete[] jl_old; - delete[] jl_new; +} + +TEST_F(Sphbes, SphericalBesselPrecisionNearZero) +{ + // This test checks whether sphbesj agrees with the Octave implementation + // when x is near zero point for a range of l. + const int n = 16; + const int l_lo = 0, l_hi = 11; + double* x = new double[n + 10]; + double* Y = new double[n * (l_hi - l_lo + 1) + 10]; + + // read reference data + double y; + std::ifstream fin("data/bjxo.bin", std::ios::binary); + int i = 0; + while (fin.read(reinterpret_cast(&y), sizeof(double))){ + Y[i] = y; ++i; + } + fin.close(); + // generate x + x[0] = 1.0 / (1 << 5); + for(int i = 1; i < n; i++){ + x[i] = x[i-1] / 2; + } + + // test for sphbesj near zero + for (int l = l_lo; l <= l_hi; ++l) + { + for(int i = 0; i< n; ++i){ + EXPECT_NEAR(ModuleBase::Sphbes::sphbesj(l, x[i]), Y[l * n + i], 1e-12); + } + } + delete[] x; + delete[] Y; } int main(int argc, char **argv) @@ -266,4 +358,4 @@ TEST_F(Sphbes,SphericalBesselsjp) EXPECT_EQ(sjp[iii], 1.0); } delete [] sjp; -} +} \ No newline at end of file