! File: submodule_bnsfuka_constructor.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 (bns_fuka) constructor !********************************************************* ! !# Implementation of the constructor and ! destructor of TYPE [[bnsfuka]], and of the ! [[bnsfuka]]-member ! PROCEDURES that call the C-bound PROCEDURES ! constructig and destructing the |fuka| ! |bnsexp| object ! ! FT 23.10.2020 ! !********************************************************* IMPLICIT NONE CONTAINS ! !-- Implementation of the constructor of the bns object ! MODULE PROCEDURE construct_bnsfuka !**************************************************** ! !# Constructs an object of TYPE [[bnsfuka]] ! ! Created: FT 09.02.2022 ! Last update: FT 27.05.2022 ! !**************************************************** USE utility, ONLY: ten, Msun_geo, flag$tpo IMPLICIT NONE INTEGER, SAVE:: bns_counter= 1 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: length_scale_pressure IF(max$tpo < flag$tpo)THEN PRINT *, "** ERROR in MODULE bns_fuka! It must hold max$tpo > flag$tpo,",& "where both variables are negative integers. Instead:" PRINT *, " * max$tpo=", max$tpo PRINT *, " * flag$tpo=", flag$tpo PRINT *, " * Stopping..." PRINT * STOP ENDIF CALL derived_type% set_n_matter(2) CALL derived_type% set_cold_system(.TRUE.) derived_type% construction_timer= timer("binary_construction_timer") ! Construct |fuka| |bnsexp| object CALL derived_type% construct_binary(filename) derived_type% filename= filename ! Read the properties of the binary system CALL read_bns_properties(derived_type) ! Assign a unique identifier to the bns_fuka object derived_type% bns_identifier= bns_counter bns_counter= bns_counter + 1 ! Do not use the geodesic gauge by default CALL derived_type% set_one_lapse (.FALSE.) CALL derived_type% set_zero_shift(.FALSE.) ! Compute typical length scales of the system using the pressure IF( derived_type% get_estimate_length_scale() )THEN ALLOCATE( length_scale_pressure(derived_type% get_n_matter()) ) length_scale_pressure= derived_type% estimate_lengthscale_field( & get_pressure, & derived_type% get_n_matter() ) PRINT *, " * Minimum length scale to resolve on star 1, based on ", & "pressure= ", length_scale_pressure(1)*Msun_geo*ten*ten*ten, "m" PRINT *, " * Minimum length scale to resolve on star 2, based on ", & "pressure= ", length_scale_pressure(1)*Msun_geo*ten*ten*ten, "m" PRINT * ENDIF ! Assign PROCEDURE POINTER to the desired PROCEDURE derived_type% finalize_sph_id_ptr => finalize CONTAINS FUNCTION get_pressure( x, y, z ) RESULT( val ) !! Returns the value of the pressure at the desired point DOUBLE PRECISION, INTENT(IN):: x !! \(x\) coordinate of the desired point DOUBLE PRECISION, INTENT(IN):: y !! \(y\) coordinate of the desired point DOUBLE PRECISION, INTENT(IN):: z !! \(z\) coordinate of the desired point DOUBLE PRECISION:: val !! Pressure at \((x,y,z)\) val= derived_type% read_pressure( x, y, z ) END FUNCTION get_pressure END PROCEDURE construct_bnsfuka MODULE PROCEDURE finalize !*********************************************** ! !# This SUBROUTINE is curretly just a placeholder. ! It could be used, for exmaple, to correct ! the ADM momentum at the end of the execution, ! or correct for residual eccentricity, etc. ! ! @note Temporary implementation, to avoid warnings ! about unused variables. ! ! FT 14.04.2022 ! !*********************************************** IMPLICIT NONE pos = pos nlrf = nlrf nu = nu pr = pr vel_u= vel_u theta= theta nstar= nstar u = u END PROCEDURE finalize MODULE PROCEDURE initialize_id_bnsfuka !*********************************************** ! !# Initialize the |fuka| |bns| |id|. ! ! - If `flag`= [[utility:flag$sph]], set up the ! lattices around the stars for the |bns| ! produced with |fuka|. ! - If `flag`= [[utility:flag$tpo]], allocate ! memory for the hydro grid functions. ! - If `flag` > 0, assign its value to [[bnsfuka:l_curr]]. ! - If [[utility:flag$tpo]] < `flag` < 0, ! assign its value to [[bnsfuka:tpo_curr]]. ! ! FT 16.09.2022 ! !*********************************************** USE utility, ONLY: flag$sph, flag$tpo USE mesh_refinement, ONLY: allocate_grid_function IMPLICIT NONE INTEGER, SAVE:: tpo_counter= 1 !! Counts how many times the PROCEDURE construct_particles_idase is called INTEGER:: i CHARACTER(LEN= 3):: cnt_i CHARACTER(LEN= :), ALLOCATABLE:: name_mass_density CHARACTER(LEN= :), ALLOCATABLE:: name_specific_energy CHARACTER(LEN= :), ALLOCATABLE:: name_pressure CHARACTER(LEN= :), ALLOCATABLE:: name_v_euler_x CHARACTER(LEN= :), ALLOCATABLE:: name_v_euler_y CHARACTER(LEN= :), ALLOCATABLE:: name_v_euler_z LOGICAL:: wanted_tpo IF( PRESENT(switch) )THEN IF( switch .EQV. .TRUE. )THEN this% tpo$log(tpo_counter)= flag this% tpo_curr= tpo_counter tpo_counter= tpo_counter + 1 CALL initialize_id_bnsfuka(this, flag$tpo) ENDIF ENDIF IF( flag /= flag$sph .AND. flag /= flag$tpo .AND. flag < -max$tpo )THEN PRINT *, "** ERROR in SUBROUTINE initialize_id_bnsfuka! The INTEGER ", & "argument 'flag' should be in the set [1,nlevels], or ", & "it should be equal to either flag$sph or flag$tpo, defined ", & "in MODULE utility." PRINT *, " * Stopping..." PRINT * STOP ENDIF IF( flag == flag$sph )THEN ! Since Kadath is not thread-safe, we cannot parallelize it using OMP ! within SPHINCS_ID. Hence, we chose to make a system call to a program ! within Kadath that reads the ID from the FUKA output file and prints it ! on a lattice. The ID on the particles will be interplated from this fine ! lattice. CALL this% set_up_lattices_around_stars() IF(ALLOCATED(this% surfaces)) DEALLOCATE(this% surfaces) !PRINT *, ALLOCATED(this% surfaces(1)% points) !PRINT *, ALLOCATED(this% surfaces(2)% points) ! Find the surfaces of the stars and print them to a formatted file CALL this% find_print_surfaces() ELSEIF( flag >= flag$tpo .AND. flag <= -1 )THEN wanted_tpo= .FALSE. DO i= 1, max$tpo, 1 IF(this% tpo$log(i) == flag)THEN wanted_tpo= .TRUE. this% tpo_curr= i ENDIF ENDDO IF(flag == flag$tpo)THEN IF( this% tpo_curr <= 9 ) WRITE( cnt_i, "(I1)" ) this% tpo_curr IF( this% tpo_curr >= 10 .AND. flag <= 99 ) WRITE( cnt_i, "(I2)" ) & this% tpo_curr IF( this% tpo_curr >= 100 .AND. flag <= 999 ) WRITE( cnt_i, "(I3)" ) & this% tpo_curr !PRINT * !PRINT *, "cnt_i=", TRIM(cnt_i) !PRINT * name_mass_density = "mass_density_fuka-"//TRIM(cnt_i) name_specific_energy= "specific_energy_fuka-"//TRIM(cnt_i) name_pressure = "pressure_fuka-"//TRIM(cnt_i) name_v_euler_x = "v_euler_x_fuka-"//TRIM(cnt_i) name_v_euler_y = "v_euler_y_fuka-"//TRIM(cnt_i) name_v_euler_z = "v_euler_z_fuka-"//TRIM(cnt_i) !PRINT * !PRINT *, "name_mass_density=", name_mass_density !PRINT * CALL allocate_grid_function( this% mass_density(this% tpo_curr), & TRIM(name_mass_density), 1 ) CALL allocate_grid_function( this% specific_energy(this% tpo_curr), & TRIM(name_specific_energy), 1 ) CALL allocate_grid_function( this% pressure(this% tpo_curr), & TRIM(name_pressure), 1 ) CALL allocate_grid_function( this% v_euler_x(this% tpo_curr), & TRIM(name_v_euler_x), 1 ) CALL allocate_grid_function( this% v_euler_y(this% tpo_curr), & TRIM(name_v_euler_y), 1 ) CALL allocate_grid_function( this% v_euler_z(this% tpo_curr), & TRIM(name_v_euler_z), 1 ) ELSE IF(.NOT.wanted_tpo)THEN PRINT *, "** ERROR! Mismatch between bns_fuka and bssn objects!" PRINT *, " * This should never happen: there is most likely a bug ", & "in SUBROUTINE initialize_id_bnsfuka, or a bug in the ", & "places where it is called." PRINT *, " * Stopping..." PRINT * STOP ENDIF ENDIF ELSE this% l_curr= flag ENDIF END PROCEDURE initialize_id_bnsfuka ! !-- Implementation of the destructor of the bns object ! MODULE PROCEDURE destruct_bnsfuka !*********************************************** ! !# Destructs an object of TYPE [[bnsfuka]] ! ! Created: FT 09.02.2022 ! Last update: FT 30.06.2022 ! !*********************************************** IMPLICIT NONE INTEGER:: i_star ! Deallocate memory CALL this% deallocate_bnsfuka_memory() DO i_star=1, 2, 1 CALL this% star_lattice(i_star)% deallocate_lattice_memory() ENDDO END PROCEDURE destruct_bnsfuka MODULE PROCEDURE construct_binary !*********************************************** ! !# Construct the |fuka| |bnsexp| object ! ! Created: FT 09.02.2022 ! Last update: FT 27.05.2022 ! !*********************************************** IMPLICIT NONE LOGICAL:: exist #ifdef __INTEL_COMPILER IF ( C_ASSOCIATED( this% bns_ptr ) ) THEN CALL destruct_bns_fuka( this% bns_ptr ) ENDIF #endif INQUIRE( FILE= fukafile, EXIST= exist ) IF( exist )THEN CALL this% construction_timer% start_timer() this% bns_ptr = construct_bns_fuka( fukafile//C_NULL_CHAR ) CALL this% construction_timer% stop_timer() ELSE PRINT *, "** ERROR in SUBROUTINE construct_binary: file ", & fukafile, " cannot be found!" PRINT * STOP ENDIF END PROCEDURE construct_binary MODULE PROCEDURE destruct_binary !************************************************ ! !# Destructs the |fuka| |bnsexp| object and frees ! the pointer [[bnsfuka:bns_ptr]] pointing to it ! ! Created: FT 09.02.2022 ! Last update: FT 27.05.2022 ! !************************************************ IMPLICIT NONE IF ( C_ASSOCIATED( this% bns_ptr ) ) THEN CALL destruct_bns_fuka( this% bns_ptr ) this% bns_ptr = C_NULL_PTR ENDIF END PROCEDURE destruct_binary END SUBMODULE constructor