Skip to content

Commit

Permalink
Multipole: beautify interpolate.cc
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwJi
lwJi committed Jun 23, 2024
1 parent 3cc5f29 commit 6807bc4
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 260 deletions.
170 changes: 0 additions & 170 deletions Multipole/src/interp_exp.cc
View file

This file was deleted.

137 changes: 47 additions & 90 deletions Multipole/src/interpolate.cc
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#include <stdio.h>

#include "interpolate.hh"
#include <loop_device.hxx>

static void report_interp_error(int ierr) {
if (ierr < 0) {
Expand All @@ -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);

Expand All @@ -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);
// }

0 comments on commit 6807bc4

Please sign in to comment.