! File: submodule_sph_particles_compose.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'. * !************************************************************************ SUBMODULE (sph_particles) compose !*************************************************** ! !# This SUBMODULE contains the implementation of ! the methods of TYPE [[particles]] ! that compute \(Y_e\) on the particles, using the ! data from the |compose| database ! ! <https://compose.obspm.fr/> ! ! FT 12.07.2021 ! !*************************************************** IMPLICIT NONE CONTAINS !-------------------! !-- SUBROUTINES --! !-------------------! MODULE PROCEDURE read_compose_composition !************************************************ ! !# Reads the electron fraction \(Y_e\) = \(n_e/n_b\), ! with \(n_e\) electron number density and \(n_b\) ! baryon number density, from the .compo file ! taken from the |compose| database of EoS. ! \(Y_e\) is given as a function of \(T\), ! \(n_b\), \(Y_q\) on ! a grid; the computation of \(Y_e\) on the stars is ! done by the SUBROUTINE compute_Ye_on_stars. ! ! FT 1.03.2021 ! ! ! @warning DEPRECATED ! ! FT 13.03.2023 ! !************************************************ USE constants, ONLY: fm2cm, cm2km USE utility, ONLY: km2Msun_geo, zero IMPLICIT NONE ! The commented variables might be useful in the future INTEGER:: itr, cntr!, & !i_t, i_nb, i_phase, n_pairs, i_e, i_n, Y_n, & !n_quad, & !i_i, A_i, Z_i, Y_i, i_leptons, i_ns INTEGER, PARAMETER:: unit_compose= 56 INTEGER, PARAMETER:: max_length_eos= 10000 !DOUBLE PRECISION:: m_n, m_p, & ! p_nb, s_nb, mub_mn, muq_mn, mul_mn, f_nbmn, e_nbmn, h !DOUBLE PRECISION, DIMENSION( : ), ALLOCATABLE:: n_b, Y_e LOGICAL:: exist CHARACTER(LEN=:), ALLOCATABLE:: finalnamefile PRINT *, "** Executing the read_compose_composition subroutine..." ALLOCATE( this% nb_table( max_length_eos ) ) ALLOCATE( this% Ye_table( max_length_eos ) ) this% nb_table= zero this% Ye_table= zero IF( PRESENT(namefile) )THEN finalnamefile= TRIM(namefile)//".beta" ELSE finalnamefile= "../../CompOSE_EOS/SFHO_with_electrons/eos.beta" ENDIF 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 STOP ENDIF !CALL test_status( ios, err_msg, "...error when opening " & ! // TRIM(finalnamefile) ) ELSE PRINT *, "** ERROR! Unable to find file " // TRIM(finalnamefile) STOP ENDIF !READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! m_n, m_p, i_leptons PRINT *, " * Reading file " // TRIM(finalnamefile) // "..." cntr= 1 read_compose_beta: DO itr= 1, max_length_eos, 1 READ( UNIT= unit_compose, FMT= *, IOSTAT = ios, IOMSG= err_msg ) & ! Variables for .thermo.ns file !i_T, i_nb, i_yq, & !p_nb, s_nb, mub_mn, muq_mn, mul_mn, f_nbmn, e_nbmn, & !i_ns, this% Ye_table(itr), h ! Variables for .beta file this% nb_table(itr), & this% Ye_table(itr) ! Variables for .compo file !i_T, i_nb, i_yq, & !i_phase, n_pairs, & !i_e, Y_e(itr)!, & !i_n, Y_n, & !n_quad, & !i_i, A_i, Z_i, Y_i IF( ios > 0 )THEN PRINT *, "...error when reading " // TRIM(finalnamefile), & ". The error message is", err_msg STOP ENDIF IF( ios < 0 )THEN PRINT *, " * Reached end of file " // TRIM(finalnamefile) EXIT ENDIF cntr= cntr + 1 ENDDO read_compose_beta !PRINT *, "cntr= ", cntr ! Reallocate the arrays to delete the empty elements this% nb_table= this% nb_table( 1:cntr )/(fm2cm**3*cm2km**3*km2Msun_geo**3) this% Ye_table= this% Ye_table( 1:cntr ) !PRINT *, i_T, i_nb, i_yq, i_phase, n_pairs, i_e !PRINT *, "SIZE(n_b)= ", SIZE(this% n_b), "SIZE(Y_e)= ", SIZE(this% Y_e) !PRINT *, "n_b(1)= ", this% n_b(1), "Y_e(1)= ", this% Y_e(1) !PRINT *, "n_b(cntr)= ", this% n_b(cntr), "Y_e(cntr)= ", this% Y_e(cntr) !STOP CLOSE( unit_compose ) PRINT *, "** Subroutine read_compose_composition executed." PRINT * END PROCEDURE read_compose_composition MODULE PROCEDURE compute_Ye !************************************************ ! !# Interpolates the electron fraction ! \(Y_e\) = \(n_e/n_b\) ! at the particle positions, using the data ! read by read_compose_composition. ! ! FT 3.03.2021 ! ! Uses new infrastructure ! ! FT 13.03.2023 ! !************************************************ USE numerics, ONLY: linear_interpolation USE units, ONLY: m0c2_cu IMPLICIT NONE INTEGER:: dim_table, a, i_matter DOUBLE PRECISION:: min_nb_table, max_nb_table DO i_matter= 1, this% n_matter, 1 !dim_table= SIZE(this% nb_table) !min_nb_table= MINVAL(this% nb_table, DIM= 1) !max_nb_table= MAXVAL(this% nb_table, DIM= 1) dim_table= SIZE(this% all_eos(i_matter)% table_eos(1,:)) min_nb_table= MINVAL(this% all_eos(i_matter)% table_eos(14,:), DIM= 1) max_nb_table= MAXVAL(this% all_eos(i_matter)% table_eos(14,:), DIM= 1) !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( this, min_nb_table, max_nb_table, dim_table, & !$OMP i_matter, m0c2_cu ) & !$OMP PRIVATE( a ) particle_loop: DO a= 1, this% npart, 1 IF( this% nlrf(a)*m0c2_cu < min_nb_table )THEN PRINT *, "** ERROR! The value of nlrf(", a, ")=", & this% nlrf(a)*m0c2_cu, & "is lower than the minimum value in the table =", min_nb_table PRINT *, " * Is nlrf computed when you call this SUBROUTINE? " // & "If yes, please select a table with a wider range." PRINT * STOP ELSEIF( this% nlrf(a)*m0c2_cu > max_nb_table )THEN PRINT *, "** ERROR! The value of nlrf(", a, ")=", & this% nlrf(a)*m0c2_cu, & "is larger than the maximum value in the table =", max_nb_table PRINT *, " * Is nlrf computed when you call this SUBROUTINE? " // & "If yes, please select a table with a wider range." PRINT * STOP ENDIF !Ye_linear_interpolation_loop: DO i_table= 1, dim_table - 1, 1 ! ! IF( this% nb_table(i_table) < this% nlrf(a) .AND. & ! this% nlrf(a) < this% nb_table(i_table + 1) )THEN ! ! this% Ye(a)= this% Ye_table(i_table) & ! + (this% Ye_table(i_table + 1) - this% Ye_table(i_table)) & ! /(this% nb_table(i_table + 1) - this% nb_table(i_table)) & ! *(this% nlrf(a) - this% nb_table(i_table)) ! EXIT ! ! ENDIF ! !ENDDO Ye_linear_interpolation_loop this% Ye(a)= linear_interpolation & (this% nlrf(a)*m0c2_cu, dim_table, & this% all_eos(i_matter)% table_eos(14,:), & this% all_eos(i_matter)% table_eos(12,:)) ENDDO particle_loop !$OMP END PARALLEL DO ENDDO END PROCEDURE compute_Ye END SUBMODULE compose