Skip to content

Weyl inv 3 #4

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

Open
cjoana wants to merge 3 commits into
base: WeylInvariants
Choose a base branch
from
Open

Weyl inv 3 #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
11f9a2a
evo working, Weyl4 needed
cjoana Jul 4, 2023
1240a2c
working evo and weyl
cjoana Jul 4, 2023
dc332f5
cleaning oSorce and del newGeometry
cjoana Jul 4, 2023
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
7 changes: 2 additions & 5 deletions Source/AMRInterpolator/AMRInterpolator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@

// Chombo includes
#include "AMRLevel.H"
#include "LoHiSide.H"

// Our includes
#include "BoundaryConditions.hpp"
Expand Down Expand Up @@ -57,10 +56,8 @@ 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;
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;
const std::array<double, CH_SPACEDIM> &get_coarsest_dx();
const std::array<double, CH_SPACEDIM> &get_coarsest_origin();

private:
void computeLevelLayouts();
Expand Down
71 changes: 34 additions & 37 deletions Source/AMRInterpolator/AMRInterpolator.impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,6 @@
* 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 @@ -77,38 +73,23 @@ const AMR &AMRInterpolator<InterpAlgo>::getAMR() const

template <typename InterpAlgo>
const std::array<double, CH_SPACEDIM> &
AMRInterpolator<InterpAlgo>::get_coarsest_dx() const
AMRInterpolator<InterpAlgo>::get_coarsest_dx()
{
return m_coarsest_dx;
}
template <typename InterpAlgo>
const std::array<double, CH_SPACEDIM> &
AMRInterpolator<InterpAlgo>::get_coarsest_origin() const
AMRInterpolator<InterpAlgo>::get_coarsest_origin()
{
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)
{
CH_TIME("AMRInterpolator::limit_num_levels");

int max_num_levels = const_cast<AMR &>(m_gr_amr).getAMRLevels().size();
int max_num_levels = const_cast<GRAMR &>(m_gr_amr).getAMRLevels().size();
if (num_levels > max_num_levels || num_levels == 0)
{
m_num_levels = max_num_levels;
Expand Down Expand Up @@ -258,14 +239,10 @@ void AMRInterpolator<InterpAlgo>::computeLevelLayouts()
if (m_verbosity >= 2)
{
_pout << " Level " << level_idx << "\t"
<< "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])
<< "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]
<< ")" << endl;
}

Expand Down Expand Up @@ -316,7 +293,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 @@ -404,7 +381,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 @@ -654,15 +631,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 @@ -680,6 +657,10 @@ 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 @@ -715,6 +696,22 @@ 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 @@ -740,7 +737,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);
algo.setup(deriv, m_dx[level_idx], grid_coord, nearest);

for (typename comps_t::iterator it = comps.begin();
it != comps.end(); ++it)
Expand Down Expand Up @@ -809,7 +806,7 @@ void AMRInterpolator<InterpAlgo>::set_reflective_BC()
.domainBox()
.bigEnd();

FOR1(i)
FOR(i)
{
m_upper_corner[i] = (big_end[i] + 1) * m_coarsest_dx[i];

Expand Down Expand Up @@ -837,7 +834,7 @@ int AMRInterpolator<InterpAlgo>::get_var_parity(int comp,
"extracting diagnostic variables with reflective BC");

int parity = 1;
FOR1(dir)
FOR(dir)
{
double coord = query.m_coords[dir][point_idx];
if ((m_lo_boundary_reflective[dir] && coord < 0.) ||
Expand Down
38 changes: 15 additions & 23 deletions Source/AMRInterpolator/CylindricalGeometry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,43 +32,36 @@ class CylindricalGeometry
{
}

ALWAYS_INLINE double get_domain_u_min() const { return 0.; }
ALWAYS_INLINE double get_domain_u_max() const { return m_z_length; }
ALWAYS_INLINE double get_domain_v_min() const { return 0.; }
ALWAYS_INLINE double get_domain_v_max() const { return 2. * M_PI; }

//! returns the grid spacing in z
ALWAYS_INLINE double du(int a_num_points_z) const
inline double du(int a_num_points_z) const
{
return (get_domain_u_max() - get_domain_u_min()) / (a_num_points_z - 1);
return m_z_length / (a_num_points_z - 1);
}

//! returns the grid spacing in phi
ALWAYS_INLINE double dv(int a_num_points_phi) const
inline double dv(int a_num_points_phi) const
{
return (get_domain_v_max() - get_domain_v_min()) /
((double)a_num_points_phi);
return 2.0 * M_PI / ((double)a_num_points_phi);
}

//! returns the z coordinate associated to the theta/u index
ALWAYS_INLINE double u(int a_iz, int a_num_points_z) const
inline double u(int a_iz, int a_num_points_z) const
{
return a_iz * du(a_num_points_z) -
((get_domain_v_max() - get_domain_v_min()) / 2.0);
return a_iz * du(a_num_points_z) - (m_z_length / 2.0);
}

//! returns the phi coordinate associated to the phi/v index
ALWAYS_INLINE double v(int a_iphi, int a_num_points_phi) const
inline double v(int a_iphi, int a_num_points_phi) const
{
return a_iphi * dv(a_num_points_phi);
}

ALWAYS_INLINE bool is_u_periodic() const { return false; }
ALWAYS_INLINE bool is_v_periodic() const { return true; }
inline bool is_u_periodic() const { return false; }
inline bool is_v_periodic() const { return true; }
//! returns the Cartesian coordinate in direction a_dir with specified
//! radius, z and phi.
ALWAYS_INLINE double get_grid_coord(int a_dir, double a_radius, double a_z,
double a_phi) const
inline double get_grid_coord(int a_dir, double a_radius, double a_z,
double a_phi) const
{
switch (a_dir)
{
Expand All @@ -85,17 +78,16 @@ class CylindricalGeometry

//! returns the area element on a cylinder with radius a_radius at the point
//! (a_z, a_phi)
ALWAYS_INLINE double area_element(double a_radius, double a_z,
double a_phi) const
inline double area_element(double a_radius, double a_z, double a_phi) const
{
return a_radius;
}

ALWAYS_INLINE std::string param_name() const { return "r"; }
inline std::string param_name() const { return "r"; }

ALWAYS_INLINE std::string u_name() const { return "z"; }
inline std::string u_name() const { return "z"; }

ALWAYS_INLINE std::string v_name() const { return "phi"; }
inline std::string v_name() const { return "phi"; }
};

#endif /* CYLINDRICALGEOMETRY_HPP_ */
18 changes: 7 additions & 11 deletions Source/AMRInterpolator/Derivative.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,40 +108,36 @@ class Derivative : public std::array<int, CH_SPACEDIM>

static const Derivative dx;
static const Derivative dy;
static const Derivative dz;

static const Derivative dxdx;
static const Derivative dxdy;
static const Derivative dydy;
static const Derivative dzdz;

#if CH_SPACEDIM == 3
static const Derivative dz;
static const Derivative dxdy;
static const Derivative dxdz;
static const Derivative dydz;
static const Derivative dzdz;
#endif

static std::string name(const Derivative &deriv)
{
if (deriv == dx)
return "dx";
else if (deriv == dy)
return "dy";
else if (deriv == dz)
return "dz";
else if (deriv == dxdx)
return "dxdx";
else if (deriv == dxdy)
return "dxdy";
else if (deriv == dydy)
return "dydy";
#if CH_SPACEDIM == 3
else if (deriv == dz)
return "dz";
else if (deriv == dzdz)
return "dzdz";
else if (deriv == dxdy)
return "dxdy";
else if (deriv == dxdz)
return "dxdz";
else if (deriv == dydz)
return "dydz";
#endif
else
return "";
}
Expand Down
8 changes: 3 additions & 5 deletions Source/AMRInterpolator/DerivativeSetup.hpp
100755 → 100644
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,14 @@ const Derivative Derivative::LOCAL;

const Derivative Derivative::dx(0);
const Derivative Derivative::dy(1);
const Derivative Derivative::dz(2);

const Derivative Derivative::dxdx(0, 0);
const Derivative Derivative::dxdy(0, 1);
const Derivative Derivative::dydy(1, 1);
const Derivative Derivative::dzdz(2, 2);

#if CH_SPACEDIM == 3
const Derivative Derivative::dz(2);
const Derivative Derivative::dxdy(0, 1);
const Derivative Derivative::dxdz(0, 2);
const Derivative Derivative::dydz(1, 2);
const Derivative Derivative::dzdz(2, 2);
#endif

#endif /* DERIVATIVESETUP_HPP_ */
5 changes: 3 additions & 2 deletions Source/AMRInterpolator/InterpSource.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,13 @@
#include "UsingNamespace.H"

// Abstrace base class to get the FABs out of an AMRLevel
template <int N_DIMS = CH_SPACEDIM> class InterpSource
class InterpSource
{
public:
virtual const LevelData<FArrayBox> &getLevelData(
const VariableType var_type = VariableType::evolution) const = 0;
virtual bool contains(const std::array<double, N_DIMS> &point) const = 0;
virtual bool
contains(const std::array<double, CH_SPACEDIM> &point) const = 0;
};

#endif /* INTERPSOURCE_H_ */
Loading