From 86fa90915444a7c4b9571a27978f5e11724ba2a0 Mon Sep 17 00:00:00 2001 From: hinder Date: Wed, 9 Oct 2013 17:56:00 +0000 Subject: Add convergence test for Simpson integration method git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@94 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843 --- interface.ccl | 5 ++ schedule.ccl | 9 ++ src/make.code.defn | 2 +- src/tests.cc | 97 ++++++++++++++++++++++ test/test_simpson.par | 5 +- .../test_simpson_convergence_order..asc | 4 + 6 files changed, 120 insertions(+), 2 deletions(-) create mode 100644 src/tests.cc create mode 100644 test/test_simpson/test_simpson_convergence_order..asc diff --git a/interface.ccl b/interface.ccl index a8eda59..6c12ff7 100644 --- a/interface.ccl +++ b/interface.ccl @@ -13,3 +13,8 @@ CCTK_REAL harmonics type=GF timelevels=1 { harmonic_re, harmonic_im } "Spherical harmonics" + +CCTK_REAL test_integration_convergence_orders type=SCALAR +{ + test_simpson_convergence_order +} "Test integration convergence orders" diff --git a/schedule.ccl b/schedule.ccl index 9c94f21..d460e0f 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,5 +1,7 @@ #schedule.ccl for thorn Multipole +STORAGE: test_integration_convergence_orders + if (enable_test) { STORAGE: harmonics[1] @@ -24,3 +26,10 @@ schedule Multipole_ParamCheck at CCTK_PARAMCHECK LANG: C OPTIONS: GLOBAL } "Check Multipole parameters" + +# # Tests + +schedule Multipole_TestSimpson at CCTK_PARAMCHECK +{ + LANG: C +} "Test simpson integration method" diff --git a/src/make.code.defn b/src/make.code.defn index cccefbd..4d49656 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -1 +1 @@ -SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc +SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc tests.cc diff --git a/src/tests.cc b/src/tests.cc new file mode 100644 index 0000000..b63e18b --- /dev/null +++ b/src/tests.cc @@ -0,0 +1,97 @@ +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "utils.hh" +#include "integrate.hh" + +static CCTK_REAL test_simpson_integral(int n) +{ + const int nx = n; + const int ny = n; + const int array_size=(nx+1)*(ny+1); + + CCTK_REAL *f = new CCTK_REAL[array_size]; + + const CCTK_REAL dx = 1./nx; + const CCTK_REAL dy = 1./ny; + + for (int ix = 0; ix <= nx; ix++) + { + for (int iy = 0; iy <= ny; iy++) + { + const int i = Multipole_Index(ix, iy, nx); + + const CCTK_REAL x = ix*dx; + const CCTK_REAL y = iy*dy; + const CCTK_REAL PI = acos(-1.0); + + f[i] = x*pow(y,2)*pow(cos(2*PI*y),2)*pow(sin(2*PI*x),2); + } + } + + const CCTK_REAL result = Simpson2DIntegral(f, nx, ny, dx, dy); + delete [] f; + return result; +} + + +void Multipole_TestSimpson(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + + const int n1 = 100; + const int n2 = 200; + const CCTK_REAL PI = acos(-1.0); + const CCTK_REAL result1 = test_simpson_integral(100); + const CCTK_REAL result2 = test_simpson_integral(200); + const CCTK_REAL exact = 1./24 + 1./(64 * pow(PI,2)); + const CCTK_REAL error1 = fabs(result1 - exact); + const CCTK_REAL error2 = fabs(result2 - exact); + *test_simpson_convergence_order = log10(error1/error2) / log10((CCTK_REAL) n2/n1); +} + +// void Multipole_TestIntegrate(CCTK_ARGUMENTS) +// { +// const int n = 100; + +// const int nth = n; +// const int nph = n; +// const int array_size=(nth+1)*(nph+1); + +// CCTK_REAL *f = new CCTK_REAL[array_size]; + +// const CCTK_REAL dth = 1/nx; +// const CCTK_REAL dph = 1/ny; + +// for (int ix = 0; ix <= nx; ix++) +// { +// for (int iy = 0; iy <= ny; iy++) +// { +// const int i = Multipole_Index(ix, iy, nx); + +// const CCTK_REAL x = ix*dx; +// const CCTK_REAL y = iy*dy; + +// f[i] = sin(2*PI*x)*cos(2*PI*y); +// } +// } + + + + +// CCTK_REAL result = Multipole_Integrate(int array_size, int nthetap, +// CCTK_REAL const array1r[], CCTK_REAL const array1i[], +// CCTK_REAL const array2r[], CCTK_REAL const array2i[], +// CCTK_REAL const th[], CCTK_REAL const ph[], +// CCTK_REAL *outre, CCTK_REAL *outim) + + + +// printf("Integration result: %.19g\n", result); +// } diff --git a/test/test_simpson.par b/test/test_simpson.par index 78daaf4..772d2a9 100644 --- a/test/test_simpson.par +++ b/test/test_simpson.par @@ -1,5 +1,5 @@ -ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib CarpetInterp AEILocalInterp InitBase Multipole LoopControl" +ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib CarpetInterp AEILocalInterp InitBase Multipole LoopControl CarpetIOASCII" ############################################################# # Grid @@ -87,3 +87,6 @@ IO::out_fileinfo = "none" # CarpetIOASCII::out1d_x = yes # CarpetIOASCII::out1d_y = yes # CarpetIOASCII::out1d_z = yes + +CarpetIOASCII::out0d_vars = "Multipole::test_simpson_convergence_order" +CarpetIOASCII::out0d_every = 1 diff --git a/test/test_simpson/test_simpson_convergence_order..asc b/test/test_simpson/test_simpson_convergence_order..asc new file mode 100644 index 0000000..446999a --- /dev/null +++ b/test/test_simpson/test_simpson_convergence_order..asc @@ -0,0 +1,4 @@ +# 0D ASCII output created by CarpetIOASCII +# +0 0 0 0 0 0 0 0 0 0 0 0 4.00339319061966 + -- cgit v1.2.3