Skip to content

Commit

Permalink
Multipole: add memember func Surface::integrate
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwJi
lwJi committed Jun 26, 2024
1 parent 6d77a9f commit 5721d3b
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 19 deletions.
70 changes: 56 additions & 14 deletions Multipole/src/surface.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,34 +9,30 @@

namespace Multipole {

void Surface::interp(CCTK_ARGUMENTS, int realIdx, int imagIdx) {
void Surface::interpolate(CCTK_ARGUMENTS, int realIdx, int imagIdx) {
DECLARE_CCTK_PARAMETERS;

const CCTK_INT nPoints =
CCTK_MyProc(cctkGH) == 0 ? (nTheta_ + 1) * (nPhi_ + 1) : 0;

const void *interpCoords[Loop::dim] = {(const void *)x_.data(),
(const void *)y_.data(),
(const void *)z_.data()};
const void *interpCoords[Loop::dim] = {x_.data(), y_.data(), z_.data()};

CCTK_INT nInputArrays = imagIdx == -1 ? 1 : 2;
CCTK_INT nOutputArrays = imagIdx == -1 ? 1 : 2;

const CCTK_INT inputArrayIndices[2] = {realIdx, imagIdx};

// Interpolation result
CCTK_POINTER outputArrays[2];
outputArrays[0] = real_.data();
outputArrays[1] = imag_.data();
CCTK_POINTER outputArrays[2] = {real_.data(), imag_.data()};

/* DriverInterpolate arguments that aren't currently used */
const int coordSystemHandle = 0;
CCTK_INT const interpCoords_type_code = 0;
CCTK_INT const outputArrayTypes[1] = {0};
const CCTK_INT interpCoordsTypeCode = 0;
const CCTK_INT outputArrayTypes[1] = {0};

int interpHandle = CCTK_InterpHandle("CarpetX");
if (interpHandle < 0) {
CCTK_VERROR("Could not obtain inteprolator handle for built-in 'CarpetX' "
CCTK_VERROR("Could not obtain interpolator handle for built-in 'CarpetX' "
"interpolator: %d",
interpHandle);
}
Expand All @@ -47,7 +43,7 @@ void Surface::interp(CCTK_ARGUMENTS, int realIdx, int imagIdx) {
int ierr = Util_TableSetFromString(paramTableHandle, interpolator_pars);

ierr = DriverInterpolate(cctkGH, Loop::dim, interpHandle, paramTableHandle,
coordSystemHandle, nPoints, interpCoords_type_code,
coordSystemHandle, nPoints, interpCoordsTypeCode,
interpCoords, nInputArrays, inputArrayIndices,
nOutputArrays, outputArrayTypes, outputArrays);

Expand All @@ -57,12 +53,58 @@ void Surface::interp(CCTK_ARGUMENTS, int realIdx, int imagIdx) {
}

if (imagIdx == -1) {
for (int i = 0; i < (nTheta_ + 1) * (nPhi_ + 1); i++) {
imag_[i] = 0;
}
std::fill(imag_.begin(), imag_.end(), 0);
}

Util_TableDestroy(paramTableHandle);
}

// Take the integral of conj(array1)*array2*sin(th)
void Surface::integrate(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 *outRe,
CCTK_REAL *outIm) {
DECLARE_CCTK_PARAMETERS;

std::vector<CCTK_REAL> fReal(theta_.size());
std::vector<CCTK_REAL> fImag(theta_.size());

// integrand: conj(array1)*array2*sin(th)
for (size_t i = 0; i < theta_.size(); ++i) {
fReal[i] = (array1r[i] * array2r[i] + array1i[i] * array2i[i]) *
std::sin(theta_[i]);
fImag[i] = (array1r[i] * array2i[i] - array1i[i] * array2r[i]) *
std::sin(theta_[i]);
}

if (CCTK_Equals(integration_method, "midpoint")) {
*outRe = Midpoint2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
*outIm = Midpoint2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
} else if (CCTK_Equals(integration_method, "trapezoidal")) {
*outRe =
Trapezoidal2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
*outIm =
Trapezoidal2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
} else if (CCTK_Equals(integration_method, "Simpson")) {
if (nPhi_ % 2 != 0 || nTheta_ % 2 != 0) {
CCTK_WARN(CCTK_WARN_ABORT, "The Simpson integration method requires even "
"nTheta_ and even nPhi_");
}
*outRe = Simpson2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
*outIm = Simpson2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
} else if (CCTK_Equals(integration_method, "DriscollHealy")) {
if (nTheta_ % 2 != 0) {
CCTK_WARN(CCTK_WARN_ABORT,
"The Driscoll&Healy integration method requires even nTheta_");
}
*outRe =
DriscollHealy2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
*outIm =
DriscollHealy2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
} else {
CCTK_WARN(CCTK_WARN_ABORT, "internal error");
}
}

} // namespace Multipole
18 changes: 13 additions & 5 deletions Multipole/src/surface.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#ifndef MULTIPOLE_SURFACE_HXX
#define MULTIPOLE_SURFACE_HXX

#include "integrate.hxx"
#include "sphericalharmonic.hxx"

#include <cctk.h>
Expand All @@ -29,14 +30,14 @@ public:

// Add an offset for midpoint integration.
const CCTK_REAL is_midpoint = static_cast<CCTK_REAL>(isMidPoint);
const CCTK_REAL dTheta = PI / (nTheta + is_midpoint);
const CCTK_REAL dPhi = 2 * PI / (nPhi + is_midpoint);
dTheta_ = PI / (nTheta + is_midpoint);
dPhi_ = 2 * PI / (nPhi + is_midpoint);

for (int it = 0; it <= nTheta; ++it) {
for (int ip = 0; ip <= nPhi; ++ip) {
const int i = index2D(it, ip);
theta_[i] = it * dTheta + 0.5 * dTheta * is_midpoint;
phi_[i] = ip * dPhi + 0.5 * dPhi * is_midpoint;
theta_[i] = it * dTheta_ + 0.5 * dTheta_ * is_midpoint;
phi_[i] = ip * dPhi_ + 0.5 * dPhi_ * is_midpoint;
}
}
}
Expand All @@ -46,12 +47,19 @@ public:
const std::vector<CCTK_REAL> &getTheta() const { return theta_; }
const std::vector<CCTK_REAL> &getPhi() const { return phi_; }

void interp(CCTK_ARGUMENTS, int real_idx, int imag_idx);
void interpolate(CCTK_ARGUMENTS, int real_idx, int imag_idx);

void integrate(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 *outre,
CCTK_REAL *outim);

protected:
inline int index2D(int it, int ip) const { return it + (nTheta_ + 1) * ip; }

const int nTheta_, nPhi_;
CCTK_REAL dTheta_, dPhi_;

std::vector<CCTK_REAL> theta_, phi_;
std::vector<CCTK_REAL> x_, y_, z_; // embedding map
Expand Down

0 comments on commit 5721d3b

Please sign in to comment.