From e7337578ea9b1751bde6d9df8dd2ee263f391597 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Wed, 6 May 2020 20:06:34 +0100 Subject: [PATCH 1/8] Add ADM Mass and Angular Momentum Extraction Add Richardson Extrapolation to Infinity --- Examples/KerrBH/DiagnosticVariables.hpp | 7 +- Examples/KerrBH/KerrBHLevel.cpp | 32 ++++- Examples/KerrBH/KerrBHLevel.hpp | 3 + Examples/KerrBH/Main_KerrBH.cpp | 9 ++ Examples/KerrBH/SimulationParameters.hpp | 3 + Examples/KerrBH/params_cheap.txt | 7 +- Source/AMRInterpolator/SurfaceExtraction.hpp | 5 + .../SurfaceExtraction.impl.hpp | 101 ++++++++++++-- Source/CCZ4/ADMMass.hpp | 126 ++++++++++++++++++ Source/utils/ADMMassExtraction.hpp | 77 +++++++++++ 10 files changed, 349 insertions(+), 21 deletions(-) create mode 100644 Source/CCZ4/ADMMass.hpp create mode 100644 Source/utils/ADMMassExtraction.hpp diff --git a/Examples/KerrBH/DiagnosticVariables.hpp b/Examples/KerrBH/DiagnosticVariables.hpp index 88ff5b7dc..b8cea3394 100644 --- a/Examples/KerrBH/DiagnosticVariables.hpp +++ b/Examples/KerrBH/DiagnosticVariables.hpp @@ -15,6 +15,9 @@ enum c_Mom2, c_Mom3, + c_Madm, + c_Jadm, + NUM_DIAGNOSTIC_VARS }; @@ -23,7 +26,9 @@ namespace DiagnosticVariables static const std::array variable_names = { "Ham", - "Mom1", "Mom2", "Mom3"}; + "Mom1", "Mom2", "Mom3", + + "M_adm", "J_adm"}; } #endif /* DIAGNOSTICVARIABLES_HPP */ diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp index 647a80a7b..2e938f887 100644 --- a/Examples/KerrBH/KerrBHLevel.cpp +++ b/Examples/KerrBH/KerrBHLevel.cpp @@ -20,6 +20,9 @@ #include "GammaCalculator.hpp" #include "KerrBH.hpp" +#include "ADMMass.hpp" +#include "ADMMassExtraction.hpp" + void KerrBHLevel::specificAdvance() { // Enforce the trace free A_ij condition and positive chi and alpha @@ -38,9 +41,9 @@ void KerrBHLevel::initialData() if (m_verbosity) pout() << "KerrBHLevel::initialData " << m_level << endl; - // First set everything to zero then calculate initial data Get the Kerr - // solution in the variables, then calculate the \tilde\Gamma^i numerically - // as these are non zero and not calculated in the Kerr ICs + // First set everything to zero then calculate initial data + // Get the Kerr solution in the variables, then calculate the \tilde\Gamma^i + // numerically as these are non zero and not calculated in the Kerr ICs BoxLoops::loop( make_compute_pack(SetValue(0.), KerrBH(m_p.kerr_params, m_dx)), m_state_new, m_state_new, INCLUDE_GHOST_CELLS); @@ -99,3 +102,26 @@ void KerrBHLevel::computeTaggingCriterion(FArrayBox &tagging_criterion, { BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion); } + +void KerrBHLevel::specificPostTimeStep() +{ + CH_TIME("KerrBHLevel::specificPostTimeStep"); + if (m_p.activate_extraction == 1) + { + // Populate the Weyl Scalar values on the grid + fillAllGhosts(); + BoxLoops::loop(ADMMass(m_p.extraction_params.center, m_dx), m_state_new, + m_state_new, EXCLUDE_GHOST_CELLS); + + // Do the extraction on the min extraction level + if (m_level == m_p.extraction_params.min_extraction_level()) + { + CH_TIME("ADMExtraction"); + // Now refresh the interpolator and do the interpolation + m_gr_amr.m_interpolator->refresh(); + ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, m_time, + m_restart_time); + my_extraction.execute_query(m_gr_amr.m_interpolator); + } + } +} \ No newline at end of file diff --git a/Examples/KerrBH/KerrBHLevel.hpp b/Examples/KerrBH/KerrBHLevel.hpp index cf19cc958..5f8753bfa 100644 --- a/Examples/KerrBH/KerrBHLevel.hpp +++ b/Examples/KerrBH/KerrBHLevel.hpp @@ -39,6 +39,9 @@ class KerrBHLevel : public GRAMRLevel /// Things to do before tagging cells (i.e. filling ghosts) virtual void preTagCells() override; + // to do post each time step on every level + virtual void specificPostTimeStep() override; + virtual void computeTaggingCriterion(FArrayBox &tagging_criterion, const FArrayBox ¤t_state) override; diff --git a/Examples/KerrBH/Main_KerrBH.cpp b/Examples/KerrBH/Main_KerrBH.cpp index 46f8c28f0..0d7944705 100644 --- a/Examples/KerrBH/Main_KerrBH.cpp +++ b/Examples/KerrBH/Main_KerrBH.cpp @@ -40,6 +40,15 @@ int runGRChombo(int argc, char *argv[]) DefaultLevelFactory kerr_bh_level_fact(gr_amr, sim_params); setupAMRObject(gr_amr, kerr_bh_level_fact); + // Set up interpolator: + // call this after amr object setup so grids known + // and need it to stay in scope throughout run + // 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.set_interpolator(&interpolator); + using Clock = std::chrono::steady_clock; using Minutes = std::chrono::duration>; diff --git a/Examples/KerrBH/SimulationParameters.hpp b/Examples/KerrBH/SimulationParameters.hpp index 0192e79b7..f5c177343 100644 --- a/Examples/KerrBH/SimulationParameters.hpp +++ b/Examples/KerrBH/SimulationParameters.hpp @@ -29,6 +29,8 @@ class SimulationParameters : public SimulationParametersBase pp.load("kerr_mass", kerr_params.mass); pp.load("kerr_spin", kerr_params.spin); pp.load("kerr_center", kerr_params.center, center); + + pp.load("activate_extraction", activate_extraction, false); } void check_params() @@ -51,6 +53,7 @@ class SimulationParameters : public SimulationParametersBase } KerrBH::params_t kerr_params; + bool activate_extraction; }; #endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Examples/KerrBH/params_cheap.txt b/Examples/KerrBH/params_cheap.txt index 838137673..ed2cd0a22 100644 --- a/Examples/KerrBH/params_cheap.txt +++ b/Examples/KerrBH/params_cheap.txt @@ -20,7 +20,7 @@ max_spatial_derivative_order = 4 # can be 4 or 6 # Params for Kerr BH kerr_mass = 1.0 -kerr_spin = 0.5 +kerr_spin = 0.3 # This is now defaulted to center of the grid so it does not need to be updated # when changing L (unless one wants to put the BH somewhere off center) #kerr_center = 64 64 64 @@ -29,7 +29,7 @@ kerr_spin = 0.5 # Level data # Maximum number of times you can regrid above coarsest level -max_level = 4 #4 There are (max_level+1) grids, so min is zero +max_level = 3 #4 There are (max_level+1) grids, so min is zero # Frequency of regridding at each level and thresholds on the tagging # Need one for each level, ie max_level+1 items @@ -65,6 +65,7 @@ vars_parity = 0 0 4 6 0 5 0 #chi and hij 0 1 2 3 #Theta and Gamma 0 1 2 3 1 2 3 #lapse shift and B vars_parity_diagnostic = 0 1 2 3 #Ham and Mom + 0 0 #ADM M and J # if sommerfeld boundaries selected, must select # asymptotic values (in order given by UserVariables.hpp) @@ -106,7 +107,7 @@ covariantZ4 = 1 # 0: keep kappa1; 1 [default]: replace kappa1 -> kappa1/lapse sigma = 0.3 #extraction params -activate_extraction = 0 +activate_extraction = 1 num_extraction_radii = 3 extraction_radii = 30. 40. 50. extraction_levels = 1 1 0 diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index c46dbf636..d5fae3eb2 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -187,6 +187,11 @@ template class SurfaceExtraction void write_integral(const std::string &a_filename, const std::vector a_integrals, const std::string a_label = "") const; + +protected: + std::vector + richardson_extrapolation(const std::vector> &integrals, + int extrapolation_order) const; }; #include "SurfaceExtraction.impl.hpp" diff --git a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp index 16eb7e951..4d2796000 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -397,11 +397,21 @@ template void SurfaceExtraction::write_integrals( const std::string &a_filename, const std::vector> &a_integrals, - const std::vector &a_labels) const + const std::vector &a_labels, int extrapolation_order) const { if (procID() == 0) { + std::vector extrapolations; + if (extrapolation_order) + { + extrapolations = + richardson_extrapolation(a_integrals, extrapolation_order); + } + const int num_integrals_per_surface = a_integrals.size(); + const int total_num_surfaces = + m_params.num_surfaces + (extrapolation_order > 0 ? 1 : 0); + // if labels are provided there must be the same number of labels as // there are integrals if (!a_labels.empty()) @@ -425,13 +435,13 @@ void SurfaceExtraction::write_integrals( { // make header strings std::vector header1_strings(num_integrals_per_surface * - m_params.num_surfaces); + total_num_surfaces); std::vector header2_strings(num_integrals_per_surface * - m_params.num_surfaces); - for (int isurface = 0; isurface < m_params.num_surfaces; ++isurface) + 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; if (a_labels.empty()) @@ -441,6 +451,13 @@ void SurfaceExtraction::write_integrals( header2_strings[idx] = std::to_string(m_params.surface_param_values[isurface]); } + if (extrapolation_order) + { + int idx = m_params.num_surfaces * num_integrals_per_surface + + iintegral; + header1_strings[idx] = a_labels[iintegral]; + header2_strings[idx] = "Infinity"; + } } std::string pre_header2_string = m_geom.param_name() + " = "; @@ -452,15 +469,20 @@ void SurfaceExtraction::write_integrals( // make vector of data for writing std::vector data_for_writing(num_integrals_per_surface * - m_params.num_surfaces); - for (int isurface = 0; isurface < m_params.num_surfaces; ++isurface) - { - for (int iintegral = 0; iintegral < num_integrals_per_surface; - ++iintegral) + total_num_surfaces); + 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; data_for_writing[idx] = a_integrals[iintegral][isurface]; } + if (extrapolation_order) + { + int idx = + m_params.num_surfaces * num_integrals_per_surface + iintegral; + data_for_writing[idx] = extrapolations[iintegral]; + } } // write data @@ -473,16 +495,67 @@ void SurfaceExtraction::write_integrals( template void SurfaceExtraction::write_integral( const std::string &a_filename, const std::vector a_integrals, - const std::string a_label) const + const std::string a_label, int extrapolation_order) const { std::vector> integrals(1, a_integrals); if (!a_label.empty()) { std::vector labels(1, a_label); - write_integrals(a_filename, integrals, labels); + write_integrals(a_filename, integrals, labels, extrapolation_order); } else - write_integrals(a_filename, integrals); + write_integrals(a_filename, integrals, extrapolation_order); +} + +// function to apply richardson extrapolation of 1st or 2nd order to calculated +// integrals +template +std::vector +SurfaceExtraction::richardson_extrapolation( + const std::vector> &integrals, + int extrapolation_order) const +{ + CH_TIME("SurfaceExtraction::richardson_extrapolation"); + int num_comps = integrals.size(); + int num_surfaces = m_params.num_surfaces; + CH_assert(extrapolation_order <= 2 && extrapolation_order >= 0); + + std::vector extrapolations(num_comps); + + if (num_surfaces >= 3 && extrapolation_order == 2) + { + for (int icomp = 0; icomp < num_comps; ++icomp) + { + int i1 = num_surfaces - 1; + int i2 = num_surfaces - 2; + int i3 = num_surfaces - 3; + double r1 = m_params.surface_param_values[i1]; + double r2 = m_params.surface_param_values[i2]; + double r3 = m_params.surface_param_values[i3]; + double c2 = r2 / r1 * (1. - r3 / r1) / (1. - r3 / r2); + double c3 = -r3 / r1 * (1. - r2 / r1) / (r2 / r3 - 1.); + double comp_inf = + (integrals[icomp][i1] - c2 * integrals[icomp][i2] - + c3 * integrals[icomp][i3]) / + (1. - c2 - c3); + extrapolations[icomp] = comp_inf; + } + } + else if (num_surfaces >= 2 && extrapolation_order >= 1) + { + for (int icomp = 0; icomp < num_comps; ++icomp) + { + int i1 = num_surfaces - 1; + int i2 = num_surfaces - 2; + double r1 = m_params.surface_param_values[i1]; + double r2 = m_params.surface_param_values[i2]; + double comp_inf = + (integrals[icomp][i1] - r2 / r1 * integrals[icomp][i2]) / + (1. - r2 / r1); + extrapolations[icomp] = comp_inf; + } + } + return extrapolations; } #endif /* SURFACEEXTRACTION_IMPL_HPP_ */ diff --git a/Source/CCZ4/ADMMass.hpp b/Source/CCZ4/ADMMass.hpp new file mode 100644 index 000000000..3ce1b77f8 --- /dev/null +++ b/Source/CCZ4/ADMMass.hpp @@ -0,0 +1,126 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ADMMASS_HPP_ +#define ADMMASS_HPP_ + +#include "ADMConformalVars.hpp" +#include "CCZ4Geometry.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "FourthOrderDerivatives.hpp" +#include "GRInterval.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "VarsTools.hpp" +#include "simd.hpp" + +//! Calculates the ADM massa +class ADMMass +{ + // Use the variable definition in ADMConformalVars - only require the key + // vars + template using Vars = ADMConformalVars::VarsNoGauge; + + template + using Diff1Vars = ADMConformalVars::Diff2VarsNoGauge; + + public: + enum DIR + { + NONE, + X, + Y, + Z + }; + + ADMMass(const std::array &a_center, double a_dx, + DIR spin_direction = Z, double a_G_Newton = 1.0) + : m_deriv(a_dx), m_center(a_center), m_G_Newton(a_G_Newton), + dir(spin_direction) + { + } + + template void compute(Cell current_cell) const + { + // copy data from chombo gridpoint into local variables, and calc 1st + // derivs + const auto vars = current_cell.template load_vars(); + const auto d1 = m_deriv.template diff1(current_cell); + + using namespace TensorAlgebra; + const auto h_UU = compute_inverse_sym(vars.h); + + // 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.; + FOR2(i, j) { dS_norm += vars.h[i][j] / vars.chi * dS_U[i] * dS_U[j]; } + dS_norm = sqrt(dS_norm); + FOR1(i) { dS_U[i] /= dS_norm; } + + Tensor<1, data_t> dS_L; + FOR1(i) + { + // dS_L[i] = dS_U[i]; + dS_L[i] = 0.; + FOR1(j) { dS_L[i] += vars.h[i][j] / vars.chi * dS_U[j]; } + } + + data_t Madm = 0.0; + FOR4(i, j, k, l) + { + Madm += dS_L[i] / (16. * M_PI * m_G_Newton) * pow(vars.chi, -1.5) * + h_UU[j][k] * h_UU[i][l] * + (vars.chi * (d1.h[l][k][j] - d1.h[j][k][l]) - + (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l])); + } + + // assign values of ADMMass in output box + current_cell.store_vars(Madm, c_Madm); + + if (dir == NONE) + return; + + // spin about z axis + data_t Jadm = 0.0; + + // note this is the levi civita symbol, + // not tensor (eps_tensor = eps_symbol * chi^-1.5) + const Tensor<3, double> epsilon = TensorAlgebra::epsilon(); + + FOR3(i, j, k) + { + Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) * + epsilon[dir - 1][j][k] * x[j] * vars.K * + TensorAlgebra::delta(i, k); + + FOR2(l, m) + { + Jadm += dS_L[i] / (8. * M_PI * m_G_Newton) * + epsilon[dir - 1][j][k] * x[j] * h_UU[i][l] * + h_UU[k][m] * vars.chi * + (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); + } + } + current_cell.store_vars(Jadm, c_Jadm); + } + + protected: + const FourthOrderDerivatives + m_deriv; //!< An object for calculating derivatives of the variables + const std::array &m_center; + const double m_G_Newton; //!< Newton's constant + + DIR dir; +}; + +#endif /* ADMMASS_HPP_ */ diff --git a/Source/utils/ADMMassExtraction.hpp b/Source/utils/ADMMassExtraction.hpp new file mode 100644 index 000000000..d1586138c --- /dev/null +++ b/Source/utils/ADMMassExtraction.hpp @@ -0,0 +1,77 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ADMMASSEXTRACTION_HPP_ +#define ADMMASSEXTRACTION_HPP_ + +#include "SphericalExtraction.hpp" +//! The class allows extraction of the values of the ADM mass and angular +//! momentum over spherical shells at specified radii, and integration over +//! those shells +/*! + The class allows the user to extract data from the grid for the ADM mass and + angular momentum over spherical shells at specified radii. The values may + then be written to an output file, or integrated across the surfaces. +*/ +class ADMMassExtraction : public SphericalExtraction +{ + public: + //! The constructor + ADMMassExtraction(SphericalExtraction::params_t &a_params, double a_dt, + double a_time, bool a_first_step, + double a_restart_time = 0.0, bool a_extract_J = true) + : SphericalExtraction(a_params, a_dt, a_time, a_first_step, + a_restart_time), + extract_J(a_extract_J) + { + add_var(c_Madm); + if (extract_J) + add_var(c_Jadm); + } + + //! The old constructor which assumes it is called in specificPostTimeStep + //! so the first time step is when m_time == m_dt + ADMMassExtraction(SphericalExtraction::params_t a_params, double a_dt, + double a_time, double a_restart_time = 0.0, + bool a_extract_J = true) + : ADMMassExtraction(a_params, a_dt, a_time, (a_dt == a_time), + a_restart_time, a_extract_J) + { + } + + //! Execute the query + void execute_query(AMRInterpolator> *a_interpolator) + { + // extract the values of the Weyl scalars on the spheres + extract(a_interpolator); + + if (m_params.write_extraction) + { + write_extraction("Weyl4ExtractionOut_"); + } + + // now calculate and write the requested spherical harmonic modes + std::vector> out_integrals(1 + extract_J); + + add_var_integrand(0, out_integrals[0], IntegrationMethod::simpson); + if (extract_J) + add_var_integrand(1, out_integrals[1], IntegrationMethod::simpson); + + // do the integration over the surface + integrate(); + int extrapolation_order = 2; + std::vector labels(1 + extract_J); + labels[0] = "M_adm"; + if (extract_J) + labels[1] = "J_adm"; + write_integrals("IntegralADMmass", out_integrals, labels, + extrapolation_order); + } + + private: + bool extract_J; +}; + +#endif /* ADMMASSEXTRACTION_HPP_ */ From 6e165f7dc76238d9fe749fbf06a564f9565818ef Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Wed, 13 May 2020 17:30:15 +0100 Subject: [PATCH 2/8] Resolve Katy's comments for ADM Mass PR --- Examples/KerrBH/KerrBHLevel.cpp | 22 +++++++++++++++++----- Examples/KerrBH/params.txt | 3 ++- Source/CCZ4/ADMMass.hpp | 10 +++++----- Source/utils/ADMMassExtraction.hpp | 4 ++-- 4 files changed, 26 insertions(+), 13 deletions(-) diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp index 2e938f887..5984b16e9 100644 --- a/Examples/KerrBH/KerrBHLevel.cpp +++ b/Examples/KerrBH/KerrBHLevel.cpp @@ -6,7 +6,7 @@ #include "KerrBHLevel.hpp" #include "BoxLoops.hpp" #include "CCZ4RHS.hpp" -#include "ChiTaggingCriterion.hpp" +#include "ChiExtractionTaggingCriterion.hpp" #include "ComputePack.hpp" #include "KerrBHLevel.hpp" #include "NanCheck.hpp" @@ -100,23 +100,35 @@ void KerrBHLevel::preTagCells() void KerrBHLevel::computeTaggingCriterion(FArrayBox &tagging_criterion, const FArrayBox ¤t_state) { - BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion); + BoxLoops::loop(ChiExtractionTaggingCriterion(m_dx, m_level, m_p.max_level, + m_p.extraction_params, + m_p.activate_extraction), + current_state, tagging_criterion); } void KerrBHLevel::specificPostTimeStep() { CH_TIME("KerrBHLevel::specificPostTimeStep"); + // Do the extraction on the min extraction level if (m_p.activate_extraction == 1) { - // Populate the Weyl Scalar values on the grid + auto min_level = m_p.extraction_params.min_extraction_level(); + + // ensure ADM Mass is only calculated on full cycles of 'min_level' + double extraction_dt = m_dt * pow(2., m_level - min_level); + double finest_dt = m_dt * pow(2., m_level - m_p.max_level); + if (fabs(remainder(m_time, extraction_dt)) > finest_dt / 2.) + return; + + // Populate the ADM Mass and Spin values on the grid fillAllGhosts(); BoxLoops::loop(ADMMass(m_p.extraction_params.center, m_dx), m_state_new, m_state_new, EXCLUDE_GHOST_CELLS); - // Do the extraction on the min extraction level - if (m_level == m_p.extraction_params.min_extraction_level()) + if (m_level == min_level) { CH_TIME("ADMExtraction"); + // Now refresh the interpolator and do the interpolation m_gr_amr.m_interpolator->refresh(); ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, m_time, diff --git a/Examples/KerrBH/params.txt b/Examples/KerrBH/params.txt index 81f66a6ca..a0ce00037 100644 --- a/Examples/KerrBH/params.txt +++ b/Examples/KerrBH/params.txt @@ -65,6 +65,7 @@ vars_parity = 0 0 4 6 0 5 0 #chi and hij 0 1 2 3 #Theta and Gamma 0 1 2 3 1 2 3 #lapse shift and B vars_parity_diagnostic = 0 1 2 3 #Ham and Mom + 0 0 #ADM M and J # if sommerfeld boundaries selected, must select # non zero asymptotic values @@ -104,7 +105,7 @@ covariantZ4 = 1 # 0: keep kappa1; 1 [default]: replace kappa1 -> kappa1/lapse sigma = 0.3 #extraction params -activate_extraction = 0 +activate_extraction = 1 num_extraction_radii = 3 extraction_radii = 30. 40. 50. extraction_levels = 2 1 0 diff --git a/Source/CCZ4/ADMMass.hpp b/Source/CCZ4/ADMMass.hpp index 3ce1b77f8..9b3599720 100644 --- a/Source/CCZ4/ADMMass.hpp +++ b/Source/CCZ4/ADMMass.hpp @@ -40,7 +40,7 @@ class ADMMass ADMMass(const std::array &a_center, double a_dx, DIR spin_direction = Z, double a_G_Newton = 1.0) : m_deriv(a_dx), m_center(a_center), m_G_Newton(a_G_Newton), - dir(spin_direction) + m_dir(spin_direction) { } @@ -87,7 +87,7 @@ class ADMMass // assign values of ADMMass in output box current_cell.store_vars(Madm, c_Madm); - if (dir == NONE) + if (m_dir == NONE) return; // spin about z axis @@ -100,13 +100,13 @@ class ADMMass FOR3(i, j, k) { Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) * - epsilon[dir - 1][j][k] * x[j] * vars.K * + epsilon[m_dir - 1][j][k] * x[j] * vars.K * TensorAlgebra::delta(i, k); FOR2(l, m) { Jadm += dS_L[i] / (8. * M_PI * m_G_Newton) * - epsilon[dir - 1][j][k] * x[j] * h_UU[i][l] * + epsilon[m_dir - 1][j][k] * x[j] * h_UU[i][l] * h_UU[k][m] * vars.chi * (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); } @@ -120,7 +120,7 @@ class ADMMass const std::array &m_center; const double m_G_Newton; //!< Newton's constant - DIR dir; + DIR m_dir; }; #endif /* ADMMASS_HPP_ */ diff --git a/Source/utils/ADMMassExtraction.hpp b/Source/utils/ADMMassExtraction.hpp index d1586138c..a6171b56e 100644 --- a/Source/utils/ADMMassExtraction.hpp +++ b/Source/utils/ADMMassExtraction.hpp @@ -44,12 +44,12 @@ class ADMMassExtraction : public SphericalExtraction //! Execute the query void execute_query(AMRInterpolator> *a_interpolator) { - // extract the values of the Weyl scalars on the spheres + // extract the values of the ADM mass and spin on the spheres extract(a_interpolator); if (m_params.write_extraction) { - write_extraction("Weyl4ExtractionOut_"); + write_extraction("MadmExtractionOut_"); } // now calculate and write the requested spherical harmonic modes From 23e2bfb37dba812735c898972ec3fecb07050174 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Thu, 14 May 2020 17:41:16 +0100 Subject: [PATCH 3/8] Extraction improvements for PR In BinaryBH Example, calculate Weyl only on the minimum extraction level full cycles Allow the user to extract only the ADM mass without adding the spin as a variable --- Examples/KerrBH/KerrBHLevel.cpp | 1 - Source/CCZ4/ADMMass.hpp | 8 +++++++- Source/utils/ADMMassExtraction.hpp | 10 +++++++++- 3 files changed, 16 insertions(+), 3 deletions(-) diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp index 5984b16e9..813334120 100644 --- a/Examples/KerrBH/KerrBHLevel.cpp +++ b/Examples/KerrBH/KerrBHLevel.cpp @@ -128,7 +128,6 @@ void KerrBHLevel::specificPostTimeStep() if (m_level == min_level) { CH_TIME("ADMExtraction"); - // Now refresh the interpolator and do the interpolation m_gr_amr.m_interpolator->refresh(); ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, m_time, diff --git a/Source/CCZ4/ADMMass.hpp b/Source/CCZ4/ADMMass.hpp index 9b3599720..6f99e368d 100644 --- a/Source/CCZ4/ADMMass.hpp +++ b/Source/CCZ4/ADMMass.hpp @@ -90,6 +90,12 @@ class ADMMass if (m_dir == NONE) return; + // user should be able to run this just for the ADM mass + auto int_Jadm = ChomboParameters::variable_name_to_enum("Jadm"); + if (int_Jadm < 0) + MayDay::Error("Please include 'c_Jadm' in UserVariables with name " + "'Jadm'"); + // spin about z axis data_t Jadm = 0.0; @@ -111,7 +117,7 @@ class ADMMass (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); } } - current_cell.store_vars(Jadm, c_Jadm); + current_cell.store_vars(Jadm, int_Jadm); } protected: diff --git a/Source/utils/ADMMassExtraction.hpp b/Source/utils/ADMMassExtraction.hpp index a6171b56e..4480af6af 100644 --- a/Source/utils/ADMMassExtraction.hpp +++ b/Source/utils/ADMMassExtraction.hpp @@ -28,7 +28,15 @@ class ADMMassExtraction : public SphericalExtraction { add_var(c_Madm); if (extract_J) - add_var(c_Jadm); + { + // user should be able to run this just for the ADM mass + auto int_Jadm = ChomboParameters::variable_name_to_enum("Jadm"); + if (int_Jadm < 0) + MayDay::Error( + "Please include 'c_Jadm' in UserVariables with name " + "'Jadm'"); + add_var(int_Jadm); + } } //! The old constructor which assumes it is called in specificPostTimeStep From 229908c40e99deea6520f10cbf080f959ca3e125 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Fri, 31 Jul 2020 15:23:48 +0100 Subject: [PATCH 4/8] Add arguments for user to specify which DiagnosticVariables to use for storage (more like the style of Miren's NewConstraints class) --- Examples/KerrBH/KerrBHLevel.cpp | 39 +++++++-------- Source/CCZ4/ADMMass.hpp | 79 +++++++++++++++--------------- Source/utils/ADMMassExtraction.hpp | 53 ++++++++++---------- 3 files changed, 86 insertions(+), 85 deletions(-) diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp index 813334120..6e665c356 100644 --- a/Examples/KerrBH/KerrBHLevel.cpp +++ b/Examples/KerrBH/KerrBHLevel.cpp @@ -112,27 +112,26 @@ void KerrBHLevel::specificPostTimeStep() // Do the extraction on the min extraction level if (m_p.activate_extraction == 1) { - auto min_level = m_p.extraction_params.min_extraction_level(); - - // ensure ADM Mass is only calculated on full cycles of 'min_level' - double extraction_dt = m_dt * pow(2., m_level - min_level); - double finest_dt = m_dt * pow(2., m_level - m_p.max_level); - if (fabs(remainder(m_time, extraction_dt)) > finest_dt / 2.) - return; - - // Populate the ADM Mass and Spin values on the grid - fillAllGhosts(); - BoxLoops::loop(ADMMass(m_p.extraction_params.center, m_dx), m_state_new, - m_state_new, EXCLUDE_GHOST_CELLS); - - if (m_level == min_level) + int min_level = m_p.extraction_params.min_extraction_level(); + bool calculate_adm = at_level_timestep_multiple(min_level); + if (calculate_adm) { - CH_TIME("ADMExtraction"); - // Now refresh the interpolator and do the interpolation - m_gr_amr.m_interpolator->refresh(); - ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, m_time, - m_restart_time); - my_extraction.execute_query(m_gr_amr.m_interpolator); + // Populate the ADM Mass and Spin values on the grid + fillAllGhosts(); + BoxLoops::loop( + ADMMass(m_p.extraction_params.center, m_dx, c_Madm, c_Jadm), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); + + if (m_level == min_level) + { + CH_TIME("ADMExtraction"); + // Now refresh the interpolator and do the interpolation + m_gr_amr.m_interpolator->refresh(); + ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, + m_time, m_restart_time, c_Madm, + c_Jadm); + my_extraction.execute_query(m_gr_amr.m_interpolator); + } } } } \ No newline at end of file diff --git a/Source/CCZ4/ADMMass.hpp b/Source/CCZ4/ADMMass.hpp index 6f99e368d..12ab6cb9f 100644 --- a/Source/CCZ4/ADMMass.hpp +++ b/Source/CCZ4/ADMMass.hpp @@ -31,19 +31,22 @@ class ADMMass public: enum DIR { - NONE, X, Y, Z }; ADMMass(const std::array &a_center, double a_dx, - DIR spin_direction = Z, double a_G_Newton = 1.0) - : m_deriv(a_dx), m_center(a_center), m_G_Newton(a_G_Newton), - m_dir(spin_direction) + 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) { } + // in case user wants to change direction of spin calculation to something + // other than Z + void set_spin_dir(DIR spin_direction) { m_dir = spin_direction; } + template void compute(Cell current_cell) const { // copy data from chombo gridpoint into local variables, and calc 1st @@ -75,49 +78,46 @@ class ADMMass FOR1(j) { dS_L[i] += vars.h[i][j] / vars.chi * dS_U[j]; } } - data_t Madm = 0.0; - FOR4(i, j, k, l) + if (m_c_Madm >= 0) { - Madm += dS_L[i] / (16. * M_PI * m_G_Newton) * pow(vars.chi, -1.5) * - h_UU[j][k] * h_UU[i][l] * - (vars.chi * (d1.h[l][k][j] - d1.h[j][k][l]) - - (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l])); - } - - // assign values of ADMMass in output box - current_cell.store_vars(Madm, c_Madm); - - if (m_dir == NONE) - return; - - // user should be able to run this just for the ADM mass - auto int_Jadm = ChomboParameters::variable_name_to_enum("Jadm"); - if (int_Jadm < 0) - MayDay::Error("Please include 'c_Jadm' in UserVariables with name " - "'Jadm'"); - - // spin about z axis - data_t Jadm = 0.0; + data_t Madm = 0.0; + FOR4(i, j, k, l) + { + Madm += dS_L[i] / (16. * M_PI * m_G_Newton) * + pow(vars.chi, -1.5) * h_UU[j][k] * h_UU[i][l] * + (vars.chi * (d1.h[l][k][j] - d1.h[j][k][l]) - + (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l])); + } - // note this is the levi civita symbol, - // not tensor (eps_tensor = eps_symbol * chi^-1.5) - const Tensor<3, double> epsilon = TensorAlgebra::epsilon(); + // assign values of ADMMass in output box + current_cell.store_vars(Madm, m_c_Madm); + } - FOR3(i, j, k) + if (m_c_Jadm >= 0) { - Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) * - epsilon[m_dir - 1][j][k] * x[j] * vars.K * - TensorAlgebra::delta(i, k); + // spin about m_dir axis (x, y or z) + data_t Jadm = 0.0; + + // note this is the levi civita symbol, + // not tensor (eps_tensor = eps_symbol * chi^-1.5) + const Tensor<3, double> epsilon = TensorAlgebra::epsilon(); - FOR2(l, m) + FOR3(i, j, k) { - Jadm += dS_L[i] / (8. * M_PI * m_G_Newton) * - epsilon[m_dir - 1][j][k] * x[j] * h_UU[i][l] * - h_UU[k][m] * vars.chi * - (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); + Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) * + epsilon[m_dir][j][k] * x[j] * vars.K * + TensorAlgebra::delta(i, k); + + FOR2(l, m) + { + Jadm += dS_L[i] / (8. * M_PI * m_G_Newton) * + epsilon[m_dir][j][k] * x[j] * h_UU[i][l] * + h_UU[k][m] * vars.chi * + (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); + } } + current_cell.store_vars(Jadm, m_c_Jadm); } - current_cell.store_vars(Jadm, int_Jadm); } protected: @@ -125,6 +125,7 @@ class ADMMass m_deriv; //!< An object for calculating derivatives of the variables const std::array &m_center; const double m_G_Newton; //!< Newton's constant + const int m_c_Madm, m_c_Jadm; DIR m_dir; }; diff --git a/Source/utils/ADMMassExtraction.hpp b/Source/utils/ADMMassExtraction.hpp index 4480af6af..7ed807242 100644 --- a/Source/utils/ADMMassExtraction.hpp +++ b/Source/utils/ADMMassExtraction.hpp @@ -21,31 +21,25 @@ class ADMMassExtraction : public SphericalExtraction //! The constructor ADMMassExtraction(SphericalExtraction::params_t &a_params, double a_dt, double a_time, bool a_first_step, - double a_restart_time = 0.0, bool a_extract_J = true) + double a_restart_time = 0.0, int a_c_Madm = -1, + int a_c_Jadm = -1) : SphericalExtraction(a_params, a_dt, a_time, a_first_step, a_restart_time), - extract_J(a_extract_J) + m_c_Madm(a_c_Madm), m_c_Jadm(a_c_Jadm) { - add_var(c_Madm); - if (extract_J) - { - // user should be able to run this just for the ADM mass - auto int_Jadm = ChomboParameters::variable_name_to_enum("Jadm"); - if (int_Jadm < 0) - MayDay::Error( - "Please include 'c_Jadm' in UserVariables with name " - "'Jadm'"); - add_var(int_Jadm); - } + if (m_c_Madm >= 0) + add_var(m_c_Madm, VariableType::diagnostic); + if (m_c_Jadm >= 0) + add_var(m_c_Jadm, VariableType::diagnostic); } //! The old constructor which assumes it is called in specificPostTimeStep //! so the first time step is when m_time == m_dt ADMMassExtraction(SphericalExtraction::params_t a_params, double a_dt, double a_time, double a_restart_time = 0.0, - bool a_extract_J = true) + int a_c_Madm = -1, int a_c_Jadm = -1) : ADMMassExtraction(a_params, a_dt, a_time, (a_dt == a_time), - a_restart_time, a_extract_J) + a_restart_time, a_c_Madm, a_c_Jadm) { } @@ -56,30 +50,37 @@ class ADMMassExtraction : public SphericalExtraction extract(a_interpolator); if (m_params.write_extraction) - { write_extraction("MadmExtractionOut_"); - } + + int num_integrals = (m_c_Madm >= 0) + (m_c_Jadm >= 0); + if (!num_integrals) + return; + + int J_index = m_c_Madm >= 0 ? 1 : 0; // now calculate and write the requested spherical harmonic modes - std::vector> out_integrals(1 + extract_J); + std::vector> out_integrals(num_integrals); - add_var_integrand(0, out_integrals[0], IntegrationMethod::simpson); - if (extract_J) - add_var_integrand(1, out_integrals[1], IntegrationMethod::simpson); + if (m_c_Madm >= 0) + add_var_integrand(0, out_integrals[0], IntegrationMethod::simpson); + if (m_c_Jadm >= 0) + add_var_integrand(J_index, out_integrals[J_index], + IntegrationMethod::simpson); // do the integration over the surface integrate(); int extrapolation_order = 2; - std::vector labels(1 + extract_J); - labels[0] = "M_adm"; - if (extract_J) - labels[1] = "J_adm"; + std::vector labels(num_integrals); + if (m_c_Madm >= 0) + labels[0] = "M_adm"; + if (m_c_Jadm >= 0) + labels[J_index] = "J_adm"; write_integrals("IntegralADMmass", out_integrals, labels, extrapolation_order); } private: - bool extract_J; + const int m_c_Madm, m_c_Jadm; }; #endif /* ADMMASSEXTRACTION_HPP_ */ From b77e7dda938b47ab3f01ba75d8f7607c917ff5e5 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Tue, 12 Jan 2021 19:38:52 +0100 Subject: [PATCH 5/8] Change ADMMass name to ADMQuantities --- Examples/KerrBH/KerrBHLevel.cpp | 17 ++++++----- .../CCZ4/{ADMMass.hpp => ADMQuantities.hpp} | 16 +++++----- ...action.hpp => ADMQuantitiesExtraction.hpp} | 30 +++++++++---------- 3 files changed, 33 insertions(+), 30 deletions(-) rename Source/CCZ4/{ADMMass.hpp => ADMQuantities.hpp} (91%) rename Source/utils/{ADMMassExtraction.hpp => ADMQuantitiesExtraction.hpp} (72%) diff --git a/Examples/KerrBH/KerrBHLevel.cpp b/Examples/KerrBH/KerrBHLevel.cpp index 6e665c356..d3381518a 100644 --- a/Examples/KerrBH/KerrBHLevel.cpp +++ b/Examples/KerrBH/KerrBHLevel.cpp @@ -20,8 +20,8 @@ #include "GammaCalculator.hpp" #include "KerrBH.hpp" -#include "ADMMass.hpp" -#include "ADMMassExtraction.hpp" +#include "ADMQuantities.hpp" +#include "ADMQuantitiesExtraction.hpp" void KerrBHLevel::specificAdvance() { @@ -118,18 +118,19 @@ void KerrBHLevel::specificPostTimeStep() { // Populate the ADM Mass and Spin values on the grid fillAllGhosts(); - BoxLoops::loop( - ADMMass(m_p.extraction_params.center, m_dx, c_Madm, c_Jadm), - m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); + BoxLoops::loop(ADMQuantities(m_p.extraction_params.center, m_dx, + c_Madm, c_Jadm), + m_state_new, m_state_diagnostics, + EXCLUDE_GHOST_CELLS); if (m_level == min_level) { CH_TIME("ADMExtraction"); // Now refresh the interpolator and do the interpolation m_gr_amr.m_interpolator->refresh(); - ADMMassExtraction my_extraction(m_p.extraction_params, m_dt, - m_time, m_restart_time, c_Madm, - c_Jadm); + ADMQuantitiesExtraction my_extraction( + m_p.extraction_params, m_dt, m_time, m_restart_time, c_Madm, + c_Jadm); my_extraction.execute_query(m_gr_amr.m_interpolator); } } diff --git a/Source/CCZ4/ADMMass.hpp b/Source/CCZ4/ADMQuantities.hpp similarity index 91% rename from Source/CCZ4/ADMMass.hpp rename to Source/CCZ4/ADMQuantities.hpp index 12ab6cb9f..81f876d45 100644 --- a/Source/CCZ4/ADMMass.hpp +++ b/Source/CCZ4/ADMQuantities.hpp @@ -3,8 +3,8 @@ * Please refer to LICENSE in GRChombo's root directory. */ -#ifndef ADMMASS_HPP_ -#define ADMMASS_HPP_ +#ifndef ADMQUANTITIES_HPP_ +#define ADMQUANTITIES_HPP_ #include "ADMConformalVars.hpp" #include "CCZ4Geometry.hpp" @@ -19,7 +19,7 @@ #include "simd.hpp" //! Calculates the ADM massa -class ADMMass +class ADMQuantities { // Use the variable definition in ADMConformalVars - only require the key // vars @@ -36,8 +36,8 @@ class ADMMass Z }; - ADMMass(const std::array &a_center, double a_dx, - int a_c_Madm = -1, int a_c_Jadm = -1, double a_G_Newton = 1.0) + 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) { @@ -89,7 +89,7 @@ class ADMMass (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l])); } - // assign values of ADMMass in output box + // assign values of ADM Mass in output box current_cell.store_vars(Madm, m_c_Madm); } @@ -116,6 +116,8 @@ class ADMMass (vars.A[l][m] + vars.K * vars.h[l][m] / 3.); } } + + // assign values of ADM Momentum in output box current_cell.store_vars(Jadm, m_c_Jadm); } } @@ -130,4 +132,4 @@ class ADMMass DIR m_dir; }; -#endif /* ADMMASS_HPP_ */ +#endif /* ADMQUANTITIES_HPP_ */ diff --git a/Source/utils/ADMMassExtraction.hpp b/Source/utils/ADMQuantitiesExtraction.hpp similarity index 72% rename from Source/utils/ADMMassExtraction.hpp rename to Source/utils/ADMQuantitiesExtraction.hpp index 7ed807242..c1cf8c48a 100644 --- a/Source/utils/ADMMassExtraction.hpp +++ b/Source/utils/ADMQuantitiesExtraction.hpp @@ -3,8 +3,8 @@ * Please refer to LICENSE in GRChombo's root directory. */ -#ifndef ADMMASSEXTRACTION_HPP_ -#define ADMMASSEXTRACTION_HPP_ +#ifndef ADMQUANTITIESEXTRACTION_HPP_ +#define ADMQUANTITIESEXTRACTION_HPP_ #include "SphericalExtraction.hpp" //! The class allows extraction of the values of the ADM mass and angular @@ -15,14 +15,14 @@ angular momentum over spherical shells at specified radii. The values may then be written to an output file, or integrated across the surfaces. */ -class ADMMassExtraction : public SphericalExtraction +class ADMQuantitiesExtraction : public SphericalExtraction { public: //! The constructor - ADMMassExtraction(SphericalExtraction::params_t &a_params, double a_dt, - double a_time, bool a_first_step, - double a_restart_time = 0.0, int a_c_Madm = -1, - int a_c_Jadm = -1) + ADMQuantitiesExtraction(SphericalExtraction::params_t &a_params, + double a_dt, double a_time, bool a_first_step, + double a_restart_time = 0.0, int a_c_Madm = -1, + int a_c_Jadm = -1) : SphericalExtraction(a_params, a_dt, a_time, a_first_step, a_restart_time), m_c_Madm(a_c_Madm), m_c_Jadm(a_c_Jadm) @@ -35,11 +35,11 @@ class ADMMassExtraction : public SphericalExtraction //! The old constructor which assumes it is called in specificPostTimeStep //! so the first time step is when m_time == m_dt - ADMMassExtraction(SphericalExtraction::params_t a_params, double a_dt, - double a_time, double a_restart_time = 0.0, - int a_c_Madm = -1, int a_c_Jadm = -1) - : ADMMassExtraction(a_params, a_dt, a_time, (a_dt == a_time), - a_restart_time, a_c_Madm, a_c_Jadm) + ADMQuantitiesExtraction(SphericalExtraction::params_t a_params, double a_dt, + double a_time, double a_restart_time = 0.0, + int a_c_Madm = -1, int a_c_Jadm = -1) + : ADMQuantitiesExtraction(a_params, a_dt, a_time, (a_dt == a_time), + a_restart_time, a_c_Madm, a_c_Jadm) { } @@ -50,7 +50,7 @@ class ADMMassExtraction : public SphericalExtraction extract(a_interpolator); if (m_params.write_extraction) - write_extraction("MadmExtractionOut_"); + write_extraction("ADMQuantitiesExtractionOut_"); int num_integrals = (m_c_Madm >= 0) + (m_c_Jadm >= 0); if (!num_integrals) @@ -75,7 +75,7 @@ class ADMMassExtraction : public SphericalExtraction labels[0] = "M_adm"; if (m_c_Jadm >= 0) labels[J_index] = "J_adm"; - write_integrals("IntegralADMmass", out_integrals, labels, + write_integrals("IntegralADMQuantities", out_integrals, labels, extrapolation_order); } @@ -83,4 +83,4 @@ class ADMMassExtraction : public SphericalExtraction const int m_c_Madm, m_c_Jadm; }; -#endif /* ADMMASSEXTRACTION_HPP_ */ +#endif /* ADMQUANTITIESEXTRACTION_HPP_ */ From 8d2f6821fb72600433d9a5a1d22fe389ecebad66 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Tue, 19 Jan 2021 16:49:16 +0100 Subject: [PATCH 6/8] Add extrapolation radii as parameter --- Examples/KerrBH/params_cheap.txt | 2 + Source/AMRInterpolator/SurfaceExtraction.hpp | 20 ++-- .../SurfaceExtraction.impl.hpp | 99 ++++++++++++------- .../GRChomboCore/SimulationParametersBase.hpp | 19 ++++ Source/utils/ADMQuantitiesExtraction.hpp | 4 +- Source/utils/WeylExtraction.hpp | 2 + 6 files changed, 100 insertions(+), 46 deletions(-) diff --git a/Examples/KerrBH/params_cheap.txt b/Examples/KerrBH/params_cheap.txt index ed2cd0a22..52a83e275 100644 --- a/Examples/KerrBH/params_cheap.txt +++ b/Examples/KerrBH/params_cheap.txt @@ -113,3 +113,5 @@ extraction_radii = 30. 40. 50. 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 diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index d5fae3eb2..1e5084261 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -50,11 +50,17 @@ template class SurfaceExtraction //!< extraction for each surface bool write_extraction; //!< whether or not to write the extracted data - int min_extraction_level() + std::vector + radii_idxs_for_extrapolation; // 2 for 2nd order, 3 for 3rd order + + int min_extraction_level() const { return *(std::min_element(extraction_levels.begin(), extraction_levels.end())); } + + //! deletes invalid or repeated + void validate_extrapolation_radii(); }; using vars_t = std::tuple; @@ -62,7 +68,7 @@ template class SurfaceExtraction protected: const SurfaceGeometry m_geom; //!< the geometry class which knows about //!< the particular surface - const params_t m_params; + params_t m_params; std::vector> m_vars; //!< the vector of pairs of //!< variables and derivatives to extract @@ -185,13 +191,13 @@ template class SurfaceExtraction //! convenience caller for write_integrals in the case of just integral per //! surface void write_integral(const std::string &a_filename, - const std::vector a_integrals, - const std::string a_label = "") const; + const std::vector &a_integrals, + const std::string &a_label = "") const; protected: - std::vector - richardson_extrapolation(const std::vector> &integrals, - int extrapolation_order) const; + // returns empty vector if no extrapolation is done + std::vector richardson_extrapolation( + const std::vector> &integrals) const; }; #include "SurfaceExtraction.impl.hpp" diff --git a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp index 4d2796000..6ba7ab4ab 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -50,6 +50,8 @@ SurfaceExtraction::SurfaceExtraction( } } } + + m_params.validate_extrapolation_radii(); } //! add a single variable or derivative of variable @@ -397,20 +399,16 @@ template void SurfaceExtraction::write_integrals( const std::string &a_filename, const std::vector> &a_integrals, - const std::vector &a_labels, int extrapolation_order) const + const std::vector &a_labels) const { if (procID() == 0) { - std::vector extrapolations; - if (extrapolation_order) - { - extrapolations = - richardson_extrapolation(a_integrals, extrapolation_order); - } + std::vector extrapolations = richardson_extrapolation(a_integrals); + bool extrapolation_done = extrapolations.size() > 0; const int num_integrals_per_surface = a_integrals.size(); const int total_num_surfaces = - m_params.num_surfaces + (extrapolation_order > 0 ? 1 : 0); + m_params.num_surfaces + (extrapolation_done > 0 ? 1 : 0); // if labels are provided there must be the same number of labels as // there are integrals @@ -451,7 +449,7 @@ void SurfaceExtraction::write_integrals( header2_strings[idx] = std::to_string(m_params.surface_param_values[isurface]); } - if (extrapolation_order) + if (extrapolation_done) { int idx = m_params.num_surfaces * num_integrals_per_surface + iintegral; @@ -477,7 +475,7 @@ void SurfaceExtraction::write_integrals( int idx = isurface * num_integrals_per_surface + iintegral; data_for_writing[idx] = a_integrals[iintegral][isurface]; } - if (extrapolation_order) + if (extrapolation_done) { int idx = m_params.num_surfaces * num_integrals_per_surface + iintegral; @@ -494,46 +492,74 @@ void SurfaceExtraction::write_integrals( //! surface template void SurfaceExtraction::write_integral( - const std::string &a_filename, const std::vector a_integrals, - const std::string a_label, int extrapolation_order) const + const std::string &a_filename, const std::vector &a_integrals, + const std::string &a_label) const +{ + write_integrals(a_filename, {a_integrals}, {a_label}); +} + +template +void SurfaceExtraction< + SurfaceGeometry>::params_t::validate_extrapolation_radii() { - std::vector> integrals(1, a_integrals); - if (!a_label.empty()) + std::vector valid_radii; + for (int i = 0; i < radii_idxs_for_extrapolation.size(); ++i) { - std::vector labels(1, a_label); - write_integrals(a_filename, integrals, labels, extrapolation_order); + // if valid + if (radii_idxs_for_extrapolation[i] < num_surfaces && + radii_idxs_for_extrapolation[i] > -num_surfaces) + { + // allow negative indices, such that '-1' is 'last', '-2' is 'one + // before last' + if (radii_idxs_for_extrapolation[i] < 0) + radii_idxs_for_extrapolation[i] += num_surfaces; + + // if not repeated already + if (i == 0 || std::find(radii_idxs_for_extrapolation.begin(), + radii_idxs_for_extrapolation.begin() + i, + radii_idxs_for_extrapolation[i]) == + radii_idxs_for_extrapolation.begin() + i) + { + valid_radii.push_back(radii_idxs_for_extrapolation[i]); + } + } } - else - write_integrals(a_filename, integrals, extrapolation_order); + radii_idxs_for_extrapolation = valid_radii; } -// function to apply richardson extrapolation of 1st or 2nd order to calculated +// function to apply richardson extrapolation of 2nd or 3rd order to calculated // integrals template std::vector SurfaceExtraction::richardson_extrapolation( - const std::vector> &integrals, - int extrapolation_order) const + const std::vector> &integrals) const { CH_TIME("SurfaceExtraction::richardson_extrapolation"); int num_comps = integrals.size(); int num_surfaces = m_params.num_surfaces; - CH_assert(extrapolation_order <= 2 && extrapolation_order >= 0); + + // already validated the radii + int extrapolation_order = m_params.radii_idxs_for_extrapolation.size(); + + if (extrapolation_order < 2) + { + return std::vector(); + } std::vector extrapolations(num_comps); - if (num_surfaces >= 3 && extrapolation_order == 2) + if (extrapolation_order == 3) { + int i1 = m_params.radii_idxs_for_extrapolation[0]; + int i2 = m_params.radii_idxs_for_extrapolation[1]; + int i3 = m_params.radii_idxs_for_extrapolation[2]; + double r1 = m_params.surface_param_values[i1]; + double r2 = m_params.surface_param_values[i2]; + double r3 = m_params.surface_param_values[i3]; + double c2 = r2 / r1 * (1. - r3 / r1) / (1. - r3 / r2); + double c3 = -r3 / r1 * (1. - r2 / r1) / (r2 / r3 - 1.); for (int icomp = 0; icomp < num_comps; ++icomp) { - int i1 = num_surfaces - 1; - int i2 = num_surfaces - 2; - int i3 = num_surfaces - 3; - double r1 = m_params.surface_param_values[i1]; - double r2 = m_params.surface_param_values[i2]; - double r3 = m_params.surface_param_values[i3]; - double c2 = r2 / r1 * (1. - r3 / r1) / (1. - r3 / r2); - double c3 = -r3 / r1 * (1. - r2 / r1) / (r2 / r3 - 1.); double comp_inf = (integrals[icomp][i1] - c2 * integrals[icomp][i2] - c3 * integrals[icomp][i3]) / @@ -541,20 +567,21 @@ SurfaceExtraction::richardson_extrapolation( extrapolations[icomp] = comp_inf; } } - else if (num_surfaces >= 2 && extrapolation_order >= 1) + else if (extrapolation_order == 2) { + int i1 = m_params.radii_idxs_for_extrapolation[0]; + int i2 = m_params.radii_idxs_for_extrapolation[1]; + double r1 = m_params.surface_param_values[i1]; + double r2 = m_params.surface_param_values[i2]; for (int icomp = 0; icomp < num_comps; ++icomp) { - int i1 = num_surfaces - 1; - int i2 = num_surfaces - 2; - double r1 = m_params.surface_param_values[i1]; - double r2 = m_params.surface_param_values[i2]; double comp_inf = (integrals[icomp][i1] - r2 / r1 * integrals[icomp][i2]) / (1. - r2 / r1); extrapolations[icomp] = comp_inf; } } + return extrapolations; } diff --git a/Source/GRChomboCore/SimulationParametersBase.hpp b/Source/GRChomboCore/SimulationParametersBase.hpp index 8be50660d..0a0cbdc8f 100644 --- a/Source/GRChomboCore/SimulationParametersBase.hpp +++ b/Source/GRChomboCore/SimulationParametersBase.hpp @@ -122,6 +122,16 @@ class SimulationParametersBase : public ChomboParameters } pp.load("write_extraction", extraction_params.write_extraction, false); + + if (pp.contains("extraction_extrapolation_order")) + { + int extrapolation_order; + pp.load("extraction_extrapolation_order", extrapolation_order); + pp.load("extraction_extrapolation_radii", + extraction_params.radii_idxs_for_extrapolation, + extrapolation_order); + extraction_params.validate_extrapolation_radii(); + } } void check_params() @@ -242,6 +252,15 @@ class SimulationParametersBase : public ChomboParameters check_parameter(mode_name, value_str, (l >= 2) && (abs(m) <= l), "l must be >= 2 and m must satisfy -l <= m <= l"); } + + int extrapolation_order = + extraction_params.radii_idxs_for_extrapolation.size(); + if (extrapolation_order > 1) + pout() << "Extraction with extrapolation order " + << extrapolation_order << std::endl; + else if (extrapolation_order == + 1) // 1 radius only because other were invalid + pout() << "Extraction radii not valid." << std::endl; } protected: diff --git a/Source/utils/ADMQuantitiesExtraction.hpp b/Source/utils/ADMQuantitiesExtraction.hpp index c1cf8c48a..c4459692c 100644 --- a/Source/utils/ADMQuantitiesExtraction.hpp +++ b/Source/utils/ADMQuantitiesExtraction.hpp @@ -69,14 +69,12 @@ class ADMQuantitiesExtraction : public SphericalExtraction // do the integration over the surface integrate(); - int extrapolation_order = 2; std::vector labels(num_integrals); if (m_c_Madm >= 0) labels[0] = "M_adm"; if (m_c_Jadm >= 0) labels[J_index] = "J_adm"; - write_integrals("IntegralADMQuantities", out_integrals, labels, - extrapolation_order); + write_integrals("IntegralADMQuantities", out_integrals, labels); } private: diff --git a/Source/utils/WeylExtraction.hpp b/Source/utils/WeylExtraction.hpp index 64b809e69..b86aecf18 100644 --- a/Source/utils/WeylExtraction.hpp +++ b/Source/utils/WeylExtraction.hpp @@ -87,6 +87,8 @@ class WeylExtraction : public SphericalExtraction std::move(mode_integrals[imode].second)}; std::vector labels = {"integral Re", "integral Im"}; write_integrals(integrals_filename, integrals_for_writing, labels); + write_integral(integrals_filename, integrals_for_writing[0], + labels[0]); } } }; From 613b4b9f0e5afa85c2ec01752fd6b7e7c41b0474 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Tue, 19 Jan 2021 18:27:01 +0100 Subject: [PATCH 7/8] Fixes after rebase --- Examples/KerrBH/Main_KerrBH.cpp | 3 ++- Examples/KerrBH/params_cheap.txt | 4 ++-- Source/AMRInterpolator/SurfaceExtraction.hpp | 2 +- .../SurfaceExtraction.impl.hpp | 23 +++++++++++-------- Source/CCZ4/ADMQuantities.hpp | 7 ++---- 5 files changed, 20 insertions(+), 19 deletions(-) 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.; From 03f845589ca16999cc6b65e981914bf3ff8f1067 Mon Sep 17 00:00:00 2001 From: TaigoFr Date: Fri, 26 Mar 2021 20:23:09 +0000 Subject: [PATCH 8/8] Add ADMQuantitiesTest (<1% error on M and J on level 0 after extrapolation) --- Examples/KerrBH/params.txt | 2 +- Examples/KerrBH/params_cheap.txt | 2 +- Source/AMRInterpolator/SurfaceExtraction.hpp | 1 - .../SurfaceExtraction.impl.hpp | 5 +- Source/CCZ4/ADMQuantities.hpp | 4 +- .../GRChomboCore/SimulationParametersBase.hpp | 35 +++--- Source/utils/ADMQuantitiesExtraction.hpp | 1 - Tests/ADMQuantitiesTest/ADMQuantitiesTest.cpp | 117 ++++++++++++++++++ .../ADMQuantitiesTest.inputs | 48 +++++++ .../ADMQuantitiesTestLevel.hpp | 57 +++++++++ .../ADMQuantitiesTest/DiagnosticVariables.hpp | 24 ++++ Tests/ADMQuantitiesTest/GNUmakefile | 23 ++++ .../SimulationParameters.hpp | 76 ++++++++++++ Tests/ADMQuantitiesTest/UserVariables.hpp | 29 +++++ .../InterpolatorTestLevel.hpp | 5 +- Tests/SphericalExtractionTest/GNUmakefile | 44 +++---- .../SphericalExtractionTestLevel.hpp | 7 +- 17 files changed, 429 insertions(+), 51 deletions(-) create mode 100644 Tests/ADMQuantitiesTest/ADMQuantitiesTest.cpp create mode 100644 Tests/ADMQuantitiesTest/ADMQuantitiesTest.inputs create mode 100644 Tests/ADMQuantitiesTest/ADMQuantitiesTestLevel.hpp create mode 100644 Tests/ADMQuantitiesTest/DiagnosticVariables.hpp create mode 100644 Tests/ADMQuantitiesTest/GNUmakefile create mode 100644 Tests/ADMQuantitiesTest/SimulationParameters.hpp create mode 100644 Tests/ADMQuantitiesTest/UserVariables.hpp diff --git a/Examples/KerrBH/params.txt b/Examples/KerrBH/params.txt index 37ba620d8..c97a17e45 100644 --- a/Examples/KerrBH/params.txt +++ b/Examples/KerrBH/params.txt @@ -153,7 +153,7 @@ num_extraction_radii = 3 extraction_radii = 30. 40. 50. extraction_levels = 2 1 0 num_points_phi = 50 -num_points_theta = 52 +num_points_theta = 51 # num_modes = 3 # modes = 2 0 # l m for spherical harmonics # 2 1 diff --git a/Examples/KerrBH/params_cheap.txt b/Examples/KerrBH/params_cheap.txt index 3c5bbe856..2a96c2923 100644 --- a/Examples/KerrBH/params_cheap.txt +++ b/Examples/KerrBH/params_cheap.txt @@ -153,7 +153,7 @@ num_extraction_radii = 3 extraction_radii = 30. 40. 50. extraction_levels = 1 1 0 num_points_phi = 50 -num_points_theta = 52 +num_points_theta = 51 # num_modes = 3 # modes = 2 0 # l m for spherical harmonics # 2 1 diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index 3f119e1b0..24b5a3295 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -198,7 +198,6 @@ template class SurfaceExtraction const std::vector &a_integrals, const std::string &a_label = "") const; - 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 3feb8fe48..72f2cb892 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -556,13 +556,14 @@ SurfaceExtraction::richardson_extrapolation( // already validated the radii int extrapolation_order = m_params.radii_idxs_for_extrapolation.size(); + int num_comps = integrals.size(); - if (extrapolation_order < 2) + if (extrapolation_order < 2 || num_comps == 0 || + integrals[0].size() < extrapolation_order) { 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 27b7a8327..b0f832aca 100644 --- a/Source/CCZ4/ADMQuantities.hpp +++ b/Source/CCZ4/ADMQuantities.hpp @@ -80,8 +80,8 @@ class ADMQuantities data_t Madm = 0.0; FOR4(i, j, k, l) { - Madm += dS_L[i] / (16. * M_PI * m_G_Newton) * - pow(vars.chi, -1.5) * h_UU[j][k] * h_UU[i][l] * + Madm += dS_L[i] / (16. * M_PI * m_G_Newton) / + (vars.chi * sqrt(vars.chi)) * h_UU[j][k] * h_UU[i][l] * (vars.chi * (d1.h[l][k][j] - d1.h[j][k][l]) - (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l])); } diff --git a/Source/GRChomboCore/SimulationParametersBase.hpp b/Source/GRChomboCore/SimulationParametersBase.hpp index 07e7e663f..3c0a75a02 100644 --- a/Source/GRChomboCore/SimulationParametersBase.hpp +++ b/Source/GRChomboCore/SimulationParametersBase.hpp @@ -158,15 +158,15 @@ class SimulationParametersBase : public ChomboParameters extraction_params.integral_file_prefix, std::string("Weyl4_mode_")); - if (pp.contains("extraction_extrapolation_order")) - { - int extrapolation_order; - pp.load("extraction_extrapolation_order", extrapolation_order); - pp.load("extraction_extrapolation_radii", - extraction_params.radii_idxs_for_extrapolation, - extrapolation_order); - extraction_params.validate_extrapolation_radii(); - } + if (pp.contains("extraction_extrapolation_order")) + { + int extrapolation_order; + pp.load("extraction_extrapolation_order", extrapolation_order); + pp.load("extraction_extrapolation_radii", + extraction_params.radii_idxs_for_extrapolation, + extrapolation_order); + extraction_params.validate_extrapolation_radii(); + } } } @@ -301,16 +301,15 @@ class SimulationParametersBase : public ChomboParameters mode_name, value_str, (l >= 2) && (abs(m) <= l), "l must be >= 2 and m must satisfy -l <= m <= l"); } - int extrapolation_order = - extraction_params.radii_idxs_for_extrapolation.size(); - if (extrapolation_order > 1) - pout() << "Extraction with extrapolation order " - << extrapolation_order << std::endl; - else if (extrapolation_order == - 1) // 1 radius only because other were invalid - pout() << "Extraction radii not valid." << std::endl; + int extrapolation_order = + extraction_params.radii_idxs_for_extrapolation.size(); + if (extrapolation_order > 1) + pout() << "Extraction with extrapolation order " + << extrapolation_order << std::endl; + else if (extrapolation_order == + 1) // 1 radius only because other were invalid + pout() << "Extraction radii not valid." << std::endl; } - } protected: diff --git a/Source/utils/ADMQuantitiesExtraction.hpp b/Source/utils/ADMQuantitiesExtraction.hpp index c4459692c..da06f5faf 100644 --- a/Source/utils/ADMQuantitiesExtraction.hpp +++ b/Source/utils/ADMQuantitiesExtraction.hpp @@ -58,7 +58,6 @@ class ADMQuantitiesExtraction : public SphericalExtraction int J_index = m_c_Madm >= 0 ? 1 : 0; - // now calculate and write the requested spherical harmonic modes std::vector> out_integrals(num_integrals); if (m_c_Madm >= 0) diff --git a/Tests/ADMQuantitiesTest/ADMQuantitiesTest.cpp b/Tests/ADMQuantitiesTest/ADMQuantitiesTest.cpp new file mode 100644 index 000000000..2a41cb11b --- /dev/null +++ b/Tests/ADMQuantitiesTest/ADMQuantitiesTest.cpp @@ -0,0 +1,117 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifdef CH_LANG_CC +/* + * _______ __ + * / ___/ / ___ __ _ / / ___ + * / /__/ _ \/ _ \/ V \/ _ \/ _ \ + * \___/_//_/\___/_/_/_/_.__/\___/ + * Please refer to LICENSE, in Chombo's root directory. + */ +#endif + +// Chombo includes +#include "parstream.H" //Gives us pout() + +#include "DefaultLevelFactory.hpp" +#include "GRAMR.hpp" + +#include "GRParmParse.hpp" +#include "SetupFunctions.hpp" +#include "SimulationParameters.hpp" + +// Problem specific includes: +#include "ADMQuantitiesExtraction.hpp" +#include "ADMQuantitiesTestLevel.hpp" +#include "AMRInterpolator.hpp" +#include "Lagrange.hpp" + +// Chombo namespace +#include "UsingNamespace.H" + +int runADMQuantitiesTest(int argc, char *argv[]) +{ + // Load the parameter file and construct the SimulationParameter class + // To add more parameters edit the SimulationParameters file. + std::string in_string = argv[argc - 1]; + char const *in_file = argv[argc - 1]; + GRParmParse pp(0, argv + argc, NULL, in_file); + SimulationParameters sim_params(pp); + + GRAMR gr_amr; + DefaultLevelFactory test_level_fact(gr_amr, + sim_params); + setupAMRObject(gr_amr, test_level_fact); + + AMRInterpolator> interpolator( + gr_amr, sim_params.origin, sim_params.dx, sim_params.boundary_params, + sim_params.verbosity); + + // one integral for ADM mass and one for momentum + std::vector> out_integrals(2); + + // perform extraction + ADMQuantitiesExtraction adm_extraction(sim_params.extraction_params, + sim_params.coarsest_dt, 0.0, true, + 0.0, c_Madm, c_Jadm); + adm_extraction.extract(&interpolator); + bool broadcast_integral = true; + adm_extraction.add_var_integrand( + 0, out_integrals[0], IntegrationMethod::simpson, + IntegrationMethod::trapezium, broadcast_integral); + adm_extraction.add_var_integrand( + 1, out_integrals[1], IntegrationMethod::simpson, + IntegrationMethod::trapezium, broadcast_integral); + adm_extraction.integrate(); + + // extrapolate results to infinity + std::vector extrapolations = + adm_extraction.richardson_extrapolation(out_integrals); + + CH_assert(extrapolations.size() == 2); + + int status = 0; + + // spin = J / M + double spin = extrapolations[1] / extrapolations[0]; + double mass_rel_error_percentage = + abs(extrapolations[0] / sim_params.kerr_params.mass - 1.) * 100.; + double spin_rel_error_percentage = + abs(spin / sim_params.kerr_params.spin - 1.) * 100.; + + pout() << std::setprecision(10); + pout() << "\nExtrapolation to infinity results:" << std::endl; + pout() << "ADM Mass = " << extrapolations[0] << std::endl; + pout() << "ADM Spin = " << spin << std::endl; + pout() << "Expected results:" << std::endl; + pout() << "ADM Mass = " << sim_params.kerr_params.mass << std::endl; + pout() << "ADM Spin = " << sim_params.kerr_params.spin << std::endl; + pout() << std::setprecision(3); + pout() << "Relative Error:" << std::endl; + pout() << "ADM Mass = " << mass_rel_error_percentage << "%" << std::endl; + pout() << "ADM Spin = " << spin_rel_error_percentage << "%" << std::endl; + + status |= (mass_rel_error_percentage > 1.); + status |= (spin_rel_error_percentage > 1.); + + return status; +} + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); + + int status = runADMQuantitiesTest(argc, argv); + + if (status == 0) + pout() << "ADMQuantitiesTest test passed." << endl; + else + pout() << "ADMQuantitiesTest test failed with return code " << status + << endl; + + mainFinalize(); + return status; +} diff --git a/Tests/ADMQuantitiesTest/ADMQuantitiesTest.inputs b/Tests/ADMQuantitiesTest/ADMQuantitiesTest.inputs new file mode 100644 index 000000000..c899df93e --- /dev/null +++ b/Tests/ADMQuantitiesTest/ADMQuantitiesTest.inputs @@ -0,0 +1,48 @@ + +################################################# +# Filesystem parameters + +chk_prefix = KerrBH_ +plot_prefix = KerrBHp_ +checkpoint_interval = 1 + +################################################# +# Initial Data parameters + +kerr_mass = 1. +kerr_spin = 0.5 + +################################################# +# Grid parameters + +N_full = 64 +L_full = 128 + +max_level = 0 +regrid_interval = 0 +regrid_threshold = 0.3 + +# Max and min box sizes +max_box_size = 16 +min_box_size = 16 + +################################################# +# Boundary Conditions parameters + +isPeriodic = 1 1 1 + +################################################# +# Extraction parameters + +num_extraction_radii = 3 +extraction_radii = 30. 40. 50. +extraction_levels = 0 0 0 +num_points_phi = 50 +num_points_theta = 51 + +extraction_extrapolation_order = 3 +extraction_extrapolation_radii = 0 1 2 + +################################################# + + diff --git a/Tests/ADMQuantitiesTest/ADMQuantitiesTestLevel.hpp b/Tests/ADMQuantitiesTest/ADMQuantitiesTestLevel.hpp new file mode 100644 index 000000000..42fe47939 --- /dev/null +++ b/Tests/ADMQuantitiesTest/ADMQuantitiesTestLevel.hpp @@ -0,0 +1,57 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ADMQUANTITIESLEVEL_HPP_ +#define ADMQUANTITIESLEVEL_HPP_ + +#include "ADMQuantities.hpp" +#include "BoxLoops.hpp" +#include "DiagnosticVariables.hpp" +#include "GRAMRLevel.hpp" +#include "SetValue.hpp" + +class ADMQuantitiesTestLevel : public GRAMRLevel +{ + friend class DefaultLevelFactory; + // Inherit the contructors from GRAMRLevel + using GRAMRLevel::GRAMRLevel; + + // initialize data + virtual void initialData() + { + // First set everything to zero then calculate initial data + // Get the Kerr solution in the variables, then calculate the + // \tilde\Gamma^i numerically as these are non zero and not calculated + // in the Kerr ICs + BoxLoops::loop( + make_compute_pack(SetValue(0.), KerrBH(m_p.kerr_params, m_dx)), + m_state_new, m_state_new, INCLUDE_GHOST_CELLS); + + fillAllGhosts(); + // we don't need the Gamma's for the ADM mass calculation + // BoxLoops::loop(GammaCalculator(m_dx), m_state_new, m_state_new, + // EXCLUDE_GHOST_CELLS); + + // set ghosts to 0 for extraction + BoxLoops::loop(SetValue(0.), m_state_diagnostics, m_state_diagnostics, + INCLUDE_GHOST_CELLS); + BoxLoops::loop( + ADMQuantities(m_p.extraction_params.center, m_dx, c_Madm, c_Jadm), + m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS); + } + + virtual void specificEvalRHS(GRLevelData &a_soln, GRLevelData &a_rhs, + const double a_time) + { + } + + virtual void computeTaggingCriterion(FArrayBox &tagging_criterion, + const FArrayBox ¤t_state) + { + tagging_criterion.setVal(0.); + }; +}; + +#endif /* ADMQUANTITIESLEVEL_HPP_ */ diff --git a/Tests/ADMQuantitiesTest/DiagnosticVariables.hpp b/Tests/ADMQuantitiesTest/DiagnosticVariables.hpp new file mode 100644 index 000000000..c03d07a7b --- /dev/null +++ b/Tests/ADMQuantitiesTest/DiagnosticVariables.hpp @@ -0,0 +1,24 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef DIAGNOSTICVARIABLES_HPP +#define DIAGNOSTICVARIABLES_HPP + +// assign an enum to each variable +enum +{ + c_Madm, + c_Jadm, + + NUM_DIAGNOSTIC_VARS +}; + +namespace DiagnosticVariables +{ +static const std::array variable_names = { + "M_adm", "J_adm"}; +} + +#endif /* DIAGNOSTICVARIABLES_HPP */ diff --git a/Tests/ADMQuantitiesTest/GNUmakefile b/Tests/ADMQuantitiesTest/GNUmakefile new file mode 100644 index 000000000..5506780ea --- /dev/null +++ b/Tests/ADMQuantitiesTest/GNUmakefile @@ -0,0 +1,23 @@ +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = ../../Source + +ebase := ADMQuantitiesTest + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/CCZ4 \ + $(GRCHOMBO_SOURCE)/BoxUtils \ + $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/AMRInterpolator \ + $(GRCHOMBO_SOURCE)/InitialConditions/BlackHoles + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Tests/ADMQuantitiesTest/SimulationParameters.hpp b/Tests/ADMQuantitiesTest/SimulationParameters.hpp new file mode 100644 index 000000000..14cfba215 --- /dev/null +++ b/Tests/ADMQuantitiesTest/SimulationParameters.hpp @@ -0,0 +1,76 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef SIMULATIONPARAMETERS_HPP_ +#define SIMULATIONPARAMETERS_HPP_ + +// General includes +#include "ChomboParameters.hpp" +#include "GRParmParse.hpp" + +// Problem specific includes +#include "KerrBH.hpp" +#include "SphericalExtraction.hpp" + +class SimulationParameters : public ChomboParameters +{ + public: + SimulationParameters(GRParmParse &pp) : ChomboParameters(pp) + { + read_params(pp); + } + + private: + void read_params(GRParmParse &pp) + { + // Initial Kerr data + pp.load("kerr_mass", kerr_params.mass); + pp.load("kerr_spin", kerr_params.spin); + pp.load("kerr_center", kerr_params.center, center); + + // Extraction params + pp.load("num_extraction_radii", extraction_params.num_extraction_radii, + 1); + + if (pp.contains("extraction_radii")) + { + pp.load("extraction_radii", extraction_params.extraction_radii, + extraction_params.num_extraction_radii); + } + else + { + pp.load("extraction_radius", extraction_params.extraction_radii, 1, + 0.1); + } + + pp.load("num_points_phi", extraction_params.num_points_phi, 2); + pp.load("num_points_theta", extraction_params.num_points_theta, 5); + if (extraction_params.num_points_theta % 2 == 0) + { + extraction_params.num_points_theta += 1; + pout() << "Parameter: num_points_theta incompatible with " + "Simpson's " + << "rule so increased by 1.\n"; + } + pp.load("extraction_center", extraction_params.center, center); + pp.load("write_extraction", extraction_params.write_extraction, false); + + if (pp.contains("extraction_extrapolation_order")) + { + int extrapolation_order; + pp.load("extraction_extrapolation_order", extrapolation_order); + pp.load("extraction_extrapolation_radii", + extraction_params.radii_idxs_for_extrapolation, + extrapolation_order); + extraction_params.validate_extrapolation_radii(); + } + } + + public: + KerrBH::params_t kerr_params; + SphericalExtraction::params_t extraction_params; +}; + +#endif /* SIMULATIONPARAMETERS_HPP_ */ diff --git a/Tests/ADMQuantitiesTest/UserVariables.hpp b/Tests/ADMQuantitiesTest/UserVariables.hpp new file mode 100644 index 000000000..b3ca4c817 --- /dev/null +++ b/Tests/ADMQuantitiesTest/UserVariables.hpp @@ -0,0 +1,29 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "ArrayTools.hpp" +#include "CCZ4UserVariables.hpp" +#include "DiagnosticVariables.hpp" + +/// This enum gives the index of every variable stored in the grid +enum +{ + // Note that it is important that the first enum value is set to 1 more than + // the last CCZ4 var enum + NUM_VARS = NUM_CCZ4_VARS, +}; + +namespace UserVariables +{ +static const std::array variable_names = + ccz4_variable_names; +} // namespace UserVariables + +#include "UserVariables.inc.hpp" + +#endif /* USERVARIABLES_HPP */ diff --git a/Tests/AMRInterpolatorTest/InterpolatorTestLevel.hpp b/Tests/AMRInterpolatorTest/InterpolatorTestLevel.hpp index b4f9fa5fd..a54bb0b81 100644 --- a/Tests/AMRInterpolatorTest/InterpolatorTestLevel.hpp +++ b/Tests/AMRInterpolatorTest/InterpolatorTestLevel.hpp @@ -30,7 +30,10 @@ class InterpolatorTestLevel : public GRAMRLevel } virtual void computeTaggingCriterion(FArrayBox &tagging_criterion, - const FArrayBox ¤t_state){}; + const FArrayBox ¤t_state) + { + tagging_criterion.setVal(0.); + }; }; #endif /* INTERPOLATORTESTLEVEL_HPP_ */ diff --git a/Tests/SphericalExtractionTest/GNUmakefile b/Tests/SphericalExtractionTest/GNUmakefile index 8b72a8165..f224a99d5 100644 --- a/Tests/SphericalExtractionTest/GNUmakefile +++ b/Tests/SphericalExtractionTest/GNUmakefile @@ -1,22 +1,22 @@ -# -*- Mode: Makefile -*- - -### This makefile produces an executable for each name in the `ebase' -### variable using the libraries named in the `LibNames' variable. - -# Included makefiles need an absolute path to the Chombo installation -# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) - -GRCHOMBO_SOURCE = $(shell pwd)/../../Source - -ebase := SphericalExtractionTest - -LibNames := AMRTimeDependent AMRTools BoxTools - -src_dirs := $(GRCHOMBO_SOURCE)/utils \ - $(GRCHOMBO_SOURCE)/simd \ - $(GRCHOMBO_SOURCE)/BoxUtils \ - $(GRCHOMBO_SOURCE)/CCZ4 \ - $(GRCHOMBO_SOURCE)/GRChomboCore \ - $(GRCHOMBO_SOURCE)/AMRInterpolator - -include $(CHOMBO_HOME)/mk/Make.test +# -*- Mode: Makefile -*- + +### This makefile produces an executable for each name in the `ebase' +### variable using the libraries named in the `LibNames' variable. + +# Included makefiles need an absolute path to the Chombo installation +# CHOMBO_HOME := Please set the CHOMBO_HOME locally (e.g. export CHOMBO_HOME=... in bash) + +GRCHOMBO_SOURCE = $(shell pwd)/../../Source + +ebase := SphericalExtractionTest + +LibNames := AMRTimeDependent AMRTools BoxTools + +src_dirs := $(GRCHOMBO_SOURCE)/utils \ + $(GRCHOMBO_SOURCE)/simd \ + $(GRCHOMBO_SOURCE)/BoxUtils \ + $(GRCHOMBO_SOURCE)/CCZ4 \ + $(GRCHOMBO_SOURCE)/GRChomboCore \ + $(GRCHOMBO_SOURCE)/AMRInterpolator + +include $(CHOMBO_HOME)/mk/Make.test diff --git a/Tests/SphericalExtractionTest/SphericalExtractionTestLevel.hpp b/Tests/SphericalExtractionTest/SphericalExtractionTestLevel.hpp index def66a1bc..e1a010127 100644 --- a/Tests/SphericalExtractionTest/SphericalExtractionTestLevel.hpp +++ b/Tests/SphericalExtractionTest/SphericalExtractionTestLevel.hpp @@ -32,7 +32,10 @@ class SphericalExtractionTestLevel : public GRAMRLevel } virtual void computeTaggingCriterion(FArrayBox &tagging_criterion, - const FArrayBox ¤t_state){}; + const FArrayBox ¤t_state) + { + tagging_criterion.setVal(0.); + }; }; -#endif /* SURFACEEXTRACTIONTESTLEVEL_HPP_ */ +#endif /* SPHERICALEXTRACTIONTESTLEVEL_HPP_ */