From 2ff47554b5b591090e770613974a0c5e58915f0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jo=C3=A3o=20Viana?= <57032457+vollous@users.noreply.github.com> Date: Thu, 25 Apr 2024 10:09:57 +0100 Subject: [PATCH] Improve names of some variables, added support for T = 0 bounce solver and added unit tests for it. --- .../bounce_solution/action_calculation.h | 8 +- src/bounce_solution/action_calculation.cpp | 105 +++++++++++++----- tests/unittests/Test-gw.cpp | 58 ++++++++++ 3 files changed, 144 insertions(+), 27 deletions(-) diff --git a/include/BSMPT/bounce_solution/action_calculation.h b/include/BSMPT/bounce_solution/action_calculation.h index d375bc73..9b74d76a 100644 --- a/include/BSMPT/bounce_solution/action_calculation.h +++ b/include/BSMPT/bounce_solution/action_calculation.h @@ -73,7 +73,7 @@ class BounceActionInt * @brief l0 - Initial_lmin for solutions starting very near the true vacuum * */ - double l0Initiallmin; + double l0_minus_lmin; /** * @brief True if path deformation reached the desired results without solving @@ -276,6 +276,12 @@ class BounceActionInt */ tk::spline RasterizeddVdl; // RasterizeddVdl + /** + * @brief Default constructor (unit tests) + * + */ + BounceActionInt(); + /** * @brief Construct a new Bounce Action Int object * diff --git a/src/bounce_solution/action_calculation.cpp b/src/bounce_solution/action_calculation.cpp index 5e4df4b4..5b2c1263 100644 --- a/src/bounce_solution/action_calculation.cpp +++ b/src/bounce_solution/action_calculation.cpp @@ -12,6 +12,8 @@ namespace BSMPT { +BounceActionInt::BounceActionInt(){}; + BounceActionInt::BounceActionInt( std::vector> InitPath_In, std::vector TrueVacuum_In, @@ -418,7 +420,7 @@ void BounceActionInt::RK5_step(std::vector y, return; } -double BounceActionInt::BesselI(double Alpha, double x, int terms) +double BounceActionInt::BesselI(double alpha, double x, int terms) { // This implementation seems to converge quite quicly // https://en.wikipedia.org/wiki/Bessel_function#:~:text=Modified%20Bessel%20functions%3A%20I%CE%B1%2C%20K%CE%B1%5B,first%20and%20second%20kind%20and%20are%20defined%20as%5B19%5D @@ -428,8 +430,8 @@ double BounceActionInt::BesselI(double Alpha, double x, int terms) while ((m < terms) && (abs((r - r0) / r) > 1e-15)) { r0 = r; // Save step - r += 1 / (tgamma(m + Alpha + 1) * tgamma(m + 1)) * - std::pow(x / 2.0, 2.0 * m + Alpha); // One more term + r += 1 / (tgamma(m + alpha + 1) * tgamma(m + 1)) * + std::pow(x / 2.0, 2.0 * m + alpha); // One more term m++; // Update later to not mess up summation } return r; @@ -479,19 +481,37 @@ std::vector BounceActionInt::ExactSolutionFromMinimum(double l0, { LinearSolution = [=](const double rho_in) { - return l - (Initial_lmin + l0Initiallmin * + return l - (Initial_lmin + l0_minus_lmin * sinh(sqrt(TrueVacuumHessian) * rho_in) / (sqrt(TrueVacuumHessian) * rho_in)); }; LinearSolutionDerivative = [=](const double rho_in) { - return ((l0Initiallmin * cosh(rho_in * std::sqrt(TrueVacuumHessian))) / + return ((l0_minus_lmin * cosh(rho_in * std::sqrt(TrueVacuumHessian))) / rho_in) - - (l0Initiallmin * sinh(rho_in * std::sqrt(TrueVacuumHessian))) / + (l0_minus_lmin * sinh(rho_in * std::sqrt(TrueVacuumHessian))) / (pow(rho_in, 2) * std::sqrt(TrueVacuumHessian)); }; } + // T = 0 + if (Alpha == 3) + { + LinearSolution = [=](const double rho_in) + { + return l - + (Initial_lmin + 2 * l0_minus_lmin * + BesselI(1, sqrt(TrueVacuumHessian) * rho_in) / + (sqrt(TrueVacuumHessian) * rho_in)); + }; + LinearSolutionDerivative = [=](const double rho_in) + { + return ( + (2 * l0_minus_lmin * BesselI(2, rho_in * sqrt(TrueVacuumHessian))) / + rho_in); + }; + } + // Check if lower limit is viable assert(LinearSolution(rho_down) > 0); @@ -1022,7 +1042,7 @@ void BounceActionInt::Solve1DBounce( if (mode == 0) { l0 = (lmax + lmin) / 2.0; // Perform binary search - l0Initiallmin = l0 - Initial_lmin; + l0_minus_lmin = l0 - Initial_lmin; IntegrateBounce( l0, conv, rho, l, dl_drho, d2l_drho2, 100000, error, error * 0.0015); if (StateOfBounceActionInt != ActionStatus::NotCalculated) return; @@ -1083,8 +1103,8 @@ void BounceActionInt::Solve1DBounce( if (mode == 1) { mu_middle = (mu_min + mu_max) / 2; - l0Initiallmin = exp(mu_middle); - l0 = Initial_lmin + l0Initiallmin; // Perform binary search in log space + l0_minus_lmin = exp(mu_middle); + l0 = Initial_lmin + l0_minus_lmin; // Perform binary search in log space IntegrateBounce( l0, conv, rho, l, dl_drho, d2l_drho2, 100000, error, error * 0.0015); @@ -1632,17 +1652,34 @@ double BounceActionInt::CalculateKineticTermAction(std::vector &rho, { double integral = 0; double int_delta = rho[rho.size() - 1] / 2000; - for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + if (Alpha == 2) { - // Simpson Integration (1 + 4 + 1)/ 6 * step - integral += r * r * (0.5 * std::pow(dl_drho_spl(r), 2)); - integral += - 4 * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2)); - integral += r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2)); + for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + { + // Simpson Integration (1 + 4 + 1)/ 6 * step + integral += r * r * (0.5 * std::pow(dl_drho_spl(r), 2)); + integral += + 4 * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2)); + integral += r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2)); + } + integral = integral * 4 * M_PI * int_delta / + 6.0; // Angular integration and Simpson step + return integral; + } + else if (Alpha == 3) + { + for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + { + // Simpson Integration (1 + 4 + 1)/ 6 * step + integral += r * r * r * (0.5 * std::pow(dl_drho_spl(r), 2)); + integral += + 4 * r * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2)); + integral += r * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2)); + } + integral = integral * 2 * M_PI * M_PI * int_delta / + 6.0; // Angular integration and Simpson step + return integral; } - integral = integral * 4 * M_PI * int_delta / - 6.0; // Angular integration and Simpson step - return integral; } double BounceActionInt::CalculatePotentialTermAction(std::vector &rho, @@ -1650,16 +1687,32 @@ double BounceActionInt::CalculatePotentialTermAction(std::vector &rho, { double integral = 0; double int_delta = rho[rho.size() - 1] / 2000; - for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + if (Alpha == 2) { - // Simpson Integration (1 + 4 + 1)/ 6 * step - integral += r * r * (V(Spline(l_rho_spl(r)))); - integral += 4 * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0)))); - integral += r * r * (V(Spline(l_rho_spl(r + int_delta)))); + for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + { + // Simpson Integration (1 + 4 + 1)/ 6 * step + integral += r * r * (V(Spline(l_rho_spl(r)))); + integral += 4 * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0)))); + integral += r * r * (V(Spline(l_rho_spl(r + int_delta)))); + } + integral = integral * 4 * M_PI * int_delta / + 6.0; // Angular integration and Simpson step + return integral; + } + else if (Alpha == 3) + { + for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta) + { + // Simpson Integration (1 + 4 + 1)/ 6 * step + integral += r * r * r * (V(Spline(l_rho_spl(r)))); + integral += 4 * r * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0)))); + integral += r * r * r * (V(Spline(l_rho_spl(r + int_delta)))); + } + integral = integral * 2 * M_PI * M_PI * int_delta / + 6.0; // Angular integration and Simpson step + return integral; } - integral = integral * 4 * M_PI * int_delta / - 6.0; // Angular integration and Simpson step - return integral; } void BounceActionInt::CalculateAction( diff --git a/tests/unittests/Test-gw.cpp b/tests/unittests/Test-gw.cpp index 3a2e73d5..f38331fb 100644 --- a/tests/unittests/Test-gw.cpp +++ b/tests/unittests/Test-gw.cpp @@ -19,6 +19,18 @@ using Approx = Catch::Approx; #include // for Logger Class #include +TEST_CASE("Test I_\alpha", "[gw]") +{ + // Tests bounce solver with analytical derivative + using namespace BSMPT; + + BounceActionInt BACalc; + + REQUIRE(BACalc.BesselI(3, 1) == Approx(0.0221684249).epsilon(5e-2)); + REQUIRE(BACalc.BesselI(1, 3) == Approx(3.953370217).epsilon(5e-2)); + REQUIRE(BACalc.BesselI(1, 1.5) == Approx(0.9816664).epsilon(5e-2)); +} + TEST_CASE("Solve bounce equation with analytical derivative", "[gw]") { // Tests bounce solver with analytical derivative @@ -61,6 +73,51 @@ TEST_CASE("Solve bounce equation with analytical derivative", "[gw]") REQUIRE(bc.Action == Approx(4.5011952256).epsilon(5e-2)); } +TEST_CASE( + "Solve bounce equation with analytical derivative and Alpha = 3 (T = 0)", + "[gw]") +{ + // Tests bounce solver with analytical derivative + using namespace BSMPT; + std::function)> V = [&](std::vector x) + { + double c = 5; + double fx = 0; + double fy = 80; + + double r1 = x[0] * x[0] + c * x[1] * x[1]; + double r2 = c * pow(x[0] - 1, 2) + pow(x[1] - 1, 2); + double r3 = fx * (0.25 * pow(x[0], 4) - pow(x[0], 3) / 3.); + r3 += fy * (0.25 * pow(x[1], 4) - pow(x[1], 3) / 3.); + + return (r1 * r2 + r3); + }; + + std::function(std::vector)> dV = + [&](std::vector l0) + { + int dim = 2; + std::vector result(dim); + result = {2 * l0[0] * (5 * pow(-1 + l0[0], 2) + pow(-1 + l0[1], 2)) + + 10 * (-1 + l0[0]) * (pow(l0[0], 2) + 5 * pow(l0[1], 2)), + 10 * (5 * pow(-1 + l0[0], 2) + pow(-1 + l0[1], 2)) * l0[1] + + 2 * (-1 + l0[1]) * (pow(l0[0], 2) + 5 * pow(l0[1], 2)) + + 80 * (-1. * pow(l0[1], 2) + 1. * pow(l0[1], 3))}; + return result; + }; + + std::vector FalseVacuum = {0, 0}; + std::vector TrueVacuum = {1, 1}; + + std::vector> path = {TrueVacuum, FalseVacuum}; + + BounceActionInt bc(path, TrueVacuum, FalseVacuum, V, 0, 6); + bc.Alpha = 3; + bc.CalculateAction(); + + REQUIRE(bc.Action == Approx(13.3129767888).epsilon(5e-2)); +} + TEST_CASE("Solve bounce equation with numerical derivative", "[gw]") { // Tests bounce solver with numerical derivative @@ -226,6 +283,7 @@ TEST_CASE("Solve bounce equation with numerical derivative thin walled", "[gw]") { // Tests bounce solver with numerical derivative using namespace BSMPT; + SetLogger({"--logginglevel::complete"}); std::function)> V = [&](std::vector x) {