From 2dfd2caa9018cd2988e0d44ba3fd12c709da50a6 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Wed, 1 Jul 2015 18:40:12 +0200 Subject: [PATCH] create EWSB solver interface and put EWSB routines into an implementation --- meta/EWSB.m | 77 +++++----- meta/FlexibleSUSY.m | 14 +- templates/ewsb_solver.cpp.in | 162 ++++++++------------ templates/ewsb_solver.hpp.in | 73 ++++++---- templates/mass_eigenstates.cpp.in | 235 ++---------------------------- templates/mass_eigenstates.hpp.in | 12 +- templates/module.mk.in | 2 + 7 files changed, 171 insertions(+), 404 deletions(-) diff --git a/meta/EWSB.m b/meta/EWSB.m index 39a34a572..b93eba5e5 100644 --- a/meta/EWSB.m +++ b/meta/EWSB.m @@ -178,7 +178,8 @@ CConversion`RValueToCFormString[InitialGuessFor[parametersFixedByEWSB[[i]]]] <> ";\n"; ]; - Return[result]; + Parameters`CreateLocalConstRefs[parametersFixedByEWSB] <> "\n" <> + result ]; MakeParameterUnique[(Re|Im)[par_]] := @@ -573,6 +574,18 @@ Return[{Flatten[reducedSolution], freePhases}]; ]; +GetFreeParameter[Rule[a_,b_]] := b; +GetFreeParameter[RuleDelayed[a_,b_]] := b; +GetFreeParameter[p_] := Print["Error: not a rule: ", p]; + +GetFixedParameter[Rule[a_,b_]] := a; +GetFixedParameter[RuleDelayed[a_,b_]] := a; +GetFixedParameter[p_] := Print["Error: not a rule: ", p]; + +(* returns list rhs of rules where the lhs symbols are replaced by unique symbols *) +RemoveFixedParameters[rules_List] := + (GetFreeParameter /@ rules) /. (Rule[#, Unique[CConversion`ToValidCSymbolString[#]]]& /@ (GetFixedParameter /@ rules)); + CreateTreeLevelEwsbSolver[solution_List] := Module[{result = "", body = "", i, par, expr, parStr, oldParStr, reducedSolution, @@ -585,24 +598,15 @@ FlexibleSUSY`Sign[p_] :> Global`LOCALINPUT[CConversion`ToValidCSymbol[FlexibleSUSY`Sign[p]]], FlexibleSUSY`Phase[p_] :> Global`LOCALINPUT[CConversion`ToValidCSymbol[FlexibleSUSY`Phase[p]]] }; - result = Parameters`CreateLocalConstRefsForInputParameters[reducedSolution, "LOCALINPUT"] <> "\n"; - (* save old parameters *) - For[i = 1, i <= Length[reducedSolution], i++, - par = reducedSolution[[i,1]]; - expr = reducedSolution[[i,2]]; - type = CConversion`CreateCType[CConversion`GetScalarElementType[Parameters`GetType[par]]]; - parStr = CConversion`RValueToCFormString[par]; - oldParStr = "old_" <> CConversion`ToValidCSymbolString[par]; - result = result <> - "const " <> type <>" " <> oldParStr <> " = " <> parStr <> ";\n"; - ]; - result = result <> "\n"; + result = Parameters`CreateLocalConstRefs[RemoveFixedParameters[reducedSolution]] <> "\n"; (* write solution *) For[i = 1, i <= Length[reducedSolution], i++, par = reducedSolution[[i,1]]; expr = reducedSolution[[i,2]]; type = CConversion`GetScalarElementType[Parameters`GetType[par]]; - result = result <> Parameters`SetParameter[par, expr, type]; + result = result <> + "const " <> CConversion`CreateCType[type] <> " " <> + Parameters`SetParameter[par, expr, type]; ]; result = result <> "\n"; (* check for errors *) @@ -619,13 +623,14 @@ body = ""; For[i = 1, i <= Length[reducedSolution], i++, par = reducedSolution[[i,1]]; - oldParStr = "old_" <> CConversion`ToValidCSymbolString[par]; - body = body <> Parameters`SetParameter[par, oldParStr, None]; + parStr = CConversion`ToValidCSymbolString[par]; + body = body <> Parameters`SetParameter[par, parStr, "MODEL", None]; ]; - body = body <> "error = 1;\n"; result = result <> - "if (!is_finite) {\n" <> + "if (is_finite) {\n" <> IndentText[body] <> + "} else {\n" <> + IndentText["error = 1;"] <> "\n" <> "}"; , result = "error = solve_ewsb_iteratively(0);\n"; @@ -688,16 +693,16 @@ ]; CreateNewEWSBRootFinder[] := - "new Root_finder(CLASSNAME::tadpole_equations, ¶ms, number_of_ewsb_iterations, ewsb_iteration_precision, "; + "new Root_finder(CLASSNAME::tadpole_equations, ¶ms, get_number_of_iterations(), get_precision(), "; CreateEWSBRootFinder[rootFinder_ /; rootFinder === FlexibleSUSY`FPIRelative] := - "new Fixed_point_iterator(CLASSNAME::ewsb_step, ¶ms, number_of_ewsb_iterations, ewsb_iteration_precision)"; + "new Fixed_point_iterator(CLASSNAME::step, ¶ms, get_number_of_iterations(), get_precision())"; CreateEWSBRootFinder[rootFinder_ /; rootFinder === FlexibleSUSY`FPIAbsolute] := - "new Fixed_point_iterator(CLASSNAME::ewsb_step, ¶ms, number_of_ewsb_iterations, ewsb_iteration_precision)"; + "new Fixed_point_iterator(CLASSNAME::step, ¶ms, get_number_of_iterations(), get_precision())"; CreateEWSBRootFinder[rootFinder_ /; rootFinder === FlexibleSUSY`FPITadpole] := - "new Fixed_point_iterator(CLASSNAME::ewsb_step, ¶ms, number_of_ewsb_iterations, fixed_point_iterator::Convergence_tester_tadpole(ewsb_iteration_precision, CLASSNAME::tadpole_equations, ¶ms))"; + "new Fixed_point_iterator(CLASSNAME::step, ¶ms, get_number_of_iterations(), fixed_point_iterator::Convergence_tester_tadpole(get_precision(), CLASSNAME::tadpole_equations, ¶ms))"; CreateEWSBRootFinder[rootFinder_ /; rootFinder === FlexibleSUSY`GSLHybrid] := CreateNewEWSBRootFinder[] <> "gsl_multiroot_fsolver_hybrid)"; @@ -726,14 +731,17 @@ WrapPhase[phase_, str_String] := "LOCALINPUT(" <> CConversion`ToValidCSymbolString[phase] <> ")*Abs(" <> str <> ")"; -SetEWSBSolution[par_, idx_, phase_, func_String] := +SetEWSBSolution[par_, idx_, phase_, model_String /; model === "", func_String] := Parameters`SetParameter[par, WrapPhase[phase, func <> "(" <> ToString[idx-1] <> ")"], None]; -SetEWSBSolution[parametersFixedByEWSB_List, freePhases_List, func_String] := +SetEWSBSolution[par_, idx_, phase_, model_String, func_String] := + Parameters`SetParameter[par, WrapPhase[phase, func <> "(" <> ToString[idx-1] <> ")"], model, None]; + +SetEWSBSolution[parametersFixedByEWSB_List, freePhases_List, model_String, func_String] := Module[{result = "", i, phase}, For[i = 1, i <= Length[parametersFixedByEWSB], i++, phase = FindFreePhase[parametersFixedByEWSB[[i]], freePhases]; - result = result <> SetEWSBSolution[parametersFixedByEWSB[[i]], i, phase, func]; + result = result <> SetEWSBSolution[parametersFixedByEWSB[[i]], i, phase, model, func]; ]; result ]; @@ -768,7 +776,7 @@ ]; CreateEwsbSolverWithTadpoles[solution_List, softHiggsMassToTadpoleAssociation_List] := - Module[{result = "", i, par, expr, parStr, reducedSolution, rules, type}, + Module[{result = "", i, par, expr, parStr, reducedSolution, rules, type, ctype}, reducedSolution = solution /. FlexibleSUSY`tadpole[p_] :> CConversion`ReleaseHoldAt[HoldForm[FlexibleSUSY`tadpole[[p-1]]], {1,2}]; If[reducedSolution =!= {}, @@ -778,28 +786,19 @@ FlexibleSUSY`Sign[p_] :> Global`LOCALINPUT[CConversion`ToValidCSymbol[FlexibleSUSY`Sign[p]]], FlexibleSUSY`Phase[p_] :> Global`LOCALINPUT[CConversion`ToValidCSymbol[FlexibleSUSY`Phase[p]]] }; - result = Parameters`CreateLocalConstRefsForInputParameters[reducedSolution, "LOCALINPUT"]; - (* define variables for new parameters *) - For[i = 1, i <= Length[reducedSolution], i++, - par = reducedSolution[[i,1]]; - expr = reducedSolution[[i,2]]; - type = CConversion`CreateCType[CConversion`GetScalarElementType[Parameters`GetType[par]]]; - parStr = CConversion`ToValidCSymbolString[par]; - result = result <> type <> " " <> parStr <> ";\n"; - ]; - result = result <> "\n"; + result = Parameters`CreateLocalConstRefs[RemoveFixedParameters[reducedSolution]] <> "\n"; (* write solution *) For[i = 1, i <= Length[reducedSolution], i++, par = reducedSolution[[i,1]]; expr = reducedSolution[[i,2]]; type = CConversion`GetScalarElementType[Parameters`GetType[par]]; + ctype = CConversion`CreateCType[type]; parStr = CConversion`ToValidCSymbolString[par]; - result = result <> parStr <> " = " <> + result = result <> "const " <> ctype <> " " <> parStr <> " = " <> CConversion`CastTo[CConversion`RValueToCFormString[expr],type] <> ";\n"; ]; - result = result <> "\n"; (* check for errors *) - result = result <> "const bool is_finite = "; + result = result <> "\n" <> "const bool is_finite = "; For[i = 1, i <= Length[reducedSolution], i++, par = reducedSolution[[i,1]]; parStr = CConversion`ToValidCSymbolString[par]; diff --git a/meta/FlexibleSUSY.m b/meta/FlexibleSUSY.m index b5a77474c..28f987d08 100644 --- a/meta/FlexibleSUSY.m +++ b/meta/FlexibleSUSY.m @@ -752,7 +752,7 @@ corresponding tadpole is real or imaginary (only in models with CP massCalculationPrototypes = "", massCalculationFunctions = "", calculateAllMasses = "", calculateOneLoopTadpoles = "", calculateTwoLoopTadpoles = "", - calculateOneLoopTadpolesNoStruct = "", calculateTwoLoopTadpolesNoStruct = "", + calculateOneLoopTadpolesFromModel = "", calculateTwoLoopTadpolesFromModel = "", selfEnergyPrototypes = "", selfEnergyFunctions = "", twoLoopTadpolePrototypes = "", twoLoopTadpoleFunctions = "", twoLoopSelfEnergyPrototypes = "", twoLoopSelfEnergyFunctions = "", @@ -830,11 +830,11 @@ corresponding tadpole is real or imaginary (only in models with CP higgsToEWSBEqAssociation = CreateHiggsToEWSBEqAssociation[]; oneLoopTadpoles = Cases[nPointFunctions, SelfEnergies`Tadpole[___]]; calculateOneLoopTadpoles = SelfEnergies`FillArrayWithOneLoopTadpoles[higgsToEWSBEqAssociation, "tadpole", "-"]; - calculateOneLoopTadpolesNoStruct = SelfEnergies`FillArrayWithOneLoopTadpoles[higgsToEWSBEqAssociation, "tadpole", "+"]; + calculateOneLoopTadpolesFromModel = SelfEnergies`FillArrayWithOneLoopTadpoles[higgsToEWSBEqAssociation, "tadpole", "+", "MODEL->"]; If[SARAH`UseHiggs2LoopMSSM === True || FlexibleSUSY`UseHiggs2LoopNMSSM === True, calculateTwoLoopTadpoles = SelfEnergies`FillArrayWithTwoLoopTadpoles[SARAH`HiggsBoson, "tadpole", "-"]; - calculateTwoLoopTadpolesNoStruct = SelfEnergies`FillArrayWithTwoLoopTadpoles[SARAH`HiggsBoson, "tadpole", "+"]; + calculateTwoLoopTadpolesFromModel = SelfEnergies`FillArrayWithTwoLoopTadpoles[SARAH`HiggsBoson, "tadpole", "+", "MODEL->"]; {thirdGenerationHelperPrototypes, thirdGenerationHelperFunctions} = TreeMasses`CreateThirdGenerationHelpers[]; ]; If[SARAH`UseHiggs2LoopSM === True, @@ -903,8 +903,8 @@ corresponding tadpole is real or imaginary (only in models with CP solveEWSBTemporarily = "solve_ewsb_tree_level();"; ]; EWSBSolvers = EWSB`CreateEWSBRootFinders[FlexibleSUSY`FSEWSBSolvers]; - setEWSBSolution = EWSB`SetEWSBSolution[parametersFixedByEWSB, freePhases, "solver->get_solution"]; - fillArrayWithEWSBParameters = EWSB`FillArrayWithParameters["ewsb_parameters", parametersFixedByEWSB]; + setEWSBSolution = EWSB`SetEWSBSolution[parametersFixedByEWSB, freePhases, "MODEL", "solver->get_solution"]; + fillArrayWithEWSBParameters = EWSB`FillArrayWithParameters["parameters", parametersFixedByEWSB]; solveEwsbWithTadpoles = EWSB`CreateEwsbSolverWithTadpoles[ewsbSolution, softHiggsMasses]; getEWSBParametersFromGSLVector = EWSB`GetEWSBParametersFromGSLVector[parametersFixedByEWSB, freePhases, "x"]; setEWSBParametersFromLocalCopies = EWSB`SetEWSBParametersFromLocalCopies[parametersFixedByEWSB, "model"]; @@ -927,8 +927,8 @@ corresponding tadpole is real or imaginary (only in models with CP "@calculateTreeLevelTadpoles@" -> IndentText[calculateTreeLevelTadpoles], "@calculateOneLoopTadpoles@" -> IndentText[calculateOneLoopTadpoles], "@calculateTwoLoopTadpoles@" -> IndentText[calculateTwoLoopTadpoles], - "@calculateOneLoopTadpolesNoStruct@" -> IndentText[calculateOneLoopTadpolesNoStruct], - "@calculateTwoLoopTadpolesNoStruct@" -> IndentText[calculateTwoLoopTadpolesNoStruct], + "@calculateOneLoopTadpolesFromModel@" -> IndentText[calculateOneLoopTadpolesFromModel], + "@calculateTwoLoopTadpolesFromModel@" -> IndentText[calculateTwoLoopTadpolesFromModel], "@clearOutputParameters@" -> IndentText[clearOutputParameters], "@clearPhases@" -> IndentText[clearPhases], "@copyDRbarMassesToPoleMasses@" -> IndentText[copyDRbarMassesToPoleMasses], diff --git a/templates/ewsb_solver.cpp.in b/templates/ewsb_solver.cpp.in index 0509c45f9..83a7748fe 100644 --- a/templates/ewsb_solver.cpp.in +++ b/templates/ewsb_solver.cpp.in @@ -21,45 +21,33 @@ #include "@ModelName@_ewsb_solver.hpp" #include "@ModelName@_mass_eigenstates.hpp" +#include "ewsb_solver.hpp" +#include "fixed_point_iterator.hpp" +#include "functors.hpp" +#include "logger.hpp" +#include "minimizer.hpp" +#include "numerics2.hpp" +#include "root_finder.hpp" +#include "wrappers.hpp" #include "gsl_utils.hpp" #include -namespace flexiblesusy { - -@ModelName@_ewsb_solver::@ModelName@_ewsb_solver() - : number_of_ewsb_iterations(100) - , ewsb_loop_order(2) - , ewsb_iteration_precision(1.0e-5) -{ -} - -~@ModelName@_ewsb_solver::@ModelName@_ewsb_solver() -{ -} - -void @ModelName@_ewsb_solver::set_ewsb_loop_order(unsigned loop_order) -{ - ewsb_loop_order = loop_order; -} - -void @ModelName@_ewsb_solver::set_number_of_ewsb_iterations(std::size_t iterations) -{ - number_of_ewsb_iterations = iterations; -} +#define CLASSNAME @ModelName@_ewsb_solver +#define MODEL get_model() +#define MODELPARAMETER(parameter) MODEL->get_##parameter() +#define INPUT(parameter) model->get_input().parameter +#define LOCALINPUT(parameter) MODEL->get_input().parameter +#define INPUTPARAMETER(parameter) LOCALINPUT(parameter) -void @ModelName@_ewsb_solver::set_ewsb_iteration_precision(double precision) -{ - ewsb_iteration_precision = precision; -} +namespace flexiblesusy { -double @ModelName@_ewsb_solver::get_ewsb_iteration_precision() const +@ModelName@_ewsb_solver::@ModelName@_ewsb_solver(@ModelName@_mass_eigenstates* model_) + : @ModelName@_ewsb_solver_interface(model_) { - return ewsb_iteration_precision; } -double @ModelName@_ewsb_solver::get_ewsb_loop_order() const +@ModelName@_ewsb_solver::~@ModelName@_ewsb_solver() { - return ewsb_loop_order; } /** @@ -79,14 +67,14 @@ int @ModelName@_ewsb_solver::tadpole_equations(const gsl_vector* x, void* params return GSL_EDOM; } - const @ModelName@_ewsb_solver::EWSB_args* ewsb_args + const @ModelName@_ewsb_solver::EWSB_args* args = static_cast<@ModelName@_ewsb_solver::EWSB_args*>(params); - @ModelName@_mass_eigenstates* model = ewsb_args->model; - const unsigned ewsb_loop_order = ewsb_args->ewsb_loop_order; + @ModelName@_mass_eigenstates* model = args->model; + const unsigned loop_order = args->loop_order; @setEWSBParametersFromGSLVector@ - if (ewsb_loop_order > 0) + if (loop_order > 0) model->calculate_DRbar_masses(); double tadpole[number_of_ewsb_equations] = { 0. }; @@ -99,29 +87,13 @@ int @ModelName@_ewsb_solver::tadpole_equations(const gsl_vector* x, void* params return is_finite(tadpole) ? GSL_SUCCESS : GSL_EDOM; } -/** - * Method which calculates the tadpoles at the current loop order. - * - * @param tadpole array of tadpole - */ -void @ModelName@_ewsb_solver::tadpole_equations(double tadpole[number_of_ewsb_equations]) const -{ -@calculateTreeLevelTadpoles@ - if (ewsb_loop_order > 0) { -@calculateOneLoopTadpoles@ - if (ewsb_loop_order > 1) { -@calculateTwoLoopTadpoles@ - } - } -} - /** * This method solves the EWSB conditions iteratively, trying several * root finding methods until a solution is found. */ -int @ModelName@_ewsb_solver::solve_ewsb_iteratively() +int @ModelName@_ewsb_solver::solve_iteratively() { - EWSB_args params = {this, ewsb_loop_order}; + EWSB_args params = {get_model(), this, get_loop_order()}; EWSB_solver* solvers[] = { @EWSBSolvers@ @@ -129,7 +101,7 @@ int @ModelName@_ewsb_solver::solve_ewsb_iteratively() const std::size_t number_of_solvers = sizeof(solvers)/sizeof(*solvers); double x_init[number_of_ewsb_equations] = { 0. }; - ewsb_initial_guess(x_init); + initial_guess(x_init); #ifdef ENABLE_VERBOSE std::cout << "Solving EWSB equations ...\n" @@ -142,7 +114,7 @@ int @ModelName@_ewsb_solver::solve_ewsb_iteratively() int status; for (std::size_t i = 0; i < number_of_solvers; ++i) { VERBOSE_MSG("\tStarting EWSB iteration using solver " << i); - status = solve_ewsb_iteratively_with(solvers[i], x_init); + status = solve_iteratively_with(solvers[i], x_init); if (status == EWSB_solver::SUCCESS) { VERBOSE_MSG("\tSolver " << i << " finished successfully!"); break; @@ -150,18 +122,18 @@ int @ModelName@_ewsb_solver::solve_ewsb_iteratively() #ifdef ENABLE_VERBOSE else { WARNING("\tSolver " << i << " could not find a solution!" - " (requested precision: " << ewsb_iteration_precision << ")"); + " (requested precision: " << iteration_precision << ")"); } #endif } if (status == EWSB_solver::SUCCESS) { - problems.unflag_no_ewsb(); + MODEL->get_problems().unflag_no_ewsb(); } else { - problems.flag_no_ewsb(); + MODEL->get_problems().flag_no_ewsb(); #ifdef ENABLE_VERBOSE WARNING("\tCould not find a solution to the EWSB equations!" - " (requested precision: " << ewsb_iteration_precision << ")"); + " (requested precision: " << iteration_precision << ")"); #endif } @@ -178,7 +150,7 @@ int @ModelName@_ewsb_solver::solve_ewsb_iteratively() * * @return status of the EWSB solver */ -int @ModelName@_ewsb_solver::solve_ewsb_iteratively_with( +int @ModelName@_ewsb_solver::solve_iteratively_with( EWSB_solver* solver, const double x_init[number_of_ewsb_equations] ) @@ -190,19 +162,19 @@ int @ModelName@_ewsb_solver::solve_ewsb_iteratively_with( return status; } -int @ModelName@_ewsb_solver::solve_ewsb_iteratively(unsigned loop_order) +int @ModelName@_ewsb_solver::solve_iteratively(unsigned loop_order) { - // temporarily set `ewsb_loop_order' to `loop_order' and do + // temporarily set `loop_order' to `loop_order' and do // iteration - const unsigned old_loop_order = ewsb_loop_order; - ewsb_loop_order = loop_order; - const int status = solve_ewsb_iteratively(); - ewsb_loop_order = old_loop_order; + const unsigned old_loop_order = loop_order; + loop_order = loop_order; + const int status = solve_iteratively(); + loop_order = old_loop_order; return status; } -int @ModelName@_ewsb_solver::solve_ewsb_tree_level() +int @ModelName@_ewsb_solver::solve_tree_level() { int error = 0; @@ -211,31 +183,22 @@ int @ModelName@_ewsb_solver::solve_ewsb_tree_level() return error; } -int @ModelName@_ewsb_solver::solve_ewsb_tree_level_via_soft_higgs_masses() -{ - int error = 0; - -@solveTreeLevelEWSBviaSoftHiggsMasses@ - - return error; -} - -int @ModelName@_ewsb_solver::solve_ewsb_one_loop() +int @ModelName@_ewsb_solver::solve_one_loop() { - return solve_ewsb_iteratively(1); + return solve_iteratively(1); } -int @ModelName@_ewsb_solver::solve_ewsb() +int @ModelName@_ewsb_solver::solve() { - VERBOSE_MSG("\tSolving EWSB at " << ewsb_loop_order << "-loop order"); + VERBOSE_MSG("\tSolving EWSB at " << get_loop_order() << "-loop order"); - if (ewsb_loop_order == 0) - return solve_ewsb_tree_level(); + if (get_loop_order() == 0) + return solve_tree_level(); - return solve_ewsb_iteratively(ewsb_loop_order); + return solve_iteratively(get_loop_order()); } -void @ModelName@_ewsb_solver::ewsb_initial_guess(double x_init[number_of_ewsb_equations]) +void @ModelName@_ewsb_solver::initial_guess(double x_init[number_of_ewsb_equations]) { @ewsbInitialGuess@ } @@ -243,21 +206,21 @@ void @ModelName@_ewsb_solver::ewsb_initial_guess(double x_init[number_of_ewsb_eq /** * Calculates EWSB output parameters including loop-corrections. * - * @param ewsb_parameters new EWSB output parameters. \a - * ewsb_parameters is only modified if all new parameters are finite. + * @param parameters new EWSB output parameters. \a + * parameters is only modified if all new parameters are finite. * * @return GSL_SUCCESS if new EWSB output parameters are finite, * GSL_EDOM otherwise. */ -int @ModelName@_ewsb_solver::ewsb_step(double ewsb_parameters[number_of_ewsb_equations]) const +int @ModelName@_ewsb_solver::step(double parameters[number_of_ewsb_equations]) const { int error; double tadpole[number_of_ewsb_equations] = { 0. }; - if (ewsb_loop_order > 0) { -@calculateOneLoopTadpolesNoStruct@ - if (ewsb_loop_order > 1) { -@calculateTwoLoopTadpolesNoStruct@ + if (get_loop_order() > 0) { +@calculateOneLoopTadpolesFromModel@ + if (get_loop_order() > 1) { +@calculateTwoLoopTadpolesFromModel@ } } @@ -280,33 +243,34 @@ int @ModelName@_ewsb_solver::ewsb_step(double ewsb_parameters[number_of_ewsb_equ * @param params further function parameters * @param f new EWSB output parameters * - * @return Returns status of @ModelName@_ewsb_solver::ewsb_step + * @return Returns status of @ModelName@_ewsb_solver::step */ -int @ModelName@_ewsb_solver::ewsb_step(const gsl_vector* x, void* params, gsl_vector* f) +int @ModelName@_ewsb_solver::step(const gsl_vector* x, void* params, gsl_vector* f) { if (!is_finite(x)) { gsl_vector_set_all(f, std::numeric_limits::max()); return GSL_EDOM; } - const @ModelName@_ewsb_solver::EWSB_args* ewsb_args + const @ModelName@_ewsb_solver::EWSB_args* args = static_cast<@ModelName@_ewsb_solver::EWSB_args*>(params); - @ModelName@_mass_eigenstates* model = ewsb_args->model; - const unsigned ewsb_loop_order = ewsb_args->ewsb_loop_order; + @ModelName@_mass_eigenstates* model = args->model; + @ModelName@_ewsb_solver* solver = args->solver; + const unsigned loop_order = args->loop_order; @getEWSBParametersFromGSLVector@ @setEWSBParametersFromLocalCopies@ - if (ewsb_loop_order > 0) + if (loop_order > 0) model->calculate_DRbar_masses(); - double ewsb_parameters[number_of_ewsb_equations] = + double parameters[number_of_ewsb_equations] = @ewsbParametersInitializationList@; - const int status = model->ewsb_step(ewsb_parameters); + const int status = solver->step(parameters); for (std::size_t i = 0; i < number_of_ewsb_equations; ++i) - gsl_vector_set(f, i, ewsb_parameters[i]); + gsl_vector_set(f, i, parameters[i]); return status; } diff --git a/templates/ewsb_solver.hpp.in b/templates/ewsb_solver.hpp.in index 7c8f06a0d..497ff03ad 100644 --- a/templates/ewsb_solver.hpp.in +++ b/templates/ewsb_solver.hpp.in @@ -38,48 +38,65 @@ namespace flexiblesusy { class EWSB_solver; class @ModelName@_mass_eigenstates; -class @ModelName@_ewsb_solver { +class @ModelName@_ewsb_solver_interface { +public: + @ModelName@_ewsb_solver_interface(@ModelName@_mass_eigenstates* model_ = NULL) + : number_of_iterations(100) + , loop_order(2) + , precision(1.0e-5) + , model(model_) + {} + virtual ~@ModelName@_ewsb_solver_interface() {} + + void set_loop_order(unsigned n) { loop_order = n; } + void set_model(@ModelName@_mass_eigenstates* m) { model = m; } + void set_number_of_iterations(std::size_t n) { number_of_iterations = n; } + void set_precision(double p) { precision = p; } + + unsigned get_loop_order() const { return loop_order; } + @ModelName@_mass_eigenstates* get_model() const { return model; } + unsigned get_number_of_iterations() const { return number_of_iterations; } + double get_precision() const { return precision; } + + virtual int solve_tree_level() = 0; + virtual int solve_one_loop() = 0; + virtual int solve() = 0; + +private: + std::size_t number_of_iterations; + unsigned loop_order; + double precision; + @ModelName@_mass_eigenstates* model; +}; + +class @ModelName@_ewsb_solver : public @ModelName@_ewsb_solver_interface { public: /// number of EWSB equations static const std::size_t number_of_ewsb_equations = @numberOfEWSBEquations@; - @ModelName@_ewsb_solver(@ModelName@_mass_eigenstates*); + @ModelName@_ewsb_solver(@ModelName@_mass_eigenstates* model = NULL); ~@ModelName@_ewsb_solver(); - void set_ewsb_iteration_precision(double); - void set_ewsb_loop_order(unsigned); - void set_number_of_ewsb_iterations(std::size_t); - - double get_ewsb_iteration_precision() const; - double get_ewsb_loop_order() const; - - int solve_ewsb_tree_level(); - int solve_ewsb_one_loop(); - int solve_ewsb(); ///< solve EWSB at ewsb_loop_order level - - /// calculates the tadpoles at current loop order - void tadpole_equations(double[number_of_ewsb_equations]) const; + virtual int solve_tree_level(); + virtual int solve_one_loop(); + virtual int solve(); ///< solve EWSB at loop_order level private: struct EWSB_args { @ModelName@_mass_eigenstates* model; - unsigned ewsb_loop_order; + @ModelName@_ewsb_solver* solver; + unsigned loop_order; }; - std::size_t number_of_ewsb_iterations; - unsigned ewsb_loop_order; - double ewsb_iteration_precision; - - int solve_ewsb_iteratively(); - int solve_ewsb_iteratively(unsigned); - int solve_ewsb_iteratively_with(EWSB_solver*, const double[number_of_ewsb_equations]); - int solve_ewsb_tree_level_via_soft_higgs_masses(); - void ewsb_initial_guess(double[number_of_ewsb_equations]); - int ewsb_step(double[number_of_ewsb_equations]) const; - static int ewsb_step(const gsl_vector*, void*, gsl_vector*); + int solve_iteratively(); + int solve_iteratively(unsigned); + int solve_iteratively_with(EWSB_solver*, const double[number_of_ewsb_equations]); + int solve_tree_level_via_soft_higgs_masses(); + void initial_guess(double[number_of_ewsb_equations]); + int step(double[number_of_ewsb_equations]) const; + static int step(const gsl_vector*, void*, gsl_vector*); static int tadpole_equations(const gsl_vector*, void*, gsl_vector*); - }; } // namespace flexiblesusy diff --git a/templates/mass_eigenstates.cpp.in b/templates/mass_eigenstates.cpp.in index 184c83022..eea167530 100644 --- a/templates/mass_eigenstates.cpp.in +++ b/templates/mass_eigenstates.cpp.in @@ -31,6 +31,7 @@ */ #include "@ModelName@_mass_eigenstates.hpp" +#include "@ModelName@_ewsb_solver.hpp" #include "eigen_utils.hpp" #include "wrappers.hpp" #include "linalg2.hpp" @@ -84,14 +85,12 @@ using namespace @ModelName@_info; CLASSNAME::@ModelName@_mass_eigenstates(const @ModelName@_input_parameters& input_) : @ModelName@_soft_parameters(input_) - , number_of_ewsb_iterations(100) , number_of_mass_iterations(20) - , ewsb_loop_order(2) , pole_mass_loop_order(2) , calculate_sm_pole_masses(false) , force_output(false) , precision(1.0e-3) - , ewsb_iteration_precision(1.0e-5) + , ewsb_solver(new @ModelName@_ewsb_solver(this)) , physical() , problems(@ModelName@_info::particle_names) , two_loop_corrections() @@ -106,6 +105,7 @@ CLASSNAME::@ModelName@_mass_eigenstates(const @ModelName@_input_parameters& inpu CLASSNAME::~@ModelName@_mass_eigenstates() { + delete ewsb_solver; } void CLASSNAME::do_calculate_sm_pole_masses(bool flag) @@ -130,7 +130,7 @@ bool CLASSNAME::do_force_output() const void CLASSNAME::set_ewsb_loop_order(unsigned loop_order) { - ewsb_loop_order = loop_order; + ewsb_solver->set_loop_order(loop_order); } void CLASSNAME::set_two_loop_corrections(const Two_loop_corrections& two_loop_corrections_) @@ -140,7 +140,7 @@ void CLASSNAME::set_two_loop_corrections(const Two_loop_corrections& two_loop_co void CLASSNAME::set_number_of_ewsb_iterations(std::size_t iterations) { - number_of_ewsb_iterations = iterations; + ewsb_solver->set_number_of_iterations(iterations); } void CLASSNAME::set_number_of_mass_iterations(std::size_t iterations) @@ -151,7 +151,7 @@ void CLASSNAME::set_number_of_mass_iterations(std::size_t iterations) void CLASSNAME::set_precision(double precision_) { precision = precision_; - ewsb_iteration_precision = precision_; + ewsb_solver->set_precision(precision_); } void CLASSNAME::set_pole_mass_loop_order(unsigned loop_order) @@ -161,17 +161,17 @@ void CLASSNAME::set_pole_mass_loop_order(unsigned loop_order) void CLASSNAME::set_ewsb_iteration_precision(double precision) { - ewsb_iteration_precision = precision; + ewsb_solver->set_precision(precision); } double CLASSNAME::get_ewsb_iteration_precision() const { - return ewsb_iteration_precision; + return ewsb_solver->get_precision(); } double CLASSNAME::get_ewsb_loop_order() const { - return ewsb_loop_order; + return ewsb_solver->get_loop_order(); } const @ModelName@_physical& CLASSNAME::get_physical() const @@ -206,6 +206,8 @@ Problems<@ModelName@_info::NUMBER_OF_PARTICLES>& CLASSNAME::get_problems() */ void CLASSNAME::tadpole_equations(double tadpole[number_of_ewsb_equations]) const { + const unsigned ewsb_loop_order = ewsb_solver->get_loop_order(); + @calculateTreeLevelTadpoles@ if (ewsb_loop_order > 0) { @calculateOneLoopTadpoles@ @@ -215,137 +217,9 @@ void CLASSNAME::tadpole_equations(double tadpole[number_of_ewsb_equations]) cons } } -/** - * Method which calculates the tadpoles at loop order specified in the - * pointer to the CLASSNAME::EWSB_args struct. - * - * @param x GSL vector of EWSB output parameters - * @param params pointer to CLASSNAME::EWSB_args struct - * @param f GSL vector with tadpoles - * - * @return GSL_EDOM if x contains Nans, GSL_SUCCESS otherwise. - */ -int CLASSNAME::tadpole_equations(const gsl_vector* x, void* params, gsl_vector* f) -{ - if (!is_finite(x)) { - gsl_vector_set_all(f, std::numeric_limits::max()); - return GSL_EDOM; - } - - const CLASSNAME::EWSB_args* ewsb_args - = static_cast(params); - @ModelName@_mass_eigenstates* model = ewsb_args->model; - const unsigned ewsb_loop_order = ewsb_args->ewsb_loop_order; - -@setEWSBParametersFromGSLVector@ - - if (ewsb_loop_order > 0) - model->calculate_DRbar_masses(); - - double tadpole[number_of_ewsb_equations] = { 0. }; - - model->tadpole_equations(tadpole); - - for (std::size_t i = 0; i < number_of_ewsb_equations; ++i) - gsl_vector_set(f, i, tadpole[i]); - - return is_finite(tadpole) ? GSL_SUCCESS : GSL_EDOM; -} - -/** - * This method solves the EWSB conditions iteratively, trying several - * root finding methods until a solution is found. - */ -int CLASSNAME::solve_ewsb_iteratively() -{ - EWSB_args params = {this, ewsb_loop_order}; - - EWSB_solver* solvers[] = { -@EWSBSolvers@ - }; - - const std::size_t number_of_solvers = sizeof(solvers)/sizeof(*solvers); - double x_init[number_of_ewsb_equations] = { 0. }; - ewsb_initial_guess(x_init); - -#ifdef ENABLE_VERBOSE - std::cout << "Solving EWSB equations ...\n" - "\tInitial guess: x_init ="; - for (std::size_t i = 0; i < number_of_ewsb_equations; ++i) - std::cout << ' ' << x_init[i]; - std::cout << '\n'; -#endif - - int status; - for (std::size_t i = 0; i < number_of_solvers; ++i) { - VERBOSE_MSG("\tStarting EWSB iteration using solver " << i); - status = solve_ewsb_iteratively_with(solvers[i], x_init); - if (status == EWSB_solver::SUCCESS) { - VERBOSE_MSG("\tSolver " << i << " finished successfully!"); - break; - } -#ifdef ENABLE_VERBOSE - else { - WARNING("\tSolver " << i << " could not find a solution!" - " (requested precision: " << ewsb_iteration_precision << ")"); - } -#endif - } - - if (status == EWSB_solver::SUCCESS) { - problems.unflag_no_ewsb(); - } else { - problems.flag_no_ewsb(); -#ifdef ENABLE_VERBOSE - WARNING("\tCould not find a solution to the EWSB equations!" - " (requested precision: " << ewsb_iteration_precision << ")"); -#endif - } - - for_each(solvers, solvers + number_of_solvers, Delete_object()); - - return status; -} - -/** - * Solves EWSB equations with given EWSB solver - * - * @param solver EWSB solver - * @param x_init initial values - * - * @return status of the EWSB solver - */ -int CLASSNAME::solve_ewsb_iteratively_with( - EWSB_solver* solver, - const double x_init[number_of_ewsb_equations] -) -{ - const int status = solver->solve(x_init); - -@setEWSBSolution@ - - return status; -} - -int CLASSNAME::solve_ewsb_iteratively(unsigned loop_order) -{ - // temporarily set `ewsb_loop_order' to `loop_order' and do - // iteration - const unsigned old_loop_order = ewsb_loop_order; - ewsb_loop_order = loop_order; - const int status = solve_ewsb_iteratively(); - ewsb_loop_order = old_loop_order; - return status; -} - - int CLASSNAME::solve_ewsb_tree_level() { - int error = 0; - -@solveEwsbTreeLevel@ - - return error; + return ewsb_solver->solve_tree_level(); } int CLASSNAME::solve_ewsb_tree_level_via_soft_higgs_masses() @@ -359,93 +233,12 @@ int CLASSNAME::solve_ewsb_tree_level_via_soft_higgs_masses() int CLASSNAME::solve_ewsb_one_loop() { - return solve_ewsb_iteratively(1); + return ewsb_solver->solve_one_loop(); } int CLASSNAME::solve_ewsb() { - VERBOSE_MSG("\tSolving EWSB at " << ewsb_loop_order << "-loop order"); - - if (ewsb_loop_order == 0) - return solve_ewsb_tree_level(); - - return solve_ewsb_iteratively(ewsb_loop_order); -} - -void CLASSNAME::ewsb_initial_guess(double x_init[number_of_ewsb_equations]) -{ -@ewsbInitialGuess@ -} - -/** - * Calculates EWSB output parameters including loop-corrections. - * - * @param ewsb_parameters new EWSB output parameters. \a - * ewsb_parameters is only modified if all new parameters are finite. - * - * @return GSL_SUCCESS if new EWSB output parameters are finite, - * GSL_EDOM otherwise. - */ -int CLASSNAME::ewsb_step(double ewsb_parameters[number_of_ewsb_equations]) const -{ - int error; - double tadpole[number_of_ewsb_equations] = { 0. }; - - if (ewsb_loop_order > 0) { -@calculateOneLoopTadpolesNoStruct@ - if (ewsb_loop_order > 1) { -@calculateTwoLoopTadpolesNoStruct@ - } - } - -@solveEwsbWithTadpoles@ - - if (is_finite) { - error = GSL_SUCCESS; -@fillArrayWithEWSBParameters@ - } else { - error = GSL_EDOM; - } - - return error; -} - -/** - * Calculates EWSB output parameters including loop-corrections. - * - * @param x old EWSB output parameters - * @param params further function parameters - * @param f new EWSB output parameters - * - * @return Returns status of CLASSNAME::ewsb_step - */ -int CLASSNAME::ewsb_step(const gsl_vector* x, void* params, gsl_vector* f) -{ - if (!is_finite(x)) { - gsl_vector_set_all(f, std::numeric_limits::max()); - return GSL_EDOM; - } - - const CLASSNAME::EWSB_args* ewsb_args - = static_cast(params); - @ModelName@_mass_eigenstates* model = ewsb_args->model; - const unsigned ewsb_loop_order = ewsb_args->ewsb_loop_order; - -@getEWSBParametersFromGSLVector@ -@setEWSBParametersFromLocalCopies@ - - if (ewsb_loop_order > 0) - model->calculate_DRbar_masses(); - - double ewsb_parameters[number_of_ewsb_equations] = - @ewsbParametersInitializationList@; - - const int status = model->ewsb_step(ewsb_parameters); - - for (std::size_t i = 0; i < number_of_ewsb_equations; ++i) - gsl_vector_set(f, i, ewsb_parameters[i]); - - return status; + return ewsb_solver->solve(); } void CLASSNAME::print(std::ostream& ostr) const diff --git a/templates/mass_eigenstates.hpp.in b/templates/mass_eigenstates.hpp.in index 635654de1..37edd0e1f 100644 --- a/templates/mass_eigenstates.hpp.in +++ b/templates/mass_eigenstates.hpp.in @@ -52,6 +52,7 @@ namespace flexiblesusy { class EWSB_solver; +class @ModelName@_ewsb_solver_interface; /** * @class @ModelName@_mass_eigenstates * @brief model class with routines for determing masses and mixinga and EWSB @@ -143,14 +144,12 @@ private: }; #endif - std::size_t number_of_ewsb_iterations; std::size_t number_of_mass_iterations; - unsigned ewsb_loop_order; unsigned pole_mass_loop_order; bool calculate_sm_pole_masses; ///< switch to calculate the pole masses of the Standard Model particles bool force_output; ///< switch to force output of pole masses double precision; ///< RG running precision - double ewsb_iteration_precision; + @ModelName@_ewsb_solver_interface* ewsb_solver; @ModelName@_physical physical; ///< contains the pole masses and mixings Problems<@ModelName@_info::NUMBER_OF_PARTICLES> problems; Two_loop_corrections two_loop_corrections; ///< used 2-loop corrections @@ -159,14 +158,7 @@ private: static std::mutex mtx_fortran; /// locks fortran functions #endif - int solve_ewsb_iteratively(); - int solve_ewsb_iteratively(unsigned); - int solve_ewsb_iteratively_with(EWSB_solver*, const double[number_of_ewsb_equations]); int solve_ewsb_tree_level_via_soft_higgs_masses(); - void ewsb_initial_guess(double[number_of_ewsb_equations]); - int ewsb_step(double[number_of_ewsb_equations]) const; - static int ewsb_step(const gsl_vector*, void*, gsl_vector*); - static int tadpole_equations(const gsl_vector*, void*, gsl_vector*); void copy_DRbar_masses_to_pole_masses(); // Passarino-Veltman loop functions diff --git a/templates/module.mk.in b/templates/module.mk.in index 826f6e986..3514a5dff 100644 --- a/templates/module.mk.in +++ b/templates/module.mk.in @@ -34,6 +34,7 @@ LIB@CLASSNAME@_HDR := ifneq ($(findstring two_scale,$(ALGORITHMS)),) LIB@CLASSNAME@_SRC += \ + $(DIR)/@CLASSNAME@_ewsb_solver.cpp \ $(DIR)/@CLASSNAME@_mass_eigenstates.cpp \ $(DIR)/@CLASSNAME@_info.cpp \ $(DIR)/@CLASSNAME@_input_parameters.cpp \ @@ -55,6 +56,7 @@ EXE@CLASSNAME@_SRC += \ $(DIR)/scan_@CLASSNAME@.cpp LIB@CLASSNAME@_HDR += \ $(DIR)/@CLASSNAME@_convergence_tester.hpp \ + $(DIR)/@CLASSNAME@_ewsb_solver.hpp \ $(DIR)/@CLASSNAME@_high_scale_constraint.hpp \ $(DIR)/@CLASSNAME@_mass_eigenstates.hpp \ $(DIR)/@CLASSNAME@_info.hpp \