diff --git a/Source/AMRInterpolator/AMRInterpolator.hpp b/Source/AMRInterpolator/AMRInterpolator.hpp index f68a0a143..b8b63c97d 100644 --- a/Source/AMRInterpolator/AMRInterpolator.hpp +++ b/Source/AMRInterpolator/AMRInterpolator.hpp @@ -12,7 +12,6 @@ // Chombo includes #include "AMRLevel.H" -#include "LoHiSide.H" // Our includes #include "BoundaryConditions.hpp" @@ -57,10 +56,8 @@ 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; - 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; + const std::array &get_coarsest_dx(); + const std::array &get_coarsest_origin(); private: void computeLevelLayouts(); diff --git a/Source/AMRInterpolator/AMRInterpolator.impl.hpp b/Source/AMRInterpolator/AMRInterpolator.impl.hpp index 94d5be946..22d87d36d 100644 --- a/Source/AMRInterpolator/AMRInterpolator.impl.hpp +++ b/Source/AMRInterpolator/AMRInterpolator.impl.hpp @@ -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_ @@ -77,38 +73,23 @@ const AMR &AMRInterpolator::getAMR() const template const std::array & -AMRInterpolator::get_coarsest_dx() const +AMRInterpolator::get_coarsest_dx() { return m_coarsest_dx; } template const std::array & -AMRInterpolator::get_coarsest_origin() const +AMRInterpolator::get_coarsest_origin() { 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) { CH_TIME("AMRInterpolator::limit_num_levels"); - int max_num_levels = const_cast(m_gr_amr).getAMRLevels().size(); + int max_num_levels = const_cast(m_gr_amr).getAMRLevels().size(); if (num_levels > max_num_levels || num_levels == 0) { m_num_levels = max_num_levels; @@ -258,14 +239,10 @@ void AMRInterpolator::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; } @@ -316,7 +293,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(); @@ -404,7 +381,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(); @@ -654,6 +631,7 @@ 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) { @@ -661,8 +639,7 @@ 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; @@ -680,6 +657,10 @@ 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]); @@ -715,6 +696,22 @@ 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) @@ -740,7 +737,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); + algo.setup(deriv, m_dx[level_idx], grid_coord, nearest); for (typename comps_t::iterator it = comps.begin(); it != comps.end(); ++it) @@ -809,7 +806,7 @@ void AMRInterpolator::set_reflective_BC() .domainBox() .bigEnd(); - FOR1(i) + FOR(i) { m_upper_corner[i] = (big_end[i] + 1) * m_coarsest_dx[i]; @@ -837,7 +834,7 @@ int AMRInterpolator::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.) || diff --git a/Source/AMRInterpolator/CylindricalGeometry.hpp b/Source/AMRInterpolator/CylindricalGeometry.hpp index 8e1195580..14b34ae1f 100644 --- a/Source/AMRInterpolator/CylindricalGeometry.hpp +++ b/Source/AMRInterpolator/CylindricalGeometry.hpp @@ -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) { @@ -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_ */ diff --git a/Source/AMRInterpolator/Derivative.hpp b/Source/AMRInterpolator/Derivative.hpp index 2e9808224..f9b8943e0 100644 --- a/Source/AMRInterpolator/Derivative.hpp +++ b/Source/AMRInterpolator/Derivative.hpp @@ -108,17 +108,15 @@ class Derivative : public std::array 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) { @@ -126,22 +124,20 @@ class Derivative : public std::array 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 ""; } diff --git a/Source/AMRInterpolator/DerivativeSetup.hpp b/Source/AMRInterpolator/DerivativeSetup.hpp old mode 100755 new mode 100644 index 6b5cf0b4c..9d22c1c65 --- a/Source/AMRInterpolator/DerivativeSetup.hpp +++ b/Source/AMRInterpolator/DerivativeSetup.hpp @@ -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_ */ diff --git a/Source/AMRInterpolator/InterpSource.hpp b/Source/AMRInterpolator/InterpSource.hpp index 2a53c5848..90c35a1c8 100644 --- a/Source/AMRInterpolator/InterpSource.hpp +++ b/Source/AMRInterpolator/InterpSource.hpp @@ -17,12 +17,13 @@ #include "UsingNamespace.H" // Abstrace base class to get the FABs out of an AMRLevel -template class InterpSource +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; }; #endif /* INTERPSOURCE_H_ */ diff --git a/Source/AMRInterpolator/Lagrange.hpp b/Source/AMRInterpolator/Lagrange.hpp index 199ada2ea..c54aa13a6 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, - int dim = N_DIMS - 1); + generateStencil(const std::array &deriv, + const std::array &dx, + const std::array &evalCoord, + const IntVect &nearest, int dim = CH_SPACEDIM - 1); std::vector m_interp_points; std::vector m_interp_weights; @@ -56,17 +56,13 @@ template class Lagrange multiset m_interp_pos; public: - Lagrange(const InterpSource &source, bool verbosity = false); + Lagrange(const InterpSource &source, bool verbosity = false); - // evalCoord is in 'index' coordinates, not physical coordinates - void setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord); - - // any class with a method: - // Real get(const IntVect &a_iv, int a_comp) const - template - double interpData(const GeneralArrayBox &fab, int comp = 0); + void setup(const std::array &deriv, + const std::array &dx, + const std::array &evalCoord, + const IntVect &nearest); + double interpData(const FArrayBox &fab, int comp); const static string TAG; }; diff --git a/Source/AMRInterpolator/Lagrange.impl.hpp b/Source/AMRInterpolator/Lagrange.impl.hpp index 8b0cb99ba..48bc7ec68 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,26 +151,26 @@ 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) +template +void Lagrange::setup(const std::array &deriv, + const std::array &dx, + const std::array &evalCoord, + const IntVect &nearest) { pair, std::vector> result = - generateStencil(deriv, dx, evalCoord); + generateStencil(deriv, dx, evalCoord, nearest); m_interp_points = result.first; m_interp_weights = result.second; /* pout() << TAG << "Stencil: coord = { "; - for (int i = 0; i < N_DIMS; ++i) + for (int i = 0; i < CH_SPACEDIM; ++i) { pout() << evalCoord[i] << " "; } @@ -183,9 +183,8 @@ void Lagrange::setup(const std::array &deriv, */ } -template -template -double Lagrange::interpData(const GeneralArrayBox &fab, int comp) +template +double Lagrange::interpData(const FArrayBox &fab, int comp) { /* m_interp_neg.clear(); @@ -237,11 +236,13 @@ double Lagrange::interpData(const GeneralArrayBox &fab, int comp) return accum; } -template +template pair, std::vector> -Lagrange::generateStencil( - const std::array &deriv, const std::array &dx, - const std::array &evalCoord, int dim) +Lagrange::generateStencil( + const std::array &deriv, + const std::array &dx, + const std::array &evalCoord, const IntVect &nearest, + int dim) { std::vector out_points; std::vector out_weights; @@ -265,13 +266,9 @@ Lagrange::generateStencil( int points_min = Order + deriv[dim]; int points_max = Order + deriv[dim]; - std::array interp_coord = evalCoord; - - // TF: assume nearest point is at 'std::round(evalCoord[dim])' - // this used to be an argument, but I think this assumption should always - // hold - int candidate = std::round(evalCoord[dim]); - int grown_direction = (candidate - evalCoord[dim] < 0) ? DOWN : UP; + std::array interp_coord = evalCoord; + int candidate = nearest[dim]; + int grown_direction = (nearest[dim] - evalCoord[dim] < 0) ? DOWN : UP; while ((can_grow[DOWN] || can_grow[UP]) && (points_max - points_min < Order + deriv[dim])) @@ -333,7 +330,7 @@ Lagrange::generateStencil( { // Descend to the next dimension pair, std::vector> sub_result = - generateStencil(deriv, dx, interp_coord, dim - 1); + generateStencil(deriv, dx, interp_coord, nearest, dim - 1); std::vector &sub_points = sub_result.first; std::vector &sub_weights = sub_result.second; diff --git a/Source/AMRInterpolator/MPIContext.impl.hpp b/Source/AMRInterpolator/MPIContext.impl.hpp index b63cf267a..419a93b9d 100644 --- a/Source/AMRInterpolator/MPIContext.impl.hpp +++ b/Source/AMRInterpolator/MPIContext.impl.hpp @@ -68,7 +68,7 @@ inline void MPIContext::asyncExchangeQuery(void *sendbuf, void *recvbuf, MPI_Request req; m_mpi_requests.push_back(req); -#if MPI_VERSION >= 3 +#if MPI_VERSION >= 3 && !defined(OPEN_MPI) MPI_Ialltoallv(sendbuf, m_query.countsPtr(), m_query.displsPtr(), type, recvbuf, m_answer.countsPtr(), m_answer.displsPtr(), type, Chombo_MPI::comm, &m_mpi_requests.back()); @@ -86,7 +86,7 @@ inline void MPIContext::asyncExchangeAnswer(void *sendbuf, void *recvbuf, MPI_Request req; m_mpi_requests.push_back(req); -#if MPI_VERSION >= 3 +#if MPI_VERSION >= 3 && !defined(OPEN_MPI) MPI_Ialltoallv(sendbuf, m_answer.countsPtr(), m_answer.displsPtr(), type, recvbuf, m_query.countsPtr(), m_query.displsPtr(), type, Chombo_MPI::comm, &m_mpi_requests.back()); @@ -102,7 +102,7 @@ inline void MPIContext::asyncEnd() CH_assert(m_async_active); m_async_active = false; -#if MPI_VERSION >= 3 +#if MPI_VERSION >= 3 && !defined(OPEN_MPI) MPI_Waitall(m_mpi_requests.size(), &m_mpi_requests[0], MPI_STATUSES_IGNORE); #endif diff --git a/Source/AMRInterpolator/QuinticConvolution.hpp b/Source/AMRInterpolator/QuinticConvolution.hpp index ea92d574d..c01a6337a 100644 --- a/Source/AMRInterpolator/QuinticConvolution.hpp +++ b/Source/AMRInterpolator/QuinticConvolution.hpp @@ -9,27 +9,22 @@ #include "InterpSource.hpp" #include -template class QuinticConvolution +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); - - // evalCoord is in 'index' coordinates, not physical coordinates - void setup(const std::array &deriv, - const std::array &dx, - const std::array &evalCoord); - - // any class with a method: - // Real get(const IntVect &a_iv, int a_comp) const - template - double interpData(const GeneralArrayBox &fab, int comp = 0); + 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); const static string TAG; }; diff --git a/Source/AMRInterpolator/QuinticConvolution.impl.hpp b/Source/AMRInterpolator/QuinticConvolution.impl.hpp index b9c57a8a6..38c9783fa 100644 --- a/Source/AMRInterpolator/QuinticConvolution.impl.hpp +++ b/Source/AMRInterpolator/QuinticConvolution.impl.hpp @@ -6,29 +6,26 @@ #ifndef QUINTICCONVOLUTION_IMPL_HPP_ #define QUINTICCONVOLUTION_IMPL_HPP_ -template -const string QuinticConvolution::TAG = - "\x1b[36;1m[QuinticConvolution]\x1b[0m "; +const string QuinticConvolution::TAG = "\x1b[36;1m[QuinticConvolution]\x1b[0m "; -template -QuinticConvolution::QuinticConvolution( - const InterpSource &source, bool verbosity) +QuinticConvolution::QuinticConvolution(const InterpSource &source, + bool verbosity) : m_source(source), m_verbosity(verbosity) { - CH_assert(N_DIMS <= 3); + CH_assert(CH_SPACEDIM <= 3); } -template -void QuinticConvolution::setup( - const std::array &deriv, const std::array &dx, - const std::array &evalCoord) +void QuinticConvolution::setup(const std::array &deriv, + const std::array &dx, + const std::array &evalCoord, + const IntVect &nearest) { m_interp_points.clear(); m_interp_weights.clear(); - double weights_1d[N_DIMS][6]; + double weights_1d[CH_SPACEDIM][6]; - for (int dim = 0; dim < N_DIMS; ++dim) + for (int dim = 0; dim < CH_SPACEDIM; ++dim) { double s = evalCoord[dim] - floor(evalCoord[dim]); @@ -101,18 +98,16 @@ void QuinticConvolution::setup( } } - std::array interp_coord; + std::array interp_coord; -#if N_DIMS >= 3 +#if CH_SPACEDIM >= 3 for (int z = 0; z < 6; ++z) { interp_coord[2] = floor(evalCoord[2]) + z - 2; #endif -#if N_DIMS >= 2 for (int y = 0; y < 6; ++y) { interp_coord[1] = floor(evalCoord[1]) + y - 2; -#endif for (int x = 0; x < 6; ++x) { @@ -122,27 +117,17 @@ void QuinticConvolution::setup( 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]))); -#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 + m_interp_weights.push_back(D_TERM6(weights_1d[0][x], + *weights_1d[1][y], + *weights_1d[2][z], , , )); } -#if N_DIMS >= 2 } -#endif -#if N_DIMS >= 3 +#if CH_SPACEDIM >= 3 } #endif } -template -template -double QuinticConvolution::interpData(const GeneralArrayBox &fab, - int comp) +double QuinticConvolution::interpData(const FArrayBox &fab, int comp) { long double accum = 0.0; diff --git a/Source/AMRInterpolator/SimpleArrayBox.hpp b/Source/AMRInterpolator/SimpleArrayBox.hpp deleted file mode 100644 index 86ea592f8..000000000 --- a/Source/AMRInterpolator/SimpleArrayBox.hpp +++ /dev/null @@ -1,75 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef SIMPLEARRAYBOX -#define SIMPLEARRAYBOX - -#include "CH_assert.H" -#include "MayDay.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(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 deleted file mode 100644 index 07f8f8d43..000000000 --- a/Source/AMRInterpolator/SimpleInterpSource.hpp +++ /dev/null @@ -1,46 +0,0 @@ -/* 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 -{ - LevelData m_fake; - - std::array m_points_per_dir; - std::array m_is_periodic; - - public: - SimpleInterpSource(std::array a_points_per_dir, - std::array a_is_periodic = {false}) - : m_points_per_dir(a_points_per_dir), m_is_periodic(a_is_periodic) - { - } - - const LevelData & - getLevelData(const VariableType var_type = VariableType::evolution) const - { - 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; - }; - void fillAllGhosts(const VariableType var_type = VariableType::evolution){}; -}; - -#endif /* SIMPLEINTERPSOURCE_H_ */ diff --git a/Source/AMRInterpolator/SphericalExtraction.hpp b/Source/AMRInterpolator/SphericalExtraction.hpp index 348317cc6..8d8817eaf 100644 --- a/Source/AMRInterpolator/SphericalExtraction.hpp +++ b/Source/AMRInterpolator/SphericalExtraction.hpp @@ -85,21 +85,15 @@ class SphericalExtraction : public SurfaceExtraction const IntegrationMethod &a_method_phi = IntegrationMethod::trapezium, const bool a_broadcast_integral = false) { - // {x, y, z} - SphericalGeometry::UP_DIR up_dir = m_geom.get_up_dir(); - std::array dirs = {(up_dir + 1) % 3, (up_dir + 2) % 3, up_dir}; - std::array center; - for (int i = 0; i < 3; ++i) - center[i] = (dirs[i] < CH_SPACEDIM ? m_center[dirs[i]] : 0.); - - auto integrand_re = [dirs, center, &geom = m_geom, es, el, em, + auto integrand_re = [center = m_center, &geom = m_geom, es, el, em, &a_function](std::vector &a_data_here, - double r, double theta, double phi) { + double r, double theta, double phi) + { // note that spin_Y_lm requires the coordinates with the center // at the origin - double x = geom.get_grid_coord(dirs[0], r, theta, phi) - center[0]; - double y = geom.get_grid_coord(dirs[1], r, theta, phi) - center[1]; - double z = geom.get_grid_coord(dirs[2], r, theta, phi) - center[2]; + double x = geom.get_grid_coord(0, r, theta, phi) - center[0]; + double y = geom.get_grid_coord(1, r, theta, phi) - center[1]; + double z = geom.get_grid_coord(2, r, theta, phi) - center[2]; SphericalHarmonics::Y_lm_t Y_lm = SphericalHarmonics::spin_Y_lm(x, y, z, es, el, em); auto function_here = a_function(a_data_here, r, theta, phi); @@ -110,14 +104,15 @@ class SphericalExtraction : public SurfaceExtraction add_integrand(integrand_re, out_integrals.first, a_method_theta, a_method_phi, a_broadcast_integral); - auto integrand_im = [dirs, center, &geom = m_geom, es, el, em, + auto integrand_im = [center = m_center, &geom = m_geom, es, el, em, &a_function](std::vector &a_data_here, - double r, double theta, double phi) { + double r, double theta, double phi) + { // note that spin_Y_lm requires the coordinates with the center // at the origin - double x = geom.get_grid_coord(dirs[0], r, theta, phi) - center[0]; - double y = geom.get_grid_coord(dirs[0], r, theta, phi) - center[1]; - double z = geom.get_grid_coord(dirs[0], r, theta, phi) - center[2]; + double x = geom.get_grid_coord(0, r, theta, phi) - center[0]; + double y = geom.get_grid_coord(1, r, theta, phi) - center[1]; + double z = geom.get_grid_coord(2, r, theta, phi) - center[2]; SphericalHarmonics::Y_lm_t Y_lm = SphericalHarmonics::spin_Y_lm(x, y, z, es, el, em); auto function_here = a_function(a_data_here, r, theta, phi); diff --git a/Source/AMRInterpolator/SphericalGeometry.hpp b/Source/AMRInterpolator/SphericalGeometry.hpp index aa35d4ab6..2645fe8d9 100644 --- a/Source/AMRInterpolator/SphericalGeometry.hpp +++ b/Source/AMRInterpolator/SphericalGeometry.hpp @@ -7,12 +7,9 @@ #define SPHERICALGEOMETRY_HPP_ // Chombo includes -#include "AlwaysInline.hpp" #include "MayDay.H" // Other includes -#include "AlwaysInline.hpp" -#include "MayDay.H" #include #include #include @@ -25,95 +22,73 @@ //! u = theta, v = phi class SphericalGeometry { - public: - enum UP_DIR - { - X, - Y, - Z - }; + private: + std::array m_center; + public: SphericalGeometry(const std::array &a_center) : m_center(a_center) { -#if CH_SPACEDIM == 3 // for now force 'Z', could be an argument later on - m_up_dir = Z; -#elif CH_SPACEDIM == 2 // in Cartoon method, assume 'X' as 'z' axis - m_up_dir = X; -#endif } - ALWAYS_INLINE UP_DIR get_up_dir() const { return m_up_dir; } - ALWAYS_INLINE double get_domain_u_min() const { return 0.; } - ALWAYS_INLINE double get_domain_u_max() const { return M_PI; } - 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 theta - ALWAYS_INLINE double du(int a_num_points_theta) const + inline double du(int a_num_points_theta) const { - return (get_domain_u_max() - get_domain_u_min()) / - (double)(a_num_points_theta - 1); + return M_PI / (double)(a_num_points_theta - 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 theta coordinate associated to the theta/u index - ALWAYS_INLINE double u(int a_itheta, int a_num_points_theta) const + inline double u(int a_itheta, int a_num_points_theta) const { return a_itheta * du(a_num_points_theta); } //! 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; } - - ALWAYS_INLINE std::string param_name() const { return "r"; } - - ALWAYS_INLINE std::string u_name() const { return "theta"; } - - ALWAYS_INLINE std::string v_name() const { return "phi"; } - - //! returns the area element on a sphere with radius a_radius at the point - //! (a_theta, a_phi) - ALWAYS_INLINE double area_element(double a_radius, double a_theta, - double a_phi) const - { - return a_radius * a_radius * sin(a_theta); - } + 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, theta and phi. - ALWAYS_INLINE double get_grid_coord(int a_dir, double a_radius, - double a_theta, double a_phi = 0.) const + inline double get_grid_coord(int a_dir, double a_radius, double a_theta, + double a_phi) const { - double center = (a_dir < CH_SPACEDIM ? m_center[a_dir] : 0.); - switch ((a_dir + 2 - m_up_dir) % 3) + switch (a_dir) { case (0): - return center + a_radius * sin(a_theta) * cos(a_phi); + return m_center[0] + a_radius * sin(a_theta) * cos(a_phi); case (1): - return center + a_radius * sin(a_theta) * sin(a_phi); + return m_center[1] + a_radius * sin(a_theta) * sin(a_phi); case (2): - return center + a_radius * cos(a_theta); + return m_center[2] + a_radius * cos(a_theta); default: MayDay::Error("SphericalGeometry: Direction not supported"); } } - protected: - std::array m_center; - UP_DIR m_up_dir; + //! returns the area element on a sphere with radius a_radius at the point + //! (a_theta, a_phi) + inline double area_element(double a_radius, double a_theta, + double a_phi) const + { + return a_radius * a_radius * sin(a_theta); + } + + inline std::string param_name() const { return "r"; } + + inline std::string u_name() const { return "theta"; } + + inline std::string v_name() const { return "phi"; } }; #endif /* SPHERICALGEOMETRY_HPP_ */ diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index eb444f130..e060225a7 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -13,6 +13,7 @@ // Other inclues #include "AMRInterpolator.hpp" #include "DimensionDefinitions.hpp" +#include "FilesystemTools.hpp" #include "IntegrationMethod.hpp" #include "InterpolationQuery.hpp" #include "Lagrange.hpp" @@ -50,7 +51,8 @@ template class SurfaceExtraction //!< extraction for each surface bool write_extraction; //!< whether or not to write the extracted data - std::string extraction_prefix; + std::string data_path, integral_file_prefix; + std::string extraction_path, extraction_file_prefix; int min_extraction_level() { @@ -94,18 +96,11 @@ template class SurfaceExtraction //! returns the flattened index for m_interp_data and m_interp_coords //! associated to given surface, u and v indices -#if CH_SPACEDIM == 3 int index(int a_isurface, int a_iu, int a_iv) const { return a_isurface * m_params.num_points_u * m_params.num_points_v + a_iu * m_params.num_points_v + a_iv; } -#elif CH_SPACEDIM == 2 - int index(int a_isurface, int a_iu) const - { - return a_isurface * m_params.num_points_u + a_iu; - } -#endif public: //! Normal constructor which requires vars to be added after construction diff --git a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp index ccfe05ffb..9b97a436b 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -19,18 +19,28 @@ SurfaceExtraction::SurfaceExtraction( : m_geom(a_geom), m_params(a_params), m_dt(a_dt), m_time(a_time), m_first_step(a_first_step), m_restart_time(a_restart_time), m_num_interp_points((procID() == 0) - ? m_params.num_surfaces * m_params.num_points_u -#if CH_SPACEDIM == 3 - * m_params.num_points_v -#endif + ? m_params.num_surfaces * m_params.num_points_u * + m_params.num_points_v : 0), m_du(m_geom.du(m_params.num_points_u)), m_dv(m_geom.dv(m_params.num_points_v)), m_done_extraction(false) { + // check folders only in first two timesteps + // (or at m_first_step if this is not the first two timesteps) + if (m_time < m_restart_time + 1.5 * m_dt || m_first_step) + { + if (!FilesystemTools::directory_exists(m_params.data_path)) + FilesystemTools::mkdir_recursive(m_params.data_path); + + if (m_params.write_extraction && + !FilesystemTools::directory_exists(m_params.extraction_path)) + FilesystemTools::mkdir_recursive(m_params.extraction_path); + } + // only interp points on rank 0 if (procID() == 0) { - FOR1(idir) { m_interp_coords[idir].resize(m_num_interp_points); } + FOR(idir) { m_interp_coords[idir].resize(m_num_interp_points); } for (int isurface = 0; isurface < m_params.num_surfaces; ++isurface) { @@ -39,25 +49,16 @@ SurfaceExtraction::SurfaceExtraction( for (int iu = 0; iu < m_params.num_points_u; ++iu) { double u = m_geom.u(iu, m_params.num_points_u); -#if CH_SPACEDIM == 3 for (int iv = 0; iv < m_params.num_points_v; ++iv) { double v = m_geom.v(iv, m_params.num_points_v); - FOR1(idir) + FOR(idir) { int idx = index(isurface, iu, iv); m_interp_coords[idir][idx] = m_geom.get_grid_coord( idir, surface_param_value, u, v); } } -#elif CH_SPACEDIM == 2 - FOR1(idir) - { - int idx = index(isurface, iu); - m_interp_coords[idir][idx] = - m_geom.get_grid_coord(idir, surface_param_value, u); - } -#endif } } } @@ -146,7 +147,7 @@ void SurfaceExtraction::extract( } // m_num_interp_points is 0 on ranks > 0 InterpolationQuery query(m_num_interp_points); - FOR1(idir) { query.setCoords(idir, m_interp_coords[idir].data()); } + FOR(idir) { query.setCoords(idir, m_interp_coords[idir].data()); } for (int ivar = 0; ivar < m_vars.size(); ++ivar) { // note the difference in order between the m_vars tuple in this class @@ -232,9 +233,9 @@ void SurfaceExtraction::add_var_integrand( const bool a_broadcast_integral) { CH_assert(a_var >= 0 && a_var < m_vars.size()); - integrand_t var_integrand = [var = a_var](std::vector &data, double, - double, - double) { return data[var]; }; + integrand_t var_integrand = + [var = a_var](std::vector &data, double, double, double) + { return data[var]; }; add_integrand(var_integrand, out_integrals, a_method_u, a_method_v, a_broadcast_integral); } @@ -265,11 +266,7 @@ void SurfaceExtraction::integrate() for (int ivar = 0; ivar < m_vars.size(); ++ivar) { data_here[ivar] = -#if CH_SPACEDIM == 3 m_interp_data[ivar][index(isurface, iu, iv)]; -#elif CH_SPACEDIM == 2 - m_interp_data[ivar][index(isurface, iu)]; -#endif } for (int iintegral = 0; iintegral < num_integrals; ++iintegral) @@ -345,7 +342,7 @@ void SurfaceExtraction::write_extraction( CH_assert(m_done_extraction); if (procID() == 0) { - SmallDataIO extraction_file(m_params.extraction_prefix + a_file_prefix, + SmallDataIO extraction_file(m_params.extraction_path + a_file_prefix, m_dt, m_time, m_restart_time, SmallDataIO::NEW, m_first_step); @@ -390,26 +387,17 @@ void SurfaceExtraction::write_extraction( for (int iu = 0; iu < m_params.num_points_u; ++iu) { double u = m_geom.u(iu, m_params.num_points_u); -#if CH_SPACEDIM == 3 for (int iv = 0; iv < m_params.num_points_v; ++iv) { double v = m_geom.v(iv, m_params.num_points_v); int idx = index(isurface, iu, iv); -#elif CH_SPACEDIM == 2 - { - int idx = index(isurface, iu); -#endif std::vector data(m_vars.size()); for (int ivar = 0; ivar < m_vars.size(); ++ivar) { data[ivar] = m_interp_data[ivar][idx]; } -#if CH_SPACEDIM == 3 extraction_file.write_data_line(data, {u, v}); -#elif CH_SPACEDIM == 2 - extraction_file.write_data_line(data, {u}); -#endif } } extraction_file.line_break(); @@ -440,8 +428,8 @@ void SurfaceExtraction::write_integrals( CH_assert(vect.size() == m_params.num_surfaces); } // open file for writing - SmallDataIO integral_file(m_params.extraction_prefix + a_filename, m_dt, - m_time, m_restart_time, SmallDataIO::APPEND, + SmallDataIO integral_file(m_params.data_path + a_filename, m_dt, m_time, + m_restart_time, SmallDataIO::APPEND, m_first_step); // remove any duplicate data if this is a restart diff --git a/Source/ApparentHorizonFinder/AHData.hpp b/Source/ApparentHorizonFinder/AHData.hpp deleted file mode 100644 index 8585bfa2e..000000000 --- a/Source/ApparentHorizonFinder/AHData.hpp +++ /dev/null @@ -1,88 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHDATA_HPP_ -#define _AHDATA_HPP_ - -#include "Derivative.hpp" -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists -#include "InterpolationQuery.hpp" -#include "Tensor.hpp" - -#include // memset - -// The d prefix refers to partial derivatives wrt Cartesian coordinates - -template struct AHData -{ - // vars user might want to export - std::map vars; - std::map> d1; - std::map> d2; - - void set_vars(InterpolationQuery &query, int var, key k, VariableType type, - int total) - { - vars[k].resize(total); - query.addComp(var, &vars[k][0], Derivative::LOCAL, type); - } - void set_d1(InterpolationQuery &query, int var, key k, VariableType type, - int total) - { - for (int j = 0; j < CH_SPACEDIM; ++j) - d1[k][j].resize(total); - - query.addComp(var, &d1[k][0][0], Derivative::dx, type) - .addComp(var, &d1[k][1][0], Derivative::dy, type); -#if CH_SPACEDIM == 3 - query.addComp(var, &d1[k][2][0], Derivative::dz, type); -#endif - } - void set_d2(InterpolationQuery &query, int var, key k, VariableType type, - int total) - { - for (int j = 0; j < CH_SPACEDIM * (CH_SPACEDIM + 1) / 2; ++j) - d2[k][j].resize(total); - - query.addComp(var, &d2[k][0][0], Derivative::dxdx, type) - .addComp(var, &d2[k][1][0], Derivative::dxdy, type); -#if CH_SPACEDIM == 3 - query.addComp(var, &d2[k][2][0], Derivative::dxdz, type) - .addComp(var, &d2[k][3][0], Derivative::dydy, type) - .addComp(var, &d2[k][4][0], Derivative::dydz, type) - .addComp(var, &d2[k][5][0], Derivative::dzdz, type); -#elif CH_SPACEDIM == 2 - query.addComp(var, &d2[k][2][0], Derivative::dydy, type); -#endif - } -}; - -template -AHData -get_AHData_idx(int idx, const AHData> &a_data) -{ - AHData data; - - for (auto &vec : a_data.vars) - data.vars[vec.first] = (vec.second[idx]); - - for (auto &vec : a_data.d1) - { - data.d1[vec.first] = {0.}; - for (int i = 0; i < CH_SPACEDIM; ++i) - data.d1[vec.first][i] = vec.second[i][idx]; - } - - for (auto &vec : a_data.d2) - { - data.d2[vec.first] = {0.}; - for (int i = 0; i < CH_SPACEDIM * (CH_SPACEDIM + 1) / 2; ++i) - data.d2[vec.first][i] = vec.second[i][idx]; - } - - return data; -} - -#endif /* _AHDATA_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHDeriv.hpp b/Source/ApparentHorizonFinder/AHDeriv.hpp deleted file mode 100644 index 4c0fe5883..000000000 --- a/Source/ApparentHorizonFinder/AHDeriv.hpp +++ /dev/null @@ -1,48 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHDERIV_HPP_ -#define _AHDERIV_HPP_ - -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists - -#include // memset - -#define DWIDTH 5 -#define DDWIDTH 6 - -//! class to store 1st and 2nd derivatives of 'F' -struct AHDeriv -{ - double duF; - double duduF; -#if CH_SPACEDIM == 3 - double dvF; - double dvdvF; - double dudvF; -#endif - - int du_stencil_start; - int dudu_stencil_start; -#if CH_SPACEDIM == 3 - int dv_stencil_start; - int dvdv_stencil_start; -#endif - - double du_weights[DWIDTH]; - double dudu_weights[DDWIDTH]; -#if CH_SPACEDIM == 3 - double dv_weights[DWIDTH]; - double dvdv_weights[DDWIDTH]; -#endif - - AHDeriv() - { - // force all (double) elements of AHDeriv to be 0 - memset(this, 0, sizeof(AHDeriv)); - } -}; - -#endif /* _AHDERIV_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHFinder.cpp b/Source/ApparentHorizonFinder/AHFinder.cpp deleted file mode 100644 index 4b396afa8..000000000 --- a/Source/ApparentHorizonFinder/AHFinder.cpp +++ /dev/null @@ -1,634 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifdef USE_AHFINDER - -#include "AHFinder.hpp" - -#include "AHInterpolation.hpp" -#include "ApparentHorizon.hpp" - -AHFinder::~AHFinder() -{ - // destroy horizon pointers and finalize PETSc - for (auto &ah : m_apparent_horizons) - delete ah; - PETSc_finalize(); -} - -int AHFinder::add_ah(const AHSurfaceGeometry &a_coord_system, - double a_initial_guess, const AHFinder::params &a_params, - const std::string &a_stats, const std::string &a_coords, - bool solve_first_step) -{ - return add_ah(a_coord_system, a_initial_guess, a_params, - AHFunction::params(), a_stats, a_coords, solve_first_step); -} - -int AHFinder::add_ah(const AHSurfaceGeometry &a_coord_system, - double a_initial_guess, const AHFinder::params &a_params, - const typename AHFunction::params &a_func_params, - const std::string &a_stats, const std::string &a_coords, - bool solve_first_step) -{ - if (!AHFinder::m_initialized) - AHFinder::PETSc_initialize(a_params.num_ranks); - - // determine how many AH there are already - const int num_ah = m_apparent_horizons.size(); - - AHInterpolation interp(a_coord_system, - m_interpolator); - - m_apparent_horizons.push_back( - new ApparentHorizon( - interp, a_initial_guess, a_params, a_func_params, - a_stats + std::to_string(num_ah + 1), - a_coords + std::to_string(num_ah + 1) + "_", solve_first_step)); - m_merger_pairs.push_back({-1, -1}); - - return num_ah; -} - -int AHFinder::add_ah_merger(int ah1, int ah2, const params &a_params) -{ - CH_assert(a_params.track_center); // center tracking must be on for mergers - - int num_ah = m_apparent_horizons.size(); - CH_assert(ah1 >= 0 && ah1 < num_ah); - CH_assert(ah2 >= 0 && ah2 < num_ah); - CH_assert(ah1 != ah2); - - double initial_guess_merger; - std::array origin_merger; - bool do_solve = solve_merger(ah1, ah2, initial_guess_merger, origin_merger); - - // assume same coordiante system between merging horizons - AHSurfaceGeometry coord_system = - m_apparent_horizons[ah1]->get_ah_interp().get_coord_system(); - coord_system.set_origin(origin_merger); - - std::string stats = m_apparent_horizons[ah1]->m_stats; - std::string coords = m_apparent_horizons[ah1]->m_coords; - // remove number of AH1 from file names - stats = stats.substr(0, stats.size() - std::to_string(ah1).size()); - coords = coords.substr(0, coords.size() - std::to_string(ah1).size() - 1); - - auto &function_to_optimize_params = m_apparent_horizons[ah1]->m_func_params; - - int num = add_ah(coord_system, initial_guess_merger, a_params, - function_to_optimize_params, stats, coords, do_solve); - m_merger_pairs[num] = {ah1, ah2}; - - return num; -} - -void AHFinder::solve(double a_dt, double a_time, double a_restart_time) -{ - CH_TIME("AHFinder::solve"); - - enum - { - NOT_SOLVED, - TOO_FAR, - SKIPPED, - SOLVED - }; - - // in case this was called at specificPostTimeStep at t=0 due to - // MultiLevelTask (solve already happened when AH was created) - if (a_time == 0.) - return; - - // solve if it has never converged before OR if has already converged in the - // past this means it will stop solving as soon as it stops converging but - // keeps trying (with same initial guess) if never found one - // (e.g. in the BinaryBH case, AH1 and AH2 stop solving once they - // disappear, leaving the merged AH alone from then on) - // (e.g. in the ScalarField case, only after a few time steps does an AH - // form, and from then on it keeps solving forever) - - const int num_ah = m_apparent_horizons.size(); - - // skip if AH 0 is not supposed to be solved - // assumes all AHs have same 'solve_interval' - if (num_ah == 0 || !m_apparent_horizons[0]->do_solve(a_dt, a_time)) - return; - - // the mess below is done so that AH are searched for iteratively - // taking into account merger priorities - // (a merger will only be solved if it's 'parents' have already been solved) - int ahs_solved = 0; - // sets whether or not AH[i] has already been solved - // 0 - not solved; 1 - too far away; 2 - skipped; 3 - solved - int *ah_solved = - new int[num_ah]; // with '-g' compilation flag was complaining that int - // ah_solved[num_ah] was invalid - for (int i = 0; i < num_ah; ++i) - ah_solved[i] = NOT_SOLVED; - - int max_attempts = 100, attempt = 0; - while (ahs_solved < num_ah) - { - if (max_attempts < attempt++) - { - pout() << "Reached maximum number of attempts in AH solving loop." - << std::endl; - break; - } - - for (int i = 0; i < num_ah; ++i) - { - if (ah_solved[i] != NOT_SOLVED) - continue; // continue if already solved - - // if is a merger and has not yet converged - if (m_merger_pairs[i].first >= 0 && - !m_apparent_horizons[i]->has_been_found()) - { - auto pair = m_merger_pairs[i]; - - if (ah_solved[pair.first] == NOT_SOLVED || - ah_solved[pair.second] == NOT_SOLVED) - continue; // continue if one of the 'parents' was not - // dealt with yet - - if (!(m_apparent_horizons[pair.first]->has_been_found() && - m_apparent_horizons[pair.second]->has_been_found())) - { - if (m_apparent_horizons[i]->m_params.verbose > NONE) - { - pout() << "Skipping BH #" << i << std::endl; - } - ++ahs_solved; - ah_solved[i] = SKIPPED; - continue; // skip if one of the parents was never found yet - } - - if (ah_solved[pair.first] == TOO_FAR || - ah_solved[pair.second] == TOO_FAR) - { // if one of the 'parents' is a merger too far away to be - // solved, skip children too - if (m_apparent_horizons[i]->m_params.verbose > NONE) - { - pout() << "Skipping BH #" << i << std::endl; - } - ++ahs_solved; - ah_solved[i] = TOO_FAR; - continue; - } - - double initial_guess_merger; - std::array center_merger; - bool do_solve = - solve_merger(pair.first, pair.second, initial_guess_merger, - center_merger); - - // if distance is too large, ignore this one and move on - if (!do_solve) - { - // already printed 'Skipping' when calling 'solve_merger' - ++ahs_solved; - ah_solved[i] = TOO_FAR; - continue; - } - - // if it already converged, don't do anything; - // no need to change the initial guess, as it is the same since - // the beginning; - // if one was skipped don't update the center (because it will - // not be central anymore) - if (!m_apparent_horizons[i]->get_converged() && - ah_solved[pair.first] == SOLVED && - ah_solved[pair.second] == SOLVED) - { - m_apparent_horizons[i]->set_origin(center_merger); - // m_apparent_horizons[i]->set_initial_guess( - // initial_guess_merger); // should still be the same - } - } - - // solve if it converged last time or if is has never been found in - // the past - if (m_apparent_horizons[i]->good_to_go(a_dt, a_time)) - { - if (m_apparent_horizons[i]->m_params.verbose > NONE) - pout() << "Solving AH #" << i << std::endl; - - m_apparent_horizons[i]->solve(a_dt, a_time, a_restart_time); - - ++ahs_solved; - ah_solved[i] = SOLVED; - } - else - { - if (m_apparent_horizons[i]->m_params.verbose > NONE) - { - pout() << "Skipping BH #" << i << ". Not good to go." - << std::endl; - } - ++ahs_solved; - ah_solved[i] = SKIPPED; - } - } - } - - delete ah_solved; -} - -bool AHFinder::need_diagnostics(double a_dt, double a_time) const -{ - bool out = false; - for (auto &ah : m_apparent_horizons) - out |= - ah->m_params.extra_contain_diagnostic && ah->do_print(a_dt, a_time); - return out; -} - -bool AHFinder::solve_merger(int ah1, int ah2, double &initial_guess_merger, - std::array ¢er_merger) -{ - // SKIP if 'parents' not yet close enough - auto AH1 = m_apparent_horizons[ah1]; - auto AH2 = m_apparent_horizons[ah2]; - auto center1 = AH1->get_center(); - auto center2 = AH2->get_center(); - - auto initial_guess1 = AH1->get_initial_guess(); - auto initial_guess2 = AH2->get_initial_guess(); - - if (m_merger_pairs[ah1].first >= 0) - { - double merger_pre_factor = - std::max(m_apparent_horizons[m_merger_pairs[ah1].first] - ->m_params.merger_pre_factor, - m_apparent_horizons[m_merger_pairs[ah1].second] - ->m_params.merger_pre_factor); - if (merger_pre_factor > - 0.) // undo the factor is this a merger from a merger - initial_guess1 /= (merger_pre_factor * 4.); - } - if (m_merger_pairs[ah2].first >= 0) - { - double merger_pre_factor = - std::max(m_apparent_horizons[m_merger_pairs[ah2].first] - ->m_params.merger_pre_factor, - m_apparent_horizons[m_merger_pairs[ah2].second] - ->m_params.merger_pre_factor); - if (merger_pre_factor > - 0.) // undo the factor is this a merger from a merger - initial_guess2 /= (merger_pre_factor * 4.); - } - - double initial_guess_sum = (initial_guess1 + initial_guess2); - - // some tests revealed it was about 1.2 * initial_guess_sum, but 1.5 to - // ensure we don't catch the inner AH - double merger_pre_factor = std::max(AH1->m_params.merger_pre_factor, - AH2->m_params.merger_pre_factor); - initial_guess_merger = merger_pre_factor * 4. * initial_guess_sum; - - // update center of merged, otherwise it does it by - // itself in solve - FOR1(a) - { - if (!AH1->get_ah_interp().get_interpolator()->get_boundary_reflective( - Side::Lo, a)) - center_merger[a] = (center1[a] + center2[a]) / 2.; - else - center_merger[a] = 0.; - } - - // check if distance is too big - double distance = - sqrt((center1[0] - center2[0]) * (center1[0] - center2[0]) + - (center1[1] - center2[1]) * (center1[1] - center2[1]) -#if CH_SPACEDIM == 3 - + (center1[2] - center2[2]) * (center1[2] - center2[2]) -#endif - ); - - double merger_search_factor = std::max(AH1->m_params.merger_search_factor, - AH2->m_params.merger_search_factor); - - double min_distance = merger_search_factor * 4. * initial_guess_sum; - - bool do_solve = false; - - if (AH1->has_been_found() && AH2->has_been_found()) - { - if (AH1->get_converged() && AH2->get_converged()) - { - do_solve = merger_search_factor <= 0. || distance <= min_distance; - - if (do_solve) - // make sure twice the radius of the guess is bigger than AH - // distance - CH_assert(min_distance <= 2. * initial_guess_merger); - - if (AH1->m_params.verbose > NONE) - { - pout() << "BHs #" << ah1 << " and #" << ah2 - << " at distance = " << distance; - // if distance is too large, ignore this one and move on - if (!do_solve) - pout() << " > minimum distance = " << min_distance - << ". Skipping solve for merger..."; - - pout() << std::endl; - } - } - else - { - // if AH1 or AH2 was lost, probably merger is already around - do_solve = true; - } - } - else - { - if (AH1->m_params.verbose > NONE) - { - pout() << "Original BH not yet found. Skipping solve for merger." - << std::endl; - } - } - - return do_solve; -} - -void AHFinder::set_origins( - const std::vector> &origins, - bool includes_mergers) -{ - const int num_ah = m_apparent_horizons.size(); - - int non_merger_counter = 0; // for 'includes_mergers == false' - for (int i = 0; i < num_ah; ++i) - { - if (includes_mergers) - { - m_apparent_horizons[i]->set_origin(origins[i]); - } - else - { - // if it is a merger that has not yet been found - if (m_merger_pairs[i].first >= 0) - { - if (!m_apparent_horizons[i]->has_been_found()) - { - auto origin1 = m_apparent_horizons[m_merger_pairs[i].first] - ->get_origin(); - auto origin2 = m_apparent_horizons[m_merger_pairs[i].second] - ->get_origin(); - std::array origin; - FOR1(a) { origin[a] = (origin1[a] + origin2[a]) / 2.; } - m_apparent_horizons[i]->set_origin(origin); - } - } - else - { - m_apparent_horizons[i]->set_origin(origins[non_merger_counter]); - ++non_merger_counter; - } - } - } -} - -void AHFinder::params::read_params(GRParmParse &pp, const ChomboParameters &a_p) -{ - pp.load("AH_num_ranks", num_ranks, 0); // 0 means "all" - pp.load("AH_num_points_u", num_points_u); -#if CH_SPACEDIM == 3 - pp.load("AH_num_points_v", num_points_v); -#endif - pp.load("AH_solve_interval", solve_interval, 1); - pp.load("AH_print_interval", print_interval, 1); - pp.load("AH_track_center", track_center, true); - pp.load("AH_predict_origin", predict_origin, track_center); - // can't predict if center is not being tracked - CH_assert(!(predict_origin && !track_center)); - - pp.load("AH_level_to_run", level_to_run, 0); - CH_assert(level_to_run <= a_p.max_level && - level_to_run > -(a_p.max_level + 1)); - if (level_to_run < 0) // if negative, count backwards - level_to_run += a_p.max_level + 1; - - pp.load("AH_start_time", start_time, 0.0); - pp.load("AH_give_up_time", give_up_time, -1.0); - - pp.load("AH_merger_search_factor", merger_search_factor, 1.); - pp.load("AH_merger_pre_factor", merger_pre_factor, 1.); - - pp.load("AH_allow_re_attempt", allow_re_attempt, false); - pp.load("AH_max_fails_after_lost", max_fails_after_lost, 0); - pp.load("AH_stop_if_max_fails", stop_if_max_fails, false); - - pp.load("AH_verbose", verbose, (int)AHFinder::MIN); - pp.load("AH_print_geometry_data", print_geometry_data, false); - pp.load("AH_re_solve_at_restart", re_solve_at_restart, false); - - // sanity checks - CH_assert(solve_interval > 0 && print_interval > 0); - CH_assert(level_to_run >= 0 && level_to_run <= a_p.max_level); - - // if box division ends up with size less than 3, PETSc will - // complain (this only gives an estimate of the box side) - int size = 1; -#if CH_MPI - MPI_Comm_size(Chombo_MPI::comm, &size); -#endif - size = std::min(num_ranks, size); -#if CH_SPACEDIM == 3 - CH_assert(num_points_u > 0 && num_points_v > 0); - CH_assert(num_points_u / sqrt(size) >= 3); // make sure for size 'u' - CH_assert(num_points_v / sqrt(size) >= 3); // make sure for size 'v' -#elif CH_SPACEDIM == 2 - CH_assert(num_points_u > 0); - CH_assert(num_points_u / size >= 3); // make sure for size 'u' -#endif - - // load vars to write to coord files - num_extra_vars = 0; - extra_contain_diagnostic = 0; - - int AH_num_write_vars; - pp.load("AH_num_write_vars", AH_num_write_vars, 0); - if (AH_num_write_vars > 0) - { - std::vector AH_write_var_names(AH_num_write_vars, ""); - pp.load("AH_write_vars", AH_write_var_names, AH_num_write_vars); - for (const std::string &full_name : AH_write_var_names) - { - std::string var_name = full_name; - - // variable names might start with "d1_" or "d2_" to indicate - // the user wants derivatives - int der_type = 0; - std::string der = var_name.substr(0, 3); - if (der == "d1_") - { - der_type = 1; - var_name = var_name.substr(3); - } - else if (der == "d2_") - { - der_type = 2; - var_name = var_name.substr(3); - } - - // first assume write_var is a normal evolution var - int var = UserVariables::variable_name_to_enum(var_name); - VariableType var_type = VariableType::evolution; - if (var < 0) - { - // if not an evolution var check if it's a diagnostic var - var = DiagnosticVariables::variable_name_to_enum(var_name); - if (var < 0) - { - // it's neither :( - pout() << "Variable with name " << var_name - << " not found.\n"; - } - else - { - var_type = VariableType::diagnostic; - ++extra_contain_diagnostic; - } - } - if (var >= 0) - { - extra_vars[full_name] = - std::tuple(var, var_type, der_type); - if (der_type == 0) - num_extra_vars += 1; - else if (der_type == 1) - num_extra_vars += CH_SPACEDIM; - else // if (der_type == 2) - num_extra_vars += CH_SPACEDIM * (CH_SPACEDIM + 1) / 2; - } - } - } -} - -///////////////////////////////////////////////////////// -// PETSc control methods -///////////////////////////////////////////////////////// - -void AHFinder::set_num_ranks(int a_num_ranks) -{ - if (m_initialized) - return; - -#ifdef CH_MPI - if (m_mpi_group != MPI_GROUP_NULL) - MPI_Group_free(&m_mpi_group); - - if (m_mpi_comm != MPI_COMM_NULL) - MPI_Comm_free(&m_mpi_comm); - - if (a_num_ranks > 0) // otherwise use the whole Chombo communicator - { - // don't make more ranks than there are processes - int size; - MPI_Comm_size(Chombo_MPI::comm, &size); - a_num_ranks = std::min(a_num_ranks, size); - - int rank; - MPI_Comm_rank(Chombo_MPI::comm, &rank); - if (rank == 0) - std::cout << "Using PETSc with " << a_num_ranks << " ranks" - << std::endl; - - MPI_Group MPI_GROUP_WORLD; - MPI_Comm_group(Chombo_MPI::comm, &MPI_GROUP_WORLD); - - // (TF): for things like 'SmallDataIO', rank 0 of Chombo_MPI::comm - // is the one that prints (not rank 0 from PETSC_COMM_WORLD) so pay - // attention before taking rank 0 of Chombo out of PETSc. As of - // 26/05/2019, I tested with 'int range[3] = { 1, a_num_ranks , 1 }' - // and all the AH stuff was compatible with doing so - // Not true anymore (19/08/2020) due to how - // ApparentHorizon::write_coords_file uses rank 0 of PETSc to transfer - // data and to write, but SmallDataIO uses rank 0 of Chombo to write. - int range[3] = {0, a_num_ranks - 1, 1}; - MPI_Group_range_incl(MPI_GROUP_WORLD, 1, - reinterpret_cast(&range), &m_mpi_group); - MPI_Group_free(&MPI_GROUP_WORLD); - - MPI_Comm_create(Chombo_MPI::comm, m_mpi_group, &m_mpi_comm); - } -#endif -} - -//! initialize PETSc and its MPI sub-communicator -PetscErrorCode AHFinder::PETSc_initialize(int a_num_ranks) -{ - set_num_ranks(a_num_ranks); - - PetscErrorCode err = 0; - if (!m_initialized) - { -#ifdef CH_MPI - if (m_mpi_comm != MPI_COMM_NULL) - PETSC_COMM_WORLD = m_mpi_comm; - else // use Chombo communicator if no 'set_num_ranks' was called (or - // if called with 'a_num_ranks'<=0) - PETSC_COMM_WORLD = Chombo_MPI::comm; -#endif - - if (AHFinder::is_rank_active()) - err = PetscInitializeNoArguments(); - - if (!err) - m_initialized = true; - } - return err; -} - -//! finalize PETSc -PetscErrorCode AHFinder::PETSc_finalize() -{ - PetscErrorCode err = 0; - if (m_initialized) - { - if (AHFinder::is_rank_active()) - err = PetscFinalize(); - - if (!err) - m_initialized = false; - } - return err; -} - -//! true if part of PETSc MPI sub-communicator -bool AHFinder::is_rank_active() -{ -#ifdef CH_MPI - if (m_mpi_group == MPI_GROUP_NULL) - return true; - else - { - int rank; - MPI_Group_rank(m_mpi_group, &rank); - return rank != MPI_UNDEFINED; - } -#else - return true; -#endif -} - -///////////////////////////////////////////////////////// -// initialize some "static" variables of AHFinder class -///////////////////////////////////////////////////////// - -bool AHFinder::m_initialized = false; - -#ifdef CH_MPI -MPI_Group AHFinder::m_mpi_group = MPI_GROUP_NULL; -MPI_Comm AHFinder::m_mpi_comm = MPI_COMM_NULL; -#endif - -#endif diff --git a/Source/ApparentHorizonFinder/AHFinder.hpp b/Source/ApparentHorizonFinder/AHFinder.hpp deleted file mode 100644 index 6c48d4305..000000000 --- a/Source/ApparentHorizonFinder/AHFinder.hpp +++ /dev/null @@ -1,265 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHFINDER_HPP_ -#define _AHFINDER_HPP_ - -// Chombo includes -#include "SPMD.H" //Chombo_MPI::comm - -// Other includes: - -// included so that user can re-define default AHSurfaceGeometry or AHFunction -#include "UserVariables.hpp" - -// define SurfaceGeometry of AHFinder -#ifndef AHSurfaceGeometry -#include "AHSphericalGeometry.hpp" -#define AHSurfaceGeometry AHSphericalGeometry -#endif - -// default to expansion -#ifndef AHFunction -#include "AHFunctions.hpp" -#define AHFunction ExpansionFunction -#endif - -#ifdef CH_MPI -#include -#endif - -#include - -#include "AHData.hpp" -#include "AHDeriv.hpp" -#include "AHGeometryData.hpp" -#include "AMRInterpolator.hpp" -#include "BoundaryConditions.hpp" -#include "ChomboParameters.hpp" -#include "Lagrange.hpp" -#include "VariableType.hpp" - -// Chombo namespace -#include "UsingNamespace.H" - -// declare classes, define them later -template class AHInterpolation; -template class ApparentHorizon; - -/* -TF: -A general note on the radius and the value of chi on the AH surface. -For a Kerr BH initial data with dimensionless spin 's' and mass 'M': - r_AH = M / 4 (1 + sqrt{1 - s^2}) - chi_AH = [(1-s^2)^(1/6) / 16 * ((1 + sqrt{1-s^2})/2)^{2/3} for theta=0, - (1-s^2)^{1/6} / 16 * ((1 + sqrt{1-s^2})/2)^{1/3} for theta=pi/2] - -After puncture gauge settled: - r_AH ~ M * (0.275 + 0.76 * sqrt{1-s^2}) - chi_AH ~ 0.26 * sqrt{1-s^2} - -Take it with a grain of salt, but these results are from a numerical fit I did -for numerical Kerr BH runs, approximately settled at around t~50M, from spin=0 -to spin=0.99. You can use it generically for other things, like realizing that -for spin=0 you can plot in VisIt the contour chi=0.26 to get the AH, or that for -spin=0 the AH goes from r=M/2 to r~M - -Binaries: -After some inspiral cases I've been through, I realized that merger happened -when the distance between the punctures was about (roughtly!) ~[0.5,2]*(M1+M2) -and with a radius of ~[0.7,1.1]*(M1+M2) (wobbly of course, and bigger values for -smaller final spin). It makes sense to start looking for the merger before this, -but remember that for PETSc diverging is significantly slower than converging, -so the closer the better. Hence, as default, start looking for a merger when the -separation is about ~<2*(M1+M2) and with a initial guess of ~2*(M1+M2) (twice as -big as what we'll probably get to make sure it converges to the outer AH; even -then it might converge to the inner one, but it is much easer for the AH to -converge if the initial guess is bigger than the actual radius, but the opposite -is more senstitive). -*/ - -//! Class to manage AHs and its mergers + control PETSc MPI sub-communicator -class AHFinder -{ - public: - struct params - { - int num_ranks; //!< number of ranks for PETSc sub-communicator (default - //!< 0, which is 'all') - int num_points_u; //!< number of points for 2D coordinate grid -#if CH_SPACEDIM == 3 - int num_points_v; //!< number of points for 2D coordinate grid -#endif - int solve_interval; //!< same as checkpoint_interval, for - //!< ApparentHorizon::solve (default 1) - int print_interval; //!< same as solve_interval, but for prints (default - //!< 1) - bool track_center; //!< whether or not to update the center - //!< (set to false if you know it won't move) - //!< (default true) - bool predict_origin; //!< whether or not to estimate where the next - //!< origin will be at (default = track_center) - - int level_to_run; // if negative, it will count backwards (e.g. -1 is - // 'last level') (default 0) - - double start_time; //!< time after which AH can start (default 0.) - //!< Useful for ScalarField collapse - double give_up_time; //!< stop if at this time nothing was found - //!< (<0 to never, which is default) - //!< Useful for ScalarField collapse - - //! mergers will be searched when distance between 'parent' BHs is - //! distance < merger_search_factor * 4. * (AH_initial_guess_1 + - //! AH_initial_guess_2) should be roughly '2M=2(m1+m2)' for initial - //! guess at m/2 (set to non-positive to 'always search') - double merger_search_factor; // see note above (default is 1) - //! initial guess for merger is 'merger_pre_factor * 4. * - //! (AH_initial_guess_1 + AH_initial_guess_2)' - //! set to somethig bigger to avoid finding the inner AH - double merger_pre_factor; // see note above (default to 1.) - - bool allow_re_attempt; //!< re-attempt with initial guess if - //!< previous convergence failed (default false) - int max_fails_after_lost; //!< number of time steps to try again after - //!< (-1 to never) the AH was lost - //!< (default is 0) - - int verbose; //!< print information during execution (default is 1) - - bool print_geometry_data; //!< print metric and extrinsic - //!< curvature of the AHs (default false) - - bool re_solve_at_restart; //!< whether or not to re-run even if AH - //!< already exists (useful in order to be - //!< able to provide an initial guess and - //!< re-run the AHFinder on that time step) - //!< (default false) - - bool stop_if_max_fails; //! breaks the run if AH doesn't converge - //! 'max_fails_after_lost' times or if - //! 'give_up_time' is reached without - //! convergence (default is 'false') - - std::map> - extra_vars; //! extra vars to write to coordinates file () - int num_extra_vars; // total number of extra vars (!=extra_vars.size() - // as derivative count for multiple vars) - - int extra_contain_diagnostic; // not a parameter (set internally); - // counts how many - - void read_params(GRParmParse &pp, const ChomboParameters &a_p); - }; - - enum verbose_level - { - NONE, - MIN, // minimal - SOME, // a bit technical - MORE // debug - }; - - private: - //! if this AH is supposed to track the formation of a merger - //! this pair indicates the indices of the 2 AHs in the vector - //! 'm_apparent_horizons' - std::vector> m_merger_pairs; - AMRInterpolator> *m_interpolator; //!< The interpolator pointer - - std::vector *> - m_apparent_horizons; //!< public in case user wants to solve by himself - - public: - AHFinder(){}; - ~AHFinder(); - - ALWAYS_INLINE ApparentHorizon * - get(unsigned AH_i) - { - CH_assert(AH_i < m_apparent_horizons.size()); - return m_apparent_horizons[AH_i]; - } - - //! returns the index of the AH in m_apparent_horizons - int - add_ah(const AHSurfaceGeometry &a_coord_system, - double a_initial_guess, //!< Initial guess for radius (or whatever - //!< coordinate you're solving for) - const params &a_params, //!< set of AH parameters - const std::string &a_stats = - "stats_AH", //!< name for stats file with - //!< area, spin and AH origin/center - const std::string &a_coords = - "coords_AH", //!< name for coords file with AH - //!< coordinates at each time step) - bool solve_first_step = true //!< whether or not to solve if t=0 - ); - - //! returns the index of the AH in m_apparent_horizons - //! allows for personalized optimizer that finds zero of function - //! 'AHFunction' (that can have some ::params) - int add_ah(const AHSurfaceGeometry &a_coord_system, double a_initial_guess, - const params &a_params, - const typename AHFunction::params &a_func_params, - const std::string &a_stats = "stats_AH", - const std::string &a_coords = "coords_AH", - bool solve_first_step = true); - - // returns the index of the AH in m_apparent_horizons - int add_ah_merger(int ah1, int ah2, const params &a_params); - - //! Find AH; Calculate area and spin; Update center; Print outputs - void solve(double a_dt, double a_time, double a_restart_time); - - bool need_diagnostics(double a_dt, double a_time) - const; //!< is any AH printing diagnostics? Useful if Diagnostics need - //!< to be calculated - - ALWAYS_INLINE void - set_interpolator(AMRInterpolator> *a_interpolator) - { - m_interpolator = a_interpolator; - } - - void - set_origins(const std::vector> &origins, - bool includes_mergers = false); - - private: - //! returns false if 'parent' AHs are too far - //! sets the initial guess and the origin for the merger - bool solve_merger(int ah1, int ah2, double &initial_guess_merger, - std::array &origin_merged); - - ///////////////////////////////////////////////////////// - ///////////////////////////////////////////////////////// - ///////////////// PETSc control methods ///////////////// - ///////////////////////////////////////////////////////// - ///////////////////////////////////////////////////////// - - static bool m_initialized; //!< is initialized? - -#ifdef CH_MPI - static MPI_Group m_mpi_group; //!< set of MPI ranks for PETSc - static MPI_Comm m_mpi_comm; //!< MPI sub-communicator -#endif - - //! define number of ranks of PETSc sub-communicator - static void set_num_ranks(int a_num_ranks); - - //! initialize PETSc and its MPI sub-communicator - static PetscErrorCode PETSc_initialize(int a_num_ranks); - - //! finalize PETSc - static PetscErrorCode PETSc_finalize(); - - public: - //! true if part of PETSc MPI sub-communicator - static bool is_rank_active(); -}; // namespace AHFinder - -#endif /* _AHFINDER_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHFunctionDefault.hpp b/Source/ApparentHorizonFinder/AHFunctionDefault.hpp deleted file mode 100644 index cd7b599f2..000000000 --- a/Source/ApparentHorizonFinder/AHFunctionDefault.hpp +++ /dev/null @@ -1,92 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHFUNCTIONDEFAULT_HPP_ -#define _AHFUNCTIONDEFAULT_HPP_ - -#include "AHData.hpp" -#include "AHDeriv.hpp" -#include "AHGeometryData.hpp" -#include "AlwaysInline.hpp" -#include "Tensor.hpp" - -///////////////////////////////////////////////////////// -// Predefined optimization functions -///////////////////////////////////////////////////////// - -struct AHFunctionDefault -{ - // set the minimum/maximum variable needed (and same for 1st and 2nd - // derivatives) - static ALWAYS_INLINE int vars_min() { return 0; } - static ALWAYS_INLINE int vars_max() { return -1; } - static ALWAYS_INLINE int d1_vars_min() { return 0; } - static ALWAYS_INLINE int d1_vars_max() { return -1; } - static ALWAYS_INLINE int d2_vars_min() { return 0; } - static ALWAYS_INLINE int d2_vars_max() { return -1; } - - // set how many and what variables to print for the geometry when - // 'AH_print_geometry_data' is on - // (useful to print more complex variables as the non-conformal metric or - // the extrinsic curvature, which are not elementary variables in - // BSSN/CCZ4) - static ALWAYS_INLINE int num_write_vars() { return 0; } - static ALWAYS_INLINE void write_headers(std::string *vec) {} - ALWAYS_INLINE void write_vars(double *vec) {} - - // metric to use when calculating the spin and area of the AH - // (override with spacetime metric) - ALWAYS_INLINE const Tensor<2, double> get_metric() const - { - // cartesian flat metric - Tensor<2, double> g = {0.}; - FOR1(i) { g[i][i] = 1.; } - return g; - } - ALWAYS_INLINE const Tensor<2, double> get_extrinsic_curvature() const - { - return {0.}; - } - - // case of reduced Cartoon methods that have an extra metric component -#if GR_SPACEDIM != CH_SPACEDIM // hd - higher dimensions - ALWAYS_INLINE double get_metric_hd() const { return 1.; } -#endif - - struct params // no params needed - { - }; - - // not defined by default - Tensor<1, double> - get_level_function_derivative(const AHGeometryData &geo_data, - const AHDeriv &deriv) const - { - return {0.}; - } - Tensor<1, double> get_spatial_normal_U(const Tensor<1, double> &s_L) const - { - return {0.}; - } - // not defined by default - Tensor<2, double> get_level_function_2nd_covariant_derivative( - const AHGeometryData &geo_data, const AHDeriv &deriv, - const Tensor<1, double> &s_L) const - { - return {0.}; - } - - // WHAT TO ADD TO YOUR OWN FUNCTIONS: - // some constructor with these arguments: - // AHFunctionDefault(const AHData &a_data, - // const Tensor<1, double> &a_coords, - // const Tensor<1, double> &a_coords_cartesian); - - // and some 'get' - // double get(const AHGeometryData &geo_data, const AHDeriv &deriv, - // const params &a_params) const; -}; - -#endif /* _AHFUNCTIONDEFAULT_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHFunctions.hpp b/Source/ApparentHorizonFinder/AHFunctions.hpp deleted file mode 100644 index 30339f5b8..000000000 --- a/Source/ApparentHorizonFinder/AHFunctions.hpp +++ /dev/null @@ -1,374 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHFUNCTIONS_HPP_ -#define _AHFUNCTIONS_HPP_ - -// Chombo includes -#include "CH_Timer.H" - -// Other includes -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists -#include "TensorAlgebra.hpp" -#include "UserVariables.hpp" -#include // sqrt - -#include "AHFunctionDefault.hpp" - -// Chombo namespace -#include "UsingNamespace.H" - -///////////////////////////////////////////////////////// -// Predefined optimization functions -///////////////////////////////////////////////////////// - -struct ExpansionFunction : AHFunctionDefault -{ - ////////////////////////////////// - // needed to calculate expansion: - - // induced metric - Tensor<2, double> g; - Tensor<2, double> g_UU; - Tensor<3, double> dg; - - // extrinsic curvature - Tensor<2, double> K; - double trK; - - // target coordinate (e.g. radius) - double f; - - // case of reduced Cartoon methods that have an extra metric component and - // require the coordinates to calculate the expansion -#if GR_SPACEDIM != CH_SPACEDIM - // hd - higher dimensions - Tensor<1, double> coords; // cartesian coordinates - double g_hd; - Tensor<1, double> dg_hd; - - double Kww; // not needed to calculate expansion, but may be needed for ut - ALWAYS_INLINE double get_metric_hd() const { return g_hd; } -#endif - ////////////////////////////////// - - // only require variables up to Aij (chi, hij, K, Aij) - // this assumes c_Theta comes after the last of Aij - // (and like this the code is generic for 2D and 3D) - static ALWAYS_INLINE int vars_min() { return c_chi; } - static ALWAYS_INLINE int vars_max() { return c_Theta - 1; } - // Derivatives required only for the metric component (chi and hij) - // this assumes c_K comes after the last hij - // (and like this the code is generic for 2D and 3D) - static ALWAYS_INLINE int d1_vars_min() { return c_chi; } - static ALWAYS_INLINE int d1_vars_max() { return c_K - 1; } - - ALWAYS_INLINE const Tensor<2, double> get_metric() const { return g; } - ALWAYS_INLINE const Tensor<2, double> get_extrinsic_curvature() const - { - return K; - } - - static int num_write_vars() - { - int num_components = CH_SPACEDIM * (CH_SPACEDIM + 1) / 2 * 2; -#if GR_SPACEDIM != CH_SPACEDIM - num_components += 2; // hww and Kww -#endif - return num_components; - } - static void write_headers(std::string *vec) - { - int el = 0; - for (int i = 0; i < CH_SPACEDIM; ++i) - for (int j = i; j < CH_SPACEDIM; ++j) - vec[el++] = "g" + std::to_string(i + 1) + std::to_string(j + 1); -#if GR_SPACEDIM != CH_SPACEDIM - vec[el++] = "gww"; -#endif - - for (int i = 0; i < CH_SPACEDIM; ++i) - for (int j = i; j < CH_SPACEDIM; ++j) - vec[el++] = "K" + std::to_string(i + 1) + std::to_string(j + 1); -#if GR_SPACEDIM != CH_SPACEDIM - vec[el++] = "Kww"; -#endif - } - void write_vars(double *vec) - { - int el = 0; - for (int i = 0; i < CH_SPACEDIM; ++i) - for (int j = i; j < CH_SPACEDIM; ++j) - vec[el++] = g[i][j]; -#if GR_SPACEDIM != CH_SPACEDIM - vec[el++] = g_hd; -#endif - for (int i = 0; i < CH_SPACEDIM; ++i) - for (int j = i; j < CH_SPACEDIM; ++j) - vec[el++] = K[i][j]; -#if GR_SPACEDIM != CH_SPACEDIM - vec[el++] = Kww; -#endif - } - - ExpansionFunction(const AHData &a_data, - const Tensor<1, double> &a_coords, - const Tensor<1, double> &a_coords_cartesian) - { - CH_TIME("ExpansionFunction::calculate_data"); - - f = a_coords[CH_SPACEDIM - 1]; - - // * --------------------------- - // * GR-RELATED DATA - // * --------------------------- - - int h[CH_SPACEDIM][CH_SPACEDIM]; - int A[CH_SPACEDIM][CH_SPACEDIM]; - - int comp_h = c_h11; - int comp_A = c_A11; - for (int i = 0; i < CH_SPACEDIM; ++i) - { - for (int j = i; j < CH_SPACEDIM; ++j) - { - h[i][j] = comp_h; - A[i][j] = comp_A; - if (i != j) - { - h[j][i] = comp_h; - A[j][i] = comp_A; - } - ++comp_h; - ++comp_A; - } - } - - const double chi = a_data.vars.at(c_chi); - trK = a_data.vars.at(c_K); - - // INVERSE METRIC - Tensor<2, double, CH_SPACEDIM> h_DD; - FOR2(i, j) { h_DD[i][j] = a_data.vars.at(h[i][j]); } - - Tensor<2, double, CH_SPACEDIM> h_UU = - TensorAlgebra::compute_inverse_sym(h_DD); - - // Reconstructing ADM variables - Tensor<1, double, CH_SPACEDIM> dchi; - - FOR1(i) { dchi[i] = a_data.d1.at(c_chi)[i]; } - - for (int i = 0; i < CH_SPACEDIM; ++i) - { - for (int j = i; j < CH_SPACEDIM; ++j) - { - { - const double gij = h_DD[i][j] / chi; - g[i][j] = gij; - g[j][i] = gij; - } - { - const double g_UUij = h_UU[i][j] * chi; - g_UU[i][j] = g_UUij; - g_UU[j][i] = g_UUij; - } - - { - const double Aij = a_data.vars.at(A[i][j]); - const double Kij = Aij / chi + trK * g[i][j] / GR_SPACEDIM; - K[i][j] = Kij; - K[j][i] = Kij; - } - - for (int k = 0; k < CH_SPACEDIM; ++k) - { - { - const double dhij = a_data.d1.at(h[i][j])[k]; - const double dgijk = - (dhij - (h_DD[i][j] * dchi[k]) / chi) / chi; - dg[i][j][k] = dgijk; - dg[j][i][k] = dgijk; - } - } - } - } - - // part for higher dimensions that use Cartoon Method - // Uli's paper does it for 3+1D (from 5D) - arxiv 1808.05834 -#if GR_SPACEDIM != CH_SPACEDIM - int comp_hww = - c_K - 1; // c_K-1 expected to be c_hww (is c_hww direct better?) - int comp_Aww = c_Aww; // is (c_Theta-1) more general? - Tensor<1, double, CH_SPACEDIM> dhww; - - coords = a_coords_cartesian; - - double hww = a_data.vars.at(comp_hww); - FOR1(a) { dhww[a] = a_data.d1.at(comp_hww)[a]; } - double Aww = a_data.vars.at(comp_Aww); - - g_hd = hww / chi; - FOR1(a) { dg_hd[a] = (dhww[a] - (hww * dchi[a]) / chi) / chi; } - - Kww = Aww / chi + trK * g_hd / GR_SPACEDIM; -#endif - } - - struct params - { - double expansion_radius_power = 1.; - }; - double get(const AHGeometryData &geo_data, const AHDeriv &deriv, - const params &a_params) const - { - Tensor<1, double> s_L = get_level_function_derivative(geo_data, deriv); - Tensor<2, double> Ds = - get_level_function_2nd_covariant_derivative(geo_data, deriv, s_L); - double norm_s; - Tensor<1, double> S_U; - get_spatial_normal_U_and_norm(S_U, norm_s, s_L); - - // calculate D_i S^i and S^i S^j K_ij - double DiSi = 0.; - double Kij_dot_Si_Sj = 0.; - FOR2(a, b) - { - DiSi += (g_UU[a][b] - S_U[a] * S_U[b]) * Ds[a][b] / norm_s; - Kij_dot_Si_Sj += S_U[a] * S_U[b] * K[a][b]; - } - - // Calculation of expansion - as in (6.7.9 / 6.7.13) of Alcubierre - double expansion = DiSi - trK + Kij_dot_Si_Sj; - - // part from extra dimensions in the case of Cartoon methods -#if GR_SPACEDIM != CH_SPACEDIM - expansion += (GR_SPACEDIM - CH_SPACEDIM) * S_U[CH_SPACEDIM - 1] / - coords[CH_SPACEDIM - 1]; - FOR1(a) - { - expansion += - (GR_SPACEDIM - CH_SPACEDIM) * 0.5 * dg_hd[a] / g_hd * S_U[a]; - } -#endif - - // using "r * Expansion" significantly improves the convergence - // (making a Schw. BH converge for any radius >~ 0.5*r_AH instead of - // only up to ~ 3 * r_AH as it happens just with the expansion) - // (see arXiv:gr-qc/0512169, around Fig 3.4) - return expansion * pow(f, a_params.expansion_radius_power); - } - - // extra stuff: - Tensor<1, double> - get_level_function_derivative(const AHGeometryData &geo_data, - const AHDeriv &deriv) const - { - // calculate D_a L of 6.7.12 of Alcubierre for some level function - // L picking L = f - F(u,v) - // D_a L = d_a f - dF/du * du/dx^a - dF/dv * dv/dx^a - - Tensor<1, double> s_L = {0.}; // not normalized, just D_a L - FOR1(a) - { - s_L[a] = geo_data.df[a] - (deriv.duF * geo_data.du[a]) -#if CH_SPACEDIM == 3 - - (deriv.dvF * geo_data.dv[a]) -#endif - ; - } - - return s_L; - } - Tensor<1, double> get_spatial_normal_U(const Tensor<1, double> &s_L) const - { - Tensor<1, double> S_U; - double norm_s; - get_spatial_normal_U_and_norm(S_U, norm_s, s_L); - return S_U; - } - void get_spatial_normal_U_and_norm(Tensor<1, double> &S_U, double &norm_s, - const Tensor<1, double> &s_L) const - { - // calculate S_U = the real 's' of Alcubierre - - // norm of s_L = | D_a L| (the 'u' in 6.7.12 of Alcubierre) - norm_s = 0.0; - FOR2(a, b) { norm_s += g_UU[a][b] * s_L[a] * s_L[b]; } - norm_s = sqrt(norm_s); - - // raise s_L - Tensor<1, double> s_U = {0.}; - FOR2(a, b) { s_U[a] += g_UU[a][b] * s_L[b]; } - - FOR1(a) { S_U[a] = s_U[a] / norm_s; } - } - - Tensor<2, double> get_level_function_2nd_covariant_derivative( - const AHGeometryData &geo_data, const AHDeriv &deriv, - const Tensor<1, double> &s_L) const - { - // calculates D_a D_b L, required for 6.7.13 of Alcubierre - // for the level function L = f - F(u,v) - - Tensor<2, double> ds = {0.}; - FOR2(a, b) - { - ds[a][b] = geo_data.ddf[a][b] - (deriv.duF * geo_data.ddu[a][b]) - - (deriv.duduF * geo_data.du[a] * geo_data.du[b]) -#if CH_SPACEDIM == 3 - - (deriv.dvF * geo_data.ddv[a][b]) - - (deriv.dvdvF * geo_data.dv[a] * geo_data.dv[b]) - - (deriv.dudvF * (geo_data.du[a] * geo_data.dv[b] + - geo_data.du[b] * geo_data.dv[a])) -#endif - ; - } - - // calculate Christoffels on this point (u,v,f) - Tensor<3, double> chris = {0.}; - FOR4(a, b, c, d) - { - chris[a][b][c] += - 0.5 * g_UU[a][d] * (dg[b][d][c] + dg[c][d][b] - dg[b][c][d]); - } - - // covariant derivatrive of s_a to use for DS - Tensor<2, double> Ds = {0.}; - FOR2(a, b) - { - Ds[a][b] = ds[a][b]; - FOR1(c) { Ds[a][b] -= chris[c][a][b] * s_L[c]; } - } - - return Ds; - } -}; - -// to look for chi contours -struct ChiContourFunction : AHFunctionDefault -{ - double chi; - - static int vars_min() { return c_chi; } - static int vars_max() { return c_chi; } - - ChiContourFunction(const AHData &a_data, - const Tensor<1, double> &a_coords, - const Tensor<1, double> &a_coords_cartesian) - { - chi = a_data.vars.at(c_chi); - } - - using params = double; // chi contour to look for - double get(const AHGeometryData &geo_data, const AHDeriv &deriv, - const params &a_params) const - { - double look_for_chi_contour = a_params; - return chi - look_for_chi_contour; - } -}; - -#endif /* _AHFUNCTIONS_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHGeometryData.hpp b/Source/ApparentHorizonFinder/AHGeometryData.hpp deleted file mode 100644 index 67ae91ba7..000000000 --- a/Source/ApparentHorizonFinder/AHGeometryData.hpp +++ /dev/null @@ -1,51 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHGEOMETRYDATA_HPP_ -#define _AHGEOMETRYDATA_HPP_ - -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists -#include "Tensor.hpp" - -#include // memset - -// The d prefix refers to partial derivatives wrt Cartesian coordinates - -struct AHGeometryData -{ - - // jacobian - Tensor<1, double> du; -#if CH_SPACEDIM == 3 - Tensor<1, double> dv; -#endif - Tensor<1, double> df; - - // hessian - Tensor<2, double> ddu; -#if CH_SPACEDIM == 3 - Tensor<2, double> ddv; -#endif - Tensor<2, double> ddf; - - // inverse jacobian - Tensor<1, double> dxdu; -#if CH_SPACEDIM == 3 - Tensor<1, double> dxdv; -#endif - Tensor<1, double> dxdf; - - // inverse hessian - // never needed - - AHGeometryData() - { - // TF: no need to set to 0 - // force all (double) elements of AHGeometryData to be 0 - // memset(this, 0, sizeof(AHGeometryData)); - } -}; - -#endif /* _AHGEOMETRYDATA_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHInterpolation.hpp b/Source/ApparentHorizonFinder/AHInterpolation.hpp deleted file mode 100644 index a6ae0576c..000000000 --- a/Source/ApparentHorizonFinder/AHInterpolation.hpp +++ /dev/null @@ -1,143 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHINTERPOLATION_HPP_ -#define _AHINTERPOLATION_HPP_ - -#include "AHData.hpp" -#include "AHGeometryData.hpp" -#include "AMRInterpolator.hpp" -#include "Lagrange.hpp" - -#include -#include - -//! Class used for interpolation of the variables needed to calculate expansion -//! with the data from a given 'SurfaceGeometry' -template class AHInterpolation -{ - private: - SurfaceGeometry m_coord_system; - AMRInterpolator> *m_interpolator; - - // variables of AH in 'SurfaceGeometry' - // (in Spherical Coordinates Class 'AHSphericalCoords', - // they correspond to 'theta', 'phi' and 'sqrt(radius)') - std::vector m_u; -#if CH_SPACEDIM == 3 - std::vector m_v; -#endif - std::vector m_f; - - // variables of AH in cartesian (for AMR interpolation) - std::vector m_x; - std::vector m_y; -#if CH_SPACEDIM == 3 - std::vector m_z; -#endif - - AHData> m_data; - AHData> m_extra; - - std::array m_coord_min, - m_coord_max; //!< maximum and minimum of level 0 box, used in - //!< 'fit_in_grid' - - //! when PETSc tried to diverge out of the grid, this doesn't let him do so - //! it forces him to stay on the grid, causing non-convergence - bool fit_in_grid(double &x, double &y -#if CH_SPACEDIM == 3 - , - double &z -#endif - ); - - public: - AHInterpolation(const SurfaceGeometry &a_coordSystem, - AMRInterpolator> *a_interpolator); - - const AMRInterpolator> *get_interpolator() const; - const SurfaceGeometry &get_coord_system() const; - - // several of the methods below just call the correspondent method of the - // CoordSystem, as these are needed in the 'ApparentHorizon' class - bool is_u_periodic() const; - double get_domain_u_min() const; //!< lower bound of (u,v) coordinates - double get_domain_u_max() const; //!< upper bound of (u,v) coordinates -#if CH_SPACEDIM == 3 - bool is_v_periodic() const; - double get_domain_v_min() const; //!< lower bound of (u,v) coordinates - double get_domain_v_max() const; //!< upper bound of (u,v) coordinates -#endif - - std::vector get_labels() const; //!< get all names (u, v, f) - - void set_origin( - const std::array &); //!< set origin of CoordSystem - const std::array & - get_origin() const; //!< get origin of CoordSystem - - void refresh_interpolator( - bool printing_step, - const std::map> - &extra_vars); //!< refresh AMRInterpolator 'm_interpolator' - - double get_grid_coord(int a_dir, double f, double u -#if CH_SPACEDIM == 3 - , - double v -#endif - ) const; //!< transform from spherical to cartesian - - //! returns whether any pointis outside of grid (==> diverging) - bool set_coordinates(const std::vector &f, - const std::vector &u, -#if CH_SPACEDIM == 3 - const std::vector &v, -#endif - double add_epsilon = 0.); - const AHGeometryData get_geometry_data(int idx) const; - const Tensor<1, double> get_cartesian_coords(int idx) const; - const Tensor<1, double> get_coords(int idx) const; - const AHData get_data(int idx) const; - - //! verify origin +- initial guess is inside the grid - bool is_in_grid(const std::array &a_origin, - double a_initial_guess); - - //! triplet of functions to be used together in blocks of code that require - //! PETSc AND AMRInterpolator to interpolate - //! 'keep_interpolating_if_inactive' returns 'true' immediately for PETSc - //! ranks for non-PETSc ranks, it loops in a sequence of 'interpolate()' - //! waiting for PETSc ranks to call 'interpolate()'' as well (as this needs - //! to be called by all Chombo processes) can be aborted by calling - //! 'break_interpolation_loop()' for PETSc ranks - /** Example of usage: - * if(m_geom.keep_interpolating_if_inactive()) - * { - * (... PETSc code that calls 'set_coordinates()'...) - * m_geom.break_interpolation_loop(); - * } - */ - bool keep_interpolating_if_inactive(); - void - break_interpolation_loop() const; //!< see 'keep_interpolating_if_inactive' - //!< (this one could be static) - //! returns whether or not to keep interpolating. See - //! 'keep_interpolating_if_inactive' - //! Again, all Chombo_MPI:comm need to run 'interpolate' for it to work (if - //! some are not part of PETSc, they still need to run 'interpolate' on the - //! side) - int interpolate(); - - void interpolate_extra_vars( - const std::map> - &extra_vars); - const AHData get_extra_data(int idx) const; -}; - -#include "AHInterpolation.impl.hpp" - -#endif /* _AHINTERPOLATION_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHInterpolation.impl.hpp b/Source/ApparentHorizonFinder/AHInterpolation.impl.hpp deleted file mode 100644 index fa5b2b558..000000000 --- a/Source/ApparentHorizonFinder/AHInterpolation.impl.hpp +++ /dev/null @@ -1,588 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHINTERPOLATION_HPP_ -#error "This file should only be included through AHInterpolation.hpp" -#endif - -#ifndef _AHINTERPOLATION_IMPL_HPP_ -#define _AHINTERPOLATION_IMPL_HPP_ - -#include "AHFinder.hpp" -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists -#include "TensorAlgebra.hpp" - -template -AHInterpolation::AHInterpolation( - const SurfaceGeometry &a_coord_system, - AMRInterpolator> *a_interpolator) - : m_coord_system(a_coord_system), m_interpolator(a_interpolator) -{ - // below: - // determine maximum and minimum physical coordinate of the grid - // (so that 'fit_in_grid' knows when PETSc has diverged out of the grid and - // doesn't let him do so, - // as it would cause an error in the AMRInterpolator) - - const Box &domainBox = const_cast(m_interpolator->getAMR()) - .getAMRLevels()[0] - ->problemDomain() - .domainBox(); - const IntVect &small_end = domainBox.smallEnd(); - const IntVect &big_end = domainBox.bigEnd(); - - const std::array &coarsest_dx = - m_interpolator->get_coarsest_dx(); - const std::array &coarsest_origin = - m_interpolator->get_coarsest_origin(); - - // set coordinate minimum and maximum as ONE cells before the boundary - for (unsigned i = 0; i < CH_SPACEDIM; ++i) - { - if (m_interpolator->get_boundary_reflective(Side::Lo, i)) - { - m_coord_min[i] = - -(big_end[i] - 1) * coarsest_dx[i] - coarsest_origin[i]; - } - else if (m_interpolator->get_boundary_periodic(i)) - { - m_coord_min[i] = small_end[i] * coarsest_dx[i]; - } - else - { - m_coord_min[i] = - (small_end[i] + 1) * coarsest_dx[i] + coarsest_origin[i]; - } - - if (m_interpolator->get_boundary_reflective(Side::Hi, i)) - { - m_coord_max[i] = - (2 * big_end[i] - 1) * coarsest_dx[i] + coarsest_origin[i]; - } - else if (m_interpolator->get_boundary_periodic(i)) - { - m_coord_max[i] = (big_end[i] + 1) * coarsest_dx[i]; - } - else - { - m_coord_max[i] = - (big_end[i] - 1) * coarsest_dx[i] + coarsest_origin[i]; - } - // pout() << "m_coord = " << m_coord_min[i] << "\t" << m_coord_max[i] - // << " with coarse origin " << coarsest_origin[i] << " and dx " - // << coarsest_dx[i] << std::endl; - } -} - -template -const AMRInterpolator> * -AHInterpolation::get_interpolator() const -{ - return m_interpolator; -} - -template -const SurfaceGeometry & -AHInterpolation::get_coord_system() const -{ - return m_coord_system; -} - -template -bool AHInterpolation::is_u_periodic() const -{ - return m_coord_system.is_u_periodic(); -} -template -double AHInterpolation::get_domain_u_min() const -{ - return m_coord_system.get_domain_u_min(); -} -template -double AHInterpolation::get_domain_u_max() const -{ - return m_coord_system.get_domain_u_max(); -} - -#if CH_SPACEDIM == 3 -template -bool AHInterpolation::is_v_periodic() const -{ - return m_coord_system.is_v_periodic(); -} -template -double AHInterpolation::get_domain_v_min() const -{ - return m_coord_system.get_domain_v_min(); -} -template -double AHInterpolation::get_domain_v_max() const -{ - return m_coord_system.get_domain_v_max(); -} -#endif - -template -std::vector -AHInterpolation::get_labels() const -{ - return - { - m_coord_system.u_name(), -#if CH_SPACEDIM == 3 - m_coord_system.v_name(), -#endif - m_coord_system.param_name() - }; -} - -template -void AHInterpolation::set_origin( - const std::array &origin) -{ - m_coord_system.set_origin(origin); -} - -template -const std::array & -AHInterpolation::get_origin() const -{ - return m_coord_system.get_origin(); -} - -template -bool AHInterpolation::is_in_grid( - const std::array &a_origin, double a_initial_guess) -{ - bool out_of_grid = false; - - double x = a_origin[0]; - double y = a_origin[1]; -#if CH_SPACEDIM == 3 - double z = a_origin[2]; -#endif - -#if CH_SPACEDIM == 3 - out_of_grid |= fit_in_grid(x, y, z); - - x = a_origin[0] - a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); - - x = a_origin[0] + a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); - - x = a_origin[0]; - y = a_origin[1] - a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); - - y = a_origin[1] + a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); - - y = a_origin[1]; - z = a_origin[2] - a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); - - z = a_origin[2] + a_initial_guess; - out_of_grid |= fit_in_grid(x, y, z); -#elif CH_SPACEDIM == 2 - out_of_grid |= fit_in_grid(x, y); - /* - x = a_origin[0] - a_initial_guess; - out_of_grid |= fit_in_grid(x, y); - - x = a_origin[0] + a_initial_guess; - out_of_grid |= fit_in_grid(x, y); - */ - y = a_origin[1] - a_initial_guess; - out_of_grid |= fit_in_grid(x, y); - - y = a_origin[1] + a_initial_guess; - out_of_grid |= fit_in_grid(x, y); -#endif - - return out_of_grid; -} -template -bool AHInterpolation::fit_in_grid(double &x, - double &y -#if CH_SPACEDIM == 3 - , - double &z -#endif -) -{ - CH_TIME("AHInterpolation::fit_in_grid"); - - bool out_of_grid = false; - - // if out of bounds, put back in grid - if (x < m_coord_min[0]) - { - out_of_grid = true; - x = m_coord_min[0]; - } - if (x > m_coord_max[0]) - { - out_of_grid = true; - x = m_coord_max[0]; - } - - if (y < m_coord_min[1]) - { - out_of_grid = true; - y = m_coord_min[1]; - } - if (y > m_coord_max[1]) - { - out_of_grid = true; - y = m_coord_max[1]; - } - -#if CH_SPACEDIM == 3 - if (z < m_coord_min[2]) - { - out_of_grid = true; - z = m_coord_min[2]; - } - if (z > m_coord_max[2]) - { - out_of_grid = true; - z = m_coord_max[2]; - } -#endif - - return out_of_grid; -} - -template -double AHInterpolation::get_grid_coord(int a_dir, - double f, - double u -#if CH_SPACEDIM == 3 - , - double v -#endif -) const -{ - CH_TIME("AHInterpolation::get_grid_coord"); - - return m_coord_system.get_grid_coord(a_dir, f, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); -} - -//! triplet of functions to be used together in blocks of code that require -//! PETSc AND AMRInterpolator to interpolate 'keep_interpolating_if_inactive' -//! returns 'true' immediately for PETSc ranks for non-PETSc ranks, it loops in -//! a sequence of 'interpolate()' waiting for PETSc ranks to call -//! 'interpolate()'' as well (as this needs to be called by all Chombo -//! processes) can be aborted by calling 'break_interpolation_loop()' for PETSc -//! ranks -/** Example of usage: - * if(m_geom.keep_interpolating_if_inactive()) - * { - * (... PETSc code that calls 'set_coordinates()'...) - * m_geom.break_interpolation_loop(); - * } - */ -template -bool AHInterpolation::keep_interpolating_if_inactive() -{ - CH_TIME("AMRInterpolator::keep_interpolating_if_inactive"); - - if (!AHFinder::is_rank_active()) - { - int keep_interpolating = 1; - while (keep_interpolating) - keep_interpolating = interpolate(); - return false; - } - else - return true; -} - -template -void AHInterpolation::break_interpolation_loop() - const -{ - CH_TIME("AMRInterpolator::break_interpolation_loop"); - -#ifdef CH_MPI - // break "keep_interpolating" loop - int ZERO = 0; - int keep_interpolating = 0; - // this is also called in 'interpolate()', breaking the loop - MPI_Allreduce(&ZERO, &keep_interpolating, 1, MPI_INT, MPI_LAND, - Chombo_MPI::comm); -#endif -} -template -int AHInterpolation::interpolate() -{ - CH_TIME("AHInterpolation::interpolate"); - - // Code below used to allow PETSc to run in a sub-communicator. - // For that, PETSc processes run 'set_coordinates' and 'interpolate', while - // non-PETSc processes loop through 'interpolate'. If all cores run - // interpolate, 'MPI_Allreduce' will return 'keep_interpolating' as true. To - // exit the loop for non-PETSc cores, simply do on them a call to: - // MPI_Allreduce(&ZERO, &keep_interpolating, 1, MPI_INT, MPI_LAND, - // Chombo_MPI::comm); (this is done in 'break_interpolation_loop()') - - int keep_interpolating = 1; -#ifdef CH_MPI - int ONE = 1; - MPI_Allreduce(&ONE, &keep_interpolating, 1, MPI_INT, MPI_LAND, - Chombo_MPI::comm); -#endif - - if (!keep_interpolating) - return 0; - - const int n = m_x.size(); - - InterpolationQuery query(n); - query.setCoords(0, m_x.data()).setCoords(1, m_y.data()); -#if CH_SPACEDIM == 3 - query.setCoords(2, m_z.data()); -#endif - - for (int i = AHFunction::vars_min(); i <= AHFunction::vars_max(); ++i) - m_data.set_vars(query, i, i, VariableType::evolution, n); - for (int i = AHFunction::d1_vars_min(); i <= AHFunction::d1_vars_max(); ++i) - m_data.set_d1(query, i, i, VariableType::evolution, n); - for (int i = AHFunction::d2_vars_min(); i <= AHFunction::d2_vars_max(); ++i) - m_data.set_d2(query, i, i, VariableType::evolution, n); - - m_interpolator->interp(query); - - return 1; -} - -template -void AHInterpolation::refresh_interpolator( - bool printing_step, - const std::map> &extra_vars) -{ - CH_TIME("AHInterpolation::refresh_interpolator"); - - // brute force filling of all ghosts: - // m_interpolator->refresh(); - - // OR try to only interpolate the variables needed: - - int min_evolution_var = - std::min(AHFunction::vars_min(), std::min(AHFunction::d1_vars_min(), - AHFunction::d2_vars_min())); - int max_evolution_var = - std::max(AHFunction::vars_max(), std::max(AHFunction::d1_vars_max(), - AHFunction::d2_vars_max())); - - int min_diagnostic_var = -1; - int max_diagnostic_var = -1; - - for (auto &var : extra_vars) - { - int var_enum = std::get<0>(var.second); - VariableType var_type = std::get<1>(var.second); - - if (var_type == VariableType::evolution) - { - if (var_enum < min_evolution_var) - min_evolution_var = var_enum; - else if (var_enum > max_evolution_var) - max_evolution_var = var_enum; - } - else if (printing_step) // diagnostic - { - if (min_diagnostic_var == -1) - { - min_diagnostic_var = var_enum; - max_diagnostic_var = var_enum; - } - else if (var_enum < min_diagnostic_var) - min_diagnostic_var = var_enum; - else if (var_enum > max_diagnostic_var) - max_diagnostic_var = var_enum; - } - } - - // fill ghosts manually to minimise communication - bool fill_ghosts = false; - m_interpolator->refresh(fill_ghosts); - m_interpolator->fill_multilevel_ghosts( - VariableType::evolution, - Interval(min_evolution_var, max_evolution_var)); - if (printing_step && min_diagnostic_var != -1) - { - m_interpolator->fill_multilevel_ghosts( - VariableType::diagnostic, - Interval(min_diagnostic_var, max_diagnostic_var)); - } -} - -// 'set_coordinates' calls 'interpolate'. All Chombo_MPI:comm need to run -// 'interpolate' for it to work -template -bool AHInterpolation::set_coordinates( - const vector &f, const vector &u, -#if CH_SPACEDIM == 3 - const vector &v, -#endif - double add_epsilon) -{ - CH_TIME("AHInterpolation::set_coordinates"); - - CH_assert(f.size() == u.size()); -#if CH_SPACEDIM == 3 - CH_assert(u.size() == v.size()); -#endif - - const int n = f.size(); - - // resize vectors - m_f.resize(n); - m_u.resize(n); -#if CH_SPACEDIM == 3 - m_v.resize(n); -#endif - - m_x.resize(n); - m_y.resize(n); -#if CH_SPACEDIM == 3 - m_z.resize(n); -#endif - - bool out_of_grid = false; - - // Transform to Cartesian - for (int i = 0; i < n; ++i) - { - m_f[i] = f[i] + add_epsilon; - m_u[i] = u[i]; -#if CH_SPACEDIM == 3 - m_v[i] = v[i]; -#endif - -#if CH_SPACEDIM == 3 - m_x[i] = m_coord_system.get_grid_coord(0, m_f[i], m_u[i], m_v[i]); - m_y[i] = m_coord_system.get_grid_coord(1, m_f[i], m_u[i], m_v[i]); - m_z[i] = m_coord_system.get_grid_coord(2, m_f[i], m_u[i], m_v[i]); -#elif CH_SPACEDIM == 2 - m_x[i] = m_coord_system.get_grid_coord(0, m_f[i], m_u[i]); - m_y[i] = m_coord_system.get_grid_coord(1, m_f[i], m_u[i]); -#endif - - // don't let PETSc diverge to outside of the grid (this can happen if - // there is no BH) - out_of_grid |= fit_in_grid(m_x[i], m_y[i] -#if CH_SPACEDIM == 3 - , - m_z[i] -#endif - ); - } - - return out_of_grid; -} - -template -const AHGeometryData -AHInterpolation::get_geometry_data(int idx) const -{ - CH_TIME("AHInterpolation::get_geometry_data"); - - return m_coord_system.get_geometry_data(m_f[idx], m_u[idx] -#if CH_SPACEDIM == 3 - , - m_v[idx] -#endif - ); -} - -template -const Tensor<1, double> -AHInterpolation::get_cartesian_coords( - int idx) const -{ - return - { - m_x[idx], m_y[idx] -#if CH_SPACEDIM == 3 - , - m_z[idx] -#endif - }; -} - -template -const Tensor<1, double> -AHInterpolation::get_coords(int idx) const -{ - return - { - m_u[idx], -#if CH_SPACEDIM == 3 - m_v[idx], -#endif - m_f[idx] - }; -} - -template -const AHData -AHInterpolation::get_data(int idx) const -{ - return get_AHData_idx(idx, m_data); -} - -template -void AHInterpolation::interpolate_extra_vars( - const std::map> &extra_vars) -{ - CH_TIME("AHInterpolation::interpolate_extra_vars"); - - if (extra_vars.size() == 0) - return; - - int n = m_x.size(); - - InterpolationQuery query(n); - query.setCoords(0, m_x.data()).setCoords(1, m_y.data()); -#if CH_SPACEDIM == 3 - query.setCoords(2, m_z.data()); -#endif - - for (auto &var : extra_vars) - { - int var_enum = std::get<0>(var.second); - VariableType var_type = std::get<1>(var.second); - int der_type = std::get<2>(var.second); - - if (der_type == 0) - m_extra.set_vars(query, var_enum, var.first, var_type, n); - else if (der_type == 1) - m_extra.set_d1(query, var_enum, var.first, var_type, n); - else // if(der_type == 2) - m_extra.set_d2(query, var_enum, var.first, var_type, n); - } - - // ghosts for extra evolution or diagnostic vars already filled in - // 'refresh_interpolator' - m_interpolator->interp(query); -} - -template -const AHData -AHInterpolation::get_extra_data(int idx) const -{ - return get_AHData_idx(idx, m_extra); -} - -#endif // _AHINTERPOLATION_IMPL_HPP_ \ No newline at end of file diff --git a/Source/ApparentHorizonFinder/AHSphericalGeometry.hpp b/Source/ApparentHorizonFinder/AHSphericalGeometry.hpp deleted file mode 100644 index a3d0c4ce3..000000000 --- a/Source/ApparentHorizonFinder/AHSphericalGeometry.hpp +++ /dev/null @@ -1,179 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef AHSPHERICALGEOMETRY_HPP_ -#define AHSPHERICALGEOMETRY_HPP_ - -#include "DimensionDefinitions.hpp" // make sure GR_SPACEDIM exists - -#if GR_SPACEDIM != 3 -#error "This file should only be included for GR_SPACEDIM == 3." -#endif - -// Chombo includes -#include "CH_Timer.H" - -// Other includes -#include "AHGeometryData.hpp" -#include "AlwaysInline.hpp" -#include "SphericalGeometry.hpp" -#include -#include - -// Chombo namespace -#include "UsingNamespace.H" - -//! This SurfaceGeometry template class provides spherical shell geometry -//! implementation for the SurfaceExtraction class -//! u = theta, v = phi -class AHSphericalGeometry : public SphericalGeometry -{ - public: - AHSphericalGeometry(const std::array &a_origin) - : SphericalGeometry(a_origin) - { - } - - ALWAYS_INLINE void - set_origin(const std::array &a_origin) - { - m_center = a_origin; - } - ALWAYS_INLINE const std::array &get_origin() const - { - return m_center; - } - - AHGeometryData get_geometry_data(double r, double theta, - double phi = 0.) const - { - // for 2D Cartoon spherical geometry (3D->2D), this should always be - // called with no 3rd argument (so phi=0), which means the formulas - // below are effectively polar coordinates - - CH_TIME("AHSphericalGeometry::get_geometry_data"); - - double costheta = cos(theta); - double sintheta = sin(theta); - double tantheta = tan(theta); - - double cosphi = cos(phi); - double sinphi = sin(phi); - - double cos2theta = cos(2. * theta); - double sin2theta = sin(2. * theta); - double cos2phi = cos(2. * phi); - double sin2phi = sin(2. * phi); - - // in 3D, (x,y,z) are normal - // in 2D Cartoon method, the spherical 'z' axis is the grid 'x' axis, so - // the 2D coordinates (x,y) are (z,x) spherical coordinates - // (that means commenting out everything related to 'y' and to 'v' - // [=phi]) - int iz = m_up_dir; - int ix = (m_up_dir + 1) % 3; -#if CH_SPACEDIM == 3 - int iy = (m_up_dir + 2) % 3; -#endif - - AHGeometryData out; - - out.du[ix] = (costheta * cosphi) / r; -#if CH_SPACEDIM == 3 - out.du[iy] = (costheta * sinphi) / r; -#endif - out.du[iz] = -sintheta / r; - -#if CH_SPACEDIM == 3 - out.dv[ix] = (-sinphi) / (r * sintheta); - out.dv[iy] = (cosphi) / (r * sintheta); - out.dv[iz] = 0; -#endif - - double dfdx = sintheta * cosphi; -#if CH_SPACEDIM == 3 - double dfdy = sintheta * sinphi; -#endif - double dfdz = costheta; - - out.ddu[ix][ix] = - (cos2theta * cosphi * cosphi - cos2phi) / (r * r * tantheta); - out.ddu[ix][iz] = -(cos2theta * cosphi) / (r * r); - out.ddu[iz][ix] = out.ddu[ix][iz]; - out.ddu[iz][iz] = (sin2theta) / (r * r); -#if CH_SPACEDIM == 3 - out.ddu[ix][iy] = - ((cos2theta - 2.) * sin2phi) / (2. * r * r * tantheta); - out.ddu[iy][iy] = - (cos2theta * sinphi * sinphi + cos2phi) / (r * r * tantheta); - out.ddu[iy][iz] = -(cos2theta * sinphi) / (r * r); - out.ddu[iy][ix] = out.ddu[ix][iy]; - out.ddu[iz][iy] = out.ddu[iy][iz]; -#endif - -#if CH_SPACEDIM == 3 - out.ddv[ix][ix] = sin2phi / (r * r * sintheta * sintheta); - out.ddv[ix][iz] = 0.; - out.ddv[iz][ix] = out.ddv[ix][iz]; - out.ddv[iz][iz] = 0.; - out.ddv[ix][iy] = -cos2phi / (r * r * sintheta * sintheta); - out.ddv[iy][iy] = -sin2phi / (r * r * sintheta * sintheta); - out.ddv[iy][iz] = 0.; - out.ddv[iy][ix] = out.ddv[ix][iy]; - out.ddv[iz][iy] = out.ddv[iy][iz]; -#endif - - double ddfdxdx = - (costheta * costheta + sintheta * sintheta * sinphi * sinphi) / r; - double ddfdxdz = -(sintheta * costheta * cosphi) / r; - double ddfdzdz = (sintheta * sintheta) / r; -#if CH_SPACEDIM == 3 - double ddfdxdy = -(sintheta * sintheta * sinphi * cosphi) / r; - double ddfdydy = - (costheta * costheta + sintheta * sintheta * cosphi * cosphi) / r; - double ddfdydz = -(sintheta * costheta * sinphi) / r; -#endif - - out.df[ix] = dfdx; -#if CH_SPACEDIM == 3 - out.df[iy] = dfdy; -#endif - out.df[iz] = dfdz; - - out.ddf[ix][ix] = ddfdxdx; - out.ddf[ix][iz] = ddfdxdz; - out.ddf[iz][ix] = out.ddf[ix][iz]; - out.ddf[iz][iz] = ddfdzdz; -#if CH_SPACEDIM == 3 - out.ddf[ix][iy] = ddfdxdy; - out.ddf[iy][iy] = ddfdydy; - out.ddf[iy][iz] = ddfdydz; - out.ddf[iy][ix] = out.ddf[ix][iy]; - out.ddf[iz][iy] = out.ddf[iy][iz]; -#endif - - out.dxdu[ix] = r * costheta * cosphi; -#if CH_SPACEDIM == 3 - out.dxdu[iy] = r * costheta * sinphi; -#endif - out.dxdu[iz] = -r * sintheta; - -#if CH_SPACEDIM == 3 - out.dxdv[ix] = -r * sintheta * sinphi; - out.dxdv[iy] = r * sintheta * cosphi; - out.dxdv[iz] = 0.; -#endif - - out.dxdf[ix] = sintheta * cosphi; -#if CH_SPACEDIM == 3 - out.dxdf[iy] = sintheta * sinphi; -#endif - out.dxdf[iz] = costheta; - - return out; - } -}; - -#endif /* AHSPHERICALGEOMETRY_HPP_ */ diff --git a/Source/ApparentHorizonFinder/AHStringGeometry.hpp b/Source/ApparentHorizonFinder/AHStringGeometry.hpp deleted file mode 100644 index 4cb410733..000000000 --- a/Source/ApparentHorizonFinder/AHStringGeometry.hpp +++ /dev/null @@ -1,109 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _AHSPHERICALGEOMETRY_HPP_ -#define _AHSPHERICALGEOMETRY_HPP_ - -#if CH_SPACEDIM != 2 -#error "This file should only be included for 2+1D simulations" -#endif - -// Chombo includes -#include "CH_Timer.H" -#include "MayDay.H" - -// Other includes -#include "AHGeometryData.hpp" -#include "AlwaysInline.hpp" -#include - -// Chombo namespace -#include "UsingNamespace.H" - -//! Class for coordinate system to use -//! String assumed to be along in 'x' direction -//! Periodic in x coordinate, y direction may be reflective -class AHStringGeometry -{ - private: - double m_string_length; - std::array fake_origin = {0., 0.}; - - public: - AHStringGeometry(double a_string_length) : m_string_length(a_string_length) - { - } - - ALWAYS_INLINE void - set_origin(const std::array &a_origin) - { - } - ALWAYS_INLINE const std::array &get_origin() const - { - return fake_origin; - } - - ALWAYS_INLINE double get_domain_u_min() const { return 0.; } - ALWAYS_INLINE double get_domain_u_max() const { return m_string_length; } - - ALWAYS_INLINE bool is_u_periodic() const { return true; } - - ALWAYS_INLINE std::string param_name() const { return "y"; } - - ALWAYS_INLINE std::string u_name() const { return "x"; } - - ALWAYS_INLINE double get_grid_coord(int a_dir, double a_y, double a_x) const - { - switch (a_dir) - { - case (0): - return a_x; - case (1): - return a_y; - default: - MayDay::Error("AHStringGeometry: Direction not supported"); - } - } - - AHGeometryData get_geometry_data(double y, double x) const - { - CH_TIME("AHStringGeometry::get_geometry_data"); - - AHGeometryData out; - - out.du[0] = 1.; - out.du[1] = 0.; - - double dfdx = 0.; - double dfdy = 1.; - - out.ddu[0][0] = 0.; - out.ddu[0][1] = 0.; - out.ddu[1][1] = 0.; - out.ddu[1][0] = out.ddu[0][1]; - - double ddfdxdx = 0.; - double ddfdxdy = 0.; - double ddfdydy = 0.; - - out.df[0] = dfdx; - out.df[1] = dfdy; - - out.ddf[0][0] = ddfdxdx; - out.ddf[0][1] = ddfdxdy; - out.ddf[1][1] = ddfdydy; - out.ddf[1][0] = out.ddf[0][1]; - - out.dxdu[0] = 1.; - out.dxdu[1] = 0.; - - out.dxdf[0] = 0.; - out.dxdf[1] = 1.; - - return out; - } -}; - -#endif /* _AHSPHERICALGEOMETRY_HPP_ */ diff --git a/Source/ApparentHorizonFinder/ApparentHorizon.hpp b/Source/ApparentHorizonFinder/ApparentHorizon.hpp deleted file mode 100644 index 233653ef6..000000000 --- a/Source/ApparentHorizonFinder/ApparentHorizon.hpp +++ /dev/null @@ -1,307 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _APPARENTHORIZON_HPP_ -#define _APPARENTHORIZON_HPP_ - -#include -// #include - -#include "AHData.hpp" -#include "AHDeriv.hpp" -#include "AHFinder.hpp" -#include "AHGeometryData.hpp" -#include "AHInterpolation.hpp" -#include "IntegrationMethod.hpp" - -// Class to manage ApparentHorizon for 2+1D and 3+1D simulations -// Class AHFinder manages it -//! AHFunction defines the optimizing function (see AHFunction.hpp for -//! expansion example calculation) -template class ApparentHorizon -{ - using Interpolation = AHInterpolation; - - public: - //! AH that finds the zero of expansion - ApparentHorizon( - const Interpolation &a_interp, //!< Geometry class to exchange data - double a_initial_guess, //!< Initial guess for radius (or whatever - //!< coordinate you're solving for) - const AHFinder::params &a_params, //!< set of AH parameters - const std::string &a_stats = - "stats.dat", //!< name for output file with area, spin and AH origin - const std::string &a_coords = - "coords_", //!< name for output file with AH coordinates at each - //!< time step - bool solve_first_step = true //!< whether or not to solve if t=0 - ); - //! personalized optimizer that finds zero of function - //! 'a_function_to_optimize' (a void* 'a_function_to_optimize_params' can be - //! passed for auxiliary parameters passed to 'a_function_to_optimize') - ApparentHorizon(const Interpolation &a_interp, double a_initial_guess, - const AHFinder::params &a_params, - const typename AHFunction::params &a_func_params, - const std::string &a_stats = "stats", - const std::string &a_coords = "coords", - bool solve_first_step = true); - ~ApparentHorizon(); - - void solve(double a_dt, double a_time, double a_restart_time); - - bool good_to_go(double a_dt, double a_time) const; - bool get_converged() - const; //!< PETSc didn't converge last time solve() was called - int get_failed_convergences() - const; //!< PETSc didn't converge last time solve() was called - bool has_been_found() const; //!< PETSc converged once and then - //!< stopped converging (AH collapsed) - - const std::array &get_origin() const; - const std::array &get_center() const; - double get_initial_guess() const; - void set_initial_guess(double a_initial_guess); - const Interpolation &get_ah_interp() const; - - // set origin to whatever you want (e.g. punctures) before solving if you - // want, otherwise we use the last center or the estimate next center - void set_origin(const std::array &a_origin); - - double get_max_F() const; - double get_min_F() const; - double get_ave_F() const; - double get_std_F() const; - - bool do_solve(double a_dt, double a_time) - const; //!< decide (based times passed to 'solve') whether or not - //!< to print (uses AHFinder::params::solve_interval) - bool do_print(double a_dt, double a_time) - const; //!< decide when to print (only AHFinder::params::print_interval - //!< out of all 'solve's) - - const AHFinder::params &m_params; //!< set of AH parameters - const std::string m_stats, - m_coords; //!< public base names for output files (no need for a set as - //!< they are const) - - // any parameters that want to be saved to be passed to optimization - // function - typename AHFunction::params m_func_params; - - private: - void write_outputs(double a_dt, double a_time, double a_restart_time); - - void - restart(bool solve_first_step); //!< restart AH, updating the coordinates - //!< based on the last output file and the - //!< origin based on the stats file - - double calculate_area(); //!< calculate AH area - -#if CH_SPACEDIM == 3 - Tensor<1, double> calculate_angular_momentum_J(); //!< calculate spin, ONLY - //!< FOR 3D - double calculate_spin_dimensionless( - double a_area); //!< calculate spin with 'z' direction, ONLY FOR 3D -#endif - // estimate based on area and angular momentum J - ALWAYS_INLINE double calculate_mass(double area, double J_norm) - { - return sqrt(area / (16. * M_PI) + J_norm * J_norm * 4. * M_PI / area); - } - ALWAYS_INLINE double calculate_irreducible_mass(double area) - { - return calculate_mass(area, 0.); - } - - void calculate_minmax_F() const; - void calculate_average_F() const; - - void reset_initial_guess(); - - void update_old_centers(std::array); - void predict_next_origin(); - std::array - calculate_center(); //!< update location of center by calculating the - //!< centroid of the AH - - void write_coords_file(double a_dt, double a_time, double a_restart_time, - const std::string &filename, - bool write_geometry_data = false) - const; //!< write coords in (m_u, m_v, m_F) to 'filename' - - //! write metric and extrinsic curvature for each point of AH - void write_geometry_data(const std::string &a_file_prefix, double a_dt, - double a_time, double a_restart_time) const; - - void check_convergence(); //!< check if PETSc has converged and share that - //!< info across all Chombo processes - - // default to simpson rule and change if invalid - void check_integration_methods(); - - void initialise_PETSc(); //!< initialise automatically done in constructor - void finalise_PETSc(); //!< finalise automatically done in destructor - - // variables - private: - bool m_printed_once; - - int m_converged; //!< flag saying if PETSc has converged the last 'N' times - //!<(read using 'get_converged()') - bool m_has_been_found; //!< flag saying if an horizon has ever been found - //!< (read using 'has_been_found()') - int m_num_failed_convergences; //!< the number of failed consecutive - //!< convergences - - double m_initial_guess; //!< initial guess for AH (saved so that it can be - //!< re-used when atempting to solve again) - - //! used to estimate where the new center will be before solving - //! (currently last 2 centers saved to make a parabolic regression) - //! (elements ordered in reverse time - older with higher index) - std::vector> m_old_centers; - - // mutable so that 'const' objects can still change them - //! stores the global maximum and minimum of F - calculate with - //! calculate_minmax_F - mutable double m_max_F, m_min_F, m_ave_F, m_std_F; - - // Integration method to calculate area and spin of AH. Defaults to Simpson - // method - std::array m_integration_methods; - - // just to save the result temporarily at each iteration - double m_area, m_spin, m_mass, m_irreducible_mass, m_spin_z_alt; - Tensor<1, double> m_dimensionless_spin_vector; - - // prevents resetting the origin when the user externally did 'set_origin' - // before 'solve' - bool origin_already_updated; - - ///////////////////////////////////////////////////////////////////////// - /////////////////////////// PETSc stuff below /////////////////////////// - ///////////////////////////////////////////////////////////////////////// - -#if CH_SPACEDIM == 3 - typedef PetscScalar **dmda_arr_t; -#elif CH_SPACEDIM == 2 - typedef PetscScalar *dmda_arr_t; -#endif - - //! Geometries of the AH - //! 'm_geom_plus' and 'm_geom_minus' are used to calculate the - //! jacobian of the expansion using a 'delta' numerical differentiation - Interpolation m_interp; - Interpolation m_interp_plus; - Interpolation m_interp_minus; - - //!< used to compute jacobian of expansion (numerical differentiation) - static constexpr double eps = 1e-7; - - const bool m_periodic_u; //!< is 'u' periodic? - PetscInt m_num_global_u; //!< total number of grid points in 'u' coordinate - double m_du; //!< physical 'delta' in 'u' coordinate - -#if CH_SPACEDIM == 3 - const bool m_periodic_v; //!< is 'v' periodic? - PetscInt m_num_global_v; //!< total number of grid points in 'u' coordinate - double m_dv; //!< physical 'delta' in 'v' coordinate -#endif - - //! vectors to store and manipulate 'F', u', and 'v' - //! internally, in interaction with the 'Interpolation's - std::vector m_F; - std::vector m_u; -#if CH_SPACEDIM == 3 - std::vector m_v; -#endif - - //! minimums and maximums of coordinates 'u' and 'v' - //! of the PETSc grid specific to the current rank - PetscInt m_umin; - PetscInt m_umax; - -#if CH_SPACEDIM == 3 - PetscInt m_vmin; - PetscInt m_vmax; -#endif - - //! number of points in 'u' and 'v' direction - //! (m_nu = m_umax - m_umin) - PetscInt m_nu; -#if CH_SPACEDIM == 3 - PetscInt m_nv; -#endif - - // PETSc main object - DM m_dmda; - //! Scalable Nonlinear Equations Solvers - SNES m_snes; - - Vec m_snes_soln; - Vec m_snes_rhs; - Mat m_snes_jac; - - //! interpolate (u,v) 2D grid points at restart if number of points in - //! either direction changed - //! This only interpolates the points that the PETSc rank that called it has - //! returns whether or not points were interpolated - bool interpolate_ah(dmda_arr_t f, - const std::vector> &old_coords); - - //! set the default stencils of AHDeriv at position {u,v} - void set_stencils(AHDeriv &out, int u -#if CH_SPACEDIM == 3 - , - int v -#endif - ); - - //! function to calculate 1st and 2nd derivatives of 'in' - //! (tipically corresponds to our 'f' function) - //! in the 'u' and 'v' directions - AHDeriv diff(const dmda_arr_t in, int u -#if CH_SPACEDIM == 3 - , - int v -#endif - ); - - //! private functions used to compute the RHS (the expansion) and it's - //! jacobian - void form_function(Vec F, Vec Rhs); - void form_jacobian(Vec F, Mat J); - - //! helper for 'form_jacobian' - double point_jacobian(int u, int u_stencil, -#if CH_SPACEDIM == 3 - int v, int v_stencil, -#endif - dmda_arr_t in, int idx, - const Interpolation &interp_plus, - const Interpolation &interp_minus); - - //! functions used by PETSc based on 'form_function' and 'form_jacobian' - static PetscErrorCode Petsc_form_function(SNES snes, Vec F, Vec Rhs, - void *ptr); - - static PetscErrorCode -#if PETSC_VERSION_GE(3, 5, 0) - Petsc_form_jacobian(SNES snes, Vec F, Mat Amat, Mat Pmat, void *ptr); -#else - Petsc_form_jacobian(SNES snes, Vec F, Mat *Amat, Mat *Pmat, - MatStructure *flag, void *ptr); -#endif - - // monitor function required for SNES - static PetscErrorCode Petsc_SNES_monitor(SNES snes, PetscInt its, - PetscReal norm, void *ptr); -}; - -#include "ApparentHorizon.impl.hpp" -#include "ApparentHorizon_petsc.impl.hpp" - -#endif /* _APPARENTHORIZON_HPP_ */ diff --git a/Source/ApparentHorizonFinder/ApparentHorizon.impl.hpp b/Source/ApparentHorizonFinder/ApparentHorizon.impl.hpp deleted file mode 100644 index 4830ee4d6..000000000 --- a/Source/ApparentHorizonFinder/ApparentHorizon.impl.hpp +++ /dev/null @@ -1,1951 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _APPARENTHORIZON_HPP_ -#error "This file should only be included through ApparentHorizon.hpp" -#endif - -#ifndef _APPARENTHORIZON_IMPL_HPP_ -#define _APPARENTHORIZON_IMPL_HPP_ - -#include "GRAMR.hpp" -#include "SmallDataIO.hpp" -#include "TensorAlgebra.hpp" -#include "UserVariables.hpp" - -// for 1D interpolation of restart centers -#include "Lagrange.hpp" -#include "SimpleArrayBox.hpp" -#include "SimpleInterpSource.hpp" - -#include // for 'setw' in stringstream -#include // infinity for min and max calculation - -template -ApparentHorizon::ApparentHorizon( - const AHInterpolation &a_interp, - double a_initial_guess, const AHFinder::params &a_params, - const std::string &a_stats, const std::string &a_coords, - bool solve_first_step) - : ApparentHorizon(a_interp, a_initial_guess, a_params, AHFunction::params(), - a_stats, a_coords, solve_first_step) -{ -} - -template -ApparentHorizon::ApparentHorizon( - const AHInterpolation &a_interp, - double a_initial_guess, const AHFinder::params &a_params, - const typename AHFunction::params &a_func_params, - const std::string &a_stats, const std::string &a_coords, - bool solve_first_step) - : m_params(a_params), - - m_stats(a_stats), m_coords(a_coords), - - m_func_params(a_func_params), - - m_printed_once(false), - - m_converged(0), - - m_has_been_found(false), - - m_num_failed_convergences(0), - - m_initial_guess(a_initial_guess), - - m_old_centers(3, a_interp.get_coord_system().get_origin()), - - m_max_F(0.), m_min_F(0.), m_ave_F(0.), m_std_F(0.), - - m_integration_methods({ -#if CH_SPACEDIM == 3 - IntegrationMethod::simpson, -#endif - IntegrationMethod::simpson - }), - - m_area(NAN), m_spin(NAN), m_mass(NAN), m_irreducible_mass(NAN), - m_spin_z_alt(NAN), m_dimensionless_spin_vector({NAN}), - - origin_already_updated(false), - - m_interp(a_interp), m_interp_plus(a_interp), m_interp_minus(a_interp), - - m_periodic_u(a_interp.is_u_periodic()), - m_num_global_u(a_params.num_points_u) - -#if CH_SPACEDIM == 3 - , - m_periodic_v(a_interp.is_v_periodic()), - m_num_global_v(a_params.num_points_v) -#endif -{ - initialise_PETSc(); - set_origin(a_interp.get_coord_system().get_origin()); - check_integration_methods(); - restart(solve_first_step); -} -template -ApparentHorizon::~ApparentHorizon() -{ - finalise_PETSc(); -} - -template -bool ApparentHorizon::good_to_go( - double a_dt, double a_time) const -{ - if (a_time < m_params.start_time - 1.e-7 /*just to be safe*/) - return false; - if (m_params.give_up_time >= 0. && a_time >= m_params.give_up_time && - !get_converged()) - return false; - - bool is_lost = - (m_params.max_fails_after_lost >= 0 - ? m_num_failed_convergences > m_params.max_fails_after_lost - : false); - - // stop if it has been found but was lost - return do_solve(a_dt, a_time) && !(has_been_found() && is_lost); -} - -template -bool ApparentHorizon::do_solve(double a_dt, - double a_time) const -{ - CH_assert(a_dt != 0); // Check if time was set! - return !(((int)(std::round(a_time / a_dt))) % m_params.solve_interval); -} -template -bool ApparentHorizon::do_print(double a_dt, - double a_time) const -{ - CH_assert(a_dt != 0); // Check if time was set! - return (get_converged() || !has_been_found()) && - !(((int)(std::round(a_time / a_dt))) % - (m_params.solve_interval * m_params.print_interval)); -} - -template -const std::array & -ApparentHorizon::get_origin() const -{ - return m_interp.get_origin(); -} - -template -const AHInterpolation & -ApparentHorizon::get_ah_interp() const -{ - return m_interp; -} - -template -void ApparentHorizon::set_origin( - const std::array &a_origin) -{ - m_interp.set_origin(a_origin); - m_interp_plus.set_origin(a_origin); - m_interp_minus.set_origin(a_origin); - origin_already_updated = true; - - if (m_params.verbose > AHFinder::SOME) - { - pout() << "Setting origin to (" << a_origin[0] << "," << a_origin[1] -#if CH_SPACEDIM == 3 - << "," << a_origin[2] -#endif - << ")" << std::endl; - } -} - -template -const std::array & -ApparentHorizon::get_center() const -{ - return m_old_centers[0]; -} -template -bool ApparentHorizon::get_converged() const -{ - return m_converged; -} - -template -int ApparentHorizon::get_failed_convergences() - const -{ - return m_num_failed_convergences; -} -template -bool ApparentHorizon::has_been_found() const -{ - return m_has_been_found; -} -template -double ApparentHorizon::get_max_F() const -{ - if (m_max_F <= 0.) - calculate_minmax_F(); - return m_max_F; -} -template -double ApparentHorizon::get_min_F() const -{ - if (m_min_F <= 0.) - calculate_minmax_F(); - return m_min_F; -} -template -double ApparentHorizon::get_ave_F() const -{ - if (m_ave_F <= 0.) - calculate_average_F(); - return m_ave_F; -} -template -double ApparentHorizon::get_std_F() const -{ - if (m_std_F <= 0.) - calculate_average_F(); - return m_std_F; -} - -template -double ApparentHorizon::get_initial_guess() const -{ - return m_initial_guess; -} - -template -void ApparentHorizon::set_initial_guess( - double a_initial_guess) -{ - m_initial_guess = a_initial_guess; -} -template -void ApparentHorizon::reset_initial_guess() -{ - CH_TIME("ApparentHorizon::reset_initial_guess"); - - if (!AHFinder::is_rank_active()) - return; - - auto origin = m_interp.get_origin(); - - // verify origin +- initial guess is inside the grid - bool out_of_grid = m_interp.is_in_grid(origin, m_initial_guess); - CH_assert(!out_of_grid); - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Setting Initial Guess to f=" << m_initial_guess - << " centered at (" << origin[0] << "," << origin[1] -#if CH_SPACEDIM == 3 - << "," << origin[2] -#endif - << ")" << std::endl; - } - - // read PETSc array to 'f' - dmda_arr_t f; - DMDAVecGetArray(m_dmda, m_snes_soln, &f); - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { -#if CH_SPACEDIM == 3 - double &f_point = f[v][u]; -#else - double &f_point = f[u]; -#endif - - f_point = m_initial_guess; - } - } - - // write PETSc array back - DMDAVecRestoreArray(m_dmda, m_snes_soln, &f); -} -template -void ApparentHorizon::predict_next_origin() -{ - std::array new_center = m_old_centers[0]; - if (m_converged >= 3) // add 2nd derivative - { - FOR1(a) - { - new_center[a] += (m_old_centers[0][a] + m_old_centers[2][a] - - 2. * m_old_centers[1][a]); - } - if (m_params.verbose > AHFinder::SOME) - { - pout() << "OLD[-2]: (" << m_old_centers[2][0] << "," - << m_old_centers[2][1] -#if CH_SPACEDIM == 3 - << "," << m_old_centers[2][2] -#endif - << ")" << std::endl; - } - } - if (m_converged >= 2) // add 1st derivative - { - FOR1(a) - { - new_center[a] += (m_old_centers[0][a] - m_old_centers[1][a]); - } - - if (m_params.verbose > AHFinder::SOME) - { - pout() << "OLD[-1]: (" << m_old_centers[1][0] << "," - << m_old_centers[1][1] -#if CH_SPACEDIM == 3 - << "," << m_old_centers[1][2] -#endif - << ")" << std::endl; - pout() << "OLD[0]: (" << m_old_centers[0][0] << "," - << m_old_centers[0][1] -#if CH_SPACEDIM == 3 - << "," << m_old_centers[0][2] -#endif - << ")" << std::endl; - } - if (m_params.verbose > AHFinder::MIN) - { - pout() << "Estimated center: (" << new_center[0] << "," - << new_center[1] -#if CH_SPACEDIM == 3 - << "," << new_center[2] -#endif - << ")" << std::endl; - } - } - - set_origin(new_center); -} - -template -void ApparentHorizon::update_old_centers( - std::array new_center) -{ - FOR1(a) - { - if (!m_interp.get_interpolator()->get_boundary_reflective(Side::Lo, - a) && - !m_interp.get_interpolator()->get_boundary_reflective(Side::Hi, a)) - { - m_old_centers[2][a] = m_old_centers[1][a]; - m_old_centers[1][a] = m_old_centers[0][a]; - m_old_centers[0][a] = new_center[a]; - } - } -} - -template -void ApparentHorizon::solve(double a_dt, - double a_time, - double a_restart_time) -{ - CH_TIME("ApparentHorizon::solve"); - if (!good_to_go(a_dt, a_time)) - return; - - m_interp.refresh_interpolator( - do_print(a_dt, a_time), - m_params.extra_vars); // (ALL CHOMBO ranks do it!!) - - // estimate the next position where origin will be - if (!origin_already_updated && get_converged()) - { - if (m_params.predict_origin) - predict_next_origin(); - else if (m_params.track_center) - set_origin(get_center()); - } - - // PETSc processes go inside 'if', others "wait" until 'if' gets to - // 'm_interp.break_interpolation_loop()' - if (m_interp.keep_interpolating_if_inactive()) - { - CH_TIME("ApparentHorizon::solve::solving"); - - if (!get_converged()) - reset_initial_guess(); // reset initial guess if diverged (or in - // first attempt) - - // actual solve happens here! - SNESSolve(m_snes, NULL, m_snes_soln); - - PetscInt its; - SNESGetIterationNumber(m_snes, &its); - if (m_params.verbose > AHFinder::MIN) - { - pout() << "SNES Iteration number " << its << endl; - } - SNESGetLinearSolveIterations(m_snes, &its); - if (m_params.verbose > AHFinder::MIN) - { - pout() << "KSP Iteration number " << its << endl; - } - m_interp.break_interpolation_loop(); - } - - { - if (m_params.verbose > AHFinder::SOME) - { - pout() << "In [ApparentHorizon::solve::post-solving]" << endl; - } - - CH_TIME("ApparentHorizon::solve::post-solving"); - - bool save_converged = get_converged(); - - // ask PETSc if it converged - check_convergence(); - - // retry with initial guess once if failed (and everything was ok - // before) - // on the 2nd iteration, 'save_converged' will be false and it will skip - // the 'if' - if (save_converged && !get_converged() && m_params.allow_re_attempt) - { - --m_num_failed_convergences; // reduce failed and try again - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Re-attempting to solve using initial guess." - << std::endl; - } - solve(a_dt, a_time, a_restart_time); - return; // do nothing else - } - - // update center - // note that after this, the point that corresponds to the (u,v,r) - // coordinates is 'get_origin()', not 'get_center()' - if (m_params.track_center) - calculate_center(); - - origin_already_updated = false; - - m_area = calculate_area(); -#if GR_SPACEDIM == 3 // GR_SPACEDIM, not CH_SPACEDIM !!! -#if CH_SPACEDIM == 3 - Tensor<1, double> J = calculate_angular_momentum_J(); - double J_norm = sqrt(J[0] * J[0] + J[1] * J[1] + J[2] * J[2]); -#elif CH_SPACEDIM == 2 - double J_norm = 0.; -#endif - m_mass = calculate_mass(m_area, J_norm); - m_irreducible_mass = calculate_irreducible_mass(m_area); - -#if CH_SPACEDIM == 3 - m_spin = J_norm / m_mass; - FOR1(a) { m_dimensionless_spin_vector[a] = J[a] / (m_mass * m_mass); } - - m_spin_z_alt = calculate_spin_dimensionless(m_area); -#endif - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "mass = " << m_mass << endl; -#if CH_SPACEDIM == 3 - pout() << "spin = " << m_spin << endl; -#endif - if (m_params.verbose > AHFinder::MIN) - { - pout() << "irreducible mass = " << m_irreducible_mass << endl; -#if CH_SPACEDIM == 3 - pout() << "dimensionless spin vector = (" - << m_dimensionless_spin_vector[0] << ", " - << m_dimensionless_spin_vector[1] << ", " - << m_dimensionless_spin_vector[2] << ")" << endl; - pout() << "dimensionless spin in z (from equator-length " - "integral) = " - << m_spin_z_alt << endl; -#endif - } - } -#endif - - // reset min and max F, to force re-calculation - m_max_F = 0.; - m_min_F = 0.; - m_ave_F = 0.; - m_std_F = 0.; - } - - write_outputs(a_dt, a_time, a_restart_time); - - // break if 'AH_stop_if_max_fails == true' - // only break after re-attempting - // (if max fails reached or give up time reached) - if (m_params.stop_if_max_fails && - ((m_params.max_fails_after_lost >= 0 && - m_num_failed_convergences > m_params.max_fails_after_lost) || - (m_params.give_up_time >= 0. && a_time >= m_params.give_up_time && - !get_converged()))) - MayDay::Error("Reached max fails. Stopping. Parameter " - "'stop_if_max_fails' is set to true."); - - if (m_params.verbose > AHFinder::SOME) - { - pout() << "ApparentHorizon::solve finished successfully!" << endl; - } -} - -template -void ApparentHorizon::write_outputs( - double a_dt, double a_time, double a_restart_time) -{ - // print step. Printing the step allows the user to modify params as - // 'solve_interval', 'print_interval' or even 'dt_multiplier' (that change - // the frequency of writing outputs) and the class will still be able to - // restart from the correct file (this only applies for frequency increase, - // for which the file number will suddenly jump, as for frequency decrease - // the AHFinder may override old coord files) - - // print stats (area, spin, origin, center) and coordinates - // stop printing if it stopped converging - // then in 'write_coords_file' nothing is printed if not converged - if (do_print(a_dt, a_time)) - { - CH_TIME("ApparentHorizon::solve::printing"); - if (m_params.verbose > AHFinder::MIN) - { - pout() << "Printing statistics and coordinates." << std::endl; - } - - // write stats - double fake_dt = - a_dt * m_params.solve_interval * m_params.print_interval; - SmallDataIO file(m_stats, fake_dt, a_time, a_restart_time, - SmallDataIO::APPEND, !m_printed_once); - - file.remove_duplicate_time_data(); - - // std::string coords_filename = file.get_new_file_number(fake_dt, - // a_time); - - CH_assert(CH_SPACEDIM == 3 || CH_SPACEDIM == 2); - // first '1' corresponds to 'step' - // area+mass+irred.mass+spin+spin_vec+spin_alt OR area+mass+irred.mass - // in 2D Cartoon OR only area for other cases (e.g. 4D -> 2D cartoon) - int stats = (GR_SPACEDIM == 3 ? (CH_SPACEDIM == 3 ? 8 : 3) : 1); - std::vector values(1 + stats + - CH_SPACEDIM * (1 + m_params.track_center)); - - auto origin = get_origin(); - - int step = std::round(a_time / fake_dt); - - int idx = 0; - values[idx++] = step; - values[idx++] = m_area; -#if GR_SPACEDIM == 3 // GR_SPACEDIM, not CH_SPACEDIM !!! - values[idx++] = m_mass; - values[idx++] = m_irreducible_mass; -#if CH_SPACEDIM == 3 - values[idx++] = m_spin; - FOR1(a) { values[idx++] = m_dimensionless_spin_vector[a]; } - values[idx++] = m_spin_z_alt; -#endif -#endif - for (int i = 0; i < CH_SPACEDIM; ++i) - values[idx++] = origin[i]; - if (m_params.track_center) - for (int i = 0; i < CH_SPACEDIM; ++i) - values[idx++] = m_old_centers[0][i]; - - // print headers to stats file in the beginning of evolution - if (!m_printed_once) - { - std::vector headers( - 1 + stats + CH_SPACEDIM * (1 + m_params.track_center)); - - idx = 0; - headers[idx++] = "file"; - headers[idx++] = "area"; -#if GR_SPACEDIM == 3 // GR_SPACEDIM, not CH_SPACEDIM !!! - headers[idx++] = "mass"; - headers[idx++] = "irreducible mass"; -#if CH_SPACEDIM == 3 - headers[idx++] = "spin"; - headers[idx++] = "dimless spin-x"; - headers[idx++] = "dimless spin-y"; - headers[idx++] = "dimless spin-z"; - headers[idx++] = "dimless spin-z-alt"; -#endif -#endif - headers[idx++] = "origin_x"; - headers[idx++] = "origin_y"; -#if CH_SPACEDIM == 3 - headers[idx++] = "origin_z"; -#endif - if (m_params.track_center) - { - headers[idx++] = "center_x"; - headers[idx++] = "center_y"; -#if CH_SPACEDIM == 3 - headers[idx++] = "center_z"; -#endif - } - - file.write_header_line(headers); - - m_printed_once = true; - } - - file.write_time_data_line(values); - - // write coordinates - m_interp.interpolate_extra_vars(m_params.extra_vars); - write_coords_file(a_dt, a_time, a_restart_time, m_coords, - m_params.print_geometry_data); - } -} -template -void ApparentHorizon::check_convergence() -{ - CH_TIME("ApparentHorizon::check_convergence"); - - // Check SNES convergence reasons in: - // https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESConvergedReason.html#SNESConvergedReason - - int result; - if (AHFinder::is_rank_active()) - { - SNESConvergedReason reason; - SNESGetConvergedReason(m_snes, &reason); - - // result will be 0 if any of the PETSc ranks says 'reason <=0' (PETSc - // convergence error) - int tmp = reason; // convert just to make sure it's 'int' -#ifdef CH_MPI - MPI_Allreduce(&tmp, &result, 1, MPI_INT, MPI_MIN, Chombo_MPI::comm); -#else - result = tmp; -#endif - } - else - { - int ALOT = 100; // smaller than all SNESConvergedReasons -#ifdef CH_MPI - MPI_Allreduce(&ALOT, &result, 1, MPI_INT, MPI_MIN, Chombo_MPI::comm); -#else - result = ALOT; -#endif - } - - bool converged = (result > 0); - - if ((bool)converged) - { - ++m_converged; - m_num_failed_convergences = 0; - } - else - { - m_converged = 0; // reset - ++m_num_failed_convergences; - } - - if (m_params.verbose > AHFinder::NONE) - { - pout() << (m_converged ? "Solver converged. Horizon FOUND." - : "Solver diverged. Horizon NOT found.") - << std::endl; - - if (m_params.verbose > AHFinder::MIN) - { - pout() << "SNESConvergedReason = " << result << std::endl; - } - } - - if (!m_has_been_found && m_converged) - m_has_been_found = true; // finally found :D -} - -template -void ApparentHorizon::restart( - bool solve_first_step) -{ - CH_TIME("ApparentHorizon::restart"); - - static double eps = 1.e-7; - - // GUIDING PRINCIPLE : restart should not depend on m_params, as these might - // have changed! We should be able to figure everything out without them - - const GRAMR &gramr = - (static_cast(m_interp.get_interpolator()->getAMR())); - - int current_step = AMR::s_step; - int restart_step = gramr.get_restart_step(); - - bool solve_time_0_now = (!m_params.re_solve_at_restart && - current_step == 0 && restart_step < 0); - - // solve if at t=0 and not a restart - if (solve_time_0_now) - { - // 'dt' doesn't matter for 1st time step, so set to 1 - if (solve_first_step) - solve(1., 0., 0.); - return; - } - - double current_time = gramr.getCurrentTime(); - double coarse_dt = (current_step == 0 ? 0. : current_time / current_step); - double level_dt = coarse_dt * pow(2., -m_params.level_to_run); - double current_solve_dt = level_dt * m_params.solve_interval; - - // READ STATS - - // get centre from stats file - std::string file = m_stats + ".dat"; - auto stats = SmallDataIO::read(file); - - int idx = 0; - double old_print_dt = 0.; - if (stats.size() == 0) - { // case when it never ran the AH or the file doesn't exist - if (m_params.verbose > AHFinder::NONE && current_step != 0) - { - pout() << "Empty stats file '" << file - << "'. Assuming AHFinder was not run yet for this AH." - << std::endl; - } - m_printed_once = false; - m_has_been_found = false; - m_converged = 0; - m_num_failed_convergences = 0; - } - else - { - // look for stats line with time 'current_time', as there may be - // further output messed up in the file after the last checkpoint was - // made - idx = stats[0].size(); - while ((--idx) >= 0 && stats[0][idx] - current_time > eps) - ; - - // figure out what was the old 'dt' used - // this may have changed if 'level_to_run' changed of if 'dt_multiplier' - // changed, or 'solve_interval', or 'print_interval' - if (idx > 0) - old_print_dt = stats[0][idx] - stats[0][idx - 1]; - else if (idx == 0 && stats[0].size() > 1) - old_print_dt = stats[0][idx + 1] - - stats[0][idx]; // there's still time steps forward - else // if we can't know, just default to the current one - old_print_dt = current_solve_dt * m_params.print_interval; - - if (m_params.verbose > AHFinder::SOME) - { - if (idx >= 0) - pout() << "stats[0][idx] = " << stats[0][idx] << std::endl; - pout() << "stats[0].size() = " << stats[0].size() << std::endl; - pout() << "current_time = " << current_time << std::endl; - pout() << "old_print_dt = " << old_print_dt << std::endl; - } - - m_printed_once = true; // don't delete old - - if (idx < 0) - { // case when job was restarted at a point before the AH was first ever - // found - if (m_params.verbose > AHFinder::NONE) - { - pout() - << "First time step is after restart in '" << file - << "'. Assuming AHFinder started after current step, but " - "before " - "PETSc had never converged." - << std::endl; - } - m_has_been_found = false; - m_converged = 0; - m_num_failed_convergences = 0; - } - else if (stats[0][idx] - current_time < - -(old_print_dt == 0. ? eps : old_print_dt - eps)) - { // case when the PETSc stopped converging and stopped printing - // so there is a mismatch with the last time in the file - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Last time step not found in '" << file - << "'. Assuming AH stopped converging." << std::endl; - pout() << "(number of failed convergences will be set to 1)" - << std::endl; - } - m_has_been_found = true; - m_converged = 0; - m_num_failed_convergences = 1; - } - else if (std::isnan(stats[2][idx])) - { // case when the simulation was still going without ever - // having converged - // OR when job was restarted at a point before the AH was first ever - // found - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Last time step is NAN in '" << file - << "'. Assuming AH wasn't found and PETSc didn't " - "converged." - << std::endl; - pout() << "(number of failed convergences will be set to 1)" - << std::endl; - } - m_has_been_found = false; - m_converged = 0; - m_num_failed_convergences = 1; - } - else - { // case when the AH was found and there is a coordinate file - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Last time step found in '" << file - << "'. Reading coordinates file." << std::endl; - } - m_has_been_found = true; - m_converged = 1; - m_num_failed_convergences = 0; - } - } - - // force restart if interpolation was needed (because #points changed) - // or if time step found is not the current time - // 'int' in order to use MPI_INT later - int force_restart = false; - - // if not converged, leave - // origin as the initial guessed origin - // and 'solve' will reset coords to initial guess - if (m_converged) - { - //////////////////// - // READ ORIGIN - //////////////////// - - int cols = stats.size(); - - std::array origin; -#if GR_SPACEDIM == 3 // GR_SPACEDIM, not CH_SPACEDIM !!! -#if CH_SPACEDIM == 3 - bool was_center_tracked = - (cols > - 10 + - CH_SPACEDIM); // 10 for time + file + area + mass + irreducible - // + spin + spin_vector + spin_z_alt -#else - // 3D -> 2D Cartoon method - bool was_center_tracked = - (cols > 5 + CH_SPACEDIM); // 5 for time + file + area + mass + - // irreducible mass -#endif -#else - bool was_center_tracked = - (cols > 3 + CH_SPACEDIM); // 3 for time + file + area -#endif - - int offset = CH_SPACEDIM + was_center_tracked * CH_SPACEDIM; - FOR1(a) { origin[a] = stats[cols - offset + a][idx]; } - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Setting origin from stats file '" << m_stats << "' at (" - << origin[0] << "," << origin[1] -#if CH_SPACEDIM == 3 - << "," << origin[2] -#endif - << ")" << std::endl; - } - - set_origin(origin); - - //////////////////// - // READ CENTERS - //////////////////// - - // this is only relevant if track_center was on - if (was_center_tracked) - { - // SET OLD CENTERS directly if old 'dt' matches current 'dt' - if (std::abs(old_print_dt - current_solve_dt) < eps) - { - if (m_params.verbose > AHFinder::NONE && - m_params.predict_origin) - { - pout() << "Reading old centers to predict next origin." - << std::endl; - } - - // skip i=0, this one is directly in the file and is done after - for (int i = m_old_centers.size() - 1; i > 0; --i) - { - if (idx >= i && !std::isnan(stats[2][idx - 1])) - { - m_converged++; - update_old_centers({ -#if CH_SPACEDIM == 3 - stats[cols - 3][idx - i], -#endif - stats[cols - 2][idx - i], - stats[cols - 1][idx - i] - }); - } - } - } - else - { - // calculate last NAN to make sure we don't calculate points of - // when AH didn't converge - int last_nan_idx = stats[0].size(); - while ((--last_nan_idx) >= 0 && - !std::isnan(stats[2][last_nan_idx])) - ; - - // interpolate old centers when dt has changed - std::vector old_centers_time_index; - // skip i=0, this one is directly in the file and is done after - for (int i = 1; i < m_old_centers.size(); ++i) - { - double time = current_time - current_solve_dt * i; - double index = idx - (current_time - time) / old_print_dt; - - if (index >= last_nan_idx) - old_centers_time_index.push_back(index); - } - - if (m_params.verbose > AHFinder::NONE && - m_params.predict_origin) - { - if (old_centers_time_index.size() == 0) - pout() << "AH time step changed and no old centers " - "able to " - "use for prediction." - << std::endl; - else - pout() - << "Old AH time step is different than current one " - "(this might happen if 'AH_print_interval != " - "1').\nRecovering " - << old_centers_time_index.size() - << " old centers from interpolation, for accurate " - "prediction of next origin." - << std::endl; - } - - int rows = stats[0].size(); - SimpleInterpSource<1> source({rows}); - SimpleArrayBox<1> box_x({rows}, stats[cols - CH_SPACEDIM]); - SimpleArrayBox<1> box_y({rows}, stats[cols - CH_SPACEDIM + 1]); -#if CH_SPACEDIM == 3 - SimpleArrayBox<1> box_z({rows}, stats[cols - CH_SPACEDIM + 2]); -#endif - - const int Order = 2; - Lagrange interpolator(source); - for (int i = old_centers_time_index.size() - 1; i >= 0; --i) - { - std::array old_center; - - interpolator.setup({0}, {old_print_dt}, - {old_centers_time_index[i]}); - old_center[0] = interpolator.interpData(box_x); - old_center[1] = interpolator.interpData(box_y); -#if CH_SPACEDIM == 3 - old_center[2] = interpolator.interpData(box_z); -#endif - - m_converged++; - update_old_centers(old_center); - - if (m_params.verbose > AHFinder::SOME) - { - pout() << "OLD[-" << i + 1 << "] = (" << old_center[0] - << "," << old_center[1] -#if CH_SPACEDIM == 3 - << "," << old_center[2] -#endif - << ")" << std::endl; - } - } - } - } - - // save current center no matter what - // place back OLD[0] - THE CENTER - update_old_centers({ -#if CH_SPACEDIM == 3 - stats[cols - 3][idx], -#endif - stats[cols - 2][idx], stats[cols - 1][idx] - }); - - //////////////////// - // READ COORDINATES - //////////////////// - - // determine last step in which there was an AH output - - // the lines commented below also work, but not if there's only one - // stats, as then 'old_print_dt' can't be defined std::string - // coords_filename = SmallDataIO::get_new_filename(m_coords, - // old_print_dt, current_time); - - int coords_file_number = stats[1][idx]; - std::string coords_filename = SmallDataIO::get_new_filename( - m_coords, 1. /*fake dt*/, coords_file_number); - - auto coords = SmallDataIO::read(coords_filename); - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Setting Initial Guess to previous file '" - << coords_filename << "' found." << std::endl; - } - - // now write local points based on file (or interpolated file) - if (AHFinder::is_rank_active()) - { - // read PETSc array to 'f' - dmda_arr_t f; - DMDAVecGetArray(m_dmda, m_snes_soln, &f); - - // doesn't change anything if number of points remained the same - force_restart |= interpolate_ah(f, coords); - - // write PETSc array back - DMDAVecRestoreArray(m_dmda, m_snes_soln, &f); - } - -#ifdef CH_MPI - int tmp = force_restart; - MPI_Allreduce(&tmp, &force_restart, 1, MPI_INT, MPI_LOR, - Chombo_MPI::comm); -#endif - } - - // if t=0, solve only if it's a restart and if solve_first_step = true (it - // may be false for example for a merger of a binary, for which we don't - // want to solve!) - bool solve_no_matter_what = - (m_params.re_solve_at_restart && - (current_step != 0 || - (current_step == 0 && solve_first_step && restart_step >= 0))); - - if (force_restart || solve_no_matter_what) - solve(current_time == 0. ? 1. : level_dt, current_time, current_time); -} - -template -void ApparentHorizon::write_coords_file( - double a_dt, double a_time, double a_restart_time, - const std::string &filename, bool write_geometry_data) const -{ - CH_TIME("ApparentHorizon::write_coords_file"); - - if (!get_converged()) - { - // old way: writing an empty file - // new way: just don't write any coords file if no AH was founds - // coords_file.write_header_line({"Horizon NOT found."}, ""); - return; - } - - if (m_params.verbose > AHFinder::NONE && write_geometry_data) - pout() << "Writing geometry data." << std::endl; - - CH_assert(a_dt != 0); // Check if time was set!! - - double fake_dt = a_dt * m_params.solve_interval * m_params.print_interval; - SmallDataIO file(filename, fake_dt, a_time, a_restart_time, - SmallDataIO::NEW, true); - - unsigned num_components; // only related to 'write_geometry_data' - if (write_geometry_data) - num_components = AHFunction::num_write_vars(); - else - num_components = 0; - - // Write headers - std::vector components(num_components + - m_params.num_extra_vars); - int el = 0; - - // extra vars headers - for (auto &var : m_params.extra_vars) - { - int var_enum = std::get<0>(var.second); - VariableType var_type = std::get<1>(var.second); - int der_type = std::get<2>(var.second); - - std::string var_name; - if (var_type == VariableType::evolution) - var_name = UserVariables::variable_names[var_enum]; - else - var_name = DiagnosticVariables::variable_names[var_enum]; - - static std::string xyz[CH_SPACEDIM] = { - "x", - "y" -#if CH_SPACEDIM == 3 - , - "z" -#endif - }; - - if (der_type == 0) - components[el++] = var_name; - else if (der_type == 1) - for (int i = 0; i < CH_SPACEDIM; ++i) - components[el++] = "d" + xyz[i] + "_" + var_name; - else // if (der_type == 2) - for (int i = 0; i < CH_SPACEDIM; ++i) - for (int j = i; j < CH_SPACEDIM; ++j) - components[el++] = - "d" + xyz[i] + "d" + xyz[j] + "_" + var_name; - } - - // geometry headers - if (write_geometry_data) - AHFunction::write_headers(&components[el]); - file.write_header_line(components, m_interp.get_labels()); - - ////////////////// - // method: - // 1) all PESc ranks send their coordinates to rank 0 - // (non-PETSc processes will do nothing) - // 2) rank 0 receives and stores all coordinates - // 3) non-blocking sends finish the communication ('MPI_Wait') - // 4) rank 0 writes all the information - // Note: this does NOT assume that rank 0 is part of the PETSc communicator! - // (even though for now the PETSc processes always include rank 0) - ////////////////// - -#ifdef CH_MPI - // not needed for non-PETSc ranks, but they will be waiting anyway - MPI_Request *requests = nullptr; // needed for non-blocking communication -#endif - double *output = nullptr; // needed for non-blocking communication - unsigned num_components_total = - num_components + m_params.num_extra_vars + CH_SPACEDIM; - - // step 1 - if (AHFinder::is_rank_active()) - { - -#ifdef CH_MPI - int rank_petsc; - MPI_Comm_rank(PETSC_COMM_WORLD, &rank_petsc); - - // make sure rank 0 of Chombo is rank 0 of PETSc - // this must be the case just because it is rank 0 of PETSc - // that will be writing in the file, and SmallDataIO assumes it's rank 0 - // of Chombo (!) the one writing - // Alternative 1: if rank 0 of PETSc is different than rank 0 of Chombo, - // transfer all the data to Chombo's rank 0 before writing - // Alternative 2: find a way to tell SmallDataIO to use another rank for - // writing - int rank = 0; - MPI_Comm_rank(Chombo_MPI::comm, &rank); - if ((rank != 0 && rank_petsc == 0) || (rank == 0 && rank_petsc != 0)) - MayDay::Error("PETSc's rank 0 should be Chombo's rank 0"); -#endif - -#if CH_SPACEDIM == 3 - int local_total = (m_vmax - m_vmin) * (m_umax - m_umin); -#elif CH_SPACEDIM == 2 - int local_total = (m_umax - m_umin); -#endif - -#ifdef CH_MPI - requests = new MPI_Request[local_total]; -#endif - output = new double[local_total * num_components_total]; - - int idx = 0; -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { -#if CH_SPACEDIM == 3 - output[idx * num_components_total] = m_u[idx]; - output[idx * num_components_total + 1] = m_v[idx]; - output[idx * num_components_total + 2] = m_F[idx]; -#elif CH_SPACEDIM == 2 - output[idx * num_components_total] = m_u[idx]; - output[idx * num_components_total + 1] = m_F[idx]; -#endif - - auto extra = m_interp.get_extra_data(idx); - - int el = 0; - - // extra vars - for (auto &var : m_params.extra_vars) - { - // int var_enum = std::get<0>(var.second); - // VariableType var_type = std::get<1>(var.second); - int der_type = std::get<2>(var.second); - - if (der_type == 0) - { - output[idx * num_components_total + CH_SPACEDIM + - el++] = extra.vars.at(var.first); - } - else if (der_type == 1) - { - for (int i = 0; i < CH_SPACEDIM; ++i) - output[idx * num_components_total + CH_SPACEDIM + - el++] = extra.d1.at(var.first)[i]; - } - else // if (der_type == 2) - { - for (int i = 0; i < CH_SPACEDIM * (CH_SPACEDIM + 1) / 2; - ++i) - output[idx * num_components_total + CH_SPACEDIM + - el++] = extra.d2.at(var.first)[i]; - } - } - - // geometry vars - if (write_geometry_data) - { - const auto data = m_interp.get_data(idx); - const auto coords = m_interp.get_coords(idx); - const auto coords_cart = m_interp.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - func.write_vars( - &output[idx * num_components_total + CH_SPACEDIM + el]); - } - -#ifdef CH_MPI - // all processes send their 'output' to rank 0, who receives - // and writes everything (to simplify, rank 0 also sends it - // to itself, otherwise we wouldn't know which tags rank 0 - // has) sends are tagged by global index, so that receives - // can be unique and writes indexed correctly -#if CH_SPACEDIM == 3 - int idx_global = v * m_num_global_u + u; -#elif CH_SPACEDIM == 2 - int idx_global = u; -#endif - MPI_Issend(&output[idx * num_components_total], - num_components_total, MPI_DOUBLE, 0, idx_global, - PETSC_COMM_WORLD, &(requests[idx])); -#else - // only rank 0 is active -> serial - file.write_data_line( - std::vector(output + idx * num_components_total + - CH_SPACEDIM, - output + idx * num_components_total + - num_components_total), - std::vector(output + idx * num_components_total, - output + idx * num_components_total + - CH_SPACEDIM)); -#endif - - idx++; - } - } - - // step 2 -#ifdef CH_MPI -#if CH_SPACEDIM == 3 - int total = m_num_global_u * m_num_global_v; -#elif CH_SPACEDIM == 2 - int total = m_num_global_u; -#endif - double temp[total * num_components_total]; - if (rank_petsc == 0) - { - for (unsigned i = 0; i < total; ++i) - { - MPI_Status status; - MPI_Recv(&temp[i * num_components_total], num_components_total, - MPI_DOUBLE, MPI_ANY_SOURCE, - i /* use i as tag to make them ordered*/, - PETSC_COMM_WORLD, &status); - } - } - - // step 3 - idx = 0; -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { - MPI_Wait(&(requests[idx]), MPI_STATUS_IGNORE); - idx++; - } - } - - delete[] requests; // free allocated memory - delete[] output; - - // step 4 - if (rank_petsc == 0) - { - for (unsigned i = 0; i < total; ++i) - { - file.write_data_line( - std::vector( - &temp[i * num_components_total] + CH_SPACEDIM, - &temp[i * num_components_total] + num_components_total), - std::vector(&temp[i * num_components_total], - &temp[i * num_components_total] + - CH_SPACEDIM)); - } - } -#endif - } -} - -template -void ApparentHorizon::check_integration_methods() -{ - // check if integration methods are valid given periodicity and number of - // points - bool valid_u = m_integration_methods[0].is_valid(m_params.num_points_u, - m_interp.is_u_periodic()); - - IntegrationMethod method_default = IntegrationMethod::trapezium; - if (!valid_u) - { - std::string warn = - "ApparentHorizon: Simpson IntegrationMethod for u is " - "not valid with this num_points_u.\n" - "Reverting to trapezium rule."; - MayDay::Warning(warn.c_str()); - pout() << warn << std::endl; - m_integration_methods[0] = method_default; - } - -#if CH_SPACEDIM == 3 - bool valid_v = m_integration_methods[1].is_valid(m_params.num_points_v, - m_interp.is_v_periodic()); - if (!valid_v) - { - std::string warn = - "ApparentHorizon: Simpson IntegrationMethod for v is " - "not valid with this num_points_v.\n" - "Reverting to trapezium rule."; - MayDay::Warning(warn.c_str()); - pout() << warn << std::endl; - m_integration_methods[1] = method_default; - } -#endif -} - -#if CH_SPACEDIM == 3 -// ONLY FOR 3D -// OLD method to calculate the spin using the 1D equator length -// Adopted area integral 'calculate_angular_momentum_J' that allows to get the -// direction of the spin as well -template -double -ApparentHorizon::calculate_spin_dimensionless( - double a_area) -{ - CH_assert(CH_SPACEDIM == 3); - CH_TIME("ApparentHorizon::calculate_spin_dimensionless"); - - if (!get_converged()) - return NAN; - - double equator_length; - double integral = 0.; // temporary, but defined outside to be in scope for - // non-PETSc processes - - // PETSc processes go inside 'if', others "wait" until 'if' gets to - // 'm_interp.break_interpolation_loop()' - if (m_interp.keep_interpolating_if_inactive()) - { - - int idx = 0; - - Vec localF; - - DMGetLocalVector(m_dmda, &localF); - DMGlobalToLocalBegin(m_dmda, m_snes_soln, INSERT_VALUES, localF); - DMGlobalToLocalEnd(m_dmda, m_snes_soln, INSERT_VALUES, localF); - - dmda_arr_t in; - DMDAVecGetArray(m_dmda, localF, &in); - - // m_F is already set from solve - - int u_equator = std::round(M_PI / 2.0 / m_du); - - for (int v = m_vmin; v < m_vmax; ++v) - { - for (int u = m_umin; u < m_umax; ++u) - { - if (u_equator == u) - { - AHDeriv deriv = diff(in, u, v); - const auto geometry_data = m_interp.get_geometry_data(idx); - const auto data = m_interp.get_data(idx); - const auto coords = m_interp.get_coords(idx); - const auto coords_cart = m_interp.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - auto &g = func.get_metric(); - - double dxdv[3]; - dxdv[0] = geometry_data.dxdv[0] + - deriv.dvF * geometry_data.dxdf[0]; - dxdv[1] = geometry_data.dxdv[1] + - deriv.dvF * geometry_data.dxdf[1]; - dxdv[2] = geometry_data.dxdv[2] + - deriv.dvF * geometry_data.dxdf[2]; - - double norm2 = 0.; - FOR2(i, j) norm2 += g[i][j] * dxdv[i] * dxdv[j]; - CH_assert(norm2 >= 0.); - - double weight = m_integration_methods[1].weight( - v, m_params.num_points_v, m_interp.is_v_periodic()); - - integral += sqrt(norm2) * weight * m_dv; - } - - idx++; - } - } - - DMDAVecRestoreArray(m_dmda, localF, &in); - DMRestoreLocalVector(m_dmda, &localF); - - m_interp.break_interpolation_loop(); - } - - // reduction across all Chombo processes (note that 'integral' is 0 for - // non-PETSc processes) because SmallDataIO uses rank 0 to write this - // ensures rank 0 will have 'integral' (even though for now the PETSc - // processes always include rank 0) -#ifdef CH_MPI - MPI_Allreduce(&integral, &equator_length, 1, MPI_DOUBLE, MPI_SUM, - Chombo_MPI::comm); -#else // serial - equator_length = integral; -#endif - - // mass could be estimated based on "equator_length / (4. * M_PI)" - // but calculation using area as well is probably more precise - - double factor = - ((2. * M_PI * a_area / (equator_length * equator_length)) - 1.); - double spin_dimensionless = - (factor > 1. ? 0 - : sqrt(1. - factor * factor)); // factor>1 means numerical - // error with spin as 0 - - return spin_dimensionless; -} -#endif - -#if CH_SPACEDIM == 3 -// ONLY FOR 3D -template -Tensor<1, double> -ApparentHorizon::calculate_angular_momentum_J() -{ - CH_assert(CH_SPACEDIM == 3); - CH_TIME("ApparentHorizon::calculate_angular_momentum_J"); - - if (!get_converged()) - return {NAN}; - - Tensor<1, double> J; - Tensor<1, double> integrals = {0.}; // temporary, but defined outside to be - // in scope for non-PETSc processes - - // PETSc processes go inside 'if', others "wait" until 'if' gets to - // 'm_interp.break_interpolation_loop()' - if (m_interp.keep_interpolating_if_inactive()) - { - - int idx = 0; - - Vec localF; - - DMGetLocalVector(m_dmda, &localF); - DMGlobalToLocalBegin(m_dmda, m_snes_soln, INSERT_VALUES, localF); - DMGlobalToLocalEnd(m_dmda, m_snes_soln, INSERT_VALUES, localF); - - dmda_arr_t in; - DMDAVecGetArray(m_dmda, localF, &in); - - // m_F is already set from solve - - const std::array ¢er = get_center(); - - for (int v = m_vmin; v < m_vmax; ++v) - { - Tensor<1, double> inner_integral = {0.}; - for (int u = m_umin; u < m_umax; ++u) - { - AHDeriv deriv = diff(in, u, v); - const auto geometric_data = m_interp.get_geometry_data(idx); - const auto data = m_interp.get_data(idx); - const auto coords = m_interp.get_coords(idx); - const auto coords_cart = m_interp.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - Tensor<2, double> g = func.get_metric(); - Tensor<2, double> K = func.get_extrinsic_curvature(); - Tensor<1, double> s_L = - func.get_level_function_derivative(geometric_data, deriv); - Tensor<1, double> S_U = func.get_spatial_normal_U(s_L); - - Tensor<1, double> coords_cart_centered; - FOR1(i) - { - coords_cart_centered[i] = coords_cart[i] - center[i]; - } - - Tensor<1, double> spin_integrand = {0.}; - double directions[3][3] = { - {0, -coords_cart_centered[2], coords_cart_centered[1]}, - {coords_cart_centered[2], 0., -coords_cart_centered[0]}, - {-coords_cart_centered[1], coords_cart_centered[0], 0.}}; - FOR3(a, b, c) - { - // as in arXiv:gr-qc/0206008 eq. 25 - // but computing with 'killing vector' in directions 'x,y,z' - // ('x,y,z' killing vector are directions[c][a] = - // epsilon[c][j][a] * x[j], just like for ADM Momentum) - spin_integrand[c] += directions[c][a] * S_U[b] * K[a][b]; - } - - // now calculate sqrt(-g) area element - - // Calculate Jacobian matrix for transformation from Cartesian - // to (f,u,v) coords - Tensor<2, double> Jac; - FOR1(k) - { - Jac[0][k] = geometric_data.dxdf[k]; - Jac[1][k] = geometric_data.dxdu[k]; - Jac[2][k] = geometric_data.dxdv[k]; - } - - // Now do the coordinate transformation - Tensor<2, double> g_spherical = {0.}; - FOR4(i, j, k, l) - { - g_spherical[i][j] += Jac[i][k] * Jac[j][l] * g[k][l]; - } - - // Construct the 2-metric on the horizon in (u,v) coords - // i.e. substitute df = (df/du)du + (df/dv)dv - // into the spherical metric - Tensor<2, double, CH_SPACEDIM - 1> g_horizon = {0.}; - g_horizon[0][0] = g_spherical[1][1] + - g_spherical[0][0] * deriv.duF * deriv.duF + - 2.0 * g_spherical[0][1] * deriv.duF; - g_horizon[1][1] = g_spherical[2][2] + - g_spherical[0][0] * deriv.dvF * deriv.dvF + - 2.0 * g_spherical[0][2] * deriv.dvF; - g_horizon[0][1] = g_horizon[1][0] = - g_spherical[1][2] + - g_spherical[0][0] * deriv.duF * deriv.dvF + - g_spherical[0][1] * deriv.dvF + - g_spherical[0][2] * deriv.duF; - - double det = TensorAlgebra::compute_determinant(g_horizon); - - double weight = m_integration_methods[0].weight( - u, m_params.num_points_u, m_interp.is_u_periodic()); - - FOR1(a) - { - double element = spin_integrand[a] / (8. * M_PI) * - sqrt(det) * weight * m_du; - - // hack for the poles (theta=0,\pi) - // where nans appear in 'spin_integrand', but - // 'det' should be 0 - if (std::isnan(element)) - element = 0.; - - inner_integral[a] += element; - } - idx++; - } - double weight = m_integration_methods[1].weight( - v, m_params.num_points_v, m_interp.is_v_periodic()); - FOR1(a) { integrals[a] += weight * m_dv * inner_integral[a]; } - } - - DMDAVecRestoreArray(m_dmda, localF, &in); - DMRestoreLocalVector(m_dmda, &localF); - - m_interp.break_interpolation_loop(); - } - - // reduction across all Chombo processes (note that 'integral' is 0 for - // non-PETSc processes) because SmallDataIO uses rank 0 to write this - // ensures rank 0 will have 'integral' (even though for now the PETSc - // processes always include rank 0) -#ifdef CH_MPI - MPI_Allreduce(&integrals, &J, GR_SPACEDIM, MPI_DOUBLE, MPI_SUM, - Chombo_MPI::comm); -#else // serial - FOR1(a) { J[a] = integrals[a]; } -#endif - - return J; -} -#endif - -template -double ApparentHorizon::calculate_area() -{ - CH_TIME("ApparentHorizon::calculate_area"); - - if (!get_converged()) - return NAN; - - double area; - double integral = - 0.; // temporary, but defined outside to exist for non-PETSc processes - - // PETSc processes go inside 'if', others "wait" until 'if' gets to - // 'm_interp.break_interpolation_loop()' - if (m_interp.keep_interpolating_if_inactive()) - { - - int idx = 0; - - Vec localF; - - DMGetLocalVector(m_dmda, &localF); - DMGlobalToLocalBegin(m_dmda, m_snes_soln, INSERT_VALUES, localF); - DMGlobalToLocalEnd(m_dmda, m_snes_soln, INSERT_VALUES, localF); - - dmda_arr_t in; - DMDAVecGetArray(m_dmda, localF, &in); - - // m_F is already set from solve - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - double inner_integral = 0.; - for (int u = m_umin; u < m_umax; ++u) - { - const auto geometric_data = m_interp.get_geometry_data(idx); - const auto data = m_interp.get_data(idx); - const auto coords = m_interp.get_coords(idx); - const auto coords_cart = m_interp.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - Tensor<2, double> g = func.get_metric(); - - // Calculate Jacobian matrix for transformation from Cartesian - // to (f,u,v) coords - Tensor<2, double> Jac; - FOR1(k) - { - Jac[0][k] = geometric_data.dxdf[k]; - Jac[1][k] = geometric_data.dxdu[k]; -#if CH_SPACEDIM == 3 - Jac[2][k] = geometric_data.dxdv[k]; -#endif - } - - AHDeriv deriv = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - // Now do the coordinate transformation - Tensor<2, double> g_spherical = {0.}; - FOR4(i, j, k, l) - { - g_spherical[i][j] += Jac[i][k] * Jac[j][l] * g[k][l]; - } - - // Construct the 2-metric on the horizon in (u,v) coords - // i.e. substitute df = (df/du)du + (df/dv)dv - // into the spherical metric - Tensor<2, double, CH_SPACEDIM - 1> g_horizon = {0.}; - g_horizon[0][0] = g_spherical[1][1] + - g_spherical[0][0] * deriv.duF * deriv.duF + - 2.0 * g_spherical[0][1] * deriv.duF; -#if CH_SPACEDIM == 3 - g_horizon[1][1] = g_spherical[2][2] + - g_spherical[0][0] * deriv.dvF * deriv.dvF + - 2.0 * g_spherical[0][2] * deriv.dvF; - g_horizon[0][1] = g_horizon[1][0] = - g_spherical[1][2] + - g_spherical[0][0] * deriv.duF * deriv.dvF + - g_spherical[0][1] * deriv.dvF + - g_spherical[0][2] * deriv.duF; -#endif - - double det = TensorAlgebra::compute_determinant(g_horizon); - - double weight = m_integration_methods[0].weight( - u, m_params.num_points_u, m_interp.is_u_periodic()); - - double element = sqrt(det) * weight * m_du; - -// assume a (GR_SPACEDIM - CH_SPACEDIM)-sphere leftover -#if GR_SPACEDIM != CH_SPACEDIM - double n_sphere = (GR_SPACEDIM - CH_SPACEDIM); - element *= pow(sqrt(func.get_metric_hd()) * - coords_cart[CH_SPACEDIM - 1], - n_sphere); -#endif - - inner_integral += element; - idx++; - } -#if CH_SPACEDIM == 3 - double weight = m_integration_methods[1].weight( - v, m_params.num_points_v, m_interp.is_v_periodic()); - integral += weight * m_dv * inner_integral; -#elif CH_SPACEDIM == 2 - integral += inner_integral; -#endif - } - -// do this separate so that the gamma is only calculated once -#if GR_SPACEDIM != CH_SPACEDIM - double n_sphere = (GR_SPACEDIM - CH_SPACEDIM); - // this is 2pi for n=1, 4pi for n=2, 2pi^2 for n=3, ... - double n_sphere_coeff = 2. * std::pow(M_PI, (n_sphere + 1.) / 2.) / - std::tgamma((n_sphere + 1.) / 2.); - integral *= n_sphere_coeff; -#endif - - DMDAVecRestoreArray(m_dmda, localF, &in); - DMRestoreLocalVector(m_dmda, &localF); - - m_interp.break_interpolation_loop(); - } - -// reduction across all Chombo processes (note that 'area' is 0 for non-PETSc -// processes) because SmallDataIO uses rank 0 to write this ensures rank 0 will -// have the area (even though for now the PETSc processes always include rank 0) -#ifdef CH_MPI - // we want all the processes to have the area so that all use it in the - // spin calculation (which in reality is needed not for the output files, - // but only for the pout() prints) - MPI_Allreduce(&integral, &area, 1, MPI_DOUBLE, MPI_SUM, Chombo_MPI::comm); -#else // serial - area = integral; -#endif - - if (m_params.verbose > AHFinder::MIN) - { - pout() << "area = " << area << endl; - } - - return area; -} - -template -std::array -ApparentHorizon::calculate_center() -{ - CH_TIME("ApparentHorizon::calculate_center"); - - if (!get_converged()) - { -#if CH_SPACEDIM == 3 - return {NAN, NAN, NAN}; -#elif CH_SPACEDIM == 2 - return {NAN, NAN}; -#endif - } - - // [OLD] Method: - // Calculate centroid of all the points, by summing the coordinates {x,y,z} - // of all the points, which are spread across processors, and then reducing - // them all with MPI. Finally, divide by total number of points summed, to - // get the centroid - // Above is not good. Replace by calculating the maximum and minimum point - // of the surface and estimating the center by their average - - // std::array temp = {0., 0., 0.}; // old method - std::array max_temp = get_origin(); - std::array min_temp = max_temp; - - int idx = 0; - if (AHFinder::is_rank_active()) // in principle this wouldn't be needed, as - // non-PETSc processes wouldn't even enter - // the loop anyways - { -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { - for (unsigned i = 0; i < CH_SPACEDIM; ++i) - { - double coord_i = - m_interp.get_grid_coord(i, m_F[idx], m_u[idx] -#if CH_SPACEDIM == 3 - , - m_v[idx] -#endif - ); - - // temp[i] += point[i]; // old method - if (coord_i > max_temp[i]) - max_temp[i] = coord_i; - if (coord_i < min_temp[i]) - min_temp[i] = coord_i; - } - - idx++; - } - } - } - - std::array max_point; - std::array min_point; - std::array center; - // int idx_sum; - -#ifdef CH_MPI - MPI_Allreduce(&max_temp, &max_point, CH_SPACEDIM, MPI_DOUBLE, MPI_MAX, - Chombo_MPI::comm); - MPI_Allreduce(&min_temp, &min_point, CH_SPACEDIM, MPI_DOUBLE, MPI_MIN, - Chombo_MPI::comm); - // MPI_Allreduce(¢er, ¢er, CH_SPACEDIM, MPI_DOUBLE, MPI_SUM, - // Chombo_MPI::comm); // old method - // MPI_Allreduce(&idx, &idx_sum, 1, MPI_INT, MPI_MAX, Chombo_MPI::comm); -#else // serial - max_point = max_temp; - min_point = min_temp; - // idx_sum = idx; -#endif - - // CH_assert(idx_sum != 0); // old method - - for (unsigned i = 0; i < CH_SPACEDIM; ++i) - { - // center[i] /= idx_sum; // old method - center[i] = (max_point[i] + min_point[i]) / 2.; - } - - if (m_params.verbose > AHFinder::SOME) - { - pout() << "max_point: (" << max_point[0] << ", " << max_point[1] -#if CH_SPACEDIM == 3 - << "," << max_point[2] -#endif - << ")" << std::endl; - pout() << "min_point: (" << min_point[0] << ", " << min_point[1] -#if CH_SPACEDIM == 3 - << ", " << min_point[2] -#endif - << ")" << std::endl; - } - - update_old_centers(center); // set new - - if (m_params.verbose > AHFinder::MIN) - { - pout() << "center: (" << center[0] << ", " << center[1] -#if CH_SPACEDIM == 3 - << ", " << center[2] -#endif - << ")" << std::endl; - } - - return center; -} - -template -void ApparentHorizon::calculate_minmax_F() const -{ - if (!get_converged()) - return; - - double local_max = std::numeric_limits::min(); - double local_min = std::numeric_limits::max(); - if (AHFinder::is_rank_active()) - { - auto local_minmax = (std::minmax_element(m_F.begin(), m_F.end())); - local_min = *(local_minmax.first); - local_max = *(local_minmax.second); - } - double global_max = local_max; - double global_min = local_min; - -#ifdef CH_MPI - MPI_Allreduce(&local_min, &global_min, 1, MPI_DOUBLE, MPI_MIN, - Chombo_MPI::comm); - MPI_Allreduce(&local_max, &global_max, 1, MPI_DOUBLE, MPI_MAX, - Chombo_MPI::comm); -#endif - - m_max_F = global_max; - m_min_F = global_min; -} - -template -void ApparentHorizon::calculate_average_F() const -{ - if (!get_converged()) - return; - - int local_size = 0; - double local_sum = 0.; - double local_sum_sq = 0.; - if (AHFinder::is_rank_active()) - { - std::pair sums(0., 0.); - auto lambda_sum = [](std::pair sums, double r) { - sums.first += r; // radius - sums.second += r * r; // radius^2 - return std::move(sums); - }; - - auto local_sums = - std::accumulate(m_F.begin(), m_F.end(), sums, lambda_sum); - - local_size = m_F.size(); - local_sum = local_sums.first; - local_sum_sq = local_sums.second; - } - int global_size = local_size; - double global_sum = local_sum; - double global_sum_sq = local_sum_sq; - -#ifdef CH_MPI - MPI_Allreduce(&local_size, &global_size, 1, MPI_INT, MPI_SUM, - Chombo_MPI::comm); - MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, - Chombo_MPI::comm); - MPI_Allreduce(&global_sum_sq, &global_sum_sq, 1, MPI_DOUBLE, MPI_SUM, - Chombo_MPI::comm); -#endif - - m_ave_F = global_sum / global_size; - m_std_F = sqrt(global_sum_sq / global_size - - m_ave_F * m_ave_F); // sqrt( variance ) -} - -#endif // _APPARENTHORIZON_IMPL_HPP_ \ No newline at end of file diff --git a/Source/ApparentHorizonFinder/ApparentHorizon_petsc.impl.hpp b/Source/ApparentHorizonFinder/ApparentHorizon_petsc.impl.hpp deleted file mode 100644 index 91585d733..000000000 --- a/Source/ApparentHorizonFinder/ApparentHorizon_petsc.impl.hpp +++ /dev/null @@ -1,1001 +0,0 @@ -/* GRChombo - * Copyright 2012 The GRChombo collaboration. - * Please refer to LICENSE in GRChombo's root directory. - */ - -#ifndef _APPARENTHORIZON_HPP_ -#error "This file should only be included through ApparentHorizon.hpp" -#endif - -#ifndef _APPARENTHORIZON_PETSC_IMPL_HPP_ -#define _APPARENTHORIZON_PETSC_IMPL_HPP_ - -// petsc libraries -#include "petscsys.h" -#include "petscviewerhdf5.h" - -// for interpolation at restart: -#include "Lagrange.hpp" -#include "SimpleArrayBox.hpp" -#include "SimpleInterpSource.hpp" - -///////////////////////////////////////////////////////////////////////// -///////////////////////////////////////////////////////////////////////// -/////////////////////////// PETSc stuff below /////////////////////////// -///////////////////////////////////////////////////////////////////////// -///////////////////////////////////////////////////////////////////////// - -template -void ApparentHorizon::initialise_PETSc() -{ - CH_TIME("ApparentHorizon::initialise_PETSc"); - - if (!AHFinder::is_rank_active()) - return; - - const double lo_u = m_interp.get_domain_u_min(); - const double hi_u = m_interp.get_domain_u_max(); -#if CH_SPACEDIM == 3 // lo_v, hi_v not used otherwise - const double lo_v = m_interp.get_domain_v_min(); - const double hi_v = m_interp.get_domain_v_max(); -#endif - -#if PETSC_VERSION_LT(3, 5, 0) -#define DM_BOUNDARY_PERIODIC DMDA_BOUNDARY_PERIODIC -#define DM_BOUNDARY_GHOSTED DMDA_BOUNDARY_PERIODIC -#endif - -#if CH_SPACEDIM == 3 - DMDACreate2d(PETSC_COMM_WORLD, - m_periodic_u ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_GHOSTED, - m_periodic_v ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_GHOSTED, - DMDA_STENCIL_BOX, m_num_global_u, - m_num_global_v, /* grid size (negative means that we can - override it using command line) */ - PETSC_DECIDE, PETSC_DECIDE, /* distribution */ - 1, /* number of degrees of freedom */ - 3, /* stencil width (each side from central point) */ - NULL, NULL, &m_dmda); -#elif CH_SPACEDIM == 2 - DMDACreate1d(PETSC_COMM_WORLD, - m_periodic_u ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_GHOSTED, - m_num_global_u, /* grid size (negative means that we can - override it using command line) */ - 1, /* number of degrees of freedom */ - 3, /* stencil width (each side from central point) */ - NULL, &m_dmda); -#endif - - DMSetUp(m_dmda); - - // reload -> is it the same value? -#if CH_SPACEDIM == 3 // lo_v, hi_v not used otherwise - DMDAGetInfo(m_dmda, NULL, &m_num_global_u, &m_num_global_v, NULL, NULL, - NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); -#elif CH_SPACEDIM == 2 - DMDAGetInfo(m_dmda, NULL, &m_num_global_u, NULL, NULL, NULL, NULL, NULL, - NULL, NULL, NULL, NULL, NULL, NULL); -#endif - - m_du = - (hi_u - lo_u) / (m_periodic_u ? m_num_global_u : (m_num_global_u - 1)); - -#if CH_SPACEDIM == 3 - m_dv = - (hi_v - lo_v) / (m_periodic_v ? m_num_global_v : (m_num_global_v - 1)); -#endif - -#if CH_SPACEDIM == 3 - DMDAGetCorners(m_dmda, &m_umin, &m_vmin, NULL, &m_nu, &m_nv, NULL); - m_vmax = m_vmin + m_nv; -#elif CH_SPACEDIM == 2 - DMDAGetCorners(m_dmda, &m_umin, NULL, NULL, &m_nu, NULL, NULL); -#endif - m_umax = m_umin + m_nu; - -#if CH_SPACEDIM == 3 - int vec_size = m_nu * m_nv; -#elif CH_SPACEDIM == 2 - int vec_size = m_nu; -#endif - - m_F.resize(vec_size); - m_u.resize(vec_size); -#if CH_SPACEDIM == 3 - m_v.resize(vec_size); -#endif - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { -#if CH_SPACEDIM == 3 - m_u[(v - m_vmin) * m_nu + (u - m_umin)] = lo_u + u * m_du; - m_v[(v - m_vmin) * m_nu + (u - m_umin)] = lo_v + v * m_dv; -#elif CH_SPACEDIM == 2 - m_u[u - m_umin] = lo_u + u * m_du; -#endif - } - } - - DMCreateGlobalVector(m_dmda, &m_snes_soln); - PetscObjectSetName((PetscObject)m_snes_soln, "F"); - - VecDuplicate(m_snes_soln, &m_snes_rhs); - PetscObjectSetName((PetscObject)m_snes_rhs, "expansion"); - -#if PETSC_VERSION_GE(3, 5, 0) - DMSetMatType(m_dmda, MATAIJ); - DMCreateMatrix(m_dmda, &m_snes_jac); - MatSetOption(m_snes_jac, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); -#else - DMCreateMatrix(m_dmda, MATAIJ, &m_snes_jac); - MatSetOption(m_snes_jac, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); -#endif - - SNESCreate(PETSC_COMM_WORLD, &m_snes); - SNESSetFunction(m_snes, m_snes_rhs, &Petsc_form_function, this); - SNESSetJacobian(m_snes, m_snes_jac, m_snes_jac, &Petsc_form_jacobian, this); - SNESMonitorSet(m_snes, &Petsc_SNES_monitor, this, NULL); - - SNESSetFromOptions(m_snes); - - if (m_params.verbose > AHFinder::MIN) - { - SNESType snes_type; - SNESGetType(m_snes, &snes_type); - PetscReal snes_atol, snes_rtol, snes_stol; - PetscInt snes_maxit, snes_maxf; - SNESGetTolerances(m_snes, &snes_atol, &snes_rtol, &snes_stol, - &snes_maxit, &snes_maxf); - - PetscReal snes_divtol; - SNESGetDivergenceTolerance(m_snes, &snes_divtol); - - KSP snes_ksp; - SNESGetKSP(m_snes, &snes_ksp); - KSPType ksp_type; - KSPGetType(snes_ksp, &ksp_type); - PetscReal ksp_rtol, ksp_abstol, ksp_dtol; - PetscInt ksp_maxits; - KSPGetTolerances(snes_ksp, &ksp_rtol, &ksp_abstol, &ksp_dtol, - &ksp_maxits); - - pout() << "-------------------------------------\n"; - pout() << "ApparentHorizon Options:\n"; - pout() << "-------------------------------------\n"; - pout() << "PETSc SNES Options:\n"; - pout() << "Type: " << snes_type << "\n"; - pout() << "atol = " << snes_atol << ", rtol = " << snes_rtol - << ", stol = " << snes_stol << ",\n"; - pout() << "maxit = " << snes_maxit << ", maxf = " << snes_maxf << "\n"; - pout() << "divtol = " << snes_divtol << "\n"; - pout() << "-------------------------------------\n"; - pout() << "PETSc KSP Options:\n"; - pout() << "Type: " << ksp_type << "\n"; - pout() << "rtol = " << ksp_rtol << ", abstol = " << ksp_abstol - << ", dtol = " << ksp_dtol << ", maxits = " << ksp_maxits - << "\n"; - pout() << "-------------------------------------" << std::endl; - } -} - -template -void ApparentHorizon::finalise_PETSc() -{ - if (!AHFinder::is_rank_active()) - return; - - SNESDestroy(&m_snes); - VecDestroy(&m_snes_soln); - VecDestroy(&m_snes_rhs); - MatDestroy(&m_snes_jac); - DMDestroy(&m_dmda); -} - -template -bool ApparentHorizon::interpolate_ah( - dmda_arr_t f, const std::vector> &old_coords) -{ - // check if number of points changed - bool points_changed = false; - const double du = old_coords[0][1] - old_coords[0][0]; - const double lo_u = m_interp.get_domain_u_min(); - -#if CH_SPACEDIM == 3 - // assumes ordering in file is done first by fixed 'v', varying 'u' - // such that number of equal 'v's will be the number of 'u' points - int old_number_of_u = 1; - while (abs(old_coords[1][old_number_of_u - 1] - - old_coords[1][old_number_of_u]) < 1e-7) - ++old_number_of_u; - int old_number_of_v = old_coords[0].size() / old_number_of_u; - const double dv = - old_coords[1][old_number_of_u] - old_coords[1][old_number_of_u - 1]; - const double lo_v = m_interp.get_domain_v_min(); - // int total_new_points = m_num_global_u * m_num_global_v; - - points_changed |= (old_number_of_v != m_num_global_v); -#elif CH_SPACEDIM == 2 - int old_number_of_u = old_coords[0].size(); - // int total_new_points = m_num_global_u; -#endif - points_changed |= (old_number_of_u != m_num_global_u); - - if (!points_changed) - { - // these should be equal if all is good - // remember that only PETSc ranks have m_du and m_dv defined - // but only those should enter here, so all good -#if CH_SPACEDIM == 3 - CH_assert(abs(m_dv - dv) < 1e-7); - CH_assert(old_coords[0].size() == m_num_global_u * m_num_global_v); -#elif CH_SPACEDIM == 2 - CH_assert(old_coords[0].size() == m_num_global_u); -#endif - CH_assert(abs(m_du - du) < 1e-7); - - // select that F's that matter to this ranking -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { -#if CH_SPACEDIM == 3 - int idx_global = v * m_num_global_u + u; - int idx_local = (v - m_vmin) * m_nu + (u - m_umin); - if (std::abs((lo_v + v * m_dv) - old_coords[1][idx_global]) > - 1.e-7) - MayDay::Error( - "'v' coordinates incompatible with restart data"); -#elif CH_SPACEDIM == 2 - int idx_global = u; - int idx_local = (u - m_umin); -#endif - - if (std::abs((lo_u + u * m_du) - old_coords[0][idx_global]) > - 1.e-7) - MayDay::Error( - "'u' coordinates incompatible with restart data"); - - m_F[idx_local] = old_coords[CH_SPACEDIM - 1][idx_global]; - -#if CH_SPACEDIM == 3 - // as 'write_coords_file' preserves order - f[v][u] = m_F[idx_local]; -#elif CH_SPACEDIM == 2 - f[u] = m_F[idx_local]; -#endif - } - } - - return false; - } - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Number of AH points changed. Interpolating old ones." - << std::endl; - } - - // start interpolating - -#if CH_SPACEDIM == 3 - SimpleArrayBox box( - {old_number_of_u, old_number_of_v}, old_coords[CH_SPACEDIM - 1], - {m_interp.is_u_periodic(), m_interp.is_v_periodic()}); - SimpleInterpSource source( - {old_number_of_u, old_number_of_v}, - {m_interp.is_u_periodic(), m_interp.is_v_periodic()}); -#elif CH_SPACEDIM == 2 - SimpleArrayBox box({old_number_of_u}, - old_coords[CH_SPACEDIM - 1], - {m_interp.is_u_periodic()}); - SimpleInterpSource source({old_number_of_u}, - {m_interp.is_u_periodic()}); -#endif - - const int Order = 4; - bool verbose = false; - Lagrange interpolator4(source, verbose); - -#if CH_SPACEDIM == 3 - // local, no derivative - std::array derivs = {0, 0}; - std::array dxs = {du, dv}; -#elif CH_SPACEDIM == 2 - std::array derivs = {0}; - std::array dxs = {du}; -#endif - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) - { - double v_val = lo_v + v * m_dv; - double v_old_idx = v_val / dv; - // round such that we don't get AMRInterpolator errors - // (e.g. the last point idx=N-1 should be exact, and not being exact was - // generating this error in the Lagrange intrepolator trying to access - // cells beyond the last cell, and hence error... this fixed it) - if (v_old_idx - int(v_old_idx) < 1.e-7) - v_old_idx = int(v_old_idx); -#else - { -#endif - for (int u = m_umin; u < m_umax; ++u) - { - double u_val = lo_u + u * m_du; - double u_old_idx = u_val / du; - // round such that we don't get AMRInterpolator errors (same as for - // 'v') - if (u_old_idx - int(u_old_idx) < 1.e-7) - u_old_idx = int(u_old_idx); - -#if CH_SPACEDIM == 3 - // int idx_global = v * m_num_global_u + u; - int idx_local = (v - m_vmin) * m_nu + (u - m_umin); - - std::array evalCoord = {u_old_idx, - v_old_idx}; -#elif CH_SPACEDIM == 2 - // int idx_global = u; - int idx_local = (u - m_umin); - std::array evalCoord = {u_old_idx}; -#endif - - interpolator4.setup(derivs, dxs, evalCoord); - m_F[idx_local] = interpolator4.interpData(box); - -#if CH_SPACEDIM == 3 - // as 'write_coords_file' preserves order - f[v][u] = m_F[idx_local]; -#elif CH_SPACEDIM == 2 - f[u] = m_F[idx_local]; -#endif - } - } - - if (m_params.verbose > AHFinder::NONE) - { - pout() << "Interpolation Successfull." << std::endl; - } - - return true; -} - -template -void ApparentHorizon::set_stencils(AHDeriv &out, - int u -#if CH_SPACEDIM == 3 - , - int v -#endif -) -{ -#if CH_SPACEDIM == 3 - const bool periodic[CH_SPACEDIM - 1] = {m_periodic_u, m_periodic_v}; - - int coord[CH_SPACEDIM - 1] = {u, v}; - const int N[CH_SPACEDIM - 1] = {m_num_global_u, m_num_global_v}; - - int d_start[CH_SPACEDIM - 1]; - double *const d_weights[CH_SPACEDIM - 1] = {out.du_weights, out.dv_weights}; - - int dd_start[CH_SPACEDIM - 1]; - double *const dd_weights[CH_SPACEDIM - 1] = {out.dudu_weights, - out.dvdv_weights}; -#elif CH_SPACEDIM == 2 - const bool periodic[CH_SPACEDIM - 1] = {m_periodic_u}; - - int coord[CH_SPACEDIM - 1] = {u}; - const int N[CH_SPACEDIM - 1] = {m_num_global_u}; - - int d_start[CH_SPACEDIM - 1]; - double *const d_weights[CH_SPACEDIM - 1] = {out.du_weights}; - - int dd_start[CH_SPACEDIM - 1]; - double *const dd_weights[CH_SPACEDIM - 1] = {out.dudu_weights}; -#endif - - for (int i = 0; i < CH_SPACEDIM - 1; ++i) - { - - // Non-periodic boundaries - if (!periodic[i] && (coord[i] == 0 || coord[i] == 1 || - coord[i] == N[i] - 2 || coord[i] == N[i] - 1)) - { - - if (coord[i] == 0) - { - - d_start[i] = 0; - - d_weights[i][0] = -2.083333333333333333333; - d_weights[i][1] = +4.000000000000000000000; - d_weights[i][2] = -3.000000000000000000000; - d_weights[i][3] = +1.333333333333333333333; - d_weights[i][4] = -0.250000000000000000000; - - dd_start[i] = 0; - - dd_weights[i][0] = +3.75000000000000000000; - dd_weights[i][1] = -12.8333333333333333333; - dd_weights[i][2] = +17.8333333333333333333; - dd_weights[i][3] = -13.0000000000000000000; - dd_weights[i][4] = +5.08333333333333333333; - dd_weights[i][5] = -0.83333333333333333333; - } - else if (coord[i] == 1) - { - - d_start[i] = -1; - - d_weights[i][0] = -0.250000000000000000000; - d_weights[i][1] = -0.833333333333333333333; - d_weights[i][2] = +1.500000000000000000000; - d_weights[i][3] = -0.500000000000000000000; - d_weights[i][4] = +0.083333333333333333333; - - dd_start[i] = -1; - - dd_weights[i][0] = +0.83333333333333333333; - dd_weights[i][1] = -1.25000000000000000000; - dd_weights[i][2] = -0.33333333333333333333; - dd_weights[i][3] = +1.16666666666666666667; - dd_weights[i][4] = -0.50000000000000000000; - dd_weights[i][5] = +0.08333333333333333333; - } - else if (coord[i] == N[i] - 2) - { - - d_start[i] = -3; - - d_weights[i][0] = -0.083333333333333333333; - d_weights[i][1] = +0.500000000000000000000; - d_weights[i][2] = -1.500000000000000000000; - d_weights[i][3] = +0.833333333333333333333; - d_weights[i][4] = +0.250000000000000000000; - - dd_start[i] = -4; - - dd_weights[i][0] = +0.08333333333333333333; - dd_weights[i][1] = -0.50000000000000000000; - dd_weights[i][2] = +1.16666666666666666667; - dd_weights[i][3] = -0.33333333333333333333; - dd_weights[i][4] = -1.25000000000000000000; - dd_weights[i][5] = +0.83333333333333333333; - } - else if (coord[i] == N[i] - 1) - { - - d_start[i] = -4; - - d_weights[i][0] = +0.250000000000000000000; - d_weights[i][1] = -1.333333333333333333333; - d_weights[i][2] = +3.000000000000000000000; - d_weights[i][3] = -4.000000000000000000000; - d_weights[i][4] = +2.083333333333333333333; - - dd_start[i] = -5; - - dd_weights[i][0] = -0.83333333333333333333; - dd_weights[i][1] = +5.08333333333333333333; - dd_weights[i][2] = -13.0000000000000000000; - dd_weights[i][3] = +17.8333333333333333333; - dd_weights[i][4] = -12.8333333333333333333; - dd_weights[i][5] = +3.75000000000000000000; - } - } - // Standard points - else - { - - d_start[i] = -2; - - d_weights[i][0] = +0.083333333333333333333; - d_weights[i][1] = -0.666666666666666666666; - d_weights[i][2] = 0; - d_weights[i][3] = +0.666666666666666666666; - d_weights[i][4] = -0.083333333333333333333; - - if (coord[i] == N[i] - 3) - { - dd_start[i] = -3; - - dd_weights[i][0] = 0; - dd_weights[i][1] = -0.08333333333333333333; - dd_weights[i][2] = +1.33333333333333333333; - dd_weights[i][3] = -2.50000000000000000000; - dd_weights[i][4] = +1.33333333333333333333; - dd_weights[i][5] = -0.08333333333333333333; - } - else - { - dd_start[i] = -2; - - dd_weights[i][0] = -0.08333333333333333333; - dd_weights[i][1] = +1.33333333333333333333; - dd_weights[i][2] = -2.50000000000000000000; - dd_weights[i][3] = +1.33333333333333333333; - dd_weights[i][4] = -0.08333333333333333333; - dd_weights[i][5] = 0; - } - } - } - - out.du_stencil_start = d_start[0]; - out.dudu_stencil_start = dd_start[0]; - -#if CH_SPACEDIM == 3 - out.dv_stencil_start = d_start[1]; - out.dvdv_stencil_start = dd_start[1]; -#endif -} - -template -AHDeriv ApparentHorizon::diff(const dmda_arr_t in, - int u -#if CH_SPACEDIM == 3 - , - int v -#endif -) -{ - CH_TIME("ApparentHorizon::diff"); - - AHDeriv out; - set_stencils(out, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - -#if CH_SPACEDIM == 3 - // d/du - for (int j = 0; j < DWIDTH; ++j) - { - out.duF += - out.du_weights[j] * in[v][u + out.du_stencil_start + j] / m_du; - - // d2/dudv - for (int k = 0; k < DWIDTH; ++k) - { - out.dudvF += - out.du_weights[j] * out.dv_weights[k] * - in[v + out.dv_stencil_start + k][u + out.du_stencil_start + j] / - (m_du * m_dv); - } - } - - // d/dv - for (int j = 0; j < DWIDTH; ++j) - { - out.dvF += - out.dv_weights[j] * in[v + out.dv_stencil_start + j][u] / m_dv; - } - - // d2/du2 - for (int j = 0; j < DDWIDTH; ++j) - { - out.duduF += out.dudu_weights[j] * - in[v][u + out.dudu_stencil_start + j] / (m_du * m_du); - } - - // d2/dv2 - for (int j = 0; j < DDWIDTH; ++j) - { - out.dvdvF += out.dvdv_weights[j] * - in[v + out.dvdv_stencil_start + j][u] / (m_dv * m_dv); - } - -#elif CH_SPACEDIM == 2 - - // d/du - for (int j = 0; j < DWIDTH; ++j) - { - out.duF += out.du_weights[j] * in[u + out.du_stencil_start + j] / m_du; - } - - // d2/du2 - for (int j = 0; j < DDWIDTH; ++j) - { - out.duduF += out.dudu_weights[j] * in[u + out.dudu_stencil_start + j] / - (m_du * m_du); - } - -#endif - - return out; -} - -template -void ApparentHorizon::form_function(Vec F, Vec Rhs) -{ - CH_TIME("ApparentHorizon::form_function"); - - // Scatter ghost cells - Vec localF; - DMGetLocalVector(m_dmda, &localF); - DMGlobalToLocalBegin(m_dmda, F, INSERT_VALUES, localF); - DMGlobalToLocalEnd(m_dmda, F, INSERT_VALUES, localF); - - dmda_arr_t in; - DMDAVecGetArray(m_dmda, localF, &in); - - dmda_arr_t out; - DMDAVecGetArray(m_dmda, Rhs, &out); - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { -#if CH_SPACEDIM == 3 - m_F[(v - m_vmin) * m_nu + (u - m_umin)] = in[v][u]; -#elif CH_SPACEDIM == 2 - m_F[u - m_umin] = in[u]; -#endif - } - } - - bool out_of_grid = m_interp.set_coordinates(m_F, m_u -#if CH_SPACEDIM == 3 - , - m_v -#endif - ); - // abort if out of grid - reduces the time in divergence dramatically! - if (out_of_grid) - SNESSetFunctionDomainError(m_snes); - - // must interpolate anyway, as other PETSc ranks do not know if this one is - // "out_of_grid" - m_interp.interpolate(); - - int idx = 0; - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { - -#if CH_SPACEDIM == 3 - double &_out = out[v][u]; -#elif CH_SPACEDIM == 2 - double &_out = out[u]; -#endif - - AHDeriv deriv = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - if (out_of_grid) // no need to calculate expansion - _out = 0.; - else if (!m_periodic_u && (u == 0 || u == m_num_global_u - 1)) - _out = deriv.duF; - -#if CH_SPACEDIM == 3 - else if (!m_periodic_v && (v == 0 || v == m_num_global_v - 1)) - _out = deriv.dvF; -#endif - - else - { - const auto geometry_data = m_interp.get_geometry_data(idx); - const auto data = m_interp.get_data(idx); - const auto coords = m_interp.get_coords(idx); - const auto coords_cart = m_interp.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - _out = func.get(geometry_data, deriv, m_func_params); - } - - ++idx; - } - } - - DMDAVecRestoreArray(m_dmda, localF, &in); - DMDAVecRestoreArray(m_dmda, Rhs, &out); - DMRestoreLocalVector(m_dmda, &localF); -} - -template -void ApparentHorizon::form_jacobian(Vec F, Mat J) -{ - CH_TIME("ApparentHorizon::form_jacobian"); - - // Scatter ghost cells - Vec localF; - DMGetLocalVector(m_dmda, &localF); - DMGlobalToLocalBegin(m_dmda, F, INSERT_VALUES, localF); - DMGlobalToLocalEnd(m_dmda, F, INSERT_VALUES, localF); - - dmda_arr_t in; - DMDAVecGetArray(m_dmda, localF, &in); - - bool out_of_grid = m_interp_plus.set_coordinates(m_F, m_u, -#if CH_SPACEDIM == 3 - m_v, -#endif - eps); - out_of_grid |= m_interp_minus.set_coordinates(m_F, m_u, -#if CH_SPACEDIM == 3 - m_v, -#endif - -eps); - // abort if out of grid - reduces the time in divergence dramatically! - if (out_of_grid) - SNESSetFunctionDomainError(m_snes); - - // must interpolate anyway, as other PETSc ranks do not know if this one is - // "out_of_grid" - m_interp_plus.interpolate(); - m_interp_minus.interpolate(); - - int idx = 0; - -#if CH_SPACEDIM == 3 - for (int v = m_vmin; v < m_vmax; ++v) -#endif - { - for (int u = m_umin; u < m_umax; ++u) - { - MatStencil row[1] = {0}; - row[0].i = u; -#if CH_SPACEDIM == 3 - row[0].j = v; -#endif - - if (out_of_grid) // no need to calculate expansion - { - MatStencil col[DWIDTH] = {0}; - double val[DWIDTH] = {0}; - MatSetValuesStencil(J, 1, row, DWIDTH, col, val, INSERT_VALUES); - } - else if (!m_periodic_u && (u == 0 || u == m_num_global_u - 1)) - { - MatStencil col[DWIDTH] = {0}; - double val[DWIDTH] = {0}; - - const AHDeriv deriv = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - for (int a = 0; a < DWIDTH; ++a) - { - col[a].i = u + deriv.du_stencil_start + a; -#if CH_SPACEDIM == 3 - col[a].j = v; -#endif - val[a] = deriv.du_weights[a]; - } - - MatSetValuesStencil(J, 1, row, DWIDTH, col, val, INSERT_VALUES); - } - -#if CH_SPACEDIM == 3 - else if (!m_periodic_v && (v == 0 || v == m_num_global_v - 1)) - { - MatStencil col[DWIDTH] = {0}; - double val[DWIDTH] = {0}; - - const AHDeriv deriv = diff(in, u, v); - - for (int b = 0; b < DWIDTH; ++b) - { - col[b].i = u; - col[b].j = v + deriv.dv_stencil_start + b; - val[b] = deriv.dv_weights[b]; - } - - MatSetValuesStencil(J, 1, row, DWIDTH, col, val, INSERT_VALUES); - } -#endif - else - { - -#if CH_SPACEDIM == 3 - const int NVAL = DDWIDTH * DDWIDTH; -#elif CH_SPACEDIM == 2 - const int NVAL = DDWIDTH; -#endif - MatStencil col[NVAL] = {0}; - double val[NVAL] = {0}; - - const AHDeriv deriv_default = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - // "local" and "stencil" jacobian terms - // "local" (point {u,v}) corresponds to (a == - // -deriv_default.dudu_stencil_start - // && b == -deriv_default.dvdv_stencil_start) and has to use - // 'm_interp_plus' and 'm_interp_minus' - { -#if CH_SPACEDIM == 3 - for (int b = 0; b < DDWIDTH; ++b) - { -#elif CH_SPACEDIM == 2 - { - int b = 0; -#endif - for (int a = 0; a < DDWIDTH; ++a) - { - col[b * DDWIDTH + a].i = - u + deriv_default.dudu_stencil_start + a; - const int uu = col[b * DDWIDTH + a].i; - const int u_start = - -deriv_default.dudu_stencil_start; -#if CH_SPACEDIM == 3 - col[b * DDWIDTH + a].j = - v + deriv_default.dvdv_stencil_start + b; - const int vv = col[b * DDWIDTH + a].j; - const int v_start = - -deriv_default.dvdv_stencil_start; -#elif CH_SPACEDIM == 2 - const int v_start = 0; -#endif - - if (a == u_start && b == v_start) // <=> uu=u, vv=v - { - val[b * DDWIDTH + a] = point_jacobian( - u, uu, -#if CH_SPACEDIM == 3 - v, vv, -#endif - in, idx, m_interp_plus, m_interp_minus); - } - else - { - val[b * DDWIDTH + a] = - point_jacobian(u, uu, -#if CH_SPACEDIM == 3 - v, vv, -#endif - in, idx, m_interp, m_interp); - } - } - } - } - - MatSetValuesStencil(J, 1, row, NVAL, col, val, INSERT_VALUES); - } - - ++idx; - } - } - - DMDAVecRestoreArray(m_dmda, localF, &in); - DMRestoreLocalVector(m_dmda, &localF); - - MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); -} - -template -double ApparentHorizon::point_jacobian( - int u, int u_stencil, -#if CH_SPACEDIM == 3 - int v, int v_stencil, -#endif - dmda_arr_t in, int idx, - const AHInterpolation &interp_plus, - const AHInterpolation &interp_minus) -{ -#if CH_SPACEDIM == 3 - double &_in = in[v_stencil][u_stencil]; -#elif CH_SPACEDIM == 2 - double &_in = in[u_stencil]; -#endif - - // "plus" perturbation - double expansionPlus; - { - double in_old = _in; - _in += eps; // perturb just the point {u,v} - - const AHDeriv deriv = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - const auto geometry_data = interp_plus.get_geometry_data(idx); - const auto data = interp_plus.get_data(idx); - const auto coords = interp_plus.get_coords(idx); - const auto coords_cart = interp_plus.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - expansionPlus = func.get(geometry_data, deriv, m_func_params); - - _in = in_old; - } - - // "minus" perturbation - double expansionMinus; - { - double in_old = _in; - _in -= eps; // perturb just the point {u,v} - - const AHDeriv deriv = diff(in, u -#if CH_SPACEDIM == 3 - , - v -#endif - ); - - const auto geometry_data = interp_minus.get_geometry_data(idx); - const auto data = interp_minus.get_data(idx); - const auto coords = interp_minus.get_coords(idx); - const auto coords_cart = interp_minus.get_cartesian_coords(idx); - AHFunction func(data, coords, coords_cart); - expansionMinus = func.get(geometry_data, deriv, m_func_params); - - _in = in_old; - } - - return (expansionPlus - expansionMinus) / (2. * eps); -} - -//! functions used by PETSc based on 'form_function' and 'form_jacobian' -template -PetscErrorCode -ApparentHorizon::Petsc_form_function(SNES snes, - Vec F, - Vec Rhs, - void *ptr) -{ - ApparentHorizon &ah = *reinterpret_cast(ptr); - CH_assert(ah.m_snes == snes); - ah.form_function(F, Rhs); - return 0; -} - -template -PetscErrorCode -#if PETSC_VERSION_GE(3, 5, 0) -ApparentHorizon::Petsc_form_jacobian( - SNES snes, Vec F, Mat Amat, Mat Pmat, void *ptr) -#else -ApparentHorizon::Petsc_form_jacobian( - SNES snes, Vec F, Mat *Amat, Mat *Pmat, MatStructure *flag, void *ptr) -#endif -{ - ApparentHorizon &ah = *reinterpret_cast(ptr); - -#if PETSC_VERSION_GE(3, 5, 0) - CH_assert(ah.m_snes == snes && Amat == Pmat); - ah.form_jacobian(F, Amat); -#else - CH_assert(ah.m_snes == snes && *Amat == *Pmat); - ah.form_jacobian(F, *Amat); -#endif - - return 0; -} - -template -PetscErrorCode ApparentHorizon::Petsc_SNES_monitor( - SNES snes, PetscInt its, PetscReal norm, void *ptr) -{ - ApparentHorizon &ah = *reinterpret_cast(ptr); - CH_assert(ah.m_snes == snes); - return 0; -} - -#endif // _APPARENTHORIZON_PETSC_IMPL_HPP_ diff --git a/Source/BlackHoles/BHAMR.hpp b/Source/BlackHoles/BHAMR.hpp index 9febe1ddc..f5f33e79e 100644 --- a/Source/BlackHoles/BHAMR.hpp +++ b/Source/BlackHoles/BHAMR.hpp @@ -7,13 +7,7 @@ #define BHAMR_HPP_ #include "GRAMR.hpp" -#if CH_SPACEDIM == 3 #include "PunctureTracker.hpp" -#endif - -#ifdef USE_AHFINDER -#include "AHFinder.hpp" -#endif /// A child of Chombo's AMR class to interface with tools which require /// access to the whole AMR hierarchy, and those of GRAMR @@ -23,25 +17,14 @@ class BHAMR : public GRAMR { public: -#if CH_SPACEDIM == 3 PunctureTracker m_puncture_tracker; -#endif - -#ifdef USE_AHFINDER - AHFinder m_ah_finder; -#endif BHAMR() {} void set_interpolator(AMRInterpolator> *a_interpolator) override { GRAMR::set_interpolator(a_interpolator); -#if CH_SPACEDIM == 3 m_puncture_tracker.set_interpolator(a_interpolator); -#endif -#ifdef USE_AHFINDER - m_ah_finder.set_interpolator(a_interpolator); -#endif } }; diff --git a/Source/BlackHoles/PunctureTracker.cpp b/Source/BlackHoles/PunctureTracker.cpp index 5ad520cc6..790e76d68 100644 --- a/Source/BlackHoles/PunctureTracker.cpp +++ b/Source/BlackHoles/PunctureTracker.cpp @@ -3,8 +3,6 @@ * Please refer to LICENSE in GRChombo's root directory. */ -#if CH_SPACEDIM == 3 - #include "PunctureTracker.hpp" #include "ChomboParameters.hpp" // for writing data #include "DimensionDefinitions.hpp" @@ -15,9 +13,13 @@ //! Set punctures post restart void PunctureTracker::initial_setup( const std::vector> &initial_puncture_coords, - const std::string &a_checkpoint_prefix, const int a_min_level) + const std::string &a_filename, const std::string &a_output_path, + const int a_min_level) { - m_punctures_filename = a_checkpoint_prefix + "Punctures"; + if (!FilesystemTools::directory_exists(a_output_path)) + FilesystemTools::mkdir_recursive(a_output_path); + + m_punctures_filename = a_output_path + a_filename; // first set the puncture data // m_num_punctures is only set later @@ -56,7 +58,7 @@ void PunctureTracker::set_initial_punctures() for (int ipuncture = 0; ipuncture < m_num_punctures; ipuncture++) { // assume initial shift is always zero - FOR1(i) { m_puncture_shift[ipuncture][i] = 0.0; } + FOR(i) { m_puncture_shift[ipuncture][i] = 0.0; } } // now the write out to a new file @@ -151,7 +153,7 @@ void PunctureTracker::execute_tracking(double a_time, double a_restart_time, // update puncture locations using second order update for (int ipuncture = 0; ipuncture < m_num_punctures; ipuncture++) { - FOR1(i) + FOR(i) { m_puncture_coords[ipuncture][i] += -0.5 * a_dt * @@ -239,5 +241,3 @@ std::vector PunctureTracker::get_puncture_vector() const } return puncture_vector; } - -#endif // #if CH_SPACEDIM == 3 diff --git a/Source/BlackHoles/PunctureTracker.hpp b/Source/BlackHoles/PunctureTracker.hpp index ca5d175c1..8d79b2c40 100644 --- a/Source/BlackHoles/PunctureTracker.hpp +++ b/Source/BlackHoles/PunctureTracker.hpp @@ -36,7 +36,8 @@ class PunctureTracker //! if the puncture locations are required for Tagging Criteria void initial_setup(const std::vector> &initial_puncture_coords, - const std::string &a_checkpoint_prefix, + const std::string &a_filename = "punctures", + const std::string &a_output_path = "", const int a_min_level = 0); //! set puncture locations on start (or restart) diff --git a/Source/BoxUtils/BoxLoops.impl.hpp b/Source/BoxUtils/BoxLoops.impl.hpp index e8767df61..2d9e8432c 100644 --- a/Source/BoxUtils/BoxLoops.impl.hpp +++ b/Source/BoxUtils/BoxLoops.impl.hpp @@ -28,20 +28,20 @@ void BoxLoops::innermost_loop(const ComputePack &compute_pack, // SIMD LOOP #ifdef __INTEL_COMPILER #pragma novector -#else +#elif !defined(__clang__) #pragma omp simd safelen(1) #endif /* __INTEL_COMPILER */ for (int ix = loop_lo_x; ix <= x_simd_max; ix += simd_width) { compute_pack.call_compute( - Cell>(IntVect(D_DECL(ix, iy, iz)), box_pointers)); + Cell>(IntVect(ix, iy, iz), box_pointers)); } // REMAINDER LOOP for (int ix = x_simd_max + simd::simd_len; ix <= loop_hi_x; ++ix) { compute_pack.call_compute( - Cell(IntVect(D_DECL(ix, iy, iz)), box_pointers)); + Cell(IntVect(ix, iy, iz), box_pointers)); } } @@ -54,7 +54,7 @@ void BoxLoops::innermost_loop(const ComputePack &compute_pack, for (int ix = loop_lo_x; ix <= loop_hi_x; ++ix) { compute_pack.call_compute( - Cell(IntVect(D_DECL(ix, iy, iz)), box_pointers)); + Cell(IntVect(ix, iy, iz), box_pointers)); } } @@ -71,9 +71,6 @@ void BoxLoops::loop(const ComputePack &compute_pack, const int *loop_lo = loop_box.loVect(); const int *loop_hi = loop_box.hiVect(); -#if CH_SPACEDIM == 2 - int iz = 0; -#endif #pragma omp parallel for default(shared) collapse(CH_SPACEDIM - 1) #if CH_SPACEDIM >= 3 for (int iz = loop_lo[2]; iz <= loop_hi[2]; ++iz) diff --git a/Source/BoxUtils/BoxPointers.hpp b/Source/BoxUtils/BoxPointers.hpp index d2c2cdf36..120b156c2 100644 --- a/Source/BoxUtils/BoxPointers.hpp +++ b/Source/BoxUtils/BoxPointers.hpp @@ -74,16 +74,16 @@ class BoxPointers CellIndexIn get_in_index(IntVect integer_coords) const { - return (D_TERM((integer_coords[0] - m_in_lo[0]), - +m_in_stride[1] * (integer_coords[1] - m_in_lo[1]), - +m_in_stride[2] * (integer_coords[2] - m_in_lo[2]))); + return (m_in_stride[2] * (integer_coords[2] - m_in_lo[2]) + + m_in_stride[1] * (integer_coords[1] - m_in_lo[1]) + + (integer_coords[0] - m_in_lo[0])); } CellIndexOut get_out_index(IntVect integer_coords) const { - return (D_TERM((integer_coords[0] - m_out_lo[0]), - +m_out_stride[1] * (integer_coords[1] - m_out_lo[1]), - +m_out_stride[2] * (integer_coords[2] - m_out_lo[2]))); + return (m_out_stride[2] * (integer_coords[2] - m_out_lo[2]) + + m_out_stride[1] * (integer_coords[1] - m_out_lo[1]) + + (integer_coords[0] - m_out_lo[0])); } }; diff --git a/Source/BoxUtils/Cell.impl.hpp b/Source/BoxUtils/Cell.impl.hpp index 020bd00dc..ff9eef4f2 100644 --- a/Source/BoxUtils/Cell.impl.hpp +++ b/Source/BoxUtils/Cell.impl.hpp @@ -40,9 +40,9 @@ template template