Merge branch 'release-v0.6.0'

CMakeLists.txt

@@ -4,7 +4,7 @@ if (POLICY CMP0042)
 endif ()
 PROJECT(EVOEMD CXX)
 
-SET(EVOEMD_VERSION "0.5.0")
+SET(EVOEMD_VERSION "0.6.0")
 
 FIND_PACKAGE(GSL 2.1 REQUIRED)
 FIND_PACKAGE(Eigen3 REQUIRED)
@@ -31,4 +31,5 @@ TARGET_LINK_LIBRARIES(EVOEMD_shared ${GSL_LIBRARIES} cuba)
 SET(EVOEMD_LIBRARY EVOEMD_shared)
 
 add_subdirectory(src/Models/SeesawFreezeIn)
+add_subdirectory(src/Models/ToyDM)
 # add_subdirectory(test)

README.md

@@ -2,10 +2,18 @@
 
 A framework to evaluate the evolution w/ or w/o an early matter dominated (EMD) era.
 
-## Feature
+## Features
 
 - Global parameter system
 - Global particle system
 - Global process system
-- Automatically build up the Boltzmann equation according to the user's definition of particle and process (working in progress)
-- Solving the Boltzmann equation with RK to 5th order
+- Automatically build up the Boltzmann equation according to the user's definition of particle and process
+- Solving the Boltzmann equation using 4th order Runge-Kutta method with adaptive steps tailored to Cosmology application
+
+## Document
+
+In preparation...
+
+## Examples
+
+Simple examples can be found in corresponding folders under `Models`

include/EvoEMD/BoltzmannEquation.h

@@ -0,0 +1,34 @@
+#ifndef _BOLTZMANN_EQUATION_H_
+#define _BOLTZMANN_EQUATION_H_
+
+#include <string>
+#include <vector>
+
+#include "EvoEMD/Hubble.h"
+#include "EvoEMD/ParticleBase.h"
+#include "EvoEMD/RungeKutta.h"
+
+namespace EvoEMD {
+
+class BoltzmannEquation : public ODE_FUNCS {
+private:
+    Parameter_Base *ptr_scale;
+    REAL scale;
+    std::vector<int> poi_pids;
+    std::vector<std::string> poi_names;
+    std::vector<Pseudo_Particle *> poi_ptrs;
+    Particle_Factory &pf;
+    Hubble_History &hh;
+    void Setup_Scale();
+
+public:
+    BoltzmannEquation(Parameter_Base *scale = nullptr);
+    ~BoltzmannEquation(){};
+
+    virtual VD dYdX(REAL x, VD y) override;
+    virtual VD Yeq(REAL x) override;
+};
+
+}  // namespace EvoEMD
+
+#endif  //_BOLTZMANN_EQUATION_H_

include/EvoEMD/CollisionRate.h

@@ -11,7 +11,6 @@ struct Process_Amp {
     // * used to store how many squared diagrams we have for corresponding amplitude calculation
     unsigned n_diag;
 
-    // * bool determine whether the corresponding REAL is actually exactly zero
     typedef std::vector<REAL> CTH_RES_FULL;
     typedef CTH_RES_FULL NUMERATOR_STRUCTURE;
     typedef int Propagator_ID;
@@ -40,14 +39,14 @@ class Amplitude {
     typedef std::vector<Pseudo_Particle *> INITIAL_STATES;
     typedef std::vector<Pseudo_Particle *> FINAL_STATES;
 
-    Amplitude();
+    Amplitude(){};
     virtual ~Amplitude(){};
     virtual void Update_Amp(REAL sqrt_shat) = 0;
     virtual const Process_Amp &Get_Amp(REAL sqrt_shat) {
         Update_Amp(sqrt_shat);
         return amp_res;
     };
-    virtual REAL Get_Coeff(REAL T) = 0;
+    virtual REAL Get_Coeff(REAL T, int PID) = 0;
 
     unsigned N_INITIAL;
     INITIAL_STATES INITIAL;

include/EvoEMD/Constants.h

@@ -5,5 +5,6 @@
 #define GeV 1
 #define TeV 1000
 #define eV 1e-9
+#define BESSEL_Z_MAX 200
 
 #endif  //_CONSTANTS_H_

include/EvoEMD/Hubble.h

@@ -1,6 +1,7 @@
 #ifndef _HUBBLE_H_
 #define _HUBBLE_H_
 
+#include "EvoEMD/ParameterBase.h"
 #include "EvoEMD/RealTypes.h"
 
 namespace EvoEMD {
@@ -41,6 +42,7 @@ class Hubble_For_Single_Period {
     double Get_beta_T() const { return beta_T; }
     double Get_beta_s() const { return beta_s; }
     bool Is_Isentropic() const { return Isentropic; }
+    void Print() const;
 };
 
 class Hubble_RD : public Hubble_For_Single_Period {
@@ -74,7 +76,7 @@ class Hubble_EP : public Hubble_For_Single_Period {
     virtual REAL Get_Hubble_at_T(const REAL T) override;
 };
 
-class Hubble_History {
+class Hubble_History : public Parameter_Base {
 private:
     std::vector<Hubble_For_Single_Period *> Periods;
     std::vector<REAL> Temperatures;
@@ -86,17 +88,22 @@ class Hubble_History {
 
     void Solve_Te();
 
-public:
-    Hubble_History(const REAL Ti, const REAL Tr);
+    Hubble_History();
     Hubble_History(const Hubble_History &HH);
-    ~Hubble_History();
-
     Hubble_History &operator=(const Hubble_History &HH);
+    ~Hubble_History();
 
+public:
+    static Hubble_History &Get_Hubble_History();
     int Get_N_Period() const { return Periods.size(); }
-    int Get_Period_ID_at_T(const REAL T) const;
+    int Get_Period_ID_at_T(const REAL T);
     REAL Get_Hubble_at_T(const REAL T);
+    double Get_beta_T_at_T(const REAL T);
     Hubble_For_Single_Period *operator[](const int pid);
+    Hubble_For_Single_Period *at(const int pid);
+
+    virtual void Update_Value(REAL input) override;
+    void Print_History();
 };
 
 }  // namespace EvoEMD

include/EvoEMD/ParameterBase.h

@@ -24,6 +24,7 @@ class Parameter_Base {
     void Claim_Dependence(Parameter_Base *par) { underlying_parameters.insert(par); }
     void Notify() {
         updated = false;
+        // std::cout << "Parameter: " << name << " Need to be Updated" << std::endl;
         for (auto &&bp : underlying_parameters) {
             bp->Notify();
         }
@@ -36,6 +37,7 @@ class Parameter_Base {
         Notify();
         Update_Value(input);
         updated = true;
+        // std::cout << "Parameter: " << name << " has been updated" << std::endl;
     }
     REAL Get_Value() {
         if (updated) {
@@ -84,7 +86,7 @@ class Register_Parameter {
     };
 };
 
-#define REGISTER_PARAMETER(paramName) Register_Parameter g_register_parameter_##paramName(new paramName)
+#define REGISTER_PARAMETER(paramName) const Register_Parameter g_register_parameter_##paramName(new paramName)
 #define DECLARE_FREE_PARAMETER(paramName, value)                   \
     class param_##paramName : public Free_Parameter {              \
     public:                                                        \

include/EvoEMD/ParticleBase.h

@@ -37,6 +37,7 @@ class Particle_Base {
     void Register_Process(Process *);
 
     void Notify_Client();
+    std::set<Process *> Get_Process() const { return Process_Set; }
 };
 
 class Pseudo_Particle : public Particle_Base {
@@ -47,6 +48,7 @@ class Pseudo_Particle : public Particle_Base {
     // * DOF: degree of freedom, (particle and antiparticle count seperately)
     // * Calculate Number Density or Yield at Equilibrium;
 protected:
+    std::string name;
     bool selfconjugate;
     bool massless;
     bool thermalized;
@@ -72,13 +74,20 @@ class Pseudo_Particle : public Particle_Base {
      * @param  selfconjugate: Whether
      * @retval
      */
-    Pseudo_Particle(int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
+    Pseudo_Particle(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
                     bool selfconjugate = false);
     virtual ~Pseudo_Particle(){};
 
     int Get_PID() const { return PID; }
     int Get_DOF() const { return DOF; }
-    double Get_Mass() const { return p_mass->Get_Value(); }
+    std::string Get_Name() const { return name; }
+    double Get_Mass() const {
+        if (massless) {
+            return 0;
+        } else {
+            return p_mass->Get_Value();
+        }
+    }
     bool Is_Massless() const { return massless; }
     bool Is_Selfconjugate() const { return selfconjugate; }
 
@@ -87,8 +96,10 @@ class Pseudo_Particle : public Particle_Base {
     }
     REAL Get_Equilibrium_Yield_at_T(const REAL T) const { return DOF * Get_Equilibrium_Yield_per_DOF(T); };
 
+    bool start_with_thermal;
     REAL Yield;
     REAL Numer_Density;
+
     void Set_Mass(double mass);
 };
 
@@ -98,7 +109,7 @@ class Fermion : public Pseudo_Particle {
     virtual REAL Get_Equilibrium_Yield_per_DOF(const REAL T) const override;
 
 public:
-    Fermion(int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
+    Fermion(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
             bool selfconjugate = false);
     ~Fermion(){};
 };
@@ -109,7 +120,7 @@ class Boson : public Pseudo_Particle {
     virtual REAL Get_Equilibrium_Yield_per_DOF(const REAL T) const override;
 
 public:
-    Boson(int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
+    Boson(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
           bool selfconjugate = false);
     ~Boson(){};
 };
@@ -124,6 +135,7 @@ class Particle_Factory {
     bool Register_Particle(Pseudo_Particle *);
     bool Register_POI(int PID);
     bool Set_Mass(int PID, double mass);
+    std::set<int> Get_POI() const { return POI; }
 
 private:
     Particle_Factory();
@@ -144,11 +156,23 @@ class Register_Particle {
 
 // #define REGISTER_PARTICLE(partName) Register_Particle g_register_particle_##partName(new partName)
 
-#define REGISTER_PARTICLE(className, partName, PID, DOF, MASS, WIDTH, C) \
-    class part_##partName : public className {                           \
-    public:                                                              \
-        part_##partName() : className(PID, DOF, MASS, WIDTH, C){};       \
-    };                                                                   \
+#define REGISTER_PARTICLE(className, partName, PID, DOF, MASS, WIDTH, C)      \
+    class part_##partName : public className {                                \
+    public:                                                                   \
+        part_##partName() : className(#partName, PID, DOF, MASS, WIDTH, C){}; \
+    };                                                                        \
     Register_Particle g_register_particle_##partName(new part_##partName)
 
+class Register_POI {
+public:
+    Register_POI(int PID, bool start_with_thermal = true) {
+        EvoEMD::Particle_Factory &pf = EvoEMD::Particle_Factory::Get_Particle_Factory();
+        pf.Register_POI(PID);
+        EvoEMD::Pseudo_Particle *pp = pf.Get_Particle(PID);
+        pp->start_with_thermal = start_with_thermal;
+    }
+};
+
+#define REGISTER_POI(PID, THERMAL) Register_POI g_register_poi_##PID(PID, THERMAL)
+
 #endif  //_PARTICLE_H_

include/EvoEMD/ProcessBase.h

@@ -13,7 +13,7 @@ class Process {
     ~Process();
 
     REAL Get_Collision_Rate(REAL T);
-    REAL Get_Yield_Coeff(REAL T);
+    REAL Get_Yield_Coeff(REAL T, int PID);
 
 protected:
     INITIAL_STATES INIT;

include/EvoEMD/RungeKutta.h

@@ -32,6 +32,7 @@ class ODE_FUNCS {
     VD Get_BOUNDARY_CONDITION() const { return BOUNDARY_CONDITION; }
     VD operator()(REAL x, VD y) { return dYdX(x, y); }
     virtual VD dYdX(REAL x, VD y) = 0;
+    virtual VD Yeq(REAL x) = 0;
 };
 
 /*
@@ -83,6 +84,7 @@ class RungeKutta {
 
     VD _X;      // * The vector storing the points in x
     VVD _Y;     // * The vector storing the points in y, len(_Y) = len(_X) and the second dimension is DOF
+    VVD _Yeq;   // * The vector storing the equilibrium value of Y;
     VVD _dYdX;  // * Similar to _Y, but storing dy/dx
 
     VD BOUNDARY_AT_BEGIN;  // * The boundary condition at X_BEGIN
@@ -110,7 +112,7 @@ class RungeKutta {
      * @param step_size, step size, we would like to go forward
      * @param dy_next, the dy at next step (output)
      */
-    void RK4_SingleStep(const REAL x_cur, const VD &y_cur, const VD &dy_cur, const REAL step_size, VD &y_next);
+    bool RK4_SingleStep(const REAL x_cur, const VD &y_cur, const VD &dy_cur, const REAL step_size, VD &y_next);
 
     /**
      * @brief The Runge-Kutta method taking one step forward
@@ -126,7 +128,7 @@ class RungeKutta {
      * @param step_size_did, output: actual step size we take
      * @param step_size_further, output: based on what we did, the prospect step size for next step
      */
-    void RKQC_SingleStep(REAL &x, VD &y, VD &dy, const REAL step_size_guess, const REAL eps, const VD &Y_Scale,
+    bool RKQC_SingleStep(REAL &x, VD &y, VD &dy, const REAL step_size_guess, const REAL eps, const VD &Y_Scale,
                          REAL &step_size_did, REAL &step_size_further);
 };
 }  // namespace EvoEMD

include/Models/SeesawFreezeIn/Particles.h

@@ -6,7 +6,10 @@
 using namespace EvoEMD;
 
 REGISTER_PARTICLE(Fermion, N1, 900001, 2, RETRIVE_PARAMETER(MN1), nullptr, true);
+REGISTER_POI(900001, 1);
 
 REGISTER_PARTICLE(Fermion, N2, 900002, 2, RETRIVE_PARAMETER(MN2), nullptr, true);
 
+REGISTER_PARTICLE(Fermion, N3, 900003, 2, RETRIVE_PARAMETER(MN3), nullptr, true);
+
 #endif  //_SEESAW_FI_PARTICLES_H_

include/Models/SeesawFreezeIn/SeesawTypeI.h

@@ -34,7 +34,7 @@ class SeesawTypeI : public EvoEMD::Parameter_Base {
 
     std::complex<double> Get_Yij(int i, int j);
     std::complex<double> Get_YdagYij(int i, int j);
-    virtual void Update_Value(REAL input);
+    virtual void Update_Value(REAL input) override;
 };
 
 #endif  //_SEESAW_TYPE_I_H_

include/Models/ToyDM/Amplitudes.h

@@ -0,0 +1,15 @@
+#ifndef _TOY_DM_PROCESSES_H
+#define _TOY_DM_PROCESSES_H
+
+#include "EvoEMD/CollisionRate.h"
+
+class XX_SS_Amp : public EvoEMD::Amplitude {
+public:
+    XX_SS_Amp();
+    ~XX_SS_Amp(){};
+
+    virtual void Update_Amp(REAL sqrt_shat) override;
+    virtual REAL Get_Coeff(REAL T, int PID) override;
+};
+
+#endif  //_TOY_DM_PROCESSES_H

include/Models/ToyDM/Parameters.h

@@ -0,0 +1,9 @@
+#ifndef _TOY_DM_PARAMETERS_H_
+#define _TOY_DM_PARAMETERS_H_
+
+#include "EvoEMD/ParameterBase.h"
+using namespace EvoEMD;
+DECLARE_FREE_PARAMETER(Lam, 1e-2);
+DECLARE_FREE_PARAMETER(MX, 100);
+
+#endif  //_TOY_DM_PARAMETERS_H_

include/Models/ToyDM/Particles.h

@@ -0,0 +1,12 @@
+#ifndef _TOY_DM_PARTICLES_H_
+#define _TOY_DM_PARTICLES_H_
+
+#include "EvoEMD/ParticleBase.h"
+#include "Parameters.h"
+using namespace EvoEMD;
+REGISTER_PARTICLE(Boson, X, 900001, 2, RETRIVE_PARAMETER(MX), nullptr, true);
+REGISTER_POI(900001, 1);
+
+REGISTER_PARTICLE(Boson, S, 900011, 1, nullptr, nullptr, true);
+
+#endif  //_TOY_DM_PARTICLES_H_

include/Models/ToyDM/Processes.h

@@ -0,0 +1,4 @@
+#include "Amplitudes.h"
+#include "EvoEMD/ProcessBase.h"
+
+EvoEMD::Process *g_ss_xx = new EvoEMD::Process(new XX_SS_Amp);