diff --git a/Examples/BinaryBH/params.txt b/Examples/BinaryBH/params.txt index 5c8699efe..73aa93fda 100644 --- a/Examples/BinaryBH/params.txt +++ b/Examples/BinaryBH/params.txt @@ -180,6 +180,9 @@ modes = 2 0 # l m for spherical harmonics 4 3 4 4 +extraction_extrapolation_order = 3 +extraction_extrapolation_radii = 0 1 2 + # integral_file_prefix = "Weyl4_mode_" # write_extraction = 0 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..d3381518a 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" @@ -20,6 +20,9 @@ #include "GammaCalculator.hpp" #include "KerrBH.hpp" +#include "ADMQuantities.hpp" +#include "ADMQuantitiesExtraction.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); @@ -97,5 +100,39 @@ 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) + { + int min_level = m_p.extraction_params.min_extraction_level(); + bool calculate_adm = at_level_timestep_multiple(min_level); + if (calculate_adm) + { + // Populate the ADM Mass and Spin values on the grid + fillAllGhosts(); + 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(); + 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); + } + } + } +} \ 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..15496d241 100644 --- a/Examples/KerrBH/Main_KerrBH.cpp +++ b/Examples/KerrBH/Main_KerrBH.cpp @@ -40,6 +40,16 @@ 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.boundary_params, + 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 818942c09..1e3cd9b3d 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.txt b/Examples/KerrBH/params.txt index 1ce6a12d4..c97a17e45 100644 --- a/Examples/KerrBH/params.txt +++ b/Examples/KerrBH/params.txt @@ -95,6 +95,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 @@ -144,4 +145,26 @@ sigma = 0.3 # min_lapse = 1.e-4 ################################################# +# Extraction parameters + +# extraction_center = 256 256 256 # defaults to center +activate_extraction = 1 +num_extraction_radii = 3 +extraction_radii = 30. 40. 50. +extraction_levels = 2 1 0 +num_points_phi = 50 +num_points_theta = 51 +# num_modes = 3 +# modes = 2 0 # l m for spherical harmonics +# 2 1 +# 2 2 + +# integral_file_prefix = "Weyl4_mode_" + +# write_extraction = 0 +# extraction_subpath = "data/extraction" # directory for 'write_extraction = 1' +# extraction_file_prefix = "Weyl4_extraction_" + +################################################# + diff --git a/Examples/KerrBH/params_cheap.txt b/Examples/KerrBH/params_cheap.txt index 10135d106..2a96c2923 100644 --- a/Examples/KerrBH/params_cheap.txt +++ b/Examples/KerrBH/params_cheap.txt @@ -95,6 +95,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 @@ -118,7 +119,7 @@ max_steps = 4 # Spatial derivative order (only affects CCZ4 RHS) max_spatial_derivative_order = 4 # can be 4 or 6 -nan_check = 1 +# nan_check = 1 # Lapse evolution lapse_advec_coeff = 1.0 @@ -144,3 +145,29 @@ sigma = 0.3 # min_lapse = 1.e-4 ################################################# +# Extraction parameters + +# extraction_center = 256 256 256 # defaults to center +activate_extraction = 1 +num_extraction_radii = 3 +extraction_radii = 30. 40. 50. +extraction_levels = 1 1 0 +num_points_phi = 50 +num_points_theta = 51 +# num_modes = 3 +# modes = 2 0 # l m for spherical harmonics +# 2 1 +# 2 2 + +extraction_extrapolation_order = 3 +extraction_extrapolation_radii = 0 1 2 + +# integral_file_prefix = "Weyl4_mode_" + +# write_extraction = 0 +# extraction_subpath = "data/extraction" # directory for 'write_extraction = 1' +# extraction_file_prefix = "Weyl4_extraction_" + +################################################# + + diff --git a/Source/AMRInterpolator/SurfaceExtraction.hpp b/Source/AMRInterpolator/SurfaceExtraction.hpp index e060225a7..24b5a3295 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.hpp @@ -54,11 +54,17 @@ template class SurfaceExtraction std::string data_path, integral_file_prefix; std::string extraction_path, extraction_file_prefix; + std::vector + radii_idxs_for_extrapolation; // 2 for 2nd order, 3 for 3rd order + int min_extraction_level() { return *(std::min_element(extraction_levels.begin(), extraction_levels.end())); } + + //! deletes invalid or repeated + void validate_extrapolation_radii(); }; using vars_t = std::tuple; @@ -66,7 +72,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 @@ -189,8 +195,12 @@ 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; + + // 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 93d186c75..18305b612 100644 --- a/Source/AMRInterpolator/SurfaceExtraction.impl.hpp +++ b/Source/AMRInterpolator/SurfaceExtraction.impl.hpp @@ -62,6 +62,8 @@ SurfaceExtraction::SurfaceExtraction( } } } + + m_params.validate_extrapolation_radii(); } //! add a single variable or derivative of variable @@ -414,7 +416,14 @@ void SurfaceExtraction::write_integrals( { if (procID() == 0) { + 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_done > 0 ? 1 : 0); + // if labels are provided there must be the same number of labels as // there are integrals if (!a_labels.empty()) @@ -439,13 +448,14 @@ 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()) @@ -455,6 +465,14 @@ void SurfaceExtraction::write_integrals( header2_strings[idx] = std::to_string(m_params.surface_param_values[isurface]); } + if (extrapolation_done) + { + 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() + " = "; @@ -466,15 +484,21 @@ 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) + 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; data_for_writing[idx] = a_integrals[iintegral][isurface]; } + if (extrapolation_done) + { + int idx = m_params.num_surfaces * num_integrals_per_surface + + iintegral; + data_for_writing[idx] = extrapolations[iintegral]; + } } // write data @@ -486,17 +510,97 @@ 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) 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); + // 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]); + } + } + } + radii_idxs_for_extrapolation = valid_radii; +} + +// function to apply richardson extrapolation of 2nd or 3rd order to calculated +// integrals +template +std::vector +SurfaceExtraction::richardson_extrapolation( + const std::vector> &integrals) const +{ + CH_TIME("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 || num_comps == 0 || + integrals[0].size() < extrapolation_order) + { + return std::vector(); + } + + std::vector extrapolations(num_comps); + + 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) + { + double comp_inf = + (integrals[icomp][i1] - c2 * integrals[icomp][i2] - + c3 * integrals[icomp][i3]) / + (1. - c2 - c3); + extrapolations[icomp] = comp_inf; + } } - else - write_integrals(a_filename, integrals); + 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) + { + 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/ADMQuantities.hpp b/Source/CCZ4/ADMQuantities.hpp new file mode 100644 index 000000000..52dc72c9e --- /dev/null +++ b/Source/CCZ4/ADMQuantities.hpp @@ -0,0 +1,132 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ADMQUANTITIES_HPP_ +#define ADMQUANTITIES_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 ADMQuantities +{ + // Use the variable definition in ADMConformalVars - only require the key + // vars + template using Vars = ADMConformalVars::VarsNoGauge; + + template + using Diff1Vars = ADMConformalVars::Diff2VarsNoGauge; + + public: + enum DIR + { + X, + Y, + Z + }; + + 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_c_Madm(a_c_Madm), m_c_Jadm(a_c_Jadm), m_dir(Z) + { + } + + // 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 + // 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); + Tensor<1, data_t> x = {coords.x, coords.y, coords.z}; + Tensor<1, data_t> dS_U = x; + + data_t dS_norm = 0.; + FOR(i, j) { dS_norm += vars.h[i][j] / vars.chi * dS_U[i] * dS_U[j]; } + dS_norm = sqrt(dS_norm); + FOR(i) { dS_U[i] /= dS_norm; } + + Tensor<1, data_t> dS_L; + FOR(i) + { + // dS_L[i] = dS_U[i]; + dS_L[i] = 0.; + FOR(j) { dS_L[i] += vars.h[i][j] / vars.chi * dS_U[j]; } + } + + if (m_c_Madm >= 0) + { + data_t Madm = 0.0; + FOR(i, j, k, 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])); + } + + // assign values of ADM Mass in output box + current_cell.store_vars(Madm, m_c_Madm); + } + + if (m_c_Jadm >= 0) + { + // 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(); + + FOR(i, j, k) + { + Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) * + epsilon[m_dir][j][k] * x[j] * vars.K * + TensorAlgebra::delta(i, k); + + FOR(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.); + } + } + + // assign values of ADM Momentum in output box + current_cell.store_vars(Jadm, m_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 + const int m_c_Madm, m_c_Jadm; + + DIR m_dir; +}; + +#endif /* ADMQUANTITIES_HPP_ */ diff --git a/Source/GRChomboCore/ChomboParameters.hpp b/Source/GRChomboCore/ChomboParameters.hpp index 49349a974..8b3a7b3ec 100644 --- a/Source/GRChomboCore/ChomboParameters.hpp +++ b/Source/GRChomboCore/ChomboParameters.hpp @@ -93,7 +93,7 @@ class ChomboParameters // time stepping outputs and regrid data pp.load("checkpoint_interval", checkpoint_interval, 1); pp.load("plot_interval", plot_interval, 0); - pp.load("stop_time", stop_time, 1.0); + pp.load("stop_time", stop_time, 1000000.0); pp.load("max_steps", max_steps, 1000000); #ifdef CH_USE_HDF5 pp.load("write_plot_ghosts", write_plot_ghosts, false); diff --git a/Source/GRChomboCore/SimulationParametersBase.hpp b/Source/GRChomboCore/SimulationParametersBase.hpp index e0ad7e802..762c2a780 100644 --- a/Source/GRChomboCore/SimulationParametersBase.hpp +++ b/Source/GRChomboCore/SimulationParametersBase.hpp @@ -157,6 +157,16 @@ class SimulationParametersBase : public ChomboParameters pp.load("integral_file_prefix", 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(); + } } } @@ -291,6 +301,14 @@ 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; } } diff --git a/Source/utils/ADMQuantitiesExtraction.hpp b/Source/utils/ADMQuantitiesExtraction.hpp new file mode 100644 index 000000000..da06f5faf --- /dev/null +++ b/Source/utils/ADMQuantitiesExtraction.hpp @@ -0,0 +1,83 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ADMQUANTITIESEXTRACTION_HPP_ +#define ADMQUANTITIESEXTRACTION_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 ADMQuantitiesExtraction : public SphericalExtraction +{ + public: + //! The constructor + 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) + { + 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 + 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) + { + } + + //! Execute the query + void execute_query(AMRInterpolator> *a_interpolator) + { + // extract the values of the ADM mass and spin on the spheres + extract(a_interpolator); + + if (m_params.write_extraction) + write_extraction("ADMQuantitiesExtractionOut_"); + + int num_integrals = (m_c_Madm >= 0) + (m_c_Jadm >= 0); + if (!num_integrals) + return; + + int J_index = m_c_Madm >= 0 ? 1 : 0; + + std::vector> out_integrals(num_integrals); + + 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(); + 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); + } + + private: + const int m_c_Madm, m_c_Jadm; +}; + +#endif /* ADMQUANTITIESEXTRACTION_HPP_ */ diff --git a/Source/utils/WeylExtraction.hpp b/Source/utils/WeylExtraction.hpp index 670effcae..b0d5dadc7 100644 --- a/Source/utils/WeylExtraction.hpp +++ b/Source/utils/WeylExtraction.hpp @@ -85,6 +85,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]); } } }; 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_ */