diff options
Diffstat (limited to 'Carpet/CarpetInterp/src/interp.cc')
-rw-r--r-- | Carpet/CarpetInterp/src/interp.cc | 681 |
1 files changed, 58 insertions, 623 deletions
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index 8b48e9c54..47d5b9cb8 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -1,30 +1,22 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.36 2004/09/13 14:36:13 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.1 2003/04/30 12:37:31 schnetter Exp $ #include <assert.h> -#include <math.h> -#include <algorithm> +#include <complex> #include <vector> #include <mpi.h> #include "cctk.h" -#include "util_ErrorCodes.h" -#include "util_Table.h" +#include "Carpet/CarpetLib/src/vect.hh" -#include "bbox.hh" -#include "data.hh" -#include "defs.hh" -#include "ggf.hh" -#include "vect.hh" - -#include "carpet.hh" +#include "Carpet/Carpet/src/carpet.hh" #include "interp.hh" extern "C" { - static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.36 2004/09/13 14:36:13 schnetter Exp $"; + static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.1 2003/04/30 12:37:31 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc); } @@ -39,84 +31,6 @@ namespace CarpetInterp { typedef vect<int,dim> ivect; typedef vect<CCTK_REAL,dim> rvect; - typedef bbox<int,dim> ibbox; - - - -#define ind_rc(m,rl,c) ind_rc_(m,rl,minrl,maxrl,c,maxncomps,vhh) - static inline int ind_rc_(int const m, - int const rl, int const minrl, int const maxrl, - int const c, int const maxncomps, - vector<gh<dim> *> const hh) - { - assert (m>=0 && m<maps); - assert (rl>=minrl && rl<maxrl); - assert (minrl>=0 && maxrl<=hh.at(m)->reflevels()); - assert (c>=0 && c<maxncomps); - assert (maxncomps>=0 && maxncomps<=hh.at(m)->components(rl)); - int const ind = (rl-minrl) * maxncomps + c; - assert (ind>=0 && ind < (maxrl-minrl) * maxncomps); - return ind; - } - -#define ind_prc(p,m,rl,c) \ - ind_prc_(p,nprocs,m,rl,minrl,maxrl,c,maxncomps,cgh,vhh) - static inline int ind_prc_(int const p, int const nprocs, - int const m, - int const rl, int const minrl, int const maxrl, - int const c, int const maxncomps, - cGH const * const cgh, - vector<gh<dim> *> const hh) - { - assert (p>=0 && p<nprocs); - assert (nprocs==CCTK_nProcs(cgh)); - assert (m>=0 && m<maps); - assert (rl>=minrl && rl<maxrl); - assert (minrl>=0 && maxrl<=hh.at(m)->reflevels()); - assert (c>=0 && c<maxncomps); - assert (maxncomps>=0 && maxncomps<=hh.at(m)->components(rl)); - int const ind = (p * (maxrl-minrl) + (rl-minrl)) * maxncomps + c; - assert (ind>=0 && ind < nprocs * (maxrl-minrl) * maxncomps); - return ind; - } - - - static int GetInterpNumTimelevels(const cGH * const cgh, - const int group) - { - assert (group>=0 && group<CCTK_NumGroups()); - - int ierr; - - cGroup gp; - ierr = CCTK_GroupData (group, &gp); - assert (!ierr); - - int interp_ntl = -1; - - const int length = Util_TableGetInt - (gp.tagstable, &interp_ntl, "InterpNumTimelevels"); - if (length == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - interp_ntl = -1; - } else if (length >= 0) { - // all is well - just check they're not asking too much - if (interp_ntl > CCTK_ActiveTimeLevelsGI(cgh, group)) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The tags table entry \"InterpNumTimelevels\" " - "for group \"%s\" says that %d timelevels should " - "be used for time interpolation.\n However, only %d " - "are active.", - groupname, interp_ntl, - CCTK_ActiveTimeLevelsGI(cgh, group)); - } - } else { - assert(0); - } - - return interp_ntl; - - } void CarpetInterpStartup () @@ -126,7 +40,7 @@ namespace CarpetInterp { - int InterpGridArrays (cGH const * const cgh, + int InterpGridArrays (cGH const * const cGH, int const N_dims, int const local_interp_handle, int const param_table_handle, @@ -140,556 +54,77 @@ namespace CarpetInterp { CCTK_INT const output_array_type_codes [], void * const output_arrays []) { - if (CCTK_IsFunctionAliased ("SymmetryInterpolate")) { - return SymmetryInterpolate - (cgh, N_dims, - local_interp_handle, param_table_handle, coord_system_handle, - N_interp_points, interp_coords_type_code, interp_coords, - N_input_arrays, input_array_variable_indices, - N_output_arrays, output_array_type_codes, output_arrays); - } else { - return Carpet_DriverInterpolate - (cgh, N_dims, - local_interp_handle, param_table_handle, coord_system_handle, - N_interp_points, interp_coords_type_code, interp_coords, - N_input_arrays, input_array_variable_indices, - N_output_arrays, output_array_type_codes, output_arrays); - } - } - - - - extern "C" CCTK_INT - Carpet_DriverInterpolate (CCTK_POINTER_TO_CONST const cgh_, - CCTK_INT const N_dims, - CCTK_INT const local_interp_handle, - CCTK_INT const param_table_handle, - CCTK_INT const coord_system_handle, - CCTK_INT const N_interp_points, - CCTK_INT const interp_coords_type_code, - CCTK_POINTER_TO_CONST const interp_coords [], - CCTK_INT const N_input_arrays, - CCTK_INT const input_array_variable_indices [], - CCTK_INT const N_output_arrays, - CCTK_INT const output_array_type_codes [], - CCTK_POINTER const output_arrays []) - { - cGH const * const cgh = static_cast<cGH const *> (cgh_); int ierr; - assert (cgh); + assert (cGH); assert (N_dims == dim); - - - if (is_meta_mode()) { - CCTK_WARN (0, "It is not possible to interpolate in meta mode"); - } - - bool const want_global_mode = is_global_mode(); - - - - // Get some information - MPI_Comm const comm = CarpetMPIComm (); - int const myproc = CCTK_MyProc (cgh); - int const nprocs = CCTK_nProcs (cgh); - assert (myproc>=0 && myproc<nprocs); - - // Multiple maps are not supported - // (because we don't know how to select a map) - assert (maps == 1); - const int m = 0; - - int const minrl = want_global_mode ? 0 : reflevel; - int const maxrl = want_global_mode ? vhh.at(m)->reflevels() : reflevel+1; - - // Multiple convergence levels are not supported - assert (mglevels == 1); - int const ml = 0; - - int maxncomps = 0; - for (int rl=minrl; rl<maxrl; ++rl) { - maxncomps = max(maxncomps, vhh.at(m)->components(rl)); - } - - - - // Find the time interpolation order - int partype; - void const * const parptr - = CCTK_ParameterGet ("prolongation_order_time", "Carpet", &partype); - assert (parptr); - assert (partype == PARAMETER_INTEGER); - int const prolongation_order_time = * (CCTK_INT const *) parptr; - - CCTK_REAL const current_time = cgh->cctk_time / cgh->cctk_delta_time; - - - - // Find out about the coordinates: origin and delta for the Carpet - // grid indices -#if 0 const char * coord_system_name = CCTK_CoordSystemName (coord_system_handle); assert (coord_system_name); - rvect lower, upper, delta; + rvect lower, upper; for (int d=0; d<dim; ++d) { ierr = CCTK_CoordRange - (cgh, &lower[d], &upper[d], d+1, 0, coord_system_name); + (cGH, &lower[d], &upper[d], d+1, 0, coord_system_name); assert (!ierr); - const ibbox & baseext = vhh.at(m)->baseextent; - delta[d] - = (upper[d] - lower[d]) / (baseext.upper()[d] - baseext.lower()[d]); - } -#else - rvect const lower = rvect::ref(cgh->cctk_origin_space); - rvect const delta = rvect::ref(cgh->cctk_delta_space) / maxreflevelfact; - rvect const upper = lower + delta * (vhh.at(m)->baseextent.upper() - vhh.at(m)->baseextent.lower()); -#endif - - assert (N_interp_points >= 0); - assert (interp_coords); - for (int d=0; d<dim; ++d) { - assert (N_interp_points==0 || interp_coords[d]); } - assert (N_output_arrays >= 0); - if (N_interp_points > 0) { - assert (output_arrays); - for (int j=0; j<N_output_arrays; ++j) { - assert (output_arrays[j]); - } - } - - + const ivect lbnd (cGH->cctk_lbnd); + const ivect lsh (cGH->cctk_lsh ); + const ivect gsh (cGH->cctk_gsh ); - // Assign interpolation points to components - vector<int> rlev (N_interp_points); // refinement level of point n - vector<int> home (N_interp_points); // component of point n - vector<int> homecnts ((maxrl-minrl) * maxncomps); // number of points in component c - - for (int n=0; n<N_interp_points; ++n) { - - // Find the grid point closest to the interpolation point - assert (interp_coords_type_code == CCTK_VARIABLE_REAL); - rvect pos; - for (int d=0; d<dim; ++d) { - pos[d] = static_cast<CCTK_REAL const *>(interp_coords[d])[n]; - } - - // Find the component that this grid point belongs to - rlev.at(n) = -1; - home.at(n) = -1; - if (all(pos>=lower && pos<=upper)) { - for (int rl=maxrl-1; rl>=minrl; --rl) { - - const int fact = maxreflevelfact / ipow(reffact, rl) * ipow(mgfact, mglevel); - ivect const ipos = ivect(floor((pos - lower) / (delta * fact) + 0.5)) * fact; - assert (all(ipos % vhh.at(m)->bases.at(rl).at(ml).stride() == 0)); - - // TODO: use something faster than a linear search - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - if (vhh.at(m)->extents.at(rl).at(c).at(ml).contains(ipos)) { - rlev.at(n) = rl; - home.at(n) = c; - goto found; - } - } - } - } - CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation point #%d at [%g,%g,%g] is not on any grid patch", - n, pos[0], pos[1], pos[2]); - // Map the point (arbitrarily) to the first component of the - // coarsest grid - rlev.at(n) = minrl; - home.at(n) = 0; - found: - assert (rlev.at(n)>=minrl && rlev.at(n)<maxrl); - assert (home.at(n)>=0 && home.at(n)<vhh.at(m)->components(rlev.at(n))); - ++ homecnts.at(ind_rc(m, rlev.at(n), home.at(n))); - - } // for n - - // Communicate counts - vector<int> allhomecnts(nprocs * (maxrl-minrl) * maxncomps); - MPI_Allgather - (&homecnts .at(0), (maxrl-minrl) * maxncomps, MPI_INT, - &allhomecnts.at(0), (maxrl-minrl) * maxncomps, MPI_INT, comm); - - - - // Create coordinate patches - vector<data<CCTK_REAL,dim> > allcoords - (nprocs * (maxrl-minrl) * maxncomps); - for (int p=0; p<nprocs; ++p) { - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - ivect lo (0); - ivect up (1); - up[0] = allhomecnts.at(ind_prc(p,m,rl,c)); - up[1] = dim; - ivect str (1); - ibbox extent (lo, up-str, str); - allcoords.at(ind_prc(p,m,rl,c)).allocate (extent, p); - } - } - } - - // Fill in local coordinate patches - { - vector<int> tmpcnts ((maxrl-minrl) * maxncomps); - for (int n=0; n<N_interp_points; ++n) { - int const rl = rlev.at(n); - int const c = home.at(n); - assert (rl>=minrl && rl<maxrl); - assert (c>=0 && c<vhh.at(m)->components(rl)); - assert (tmpcnts.at(ind_rc(m,rl,c)) >= 0); - assert (tmpcnts.at(ind_rc(m,rl,c)) < homecnts.at(ind_rc(m,rl,c))); - assert (dim==3); - for (int d=0; d<dim; ++d) { - allcoords.at(ind_prc(myproc,m,rl,c))[ivect(tmpcnts.at(ind_rc(m,rl,c)),d,0)] - = static_cast<CCTK_REAL const *>(interp_coords[d])[n]; - } - ++ tmpcnts.at(c + (rl-minrl)*maxncomps); - } - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - assert (tmpcnts.at(ind_rc(m,rl,c)) == homecnts.at(ind_rc(m,rl,c))); - } - } - } - - // Transfer coordinate patches - for (comm_state<dim> state; !state.done(); state.step()) { - for (int p=0; p<nprocs; ++p) { - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - allcoords.at(ind_prc(p,m,rl,c)).change_processor - (state, vhh.at(m)->processors.at(rl).at(c)); - } - } - } - } - - - - // Create output patches - vector<data<CCTK_REAL,dim> > alloutputs - (nprocs * (maxrl-minrl) * maxncomps, -1); - for (int p=0; p<nprocs; ++p) { - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - ivect lo (0); - ivect up (1); - up[0] = allhomecnts.at(ind_prc(p,m,rl,c)); - up[1] = N_output_arrays; - ivect str (1); - ibbox extent (lo, up-str, str); - alloutputs.at(ind_prc(p,m,rl,c)).allocate - (extent, vhh.at(m)->processors.at(rl).at(c)); - } - } - } + const rvect coord_delta + = (upper - lower) / rvect((gsh - 1) * reflevelfact); + const rvect coord_origin = lower + rvect(lbnd) * coord_delta; - - - // - // Do the local interpolation - // - int overall_ierr = 0; - BEGIN_GLOBAL_MODE(cgh) { - for (int rl=minrl; rl<maxrl; ++rl) { - enter_level_mode (const_cast<cGH*>(cgh), rl); + vector<CCTK_INT> input_array_type_codes(N_input_arrays); + vector<const void *> input_arrays(N_input_arrays); + for (int n=0; n<N_input_arrays; ++n) { + if (input_array_variable_indices[n] == -1) { - // Number of necessary time levels - CCTK_REAL const level_time = cgh->cctk_time / cgh->cctk_delta_time; - bool const need_time_interp - = (fabs(current_time - level_time) - > 1e-12 * (fabs(level_time) + fabs(current_time) - + fabs(cgh->cctk_delta_time))); - assert (! (! want_global_mode && need_time_interp)); - int const num_tl = need_time_interp ? prolongation_order_time + 1 : 1; + // Ignore this entry + input_array_type_codes[n] = -1; + input_arrays[n] = 0; + + } else { + + int const vi = input_array_variable_indices[n]; + assert (vi>=0 && vi<CCTK_NumVars()); + + int const gi = CCTK_GroupIndexFromVarI (vi); + assert (gi>=0 && gi<CCTK_NumGroups()); + + cGroup group; + ierr = CCTK_GroupData (gi, &group); + assert (!ierr); + assert (group.grouptype == CCTK_GF); + assert (group.dim == dim); + assert (group.disttype == CCTK_DISTRIB_DEFAULT); + assert (group.stagtype == 0); // not staggered + + input_array_type_codes[n] = group.vartype; + input_arrays[n] = CCTK_VarDataPtrI (cGH, 0, vi); - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - - // Find out about the local geometry - ivect lsh; - rvect coord_origin, coord_delta; - for (int d=0; d<dim; ++d) { - lsh[d] = cgh->cctk_lsh[d]; - coord_delta[d] = cgh->cctk_delta_space[d] / cgh->cctk_levfac[d]; - coord_origin[d] = cgh->cctk_origin_space[d] + (cgh->cctk_lbnd[d] + 1.0 * cgh->cctk_levoff[d] / cgh->cctk_levoffdenom[d]) * coord_delta[d]; - } - - - - // Find out about the grid functions - vector<CCTK_INT> input_array_type_codes(N_input_arrays); - vector<const void *> input_arrays(N_input_arrays); - for (int n=0; n<N_input_arrays; ++n) { - if (input_array_variable_indices[n] == -1) { - - // Ignore this entry - input_array_type_codes.at(n) = -1; - - } else { - - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - - int const gi = CCTK_GroupIndexFromVarI (vi); - assert (gi>=0 && gi<CCTK_NumGroups()); - - cGroup group; - ierr = CCTK_GroupData (gi, &group); - assert (!ierr); - assert (group.grouptype == CCTK_GF); - assert (group.dim == dim); - assert (group.disttype == CCTK_DISTRIB_DEFAULT); - assert (group.stagtype == 0); // not staggered - - int interp_ntls = GetInterpNumTimelevels(cgh, gi); - // If -ve then key wasn't set; use all timelevels. - if (interp_ntls < 0) interp_ntls = group.numtimelevels; - - // TODO: emit better warning - if ((num_tl > group.numtimelevels)&& - (interp_ntls > group.numtimelevels)) { - CCTK_VWarn(0, __LINE__,__FILE__,"CarpetInterp", - "Tried to interpolate variable '%s' " - "in time.\nIt has insufficient timelevels " - "(%d are required).", - CCTK_FullName(vi), - min(num_tl,interp_ntls)); - } - assert (group.numtimelevels >= min(num_tl,interp_ntls)); - - input_array_type_codes.at(n) = group.vartype; - - } - } // for input arrays - - - - // Work on the data from all processors - for (int p=0; p<nprocs; ++p) { - assert (allcoords.at(ind_prc(p,m,reflevel,component)).owns_storage()); - assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == allcoords.at(ind_prc(p,m,reflevel,component)).shape()[0]); - assert (allhomecnts.at(ind_prc(p,m,reflevel,component)) == alloutputs.at(ind_prc(p,m,reflevel,component)).shape()[0]); - - int const npoints = allhomecnts.at(ind_prc(p,m,reflevel,component)); - - // Do the processor-local interpolation - vector<const void *> tmp_interp_coords (dim); - for (int d=0; d<dim; ++d) { - tmp_interp_coords.at(d) = &allcoords.at(ind_prc(p,m,reflevel,component))[ivect(0,d,0)]; - } - - - - vector<vector<void *> > tmp_output_arrays (num_tl); - - for (int tl=0; tl<num_tl; ++tl) { - - for (int n=0; n<N_input_arrays; ++n) { - - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - - int const gi = CCTK_GroupIndexFromVarI (vi); - assert (gi>=0 && gi<CCTK_NumGroups()); - - int interp_ntls = GetInterpNumTimelevels(cgh, gi); - // If -ve then key wasn't set; use all timelevels. - if (interp_ntls < 0) interp_ntls = num_tl; - - if (input_array_variable_indices[n] == -1) { - input_arrays.at(n) = 0; - } else if (tl > interp_ntls-1) { - // Do the interpolation anyway, just from - // an earlier timelevel. - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - input_arrays.at(n) = - CCTK_VarDataPtrI (cgh, interp_ntls-1, vi); - } else { - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - input_arrays.at(n) = CCTK_VarDataPtrI (cgh, tl, vi); - } - } // for input arrays - - tmp_output_arrays.at(tl).resize (N_output_arrays); - for (int j=0; j<N_output_arrays; ++j) { - assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); - if (need_time_interp) { - tmp_output_arrays.at(tl).at(j) = new CCTK_REAL [npoints]; - } else { - tmp_output_arrays.at(tl).at(j) = &alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(0,j,0)]; - } - } - - ierr = CCTK_InterpLocalUniform - (N_dims, local_interp_handle, param_table_handle, - &coord_origin[0], &coord_delta[0], - npoints, - interp_coords_type_code, &tmp_interp_coords[0], - N_input_arrays, &lsh[0], - &input_array_type_codes[0], &input_arrays[0], - N_output_arrays, - output_array_type_codes, &tmp_output_arrays.at(tl)[0]); - if (ierr) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "The local interpolator returned the error code %d", ierr); - } - overall_ierr = min(overall_ierr, ierr); - - } // for tl - - // Interpolate in time, if necessary - if (need_time_interp) { - - for (int j=0; j<N_output_arrays; ++j) { - - // Find output variable indices - vector<CCTK_INT> operand_indices (N_output_arrays); - ierr = Util_TableGetIntArray - (param_table_handle, N_output_arrays, - &operand_indices.front(), "operand_indices"); - if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - assert (N_output_arrays == N_input_arrays); - for (int m=0; m<N_output_arrays; ++m) { - operand_indices.at(m) = m; - } - } else { - assert (ierr == N_output_arrays); - } - - // Find input array for this output array - int const n = operand_indices.at(j); - - int const vi = input_array_variable_indices[n]; - assert (vi>=0 && vi<CCTK_NumVars()); - - int const gi = CCTK_GroupIndexFromVarI (vi); - assert (gi>=0 && gi<CCTK_NumGroups()); - - int interp_ntls = GetInterpNumTimelevels(cgh, gi); - // If -ve then key wasn't set; use all timelevels. - if (interp_ntls < 0) interp_ntls = num_tl; - - // Get interpolation times - CCTK_REAL const time = current_time; - vector<CCTK_REAL> times(interp_ntls); - for (int tl=0; tl<interp_ntls; ++tl) { - times.at(tl) = vtt.at(m)->time (-tl, reflevel, mglevel); - } - - // Calculate interpolation weights - vector<CCTK_REAL> tfacs(interp_ntls); - switch (interp_ntls) { - case 1: - // no interpolation - // We have to assume that any GF with one timelevel - // is constant in time!!! - // assert (fabs((time - times.at(0)) / fabs(time + times.at(0) + cgh->cctk_delta_time)) < 1e-12); - tfacs.at(0) = 1.0; - break; - case 2: - // linear (2-point) interpolation - tfacs.at(0) = (time - times.at(1)) / (times.at(0) - times.at(1)); - tfacs.at(1) = (time - times.at(0)) / (times.at(1) - times.at(0)); - break; - case 3: - // quadratic (3-point) interpolation - tfacs.at(0) = (time - times.at(1)) * (time - times.at(2)) / ((times.at(0) - times.at(1)) * (times.at(0) - times.at(2))); - tfacs.at(1) = (time - times.at(0)) * (time - times.at(2)) / ((times.at(1) - times.at(0)) * (times.at(1) - times.at(2))); - tfacs.at(2) = (time - times.at(0)) * (time - times.at(1)) / ((times.at(2) - times.at(0)) * (times.at(2) - times.at(1))); - break; - default: - assert (0); - } - - // Interpolate - assert (output_array_type_codes[j] == CCTK_VARIABLE_REAL); - for (int k=0; k<npoints; ++k) { - CCTK_REAL & dest = alloutputs.at(ind_prc(p,m,reflevel,component))[ivect(k,j,0)]; - dest = 0; - for (int tl=0; tl<interp_ntls; ++tl) { - CCTK_REAL const src = ((CCTK_REAL const *)tmp_output_arrays.at(tl).at(j))[k]; - dest += tfacs[tl] * src; - } - } - } - - for (int tl=0; tl<num_tl; ++tl) { - for (int j=0; j<N_output_arrays; ++j) { - delete [] (CCTK_REAL *)tmp_output_arrays.at(tl).at(j); - } - } - - } // if need_time_interp - - } // for processors - - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - - leave_level_mode (const_cast<cGH*>(cgh)); - } // for rl - - } END_GLOBAL_MODE; - - - - // Transfer output patches back - for (comm_state<dim> state; !state.done(); state.step()) { - for (int p=0; p<nprocs; ++p) { - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - alloutputs.at(ind_prc(p,m,rl,c)).change_processor (state, p); - } - } - } - } - - - - // Read out local output patches - { - vector<int> tmpcnts ((maxrl-minrl) * maxncomps); - for (int n=0; n<N_interp_points; ++n) { - int const rl = rlev.at(n); - int const c = home.at(n); - for (int j=0; j<N_output_arrays; ++j) { - assert (interp_coords_type_code == CCTK_VARIABLE_REAL); - assert (alloutputs.at(ind_prc(myproc,m,rl,c)).owns_storage()); - assert (output_arrays); - assert (output_arrays[j]); - static_cast<CCTK_REAL *>(output_arrays[j])[n] = alloutputs.at(ind_prc(myproc,m,rl,c))[ivect(tmpcnts.at(ind_rc(m,rl,c)),j,0)]; - } - ++ tmpcnts.at(ind_rc(m,rl,c)); - } - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<vhh.at(m)->components(rl); ++c) { - assert (tmpcnts.at(ind_rc(m,rl,c)) == homecnts.at(ind_rc(m,rl,c))); - } } } - - - int global_overall_ierr; - MPI_Allreduce - (&overall_ierr, &global_overall_ierr, 1, MPI_INT, MPI_MIN, comm); - - - - // Done. - return global_overall_ierr; + ierr = CCTK_InterpLocalUniform (N_dims, + local_interp_handle, + param_table_handle, + &coord_origin[0], + &coord_delta[0], + N_interp_points, + interp_coords_type_code, + interp_coords, + N_input_arrays, + &lsh[0], + &input_array_type_codes[0], + &input_arrays[0], + N_output_arrays, + output_array_type_codes, + output_arrays); + + return ierr; } } // namespace CarpetInterp |