From 11f9a2a206d01ff7605f5d98995a7b7a74ece211 Mon Sep 17 00:00:00 2001 From: cjoana Date: Tue, 4 Jul 2023 23:50:43 +0800 Subject: [PATCH 1/3] evo working, Weyl4 needed --- Source/CCZ4/Weyl4.hpp | 53 +++++-- Source/CCZ4/Weyl4.hpp_old | 102 ++++++++++++ Source/CCZ4/Weyl4.impl.hpp | 207 +++++++++++++++++------- Source/CCZ4/Weyl4.impl.hpp_old | 247 +++++++++++++++++++++++++++++ Source/Matter/MatterWeyl4.hpp | 54 +++++++ Source/Matter/MatterWeyl4.impl.hpp | 78 +++++++++ 6 files changed, 672 insertions(+), 69 deletions(-) create mode 100644 Source/CCZ4/Weyl4.hpp_old create mode 100644 Source/CCZ4/Weyl4.impl.hpp_old create mode 100644 Source/Matter/MatterWeyl4.hpp create mode 100644 Source/Matter/MatterWeyl4.impl.hpp diff --git a/Source/CCZ4/Weyl4.hpp b/Source/CCZ4/Weyl4.hpp index 46f2de58b..4187f89b4 100644 --- a/Source/CCZ4/Weyl4.hpp +++ b/Source/CCZ4/Weyl4.hpp @@ -6,8 +6,7 @@ #ifndef WEYL4_HPP_ #define WEYL4_HPP_ -#include "BSSNVars.hpp" -#include "CCZ4Geometry.hpp" +#include "CCZ4.hpp" #include "Cell.hpp" #include "Coordinates.hpp" #include "FourthOrderDerivatives.hpp" @@ -17,6 +16,16 @@ #include "simd.hpp" #include +// CJ addition +//! Struct the Weyl curvature and Chern-Pontryagin invariant. +template struct WeylScalars_t +{ + data_t Weyl_curv; // Weyl curvature + data_t ChP_inv; // Chern-Pontryagin invariant +}; + + + //! Struct for the E and B fields template struct EBFields_t { @@ -51,16 +60,19 @@ class Weyl4 { public: // Use the variable definitions containing the needed quantities - template using Vars = BSSNVars::VarsWithGauge; + template using Vars = CCZ4Vars::VarsWithGauge; template using Diff2Vars = ADMConformalVars::Diff2VarsNoGauge; //! Constructor of class Weyl4 /*! - Takes in the centre for the calculation of the tetrads, and grid spacing + Takes in the centre for the calculation of the tetrads, grid spacing and + the formulation. */ - Weyl4(const std::array a_center, const double a_dx) - : m_center(a_center), m_dx(a_dx), m_deriv(a_dx) + Weyl4(const std::array a_center, const double a_dx, + const int a_formulation = CCZ4RHS<>::USE_CCZ4) + : m_center(a_center), m_dx(a_dx), m_deriv(a_dx), + m_formulation(a_formulation) { } @@ -72,6 +84,12 @@ class Weyl4 const std::array m_center; //!< The grid center const double m_dx; //!< the grid spacing const FourthOrderDerivatives m_deriv; //!< for calculating derivs of vars + const int m_formulation; //!< CCZ4 or BSSN? + + //! Compute spatial volume element + template + Tensor<3, data_t> compute_epsilon3_LUU(const Vars &vars, + const Tensor<2, data_t> &h_UU) const; //! Calculation of Weyl_4 scalar template @@ -79,22 +97,33 @@ class Weyl4 const Vars &vars, const Vars> &d1, const Diff2Vars> &d2, + const Tensor<2, data_t> &h_UU, const Coordinates &coords) const; + //! Calculation of the Weyl Scalars + template + WeylScalars_t compute_WeylCurvature(const EBFields_t &ebfields, + const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Tensor<2, data_t> &h_UU, + const Coordinates &coords) const; + //! Calculation of the tetrads template Tetrad_t - compute_null_tetrad(const Vars &vars, + compute_null_tetrad(const Vars &vars, const Tensor<2, data_t> &h_UU, const Coordinates &coords) const; //! Calulation of the decomposition of the Weyl tensor in Electric and //! Magnetic fields template - EBFields_t - compute_EB_fields(const Vars &vars, - const Vars> &d1, - const Diff2Vars> &d2, - const Coordinates &coords) const; + EBFields_t compute_EB_fields(const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Tensor<3, data_t> &epsilon3_LUU, + const Tensor<2, data_t> &h_UU, + const chris_t &chris) const; }; #include "Weyl4.impl.hpp" diff --git a/Source/CCZ4/Weyl4.hpp_old b/Source/CCZ4/Weyl4.hpp_old new file mode 100644 index 000000000..46f2de58b --- /dev/null +++ b/Source/CCZ4/Weyl4.hpp_old @@ -0,0 +1,102 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef WEYL4_HPP_ +#define WEYL4_HPP_ + +#include "BSSNVars.hpp" +#include "CCZ4Geometry.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "FourthOrderDerivatives.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs c_NUM - total number of components +#include "simd.hpp" +#include + +//! Struct for the E and B fields +template struct EBFields_t +{ + Tensor<2, data_t> E; //!< Electric component of Weyltensor + Tensor<2, data_t> B; //!< Magnetic component of Weyltensor +}; + +//! Struct for the null tetrad +template struct Tetrad_t +{ + Tensor<1, data_t> u; //!< the vector u^i + Tensor<1, data_t> v; //!< the vector v^i + Tensor<1, data_t> w; //!< the vector w^i +}; + +//! Struct for the Newman Penrose scalar +template struct NPScalar_t +{ + data_t Real; // Real component + data_t Im; // Imaginary component +}; + +//! Calculates the Weyl4 scalar for spacetimes without matter content +/*! + This class calculates the Weyl4 scalar real and im parts using definitions + from Alcubierres book "Introduction to 3+1 Numerical Relativity". We use a + decomposition of the Weyl tensor in electric and magnetic parts, which then + is used together with the tetrads defined in "gr-qc/0104063" to calculate the + Weyl4 scalar. +*/ +class Weyl4 +{ + public: + // Use the variable definitions containing the needed quantities + template using Vars = BSSNVars::VarsWithGauge; + template + using Diff2Vars = ADMConformalVars::Diff2VarsNoGauge; + + //! Constructor of class Weyl4 + /*! + Takes in the centre for the calculation of the tetrads, and grid spacing + */ + Weyl4(const std::array a_center, const double a_dx) + : m_center(a_center), m_dx(a_dx), m_deriv(a_dx) + { + } + + //! The compute member which calculates the wave quantities at each point on + //! the grid + template void compute(Cell current_cell) const; + + protected: + const std::array m_center; //!< The grid center + const double m_dx; //!< the grid spacing + const FourthOrderDerivatives m_deriv; //!< for calculating derivs of vars + + //! Calculation of Weyl_4 scalar + template + NPScalar_t compute_Weyl4(const EBFields_t &ebfields, + const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Coordinates &coords) const; + + //! Calculation of the tetrads + template + Tetrad_t + compute_null_tetrad(const Vars &vars, + const Coordinates &coords) const; + + //! Calulation of the decomposition of the Weyl tensor in Electric and + //! Magnetic fields + template + EBFields_t + compute_EB_fields(const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Coordinates &coords) const; +}; + +#include "Weyl4.impl.hpp" + +#endif /* WEYL4_HPP_ */ diff --git a/Source/CCZ4/Weyl4.impl.hpp b/Source/CCZ4/Weyl4.impl.hpp index 23a94ffb5..1c0b96eb7 100644 --- a/Source/CCZ4/Weyl4.impl.hpp +++ b/Source/CCZ4/Weyl4.impl.hpp @@ -21,88 +21,117 @@ template void Weyl4::compute(Cell current_cell) const // Get the coordinates const Coordinates coords(current_cell, m_dx, m_center); + // Compute the inverse metric and Christoffel symbols + using namespace TensorAlgebra; + const auto h_UU = compute_inverse_sym(vars.h); + const auto chris = compute_christoffel(d1.h, h_UU); + + // Compute the spatial volume element + const auto epsilon3_LUU = compute_epsilon3_LUU(vars, h_UU); + // Compute the E and B fields - EBFields_t ebfields = compute_EB_fields(vars, d1, d2, coords); + EBFields_t ebfields = + compute_EB_fields(vars, d1, d2, epsilon3_LUU, h_UU, chris); // work out the Newman Penrose scalar - NPScalar_t out = compute_Weyl4(ebfields, vars, d1, d2, coords); + NPScalar_t out = + compute_Weyl4(ebfields, vars, d1, d2, h_UU, coords); // Write the rhs into the output FArrayBox current_cell.store_vars(out.Real, c_Weyl4_Re); current_cell.store_vars(out.Im, c_Weyl4_Im); + + // work out the Newman Penrose scalar + WeylScalars_t out = + compute_WeylCurvature(ebfields, vars, d1, d2, h_UU, coords); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(out.Weyl_curv, c_Weyl_curv); + current_cell.store_vars(out.ChP_inv, c_ChP_inv); } -// Calculation of E and B fields, using tetrads from gr-qc/0104063 -// Formalism from Alcubierre book template -EBFields_t -Weyl4::compute_EB_fields(const Vars &vars, - const Vars> &d1, - const Diff2Vars> &d2, - const Coordinates &coords) const +Tensor<3, data_t> +Weyl4::compute_epsilon3_LUU(const Vars &vars, + const Tensor<2, data_t> &h_UU) const { - EBFields_t out; - // raised normal vector, NB index 3 is time data_t n_U[4]; n_U[3] = 1. / vars.lapse; - FOR1(i) { n_U[i] = -vars.shift[i] / vars.lapse; } + FOR(i) { n_U[i] = -vars.shift[i] / vars.lapse; } // 4D levi civita symbol and 3D levi civita tensor in LLL and LUU form - const std::array, 4>, 4>, 4> - epsilon4 = TensorAlgebra::epsilon4D(); + const auto epsilon4 = TensorAlgebra::epsilon4D(); Tensor<3, data_t> epsilon3_LLL; Tensor<3, data_t> epsilon3_LUU; // Projection of antisymmentric Tensor onto hypersurface - see 8.3.17, // Alcubierre - FOR3(i, j, k) + FOR(i, j, k) { epsilon3_LLL[i][j][k] = 0.0; epsilon3_LUU[i][j][k] = 0.0; } // projection of 4-antisymetric tensor to 3-tensor on hypersurface // note last index contracted as per footnote 86 pg 290 Alcubierre - FOR3(i, j, k) + FOR(i, j, k) { for (int l = 0; l < 4; ++l) { epsilon3_LLL[i][j][k] += n_U[l] * epsilon4[i][j][k][l] * - vars.lapse * pow(vars.chi, -1.5); + vars.lapse / (vars.chi * sqrt(vars.chi)); } } // rasing indices - const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); - FOR3(i, j, k) + FOR(i, j, k) { - FOR2(m, n) + FOR(m, n) { epsilon3_LUU[i][j][k] += epsilon3_LLL[i][m][n] * h_UU[m][j] * vars.chi * h_UU[n][k] * vars.chi; } } + return epsilon3_LUU; +} + +// Calculation of E and B fields, using tetrads from gr-qc/0104063 +// BSSN expressions from Alcubierre book +// CCZ4 expressions calculated by MR and checked with TF see: +// https://www.overleaf.com/read/tvqjbyhvqqtp +template +EBFields_t Weyl4::compute_EB_fields( + const Vars &vars, const Vars> &d1, + const Diff2Vars> &d2, + const Tensor<3, data_t> &epsilon3_LUU, const Tensor<2, data_t> &h_UU, + const chris_t &chris) const +{ + EBFields_t out; + // Extrinsic curvature Tensor<2, data_t> K_tensor; Tensor<3, data_t> d1_K_tensor; Tensor<3, data_t> covariant_deriv_K_tensor; - // Compute inverse, Christoffel symbols and Ricci Tensor - using namespace TensorAlgebra; - const auto chris = compute_christoffel(d1.h, h_UU); - const auto ricci = CCZ4Geometry::compute_ricci(vars, d1, d2, h_UU, chris); + // Compute inverse, Christoffel symbols, Ricci tensor and Z terms + // Note that unlike in CCZ4 equations we want R_ij + 0.5(D_iZ_j + D_jZ_i) + // rather than R_ij + D_iZ_j + D_jZ_i hence use compute_ricci_Z_general + double dZ_coeff = (m_formulation == CCZ4RHS<>::USE_CCZ4) ? 1. : 0.; + auto ricci_and_Z_terms = CCZ4Geometry::compute_ricci_Z_general( + vars, d1, d2, h_UU, chris, dZ_coeff); // Compute full spatial Christoffel symbols + using namespace TensorAlgebra; const Tensor<3, data_t> chris_phys = compute_phys_chris(d1.chi, vars.chi, vars.h, h_UU, chris.ULL); // Extrinsic curvature and corresponding covariant and partial derivatives - FOR2(i, j) + FOR(i, j) { K_tensor[i][j] = vars.A[i][j] / vars.chi + 1. / 3. * (vars.h[i][j] * vars.K) / vars.chi; - FOR1(k) + FOR(k) { d1_K_tensor[i][j][k] = d1.A[i][j][k] / vars.chi - d1.chi[k] / vars.chi * K_tensor[i][j] + @@ -111,33 +140,96 @@ Weyl4::compute_EB_fields(const Vars &vars, } } // covariant derivative of K - FOR3(i, j, k) { covariant_deriv_K_tensor[i][j][k] = d1_K_tensor[i][j][k]; } - - FOR4(i, j, k, l) + FOR(i, j, k) { - covariant_deriv_K_tensor[i][j][k] += - -chris_phys[l][k][i] * K_tensor[l][j] - - chris_phys[l][k][j] * K_tensor[i][l]; + covariant_deriv_K_tensor[i][j][k] = d1_K_tensor[i][j][k]; + + FOR(l) + { + covariant_deriv_K_tensor[i][j][k] += + -chris_phys[l][k][i] * K_tensor[l][j] - + chris_phys[l][k][j] * K_tensor[i][l]; + } } + // Use 'K-Theta' in CCZ4. Just 'K' in BSSN. Not a mistake, this is not to + // confuse with the typical 'K-2*Theta' that appears in the CCZ4 equations + data_t K_minus_theta = vars.K; + if (m_formulation == CCZ4RHS<>::USE_CCZ4) + K_minus_theta -= vars.Theta; + // Calculate electric and magnetic fields - FOR2(i, j) + FOR(i, j) + { + out.E[i][j] = 0.0; + out.B[i][j] = 0.0; + } + + FOR(i, j) { - out.E[i][j] = 0; - out.B[i][j] = 0; + out.E[i][j] += + ricci_and_Z_terms.LL[i][j] + K_minus_theta * K_tensor[i][j]; + + FOR(k, l) + { + out.E[i][j] += + -K_tensor[i][k] * K_tensor[l][j] * h_UU[k][l] * vars.chi; + + out.B[i][j] += + epsilon3_LUU[i][k][l] * covariant_deriv_K_tensor[l][j][k]; + } } - FOR4(i, j, k, l) + // For CCZ4, explicit symmetrization appears naturally; + // For BSSN, only extra matter terms appear in the original expression, but + // assuming the Momentum constraints are satisfied, we can make the + // expression explicitly symmetric, which we enforce below + // (see Alcubierre chapter 8.3, from eq. 8.3.15 onwards) + TensorAlgebra::make_symmetric(out.B); + + // For CCZ4, Eij is explicitly trace-free; + // For BSSN, only extra matter terms appear in the original expression, but + // assuming the Hamiltonian constraint is satisfied, we can make the + // expression explicitly trace free, which we enforce below + // (see Alcubierre chapter 8.3, from eq. 8.3.15 onwards) + TensorAlgebra::make_trace_free(out.E, vars.h, h_UU); + + return out; +} + +// Calculation of the Weyl Scalars +template +WeylScalars_t Weyl4::compute_WeylCurvature(const EBFields_t &ebfields, + const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Tensor<2, data_t> &h_UU, + const Coordinates &coords) const +{ + WeylScalars_t out; + using namespace TensorAlgebra; + + + + // Projection of Electric and magnetic field components using tetrads + out.Weyl_curv = 0.0; + out.ChP_inv = 0.0; + + Tensor<2, data_t> chi_h_UU; + + FOR(i, j) { - out.B[i][j] += - epsilon3_LUU[i][k][l] * (covariant_deriv_K_tensor[l][j][k]); + chi_h_UU[i][j] = h_UU[i][j] * vars.chi; } - FOR2(i, j) { out.E[i][j] += ricci.LL[i][j] + vars.K * K_tensor[i][j]; } + // h_UU now is chi * h_UU ! + Tensor<2, data_t> E_UU = raise_all(ebfields.E, h_UU); + Tensor<2, data_t> B_UU = raise_all(ebfields.B, h_UU); - FOR4(i, j, k, l) + FOR(i, j) { - out.E[i][j] += -K_tensor[i][k] * K_tensor[l][j] * h_UU[k][l] * vars.chi; + out.Weyl_curv += 8.0 * (ebfields.E[i][j] * E_UU[i][j] - ebfields.B[i][j] * B_UU[i][j]); + out.ChP_inv += 16.0 * (ebfields.E[i][j] * B_UU[i][j]); } return out; @@ -149,17 +241,18 @@ NPScalar_t Weyl4::compute_Weyl4(const EBFields_t &ebfields, const Vars &vars, const Vars> &d1, const Diff2Vars> &d2, + const Tensor<2, data_t> &h_UU, const Coordinates &coords) const { NPScalar_t out; // Calculate the tetrads - const Tetrad_t tetrad = compute_null_tetrad(vars, coords); + const Tetrad_t tetrad = compute_null_tetrad(vars, h_UU, coords); // Projection of Electric and magnetic field components using tetrads out.Real = 0.0; out.Im = 0.0; - FOR2(i, j) + FOR(i, j) { out.Real += 0.5 * (ebfields.E[i][j] * (tetrad.w[i] * tetrad.w[j] - tetrad.v[i] * tetrad.v[j]) - @@ -179,6 +272,7 @@ NPScalar_t Weyl4::compute_Weyl4(const EBFields_t &ebfields, template Tetrad_t Weyl4::compute_null_tetrad(const Vars &vars, + const Tensor<2, data_t> &h_UU, const Coordinates &coords) const { Tetrad_t out; @@ -188,8 +282,7 @@ Weyl4::compute_null_tetrad(const Vars &vars, const double y = coords.y; const double z = coords.z; - // the inverse metric and alternating levi civita symbol - const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + // the alternating levi civita symbol const Tensor<3, double> epsilon = TensorAlgebra::epsilon(); // calculate the tetrad @@ -208,38 +301,38 @@ Weyl4::compute_null_tetrad(const Vars &vars, // floor on chi const data_t chi = simd_max(vars.chi, 1e-4); - FOR4(i, j, k, m) + FOR(i, j, k, m) { - out.w[i] += pow(chi, -0.5) * h_UU[i][j] * epsilon[j][k][m] * out.v[k] * + out.w[i] += 1. / sqrt(chi) * h_UU[i][j] * epsilon[j][k][m] * out.v[k] * out.u[m]; } // Gram Schmitt orthonormalisation // Choice of orthonormalisaion to avoid frame-dragging data_t omega_11 = 0.0; - FOR2(i, j) { omega_11 += out.v[i] * out.v[j] * vars.h[i][j] / chi; } - FOR1(i) { out.v[i] = out.v[i] / sqrt(omega_11); } + FOR(i, j) { omega_11 += out.v[i] * out.v[j] * vars.h[i][j] / chi; } + FOR(i) { out.v[i] = out.v[i] / sqrt(omega_11); } data_t omega_12 = 0.0; - FOR2(i, j) { omega_12 += out.v[i] * out.u[j] * vars.h[i][j] / chi; } - FOR1(i) { out.u[i] += -omega_12 * out.v[i]; } + FOR(i, j) { omega_12 += out.v[i] * out.u[j] * vars.h[i][j] / chi; } + FOR(i) { out.u[i] += -omega_12 * out.v[i]; } data_t omega_22 = 0.0; - FOR2(i, j) { omega_22 += out.u[i] * out.u[j] * vars.h[i][j] / chi; } - FOR1(i) { out.u[i] = out.u[i] / sqrt(omega_22); } + FOR(i, j) { omega_22 += out.u[i] * out.u[j] * vars.h[i][j] / chi; } + FOR(i) { out.u[i] = out.u[i] / sqrt(omega_22); } data_t omega_13 = 0.0; data_t omega_23 = 0.0; - FOR2(i, j) + FOR(i, j) { omega_13 += out.v[i] * out.w[j] * vars.h[i][j] / chi; omega_23 += out.u[i] * out.w[j] * vars.h[i][j] / chi; } - FOR1(i) { out.w[i] += -(omega_13 * out.v[i] + omega_23 * out.u[i]); } + FOR(i) { out.w[i] += -(omega_13 * out.v[i] + omega_23 * out.u[i]); } data_t omega_33 = 0.0; - FOR2(i, j) { omega_33 += out.w[i] * out.w[j] * vars.h[i][j] / chi; } - FOR1(i) { out.w[i] = out.w[i] / sqrt(omega_33); } + FOR(i, j) { omega_33 += out.w[i] * out.w[j] * vars.h[i][j] / chi; } + FOR(i) { out.w[i] = out.w[i] / sqrt(omega_33); } return out; } diff --git a/Source/CCZ4/Weyl4.impl.hpp_old b/Source/CCZ4/Weyl4.impl.hpp_old new file mode 100644 index 000000000..23a94ffb5 --- /dev/null +++ b/Source/CCZ4/Weyl4.impl.hpp_old @@ -0,0 +1,247 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(WEYL4_HPP_) +#error "This file should only be included through Weyl4.hpp" +#endif + +#ifndef WEYL4_IMPL_HPP_ +#define WEYL4_IMPL_HPP_ + +template void Weyl4::compute(Cell current_cell) const +{ + + // copy data from chombo gridpoint into local variables + const auto vars = current_cell.template load_vars(); + const auto d1 = m_deriv.template diff1(current_cell); + const auto d2 = m_deriv.template diff2(current_cell); + + // Get the coordinates + const Coordinates coords(current_cell, m_dx, m_center); + + // Compute the E and B fields + EBFields_t ebfields = compute_EB_fields(vars, d1, d2, coords); + + // work out the Newman Penrose scalar + NPScalar_t out = compute_Weyl4(ebfields, vars, d1, d2, coords); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(out.Real, c_Weyl4_Re); + current_cell.store_vars(out.Im, c_Weyl4_Im); +} + +// Calculation of E and B fields, using tetrads from gr-qc/0104063 +// Formalism from Alcubierre book +template +EBFields_t +Weyl4::compute_EB_fields(const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Coordinates &coords) const +{ + EBFields_t out; + + // raised normal vector, NB index 3 is time + data_t n_U[4]; + n_U[3] = 1. / vars.lapse; + FOR1(i) { n_U[i] = -vars.shift[i] / vars.lapse; } + + // 4D levi civita symbol and 3D levi civita tensor in LLL and LUU form + const std::array, 4>, 4>, 4> + epsilon4 = TensorAlgebra::epsilon4D(); + Tensor<3, data_t> epsilon3_LLL; + Tensor<3, data_t> epsilon3_LUU; + + // Projection of antisymmentric Tensor onto hypersurface - see 8.3.17, + // Alcubierre + FOR3(i, j, k) + { + epsilon3_LLL[i][j][k] = 0.0; + epsilon3_LUU[i][j][k] = 0.0; + } + // projection of 4-antisymetric tensor to 3-tensor on hypersurface + // note last index contracted as per footnote 86 pg 290 Alcubierre + FOR3(i, j, k) + { + for (int l = 0; l < 4; ++l) + { + epsilon3_LLL[i][j][k] += n_U[l] * epsilon4[i][j][k][l] * + vars.lapse * pow(vars.chi, -1.5); + } + } + // rasing indices + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + FOR3(i, j, k) + { + FOR2(m, n) + { + epsilon3_LUU[i][j][k] += epsilon3_LLL[i][m][n] * h_UU[m][j] * + vars.chi * h_UU[n][k] * vars.chi; + } + } + + // Extrinsic curvature + Tensor<2, data_t> K_tensor; + Tensor<3, data_t> d1_K_tensor; + Tensor<3, data_t> covariant_deriv_K_tensor; + + // Compute inverse, Christoffel symbols and Ricci Tensor + using namespace TensorAlgebra; + const auto chris = compute_christoffel(d1.h, h_UU); + const auto ricci = CCZ4Geometry::compute_ricci(vars, d1, d2, h_UU, chris); + + // Compute full spatial Christoffel symbols + const Tensor<3, data_t> chris_phys = + compute_phys_chris(d1.chi, vars.chi, vars.h, h_UU, chris.ULL); + + // Extrinsic curvature and corresponding covariant and partial derivatives + FOR2(i, j) + { + K_tensor[i][j] = vars.A[i][j] / vars.chi + + 1. / 3. * (vars.h[i][j] * vars.K) / vars.chi; + + FOR1(k) + { + d1_K_tensor[i][j][k] = d1.A[i][j][k] / vars.chi - + d1.chi[k] / vars.chi * K_tensor[i][j] + + 1. / 3. * d1.h[i][j][k] * vars.K / vars.chi + + 1. / 3. * vars.h[i][j] * d1.K[k] / vars.chi; + } + } + // covariant derivative of K + FOR3(i, j, k) { covariant_deriv_K_tensor[i][j][k] = d1_K_tensor[i][j][k]; } + + FOR4(i, j, k, l) + { + covariant_deriv_K_tensor[i][j][k] += + -chris_phys[l][k][i] * K_tensor[l][j] - + chris_phys[l][k][j] * K_tensor[i][l]; + } + + // Calculate electric and magnetic fields + FOR2(i, j) + { + out.E[i][j] = 0; + out.B[i][j] = 0; + } + + FOR4(i, j, k, l) + { + out.B[i][j] += + epsilon3_LUU[i][k][l] * (covariant_deriv_K_tensor[l][j][k]); + } + + FOR2(i, j) { out.E[i][j] += ricci.LL[i][j] + vars.K * K_tensor[i][j]; } + + FOR4(i, j, k, l) + { + out.E[i][j] += -K_tensor[i][k] * K_tensor[l][j] * h_UU[k][l] * vars.chi; + } + + return out; +} + +// Calculation of the Weyl4 scalar +template +NPScalar_t Weyl4::compute_Weyl4(const EBFields_t &ebfields, + const Vars &vars, + const Vars> &d1, + const Diff2Vars> &d2, + const Coordinates &coords) const +{ + NPScalar_t out; + + // Calculate the tetrads + const Tetrad_t tetrad = compute_null_tetrad(vars, coords); + + // Projection of Electric and magnetic field components using tetrads + out.Real = 0.0; + out.Im = 0.0; + FOR2(i, j) + { + out.Real += 0.5 * (ebfields.E[i][j] * (tetrad.w[i] * tetrad.w[j] - + tetrad.v[i] * tetrad.v[j]) - + 2.0 * ebfields.B[i][j] * tetrad.w[i] * tetrad.v[j]); + out.Im += 0.5 * (ebfields.B[i][j] * (-tetrad.w[i] * tetrad.w[j] + + tetrad.v[i] * tetrad.v[j]) - + 2.0 * ebfields.E[i][j] * tetrad.w[i] * tetrad.v[j]); + } + + return out; +} + +// Calculation of the null tetrad +// Defintions from gr-qc/0104063 +// "The Lazarus project: A pragmatic approach to binary black hole evolutions", +// Baker et al. +template +Tetrad_t +Weyl4::compute_null_tetrad(const Vars &vars, + const Coordinates &coords) const +{ + Tetrad_t out; + + // compute coords + const data_t x = coords.x; + const double y = coords.y; + const double z = coords.z; + + // the inverse metric and alternating levi civita symbol + const auto h_UU = TensorAlgebra::compute_inverse_sym(vars.h); + const Tensor<3, double> epsilon = TensorAlgebra::epsilon(); + + // calculate the tetrad + out.u[0] = x; + out.u[1] = y; + out.u[2] = z; + + out.v[0] = -y; + out.v[1] = x; + out.v[2] = 0.0; + + out.w[0] = 0.0; + out.w[1] = 0.0; + out.w[2] = 0.0; + + // floor on chi + const data_t chi = simd_max(vars.chi, 1e-4); + + FOR4(i, j, k, m) + { + out.w[i] += pow(chi, -0.5) * h_UU[i][j] * epsilon[j][k][m] * out.v[k] * + out.u[m]; + } + + // Gram Schmitt orthonormalisation + // Choice of orthonormalisaion to avoid frame-dragging + data_t omega_11 = 0.0; + FOR2(i, j) { omega_11 += out.v[i] * out.v[j] * vars.h[i][j] / chi; } + FOR1(i) { out.v[i] = out.v[i] / sqrt(omega_11); } + + data_t omega_12 = 0.0; + FOR2(i, j) { omega_12 += out.v[i] * out.u[j] * vars.h[i][j] / chi; } + FOR1(i) { out.u[i] += -omega_12 * out.v[i]; } + + data_t omega_22 = 0.0; + FOR2(i, j) { omega_22 += out.u[i] * out.u[j] * vars.h[i][j] / chi; } + FOR1(i) { out.u[i] = out.u[i] / sqrt(omega_22); } + + data_t omega_13 = 0.0; + data_t omega_23 = 0.0; + FOR2(i, j) + { + omega_13 += out.v[i] * out.w[j] * vars.h[i][j] / chi; + omega_23 += out.u[i] * out.w[j] * vars.h[i][j] / chi; + } + FOR1(i) { out.w[i] += -(omega_13 * out.v[i] + omega_23 * out.u[i]); } + + data_t omega_33 = 0.0; + FOR2(i, j) { omega_33 += out.w[i] * out.w[j] * vars.h[i][j] / chi; } + FOR1(i) { out.w[i] = out.w[i] / sqrt(omega_33); } + + return out; +} + +#endif /* WEYL4_HPP_ */ diff --git a/Source/Matter/MatterWeyl4.hpp b/Source/Matter/MatterWeyl4.hpp new file mode 100644 index 000000000..06df1de4f --- /dev/null +++ b/Source/Matter/MatterWeyl4.hpp @@ -0,0 +1,54 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef MATTERWEYL4_HPP_ +#define MATTERWEYL4_HPP_ + +#include "MatterCCZ4.hpp" +#include "Weyl4.hpp" + +//! Calculates the Weyl4 scalar for spacetimes with matter content +/*! + This class calculates the Weyl4 scalar real and im parts. It inherits from + the Weyl4 class and adds in the matter terms as appropriate depending on the + formulation +*/ +template class MatterWeyl4 : public Weyl4 +{ + public: + template + using Vars = typename MatterCCZ4::template Vars; + + //! Constructor + MatterWeyl4(matter_t a_matter, + const std::array a_center, + const double a_dx, + const int a_formulation = CCZ4RHS<>::USE_CCZ4, + double a_G_Newton = 1.0) + : Weyl4(a_center, a_dx, a_formulation), m_matter(a_matter), + m_G_Newton(a_G_Newton) + { + } + + //! The compute member which calculates the wave quantities at each point on + //! the grid + template void compute(Cell current_cell) const; + + protected: + matter_t m_matter; //!< The matter object, e.g. a scalar field + const double m_G_Newton; //!< Newton's constant, set to one by default + + //! Add matter terms to electric and magnetic parts + template + void add_matter_EB(EBFields_t &eb_fields, const Vars &vars, + const Vars> &d1, + const Tensor<3, data_t> &epsilon3_LUU, + const Tensor<2, data_t> &h_UU, + const chris_t &chris) const; +}; + +#include "MatterWeyl4.impl.hpp" + +#endif /* MATTERWEYL4_HPP_ */ diff --git a/Source/Matter/MatterWeyl4.impl.hpp b/Source/Matter/MatterWeyl4.impl.hpp new file mode 100644 index 000000000..d0bcb7b58 --- /dev/null +++ b/Source/Matter/MatterWeyl4.impl.hpp @@ -0,0 +1,78 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(MATTERWEYL4_HPP_) +#error "This file should only be included through MatterWeyl4.hpp" +#endif + +#ifndef MATTERWEYL4_IMPL_HPP_ +#define MATTERWEYL4_IMPL_HPP_ + +template +template +void MatterWeyl4::compute(Cell current_cell) const +{ + + // copy data from chombo gridpoint into local variables + const auto vars = current_cell.template load_vars(); + const auto d1 = m_deriv.template diff1(current_cell); + const auto d2 = m_deriv.template diff2(current_cell); + + // Get the coordinates + const Coordinates coords(current_cell, m_dx, m_center); + + // Compute the inverse metric + using namespace TensorAlgebra; + const auto h_UU = compute_inverse_sym(vars.h); + const auto chris = compute_christoffel(d1.h, h_UU); + + // Compute the spatial volume element + const auto epsilon3_LUU = compute_epsilon3_LUU(vars, h_UU); + + // Compute the E and B fields + EBFields_t ebfields = + compute_EB_fields(vars, d1, d2, epsilon3_LUU, h_UU, chris); + + // Add in matter terms to E and B fields + add_matter_EB(ebfields, vars, d1, epsilon3_LUU, h_UU, chris); + + // work out the Newman Penrose scalar + NPScalar_t out = + compute_Weyl4(ebfields, vars, d1, d2, h_UU, coords); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(out.Real, c_Weyl4_Re); + current_cell.store_vars(out.Im, c_Weyl4_Im); + + // work out the Newman Penrose scalar + WeylScalars_t out2 = + compute_WeylCurvature(ebfields, vars, d1, d2, h_UU, coords); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(out2.Weyl_curv, c_Weyl_curv); + current_cell.store_vars(out2.ChP_inv, c_ChP_inv); +} + +template +template +void MatterWeyl4::add_matter_EB(EBFields_t &ebfields, + const Vars &vars, + const Vars> &d1, + const Tensor<3, data_t> &epsilon3_LUU, + const Tensor<2, data_t> &h_UU, + const chris_t &chris) const +{ + // Calculate decomposed energy momentum tensor components + const auto emtensor = m_matter.compute_emtensor(vars, d1, h_UU, chris.ULL); + + Tensor<2, data_t> Sij_TF = emtensor.Sij; + TensorAlgebra::make_trace_free(Sij_TF, vars.h, h_UU); + + // as we made the vacuum expression of Bij explictly symmetric and Eij + // explictly trace-free, only Eij has matter terms + FOR(i, j) { ebfields.E[i][j] += -4.0 * M_PI * m_G_Newton * Sij_TF[i][j]; } +} + +#endif /* MATTERWEYL4_IMPL_HPP_ */ From 1240a2c47070dd93ac44ab626da7cc4ff513e393 Mon Sep 17 00:00:00 2001 From: cjoana Date: Wed, 5 Jul 2023 00:41:48 +0800 Subject: [PATCH 2/3] working evo and weyl --- Source/AMRInterpolator/AMRInterpolator.hpp | 7 +- .../AMRInterpolator/AMRInterpolator.impl.hpp | 71 +- .../AMRInterpolator/CylindricalGeometry.hpp | 38 +- Source/AMRInterpolator/Derivative.hpp | 18 +- Source/AMRInterpolator/DerivativeSetup.hpp | 8 +- Source/AMRInterpolator/InterpSource.hpp | 5 +- Source/AMRInterpolator/Lagrange.hpp | 28 +- Source/AMRInterpolator/Lagrange.impl.hpp | 79 +- Source/AMRInterpolator/MPIContext.impl.hpp | 6 +- Source/AMRInterpolator/QuinticConvolution.hpp | 23 +- .../QuinticConvolution.impl.hpp | 49 +- .../AMRInterpolator/SphericalExtraction.hpp | 29 +- Source/AMRInterpolator/SphericalGeometry.hpp | 85 +- Source/AMRInterpolator/SurfaceExtraction.hpp | 11 +- .../SurfaceExtraction.impl.hpp | 58 +- Source/BlackHoles/BHAMR.hpp | 17 - Source/BlackHoles/PunctureTracker.cpp | 16 +- Source/BlackHoles/PunctureTracker.hpp | 3 +- Source/BoxUtils/BoxLoops.impl.hpp | 11 +- Source/BoxUtils/BoxPointers.hpp | 12 +- Source/BoxUtils/Cell.impl.hpp | 13 +- Source/BoxUtils/ComputeModGrad.hpp | 4 +- Source/BoxUtils/FourthOrderDerivatives.hpp | 150 ++- Source/BoxUtils/NanCheck.hpp | 28 +- Source/BoxUtils/SixthOrderDerivatives.hpp | 155 ++- Source/CCZ4/CCZ4Geometry.hpp | 107 ++ Source/CCZ4/newGeometry.hpp | 193 +++ Source/GRChomboCore/BoundaryConditions.cpp | 63 +- Source/GRChomboCore/ChomboParameters.hpp | 149 ++- Source/GRChomboCore/GRAMR.cpp | 12 +- Source/GRChomboCore/GRAMR.hpp | 3 - Source/GRChomboCore/GRAMRLevel.cpp | 97 +- Source/GRChomboCore/GRAMRLevel.hpp | 18 +- Source/GRChomboCore/SetupFunctions.hpp | 21 +- .../GRChomboCore/SimulationParametersBase.hpp | 354 +++-- .../BlackHoles/BinaryBH.impl.hpp | 4 +- .../BlackHoles/BoostedBH.impl.hpp | 2 +- .../InitialConditions/BlackHoles/KerrBH.hpp | 4 +- .../BlackHoles/KerrBH.impl.hpp | 58 +- .../TwoPuncturesInitialData.impl.hpp | 4 +- .../ScalarFields/ScalarBubble.impl.hpp | 2 +- Source/Matter/MatterCCZ4.hpp | 107 +- .../ChiAndPhiTaggingCriterion.hpp | 2 +- .../ChiExtractionTaggingCriterion.hpp | 2 +- .../ChiPunctureExtractionTaggingCriterion.hpp | 2 +- .../TaggingCriteria/ChiTaggingCriterion.hpp | 4 +- .../PhiAndKTaggingCriterion.hpp | 6 +- Source/simd/arm/arm.hpp | 28 + Source/simd/arm/neon.hpp | 178 +++ Source/simd/arm/sve.hpp | 217 ++++ Source/simd/simd.hpp | 79 +- Source/simd/simd_base.hpp | 4 - Source/simd/x64/avx.hpp | 1 + Source/utils/CoordinateTransformations.hpp | 320 +++-- Source/utils/Coordinates.hpp | 25 +- Source/utils/DimensionDefinitions.hpp | 11 + Source/utils/FilesystemTools.hpp | 90 ++ Source/utils/SmallDataIO.cpp | 130 +- Source/utils/SmallDataIO.hpp | 8 +- Source/utils/Tensor.hpp | 2 +- Source/utils/TensorAlgebra.hpp | 196 +-- Source/utils/VarsTools.hpp | 41 +- Source/utils/WeylExtraction.hpp | 11 +- oSource/AMRInterpolator/AMRInterpolator.hpp | 134 ++ .../AMRInterpolator/AMRInterpolator.impl.hpp | 868 +++++++++++++ .../AMRInterpolator/CylindricalExtraction.hpp | 69 + .../AMRInterpolator/CylindricalGeometry.hpp | 101 ++ oSource/AMRInterpolator/Derivative.hpp | 167 +++ oSource/AMRInterpolator/DerivativeSetup.hpp | 27 + oSource/AMRInterpolator/IntegrationMethod.hpp | 79 ++ .../IntegrationMethodSetup.hpp | 20 + oSource/AMRInterpolator/InterpSource.hpp | 28 + .../InterpolationAlgorithm.hpp | 35 + .../AMRInterpolator/InterpolationLayout.hpp | 27 + .../AMRInterpolator/InterpolationQuery.hpp | 95 ++ oSource/AMRInterpolator/Lagrange.hpp | 76 ++ oSource/AMRInterpolator/Lagrange.impl.hpp | 367 ++++++ oSource/AMRInterpolator/MPIContext.hpp | 81 ++ oSource/AMRInterpolator/MPIContext.impl.hpp | 113 ++ oSource/AMRInterpolator/MPILayout.hpp | 44 + oSource/AMRInterpolator/MPILayout.impl.hpp | 71 ++ .../AMRInterpolator/QuinticConvolution.hpp | 39 + .../QuinticConvolution.impl.hpp | 196 +++ .../AMRInterpolator/SimpleArrayBox.hpp | 0 .../AMRInterpolator/SimpleInterpSource.hpp | 0 .../AMRInterpolator/SphericalExtraction.hpp | 151 +++ oSource/AMRInterpolator/SphericalGeometry.hpp | 119 ++ oSource/AMRInterpolator/SurfaceExtraction.hpp | 203 +++ .../SurfaceExtraction.impl.hpp | 514 ++++++++ .../ApparentHorizonFinder/AHData.hpp | 0 .../ApparentHorizonFinder/AHDeriv.hpp | 0 .../ApparentHorizonFinder/AHFinder.cpp | 0 .../ApparentHorizonFinder/AHFinder.hpp | 0 .../AHFunctionDefault.hpp | 0 .../ApparentHorizonFinder/AHFunctions.hpp | 0 .../ApparentHorizonFinder/AHGeometryData.hpp | 0 .../ApparentHorizonFinder/AHInterpolation.hpp | 0 .../AHInterpolation.impl.hpp | 0 .../AHSphericalGeometry.hpp | 0 .../AHStringGeometry.hpp | 0 .../ApparentHorizonFinder/ApparentHorizon.hpp | 0 .../ApparentHorizon.impl.hpp | 0 .../ApparentHorizon_petsc.impl.hpp | 0 oSource/BlackHoles/BHAMR.hpp | 48 + oSource/BlackHoles/PunctureTracker.cpp | 243 ++++ oSource/BlackHoles/PunctureTracker.hpp | 77 ++ oSource/BoxUtils/AMRReductions.hpp | 85 ++ oSource/BoxUtils/AMRReductions.impl.hpp | 143 +++ oSource/BoxUtils/BoxLoops.hpp | 90 ++ oSource/BoxUtils/BoxLoops.impl.hpp | 152 +++ oSource/BoxUtils/BoxPointers.hpp | 90 ++ oSource/BoxUtils/Cell.hpp | 147 +++ oSource/BoxUtils/Cell.impl.hpp | 94 ++ oSource/BoxUtils/CellIndex.hpp | 35 + oSource/BoxUtils/ComputeModGrad.hpp | 48 + oSource/BoxUtils/ComputePack.hpp | 78 ++ oSource/BoxUtils/FourthOrderDerivatives.hpp | 415 ++++++ oSource/BoxUtils/NanCheck.hpp | 73 ++ oSource/BoxUtils/SetValue.hpp | 51 + oSource/BoxUtils/SixthOrderDerivatives.hpp | 484 +++++++ oSource/CCZ4/ADMConformalVars.hpp | 143 +++ oSource/CCZ4/BSSNVars.hpp | 93 ++ oSource/CCZ4/CCZ4.hpp | 1 + oSource/CCZ4/CCZ4Geometry.hpp | 154 +++ oSource/CCZ4/CCZ4RHS.hpp | 130 ++ oSource/CCZ4/CCZ4RHS.impl.hpp | 229 ++++ oSource/CCZ4/CCZ4UserVariables.hpp | 76 ++ oSource/CCZ4/CCZ4Vars.hpp | 69 + oSource/CCZ4/Constraints.hpp | 66 + oSource/CCZ4/Constraints.impl.hpp | 80 ++ oSource/CCZ4/GammaCalculator.hpp | 59 + .../CCZ4/IntegratedMovingPunctureGauge.hpp | 96 ++ oSource/CCZ4/MovingPunctureGauge.hpp | 69 + oSource/CCZ4/MyGauge.hpp | 80 ++ oSource/CCZ4/MyNewConstraints.hpp | 99 ++ oSource/CCZ4/MyNewConstraints.impl.hpp | 202 +++ oSource/CCZ4/NewConstraints.hpp | 77 ++ oSource/CCZ4/NewConstraints.impl.hpp | 148 +++ oSource/CCZ4/PositiveChiAndAlpha.hpp | 41 + oSource/CCZ4/TraceARemoval.hpp | 51 + oSource/CCZ4/Weyl4.hpp | 131 ++ oSource/CCZ4/Weyl4.hpp_old | 102 ++ oSource/CCZ4/Weyl4.impl.hpp | 340 +++++ oSource/CCZ4/Weyl4.impl.hpp_old | 247 ++++ oSource/CCZ4/WeylCurvature.hpp | 108 ++ oSource/CCZ4/WeylCurvature.impl.hpp | 251 ++++ oSource/CCZ4/__CCZ4.hpp | 117 ++ oSource/CCZ4/__CCZ4.impl.hpp | 248 ++++ oSource/GRChomboCore/BoundaryConditions.cpp | 1135 +++++++++++++++++ oSource/GRChomboCore/BoundaryConditions.hpp | 225 ++++ oSource/GRChomboCore/ChomboParameters.hpp | 493 +++++++ oSource/GRChomboCore/DefaultLevelFactory.hpp | 40 + .../GRChomboCore/EmptyDiagnosticVariables.hpp | 23 + oSource/GRChomboCore/GRAMR.cpp | 57 + oSource/GRChomboCore/GRAMR.hpp | 79 ++ oSource/GRChomboCore/GRAMRLevel.cpp | 1101 ++++++++++++++++ oSource/GRChomboCore/GRAMRLevel.hpp | 243 ++++ oSource/GRChomboCore/GRLevelData.cpp | 127 ++ oSource/GRChomboCore/GRLevelData.hpp | 32 + oSource/GRChomboCore/SetupFunctions.hpp | 177 +++ .../GRChomboCore/SimulationParametersBase.hpp | 197 +++ oSource/GRChomboCore/UserVariables.inc.hpp | 145 +++ oSource/GRChomboCore/VariableType.hpp | 16 + .../__SimulationParametersBase.hpp | 0 .../InitialConditions/BlackHoles/BinaryBH.hpp | 52 + .../BlackHoles/BinaryBH.impl.hpp | 73 ++ .../BlackHoles/BoostedBH.hpp | 56 + .../BlackHoles/BoostedBH.impl.hpp | 96 ++ .../InitialConditions/BlackHoles/KerrBH.hpp | 71 ++ .../BlackHoles/KerrBH.impl.hpp | 173 +++ .../InitialConditions/BlackHoles/SingleBH.hpp | 0 .../BlackHoles/SingleBH.impl.hpp | 0 .../BlackHoles/TwoPuncturesInitialData.hpp | 55 + .../TwoPuncturesInitialData.impl.hpp | 93 ++ .../ScalarFields/ScalarBubble.hpp | 51 + .../ScalarFields/ScalarBubble.impl.hpp | 53 + oSource/Matter/BiScalarField.hpp | 156 +++ oSource/Matter/BiScalarField.impl.hpp | 264 ++++ oSource/Matter/ChiRelaxation.hpp | 83 ++ oSource/Matter/ChiRelaxation.impl.hpp | 74 ++ oSource/Matter/DefaultEOS.hpp | 33 + oSource/Matter/DefaultPotential.hpp | 47 + oSource/Matter/ExtendedMatterConstraints.hpp | 67 + .../Matter/ExtendedMatterConstraints.impl.hpp | 67 + oSource/Matter/FluidInitData.hpp | 71 ++ oSource/Matter/FluidInitData.impl.hpp | 69 + oSource/Matter/MatterCCZ4.hpp | 106 ++ oSource/Matter/MatterCCZ4.impl.hpp | 107 ++ oSource/Matter/MatterCCZ4RHS.hpp | 123 ++ oSource/Matter/MatterCCZ4RHS.impl.hpp | 108 ++ oSource/Matter/MatterConstraints.hpp | 68 + oSource/Matter/MatterConstraints.impl.hpp | 51 + oSource/Matter/MatterWeyl4.hpp | 54 + oSource/Matter/MatterWeyl4.impl.hpp | 78 ++ oSource/Matter/MyNewMatterConstraints.hpp | 75 ++ .../Matter/MyNewMatterConstraints.impl.hpp | 77 ++ oSource/Matter/NewMatterConstraints.hpp | 69 + oSource/Matter/NewMatterConstraints.impl.hpp | 65 + oSource/Matter/PerfectFluid.hpp | 203 +++ oSource/Matter/PerfectFluid.impl.bu.hpp | 534 ++++++++ oSource/Matter/PerfectFluid.impl.hpp | 707 ++++++++++ oSource/Matter/ScalarField.hpp | 131 ++ oSource/Matter/ScalarField.impl.hpp | 148 +++ .../ChiAndBiPhiTaggingCriterion.hpp | 0 .../ChiAndPhiTaggingCriterion.hpp | 74 ++ .../ChiExtractionTaggingCriterion.hpp | 94 ++ .../ChiPunctureExtractionTaggingCriterion.hpp | 130 ++ .../TaggingCriteria/ChiTaggingCriterion.hpp | 38 + .../FixedGridsTaggingCriterion.hpp | 46 + .../TaggingCriteria/KTaggingCriterion.hpp | 0 .../TaggingCriteria/KandWTaggingCriterion.hpp | 0 .../PhiAndKTaggingCriterion.hpp | 51 + oSource/simd/simd.hpp | 144 +++ oSource/simd/simd_base.hpp | 188 +++ oSource/simd/simdify.hpp | 199 +++ oSource/simd/x64/avx.hpp | 203 +++ oSource/simd/x64/avx512.hpp | 223 ++++ oSource/simd/x64/sse.hpp | 226 ++++ oSource/simd/x64/x64.hpp | 31 + oSource/utils/AlwaysInline.hpp | 15 + oSource/utils/ArrayTools.hpp | 71 ++ oSource/utils/Combinatorics.hpp | 38 + oSource/utils/CoordinateTransformations.hpp | 298 +++++ oSource/utils/Coordinates.hpp | 118 ++ oSource/utils/DebuggingTools.hpp | 84 ++ oSource/utils/DimensionDefinitions.hpp | 22 + oSource/utils/GRInterval.hpp | 30 + oSource/utils/GRParmParse.hpp | 159 +++ oSource/utils/MultiLevelTask.hpp | 93 ++ oSource/utils/SmallDataIO.cpp | 472 +++++++ oSource/utils/SmallDataIO.hpp | 163 +++ oSource/utils/SphericalHarmonics.hpp | 64 + oSource/utils/Tensor.hpp | 54 + oSource/utils/TensorAlgebra.hpp | 393 ++++++ oSource/utils/VarsTools.hpp | 98 ++ oSource/utils/WeylExtraction.hpp | 94 ++ 236 files changed, 23979 insertions(+), 1395 deletions(-) mode change 100755 => 100644 Source/AMRInterpolator/DerivativeSetup.hpp create mode 100644 Source/CCZ4/newGeometry.hpp mode change 100644 => 120000 Source/Matter/MatterCCZ4.hpp create mode 100644 Source/simd/arm/arm.hpp create mode 100644 Source/simd/arm/neon.hpp create mode 100644 Source/simd/arm/sve.hpp create mode 100644 Source/utils/FilesystemTools.hpp create mode 100644 oSource/AMRInterpolator/AMRInterpolator.hpp create mode 100644 oSource/AMRInterpolator/AMRInterpolator.impl.hpp create mode 100644 oSource/AMRInterpolator/CylindricalExtraction.hpp create mode 100644 oSource/AMRInterpolator/CylindricalGeometry.hpp create mode 100644 oSource/AMRInterpolator/Derivative.hpp create mode 100755 oSource/AMRInterpolator/DerivativeSetup.hpp create mode 100644 oSource/AMRInterpolator/IntegrationMethod.hpp create mode 100644 oSource/AMRInterpolator/IntegrationMethodSetup.hpp create mode 100644 oSource/AMRInterpolator/InterpSource.hpp create mode 100644 oSource/AMRInterpolator/InterpolationAlgorithm.hpp create mode 100644 oSource/AMRInterpolator/InterpolationLayout.hpp create mode 100644 oSource/AMRInterpolator/InterpolationQuery.hpp create mode 100644 oSource/AMRInterpolator/Lagrange.hpp create mode 100644 oSource/AMRInterpolator/Lagrange.impl.hpp create mode 100644 oSource/AMRInterpolator/MPIContext.hpp create mode 100644 oSource/AMRInterpolator/MPIContext.impl.hpp create mode 100644 oSource/AMRInterpolator/MPILayout.hpp create mode 100644 oSource/AMRInterpolator/MPILayout.impl.hpp create mode 100644 oSource/AMRInterpolator/QuinticConvolution.hpp create mode 100644 oSource/AMRInterpolator/QuinticConvolution.impl.hpp rename {Source => oSource}/AMRInterpolator/SimpleArrayBox.hpp (100%) rename {Source => oSource}/AMRInterpolator/SimpleInterpSource.hpp (100%) create mode 100644 oSource/AMRInterpolator/SphericalExtraction.hpp create mode 100644 oSource/AMRInterpolator/SphericalGeometry.hpp create mode 100644 oSource/AMRInterpolator/SurfaceExtraction.hpp create mode 100644 oSource/AMRInterpolator/SurfaceExtraction.impl.hpp rename {Source => oSource}/ApparentHorizonFinder/AHData.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHDeriv.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHFinder.cpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHFinder.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHFunctionDefault.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHFunctions.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHGeometryData.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHInterpolation.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHInterpolation.impl.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHSphericalGeometry.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/AHStringGeometry.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/ApparentHorizon.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/ApparentHorizon.impl.hpp (100%) rename {Source => oSource}/ApparentHorizonFinder/ApparentHorizon_petsc.impl.hpp (100%) create mode 100644 oSource/BlackHoles/BHAMR.hpp create mode 100644 oSource/BlackHoles/PunctureTracker.cpp create mode 100644 oSource/BlackHoles/PunctureTracker.hpp create mode 100644 oSource/BoxUtils/AMRReductions.hpp create mode 100644 oSource/BoxUtils/AMRReductions.impl.hpp create mode 100644 oSource/BoxUtils/BoxLoops.hpp create mode 100644 oSource/BoxUtils/BoxLoops.impl.hpp create mode 100644 oSource/BoxUtils/BoxPointers.hpp create mode 100644 oSource/BoxUtils/Cell.hpp create mode 100644 oSource/BoxUtils/Cell.impl.hpp create mode 100644 oSource/BoxUtils/CellIndex.hpp create mode 100644 oSource/BoxUtils/ComputeModGrad.hpp create mode 100644 oSource/BoxUtils/ComputePack.hpp create mode 100644 oSource/BoxUtils/FourthOrderDerivatives.hpp create mode 100644 oSource/BoxUtils/NanCheck.hpp create mode 100644 oSource/BoxUtils/SetValue.hpp create mode 100644 oSource/BoxUtils/SixthOrderDerivatives.hpp create mode 100644 oSource/CCZ4/ADMConformalVars.hpp create mode 100644 oSource/CCZ4/BSSNVars.hpp create mode 120000 oSource/CCZ4/CCZ4.hpp create mode 100644 oSource/CCZ4/CCZ4Geometry.hpp create mode 100644 oSource/CCZ4/CCZ4RHS.hpp create mode 100644 oSource/CCZ4/CCZ4RHS.impl.hpp create mode 100644 oSource/CCZ4/CCZ4UserVariables.hpp create mode 100644 oSource/CCZ4/CCZ4Vars.hpp create mode 100644 oSource/CCZ4/Constraints.hpp create mode 100644 oSource/CCZ4/Constraints.impl.hpp create mode 100644 oSource/CCZ4/GammaCalculator.hpp create mode 100644 oSource/CCZ4/IntegratedMovingPunctureGauge.hpp create mode 100644 oSource/CCZ4/MovingPunctureGauge.hpp create mode 100644 oSource/CCZ4/MyGauge.hpp create mode 100644 oSource/CCZ4/MyNewConstraints.hpp create mode 100644 oSource/CCZ4/MyNewConstraints.impl.hpp create mode 100644 oSource/CCZ4/NewConstraints.hpp create mode 100644 oSource/CCZ4/NewConstraints.impl.hpp create mode 100644 oSource/CCZ4/PositiveChiAndAlpha.hpp create mode 100644 oSource/CCZ4/TraceARemoval.hpp create mode 100644 oSource/CCZ4/Weyl4.hpp create mode 100644 oSource/CCZ4/Weyl4.hpp_old create mode 100644 oSource/CCZ4/Weyl4.impl.hpp create mode 100644 oSource/CCZ4/Weyl4.impl.hpp_old create mode 100644 oSource/CCZ4/WeylCurvature.hpp create mode 100644 oSource/CCZ4/WeylCurvature.impl.hpp create mode 100644 oSource/CCZ4/__CCZ4.hpp create mode 100644 oSource/CCZ4/__CCZ4.impl.hpp create mode 100644 oSource/GRChomboCore/BoundaryConditions.cpp create mode 100644 oSource/GRChomboCore/BoundaryConditions.hpp create mode 100644 oSource/GRChomboCore/ChomboParameters.hpp create mode 100644 oSource/GRChomboCore/DefaultLevelFactory.hpp create mode 100644 oSource/GRChomboCore/EmptyDiagnosticVariables.hpp create mode 100644 oSource/GRChomboCore/GRAMR.cpp create mode 100644 oSource/GRChomboCore/GRAMR.hpp create mode 100644 oSource/GRChomboCore/GRAMRLevel.cpp create mode 100644 oSource/GRChomboCore/GRAMRLevel.hpp create mode 100644 oSource/GRChomboCore/GRLevelData.cpp create mode 100644 oSource/GRChomboCore/GRLevelData.hpp create mode 100644 oSource/GRChomboCore/SetupFunctions.hpp create mode 100644 oSource/GRChomboCore/SimulationParametersBase.hpp create mode 100644 oSource/GRChomboCore/UserVariables.inc.hpp create mode 100644 oSource/GRChomboCore/VariableType.hpp rename {Source => oSource}/GRChomboCore/__SimulationParametersBase.hpp (100%) create mode 100644 oSource/InitialConditions/BlackHoles/BinaryBH.hpp create mode 100644 oSource/InitialConditions/BlackHoles/BinaryBH.impl.hpp create mode 100644 oSource/InitialConditions/BlackHoles/BoostedBH.hpp create mode 100644 oSource/InitialConditions/BlackHoles/BoostedBH.impl.hpp create mode 100644 oSource/InitialConditions/BlackHoles/KerrBH.hpp create mode 100644 oSource/InitialConditions/BlackHoles/KerrBH.impl.hpp rename {Source => oSource}/InitialConditions/BlackHoles/SingleBH.hpp (100%) rename {Source => oSource}/InitialConditions/BlackHoles/SingleBH.impl.hpp (100%) create mode 100644 oSource/InitialConditions/BlackHoles/TwoPuncturesInitialData.hpp create mode 100644 oSource/InitialConditions/BlackHoles/TwoPuncturesInitialData.impl.hpp create mode 100644 oSource/InitialConditions/ScalarFields/ScalarBubble.hpp create mode 100644 oSource/InitialConditions/ScalarFields/ScalarBubble.impl.hpp create mode 100644 oSource/Matter/BiScalarField.hpp create mode 100644 oSource/Matter/BiScalarField.impl.hpp create mode 100644 oSource/Matter/ChiRelaxation.hpp create mode 100644 oSource/Matter/ChiRelaxation.impl.hpp create mode 100644 oSource/Matter/DefaultEOS.hpp create mode 100644 oSource/Matter/DefaultPotential.hpp create mode 100644 oSource/Matter/ExtendedMatterConstraints.hpp create mode 100644 oSource/Matter/ExtendedMatterConstraints.impl.hpp create mode 100644 oSource/Matter/FluidInitData.hpp create mode 100644 oSource/Matter/FluidInitData.impl.hpp create mode 100644 oSource/Matter/MatterCCZ4.hpp create mode 100644 oSource/Matter/MatterCCZ4.impl.hpp create mode 100644 oSource/Matter/MatterCCZ4RHS.hpp create mode 100644 oSource/Matter/MatterCCZ4RHS.impl.hpp create mode 100644 oSource/Matter/MatterConstraints.hpp create mode 100644 oSource/Matter/MatterConstraints.impl.hpp create mode 100644 oSource/Matter/MatterWeyl4.hpp create mode 100644 oSource/Matter/MatterWeyl4.impl.hpp create mode 100644 oSource/Matter/MyNewMatterConstraints.hpp create mode 100644 oSource/Matter/MyNewMatterConstraints.impl.hpp create mode 100644 oSource/Matter/NewMatterConstraints.hpp create mode 100644 oSource/Matter/NewMatterConstraints.impl.hpp create mode 100644 oSource/Matter/PerfectFluid.hpp create mode 100644 oSource/Matter/PerfectFluid.impl.bu.hpp create mode 100644 oSource/Matter/PerfectFluid.impl.hpp create mode 100644 oSource/Matter/ScalarField.hpp create mode 100644 oSource/Matter/ScalarField.impl.hpp rename {Source => oSource}/TaggingCriteria/ChiAndBiPhiTaggingCriterion.hpp (100%) create mode 100644 oSource/TaggingCriteria/ChiAndPhiTaggingCriterion.hpp create mode 100644 oSource/TaggingCriteria/ChiExtractionTaggingCriterion.hpp create mode 100644 oSource/TaggingCriteria/ChiPunctureExtractionTaggingCriterion.hpp create mode 100644 oSource/TaggingCriteria/ChiTaggingCriterion.hpp create mode 100644 oSource/TaggingCriteria/FixedGridsTaggingCriterion.hpp rename {Source => oSource}/TaggingCriteria/KTaggingCriterion.hpp (100%) rename {Source => oSource}/TaggingCriteria/KandWTaggingCriterion.hpp (100%) create mode 100644 oSource/TaggingCriteria/PhiAndKTaggingCriterion.hpp create mode 100644 oSource/simd/simd.hpp create mode 100644 oSource/simd/simd_base.hpp create mode 100644 oSource/simd/simdify.hpp create mode 100644 oSource/simd/x64/avx.hpp create mode 100644 oSource/simd/x64/avx512.hpp create mode 100644 oSource/simd/x64/sse.hpp create mode 100644 oSource/simd/x64/x64.hpp create mode 100644 oSource/utils/AlwaysInline.hpp create mode 100644 oSource/utils/ArrayTools.hpp create mode 100644 oSource/utils/Combinatorics.hpp create mode 100644 oSource/utils/CoordinateTransformations.hpp create mode 100644 oSource/utils/Coordinates.hpp create mode 100644 oSource/utils/DebuggingTools.hpp create mode 100644 oSource/utils/DimensionDefinitions.hpp create mode 100644 oSource/utils/GRInterval.hpp create mode 100644 oSource/utils/GRParmParse.hpp create mode 100644 oSource/utils/MultiLevelTask.hpp create mode 100644 oSource/utils/SmallDataIO.cpp create mode 100644 oSource/utils/SmallDataIO.hpp create mode 100644 oSource/utils/SphericalHarmonics.hpp create mode 100644 oSource/utils/Tensor.hpp create mode 100644 oSource/utils/TensorAlgebra.hpp create mode 100644 oSource/utils/VarsTools.hpp create mode 100644 oSource/utils/WeylExtraction.hpp 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/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/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