From 4936705b6751a4b63ada129ae8562e26c1850d67 Mon Sep 17 00:00:00 2001 From: ShaunFell Date: Thu, 18 Apr 2024 09:40:27 +0200 Subject: [PATCH] Revamped data reader test. Redoing KdS tests --- tests/KdS_test/ADMFixedBGVars.hpp | 42 ++ tests/KdS_test/AssignFixedBGtoBSSNVars.hpp | 64 +++ tests/KdS_test/BaseProcaFieldTest.hpp | 78 +++ tests/KdS_test/BaseProcaFieldTest.impl.hpp | 141 +++++ tests/KdS_test/EmptyDiagnosticVariables.hpp | 20 + tests/KdS_test/ExcisionTest.hpp | 69 +++ tests/KdS_test/GNUmakefile | 3 +- tests/KdS_test/Main_Tests.cpp | 519 +++++++++++++----- tests/KdS_test/MatterCCZ4.hpp | 123 +++++ tests/KdS_test/MatterCCZ4RHS.hpp | 123 +++++ tests/KdS_test/MatterCCZ4RHS.impl.hpp | 108 ++++ tests/KdS_test/MatterEvolutionTest.hpp | 90 +++ tests/KdS_test/MetricVariablesInterface.hpp | 73 +++ tests/KdS_test/ProcaField.hpp | 132 +++++ tests/KdS_test/ProcaFieldTest.hpp | 131 +++++ tests/KdS_test/UserVariables.hpp | 51 ++ .../GNUmakefile | 0 .../SimpleDataReader_test/Main_datareader.cpp | 83 +++ .../ProcaField.hpp | 0 .../SimulationParameters.hpp | 0 .../UserVariables.hpp | 0 .../TestSimpleDataReader/Main_datareader.cpp | 82 --- 22 files changed, 1722 insertions(+), 210 deletions(-) create mode 100644 tests/KdS_test/ADMFixedBGVars.hpp create mode 100644 tests/KdS_test/AssignFixedBGtoBSSNVars.hpp create mode 100644 tests/KdS_test/BaseProcaFieldTest.hpp create mode 100644 tests/KdS_test/BaseProcaFieldTest.impl.hpp create mode 100644 tests/KdS_test/EmptyDiagnosticVariables.hpp create mode 100644 tests/KdS_test/ExcisionTest.hpp create mode 100644 tests/KdS_test/MatterCCZ4.hpp create mode 100644 tests/KdS_test/MatterCCZ4RHS.hpp create mode 100644 tests/KdS_test/MatterCCZ4RHS.impl.hpp create mode 100644 tests/KdS_test/MatterEvolutionTest.hpp create mode 100644 tests/KdS_test/MetricVariablesInterface.hpp create mode 100644 tests/KdS_test/ProcaField.hpp create mode 100644 tests/KdS_test/ProcaFieldTest.hpp create mode 100644 tests/KdS_test/UserVariables.hpp rename tests/{TestSimpleDataReader => SimpleDataReader_test}/GNUmakefile (100%) create mode 100644 tests/SimpleDataReader_test/Main_datareader.cpp rename tests/{TestSimpleDataReader => SimpleDataReader_test}/ProcaField.hpp (100%) rename tests/{TestSimpleDataReader => SimpleDataReader_test}/SimulationParameters.hpp (100%) rename tests/{TestSimpleDataReader => SimpleDataReader_test}/UserVariables.hpp (100%) delete mode 100644 tests/TestSimpleDataReader/Main_datareader.cpp diff --git a/tests/KdS_test/ADMFixedBGVars.hpp b/tests/KdS_test/ADMFixedBGVars.hpp new file mode 100644 index 0000000..625aaba --- /dev/null +++ b/tests/KdS_test/ADMFixedBGVars.hpp @@ -0,0 +1,42 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + * + * Modified by GRBoondi 2024 + * Please refer to LICENSE in GRBoondi's root directory. + */ + +#ifndef ADMFIXEDBGVARS_HPP_ +#define ADMFIXEDBGVARS_HPP_ + +#include "Tensor.hpp" +#include "UserVariables.hpp" +#include "VarsTools.hpp" + +/// Namespace for ADM vars for fixed BG evolution +namespace ADMFixedBGVars +{ +/// Vars object for ADM vars used in FixedBG evolution +template struct Vars +{ + // ADM vars needed in matter only rhs (ok for Proca and SF) + Tensor<2, data_t> K_tensor; + data_t K; + + data_t lapse; + Tensor<1, data_t> shift; + Tensor<2, data_t> gamma; + + Tensor<1, data_t> d1_lapse; + Tensor<1, Tensor<1, data_t>> d1_shift; + Tensor<2, Tensor<1, data_t>> d1_gamma; + + // Optional second derivatives of the vars + Tensor<2, data_t> d2_lapse; + Tensor<1, Tensor<2, data_t>> d2_shift; + Tensor<2, Tensor<2, data_t>> d2_gamma; +}; + +}; // namespace ADMFixedBGVars + +#endif /* ADMFIXEDBGVARS_HPP_ */ diff --git a/tests/KdS_test/AssignFixedBGtoBSSNVars.hpp b/tests/KdS_test/AssignFixedBGtoBSSNVars.hpp new file mode 100644 index 0000000..3a719d0 --- /dev/null +++ b/tests/KdS_test/AssignFixedBGtoBSSNVars.hpp @@ -0,0 +1,64 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef ASSIGNFIXEDBGTOBSSNVARS_HPP_ +#define ASSIGNFIXEDBGTOBSSNVARS_HPP_ + +#include "ADMFixedBGVars.hpp" +#include "BoxLoops.hpp" +#include "CCZ4RHS.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "MetricVariablesInterface.hpp" +#include "Tensor.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "simd.hpp" + +template class AssignFixedBGtoBSSNVars +{ + public: + AssignFixedBGtoBSSNVars(const background_t a_background, const double a_dx, + const std::array a_center) + : m_dx(a_dx), m_center(a_center), m_background(a_background) + { + } + + template void compute(Cell current_cell) const + { + // compute the ADM vars + const Coordinates coords(current_cell, m_dx, m_center); + ADMFixedBGVars::Vars adm_vars; + m_background.compute_metric_background(adm_vars, coords); + + // assign values to the BSSN vars + using namespace TensorAlgebra; + CCZ4RHS<>::Vars bssn_vars; + + bssn_vars.K = adm_vars.K; + bssn_vars.lapse = adm_vars.lapse; + FOR1(i) { bssn_vars.shift[i] = adm_vars.shift[i]; } + bssn_vars.chi = compute_determinant_sym(adm_vars.gamma); + bssn_vars.chi = pow(bssn_vars.chi, -1.0 / 3.0); + FOR2(i, j) + { + bssn_vars.h[i][j] = adm_vars.gamma[i][j] * bssn_vars.chi; + bssn_vars.A[i][j] = adm_vars.K_tensor[i][j] * bssn_vars.chi - + 1.0 / 3.0 * bssn_vars.K * bssn_vars.h[i][j]; + } + + // set the field to something arbitrary + data_t r = coords.get_radius(); + current_cell.store_vars(0.1 * sin(coords.x) * exp(-r) * coords.z, + c_Avec1); + current_cell.store_vars(bssn_vars); + } + + protected: + const double m_dx; + const std::array m_center; + const background_t m_background; +}; + +#endif /* ASSIGNFIXEDBGTOBSSNVARS_HPP_ */ diff --git a/tests/KdS_test/BaseProcaFieldTest.hpp b/tests/KdS_test/BaseProcaFieldTest.hpp new file mode 100644 index 0000000..a69e814 --- /dev/null +++ b/tests/KdS_test/BaseProcaFieldTest.hpp @@ -0,0 +1,78 @@ +/* GRBoondi 2024 + * Please refer to LICENSE in GRBoondi's root directory. + */ + +/* + Base class for the Proca field. Must be inherited by the Proca field class +defined by the user + +NOTE: + We use the method of 'Curiously Recurring Template Pattern' to allow +arbitrary templated modifications to the theory + https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern + + This allows us to use virtual functions that can be specified by derived +classes. e.g. see ProcaKerrBH Example +*/ +#ifndef BASEPROCAFIELDTEST_H_INCLUDED +#define BASEPROCAFIELDTEST_H_INCLUDED + +#include "ADMFixedBGVars.hpp" //For metric variables +#include "ADMProcaVars.hpp" //For matter variables +#include "CCZ4RHS.hpp" +#include "DefaultBackground.hpp" //Minkowski background as default +#include "FourthOrderDerivatives.hpp" //For calculating derivatives +#include "Tensor.hpp" //For performing tensorial operations +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //For user-defined variables (e.g. see EMKerrBH) +#include "VarsTools.hpp" +#include "simd.hpp" //for SIMD operations + +#include "MetricVariablesInterface.hpp" +#include "ProcaField.hpp" + +template class BaseProcaFieldTest +{ + protected: + const background_t m_background; + + struct params_t + { + }; + + public: + // constructor, inputs are matter params + BaseProcaFieldTest(background_t a_background) + : m_background{a_background} {}; + + template using Vars = ADMProcaVars::MatterVars; + + template + using Diff2Vars = ADMProcaVars::Diff2MatterVars; + + // we set G=0, so this doenst really matter, need it just for using + // MatterCCZ4 + template class vars_t> + emtensor_t compute_emtensor( + const vars_t &vars, //!< the value of the variables + const vars_t> &d1, //!< the value of the 1st derivs + const Tensor<2, data_t> &h_UU, //!< the inverse metric (raised indices) + const Tensor<3, data_t> &chris_ULL) const + { + return emtensor_t{}; + }; //!< the conformal christoffel symbol + + template class vars_t, + template class diff2_vars_t, + template class rhs_vars_t> + void add_matter_rhs( + rhs_vars_t &total_rhs, // RHS terms for all vars + const vars_t &matter_vars, // the value fo the variables + const vars_t> &d1, // value of 1st derivs + const diff2_vars_t> &d2, // 2nd derivs + const vars_t &advec // value of the beta^i d_i(var) terms + ) const; +}; + +#include "BaseProcaFieldTest.impl.hpp" +#endif // BASEPROCAFIELDTEST_H_INCLUDED diff --git a/tests/KdS_test/BaseProcaFieldTest.impl.hpp b/tests/KdS_test/BaseProcaFieldTest.impl.hpp new file mode 100644 index 0000000..692503a --- /dev/null +++ b/tests/KdS_test/BaseProcaFieldTest.impl.hpp @@ -0,0 +1,141 @@ +/* GRBoondi 2024 + * Please refer to LICENSE in GRBoondi's root directory. + */ + +/* +implementation file for ProcaField.hpp +*/ + +#if !defined(BASEPROCAFIELDTEST_H_INCLUDED) +#error "This file should only be included through GeneralizedProcaField.hpp" +#endif + +#ifndef BASEPROCAFIELDTEST_IMPL_H_INCLUDED +#define BASEPROCAFIELDTEST_IMPL_H_INCLUDED + +template +template class vars_t, + template class diff2_vars_t, + template class rhs_vars_t> +void BaseProcaFieldTest::add_matter_rhs( + rhs_vars_t &total_rhs, // RHS terms for all vars + const vars_t &matter_vars, // the value fo the variables + const vars_t> &d1, // value of 1st derivs + const diff2_vars_t> &d2, // 2nd derivs + const vars_t &advec // value of the beta^i d_i(var) terms +) const +{ + + // compute non-conformal metric + CCZ4RHS<>::Vars ccz4_vars; + ccz4_vars.lapse = matter_vars.lapse; + ccz4_vars.chi = matter_vars.chi; + ccz4_vars.K = matter_vars.K; + FOR1(i) { ccz4_vars.shift[i] = matter_vars.shift[i]; }; + FOR2(i, j) + { + ccz4_vars.h[i][j] = matter_vars.h[i][j]; + ccz4_vars.A[i][j] = matter_vars.A[i][j]; + } + + /* ADMFixedBGVars::Vars fixed_bg_vars { + * MetricVarsInterface(ccz4_vars) }; + */ + + ADMFixedBGVars::Vars fixed_bg_vars; + fixed_bg_vars.lapse = matter_vars.lapse; + fixed_bg_vars.K = matter_vars.K; + FOR1(i) + { + fixed_bg_vars.shift[i] = matter_vars.shift[i]; + fixed_bg_vars.d1_lapse[i] = d1.lapse[i]; + }; + FOR2(i, j) + { + fixed_bg_vars.gamma[i][j] = matter_vars.h[i][j] / matter_vars.chi; + fixed_bg_vars.K_tensor[i][j] = + (matter_vars.A[i][j] + + 1.0 / 3.0 * matter_vars.K * matter_vars.h[i][j]) / + matter_vars.chi; + fixed_bg_vars.d1_shift[i][j] = d1.shift[i][j]; + } + FOR3(i, j, k) + { + fixed_bg_vars.d1_gamma[i][j][k] = + d1.h[i][j][k] / matter_vars.chi - + matter_vars.h[i][j] / matter_vars.chi / matter_vars.chi * d1.chi[k]; + } + + auto h_UU = TensorAlgebra::compute_inverse_sym(ccz4_vars.h); + const auto non_phys_chris = + TensorAlgebra::compute_christoffel(d1.h, h_UU).ULL; + const Tensor<3, data_t> chris_phys = TensorAlgebra::compute_phys_chris( + d1.chi, ccz4_vars.chi, ccz4_vars.h, h_UU, non_phys_chris); + + // calculate conformal contravariant metric and conformal christoffel + // symbols + const Tensor<2, data_t> gamma_UU = + TensorAlgebra::compute_inverse(fixed_bg_vars.gamma); + + // covariant derivative of spatial part of Proca field + Tensor<2, data_t> DA; + FOR2(i, j) + { + DA[i][j] = d1.Avec[j][i]; + FOR1(k) { DA[i][j] -= chris_phys[k][i][j] * matter_vars.Avec[k]; } + } + + // evolution equations for spatial part of vector field (index down) + FOR1(i) + { + total_rhs.Avec[i] = -fixed_bg_vars.lapse * d1.phi[i] - + matter_vars.phi * d1.lapse[i] + advec.Avec[i]; + + FOR1(j) + { + total_rhs.Avec[i] += -fixed_bg_vars.lapse * + fixed_bg_vars.gamma[i][j] * + matter_vars.Evec[j] + + matter_vars.Avec[j] * d1.shift[j][i]; + }; + }; + + // evolution equations for Electric vector field (index up) + FOR1(i) + { + total_rhs.Evec[i] = + fixed_bg_vars.lapse * fixed_bg_vars.K * matter_vars.Evec[i] + + advec.Evec[i]; + + FOR1(j) { total_rhs.Evec[i] += -matter_vars.Evec[j] * d1.shift[i][j]; } + + FOR3(j, k, l) + { + total_rhs.Evec[i] += + gamma_UU[j][k] * gamma_UU[i][l] * + (d1.lapse[j] * (d1.Avec[k][l] - d1.Avec[l][k]) + + fixed_bg_vars.lapse * (d2.Avec[k][l][j] - d2.Avec[l][k][j])); + + FOR1(m) + { + total_rhs.Evec[i] += + -fixed_bg_vars.lapse * gamma_UU[j][k] * gamma_UU[i][l] * + (chris_phys[m][j][l] * (d1.Avec[k][m] - d1.Avec[m][k]) + + chris_phys[m][j][k] * (d1.Avec[m][l] - d1.Avec[l][m])); + }; + }; + }; + + // Evolution equations for phi field totally depend on the theory, so we + // leave is up to the user to specify them for their model + total_rhs.phi = 0.; + + // Evolution for auxiliary Z field is also left up to the user in how they + // want to add damping terms + total_rhs.Z = 0.; + + static_cast(this)->matter_rhs_modification( + total_rhs, matter_vars, fixed_bg_vars, d1, d2, advec); +}; + +#endif // BASEPROCAFIELDTEST_IMPL_H_INCLUDED diff --git a/tests/KdS_test/EmptyDiagnosticVariables.hpp b/tests/KdS_test/EmptyDiagnosticVariables.hpp new file mode 100644 index 0000000..377f327 --- /dev/null +++ b/tests/KdS_test/EmptyDiagnosticVariables.hpp @@ -0,0 +1,20 @@ +/* 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 +{ + NUM_DIAGNOSTIC_VARS +}; + +namespace DiagnosticVariables +{ +static const std::array variable_names = {}; +} + +#endif /* DIAGNOSTICVARIABLES_HPP */ diff --git a/tests/KdS_test/ExcisionTest.hpp b/tests/KdS_test/ExcisionTest.hpp new file mode 100644 index 0000000..d2be8a3 --- /dev/null +++ b/tests/KdS_test/ExcisionTest.hpp @@ -0,0 +1,69 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef EXCISIONTEST_HPP_ +#define EXCISIONTEST_HPP_ + +#include "CCZ4.hpp" +// #include "CCZ4Geometry.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "GRInterval.hpp" +#include "Tensor.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "VarsTools.hpp" +#include "simd.hpp" + +//! Does excision for fixed BG BH solutions +//! Note that it is does not using simd so one must set disable_simd() +template class ExcisionTest +{ + // Use matter_t class + using Vars = typename matter_t::template MatterVars; + + protected: + const double m_dx; //!< The grid spacing + const std::array m_center; //!< The BH center + const background_t m_background; + + public: + ExcisionTest(const double a_dx, + const std::array a_center, + background_t a_background) + : m_dx(a_dx), m_center(a_center), m_background(a_background) + { + } + + void compute(Cell current_cell) const + { + const Coordinates coords(current_cell, m_dx, m_center); + bool do_I_excise = m_background.check_if_excised(coords); + if (do_I_excise) + { + // the matter rhs vars within the excision zone + // recalculate them - for now set to zero + Vars matter_vars; + matter_vars.enum_mapping( + [this, ¤t_cell](const int &ivar, double &var) + { var = 0.0; }); + + CCZ4::Vars ccz4_vars; + ccz4_vars.enum_mapping( + [this, ¤t_cell](const int &ivar, double &var) + { var = 0.0; }); + + // assign values of rhs in output box + // also zero the constraints + current_cell.store_vars(matter_vars); + current_cell.store_vars(ccz4_vars); + current_cell.store_vars(0.0, c_Ham); + current_cell.store_vars(0.0, c_Mom1); + current_cell.store_vars(0.0, c_Mom2); + current_cell.store_vars(0.0, c_Mom3); + } // else do nothing + } +}; + +#endif /* EXCISIONTEST_HPP_ */ diff --git a/tests/KdS_test/GNUmakefile b/tests/KdS_test/GNUmakefile index 11591f7..53df05f 100644 --- a/tests/KdS_test/GNUmakefile +++ b/tests/KdS_test/GNUmakefile @@ -16,9 +16,8 @@ LibNames := AMRTimeDependent AMRTools BoxTools src_dirs := $(GRCHOMBO_SOURCE)/utils \ $(GRCHOMBO_SOURCE)/simd \ $(GRCHOMBO_SOURCE)/BoxUtils \ - $(GRCHOMBO_SOURCE)/TaggingCriteria \ - $(GRCHOMBO_SOURCE)/GRChomboCore \ $(GRCHOMBO_SOURCE)/AMRInterpolator \ + $(GRCHOMBO_SOURCE)/CCZ4 \ $(GRBoondi_SOURCE)/Background \ $(GRBoondi_SOURCE)/Matter \ $(GRBoondi_SOURCE)/Diagnostics \ diff --git a/tests/KdS_test/Main_Tests.cpp b/tests/KdS_test/Main_Tests.cpp index b81c27c..a233e39 100644 --- a/tests/KdS_test/Main_Tests.cpp +++ b/tests/KdS_test/Main_Tests.cpp @@ -1,163 +1,430 @@ +/* +GRBoondi -// background includes -#include "ADMFixedBGVars.hpp" -#include "KerrdeSitter.hpp" +Some of these tests were borrowed from the GRDzhadzha test suite. Please see +https://github.com/GRTLCollaboration/GRDzhadzha/tree/main/Tests/KerrBHScalarTest +for more details +*/ + +#ifdef _OPENMP +#include +#endif + +#define NUM_ERR_LIM 0.1 + +// misc includes +#include "AssignFixedBGtoBSSNVars.hpp" +#include "ExcisionTest.hpp" // GRChombo includes #include "Cell.hpp" #include "Coordinates.hpp" +#include "GammaCalculator.hpp" +#include "MatterCCZ4.hpp" +#include "MatterCCZ4RHS.hpp" +#include "NewConstraints.hpp" // Chombo includes +#include "BoxLoops.hpp" +#include "Cell.hpp" +#include "ComputePack.hpp" +#include "DebuggingTools.hpp" +#include "FArrayBox.H" #include "IntVect.H" +#include "SetValue.hpp" #include "UsingNamespace.H" +#include +#include +#include +#include + +// background includes +#include "ADMFixedBGVars.hpp" +#include "BaseProcaFieldTest.hpp" +#include "BoostedBH.hpp" +#include "KerrSchild.hpp" +#include "KerrSchildFixedBG.hpp" +#include "MatterEvolution.hpp" +#include "ProcaField.hpp" + +#include "BaseProcaField.hpp" +#include "ProcaField.hpp" +#include "ProcaFieldTest.hpp" // c++ includes #include #include #include -#define ERR 1e-11 +std::string tab{"\t"}; -int main(int argc, char *argv[]) +void test1() { - std::cout << std::string(10, '#') << std::endl; - std::cout << "Running GRBoondi tests" << std::endl; - std::cout << "current tests: KerrdeSitter" << std::endl; - std::cout << std::string(10, '#') << std::endl; + std::cout << "\n\n KerrSchild test:" << std::endl; - std::string tab{"\t"}; + double dx{1e-3}; - // Here we test the outputs of the background functions + // initialize an IntVect object + IntVect intvect; + intvect[0] = 0.5; + intvect[1] = 0.5; + intvect[2] = 0.5; + + // create a coordinate object + // We set the coordinates to some arbitrary value, but which + // correspond to the + // values in the Mathematica notebook + Coordinates coords( + intvect, dx, + {0, 0, 0}); // should correspond to the point (x,y,z) = (1,1,1) + const double x = coords.x; + const double y = coords.y; + const double z = coords.z; + std::cout << tab << "(x,y,z) = (" << x << "," << y << "," << z << ")" + << std::endl; + + // initialize the kerr params and the kerr object + KerrSchild::params_t kerr_params; + kerr_params.mass = 1.0; + kerr_params.spin = 0.5; + // kerr_params.center = {0.0, 0.0, 0.0}; + KerrSchild kerr_init(kerr_params, dx); + std::cout << tab << " kerr mass = " << kerr_params.mass << std::endl; + std::cout << tab << " kerr spin = " << kerr_params.spin << std::endl; + + // instantiate the metric variables + ADMFixedBGVars::Vars metric_vars; + + // the real work + // init reference class + KerrSchildFixedBG::params_t ref_params = { + kerr_params.mass, kerr_params.center, kerr_params.spin}; + KerrSchildFixedBG kerr_ref(ref_params, dx); + + // Excision test { + // setup container objects for testing excision code + std::vector ref_excision; + std::vector my_excision; + IntVect excision_intvect(1., 1., 1.); - // KerrdeSitter + // First check BH params match + if (kerr_init.m_params.mass != kerr_ref.m_params.mass) { - std::cout << "\n\n KerrdeSitter test:" << std::endl; - - double dx{1e-3}; - - // initialize an IntVect object - IntVect intvect; - intvect[0] = 0.5; - intvect[1] = 0.5; - intvect[2] = 0.5; - - // create a coordinate object - // We set the coordinates to some arbitrary value, but which - // correspond to the - // values in the Mathematica notebook - Coordinates coords( - intvect, dx, - {0, 0, 0}); // should correspond to the point (x,y,z) = (1,1,1) - const double x = coords.x; - const double y = coords.y; - const double z = coords.z; - std::cout << tab << "(x,y,z) = (" << x << "," << y << "," << z - << ")" << std::endl; - - // initialize the kerr params and the kerr object - KerrdeSitter::params_t kerr_params; - kerr_params.mass = 1.4; - kerr_params.spin = 0.9; - kerr_params.cosmo_constant = 0.003; - // kerr_params.center = {0.0, 0.0, 0.0}; - KerrdeSitter kerr_init(kerr_params, dx); - std::cout << tab << " kerr mass = " << kerr_params.mass - << std::endl; - std::cout << tab << " kerr spin = " << kerr_params.spin - << std::endl; - std::cout << tab - << " kerr cosmo_constant = " << kerr_params.cosmo_constant - << std::endl; - - // Discriminant test - double Q{KerrdeSitter::discriminant(kerr_params)}; - std::cout << tab << "Q = " << Q << std::endl; - double mathematica_result_Q{0.0175359611526745}; - if (abs(Q - mathematica_result_Q) > ERR) - { - std::cout << tab << "Discriminant test failed" << std::endl; - } - else - { - std::cout << tab << "Discriminant test passed" << std::endl; - } - - // Metric variable tests - std::cout << "Metric variable tests: " << std::endl; - // instantiate the metric variables - ADMFixedBGVars::Vars metric_vars; + std::cout << "Mass test failed" << std::endl; + } + if (kerr_init.m_params.spin != kerr_ref.m_params.spin) + { + std::cout << "Spin test failed" << std::endl; + } - IntVect mv_intvec; - mv_intvec[0] = 10000; - mv_intvec[1] = 10000; - mv_intvec[2] = 10000; + bool ExcisionPass{true}; + double buffer{0.9554}; + for (int i{0}; i < 3 / dx; i++) + { + excision_intvect[0] = 1.0; + excision_intvect[1] = 1.0; + excision_intvect[2] += i; Box box(IntVect(0, 0, 0), IntVect(8, 8, 8)); FArrayBox fab_in(box, 3); FArrayBox fab_out(box, 3); auto box_pointers = BoxPointers{fab_in, fab_out}; - Cell current_cell(mv_intvec, box_pointers); - Coordinates coords_metricvars(current_cell, dx, - kerr_params.center); - - std::cout << "coords: " << coords_metricvars.x << " " - << coords_metricvars.y << " " << coords_metricvars.z - << std::endl; - - kerr_init.compute_metric_background(metric_vars, coords_metricvars); - - std::cout << "lapse: " << metric_vars.lapse << std::endl; - - std::cout << "shift: " << metric_vars.shift[0] << " " - << metric_vars.shift[1] << " " << metric_vars.shift[2] - << std::endl; - - std::cout << "gamma: " << metric_vars.gamma[0][0] << " " - << metric_vars.gamma[0][1] << " " - << metric_vars.gamma[0][2] << "\n " - << metric_vars.gamma[1][0] << " " - << metric_vars.gamma[1][1] << " " - << metric_vars.gamma[1][2] << "\n " - << metric_vars.gamma[2][0] << " " - << metric_vars.gamma[2][1] << " " - << metric_vars.gamma[2][2] << std::endl; - - std::cout << "metric_vars.d1_lapse: " << metric_vars.d1_lapse[0] - << " " << metric_vars.d1_lapse[1] << " " - << metric_vars.d1_lapse[2] << std::endl; - - std::cout << "metric_vars.d1_shift: " << metric_vars.d1_shift[0][0] - << " " << metric_vars.d1_shift[0][1] << " " - << metric_vars.d1_shift[0][2] << std::endl - << metric_vars.d1_shift[1][0] << " " - << metric_vars.d1_shift[1][1] << " " - << metric_vars.d1_shift[1][2] << std::endl - << metric_vars.d1_shift[2][0] << " " - << metric_vars.d1_shift[2][1] << " " - << metric_vars.d1_shift[2][2] << std::endl; - - std::cout << "metric_vars.d1_gamma: " << std::endl; - FOR2(i, j) + Cell current_cell(excision_intvect, box_pointers); + const Coordinates coords(current_cell, dx, + kerr_params.center); + bool kerr_ref_excise_result{kerr_ref.excise(current_cell) < buffer}; + bool kerr_init_excise_result{ + kerr_init.check_if_excised(coords, buffer)}; + if (kerr_ref_excise_result != kerr_init_excise_result) + { + + std::cout << "Excision test failed" << std::endl; + std::cout << "Coordinate location: " << coords << std::endl; + std::cout << "kerr reference output: " << kerr_ref_excise_result + << std::endl; + std::cout << "kerr my output: " << kerr_init_excise_result + << std::endl; + ExcisionPass = false; + } + } + std::cout << "Excision passed: " << std::boolalpha << ExcisionPass + << std::endl + << std::endl + << std::endl; + } +}; + +// This struct should be a copy of the background param_t struct, with specified +// values +struct +{ + double mass = 1.0; + std::array + center; // automatically filled in during test2 eval + double velocity = 0.5; +} BackgroundParams; + +template int test2() +{ + +#ifdef _OPENMP + std::cout << " Number of threads: " << omp_get_max_threads() << std::endl; +#endif + + int failed{0}; // flag for failed + const bool debug_plots_on{false}; // export data for plotting + const int num_resolutions{2}; // number of resolutions to run at + + // vector of norms for convergence checking + std::array, num_resolutions> error_norms; + + for (int ires{0}; ires < num_resolutions; ++ires) + { + std::cout << std::endl << std::endl; + std::cout << "*******************" << std::endl; + std::cout << "Resolution: " << ires << std::endl; + std::cout << "*******************" << std::endl; + + // fill current error norm array with zeroes + error_norms[ires].fill(0.0); + + // setup the array boxes for various inputs and outputs + const int N_GRID{ + 128 * (int)pow(2, ires)}; // number of cells on each side of box + + std::cout << "Creating grid (" << N_GRID << " x " << N_GRID << " x " + << N_GRID << ")" << std::endl; + // setup boxes + Box box(IntVect(0, 0, 0), IntVect(N_GRID - 1, N_GRID - 1, + N_GRID - 1)); // The computational box + Box ghosted_box( + IntVect(-3, -3, -3), + IntVect(N_GRID + 2, N_GRID + 2, + N_GRID + + 2)); // The computational box with added ghost cells + Box double_ghosted_box( + IntVect(-6, -6, -6), + IntVect( + N_GRID + 5, N_GRID + 5, + N_GRID + + 5)); // The computational box with added double ghost cells + + // setup containers for data + FArrayBox fixedbg_FAB(double_ghosted_box, NUM_VARS); + FArrayBox deriv_fixedbg_FAB(ghosted_box, NUM_VARS); + FArrayBox rhs_FAB(box, NUM_VARS); + FArrayBox fixedbg_rhs_FAB(box, NUM_VARS); + + // loop over box cells and set everything to zero + BoxLoops::loop(make_compute_pack(SetValue(0.0)), fixedbg_FAB, + fixedbg_FAB); + BoxLoops::loop(make_compute_pack(SetValue(0.0)), deriv_fixedbg_FAB, + deriv_fixedbg_FAB); + BoxLoops::loop(make_compute_pack(SetValue(0.0)), rhs_FAB, rhs_FAB); + BoxLoops::loop(make_compute_pack(SetValue(0.0)), fixedbg_rhs_FAB, + fixedbg_rhs_FAB); + + // setup various aspects of the grid + const double length = 16.0; // coordinate size of the computational box + const double dx = length / (N_GRID); // grid spacing + const double center = length / 2.0; // center of the computational box + const std::array center_vector = { + center, center, center}; // center of the computational box + BackgroundParams.center = center_vector; + + // create background + typename background_t::params_t bg_params; + bg_params.mass = BackgroundParams.mass; + bg_params.velocity = BackgroundParams.velocity; + bg_params.center = center_vector; + std::cout << "Mass: " << bg_params.mass << std::endl; + std::cout << "Velocity: " << bg_params.velocity << std::endl; + std::cout << "Center: " << bg_params.center[0] << " " + << bg_params.center[1] << " " << bg_params.center[2] + << std::endl; + background_t background_init(bg_params, dx); + + std::cout << "Computing fixed background..." << std::endl; + // assign background variables to grid + BoxLoops::loop(AssignFixedBGtoBSSNVars(background_init, + dx, center_vector), + fixedbg_FAB, fixedbg_FAB); + GammaCalculator gamamcalc(dx); + BoxLoops::loop(gamamcalc, fixedbg_FAB, deriv_fixedbg_FAB); + + // add derivatives to variables + fixedbg_FAB += deriv_fixedbg_FAB; + + std::cout << "computing constraints..." << std::endl; + + // compute the hamiltonian and momentum constraints and put them in the + // RHS fab + BoxLoops::loop(Constraints(dx, c_Ham, Interval(c_Mom1, c_Mom3)), + fixedbg_FAB, rhs_FAB, box); + + // set CCZ4 parameters + const double G_Newt{0.}; // dont evolve the background + const double sigma{0.}; // turn off Kreiss-Oliger dissipation + CCZ4RHS<>::params_t ccz4_params; + ccz4_params.kappa1 = 0.0; + ccz4_params.kappa2 = 0.0; + ccz4_params.kappa3 = 0.0; + ccz4_params.lapse_coeff = 0.0; // no evolution lapse or shift + ccz4_params.shift_Gamma_coeff = 0.0; + + // setup matter class + ProcaFieldTest::params_t proca_params_test = {1, 1, 1}; + ProcaFieldTest matter(background_init, proca_params_test); + + // setup matterccz4 rhs with matter class + MatterCCZ4RHS matter_ccz4_rhs( + matter, ccz4_params, dx, sigma, CCZ4RHS<>::USE_BSSN, G_Newt); + + std::cout << "Numerically computing rhs..." << std::endl; + + // numerically compute RHS + BoxLoops::loop(matter_ccz4_rhs, fixedbg_FAB, rhs_FAB); + + std::cout << "Analytic computing rhs..." << std::endl; + + // setup analytic Proca field, for computing analytic derivatives + ProcaField::params_t proca_params = {1, 1, 1}; + ProcaField analytic_matter(background_init, proca_params); + + // compute RHS using analytic expressions + MatterEvolution my_an_evolution( + analytic_matter, background_init, sigma, dx, center_vector); + BoxLoops::loop(my_an_evolution, fixedbg_FAB, fixedbg_rhs_FAB); + + // take the difference between the numerically computed RHS and the + // analytically computed RHS. should converge to zero for increasing + // resolution + rhs_FAB -= fixedbg_rhs_FAB; + + std::cout << "Excising..." << std::endl; + + // excise the center where values are always large + ExcisionTest, background_t> + excision(dx, center_vector, background_init); + BoxLoops::loop(excision, rhs_FAB, rhs_FAB, disable_simd()); + + std::cout << "Checking results..." << std::endl; + + // Now we check the results + const int max_norm = 0; + const int L1_norm = 1; + const int num_comps = 1; + const double error_limit = NUM_ERR_LIM; + + // check that the initial data satisfies the constraints + for (int i{c_Ham}; i <= c_Mom3; i++) + { + double max_err = rhs_FAB.norm(max_norm, i, num_comps); + if (max_err > error_limit) + { + std::cout << "CONSTRAINT " << UserVariables::variable_names[i] + << " IS NON ZERO: MAX ERROR = " << max_err + << std::endl; + failed = -1; + } + + // save the L1 norm for convergence checking + error_norms[ires][i] = + rhs_FAB.norm(L1_norm, i, num_comps) * pow(N_GRID, -3); + } + + // compare the numerically calculated and the analytically calculated + // RHS. This checks the expressions for the derivatives of the analytic + // background variables (d1_gamma, d1_lapse, etc.) + for (int i{c_phi}; i <= c_Z; i++) + { + double max_err = rhs_FAB.norm(max_norm, i, num_comps); + if (max_err > error_limit) { - FOR1(k) - { - std::cout << " d1_gamma[" << i << "][" << j << "][" << k - << "] = " << metric_vars.d1_gamma[i][j][k]; - } - std::cout << std::endl; + std::cout + << "ANALYTIC MATTER VARS RHS FOR " + << UserVariables::variable_names[i] + << " DOES NOT MATCH FINITE DIFFERENCE RHS: MAX ERROR = " + << max_err << std::endl; + failed = -1; } - std::cout << "metric_vars.k_tensor: " << std::endl; - FOR2(i, j) + // save for convergence checking + error_norms[ires][i] = + rhs_FAB.norm(L1_norm, i, num_comps) * pow(N_GRID, -3); + } + + } // end of resolution loop + + // check convergence for increasing resolution + double min_convergence_factor = 16.0; + for (int i{0}; i < NUM_VARS; i++) + { + for (int ires{0}; ires < num_resolutions - 1; ires++) + { + double hi_res_norm = error_norms[ires + 1][i]; + double lo_res_norm = error_norms[ires][i]; + // ignore exact zero values + if (abs(hi_res_norm) < 1e-16 && abs(lo_res_norm) < 1e-16) + { + lo_res_norm = 1e-8; + hi_res_norm = 1e-10; + } + double convergence_factor = lo_res_norm / hi_res_norm; + if (convergence_factor < min_convergence_factor) { - std::cout << " k_tensor[" << i << "][" << j - << "]= " << metric_vars.K_tensor[i][j]; - std::cout << std::endl; + min_convergence_factor = convergence_factor; } + if (convergence_factor < 11) + { + failed = -1; + std::cout << "CONVERGENCE FACTOR FOR COMPONENT " + << UserVariables::variable_names[i] << " ON LEVEL " + << ires << " IS LOW: VALUE = " << convergence_factor + << " " << hi_res_norm << " " << lo_res_norm + << std::endl; + } + } + } - } // end of kerr tests + if (failed == 0) + { + std::cout << "The minimum convergence factor was " + << min_convergence_factor << std::endl; + std::cout << "Fixed Background test ........ passed" << std::endl; } + else + { + std::cout << "The minimum convergence factor was " + << min_convergence_factor << std::endl; + std::cout << "Fixed Background test ........ failed!!!" << std::endl; + } + + return failed; +} + +int main(int argc, char *argv[]) +{ + std::cout << std::string(10, '#') << std::endl; + std::cout << "Running GRBoondi tests" << std::endl; + std::cout << std::string(10, '#') << std::endl; + + // Here we test the outputs of the background functions + + std::cout << std::string(10, '#') << std::endl; + std::cout << std::string(10, '#') << std::endl; + std::cout << "current tests: KerrSchild" << std::endl; + std::cout << std::string(10, '#') << std::endl; + std::cout << std::string(10, '#') << std::endl; + test1(); + + std::cout << std::string(10, '#') << std::endl; + std::cout << std::string(10, '#') << std::endl; + std::cout << "current tests: KerrSchild Extended Test" << std::endl; + std::cout << std::string(10, '#') << std::endl; + std::cout << std::string(10, '#') << std::endl; + int test2res = test2(); - return 0; + return test2res; } diff --git a/tests/KdS_test/MatterCCZ4.hpp b/tests/KdS_test/MatterCCZ4.hpp new file mode 100644 index 0000000..72bf3cd --- /dev/null +++ b/tests/KdS_test/MatterCCZ4.hpp @@ -0,0 +1,123 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef MATTERCCZ4RHS_HPP_ +#define MATTERCCZ4RHS_HPP_ + +// #include "CCZ4Geometry.hpp" +#include "CCZ4RHS.hpp" +#include "Cell.hpp" +#include "FourthOrderDerivatives.hpp" +#include "MovingPunctureGauge.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 RHS using CCZ4 including matter terms, and matter variable +//! evolution +/*! + The class calculates the RHS evolution for all the variables. It inherits + from the CCZ4RHS class, which it uses to do the non matter evolution of + variables. It then adds in the additional matter terms to the CCZ4 evolution + (those including the stress energy tensor), and calculates the evolution of + the matter variables. It does not assume a specific form of matter but is + templated over a matter class matter_t. Please see the class ScalarField as + an example of a matter_t. \sa CCZ4RHS(), ScalarField() +*/ + +template +class MatterCCZ4RHS : public CCZ4RHS +{ + public: + // Use this alias for the same template instantiation as this class + using CCZ4 = CCZ4RHS; + + using params_t = CCZ4_params_t; + + template + using MatterVars = typename matter_t::template Vars; + + template + using MatterDiff2Vars = typename matter_t::template Diff2Vars; + + template + using CCZ4Vars = typename CCZ4::template Vars; + + template + using CCZ4Diff2Vars = typename CCZ4::template Diff2Vars; + + // Inherit the variable definitions from CCZ4RHS + matter_t + template + struct Vars : public CCZ4Vars, public MatterVars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Vars::enum_mapping(mapping_function); + MatterVars::enum_mapping(mapping_function); + } + }; + + template + struct Diff2Vars : public CCZ4Diff2Vars, + public MatterDiff2Vars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Diff2Vars::enum_mapping(mapping_function); + MatterDiff2Vars::enum_mapping(mapping_function); + } + }; + + //! Constructor of class MatterCCZ4 + /*! + Inputs are the grid spacing, plus the CCZ4 evolution parameters and a + matter object. It also takes the dissipation parameter sigma, and allows + the formulation to be toggled between CCZ4 and BSSN. The default is CCZ4. + It allows the user to set the value of Newton's constant, which is set to + one by default. + */ + MatterCCZ4RHS(matter_t a_matter, params_t a_params, double a_dx, + double a_sigma, int a_formulation = CCZ4RHS<>::USE_CCZ4, + double a_G_Newton = 1.0); + + //! The compute member which calculates the RHS at each point in the box + //! \sa matter_rhs_equation() + template void compute(Cell current_cell) const; + + protected: + //! The function which adds in the EM Tensor terms to the CCZ4 rhs \sa + //! compute() + template + void add_emtensor_rhs( + Vars + &matter_rhs, //!< the RHS data for each variable at that point. + const Vars &vars, //!< the value of the variables at the point. + const Vars> + &d1 //!< the value of the first derivatives of the variables. + ) const; + + // Class members + matter_t my_matter; //!< The matter object, e.g. a scalar field. + const double m_G_Newton; //!< Newton's constant, set to one by default. +}; + +#include "MatterCCZ4RHS.impl.hpp" + +// This is here for backwards compatibility though the MatterCCZ4RHS +// class should be used in future hence mark as deprecated +template +using MatterCCZ4 [[deprecated("Use MatterCCZ4RHS instead")]] = + MatterCCZ4RHS; + +#endif /* MATTERCCZ4RHS_HPP_ */ diff --git a/tests/KdS_test/MatterCCZ4RHS.hpp b/tests/KdS_test/MatterCCZ4RHS.hpp new file mode 100644 index 0000000..72bf3cd --- /dev/null +++ b/tests/KdS_test/MatterCCZ4RHS.hpp @@ -0,0 +1,123 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef MATTERCCZ4RHS_HPP_ +#define MATTERCCZ4RHS_HPP_ + +// #include "CCZ4Geometry.hpp" +#include "CCZ4RHS.hpp" +#include "Cell.hpp" +#include "FourthOrderDerivatives.hpp" +#include "MovingPunctureGauge.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 RHS using CCZ4 including matter terms, and matter variable +//! evolution +/*! + The class calculates the RHS evolution for all the variables. It inherits + from the CCZ4RHS class, which it uses to do the non matter evolution of + variables. It then adds in the additional matter terms to the CCZ4 evolution + (those including the stress energy tensor), and calculates the evolution of + the matter variables. It does not assume a specific form of matter but is + templated over a matter class matter_t. Please see the class ScalarField as + an example of a matter_t. \sa CCZ4RHS(), ScalarField() +*/ + +template +class MatterCCZ4RHS : public CCZ4RHS +{ + public: + // Use this alias for the same template instantiation as this class + using CCZ4 = CCZ4RHS; + + using params_t = CCZ4_params_t; + + template + using MatterVars = typename matter_t::template Vars; + + template + using MatterDiff2Vars = typename matter_t::template Diff2Vars; + + template + using CCZ4Vars = typename CCZ4::template Vars; + + template + using CCZ4Diff2Vars = typename CCZ4::template Diff2Vars; + + // Inherit the variable definitions from CCZ4RHS + matter_t + template + struct Vars : public CCZ4Vars, public MatterVars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Vars::enum_mapping(mapping_function); + MatterVars::enum_mapping(mapping_function); + } + }; + + template + struct Diff2Vars : public CCZ4Diff2Vars, + public MatterDiff2Vars + { + /// Defines the mapping between members of Vars and Chombo grid + /// variables (enum in User_Variables) + template + void enum_mapping(mapping_function_t mapping_function) + { + CCZ4Diff2Vars::enum_mapping(mapping_function); + MatterDiff2Vars::enum_mapping(mapping_function); + } + }; + + //! Constructor of class MatterCCZ4 + /*! + Inputs are the grid spacing, plus the CCZ4 evolution parameters and a + matter object. It also takes the dissipation parameter sigma, and allows + the formulation to be toggled between CCZ4 and BSSN. The default is CCZ4. + It allows the user to set the value of Newton's constant, which is set to + one by default. + */ + MatterCCZ4RHS(matter_t a_matter, params_t a_params, double a_dx, + double a_sigma, int a_formulation = CCZ4RHS<>::USE_CCZ4, + double a_G_Newton = 1.0); + + //! The compute member which calculates the RHS at each point in the box + //! \sa matter_rhs_equation() + template void compute(Cell current_cell) const; + + protected: + //! The function which adds in the EM Tensor terms to the CCZ4 rhs \sa + //! compute() + template + void add_emtensor_rhs( + Vars + &matter_rhs, //!< the RHS data for each variable at that point. + const Vars &vars, //!< the value of the variables at the point. + const Vars> + &d1 //!< the value of the first derivatives of the variables. + ) const; + + // Class members + matter_t my_matter; //!< The matter object, e.g. a scalar field. + const double m_G_Newton; //!< Newton's constant, set to one by default. +}; + +#include "MatterCCZ4RHS.impl.hpp" + +// This is here for backwards compatibility though the MatterCCZ4RHS +// class should be used in future hence mark as deprecated +template +using MatterCCZ4 [[deprecated("Use MatterCCZ4RHS instead")]] = + MatterCCZ4RHS; + +#endif /* MATTERCCZ4RHS_HPP_ */ diff --git a/tests/KdS_test/MatterCCZ4RHS.impl.hpp b/tests/KdS_test/MatterCCZ4RHS.impl.hpp new file mode 100644 index 0000000..acc69b1 --- /dev/null +++ b/tests/KdS_test/MatterCCZ4RHS.impl.hpp @@ -0,0 +1,108 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#if !defined(MATTERCCZ4RHS_HPP_) +#error "This file should only be included through MatterCCZ4RHS.hpp" +#endif + +#ifndef MATTERCCZ4RHS_IMPL_HPP_ +#define MATTERCCZ4RHS_IMPL_HPP_ +#include "DimensionDefinitions.hpp" + +template +MatterCCZ4RHS::MatterCCZ4RHS( + matter_t a_matter, CCZ4_params_t a_params, + double a_dx, double a_sigma, int a_formulation, double a_G_Newton) + : CCZ4RHS(a_params, a_dx, a_sigma, a_formulation, + 0.0 /*No cosmological constant*/), + my_matter(a_matter), m_G_Newton(a_G_Newton) +{ +} + +template +template +void MatterCCZ4RHS::compute( + Cell current_cell) const +{ + // copy data from chombo gridpoint into local variables + const auto matter_vars = current_cell.template load_vars(); + const auto d1 = this->m_deriv.template diff1(current_cell); + const auto d2 = this->m_deriv.template diff2(current_cell); + const auto advec = + this->m_deriv.template advection(current_cell, matter_vars.shift); + + // Call CCZ4 RHS - work out RHS without matter, no dissipation + Vars matter_rhs; + this->rhs_equation(matter_rhs, matter_vars, d1, d2, advec); + + // add RHS matter terms from EM Tensor + add_emtensor_rhs(matter_rhs, matter_vars, d1); + + // add evolution of matter fields themselves + my_matter.add_matter_rhs(matter_rhs, matter_vars, d1, d2, advec); + + // Add dissipation to all terms + this->m_deriv.add_dissipation(matter_rhs, current_cell, this->m_sigma); + + // Write the rhs into the output FArrayBox + current_cell.store_vars(matter_rhs); +} + +// Function to add in EM Tensor matter terms to CCZ4 rhs +template +template +void MatterCCZ4RHS::add_emtensor_rhs( + Vars &matter_rhs, const Vars &matter_vars, + const Vars> &d1) const +{ + using namespace TensorAlgebra; + + const auto h_UU = compute_inverse_sym(matter_vars.h); + const auto chris = compute_christoffel(d1.h, h_UU); + + // Calculate elements of the decomposed stress energy tensor + const auto emtensor = + my_matter.compute_emtensor(matter_vars, d1, h_UU, chris.ULL); + + // Update RHS for K and Theta depending on formulation + if (this->m_formulation == CCZ4RHS<>::USE_BSSN) + { + matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse * + (emtensor.S + emtensor.rho); + matter_rhs.Theta += 0.0; + } + else + { + matter_rhs.K += 4.0 * M_PI * m_G_Newton * matter_vars.lapse * + (emtensor.S - 3 * emtensor.rho); + matter_rhs.Theta += + -8.0 * M_PI * m_G_Newton * matter_vars.lapse * emtensor.rho; + } + + // Update RHS for other variables + Tensor<2, data_t> Sij_TF = emtensor.Sij; + make_trace_free(Sij_TF, matter_vars.h, h_UU); + + FOR(i, j) + { + matter_rhs.A[i][j] += -8.0 * M_PI * m_G_Newton * matter_vars.chi * + matter_vars.lapse * Sij_TF[i][j]; + } + + FOR(i) + { + data_t matter_term_Gamma = 0.0; + FOR(j) + { + matter_term_Gamma += -16.0 * M_PI * m_G_Newton * matter_vars.lapse * + h_UU[i][j] * emtensor.Si[j]; + } + + matter_rhs.Gamma[i] += matter_term_Gamma; + matter_rhs.B[i] += matter_term_Gamma; + } +} + +#endif /* MATTERCCZ4RHS_IMPL_HPP_ */ diff --git a/tests/KdS_test/MatterEvolutionTest.hpp b/tests/KdS_test/MatterEvolutionTest.hpp new file mode 100644 index 0000000..72ad44e --- /dev/null +++ b/tests/KdS_test/MatterEvolutionTest.hpp @@ -0,0 +1,90 @@ +/* GRChombo + * Copyright 2012 The GRChombo collaboration. + * Please refer to LICENSE in GRChombo's root directory. + */ + +#ifndef MATTEREVOLUTION_HPP_ +#define MATTEREVOLUTION_HPP_ + +#include "ADMFixedBGVars.hpp" +#include "ADMProcaVars.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "FourthOrderDerivatives.hpp" +#include "Tensor.hpp" +#include "TensorAlgebra.hpp" +#include "UserVariables.hpp" //This files needs c_NUM - total number of components +#include "VarsTools.hpp" +#include "simd.hpp" + +//! Calculates RHS of matter variables only, metric vars assumed analytic +/*! + The class calculates the RHS evolution for the matter variables. + It does not assume a specific form of matter or background - + it is templated over a matter class matter_t, and over a background metric, + background_t +*/ + +template class MatterEvolution +{ + public: + //! Inherit the variable definitions from the Matter vars + template using MatterVars = ADMProcaVars::MatterVars; + + // Need d2 of certain matter vars + template + using MatterDiff2Vars = ADMProcaVars::Diff2MatterVars; + + // This is used for the non evolved ADM vars + template using MetricVars = ADMFixedBGVars::Vars; + + //! Constructor of class MatterEvolution + MatterEvolution(matter_t a_matter, background_t a_background, double sigma, + double dx, std::array a_center) + : m_matter(a_matter), m_background(a_background), m_sigma(sigma), + m_deriv(dx), m_dx(dx), m_center(a_center) + { + } + + //! The compute member which calculates the RHS at each point in the box + //! \sa matter_rhs_equation() + template void compute(Cell current_cell) const + { + // copy matter data from chombo gridpoint into local variable + const auto matter_vars = current_cell.template load_vars(); + + // compute the background metric vars + MetricVars metric_vars; + const Coordinates coords(current_cell, m_dx, m_center); + m_background.compute_metric_background(metric_vars, coords); + + // compute derivs for matter grid vars + const auto d1 = m_deriv.template diff1(current_cell); + const auto d2 = m_deriv.template diff2(current_cell); + const auto advec = m_deriv.template advection( + current_cell, metric_vars.shift); + + // the RHS + MatterVars matter_rhs; + VarsTools::assign(matter_rhs, 0.); // All components set to zero + + // add evolution of matter fields and dissipation + m_matter.matter_rhs(matter_rhs, matter_vars, metric_vars, d1, d2, + advec); + m_deriv.add_dissipation(matter_rhs, current_cell, m_sigma); + + // Write the rhs into the output vars for this cell + current_cell.store_vars(matter_rhs); + } + + protected: + const matter_t m_matter; //!< The matter object + const background_t m_background; //!< The metric background + const FourthOrderDerivatives + m_deriv; //!< An object for calculating derivatives of the vars + const double m_sigma; //!< Sigma for dissipation + const double m_dx; //!< Grid spacing + const std::array m_center; //!< Grid center +}; + +#endif /* MATTEREVOLUTION_HPP_ */ diff --git a/tests/KdS_test/MetricVariablesInterface.hpp b/tests/KdS_test/MetricVariablesInterface.hpp new file mode 100644 index 0000000..c35fb61 --- /dev/null +++ b/tests/KdS_test/MetricVariablesInterface.hpp @@ -0,0 +1,73 @@ +/* +GRBoondi +Copyright 2024, Shaun Fell +Please refer to LICENSE in GRBoondi's root directory +*/ + +/* +This function provides a unifying interface between the fixed BG vars and the +CCZ4 ones +*/ + +#ifndef METRICVARIABLESINTERFACE_HPP_ +#define METRICVARIABLESINTERFACE_HPP_ + +#include "ADMFixedBGVars.hpp" +#include "ADMProcaVars.hpp" +#include "BoxLoops.hpp" +#include "CCZ4RHS.hpp" +#include "Cell.hpp" +#include "Coordinates.hpp" +#include "Tensor.hpp" +#include "UserVariables.hpp" //This files needs NUM_VARS - total number of components +#include "simd.hpp" + +template using CCZ4_Vars = CCZ4RHS<>::Vars; + +template using FixedBGVars = ADMFixedBGVars::Vars; + +template +FixedBGVars MetricVarsInterface(CCZ4_Vars ccz4_metric_vars) +{ + FixedBGVars fixed_bg_vars; + + fixed_bg_vars.K = ccz4_metric_vars.K; + fixed_bg_vars.lapse = ccz4_metric_vars.lapse; + FOR1(i) { fixed_bg_vars.shift[i] = ccz4_metric_vars.shift[i]; }; + + FOR2(i, j) + { + fixed_bg_vars.gamma[i][j] = + ccz4_metric_vars.h[i][j] / ccz4_metric_vars.chi; + fixed_bg_vars.K_tensor[i][j] = + (1. / ccz4_metric_vars.chi) * + (ccz4_metric_vars.A[i][j] + + 1. / 3. * ccz4_metric_vars.K * ccz4_metric_vars.h[i][j]); + } + + return fixed_bg_vars; +} + +template +CCZ4_Vars MetricVarsInterface(FixedBGVars fixed_bg_vars) +{ + CCZ4_Vars ccz4_metric_vars; + ccz4_metric_vars.K = fixed_bg_vars.K; + ccz4_metric_vars.lapse = fixed_bg_vars.lapse; + FOR1(i) { ccz4_metric_vars.shift[i] = fixed_bg_vars.shift[i]; }; + auto detgamma = TensorAlgebra::compute_determinant_sym(fixed_bg_vars.gamma); + ccz4_metric_vars.chi = pow(detgamma, -1.0 / 3.0); + + FOR2(i, j) + { + ccz4_metric_vars.h[i][j] = + fixed_bg_vars.gamma[i][j] * ccz4_metric_vars.chi; + ccz4_metric_vars.A[i][j] = + fixed_bg_vars.K_tensor[i][j] * ccz4_metric_vars.chi - + 1.0 / 3.0 * ccz4_metric_vars.K * ccz4_metric_vars.h[i][j]; + } + + return ccz4_metric_vars; +} + +#endif // METRICVARIABLESINTERFACE_HPP_ \ No newline at end of file diff --git a/tests/KdS_test/ProcaField.hpp b/tests/KdS_test/ProcaField.hpp new file mode 100644 index 0000000..b82ca55 --- /dev/null +++ b/tests/KdS_test/ProcaField.hpp @@ -0,0 +1,132 @@ +#ifndef PROCAFIELD_H_INCLUDED +#define PROCAFIELD_H_INCLUDED + +/* +This class adds the simplest L2 lagrangian to the base equations of motion +*/ +#include "ADMFixedBGVars.hpp" +#include "ADMProcaVars.hpp" +#include "BaseProcaField.hpp" +#include "BoostedBH.hpp" +#include "DefaultG.hpp" +#include "KerrSchild.hpp" +#include "L2_simp.hpp" + +// Note: base class BaseProcaField uses CRTP, so pass ProcaField itself as +// template argument +class ProcaField : public BaseProcaField +{ + + protected: + template + using MatterVars = typename ADMProcaVars::MatterVars; + + template + using MatterVarsD2 = typename ADMProcaVars::Diff2MatterVars; + + template + using MetricVars = typename ADMFixedBGVars::template Vars; + + using L2_t = L2; + + public: + struct params_t + { + double mass; + double alpha2; + double vector_damping; + }; + + BoostedBH m_background; + params_t m_params; + L2_t m_L2; + DefaultG m_G2; + + ProcaField(BoostedBH a_background, params_t a_params) + : BaseProcaField(a_background), + m_background(a_background), m_params(a_params) + { + // set up the L2 lagrangian + + DefaultG::params_t G2_params{ + m_params.mass}; // Initialize G2 function parameters + L2_t::params_t L2_params{ + m_params.alpha2}; // Initialize L2 Lagrangian parameters + + DefaultG a_G2(G2_params); + this->m_L2 = L2_t(a_G2, L2_params); + this->m_G2 = a_G2; + }; + + template + void compute_emtensor_modification( + emtensor_t + &base_emtensor, // pass by reference to allow modifications + const MatterVars &vars, const MetricVars &metric_vars, + const MatterVars> &d1, + const MatterVarsD2> &d2, // 2nd derivs + const MatterVars &advec // value of the beta^i d_i(var) terms + ) const + { + // add modifications coming from L2 lagrangian + m_L2.compute_emtensor_modification(base_emtensor, vars, metric_vars, d1, + d2, advec); + }; + + template class rhs_vars_t> + void matter_rhs_modification( + rhs_vars_t &total_rhs, // RHS terms for all vars + const MatterVars &matter_vars, // the value fo the variables + const MetricVars &metric_vars, + const MatterVars> &d1, // value of 1st derivs + const MatterVarsD2> &d2, // 2nd derivs + const MatterVars &advec // value of the beta^i d_i(var) terms + ) const + { + // add modifications coming from L2 lagrangian + m_L2.matter_rhs_modification(total_rhs, matter_vars, metric_vars, d1, + d2, advec); + + // add auxiliary field damping to minimize violation of gauss constraint + + Tensor<2, data_t> gamma_UU{ + TensorAlgebra::compute_inverse_sym(metric_vars.gamma)}; + auto chris_phys{ + TensorAlgebra::compute_christoffel(metric_vars.d1_gamma, gamma_UU) + .ULL}; + + // compute G2 function and its derivatives + data_t V{0}; + data_t dVdA{0}; + data_t dVddA{0}; + m_G2.compute_function(V, dVdA, dVddA, matter_vars, metric_vars); + + // Compute constraint algebra term + data_t gnn{dVdA - 2.0 * dVddA * matter_vars.phi * matter_vars.phi}; + + // Add evolution of auxiliary scalar field + total_rhs.Z += + 2 * metric_vars.lapse * dVdA * matter_vars.phi - + m_params.vector_damping * metric_vars.lapse * matter_vars.Z + + advec.Z; + FOR1(i) + { + total_rhs.Z += metric_vars.lapse * d1.Evec[i][i]; + FOR1(j) + { + total_rhs.Z += metric_vars.lapse * chris_phys[i][i][j] * + matter_vars.Evec[j]; + + // Add damping term to electric field evolution + total_rhs.Evec[i] += + metric_vars.lapse * gamma_UU[i][j] * d1.Z[j]; + } + } + + // Add damping term to scalar field evolution + total_rhs.phi += -metric_vars.lapse * matter_vars.Z * m_params.mass * + m_params.mass / (2 * gnn); + } +}; + +#endif // PROCAFIELD_H_INCLUDED \ No newline at end of file diff --git a/tests/KdS_test/ProcaFieldTest.hpp b/tests/KdS_test/ProcaFieldTest.hpp new file mode 100644 index 0000000..68543a5 --- /dev/null +++ b/tests/KdS_test/ProcaFieldTest.hpp @@ -0,0 +1,131 @@ +#ifndef PROCAFIELDTEST_H_INCLUDED +#define PROCAFIELDTEST_H_INCLUDED + +/* +This class adds the simplest L2 lagrangian to the base equations of motion +*/ +#include "ADMFixedBGVars.hpp" +#include "ADMProcaVars.hpp" +#include "BaseProcaField.hpp" +#include "BoostedBH.hpp" +#include "DefaultG.hpp" +#include "L2_simp.hpp" + +// Note: base class BaseProcaField uses CRTP, so pass ProcaField itself as +// template argument +class ProcaFieldTest : public BaseProcaFieldTest +{ + + protected: + template + using MatterVars = typename ADMProcaVars::MatterVars; + + template + using MatterVarsD2 = typename ADMProcaVars::Diff2MatterVars; + + template + using MetricVars = typename ADMFixedBGVars::template Vars; + + using L2_t = L2; + + public: + struct params_t + { + double mass; + double alpha2; + double vector_damping; + }; + + BoostedBH m_background; + params_t m_params; + L2_t m_L2; + DefaultG m_G2; + + ProcaFieldTest(BoostedBH a_background, params_t a_params) + : BaseProcaFieldTest(a_background), + m_background(a_background), m_params(a_params) + { + // set up the L2 lagrangian + + DefaultG::params_t G2_params{ + m_params.mass}; // Initialize G2 function parameters + L2_t::params_t L2_params{ + m_params.alpha2}; // Initialize L2 Lagrangian parameters + + DefaultG a_G2(G2_params); + this->m_L2 = L2_t(a_G2, L2_params); + this->m_G2 = a_G2; + }; + + template + void compute_emtensor_modification( + emtensor_t + &base_emtensor, // pass by reference to allow modifications + const MatterVars &vars, const MetricVars &metric_vars, + const MatterVars> &d1, + const MatterVarsD2> &d2, // 2nd derivs + const MatterVars &advec // value of the beta^i d_i(var) terms + ) const + { + // add modifications coming from L2 lagrangian + m_L2.compute_emtensor_modification(base_emtensor, vars, metric_vars, d1, + d2, advec); + }; + + template class rhs_vars_t> + void matter_rhs_modification( + rhs_vars_t &total_rhs, // RHS terms for all vars + const MatterVars &matter_vars, // the value fo the variables + const MetricVars &metric_vars, + const MatterVars> &d1, // value of 1st derivs + const MatterVarsD2> &d2, // 2nd derivs + const MatterVars &advec // value of the beta^i d_i(var) terms + ) const + { + // add modifications coming from L2 lagrangian + m_L2.matter_rhs_modification(total_rhs, matter_vars, metric_vars, d1, + d2, advec); + + // add auxiliary field damping to minimize violation of gauss constraint + + Tensor<2, data_t> gamma_UU{ + TensorAlgebra::compute_inverse_sym(metric_vars.gamma)}; + auto chris_phys{ + TensorAlgebra::compute_christoffel(metric_vars.d1_gamma, gamma_UU) + .ULL}; + + // compute G2 function and its derivatives + data_t V{0}; + data_t dVdA{0}; + data_t dVddA{0}; + m_G2.compute_function(V, dVdA, dVddA, matter_vars, metric_vars); + + // Compute constraint algebra term + data_t gnn{dVdA - 2.0 * dVddA * matter_vars.phi * matter_vars.phi}; + + // Add evolution of auxiliary scalar field + total_rhs.Z += + 2 * metric_vars.lapse * dVdA * matter_vars.phi - + m_params.vector_damping * metric_vars.lapse * matter_vars.Z + + advec.Z; + FOR1(i) + { + total_rhs.Z += metric_vars.lapse * d1.Evec[i][i]; + FOR1(j) + { + total_rhs.Z += metric_vars.lapse * chris_phys[i][i][j] * + matter_vars.Evec[j]; + + // Add damping term to electric field evolution + total_rhs.Evec[i] += + metric_vars.lapse * gamma_UU[i][j] * d1.Z[j]; + } + } + + // Add damping term to scalar field evolution + total_rhs.phi += -metric_vars.lapse * matter_vars.Z * m_params.mass * + m_params.mass / (2 * gnn); + } +}; + +#endif // PROCAFIELDTEST_H_INCLUDED \ No newline at end of file diff --git a/tests/KdS_test/UserVariables.hpp b/tests/KdS_test/UserVariables.hpp new file mode 100644 index 0000000..c01379e --- /dev/null +++ b/tests/KdS_test/UserVariables.hpp @@ -0,0 +1,51 @@ +/* GRBoondi 2024 + * Please refer to LICENSE in GRBoondi's root directory. + */ + +#ifndef USERVARIABLES_HPP +#define USERVARIABLES_HPP + +#include "ArrayTools.hpp" +#include "CCZ4UserVariables.hpp" +#include "EmptyDiagnosticVariables.hpp" + +/// This enum gives the index of every variable stored in the grid +enum +{ + c_phi = NUM_CCZ4_VARS, // scalar part of proca field + + c_Avec1, // spatial part of proca field + c_Avec2, + c_Avec3, + + c_Evec1, // conjugate momentum of proca field + c_Evec2, + c_Evec3, + + c_Z, // auxiliary scalar field + + c_Ham, + c_Mom1, + c_Mom2, + c_Mom3, + + NUM_VARS +}; + +namespace UserVariables +{ +static const std::array + user_variable_names = {"phi", + + "Avec1", "Avec2", "Avec3", + + "Evec1", "Evec2", "Evec3", + + "Z"}; + +static const std::array variable_names = + ArrayTools::concatenate(ccz4_variable_names, user_variable_names); + +} // namespace UserVariables + +#endif /* USERVARIABLES_HPP */ diff --git a/tests/TestSimpleDataReader/GNUmakefile b/tests/SimpleDataReader_test/GNUmakefile similarity index 100% rename from tests/TestSimpleDataReader/GNUmakefile rename to tests/SimpleDataReader_test/GNUmakefile diff --git a/tests/SimpleDataReader_test/Main_datareader.cpp b/tests/SimpleDataReader_test/Main_datareader.cpp new file mode 100644 index 0000000..9b707b3 --- /dev/null +++ b/tests/SimpleDataReader_test/Main_datareader.cpp @@ -0,0 +1,83 @@ + +// GRBoondi +#include "DataContainer.hpp" +#include "DataManipulation.hpp" +#include "SetupFunctions.hpp" //For setting up MPI processes +#include "SimpleDataReader.hpp" + +#define SMALL_ERR 1e-4 + +int main(int argc, char *argv[]) +{ + mainSetup(argc, argv); // setup MPI processes + + int reader_failed { 0 }; + int interp_failed { 0 }; + + SimpleDataReader reader{"data_test.dat"}; + DataContainer data = reader.get_data(); + + // Test the nearest neighbor algorithm + std::cout << "###########################################" << std::endl; + std::cout << "Testing nearest neighbor algorithm" << std::endl; + std::cout << "###########################################" << std::endl; + + //get the coordinates and data values + std::vector> coords{data.get_coords()}; + std::vector vals{data.get_data()}; + + //set the query point + std::vector query_point{2.5, 2.5}; + + //find the nearest neighbors + std::pair, std::vector>> + nearest_neighbors{ + DataManipulation::find_nearest_neighbors(coords, query_point, 4)}; + + //extract nearest neighbor coordinates + std::vector> nearest_neighbor_coords{nearest_neighbors.second}; + + //compare against known neighbors + std::vector> ref_data { {2.,2.},{2.,3.},{3.,2.},{3.,3.} }; + if ( nearest_neighbor_coords == ref_data ) + { + std::cout << "Data reader test ..... PASSED" << std::endl; + } else + { + std::cout << "Data reader test ..... FAILED" << std::endl; + reader_failed = 1; + } + + // Test the interpolation algorithm + std::cout << "\n\n###########################################" << std::endl; + std::cout << "Testing interpolation algorithm" << std::endl; + std::cout << "###########################################" << std::endl; + + // create 3D point + std::vector point1{nearest_neighbors.second[0][0], + nearest_neighbors.second[0][1], + vals[nearest_neighbors.first[0]]}; + std::vector point2{nearest_neighbors.second[1][0], + nearest_neighbors.second[1][1], + vals[nearest_neighbors.first[1]]}; + std::vector point3{nearest_neighbors.second[2][0], + nearest_neighbors.second[2][1], + vals[nearest_neighbors.first[2]]}; + + //get interpolated value + double interp_data { DataManipulation::lin_interp_2d(point1, point2, point3, query_point) }; + + //compare against known value + double ref_interp_data { 4. }; + if (interp_data == ref_interp_data) + { + std::cout << "Interpolation test ..... PASSED" << std::endl; + } else + { + std::cout << "Interpolation test ..... FAILED" << std::endl; + interp_failed = 1; + } + + mainFinalize(); // cleanup MPI processes + return interp_failed && reader_failed; +} \ No newline at end of file diff --git a/tests/TestSimpleDataReader/ProcaField.hpp b/tests/SimpleDataReader_test/ProcaField.hpp similarity index 100% rename from tests/TestSimpleDataReader/ProcaField.hpp rename to tests/SimpleDataReader_test/ProcaField.hpp diff --git a/tests/TestSimpleDataReader/SimulationParameters.hpp b/tests/SimpleDataReader_test/SimulationParameters.hpp similarity index 100% rename from tests/TestSimpleDataReader/SimulationParameters.hpp rename to tests/SimpleDataReader_test/SimulationParameters.hpp diff --git a/tests/TestSimpleDataReader/UserVariables.hpp b/tests/SimpleDataReader_test/UserVariables.hpp similarity index 100% rename from tests/TestSimpleDataReader/UserVariables.hpp rename to tests/SimpleDataReader_test/UserVariables.hpp diff --git a/tests/TestSimpleDataReader/Main_datareader.cpp b/tests/TestSimpleDataReader/Main_datareader.cpp deleted file mode 100644 index 0924479..0000000 --- a/tests/TestSimpleDataReader/Main_datareader.cpp +++ /dev/null @@ -1,82 +0,0 @@ - -// GRBoondi -#include "DataContainer.hpp" -#include "DataManipulation.hpp" -#include "SetupFunctions.hpp" //For setting up MPI processes -#include "SimpleDataReader.hpp" - -#define SMALL_ERR 1e-4 - -int main(int argc, char *argv[]) -{ - mainSetup(argc, argv); // setup MPI processes - - SimpleDataReader reader{"KerrdeSitter_rPlus.dat"}; - DataContainer data = reader.get_data(); - - // Test the nearest neighbor algorithm - std::cout << "###########################################" << std::endl; - std::cout << "Testing nearest neighbor algorithm" << std::endl; - std::cout << "###########################################" << std::endl; - - std::vector> coords{data.get_coords()}; - std::vector vals{data.get_data()}; - - std::vector query_point{0.0526, 0.055}; - - std::cout << std::setprecision(15); - - std::pair, std::vector>> - nearest_3_neighbors{ - DataManipulation::find_nearest_neighbors(coords, query_point, 3)}; - std::cout << "Indicies of neighbors: " << std::endl; - for (int i{0}; i < nearest_3_neighbors.first.size(); i++) - { - std::cout << nearest_3_neighbors.first[i] << " "; - } - std::cout << std::endl; - - std::cout << "Coords of neighbors: " << std::endl; - for (int i{0}; i < nearest_3_neighbors.second.size(); i++) - { - for (int j{0}; j < nearest_3_neighbors.second[i].size(); j++) - { - std::cout << nearest_3_neighbors.second[i][j] << " "; - } - std::cout << std::endl; - } - - // Test the interpolation algorithm - std::cout << "###########################################" << std::endl; - std::cout << "Testing interpolation algorithm" << std::endl; - std::cout << "###########################################" << std::endl; - - // create 3D point - std::vector point1{nearest_3_neighbors.second[0][0], - nearest_3_neighbors.second[0][1], - vals[nearest_3_neighbors.first[0]]}; - std::vector point2{nearest_3_neighbors.second[1][0], - nearest_3_neighbors.second[1][1], - vals[nearest_3_neighbors.first[1]]}; - std::vector point3{nearest_3_neighbors.second[2][0], - nearest_3_neighbors.second[2][1], - vals[nearest_3_neighbors.first[2]]}; - - std::cout << "Point 1: " << point1[0] << " " << point1[1] << " " - << point1[2] << std::endl; - std::cout << "Point 2: " << point2[0] << " " << point2[1] << " " - << point2[2] << std::endl; - std::cout << "Point 3: " << point3[0] << " " << point3[1] << " " - << point3[2] << std::endl; - - std::cout << "Query Point: " << query_point[0] << " " << query_point[1] - << std::endl; - - std::cout << "Interpolated value: " - << DataManipulation::lin_interp_2d(point1, point2, point3, - query_point) - << std::endl; - - mainFinalize(); // cleanup MPI processes - return 0.; -} \ No newline at end of file