[Commits] [svn:einsteintoolkit] Multipole/trunk/ (Rev. 94)

[ian.hinder at aei.mpg.de]

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

Added:
 /trunk/src/
  tests.cc
 /trunk/test/test_simpson/
  test_simpson_convergence_order..asc

Modified:
 /trunk/
  interface.ccl, schedule.ccl
 /trunk/src/
  make.code.defn
 /trunk/test/
  test_simpson.par

Log:
 Add convergence test for Simpson integration method

File Changes:

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

File [modified]: make.code.defn
Delta lines: +1 -1
===================================================================
--- trunk/src/make.code.defn	2013-10-09 17:55:48 UTC (rev 93)
+++ trunk/src/make.code.defn	2013-10-09 17:56:00 UTC (rev 94)
@@ -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

File [added]: tests.cc
Delta lines: +97 -0
===================================================================
--- trunk/src/tests.cc	                        (rev 0)
+++ trunk/src/tests.cc	2013-10-09 17:56:00 UTC (rev 94)
@@ -0,0 +1,97 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#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);
+// }

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

File [modified]: test_simpson.par
Delta lines: +4 -1
===================================================================
--- trunk/test/test_simpson.par	2013-10-09 17:55:48 UTC (rev 93)
+++ trunk/test/test_simpson.par	2013-10-09 17:56:00 UTC (rev 94)
@@ -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 @@
 # CarpetIOASCII::out1d_x = yes
 # CarpetIOASCII::out1d_y = yes
 # CarpetIOASCII::out1d_z = yes
+
+CarpetIOASCII::out0d_vars = "Multipole::test_simpson_convergence_order"
+CarpetIOASCII::out0d_every = 1

Directory: /trunk/test/test_simpson/
====================================

File [added]: test_simpson_convergence_order..asc
Delta lines: +4 -0
===================================================================
--- trunk/test/test_simpson/test_simpson_convergence_order..asc	                        (rev 0)
+++ trunk/test/test_simpson/test_simpson_convergence_order..asc	2013-10-09 17:56:00 UTC (rev 94)
@@ -0,0 +1,4 @@
+# 0D ASCII output created by CarpetIOASCII
+#
+0	0	0 0 0	0 0 0	0	0 0 0	4.00339319061966
+

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

File [modified]: interface.ccl
Delta lines: +5 -0
===================================================================
--- trunk/interface.ccl	2013-10-09 17:55:48 UTC (rev 93)
+++ trunk/interface.ccl	2013-10-09 17:56:00 UTC (rev 94)
@@ -13,3 +13,8 @@
 {
   harmonic_re, harmonic_im
 } "Spherical harmonics"
+
+CCTK_REAL test_integration_convergence_orders type=SCALAR
+{
+  test_simpson_convergence_order
+} "Test integration convergence orders"

File [modified]: schedule.ccl
Delta lines: +9 -0
===================================================================
--- trunk/schedule.ccl	2013-10-09 17:55:48 UTC (rev 93)
+++ trunk/schedule.ccl	2013-10-09 17:56:00 UTC (rev 94)
@@ -1,5 +1,7 @@
 #schedule.ccl for thorn Multipole
 
+STORAGE: test_integration_convergence_orders
+
 if (enable_test)
 {
   STORAGE: harmonics[1]
@@ -24,3 +26,10 @@
   LANG: C
   OPTIONS: GLOBAL
 } "Check Multipole parameters"
+
+# # Tests
+
+schedule Multipole_TestSimpson at CCTK_PARAMCHECK
+{
+  LANG: C
+} "Test simpson integration method"