semi_analytic_ewsb_solver.cpp.in

// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY 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.
//
// FlexibleSUSY 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 FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================

// File generated at @DateAndTime@

/**
 * @file @ModelName@_semi_analytic_ewsb_solver.cpp
 *
 * @brief implementation of EWSB solver for semi-analytic iteration
 *
 * This file was generated at @DateAndTime@ with FlexibleSUSY
 * @FlexibleSUSYVersion@ (git commit: @FlexibleSUSYGitCommit@) and SARAH @SARAHVersion@ .
 */

#include "@ModelName@_semi_analytic_ewsb_solver.hpp"
#include "@ModelName@_mass_eigenstates.hpp"
#include "@ModelName@_semi_analytic_solutions.hpp"
#include "error.hpp"
#include "ew_input.hpp"
#include "logger.hpp"
#include "root_finder.hpp"
#include "fixed_point_iterator.hpp"
#include "raii.hpp"

#include <memory>

namespace flexiblesusy {

#define INPUT(parameter) model.get_input().parameter
#define LOCALINPUT(parameter) INPUT(parameter)
#define INPUTPARAMETER(parameter) LOCALINPUT(parameter)
#define MODELPARAMETER(parameter) model.get_##parameter()
#define EXTRAPARAMETER(parameter) model.get_##parameter()
#define DERIVEDPARAMETER(parameter) model.parameter()
#define PHASE(p) MODELPARAMETER(p)
#define SEMIANALYTICPARAMETER(p) solutions->get_##p()
#define LowEnergyConstant(p) Electroweak_constants::p
#define CLASSNAME @ModelName@_ewsb_solver<Semi_analytic>

/**
 * This method solves the EWSB conditions iteratively, trying several
 * root finding methods until a solution is found.
 */
int CLASSNAME::solve_iteratively(@ModelName@_mass_eigenstates& model_to_solve)
{
   auto model = model_to_solve;
   model.set_ewsb_loop_order(loop_order);

   auto ewsb_stepper = [this, model](const EWSB_vector_t& ewsb_pars) mutable -> EWSB_vector_t {
@getEWSBParametersFromVector@
@setEWSBParametersFromLocalCopies@
@setBoundaryValueParametersFromLocalCopies@
      solutions->evaluate_solutions(model);
@applyEWSBSubstitutions@
      if (this->loop_order > 0)
         model.calculate_DRbar_masses();

      return this->ewsb_step(model);
   };

   auto tadpole_stepper = [this, model](const EWSB_vector_t& ewsb_pars) mutable -> EWSB_vector_t {
@getEWSBParametersFromVector@
@setEWSBParametersFromLocalCopies@
@setBoundaryValueParametersFromLocalCopies@
      solutions->evaluate_solutions(model);
@applyEWSBSubstitutions@
      if (this->loop_order > 0)
         model.calculate_DRbar_masses();

      return this->tadpole_equations(model);
   };

   std::unique_ptr<EWSB_solver> solvers[] = {
@EWSBSolvers@
   };

   const auto x_init(initial_guess(model_to_solve));

   VERBOSE_MSG("\t\tSolving EWSB equations ...");
   VERBOSE_MSG("\t\tInitial guess: x_init = " << x_init.transpose());

   int status;
   for (auto& solver: solvers) {
      VERBOSE_MSG("\t\t\tStarting EWSB iteration using " << solver->name());
      status = solve_iteratively_with(model_to_solve, solver.get(), x_init);
      if (status == EWSB_solver::SUCCESS) {
         VERBOSE_MSG("\t\t\t" << solver->name() << " finished successfully!");
         break;
      }
#ifdef ENABLE_VERBOSE
      else {
         WARNING("\t\t\t" << solver->name() << " could not find a solution!"
                 " (requested precision: " << precision << ")");
      }
#endif
   }

   if (status == EWSB_solver::SUCCESS) {
      model_to_solve.get_problems().unflag_no_ewsb();
   } else {
      set_best_ewsb_solution(model_to_solve, std::begin(solvers), std::end(solvers));
      model_to_solve.get_problems().flag_no_ewsb();
#ifdef ENABLE_VERBOSE
      WARNING("\t\tCould not find a solution to the EWSB equations!"
              " (requested precision: " << precision << ")");
#endif
   }

   return status;
}

/**
 * Solves EWSB equations with given EWSB solver
 *
 * @param model_to_solve model to solve EWSB conditions for
 * @param solver EWSB solver
 * @param x_init initial values
 *
 * @return status of the EWSB solver
 */
int CLASSNAME::solve_iteratively_with(
   @ModelName@_mass_eigenstates& model_to_solve, EWSB_solver* solver, const EWSB_vector_t& x_init)
{
   const int status = solver->solve(x_init);

   if (status == EWSB_solver::SUCCESS)
      set_ewsb_solution(model_to_solve, solver);

   return status;
}

/**
 * Sets EWSB output parameters from given solver.
 *
 * @param solver solver
 */
void CLASSNAME::set_ewsb_solution(@ModelName@_mass_eigenstates& model, const EWSB_solver* solver)
{
   const auto solution = solver->get_solution();

@setEWSBSolution@
@setBoundaryValueParametersFromSolution@
   solutions->evaluate_solutions(model);
@setModelParametersFromEWSB@
   model.calculate_DRbar_masses();
}

/**
 * Sets EWSB output parameters from the solver from the range [first,
 * last), which minimizes the tadpole equations at most.
 *
 * @param first iterator to first solver
 * @param last iterator to last solver
 */
template <typename It>
void CLASSNAME::set_best_ewsb_solution(@ModelName@_mass_eigenstates& model, It first, It last)
{
   auto ma(model), mb(model);

   const auto best_solver =
      std::min_element(first, last,
                       [this, &ma, &mb](const std::unique_ptr<EWSB_solver>& a, const std::unique_ptr<EWSB_solver>& b) {
                          this->set_ewsb_solution(ma, a.get());
                          this->set_ewsb_solution(mb, b.get());
                          return Total(Abs(Re(ma.tadpole_equations()))) < Total(Abs(Re(mb.tadpole_equations())));
                       });

   VERBOSE_MSG("\t\tUsing best solution from " << (*best_solver)->name());

   set_ewsb_solution(model, best_solver->get());
}

int CLASSNAME::solve_iteratively_at(@ModelName@_mass_eigenstates& model_to_solve, int l)
{
   // temporarily set `ewsb_loop_order' to `loop_order' and do
   // iteration
   const auto save_loop_order_raii = make_raii_save(loop_order);
   loop_order = l;

   return solve_iteratively(model_to_solve);
}

int CLASSNAME::solve(@ModelName@_mass_eigenstates& model_to_solve)
{
   if (!solutions) {
      throw SetupError("@ModelName@_ewsb_solver<Semi_analytic>:solve: "
                       "pointer to semi-analytic solutions is zero!");
   }

   if (loop_order == 0) {
      return solve_tree_level(model_to_solve);
   }
   return solve_iteratively_at(model_to_solve, loop_order);
}

int CLASSNAME::solve_tree_level(@ModelName@_mass_eigenstates& model)
{
   int error = 0;

@solveEwsbTreeLevel@
@setTreeLevelSolution@
   return error;
}

CLASSNAME::EWSB_vector_t CLASSNAME::initial_guess(
   const @ModelName@_mass_eigenstates& model) const
{
   EWSB_vector_t x_init(EWSB_vector_t::Zero());

@ewsbInitialGuess@

   return x_init;
}

CLASSNAME::EWSB_vector_t CLASSNAME::tadpole_equations(
   const @ModelName@_mass_eigenstates& model) const
{
   return model.tadpole_equations();
}

/**
 * Calculates EWSB output parameters including loop-corrections.
 *
 * Throws exception of type EEWSBStepFailed if new EWSB parameters are
 * inf or nan.
 *
 * @return new set of EWSB output parameters
 */
CLASSNAME::EWSB_vector_t CLASSNAME::ewsb_step(
   const @ModelName@_mass_eigenstates& model) const
{
   std::array<double, number_of_ewsb_equations> tadpole{};
   EWSB_vector_t ewsb_parameters(EWSB_vector_t::Zero());

   if (loop_order > 0) {
@calculateOneLoopTadpolesNoStruct@
      if (loop_order > 1) {
@calculateTwoLoopTadpolesNoStruct@
      }
   }

@solveEwsbWithTadpoles@

   if (!is_finite)
      throw EEWSBStepFailed();

@fillArrayWithEWSBParameters@

   return ewsb_parameters;
}

void CLASSNAME::set_semi_analytic_solutions(
   @ModelName@_semi_analytic_solutions* s)
{
   solutions = s;
}

} // namespace flexiblesusy