diff --git a/include/BSMPT/bounce_solution/action_calculation.h b/include/BSMPT/bounce_solution/action_calculation.h index 3d71cb23..e512aa3b 100644 --- a/include/BSMPT/bounce_solution/action_calculation.h +++ b/include/BSMPT/bounce_solution/action_calculation.h @@ -78,6 +78,20 @@ class BounceActionInt */ bool PathDeformationConvergedWithout1D = false; + /** + * @brief Lmin - L0 that splits between the two branches of the H > 0 + * analytical solution. If unset, we only use the linear approximation + * + */ + std::optional ExactSolutionThreshold; + + /** + * @brief First guess of l(rho) - l0 used to to integrate in the analytical + * solution. Can be decreased if the error is too high. + * + */ + double FractionOfThePathExact = 1e-4; + public: /** * @brief Dimension of the VEV space @@ -484,11 +498,32 @@ class BounceActionInt * linear in l, i.e. \f$ \frac{dV}{dl} \approx H (l - l_{min}) \f$. This * correspondes to a purely qudratic potential. * - * @param l0 is the integration starting point. * @param l is the integration final point. * @return std::vector */ - std::vector ExactSolutionFromMinimum(double l0, double l); + std::vector ExactSolutionFromMinimum(double l); + + /** + * @brief Logistic function with patched edges to account for numerical + * instability/nans + * + * \f$ \text{LogisticFunction}(x)=\begin{cases}1,\quad x\ge10\\0,\quad x\le + * -10\\\frac{1}{1+e^{-x}},\quad\text{otherwise}\end{cases} \f$ + * + * @param x independent variable + * @return double logistic function at at x + */ + double LogisticFunction(const double &x); + + /** + * @brief Calculate \f$l_\text{threshold}\f$ that splits the two branches of + * analytical integration. This function works recursively until the error + * **MinError** is small enough or the integration step is too small to be + * numerically stable. + * + * @param MinError Lowest error found + */ + void CalculateExactSolutionThreshold(double MinError = 1e100); /** * @brief Calculates the 1D solution by comparing the @ref diff --git a/src/bounce_solution/action_calculation.cpp b/src/bounce_solution/action_calculation.cpp index 73d2498b..a6bdccc1 100644 --- a/src/bounce_solution/action_calculation.cpp +++ b/src/bounce_solution/action_calculation.cpp @@ -263,10 +263,8 @@ double BounceActionInt::BesselJ(double x, int terms) return r; } -std::vector BounceActionInt::ExactSolutionFromMinimum(double l0, - double l) +std::vector BounceActionInt::ExactSolutionFromMinimum(double l) { - (void)l0; // TODO: if l0 is unused, then we can remove it double rho_down = 1e-100; double rho_up = 1; double rho_middle = 0; @@ -352,11 +350,6 @@ std::vector BounceActionInt::ExactSolutionLin(double l0, double rho_up = 1; std::function LinearSolution, LinearSolutionDerivative; - if (abs(dVdl) < 1e-3 or l0 - Initial_lmin < 1e-2) - { - return BounceActionInt::ExactSolutionFromMinimum(l0, l); - } - // Set maximum values that rho can have if (Alpha == 2 and d2Vdl2 > 0) { @@ -526,19 +519,64 @@ double BounceActionInt::Calc_d2Vdl2(double l) ((Hessian(Spline(l)) * Spline.dl(l)) * Spline.dl(l)); } +double BounceActionInt::LogisticFunction(const double &x) +{ + if (x >= 10) return 1; + if (x <= -10) return 0; + return 1 / (1 + exp(-x)); +} + +void BounceActionInt::CalculateExactSolutionThreshold(double MinError) +{ + double NumberOfSteps = 1000; + double Error, l0, l, inital_exponent, final_exponent; + + std::vector MinSol, LinSol; + inital_exponent = log((Spline.L - Initial_lmin) / 100.) - + 10; // Search from 2e-21% spline path + final_exponent = + log((Spline.L - Initial_lmin) / 100.); // Search until 10% spline path + + for (double exponent = inital_exponent; exponent <= final_exponent; + exponent += (final_exponent - inital_exponent) / NumberOfSteps) + { + l0_minus_lmin = exp(exponent); + l0 = Initial_lmin + l0_minus_lmin; + l = l0 + Spline.L * FractionOfThePathExact; + MinSol = ExactSolutionFromMinimum(l); + LinSol = ExactSolutionLin(l0, l, Calc_dVdl(l0), Calc_d2Vdl2(l0)); + Error = 0.5 * + abs((LinSol.at(0) - MinSol.at(0)) / (LinSol.at(0) + MinSol.at(0))); + if (Error < MinError) + { + // Found a better ExactSolutionThreshold + MinError = Error; + ExactSolutionThreshold = l0_minus_lmin; + } + } + + if (MinError > 1e-2) // Error not small enough + { + if (FractionOfThePathExact <= 1e-6) return; + FractionOfThePathExact /= 10.; + CalculateExactSolutionThreshold(MinError); + } +} std::vector BounceActionInt::ExactSolution(double l0) { - double dVdl = Calc_dVdl(l0); + // FractionOfThePathExact = 1e-5; + // Solving FractionOfThePathExact of the path + double l = l0 + Spline.L * FractionOfThePathExact; // Computing the second derivative of the potential beforehand double d2Vdl2 = Calc_d2Vdl2(l0); - // Spline length - double L = Spline.L; + // dVdl + double dVdl = Calc_dVdl(l0); std::stringstream ss; - double l_inf = l0; - double l_sup = l0 + L / 10000.0; // Solving, at most, 1% of the path - double lmiddle = (l_inf + l_sup) / 2; // Middle point + // In the case Backwards propagation failed we use only the Linear Solution + if (not ExactSolutionThreshold.has_value()) + return ExactSolutionLin(l0, l, dVdl, d2Vdl2); if (dVdl <= -1e-3) { @@ -550,10 +588,30 @@ std::vector BounceActionInt::ExactSolution(double l0) return {}; } - std::vector linS = - ExactSolutionLin(l0, lmiddle, dVdl, d2Vdl2); // Linear Solution + std::vector MinSol = ExactSolutionFromMinimum(l); + std::vector LinSol = ExactSolutionLin(l0, l, dVdl, d2Vdl2); + + // Calculate the LogisticExponent which indicates the contribution from both + // branches + double LogisticExponent = 100 * + (l0_minus_lmin - ExactSolutionThreshold.value()) / + (ExactSolutionThreshold.value() - Initial_lmin); + + // If the contribution is too close to 0 or 1 we assume only a single type of + // solution. The protects the code againts NANs from the LinSol when + // the gradient is negative due to numerical instabilities. + if (LogisticFunction(-LogisticExponent) == 1) return MinSol; + if (LogisticFunction(LogisticExponent) == 1) return LinSol; - return linS; + // Interpolate between both solution at ExactSolutionThreshold) + std::vector MinLinInterpolated = LinSol; + + MinLinInterpolated.at(0) = MinSol.at(0) * LogisticFunction(-LogisticExponent); + MinLinInterpolated.at(0) += LinSol.at(0) * LogisticFunction(LogisticExponent); + MinLinInterpolated.at(2) = MinSol.at(2) * LogisticFunction(-LogisticExponent); + MinLinInterpolated.at(2) += LinSol.at(2) * LogisticFunction(LogisticExponent); + + return MinLinInterpolated; } void BounceActionInt::IntegrateBounce(double l0, @@ -672,6 +730,9 @@ double BounceActionInt::BackwardsPropagation() if (abs((l0 - l00) / Spline.L) < 1e-8 and Calc_d2Vdl2(l0) > 0 and l0 <= Spline.L / 100) { + // Calculate the threshold between linear solution and solution from + // minimum + CalculateExactSolutionThreshold(); return l0; } } @@ -686,19 +747,20 @@ double BounceActionInt::BackwardsPropagation() { l00 = l0; l0 -= Calc_dVdl(l0) / 100; + if (abs((l0 - l00) / Spline.L) < 1e-8 and Calc_d2Vdl2(l0) > 0 and + l0 <= Spline.L / 100) + { + // Calculate the threshold between linear solution and solution from + // minimum + CalculateExactSolutionThreshold(); + return l0; + } } - if (l0 <= Spline.L / 100) - { - return l0; - BSMPT::Logger::Write(BSMPT::LoggingLevel::BounceDetailed, ss.str()); - } - ss << "Backwards propagation converged to the other minimum...\t" << l0 - << "\t using minus 0.1% Spline length as backwards propagation\n"; - return (-1 * Spline.L / 100.0); ss << "Backwards propagation not converging\t" << l0 << "\t using minus 0.1% Spline length as backwards propagation\n"; BSMPT::Logger::Write(BSMPT::LoggingLevel::BounceDetailed, ss.str()); - + ExactSolutionThreshold + .reset(); // Use always the linear solution in exact solution return (-1 * Spline.L / 1000.0); }