diff options
Diffstat (limited to 'src/multipole.cc')
| -rw-r--r-- | src/multipole.cc | 320 |
1 files changed, 320 insertions, 0 deletions
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 <stdio.h> +#include <string.h> +#include <assert.h> + +#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<nradii; i++) + { + // Compute x^i = r * \hat x^i + Multipole_ScaleCartesian(ntheta, nphi, radius[i], xhat, yhat, zhat, xs, ys, zs); + + // Interpolate Psi4r and Psi4i + Multipole_Interp(CCTK_PASS_CTOC, xs, ys, zs, vars[v].index, vars[v].imag_index, + real, imag); + for (int l=lmin; l <= lmax; l++) + { + mmin = (allmodes ? (m_mode > 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 +} |