! File: submodule_id_base_mass_profile.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 (id_base) mass_profile !******************************************** ! !# Implementation of the method of TYPE idbase ! that integrates the baryon mass density to ! extract the radial baryon mass profile. ! ! FT 12.07.2021 ! !******************************************** IMPLICIT NONE CONTAINS !-------------------! !-- SUBROUTINES --! !-------------------! MODULE PROCEDURE integrate_baryon_mass_density !************************************************ ! !# Perform 3D integration over a spherical grid ! of the baryon mass density. Output baryon ! mass and radial mass profile. ! ! FT 19.02.2021 ! ! Upgraded to ellipsoidal grid ! ! FT 15.11.2022 ! !************************************************ USE constants, ONLY: pi USE utility, ONLY: zero, one, two, three, four USE NR, ONLY: indexx USE tensor, ONLY: jxx, jxy, jxz, jyy, jyz, jzz USe utility, ONLY: determinant_sym3x3, cartesian_from_spherical USE numerics, ONLY: bilinear_interpolation IMPLICIT NONE INTEGER:: r, th, phi DOUBLE PRECISION:: rad, rad_coord, colat, long, mass_element, max_radius DOUBLE PRECISION:: sq_g, baryon_density, gamma_euler DOUBLE PRECISION:: a_x, a_y, a_z, xtemp, ytemp, ztemp DOUBLE PRECISION, DIMENSION(6):: g CHARACTER(LEN=:), ALLOCATABLE:: surface_type LOGICAL, PARAMETER:: debug= .FALSE. a_x= one a_y= one a_z= one max_radius= radius surface_type= "spherical" IF(PRESENT(surf))THEN IF(surf% is_known)THEN surface_type= "oval" ENDIF ENDIF IF(PRESENT(radii))THEN IF(.NOT.PRESENT(surf)) surface_type= "ellipsoidal" IF(PRESENT(surf))THEN IF(.NOT.(surf% is_known))THEN surface_type= "ellipsoidal" ENDIF ENDIF max_radius= MAXVAL([radius,radii(1),radii(2)]) IF(max_radius == radius)THEN a_x= one a_y= radii(1)/max_radius a_z= radii(2)/max_radius ELSEIF(max_radius == radii(1))THEN a_x= radius/max_radius a_y= one a_z= radii(2)/max_radius ELSEIF(max_radius == radii(2))THEN a_x= radius/max_radius a_y= radii(1)/max_radius a_z= one ENDIF ENDIF IF(debug) PRINT *, "inside integrate_baryon_mass_density, ", & "surface_type is:", surface_type ! STOP mass_profile( 1, 0 )= zero mass_profile( 2, 0 )= four/three*pi*dr**three*central_density mass_profile( 3, 0 )= four/three*pi*dr**three*central_density !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP SHARED(dr, dphi, dth, center, max_radius, & !$OMP mass_profile, this, a_x, a_y, a_z, surf) & !$OMP PRIVATE(r, th, phi, rad, rad_coord, long, colat, sq_g, & !$OMP gamma_euler, g, baryon_density, mass_element, & !$OMP mass, xtemp, ytemp, ztemp) radius_loop: DO r= 1, NINT(max_radius/dr), 1 mass= zero rad_coord= r*dr longitude_loop: DO phi= 1, NINT(two*pi/dphi), 1 long= phi*dphi colatitude_loop: DO th= 1, NINT(pi/two/dth), 1 colat= th*dth ! The definition of the baryon mass for the LORENE ID is in eq.(69) ! of Gourgoulhon et al., PRD 63 064029 (2001) rad= rad_coord IF(PRESENT(surf))THEN IF(surf% is_known)THEN rad= bilinear_interpolation( colat, long, & SIZE(surf% points(:,1,5)), & SIZE(surf% points(1,:,6)), & surf% points(:,:,5:6), surf% points(:,:,4) ) rad= rad*rad_coord/max_radius ENDIF ENDIF CALL cartesian_from_spherical( & a_x*(rad + dr), colat, long, & center(1), center(2), center(3), & xtemp, ytemp, ztemp, a_y/a_x, a_z/a_x ) !CALL this% read_id_mass_b( & ! center(1) + a_x*(rad_coord + dr)*SIN(colat)*COS(long), & ! center(2) + a_y*(rad_coord + dr)*SIN(colat)*SIN(long), & ! center(3) + a_z*(rad_coord + dr)*COS(colat), & ! g, baryon_density, gamma_euler ) CALL this% read_id_mass_b( xtemp, ytemp, ztemp, & g, baryon_density, gamma_euler ) IF( ISNAN( g(jxx) ) .OR. ISNAN( g(jxy) ) .OR. ISNAN( g(jxz) ) & .OR. ISNAN( g(jyy) ) .OR. ISNAN( g(jyz) ) .OR. ISNAN( g(jzz) ) & .OR. ISNAN( baryon_density ) .OR. ISNAN( gamma_euler ) ) & CYCLE ! Compute square root of the determinant of the spatial metric CALL determinant_sym3x3(g, sq_g) sq_g= SQRT(sq_g) !mass_element= a_x*a_y*a_z*(rad_coord**two)*SIN(colat)*dr*dth*dphi & ! *sq_g*gamma_euler*baryon_density mass_element= a_x*a_y*a_z*(rad**two)*SIN(colat)*dr*dth*dphi & *sq_g*gamma_euler*baryon_density mass= mass + two*mass_element ENDDO colatitude_loop ENDDO longitude_loop mass_profile(1, r)= rad_coord mass_profile(2, r)= mass ENDDO radius_loop !$OMP END PARALLEL DO DO r= 1, NINT(radius/dr), 1 mass_profile(3, r)= mass_profile(3, r - 1) + mass_profile(2, r) ENDDO mass= mass_profile(3, NINT(radius/dr)) IF( ISNAN(mass) )THEN PRINT *, "** ERROR! The integrated mass is a NaN!" PRINT * STOP ELSEIF( mass <= 0 )THEN PRINT *, "** ERROR! The integrated mass is mass=", mass, "<= 0!" PRINT * STOP ENDIF PRINT *, " * Radius covered by the integration of baryon mass density=", & MAXVAL( mass_profile( 1, : ), DIM= 1 ) PRINT *, " * Integrated baryon mass of the star=", mass PRINT * CALL indexx( NINT(radius/dr) + 1, mass_profile( 1, : ), mass_profile_idx ) END PROCEDURE integrate_baryon_mass_density END SUBMODULE mass_profile ! TODO: deprecated? It computes the relative Lorentz factor between the fluid ! and the Eulerian observer from the velocity ! ! CALL bns_obj% import_id( & ! center1 + rad_coord*SIN(lat)*COS(long), & ! rad_coord*SIN(lat)*SIN(long), & ! rad_coord*COS(lat), & ! g_xx, baryon_density, & ! gamma_euler ) ! ! ! Compute covariant spatial fluid velocity (metric is diagonal and ! ! conformally flat) ! !v_euler_x_l= g_xx*v_euler_x ! !v_euler_y_l= g_xx*v_euler_y ! !v_euler_z_l= g_xx*v_euler_z ! ! ! !! Compute the corresponding Lorentz factor ! !lorentz_factor= 1.0D0/SQRT( 1.0D0 - ( v_euler_x_l*v_euler_x & ! ! + v_euler_y_l*v_euler_y & ! ! + v_euler_z_l*v_euler_z ) ) ! ! ! !! Compute covariant fluid 4-velocity ! !u_euler_t_l= lorentz_factor *( - lapse + v_euler_x_l*shift_x & ! ! + v_euler_y_l*shift_y & ! ! + v_euler_z_l*shift_z ) ! !u_euler_x_l= lorentz_factor*v_euler_x_l ! !u_euler_y_l= lorentz_factor*v_euler_y_l ! !u_euler_z_l= lorentz_factor*v_euler_z_l ! ! ! !! Compute vector normal to spacelike hypersurface ! !! (4-velocity of the Eulerian observer) ! !n_t= 1.0D0/lapse ! !n_x= - shift_x/lapse ! !n_y= - shift_y/lapse ! !n_z= - shift_z/lapse ! ! ! !! Compute relative Lorentz factor between 4-velocity of the fluid ! !! wrt the Eulerian observer and the 4-velocity of the Eulerian observer ! !lorentz_factor_rel= - ( n_t*u_euler_t_l + n_x*u_euler_x_l & ! ! + n_y*u_euler_y_l + n_z*u_euler_z_l )