From 488d08f721aff9cd92ba4f297d54fcf1d9588135 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Wed, 19 Jun 2024 18:54:32 +0000 Subject: [PATCH] Multipole: clang-format files --- Multipole/src/integrate.cc | 126 +++++----- Multipole/src/integrate.hh | 11 +- Multipole/src/interp_exp.cc | 103 ++++---- Multipole/src/interpolate.cc | 74 +++--- Multipole/src/interpolate.hh | 11 +- Multipole/src/multipole.cc | 296 +++++++++++------------ Multipole/src/multipole.hh | 66 +++--- Multipole/src/sphericalharmonic.cc | 172 +++++++------- Multipole/src/sphericalharmonic.hh | 16 +- Multipole/src/tests.cc | 183 +++++++-------- Multipole/src/utils.cc | 361 +++++++++++++---------------- Multipole/src/utils.hh | 55 ++--- 12 files changed, 678 insertions(+), 796 deletions(-) diff --git a/Multipole/src/integrate.cc b/Multipole/src/integrate.cc index 6ef39ceb..98a30ac3 100644 --- a/Multipole/src/integrate.cc +++ b/Multipole/src/integrate.cc @@ -19,7 +19,7 @@ interval [a,b] into nx subintervals of spacing h = (b-a)/nx. These have coordinates [x_i-1, xi] where x_i = x_0 + i h. So i runs from 0 to nx. We require the function to integrate at the points F[x_i, y_i]. We have x_0 = a and x_n = b. Check: x_n = x_0 + n (b-a)/n = a -+ b - a = b. Good. ++ b - a = b. Good. If we are given these points in an array, we also need the width and height of the array. To get an actual integral, we also need the grid @@ -28,57 +28,64 @@ the integral. */ - - -#define idx(xx,yy) (assert((xx) <= nx), assert((xx) >= 0), assert((yy) <= ny), assert((yy) >= 0), ((xx) + (yy) * (nx+1))) +#define idx(xx, yy) \ + (assert((xx) <= nx), assert((xx) >= 0), assert((yy) <= ny), \ + assert((yy) >= 0), ((xx) + (yy) * (nx + 1))) // Hard coded 2D integrals -CCTK_REAL Midpoint2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy) -{ +CCTK_REAL Midpoint2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, + CCTK_REAL hy) { CCTK_REAL integrand_sum = 0.0; int ix = 0, iy = 0; - assert(nx > 0); assert(ny > 0); assert (f); + assert(nx > 0); + assert(ny > 0); + assert(f); for (iy = 0; iy <= ny; iy++) for (ix = 0; ix <= nx; ix++) - integrand_sum += f[idx(ix,iy)]; + integrand_sum += f[idx(ix, iy)]; return hx * hy * integrand_sum; } -CCTK_REAL Trapezoidal2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy) -{ +CCTK_REAL Trapezoidal2DIntegral(CCTK_REAL const *f, int nx, int ny, + CCTK_REAL hx, CCTK_REAL hy) { CCTK_REAL integrand_sum = 0.0; int ix = 0, iy = 0; - assert(nx > 0); assert(ny > 0); assert (f); + assert(nx > 0); + assert(ny > 0); + assert(f); // Corners - integrand_sum += f[idx(0,0)] + f[idx(nx,0)] + f[idx(0,ny)] + f[idx(nx,ny)]; + integrand_sum += + f[idx(0, 0)] + f[idx(nx, 0)] + f[idx(0, ny)] + f[idx(nx, ny)]; // Edges - for (ix = 1; ix <= nx-1; ix++) - integrand_sum += 2 * f[idx(ix,0)] + 2 * f[idx(ix,ny)]; + for (ix = 1; ix <= nx - 1; ix++) + integrand_sum += 2 * f[idx(ix, 0)] + 2 * f[idx(ix, ny)]; - for (iy = 1; iy <= ny-1; iy++) - integrand_sum += 2 * f[idx(0,iy)] + 2 * f[idx(nx,iy)]; + for (iy = 1; iy <= ny - 1; iy++) + integrand_sum += 2 * f[idx(0, iy)] + 2 * f[idx(nx, iy)]; // Interior - for (iy = 1; iy <= ny-1; iy++) - for (ix = 1; ix <= nx-1; ix++) - integrand_sum += 4 * f[idx(ix,iy)]; + for (iy = 1; iy <= ny - 1; iy++) + for (ix = 1; ix <= nx - 1; ix++) + integrand_sum += 4 * f[idx(ix, iy)]; - return (double) 1 / (double) 4 * hx * hy * integrand_sum; + return (double)1 / (double)4 * hx * hy * integrand_sum; } -CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy) -{ +CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, + CCTK_REAL hy) { CCTK_REAL integrand_sum = 0; int ix = 0, iy = 0; - assert(nx > 0); assert(ny > 0); assert (f); + assert(nx > 0); + assert(ny > 0); + assert(f); assert(nx % 2 == 0); assert(ny % 2 == 0); @@ -86,90 +93,86 @@ CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CC int py = ny / 2; // Corners - integrand_sum += f[idx(0,0)] + f[idx(nx,0)] + f[idx(0,ny)] + f[idx(nx,ny)]; + integrand_sum += + f[idx(0, 0)] + f[idx(nx, 0)] + f[idx(0, ny)] + f[idx(nx, ny)]; // Edges for (iy = 1; iy <= py; iy++) - integrand_sum += 4 * f[idx(0,2*iy-1)] + 4 * f[idx(nx,2*iy-1)]; + integrand_sum += 4 * f[idx(0, 2 * iy - 1)] + 4 * f[idx(nx, 2 * iy - 1)]; - for (iy = 1; iy <= py-1; iy++) - integrand_sum += 2 * f[idx(0,2*iy)] + 2 * f[idx(nx,2*iy)]; + for (iy = 1; iy <= py - 1; iy++) + integrand_sum += 2 * f[idx(0, 2 * iy)] + 2 * f[idx(nx, 2 * iy)]; for (ix = 1; ix <= px; ix++) - integrand_sum += 4 * f[idx(2*ix-1,0)] + 4 * f[idx(2*ix-1,ny)]; + integrand_sum += 4 * f[idx(2 * ix - 1, 0)] + 4 * f[idx(2 * ix - 1, ny)]; - for (ix = 1; ix <= px-1; ix++) - integrand_sum += 2 * f[idx(2*ix,0)] + 2 * f[idx(2*ix,ny)]; + for (ix = 1; ix <= px - 1; ix++) + integrand_sum += 2 * f[idx(2 * ix, 0)] + 2 * f[idx(2 * ix, ny)]; // Interior for (iy = 1; iy <= py; iy++) for (ix = 1; ix <= px; ix++) - integrand_sum += 16 * f[idx(2*ix-1,2*iy-1)]; + integrand_sum += 16 * f[idx(2 * ix - 1, 2 * iy - 1)]; - for (iy = 1; iy <= py-1; iy++) + for (iy = 1; iy <= py - 1; iy++) for (ix = 1; ix <= px; ix++) - integrand_sum += 8 * f[idx(2*ix-1,2*iy)]; + integrand_sum += 8 * f[idx(2 * ix - 1, 2 * iy)]; for (iy = 1; iy <= py; iy++) - for (ix = 1; ix <= px-1; ix++) - integrand_sum += 8 * f[idx(2*ix,2*iy-1)]; + for (ix = 1; ix <= px - 1; ix++) + integrand_sum += 8 * f[idx(2 * ix, 2 * iy - 1)]; - for (iy = 1; iy <= py-1; iy++) - for (ix = 1; ix <= px-1; ix++) - integrand_sum += 4 * f[idx(2*ix,2*iy)]; + for (iy = 1; iy <= py - 1; iy++) + for (ix = 1; ix <= px - 1; ix++) + integrand_sum += 4 * f[idx(2 * ix, 2 * iy)]; - return ((double) 1 / (double) 9) * hx * hy * integrand_sum; + return ((double)1 / (double)9) * hx * hy * integrand_sum; } // See: J.R. Driscoll and D.M. Healy Jr., Computing Fourier transforms // and convolutions on the 2-sphere. Advances in Applied Mathematics, // 15(2):202–250, 1994. -CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f, - int const nx, int const ny, - CCTK_REAL const hx, CCTK_REAL const hy) -{ +CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f, int const nx, + int const ny, CCTK_REAL const hx, + CCTK_REAL const hy) { assert(f); assert(nx >= 0); assert(ny >= 0); assert(nx % 2 == 0); - + CCTK_REAL integrand_sum = 0.0; - + // Skip the poles (ix=0 and ix=nx), since the weight there is zero // anyway -#pragma omp parallel for reduction(+: integrand_sum) - for (int ix = 1; ix < nx; ++ ix) - { - +#pragma omp parallel for reduction(+ : integrand_sum) + for (int ix = 1; ix < nx; ++ix) { + // These weights lead to an almost spectral convergence CCTK_REAL const theta = M_PI * ix / nx; CCTK_REAL weight = 0.0; - for (int l = 0; l < nx/2; ++ l) - { - weight += sin((2*l+1)*theta)/(2*l+1); + for (int l = 0; l < nx / 2; ++l) { + weight += sin((2 * l + 1) * theta) / (2 * l + 1); } weight *= 4.0 / M_PI; - + CCTK_REAL local_sum = 0.0; // Skip the last point (iy=ny), since we assume periodicity and // therefor it has the same value as the first point. We don't use // weights in this direction, which leads to spectral convergence. // (Yay periodicity!) - for (int iy = 0; iy < ny; ++ iy) - { - local_sum += f[idx(ix,iy)]; + for (int iy = 0; iy < ny; ++iy) { + local_sum += f[idx(ix, iy)]; } - + integrand_sum += weight * local_sum; - } - + return hx * hy * integrand_sum; } //// 1D integrals // -//static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h) +// static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h) //{ // CCTK_REAL integrand_sum = 0; // int i = 0; @@ -192,7 +195,8 @@ CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f, // //// 2D integral built up from 1D // -//static CCTK_REAL Composite2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy) +// static CCTK_REAL Composite2DIntegral(CCTK_REAL const *f, int nx, int ny, +// CCTK_REAL hx, CCTK_REAL hy) //{ // CCTK_REAL integrand_sum = 0; // diff --git a/Multipole/src/integrate.hh b/Multipole/src/integrate.hh index 85e2abba..d8d3e6db 100644 --- a/Multipole/src/integrate.hh +++ b/Multipole/src/integrate.hh @@ -3,13 +3,14 @@ #define __integrate_h #include "cctk.h" -CCTK_REAL Midpoint2DIntegral(CCTK_REAL const *f, int nx, int ny, - CCTK_REAL hx, CCTK_REAL hy); +CCTK_REAL Midpoint2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, + CCTK_REAL hy); -CCTK_REAL Trapezoidal2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy); +CCTK_REAL Trapezoidal2DIntegral(CCTK_REAL const *f, int nx, int ny, + CCTK_REAL hx, CCTK_REAL hy); -CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, - CCTK_REAL hx, CCTK_REAL hy); +CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, + CCTK_REAL hy); CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy); diff --git a/Multipole/src/interp_exp.cc b/Multipole/src/interp_exp.cc index cb68b0b9..faf23127 100644 --- a/Multipole/src/interp_exp.cc +++ b/Multipole/src/interp_exp.cc @@ -2,74 +2,66 @@ #include "interpolate.hh" -static -void report_interp_error(int ierr) -{ - if (ierr < 0) - { +static void report_interp_error(int ierr) { + if (ierr < 0) { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "CCTK_InterpGridArrays returned error code %d",ierr); + "CCTK_InterpGridArrays returned error code %d", ierr); } } - -void Multipole_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[]) -{ + +void Multipole_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[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; // Need parameters for the following: - // ntheta (dtheta = pi/(ntheta) + // ntheta (dtheta = pi/(ntheta) // nphi (dphi = 2pi/(nphi) // r (radius of sphere) // NOTE: depending on the interval of integration, denominator above may // need to be modified to avoid double counting - - CCTK_INT num_input_arrays = imag_idx == -1 ? 1 : 2; + CCTK_INT num_input_arrays = imag_idx == -1 ? 1 : 2; CCTK_INT num_output_arrays = imag_idx == -1 ? 1 : 2; const CCTK_INT num_dims = 3; - - const CCTK_INT num_points = CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0; // Only the 0 processor needs the points + + const CCTK_INT num_points = CCTK_MyProc(cctkGH) == 0 + ? (ntheta + 1) * (nphi + 1) + : 0; // Only the 0 processor needs the points int ierr = -1; - const void* interp_coords[num_dims] - = { (const void *) xs, - (const void *) ys, - (const void *) zs }; + const void *interp_coords[num_dims] = {(const void *)xs, (const void *)ys, + (const void *)zs}; - const CCTK_INT input_array_indices[2] - = { real_idx, - imag_idx }; + const CCTK_INT input_array_indices[2] = {real_idx, imag_idx}; // const CCTK_INT output_array_types[2] // = { CCTK_VARIABLE_REAL, // CCTK_VARIABLE_REAL }; // Interpolation result types: Not used by CarpetX DriverInterp CCTK_INT const output_array_types[1] = {0}; - + // void * output_arrays[2] // = { (void *) sphere_real, // (void *) sphere_imag }; - // Interpolation result - CCTK_POINTER output_arrays[2]; - output_arrays[0] = sphere_real; - output_arrays[1] = sphere_imag; + // Interpolation result + CCTK_POINTER output_arrays[2]; + output_arrays[0] = sphere_real; + output_arrays[1] = sphere_imag; - const int operator_handle = 0; //not used by CarpetX Interpolator + const int operator_handle = 0; // not used by CarpetX Interpolator // Interpolation parameter table CCTK_INT operations[1][num_dims]; - for (int var = 0 ; var < num_points; var++) { - operations[0][var] = 0; + for (int var = 0; var < num_points; var++) { + operations[0][var] = 0; } - - int operands[1][num_dims]; - for (int var = 0 ; var < num_points; var++) { - operands[0][var] = var; + + int operands[1][num_dims]; + for (int var = 0; var < num_points; var++) { + operands[0][var] = var; } int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); @@ -77,11 +69,13 @@ void Multipole_Interp(CCTK_ARGUMENTS, CCTK_VERROR("Can't create parameter table: %d", param_table_handle); if ((ierr = Util_TableSetInt(param_table_handle, 1, "order")) < 0) CCTK_VERROR("Can't set order in parameter table: %d", ierr); - if ((ierr = Util_TableSetIntArray(param_table_handle, num_points, (int const*const)operands, - "operand_indices")) < 0) + if ((ierr = Util_TableSetIntArray(param_table_handle, num_points, + (int const *const)operands, + "operand_indices")) < 0) CCTK_VERROR("Can't set operand_indices array in parameter table: %d", ierr); - if ((ierr = Util_TableSetIntArray(param_table_handle, num_points, (int const*const)operations, - "operation_codes")) < 0) + if ((ierr = Util_TableSetIntArray(param_table_handle, num_points, + (int const *const)operations, + "operation_codes")) < 0) CCTK_VERROR("Can't set operation_codes array in parameter table: %d", ierr); const int coord_system_handle = 0; // not used by CarpetX Interpolator @@ -92,27 +86,24 @@ void Multipole_Interp(CCTK_ARGUMENTS, // operator_handle, // param_table_handle, // coord_system_handle, - // CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0, // Only the 0 processor needs the points - // CCTK_VARIABLE_REAL, - // interp_coords, + // CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0, // Only the 0 + // processor needs the points CCTK_VARIABLE_REAL, interp_coords, // num_input_arrays, // input_array_indices, // num_output_arrays, // output_array_types, // output_arrays); - ierr = DriverInterpolate( - cctkGH, num_dims, operator_handle, param_table_handle, coord_system_handle, - num_points, CCTK_VARIABLE_REAL, interp_coords, num_input_arrays, (int const * const)input_array_indices, - num_output_arrays, output_array_types, output_arrays); - + ierr = DriverInterpolate( + cctkGH, num_dims, operator_handle, param_table_handle, + coord_system_handle, num_points, CCTK_VARIABLE_REAL, interp_coords, + num_input_arrays, (int const *const)input_array_indices, + num_output_arrays, output_array_types, output_arrays); report_interp_error(ierr); - if (imag_idx == -1) - { - for (int i = 0; i < (ntheta+1)*(nphi+1); i++) - { + if (imag_idx == -1) { + for (int i = 0; i < (ntheta + 1) * (nphi + 1); i++) { sphere_imag[i] = 0; } } @@ -122,14 +113,14 @@ void Multipole_Interp(CCTK_ARGUMENTS, // // Debugging routine // void Multipole_InterpVar(CCTK_ARGUMENTS, -// CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name, -// CCTK_REAL sphere_var[]) +// CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const +// char *var_name, CCTK_REAL sphere_var[]) // { // DECLARE_CCTK_ARGUMENTS; // DECLARE_CCTK_PARAMETERS; // // Need parameters for the following: -// // ntheta (dtheta = pi/(ntheta) +// // ntheta (dtheta = pi/(ntheta) // // nphi (dphi = 2pi/(nphi) // // r (radius of sphere) // // NOTE: depending on the interval of integration, denominator above may @@ -140,7 +131,7 @@ void Multipole_Interp(CCTK_ARGUMENTS, // const CCTK_INT num_dims = 3; // int ierr = -1; -// const void* interp_coords[num_dims] +// const void* interp_coords[num_dims] // = { (const void *) x, // (const void *) y, // (const void *) z }; diff --git a/Multipole/src/interpolate.cc b/Multipole/src/interpolate.cc index 452cb000..c495470c 100644 --- a/Multipole/src/interpolate.cc +++ b/Multipole/src/interpolate.cc @@ -3,53 +3,40 @@ #include "interpolate.hh" #include -static -void report_interp_error(int ierr) -{ - if (ierr < 0) - { +static void report_interp_error(int ierr) { + if (ierr < 0) { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "CCTK_InterpGridArrays returned error code %d",ierr); + "CCTK_InterpGridArrays returned error code %d", ierr); } } - -void Multipole_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[]) -{ + +void Multipole_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[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; // Need parameters for the following: - // ntheta (dtheta = pi/(ntheta) + // ntheta (dtheta = pi/(ntheta) // nphi (dphi = 2pi/(nphi) // r (radius of sphere) // NOTE: depending on the interval of integration, denominator above may // need to be modified to avoid double counting - - CCTK_INT num_input_arrays = imag_idx == -1 ? 1 : 2; + CCTK_INT num_input_arrays = imag_idx == -1 ? 1 : 2; CCTK_INT num_output_arrays = imag_idx == -1 ? 1 : 2; const CCTK_INT num_dims = 3; int ierr = -1; - const void* interp_coords[num_dims] - = { (const void *) xs, - (const void *) ys, - (const void *) zs }; + const void *interp_coords[num_dims] = {(const void *)xs, (const void *)ys, + (const void *)zs}; - const CCTK_INT input_array_indices[2] - = { real_idx, - imag_idx }; + const CCTK_INT input_array_indices[2] = {real_idx, imag_idx}; - const CCTK_INT output_array_types[2] - = { CCTK_VARIABLE_REAL, - CCTK_VARIABLE_REAL }; + const CCTK_INT output_array_types[2] = {CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL}; - void * output_arrays[2] - = { (void *) sphere_real, - (void *) sphere_imag }; + void *output_arrays[2] = {(void *)sphere_real, (void *)sphere_imag}; // const int operator_handle = CCTK_InterpHandle(interpolator_name); const int operator_handle = 0; // not used by CarpetX interpolator @@ -60,26 +47,17 @@ void Multipole_Interp(CCTK_ARGUMENTS, const int coord_system_handle = CCTK_CoordSystemHandle(coord_system); ierr = DriverInterpolate( - cctkGH, - num_dims, - operator_handle, - param_table_handle, + cctkGH, num_dims, operator_handle, param_table_handle, coord_system_handle, - CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0, // Only the 0 processor needs the points - CCTK_VARIABLE_REAL, - interp_coords, - num_input_arrays, - input_array_indices, - num_output_arrays, - output_array_types, - output_arrays); + CCTK_MyProc(cctkGH) == 0 ? (ntheta + 1) * (nphi + 1) + : 0, // Only the 0 processor needs the points + CCTK_VARIABLE_REAL, interp_coords, num_input_arrays, input_array_indices, + num_output_arrays, output_array_types, output_arrays); report_interp_error(ierr); - if (imag_idx == -1) - { - for (int i = 0; i < (ntheta+1)*(nphi+1); i++) - { + if (imag_idx == -1) { + for (int i = 0; i < (ntheta + 1) * (nphi + 1); i++) { sphere_imag[i] = 0; } } @@ -89,14 +67,14 @@ void Multipole_Interp(CCTK_ARGUMENTS, // // Debugging routine // void Multipole_InterpVar(CCTK_ARGUMENTS, -// CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name, -// CCTK_REAL sphere_var[]) +// CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const +// char *var_name, CCTK_REAL sphere_var[]) // { // DECLARE_CCTK_ARGUMENTS; // DECLARE_CCTK_PARAMETERS; // // Need parameters for the following: -// // ntheta (dtheta = pi/(ntheta) +// // ntheta (dtheta = pi/(ntheta) // // nphi (dphi = 2pi/(nphi) // // r (radius of sphere) // // NOTE: depending on the interval of integration, denominator above may @@ -107,7 +85,7 @@ void Multipole_Interp(CCTK_ARGUMENTS, // const CCTK_INT num_dims = 3; // int ierr = -1; -// const void* interp_coords[num_dims] +// const void* interp_coords[num_dims] // = { (const void *) x, // (const void *) y, // (const void *) z }; diff --git a/Multipole/src/interpolate.hh b/Multipole/src/interpolate.hh index 6668abc5..2a1e2b9f 100644 --- a/Multipole/src/interpolate.hh +++ b/Multipole/src/interpolate.hh @@ -7,13 +7,12 @@ #include "util_Table.h" // Multipole_Interp: -// This function interpolates psi4 onto the sphere in cartesian +// This function interpolates psi4 onto the sphere in cartesian // coordinates as created by Multipole_CoordSetup. -void Multipole_Interp(CCTK_ARGUMENTS, - CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], - int real_idx, int imag_idx, +void Multipole_Interp(CCTK_ARGUMENTS, CCTK_REAL x[], CCTK_REAL y[], + CCTK_REAL z[], int real_idx, int imag_idx, CCTK_REAL psi4r[], CCTK_REAL psi4i[]); -void Multipole_InterpVar(CCTK_ARGUMENTS, - CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name, +void Multipole_InterpVar(CCTK_ARGUMENTS, CCTK_REAL x[], CCTK_REAL y[], + CCTK_REAL z[], const char *var_name, CCTK_REAL sphere_var[]); diff --git a/Multipole/src/multipole.cc b/Multipole/src/multipole.cc index 9d3718b2..f3c0af8b 100644 --- a/Multipole/src/multipole.cc +++ b/Multipole/src/multipole.cc @@ -29,19 +29,16 @@ static const int max_spin_weights = max_vars; static const int max_l_modes = 10; static const int max_m_modes = 2 * max_l_modes + 1; -typedef struct -{ +typedef struct { int n_vars; Multipole::variable_desc *vars; -} -variables_desc; +} variables_desc; -static void fill_variable(int idx, const char *optstring, void *callback_arg) -{ +static void fill_variable(int idx, const char *optstring, void *callback_arg) { assert(idx >= 0); assert(callback_arg != 0); - variables_desc *vs = (variables_desc * ) callback_arg; + variables_desc *vs = (variables_desc *)callback_arg; assert(vs->n_vars < max_vars); // Too many variables in the variables list Multipole::variable_desc *v = &vs->vars[vs->n_vars]; @@ -54,28 +51,23 @@ static void fill_variable(int idx, const char *optstring, void *callback_arg) v->spin_weight = 0; v->name = string(CCTK_VarName(v->index)); - if (optstring != 0) - { + if (optstring != 0) { int table = Util_TableCreateFromString(optstring); - if (table >= 0) - { + if (table >= 0) { const int buffer_length = 256; char buffer[buffer_length]; - - Util_TableGetInt(table, &v->spin_weight , "sw"); - if (Util_TableGetString(table, buffer_length, buffer , "cmplx") >= 0) - { + + Util_TableGetInt(table, &v->spin_weight, "sw"); + if (Util_TableGetString(table, buffer_length, buffer, "cmplx") >= 0) { v->imag_index = CCTK_VarIndex(buffer); } - if (Util_TableGetString(table, buffer_length, buffer , "name") >= 0) - { + if (Util_TableGetString(table, buffer_length, buffer, "name") >= 0) { v->name = string(buffer); } const int ierr = Util_TableDestroy(table); - if (ierr) - { + if (ierr) { CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, "Could not destroy table: %d", ierr); } @@ -84,115 +76,106 @@ static void fill_variable(int idx, const char *optstring, void *callback_arg) vs->n_vars++; } -static void parse_variables_string(const string &var_string, Multipole::variable_desc v[max_vars], int *n_variables) -{ - variables_desc vars; +static void parse_variables_string(const string &var_string, + Multipole::variable_desc v[max_vars], + int *n_variables) { + variables_desc vars; vars.n_vars = 0; vars.vars = v; - int ierr = CCTK_TraverseString(var_string.c_str(), fill_variable, &vars, CCTK_GROUP_OR_VAR); + int ierr = CCTK_TraverseString(var_string.c_str(), fill_variable, &vars, + CCTK_GROUP_OR_VAR); assert(ierr >= 0); *n_variables = vars.n_vars; } -static void output_modes(CCTK_ARGUMENTS, const Multipole::variable_desc vars[], const CCTK_REAL radii[], - const Multipole::mode_array& modes) -{ +static void output_modes(CCTK_ARGUMENTS, const Multipole::variable_desc vars[], + const CCTK_REAL radii[], + const Multipole::mode_array &modes) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - if (output_ascii) - { - if (CCTK_MyProc(cctkGH) == 0) - { - for (int v = 0; v < modes.get_nvars(); v++) - { - for(int i=0; iindex)) << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; - Multipole_Output1D(CCTK_PASS_CTOC, real_base.str()+string(".th.asc"), array_size, th, ph, mp_theta, real); - Multipole_Output1D(CCTK_PASS_CTOC, real_base.str()+string(".ph.asc"), array_size, th, ph, mp_phi, real); - - if (v->imag_index != -1) - { + real_base << "mp_" << string(CCTK_VarName(v->index)) << "_r" + << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_Output1D(CCTK_PASS_CTOC, real_base.str() + string(".th.asc"), + array_size, th, ph, mp_theta, real); + Multipole_Output1D(CCTK_PASS_CTOC, real_base.str() + string(".ph.asc"), + array_size, th, ph, mp_phi, real); + + if (v->imag_index != -1) { ostringstream imag_base; - imag_base << "mp_" << string(CCTK_VarName(v->imag_index)) << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad; - Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str()+string(".th.asc"), array_size, th, ph, mp_theta, imag); - Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str()+string(".ph.asc"), array_size, th, ph, mp_phi, imag); + imag_base << "mp_" << string(CCTK_VarName(v->imag_index)) << "_r" + << setiosflags(ios::fixed) << setprecision(2) << rad; + Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str() + string(".th.asc"), + array_size, th, ph, mp_theta, imag); + Multipole_Output1D(CCTK_PASS_CTOC, imag_base.str() + string(".ph.asc"), + array_size, th, ph, mp_phi, imag); } } } } -bool int_in_array(int a, const int array[], int len) -{ - for (int i = 0; i < len; i++) - { +bool int_in_array(int a, const int array[], int len) { + for (int i = 0; i < len; i++) { if (array[i] == a) return true; } return false; } -int find_int_in_array(int a, const int array[], int len) -{ - for (int i = 0; i < len; i++) - { +int find_int_in_array(int a, const int array[], int len) { + for (int i = 0; i < len; i++) { if (array[i] == a) return i; } return -1; } - -static -void get_spin_weights(Multipole::variable_desc vars[], int n_vars, int spin_weights[max_spin_weights], int *n_weights) -{ +static void get_spin_weights(Multipole::variable_desc vars[], int n_vars, + int spin_weights[max_spin_weights], + int *n_weights) { int n_spin_weights = 0; - for (int i = 0; i < n_vars; i++) - { - if (!int_in_array(vars[i].spin_weight, spin_weights, n_spin_weights)) - { + for (int i = 0; i < n_vars; i++) { + if (!int_in_array(vars[i].spin_weight, spin_weights, n_spin_weights)) { assert(n_spin_weights < max_spin_weights); spin_weights[n_spin_weights] = vars[i].spin_weight; n_spin_weights++; @@ -203,81 +186,73 @@ void get_spin_weights(Multipole::variable_desc vars[], int n_vars, int spin_weig // For backward compatibility we allow the user to set l_mode instead // of l_max, but if it is left at the default of -1, l_max is used. -static int get_l_max() -{ +static int get_l_max() { DECLARE_CCTK_PARAMETERS; return l_mode == -1 ? l_max : l_mode; } -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++) - { +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; } } } - for (int si = 0; si < n_spin_weights; si++) - { + for (int si = 0; si < n_spin_weights; si++) { int sw = spin_weights[si]; - for (int l=0; l <= lmax; 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]; + for (int l = 0; l <= lmax; 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]; - 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]); } } - } + } } -extern "C" void Multipole_ParamCheck(CCTK_ARGUMENTS) -{ +extern "C" void Multipole_ParamCheck(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; DECLARE_CCTK_ARGUMENTS; - if (l_mode != -1) - { - CCTK_WARN(CCTK_WARN_ALERT, "The parameter l_mode is deprecated. Use l_max instead. For compatibility, l_max = l_mode is being used."); + if (l_mode != -1) { + CCTK_WARN(CCTK_WARN_ALERT, + "The parameter l_mode is deprecated. Use l_max instead. For " + "compatibility, l_max = l_mode is being used."); } - if (!CCTK_Equals(mode_type, "deprecated")) - { - CCTK_WARN(CCTK_WARN_ALERT, "The parameter mode_type is deprecated and is no longer used. All modes will be computed."); + if (!CCTK_Equals(mode_type, "deprecated")) { + CCTK_WARN(CCTK_WARN_ALERT, "The parameter mode_type is deprecated and is " + "no longer used. All modes will be computed."); } - if (l_min != -1) - { - CCTK_WARN(CCTK_WARN_ALERT, "The parameter l_min is deprecated and is no longer used. Modes from l = 0 will be computed."); + if (l_min != -1) { + CCTK_WARN(CCTK_WARN_ALERT, + "The parameter l_min is deprecated and is no longer used. Modes " + "from l = 0 will be computed."); } - if (m_mode != -100) - { - CCTK_WARN(CCTK_WARN_ALERT, "The parameter m_mode is deprecated. All m modes will be computed."); + if (m_mode != -100) { + CCTK_WARN( + CCTK_WARN_ALERT, + "The parameter m_mode is deprecated. All m modes will be computed."); } } -extern "C" void Multipole_Init(CCTK_ARGUMENTS) -{ - using namespace Loop; +extern "C" void Multipole_Init(CCTK_ARGUMENTS) { + using namespace Loop; DECLARE_CCTK_ARGUMENTS_Multipole_Init; DECLARE_CCTK_PARAMETERS; - const int dim = 3; + const int dim = 3; const array indextype = {0, 0, 0}; const GF3D2layout layout(cctkGH, indextype); @@ -286,18 +261,14 @@ extern "C" void Multipole_Init(CCTK_ARGUMENTS) const GF3D2 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 { - re_(p.I) = 0; - im_(p.I) = 0; - }); + grid.loop_int<0, 0, 0>(grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + re_(p.I) = 0; + im_(p.I) = 0; + }); } - - -extern "C" void Multipole_Calc(CCTK_ARGUMENTS) -{ +extern "C" void Multipole_Calc(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS_Multipole_Calc; DECLARE_CCTK_PARAMETERS; @@ -307,14 +278,14 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) 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 Multipole::variable_desc vars[max_vars]; + static Multipole::variable_desc vars[max_vars]; static int n_variables = 0; static int spin_weights[max_spin_weights]; static int n_spin_weights = 0; static bool initialized = false; - const int array_size=(ntheta+1)*(nphi+1); + const int array_size = (ntheta + 1) * (nphi + 1); if (out_every == 0 || cctk_iteration % out_every != 0) return; @@ -323,11 +294,10 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) assert(lmax + 1 <= max_l_modes); - if (!initialized) - { + if (!initialized) { real = new CCTK_REAL[array_size]; imag = new CCTK_REAL[array_size]; - th = new CCTK_REAL[array_size]; + th = new CCTK_REAL[array_size]; ph = new CCTK_REAL[array_size]; xs = new CCTK_REAL[array_size]; ys = new CCTK_REAL[array_size]; @@ -335,47 +305,45 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) xhat = new CCTK_REAL[array_size]; yhat = new CCTK_REAL[array_size]; zhat = new CCTK_REAL[array_size]; - + parse_variables_string(string(variables), vars, &n_variables); get_spin_weights(vars, n_variables, spin_weights, &n_spin_weights); Multipole_CoordSetup(xhat, yhat, zhat, th, ph); - setup_harmonics(spin_weights, n_spin_weights, lmax, th, ph, - array_size, reY, imY); + setup_harmonics(spin_weights, n_spin_weights, lmax, th, ph, array_size, reY, + imY); initialized = true; - CCTK_VINFO("initialized arrays"); + CCTK_VINFO("initialized arrays"); } Multipole::mode_array modes(n_variables, nradii, lmax); - for (int v = 0; v < n_variables; v++) - { - //assert(vars[v].spin_weight == -2); + for (int v = 0; v < n_variables; v++) { + // assert(vars[v].spin_weight == -2); - int si = find_int_in_array(vars[v].spin_weight, spin_weights, n_spin_weights); + int si = + find_int_in_array(vars[v].spin_weight, spin_weights, n_spin_weights); assert(si != -1); - for(int i=0; i= 0 && v < nvars); - assert(ri >= 0 && ri < nradii); - assert(l >= 0 && l <= lmax); - assert(m <= l && -m <= l); - return size_t(v * nradii * (lmax+1)*(lmax+1) * 2 + - ri * (lmax+1)*(lmax+1) * 2 + - (l*l + (m+l)) * 2 + is_im); - } + int get_nvars() const { return nvars; } + int get_nradii() const { return nradii; } + int get_lmax() const { return lmax; } - const int nvars; - const int nradii; - const int lmax; - std::vector modes; +private: + size_t mode_idx(int v, int ri, int l, int m, int is_im) const { + assert(v >= 0 && v < nvars); + assert(ri >= 0 && ri < nradii); + assert(l >= 0 && l <= lmax); + assert(m <= l && -m <= l); + return size_t(v * nradii * (lmax + 1) * (lmax + 1) * 2 + + ri * (lmax + 1) * (lmax + 1) * 2 + (l * l + (m + l)) * 2 + + is_im); + } + + const int nvars; + const int nradii; + const int lmax; + std::vector modes; }; -} +} // namespace Multipole diff --git a/Multipole/src/sphericalharmonic.cc b/Multipole/src/sphericalharmonic.cc index a5e89c31..e25459ce 100644 --- a/Multipole/src/sphericalharmonic.cc +++ b/Multipole/src/sphericalharmonic.cc @@ -11,103 +11,89 @@ static const CCTK_REAL PI = acos(-1.0); -static double factorial(int n) -{ +static double factorial(int n) { double returnval = 1; - for (int i = n; i >= 1; i--) - { + for (int i = n; i >= 1; i--) { returnval *= i; } return returnval; } -static inline double combination(int n, int m) -{ +static inline double combination(int n, int m) { // Binomial coefficient is undefined if these conditions do not hold assert(n >= 0); assert(m >= 0); assert(m <= n); - return factorial(n) / (factorial(m) * factorial(n-m)); + return factorial(n) / (factorial(m) * factorial(n - m)); } -static inline int imin(int a, int b) -{ - return a < b ? a : b; -} +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; -} +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 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); double all_coeff = 0, sum = 0; all_coeff = pow(-1.0, m); - all_coeff *= sqrt(factorial(l+m)*factorial(l-m)*(2*l+1) / (4.*PI*factorial(l+s)*factorial(l-s))); + all_coeff *= sqrt(factorial(l + m) * factorial(l - m) * (2 * l + 1) / + (4. * PI * factorial(l + s) * factorial(l - s))); sum = 0.; - for (int i = imax(m - s, 0); i <= imin(l + m, l - s); i++) - { - double sum_coeff = combination(l-s, i) * combination(l+s, i+s-m); - sum += sum_coeff * pow(-1.0, l-i-s) * pow(cos(th/2.), 2 * i + s - m) * - pow(sin(th/2.), 2*(l-i)+m-s); + for (int i = imax(m - s, 0); i <= imin(l + m, l - s); i++) { + double sum_coeff = combination(l - s, i) * combination(l + s, i + s - m); + sum += sum_coeff * pow(-1.0, l - i - s) * pow(cos(th / 2.), 2 * i + s - m) * + pow(sin(th / 2.), 2 * (l - i) + m - s); } - *reY = all_coeff*sum*cos(m*ph); - *imY = all_coeff*sum*sin(m*ph); + *reY = all_coeff * sum * cos(m * ph); + *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[]) -{ - for (int i = 0; i < array_size; i++) - { - Multipole_SphericalHarmonic(s,l,m,th[i],ph[i],&reY[i], &imY[i]); +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[]) { + for (int i = 0; i < array_size; i++) { + Multipole_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) -{ +extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) { using namespace Loop; using namespace std; DECLARE_CCTK_ARGUMENTS_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 { - CCTK_REAL vcoordr = sqrt(p.x*p.x + p.y*p.y + p.z*p.z); - CCTK_REAL theta = acos(p.z/vcoordr); - if (vcoordr == 0) theta = 0; - CCTK_REAL phi = atan2(p.y,p.x); - - 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_(p.I) = re * fac; - harmonic_im_(p.I) = im * fac; - }); + const GridDescBaseDevice grid(cctkGH); + grid.loop_int<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL vcoordr = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); + CCTK_REAL theta = acos(p.z / vcoordr); + if (vcoordr == 0) + theta = 0; + CCTK_REAL phi = atan2(p.y, p.x); + + 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_(p.I) = re * fac; + harmonic_im_(p.I) = im * fac; + }); // for (int k = 0; k < cctk_lsh[2]+1; k++) // { @@ -115,8 +101,10 @@ extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) // { // 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]); + // 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; @@ -138,53 +126,49 @@ extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) return; } - -extern "C" void Multipole_SetHarmonicWeyl(CCTK_ARGUMENTS) -{ +extern "C" void Multipole_SetHarmonicWeyl(CCTK_ARGUMENTS) { using namespace Loop; using namespace std; DECLARE_CCTK_ARGUMENTS_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 { - CCTK_REAL vcoordr = sqrt(p.x*p.x + p.y*p.y + p.z*p.z); - CCTK_REAL theta = acos(p.z/vcoordr); - if (vcoordr == 0) theta = 0; - CCTK_REAL phi = atan2(p.y,p.x); + grid.loop_int<0, 0, 0>( + grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL vcoordr = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); + CCTK_REAL theta = acos(p.z / vcoordr); + if (vcoordr == 0) + theta = 0; + CCTK_REAL phi = atan2(p.y, p.x); - CCTK_REAL re22 = 0; - CCTK_REAL im22 = 0; + CCTK_REAL re22 = 0; + CCTK_REAL im22 = 0; - Multipole_SphericalHarmonic(test_sw,2,2,theta,phi, - &re22, &im22); + Multipole_SphericalHarmonic(test_sw, 2, 2, theta, phi, &re22, &im22); - CCTK_REAL re2m2 = 0; - CCTK_REAL im2m2 = 0; + CCTK_REAL re2m2 = 0; + CCTK_REAL im2m2 = 0; - Multipole_SphericalHarmonic(test_sw,2,-2,theta,phi, - &re2m2, &im2m2); + Multipole_SphericalHarmonic(test_sw, 2, -2, theta, phi, &re2m2, &im2m2); - CCTK_REAL re31 = 0; - CCTK_REAL im31 = 0; + CCTK_REAL re31 = 0; + CCTK_REAL im31 = 0; - Multipole_SphericalHarmonic(test_sw,3,1,theta,phi, - &re31, &im31); + Multipole_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; - }); + 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; + }); return; } diff --git a/Multipole/src/sphericalharmonic.hh b/Multipole/src/sphericalharmonic.hh index 5c140947..987eb10e 100644 --- a/Multipole/src/sphericalharmonic.hh +++ b/Multipole/src/sphericalharmonic.hh @@ -3,16 +3,14 @@ #define __sphericalharmonic_h #include "cctk.h" -void Multipole_HarmonicSetup(int s, int l, int m, - int array_size, CCTK_REAL const th[], CCTK_REAL const 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 Multipole_SphericalHarmonic(int s, int l, int m, CCTK_REAL th, + CCTK_REAL ph, CCTK_REAL *reY, CCTK_REAL *imY); -void Multipole_SphericalHarmonic(int s, int l, int m, - CCTK_REAL th, CCTK_REAL ph, - CCTK_REAL *reY, CCTK_REAL *imY); - -void Multipole_SphericalHarmonicWeyl(int s, int l, int m, - CCTK_REAL th, CCTK_REAL ph, - CCTK_REAL *reY, CCTK_REAL *imY); +void Multipole_SphericalHarmonicWeyl(int s, int l, int m, CCTK_REAL th, + CCTK_REAL ph, CCTK_REAL *reY, + CCTK_REAL *imY); #endif diff --git a/Multipole/src/tests.cc b/Multipole/src/tests.cc index 9e2189e1..325363f1 100644 --- a/Multipole/src/tests.cc +++ b/Multipole/src/tests.cc @@ -15,69 +15,72 @@ static const int max_l_modes = 10; static const int max_m_modes = 2 * max_l_modes + 1; -static CCTK_REAL test_integral(int n, CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL),const int is_midpoint) -{ +static CCTK_REAL test_integral(int n, + CCTK_REAL (*integration_fn)(const CCTK_REAL *, + int, int, CCTK_REAL, + CCTK_REAL), + const int is_midpoint) { const int nx = n; const int ny = n; - const int array_size=(nx+1)*(ny+1); + const int array_size = (nx + 1) * (ny + 1); CCTK_REAL *f = new CCTK_REAL[array_size]; - const CCTK_REAL dx = 1./(nx + is_midpoint); - const CCTK_REAL dy = 1./(ny + is_midpoint); + const CCTK_REAL dx = 1. / (nx + is_midpoint); + const CCTK_REAL dy = 1. / (ny + is_midpoint); - for (int ix = 0; ix <= nx; ix++) - { - for (int iy = 0; iy <= ny; iy++) - { + 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 + 0.5*dx*is_midpoint; - const CCTK_REAL y = iy*dy + 0.5*dy*is_midpoint; + const CCTK_REAL x = ix * dx + 0.5 * dx * is_midpoint; + const CCTK_REAL y = iy * dy + 0.5 * dy * is_midpoint; 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); + f[i] = x * pow(y, 2) * pow(cos(2 * PI * y), 2) * pow(sin(2 * PI * x), 2); } } const CCTK_REAL result = integration_fn(f, nx, ny, dx, dy); - delete [] f; + delete[] f; return result; } -static CCTK_REAL test_pi_symmetric_sphere_integral(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL),const int is_midpoint) -{ +static CCTK_REAL test_pi_symmetric_sphere_integral( + CCTK_REAL (*integration_fn)(const CCTK_REAL *, int, int, CCTK_REAL, + CCTK_REAL), + const int is_midpoint) { const int n = 100; const int nth = n; const int nph = n; - const int array_size=(nth+1)*(nph+1); + const int array_size = (nth + 1) * (nph + 1); const CCTK_REAL PI = acos(-1.0); CCTK_REAL *f = new CCTK_REAL[array_size]; - const CCTK_REAL dth = PI/(nth + is_midpoint); - const CCTK_REAL dph = 2*PI/(nph + is_midpoint); + const CCTK_REAL dth = PI / (nth + is_midpoint); + const CCTK_REAL dph = 2 * PI / (nph + is_midpoint); - for (int ith = 0; ith <= nth; ith++) - { - for (int iph = 0; iph <= nph; iph++) - { + for (int ith = 0; ith <= nth; ith++) { + for (int iph = 0; iph <= nph; iph++) { const int i = Multipole_Index(ith, iph, nth); - const CCTK_REAL th = ith*dth + 0.5*dth*is_midpoint; - const CCTK_REAL ph = iph*dph + 0.5*dph*is_midpoint; + const CCTK_REAL th = ith * dth + 0.5 * dth * is_midpoint; + const CCTK_REAL ph = iph * dph + 0.5 * dph * is_midpoint; - f[i] = -(cos(ph)*sqrt(5/PI)*pow(cos(th/2.),3)*sin(th/2.)); + f[i] = -(cos(ph) * sqrt(5 / PI) * pow(cos(th / 2.), 3) * sin(th / 2.)); } } const CCTK_REAL result = integration_fn(f, nth, nph, dth, dph); - delete [] f; + delete[] f; return result; } -CCTK_REAL integration_convergence_order(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL), CCTK_REAL *store_low, CCTK_REAL *store_high, const int is_midpoint) -{ +CCTK_REAL integration_convergence_order( + CCTK_REAL (*integration_fn)(const CCTK_REAL *, int, int, CCTK_REAL, + CCTK_REAL), + CCTK_REAL *store_low, CCTK_REAL *store_high, const int is_midpoint) { const int n1 = 100; const int n2 = 200; const CCTK_REAL PI = acos(-1.0); @@ -85,35 +88,41 @@ CCTK_REAL integration_convergence_order(CCTK_REAL (*integration_fn) (const CCTK_ *store_low = result1; const CCTK_REAL result2 = test_integral(200, integration_fn, is_midpoint); *store_high = result2; - const CCTK_REAL exact = 1./24 + 1./(64 * pow(PI,2)); + 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); - return log10(error1/error2) / log10((CCTK_REAL) n2/n1); + return log10(error1 / error2) / log10((CCTK_REAL)n2 / n1); } -void Multipole_TestIntegrationConvergence(CCTK_ARGUMENTS) -{ +void Multipole_TestIntegrationConvergence(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; *test_simpson_convergence_order = integration_convergence_order( - &Simpson2DIntegral, test_simpson_result_low, test_simpson_result_high,0); + &Simpson2DIntegral, test_simpson_result_low, test_simpson_result_high, 0); *test_trapezoidal_convergence_order = integration_convergence_order( - &Trapezoidal2DIntegral, test_trapezoidal_result_low, test_trapezoidal_result_high,0); + &Trapezoidal2DIntegral, test_trapezoidal_result_low, + test_trapezoidal_result_high, 0); *test_midpoint_convergence_order = integration_convergence_order( - &Midpoint2DIntegral, test_midpoint_result_low, test_midpoint_result_high,1); + &Midpoint2DIntegral, test_midpoint_result_low, test_midpoint_result_high, + 1); } -void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS) -{ +void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; - *test_simpson_pi_symmetry = test_pi_symmetric_sphere_integral(&Simpson2DIntegral,0); - *test_midpoint_pi_symmetry = test_pi_symmetric_sphere_integral(&Midpoint2DIntegral,1); - *test_trapezoidal_pi_symmetry = test_pi_symmetric_sphere_integral(&Trapezoidal2DIntegral,0); - *test_driscollhealy_pi_symmetry = test_pi_symmetric_sphere_integral(&DriscollHealy2DIntegral,0); + *test_simpson_pi_symmetry = + test_pi_symmetric_sphere_integral(&Simpson2DIntegral, 0); + *test_midpoint_pi_symmetry = + test_pi_symmetric_sphere_integral(&Midpoint2DIntegral, 1); + *test_trapezoidal_pi_symmetry = + test_pi_symmetric_sphere_integral(&Trapezoidal2DIntegral, 0); + *test_driscollhealy_pi_symmetry = + test_pi_symmetric_sphere_integral(&DriscollHealy2DIntegral, 0); printf("Pi symmetry Simpson integral: %.19g\n", *test_simpson_pi_symmetry); printf("Pi symmetry midpoint integral: %.19g\n", *test_midpoint_pi_symmetry); - printf("Pi symmetry trapezoidal integral: %.19g\n", *test_trapezoidal_pi_symmetry); - printf("Pi symmetry Driscoll and Healy integral: %.19g\n", *test_driscollhealy_pi_symmetry); + printf("Pi symmetry trapezoidal integral: %.19g\n", + *test_trapezoidal_pi_symmetry); + printf("Pi symmetry Driscoll and Healy integral: %.19g\n", + *test_driscollhealy_pi_symmetry); } // void Multipole_TestIntegrate(CCTK_ARGUMENTS) @@ -142,27 +151,23 @@ void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS) // } // } - - - // 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) - - +// 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); // } -void Multipole_TestOrthonormality(CCTK_ARGUMENTS) -{ +void Multipole_TestOrthonormality(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS /* Campute Cartesian coordinates of points on the sphere */ - int array_size=(ntheta+1)*(nphi+1); + int array_size = (ntheta + 1) * (nphi + 1); CCTK_REAL *th = new CCTK_REAL[array_size]; CCTK_REAL *ph = new CCTK_REAL[array_size]; @@ -176,17 +181,14 @@ void Multipole_TestOrthonormality(CCTK_ARGUMENTS) CCTK_REAL *reY[1][max_l_modes][max_m_modes]; CCTK_REAL *imY[1][max_l_modes][max_m_modes]; - for (int sw = 0; sw <= 0; sw++) - { - for (int l = 0; l < max_l_modes; 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]; - - Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, - reY[sw][l][m+l], imY[sw][l][m+l]); + for (int sw = 0; sw <= 0; sw++) { + for (int l = 0; l < max_l_modes; 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]; + + Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[sw][l][m + l], + imY[sw][l][m + l]); } } } @@ -197,48 +199,39 @@ void Multipole_TestOrthonormality(CCTK_ARGUMENTS) // comparing each mode with each other but skipping the duplicates // gives N*(N+1)/2 // only 1 spin-weight mode is tested - const int N = max_l_modes*max_l_modes; + const int N = max_l_modes * max_l_modes; int idx = 0; - for (int sw = 0; sw <= 0; sw++) - { - for (int l = 0; l < max_l_modes; l++) - { - for (int m = -l; m <= l; m++) - { + for (int sw = 0; sw <= 0; sw++) { + for (int l = 0; l < max_l_modes; l++) { + for (int m = -l; m <= l; m++) { /* Compute scalar product of (real,imag) and all the Ylimi */ - for (int li = 0; li < max_l_modes; li++) - { - for (int mi = -li; mi <= li; mi++) - { + for (int li = 0; li < max_l_modes; li++) { + for (int mi = -li; mi <= li; mi++) { // only handle lower triangle in ((l,m),(li,mi)) space - if(l*l+l+m < li*li+li+mi) + if (l * l + l + m < li * li + li + mi) continue; CCTK_REAL real_lm = 0.0, imag_lm = 0.0; - Multipole_Integrate(array_size, ntheta, - reY[sw][li][mi+li], imY[sw][li][mi+li], - reY[sw][l][m+l], imY[sw][l][m+l], th, ph, - &real_lm, &imag_lm); + Multipole_Integrate(array_size, ntheta, reY[sw][li][mi + li], + imY[sw][li][mi + li], reY[sw][l][m + l], + imY[sw][l][m + l], th, ph, &real_lm, &imag_lm); - assert(idx < 1*N*(N+1)/2); - test_orthonormality[idx++] = sqrt(real_lm*real_lm + imag_lm*imag_lm); + assert(idx < 1 * N * (N + 1) / 2); + test_orthonormality[idx++] = + sqrt(real_lm * real_lm + imag_lm * imag_lm); } } - } } } - 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]; + 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]; } } } diff --git a/Multipole/src/utils.cc b/Multipole/src/utils.cc index 19ec6629..03edb942 100644 --- a/Multipole/src/utils.cc +++ b/Multipole/src/utils.cc @@ -28,21 +28,18 @@ #include "integrate.hh" #include "multipole.hh" -// check return code of HDF5 call abort with an error message if there was an error. -// adapted from CarpetIOHDF5/src/CarpetIOHDF5.hh. -#define HDF5_ERROR(fn_call) \ - do { \ - hid_t _error_code = fn_call; \ - \ - \ - if (_error_code < 0) \ - { \ - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, \ - "HDF5 call '%s' returned error code %d", \ - #fn_call, (int)_error_code); \ - } \ - } while (0) - +// check return code of HDF5 call abort with an error message if there was an +// error. adapted from CarpetIOHDF5/src/CarpetIOHDF5.hh. +#define HDF5_ERROR(fn_call) \ + do { \ + hid_t _error_code = fn_call; \ + \ + if (_error_code < 0) { \ + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, \ + "HDF5 call '%s' returned error code %d", #fn_call, \ + (int)_error_code); \ + } \ + } while (0) using namespace std; @@ -50,27 +47,27 @@ using namespace std; // File manipulation //////////////////////////////////////////////////////////////////////// -FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name) -{ +FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - + bool first_time = cctk_iteration == 0; const char *mode = first_time ? "w" : "a"; const char *my_out_dir = strcmp(out_dir, "") ? out_dir : io_out_dir; const int err = CCTK_CreateDirectory(0755, my_out_dir); if (err < 0) - CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, - "Multipole output directory %s could not be created (error code %d)", - my_out_dir, err); + CCTK_VWarn( + CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Multipole output directory %s could not be created (error code %d)", + my_out_dir, err); string output_name(string(my_out_dir) + "/" + name); FILE *fp = fopen(output_name.c_str(), mode); - if (fp == 0) - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "%s", (string("Could not open output file ") + output_name).c_str()); + if (fp == 0) { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "%s", + (string("Could not open output file ") + output_name).c_str()); } return fp; @@ -80,36 +77,34 @@ FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, const string &name) // Unused //////////////////////////////////////////////////////////////////////// -void Multipole_OutputArray(CCTK_ARGUMENTS, FILE *f, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], - CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[], - CCTK_REAL const data[]) -{ +void Multipole_OutputArray(CCTK_ARGUMENTS, FILE *f, int array_size, + CCTK_REAL const th[], CCTK_REAL const ph[], + CCTK_REAL const xs[], CCTK_REAL const ys[], + CCTK_REAL const zs[], CCTK_REAL const data[]) { DECLARE_CCTK_PARAMETERS; DECLARE_CCTK_ARGUMENTS; CCTK_REAL last_ph = ph[0]; - for (int i = 0; i < array_size; i++) - { + for (int i = 0; i < array_size; i++) { if (ph[i] != last_ph) // Separate blocks for gnuplot fprintf(f, "\n"); - fprintf(f, "%f %f %f %f %f %f %.19g\n", cctk_time, th[i], ph[i], xs[i], ys[i], zs[i], data[i]); + fprintf(f, "%f %f %f %f %f %f %.19g\n", cctk_time, th[i], ph[i], xs[i], + ys[i], zs[i], data[i]); last_ph = ph[i]; } } - -void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], - CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[], - CCTK_REAL const data[]) -{ +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, + int array_size, CCTK_REAL const th[], + CCTK_REAL const ph[], CCTK_REAL const xs[], + CCTK_REAL const ys[], CCTK_REAL const zs[], + CCTK_REAL const data[]) { DECLARE_CCTK_ARGUMENTS; - if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) - { - Multipole_OutputArray(CCTK_PASS_CTOC, fp, array_size, th, ph, xs, ys, zs, data); + if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) { + Multipole_OutputArray(CCTK_PASS_CTOC, fp, array_size, th, ph, xs, ys, zs, + data); fclose(fp); } } @@ -119,28 +114,21 @@ void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_s //////////////////////////////////////////////////////////////////////// void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord, - CCTK_REAL const data[]) -{ + CCTK_REAL const th[], CCTK_REAL const ph[], + mp_coord coord, CCTK_REAL const data[]) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - if (FILE *f = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) - { + if (FILE *f = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) { fprintf(f, "\"Time = %.19g\n", cctk_time); - if (coord == mp_theta) - { - for (int i = 0; i <= ntheta; i++) - { + if (coord == mp_theta) { + for (int i = 0; i <= ntheta; i++) { int idx = Multipole_Index(i, 0, ntheta); fprintf(f, "%f %.19g\n", th[idx], data[idx]); } - } - else if (coord == mp_phi) - { - for (int i = 0; i <= nphi; i++) - { + } else if (coord == mp_phi) { + for (int i = 0; i <= nphi; i++) { int idx = Multipole_Index(ntheta / 4, i, ntheta); fprintf(f, "%f %.19g\n", ph[idx], data[idx]); } @@ -154,19 +142,18 @@ void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, // Complex IO //////////////////////////////////////////////////////////////////////// -void Multipole_OutputComplex(CCTK_ARGUMENTS, FILE *fp, CCTK_REAL redata, CCTK_REAL imdata) -{ +void Multipole_OutputComplex(CCTK_ARGUMENTS, FILE *fp, CCTK_REAL redata, + CCTK_REAL imdata) { DECLARE_CCTK_PARAMETERS; DECLARE_CCTK_ARGUMENTS; fprintf(fp, "%f %.19g %.19g\n", cctk_time, redata, imdata); } -void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata) -{ +void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, + CCTK_REAL redata, CCTK_REAL imdata) { DECLARE_CCTK_ARGUMENTS; - if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) - { + if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) { Multipole_OutputComplex(CCTK_PASS_CTOC, fp, redata, imdata); fclose(fp); } @@ -178,14 +165,12 @@ void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL #ifdef HAVE_CAPABILITY_HDF5 -static bool file_exists(const string &name) -{ +static bool file_exists(const string &name) { struct stat sts; return !(stat(name.c_str(), &sts) == -1 && errno == ENOENT); } -static bool dataset_exists(hid_t file, const string &dataset_name) -{ +static bool dataset_exists(hid_t file, const string &dataset_name) { // To test whether a dataset exists, the recommended way in API 1.6 // is to use H5Gget_objinfo, but this prints an error to stderr if // the dataset does not exist. We explicitly avoid this by wrapping @@ -193,107 +178,102 @@ static bool dataset_exists(hid_t file, const string &dataset_name) // H5Gget_objinfo is deprecated, and H5Lexists does the job. See // http://www.mail-archive.com/hdf-forum@hdfgroup.org/msg00125.html - #if 1 +#if 1 bool exists; - H5E_BEGIN_TRY - { + H5E_BEGIN_TRY { exists = H5Gget_objinfo(file, dataset_name.c_str(), 1, NULL) >= 0; } H5E_END_TRY; return exists; - #else +#else return H5Lexists(file, dataset_name.c_str(), H5P_DEFAULT); - #endif +#endif } -void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const Multipole::variable_desc vars[], - const CCTK_REAL radii[], const Multipole::mode_array& modes) -{ +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, + const Multipole::variable_desc vars[], + const CCTK_REAL radii[], + const Multipole::mode_array &modes) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; const char *my_out_dir = strcmp(out_dir, "") ? out_dir : io_out_dir; if (CCTK_CreateDirectory(0755, my_out_dir) < 0) - CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, "Multipole output directory %s could not be created", my_out_dir); - static map checked; // Has the given file been checked - // for truncation? map<*,bool> - // defaults to false - for (int v = 0; v < modes.get_nvars(); v++) - { + static map checked; // Has the given file been checked + // for truncation? map<*,bool> + // defaults to false + for (int v = 0; v < modes.get_nvars(); v++) { string basename = "mp_" + vars[v].name + ".h5"; string output_name = my_out_dir + string("/") + basename; hid_t file; - if (!file_exists(output_name) || (!checked[output_name] && IO_TruncateOutputFiles(cctkGH))) - { - file = H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - } - else - { + if (!file_exists(output_name) || + (!checked[output_name] && IO_TruncateOutputFiles(cctkGH))) { + file = H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, + H5P_DEFAULT); + } else { file = H5Fopen(output_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT); } checked[output_name] = true; - for(int i=0; i ntheta = PI/dth - 1 - // th[i] = i * dth + 0.5*dth - // Therefore: - // th[0] = 0.5*dth <- GOOD. - // th[ntheta] = ntheta*dth + 0.5*dth - // = ((PI/dth)-1)*dth + 0.5*dth - // = PI - dth + 0.5*dth - // = PI - 0.5*dth <- GOOD. - // Similarly for ph. - - // Check for when midpoint disabled: - // dth = PI/ntheta -> ntheta = PI/dth - // th[i] = i * dth - // Therefore: - // th[0] = 0 <- GOOD. - // th[ntheta] = ntheta*dth - // = PI/dth*dth - // = PI <- GOOD. - // Similarly for ph. - th[i] = it * dth + 0.5*dth*is_midpoint; - ph[i] = ip * dph + 0.5*dph*is_midpoint; - xhat[i] = cos(ph[i])*sin(th[i]); - yhat[i] = sin(ph[i])*sin(th[i]); - zhat[i] = cos(th[i]); - } + const CCTK_REAL dth = PI / (ntheta + is_midpoint); + const CCTK_REAL dph = 2 * PI / (nphi + is_midpoint); + for (int it = 0; it <= ntheta; it++) { + for (int ip = 0; ip <= nphi; ip++) { + const int i = Multipole_Index(it, ip, ntheta); + // Check for when midpoint enabled: + // dth = PI/(ntheta+1) -> ntheta = PI/dth - 1 + // th[i] = i * dth + 0.5*dth + // Therefore: + // th[0] = 0.5*dth <- GOOD. + // th[ntheta] = ntheta*dth + 0.5*dth + // = ((PI/dth)-1)*dth + 0.5*dth + // = PI - dth + 0.5*dth + // = PI - 0.5*dth <- GOOD. + // Similarly for ph. + + // Check for when midpoint disabled: + // dth = PI/ntheta -> ntheta = PI/dth + // th[i] = i * dth + // Therefore: + // th[0] = 0 <- GOOD. + // th[ntheta] = ntheta*dth + // = PI/dth*dth + // = PI <- GOOD. + // Similarly for ph. + th[i] = it * dth + 0.5 * dth * is_midpoint; + ph[i] = ip * dph + 0.5 * dph * is_midpoint; + xhat[i] = cos(ph[i]) * sin(th[i]); + yhat[i] = sin(ph[i]) * sin(th[i]); + zhat[i] = cos(th[i]); } + } } void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r, - CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[], - CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[]) -{ - for (int it = 0; it <= ntheta; it++) - { - for (int ip = 0; ip <= nphi; ip++) - { + CCTK_REAL const xhat[], CCTK_REAL const yhat[], + CCTK_REAL const zhat[], CCTK_REAL x[], + CCTK_REAL y[], CCTK_REAL z[]) { + for (int it = 0; it <= ntheta; it++) { + for (int ip = 0; ip <= nphi; ip++) { const int i = Multipole_Index(it, ip, ntheta); x[i] = r * xhat[i]; @@ -389,18 +365,17 @@ void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r, // Integration //////////////////////////////////////////////////////////////////////// -void 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) -{ +void 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) { DECLARE_CCTK_PARAMETERS; - int il = Multipole_Index(0,0,ntheta); - int iu = Multipole_Index(1,0,ntheta); + int il = Multipole_Index(0, 0, ntheta); + int iu = Multipole_Index(1, 0, ntheta); CCTK_REAL dth = th[iu] - th[il]; - iu = Multipole_Index(0,1,ntheta); + iu = Multipole_Index(0, 1, ntheta); CCTK_REAL dph = ph[iu] - ph[il]; static CCTK_REAL *fr = 0; @@ -408,52 +383,40 @@ void Multipole_Integrate(int array_size, int nthetap, static bool allocated_memory = false; // Construct an array for the real integrand - if (!allocated_memory) - { + if (!allocated_memory) { fr = new CCTK_REAL[array_size]; fi = new CCTK_REAL[array_size]; allocated_memory = true; } - + // the below calculations take the integral of conj(array1)*array2*sin(th) - for (int i = 0; i < array_size; i++) - { - fr[i] = (array1r[i] * array2r[i] + - array1i[i] * array2i[i] ) * sin(th[i]); - fi[i] = (array1r[i] * array2i[i] - - array1i[i] * array2r[i] ) * sin(th[i]); + for (int i = 0; i < array_size; i++) { + fr[i] = (array1r[i] * array2r[i] + array1i[i] * array2i[i]) * sin(th[i]); + fi[i] = (array1r[i] * array2i[i] - array1i[i] * array2r[i]) * sin(th[i]); } - - if (CCTK_Equals(integration_method, "midpoint")) - { + + if (CCTK_Equals(integration_method, "midpoint")) { *outre = Midpoint2DIntegral(fr, ntheta, nphi, dth, dph); *outim = Midpoint2DIntegral(fi, ntheta, nphi, dth, dph); - } - else if (CCTK_Equals(integration_method, "trapezoidal")) - { + } else if (CCTK_Equals(integration_method, "trapezoidal")) { *outre = Trapezoidal2DIntegral(fr, ntheta, nphi, dth, dph); *outim = Trapezoidal2DIntegral(fi, ntheta, nphi, dth, dph); - } - 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"); + } 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(fr, ntheta, nphi, dth, dph); *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph); - } - else if (CCTK_Equals(integration_method, "DriscollHealy")) - { - if (ntheta % 2 != 0) - { - CCTK_WARN (CCTK_WARN_ABORT, "The Driscoll&Healy integration method requires even ntheta"); + } 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(fr, ntheta, nphi, dth, dph); *outim = DriscollHealy2DIntegral(fi, ntheta, nphi, dth, dph); - } - else - { + } else { CCTK_WARN(CCTK_WARN_ABORT, "internal error"); } -} +} diff --git a/Multipole/src/utils.hh b/Multipole/src/utils.hh index 68d59cd6..3d46ac2f 100644 --- a/Multipole/src/utils.hh +++ b/Multipole/src/utils.hh @@ -7,44 +7,47 @@ using namespace std; -enum mp_coord {mp_theta, mp_phi}; +enum mp_coord { mp_theta, mp_phi }; namespace Multipole { - class mode_array; - struct variable_desc; -} - -void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], - CCTK_REAL const x[], CCTK_REAL const y[], CCTK_REAL const z[], +class mode_array; +struct variable_desc; +} // namespace Multipole + +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, const string &name, + int array_size, CCTK_REAL const th[], + CCTK_REAL const ph[], CCTK_REAL const x[], + CCTK_REAL const y[], CCTK_REAL const z[], CCTK_REAL const data[]); void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size, - CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord, - CCTK_REAL const data[]); + CCTK_REAL const th[], CCTK_REAL const ph[], + mp_coord coord, CCTK_REAL const data[]); -void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata); +void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, + CCTK_REAL redata, CCTK_REAL imdata); -void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const Multipole::variable_desc vars[], const CCTK_REAL radii[], - const Multipole::mode_array& modes); +void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, + const Multipole::variable_desc vars[], + const CCTK_REAL radii[], + const Multipole::mode_array &modes); -void Multipole_CoordSetup(CCTK_REAL xhat[], CCTK_REAL yhat[], - CCTK_REAL zhat[], CCTK_REAL th[], - CCTK_REAL ph[]); +void Multipole_CoordSetup(CCTK_REAL xhat[], CCTK_REAL yhat[], CCTK_REAL zhat[], + CCTK_REAL th[], CCTK_REAL ph[]); void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r, - CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[], - CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[]); + CCTK_REAL const xhat[], CCTK_REAL const yhat[], + CCTK_REAL const zhat[], CCTK_REAL x[], + CCTK_REAL y[], CCTK_REAL z[]); -static inline int Multipole_Index(int it, int ip, int ntheta) -{ - return it + (ntheta+1)*ip; +static inline int Multipole_Index(int it, int ip, int ntheta) { + return it + (ntheta + 1) * ip; } -void Multipole_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 Multipole_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[]); #endif