ThermalModelVDW.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2016-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALMODELVDW_H
#define THERMALMODELVDW_H
 
#include "HRGBase/ThermalModelBase.h"
 
namespace thermalfist {

  class ThermalModelVDW : public ThermalModelBase
  {
  public:
    ThermalModelVDW(ThermalParticleSystem *TPS_, const ThermalModelParameters& params = ThermalModelParameters());

    virtual ~ThermalModelVDW(void);
    
    void FillVirialEV(const std::vector< std::vector<double> > & bij = std::vector< std::vector<double> >(0));

    void SetVirialdT(int i, int j, double dbdT) { if (i >= 0 && i < static_cast<int>(m_VirialdT.size()) && j >= 0 && j < static_cast<int>(m_VirialdT[i].size())) m_VirialdT[i][j] = dbdT; }
    
    void SetAttractiondT(int i, int j, double dadT) { if (i >= 0 && i < static_cast<int>(m_AttrdT.size()) && j >= 0 && j < static_cast<int>(m_AttrdT[i].size()))     m_AttrdT[i][j] = dadT; }

    double VirialCoefficientdT(int i, int j) const;
    
    double AttractionCoefficientdT(int i, int j) const;

    void SetTemperatureDependentAB(bool Tdep) { m_TemperatureDependentAB = Tdep; }
   
    bool TemperatureDependentAB() const { return m_TemperatureDependentAB; }
    

    virtual void SetMultipleSolutionsMode(bool search) { m_SearchMultipleSolutions = search; }

    bool UseMultipleSolutionsMode() const { return m_SearchMultipleSolutions; }

    double MuStar(int i) const { return m_MuStar[i]; }

    std::vector<double> GetMuStar() const { return m_MuStar; }
    
    void SetMuStar(const std::vector<double> & MuStar) { m_MuStar = MuStar; }
 
    // Override functions begin
 
    void FillChemicalPotentials();
 
    virtual void SetChemicalPotentials(const std::vector<double> & chem = std::vector<double>(0));
 
    virtual void FillVirial(const std::vector<double> & ri = std::vector<double>(0));
 
    virtual void FillAttraction(const std::vector< std::vector<double> > & aij = std::vector< std::vector<double> >(0));
 
    virtual void ReadInteractionParameters(const std::string &filename);
 
    virtual void WriteInteractionParameters(const std::string &filename);
 
    virtual void SetVirial(int i, int j, double b) { if (i >= 0 && i < static_cast<int>(m_Virial.size()) && j >= 0 && j < static_cast<int>(m_Virial[i].size())) m_Virial[i][j] = b; m_VDWComponentMapCalculated = false; }
    
    virtual void SetAttraction(int i, int j, double a) { if (i >= 0 && i < static_cast<int>(m_Attr.size()) && j >= 0 && j < static_cast<int>(m_Attr[i].size()))     m_Attr[i][j] = a; m_VDWComponentMapCalculated = false; }
 
    double VirialCoefficient(int i, int j) const;
 
    double AttractionCoefficient(int i, int j) const;
 
    virtual void ChangeTPS(ThermalParticleSystem *TPS);
 
    virtual void CalculatePrimordialDensities();
 
    virtual std::vector<double> CalculateChargeFluctuations(const std::vector<double> &chgs, int order = 4);
 
    virtual std::vector< std::vector<double> >  CalculateFluctuations(int order);
 
    void CalculateTwoParticleCorrelations();
 
    void CalculateFluctuations();
    
    virtual double CalculatePressure();
 
    virtual double CalculateEnergyDensity();
 
    virtual double CalculateEntropyDensity();
 
    virtual double CalculateEnergyDensityDerivativeT();
 
    virtual void CalculateTemperatureDerivatives();
 
    // Dummy
    virtual double CalculateBaryonMatterEntropyDensity() { return 0.; }
 
    virtual double CalculateMesonMatterEntropyDensity() { return 0.; }
 
    virtual double ParticleScalarDensity(int part);
 
    bool   IsLastSolutionOK() const { return m_LastBroydenSuccessFlag && m_LastCalculationSuccessFlag; }
 
    double DensityId(int index) { return m_DensitiesId[index]; }
 
    // Override functions end
 
    const std::vector< std::vector<int> >& VDWComponentIndices() const { return m_dMuStarIndices; }
    virtual double DeltaMu(int i) const { return MuShift(i); }
    const std::vector< std::vector<double> >& VirialMatrix() const { return m_Virial; }
    const std::vector< std::vector<double> >& AttractionMatrix() const { return m_Attr; }
 
  protected:
    std::vector<double> ComputeNp(const std::vector<double>& dmustar);
    
    std::vector<double> ComputeNp(const std::vector<double>& dmustar, const std::vector<double>& ns);

    void CalculateVDWComponentsMap();

    virtual std::vector<double> SearchSingleSolution(const std::vector<double> & muStarInit);

    std::vector<double> SearchMultipleSolutions(int iters = 300);

    void SolveEquations();

    virtual double MuShift(int id) const;

    std::vector<double> m_DensitiesId;

    std::vector<double> m_scaldens;

    std::vector< std::vector<double> > m_Virial;

    std::vector< std::vector<double> > m_Attr;

    std::vector< std::vector<double> > m_VirialdT;

    std::vector< std::vector<double> > m_AttrdT;

    bool   m_TemperatureDependentAB;

    bool   m_SearchMultipleSolutions;

    bool   m_LastBroydenSuccessFlag;

    bool   m_VDWComponentMapCalculated;

    std::vector<double> m_MuStar;
 
    std::vector<int> m_MapTodMuStar;
 
    std::vector<int> m_MapFromdMuStar;
 
    std::vector< std::vector<int> > m_dMuStarIndices;
 
    std::vector< std::vector<double> > m_chi;
 
    std::vector<double> m_chiarb;
  private:
    
    
    virtual void CalculatePrimordialDensitiesOld();
 
    virtual void CalculatePrimordialDensitiesNew();
 
    class BroydenEquationsVDW : public BroydenEquations
    {
    public:
      BroydenEquationsVDW(ThermalModelVDW *model) : BroydenEquations(), m_THM(model) { m_N = model->m_MapFromdMuStar.size(); }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      ThermalModelVDW *m_THM;
    };
 
    class BroydenJacobianVDW : public BroydenJacobian
    {
    public:
      BroydenJacobianVDW(ThermalModelVDW *model) : BroydenJacobian(), m_THM(model) { }
      std::vector<double> Jacobian(const std::vector<double> &x);
    private:
      ThermalModelVDW *m_THM;
    };
 
    class BroydenSolutionCriteriumVDW : public Broyden::BroydenSolutionCriterium
    {
    public:
      BroydenSolutionCriteriumVDW(ThermalModelVDW *model, double relative_error = Broyden::TOL) : Broyden::BroydenSolutionCriterium(relative_error), m_THM(model) { }
      virtual bool IsSolved(const std::vector<double>& x, const std::vector<double>& f, const std::vector<double>& xdelta = std::vector<double>()) const;
    protected:
      ThermalModelVDW *m_THM;
    };
  };

  typedef ThermalModelVDW ThermalModelVDWFull;

  void SetVDWHRGInteractionParameters(ThermalModelBase *model, double a, double b);

  void SetVDWHRGInteractionParameters(ThermalModelBase *model, const std::vector<std::vector<double>> & a, const std::vector<std::vector<double>> & b);
} // namespace thermalfist
 
#endif