From b49f04e6e2aadf21a9e98f558fd6ec5365d045ef Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Sun, 23 Jun 2024 14:30:07 +0000 Subject: [PATCH] Multipole: decorate sphericalharmonic --- Multipole/schedule.ccl | 6 +- Multipole/src/multipole.cc | 4 +- Multipole/src/sphericalharmonic.cxx | 96 +++++++---------------------- Multipole/src/sphericalharmonic.hxx | 13 ++-- Multipole/src/tests.cc | 4 +- 5 files changed, 37 insertions(+), 86 deletions(-) diff --git a/Multipole/schedule.ccl b/Multipole/schedule.ccl index 8a9fe758..b6d20b5b 100644 --- a/Multipole/schedule.ccl +++ b/Multipole/schedule.ccl @@ -34,7 +34,6 @@ if (enable_test) SCHEDULE Multipole_SetHarmonic AT CCTK_INITIAL { LANG: C - READS: CoordinatesX::vertex_coords WRITES: Multipole::harmonics(interior) SYNC: Multipole::harmonics } "Populate grid functions with spherical harmonics" @@ -51,9 +50,8 @@ if (enable_test_Weyl) SCHEDULE Multipole_SetHarmonicWeyl AT CCTK_INITIAL { LANG: C - READS: CoordinatesX::vertex_coords - WRITES: WeylScal4::Psi4i_group(Interior) - WRITES: WeylScal4::Psi4r_group(Interior) + WRITES: WeylScal4::Psi4i_group(interior) + WRITES: WeylScal4::Psi4r_group(interior) SYNC: WeylScal4::Psi4i_group SYNC: WeylScal4::Psi4r_group } "Populate Weyl grid functions with spherical harmonics" diff --git a/Multipole/src/multipole.cc b/Multipole/src/multipole.cc index 6b3b9a14..b0222e18 100644 --- a/Multipole/src/multipole.cc +++ b/Multipole/src/multipole.cc @@ -205,8 +205,8 @@ setup_harmonics(const int spin_weights[max_spin_weights], int n_spin_weights, reY[si][l][m + l] = new CCTK_REAL[array_size]; imY[si][l][m + l] = new CCTK_REAL[array_size]; - Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[si][l][m + l], - imY[si][l][m + l]); + Multipole::HarmonicSetup(sw, l, m, array_size, th, ph, + reY[si][l][m + l], imY[si][l][m + l]); } } } diff --git a/Multipole/src/sphericalharmonic.cxx b/Multipole/src/sphericalharmonic.cxx index 78b8e641..2ac85a65 100644 --- a/Multipole/src/sphericalharmonic.cxx +++ b/Multipole/src/sphericalharmonic.cxx @@ -1,15 +1,19 @@ #include "sphericalharmonic.hxx" -#include #include #include +#include #include + #include +namespace Multipole { +using namespace Loop; +using namespace std; static const CCTK_REAL PI = acos(-1.0); -static double factorial(int n) { +static inline double factorial(int n) { double returnval = 1; for (int i = n; i >= 1; i--) { returnval *= i; @@ -22,7 +26,6 @@ static inline double combination(int n, int m) { assert(n >= 0); assert(m >= 0); assert(m <= n); - return factorial(n) / (factorial(m) * factorial(n - m)); } @@ -30,11 +33,8 @@ static inline int imin(int a, int b) { return a < b ? a : b; } static inline int imax(int a, int b) { return a > b ? a : b; } -void Multipole_SphericalHarmonic(int s, int l, int m, CCTK_REAL th, - CCTK_REAL ph, CCTK_REAL *reY, CCTK_REAL *imY) { - // assert(s == -2 && l == 2 && m == 2); - // *reY = 1.0/2.0 * sqrt(5/PI) * pow(cos(th/2), 4) * cos(2*ph); - // *imY = 1.0/2.0 * sqrt(5/PI) * pow(cos(th/2), 4) * sin(2*ph); +void SphericalHarmonic(int s, int l, int m, CCTK_REAL th, CCTK_REAL ph, + CCTK_REAL *reY, CCTK_REAL *imY) { double all_coeff = 0, sum = 0; all_coeff = pow(-1.0, m); all_coeff *= sqrt(factorial(l + m) * factorial(l - m) * (2 * l + 1) / @@ -49,30 +49,18 @@ void Multipole_SphericalHarmonic(int s, int l, int m, CCTK_REAL th, *imY = all_coeff * sum * sin(m * ph); } -void Multipole_HarmonicSetup(int s, int l, int m, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], - CCTK_REAL reY[], CCTK_REAL imY[]) { +void HarmonicSetup(int s, int l, int m, int array_size, CCTK_REAL const th[], + CCTK_REAL const ph[], CCTK_REAL reY[], CCTK_REAL imY[]) { for (int i = 0; i < array_size; i++) { - Multipole_SphericalHarmonic(s, l, m, th[i], ph[i], &reY[i], &imY[i]); + SphericalHarmonic(s, l, m, th[i], ph[i], &reY[i], &imY[i]); } } // Fill a grid function with a given spherical harmonic extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) { - using namespace Loop; - using namespace std; - DECLARE_CCTK_ARGUMENTS_Multipole_SetHarmonic; + DECLARE_CCTK_ARGUMENTSX_Multipole_SetHarmonic; DECLARE_CCTK_PARAMETERS; - const int dim = 3; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout(cctkGH, indextype); - - const GF3D2 harmonic_re_(layout, harmonic_re); - const GF3D2 harmonic_im_(layout, harmonic_im); - - const GridDescBaseDevice grid(cctkGH); grid.loop_int<0, 0, 0>( grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { @@ -85,60 +73,20 @@ extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) { CCTK_REAL re = 0; CCTK_REAL im = 0; - Multipole_SphericalHarmonic(test_sw, test_l, test_m, theta, phi, &re, - &im); + SphericalHarmonic(test_sw, test_l, test_m, theta, phi, &re, &im); CCTK_REAL fac = test_mode_proportional_to_r ? vcoordr : 1.0; - harmonic_re_(p.I) = re * fac; - harmonic_im_(p.I) = im * fac; + harmonic_re(p.I) = re * fac; + harmonic_im(p.I) = im * fac; }); - // for (int k = 0; k < cctk_lsh[2]+1; k++) - // { - // for (int j = 0; j < cctk_lsh[1]+1; j++) - // { - // for (int i = 0; i < cctk_lsh[0]+1; i++) - // { - // int index = i + j * (cctk_lsh[0]+1) + k * - // (cctk_lsh[0]+1)*(cctk_lsh[1]+1) ; vcoordr = - // sqrt(vcoordx[index]*vcoordx[index] + vcoordy[index]*vcoordy[index] + - // vcoordz[index]*vcoordz[index]); - - // CCTK_REAL theta = acos(vcoordz[index]/vcoordr); - // if (vcoordr == 0) theta = 0; - // CCTK_REAL phi = atan2(vcoordy[index],vcoordx[index]); - - // CCTK_REAL re = 0; - // CCTK_REAL im = 0; - - // Multipole_SphericalHarmonic(test_sw,test_l,test_m,theta,phi, - // &re, &im); - - // CCTK_REAL fac = test_mode_proportional_to_r ? vcoordr : 1.0; - - // harmonic_re[index] = re * fac; - // harmonic_im[index] = im * fac; - // } - // } - // } return; } extern "C" void Multipole_SetHarmonicWeyl(CCTK_ARGUMENTS) { - using namespace Loop; - using namespace std; - DECLARE_CCTK_ARGUMENTS_Multipole_SetHarmonicWeyl; + DECLARE_CCTK_ARGUMENTSX_Multipole_SetHarmonicWeyl; DECLARE_CCTK_PARAMETERS; - const int dim = 3; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout(cctkGH, indextype); - - const GF3D2 Psi4r_(layout, Psi4r); - const GF3D2 Psi4i_(layout, Psi4i); - - const GridDescBaseDevice grid(cctkGH); grid.loop_int<0, 0, 0>( grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { @@ -151,22 +99,24 @@ extern "C" void Multipole_SetHarmonicWeyl(CCTK_ARGUMENTS) { CCTK_REAL re22 = 0; CCTK_REAL im22 = 0; - Multipole_SphericalHarmonic(test_sw, 2, 2, theta, phi, &re22, &im22); + SphericalHarmonic(test_sw, 2, 2, theta, phi, &re22, &im22); CCTK_REAL re2m2 = 0; CCTK_REAL im2m2 = 0; - Multipole_SphericalHarmonic(test_sw, 2, -2, theta, phi, &re2m2, &im2m2); + SphericalHarmonic(test_sw, 2, -2, theta, phi, &re2m2, &im2m2); CCTK_REAL re31 = 0; CCTK_REAL im31 = 0; - Multipole_SphericalHarmonic(test_sw, 3, 1, theta, phi, &re31, &im31); + SphericalHarmonic(test_sw, 3, 1, theta, phi, &re31, &im31); CCTK_REAL fac = test_mode_proportional_to_r ? vcoordr : 1.0; - Psi4r_(p.I) = (re22 + 0.5 * re2m2 + 0.25 * re31) * fac; - Psi4i_(p.I) = (im22 + 0.5 * im2m2 + 0.25 * im31) * fac; + Psi4r(p.I) = (re22 + 0.5 * re2m2 + 0.25 * re31) * fac; + Psi4i(p.I) = (im22 + 0.5 * im2m2 + 0.25 * im31) * fac; }); return; } + +} // namespace Multipole diff --git a/Multipole/src/sphericalharmonic.hxx b/Multipole/src/sphericalharmonic.hxx index 49e4cc99..695f7c3c 100644 --- a/Multipole/src/sphericalharmonic.hxx +++ b/Multipole/src/sphericalharmonic.hxx @@ -5,11 +5,14 @@ #include #include -void Multipole_HarmonicSetup(int s, int l, int m, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], - CCTK_REAL reY[], CCTK_REAL imY[]); +namespace Multipole { -void Multipole_SphericalHarmonic(int s, int l, int m, CCTK_REAL th, - CCTK_REAL ph, CCTK_REAL *reY, CCTK_REAL *imY); +void SphericalHarmonic(int s, int l, int m, CCTK_REAL th, CCTK_REAL ph, + CCTK_REAL *reY, CCTK_REAL *imY); + +void HarmonicSetup(int s, int l, int m, int array_size, CCTK_REAL const th[], + CCTK_REAL const ph[], CCTK_REAL reY[], CCTK_REAL imY[]); + +} // namespace Multipole #endif // #ifndef MULTIPOLE_SPHERICALHARMONIC_HXX diff --git a/Multipole/src/tests.cc b/Multipole/src/tests.cc index 20e24171..5e544f72 100644 --- a/Multipole/src/tests.cc +++ b/Multipole/src/tests.cc @@ -189,8 +189,8 @@ void Multipole_TestOrthonormality(CCTK_ARGUMENTS) { reY[sw][l][m + l] = new CCTK_REAL[array_size]; imY[sw][l][m + l] = new CCTK_REAL[array_size]; - Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[sw][l][m + l], - imY[sw][l][m + l]); + HarmonicSetup(sw, l, m, array_size, th, ph, reY[sw][l][m + l], + imY[sw][l][m + l]); } } }