Skip to content

Commit

Permalink
Rework "ExactSolution()"
Browse files Browse the repository at this point in the history
  • Loading branch information
@vollous
vollous committed Apr 30, 2024
1 parent 88c89b6 commit 65910a4
Show file tree
Hide file tree
Showing 2 changed files with 125 additions and 28 deletions.
39 changes: 37 additions & 2 deletions include/BSMPT/bounce_solution/action_calculation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> 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
Expand Down Expand Up @@ -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<double>
*/
std::vector<double> ExactSolutionFromMinimum(double l0, double l);
std::vector<double> 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
Expand Down
114 changes: 88 additions & 26 deletions src/bounce_solution/action_calculation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,10 +263,8 @@ double BounceActionInt::BesselJ(double x, int terms)
return r;
}

std::vector<double> BounceActionInt::ExactSolutionFromMinimum(double l0,
double l)
std::vector<double> 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;
Expand Down Expand Up @@ -352,11 +350,6 @@ std::vector<double> BounceActionInt::ExactSolutionLin(double l0,
double rho_up = 1;
std::function<double(double)> 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)
{
Expand Down Expand Up @@ -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<double> 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<double> 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)
{
Expand All @@ -550,10 +588,30 @@ std::vector<double> BounceActionInt::ExactSolution(double l0)
return {};
}

std::vector<double> linS =
ExactSolutionLin(l0, lmiddle, dVdl, d2Vdl2); // Linear Solution
std::vector<double> MinSol = ExactSolutionFromMinimum(l);
std::vector<double> 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<double> 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,
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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);
}

Expand Down

0 comments on commit 65910a4

Please sign in to comment.