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; }
    

    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));

    void FillVirial(const std::vector<double> & ri = std::vector<double>(0));

    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);

    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; }
    
    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; }

    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();

    // 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

  protected:
    std::vector<double> ComputeNp(const std::vector<double>& dmustar);
    
    std::vector<double> ComputeNp(const std::vector<double>& dmustar, const std::vector<double>& ns);

    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;

    std::vector<double> m_MuStar;

    std::vector<int> m_MapTodMuStar;

    std::vector<int> m_MapFromdMuStar;

    std::vector< std::vector<int> > m_dMuStarIndices;

  private:
    std::vector< std::vector<double> > m_chi;

    std::vector<double> m_chiarb;
    
    
    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;

} // namespace thermalfist

#endif