diff --git a/Examples/KerrBH/Main_KerrBH.cpp b/Examples/KerrBH/Main_KerrBH.cpp index 0d7944705..15496d241 100644 --- a/Examples/KerrBH/Main_KerrBH.cpp +++ b/Examples/KerrBH/Main_KerrBH.cpp @@ -46,7 +46,8 @@ int runGRChombo(int argc, char *argv[]) // Note: 'interpolator' needs to be in scope when gr_amr.run() is called, // otherwise pointer is lost AMRInterpolator> interpolator( - gr_amr, sim_params.origin, sim_params.dx, sim_params.verbosity); + gr_amr, sim_params.origin, sim_params.dx, sim_params.boundary_params, + sim_params.verbosity); gr_amr.set_interpolator(&interpolator); using Clock = std::chrono::steady_clock; diff --git a/Examples/KerrBH/params_cheap.txt b/Examples/KerrBH/params_cheap.txt index 52a83e275..d21801559 100644 --- a/Examples/KerrBH/params_cheap.txt +++ b/Examples/KerrBH/params_cheap.txt @@ -82,7 +82,7 @@ plot_interval = 0 num_plot_vars = 0 dt_multiplier = 0.25 stop_time = 10.0 -max_steps = 2 +max_steps = 1 #Lapse evolution lapse_power = 1.0 @@ -114,4 +114,4 @@ extraction_levels = 1 1 0 num_points_phi = 50 num_points_theta = 52 extraction_extrapolation_order = 3 -extraction_extrapolation_radii = 2 0 1 \ No newline at end of file +extraction_extrapolation_radii = 0 1 2 \ No newline at end of file diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index 1e5084261..6c37d7f49 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -194,7 +194,7 @@ template class SurfaceExtraction const std::vector &a_integrals, const std::string &a_label = "") const; -protected: + protected: // returns empty vector if no extrapolation is done std::vector richardson_extrapolation( const std::vector> &integrals) const; diff --git a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp index 6ba7ab4ab..e5249e83e 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -403,7 +403,8 @@ void SurfaceExtraction::write_integrals( { if (procID() == 0) { - std::vector extrapolations = richardson_extrapolation(a_integrals); + std::vector extrapolations = + richardson_extrapolation(a_integrals); bool extrapolation_done = extrapolations.size() > 0; const int num_integrals_per_surface = a_integrals.size(); @@ -439,7 +440,8 @@ void SurfaceExtraction::write_integrals( for (int iintegral = 0; iintegral < num_integrals_per_surface; ++iintegral) { - for (int isurface = 0; isurface < m_params.num_surfaces; ++isurface) + for (int isurface = 0; isurface < m_params.num_surfaces; + ++isurface) { int idx = isurface * num_integrals_per_surface + iintegral; if (a_labels.empty()) @@ -451,8 +453,9 @@ void SurfaceExtraction::write_integrals( } if (extrapolation_done) { - int idx = m_params.num_surfaces * num_integrals_per_surface + - iintegral; + int idx = + m_params.num_surfaces * num_integrals_per_surface + + iintegral; header1_strings[idx] = a_labels[iintegral]; header2_strings[idx] = "Infinity"; } @@ -468,8 +471,9 @@ void SurfaceExtraction::write_integrals( // make vector of data for writing std::vector data_for_writing(num_integrals_per_surface * total_num_surfaces); - for (int iintegral = 0; iintegral < num_integrals_per_surface; ++iintegral) - { + for (int iintegral = 0; iintegral < num_integrals_per_surface; + ++iintegral) + { for (int isurface = 0; isurface < m_params.num_surfaces; ++isurface) { int idx = isurface * num_integrals_per_surface + iintegral; @@ -477,8 +481,8 @@ void SurfaceExtraction::write_integrals( } if (extrapolation_done) { - int idx = - m_params.num_surfaces * num_integrals_per_surface + iintegral; + int idx = m_params.num_surfaces * num_integrals_per_surface + + iintegral; data_for_writing[idx] = extrapolations[iintegral]; } } @@ -535,8 +539,6 @@ SurfaceExtraction::richardson_extrapolation( const std::vector> &integrals) const { CH_TIME("SurfaceExtraction::richardson_extrapolation"); - int num_comps = integrals.size(); - int num_surfaces = m_params.num_surfaces; // already validated the radii int extrapolation_order = m_params.radii_idxs_for_extrapolation.size(); @@ -546,6 +548,7 @@ SurfaceExtraction::richardson_extrapolation( return std::vector(); } + int num_comps = integrals.size(); std::vector extrapolations(num_comps); if (extrapolation_order == 3) diff --git a/Source/CCZ4/ADMQuantities.hpp b/Source/CCZ4/ADMQuantities.hpp index 81f876d45..27b7a8327 100644 --- a/Source/CCZ4/ADMQuantities.hpp +++ b/Source/CCZ4/ADMQuantities.hpp @@ -38,8 +38,8 @@ class ADMQuantities ADMQuantities(const std::array &a_center, double a_dx, int a_c_Madm = -1, int a_c_Jadm = -1, double a_G_Newton = 1.0) - : m_deriv(a_dx), m_center(a_center), m_G_Newton(a_G_Newton), m_dir(Z), - m_c_Madm(a_c_Madm), m_c_Jadm(a_c_Jadm) + : m_deriv(a_dx), m_center(a_center), m_G_Newton(a_G_Newton), + m_c_Madm(a_c_Madm), m_c_Jadm(a_c_Jadm), m_dir(Z) { } @@ -59,10 +59,7 @@ class ADMQuantities // Surface element for integration Coordinates coords(current_cell, m_deriv.m_dx, m_center); - data_t r = coords.get_radius(); Tensor<1, data_t> x = {coords.x, coords.y, coords.z}; - // This is multiplied by r^2 as SphericalExtraction assumes it is - // normalised as such. Tensor<1, data_t> dS_U = x; data_t dS_norm = 0.;