Skip to content

Commit 8751d9f

Browse files
pxlxingliangpxlxingliang
authored
Merge pull request #720 from pxlxingliang/math-sphbes
test: add the unit test and comments of math_sphbes.h
2 parents 2dbcd0c + 1ce08d0 commit 8751d9f

File tree

3 files changed

+223
-11
lines changed

3 files changed

+223
-11
lines changed

‎source/module_base/math_sphbes.h‎

Lines changed: 38 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -14,26 +14,53 @@ class Sphbes
1414
Sphbes();
1515
~Sphbes();
1616

17+
/**
18+
* @brief spherical bessel
19+
*
20+
* @param msh [in] number of grid points
21+
* @param r [in] radial grid
22+
* @param q [in] k_radial
23+
* @param l [in] angular momentum
24+
* @param jl [out] jl spherical bessel function
25+
*/
1726
static void Spherical_Bessel
1827
(
19-
const int &msh, //number of grid points
20-
const double *r,//radial grid
21-
const double &q, //
22-
const int &l, //angular momentum
23-
double *jl //jl(1:msh) = j_l(q*r(i)),spherical bessel function
28+
const int &msh,
29+
const double *r,
30+
const double &q,
31+
const int &l,
32+
double *jl
2433
);
2534

35+
/**
36+
* @brief spherical bessel
37+
*
38+
* @param msh [in] number of grid points
39+
* @param r [in] radial grid
40+
* @param q [in] k_radial
41+
* @param l [in] angular momentum
42+
* @param jl [out] jl spherical bessel function
43+
* @param sjp [out] sjp[i] is assigned to be 1.0. i < msh.
44+
*/
2645
static void Spherical_Bessel
2746
(
28-
const int &msh, //number of grid points
29-
const double *r,//radial grid
30-
const double &q, //
31-
const int &l, //angular momentum
32-
double *sj, //jl(1:msh) = j_l(q*r(i)),spherical bessel function
47+
const int &msh,
48+
const double *r,
49+
const double &q,
50+
const int &l,
51+
double *sj,
3352
double *sjp
3453
);
3554

36-
55+
/**
56+
* @brief return num eigenvalues of spherical bessel function
57+
*
58+
* @param num [in] the number of eigenvalues
59+
* @param l [in] angular number
60+
* @param epsilon [in] the accuracy
61+
* @param eigenvalue [out] the calculated eigenvalues
62+
* @param rcut [in] the cutoff the radial function
63+
*/
3764
static void Spherical_Bessel_Roots
3865
(
3966
const int &num,

‎source/module_base/test/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,10 @@ AddTest(
4747
TARGET base_sph_bessel_recursive
4848
SOURCES sph_bessel_recursive_test.cpp ../sph_bessel_recursive-d1.cpp ../sph_bessel_recursive-d2.cpp
4949
)
50+
AddTest(
51+
TARGET base_math_sphbes
52+
SOURCES math_sphbes_test.cpp ../math_sphbes.cpp ../timer.cpp
53+
)
5054
AddTest(
5155
TARGET base_realarray
5256
SOURCES realarray_test.cpp ../realarray.cpp

‎source/module_base/test/math_sphbes_test.cpp‎

Lines changed: 181 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,181 @@
1+
#include"../math_sphbes.h"
2+
#include<iostream>
3+
4+
#ifdef __MPI
5+
#include"mpi.h"
6+
#endif
7+
8+
#include"gtest/gtest.h"
9+
10+
#define doublethreshold 1e-7
11+
12+
13+
/************************************************
14+
* unit test of class Integral
15+
***********************************************/
16+
17+
/**
18+
* Note: this unit test try to ensure the invariance
19+
* of the spherical Bessel produced by class Sphbes,
20+
* and the reference results are produced by ABACUS
21+
* at 2022-1-27.
22+
*
23+
* Tested function:
24+
* - Spherical_Bessel.
25+
* - Spherical_Bessel_Roots
26+
*
27+
*/
28+
29+
double mean(const double* vect, const int totN)
30+
{
31+
double meanv = 0.0;
32+
for (int i=0; i< totN; ++i) {meanv += vect[i]/totN;}
33+
return meanv;
34+
}
35+
36+
class Sphbes : public testing::Test
37+
{
38+
protected:
39+
40+
int msh = 700;
41+
int l0 = 0;
42+
int l1 = 1;
43+
int l2 = 2;
44+
int l3 = 3;
45+
int l4 = 4;
46+
int l5 = 5;
47+
int l6 = 6;
48+
int l7 = 7;
49+
double q = 1.0;
50+
double *r = new double[msh];
51+
double *jl = new double[msh];
52+
53+
void SetUp()
54+
{
55+
for(int i=0; i<msh; ++i) {r[i] = 0.01*(i);}
56+
}
57+
58+
void TearDown()
59+
{
60+
delete [] r;
61+
delete [] jl;
62+
}
63+
};
64+
65+
66+
TEST_F(Sphbes,SphericalBessel)
67+
{
68+
//int l = 0;
69+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l0,jl);
70+
//reference result is from bessel_test.cpp which is calculated by
71+
//ModuleBase::Sph_Bessel_Recursive::D1
72+
EXPECT_NEAR(mean(jl,msh)/0.2084468748396, 1.0,doublethreshold);
73+
74+
75+
//int l = 1;
76+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l1,jl);
77+
//reference result is from bessel_test.cpp which is calculated by
78+
//ModuleBase::Sph_Bessel_Recursive::D1
79+
EXPECT_NEAR(mean(jl,msh)/0.12951635180384, 1.0,doublethreshold);
80+
81+
82+
//int l = 2;
83+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l2,jl);
84+
//the result from bessel_test.cpp calculated by
85+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.124201140093879
86+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
87+
EXPECT_NEAR(mean(jl,msh)/0.12420114009271901456, 1.0,doublethreshold);
88+
89+
//int l = 3;
90+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l3,jl);
91+
//the result from bessel_test.cpp calculated by
92+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.118268654505568
93+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
94+
EXPECT_NEAR(mean(jl,msh)/0.11826865448477598408, 1.0,doublethreshold);
95+
96+
97+
//int l = 4;
98+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l4,jl);
99+
//the result from bessel_test.cpp calculated by
100+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.0933871035384385
101+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
102+
EXPECT_NEAR(mean(jl,msh)/0.093387100084701621383, 1.0,doublethreshold);
103+
104+
105+
//int l = 5;
106+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l5,jl);
107+
//the result from bessel_test.cpp calculated by
108+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.0603800487910689
109+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
110+
EXPECT_NEAR(mean(jl,msh)/0.060380048719821471925, 1.0,doublethreshold);
111+
112+
113+
//int l = 6;
114+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l6,jl);
115+
//the result from bessel_test.cpp calculated by
116+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.0327117051555907
117+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
118+
EXPECT_NEAR(mean(jl,msh)/0.032711705053977857549, 1.0,doublethreshold);
119+
120+
121+
//int l = 7;
122+
ModuleBase::Sphbes::Spherical_Bessel(msh,r,q,l7,jl);
123+
//the result from bessel_test.cpp calculated by
124+
//ModuleBase::Sph_Bessel_Recursive::D1 is 0.0152155566653926
125+
//reference result is calculated by Sphbes::Spherical_Bessel(msh,r,q,l,jl)
126+
EXPECT_NEAR(mean(jl,msh)/0.015215556095798710851, 1.0,doublethreshold);
127+
}
128+
129+
TEST_F(Sphbes,SphericalBesselRoots)
130+
{
131+
int neign = 100;
132+
double **eign = new double*[8];
133+
for(int i=0;i<8;++i)
134+
{
135+
eign[i] = new double[neign];
136+
ModuleBase::Sphbes::Spherical_Bessel_Roots(neign,i,1.0e-12,eign[i],10.0);
137+
}
138+
139+
EXPECT_NEAR(eign[0][0]/0.31415926535899563188, 1.0,doublethreshold);
140+
EXPECT_NEAR(eign[0][99]/31.415926535896932847, 1.0,doublethreshold);
141+
EXPECT_NEAR(mean(eign[0],100)/15.865042900628463229, 1.0,doublethreshold);
142+
EXPECT_NEAR(eign[1][0]/0.44934094579091843347, 1.0,doublethreshold);
143+
EXPECT_NEAR(eign[1][99]/31.572689440204385392, 1.0,doublethreshold);
144+
EXPECT_NEAR(mean(eign[1],100)/16.020655759558295017, 1.0,doublethreshold);
145+
EXPECT_NEAR(eign[2][0]/0.57634591968946913276, 1.0,doublethreshold);
146+
EXPECT_NEAR(eign[2][99]/31.729140298172534784, 1.0,doublethreshold);
147+
EXPECT_NEAR(mean(eign[2],100)/16.175128483074864505, 1.0,doublethreshold);
148+
EXPECT_NEAR(eign[3][0]/0.69879320005004752492, 1.0,doublethreshold);
149+
EXPECT_NEAR(eign[3][99]/31.885283678838447941, 1.0,doublethreshold);
150+
EXPECT_NEAR(mean(eign[3],100)/16.328616567969248763, 1.0,doublethreshold);
151+
EXPECT_NEAR(eign[4][0]/0.81825614525711076741, 1.0,doublethreshold);
152+
EXPECT_NEAR(eign[4][99]/32.041124042016576823, 1.0,doublethreshold);
153+
EXPECT_NEAR(mean(eign[4],100)/16.481221742387987206, 1.0,doublethreshold);
154+
EXPECT_NEAR(eign[5][0]/0.93558121110426506473, 1.0,doublethreshold);
155+
EXPECT_NEAR(eign[5][99]/32.196665741899131774, 1.0,doublethreshold);
156+
EXPECT_NEAR(mean(eign[5],100)/16.633019118735202113, 1.0,doublethreshold);
157+
EXPECT_NEAR(eign[6][0]/1.051283540809391015, 1.0,doublethreshold);
158+
EXPECT_NEAR(eign[6][99]/32.351913030537232885, 1.0,doublethreshold);
159+
EXPECT_NEAR(mean(eign[6],100)/16.784067905062840964, 1.0,doublethreshold);
160+
EXPECT_NEAR(eign[7][0]/1.1657032192516516567, 1.0,doublethreshold);
161+
EXPECT_NEAR(eign[7][99]/32.506870061157627561, 1.0,doublethreshold);
162+
EXPECT_NEAR(mean(eign[7],100)/16.934416735327332049, 1.0,doublethreshold);
163+
164+
for(int i=0;i<8;++i) delete [] eign[i];
165+
delete [] eign;
166+
}
167+
168+
169+
int main(int argc, char **argv)
170+
{
171+
#ifdef __MPI
172+
MPI_Init(&argc, &argv);
173+
#endif
174+
175+
testing::InitGoogleTest(&argc, argv);
176+
int result = RUN_ALL_TESTS();
177+
#ifdef __MPI
178+
MPI_Finalize();
179+
#endif
180+
return result;
181+
}

0 commit comments

Comments
 (0)