Skip to content

Commit

Permalink
Multipole: decorate sphericalharmonic
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwJi
lwJi committed Jun 23, 2024
1 parent ab2e3fe commit b49f04e
Show file tree
Hide file tree
Showing 5 changed files with 37 additions and 86 deletions.
6 changes: 2 additions & 4 deletions Multipole/schedule.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand Down
4 changes: 2 additions & 2 deletions Multipole/src/multipole.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]);
}
}
}
Expand Down
96 changes: 23 additions & 73 deletions Multipole/src/sphericalharmonic.cxx
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,19 @@
#include "sphericalharmonic.hxx"

#include <math.h>
#include <assert.h>
#include <iostream>
#include <math.h>
#include <vector>

#include <loop_device.hxx>

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;
Expand All @@ -22,19 +26,15 @@ 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));
}

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) /
Expand All @@ -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<int, dim> indextype = {0, 0, 0};
const GF3D2layout layout(cctkGH, indextype);

const GF3D2<CCTK_REAL> harmonic_re_(layout, harmonic_re);
const GF3D2<CCTK_REAL> 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 {
Expand All @@ -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<int, dim> indextype = {0, 0, 0};
const GF3D2layout layout(cctkGH, indextype);

const GF3D2<CCTK_REAL> Psi4r_(layout, Psi4r);
const GF3D2<CCTK_REAL> Psi4i_(layout, Psi4i);

const GridDescBaseDevice grid(cctkGH);
grid.loop_int<0, 0, 0>(
grid.nghostzones,
[=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
Expand All @@ -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
13 changes: 8 additions & 5 deletions Multipole/src/sphericalharmonic.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,14 @@
#include <cctk_Parameters.h>
#include <cctk_Arguments.h>

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
4 changes: 2 additions & 2 deletions Multipole/src/tests.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]);
}
}
}
Expand Down

0 comments on commit b49f04e

Please sign in to comment.