! File: module_wd_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 wd_eos !******************************************** ! !# This MODULE implements the Chandrasekhar's ! degenerate |eos| for white dwarfs. ! ! See [Benz W., Bowers R.L., Cameron A.G.W., Press W.H., 1990, APJ, 348, 647.](https://ui.adsabs.harvard.edu/abs/1990ApJ...348..647B/abstract){:target="_blank"}, eqs.(2.4)-(2.5) ! ! See also [S. Chandrasekhar, The Highly Collapsed Configurations of a Stellar Mass. (Second Paper.), Monthly Notices of the Royal Astronomical Society, Volume 95, Issue 3, January 1935, Pages 207225](https://doi.org/10.1093/mnras/95.3.207){:target="_blank"} ! ! FT 19.12.2022 ! !******************************************** USE constants, ONLY: third, c_light2, kg2g, cm2m, press_si2cgs, m2cm USE units, ONLY: m0c2_cu IMPLICIT NONE PRIVATE DOUBLE PRECISION, PARAMETER:: dens_si2cu= 1.618654158231174D-21 !! Conversion factor for the baryon mass density, from SI units to code units DOUBLE PRECISION, PARAMETER:: mu_e= 2.D0 !! Mean molecular weight per electron DOUBLE PRECISION, PARAMETER:: a_wd_cgs= 6.02D+21*press_si2cgs !! Constant with dimensions of a pressure in CGS units ! (only used in [[test_wd_eos_cgs]]) DOUBLE PRECISION, PARAMETER:: b_wd_cgs= mu_e*9.82D+8*kg2g/(m2cm**3) !! Constant with dimensions of a density in CGS units ! (only used in [[test_wd_eos_cgs]]) DOUBLE PRECISION, PARAMETER:: a_wd= 6.02D+21/(c_light2*cm2m)*dens_si2cu !# Constant with dimensions of a pressure in code units DOUBLE PRECISION, PARAMETER:: b_wd= mu_e*9.82D+8*dens_si2cu !# Constant with dimensions of a density in code units PUBLIC:: pr_wd, u_wd, rho_wd, test_wd_eos_cgs CONTAINS !--------------------------! !-- PRIVATE PROCEDURES --! !--------------------------! DOUBLE PRECISION PURE FUNCTION f_wd(x) !******************************************** ! !# Function used in the computation of ! the degenerate pressure ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: x f_wd= x*(2.D0*x**2 - 3.D0)*SQRT(x**2 + 1.D0) + 3.D0*ASINH(x) END FUNCTION f_wd DOUBLE PRECISION PURE FUNCTION g_wd(x) !******************************************** ! !# Function used in the computation of ! the degenerate specific internal energy ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: x g_wd= 2.D0*4.D0*x**3*(SQRT(x**2 + 1.D0) - 1.D0) - f_wd(x) END FUNCTION g_wd !-------------------------! !-- PUBLIC PROCEDURES --! !-------------------------! DOUBLE PRECISION PURE FUNCTION pr_wd(rho) !******************************************** ! !# Degenerate pressure as a function of density ! ! See [Benz W., Bowers R.L., Cameron A.G.W., Press W.H., 1990, APJ, 348, 647.](https://ui.adsabs.harvard.edu/abs/1990ApJ...348..647B/abstract){:target="_blank"}, eqs.(2.4) ! ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: rho DOUBLE PRECISION:: x x= (rho/b_wd)**third pr_wd= a_wd*f_wd(x) END FUNCTION pr_wd DOUBLE PRECISION PURE FUNCTION pr_wd_cgs(rho) !******************************************** ! !# Degenerate pressure as a function of density, ! all in CGS units ! ! See [Benz W., Bowers R.L., Cameron A.G.W., Press W.H., 1990, APJ, 348, 647.](https://ui.adsabs.harvard.edu/abs/1990ApJ...348..647B/abstract){:target="_blank"}, eqs.(2.4) ! ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: rho DOUBLE PRECISION:: x x= (rho/b_wd_cgs)**third pr_wd_cgs= a_wd_cgs*f_wd(x) END FUNCTION pr_wd_cgs DOUBLE PRECISION PURE FUNCTION u_wd(rho) !******************************************** ! !# Degenerate specific internal energy as a ! function of density ! ! See [Benz W., Bowers R.L., Cameron A.G.W., Press W.H., 1990, APJ, 348, 647.](https://ui.adsabs.harvard.edu/abs/1990ApJ...348..647B/abstract){:target="_blank"}, eqs.(2.5) ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: rho DOUBLE PRECISION, PARAMETER:: rho_min= 1.D-35 DOUBLE PRECISION:: x IF(rho == 0)THEN u_wd= 0.D0 RETURN ENDIF x= (rho/b_wd)**third u_wd= a_wd*g_wd(x)/rho END FUNCTION u_wd DOUBLE PRECISION FUNCTION rho_wd(pr, rho_left_bracket, rho_right_bracket) !******************************************** ! !# Degenerate density as a function of pressure ! ! Use bisection to find the value of the density, ! given the pressure ! ! FT 19.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN):: pr DOUBLE PRECISION, INTENT(IN):: rho_left_bracket DOUBLE PRECISION, INTENT(IN):: rho_right_bracket INTEGER, PARAMETER:: tolerance_magnitude= 3 DOUBLE PRECISION, PARAMETER:: tolerance_pr = 1.D-5 DOUBLE PRECISION, PARAMETER:: tolerance_rho = 1.D-10 DOUBLE PRECISION, PARAMETER:: pr_min = 1.D-30 INTEGER:: cnt DOUBLE PRECISION:: rho_left, rho_right, rho_mean, & pr_left, pr_right, pr_mean, pr_cgs LOGICAL, PARAMETER:: debug= .FALSE. IF(pr < pr_min)THEN rho_wd= 0.D0 RETURN ENDIF rho_left = FLOOR(LOG10(rho_left_bracket)) rho_right= CEILING(LOG10(rho_right_bracket)) ! Bisection in logarithmic scale, to find the approximate order of magnitude DO pr_left = pr_wd(1.D1**rho_left) pr_right= pr_wd(1.D1**rho_right) IF( pr_left <= pr .AND. pr_right > pr )THEN rho_mean= FLOOR((rho_left + rho_right)/2.D0) pr_mean = pr_wd(1.D1**rho_mean) IF(debug) PRINT *, FLOOR(ABS(LOG10(pr_right) - LOG10(pr_left))) IF(debug) PRINT *, LOG10(pr_right) IF(debug) PRINT *, LOG10(pr_left) IF( FLOOR(ABS(LOG10(pr_right) - LOG10(pr_left))) & <= tolerance_magnitude )THEN EXIT ENDIF IF(pr_mean < pr)THEN rho_left = rho_mean ELSE rho_right= rho_mean ENDIF ELSEIF( pr_left < pr .AND. pr_right < pr )THEN rho_right= rho_right + 2.D0 ELSEIF( pr_left > pr .AND. pr_right > pr )THEN rho_left= rho_left - 2.D0 ELSE PRINT * PRINT *, "** ERROR in SUBROUTINE rho_wd in MODULE wd_eos!" PRINT *, "rho_left=", rho_left PRINT *, "rho_right=", rho_right PRINT *, "pr_left=", pr_left PRINT *, "pr=", pr PRINT *, "pr_right=", pr_right PRINT * STOP ENDIF ENDDO ! Bisection in linear scale, to find the precise value pr_cgs = pr rho_left = (1.D1**rho_left ) rho_right= (1.D1**rho_right) IF(debug) PRINT *, "rho_left =", rho_left IF(debug) PRINT *, "rho_right=", rho_right IF(debug) PRINT *, "pr_cgs =", pr_cgs cnt= 0 DO pr_left = pr_wd(rho_left) pr_right= pr_wd(rho_right) IF( pr_left <= pr_cgs .AND. pr_right > pr_cgs )THEN rho_mean= (rho_left + rho_right)/2.D0 pr_mean = pr_wd(rho_mean) IF(debug) PRINT *, "pr_left =", pr_left IF(debug) PRINT *, "pr_right =", pr_right IF(debug) PRINT *, "rho_mean =", rho_mean IF(debug) PRINT *, "pr_mean =", pr_mean IF(debug .AND. cnt > 100) & PRINT *,"ABS((pr_left - pr_right)/pr_right)=", & ABS((pr_left - pr_right)/pr_right) IF(debug .AND. cnt > 100) PRINT *, "pr_left =", pr_left IF(debug .AND. cnt > 100) PRINT *, "pr_right =", pr_right IF(debug .AND. cnt > 100) PRINT *, "pr_mean =", pr_mean IF(debug .AND. cnt > 100) PRINT *, "rho_left =", rho_left IF(debug .AND. cnt > 100) PRINT *, "rho_right=", rho_right IF(debug .AND. cnt > 100) PRINT *, "rho_mean =", rho_mean IF(debug .AND. cnt > 100) PRINT *, "pr_cgs =", pr_cgs IF(debug) STOP IF( ABS((pr_left - pr_right)/pr_right) < tolerance_pr & .OR. & ABS((rho_left - rho_right)/rho_right) < tolerance_rho )THEN EXIT ENDIF IF(pr_mean < pr_cgs)THEN rho_left = rho_mean ELSE rho_right= rho_mean ENDIF ELSEIF( pr_left < pr_cgs .AND. pr_right < pr_cgs )THEN rho_right= 2.D0*rho_right ELSEIF( pr_left > pr_cgs .AND. pr_right > pr_cgs )THEN rho_left= rho_left/2.D0 ELSE PRINT * PRINT *, "** ERROR in SUBROUTINE rho_wd in MODULE wd_eos!" PRINT *, "rho_left=", rho_left PRINT *, "rho_right=", rho_right PRINT *, "pr_left=", pr_left PRINT *, "pr_cgs=", pr_cgs PRINT *, "pr_right=", pr_right PRINT * STOP ENDIF cnt= cnt + 1 IF(cnt > 150)THEN STOP ENDIF ENDDO IF( ABS(pr - pr_wd(rho_mean))/pr > tolerance_pr )THEN PRINT *, "** ERROR! The value of rho_mean found by FUNCTION rho_wd ", & "in MODULE wd_eos, does not give a value of pr_wd(rho_mean) ", & "compatible with the input pressure pr, within the required ", & "tolerance." PRINT * PRINT *, "rho_mean=", rho_mean PRINT *, "pr=", pr PRINT *, "pr_wd(rho_mean)=", pr_wd(rho_mean) PRINT *, "ABS(pr - pr_wd(rho_mean))/pr", ABS(pr - pr_wd(rho_mean))/pr PRINT *, "tolerance_pr=", tolerance_pr PRINT * STOP ENDIF rho_wd= rho_mean END FUNCTION rho_wd SUBROUTINE test_wd_eos_cgs(rho_input, rho, pr, u) !******************************************** ! !# Test that the implementation is correct. ! Takes the density in code units as input, ! and assigns CGS values to the pressure, ! and the density recomputed from the pressure. ! It also computed the dimensionless ! specific internal energy. ! ! FT 20.12.2022 ! !******************************************** IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: rho_input DOUBLE PRECISION, INTENT(OUT):: rho, pr, u PRINT *, "** Testing that the implementation of the Chandrasekhar EOS", & " for white dwarfs is correct." PRINT * PRINT *, " a_wd in code units=", a_wd PRINT *, " b_wd in code units=", b_wd PRINT *, " a_wd in CGS units=", a_wd_cgs PRINT *, " b_wd in CGS units=", b_wd_cgs PRINT * PRINT *, " rho_input in code units=", rho_input PRINT *, " rho_input in CGS units=", & rho_input/dens_si2cu*kg2g/(m2cm**3) PRINT * pr = pr_wd(rho_input) u = u_wd(rho_input) rho= rho_wd(pr,0.D0,1.D0) PRINT *, " pr in code units=", pr pr= pr/dens_si2cu*kg2g/(m2cm**3)*c_light2 PRINT *, " pr in CGS units=", pr PRINT * PRINT *, " u (dimensionless)=", u PRINT * PRINT *, " Recomputed rho in code units=", rho rho= rho/dens_si2cu*kg2g/(m2cm**3) PRINT *, " Recomputed rho in CGS units=", rho PRINT * END SUBROUTINE test_wd_eos_cgs END MODULE wd_eos