diff options
Diffstat (limited to 'Carpet/CarpetInterp/src/interp.cc')
| -rw-r--r-- | Carpet/CarpetInterp/src/interp.cc | 695 |
1 files changed, 695 insertions, 0 deletions
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc new file mode 100644 index 000000000..8b48e9c54 --- /dev/null +++ b/Carpet/CarpetInterp/src/interp.cc @@ -0,0 +1,695 @@ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.36 2004/09/13 14:36:13 schnetter Exp $ + +#include <assert.h> +#include <math.h> + +#include <algorithm> +#include <vector> + +#include <mpi.h> + +#include "cctk.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "bbox.hh" +#include "data.hh" +#include "defs.hh" +#include "ggf.hh" +#include "vect.hh" + +#include "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 $"; + CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc); +} + + + +namespace CarpetInterp { + + using namespace Carpet; + + + + 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 () + { + CCTK_OverloadInterpGridArrays (InterpGridArrays); + } + + + + int InterpGridArrays (cGH const * const cgh, + int const N_dims, + int const local_interp_handle, + int const param_table_handle, + int const coord_system_handle, + int const N_interp_points, + int const interp_coords_type_code, + void const * const interp_coords [], + int const N_input_arrays, + CCTK_INT const input_array_variable_indices [], + int const N_output_arrays, + 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 (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; + for (int d=0; d<dim; ++d) { + ierr = CCTK_CoordRange + (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]); + } + } + + + + // 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)); + } + } + } + + + + // + // 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); + + // 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; + + 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; + } + +} // namespace CarpetInterp |