From 6807bc4598ace1cb99191638827984bf5bf10aa6 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Sun, 23 Jun 2024 00:04:40 +0000 Subject: [PATCH] Multipole: beautify interpolate.cc --- Multipole/src/interp_exp.cc | 170 ----------------------------------- Multipole/src/interpolate.cc | 137 ++++++++++------------------ 2 files changed, 47 insertions(+), 260 deletions(-) delete mode 100644 Multipole/src/interp_exp.cc diff --git a/Multipole/src/interp_exp.cc b/Multipole/src/interp_exp.cc deleted file mode 100644 index faf23127..00000000 --- a/Multipole/src/interp_exp.cc +++ /dev/null @@ -1,170 +0,0 @@ -#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; - - 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 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; - - 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; - } - - 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); - if (param_table_handle < 0) - 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) - 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) - CCTK_VERROR("Can't set operation_codes array in parameter table: %d", ierr); - - const int coord_system_handle = 0; // not used by CarpetX Interpolator - - // 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); - - 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++) { - 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 - -// 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/Multipole/src/interpolate.cc b/Multipole/src/interpolate.cc index c495470c..76fefa13 100644 --- a/Multipole/src/interpolate.cc +++ b/Multipole/src/interpolate.cc @@ -1,7 +1,6 @@ #include #include "interpolate.hh" -#include static void report_interp_error(int ierr) { if (ierr < 0) { @@ -16,43 +15,59 @@ void Multipole_Interp(CCTK_ARGUMENTS, CCTK_REAL xs[], CCTK_REAL ys[], 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 N_dims = 3; + const CCTK_INT npoints = + CCTK_MyProc(cctkGH) == 0 ? (ntheta + 1) * (nphi + 1) : 0; + const void *interp_coords[N_dims] = {(const void *)xs, (const void *)ys, + (const void *)zs}; + CCTK_INT N_input_arrays = imag_idx == -1 ? 1 : 2; + CCTK_INT N_output_arrays = imag_idx == -1 ? 1 : 2; const CCTK_INT input_array_indices[2] = {real_idx, imag_idx}; + // Interpolation result + CCTK_POINTER output_arrays[2]; + output_arrays[0] = sphere_real; + output_arrays[1] = sphere_imag; - 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); - const int operator_handle = 0; // not used by CarpetX interpolator - + int ierr = -1; + /* DriverInterpolate arguments that aren't currently used */ + const int coord_system_handle = 0; + CCTK_INT const interp_coords_type_code = 0; + CCTK_INT const output_array_types[1] = {0}; + int interp_handle = CCTK_InterpHandle("CarpetX"); + + if (interp_handle < 0) + CCTK_VERROR("Could not obtain inteprolator handle for built-in 'CarpetX' " + "interpolator: %d", + interp_handle); + + // Interpolation parameter table + CCTK_INT operations[1][N_dims]; + for (int var = 0; var < npoints; var++) { + operations[0][var] = 0; + } + int operands[1][N_dims]; + for (int var = 0; var < npoints; var++) { + operands[0][var] = var; + } 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); + if (param_table_handle < 0) + 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, N_input_arrays, + (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, N_input_arrays, + (int const *const)operations, + "operation_codes")) < 0) + CCTK_VERROR("Can't set operation_codes array in parameter table: %d", ierr); ierr = DriverInterpolate( - 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); + cctkGH, N_dims, interp_handle, param_table_handle, coord_system_handle, + npoints, interp_coords_type_code, interp_coords, N_input_arrays, + input_array_indices, N_output_arrays, output_array_types, output_arrays); report_interp_error(ierr); @@ -64,61 +79,3 @@ void Multipole_Interp(CCTK_ARGUMENTS, CCTK_REAL xs[], CCTK_REAL ys[], 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 - -// 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); -// }