Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize Interpolator classes #176

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Generalize Interpolator classes #176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 19 additions & 7 deletions Source/AMRInterpolator/AMRInterpolator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand All @@ -31,6 +28,19 @@

template <typename InterpAlgo> class AMRInterpolator
{
struct InterpolationLayout
{
std::vector<int> rank;
std::vector<int> level_idx;
std::vector<int> box_idx;

InterpolationLayout(int num_points)
: rank(num_points, -1), level_idx(num_points, -1),
box_idx(num_points, -1)
{
}
};
Comment on lines +31 to +42
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I'm ok with moving this here but we should delete the InterpolationLayout.hpp file too.


public:
// constructor for backward compatibility
// (adds an artificial BC with only periodic BC)
Expand All @@ -56,15 +66,17 @@ template <typename InterpAlgo> class AMRInterpolator
void limit_num_levels(unsigned int num_levels);
void interp(InterpolationQuery &query);
const AMR &getAMR() const;
const std::array<double, CH_SPACEDIM> &get_coarsest_dx();
const std::array<double, CH_SPACEDIM> &get_coarsest_origin();
const std::array<double, CH_SPACEDIM> &get_coarsest_dx() const;
const std::array<double, CH_SPACEDIM> &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();
Expand Down
73 changes: 37 additions & 36 deletions Source/AMRInterpolator/AMRInterpolator.impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_

Expand Down Expand Up @@ -41,9 +45,7 @@ void AMRInterpolator<InterpAlgo>::refresh(const bool a_fill_ghosts)
{
CH_TIME("AMRInterpolator::refresh");

const Vector<AMRLevel *> &levels =
const_cast<GRAMR &>(m_gr_amr).getAMRLevels();
m_num_levels = levels.size();
m_num_levels = const_cast<GRAMR &>(m_gr_amr).getAMRLevels().size();

m_mem_level.clear();
m_mem_box.clear();
Expand Down Expand Up @@ -73,17 +75,32 @@ const AMR &AMRInterpolator<InterpAlgo>::getAMR() const

template <typename InterpAlgo>
const std::array<double, CH_SPACEDIM> &
AMRInterpolator<InterpAlgo>::get_coarsest_dx()
AMRInterpolator<InterpAlgo>::get_coarsest_dx() const
{
return m_coarsest_dx;
}
template <typename InterpAlgo>
const std::array<double, CH_SPACEDIM> &
AMRInterpolator<InterpAlgo>::get_coarsest_origin()
AMRInterpolator<InterpAlgo>::get_coarsest_origin() const
{
return m_coarsest_origin;
}

template <typename InterpAlgo>
bool AMRInterpolator<InterpAlgo>::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 <typename InterpAlgo>
bool AMRInterpolator<InterpAlgo>::get_boundary_periodic(int a_dir) const
{
return m_bc_params.is_periodic[a_dir];
}

template <typename InterpAlgo>
void AMRInterpolator<InterpAlgo>::limit_num_levels(unsigned int num_levels)
{
Expand Down Expand Up @@ -239,10 +256,14 @@ void AMRInterpolator<InterpAlgo>::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;
}

Expand All @@ -257,7 +278,7 @@ void AMRInterpolator<InterpAlgo>::computeLevelLayouts()
}

template <typename InterpAlgo>
InterpolationLayout
typename AMRInterpolator<InterpAlgo>::InterpolationLayout
AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
{
CH_TIME("AMRInterpolator::findBoxes");
Expand Down Expand Up @@ -293,7 +314,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
// const AMRLevel& level = *levels[level_idx];

// const LevelData<FArrayBox>& level_data = dynamic_cast<const
// InterpSource&>(level).getLevelData(); const DisjointBoxLayout&
// InterpSource<>&>(level).getLevelData(); const DisjointBoxLayout&
// box_layout = level_data.disjointBoxLayout(); const Box& domain_box =
// level.problemDomain().domainBox();

Expand Down Expand Up @@ -381,7 +402,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)
const AMRLevel &level = *levels[level_idx];

const LevelData<FArrayBox> &level_data =
dynamic_cast<const InterpSource &>(level).getLevelData(
dynamic_cast<const InterpSource<> &>(level).getLevelData(
VariableType::evolution);
const DisjointBoxLayout &box_layout = level_data.disjointBoxLayout();
const Box &domain_box = level.problemDomain().domainBox();
Expand Down Expand Up @@ -495,7 +516,7 @@ AMRInterpolator<InterpAlgo>::findBoxes(InterpolationQuery &query)

template <typename InterpAlgo>
void AMRInterpolator<InterpAlgo>::prepareMPI(InterpolationQuery &query,
const InterpolationLayout layout)
const InterpolationLayout &layout)
{
CH_TIME("AMRInterpolator::prepareMPI");

Expand Down Expand Up @@ -631,15 +652,15 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)
const int num_answers = m_mpi.totalAnswerCount();

std::array<double, CH_SPACEDIM> grid_coord;
IntVect nearest;

for (int answer_idx = 0; answer_idx < num_answers; ++answer_idx)
{
const int box_idx = m_answer_box[answer_idx];
const int level_idx = m_answer_level[answer_idx];

const AMRLevel &level = *levels[level_idx];
const InterpSource &source = dynamic_cast<const InterpSource &>(level);
const InterpSource<> &source =
dynamic_cast<const InterpSource<> &>(level);
const LevelData<FArrayBox> *const evolution_level_data_ptr =
&source.getLevelData(VariableType::evolution);
const LevelData<FArrayBox> *diagnostics_level_data_ptr;
Expand All @@ -657,10 +678,6 @@ void AMRInterpolator<InterpAlgo>::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]);
Expand Down Expand Up @@ -696,22 +713,6 @@ void AMRInterpolator<InterpAlgo>::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)
Expand All @@ -737,7 +738,7 @@ void AMRInterpolator<InterpAlgo>::calculateAnswers(InterpolationQuery &query)

typedef std::vector<typename InterpolationQuery::out_t> 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)
Expand Down
15 changes: 11 additions & 4 deletions Source/AMRInterpolator/InterpSource.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might as well fix the typo while we're at it

Suggested change
// Abstrace base class to identify points in a grid structure
// Abstract base class to identify points in a grid structure

template <int N_DIMS = CH_SPACEDIM> class InterpSource
{
public:
virtual const LevelData<FArrayBox> &getLevelData(
const VariableType var_type = VariableType::evolution) const = 0;
virtual bool
contains(const std::array<double, CH_SPACEDIM> &point) const = 0;
virtual bool contains(const std::array<double, N_DIMS> &point) const = 0;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add virtual and override specifiers to the corresponding override in GRAMRLevel?

virtual double get_dx(int dir) const = 0;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like this function has to be defined even if the specialised class represents an InterpSource with grid spacings that are different in each direction which is not ideal. Perhaps only the one below should be pure virtual and this can be removed?

virtual std::array<double, N_DIMS> get_dxs() const
{
std::array<double, N_DIMS> dxs;
for (int i = 0; i < N_DIMS; ++i)
dxs[i] = get_dx(i);
return dxs;
}
};

#endif /* INTERPSOURCE_H_ */
27 changes: 15 additions & 12 deletions Source/AMRInterpolator/Lagrange.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@
#include "InterpSource.hpp"
#include <utility>

template <int Order> class Lagrange
template <int Order, int N_DIMS = CH_SPACEDIM> class Lagrange
{
const InterpSource &m_source;
const InterpSource<N_DIMS> &m_source;
bool m_verbosity;

struct Stencil;
Expand Down Expand Up @@ -41,10 +41,10 @@ template <int Order> class Lagrange
// Helper function to generate tensor product weights
// Argument 'dim' is used for recursion over dimensions.
pair<std::vector<IntVect>, std::vector<double>>
generateStencil(const std::array<int, CH_SPACEDIM> &deriv,
const std::array<double, CH_SPACEDIM> &dx,
const std::array<double, CH_SPACEDIM> &evalCoord,
const IntVect &nearest, int dim = CH_SPACEDIM - 1);
generateStencil(const std::array<int, N_DIMS> &deriv,
const std::array<double, N_DIMS> &dx,
const std::array<double, N_DIMS> &eval_index,
int dim = N_DIMS - 1);

std::vector<IntVect> m_interp_points;
std::vector<double> m_interp_weights;
Expand All @@ -56,13 +56,16 @@ template <int Order> class Lagrange
multiset<double> m_interp_pos;

public:
Lagrange(const InterpSource &source, bool verbosity = false);
Lagrange(const InterpSource<N_DIMS> &source, bool verbosity = false);

void setup(const std::array<int, CH_SPACEDIM> &deriv,
const std::array<double, CH_SPACEDIM> &dx,
const std::array<double, CH_SPACEDIM> &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<int, N_DIMS> &deriv,
const std::array<double, N_DIMS> &eval_index);

// any class with a method:
// Real get(const IntVect &a_iv, int a_comp) const
template <class GeneralArrayBox>
double interpData(const GeneralArrayBox &fab, int comp = 0);

const static string TAG;
};
Expand Down
Loading