action_calculation.h

Go to the documentation of this file.

// Copyright (C) 2024 Lisa Biermann, Margarete Mühlleitner, Rui Santos, João
// Viana
//
// SPDX-FileCopyrightText: 2024 Lisa Biermann, Margarete Mühlleitner, Rui
// Santos, João Viana
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#pragma once
 
#include <BSMPT/utility/const_velocity_spline.h>
#include <BSMPT/utility/utility.h>
#include <Eigen/Dense>
#include <gsl/gsl_math.h>
#include <numeric>
#include <optional> // std::optional
#include <sys/stat.h>
#include <sys/types.h>
 
using Eigen::MatrixXd;
using Eigen::VectorXd;
 
namespace BSMPT
{
 
class BounceActionInt
{
private:
  // These two variables check if the method undershot and overshot at least
  // once. If only one type of convergence was achieved it is an indication that
  // no solution has been found.
  bool UndershotOnce = false;
  bool OvershotOnce = false;
 
  double TrueVacuumHessian;
 
  double Initial_lmin;
 
  double Vfalse;
 
  double l0_minus_lmin;
 
  bool PathDeformationConvergedWithout1D = false;
 
  std::optional<double> ExactSolutionThreshold;
 
  double FractionOfThePathExact = 1e-4;
 
public:
  int dim = -1;
 
  enum class ActionStatus
  {
    Success,
    NotCalculated,
    Integration1DFailed,
    PathDeformationNotConverged,
    PathDeformationCrashed,
    FalseVacuumNotMinimum,
    BackwardsPropagationFailed,
    NeverUndershootOvershoot,
    UndershootOvershootNegativeGrad,
    NotEnoughPointsForSpline
  };
 
  enum class UndershootOvershootStatus
  {
    Converged,
    Undershoot,
    Overshoot
  };
 
  enum class Integration1DStatus
  {
    Converged,
    NotConverged,
  };
 
  enum class PathDeformationStatus
  {
    Converged,
    NotConverged,
  };
 
  double Action = -1;
 
  ActionStatus StateOfBounceActionInt = ActionStatus::NotCalculated;
 
  Integration1DStatus StateOf1DIntegration = Integration1DStatus::NotConverged;
 
  PathDeformationStatus StateOfPathDeformation =
      PathDeformationStatus::NotConverged;
 
  double Alpha = 2; // alpha = 2 if T > 0 | alpha = 3 if T = 0
 
  double T = -1;
 
  int MaxPathIntegrations;
 
  int MaxSinglePathDeformations = 200;
 
  std::vector<double> rho_sol;
 
  std::vector<double> l_sol;
 
  std::vector<double> dldrho_sol;
 
  double eps = 0.01;
 
  int BernsteinDegree = 10;
 
  double NumberPathKnots = 50;
  std::vector<double> TrueVacuum;
  std::vector<double> FalseVacuum;
  std::function<double(std::vector<double>)> V; // Potential
  std::function<std::vector<double>(std::vector<double>)>
      dV; // Potential gradient
  std::function<std::vector<std::vector<double>>(std::vector<double>)> Hessian;
  std::vector<std::vector<double>> InitPath; // Initial path given to class
  std::vector<std::vector<double>> Path; // Initial path given to class
  cvspline Spline; // Constant velocity spline object
  tk::spline RasterizeddVdl; // RasterizeddVdl
 
  BounceActionInt();
 
  BounceActionInt(
      std::vector<std::vector<double>> InitPath_In,
      std::vector<double> TrueVacuum_In,
      std::vector<double> FalseVacuum_In,
      std::function<double(std::vector<double>)> &V_In,
      std::function<std::vector<double>(std::vector<double>)> &dV_In,
      double T_In,
      int MaxPathIntegrations_in);
 
  BounceActionInt(std::vector<std::vector<double>> InitPath_In,
                  std::vector<double> TrueVacuum_In,
                  std::vector<double> FalseVacuum_In,
                  std::function<double(std::vector<double>)> &V_In,
                  double T_In,
                  int MaxPathIntegrations_in);
 
  void SetPath(std::vector<std::vector<double>> InitPath_In);
 
  void RasterizedVdl(double l_start = 0);
 
  void PrintVector(std::vector<double> vec);
 
  double Calc_dVdl(double l);
 
  double Calc_d2Vdl2(double l);
 
  std::vector<double> NormalForce(const double &l,
                                  const double &dldrho,
                                  const std::vector<double> &gradient);
 
  void AuxFunctionDev(const double &rho,
                      const std::vector<double> &dvs,
                      std::vector<double> &aks);
 
  void RK5_step(const std::vector<double> &y,
                const std::vector<double> &dydx,
                int n,
                float rho,
                float h,
                std::vector<double> &yout,
                std::vector<double> &yerr);
 
  double BesselI(double alpha, double x, int terms = 100);
 
  double BesselJ(double x, int terms = 100);
 
  std::vector<double>
  ExactSolutionLin(double l0, double l, double dVdl, double d2Vdl2);
 
  std::vector<double> ExactSolutionFromMinimum(double l);
 
  double LogisticFunction(const double &x);
 
  void CalculateExactSolutionThreshold(double MinError = 1e100);
 
  std::vector<double> ExactSolution(double l0);
 
  void BackwardsPropagation();
 
  void IntegrateBounce(double l0,
                       UndershootOvershootStatus &conv,
                       std::vector<double> &rho,
                       std::vector<double> &l,
                       std::vector<double> &dl_drho,
                       std::vector<double> &d2l_drho2,
                       int maxiter,
                       double error,
                       double eps_abs,
                       double max_step = 0);
  void Solve1DBounce(std::vector<double> &rho,
                     std::vector<double> &l,
                     std::vector<double> &dl_drho,
                     std::vector<double> &d2l_drho2,
                     double error = 1e-7,
                     int maxiter  = 100);
 
  double ReductorCalculator(const double &MaximumGradient);
 
  bool PathDeformationCheck(std::vector<double> &l, tk::spline &rho_l_spl);
 
  void SinglePathDeformation(double &stepsize,
                             double &reductor,
                             std::vector<double> &l,
                             tk::spline &rho_l_spl,
                             std::vector<double> &l_fornextpath,
                             std::vector<std::vector<double>> &best_path,
                             std::vector<std::vector<double>> &next_path,
                             double &MaximumGradient,
                             double &MaximumForce,
                             double &MaximumRelativeError,
                             double &Maximum_dldrho,
                             double &PerpendicularGradient,
                             MatrixXd &inverseK,
                             std::vector<std::vector<double>> &forces);
 
  void PathDeformation(std::vector<double> &l, tk::spline &rho_l_spl);
  unsigned nChoosek(unsigned n, unsigned k);
 
  double Bernstein(int n, int nu, double x);
 
  std::vector<double> NormalForceBernstein(const double &dldrho,
                                           const std::vector<double> &gradient,
                                           const std::vector<double> &dphidl,
                                           const std::vector<double> &d2phidl2);
 
  double d2ldrho2(double l, double rho, double dldrho);
 
  double CalculateKineticTermAction(const std::vector<double> &rho,
                                    const tk::spline &dl_drho_spl);
 
  double CalculatePotentialTermAction(const std::vector<double> &rho,
                                      const tk::spline &l_rho_spl);
 
  void CalculateAction(double error = 1e-6);
};
} // namespace BSMPT