diff --git a/Source/AMRInterpolator/AMRInterpolator.hpp b/Source/AMRInterpolator/AMRInterpolator.hpp index b8b63c97d..5e14f983c 100644 --- a/Source/AMRInterpolator/AMRInterpolator.hpp +++ b/Source/AMRInterpolator/AMRInterpolator.hpp @@ -12,15 +12,12 @@ // Chombo includes #include "AMRLevel.H" +#include "LoHiSide.H" // Our includes #include "BoundaryConditions.hpp" #include "GRAMR.hpp" -#include "InterpSource.hpp" -#include "InterpolationAlgorithm.hpp" -#include "InterpolationLayout.hpp" #include "InterpolationQuery.hpp" - #include "MPIContext.hpp" #include "UserVariables.hpp" @@ -31,6 +28,19 @@ template class AMRInterpolator { + struct InterpolationLayout + { + std::vector rank; + std::vector level_idx; + std::vector box_idx; + + InterpolationLayout(int num_points) + : rank(num_points, -1), level_idx(num_points, -1), + box_idx(num_points, -1) + { + } + }; + public: // constructor for backward compatibility // (adds an artificial BC with only periodic BC) @@ -56,15 +66,17 @@ template class AMRInterpolator void limit_num_levels(unsigned int num_levels); void interp(InterpolationQuery &query); const AMR &getAMR() const; - const std::array &get_coarsest_dx(); - const std::array &get_coarsest_origin(); + const std::array &get_coarsest_dx() const; + const std::array &get_coarsest_origin() const; + bool get_boundary_reflective(Side::LoHiSide a_side, int a_dir) const; + bool get_boundary_periodic(int a_dir) const; private: void computeLevelLayouts(); InterpolationLayout findBoxes(InterpolationQuery &query); void prepareMPI(InterpolationQuery &query, - const InterpolationLayout layout); + const InterpolationLayout &layout); void exchangeMPIQuery(); void calculateAnswers(InterpolationQuery &query); void exchangeMPIAnswer(); diff --git a/Source/AMRInterpolator/AMRInterpolator.impl.hpp b/Source/AMRInterpolator/AMRInterpolator.impl.hpp index c5a970b39..75db74614 100644 --- a/Source/AMRInterpolator/AMRInterpolator.impl.hpp +++ b/Source/AMRInterpolator/AMRInterpolator.impl.hpp @@ -3,6 +3,10 @@ * Please refer to LICENSE in GRChombo's root directory. */ +#if !defined(AMRINTERPOLATOR_HPP_) +#error "This file should only be included through AMRInterpolator.hpp" +#endif + #ifndef AMRINTERPOLATOR_IMPL_HPP_ #define AMRINTERPOLATOR_IMPL_HPP_ @@ -41,9 +45,7 @@ void AMRInterpolator::refresh(const bool a_fill_ghosts) { CH_TIME("AMRInterpolator::refresh"); - const Vector &levels = - const_cast(m_gr_amr).getAMRLevels(); - m_num_levels = levels.size(); + m_num_levels = const_cast(m_gr_amr).getAMRLevels().size(); m_mem_level.clear(); m_mem_box.clear(); @@ -73,17 +75,32 @@ const AMR &AMRInterpolator::getAMR() const template const std::array & -AMRInterpolator::get_coarsest_dx() +AMRInterpolator::get_coarsest_dx() const { return m_coarsest_dx; } template const std::array & -AMRInterpolator::get_coarsest_origin() +AMRInterpolator::get_coarsest_origin() const { return m_coarsest_origin; } +template +bool AMRInterpolator::get_boundary_reflective(Side::LoHiSide a_side, + int a_dir) const +{ + if (a_side == Side::Lo) + return m_lo_boundary_reflective[a_dir]; + else + return m_hi_boundary_reflective[a_dir]; +} +template +bool AMRInterpolator::get_boundary_periodic(int a_dir) const +{ + return m_bc_params.is_periodic[a_dir]; +} + template void AMRInterpolator::limit_num_levels(unsigned int num_levels) { @@ -239,10 +256,14 @@ void AMRInterpolator::computeLevelLayouts() if (m_verbosity >= 2) { _pout << " Level " << level_idx << "\t" - << "dx=(" << m_dx[level_idx][0] << "," << m_dx[level_idx][1] - << "," << m_dx[level_idx][2] << ")\t" - << "grid_origin=(" << m_origin[level_idx][0] << "," - << m_origin[level_idx][1] << "," << m_origin[level_idx][2] + << "dx=(" + << D_TERM(m_dx[level_idx][0], << "," << m_dx[level_idx][1], + << "," << m_dx[level_idx][2]) + << ")\t" + << "grid_origin=(" + << D_TERM(m_origin[level_idx][0], + << "," << m_origin[level_idx][1], + << "," << m_origin[level_idx][2]) << ")" << endl; } @@ -257,7 +278,7 @@ void AMRInterpolator::computeLevelLayouts() } template -InterpolationLayout +typename AMRInterpolator::InterpolationLayout AMRInterpolator::findBoxes(InterpolationQuery &query) { CH_TIME("AMRInterpolator::findBoxes"); @@ -293,7 +314,7 @@ AMRInterpolator::findBoxes(InterpolationQuery &query) // const AMRLevel& level = *levels[level_idx]; // const LevelData& level_data = dynamic_cast(level).getLevelData(); const DisjointBoxLayout& + // InterpSource<>&>(level).getLevelData(); const DisjointBoxLayout& // box_layout = level_data.disjointBoxLayout(); const Box& domain_box = // level.problemDomain().domainBox(); @@ -381,7 +402,7 @@ AMRInterpolator::findBoxes(InterpolationQuery &query) const AMRLevel &level = *levels[level_idx]; const LevelData &level_data = - dynamic_cast(level).getLevelData( + dynamic_cast &>(level).getLevelData( VariableType::evolution); const DisjointBoxLayout &box_layout = level_data.disjointBoxLayout(); const Box &domain_box = level.problemDomain().domainBox(); @@ -495,7 +516,7 @@ AMRInterpolator::findBoxes(InterpolationQuery &query) template void AMRInterpolator::prepareMPI(InterpolationQuery &query, - const InterpolationLayout layout) + const InterpolationLayout &layout) { CH_TIME("AMRInterpolator::prepareMPI"); @@ -631,7 +652,6 @@ void AMRInterpolator::calculateAnswers(InterpolationQuery &query) const int num_answers = m_mpi.totalAnswerCount(); std::array grid_coord; - IntVect nearest; for (int answer_idx = 0; answer_idx < num_answers; ++answer_idx) { @@ -639,7 +659,8 @@ void AMRInterpolator::calculateAnswers(InterpolationQuery &query) const int level_idx = m_answer_level[answer_idx]; const AMRLevel &level = *levels[level_idx]; - const InterpSource &source = dynamic_cast(level); + const InterpSource<> &source = + dynamic_cast &>(level); const LevelData *const evolution_level_data_ptr = &source.getLevelData(VariableType::evolution); const LevelData *diagnostics_level_data_ptr; @@ -657,10 +678,6 @@ void AMRInterpolator::calculateAnswers(InterpolationQuery &query) &diagnostics_level_data_ptr->disjointBoxLayout(); } - const Box &domain_box = level.problemDomain().domainBox(); - const IntVect &small_end = domain_box.smallEnd(); - const IntVect &big_end = domain_box.bigEnd(); - // Convert the LayoutIndex to DataIndex const DataIndex evolution_data_idx( evolution_box_layout_ptr->layoutIterator()[box_idx]); @@ -696,22 +713,6 @@ void AMRInterpolator::calculateAnswers(InterpolationQuery &query) << (box.bigEnd()[i] + 0.5) << "]"; MayDay::Abort(s.str().c_str()); } - - // point lies beyond the "small end" of the whole domain, but still - // within the boundary cell - if (grid_coord[i] < small_end[i] && - grid_coord[i] >= small_end[i] - 0.5) - nearest[i] = small_end[i]; - - // point lies beyond the "big end" of the whole domain, but still - // within the boundary cell - else if (grid_coord[i] > big_end[i] && - grid_coord[i] <= big_end[i] + 0.5) - nearest[i] = big_end[i]; - - // otherwise we round to nearest grid point - else - nearest[i] = (int)ceil(grid_coord[i] - 0.5); } if (m_verbosity >= 2) @@ -737,7 +738,7 @@ void AMRInterpolator::calculateAnswers(InterpolationQuery &query) typedef std::vector comps_t; comps_t &comps = deriv_it->second; - algo.setup(deriv, m_dx[level_idx], grid_coord, nearest); + algo.setup(deriv, grid_coord); for (typename comps_t::iterator it = comps.begin(); it != comps.end(); ++it) diff --git a/Source/AMRInterpolator/InterpSource.hpp b/Source/AMRInterpolator/InterpSource.hpp index 90c35a1c8..581b2b8fc 100644 --- a/Source/AMRInterpolator/InterpSource.hpp +++ b/Source/AMRInterpolator/InterpSource.hpp @@ -16,14 +16,21 @@ // Chombo namespace #include "UsingNamespace.H" -// Abstrace base class to get the FABs out of an AMRLevel -class InterpSource +// Abstrace base class to identify points in a grid structure +template class InterpSource { public: virtual const LevelData &getLevelData( const VariableType var_type = VariableType::evolution) const = 0; - virtual bool - contains(const std::array &point) const = 0; + virtual bool contains(const std::array &point) const = 0; + virtual double get_dx(int dir) const = 0; + virtual std::array get_dxs() const + { + std::array dxs; + for (int i = 0; i < N_DIMS; ++i) + dxs[i] = get_dx(i); + return dxs; + } }; #endif /* INTERPSOURCE_H_ */ diff --git a/Source/AMRInterpolator/Lagrange.hpp b/Source/AMRInterpolator/Lagrange.hpp index c54aa13a6..838cf06d0 100644 --- a/Source/AMRInterpolator/Lagrange.hpp +++ b/Source/AMRInterpolator/Lagrange.hpp @@ -9,9 +9,9 @@ #include "InterpSource.hpp" #include -template class Lagrange +template class Lagrange { - const InterpSource &m_source; + const InterpSource &m_source; bool m_verbosity; struct Stencil; @@ -41,10 +41,10 @@ template class Lagrange // Helper function to generate tensor product weights // Argument 'dim' is used for recursion over dimensions. pair, std::vector> - generateStencil(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, - const IntVect &nearest, int dim = CH_SPACEDIM - 1); + generateStencil(const std::array &deriv, + const std::array &dx, + const std::array &eval_index, + int dim = N_DIMS - 1); std::vector m_interp_points; std::vector m_interp_weights; @@ -56,13 +56,16 @@ template class Lagrange multiset m_interp_pos; public: - Lagrange(const InterpSource &source, bool verbosity = false); + Lagrange(const InterpSource &source, bool verbosity = false); - void setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, - const IntVect &nearest); - double interpData(const FArrayBox &fab, int comp); + // eval_index is in 'index' coordinates, not physical coordinates + void setup(const std::array &deriv, + const std::array &eval_index); + + // any class with a method: + // Real get(const IntVect &a_iv, int a_comp) const + template + double interpData(const GeneralArrayBox &fab, int comp = 0); const static string TAG; }; diff --git a/Source/AMRInterpolator/Lagrange.impl.hpp b/Source/AMRInterpolator/Lagrange.impl.hpp index 48bc7ec68..31cfe58ee 100644 --- a/Source/AMRInterpolator/Lagrange.impl.hpp +++ b/Source/AMRInterpolator/Lagrange.impl.hpp @@ -6,8 +6,8 @@ #ifndef LAGRANGE_IMPL_HPP_ #define LAGRANGE_IMPL_HPP_ -template -const string Lagrange::TAG = "\x1b[36;1m[Lagrange]\x1b[0m "; +template +const string Lagrange::TAG = "\x1b[36;1m[Lagrange]\x1b[0m "; /* Finite difference weight generation algorithm * @@ -20,9 +20,9 @@ const string Lagrange::TAG = "\x1b[36;1m[Lagrange]\x1b[0m "; * routine: - change type of c1,c2,c3 to 'double' - replace '0', 'i', 'j' in the * annotated lines with 'grid[0]', 'grid[i]', 'grid[j]'. */ -template -Lagrange::Stencil::Stencil(int width, int deriv, double dx, - double point_offset) +template +Lagrange::Stencil::Stencil(int width, int deriv, double dx, + double point_offset) : m_width(width), m_deriv(deriv), m_dx(dx), m_point_offset(point_offset) { int c1 = 1; @@ -102,17 +102,17 @@ Lagrange::Stencil::Stencil(int width, int deriv, double dx, } } -template -bool Lagrange::Stencil::operator==( - const Lagrange::Stencil &rhs) const +template +bool Lagrange::Stencil::operator==( + const Lagrange::Stencil &rhs) const { return (rhs.m_width == m_width) && (rhs.m_deriv == m_deriv) && (rhs.m_point_offset == m_point_offset) && (rhs.dx == m_dx); } -template -bool Lagrange::Stencil::isSameAs(int width, int deriv, double dx, - double point_offset) const +template +bool Lagrange::Stencil::isSameAs(int width, int deriv, double dx, + double point_offset) const { return (width == m_width) && (deriv == m_deriv) && (dx == m_dx) && (point_offset == m_point_offset); @@ -120,10 +120,10 @@ bool Lagrange::Stencil::isSameAs(int width, int deriv, double dx, /* STENCIL ACCESSOR */ -template -typename Lagrange::Stencil -Lagrange::getStencil(int width, int deriv, double dx, - double point_offset) +template +typename Lagrange::Stencil +Lagrange::getStencil(int width, int deriv, double dx, + double point_offset) { for (typename stencil_collection_t::iterator it = m_memoized_stencils.begin(); @@ -142,8 +142,8 @@ Lagrange::getStencil(int width, int deriv, double dx, Stencil(width, deriv, dx, point_offset)); } -template -const double &Lagrange::Stencil::operator[](unsigned int i) const +template +const double &Lagrange::Stencil::operator[](unsigned int i) const { CH_assert(i < m_width); return m_weights[i]; @@ -151,28 +151,29 @@ const double &Lagrange::Stencil::operator[](unsigned int i) const /* LAGRANGE TENSOR PRODUCT LOGIC */ -template -Lagrange::Lagrange(const InterpSource &source, bool verbosity) +template +Lagrange::Lagrange(const InterpSource &source, + bool verbosity) : m_source(source), m_verbosity(verbosity) { } -template -void Lagrange::setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, - const IntVect &nearest) +template +void Lagrange::setup( + const std::array &deriv, + const std::array &eval_index) { + std::array dxs = m_source.get_dxs(); pair, std::vector> result = - generateStencil(deriv, dx, evalCoord, nearest); + generateStencil(deriv, dxs, eval_index); m_interp_points = result.first; m_interp_weights = result.second; /* - pout() << TAG << "Stencil: coord = { "; - for (int i = 0; i < CH_SPACEDIM; ++i) + pout() << TAG << "Stencil: point = { "; + for (int i = 0; i < N_DIMS; ++i) { - pout() << evalCoord[i] << " "; + pout() << eval_index[i] << " "; } pout() << "}, weights = { "; for (int i = 0; i < m_interp_weights.size(); ++i) @@ -183,8 +184,9 @@ void Lagrange::setup(const std::array &deriv, */ } -template -double Lagrange::interpData(const FArrayBox &fab, int comp) +template +template +double Lagrange::interpData(const GeneralArrayBox &fab, int comp) { /* m_interp_neg.clear(); @@ -236,13 +238,11 @@ double Lagrange::interpData(const FArrayBox &fab, int comp) return accum; } -template +template pair, std::vector> -Lagrange::generateStencil( - const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, const IntVect &nearest, - int dim) +Lagrange::generateStencil( + const std::array &deriv, const std::array &dx, + const std::array &eval_index, int dim) { std::vector out_points; std::vector out_weights; @@ -266,16 +266,20 @@ Lagrange::generateStencil( int points_min = Order + deriv[dim]; int points_max = Order + deriv[dim]; - std::array interp_coord = evalCoord; - int candidate = nearest[dim]; - int grown_direction = (nearest[dim] - evalCoord[dim] < 0) ? DOWN : UP; + std::array interp_point = eval_index; + + // TF: assume nearest point is at 'std::round(eval_index[dim])' + // this used to be an argument, but I think this assumption should always + // hold + int candidate = std::round(eval_index[dim]); + int grown_direction = (candidate - eval_index[dim] < 0) ? DOWN : UP; while ((can_grow[DOWN] || can_grow[UP]) && (points_max - points_min < Order + deriv[dim])) { - interp_coord[dim] = candidate; + interp_point[dim] = candidate; - if (m_source.contains(interp_coord)) + if (m_source.contains(interp_point)) { int idx = (grown_direction == DOWN) ? (--points_min) : (points_max++); @@ -301,12 +305,12 @@ Lagrange::generateStencil( const Stencil my_weights = getStencil(stencil_width, deriv[dim], dx[dim], - evalCoord[dim] - my_points[points_min]); + eval_index[dim] - my_points[points_min]); if (m_verbosity) { pout() << TAG << "Stencil: dim = " << dim - << ", coord = " << evalCoord[dim] << ", points = { "; + << ", point = " << eval_index[dim] << ", points = { "; for (int i = points_min; i < points_max; ++i) { pout() << my_points[i] << " "; @@ -324,13 +328,13 @@ Lagrange::generateStencil( // first. for (int i = 0; i < stencil_width; ++i) { - interp_coord[dim] = my_points[i + points_min]; + interp_point[dim] = my_points[i + points_min]; if (dim > 0) { // Descend to the next dimension pair, std::vector> sub_result = - generateStencil(deriv, dx, interp_coord, nearest, dim - 1); + generateStencil(deriv, dx, interp_point, dim - 1); std::vector &sub_points = sub_result.first; std::vector &sub_weights = sub_result.second; @@ -350,8 +354,8 @@ Lagrange::generateStencil( if (my_weights[i] != 0) { out_points.push_back(IntVect(D_DECL6( - interp_coord[0], interp_coord[1], interp_coord[2], - interp_coord[3], interp_coord[4], interp_coord[5]))); + interp_point[0], interp_point[1], interp_point[2], + interp_point[3], interp_point[4], interp_point[5]))); out_weights.push_back(my_weights[i]); } } diff --git a/Source/AMRInterpolator/QuinticConvolution.hpp b/Source/AMRInterpolator/QuinticConvolution.hpp index c01a6337a..1e5cfcd95 100644 --- a/Source/AMRInterpolator/QuinticConvolution.hpp +++ b/Source/AMRInterpolator/QuinticConvolution.hpp @@ -9,22 +9,26 @@ #include "InterpSource.hpp" #include -class QuinticConvolution +template class QuinticConvolution { - const InterpSource &m_source; + const InterpSource &m_source; bool m_verbosity; std::vector m_interp_points; std::vector m_interp_weights; public: - QuinticConvolution(const InterpSource &source, bool verbosity = false); + QuinticConvolution(const InterpSource &source, + bool verbosity = false); - void setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, - const IntVect &nearest); - double interpData(const FArrayBox &fab, int comp); + // eval_index is in 'index' coordinates, not physical coordinates + void setup(const std::array &deriv, + const std::array &eval_index); + + // any class with a method: + // Real get(const IntVect &a_iv, int a_comp) const + template + double interpData(const GeneralArrayBox &fab, int comp = 0); const static string TAG; }; diff --git a/Source/AMRInterpolator/QuinticConvolution.impl.hpp b/Source/AMRInterpolator/QuinticConvolution.impl.hpp index 38c9783fa..877be4d31 100644 --- a/Source/AMRInterpolator/QuinticConvolution.impl.hpp +++ b/Source/AMRInterpolator/QuinticConvolution.impl.hpp @@ -6,28 +6,33 @@ #ifndef QUINTICCONVOLUTION_IMPL_HPP_ #define QUINTICCONVOLUTION_IMPL_HPP_ -const string QuinticConvolution::TAG = "\x1b[36;1m[QuinticConvolution]\x1b[0m "; +template +const string QuinticConvolution::TAG = + "\x1b[36;1m[QuinticConvolution]\x1b[0m "; -QuinticConvolution::QuinticConvolution(const InterpSource &source, - bool verbosity) +template +QuinticConvolution::QuinticConvolution( + const InterpSource &source, bool verbosity) : m_source(source), m_verbosity(verbosity) { - CH_assert(CH_SPACEDIM <= 3); + CH_assert(N_DIMS <= 3); } -void QuinticConvolution::setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord, - const IntVect &nearest) +template +void QuinticConvolution::setup( + const std::array &deriv, + const std::array &eval_index) { m_interp_points.clear(); m_interp_weights.clear(); - double weights_1d[CH_SPACEDIM][6]; + std::array dxs = m_source.get_dxs(); - for (int dim = 0; dim < CH_SPACEDIM; ++dim) + double weights_1d[N_DIMS][6]; + + for (int dim = 0; dim < N_DIMS; ++dim) { - double s = evalCoord[dim] - floor(evalCoord[dim]); + double s = eval_index[dim] - floor(eval_index[dim]); if (deriv[dim] == 0) { @@ -56,40 +61,40 @@ void QuinticConvolution::setup(const std::array &deriv, weights_1d[dim][0] = (0.046875 + s * (-0.375 + s * (0.84375 + (-0.75 + 0.234375 * s) * s))) / - dx[dim]; + dxs[dim]; weights_1d[dim][1] = (-0.59375 + s * (2.5 + s * (-1.6875 + s * (-1.1875 + 1.015625 * s)))) / - dx[dim]; + dxs[dim]; weights_1d[dim][2] = - (s * (-4.25 + (7.875 - 4.21875 * s) * (s * s))) / dx[dim]; + (s * (-4.25 + (7.875 - 4.21875 * s) * (s * s))) / dxs[dim]; weights_1d[dim][3] = (0.59375 + s * (2.5 + s * (1.6875 + s * (-9. + 4.21875 * s)))) / - dx[dim]; + dxs[dim]; weights_1d[dim][4] = (-0.046875 + s * (-0.375 + s * (-0.84375 + (2.875 - 1.015625 * s) * s))) / - dx[dim]; + dxs[dim]; weights_1d[dim][5] = - ((0.1875 - 0.234375 * s) * s * (s * s)) / dx[dim]; + ((0.1875 - 0.234375 * s) * s * (s * s)) / dxs[dim]; } else if (deriv[dim] == 2) { weights_1d[dim][0] = (-0.375 + s * (1.6875 + (-2.25 + 0.9375 * s) * s)) / - (dx[dim] * dx[dim]); + (dxs[dim] * dxs[dim]); weights_1d[dim][1] = (2.5 + s * (-3.375 + s * (-3.5625 + 4.0625 * s))) / - (dx[dim] * dx[dim]); - weights_1d[dim][2] = - (-4.25 + (23.625 - 16.875 * s) * (s * s)) / (dx[dim] * dx[dim]); + (dxs[dim] * dxs[dim]); + weights_1d[dim][2] = (-4.25 + (23.625 - 16.875 * s) * (s * s)) / + (dxs[dim] * dxs[dim]); weights_1d[dim][3] = (2.5 + s * (3.375 + s * (-27. + 16.875 * s))) / - (dx[dim] * dx[dim]); + (dxs[dim] * dxs[dim]); weights_1d[dim][4] = (-0.375 + s * (-1.6875 + (8.625 - 4.0625 * s) * s)) / - (dx[dim] * dx[dim]); + (dxs[dim] * dxs[dim]); weights_1d[dim][5] = - ((0.5625 - 0.9375 * s) * (s * s)) / (dx[dim] * dx[dim]); + ((0.5625 - 0.9375 * s) * (s * s)) / (dxs[dim] * dxs[dim]); } else { @@ -98,36 +103,48 @@ void QuinticConvolution::setup(const std::array &deriv, } } - std::array interp_coord; + std::array interp_index; -#if CH_SPACEDIM >= 3 +#if N_DIMS >= 3 for (int z = 0; z < 6; ++z) { - interp_coord[2] = floor(evalCoord[2]) + z - 2; + interp_index[2] = floor(eval_index[2]) + z - 2; #endif +#if N_DIMS >= 2 for (int y = 0; y < 6; ++y) { - interp_coord[1] = floor(evalCoord[1]) + y - 2; + interp_index[1] = floor(eval_index[1]) + y - 2; +#endif for (int x = 0; x < 6; ++x) { - interp_coord[0] = floor(evalCoord[0]) + x - 2; - CH_assert(m_source.contains(interp_coord)); + interp_index[0] = floor(eval_index[0]) + x - 2; + CH_assert(m_source.contains(interp_index)); m_interp_points.push_back(IntVect(D_DECL6( - interp_coord[0], interp_coord[1], interp_coord[2], - interp_coord[3], interp_coord[4], interp_coord[5]))); - m_interp_weights.push_back(D_TERM6(weights_1d[0][x], - *weights_1d[1][y], - *weights_1d[2][z], , , )); + interp_index[0], interp_index[1], interp_index[2], + interp_index[3], interp_index[4], interp_index[5]))); +#if N_DIMS >= 3 + m_interp_weights.push_back(weights_1d[0][x] * weights_1d[1][y] * + weights_1d[2][z]); +#elif N_DIMS >= 2 + m_interp_weights.push_back(weights_1d[0][x] * weights_1d[1][y]); +#else + m_interp_weights.push_back(weights_1d[0][x]); +#endif } +#if N_DIMS >= 2 } -#if CH_SPACEDIM >= 3 +#endif +#if N_DIMS >= 3 } #endif } -double QuinticConvolution::interpData(const FArrayBox &fab, int comp) +template +template +double QuinticConvolution::interpData(const GeneralArrayBox &fab, + int comp) { long double accum = 0.0; diff --git a/Source/AMRInterpolator/SimpleArrayBox.hpp b/Source/AMRInterpolator/SimpleArrayBox.hpp new file mode 100644 index 000000000..8fef95ec0 --- /dev/null +++ b/Source/AMRInterpolator/SimpleArrayBox.hpp @@ -0,0 +1,89 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMPLEARRAYBOX +#define SIMPLEARRAYBOX + +// Chombo includes +#include "CH_assert.H" +#include "MayDay.H" + +// Chombo namespace +#include "UsingNamespace.H" + +// Global Array Box class to interpolate N_DIM dimensional data +// where 'm_f' is a flat array over the dimensions +// Created to replace an FArrayBox in interpolation ::get methods (e.g. +// Lagrange) +template class SimpleArrayBox +{ + // coordinates (u,v,f) + std::array m_points_per_dir; + std::vector m_f; + std::array m_is_periodic; + + int apply_periodic(const IntVect &a_iv, int dir) const + { + int idx = a_iv[dir]; + if (idx < 0 || idx >= m_points_per_dir[dir]) + { + if (m_is_periodic[dir]) + { + while (idx < 0) + idx += m_points_per_dir[dir]; + while (idx >= m_points_per_dir[dir]) + idx -= m_points_per_dir[dir]; + } + else + { + pout() << "Direction " << dir << " invalid in IntVect " << a_iv + << std::endl; + MayDay::Error("[SimpleArrayBox] Trying to access index out of " + "valid range."); + } + } + return idx; + } + + public: + // IntVect will be CH_SPACEDIM dimension, but ignore the lasts components if + // N_DIMS < CH_SPACEDIM + // the 'comp' argument doesn't matter, always assume 'm_f' is what matters + Real get(const IntVect &a_iv, int a_comp) const + { + int global_idx = apply_periodic(a_iv, N_DIMS - 1); + + for (int i = N_DIMS - 2; i >= 0; --i) + global_idx = + global_idx * m_points_per_dir[i] + apply_periodic(a_iv, i); + + return m_f[global_idx]; + } + + SimpleArrayBox() {} + + SimpleArrayBox(std::array a_points_per_dir, + std::vector a_f, + std::array a_is_periodic = {false}) + { + set(a_points_per_dir, a_f, a_is_periodic); + } + + void set(std::array a_points_per_dir, std::vector a_f, + std::array a_is_periodic = {false}) + { + m_points_per_dir = a_points_per_dir; + m_f = a_f; + m_is_periodic = a_is_periodic; + + CH_assert(N_DIMS <= CH_SPACEDIM); + int total = 1.; + for (int i = 0; i < N_DIMS; ++i) + total *= m_points_per_dir[i]; + CH_assert(m_f.size() == total); + } +}; + +#endif /* SIMPLEARRAYBOX */ diff --git a/Source/AMRInterpolator/SimpleInterpSource.hpp b/Source/AMRInterpolator/SimpleInterpSource.hpp new file mode 100644 index 000000000..73a25f5ae --- /dev/null +++ b/Source/AMRInterpolator/SimpleInterpSource.hpp @@ -0,0 +1,62 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMPLEINTERPSOURCE_H_ +#define SIMPLEINTERPSOURCE_H_ + +#include "InterpSource.hpp" + +// Simple interpolation source where all grid points are in the same "Box" +// Relevant to be able to use Interpolation classes by themselves +template class SimpleInterpSource : public InterpSource +{ + std::array m_points_per_dir; + std::array m_is_periodic; + std::array m_dxs; + + LevelData *m_fake; + + public: + SimpleInterpSource() {} + SimpleInterpSource(std::array a_points_per_dir, + std::array a_dxs, + std::array a_is_periodic = {false}) + { + set(a_points_per_dir, a_dxs, a_is_periodic); + } + + void set(std::array a_points_per_dir, + std::array a_dxs, + std::array a_is_periodic = {false}) + { + m_points_per_dir = a_points_per_dir; + m_dxs = a_dxs; + m_is_periodic = a_is_periodic; + } + + const LevelData & + getLevelData(const VariableType var_type = VariableType::evolution) const + { + CH_assert("This should not be needed."); + return *m_fake; + } + + bool contains(const std::array &point) const + { + bool in = true; + for (int idir = 0; idir < N_DIMS; ++idir) + { + if (!m_is_periodic[idir] && + (point[idir] < 0. || point[idir] > m_points_per_dir[idir] - 1)) + in = false; + } + return in; + }; + + double get_dx(int dir) const { return m_dxs[dir]; }; + std::array get_dxs() const { return m_dxs; } +}; + +#endif /* SIMPLEINTERPSOURCE_H_ */ diff --git a/Source/GRChomboCore/GRAMR.cpp b/Source/GRChomboCore/GRAMR.cpp index a8bf8b070..eac679ffe 100644 --- a/Source/GRChomboCore/GRAMR.cpp +++ b/Source/GRChomboCore/GRAMR.cpp @@ -11,6 +11,8 @@ GRAMR::GRAMR() : m_interpolator(nullptr) {} // Called after AMR object set up void GRAMR::set_interpolator(AMRInterpolator> *a_interpolator) { + // check AMRInterpolator was set + CH_assert(a_interpolator != nullptr); m_interpolator = a_interpolator; } diff --git a/Source/GRChomboCore/GRAMR.hpp b/Source/GRChomboCore/GRAMR.hpp index 06f51614e..4cf35c894 100644 --- a/Source/GRChomboCore/GRAMR.hpp +++ b/Source/GRChomboCore/GRAMR.hpp @@ -65,6 +65,9 @@ class GRAMR : public AMR // const version of above std::vector get_gramrlevels() const; + int get_max_level() const { return m_max_level; } + int get_restart_step() const { return m_restart_step; } + // Fill ghosts on multiple levels void fill_multilevel_ghosts( const VariableType a_var_type, diff --git a/Source/GRChomboCore/GRAMRLevel.cpp b/Source/GRChomboCore/GRAMRLevel.cpp index c0c061784..ee6df819d 100644 --- a/Source/GRChomboCore/GRAMRLevel.cpp +++ b/Source/GRChomboCore/GRAMRLevel.cpp @@ -221,30 +221,27 @@ void GRAMRLevel::tagCells(IntVectSet &a_tags) const IntVect &smallEnd = b.smallEnd(); const IntVect &bigEnd = b.bigEnd(); - const int xmin = smallEnd[0]; - const int ymin = smallEnd[1]; - const int zmin = smallEnd[2]; - - const int xmax = bigEnd[0]; - const int ymax = bigEnd[1]; - const int zmax = bigEnd[2]; - -#pragma omp parallel for collapse(3) schedule(static) default(shared) - for (int iz = zmin; iz <= zmax; ++iz) - for (int iy = ymin; iy <= ymax; ++iy) - for (int ix = xmin; ix <= xmax; ++ix) - { - IntVect iv(ix, iy, iz); - if (tagging_criterion(iv, 0) >= - m_p.regrid_thresholds[m_level]) - { + D_TERM(const int xmin = smallEnd[0];, const int ymin = smallEnd[1]; + , const int zmin = smallEnd[2];) + + D_TERM(const int xmax = bigEnd[0];, const int ymax = bigEnd[1]; + , const int zmax = bigEnd[2];) + +#pragma omp parallel for collapse(CH_SPACEDIM) schedule(static) default(shared) + D_INVTERM(for (int ix = xmin; ix <= xmax; ++ix), + for (int iy = ymin; iy <= ymax; ++iy), + for (int iz = zmin; iz <= zmax; ++iz)) + { + IntVect iv(D_DECL(ix, iy, iz)); + if (tagging_criterion(iv, 0) >= m_p.regrid_thresholds[m_level]) + { // local_tags |= is not thread safe. #pragma omp critical - { - local_tags |= iv; - } - } + { + local_tags |= iv; } + } + } } local_tags.grow(m_p.tag_buffer_size); @@ -985,8 +982,6 @@ void GRAMRLevel::copySolnData(GRLevelData &dest, const GRLevelData &src) copyBdyGhosts(src, dest); } -double GRAMRLevel::get_dx() const { return m_dx; } - bool GRAMRLevel::at_level_timestep_multiple(int a_level) const { double target_dt = m_p.coarsest_dt; diff --git a/Source/GRChomboCore/GRAMRLevel.hpp b/Source/GRChomboCore/GRAMRLevel.hpp index cc99d1d9c..5b59e537c 100644 --- a/Source/GRChomboCore/GRAMRLevel.hpp +++ b/Source/GRChomboCore/GRAMRLevel.hpp @@ -28,7 +28,7 @@ // Chombo namespace #include "UsingNamespace.H" -class GRAMRLevel : public AMRLevel, public InterpSource +class GRAMRLevel : public AMRLevel, public InterpSource<> { public: GRAMRLevel(GRAMR &gr_amr, const SimulationParameters &a_p, int a_verbosity); @@ -165,7 +165,8 @@ class GRAMRLevel : public AMRLevel, public InterpSource { } - double get_dx() const; + // direction irrelevant, but relevant for InterpSource + ALWAYS_INLINE double get_dx(int dir = 0) const { return m_dx; }; /// Returns true if m_time is the same as the time at the end of the current /// timestep on level a_level and false otherwise diff --git a/Tests/InterpolatorTest/GNUmakefile b/Tests/InterpolatorTest/GNUmakefile new file mode 100644 index 000000000..d89500efb --- /dev/null +++ b/Tests/InterpolatorTest/GNUmakefile @@ -0,0 +1,20 @@ +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = $(shell pwd)/../../Source + +ebase := InterpolatorTest + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/AMRInterpolator + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Tests/InterpolatorTest/InterpolatorTest.cpp b/Tests/InterpolatorTest/InterpolatorTest.cpp new file mode 100644 index 000000000..d986f329f --- /dev/null +++ b/Tests/InterpolatorTest/InterpolatorTest.cpp @@ -0,0 +1,123 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifdef CH_LANG_CC +/* + * _______ __ + * / ___/ / ___ __ _ / / ___ + * / /__/ _ \/ _ \/ V \/ _ \/ _ \ + * \___/_//_/\___/_/_/_/_.__/\___/ + * Please refer to LICENSE, in Chombo's root directory. + */ +#endif + +// Chombo includes +#include "parstream.H" //Gives us pout() + +// General includes: +#include +#include +#include +#include +#include + +#include "Lagrange.hpp" +#include "QuinticConvolution.hpp" +#include "SetupFunctions.hpp" +#include "SimpleArrayBox.hpp" +#include "SimpleInterpSource.hpp" + +// Chombo namespace +#include "UsingNamespace.H" + +// cell centered grid +// double get_dx(double L, int num_points_u) { return L / num_points_u; } +// double get_x(int i, double m_dx) { return m_dx * (i + 0.5); } +// double get_idx(double x, double m_dx) { return x / m_dx - 0.5; } + +// node centered grid +double get_dx(double L, int num_points_u) { return L / (num_points_u - 1.); } +double get_x(int i, double m_dx) { return m_dx * i; } +double get_idx(double x, double m_dx) { return x / m_dx; } + +double func(double x, double L) +{ + // cubic + // maximum is 1 at x = L / 6 * ( 3 +- sqrt(3) ) + // return x * (x - L / 2.) * (x - L) / (L * L * L) * 12. * sqrt(3.); + + // quadratic + // maximum is 1 at x = L / 2 + // return x * (x - L) / (L * L) * 4.; + + // non-linear + return sin(2. * M_PI * x / L); +} + +int runInterpolatorTest(int argc, char *argv[]) +{ + + const int num_points_u = 200; + const double L = 10.; + const double m_dx = get_dx(L, num_points_u); + + std::vector f(num_points_u); + + for (int i = 0; i < num_points_u; ++i) + { + double x = get_x(i, m_dx); + f[i] = func(x, L); + } + + // Note that the 2nd argument of .setup 'dx' does not matter for 0th order + // interpolation -> set to 0 + SimpleInterpSource<1> source({num_points_u}, {0.}); + SimpleArrayBox<1> box({num_points_u}, f); + + double test_point = + M_PI; // just because it's irrational, hence for sure not in the grid + double test_index = get_idx(test_point, m_dx); + + bool verbosity = true; + Lagrange<4, 1> interpolator1(source, verbosity); + interpolator1.setup({0}, {test_index}); + double val1 = interpolator1.interpData(box); + + QuinticConvolution<1> interpolator2(source, verbosity); + interpolator2.setup({0}, {test_index}); + double val2 = interpolator2.interpData(box); + + double exact = func(test_point, L); + + double err1 = abs(val1 / exact - 1.) * 100.; + double err2 = abs(val2 / exact - 1.) * 100.; + + pout() << std::setprecision(9); + pout() << "Exact value is: " << exact << std::endl; + pout() << "Lagrange Interpolator gave: " << val1 << " (" << err1 + << "% error)" << std::endl; + pout() << "Quintic Convolution Interpolator gave: " << val2 << " (" << err2 + << "% error)" << std::endl; + + bool wrong = (err1 + err2 > 1e-5); + + return wrong; +} + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); + + int status = runInterpolatorTest(argc, argv); + + if (status == 0) + std::cout << "Interpolator test passed." << endl; + else + std::cout << "Interpolator test failed with return code " << status + << endl; + + mainFinalize(); + return status; +} diff --git a/Tests/InterpolatorTest/InterpolatorTest.inputs b/Tests/InterpolatorTest/InterpolatorTest.inputs new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/Tests/InterpolatorTest/InterpolatorTest.inputs @@ -0,0 +1 @@ + diff --git a/Tests/InterpolatorTest/SimulationParameters.hpp b/Tests/InterpolatorTest/SimulationParameters.hpp new file mode 100644 index 000000000..5e00a7de6 --- /dev/null +++ b/Tests/InterpolatorTest/SimulationParameters.hpp @@ -0,0 +1,19 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMULATIONPARAMETERS_HPP_ +#define SIMULATIONPARAMETERS_HPP_ + +// General includes +#include "ChomboParameters.hpp" +#include "GRParmParse.hpp" + +class SimulationParameters : public ChomboParameters +{ + public: + SimulationParameters(GRParmParse &pp) : ChomboParameters(pp) {} +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Tests/InterpolatorTest/UserVariables.hpp b/Tests/InterpolatorTest/UserVariables.hpp new file mode 100644 index 000000000..f7642667d --- /dev/null +++ b/Tests/InterpolatorTest/UserVariables.hpp @@ -0,0 +1,26 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "EmptyDiagnosticVariables.hpp" +#include +#include + +// assign enum to each variable +enum +{ + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array variable_names = {}; +} + +#include "UserVariables.inc.hpp" + +#endif /* USERVARIABLES_HPP */