diff --git a/CarpetX/src/driver.cxx b/CarpetX/src/driver.cxx index ec535a563..a139a12fc 100644 --- a/CarpetX/src/driver.cxx +++ b/CarpetX/src/driver.cxx @@ -1332,13 +1332,13 @@ GHExt::PatchData::LevelData::LevelData(const int patch, const int level, cGH *const cctkGH = ghext->get_level_cctkGH(level); enter_patch_mode(cctkGH, leveldata); patch_cctkGH = GHExt::cctkGHptr(copy_cctkGH(cctkGH)); - int block = 0; + int component = 0; const auto mfitinfo = amrex::MFItInfo().EnableTiling(); for (amrex::MFIter mfi(*leveldata.fab, mfitinfo); mfi.isValid(); - ++mfi, ++block) { + ++mfi, ++component) { const MFPointer mfp(mfi); enter_local_mode(cctkGH, leveldata, mfp); - assert(int(local_cctkGHs.size()) == block); + assert(int(local_cctkGHs.size()) == component); local_cctkGHs.emplace_back(copy_cctkGH(cctkGH)); leave_local_mode(cctkGH, leveldata, mfp); } @@ -1494,7 +1494,7 @@ void GHExt::PatchData::LevelData::GroupData::apply_boundary_conditions( if (geom.isPeriodic(d)) gdomain.grow(d, mfab.nGrow(d)); - loop_over_blocks(mfab, [&](int index, int block) { + loop_over_components(mfab, [&](int index, int component) { amrex::FArrayBox &dest = mfab[index]; // const amrex::Box &box = mfp.fabbox(); // const amrex::Box &box = dest.box(); diff --git a/CarpetX/src/driver.hxx b/CarpetX/src/driver.hxx index 46e10fbdc..65cbfceab 100644 --- a/CarpetX/src/driver.hxx +++ b/CarpetX/src/driver.hxx @@ -239,11 +239,11 @@ struct GHExt { unique_ptr fab; cctkGHptr patch_cctkGH; - vector local_cctkGHs; // [block] + vector local_cctkGHs; // [component] cGH *get_patch_cctkGH() const { return patch_cctkGH.get(); } - cGH *get_local_cctkGH(const int block) const { - return local_cctkGHs.at(block).get(); + cGH *get_local_cctkGH(const int component) const { + return local_cctkGHs.at(component).get(); } struct GroupData : public CommonGroupData { @@ -330,10 +330,10 @@ struct GHExt { return patchdata.at(patch).leveldata.at(level).patch_cctkGH.get(); } cGH *get_local_cctkGH(const int level, const int patch, - const int block) const { + const int component) const { return patchdata.at(patch) .leveldata.at(level) - .local_cctkGHs.at(block) + .local_cctkGHs.at(component) .get(); } diff --git a/CarpetX/src/interpolate.cxx b/CarpetX/src/interpolate.cxx index 521140c22..ba04f1d18 100644 --- a/CarpetX/src/interpolate.cxx +++ b/CarpetX/src/interpolate.cxx @@ -20,7 +20,9 @@ #include #include #include +#include #include +#include #include #include @@ -31,11 +33,22 @@ namespace { // Interpolate a grid function at one point, dimensionally recursive template struct interpolator { - const amrex::Array4 &vars; + // TODO: `centering` is not used; remove it + // TODO: remove `dx` in favour of `x1`? + + static constexpr vect indextype{(centering & 0b100) != 0, + (centering & 0b010) != 0, + (centering & 0b001) != 0}; + + const GridDescBase &grid; const int vi; + const amrex::Array4 &vars; const vect &derivs; - const vect &x0; - const vect &dx; + + static constexpr T eps() { + using std::pow; + return pow(std::numeric_limits::epsilon(), T(3) / 4); + } // TODO: Check whether interpolated variables are valid @@ -43,12 +56,28 @@ template struct interpolator { template std::enable_if_t<(dir == -1), T> interpolate(const vect &i, const vect &di) const { + const amrex::IntVect j(i[0] + vars.begin.x, i[1] + vars.begin.y, + i[2] + vars.begin.z); #ifdef CCTK_DEBUG - assert(vars.contains(i[0], i[1], i[2])); + assert(vars.contains(j[0], j[1], j[2])); #endif - const auto val = vars(i[0], i[1], i[2], vi); + const T val = vars(j, vi); #ifdef CCTK_DEBUG using std::isfinite; + if (!(isfinite(val))) { + std::cerr << "!isfinite i=" << i << " di=" << di << " val=" << val + << "\n"; + for (int c = -1; c <= +1; ++c) + for (int b = -1; b <= +1; ++b) + for (int a = -1; a <= +1; ++a) + if (vars.contains(j[0] + a, j[1] + b, j[2] + c)) + // std::cerr << " val[" << a << "," << b << "," << c + // << "]=" << vars(j[0] + a, j[1] + b, j[2] + c, vi) << + // "\n"; + std::fprintf(stderr, "val[%d,%d,%d]=%.17g\n", j[0] + a, j[1] + b, + j[2] + c, + double(vars(j[0] + a, j[1] + b, j[2] + c, vi))); + } assert(isfinite(val)); #endif return val; @@ -58,17 +87,20 @@ template struct interpolator { template std::enable_if_t<(dir >= 0), T> interpolate(const vect &i, const vect &di) const { + static_assert(dir < dim); const auto DI = vect::unit(dir); // Ignore the centering for interpolation // switch ((centering >> (2 - dir)) & 1) + const T x = di[dir] - order / T(2); + // #ifdef CCTK_DEBUG + // using std::fabs; + // assert(fabs(x) <= T(0.5) + eps()); + // #endif + switch (order) { case 0: { -#ifdef CCTK_DEBUG - const T x = di[dir]; - assert(fabs(x) <= T(0.5)); -#endif const T y0 = interpolate(i, di); switch (derivs[dir]) { case 0: @@ -80,26 +112,18 @@ template struct interpolator { } } case 1: { - const T x = di[dir]; -#ifdef CCTK_DEBUG - assert(x >= T(0) && x <= T(1)); -#endif const T y0 = interpolate(i, di); const T y1 = interpolate(i + DI, di); switch (derivs[dir]) { case 0: - return (1 - x) * y0 + x * y1; + return (1 / T(2) - x) * y0 + (1 / T(2) + x) * y1; case 1: - return (-y0 + y1) / dx[dir]; + return (-y0 + y1) / grid.dx[dir]; case 2: return 0; } } case 2: { - const T x = di[dir] - order / T(2); -#ifdef CCTK_DEBUG - assert(fabs(x) <= T(0.5)); -#endif const T y0 = interpolate(i, di); const T y1 = interpolate(i + DI, di); const T y2 = interpolate(i + 2 * DI, di); @@ -110,16 +134,12 @@ template struct interpolator { (1 / T(2) * x + 1 / T(2) * pown(x, 2)) * y2; case 1: return ((-1 / T(2) + x) * y0 - 2 * x * y1 + (1 / T(2) + x) * y2) / - dx[dir]; + grid.dx[dir]; case 2: - return (y0 - 2 * y1 + y2) / pown(dx[dir], 2); + return (y0 - 2 * y1 + y2) / pown(grid.dx[dir], 2); } } case 3: { - const T x = di[dir] - order / T(2); -#ifdef CCTK_DEBUG - assert(fabs(x) <= T(0.5)); -#endif const T y0 = interpolate(i, di); const T y1 = interpolate(i + DI, di); const T y2 = interpolate(i + 2 * DI, di); @@ -143,18 +163,14 @@ template struct interpolator { (-9 / T(8) - 1 / T(2) * x + 3 / T(2) * pown(x, 2)) * y1 + (9 / T(8) - 1 / T(2) * x - 3 / T(2) * pown(x, 2)) * y2 + (-1 / T(24) + 1 / T(2) * x + 1 / T(2) * pown(x, 2)) * y3) / - dx[dir]; + grid.dx[dir]; case 2: return ((1 / T(2) - x) * y0 + (-1 / T(2) + 3 * x) * y1 + (-1 / T(2) - 3 * x) * y2 + (1 / T(2) + x) * y3) / - pown(dx[dir], 2); + pown(grid.dx[dir], 2); } } case 4: { - const T x = di[dir] - order / T(2); -#ifdef CCTK_DEBUG - assert(fabs(x) <= T(0.5)); -#endif const T y0 = interpolate(i, di); const T y1 = interpolate(i + DI, di); const T y2 = interpolate(i + 2 * DI, di); @@ -189,14 +205,14 @@ template struct interpolator { (-1 / T(12) - 1 / T(12) * x + 1 / T(4) * pown(x, 2) + 1 / T(6) * pown(x, 3)) * y4) / - dx[dir]; + grid.dx[dir]; case 2: return ((-1 / T(12) - 1 / T(2) * x + 1 / T(2) * pown(x, 2)) * y0 + (4 / T(3) + x - 2 * pown(x, 2)) * y1 + (-5 / T(2) + 3 * pown(x, 2)) * y2 + (4 / T(3) - x - 2 * pown(x, 2)) * y3 + (-1 / T(12) + 1 / T(2) * x + 1 / T(2) * pown(x, 2)) * y4) / - pown(dx[dir], 2); + pown(grid.dx[dir], 2); } } default: @@ -207,21 +223,116 @@ template struct interpolator { template void interpolate3d(const Particles &particles, std::vector &varresult) const { + const auto x0 = grid.x0 + (2 * grid.lbnd - !indextype) * grid.dx / 2; + const auto x1 = x0 + (grid.lsh - 1 - indextype) * grid.dx; + const auto dx = grid.dx; + + assert(vars.end.x - vars.begin.x == grid.lsh[0]); + assert(vars.end.y - vars.begin.y == grid.lsh[1]); + assert(vars.end.z - vars.begin.z == grid.lsh[2]); + + // We assume that the input is synchronized, i.e. that all ghost + // zones are valid, but all outer boundaries are invalid. + // TODO: Take multipatch boundary directions into account; forbid + // only interpatch boundaries but allow true outer boundaries. + vect, 2> allowed_boundaries; + for (int f = 0; f < 2; ++f) + for (int d = 0; d < dim; ++d) + allowed_boundaries[f][d] = grid.bbox[f][d]; + + // The point must lie inside the domain. At outer boundaries the + // point may be in the boundary region, but at ghost boundaries + // the point cannot be in the ghost region. We define as "ghost" + // region here the interpolation stencil size which is `order / + // 2`. + const auto x0_allowed = + x0 + + (!allowed_boundaries[0] * grid.nghostzones + (order - 1) / T(2)) * dx - + eps(); + const auto x1_allowed = + x1 - + (!allowed_boundaries[1] * grid.nghostzones + (order - 1) / T(2)) * dx + + eps(); + + // The allowed index range is [i0, i1] + const auto i0_allowed = !allowed_boundaries[0] * grid.nghostzones; + const auto i1_allowed = + grid.lsh - (!allowed_boundaries[1] * grid.nghostzones + order); + + // The interior index range is [i0, i1] + const auto i0_interior = grid.nghostzones; + const auto i1_interior = grid.lsh - (grid.nghostzones + order); + const int np = int(varresult.size()); + #pragma omp simd for (int n = 0; n < np; ++n) { - // const vect imin{vars.begin.x, vars.begin.y, vars.begin.z}; - // const vect imax{vars.end.x, vars.end.y, vars.end.z}; - vect i; - vect di; - for (int d = 0; d < dim; ++d) { - const T x = particles[n].rdata(d); - const T ri = (x - x0[d]) / dx[d]; + const vect x{particles[n].rdata(0), particles[n].rdata(1), + particles[n].rdata(2)}; + // // Ensure the point is inside the domain + // if (!(all(x >= x0_allowed && x <= x1_allowed))) + // std::cerr << "grid=" << grid + // << " allowed_boundaries=" << allowed_boundaries + // << " x0_allowed=" << x0_allowed + // << " x1_allowed=" << x1_allowed << " x=" << x << "\n"; + // assert(all(x >= x0_allowed && x <= x1_allowed)); + + // Find stencil anchor + const auto qi = (x - x0) / dx; + const auto lrint1 = [](auto a) { using std::lrint; - i[d] = lrint(ri - (order / T(2))); - di[d] = ri - i[d]; + return int(lrint(a)); + }; + const auto i = fmap(lrint1, qi - order / T(2)); + const auto di = qi - i; + // Consistency check + assert(all(i >= 0 && i + order < grid.lsh)); + + // // Push point away from boundaries + // const auto clamp1 = [](auto i, auto i0, auto i1) { + // using std::clamp; + // return clamp(i, i0, i1); + // }; + // const auto j = fmap(clamp1, i, i0_allowed, i1_allowed-1); + // const auto dj = di + (j - i); + + // Push point away from boundaries if they are just a little outside + auto j = i; + auto dj = di; + for (int d = 0; d < dim; ++d) { + if (j[d] < i0_allowed[d] && dj[d] - order / T(2) >= +T(0.5) - eps()) { + j[d] += 1; + dj[d] -= 1; + } + if (j[d] >= i1_allowed[d] && dj[d] - order / T(2) <= -T(0.5) + eps()) { + j[d] -= 1; + dj[d] += 1; + } } - varresult[n] = interpolate(i, di); + // assert(all(j >= i0_allowed && j < i1_allowed)); + + // Push points into the interior if they are just a little outside + for (int d = 0; d < dim; ++d) { + if (j[d] < i0_interior[d] && dj[d] - order / T(2) >= +T(0.5) - eps()) { + j[d] += 1; + dj[d] -= 1; + } + if (j[d] >= i1_interior[d] && dj[d] - order / T(2) <= -T(0.5) + eps()) { + j[d] -= 1; + dj[d] += 1; + } + } + + // Avoid points on boundaries + const bool is_allowed = all(j >= i0_allowed && j < i1_allowed); + const bool is_inside = all(j >= i0_interior && j < i1_interior); + assert(is_inside <= is_allowed); // `P <= Q` is `P implies Q` + + const T res = !is_allowed ? -2 + : !is_inside ? -1 + : interpolate(j, dj); + + varresult[n] = res; } } }; @@ -284,8 +395,7 @@ extern "C" CCTK_INT CarpetX_DriverInterpolate( varinds.at(i) = input_array_variable_indices[varinds.at(i)]; } } else { - // TODO: actually output the error code - CCTK_ERROR("TableGetIntArray failed."); + CCTK_VERROR("TableGetIntArray failed with error code %d", n_elems); } std::vector operations; @@ -381,6 +491,8 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const CCTK_REAL *restrict const xmax = geom.ProbHi(); const CCTK_REAL *restrict const dx = geom.CellSize(); using std::clamp; + // Push the particle at least 1/2 grid spacing into the domain + // TODO: push by less because 1/2 is too much if there are many AMR levels posx[n] = clamp(localsx[n], xmin[0] + dx[0] / 2, xmax[0] - dx[0] / 2); posy[n] = clamp(localsy[n], xmin[1] + dx[1] / 2, xmax[1] - dx[1] / 2); posz[n] = clamp(localsz[n], xmin[2] + dx[2] / 2, xmax[2] - dx[2] / 2); @@ -390,7 +502,6 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, using Container = amrex::AmrParticleContainer<3, 2>; using ParticleTile = Container::ParticleTileType; std::vector containers(ghext->num_patches()); - std::vector particle_tiles(ghext->num_patches()); for (int patch = 0; patch < ghext->num_patches(); ++patch) { const auto &restrict patchdata = ghext->patchdata.at(patch); containers.at(patch) = Container(patchdata.amrcore.get()); @@ -398,63 +509,115 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const auto &restrict leveldata = patchdata.leveldata.at(level); const amrex::MFIter mfi(*leveldata.fab); assert(mfi.isValid()); - particle_tiles.at(patch) = &containers.at(patch).GetParticles( + ParticleTile *const particle_tile = &containers.at(patch).GetParticles( level)[make_pair(mfi.index(), mfi.LocalTileIndex())]; - } - // Set particle positions - { + // Set particle positions const int proc = amrex::ParallelDescriptor::MyProc(); for (int n = 0; n < npoints; ++n) { - const int patch = patches[n]; - amrex::Particle<3, 2> p; - p.id() = Container::ParticleType::NextID(); - p.cpu() = proc; - p.pos(0) = posx[n]; // AMReX distribution position - p.pos(1) = posy[n]; - p.pos(2) = posz[n]; - p.rdata(0) = localsx[n]; // actual particle coordinate - p.rdata(1) = localsy[n]; - p.rdata(2) = localsz[n]; - p.idata(0) = proc; // source process - p.idata(1) = n; // source index - particle_tiles.at(patch)->push_back(p); + // TODO: Loop over points only once + if (patches.at(n) == patch) { + amrex::Particle<3, 2> p; + p.id() = Container::ParticleType::NextID(); + p.cpu() = proc; + p.pos(0) = posx[n]; // AMReX distribution position + p.pos(1) = posy[n]; + p.pos(2) = posz[n]; + p.rdata(0) = localsx[n]; // actual particle coordinate + p.rdata(1) = localsy[n]; + p.rdata(2) = localsz[n]; + p.idata(0) = proc; // source process + p.idata(1) = n; // source index + particle_tile->push_back(p); + } } } // Send particles to interpolation points for (auto &container : containers) { + const int patch = int(&container - containers.data()); #ifdef CCTK_DEBUG std::size_t old_nparticles = 0; + std::set oldids; { const auto &levels = container.GetParticles(); for (const auto &level : levels) { - for (const auto &[grid_tile, particle_tile] : level) { - CCTK_VINFO("old_nparticles: %zu", particle_tile.size()); - old_nparticles += particle_tile.size(); + const int lev = int(&level - levels.data()); + for (amrex::ParConstIter<3, 2> pti(container, lev); pti.isValid(); + ++pti) { + const auto &particles = pti.GetArrayOfStructs(); + const int component = MFPointer(pti).index(); + CCTK_VINFO("patch %d level %d component %d old_nparticles: %zu", + patch, lev, component, particles.size()); + for (const auto &particle : particles) { + oldids.insert(particle.id()); + // CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] " + // "proc=%d n=%d", + // int(particle.id()), int(particle.cpu()), + // particle.pos(0), particle.pos(1), particle.pos(2), + // particle.rdata(0), particle.rdata(1), + // particle.rdata(2), particle.idata(0), + // particle.idata(1)); + } + old_nparticles += particles.size(); } } + const MPI_Comm comm = amrex::ParallelDescriptor::Communicator(); + MPI_Allreduce(MPI_IN_PLACE, &old_nparticles, 1, + mpi_datatype::value, MPI_SUM, comm); } #endif + container.Redistribute(); + #ifdef CCTK_DEBUG std::size_t new_nparticles = 0; + std::set newids; { const auto &levels = container.GetParticles(); for (const auto &level : levels) { - for (const auto &[grid_tile, particle_tile] : level) { - CCTK_VINFO("new_nparticles: %zu", particle_tile.size()); - new_nparticles += particle_tile.size(); + const int lev = int(&level - levels.data()); + for (amrex::ParConstIter<3, 2> pti(container, lev); pti.isValid(); + ++pti) { + const int component = MFPointer(pti).index(); + const auto &particles = pti.GetArrayOfStructs(); + CCTK_VINFO("patch %d level %d component %d new_nparticles: %zu", + patch, lev, component, particles.size()); + for (const auto &particle : particles) { + newids.insert(particle.id()); + // CCTK_VINFO(" id=%d proc=%d pos=[%g,%g,%g] locals=[%g,%g,%g] " + // "proc=%d n=%d", + // int(particle.id()), int(particle.cpu()), + // particle.pos(0), particle.pos(1), particle.pos(2), + // particle.rdata(0), particle.rdata(1), + // particle.rdata(2), particle.idata(0), + // particle.idata(1)); + } + new_nparticles += particles.size(); } } + const MPI_Comm comm = amrex::ParallelDescriptor::Communicator(); + MPI_Allreduce(MPI_IN_PLACE, &new_nparticles, 1, + mpi_datatype::value, MPI_SUM, comm); + } + if (new_nparticles != old_nparticles) { + for (const auto oldid : oldids) + if (!newids.count(oldid)) + CCTK_VINFO("old id %d not present in new ids", oldid); + for (const auto newid : newids) + if (!oldids.count(newid)) + CCTK_VINFO("new id %d not present in old ids", newid); + CCTK_VERROR( + "We lost interpolation points on patch %d. Before redistributing: " + "%zu particles, after redistributing: %zu particles", + patch, old_nparticles, new_nparticles); } - if (new_nparticles != old_nparticles) - CCTK_ERROR("We lost interpolation points"); #endif } // Define result variables - std::map > results; + const int nprocs = amrex::ParallelDescriptor::NProcs(); + std::vector > results(nprocs); // [nprocs] // Interpolate constexpr int tl = 0; @@ -477,34 +640,30 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, const int patch = patchdata.patch; for (const auto &leveldata : patchdata.leveldata) { const int level = leveldata.level; - // CCTK_VINFO("interpolating patch %d level %d", patch, level); // TODO: use OpenMP for (amrex::ParIter<3, 2> pti(containers.at(patch), level); pti.isValid(); ++pti) { - const amrex::Geometry &geom = - ghext->patchdata.at(patch).amrcore->Geom(level); - const vect x0 = {geom.ProbLo()[0], geom.ProbLo()[1], - geom.ProbLo()[2]}; - const vect dx = {geom.CellSize()[0], geom.CellSize()[1], - geom.CellSize()[2]}; + const MFPointer mfp(pti); + const GridDesc grid(leveldata, mfp); + const int component = mfp.index(); const int np = pti.numParticles(); const auto &particles = pti.GetArrayOfStructs(); + std::cerr << "patch=" << patch << " level=" << level + << " component=" << component << " npoints=" << np << "\n"; std::vector > varresults(nvars); // TODO: Don't re-calculate interpolation coefficients for each // variable for (int v = 0; v < nvars; ++v) { - // CCTK_VINFO("interpolating level %d, variable %d", level, v); const int gi = givis.at(v).gi; const int vi = givis.at(v).vi; const auto &restrict groupdata = *leveldata.groupdata.at(gi); - const auto x0group = - x0 + vect(groupdata.indextype) * dx / 2; const int centering = groupdata.indextype[0] * 0b100 + groupdata.indextype[1] * 0b010 + groupdata.indextype[2] * 0b001; + assert(all(groupdata.nghostzones == grid.nghostzones)); const amrex::Array4 &vars = groupdata.mfab.at(tl)->array(pti); vect derivs; @@ -526,32 +685,32 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, switch (interpolation_order) { case 0: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 1: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 2: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 3: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 4: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } @@ -569,32 +728,32 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, switch (interpolation_order) { case 0: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 1: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 2: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 3: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } case 4: { - const interpolator interp{vars, vi, derivs, - x0group, dx}; + const interpolator interp{grid, vi, vars, + derivs}; interp.interpolate3d(particles, varresult); break; } @@ -618,8 +777,6 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, for (int n = 0; n < np; ++n) { const int proc = particles[n].idata(0); const int id = particles[n].idata(1); - if (!results.count(proc)) - results[proc]; auto &result = results.at(proc); result.push_back(id); for (int v = 0; v < nvars; ++v) @@ -638,48 +795,45 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_, // Collect particles back // CCTK_VINFO("collecting results"); - const int nprocs = amrex::ParallelDescriptor::NProcs(); const MPI_Comm comm = amrex::ParallelDescriptor::Communicator(); const MPI_Datatype datatype = mpi_datatype::value; - std::vector sendcounts(nprocs, 0); - for (const auto &proc_result : results) { - const int p = proc_result.first; - const auto &result = proc_result.second; - sendcounts.at(p) = result.size(); - } + // int total_npoints; + // MPI_Allreduce(&npoints, &total_npoints, 1, MPI_INT, MPI_SUM, comm); + + std::vector sendcounts(nprocs); std::vector senddispls(nprocs); - int sendcount = 0; + int total_sendcount = 0; for (int p = 0; p < nprocs; ++p) { - senddispls.at(p) = sendcount; - sendcount += sendcounts.at(p); + const auto &result = results.at(p); + sendcounts.at(p) = result.size(); + senddispls.at(p) = total_sendcount; + total_sendcount += sendcounts.at(p); } - // for (int p = 0; p < nprocs; ++p) - // CCTK_VINFO("[%d] senddispl=%d sendcount=%d", p, senddispls.at(p), - // sendcounts.at(p)); - // CCTK_VINFO("sendcount=%d", sendcount); std::vector recvcounts(nprocs); MPI_Alltoall(sendcounts.data(), 1, MPI_INT, recvcounts.data(), 1, MPI_INT, comm); std::vector recvdispls(nprocs); - int recvcount = 0; + int total_recvcount = 0; for (int p = 0; p < nprocs; ++p) { - recvdispls.at(p) = recvcount; - recvcount += recvcounts.at(p); + recvdispls.at(p) = total_recvcount; + total_recvcount += recvcounts.at(p); } - // for (int p = 0; p < nprocs; ++p) - // CCTK_VINFO("[%d] recvdispl=%d recvcount=%d", p, recvdispls.at(p), - // recvcounts.at(p)); - // CCTK_VINFO("recvcount=%d", recvcount); - // If this fails then there might be particles out of bounds - assert(recvcount == (nvars + 1) * npoints); - std::vector sendbuf(sendcount); - for (const auto &proc_result : results) { - const int p = proc_result.first; - const auto &result = proc_result.second; - copy(result.begin(), result.end(), &sendbuf.at(senddispls.at(p))); + + std::vector sendbuf(total_sendcount); + for (int p = 0; p < nprocs; ++p) { + // TODO: Don't copy, store data here right away + assert(p >= 0); + assert(p < int(results.size())); + const auto &result = results.at(p); + assert(sendcounts.at(p) == result.size()); + assert(p >= 0); + assert(p < int(senddispls.size())); + assert(senddispls.at(p) >= 0); + assert(senddispls.at(p) + sendcounts.at(p) <= int(sendbuf.size())); + std::copy(result.begin(), result.end(), sendbuf.data() + senddispls.at(p)); } - std::vector recvbuf(recvcount); + std::vector recvbuf(total_recvcount); MPI_Alltoallv(sendbuf.data(), sendcounts.data(), senddispls.data(), datatype, recvbuf.data(), recvcounts.data(), recvdispls.data(), datatype, comm); diff --git a/CarpetX/src/schedule.cxx b/CarpetX/src/schedule.cxx index ba177668d..79a623ecd 100644 --- a/CarpetX/src/schedule.cxx +++ b/CarpetX/src/schedule.cxx @@ -45,14 +45,29 @@ static inline int omp_in_parallel() { return 0; } namespace CarpetX { using namespace std; -#ifndef CCTK_HAVE_CGH_TILE +#ifndef CCTK_HAVE_CGH_LEVEL #error \ - "The Cactus flesh does not support cctk_tile_min etc. in the cGH structure. Update the flesh." + "The Cactus flesh does not support cctk_level in the cGH structure. Update the flesh." #endif #ifndef CCTK_HAVE_CGH_PATCH #error \ "The Cactus flesh does not support cctk_patch in the cGH structure. Update the flesh." #endif +#ifndef CCTK_HAVE_CGH_COMPONENT +#error \ + "The Cactus flesh does not support cctk_component in the cGH structure. Update the flesh." +#endif +#ifndef CCTK_HAVE_CGH_TILE +#error \ + "The Cactus flesh does not support cctk_tile_min etc. in the cGH structure. Update the flesh." +#endif + +#if defined _OPENMP +#if !defined AMREX_USE_OMP +#error \ + "Cactus is configured with OpenMP, but AMReX is configured without OpenMP. This does not work." +#endif +#endif namespace { double gettime() { @@ -82,20 +97,29 @@ int GroupStorageCrease(const cGH *cctkGH, int n_groups, const int *groups, GridDesc::GridDesc(const GHExt::PatchData::LevelData &leveldata, const MFPointer &mfp) { - // DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_PARAMETERS; + + // The number of ghostzones in each direction + // for (int d = 0; d < dim; ++d) + // nghostzones[d] = mfp.nGrowVect()[d]; + nghostzones = {ghost_size >= 0 ? ghost_size : ghost_size_x, + ghost_size >= 0 ? ghost_size : ghost_size_y, + ghost_size >= 0 ? ghost_size : ghost_size_z}; const auto &patchdata = ghext->patchdata.at(leveldata.patch); - const amrex::Box &fbx = mfp.fabbox(); // allocated array - const amrex::Box &vbx = mfp.validbox(); // interior region (without ghosts) - const amrex::Box &gbx = mfp.growntilebox(); // current region (with ghosts) + const amrex::IntVect ng(nghostzones[0], nghostzones[1], nghostzones[2]); const amrex::Box &domain = patchdata.amrcore->Geom(leveldata.level).Domain(); + const amrex::Box &vbx = mfp.validbox(); // interior region (without ghosts) + const amrex::Box &fbx = mfp.fabbox(ng); // allocated array + const amrex::Box &gbx = mfp.growntilebox(ng); // current region (with ghosts) for (int d = 0; d < dim; ++d) assert(domain.type(d) == amrex::IndexType::CELL); - // The number of ghostzones in each direction - for (int d = 0; d < dim; ++d) - nghostzones[d] = mfp.nGrowVect()[d]; + // Level, patch, and component + level = leveldata.level; + patch = leveldata.patch; + component = mfp.index(); // Global shape for (int d = 0; d < dim; ++d) @@ -177,22 +201,28 @@ GridDesc::GridDesc(const GHExt::PatchData::LevelData &leveldata, } GridDesc::GridDesc(const GHExt::PatchData::LevelData &leveldata, - const int block) { - // `global_block` is the global block index. + const int component) { + // `global_component` is the global component index. // There is no tiling. const auto &patchdata = ghext->patchdata.at(leveldata.patch); const amrex::FabArrayBase &fab = *leveldata.fab; - const amrex::Box &fbx = fab.fabbox(block); // allocated array - const amrex::Box &vbx = fab.box(block); // interior region (without ghosts) - const amrex::Box &gbx = fbx; // current region (with ghosts) + const amrex::Box &fbx = fab.fabbox(component); // allocated array + const amrex::Box &vbx = + fab.box(component); // interior region (without ghosts) + const amrex::Box &gbx = fbx; // current region (with ghosts) const amrex::Box &domain = patchdata.amrcore->Geom(leveldata.level).Domain(); for (int d = 0; d < dim; ++d) assert(domain.type(d) == amrex::IndexType::CELL); + // Level, patch, and component + level = leveldata.level; + patch = leveldata.patch; + this->component = component; + // The number of ghostzones in each direction for (int d = 0; d < dim; ++d) nghostzones[d] = fab.nGrowVect()[d]; @@ -281,7 +311,8 @@ GridDesc::GridDesc(const GHExt::PatchData::LevelData &leveldata, GridPtrDesc::GridPtrDesc(const GHExt::PatchData::LevelData &leveldata, const MFPointer &mfp) : GridDesc(leveldata, mfp) { - const amrex::Box &fbx = mfp.fabbox(); // allocated array + const amrex::IntVect ng(nghostzones[0], nghostzones[1], nghostzones[2]); + const amrex::Box &fbx = mfp.fabbox(ng); // allocated array cactus_offset = lbound(fbx); } @@ -290,7 +321,9 @@ GridPtrDesc1::GridPtrDesc1( const GHExt::PatchData::LevelData::GroupData &groupdata, const MFPointer &mfp) : GridDesc(leveldata, mfp) { - const amrex::Box &fbx = mfp.fabbox(); // allocated array + DECLARE_CCTK_PARAMETERS; + const amrex::IntVect ng(nghostzones[0], nghostzones[1], nghostzones[2]); + const amrex::Box &fbx = mfp.fabbox(ng); // allocated array cactus_offset = lbound(fbx); for (int d = 0; d < dim; ++d) { assert(groupdata.nghostzones.at(d) >= 0); @@ -377,9 +410,9 @@ void delete_cctkGH(cGH *cctkGH) { enum class mode_t { unknown, local, patch, level, global, meta }; mode_t current_mode(const cGH *restrict cctkGH) { - const bool have_local = cctkGH->cctk_lsh[0] != undefined; + const bool have_local = cctkGH->cctk_component != undefined; const bool have_patch = cctkGH->cctk_patch != undefined; - const bool have_level = cctkGH->cctk_levfac[0] != undefined; + const bool have_level = cctkGH->cctk_level != undefined; const bool have_global = cctkGH->cctk_nghostzones[0] != undefined; if (have_local && have_patch && have_level && have_global) return mode_t::local; @@ -447,8 +480,8 @@ void setup_cctkGH(cGH *restrict cctkGH) { cctkGH->cctk_delta_time = NAN; // init into meta mode - cctkGH->cctk_lsh[0] = undefined; - cctkGH->cctk_levfac[0] = undefined; + cctkGH->cctk_component = undefined; + cctkGH->cctk_level = undefined; cctkGH->cctk_patch = undefined; cctkGH->cctk_nghostzones[0] = undefined; assert(in_meta_mode(cctkGH)); @@ -541,6 +574,7 @@ void enter_level_mode(cGH *restrict cctkGH, const int level) { DECLARE_CCTK_PARAMETERS; assert(in_global_mode(cctkGH)); + cctkGH->cctk_level = level; for (int d = 0; d < dim; ++d) { // The refinement factor over the top level (coarsest) grid const int levfac = 1 << level; @@ -557,6 +591,7 @@ void enter_level_mode(cGH *restrict cctkGH, const int level) { } void leave_level_mode(cGH *restrict cctkGH, const int level) { assert(in_level_mode(cctkGH)); + cctkGH->cctk_level = undefined; for (int d = 0; d < dim; ++d) cctkGH->cctk_levfac[d] = undefined; for (int d = 0; d < dim; ++d) { @@ -617,6 +652,7 @@ void enter_local_mode(cGH *restrict cctkGH, assert(in_patch_mode(cctkGH)); const GridPtrDesc grid(leveldata, mfp); + cctkGH->cctk_component = mfp.index(); for (int d = 0; d < dim; ++d) cctkGH->cctk_lsh[d] = grid.lsh[d]; for (int d = 0; d < dim; ++d) @@ -685,6 +721,7 @@ void leave_local_mode(cGH *restrict cctkGH, const GHExt::PatchData::LevelData &restrict leveldata, const MFPointer &mfp) { assert(in_local_mode(cctkGH)); + cctkGH->cctk_component = undefined; for (int d = 0; d < dim; ++d) cctkGH->cctk_lsh[d] = undefined; for (int d = 0; d < dim; ++d) @@ -731,9 +768,9 @@ extern "C" CCTK_INT CarpetX_GetCallFunctionCount() { #endif } -void loop_over_blocks( +void loop_over_components( amrex::FabArrayBase &fab, - const std::function &block_kernel) { + const std::function &component_kernel) { DECLARE_CCTK_PARAMETERS; // Choose kernel launch method @@ -761,11 +798,11 @@ void loop_over_blocks( // No parallelism // Note: The amrex::MFIter uses global variables and OpenMP barriers - int block = 0; + int component = 0; const auto mfitinfo = amrex::MFItInfo().EnableTiling(); - for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++block) { + for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++component) { const MFPointer mfp(mfi); - block_kernel(mfp.index(), block); + component_kernel(mfp.index(), component); } break; } @@ -776,12 +813,12 @@ void loop_over_blocks( std::vector > tasks; // Note: The amrex::MFIter uses global variables and OpenMP barriers - int block = 0; + int component = 0; const auto mfitinfo = amrex::MFItInfo().EnableTiling(); - for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++block) { + for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++component) { const MFPointer mfp(mfi); - auto task = [&block_kernel, mfp, block]() { - block_kernel(mfp.index(), block); + auto task = [&component_kernel, mfp, component]() { + component_kernel(mfp.index(), component); }; tasks.push_back(std::move(task)); } @@ -800,11 +837,11 @@ void loop_over_blocks( // CUDA // No OpenMP parallelization when using GPUs - int block = 0; + int component = 0; const auto mfitinfo = amrex::MFItInfo().DisableDeviceSync().EnableTiling(); - for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++block) { + for (amrex::MFIter mfi(fab, mfitinfo); mfi.isValid(); ++mfi, ++component) { const MFPointer mfp(mfi); - block_kernel(mfp.index(), block); + component_kernel(mfp.index(), component); #ifdef AMREX_USE_GPU if (gpu_sync_after_every_kernel) { amrex::Gpu::streamSynchronize(); @@ -821,18 +858,20 @@ void loop_over_blocks( } } -void loop_over_blocks( +void loop_over_components( const active_levels_t &active_levels, - const std::function &block_kernel) { + const std::function &component_kernel) { DECLARE_CCTK_PARAMETERS; active_levels.loop([&](const auto &restrict leveldata) { - loop_over_blocks(*leveldata.fab, [&leveldata, &block_kernel]( - const int index, const int block) { - cGH *restrict const localGH = leveldata.get_local_cctkGH(block); - block_kernel(leveldata.patch, leveldata.level, index, block, localGH); - }); + loop_over_components( + *leveldata.fab, + [&leveldata, &component_kernel](const int index, const int component) { + cGH *restrict const localGH = leveldata.get_local_cctkGH(component); + component_kernel(leveldata.patch, leveldata.level, index, component, + localGH); + }); }); } @@ -1133,9 +1172,9 @@ int Initialise(tFleshConfig *config) { CCTK_REAL x0[dim], x1[dim], dx[dim]; for (int d = 0; d < dim; ++d) { dx[d] = patchGH->cctk_delta_space[d]; - x0[d] = patchGH->cctk_origin_space[d] - - (1 - 2 * nghostzones[d]) * dx[d] / 2; - x1[d] = x0[d] + (gsh[d] - 1 - 2 * nghostzones[d]) * dx[d]; + x0[d] = patchGH->cctk_origin_space[d] + + (2 * nghostzones[d] - 1) * dx[d] / 2; + x1[d] = x0[d] + (gsh[d] - 2 * nghostzones[d] - 1) * dx[d]; } #pragma omp critical { @@ -1975,8 +2014,9 @@ int CallFunction(void *function, cFunctionData *restrict attribute, switch (mode) { case mode_t::local: // Call function once per tile - loop_over_blocks(*active_levels, [&](int patch, int level, int index, - int block, const cGH *local_cctkGH) { + loop_over_components(*active_levels, [&](int patch, int level, int index, + int component, + const cGH *local_cctkGH) { update_cctkGH(const_cast(local_cctkGH), cctkGH); CCTK_CallFunction(function, attribute, const_cast(local_cctkGH)); }); @@ -2231,6 +2271,8 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups, CCTK_IsFunctionAliased("MultiPatch_Interpolate"); active_levels->loop([&](auto &restrict leveldata) { +#warning "TODO" + CCTK_VINFO("sync level=%d patch=%d", leveldata.level, leveldata.patch); for (const int gi : groups) { auto &restrict groupdata = *leveldata.groupdata.at(gi); const nan_handling_t nan_handling = groupdata.do_checkpoint diff --git a/CarpetX/src/schedule.hxx b/CarpetX/src/schedule.hxx index 896db43ef..2f9ed49f1 100644 --- a/CarpetX/src/schedule.hxx +++ b/CarpetX/src/schedule.hxx @@ -108,10 +108,8 @@ extern optional active_levels; // Like an MFIter, but does not support iteration, instead it can be copied struct MFPointer { int m_index; - amrex::Box m_fabbox; - amrex::Box m_growntilebox; - amrex::Box m_validbox; - amrex::IntVect m_nGrowVect; + amrex::Box m_validbox; // interior of component + amrex::Box m_tilebox; // interior of tile MFPointer() = delete; MFPointer(const MFPointer &) = default; @@ -119,15 +117,31 @@ struct MFPointer { MFPointer &operator=(const MFPointer &) = default; MFPointer &operator=(MFPointer &&) = default; MFPointer(const amrex::MFIter &mfi) - : m_index((assert(mfi.isValid()), mfi.index())), m_fabbox(mfi.fabbox()), - m_growntilebox(mfi.growntilebox()), m_validbox(mfi.validbox()), - m_nGrowVect(mfi.theFabArrayBase().nGrowVect()) {} + : m_index((assert(mfi.isValid()), mfi.index())), + m_validbox(mfi.validbox()), m_tilebox(mfi.tilebox()) {} constexpr int index() const noexcept { return m_index; } - constexpr amrex::Box fabbox() const noexcept { return m_fabbox; } - constexpr amrex::Box growntilebox() const noexcept { return m_growntilebox; } + constexpr amrex::Box validbox() const noexcept { return m_validbox; } - constexpr amrex::IntVect nGrowVect() const noexcept { return m_nGrowVect; } + + constexpr amrex::Box tilebox() const noexcept { return m_tilebox; } + + const amrex::Box fabbox(const amrex::IntVect &ng) const noexcept { + return amrex::grow(validbox(), ng); + } + + const amrex::Box growntilebox(const amrex::IntVect &ng) const noexcept { + // return m_growntilebox; + amrex::Box bx = tilebox(); + const amrex::Box &vbx = validbox(); + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + if (bx.smallEnd(d) == vbx.smallEnd(d)) + bx.growLo(d, ng[d]); + if (bx.bigEnd(d) == vbx.bigEnd(d)) + bx.growHi(d, ng[d]); + } + return bx; + } }; struct GridDesc : GridDescBase { @@ -135,7 +149,7 @@ struct GridDesc : GridDescBase { GridDesc() = delete; GridDesc(const GHExt::PatchData::LevelData &leveldata, const MFPointer &mfp); GridDesc(const GHExt::PatchData::LevelData &leveldata, - const int global_block); + const int global_component); GridDesc(const cGH *cctkGH) : GridDescBase(cctkGH) {} }; @@ -227,15 +241,15 @@ void leave_local_mode(cGH *restrict cctkGH, const GHExt::PatchData::LevelData &restrict leveldata, const MFPointer &mfp); -// Loop over all blocks of a single patch and level -void loop_over_blocks( +// Loop over all components of a single patch and level +void loop_over_components( amrex::FabArrayBase &fab, - const std::function &block_kernel); -// Loop over all blocks of several patches and levels -void loop_over_blocks( + const std::function &component_kernel); +// Loop over all components of several patches and levels +void loop_over_components( const active_levels_t &active_levels, - const std::function &block_kernel); + const std::function &component_kernel); void synchronize(); // These functions are defined in valid.cxx. These prototypes should diff --git a/CarpetX/src/valid.cxx b/CarpetX/src/valid.cxx index 8fc815253..93c745c8e 100644 --- a/CarpetX/src/valid.cxx +++ b/CarpetX/src/valid.cxx @@ -180,8 +180,8 @@ void poison_invalid(const GHExt::PatchData::LevelData &leveldata, const active_levels_t active_levels(leveldata.level, leveldata.level + 1, leveldata.patch, leveldata.patch + 1); - loop_over_blocks(active_levels, [&](int patch, int level, int index, - int block, const cGH *cctkGH) { + loop_over_components(active_levels, [&](int patch, int level, int index, + int component, const cGH *cctkGH) { const Loop::GridDescBaseDevice grid(cctkGH); const Loop::GF3D2layout layout(cctkGH, groupdata.indextype); const Loop::GF3D2 gf( @@ -279,8 +279,8 @@ void check_valid(const GHExt::PatchData::LevelData &leveldata, const active_levels_t active_levels(leveldata.level, leveldata.level + 1, leveldata.patch, leveldata.patch + 1); - loop_over_blocks(active_levels, [&](int patch, int level, int index, - int block, const cGH *cctkGH) { + loop_over_components(active_levels, [&](int patch, int level, int index, + int component, const cGH *cctkGH) { const Loop::GridDescBaseDevice grid(cctkGH); const Loop::GF3D2layout layout(cctkGH, groupdata.indextype); const Loop::GF3D2 gf( @@ -326,49 +326,68 @@ void check_valid(const GHExt::PatchData::LevelData &leveldata, struct info_t { where_t where; + int patch, level, component; vect I; vect X; CCTK_REAL val; }; std::vector infos; - loop_over_blocks(active_levels, [&](int patch, int level, int index, - int block, const cGH *cctkGH) { + loop_over_components(active_levels, [&](int patch, int level, int index, + int component, + const cGH *cctkGH) { const Loop::GridDescBaseDevice grid(cctkGH); const Loop::GF3D2layout layout(cctkGH, groupdata.indextype); const Loop::GF3D2 gf( - layout, static_cast(CCTK_VarDataPtrI( + layout, static_cast(CCTK_VarDataPtrI( cctkGH, tl, groupdata.firstvarindex + vi))); +#pragma omp critical(CarpetX_check_valid) + std::cerr << "valid grid=" << grid << "\n"; if (valid.valid_int) - grid.loop_idx( - where_t::interior, groupdata.indextype, groupdata.nghostzones, - [&](const Loop::PointDesc &p) { - if (nan_check(grid, gf, p)) + grid.loop_idx(where_t::interior, groupdata.indextype, + groupdata.nghostzones, [&](const Loop::PointDesc &p) { + if (nan_check(grid, gf, p)) #pragma omp critical(CarpetX_check_valid) - infos.push_back(info_t{where_t::interior, p.I, p.X, gf(p.I)}); - }); + infos.push_back(info_t{where_t::interior, p.patch, + p.level, p.component, p.I, + p.X, gf(p.I)}); + }); if (valid.valid_outer) - grid.loop_idx( - where_t::boundary, groupdata.indextype, groupdata.nghostzones, - [&](const Loop::PointDesc &p) { - if (nan_check(grid, gf, p)) + grid.loop_idx(where_t::boundary, groupdata.indextype, + groupdata.nghostzones, [&](const Loop::PointDesc &p) { + if (nan_check(grid, gf, p)) #pragma omp critical(CarpetX_check_valid) - infos.push_back(info_t{where_t::interior, p.I, p.X, gf(p.I)}); - }); + infos.push_back(info_t{where_t::boundary, p.patch, + p.level, p.component, p.I, + p.X, gf(p.I)}); + }); if (valid.valid_ghosts) - grid.loop_idx( - where_t::ghosts, groupdata.indextype, groupdata.nghostzones, - [&](const Loop::PointDesc &p) { - if (nan_check(grid, gf, p)) + grid.loop_idx(where_t::ghosts, groupdata.indextype, + groupdata.nghostzones, [&](const Loop::PointDesc &p) { + if (nan_check(grid, gf, p)) #pragma omp critical(CarpetX_check_valid) - infos.push_back(info_t{where_t::interior, p.I, p.X, gf(p.I)}); - }); + infos.push_back(info_t{where_t::ghosts, p.patch, + p.level, p.component, p.I, + p.X, gf(p.I)}); + }); }); synchronize(); std::sort(infos.begin(), infos.end(), [](const info_t &a, const info_t &b) { + if (a.level < b.level) + return true; + if (a.level > b.level) + return false; + if (a.patch < b.patch) + return true; + if (a.patch > b.patch) + return false; + if (a.component < b.component) + return true; + if (a.component > b.component) + return false; const std::less > lt; return lt(reversed(a.I), reversed(b.I)); }); @@ -377,7 +396,9 @@ void check_valid(const GHExt::PatchData::LevelData &leveldata, buf << setprecision(std::numeric_limits::digits10 + 1); for (const auto &info : infos) buf << "\n" - << info.where << " " << info.I << " " << info.X << " " << info.val; + << info.where << " level " << info.level << " patch " << info.patch + << " component " << info.component << " " << info.I << " " << info.X + << " " << info.val; CCTK_WARN(CCTK_WARN_ALERT, buf.str().c_str()); CCTK_VERROR("%s: Grid function \"%s\" contains nans, infinities, or " diff --git a/Loop/src/loop.cxx b/Loop/src/loop.cxx index cb119c665..f57aad2bf 100644 --- a/Loop/src/loop.cxx +++ b/Loop/src/loop.cxx @@ -23,7 +23,11 @@ std::ostream &operator<<(std::ostream &os, const where_t where) { std::ostream &operator<<(std::ostream &os, const PointDesc &p) { return os << "PointDesc{" + << "level:" << p.level << ", " + << "patch:" << p.patch << ", " + << "component:" << p.component << ", " << "I:" << p.I << ", " + << "iter:" << p.iter << ", " << "NI:" << p.NI << ", " << "I0:" << p.I0 << ", " << "BI:" << p.BI << ", " @@ -35,6 +39,10 @@ std::ostream &operator<<(std::ostream &os, const PointDesc &p) { GridDescBase::GridDescBase() {} GridDescBase::GridDescBase(const cGH *restrict cctkGH) { + level = cctkGH->cctk_level; + patch = cctkGH->cctk_patch; + component = cctkGH->cctk_component; + for (int d = 0; d < dim; ++d) { assert(cctkGH->cctk_gsh[d] != undefined); gsh[d] = cctkGH->cctk_gsh[d]; @@ -81,11 +89,17 @@ GridDescBase::GridDescBase(const cGH *restrict cctkGH) { } std::ostream &operator<<(std::ostream &os, const GridDescBase &grid) { + // Convert to vertex-centred boundaries + const auto x0 = grid.x0 - grid.dx / 2; + const auto x1 = x0 + (grid.gsh - 1) * grid.dx; return os << "GridDescBase{" - << "gsh" << grid.gsh << ",lbnd" << grid.lbnd << ",ubnd" << grid.ubnd - << ",lsh" << grid.lsh << ",bbox" << grid.bbox << ",nghostzones" - << grid.nghostzones << ",tmin" << grid.tmin << ",tmax" << grid.tmax - << "}"; + << "level=" << grid.level << ",patch=" << grid.patch + << ",component=" << grid.component << ",gsh=" << grid.gsh + << ",lbnd=" << grid.lbnd << ",ubnd=" << grid.ubnd + << ",lsh=" << grid.lsh << ",ash=" << grid.ash + << ",bbox=" << grid.bbox << ",nghostzones=" << grid.nghostzones + << ",tmin=" << grid.tmin << ",tmax=" << grid.tmax + << ",x0[vc]=" << x0 << ",x1[vc]=" << x1 << ",dx=" << grid.dx << "}"; } } // namespace Loop diff --git a/Loop/src/loop.hxx b/Loop/src/loop.hxx index 453016b67..ed3d408e9 100644 --- a/Loop/src/loop.hxx +++ b/Loop/src/loop.hxx @@ -63,6 +63,7 @@ template struct units_t { struct PointDesc { units_t DI; // direction unit vectors + int level, patch, component; vect I; // grid point int iter; // iteration // outward boundary normal (if in outer boundary), else zero @@ -72,7 +73,7 @@ struct PointDesc { vect BI; // outer boundary points for this grid function (might be outside the current - // grid function block) + // grid function component) vect bnd_min, bnd_max; vect loop_min, loop_max; // loop shape @@ -92,20 +93,23 @@ struct PointDesc { PointDesc &operator=(PointDesc &&) = default; constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST - PointDesc(const vect &I, const int iter, const vect &NI, + PointDesc(const int level, const int patch, const int component, + const vect &I, const int iter, const vect &NI, const vect &I0, const vect &BI, const vect &bnd_min, const vect &bnd_max, const vect &loop_min, const vect &loop_max, const vect &X, const vect &DX) - : I(I), iter(iter), NI(NI), I0(I0), BI(BI), bnd_min(bnd_min), - bnd_max(bnd_max), loop_min(loop_min), loop_max(loop_max), X(X), DX(DX), - imin(loop_min[0]), imax(loop_max[0]), i(I[0]), j(I[1]), k(I[2]), - x(X[0]), y(X[1]), z(X[2]), dx(DX[0]), dy(DX[1]), dz(DX[2]) {} + : level(level), patch(patch), component(component), I(I), iter(iter), + NI(NI), I0(I0), BI(BI), bnd_min(bnd_min), bnd_max(bnd_max), + loop_min(loop_min), loop_max(loop_max), X(X), DX(DX), imin(loop_min[0]), + imax(loop_max[0]), i(I[0]), j(I[1]), k(I[2]), x(X[0]), y(X[1]), z(X[2]), + dx(DX[0]), dy(DX[1]), dz(DX[2]) {} friend std::ostream &operator<<(std::ostream &os, const PointDesc &p); }; struct GridDescBase { + int level, patch, component; vect gsh; vect lbnd, ubnd; vect lsh; @@ -115,8 +119,10 @@ struct GridDescBase { vect tmin, tmax; // for current level - vect x0; - vect dx; + // TODO: these are still cell centred, and so are cctk_origin_space and + // cctk_delta_space; fix this! + vect x0; // origin_space + vect dx; // delta_space friend std::ostream &operator<<(std::ostream &os, const GridDescBase &grid); @@ -140,8 +146,8 @@ public: const vect X = x0 + (lbnd + I - vect(!CI) / 2) * dx; const vect DX = dx; - return PointDesc(I, iter, NI, I0, BI, bnd_min, bnd_max, loop_min, loop_max, - X, DX); + return PointDesc(level, patch, component, I, iter, NI, I0, BI, bnd_min, + bnd_max, loop_min, loop_max, X, DX); } // Loop over a given box @@ -181,7 +187,8 @@ public: } } - // Box for outer boundaries (might be outside the current grid function block) + // Box for outer boundaries (might be outside the current grid function + // component) template void boundary_box(const vect &group_nghostzones, vect &restrict bnd_min, @@ -194,7 +201,7 @@ public: } // Box for all points and for interior (non-ghost) points in the current grid - // function block (not restricted to a single tile) + // function component (not restricted to a single tile) template void domain_boxes(const vect &group_nghostzones, vect &restrict all_min, diff --git a/PDESolvers/src/solve.cxx b/PDESolvers/src/solve.cxx index 9a149a733..05998ef6b 100644 --- a/PDESolvers/src/solve.cxx +++ b/PDESolvers/src/solve.cxx @@ -117,10 +117,10 @@ void define_point_type() { // Initialize point type everywhere, assuming there is no synchronization, // prolongation, or restriction - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, @@ -148,10 +148,10 @@ void define_point_type() { }); // Set indicator to level - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_ind( layout1, @@ -198,10 +198,10 @@ void define_point_type() { } // Check where the indicator changed; these are the restricted points - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, @@ -220,10 +220,10 @@ void define_point_type() { }); // Set indicator to index - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_ind( layout1, @@ -253,10 +253,10 @@ void define_point_type() { // Check where the indicator changed; these are the synchronized points. // If points are both restricted and synchronized, then we count them as // synchronized. - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, @@ -276,10 +276,10 @@ void define_point_type() { }); // Set indicator to level - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_ind( layout1, @@ -319,10 +319,10 @@ void define_point_type() { // Check where the indicator changed; these are the prolongated points. If // points are both restricted and prolongated, then we count them as // prolongated. - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, @@ -345,10 +345,10 @@ void define_point_type() { }); // Invalidate indicator - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const auto &patchdata = CarpetX::ghext->patchdata.at(patch); const auto &leveldata = patchdata.leveldata.at(level); leveldata.groupdata.at(gi_ind)->valid.at(tl).at(vi_ind).set_all( @@ -358,10 +358,10 @@ void define_point_type() { // Collect some statistics Arith::vect npoints{0, 0, 0, 0, 0, 0}; - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, static_cast( @@ -393,12 +393,12 @@ void define_point_type() { void enumerate_points( int &restrict npoints_local, int &restrict npoints_global, - std::vector > &restrict block_offsets, - std::vector > &restrict block_sizes, + std::vector > &restrict component_offsets, + std::vector > &restrict component_sizes, int &restrict npoints_prolongated_local, int &restrict npoints_prolongated_global, - std::vector > &restrict block_prolongated_offsets, - std::vector > &restrict block_prolongated_sizes, + std::vector > &restrict component_prolongated_offsets, + std::vector > &restrict component_prolongated_sizes, csr_t &Jp) { int myproc; MPI_Comm_rank(MPI_COMM_WORLD, &myproc); @@ -429,41 +429,42 @@ void enumerate_points( // const int vi_vcx = vn_vcx - v0_vcx; // assert(vi_vcx >= 0); - // Determine number of blocks per level + // Determine number of components per level assert(CarpetX::ghext->num_patches() == 1); const auto &patchdata = CarpetX::ghext->patchdata.at(0); std::vector level_sizes(patchdata.leveldata.size(), 0); - std::vector level_maxblocks(patchdata.leveldata.size(), -1); - CarpetX::loop_over_blocks(*CarpetX::active_levels, - [&](const int patch, const int level, - const int index, const int block, - const cGH *restrict const cctkGH) { + std::vector level_maxcomponents(patchdata.leveldata.size(), -1); + CarpetX::loop_over_components( + *CarpetX::active_levels, + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { #pragma omp atomic - ++level_sizes.at(level); - int &maxblock = level_maxblocks.at(level); - using std::max; + ++level_sizes.at(level); + int &maxcomponent = level_maxcomponents.at(level); + using std::max; #pragma omp critical - maxblock = max(maxblock, block); - }); + maxcomponent = max(maxcomponent, component); + }); // Allocate data structure - block_offsets.resize(patchdata.leveldata.size()); // process local - block_sizes.resize(patchdata.leveldata.size()); - block_prolongated_offsets.resize(patchdata.leveldata.size()); // process local - block_prolongated_sizes.resize(patchdata.leveldata.size()); + component_offsets.resize(patchdata.leveldata.size()); // process local + component_sizes.resize(patchdata.leveldata.size()); + component_prolongated_offsets.resize( + patchdata.leveldata.size()); // process local + component_prolongated_sizes.resize(patchdata.leveldata.size()); for (std::size_t level = 0; level < patchdata.leveldata.size(); ++level) { - assert(level_maxblocks.at(level) + 1 == level_sizes.at(level)); - block_offsets.at(level).resize(level_sizes.at(level)); - block_sizes.at(level).resize(level_sizes.at(level)); - block_prolongated_offsets.at(level).resize(level_sizes.at(level)); - block_prolongated_sizes.at(level).resize(level_sizes.at(level)); + assert(level_maxcomponents.at(level) + 1 == level_sizes.at(level)); + component_offsets.at(level).resize(level_sizes.at(level)); + component_sizes.at(level).resize(level_sizes.at(level)); + component_prolongated_offsets.at(level).resize(level_sizes.at(level)); + component_prolongated_sizes.at(level).resize(level_sizes.at(level)); } // Enumerate and count points - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, static_cast( @@ -489,20 +490,23 @@ void enumerate_points( assert(0); } }); - block_sizes.at(level).at(block) = npoints; - block_prolongated_sizes.at(level).at(block) = npoints_prolongated; + component_sizes.at(level).at(component) = npoints; + component_prolongated_sizes.at(level).at(component) = + npoints_prolongated; }); // Local exclusive prefix sum int npoints = 0; int npoints_prolongated = 0; - for (std::size_t level = 0; level < block_offsets.size(); ++level) { - for (std::size_t block = 0; block < block_offsets.at(level).size(); - ++block) { - block_offsets.at(level).at(block) = npoints; - npoints += block_sizes.at(level).at(block); - block_prolongated_offsets.at(level).at(block) = npoints_prolongated; - npoints_prolongated += block_prolongated_sizes.at(level).at(block); + for (std::size_t level = 0; level < component_offsets.size(); ++level) { + for (std::size_t component = 0; + component < component_offsets.at(level).size(); ++component) { + component_offsets.at(level).at(component) = npoints; + npoints += component_sizes.at(level).at(component); + component_prolongated_offsets.at(level).at(component) = + npoints_prolongated; + npoints_prolongated += + component_prolongated_sizes.at(level).at(component); } } npoints_local = npoints; @@ -521,10 +525,10 @@ void enumerate_points( MPI_INT, MPI_SUM, MPI_COMM_WORLD); // Set indices - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_idx( layout1, @@ -533,9 +537,10 @@ void enumerate_points( layout1, static_cast( CCTK_VarDataPtrI(cctkGH, tl, vn_pt))); const CarpetX::GridDescBase grid(cctkGH); - int idx = npoints_offset + block_offsets.at(level).at(block); - int idx_prolongated = npoints_prolongated_offset + - block_prolongated_offsets.at(level).at(block); + int idx = npoints_offset + component_offsets.at(level).at(component); + int idx_prolongated = + npoints_prolongated_offset + + component_prolongated_offsets.at(level).at(component); grid.loop_all<0, 0, 0>( grid.nghostzones, [&](const Loop::PointDesc &p) ARITH_INLINE { switch (int(gf_pt(p.I))) { @@ -558,26 +563,29 @@ void enumerate_points( assert(0); } }); - assert(idx == npoints_offset + block_offsets.at(level).at(block) + - block_sizes.at(level).at(block)); + assert(idx == npoints_offset + + component_offsets.at(level).at(component) + + component_sizes.at(level).at(component)); if (!(idx_prolongated == npoints_prolongated_offset + - block_prolongated_offsets.at(level).at(block) + - block_prolongated_sizes.at(level).at(block))) { + component_prolongated_offsets.at(level).at(component) + + component_prolongated_sizes.at(level).at(component))) { std::cout << "level=" << level << "\n"; - std::cout << "block=" << block << "\n"; + std::cout << "component=" << component << "\n"; std::cout << "idx_prolongated=" << idx_prolongated << "\n"; std::cout << "npoints_prolongated_offset=" << npoints_prolongated_offset << "\n"; - std::cout << "block_prolongated_offsets=" - << block_prolongated_offsets.at(level).at(block) << "\n"; - std::cout << "block_prolongated_sizes=" - << block_prolongated_sizes.at(level).at(block) << "\n"; + std::cout << "component_prolongated_offsets=" + << component_prolongated_offsets.at(level).at(component) + << "\n"; + std::cout << "component_prolongated_sizes=" + << component_prolongated_sizes.at(level).at(component) + << "\n"; } assert(idx_prolongated == npoints_prolongated_offset + - block_prolongated_offsets.at(level).at(block) + - block_prolongated_sizes.at(level).at(block)); + component_prolongated_offsets.at(level).at(component) + + component_prolongated_sizes.at(level).at(component)); const auto &leveldata = patchdata.leveldata.at(level); leveldata.groupdata.at(gi_idx)->valid.at(tl).at(vi_idx).set_all( CarpetX::make_valid_all(), @@ -624,10 +632,10 @@ void enumerate_points( } // Check that restriction and synchronization worked - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const Loop::GF3D2 gf_pt( layout1, static_cast( @@ -669,14 +677,14 @@ void enumerate_points( locations(patchdata.leveldata.size()); for (int level = 0; level < int(patchdata.leveldata.size()); ++level) locations.at(level).resize(level_sizes.at(level)); - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { // Level 0 has no prolongated points if (level == 0) - assert(block_prolongated_sizes.at(level).at(block) == 0); - if (block_prolongated_sizes.at(level).at(block) == 0) + assert(component_prolongated_sizes.at(level).at(component) == 0); + if (component_prolongated_sizes.at(level).at(component) == 0) return; const auto &mfab = *patchdata.leveldata.at(level).groupdata.at(gi_idx)->mfab.at(tl); @@ -716,11 +724,11 @@ void enumerate_points( } } }); - locations.at(level).at(block) = std::move(locs); + locations.at(level).at(component) = std::move(locs); }); - // Find the blocks where the prolongation sources live, and determine their - // Jacobian indices + // Find the components where the prolongation sources live, and determine + // their Jacobian indices std::vector > > indices( patchdata.leveldata.size()); // std::vector > > vcoordxs( @@ -728,17 +736,17 @@ void enumerate_points( for (int level = 0; level < int(patchdata.leveldata.size()); ++level) { indices.at(level).resize(level_sizes.at(level)); // vcoordxs.at(level).resize(level_sizes.at(level)); - for (int block = 0; block < level_sizes.at(level); ++block) { - indices.at(level).at(block).resize(locations.at(level).at(block).size(), - 0x80000000U); - // vcoordxs.at(level).at(block).resize( - // locations.at(level).at(block).size(), -1.0 / 0.0); + for (int component = 0; component < level_sizes.at(level); ++component) { + indices.at(level).at(component).resize( + locations.at(level).at(component).size(), 0x80000000U); + // vcoordxs.at(level).at(component).resize( + // locations.at(level).at(component).size(), -1.0 / 0.0); } } - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const auto &patchdata = CarpetX::ghext->patchdata.at(patch); // The finest level has no prolongation sources if (level == int(patchdata.leveldata.size()) - 1) @@ -766,10 +774,11 @@ void enumerate_points( for (int level1 = 0; level1 < int(patchdata.leveldata.size()); ++level1) { - for (int block1 = 0; block1 < level_sizes.at(level1); ++block1) { - const auto &locs = locations.at(level1).at(block1); - auto &idxs = indices.at(level1).at(block1); - // auto &vcxs = vcoordxs.at(level1).at(block1); + for (int component1 = 0; component1 < level_sizes.at(level1); + ++component1) { + const auto &locs = locations.at(level1).at(component1); + auto &idxs = indices.at(level1).at(component1); + // auto &vcxs = vcoordxs.at(level1).at(component1); for (int n = 0; n < int(locs.size()); ++n) { const auto &loc = locs.at(n); const auto &cL = std::get<0>(loc); @@ -795,10 +804,10 @@ void enumerate_points( }); // Note: Prolongation boundaries do not yet work with multiple processes for (int level = 0; level < int(patchdata.leveldata.size()); ++level) { - for (int block = 0; block < level_sizes.at(level); ++block) { - for (const auto &idx : indices.at(level).at(block)) + for (int component = 0; component < level_sizes.at(level); ++component) { + for (const auto &idx : indices.at(level).at(component)) assert(idx != int(0x80000000U)); - // for (const auto &vcx : vcoordxs.at(level).at(block)) + // for (const auto &vcx : vcoordxs.at(level).at(component)) // assert(vcx != -1.0 / 0.0); } } @@ -808,14 +817,14 @@ void enumerate_points( Jpvalss(patchdata.leveldata.size()); for (int level = 0; level < int(patchdata.leveldata.size()); ++level) Jpvalss.at(level).resize(level_sizes.at(level)); - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { // Level 0 has no prolongated points if (level == 0) - assert(block_prolongated_sizes.at(level).at(block) == 0); - if (block_prolongated_sizes.at(level).at(block) == 0) + assert(component_prolongated_sizes.at(level).at(component) == 0); + if (component_prolongated_sizes.at(level).at(component) == 0) return; const auto &mfab = *CarpetX::ghext->patchdata.at(patch) .leveldata.at(level) @@ -838,9 +847,9 @@ void enumerate_points( const CarpetX::GridDescBase grid(cctkGH); const auto offset = amrex_origin; - const auto &locs = locations.at(level).at(block); - const auto &idxs = indices.at(level).at(block); - // const auto &vcxs = vcoordxs.at(level).at(block); + const auto &locs = locations.at(level).at(component); + const auto &idxs = indices.at(level).at(component); + // const auto &vcxs = vcoordxs.at(level).at(component); assert(idxs.size() == locs.size()); // assert(vcxs.size() == locs.size()); int n = 0; @@ -886,7 +895,7 @@ void enumerate_points( } }); assert(n == int(locs.size())); - Jpvalss.at(level).at(block) = std::move(Jpvals); + Jpvalss.at(level).at(component) = std::move(Jpvals); }); // Convert Jpvals into sparse matrix @@ -894,9 +903,10 @@ void enumerate_points( } } -void copy_Cactus_to_PETSc(Vec vec, const std::vector &varinds, - const std::vector > &block_offsets, - const std::vector > &block_sizes) { +void copy_Cactus_to_PETSc( + Vec vec, const std::vector &varinds, + const std::vector > &component_offsets, + const std::vector > &component_sizes) { PetscErrorCode ierr; int myproc; @@ -934,10 +944,10 @@ void copy_Cactus_to_PETSc(Vec vec, const std::vector &varinds, CCTK_REAL *restrict const petsc_ptr = vec_ptr; // Copy Cactus vector to PETSc - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const CarpetX::GridDescBase grid(cctkGH); const auto &patchdata = CarpetX::ghext->patchdata.at(patch); @@ -956,14 +966,15 @@ void copy_Cactus_to_PETSc(Vec vec, const std::vector &varinds, gfs.emplace_back(layout1, static_cast( CCTK_VarDataPtrI(cctkGH, tl, varinds.at(n)))); - const int block_offset = block_offsets.at(level).at(block); + const int component_offset = component_offsets.at(level).at(component); int nelems = 0; grid.loop_all<0, 0, 0>( grid.nghostzones, [&](const Loop::PointDesc &p) ARITH_INLINE { switch (int(gf_pt(p.I))) { case int(point_type_t::intr): for (int n = 0; n < nvars; ++n) - petsc_ptr[nvars * block_offset + nelems++] = gfs.at(n)(p.I); + petsc_ptr[nvars * component_offset + nelems++] = + gfs.at(n)(p.I); break; case int(point_type_t::bdry): case int(point_type_t::rest): @@ -975,16 +986,17 @@ void copy_Cactus_to_PETSc(Vec vec, const std::vector &varinds, assert(0); } }); - assert(nelems == nvars * block_sizes.at(level).at(block)); + assert(nelems == nvars * component_sizes.at(level).at(component)); }); ierr = VecRestoreArray(vec, &vec_ptr); assert(!ierr); } -void copy_PETSc_to_Cactus(Vec vec, const std::vector &varinds, - const std::vector > &block_offsets, - const std::vector > &block_sizes) { +void copy_PETSc_to_Cactus( + Vec vec, const std::vector &varinds, + const std::vector > &component_offsets, + const std::vector > &component_sizes) { PetscErrorCode ierr; int myproc; @@ -1022,10 +1034,10 @@ void copy_PETSc_to_Cactus(Vec vec, const std::vector &varinds, const CCTK_REAL *restrict const petsc_ptr = vec_ptr; // Copy PETSc vector to Cactus - CarpetX::loop_over_blocks( + CarpetX::loop_over_components( *CarpetX::active_levels, - [&](const int patch, const int level, const int index, const int block, - const cGH *restrict const cctkGH) { + [&](const int patch, const int level, const int index, + const int component, const cGH *restrict const cctkGH) { const Loop::GF3D2layout layout1(cctkGH, indextype); const CarpetX::GridDescBase grid(cctkGH); const auto &patchdata = CarpetX::ghext->patchdata.at(patch); @@ -1038,14 +1050,15 @@ void copy_PETSc_to_Cactus(Vec vec, const std::vector &varinds, for (int n = 0; n < nvars; ++n) gfs.emplace_back(layout1, static_cast(CCTK_VarDataPtrI( cctkGH, tl, varinds.at(n)))); - const int block_offset = block_offsets.at(level).at(block); + const int component_offset = component_offsets.at(level).at(component); int nelems = 0; grid.loop_all<0, 0, 0>( grid.nghostzones, [&](const Loop::PointDesc &p) ARITH_INLINE { switch (int(gf_pt(p.I))) { case int(point_type_t::intr): for (int n = 0; n < nvars; ++n) - gfs.at(n)(p.I) = petsc_ptr[nvars * block_offset + nelems++]; + gfs.at(n)(p.I) = + petsc_ptr[nvars * component_offset + nelems++]; break; case int(point_type_t::bdry): case int(point_type_t::rest): @@ -1057,7 +1070,7 @@ void copy_PETSc_to_Cactus(Vec vec, const std::vector &varinds, assert(0); } }); - assert(nelems == nvars * block_sizes.at(level).at(block)); + assert(nelems == nvars * component_sizes.at(level).at(component)); for (int n = 0; n < nvars; ++n) leveldata.groupdata.at(gis.at(n))->valid.at(tl).at(vis.at(n)).set_all( CarpetX::make_valid_int(), @@ -1121,14 +1134,15 @@ extern "C" void PDESolvers_Solve(CCTK_ARGUMENTS) { int npoints_local, npoints_global; int npoints_prolongated_local, npoints_prolongated_global; - // block_offsets are process local - std::vector > block_offsets, block_sizes; - std::vector > block_prolongated_offsets, - block_prolongated_sizes; + // component_offsets are process local + std::vector > component_offsets, component_sizes; + std::vector > component_prolongated_offsets, + component_prolongated_sizes; csr_t Jp; - enumerate_points(npoints_local, npoints_global, block_offsets, block_sizes, - npoints_prolongated_local, npoints_prolongated_global, - block_prolongated_offsets, block_prolongated_sizes, Jp); + enumerate_points(npoints_local, npoints_global, component_offsets, + component_sizes, npoints_prolongated_local, + npoints_prolongated_global, component_prolongated_offsets, + component_prolongated_sizes, Jp); // TODO: fix this const std::vector solinds{CCTK_VarIndex("Poisson::sol")}; @@ -1157,9 +1171,9 @@ extern "C" void PDESolvers_Solve(CCTK_ARGUMENTS) { std::function evalf = [&](SNES snes, Vec x, Vec f) { - copy_PETSc_to_Cactus(x, solinds, block_offsets, block_sizes); + copy_PETSc_to_Cactus(x, solinds, component_offsets, component_sizes); CallScheduleGroup(cctkGH, "PDESolvers_Residual"); - copy_Cactus_to_PETSc(f, resinds, block_offsets, block_sizes); + copy_Cactus_to_PETSc(f, resinds, component_offsets, component_sizes); return 0; }; ierr = SNESSetFunction(snes, r, FormFunction, &evalf); @@ -1183,7 +1197,7 @@ extern "C" void PDESolvers_Solve(CCTK_ARGUMENTS) { std::function evalJ = [&](SNES snes, Vec x, Mat J, Mat B) { - copy_PETSc_to_Cactus(x, solinds, block_offsets, block_sizes); + copy_PETSc_to_Cactus(x, solinds, component_offsets, component_sizes); CallScheduleGroup(cctkGH, "PDESolvers_Jacobian"); jacobians->define_matrix(Jp, J); jacobians->clear(); @@ -1216,7 +1230,7 @@ extern "C" void PDESolvers_Solve(CCTK_ARGUMENTS) { Vec x; ierr = VecDuplicate(r, &x); assert(!ierr); - copy_Cactus_to_PETSc(x, solinds, block_offsets, block_sizes); + copy_Cactus_to_PETSc(x, solinds, component_offsets, component_sizes); // Solve @@ -1269,7 +1283,7 @@ extern "C" void PDESolvers_Solve(CCTK_ARGUMENTS) { // Extract solution - copy_PETSc_to_Cactus(x, solinds, block_offsets, block_sizes); + copy_PETSc_to_Cactus(x, solinds, component_offsets, component_sizes); CallScheduleGroup(cctkGH, "PDESolvers_Residual"); { const int nvars = resinds.size(); diff --git a/scripts/carpetx.th b/scripts/carpetx.th index d0018a494..adedd122c 100644 --- a/scripts/carpetx.th +++ b/scripts/carpetx.th @@ -723,3 +723,15 @@ CarpetX/TestProlongate CarpetX/TestSymmetries CarpetX/TmunuBaseX CarpetX/WaveToyX + + +# # CapyrX thorns +# !TARGET = $ARR +# !TYPE = git +# # !URL = https://github.com/lucass-carneiro/CapyrX +# !URL = https://github.com/eschnett/CapyrX +# !REPO_BRANCH = eschnett/interpolate +# !REPO_PATH = $2 +# !CHECKOUT = +# CapyrX/MultiPatch +# CapyrX/TestMultiPatch