scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #953 from markdewing/spline_lapl_space 

Fix BG/Q unit tests
This commit is contained in:
Paul R. C. Kent 2018-08-03 13:46:24 -04:00 committed by GitHub
parent 5b6fb4435e c97443ff0b
commit ad0d4e52b7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 78 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
src/Particle/tests/test_particle.cpp
View File

@ -66,8 +66,18 @@ TEST_CASE("particle set lattice with vacuum", "[particle]")
Uniform3DGridLayout grid;
// PPP case
grid.BoxBConds = true;
for(int i=0; i<9; i++)
grid.R(i) = i+1;
grid.R(0) = 1.0;
grid.R(1) = 2.0;
grid.R(2) = 3.0;
grid.R(3) = 0.0;
grid.R(4) = 1.0;
grid.R(5) = 0.0;
grid.R(6) = 0.0;
grid.R(7) = 0.0;
grid.R(8) = 1.0;
grid.VacuumScale=2.0;
grid.reset();
@ -85,9 +95,9 @@ TEST_CASE("particle set lattice with vacuum", "[particle]")
source.Lattice.copy(grid);
source.createSK();
REQUIRE( source.LRBox.R(2,0) == 14.0 );
REQUIRE( source.LRBox.R(2,1) == 16.0 );
REQUIRE( source.LRBox.R(2,2) == 18.0 );
REQUIRE( source.LRBox.R(2,0) == 0.0 );
REQUIRE( source.LRBox.R(2,1) == 0.0 );
REQUIRE( source.LRBox.R(2,2) == 2.0 );
// PNN case
grid.BoxBConds[1] = false;
@ -98,12 +108,12 @@ TEST_CASE("particle set lattice with vacuum", "[particle]")
REQUIRE( source.LRBox.R(0,0) == 1.0 );
REQUIRE( source.LRBox.R(0,1) == 2.0 );
REQUIRE( source.LRBox.R(0,2) == 3.0 );
REQUIRE( source.LRBox.R(1,0) == 8.0 );
REQUIRE( source.LRBox.R(1,1) == 10.0 );
REQUIRE( source.LRBox.R(1,2) == 12.0 );
REQUIRE( source.LRBox.R(2,0) == 14.0 );
REQUIRE( source.LRBox.R(2,1) == 16.0 );
REQUIRE( source.LRBox.R(2,2) == 18.0 );
REQUIRE( source.LRBox.R(1,0) == 0.0 );
REQUIRE( source.LRBox.R(1,1) == 2.0 );
REQUIRE( source.LRBox.R(1,2) == 0.0 );
REQUIRE( source.LRBox.R(2,0) == 0.0 );
REQUIRE( source.LRBox.R(2,1) == 0.0 );
REQUIRE( source.LRBox.R(2,2) == 2.0 );
}
}

40
src/spline2/tests/gen_bspline_values.py
View File

@ -233,16 +233,16 @@ def evaluate_spline_3D():
print
print ' // symbolic value at pos = ',e.subs(pos)
print
print ' Array<T, 1> v(1);'
print ' aligned_vector<T> v(npad);'
print ' bs.evaluate(pos, v);'
# the 'expand' after substituting pos is necessary for the triples of coefficients (cx*cy*cz) to
# be formed so the substitution of coefficients works
val = e.subs(pos).expand().subs(subslist)
print ' REQUIRE(v(0) == Approx(%15.10g));'%(val)
print ' REQUIRE(v[0] == Approx(%15.10g));'%(val)
print
print ' Array<T, 2> dv(1,3); // 3 - ndim'
print ' Array<T, 2> hess(1, 6); // 6 - number of unique hessian components'
print ' VectorSoaContainer<T,3> dv(npad);'
print ' VectorSoaContainer<T,6> hess(npad); // 6 - number of unique hessian components'
print ' bs.evaluate_vgh(pos, v, dv, hess);'
@ -250,14 +250,14 @@ def evaluate_spline_3D():
ey = diff(e, ys).subs(pos).expand().subs(subslist)
ez = diff(e, zs).subs(pos).expand().subs(subslist)
print ' // Gradient'
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(2,ez)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(2,ez)
print
print ' // Hessian'
print ' for (int i = 0; i < 6; i++) {'
print ' REQUIRE(hess(0,i) == Approx(0.0));'
print ' REQUIRE(hess[0][i] == Approx(0.0));'
print ' }'
@ -269,40 +269,40 @@ def evaluate_spline_3D():
val = e.subs(pos1).expand().subs(subslist)
print ' // Value'
print ' REQUIRE(v(0) == Approx(%15.10g));'%(val)
print ' REQUIRE(v[0] == Approx(%15.10g));'%(val)
print
print ' bs.evaluate_vgh(pos, v, dv, hess);'
print ' // Value'
print ' REQUIRE(v(0) == Approx(%15.10g));'%(val)
print ' REQUIRE(v[0] == Approx(%15.10g));'%(val)
ex = diff(e, xs).subs(pos1).expand().subs(subslist)
ey = diff(e, ys).subs(pos1).expand().subs(subslist)
ez = diff(e, zs).subs(pos1).expand().subs(subslist)
print ' // Gradient'
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(2,ez)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(2,ez)
print ' // Hessian'
for idx,(d1,d2) in enumerate([(xs,xs), (xs,ys), (xs,zs), (ys, ys), (ys, zs), (zs, zs)]):
hess = diff(diff(e,d1), d2).subs(pos1).expand().subs(subslist)
print ' REQUIRE(hess(0, %d) == Approx(%15.10g));'%(idx, hess)
print ' REQUIRE(hess[0][%d] == Approx(%15.10g));'%(idx, hess)
#print d1,d2,hess
print
print
print ' Array<T, 1> lap(1);'
print ' VectorSoaContainer<T,3> lap(npad);'
print ' bs.evaluate_vgl(pos, v, dv, lap);'
lap = diff(e,xs,2) + diff(e,ys,2) + diff(e,zs,2)
lap_val = lap.subs(pos1).expand().subs(subslist)
print ' // Value'
print ' REQUIRE(v(0) == Approx(%15.10g));'%(val)
print ' REQUIRE(v[0] == Approx(%15.10g));'%(val)
print ' // Gradient'
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv(0,%d) == Approx(%15.10g));'%(2,ez)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(0,ex)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(1,ey)
print ' REQUIRE(dv[0][%d] == Approx(%15.10g));'%(2,ez)
print ' // Laplacian'
print ' REQUIRE(lap(0) == Approx(%15.10g));'%(lap_val)
print ' REQUIRE(lap[0][0] == Approx(%15.10g));'%(lap_val)
if __name__ == '__main__':

73
src/spline2/tests/test_multi_spline.cpp
View File

@ -13,7 +13,7 @@
#define CATCH_CONFIG_MAIN
#include "catch.hpp"
#include "OhmmsPETE/OhmmsArray.h"
#include <OhmmsSoA/Container.h>
#include "spline2/MultiBspline.hpp"
namespace qmcplusplus
@ -128,6 +128,7 @@ void test_splines()
MultiBspline<T> bs;
int num_splines = 1;
const int npad = getAlignedSize<T>(num_splines);
BCtype_d bc[3];
Ugrid grid[3];
@ -163,25 +164,25 @@ void test_splines()
double tpi = 2*M_PI;
Array<T, 3> data(N, N, N);
std::vector<T> data(N*N*N);
// Generate the data in double precision regardless of the target spline precision
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
for (int k = 0; k < N; k++) {
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
{
double x = delta*i;
double y = delta*j;
double z = delta*k;
data(i,j,k) = std::sin(tpi*x) + std::sin(3*tpi*y) + std::sin(4*tpi*z);
}
}
data[i*N*N+j*N+k] = std::sin(tpi*x) + std::sin(3*tpi*y) + std::sin(4*tpi*z);
}
bs.create(grid, bc, npad);
bs.create(grid, bc, num_splines);
REQUIRE(bs.num_splines() == npad);
REQUIRE(bs.num_splines() == 1);
for (int i = 0; i < num_splines; i++)
bs.set(i, data);
bs.set(0, data);
// Code from here to the end of the function is generated by gen_bspline_values.py
@ -189,55 +190,55 @@ void test_splines()
// symbolic value at pos = (cx[0]/6 + 2*cx[1]/3 + cx[2]/6)*(cy[0]/6 + 2*cy[1]/3 + cy[2]/6)*(cz[0]/6 + 2*cz[1]/3 + cz[2]/6)
Array<T, 1> v(1);
aligned_vector<T> v(npad);
bs.evaluate(pos, v);
REQUIRE(v(0) == Approx(-3.529930688e-12));
REQUIRE(v[0] == Approx(-3.529930688e-12));
Array<T, 2> dv(1,3); // 3 - ndim
Array<T, 2> hess(1, 6); // 6 - number of unique hessian components
VectorSoaContainer<T,3> dv(npad);
VectorSoaContainer<T,6> hess(npad);
bs.evaluate_vgh(pos, v, dv, hess);
// Gradient
REQUIRE(dv(0,0) == Approx( 6.178320809));
REQUIRE(dv(0,1) == Approx( -7.402942564));
REQUIRE(dv(0,2) == Approx( -6.178320809));
REQUIRE(dv[0][0] == Approx( 6.178320809));
REQUIRE(dv[0][1] == Approx( -7.402942564));
REQUIRE(dv[0][2] == Approx( -6.178320809));
// Hessian
for (int i = 0; i < 6; i++) {
REQUIRE(hess(0,i) == Approx(0.0));
REQUIRE(hess[0][i] == Approx(0.0));
}
pos = {0.1, 0.2, 0.3};
bs.evaluate(pos, v);
// Value
REQUIRE(v(0) == Approx( -0.9476393279));
REQUIRE(v[0] == Approx( -0.9476393279));
bs.evaluate_vgh(pos, v, dv, hess);
// Value
REQUIRE(v(0) == Approx( -0.9476393279));
REQUIRE(v[0] == Approx( -0.9476393279));
// Gradient
REQUIRE(dv(0,0) == Approx( 5.111042137));
REQUIRE(dv(0,1) == Approx( 5.989106342));
REQUIRE(dv(0,2) == Approx( 1.952244379));
REQUIRE(dv[0][0] == Approx( 5.111042137));
REQUIRE(dv[0][1] == Approx( 5.989106342));
REQUIRE(dv[0][2] == Approx( 1.952244379));
// Hessian
REQUIRE(hess(0, 0) == Approx( -21.34557341));
REQUIRE(hess(0, 1) == Approx(1.174505743e-09));
REQUIRE(hess(0, 2) == Approx( -1.1483271e-09));
REQUIRE(hess(0, 3) == Approx( 133.9204891));
REQUIRE(hess(0, 4) == Approx(-2.15319293e-09));
REQUIRE(hess(0, 5) == Approx( 34.53786329));
REQUIRE(hess[0][0] == Approx( -21.34557341));
REQUIRE(hess[0][1] == Approx(1.174505743e-09));
REQUIRE(hess[0][2] == Approx( -1.1483271e-09));
REQUIRE(hess[0][3] == Approx( 133.9204891));
REQUIRE(hess[0][4] == Approx(-2.15319293e-09));
REQUIRE(hess[0][5] == Approx( 34.53786329));
Array<T, 1> lap(1);
VectorSoaContainer<T,3> lap(npad);
bs.evaluate_vgl(pos, v, dv, lap);
// Value
REQUIRE(v(0) == Approx( -0.9476393279));
REQUIRE(v[0] == Approx( -0.9476393279));
// Gradient
REQUIRE(dv(0,0) == Approx( 5.111042137));
REQUIRE(dv(0,1) == Approx( 5.989106342));
REQUIRE(dv(0,2) == Approx( 1.952244379));
REQUIRE(dv[0][0] == Approx( 5.111042137));
REQUIRE(dv[0][1] == Approx( 5.989106342));
REQUIRE(dv[0][2] == Approx( 1.952244379));
// Laplacian
REQUIRE(lap(0) == Approx( 147.1127789));
REQUIRE(lap[0][0] == Approx( 147.1127789));
}