User: hinder Date: 2013/10/09 12:57 PM

Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ tests.cc

Log: Add tests for symmetry of integration

It is desirable that integrands with exact pi-antisymmetry integrate to zero; we add measures of this here. This is the case for all but the midpoint method.

File Changes:

Directory: /trunk/src/ ======================

File [modified]: tests.cc Delta lines: +46 -0 =================================================================== --- trunk/src/tests.cc 2013-10-09 17:57:06 UTC (rev 101) +++ trunk/src/tests.cc 2013-10-09 17:57:15 UTC (rev 102) @@ -40,6 +40,38 @@ return result; }

+ +static CCTK_REAL test_pi_symmetric_sphere_integral(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL)) +{ + const int n = 100; + const int nth = n; + const int nph = n; + const int array_size=(nth+1)*(nph+1); + const CCTK_REAL PI = acos(-1.0); + + CCTK_REAL *f = new CCTK_REAL[array_size]; + + const CCTK_REAL dth = PI/nth; + const CCTK_REAL dph = 2*PI/nph; + + for (int ith = 0; ith <= nth; ith++) + { + for (int iph = 0; iph <= nph; iph++) + { + const int i = Multipole_Index(ith, iph, nth); + + const CCTK_REAL th = ith*dth; + const CCTK_REAL ph = iph*dph; + + f[i] = -(cos(ph)*sqrt(5/PI)*pow(cos(th/2.),3)*sin(th/2.)); + } + } + + const CCTK_REAL result = integration_fn(f, nth, nph, dth, dph); + delete [] f; + return result; +} + CCTK_REAL integration_convergence_order(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL)) { const int n1 = 100; @@ -61,6 +93,20 @@ *test_midpoint_convergence_order = integration_convergence_order(&Midpoint2DIntegral); }

+void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + + *test_simpson_pi_symmetry = test_pi_symmetric_sphere_integral(&Simpson2DIntegral); + *test_midpoint_pi_symmetry = test_pi_symmetric_sphere_integral(&Midpoint2DIntegral); + *test_trapezoidal_pi_symmetry = test_pi_symmetric_sphere_integral(&Trapezoidal2DIntegral); + *test_driscollhealy_pi_symmetry = test_pi_symmetric_sphere_integral(&DriscollHealy2DIntegral); + printf("Pi symmetry Simpson integral: %.19g\n", *test_simpson_pi_symmetry); + printf("Pi symmetry midpoint integral: %.19g\n", *test_midpoint_pi_symmetry); + printf("Pi symmetry trapezoidal integral: %.19g\n", *test_trapezoidal_pi_symmetry); + printf("Pi symmetry Driscoll and Healy integral: %.19g\n", *test_driscollhealy_pi_symmetry); +} + // void Multipole_TestIntegrate(CCTK_ARGUMENTS) // { // const int n = 100;

Directory: /trunk/ ==================

File [modified]: interface.ccl Delta lines: +8 -0 =================================================================== --- trunk/interface.ccl 2013-10-09 17:57:06 UTC (rev 101) +++ trunk/interface.ccl 2013-10-09 17:57:15 UTC (rev 102) @@ -20,3 +20,11 @@ test_trapezoidal_convergence_order, test_simpson_convergence_order } "Test integration convergence orders" + +CCTK_REAL test_integration_symmetries type=SCALAR +{ + test_midpoint_pi_symmetry, + test_trapezoidal_pi_symmetry, + test_simpson_pi_symmetry, + test_driscollhealy_pi_symmetry +} "Test integration symmetries"

File [modified]: schedule.ccl Delta lines: +6 -0 =================================================================== --- trunk/schedule.ccl 2013-10-09 17:57:06 UTC (rev 101) +++ trunk/schedule.ccl 2013-10-09 17:57:15 UTC (rev 102) @@ -4,6 +4,7 @@ { STORAGE: harmonics[1] STORAGE: test_integration_convergence_orders + STORAGE: test_integration_symmetries }

schedule Multipole_Calc at CCTK_ANALYSIS after (calc_np,PsiKadelia,Accelerator_CopyBack) @@ -34,4 +35,9 @@ { LANG: C } "Test convergence of integration" + + schedule Multipole_TestIntegrationSymmetry at CCTK_PARAMCHECK + { + LANG: C + } "Test symmetry of integration" }