From c8294a8499cc8bbd616534aa0bfdf88c0336e731 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Thu, 27 May 2021 17:18:58 +0100 Subject: [PATCH] Generalize Interpolator classes This is an attempt to allow the interpolation algorithms (namely Lagrange) to be used outside of the AMRInterpolator to anyone that needs it and to the AHFinder. The algorithm we have is prepared for any dimension, but we only use it for CH_SPACEDIM. If we want to interpolate something else, this should be a ready to use class. That is the purpose of this commit. The algorithm is actually "ready" to even do non-uniform grids ('dx' not constant), but changing that required the structure to change too much. For these alternative interpolations, the classes 'SimpleArrayBox' and 'SimpleInterpSource' were added, mimicking 'FArrayBox' and 'GRAMRLevel'. They represent a D-dimensional grid that can be provided based on a flattened 1D array, potentially with periodic BCs and constant 'dx' (per direction). Use cases: 1) interpolating the AHFinder AH centers from the 'stats' output when frequency of solving the AHFinder changes at a restart (1D interpolation) 2) interpolating the AH surface when 'num_points_u' or 'num_points_v' change at a restart (2D interpolation) 3) interpolating spherically symmetric initial data into the 3D grid (1D interpolation) --- Source/AMRInterpolator/AMRInterpolator.hpp | 26 +++- .../AMRInterpolator/AMRInterpolator.impl.hpp | 73 ++++++----- Source/AMRInterpolator/InterpSource.hpp | 15 ++- Source/AMRInterpolator/Lagrange.hpp | 27 ++-- Source/AMRInterpolator/Lagrange.impl.hpp | 98 +++++++------- Source/AMRInterpolator/QuinticConvolution.hpp | 20 +-- .../QuinticConvolution.impl.hpp | 91 +++++++------ Source/AMRInterpolator/SimpleArrayBox.hpp | 89 +++++++++++++ Source/AMRInterpolator/SimpleInterpSource.hpp | 62 +++++++++ Source/GRChomboCore/GRAMR.cpp | 2 + Source/GRChomboCore/GRAMR.hpp | 3 + Source/GRChomboCore/GRAMRLevel.cpp | 41 +++--- Source/GRChomboCore/GRAMRLevel.hpp | 5 +- Tests/InterpolatorTest/GNUmakefile | 20 +++ Tests/InterpolatorTest/InterpolatorTest.cpp | 123 ++++++++++++++++++ .../InterpolatorTest/InterpolatorTest.inputs | 1 + .../InterpolatorTest/SimulationParameters.hpp | 19 +++ Tests/InterpolatorTest/UserVariables.hpp | 26 ++++ 18 files changed, 565 insertions(+), 176 deletions(-) create mode 100644 Source/AMRInterpolator/SimpleArrayBox.hpp create mode 100644 Source/AMRInterpolator/SimpleInterpSource.hpp create mode 100644 Tests/InterpolatorTest/GNUmakefile create mode 100644 Tests/InterpolatorTest/InterpolatorTest.cpp create mode 100644 Tests/InterpolatorTest/InterpolatorTest.inputs create mode 100644 Tests/InterpolatorTest/SimulationParameters.hpp create mode 100644 Tests/InterpolatorTest/UserVariables.hpp 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 1b2fca558..fa8ee0fca 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 */