! File: module_tabulated_eos.f90 ! Authors: Francesco Torsello (FT) !************************************************************************ ! Copyright (C) 2020-2023 Francesco Torsello * ! * ! This file is part of SPHINCS_ID * ! * ! SPHINCS_ID is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by * ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! SPHINCS_ID is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>. * ! The copy of the GNU General Public License should be in the file * ! 'COPYING'. * !************************************************************************ MODULE tabulated_eos !*********************************************************** ! !# This module contains the PROCEDURES needed to deal ! with tabulated |eos|. ! ! FT 10.03.2023 ! !*********************************************************** USE utility, ONLY: ios, err_msg, km2Msun_geo, MeV2amuc2 USE constants, ONLY: fm2cm, cm2km, MeV2erg, amu, Msun !USE units, ONLY: m0c2_cu, set_units ! needed for debugging IMPLICIT NONE DOUBLE PRECISION, PARAMETER:: fm2Msun_geo= fm2cm*cm2km*km2Msun_geo CONTAINS SUBROUTINE read_compose_beta_equilibrated_eos(namefile, table_eos) !********************************************** ! !# Read the |compose| files \*.thermo.ns, ! \*.nb.ns, containing the \(\beta\)-equilibrated ! |eos| computed by the |compose| software. ! ! See [the |compose| Manual](https://compose.obspm.fr/manual){:target="_blank"}, specifically Section 4.2.2 of ! v3.00 of the CompOSE Manual, for more details. ! ! FT 10.03.2023 ! !********************************************** IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN):: namefile !! Path of the |eos| file without '.thermo.ns' extension DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT):: table_eos !# This |eos| table contains the following quantities ! (except for the pressure and the baryon mass density, which are ! converted into code units, the other quantities are in units ! such that \(\hbar=c=k_\mathrm{B}=1\)): ! ! - First table index; temperature ! (not used for \(\beta\)-equilibrated |eos|) ! ! - Second table index; baryon number density ! ! - Third table index; charge fraction of strongly interacting particles ! (not used for \(\beta\)-equilibrated |eos|) ! ! - Pressure divided by the baryon number density: ! \(\dfrac{p}{|nb|}\mathrm{(MeV)}\). After reading it, ! it is converted into code units \(M_\odot c^2 L_\odot^{-3}\) ! ! - Entropy density divided by the baryon number density ! (entropy per baryon): \(\dfrac{s}{|nb|} \mathrm{(dimensionless)}\) ! ! - Scaled and shifted baryon chemical potential: ! \(\dfrac{\mu_\mathrm{b}}{m_\mathrm{n}} - 1 \mathrm{(dimensionless)}\) ! ! - Scaled charge chemical potential: ! \(\dfrac{\mu_\mathrm{q}}{m_\mathrm{n}} \mathrm{(dimensionless)}\) ! ! - Scaled effective lepton chemical potential: ! \(\dfrac{\mu_\mathrm{l}}{m_\mathrm{n}} \mathrm{(dimensionless)}\) ! ! - Scaled free energy per baryon: ! \(\dfrac{f}{|nb|m_\mathrm{n}} - 1 \mathrm{(dimensionless)}\) ! ! - Scaled internal energy per baryon (this is equal to the specific ! internal energy \(u\) used in |sphincsid| and |sphincsbssn|): ! \(\dfrac{e}{|nb|m_\mathrm{n}} - 1 \mathrm{(dimensionless)}\) ! ! - Number of additional quantities (this should be 2, as per CompOSE ! Manual v3.00, 13.03.2023) ! ! - Electron fraction \(Y_e\) in \(\beta\)-equilibrium ! ! - Enthalpy density without the temperature term \(Ts\): ! \(h= e + p \,(\mathrm{MeV}\,\mathrm{fm}^{-3})\) ! ! - Baryon number density \((\mathrm{fm}^{-3})\). After reading it, ! it is converted in baryon mass density in code units ! \(M_\odot L_\odot^{-3}\) ! INTEGER, PARAMETER:: unit_compose= 876 INTEGER, PARAMETER:: max_length_table= 10000 INTEGER:: itr, n_lines, tmp ! !-- Data contained in the *.thermo.ns file ! DOUBLE PRECISION:: m_n !! Mass of the neutron \(m_\mathrm{n} \mathrm{(MeV)}\) DOUBLE PRECISION:: m_p !! Mass of the proton \(m_\mathrm{p} \mathrm{(MeV)}\) INTEGER:: i_leptons !# Index that says if the |eos| considers leptons. ! ! - `i_leptons`= 1 if leptons are considered ! - `i_leptons`= anything else if leptons are not considered LOGICAL:: exist CHARACTER(LEN=:), ALLOCATABLE:: finalnamefile !PRINT * !PRINT *, "** Executing the read_compose_eos subroutine..." PRINT * ALLOCATE( table_eos(14,max_length_table) ) table_eos= 0.D0 finalnamefile= TRIM(namefile)//"eos.thermo.ns" INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) IF( exist )THEN OPEN( UNIT= unit_compose, FILE= TRIM(finalnamefile), & FORM= "FORMATTED", ACTION= "READ", IOSTAT= ios, & IOMSG= err_msg ) IF( ios > 0 )THEN PRINT *, "...error when opening " // TRIM(finalnamefile), & ". The error message is", err_msg PRINT * STOP ENDIF ELSE PRINT *, "** ERROR! Unable to find file " // TRIM(finalnamefile) PRINT * STOP ENDIF READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & m_n, m_p, i_leptons PRINT *, " * Reading file " // TRIM(finalnamefile) // "..." n_lines= 0 read_compose_table: DO itr= 1, max_length_table, 1 READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! Variables for .thermo.ns file table_eos(1,itr), table_eos(2,itr), table_eos(3,itr), & table_eos(4,itr), table_eos(5,itr), table_eos(6,itr), & table_eos(7,itr), table_eos(8,itr), table_eos(9,itr), & table_eos(10,itr), table_eos(11,itr), & table_eos(12,itr), table_eos(13,itr) IF( ios > 0 )THEN PRINT *, "...error when reading " // TRIM(finalnamefile), & ". The error message is", err_msg PRINT * STOP ENDIF IF( ios < 0 )THEN PRINT *, " * Reached end of file " // TRIM(finalnamefile) PRINT * EXIT read_compose_table ENDIF n_lines= n_lines + 1 ENDDO read_compose_table ! Reallocate arrays to delete unnecessary elements table_eos= table_eos(:,1:n_lines) CLOSE( unit_compose ) ! Read baryon number density from *.nb.ns file finalnamefile= TRIM(namefile)//"eos.nb.ns" INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) IF( exist )THEN OPEN( UNIT= unit_compose, FILE= TRIM(finalnamefile), & FORM= "FORMATTED", ACTION= "READ", IOSTAT= ios, & IOMSG= err_msg ) IF( ios > 0 )THEN PRINT *, "...error when opening " // TRIM(finalnamefile), & ". The error message is", err_msg PRINT * STOP ENDIF ELSE PRINT *, "** ERROR! Unable to find file " // TRIM(finalnamefile) PRINT * STOP ENDIF READ( UNIT= unit_compose, FMT= *, IOSTAT= ios, IOMSG= err_msg ) tmp, n_lines PRINT *, " * Reading file " // TRIM(finalnamefile) // "..." read_compose_nb: DO itr= 1, n_lines, 1 READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! Variable for .nb file ! (baryon number density, independent variable) table_eos(14,itr) IF( ios > 0 )THEN PRINT *, "...error when reading " // TRIM(finalnamefile), & ". The error message is", err_msg PRINT * STOP ENDIF ENDDO read_compose_nb CLOSE( unit_compose ) PRINT *, "...done." PRINT * ! !-- Convert (some of the) physical quantities to desired quantities !-- in the desiredunits ! ! Pressure: MeV fm^{-3} to Msun_geo c^2 Msun_geo^{-3} ! (remember that c=1, so we do not have to divide by it) table_eos(4,:)= table_eos(4,:)*table_eos(14,:) & *(MeV2amuc2*amu/Msun)/(fm2Msun_geo**3) ! Baryon mass density: MeV fm^{-3} to Msun_geo Msun_geo^{-3} ! (remember that c=1, so we do not have to divide by it) table_eos(14,:)= m_n*table_eos(14,:) & *(MeV2amuc2*amu/Msun)/(fm2Msun_geo**3) !CALL set_units('NSM') ! !PRINT *, "n_lines:", n_lines !PRINT *, "SIZE(pressure):", SIZE(table_eos(4,:)) !PRINT *, "SIZE(baryon mass density):", SIZE(table_eos(14,:)) !PRINT *, "SIZE(specific internal energy):", SIZE(table_eos(10,:)) !PRINT *, "pressure:", table_eos(4,210)/lorene2hydrobase !PRINT *, "baryon mass density:", table_eos(14,210)/lorene2hydrobase !PRINT *, "specific internal energy:", table_eos(10,210) !PRINT *, "m0c2_cu:", m0c2_cu !PRINT *, "pressure*m0c2_cu:", table_eos(4,210)/m0c2_cu !PRINT *, "baryon mass density*m0c2_cu:", table_eos(14,210)/m0c2_cu !PRINT * !STOP !PRINT *, "** Subroutine read_compose_eos executed." !PRINT * END SUBROUTINE read_compose_beta_equilibrated_eos END MODULE tabulated_eos ! These code may be useful in the future, if the non-beta-equilibrated EOS ! from CompOSE should be read from the *.thermo, *.nb, *.y, *.yq files ! ! ! Read charge fraction of strongly interacting particles from *.t file ! ! (we don't have charge fraction of strongly interacting particles for the ! ! beta equilibrated data) ! finalnamefile= TRIM(namefile)//".yq" ! ! INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) ! ! IF( exist )THEN ! OPEN( UNIT= unit_compose, FILE= TRIM(finalnamefile), & ! FORM= "FORMATTED", ACTION= "READ", IOSTAT= ios, & ! IOMSG= err_msg ) ! IF( ios > 0 )THEN ! PRINT *, "...error when opening " // TRIM(finalnamefile), & ! ". The error message is", err_msg ! PRINT * ! STOP ! ENDIF ! ELSE ! PRINT *, "** ERROR! Unable to find file " // TRIM(finalnamefile) ! PRINT * ! STOP ! ENDIF ! ! READ( UNIT= unit_compose, FMT= *, IOSTAT= ios, IOMSG= err_msg )tmp, n_lines ! ! PRINT *, " * Reading file " // TRIM(finalnamefile) // "..." ! ! read_compose_nb: DO itr= 1, n_lines, 1 ! READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! ! Variable for .nb file ! ! (charge fraction of strongly ! ! interacting particles, independent variable) ! table_eos(6,itr) ! IF( ios > 0 )THEN ! PRINT *, "...error when reading " // TRIM(finalnamefile), & ! ". The error message is", err_msg ! PRINT * ! STOP ! ENDIF ! ! ENDDO read_compose_nb ! ! CLOSE( unit_compose ) ! ! Read temperature from *.t file (we don't have temperature for the ! ! beta equilibrated data) ! finalnamefile= TRIM(namefile)//".t" ! ! INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) ! ! IF( exist )THEN ! OPEN( UNIT= unit_compose, FILE= TRIM(finalnamefile), & ! FORM= "FORMATTED", ACTION= "READ", IOSTAT= ios, & ! IOMSG= err_msg ) ! IF( ios > 0 )THEN ! PRINT *, "...error when opening " // TRIM(finalnamefile), & ! ". The error message is", err_msg ! PRINT * ! STOP ! ENDIF ! ELSE ! PRINT *, "** ERROR! Unable to find file " // TRIM(finalnamefile) ! PRINT * ! STOP ! ENDIF ! ! READ( UNIT= unit_compose, FMT= *, IOSTAT= ios, IOMSG= err_msg )tmp, n_lines ! ! PRINT *, " * Reading file " // TRIM(finalnamefile) // "..." ! ! read_compose_t: DO itr= 1, n_lines, 1 ! READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! ! Variable for .t file (temperature, independent variable) ! table_eos(4,itr) ! IF( ios > 0 )THEN ! PRINT *, "...error when reading " // TRIM(finalnamefile), & ! ". The error message is", err_msg ! PRINT * ! STOP ! ENDIF ! ! ENDDO read_compose_t ! ! CLOSE( unit_compose )