From ff2a71cdedcf7b522fe3f6f23fa02335a11f3fc3 Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Tue, 31 Jul 2018 18:06:57 -0500 Subject: [PATCH 1/5] Allocate more space for laplacian in spline test. The laplacian is stored as x, y, and z components during the spline evaluation and summed into the x component at the end. Need to allocate more space for it in the unit test. --- src/spline2/tests/gen_bspline_values.py | 2 +- src/spline2/tests/test_multi_spline.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/spline2/tests/gen_bspline_values.py b/src/spline2/tests/gen_bspline_values.py index e1b36d31c..45bdb9ef8 100644 --- a/src/spline2/tests/gen_bspline_values.py +++ b/src/spline2/tests/gen_bspline_values.py @@ -290,7 +290,7 @@ def evaluate_spline_3D(): print print - print ' Array lap(1);' + print ' Array lap(1,3);' print ' bs.evaluate_vgl(pos, v, dv, lap);' lap = diff(e,xs,2) + diff(e,ys,2) + diff(e,zs,2) diff --git a/src/spline2/tests/test_multi_spline.cpp b/src/spline2/tests/test_multi_spline.cpp index 61887d304..8a5eb6609 100644 --- a/src/spline2/tests/test_multi_spline.cpp +++ b/src/spline2/tests/test_multi_spline.cpp @@ -228,7 +228,7 @@ void test_splines() REQUIRE(hess(0, 5) == Approx( 34.53786329)); - Array lap(1); + Array lap(1,3); bs.evaluate_vgl(pos, v, dv, lap); // Value REQUIRE(v(0) == Approx( -0.9476393279)); From b23691b8ec2f1c32a68a5fb1b6c6c7acb3c57e9c Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Thu, 2 Aug 2018 11:04:00 -0500 Subject: [PATCH 2/5] Fix spline test for SoA and BGQ. Increase num_splines to ensure correct alignment on BG/Q. Rearrange the test accesses to gradients and hessians to match the SoA layout. Should fix #941 --- src/spline2/tests/gen_bspline_values.py | 30 ++++++++--------- src/spline2/tests/test_multi_spline.cpp | 43 +++++++++++++------------ 2 files changed, 38 insertions(+), 35 deletions(-) diff --git a/src/spline2/tests/gen_bspline_values.py b/src/spline2/tests/gen_bspline_values.py index 45bdb9ef8..f56f68de0 100644 --- a/src/spline2/tests/gen_bspline_values.py +++ b/src/spline2/tests/gen_bspline_values.py @@ -233,7 +233,7 @@ def evaluate_spline_3D(): print print ' // symbolic value at pos = ',e.subs(pos) print - print ' Array v(1);' + print ' Array v(num_splines);' 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 @@ -241,8 +241,8 @@ def evaluate_spline_3D(): print ' REQUIRE(v(0) == Approx(%15.10g));'%(val) print - print ' Array dv(1,3); // 3 - ndim' - print ' Array hess(1, 6); // 6 - number of unique hessian components' + print ' Array dv(3,num_splines); // 3 - ndim' + print ' Array hess(6, num_splines); // 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(%d,0) == Approx(%15.10g));'%(0,ex) + print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) + print ' REQUIRE(dv(%d,0) == 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(i,0) == Approx(0.0));' print ' }' @@ -278,19 +278,19 @@ def evaluate_spline_3D(): 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(%d,0) == Approx(%15.10g));'%(0,ex) + print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) + print ' REQUIRE(dv(%d,0) == 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(%d,0) == Approx(%15.10g));'%(idx, hess) #print d1,d2,hess print print - print ' Array lap(1,3);' + print ' Array lap(3,num_splines);' print ' bs.evaluate_vgl(pos, v, dv, lap);' lap = diff(e,xs,2) + diff(e,ys,2) + diff(e,zs,2) @@ -298,9 +298,9 @@ def evaluate_spline_3D(): print ' // Value' 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(%d,0) == Approx(%15.10g));'%(0,ex) + print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) + print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(2,ez) print ' // Laplacian' print ' REQUIRE(lap(0) == Approx(%15.10g));'%(lap_val) diff --git a/src/spline2/tests/test_multi_spline.cpp b/src/spline2/tests/test_multi_spline.cpp index 8a5eb6609..f56ed50fd 100644 --- a/src/spline2/tests/test_multi_spline.cpp +++ b/src/spline2/tests/test_multi_spline.cpp @@ -127,7 +127,7 @@ void test_splines() { MultiBspline bs; - int num_splines = 1; + int num_splines = 4; // must be multiple of 4 for alignment on BG/Q BCtype_d bc[3]; Ugrid grid[3]; @@ -179,9 +179,12 @@ void test_splines() bs.create(grid, bc, num_splines); - REQUIRE(bs.num_splines() == 1); + REQUIRE(bs.num_splines() == num_splines); + + for (int i = 0; i < num_splines; i++) { + bs.set(0, data); + } - bs.set(0, data); // Code from here to the end of the function is generated by gen_bspline_values.py @@ -189,21 +192,21 @@ 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 v(1); + Array v(num_splines); bs.evaluate(pos, v); REQUIRE(v(0) == Approx(-3.529930688e-12)); - Array dv(1,3); // 3 - ndim - Array hess(1, 6); // 6 - number of unique hessian components + Array dv(3,num_splines); // 3 - ndim + Array hess(6, num_splines); // 6 - number of unique hessian components 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(1,0) == Approx( -7.402942564)); + REQUIRE(dv(2,0) == Approx( -6.178320809)); // Hessian for (int i = 0; i < 6; i++) { - REQUIRE(hess(0,i) == Approx(0.0)); + REQUIRE(hess(i,0) == Approx(0.0)); } pos = {0.1, 0.2, 0.3}; @@ -217,25 +220,25 @@ void test_splines() 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(1,0) == Approx( 5.989106342)); + REQUIRE(dv(2,0) == 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(1,0) == Approx(1.174505743e-09)); + REQUIRE(hess(2,0) == Approx( -1.1483271e-09)); + REQUIRE(hess(3,0) == Approx( 133.9204891)); + REQUIRE(hess(4,0) == Approx(-2.15319293e-09)); + REQUIRE(hess(5,0) == Approx( 34.53786329)); - Array lap(1,3); + Array lap(3,num_splines); bs.evaluate_vgl(pos, v, dv, lap); // Value 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(1,0) == Approx( 5.989106342)); + REQUIRE(dv(2,0) == Approx( 1.952244379)); // Laplacian REQUIRE(lap(0) == Approx( 147.1127789)); From 76553cc13afe03f5270ca1b8f97de1db948c0985 Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Thu, 2 Aug 2018 14:19:30 -0500 Subject: [PATCH 3/5] Choose better box vectors for test. The previous box vectors would lead to a pathological result (rc=0), which would give an inverse cutoff of inf. On some platforms this becomes a NaN and causes test failures. --- src/Particle/tests/test_particle.cpp | 32 ++++++++++++++++++---------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/src/Particle/tests/test_particle.cpp b/src/Particle/tests/test_particle.cpp index e24ae94e9..66e74265e 100644 --- a/src/Particle/tests/test_particle.cpp +++ b/src/Particle/tests/test_particle.cpp @@ -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 ); } } From 796a3214a85fdaa8718e6ec0c87b828d55d25607 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 2 Aug 2018 22:26:01 -0500 Subject: [PATCH 4/5] Fix spline2 unit test. --- src/spline2/tests/test_multi_spline.cpp | 78 ++++++++++++------------- 1 file changed, 38 insertions(+), 40 deletions(-) diff --git a/src/spline2/tests/test_multi_spline.cpp b/src/spline2/tests/test_multi_spline.cpp index f56ed50fd..7d603265a 100644 --- a/src/spline2/tests/test_multi_spline.cpp +++ b/src/spline2/tests/test_multi_spline.cpp @@ -13,7 +13,7 @@ #define CATCH_CONFIG_MAIN #include "catch.hpp" -#include "OhmmsPETE/OhmmsArray.h" +#include #include "spline2/MultiBspline.hpp" namespace qmcplusplus @@ -127,7 +127,8 @@ void test_splines() { MultiBspline bs; - int num_splines = 4; // must be multiple of 4 for alignment on BG/Q + int num_splines = 1; + const int npad = getAlignedSize(num_splines); BCtype_d bc[3]; Ugrid grid[3]; @@ -163,27 +164,24 @@ void test_splines() double tpi = 2*M_PI; - Array data(N, N, N); + std::vector 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() == num_splines); - - for (int i = 0; i < num_splines; i++) { - bs.set(0, data); - } + for (int i = 0; i < num_splines; i++) + bs.set(i, data); // Code from here to the end of the function is generated by gen_bspline_values.py @@ -192,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 v(num_splines); + aligned_vector v(npad); bs.evaluate(pos, v); - REQUIRE(v(0) == Approx(-3.529930688e-12)); + REQUIRE(v[0] == Approx(-3.529930688e-12)); - Array dv(3,num_splines); // 3 - ndim - Array hess(6, num_splines); // 6 - number of unique hessian components + VectorSoaContainer dv(npad); + VectorSoaContainer hess(npad); bs.evaluate_vgh(pos, v, dv, hess); // Gradient - REQUIRE(dv(0,0) == Approx( 6.178320809)); - REQUIRE(dv(1,0) == Approx( -7.402942564)); - REQUIRE(dv(2,0) == 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(i,0) == 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(1,0) == Approx( 5.989106342)); - REQUIRE(dv(2,0) == 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(1,0) == Approx(1.174505743e-09)); - REQUIRE(hess(2,0) == Approx( -1.1483271e-09)); - REQUIRE(hess(3,0) == Approx( 133.9204891)); - REQUIRE(hess(4,0) == Approx(-2.15319293e-09)); - REQUIRE(hess(5,0) == 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 lap(3,num_splines); + VectorSoaContainer 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(1,0) == Approx( 5.989106342)); - REQUIRE(dv(2,0) == 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)); } From 6d30d418ca7d9782dbfde9f67cf1f0d3b0d27915 Mon Sep 17 00:00:00 2001 From: Mark Dewing Date: Fri, 3 Aug 2018 10:17:05 -0500 Subject: [PATCH 5/5] Update bspline test generator script. Move the VectorSoaContainer changes from the test to the generation script. --- src/spline2/tests/gen_bspline_values.py | 40 ++++++++++++------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/spline2/tests/gen_bspline_values.py b/src/spline2/tests/gen_bspline_values.py index f56f68de0..956f773ec 100644 --- a/src/spline2/tests/gen_bspline_values.py +++ b/src/spline2/tests/gen_bspline_values.py @@ -233,16 +233,16 @@ def evaluate_spline_3D(): print print ' // symbolic value at pos = ',e.subs(pos) print - print ' Array v(num_splines);' + print ' aligned_vector 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 dv(3,num_splines); // 3 - ndim' - print ' Array hess(6, num_splines); // 6 - number of unique hessian components' + print ' VectorSoaContainer dv(npad);' + print ' VectorSoaContainer 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(%d,0) == Approx(%15.10g));'%(0,ex) - print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) - print ' REQUIRE(dv(%d,0) == 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(i,0) == 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(%d,0) == Approx(%15.10g));'%(0,ex) - print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) - print ' REQUIRE(dv(%d,0) == 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(%d,0) == Approx(%15.10g));'%(idx, hess) + print ' REQUIRE(hess[0][%d] == Approx(%15.10g));'%(idx, hess) #print d1,d2,hess print print - print ' Array lap(3,num_splines);' + print ' VectorSoaContainer 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(%d,0) == Approx(%15.10g));'%(0,ex) - print ' REQUIRE(dv(%d,0) == Approx(%15.10g));'%(1,ey) - print ' REQUIRE(dv(%d,0) == 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__':