shooting.hpp

Go to the documentation of this file.

/*
 * This file is part of BubbleProfiler.
 *
 * BubbleProfiler is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * BubbleProfiler is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with BubbleProfiler.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef BUBBLEPROFILER_SHOOTING_HPP_INCLUDED
#define BUBBLEPROFILER_SHOOTING_HPP_INCLUDED

#include "field_profiles.hpp"

#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

#include <functional>

namespace BubbleProfiler {

class Potential;

class Shooting {
public:
   enum Solver_options : unsigned int {
      Compute_action = 0x01,
      Compute_profile = 0x02
   };

   void set_bisection_precision_bits(int b) { shoot_bisect_bits = b; }
   void set_action_arrived_rel(double tol) { action_arrived_rel = tol; }
   void set_shooting_abs_tol(double tol) { shoot_ode_abs = tol; }
   void set_shooting_rel_tol(double tol) { shoot_ode_rel = tol; }
   void set_action_abs_tol(double tol) { action_ode_abs = tol; }
   void set_action_rel_tol(double tol) { action_ode_rel = tol; }
   void set_drho_frac(double frac) { drho_frac = frac; }
   void set_evolve_change_rel(double frac) { evolve_change_rel = frac; }
   void set_bisect_lambda_max(double l) { bisect_lambda_max = l; }
   void set_max_iterations(int i) { iter_max = i; }
   void set_max_periods(double p) { periods_max = p; }
   void set_f_y_max(double f) { f_y_max = f; }
   void set_f_y_min(double f) { f_y_min = f; }
   void set_y_max(double y) { y_max = y; }
   void solve(const std::function<double(double)>& potential_,
              const std::function<double(double)>& potential_first_,
              const std::function<double(double)>& potential_second_,
              double false_min_, double true_min_, double barrier_,
              int dim_ = 3,
              unsigned int options = (Solver_options::Compute_action |
                                      Solver_options::Compute_profile));
   void solve(const Potential& potential, double false_min_,
              double true_min_, double barrier_, int dim_ = 3,
              unsigned int options = (Solver_options::Compute_action |
                                      Solver_options::Compute_profile));

   const Field_profiles& get_bubble_profile() const;
   double get_euclidean_action() const;

   double quadratic_potential_solution(double) const;
   double quadratic_potential_solution(double, double) const;

private:
   using evolve_type = boost::array<double, 3>;
   using state_type = boost::array<double, 2>;
   using error_stepper_type =
      boost::numeric::odeint::runge_kutta_dopri5<state_type>;
   using controlled_stepper_type =
      boost::numeric::odeint::controlled_runge_kutta<error_stepper_type>;

   int dim{3};

   double lambda_sol{0.};
   double euclidean_action{0.};
   Field_profiles profiles{};
   bool action_computed{false};
   bool profile_computed{false};

   // ODE steppers
   controlled_stepper_type shoot_stepper{};
   controlled_stepper_type action_stepper{};

   double drho_guess{0.};
   double scale{0.};

   double barrier{0.};
   double false_min{0.};
   double true_min{1.};
   double sign_min{1.};
   double p_barrier{0.};
   double p_false_min{0.};
   double p_true_min{0.};

   std::function<double(double)> potential{nullptr};
   std::function<double(double)> potential_first{nullptr};
   std::function<double(double)> potential_second{nullptr};

   // Tolerances and integration method settings
   int shoot_bisect_bits{5};
   double action_arrived_rel{1.e-3};
   double shoot_ode_abs{1.e-4};
   double shoot_ode_rel{1.e-4};
   double action_ode_abs{1.e-6};
   double action_ode_rel{1.e-6};
   double drho_frac{1.e-3};
   double evolve_change_rel{1.e-2};
   double bisect_lambda_max{5.};
   int iter_max{100000};
   double periods_max{1.e2};
   double f_y_max{1.e6};
   double f_y_min{1.e-3};
   double y_max{1.e1};

   void initialize_potential_parameters(
      const std::function<double(double)>& potential_,
      const std::function<double(double)>& potential_first_,
      const std::function<double(double)>& potential_second_,
      double false_min_, double true_min_, double barrier_, int dim_);

   void ode(const state_type &x, state_type &dxdt, double rho);
   double calculate_roll_time(double lambda, double dVdphi,
                                  double d2Vdphi2) const;
   double calculate_roll_velocity(double lambda, double y,
                                  double dVdphi, double d2Vdphi2) const;
   evolve_type evolve(double lambda) const;
   double bubble_scale() const;

   double energy(double phi, double dot_phi) const;

   double shoot(double lambda);

   double unmap(double lambda) const;

   double shooting();

   double integrand(double dot_phi, double rho) const;

   double action(double lambda);

   double interior(double lambda);

   Field_profiles calculate_bubble_profile(double lambda);
};

}  // namespace BubbleProfiler

#endif  // BUBBLEPROFILER_SHOOTING_HPP_INCLUDED