! File: submodule_ejecta_generic_interpolate.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 (ejecta_generic) interpolate !**************************************************** ! !# Implementation of the methods of TYPE ejecta that ! interpolate the data on a grid, to a generic point ! ! FT 19.11.2021 ! !**************************************************** USE utility, ONLY: zero, one IMPLICIT NONE CONTAINS !-------------------! !-- SUBROUTINES --! !-------------------! MODULE PROCEDURE interpolate_id_full !************************************************** ! !# Stores the |id| in non-[[ejecta]]-member arrays ! with the same shape as the [[ejecta]] member arrays ! ! FT 19.11.2021 ! !************************************************** IMPLICIT NONE END PROCEDURE interpolate_id_full MODULE PROCEDURE interpolate_id_spacetime !******************************************************* ! !# Stores the spacetime |id| in multi-dimensional arrays ! needed to compute the BSSN variables and constraints ! ! FT 19.11.2021 ! !******************************************************* USE tensor, ONLY: jxx, jxy, jxz, jyy, jyz, jzz, jx, jy, jz, n_sym4x4 IMPLICIT NONE INTEGER:: i, j, k PRINT *, "** N.B. : This SUBROUTINE is not yet implemented!" PRINT *, " * Stopping..." PRINT * STOP !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( nx, ny, nz, this, pos, & !$OMP lapse, shift, g, ek ) & !$OMP PRIVATE( i, j, k ) DO k= 1, nz, 1 DO j= 1, ny, 1 DO i= 1, nx, 1 g( i, j, k, jyy )= one g( i, j, k, jyy )= one g( i, j, k, jzz )= one g( i, j, k, jxy )= zero g( i, j, k, jxz )= zero g( i, j, k, jyz )= zero ! !- Set/unset the geodesic gauge ! IF( this% get_one_lapse() )THEN lapse( i, j, k )= one ENDIF IF( this% get_zero_shift() )THEN shift( i, j, k, jx )= zero shift( i, j, k, jy )= zero shift( i, j, k, jz )= zero ENDIF ! !-- Convert the extrinsic curvature from |lorene| units to !-- |sphincs| units ! ek( i, j, k, jxx )= zero ek( i, j, k, jxy )= zero ek( i, j, k, jxz )= zero ek( i, j, k, jyy )= zero ek( i, j, k, jyz )= zero ek( i, j, k, jzz )= zero ENDDO ENDDO ENDDO !$OMP END PARALLEL DO END PROCEDURE interpolate_id_spacetime MODULE PROCEDURE interpolate_id_hydro !******************************************************* ! !# Stores the hydro |id| in the arrays needed to compute ! the constraints on the refined mesh ! ! FT 19.11.2021 ! !******************************************************* IMPLICIT NONE END PROCEDURE interpolate_id_hydro MODULE PROCEDURE interpolate_id_particles !**************************************************** ! !# Stores the hydro |id| in the arrays needed to ! compute the |sph| |id| ! ! FT 19.11.2020 ! !**************************************************** USE constants, ONLY: MSun, amu USE utility, ONLY: two USE numerics, ONLY: trilinear_interpolation IMPLICIT NONE LOGICAL, PARAMETER:: debug= .FALSE. ! INTEGER:: i, j, k DOUBLE PRECISION:: zp, xtmp, ytmp, ztmp CHARACTER( LEN= : ), ALLOCATABLE:: finalnamefile LOGICAL:: exist DOUBLE PRECISION:: foo(n), foo_exact(n), & foo_grid(this% nx_grid, this% ny_grid, this% nz_grid), & grid_coords(this%nx_grid,this%ny_grid,this%nz_grid,3), & coords(n,3) IF( debug )THEN DO i= 1, this% nx_grid, 1 DO j= 1, this% ny_grid, 1 DO k= 1, this% nz_grid, 1 grid_coords(i,j,k,1)= DBLE(i) - DBLE(this% nx_grid)/two grid_coords(i,j,k,2)= DBLE(j) - DBLE(this% ny_grid)/two grid_coords(i,j,k,3)= DBLE(k)/two! - DBLE(this% nz_grid)/two foo_grid(i,j,k)= (grid_coords(i,j,k,3))**3.D0 ENDDO ENDDO ENDDO foo= zero DO i= 1, n, 1 !CALL RANDOM_NUMBER( xsgn ) !CALL RANDOM_NUMBER( ysgn ) !CALL RANDOM_NUMBER( zsgn ) CALL RANDOM_NUMBER( xtmp ) CALL RANDOM_NUMBER( ytmp ) CALL RANDOM_NUMBER( ztmp ) coords(i,1)= xtmp*DBLE(this% nx_grid - 2) & - DBLE(this% nx_grid)/two + two coords(i,2)= ytmp*DBLE(this% ny_grid - 2) & - DBLE(this% ny_grid)/two + two coords(i,3)= (- DBLE(this% nz_grid)/two + one)*(one-ztmp) & + (DBLE(this% nz_grid)/two - one)*ztmp foo(i)= trilinear_interpolation( coords(i,1), coords(i,2), coords(i,3),& this% nx_grid, this% ny_grid, this% nz_grid, & grid_coords, foo_grid, & equator_symmetry= .FALSE., z_parity= -one, debug= .FALSE. ) foo_exact(i)= (coords(i,3))**3.D0 ENDDO ENDIF ! The density has to be converted in units of the atomic mass unit ! TODO: CHECK THAT EVERYTHING ELSE IS CONSISTENT WITH this!! DO i= 1, n, 1 zp= z(i) baryon_density(i) = this% read_mass_density( x(i), y(i), zp )*MSun/amu u_euler_x(i) = trilinear_interpolation( x(i), y(i), zp, & this% nx_grid, this% ny_grid, this% nz_grid, & this% grid, this% vel(:,:,:,1), & equator_symmetry= .TRUE., z_parity= one, & debug= .FALSE. ) u_euler_y(i) = trilinear_interpolation( x(i), y(i), zp, & this% nx_grid, this% ny_grid, this% nz_grid, & this% grid, this% vel(:,:,:,2), & equator_symmetry= .TRUE., z_parity= one, & debug= .FALSE. ) u_euler_z(i) = trilinear_interpolation( x(i), y(i), zp, & this% nx_grid, this% ny_grid, this% nz_grid, & this% grid, this% vel(:,:,:,3), & equator_symmetry= .TRUE., z_parity= -one, & debug= .FALSE. ) specific_energy(i)= trilinear_interpolation( x(i), y(i), zp, & this% nx_grid, this% ny_grid, this% nz_grid, & this% grid, this% specific_energy, & equator_symmetry= .TRUE., z_parity= one, & debug= .FALSE. ) IF( baryon_density(i) == zero )THEN specific_energy(i)= zero u_euler_x(i) = zero u_euler_y(i) = zero u_euler_z(i) = zero ENDIF ENDDO IF( debug )THEN finalnamefile= "dbg_interpolation2.dat" INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) IF( exist )THEN OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "REPLACE", & FORM= "FORMATTED", & POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, & IOMSG= err_msg ) ELSE OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "NEW", & FORM= "FORMATTED", & ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg ) ENDIF IF( ios > 0 )THEN PRINT *, "...error when opening " // TRIM(finalnamefile), & ". The error message is", err_msg STOP ENDIF DO i= 1, this% nx_grid - 1, 1 DO j= 1, this% ny_grid - 1, 1 DO k= 1, this% nz_grid - 1, 1 WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & this% grid( i, j, k, 1 ), & this% grid( i, j, k, 2 ), & this% grid( i, j, k, 3 ), & this% baryon_mass_density( i, j, k )*Msun/amu, & this% read_mass_density( & this% grid( i, j, k, 1 ) + this% dx_grid/two, & this% grid( i, j, k, 2 ), & this% grid( i, j, k, 3 ) ), & this% grid( i, j, k, 1 ) + this% dx_grid/two, & this% specific_energy( i, j, k ) ENDDO ENDDO ENDDO CLOSE( UNIT= 2 ) finalnamefile= "dbg_interpolation.dat" INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist ) IF( exist )THEN OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "REPLACE", & FORM= "FORMATTED", & POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, & IOMSG= err_msg ) ELSE OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "NEW", & FORM= "FORMATTED", & ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg ) ENDIF IF( ios > 0 )THEN PRINT *, "...error when opening " // TRIM(finalnamefile), & ". The error message is", err_msg STOP ENDIF DO i= 1, n, 1 ! IF( coords(i,3) < 0 )THEN ! PRINT *, coords(i,3) ! STOP ! ENDIF WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & i, x(i), y(i), z(i), & baryon_density(i), & u_euler_x(i), & u_euler_y(i), & u_euler_z(i), & specific_energy(i), & coords(i,1), coords(i,2), coords(i,3), & foo(i), foo_exact(i) ENDDO CLOSE( UNIT= 2 ) !STOP ENDIF energy_density = zero pressure = zero g_xx= one g_yy= one g_zz= one g_xy= zero g_xz= zero g_yz= zero lapse= one shift_x= zero shift_y= zero shift_z= zero END PROCEDURE interpolate_id_particles MODULE PROCEDURE interpolate_id_mass_b !**************************************************** ! !# Stores the hydro |id| in the arrays needed to ! compute the baryon mass, storing it to variables ! (not arrays as the others SUBROUTINES in ! the [[ejecta_generic@interpolate]] SUBMODULE). ! ! FT 19.11.2021 ! !**************************************************** USE tensor, ONLY: jxx, jxy, jxz, jyy, jyz, jzz USE constants, ONLY: MSun, amu IMPLICIT NONE baryon_density= this% read_mass_density( x, y, z ) g(jxx)= one g(jyy)= one g(jzz)= one g(jxy)= zero g(jxz)= zero g(jyz)= zero gamma_euler= one END PROCEDURE interpolate_id_mass_b MODULE PROCEDURE interpolate_id_k !**************************************************** ! !# Stores the components of the extrinsic curvature ! in arrays ! ! @warning DEPRECATED ! ! FT 19.11.2021 ! !**************************************************** IMPLICIT NONE END PROCEDURE interpolate_id_k !-----------------! !-- FUNCTIONS --! !-----------------! MODULE PROCEDURE interpolate_mass_density !*********************************************** ! !# Returns the mass density at the point ! given as argument, in units of ! \(M_\odot/L_\odot^3\). ! ! FT 19.11.2021 ! !*********************************************** !USE timing, ONLY: timer USE constants, ONLY: pi USE numerics, ONLY: trilinear_interpolation USE utility, ONLY: spherical_from_cartesian, two IMPLICIT NONE DOUBLE PRECISION:: zp, x_ell, y_ell, z_ell, theta, phi, r !TYPE(timer):: dbg_timer !dbg_timer= timer("dbg_timer") !CALL dbg_timer% start_timer() zp= z res= trilinear_interpolation( x, y, zp, & this% nx_grid, this% ny_grid, this% nz_grid, & this% grid, this% baryon_mass_density, & equator_symmetry= .TRUE., z_parity= one, & debug= .FALSE. ) CALL spherical_from_cartesian( x, y, z, & this% centers(1,1), this% centers(1,2), this% centers(1,3), & r, theta, phi ) x_ell= this% centers(1,1) & + MAX(this% sizes(1,1),this% sizes(1,2))*COS(phi)*SIN(theta) y_ell= this% centers(1,2) & + MAX(this% sizes(1,3),this% sizes(1,4))*SIN(phi)*SIN(theta) z_ell= this% centers(1,3) & + MAX(this% sizes(1,5),this% sizes(1,6))*COS(theta) IF( r >= SQRT( ( x_ell - this% centers(1,1) )**two & + ( y_ell - this% centers(1,2) )**two & + ( z_ell - this% centers(1,3) )**two ) ) res= zero IF( res < zero ) res= zero ! IF( x > this% centers(1,1) + this% sizes(1,2) & ! .OR. x < this% centers(1,1) - this% sizes(1,1) & ! .OR. y > this% centers(1,2) + this% sizes(1,4) & ! .OR. y < this% centers(1,2) - this% sizes(1,3) & ! .OR. zp > this% centers(1,3) + this% sizes(1,6) ) res= zero ! IF( x > this% xR_grid & ! .OR. x < this% xL_grid & ! .OR. y > this% yR_grid & ! .OR. y < this% yL_grid & ! .OR. zp > this% zR_grid ) res= zero !CALL dbg_timer% stop_timer() !CALL dbg_timer% print_timer( 2 ) !STOP END PROCEDURE interpolate_mass_density MODULE PROCEDURE interpolate_pressure !*********************************************** ! !# Returns the pressure at the point ! given as argument, in units of ! \(M_\odot c^2/L_\odot^3\). ! ! FT 19.11.2021 ! !*********************************************** IMPLICIT NONE ! TODO: There is no pressure in the ejecta ID. It should be computed using ! the |eos| and the specific energy (it's a hot system) PRINT *, "** The method interpolate_pressure of TYPE ejecta_generic" PRINT *, " has not been implemented yet." PRINT * STOP END PROCEDURE interpolate_pressure MODULE PROCEDURE interpolate_spatial_metric !*********************************************** ! !# Returns the spatial metric. ! ! FT 19.11.2021 ! !*********************************************** IMPLICIT NONE END PROCEDURE interpolate_spatial_metric MODULE PROCEDURE is_hydro_positive !************************************************ ! !# Return 1 if the energy density is nonpositive ! or if the specific energy is nonpositive, ! or if the pressure is nonpositive ! at the specified point; return 0 otherwise ! ! FT 19.11.2021 ! !************************************************ IMPLICIT NONE DOUBLE PRECISION, DIMENSION(3):: center DOUBLE PRECISION, DIMENSION(6):: sizes center= this% return_center(1) sizes = this% return_spatial_extent(1) IF( this% read_mass_density( x, y, z ) <= zero & .OR. & !SQRT( ( x - center(1) )**2 + ( y - center(2) )**2 & ! + ( z - center(3) )**2 ) > 50zero x > center(1) + sizes(1) & .OR. x < center(1) - sizes(2) & .OR. y > center(2) + sizes(3) & .OR. y < center(2) - sizes(4) & .OR. ABS(z) > center(3) + sizes(5) & )THEN res= .FALSE. ELSE res= .TRUE. ENDIF END PROCEDURE is_hydro_positive END SUBMODULE interpolate