! File: submodule_sph_particles_quality_indicators.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) quality_indicators !**************************************************** ! !# This SUBMODULE contains the implementation of ! the methods of TYPE [[particles]] ! that computes the quality indicators, referring to ! ! [Daniel J. Price, Smoothed Particle Hydrodynamics and Magnetohydrodynamics. Journal of Computational Physics, 231, 3, 759-794 (2012). DOI: 10.1016/j.jcp.2010.12.011](http://arxiv.org/abs/1012.1885){:target="_blank"}, ! eqs.(64), (67) and (74-75) ! ! [Rosswog, S. SPH Methods in the Modelling of Compact Objects. Living Rev Comput Astrophys 1, 1 (2015).](https://arxiv.org/abs/1406.4224){:target="_blank"}, ! eqs.(6) and (9) ! ! FT 05.10.2022 ! !**************************************************** IMPLICIT NONE CONTAINS !-------------------! !-- SUBROUTINES --! !-------------------! MODULE PROCEDURE compute_and_print_quality_indicators !**************************************************** ! !# Compute the quality indicators, referring to ! ! [Daniel J. Price, Smoothed Particle Hydrodynamics and Magnetohydrodynamics. Journal of Computational Physics, 231, 3, 759-794 (2012). DOI: 10.1016/j.jcp.2010.12.011](http://arxiv.org/abs/1012.1885){:target="_blank"}, ! eqs.(64), (67) and (74-75) ! ! [Rosswog, S. SPH Methods in the Modelling of Compact Objects. Living Rev Comput Astrophys 1, 1 (2015).](https://arxiv.org/abs/1406.4224){:target="_blank"}, ! eqs.(6) and (9) ! ! FT 05.10.2022 ! !**************************************************** USE kernel_table, ONLY: dv_table, dv_table_1,& W_no_norm, n_tab_entry, &!dWdv_no_norm, & interp_gradW_table,interp_W_gradW_table USE sphincs_sph, ONLY: ncand, all_clists USE RCB_tree_3D, ONLY: iorig, nic, nfinal, nprev, lpart, rpart USE options, ONLY: ndes USE set_h, ONLY: exact_nei_tree_update USE utility, ONLY: zero, one, two IMPLICIT NONE DOUBLE PRECISION, DIMENSION(npart) :: qi_1 DOUBLE PRECISION, DIMENSION(3,npart) :: qi_2 DOUBLE PRECISION, DIMENSION(3,npart) :: qi_3 DOUBLE PRECISION, DIMENSION(3,3,npart):: qi_4 DOUBLE PRECISION, DIMENSION(3,npart) :: qi_ham INTEGER:: a, b, ill, l, itot, inde, index1, k DOUBLE PRECISION:: ha, ha_1, ha_3, ha_4, va, xa, ya, za, & hb, hb_1, hb_3, hb_4, xb, yb, zb, rab, rab2, rab_1, & ha2_4, hb2_4, dx, dy, dz DOUBLE PRECISION:: Wab_ha, Wi, Wi1, dvv, & grW_ha_x, grW_ha_y, grW_ha_z, & grW_hb_x, grW_hb_y, grW_hb_z, eab(3), & Wa, grW, grWa, grWb, vb, vol_b LOGICAL, PARAMETER:: debug= .TRUE. INTEGER :: ios, unit_qi LOGICAL :: exist CHARACTER(LEN=:), ALLOCATABLE:: namefile, err_msg PRINT *, " * Computing the quality indicators..." ! exact_nei_tree_update updates h CALL exact_nei_tree_update( ndes, npart, pos, nu ) qi_1 = zero qi_2 = zero qi_3 = zero qi_4 = zero qi_ham= zero !$OMP PARALLEL DO DEFAULT( NONE ) & !$OMP SHARED( nfinal, nprev, iorig, lpart, rpart, nic,nu,nstar, & !$OMP ncand, all_clists, W_no_norm, pos, h, & !$OMP qi_1, qi_2, qi_3, qi_4, qi_ham ) & !$OMP PRIVATE( ill, itot, a, b, l, & !$OMP ha, ha_1, ha_3, ha_4, ha2_4, & !$OMP hb, hb_1, hb_3, hb_4, hb2_4, & !$OMP xa, ya, za, xb, yb, zb, dx, dy, dz, rab2, & !$OMP inde, va, dv_table_1, index1, Wi, & !$OMP dvv, dv_table, Wab_ha, Wa, grW, Wi1, & !$OMP grW_ha_x, grW_ha_y, grW_ha_z, grWa, grWb, & !$OMP grW_hb_x, grW_hb_y, grW_hb_z, eab, rab, vb, & !$OMP rab_1, vol_b ) lowest_level_cell_loop: DO ill= 1, nfinal itot= nprev + ill IF( nic(itot) == 0 ) CYCLE particle_loop: DO l= lpart(itot), rpart(itot) ! Set index of the particle in the lowest level cell a = iorig(l) ! Compute smoothing lengths and its powers (used below) ha = h(a) ha_1 = one/ha ha_3 = ha_1*ha_1*ha_1 ha_4 = ha_3*ha_1 ha2_4= two*two*ha*ha xa= pos(1,a) ya= pos(2,a) za= pos(3,a) candidate_loop: DO k= 1, ncand(ill) ! Set index of candidate particle b= all_clists(ill)% list(k) ! Compute smoothing lengths and its powers (used below) hb = h(b) hb_1 = one/hb hb_3 = hb_1*hb_1*hb_1 hb_4 = hb_3*hb_1 hb2_4= two*two*hb*hb xb= pos(1,b) yb= pos(2,b) zb= pos(3,b) ! Compute distance between particles a and b dx= xa - xb dy= ya - yb dz= za - zb rab2 = dx*dx + dy*dy + dz*dz rab = SQRT(rab2) + 1.D-30 rab_1= one/rab ! If the distance between the particles a and b is larger than twice ! their smoothing lengths, move to next candidate particle ! This one is not a neighbour IF( rab2 > ha2_4 .AND. rab2 > hb2_4 ) CYCLE va= rab*ha_1 ! Read kernel values inde = MIN(INT(va*dv_table_1),n_tab_entry) index1= MIN(inde + 1,n_tab_entry) Wi = W_no_norm(inde) Wi1= W_no_norm(index1) dvv = (va - DBLE(inde)*dv_table)*dv_table_1 Wab_ha= Wi + (Wi1 - Wi)*dvv ! Compute normalized distance between the particles a and b eab(1)= dx*rab_1 eab(2)= dy*rab_1 eab(3)= dz*rab_1 ! !-- Interpolate kernel value, and its gradient, !-- at the distance between the particles ! !DIR$ INLINE CALL interp_W_gradW_table( va, Wa, grW ) Wa = Wa*ha_3 grW= grW*ha_4 grW_ha_x= grW*eab(1) grW_ha_y= grW*eab(2) grW_ha_z= grW*eab(3) vb= rab*hb_1 !DIR$ INLINE CALL interp_gradW_table(vb,grW) grW= grW*hb_4 grW_hb_x= grW*eab(1) grW_hb_y= grW*eab(2) grW_hb_z= grW*eab(3) ! Compute particle volume vol_b= nu(b)/nstar(b) ! !-- Compute quality indicators ! qi_1(a)= qi_1(a) + Wab_ha*vol_b qi_2(1,a)= qi_2(1,a) + (-dx)*Wab_ha*vol_b qi_2(2,a)= qi_2(2,a) + (-dy)*Wab_ha*vol_b qi_2(3,a)= qi_2(3,a) + (-dz)*Wab_ha*vol_b qi_3(1,a)= qi_3(1,a) + grW_ha_x*vol_b qi_3(2,a)= qi_3(2,a) + grW_ha_y*vol_b qi_3(3,a)= qi_3(3,a) + grW_ha_z*vol_b qi_4(1,1,a)= qi_4(1,1,a) + (-dx)*grW_ha_x*vol_b qi_4(1,2,a)= qi_4(1,2,a) + (-dx)*grW_ha_y*vol_b qi_4(1,3,a)= qi_4(1,3,a) + (-dx)*grW_ha_z*vol_b qi_4(2,1,a)= qi_4(2,1,a) + (-dy)*grW_ha_x*vol_b qi_4(2,2,a)= qi_4(2,2,a) + (-dy)*grW_ha_y*vol_b qi_4(2,3,a)= qi_4(2,3,a) + (-dy)*grW_ha_z*vol_b qi_4(3,1,a)= qi_4(3,1,a) + (-dz)*grW_ha_x*vol_b qi_4(3,2,a)= qi_4(3,2,a) + (-dz)*grW_ha_y*vol_b qi_4(3,3,a)= qi_4(3,3,a) + (-dz)*grW_ha_z*vol_b qi_ham(1,a)= qi_ham(1,a) & + nu(b)*grW_ha_x*(one/nstar(a) + nstar(a)/nstar(b)**2) qi_ham(2,a)= qi_ham(2,a) & + nu(b)*grW_ha_y*(one/nstar(a) + nstar(a)/nstar(b)**2) qi_ham(3,a)= qi_ham(3,a) & + nu(b)*grW_ha_z*(one/nstar(a) + nstar(a)/nstar(b)**2) ENDDO candidate_loop !qi(a)= qi(a)/(nu(a)/nstar(a)) ENDDO particle_loop ENDDO lowest_level_cell_loop !$OMP END PARALLEL DO PRINT *, " * Printing the quality indicators to file..." IF(PRESENT(path))THEN namefile= path//"quality_indicators.dat" ELSE namefile= "quality_indicators.dat" ENDIF unit_qi= 279465 INQUIRE( FILE= TRIM(namefile), EXIST= exist ) IF( exist )THEN OPEN( UNIT= unit_qi, FILE= TRIM(namefile), STATUS= "REPLACE", & FORM= "FORMATTED", & POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, & IOMSG= err_msg ) ELSE OPEN( UNIT= unit_qi, FILE= TRIM(namefile), STATUS= "NEW", & FORM= "FORMATTED", & ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg ) ENDIF IF( ios > 0 )THEN PRINT *, "...error when opening " // TRIM(namefile), & ". The error message is", err_msg STOP ENDIF WRITE( UNIT = unit_qi, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & "# Run ID [ccyymmdd-hhmmss.sss]: " // run_id WRITE( UNIT = unit_qi, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & "# Quality indicators on the particles" IF( ios > 0 )THEN PRINT *, "...error when writing line 1 in " // TRIM(namefile), & ". The error message is", err_msg STOP ENDIF WRITE( UNIT = unit_qi, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & "# column: 1 2 3 4 5", & " 6 7 8", & " 9 10 11", & " 12 13 14", & " 15 16 17", & " 18 19 20", & " 21 22 23" IF( ios > 0 )THEN PRINT *, "...error when writing line 2 in " // TRIM(namefile), & ". The error message is", err_msg STOP ENDIF WRITE( UNIT = unit_qi, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & "# particle index x [Msun_geo] ", & " y [Msun_geo] z [Msun_geo] Q1", & " Q2_x Q2_y Q2_z", & " Q3_x Q3_y Q3_z", & " Q4_xx Q4_xy Q4_xz", & " Q4_yx Q4_yy Q4_yz", & " Q4_zx Q4_zy Q4_zz", & " Qhamilonian_x Qhamilonian_y Qhamilonian_z" IF( ios > 0 )THEN PRINT *, "...error when writing line 3 in " // TRIM(namefile), & ". The error message is", err_msg STOP ENDIF DO a= 1, npart, 1 WRITE( UNIT = unit_qi, IOSTAT = ios, IOMSG = err_msg, FMT = * ) & a, & ! 1 pos(1,a), & ! 2 pos(2,a), & ! 3 pos(3,a), & ! 4 qi_1(a), & ! 5 qi_2(1,a), & ! 6 qi_2(2,a), & ! 7 qi_2(3,a), & ! 8 qi_3(1,a), & ! 9 qi_3(2,a), & ! 10 qi_3(3,a), & ! 11 qi_4(1,1,a), & ! 12 qi_4(1,2,a), & ! 13 qi_4(1,3,a), & ! 14 qi_4(2,1,a), & ! 15 qi_4(2,2,a), & ! 16 qi_4(2,3,a), & ! 17 qi_4(3,1,a), & ! 18 qi_4(3,2,a), & ! 19 qi_4(3,3,a), & ! 20 qi_ham(1,a), & ! 21 qi_ham(2,a), & ! 22 qi_ham(3,a) ! 23 ENDDO CLOSE( UNIT= unit_qi ) PRINT *, " quality indicators printed to file ", namefile PRINT * END PROCEDURE compute_and_print_quality_indicators END SUBMODULE quality_indicators