From d47c3dc7324ddc372d3eeeab1366a38064fd70db Mon Sep 17 00:00:00 2001 From: knarf Date: Fri, 26 Mar 2010 10:25:54 +0000 Subject: use trunc structure git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@53 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843 --- interface.ccl | 10 ++ param.ccl | 111 ++++++++++++++++ schedule.ccl | 20 +++ src/integrate.cc | 152 ++++++++++++++++++++++ src/integrate.hh | 9 ++ src/interpolate.cc | 146 +++++++++++++++++++++ src/interpolate.hh | 19 +++ src/make.code.defn | 1 + src/multipole.cc | 320 +++++++++++++++++++++++++++++++++++++++++++++++ src/multipole.hh | 15 +++ src/sphericalharmonic.cc | 97 ++++++++++++++ src/sphericalharmonic.hh | 11 ++ src/utils.cc | 251 +++++++++++++++++++++++++++++++++++++ src/utils.hh | 40 ++++++ 14 files changed, 1202 insertions(+) create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/integrate.cc create mode 100644 src/integrate.hh create mode 100644 src/interpolate.cc create mode 100644 src/interpolate.hh create mode 100644 src/make.code.defn create mode 100644 src/multipole.cc create mode 100644 src/multipole.hh create mode 100644 src/sphericalharmonic.cc create mode 100644 src/sphericalharmonic.hh create mode 100644 src/utils.cc create mode 100644 src/utils.hh diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..463a8ce --- /dev/null +++ b/interface.ccl @@ -0,0 +1,10 @@ +#interface.ccl for thorn Multipole +implements: multipole +inherits: Grid + +public: + +CCTK_REAL harmonics type=GF timelevels=1 +{ + harmonic_re, harmonic_im +} "Spherical harmonics" diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..e752c8d --- /dev/null +++ b/param.ccl @@ -0,0 +1,111 @@ +#param.ccl for thorn Multipole + +restricted: +#Interpolator params: +CCTK_STRING interpolator_name "Which interpolator should I use" +{ + ".+" :: "Any nonempty string" +} "Hermite polynomial interpolation" + +CCTK_STRING interpolator_pars "Parameters for the interpolator" +{ + ".*" :: "Any string that Util_TableSetFromString() will take" +} "order=3" + +CCTK_STRING coord_system "What is the coord system?" +{ + ".*" :: "Any smart string will do" +} "cart3d" + +KEYWORD integration_method "How to do surface integrals" STEERABLE=always +{ + "midpoint" :: "Midpoint rule (1st order)" + "Simpson" :: "Simpson's rule (4th order)" +} "midpoint" + +CCTK_INT out_every "How often to output" \ +STEERABLE=recover +{ + * :: "Any integer" +} 1 + +CCTK_INT out_1d_every "How often to output 1d data" \ +STEERABLE=recover +{ + * :: "Any integer" +} 0 + +#physical params: +CCTK_INT nradii "How many extraction radii?" \ +STEERABLE=recover +{ + 0:101 :: "A positive integer less than 101" +} 1 + + +CCTK_REAL radius[101] "The radii for extraction" \ +STEERABLE=recover +{ + 0.0:* :: "Please keep it in the grid" +} 0.0 + +CCTK_INT ntheta "How many points in the theta direction?" \ +STEERABLE=recover +{ + 0:* :: "Positive please" +} 50 + +CCTK_INT nphi "How many points in the phi direction?" \ +STEERABLE=recover +{ + 0:* :: "Positive please" +} 100 + +CCTK_STRING variables "What variables to decompose" +{ + ".*" :: "A list of variables" +} "" + +CCTK_BOOLEAN verbose "Output detailed information about what is happening" +{ +} "no" + +KEYWORD mode_type "Which type of mode extraction do we have" STEERABLE=always +{ + "all modes" :: "Extract all modes up to (l_mode, m_mode)." + "specific mode" :: "Select one specific (l_mode, m_mode) mode" +} "all modes" + +CCTK_INT l_min "all modes: above which l mode to calculate/ specific mode: which l mode to extract" STEERABLE=always +{ + 0:* :: "l >= 0" +} 2 + +CCTK_INT l_mode "all modes: Up to which l mode to calculate/ specific mode: which l mode to extract" STEERABLE=always +{ + 0:* :: "l >= 0" +} 2 + +CCTK_INT m_mode "all modes: Up to which m mode to calculate/ specific mode: which m mode to extract " STEERABLE=always +{ + 0:* :: "Positive Please" +} 2 + +CCTK_BOOLEAN enable_test "whether to set a spherical harmonic in the 'harmonic' grid functions" +{ +} "no" + +CCTK_INT test_l "which mode to put into the test variables" +{ + * :: "Any integer" +} 2 + +CCTK_INT test_m "which mode to put into the test variables" +{ + * :: "Any integer" +} 2 + +CCTK_INT test_sw "which spin weight to put into the test variables" +{ + * :: "Any integer" +} -2 diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..a7123ee --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,20 @@ +#schedule.ccl for thorn Multipole + +if (enable_test) +{ + STORAGE: harmonics[1] +} + +schedule Multipole_Calc at CCTK_ANALYSIS after (calc_np,PsiKadelia) +{ + LANG: C + OPTIONS: GLOBAL +} "Calculate Multipoles" + +if (enable_test) +{ + schedule Multipole_SetHarmonic at CCTK_INITIAL + { + LANG: C + } "Populate grid functions with spherical harmonics" +} diff --git a/src/integrate.cc b/src/integrate.cc new file mode 100644 index 0000000..0f1c392 --- /dev/null +++ b/src/integrate.cc @@ -0,0 +1,152 @@ +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +/* + +We will want to integrate functions F(th,ph) from th = 0 to pi, ph = 0 +to 2 pi with a weighting function sin(th). Alternatively, we might +want to use u = cos(th) as the variable, in which case we will go from +u = -1 to 1 and ph = 0 to 2 pi. For simplicity, we implement an +integration routine with a weight function of 1, and require the user +to multiply the integrand by their own weight function. We divide the +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. + +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 +spacing hx and hy, but these are just multiplied by the result to give +the integral. + +*/ + + + +#define idx(xx,yy) (assert((xx) <= nx), assert((xx) >= 0), assert((yy) <= ny), assert((yy) >= 0), ((xx) + (yy) * (nx+1))) + +// Hard coded 2D integrals + +static 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); + + // Corners + 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 (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)]; + + 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 integrand_sum = 0; + int ix = 0, iy = 0; + + assert(nx > 0); assert(ny > 0); assert (f); + assert(nx % 2 == 0); + assert(ny % 2 == 0); + + int px = nx / 2; + int py = ny / 2; + + // Corners + 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)]; + + 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)]; + + 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)]; + + for (iy = 1; iy <= py-1; iy++) + for (ix = 1; ix <= px; ix++) + 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 (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; +} + +// 1D integrals + +static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h) +{ + CCTK_REAL integrand_sum = 0; + int i = 0; + + assert(n > 0); assert(f); + assert(n % 2 == 0); + + int p = n / 2; + + integrand_sum += f[0] + f[n]; + + for (i = 1; i <= p-1; i++) + integrand_sum += 4 * f[2*i-1] + 2 * f[2*i]; + + integrand_sum += 4 * f[2*p-1]; + + return 1.0/3.0 * h * integrand_sum; +} + +// 2D integral built up from 1D + +static CCTK_REAL Composite2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy) +{ + CCTK_REAL integrand_sum = 0; + + assert(nx > 0); assert(ny > 0); assert (f); + assert(nx % 2 == 0); + assert(ny % 2 == 0); + + CCTK_REAL *g = new CCTK_REAL[ny+1]; + + for (int i = 0; i <= ny; i++) + { + g[i] = Simpson1DIntegral(&f[idx(0,i)], nx, hx); + } + + integrand_sum = Simpson1DIntegral(g, ny, hy); + delete [] g; + return integrand_sum; +} diff --git a/src/integrate.hh b/src/integrate.hh new file mode 100644 index 0000000..c45917e --- /dev/null +++ b/src/integrate.hh @@ -0,0 +1,9 @@ + +#ifndef __integrate_h +#define __integrate_h +#include "cctk.h" + +CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, + CCTK_REAL hx, CCTK_REAL hy); + +#endif diff --git a/src/interpolate.cc b/src/interpolate.cc new file mode 100644 index 0000000..5c9946c --- /dev/null +++ b/src/interpolate.cc @@ -0,0 +1,146 @@ +#include + +#include "interpolate.hh" + +static +void report_interp_error(int ierr) +{ + if (ierr < 0) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "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[]) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // Need parameters for the following: + // 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_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 CCTK_INT input_array_indices[2] + = { real_idx, + imag_idx }; + + const CCTK_INT output_array_types[2] + = { CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL }; + + void * output_arrays[2] + = { (void *) sphere_real, + (void *) sphere_imag }; + + const int operator_handle = CCTK_InterpHandle(interpolator_name); + + int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); + ierr = Util_TableSetFromString(param_table_handle, interpolator_pars); + + const int coord_system_handle = CCTK_CoordSystemHandle(coord_system); + + ierr = CCTK_InterpGridArrays( + 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); + + report_interp_error(ierr); + + if (imag_idx == -1) + { + for (int i = 0; i < (ntheta+1)*(nphi+1); i++) + { + sphere_imag[i] = 0; + } + } + + Util_TableDestroy(param_table_handle); +} + +// // Debugging routine +// void Multipole_InterpVar(CCTK_ARGUMENTS, +// 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) +// // 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 + +// // printf("Multipole: Interpolating %s\n", var_name); + +// const CCTK_INT num_input_arrays = 1; +// const CCTK_INT num_output_arrays = 1; +// const CCTK_INT num_dims = 3; +// int ierr = -1; + +// const void* interp_coords[num_dims] +// = { (const void *) x, +// (const void *) y, +// (const void *) z }; + +// const CCTK_INT input_array_indices[num_input_arrays] +// = { CCTK_VarIndex(var_name) }; + +// const CCTK_INT output_array_types[num_output_arrays] +// = { CCTK_VARIABLE_REAL }; + +// void * output_arrays[num_output_arrays] +// = { (void *) sphere_var }; + +// const int operator_handle = CCTK_InterpHandle(interpolator_name); + +// int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); +// ierr = Util_TableSetFromString(param_table_handle, interpolator_pars); + +// const int coord_system_handle = CCTK_CoordSystemHandle(coord_system); + +// ierr = CCTK_InterpGridArrays( +// cctkGH, +// num_dims, +// operator_handle, +// param_table_handle, +// coord_system_handle, +// (ntheta+1)*(nphi+1), +// CCTK_VARIABLE_REAL, +// interp_coords, +// num_input_arrays, +// input_array_indices, +// num_output_arrays, +// output_array_types, +// output_arrays); +// report_interp_error(ierr); +// } diff --git a/src/interpolate.hh b/src/interpolate.hh new file mode 100644 index 0000000..6668abc --- /dev/null +++ b/src/interpolate.hh @@ -0,0 +1,19 @@ +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" +#include "util_Table.h" + +// Multipole_Interp: +// 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, + CCTK_REAL psi4r[], CCTK_REAL psi4i[]); + +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/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..cccefbd --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1 @@ +SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc diff --git a/src/multipole.cc b/src/multipole.cc new file mode 100644 index 0000000..30a6624 --- /dev/null +++ b/src/multipole.cc @@ -0,0 +1,320 @@ +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +#include "multipole.hh" +#include "interpolate.hh" +#include "utils.hh" +#include "sphericalharmonic.hh" + +static const int var_name_length = 30; +static const int max_vars = 10; + +static const int max_spin_weights = 3; +static const int max_l_modes = 10; +static const int max_m_modes = 2 * max_l_modes + 1; + +typedef struct +{ + int index; + int imag_index; + int spin_weight; + char name[var_name_length]; +} +variable_desc; + +typedef struct +{ + int n_vars; + variable_desc *vars; +} +variables_desc; + + +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; + + assert(vs->n_vars < max_vars); // Too many variables in the variables list + variable_desc *v = &vs->vars[vs->n_vars]; + + v->index = idx; + + // Default values if there is no option string or if the options are + // not present + v->imag_index = -1; + v->spin_weight = 0; + strcpy(v->name, CCTK_VarName(v->index)); + + if (optstring != 0) + { + int table = Util_TableCreateFromString(optstring); + + if (table >= 0) + { + const int buffer_length = 256; + char imag_name[buffer_length]; + + Util_TableGetInt(table, &v->spin_weight , "sw"); +///////////////////////////////////////////////////////////// +CCTK_VInfo(CCTK_THORNSTRING,"spinweight %d", v->spin_weight); +///////////////////////////////////////////////////////////// + if (Util_TableGetString(table, buffer_length, imag_name , "cmplx") >= 0) + { + v->imag_index = CCTK_VarIndex(imag_name); + } + Util_TableGetString(table, var_name_length, v->name , "name"); + } + } + vs->n_vars++; +} + +static void parse_variables_string(const char *var_string, variable_desc v[max_vars], int *n_variables) +{ + variables_desc vars; + + vars.n_vars = 0; + vars.vars = v; + + const char *actual_var_string = strlen(var_string) == 0 ? + "WeylScal4::Psi4r{sw=-2 cmplx=\'WeylScal4::Psi4i\' name=\'psi4\'}" : + var_string; + + if (strlen(var_string) == 0) + { + static bool have_warned = false; + if (!have_warned) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "WARNING: No Multipole::variables string has been specified. For compatibility, the current default is to use \"%s\" but in future, this default will change to \"\". Please update your parameter files to include this new parameter.", actual_var_string); + have_warned = true; + } + } + + int ierr = CCTK_TraverseString(actual_var_string, fill_variable, &vars, CCTK_GROUP_OR_VAR); + assert(ierr > 0); + + *n_variables = vars.n_vars; +} + +static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, + CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + char name_tmp[1000]; + + if (CCTK_MyProc(cctkGH) == 0) + { + sprintf(name_tmp, "mp_%s_l%d_m%d_r%1.2f.asc", v->name, l, m, rad); + Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name_tmp, real_lm, imag_lm); + } +} + +static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, + CCTK_REAL *th, CCTK_REAL *ph, + CCTK_REAL *real, CCTK_REAL *imag, int array_size) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + char name_tmp[1000]; + + if (CCTK_MyProc(cctkGH) == 0) + { + if (out_1d_every != 0 && (cctk_iteration) % out_1d_every == 0) + { + const char *real_name = CCTK_VarName(v->index); + const char *imag_name = CCTK_VarName(v->imag_index); + + sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", real_name, rad); + Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, real); + sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", imag_name, rad); + Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, imag); + + sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", real_name, rad); + Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, real); + sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", imag_name, rad); + Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, 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++) + { + 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++) + { + if (array[i] == a) + return i; + } + return -1; +} + + +static +void get_spin_weights(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)) + { + assert(n_spin_weights < max_spin_weights); + spin_weights[n_spin_weights] = vars[i].spin_weight; + n_spin_weights++; + } + } + *n_weights = n_spin_weights; +} + + +static +void setup_harmonics(const int spin_weights[max_spin_weights], int n_spin_weights, int lmin, int lmax, + int m_mode, bool allmodes, 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++) + { + int sw = spin_weights[si]; + + for (int l=lmin; l <= lmax; l++) + { + int mmin = (allmodes ? (m_mode > l ? -l : -m_mode) : m_mode); + int mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode); + + for (int m=mmin; m <= mmax; m++) + { + reY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size]; + imY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size]; + + Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin]); + } + } + } +} + +extern "C" void Multipole_Calc(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + 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 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; + + int i, array_size=(ntheta+1)*(nphi+1); + CCTK_REAL real_lm = 0.0, imag_lm = 0.0; + int lmin, lmax, mmin, mmax, allmodes=0; + + if (cctk_iteration % out_every != 0) + return; + + // Determine modes to compute + if (CCTK_EQUALS(mode_type, "all modes")) + { + allmodes=1; + } + //if mode_type is set to "all modes", all modes with l <= l_mode, + //abs(m) <= min{l_mode, m_mode} will be output. Otherwise, only + //the l_mode, m_mode mode will be output. + lmin = (allmodes ? l_min : l_mode); + lmax = l_mode; + + assert(lmax + 1 <= max_l_modes); + + if (!initialized) + { + printf("Multipole: Allocating memory in Multipole_Calc\n"); + 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]; + ys = new CCTK_REAL[array_size]; + zs = new CCTK_REAL[array_size]; + xhat = new CCTK_REAL[array_size]; + yhat = new CCTK_REAL[array_size]; + zhat = new CCTK_REAL[array_size]; + + printf("Multipole: parsing variables string.\n"); + parse_variables_string(variables, vars, &n_variables); + get_spin_weights(vars, n_variables, spin_weights, &n_spin_weights); + Multipole_CoordSetup(ntheta, nphi, xhat, yhat, zhat, th, ph); + setup_harmonics(spin_weights, n_spin_weights, lmin, lmax, m_mode, allmodes, th, ph, array_size, reY, imY); + initialized = true; + } + + 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); + assert(si != -1); + + for(i=0; i l ? -l : -m_mode) : m_mode) ; + mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode); + + for (int m=mmin; m <= mmax; m++) + { + // Integrate sYlm (real + i imag) over the sphere at radius r + Multipole_Integrate(array_size, ntheta, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin], real, imag, th, ph, + &real_lm, &imag_lm); + + output_mode(CCTK_PASS_CTOC, &vars[v], radius[i], real_lm, imag_lm, l, m); + }//loop over m + }//loop over l + output_1D(CCTK_PASS_CTOC, &vars[v], radius[i], th, ph, real, imag, array_size); + }//loop over radii + }//loop over variables +} diff --git a/src/multipole.hh b/src/multipole.hh new file mode 100644 index 0000000..004d351 --- /dev/null +++ b/src/multipole.hh @@ -0,0 +1,15 @@ +#include "cctk.h" +#include "cctk_Arguments.h" + +// Multipole_Calc +// This is the main scheduling file. Because we are completely local here +// and do not use cactus arrays etc, we schedule only one function and then +// like program like one would in C, C++ with this function taking the +// place of int main(void). +// +// This function calls functions to accomplish 3 things: +// 1) Interpolate psi4 onto a sphere +// 2) Integrate psi4 with the ylm's over that sphere +// 2) Output the mode decomposed psi4 +extern "C" void Multipole_Calc(CCTK_ARGUMENTS); + diff --git a/src/sphericalharmonic.cc b/src/sphericalharmonic.cc new file mode 100644 index 0000000..116613c --- /dev/null +++ b/src/sphericalharmonic.cc @@ -0,0 +1,97 @@ +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +static const CCTK_REAL PI = acos(-1.0); + +static double factorial(int n) +{ + double returnval = 1; + for (int i = n; i >= 1; i--) + { + returnval *= i; + } + return returnval; +} + +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)); +} + +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 +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))); + 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); + } + *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]); + } +} + + +// Fill a grid function with a given spherical harmonic +extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + for (int k = 0; k < cctk_lsh[2]; k++) + { + for (int j = 0; j < cctk_lsh[1]; j++) + { + for (int i = 0; i < cctk_lsh[0]; i++) + { + int index = CCTK_GFINDEX3D(cctkGH,i,j,k) ; + + CCTK_REAL theta = acos(z[index]/r[index]); + CCTK_REAL phi = atan2(y[index],x[index]); + + Multipole_SphericalHarmonic(test_sw,test_l,test_m,theta,phi, + &harmonic_re[index], &harmonic_im[index]); + } + } + } + return; +} diff --git a/src/sphericalharmonic.hh b/src/sphericalharmonic.hh new file mode 100644 index 0000000..32b1bba --- /dev/null +++ b/src/sphericalharmonic.hh @@ -0,0 +1,11 @@ + +#ifndef __sphericalharmonic_h +#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[], + CCTK_REAL reY[], CCTK_REAL imY[]); + + +#endif diff --git a/src/utils.cc b/src/utils.cc new file mode 100644 index 0000000..00c4752 --- /dev/null +++ b/src/utils.cc @@ -0,0 +1,251 @@ +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "utils.hh" +#include "integrate.hh" + +//////////////////////////////////////////////////////////////////////// +// File manipulation +//////////////////////////////////////////////////////////////////////// + +FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, char const *name) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + const int buf_size = 1024; + bool first_time = cctk_iteration == 0; + const char *mode = first_time ? "w" : "a"; + char output_name[buf_size]; + CCTK_STRING *out_dir = (CCTK_STRING *) CCTK_ParameterGet("out_dir", "IOUtil", NULL); + + if (*out_dir == 0) + { + CCTK_WARN(1,"Parameter IOUtil::out_dir not found"); + return 0; + } + + if ((int)strlen(*out_dir)+1 > buf_size) + { + CCTK_WARN(1,"Parameter IOUtil::out_dir is too long"); + return 0; + } + + sprintf(output_name, "%s/%s", *out_dir, name); + + FILE *fp = fopen(output_name, mode); + + if (fp == 0) + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "Could not open output file %s", output_name); + } + + return fp; +} + +//////////////////////////////////////////////////////////////////////// +// 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[]) +{ + DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + + CCTK_REAL last_ph = ph[0]; + + 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]); + last_ph = ph[i]; + } +} + + +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char *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); + fclose(fp); + } +} + +//////////////////////////////////////////////////////////////////////// +// Misc +//////////////////////////////////////////////////////////////////////// + +void Multipole_Output1D(CCTK_ARGUMENTS, char const *name, int array_size, + 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)) + { + fprintf(f, "\"Time = %.19g\n", cctk_time); + + 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++) + { + int idx = Multipole_Index(ntheta / 4, i, ntheta); + fprintf(f, "%f %.19g\n", ph[idx], data[idx]); + } + } + fprintf(f, "\n\n"); + fclose(f); + } +} + +//////////////////////////////////////////////////////////////////////// +// Complex IO +//////////////////////////////////////////////////////////////////////// + +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, char const *name, CCTK_REAL redata, CCTK_REAL imdata) +{ + DECLARE_CCTK_ARGUMENTS; + + if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name)) + { + Multipole_OutputComplex(CCTK_PASS_CTOC, fp, redata, imdata); + fclose(fp); + } +} + + +//////////////////////////////////////////////////////////////////////// +// Coordinates +//////////////////////////////////////////////////////////////////////// + +void Multipole_CoordSetup(int ntheta, int nphi, + CCTK_REAL xhat[], CCTK_REAL yhat[], + CCTK_REAL zhat[], CCTK_REAL th[], + CCTK_REAL ph[]) +{ + const CCTK_REAL PI = acos(-1.0); + + for (int it = 0; it <= ntheta; it++) + { + for (int ip = 0; ip <= nphi; ip++) + { + const int i = Multipole_Index(it, ip, ntheta); + + th[i] = it * PI / (ntheta); + ph[i] = ip * 2 * PI / nphi; + 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++) + { + const int i = Multipole_Index(it, ip, ntheta); + + x[i] = r * xhat[i]; + y[i] = r * yhat[i]; + z[i] = r * zhat[i]; + } + } +} + +//////////////////////////////////////////////////////////////////////// +// 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) +{ + DECLARE_CCTK_PARAMETERS + + 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); + CCTK_REAL dph = ph[iu] - ph[il]; + + if (CCTK_Equals(integration_method, "midpoint")) + { + CCTK_REAL tempr = 0.0; + CCTK_REAL tempi = 0.0; + + for (int i=0; i < array_size; i++) { + // the below calculations take the integral of conj(array1)*array2 + tempr += ( array1r[i]*array2r[i] + array1i[i]*array2i[i] ) + *sin(th[i])*dth*dph; + tempi += ( array1r[i]*array2i[i] - array1i[i]*array2r[i] ) + *sin(th[i])*dth*dph; + + *outre = tempr; + *outim = tempi; + } + } + else if (CCTK_Equals(integration_method, "simpson")) + { + static CCTK_REAL *fr = 0; + static CCTK_REAL *fi = 0; + static bool allocated_memory = false; + + // Construct an array for the real integrand + if (!allocated_memory) + { + printf("Multipole: Allocating memory for integrand\n"); + fr = new CCTK_REAL[array_size]; + fi = new CCTK_REAL[array_size]; + allocated_memory = true; + } + + 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]); + } + + *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph); + *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph); + } +} diff --git a/src/utils.hh b/src/utils.hh new file mode 100644 index 0000000..2f7d2be --- /dev/null +++ b/src/utils.hh @@ -0,0 +1,40 @@ + +#ifndef __utils_h +#define __utils_h + +#include "cctk.h" + +enum mp_coord {mp_theta, mp_phi}; + +void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char const *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, char const *name, int array_size, + CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord, + CCTK_REAL const data[]); + +void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, char const *name, CCTK_REAL redata, CCTK_REAL imdata); + +void Multipole_CoordSetup(int ntheta, int nphi, + 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[]); + +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[]); + +#endif -- cgit v1.2.3