Skip to content

Commit

Permalink
Multipole: use vector for real and imag
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwJi
lwJi committed Jun 25, 2024
1 parent 50e238e commit af52cff
Show file tree
Hide file tree
Showing 8 changed files with 74 additions and 62 deletions.
8 changes: 5 additions & 3 deletions Multipole/src/integrate.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,9 +171,11 @@ CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f, int const nx,
return hx * hy * integrand_sum;
}

void 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[],
void Integrate(int array_size, int nthetap,
const std::vector<CCTK_REAL> array1r,
const std::vector<CCTK_REAL> array1i,
const std::vector<CCTK_REAL> array2r,
const std::vector<CCTK_REAL> array2i, CCTK_REAL const th[],
CCTK_REAL const ph[], CCTK_REAL *outre, CCTK_REAL *outim) {
DECLARE_CCTK_PARAMETERS;

Expand Down
13 changes: 8 additions & 5 deletions Multipole/src/integrate.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>

#include <vector>

namespace Multipole {

CCTK_REAL Midpoint2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx,
Expand All @@ -19,11 +21,12 @@ CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx,
CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *f, int nx, int ny,
CCTK_REAL hx, CCTK_REAL hy);

void Integrate(int array_size, int ntheta, CCTK_REAL const array1r[],
CCTK_REAL const array1i[], CCTK_REAL const array2r[],
CCTK_REAL const array2i[], CCTK_REAL const th[],
CCTK_REAL const pph[], CCTK_REAL out_arrayr[],
CCTK_REAL out_arrayi[]);
void Integrate(int array_size, int nthetap,
const std::vector<CCTK_REAL> array1r,
const std::vector<CCTK_REAL> array1i,
const std::vector<CCTK_REAL> array2r,
const std::vector<CCTK_REAL> array2i, CCTK_REAL const th[],
CCTK_REAL const ph[], CCTK_REAL *outre, CCTK_REAL *outim);

} // namespace Multipole

Expand Down
8 changes: 4 additions & 4 deletions Multipole/src/interpolate.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ static void report_interp_error(int ierr) {
}

void Interp(CCTK_ARGUMENTS, CCTK_REAL xs[], CCTK_REAL ys[], CCTK_REAL zs[],
int real_idx, int imag_idx, CCTK_REAL sphere_real[],
CCTK_REAL sphere_imag[]) {
int real_idx, int imag_idx, std::vector<CCTK_REAL> &sphere_real,
std::vector<CCTK_REAL> &sphere_imag) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;

Expand All @@ -29,8 +29,8 @@ void Interp(CCTK_ARGUMENTS, CCTK_REAL xs[], CCTK_REAL ys[], CCTK_REAL zs[],
const CCTK_INT input_array_indices[2] = {real_idx, imag_idx};
// Interpolation result
CCTK_POINTER output_arrays[2];
output_arrays[0] = sphere_real;
output_arrays[1] = sphere_imag;
output_arrays[0] = sphere_real.data();
output_arrays[1] = sphere_imag.data();

int ierr = -1;
/* DriverInterpolate arguments that aren't currently used */
Expand Down
5 changes: 4 additions & 1 deletion Multipole/src/interpolate.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,15 @@
#include <cctk_Functions.h>
#include <util_Table.h>

#include <vector>

namespace Multipole {

// This function interpolates psi4 onto the sphere in cartesian coordinates as
// created by CoordSetup.
void Interp(CCTK_ARGUMENTS, CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[],
int real_idx, int imag_idx, CCTK_REAL psi4r[], CCTK_REAL psi4i[]);
int real_idx, int imag_idx, std::vector<CCTK_REAL> &psi4r,
std::vector<CCTK_REAL> &psi4i);

} // namespace Multipole

Expand Down
67 changes: 36 additions & 31 deletions Multipole/src/multipole.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,11 @@
namespace Multipole {
using namespace std;

static vector<CCTK_REAL> real, imag;

// Y[sw][l][m][i]
static vector<vector<vector<vector<CCTK_REAL> > > > reY, imY;

static const int max_vars = 10;

// since each variable can have at most one spin weight, max_vars is a good
Expand All @@ -31,9 +36,10 @@ static void fillVariable(int idx, const char *optString, void *callbackArg) {

VariableParseArray *vs = static_cast<VariableParseArray *>(callbackArg);

assert(vs->numVars < max_vars); // Ensure we don't exceed max_vars
VariableParse *v = &vs->vars[vs->numVars++]; // Increment numVars and get
// reference to next VariableParse
assert(vs->numVars < max_vars); // Ensure we don't exceed max_vars
VariableParse *v =
&vs->vars[vs->numVars++]; // Increment numVars and get
// reference to next VariableParse

v->index = idx;

Expand Down Expand Up @@ -96,31 +102,27 @@ static void get_spin_weights(VariableParse vars[], int n_vars,
*n_weights = n_spin_weights;
}

static void
setup_harmonics(const int spin_weights[max_spin_weights], int n_spin_weights,
int lmax, CCTK_REAL th[], CCTK_REAL ph[], int array_size,
CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes],
CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes]) {
for (int si = 0; si < max_spin_weights; si++) {
for (int l = 0; l < max_l_modes; l++) {
for (int mi = 0; mi < max_m_modes; mi++) {
reY[si][l][mi] = 0;
imY[si][l][mi] = 0;
}
}
}
static void setup_harmonics(const int spin_weights[max_spin_weights],
int n_spin_weights, int lmax, CCTK_REAL th[],
CCTK_REAL ph[], int array_size,
vector<vector<vector<vector<CCTK_REAL> > > > &reY,
vector<vector<vector<vector<CCTK_REAL> > > > &imY) {
for (int si = 0; si < n_spin_weights; si++) {
int sw = spin_weights[si];

vector<vector<vector<CCTK_REAL> > > reY_s, imY_s;
for (int l = 0; l <= lmax; l++) {
vector<vector<CCTK_REAL> > reY_s_l, imY_s_l;
for (int m = -l; m <= l; m++) {
reY[si][l][m + l] = new CCTK_REAL[array_size];
imY[si][l][m + l] = new CCTK_REAL[array_size];

HarmonicSetup(sw, l, m, array_size, th, ph, reY[si][l][m + l],
imY[si][l][m + l]);
vector<CCTK_REAL> reY_s_l_m(array_size), imY_s_l_m(array_size);
HarmonicSetup(sw, l, m, array_size, th, ph, reY_s_l_m, imY_s_l_m);
reY_s_l.push_back(reY_s_l_m);
imY_s_l.push_back(imY_s_l_m);
}
reY_s.push_back(reY_s_l);
imY_s.push_back(imY_s_l);
}
reY.push_back(reY_s);
imY.push_back(imY_s);
}
}

Expand Down Expand Up @@ -212,9 +214,9 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) {
static CCTK_REAL *xs, *ys, *zs;
static CCTK_REAL *xhat, *yhat, *zhat;
static CCTK_REAL *th, *ph;
static CCTK_REAL *real = 0, *imag = 0;
static CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes];
static CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes];
// static CCTK_REAL *real = 0, *imag = 0;
// static CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes];
// static CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes];
static VariableParse vars[max_vars];
static int n_variables = 0;
static int spin_weights[max_spin_weights];
Expand All @@ -232,8 +234,11 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) {
assert(lmax + 1 <= max_l_modes);

if (!initialized) {
real = new CCTK_REAL[array_size];
imag = new CCTK_REAL[array_size];
real.resize(array_size);
imag.resize(array_size);

// real = new CCTK_REAL[array_size];
// imag = new CCTK_REAL[array_size];
th = new CCTK_REAL[array_size];
ph = new CCTK_REAL[array_size];
xs = new CCTK_REAL[array_size];
Expand Down Expand Up @@ -265,8 +270,8 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) {
ScaleCartesian(ntheta, nphi, radius[i], xhat, yhat, zhat, xs, ys, zs);

// Interpolate Psi4r and Psi4i
Interp(CCTK_PASS_CTOC, xs, ys, zs, vars[v].index, vars[v].imagIndex,
real, imag);
Interp(CCTK_PASS_CTOC, xs, ys, zs, vars[v].index, vars[v].imagIndex, real,
imag);

for (int l = 0; l <= lmax; l++) {
for (int m = -l; m <= l; m++) {
Expand All @@ -277,8 +282,8 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) {

} // loop over m
} // loop over l
output_1d(CCTK_PASS_CTOC, &vars[v], radius[i], th, ph, real, imag,
array_size);
output_1d(CCTK_PASS_CTOC, &vars[v], radius[i], th, ph, real.data(),
imag.data(), array_size);
} // loop over radii
} // loop over variables

Expand Down
4 changes: 2 additions & 2 deletions Multipole/src/sphericalharmonic.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
#include <assert.h>
#include <iostream>
#include <math.h>
#include <vector>

#include <loop_device.hxx>

Expand Down Expand Up @@ -50,7 +49,8 @@ void SphericalHarmonic(int s, int l, int m, CCTK_REAL th, CCTK_REAL ph,
}

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[]) {
CCTK_REAL const ph[], std::vector<CCTK_REAL> &reY,
std::vector<CCTK_REAL> &imY) {
for (int i = 0; i < array_size; i++) {
SphericalHarmonic(s, l, m, th[i], ph[i], &reY[i], &imY[i]);
}
Expand Down
5 changes: 4 additions & 1 deletion Multipole/src/sphericalharmonic.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,16 @@
#include <cctk_Parameters.h>
#include <cctk_Arguments.h>

#include <vector>

namespace Multipole {

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[]);
CCTK_REAL const ph[], std::vector<CCTK_REAL> &reY,
std::vector<CCTK_REAL> &imY);

} // namespace Multipole

Expand Down
26 changes: 11 additions & 15 deletions Multipole/src/tests.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,19 +143,23 @@ extern "C" void Multipole_TestOrthonormality(CCTK_ARGUMENTS) {
CoordSetup(xhat, yhat, zhat, th, ph);

/* Populate spherical-harmonic array */
CCTK_REAL *reY[1][max_l_modes][max_m_modes];
CCTK_REAL *imY[1][max_l_modes][max_m_modes];
vector<vector<vector<vector<CCTK_REAL> > > > reY, imY;

for (int sw = 0; sw <= 0; sw++) {
vector<vector<vector<CCTK_REAL> > > reY_s, imY_s;
for (int l = 0; l < max_l_modes; l++) {
vector<vector<CCTK_REAL> > reY_s_l, imY_s_l;
for (int m = -l; m <= l; m++) {
reY[sw][l][m + l] = new CCTK_REAL[array_size];
imY[sw][l][m + l] = new CCTK_REAL[array_size];

HarmonicSetup(sw, l, m, array_size, th, ph, reY[sw][l][m + l],
imY[sw][l][m + l]);
vector<CCTK_REAL> reY_s_l_m(array_size), imY_s_l_m(array_size);
HarmonicSetup(sw, l, m, array_size, th, ph, reY_s_l_m, imY_s_l_m);
reY_s_l.push_back(reY_s_l_m);
imY_s_l.push_back(imY_s_l_m);
}
reY_s.push_back(reY_s_l);
imY_s.push_back(imY_s_l);
}
reY.push_back(reY_s);
imY.push_back(imY_s);
}

/* Loop over l and m, assign Ylm to (rel,imag), and compute the scalar
Expand Down Expand Up @@ -192,14 +196,6 @@ extern "C" void Multipole_TestOrthonormality(CCTK_ARGUMENTS) {
}
assert(idx == 1 * N * (N + 1) / 2);

for (int sw = 0; sw <= 0; sw++) {
for (int l = 0; l < max_l_modes; l++) {
for (int m = -l; m <= l; m++) {
delete[] imY[0][l][m + l];
delete[] reY[0][l][m + l];
}
}
}
delete[] zhat;
delete[] yhat;
delete[] xhat;
Expand Down

0 comments on commit af52cff

Please sign in to comment.