Merge branch 'release-v0.6.1'

CMakeLists.txt

@@ -4,7 +4,7 @@ if (POLICY CMP0042)
 endif ()
 PROJECT(EVOEMD CXX)
 
-SET(EVOEMD_VERSION "0.6.0")
+SET(EVOEMD_VERSION "0.6.1")
 
 FIND_PACKAGE(GSL 2.1 REQUIRED)
 FIND_PACKAGE(Eigen3 REQUIRED)

README.md

@@ -14,6 +14,12 @@ A framework to evaluate the evolution w/ or w/o an early matter dominated (EMD)
 
 In preparation...
 
+## Dependencies
+
+- GSL
+- Eigen3
+- Cuba
+
 ## Examples
 
 Simple examples can be found in corresponding folders under `Models`

include/EvoEMD/BoltzmannEquation.h

@@ -25,8 +25,12 @@ class BoltzmannEquation : public ODE_FUNCS {
     BoltzmannEquation(Parameter_Base *scale = nullptr);
     ~BoltzmannEquation(){};
 
-    virtual VD dYdX(REAL x, VD y) override;
+    void Set_X_Range(REAL X_BEGIN, REAL X_END);
+    virtual VD dYdX(REAL x, VD y, VD delta_y_ratio) override;
     virtual VD Yeq(REAL x) override;
+    virtual VB Is_Thermalized() override;
+    virtual void Update_Thermal_Status(VB status) override;
+    virtual VB Can_be_Thermalized() override;
 };
 
 }  // namespace EvoEMD

include/EvoEMD/EvoEMD.h

@@ -0,0 +1,13 @@
+#ifndef _EVOEMD_H_
+#define _EVOEMD_H_
+
+#include "EvoEMD/BoltzmannEquation.h"
+#include "EvoEMD/Constants.h"
+#include "EvoEMD/Hubble.h"
+#include "EvoEMD/ParameterBase.h"
+#include "EvoEMD/ParticleBase.h"
+#include "EvoEMD/ProcessBase.h"
+#include "EvoEMD/RealTypes.h"
+#include "EvoEMD/spdlog_wrapper.h"
+
+#endif  //_EVOEMD_H_

include/EvoEMD/ParameterBase.h

@@ -86,14 +86,16 @@ class Register_Parameter {
     };
 };
 
-#define REGISTER_PARAMETER(paramName) const Register_Parameter g_register_parameter_##paramName(new paramName)
+#define REGISTER_PARAMETER(paramName)
 #define DECLARE_FREE_PARAMETER(paramName, value)                   \
     class param_##paramName : public Free_Parameter {              \
     public:                                                        \
         param_##paramName() : Free_Parameter(#paramName, value){}; \
     };                                                             \
-    REGISTER_PARAMETER(param_##paramName)
+    const Register_Parameter g_register_parameter_##paramName(new param_##paramName)
 
 #define RETRIVE_PARAMETER(paramName) EvoEMD::Parameter_Factory::Get_Parameter_Factory().Get_Parameter(#paramName)
 
+#define GET_PARAM_VALUE(paramName) RETRIVE_PARAMETER(paramName)->Get_Value()
+
 #endif  //_PARAMETER_BASE_H_

include/EvoEMD/ParticleBase.h

@@ -10,20 +10,9 @@
 namespace EvoEMD {
 
 class Process;
-class Particle_Client {
-public:
-    Particle_Client(){};
-    virtual ~Particle_Client(){};
-
-    virtual void Update_Particle_Info() = 0;
-};
 
 class Particle_Base {
 protected:
-    // * Objects that use particle information
-    // * In any time, particle information is updated, the client should update their own properties
-    std::set<Particle_Client *> Particle_Client_Set;
-
     // * Processes that involve current particle
     // * The process calculates its value lazily, it will acquire particle's infor when it needs
     // * This set is used to keep aware what are the processes involving current particle
@@ -33,7 +22,6 @@ class Particle_Base {
     Particle_Base(){};
     virtual ~Particle_Base(){};
 
-    void Register_Client(Particle_Client *);
     void Register_Process(Process *);
 
     void Notify_Client();
@@ -49,9 +37,7 @@ class Pseudo_Particle : public Particle_Base {
     // * Calculate Number Density or Yield at Equilibrium;
 protected:
     std::string name;
-    bool selfconjugate;
     bool massless;
-    bool thermalized;
     Parameter_Base *p_mass;
     Parameter_Base *p_width;
     int PID;
@@ -63,19 +49,19 @@ class Pseudo_Particle : public Particle_Base {
 
 public:
     /**
-     * @brief
+     * @brief  ctor for a Pseudo_Particle
      * @note
-     * @param  PID: The PID for a particle
+     * @param  name: The name for the particle
+     * @param  PID: The PID for a particle (should be unique for each particle)
      * @param  DOF: The degree of freedom of the particle, particle and antiparticle will be counted seperately
      *              On the other hand, DOF is the one used to calculate the equilibrium number density which will be
      *              used in Boltzmann equation. So be sure this DOF is consistent with the collision rate.
      * @param  mass: Pointer to mass parameter, if nullptr (default), it is assumed the particle is massless
      * @param  width: Pointer to width parameter, if nullptr (default), it is assumed the particle is stable
-     * @param  selfconjugate: Whether
      * @retval
      */
     Pseudo_Particle(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
-                    bool selfconjugate = false);
+                    bool never_thermal = false);
     virtual ~Pseudo_Particle(){};
 
     int Get_PID() const { return PID; }
@@ -89,15 +75,24 @@ class Pseudo_Particle : public Particle_Base {
         }
     }
     bool Is_Massless() const { return massless; }
-    bool Is_Selfconjugate() const { return selfconjugate; }
 
     REAL Get_Equilibrium_Number_Density_at_T(const REAL T) const {
         return DOF * Get_Equilibrium_Number_Density_per_DOF(T);
     }
     REAL Get_Equilibrium_Yield_at_T(const REAL T) const { return DOF * Get_Equilibrium_Yield_per_DOF(T); };
 
-    bool start_with_thermal;
+    const bool Never_Thermal;
+    bool Thermalized;
     REAL Yield;
+
+    // * 1 - Yield/Yield_eq
+    // * At sufficient high temperature, particle might be in thermal equilibrium
+    // * In that case, this value is actually 0. But due to numerical issue, it can be a small number
+    // * However, in calculating collision rate, it might be multiplied by a extremely large number,
+    // * thus leads to numerical issue.
+    // * So we will keep this number, especially when it is zero to avoid numerical issue.
+    // * User can set their own threshold when to use this Delta_Yield_Ratio
+    REAL Delta_Yield_Ratio;
     REAL Numer_Density;
 
     void Set_Mass(double mass);
@@ -110,7 +105,7 @@ class Fermion : public Pseudo_Particle {
 
 public:
     Fermion(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
-            bool selfconjugate = false);
+            bool never_thermal = false);
     ~Fermion(){};
 };
 
@@ -121,7 +116,7 @@ class Boson : public Pseudo_Particle {
 
 public:
     Boson(std::string name, int PID, int DOF, Parameter_Base *mass = nullptr, Parameter_Base *width = nullptr,
-          bool selfconjugate = false);
+          bool never_thermal = false);
     ~Boson(){};
 };
 
@@ -156,20 +151,22 @@ 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(#partName, PID, DOF, MASS, WIDTH, C){}; \
-    };                                                                        \
+#define REGISTER_PARTICLE(className, partName, ...)              \
+    class part_##partName : public className {                   \
+    public:                                                      \
+        part_##partName() : className(#partName, __VA_ARGS__){}; \
+    };                                                           \
     Register_Particle g_register_particle_##partName(new part_##partName)
 
+#define RETRIVE_PARTICLE(PID) EvoEMD::Particle_Factory::Get_Particle_Factory().Get_Particle(PID)
+
 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;
+        pp->Thermalized = start_with_thermal;
     }
 };
 

include/EvoEMD/ProcessBase.h

@@ -1,6 +1,8 @@
 #ifndef _PROCESS_BASE_H_
 #define _PROCESS_BASE_H_
 
+#include <set>
+
 #include "EvoEMD/CollisionRate.h"
 
 namespace EvoEMD {
@@ -21,5 +23,36 @@ class Process {
     Amplitude *amp;
     Collision_Rate *CR_Calculator;
 };
+
+class Process_Factory {
+    // * This class is used to manage all process, it is a singleton
+    // * The only object of this class will own all processes
+    // * Users won't directly use this class unless one wants to dig into the process
+
+public:
+    typedef std::set<Process *> Process_List;
+    static Process_Factory &Get_Process_Factory();
+    void Register_Process(Process *proc) { PL.insert(proc); };
+    Process_List Get_Process_List() { return PL; }
+
+private:
+    Process_List PL;
+
+    Process_Factory(){};
+    ~Process_Factory();
+};
+
 }  // namespace EvoEMD
+
+class Register_Process {
+public:
+    Register_Process(EvoEMD::Process *proc) {
+        EvoEMD::Process_Factory &pf = EvoEMD::Process_Factory::Get_Process_Factory();
+        pf.Register_Process(proc);
+    }
+};
+
+#define REGISTER_PROCESS(AmpClassName) \
+    Register_Process g_register_process_##AmpClassName(new EvoEMD::Process(new AmpClassName))
+
 #endif  //_PROCESS_BASE_H_

include/EvoEMD/RealTypes.h

@@ -4,6 +4,7 @@
 #include <iostream>
 #include <vector>
 
+#define VB std::vector<bool>
 #define REAL double
 #define VD std::vector<REAL>
 #define VVD std::vector<VD>
@@ -25,7 +26,7 @@ VD operator-(const VD &rhs);
 // * Multiply Operator
 VD operator*(const VD &lhs, const REAL &s);
 VD operator*(const REAL &s, const VD &rhs);
-REAL operator*(const VD &lhs, const VD &rhs);  // Scalar Product
+// REAL operator*(const VD &lhs, const VD &rhs);  // Scalar Product
 
 // * Divide Operator
 VD operator/(const VD &lhs, const VD &rhs);  // elementary-wise divide

include/EvoEMD/RungeKutta.h

@@ -15,24 +15,31 @@ class ODE_FUNCS {
     int DOF;
     REAL X_BEGIN;
     REAL X_END;
-    VD BOUNDARY_CONDITION;
+    VD Y_BEGIN;
+    VD Delta_Y_Ratio_BEGIN;
 
 public:
-    ODE_FUNCS(){};
+    ODE_FUNCS() : DOF(-1){};
     virtual ~ODE_FUNCS(){};
 
-    void Set_DOF(int dof) { DOF = dof; }
-    void Set_X_BEGIN(REAL xi) { X_BEGIN = xi; }
-    void Set_X_END(REAL xf) { X_END = xf; }
-    void Set_BOUNDARY_CONDITION(VD boundary) { BOUNDARY_CONDITION = boundary; }
+    void Set_DOF(int DOF) {
+        this->DOF = DOF;
+        Y_BEGIN = VD(DOF, 0);
+        Delta_Y_Ratio_BEGIN = VD(DOF, 0);
+    }
+    void Set_BOUNDARIES(REAL X_BEGIN, REAL X_END, VD boundary, bool relative_to_equilibrium = false);
 
     int Get_DOF() const { return DOF; }
     REAL Get_X_BEGIN() const { return X_BEGIN; }
     REAL Get_X_END() const { return X_END; }
-    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;
+    VD Get_Y_BEGIN() const { return Y_BEGIN; }
+    VD Get_Delta_Y_Ratio_BEGIN() const { return Delta_Y_Ratio_BEGIN; }
+    VD operator()(REAL x, VD y, VD delta_y_ratio) { return dYdX(x, y, delta_y_ratio); }
+    virtual VD dYdX(REAL x, VD y, VD delta_y_ratio) = 0;  // delta_y_ratio = 1 - Y/Yeq;
     virtual VD Yeq(REAL x) = 0;
+    virtual VB Is_Thermalized() = 0;  // Check whether the particle is starting in thermalization or not
+    virtual void Update_Thermal_Status(VB status) = 0;
+    virtual VB Can_be_Thermalized() = 0;
 };
 
 /*
@@ -77,18 +84,20 @@ class RungeKutta {
     REAL TINY_REL_THRESHOLD;
     int MAXSTEPS;
 
-    bool SOLVED;   // * Store the status whether the ODE functions is solved
-    int DOF;       // * The degree of freedon of the ODE functions
-    REAL X_BEGIN;  // * The starting point in x (the argument)
-    REAL X_END;    // * The ending point in x (the argument)
+    bool SOLVED;             // * Store the status whether the ODE functions is solved
+    int DOF;                 // * The degree of freedon of the ODE functions
+    REAL X_BEGIN;            // * The starting point in x (the argument)
+    REAL X_END;              // * The ending point in x (the argument)
+    VD Y_BEGIN;              // * The starting point of y
+    VD Delta_Y_Ratio_BEGIN;  // * 1-y/yeq at the begin;
 
-    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 _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 _Delta_Y_Ratio;  // * The vector storing the value 1 - y/yeq
+    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
-    VD ZERO_THRESHOLD;     // * The threshold to treat the corresponding value as 0;
+    VD ZERO_THRESHOLD;  // * The threshold to treat the corresponding value as 0;
 
     ODE_FUNCS *
         derivs;  // * The ODE function, dy/dx = derivs(x,y,xxxx), we don't own it, we just have one pointer pointing it.
@@ -108,28 +117,38 @@ class RungeKutta {
      * @brief The 4th order Runge-Kutta method taking one step forward
      * @param x_cur, current x value
      * @param y_cur, current y value
-     * @param dy_cur, current dy/dx value
+     * @param dydx_cur, current dy/dx value
+     * @param delta_y_ratiro_cur, current 1-Y/Yeq
+     * @param thermal_status_cur, current thermalization status
      * @param step_size, step size, we would like to go forward
-     * @param dy_next, the dy at next step (output)
+     * @param y_next, the dy at next step (output)
+     * @param delta_y_ratiro_next, 1-Y/Yeq at next step (output)
+     * @param thermal_status_next, thermalization status at next step (output)
      */
-    bool 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 &dydx_cur, const VD &delta_y_ratio_cur,
+                        const VB &thermal_status_cur, const REAL step_size, VD &y_next, VD &delta_y_ratio_next,
+                        VB &thermal_status_next);
 
     /**
      * @brief The Runge-Kutta method taking one step forward
      *        Based on 4th order Runge-Kutta and adding quality control (adaptive stepsize),
      *        such that we can kind of achieve 5th order accuracy
      * @param x, input: current value of x, output: next value of x
      * @param y, input: current value of y, output: next value of y
-     * @param dy, input: current value of dy, output: next value of dy
+     * @param yeq, input: current value of yeq, output: next value of yeq
+     * @param dydx, input: current value of dy/dx, output: next value of dy/dx
+     * @param delta_y_ratio, input: current value of 1-y/yeq, output: next value of 1-y/yeq
+     * @param thermal_status, input: current status of thermalization, output: next status of thermalization
      * @param step_size_guess, input: initial guess of current step size
      * @param eps, input: tolerance
      * @param Y_Scale, input: the scale in y to set the error (it is possible in different direction of y, the scale is
      * quite different)
      * @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
      */
-    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);
+    bool RKQC_SingleStep(REAL &x, VD &y, VD &yeq, VD &dydx, VD &delta_y_ratio, VB &thermal_status,
+                         const REAL step_size_guess, const REAL eps, const VD &Y_Scale, REAL &step_size_did,
+                         REAL &step_size_further);
 };
 }  // namespace EvoEMD
 #endif  //__RUNGEKUTTA_H__