From c5b2514255faa0e743d0061e47fbe88d9aab15b8 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 21 May 2024 14:21:05 +1000 Subject: [PATCH 01/13] Adding ability to read unstructured densities from fits file --- src/density.f90 | 146 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 103 insertions(+), 43 deletions(-) diff --git a/src/density.f90 b/src/density.f90 index ed0c99579..be61fab84 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -970,10 +970,12 @@ subroutine read_density_file() real, dimension(:,:,:), allocatable :: sph_gas_dens ! (n_rad,nz,n_az) real, dimension(:,:,:,:), allocatable :: sph_V ! (n_rad,nz,n_az,3) real, dimension(:), allocatable :: a_sph, n_a_sph, log_a_sph, log_n_a_sph ! n_a + real, dimension(:,:), allocatable :: xyz_sph real(kind=dp), dimension(:,:,:,:), allocatable :: sph_dens_dp, sph_V_dp real(kind=dp), dimension(:,:,:), allocatable :: sph_gas_dens_dp real(kind=dp), dimension(:), allocatable :: a_sph_dp + real(kind=dp), dimension(:,:), allocatable :: xyz_sph_dp real(kind=dp) :: f logical :: lread_gas_density, lread_gas_velocity @@ -1114,55 +1116,66 @@ subroutine read_density_file() ! determine the size of density file write(*,*) "Reading dust density ..." call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) - if ((nfound /= 3) .and. (nfound /= 4)) then - write(*,*) "I found", nfound, "axis instead of 3 or 4" + if ((nfound /= 1) .and. (nfound /= 3) .and. (nfound /= 4)) then + write(*,*) "I found", nfound, "axis instead of 1, 3 or 4" call error('failed to read the NAXIS keyword in HDU 1 of '//trim(density_file)//' file') endif - if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then - write(*,*) "# fits_file vs mcfost_grid" - write(*,*) naxes(1), n_rad - write(*,*) naxes(2), nz - write(*,*) naxes(3), n_az - !write(*,*) naxes(4), n_a - call error(trim(density_file)//" does not have the right dimensions in HDU 1.") - endif + if (nfound == 1) then ! Voronoi mesh + write(*,*) "Found 1D density structure, using a Voronoi mesh" + lVoronoi = .true. + l3D = .true. + n_cells = naxes(1) + npixels = naxes(1) + n_a = 1 ! for now - if (nfound == 3) then - n_a = 1 - write(*,*) "No grain size found" - npixels=naxes(1)*naxes(2)*naxes(3) - else - if (.not.lvariable_dust) then - call warning("Forcing variable dust") - lvariable_dust = .true. + allocate(sph_dens(n_cells,1,1,n_a)) + else ! Structured grid + if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then + write(*,*) "# fits_file vs mcfost_grid" + write(*,*) naxes(1), n_rad + write(*,*) naxes(2), nz + write(*,*) naxes(3), n_az + !write(*,*) naxes(4), n_a + call error(trim(density_file)//" does not have the right dimensions in HDU 1.") endif - n_a = naxes(4) - write(*,*) n_a, "grain sizes found" - npixels=naxes(1)*naxes(2)*naxes(3)*naxes(4) - endif - if (lvariable_dust) then - p_n_cells = n_cells - p_n_rad=n_rad ; p_nz = nz - else - p_n_cells = 1 - p_n_rad=1 ; p_nz=1 - endif + if (nfound == 3) then + n_a = 1 + write(*,*) "No grain size found" + npixels=naxes(1)*naxes(2)*naxes(3) + else + if (.not.lvariable_dust) then + call warning("Forcing variable dust") + lvariable_dust = .true. + endif + n_a = naxes(4) + write(*,*) n_a, "grain sizes found" + npixels=naxes(1)*naxes(2)*naxes(3)*naxes(4) + endif - if (naxes(2) == 2*nz) then - l3D_file = .true. - else - write(*,*) "The density file only has positive z, making it symmetric" - l3D_file = .false. - endif + if (lvariable_dust) then + p_n_cells = n_cells + p_n_rad=n_rad ; p_nz = nz + else + p_n_cells = 1 + p_n_rad=1 ; p_nz=1 + endif - if (l3D_file) then - allocate(sph_dens(n_rad,2*nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) - else - allocate(sph_dens(n_rad,nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) + if (naxes(2) == 2*nz) then + l3D_file = .true. + else + write(*,*) "The density file only has positive z, making it symmetric" + l3D_file = .false. + endif + + if (l3D_file) then + allocate(sph_dens(n_rad,2*nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) + else + allocate(sph_dens(n_rad,nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) + endif + sph_dens = 0.0 ; a_sph = 0.0 ; n_a_sph = 0.0 endif - sph_dens = 0.0 ; a_sph = 0.0 ; n_a_sph = 0.0 bitpix = 0 call ftgkyj(unit,"bitpix",bitpix,comment,status) @@ -1171,10 +1184,14 @@ subroutine read_density_file() if (bitpix==-32) then call ftgpve(unit,group,firstpix,npixels,nullval,sph_dens,anynull,status) else if (bitpix==-64) then - if (l3D_file) then - allocate(sph_dens_dp(n_rad,2*nz,n_az,n_a)) + if (lVoronoi) then + allocate(sph_dens_dp(n_cells,1,1,n_a)) else - allocate(sph_dens_dp(n_rad,nz,n_az,n_a)) + if (l3D_file) then + allocate(sph_dens_dp(n_rad,2*nz,n_az,n_a)) + else + allocate(sph_dens_dp(n_rad,nz,n_az,n_a)) + endif endif sph_dens_dp = 0.0_dp call ftgpvd(unit,group,firstpix,npixels,nullval,sph_dens_dp,anynull,status) @@ -1190,6 +1207,49 @@ subroutine read_density_file() sph_dens = sph_dens/maxval(sph_dens) ! normalization avant d'ajouter une constante sph_dens = max(sph_dens,1e10*tiny_real) + + if (lVoronoi) then + !--------------------------------------------------------- + ! hdu2 is particle positions + !--------------------------------------------------------- + ! move to next hdu + call ftmrhd(unit,1,hdutype,status) + nfound=1 + ! Check dimensions + naxes(:) = 0 + call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) + if (nfound /= 2) then + write(*,*) 'HDU 2 has', nfound, 'dimensions.' + call error('did not find 2 dimension in HDU 2') + endif + if ((naxes(1) /= 3)) then + write(*,*) "HDU2 dim 1 is ", naxes(1), "instead of 3" + call error("HDU 2 does not have the right dimension") + endif + if ((naxes(2) /= n_cells)) then + write(*,*) "HDU2 dim 2 is ", naxes(3), "instead of ",n_cells + call error("HDU 2 does not have the right dimension") + endif + + npixels=naxes(1)*naxes(2) + + bitpix=0 + call ftgkyj(unit,"bitpix",bitpix,comment,status) + + ! read_image + if (bitpix==-32) then + allocate(xyz_sph(3,n_cells),xyz_sph_dp(3,n_cells)) + call ftgpve(unit,group,firstpix,npixels,nullval,xyz_sph,anynull,status) + xyz_sph_dp = real(xyz_sph,kind=dp) + deallocate(xyz_sph) + else if (bitpix==-64) then + allocate(xyz_sph_dp(3,n_cells)) + call ftgpvd(unit,group,firstpix,npixels,nullval,xyz_sph_dp,anynull,status) + else + call error("cannot read bitpix in fits file") + endif + endif + if (n_a > 1) then write(*,*) "Reading grain sizes ..." read_n_a = 0 From 50f989a965f727481643b84fea62a1ab0c3f22c6 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 22 May 2024 09:42:11 +1000 Subject: [PATCH 02/13] Checking first if density file is unstructured --- src/density.f90 | 43 ++++++++++++++++++++++++++++++++++++++----- src/dust_transfer.f90 | 22 +++++++++++++--------- 2 files changed, 51 insertions(+), 14 deletions(-) diff --git a/src/density.f90 b/src/density.f90 index be61fab84..806abeedc 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -16,13 +16,11 @@ module density integer, public :: species_removed real, public :: T_rm - public :: densite_gaz, masse_gaz, surface_density, densite_gaz_midplane, densite_pouss, masse, icell_not_empty - - public :: define_density, define_density_wall3d, define_dust_density, read_density_file, & + public :: densite_gaz, masse_gaz, surface_density, densite_gaz_midplane, densite_pouss, masse, icell_not_empty, & + define_density, define_density_wall3d, define_dust_density, read_density_file, is_density_file_Voronoi, & densite_seb_charnoz2, densite_seb_charnoz, remove_species, read_sigma_file, normalize_dust_density, & reduce_density - private real(kind=dp), dimension(:), allocatable :: densite_gaz, masse_gaz ! n_rad, nz, n_az, Unites: part.m-3 et g : H2 @@ -940,6 +938,40 @@ end subroutine define_density_wall3D !******************************************************************** +subroutine is_density_file_Voronoi() + + integer :: unit, status, nfound, readwrite, blocksize + integer, dimension(4) :: naxes + + status=0 + ! Get an unused Logical Unit Number to use to open the FITS file. + call ftgiou(unit,status) + + ! Opening file + readwrite=0 + call ftopen(unit,density_file,readwrite,blocksize,status) + if (status /= 0) call error("density file needed") + + nfound=0 + call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) + ! Closing file + call ftclos(unit, status) + call ftfiou(unit, status) + + if (status /= 0) call error("error reading density file") + + if (nfound == 1) then + lVoronoi = .true. + else + lVoronoi = .false. + endif + + return + +end subroutine is_density_file_Voronoi + +!******************************************************************** + subroutine read_density_file() ! Nouvelle routine pour lire les grilles de densite ! calculees par Yorick directement a partir des donnees SPH (ou autre) @@ -1116,6 +1148,7 @@ subroutine read_density_file() ! determine the size of density file write(*,*) "Reading dust density ..." call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) + if ((nfound /= 1) .and. (nfound /= 3) .and. (nfound /= 4)) then write(*,*) "I found", nfound, "axis instead of 1, 3 or 4" call error('failed to read the NAXIS keyword in HDU 1 of '//trim(density_file)//' file') @@ -1129,7 +1162,7 @@ subroutine read_density_file() npixels = naxes(1) n_a = 1 ! for now - allocate(sph_dens(n_cells,1,1,n_a)) + allocate(sph_dens(n_cells,1,1,n_a), a_sph(n_a), n_a_sph(n_a)) else ! Structured grid if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then write(*,*) "# fits_file vs mcfost_grid" diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index 2eee4a264..a2bb7edfe 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -107,21 +107,25 @@ subroutine transfert_poussiere() ! Building the dust grain population call build_grain_size_distribution() - if (lphantom_file .or. lgadget2_file .or. lascii_SPH_file) then - call setup_SPH2mcfost(density_file, limits_file, n_SPH, extra_heating) - call setup_grid() - else if (lmhd_voronoi) then - call setup_mhd_to_mcfost() !uses sph_to_voronoi + laffichage=.true. + + if (ldensity_file) call is_density_file_Voronoi() + + write(*,*) lVoronoi + + if (lVoronoi) then + if (lmhd_voronoi) then + call setup_mhd_to_mcfost() !uses sph_to_voronoi + else + call setup_SPH2mcfost(density_file, limits_file, n_SPH, extra_heating) + endif call setup_grid() else + ! Setting up a regular grid call setup_grid() call define_grid() ! included in setup_phantom2mcfost call stars_cell_indices() - endif - - laffichage=.true. - if (.not.(lphantom_file .or. lgadget2_file .or. lascii_SPH_file .or. lmhd_voronoi)) then ! already done by setup_SPH2mcfost call allocate_densities() if (ldensity_file) then call read_density_file() From 3b149980a67cc135cb19fe924c9a2b3e9e9c74df Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 22 May 2024 10:42:46 +1000 Subject: [PATCH 03/13] Adding error messages --- src/SPH2mcfost.f90 | 2 ++ src/density.f90 | 75 ++++++++++++++++++++++++++----------------- src/dust_transfer.f90 | 2 -- 3 files changed, 48 insertions(+), 31 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 23b23829c..233adb630 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -90,6 +90,8 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) write(*,*) "Performing SPH2mcfost setup" write(*,*) "Reading SPH density file: "//trim(SPH_file) call read_ascii_SPH_file(iunit,SPH_file, x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,particle_id, ndusttypes,n_SPH,ierr) + else + call error("Unknown SPH structure.") endif write(*,*) "Done" diff --git a/src/density.f90 b/src/density.f90 index 806abeedc..0a8cdfc49 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -949,6 +949,7 @@ subroutine is_density_file_Voronoi() ! Opening file readwrite=0 + write(*,*) "Checking "//trim(density_file) call ftopen(unit,density_file,readwrite,blocksize,status) if (status /= 0) call error("density file needed") @@ -961,8 +962,12 @@ subroutine is_density_file_Voronoi() if (status /= 0) call error("error reading density file") if (nfound == 1) then + write(*,*) "Found 1D density structure, using a Voronoi mesh" lVoronoi = .true. + l3D = .true. + n_cells = naxes(1) else + write(*,*) "Found a structured mesh" lVoronoi = .false. endif @@ -1155,7 +1160,6 @@ subroutine read_density_file() endif if (nfound == 1) then ! Voronoi mesh - write(*,*) "Found 1D density structure, using a Voronoi mesh" lVoronoi = .true. l3D = .true. n_cells = naxes(1) @@ -1558,38 +1562,51 @@ subroutine read_density_file() call ftclos(unit, status) call ftfiou(unit, status) + ! Passing the densities and velocities to mcfost arrays - do k=1, n_az - do j=j_start,nz - if (l3D_file) then - jj = j - if (jj > 0) then - jj = nz + jj - else - jj = nz+1 + jj - endif - else - jj = abs(j) - endif - if (j==0) then - !densite_gaz(cell_map(i,j,k)) =0.0 - else - do i=1, n_rad - icell = cell_map(i,j,k) - if (lread_gas_density) then - densite_gaz(icell) = sph_gas_dens(i,jj,k) + if (lVoronoi) then + if (lread_gas_density) then + do icell=1,n_cells + densite_gaz(icell) = sph_gas_dens(icell,1,1) + enddo + else + do icell=1,n_cells + densite_gaz(icell) = sph_dens(icell,1,1,1) + enddo + endif + else + do k=1, n_az + do j=j_start,nz + if (l3D_file) then + jj = j + if (jj > 0) then + jj = nz + jj else - densite_gaz(icell) = sph_dens(i,jj,k,1) ! gaz = plus petites particules + jj = nz+1 + jj endif + else + jj = abs(j) + endif + if (j==0) then + !densite_gaz(cell_map(i,j,k)) =0.0 + else + do i=1, n_rad + icell = cell_map(i,j,k) + if (lread_gas_density) then + densite_gaz(icell) = sph_gas_dens(i,jj,k) + else + densite_gaz(icell) = sph_dens(i,jj,k,1) ! gaz = plus petites particules + endif - if (lread_gas_velocity) then - vfield3d(icell,:) = sph_V(i,jj,k,:) - if ((.not.l3D_file).and.(j<0)) vfield3d(icell,3) = - vfield3d(icell,3) - endif - enddo ! i - endif ! j==0 - enddo !j - enddo ! k + if (lread_gas_velocity) then + vfield3d(icell,:) = sph_V(i,jj,k,:) + if ((.not.l3D_file).and.(j<0)) vfield3d(icell,3) = - vfield3d(icell,3) + endif + enddo ! i + endif ! j==0 + enddo !j + enddo ! k + endif ! Calcul de la masse de gaz de la zone mass = 0. diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index a2bb7edfe..2ebe03a49 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -111,8 +111,6 @@ subroutine transfert_poussiere() if (ldensity_file) call is_density_file_Voronoi() - write(*,*) lVoronoi - if (lVoronoi) then if (lmhd_voronoi) then call setup_mhd_to_mcfost() !uses sph_to_voronoi From 48f17077a35a229de7dad9a478b3a35d25390c06 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 22 May 2024 14:08:39 +1000 Subject: [PATCH 04/13] Performing tesselation --- src/density.f90 | 43 ++++++++++++++++++++++++++++++++++++++++--- src/dust_transfer.f90 | 6 ++++-- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/src/density.f90 b/src/density.f90 index 0a8cdfc49..a5de763a6 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -1020,6 +1020,12 @@ subroutine read_density_file() type(star_type), dimension(:), allocatable :: etoile_old integer :: i_etoile, n_etoiles_old + real(dp), dimension(6) :: limits + real(dp), dimension(:), allocatable :: x,y,z, vx,vy,vz, h + integer, dimension(:), allocatable :: is_ghost, mask, particle_id + logical :: check_previous_tesselation + integer :: n_points + if (lplanet_az) then ! Planet is along x-axis by default in wakeflow RT_n_az = 1 RT_az_min = planet_az @@ -1563,6 +1569,8 @@ subroutine read_density_file() call ftfiou(unit, status) + write(*,*) "lVoronoi", lVoronoi + ! Passing the densities and velocities to mcfost arrays if (lVoronoi) then if (lread_gas_density) then @@ -1608,6 +1616,37 @@ subroutine read_density_file() enddo ! k endif + ! Performing tesselation + if (lVoronoi) then + n_points = n_cells + allocate(x(n_points),y(n_points),z(n_points),vx(n_points),vy(n_points),vz(n_points), & + h(n_points), is_ghost(n_points), mask(n_points), particle_id(n_points)) + vx = 0.0 ; vy = 0.0 ; vz = 0.0 ; h=0.0 + + do i=1,n_points + particle_id(i) = i + x(i) = xyz_sph_dp(1,i) + y(i) = xyz_sph_dp(2,i) + z(i) = xyz_sph_dp(3,i) + enddo + + mask = 0 + is_ghost = 0 + + ! We do not trum particles for now + limits(1) = minval(x) ; limits(2) = maxval(x) + limits(3) = minval(y) ; limits(2) = maxval(y) + limits(5) = minval(z) ; limits(2) = maxval(z) + + check_previous_tesselation = .false. + + call Voronoi_tesselation(n_points, particle_id, x,y,z,h,vx,vy,vz, is_ghost, mask, limits, check_previous_tesselation) + !deallocate(x,y,z) + write(*,*) "Using n_cells =", n_cells + endif + + + ! Calcul de la masse de gaz de la zone mass = 0. do icell=1,n_cells @@ -1631,6 +1670,7 @@ subroutine read_density_file() do icell=1,n_cells masse_gaz(icell) = densite_gaz(icell) * masse_mol_gaz * volume(icell) * AU3_to_m3 enddo ! icell + write(*,*) 'Total gas mass in model:', real(sum(masse_gaz) * g_to_Msun),' Msun' if (lvariable_dust) then write(*,*) "Differential spatial distribution" @@ -1720,12 +1760,9 @@ subroutine read_density_file() enddo ! k endif !lvariable_dust - write(*,*) 'Total gas mass in model:', real(sum(masse_gaz) * g_to_Msun),' Msun' call normalize_dust_density() deallocate(sph_dens,a_sph) - - return end subroutine read_density_file diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index 2ebe03a49..41e98a21a 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -114,12 +114,14 @@ subroutine transfert_poussiere() if (lVoronoi) then if (lmhd_voronoi) then call setup_mhd_to_mcfost() !uses sph_to_voronoi + else if (ldensity_file) then + call allocate_densities(n_cells + n_etoiles) + call read_density_file() else call setup_SPH2mcfost(density_file, limits_file, n_SPH, extra_heating) endif call setup_grid() - else - ! Setting up a regular grid + else ! Setting up a regular grid call setup_grid() call define_grid() ! included in setup_phantom2mcfost call stars_cell_indices() From e833bee1a0f07f020d9fc7832f4372b2d5d22e23 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 22 May 2024 19:40:58 +1000 Subject: [PATCH 05/13] Small update --- src/density.f90 | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/src/density.f90 b/src/density.f90 index a5de763a6..2a1aa62ee 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -940,7 +940,8 @@ end subroutine define_density_wall3D subroutine is_density_file_Voronoi() - integer :: unit, status, nfound, readwrite, blocksize + integer :: unit, status, naxis, readwrite, blocksize, nfound + character(len=8) :: comment integer, dimension(4) :: naxes status=0 @@ -953,24 +954,28 @@ subroutine is_density_file_Voronoi() call ftopen(unit,density_file,readwrite,blocksize,status) if (status /= 0) call error("density file needed") - nfound=0 - call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) - ! Closing file - call ftclos(unit, status) - call ftfiou(unit, status) - + !nfound=0 + call ftgkyj(unit,'NAXIS',naxis,comment,status) if (status /= 0) call error("error reading density file") - if (nfound == 1) then + nfound=0 + call ftgknj(unit,'NAXIS',1,4,naxes,nfound,status) + + if (naxis == 1) then write(*,*) "Found 1D density structure, using a Voronoi mesh" lVoronoi = .true. l3D = .true. + n_cells = naxes(1) else write(*,*) "Found a structured mesh" lVoronoi = .false. endif + ! Closing file + call ftclos(unit, status) + call ftfiou(unit, status) + return end subroutine is_density_file_Voronoi @@ -1621,7 +1626,7 @@ subroutine read_density_file() n_points = n_cells allocate(x(n_points),y(n_points),z(n_points),vx(n_points),vy(n_points),vz(n_points), & h(n_points), is_ghost(n_points), mask(n_points), particle_id(n_points)) - vx = 0.0 ; vy = 0.0 ; vz = 0.0 ; h=0.0 + vx = 0.0 ; vy = 0.0 ; vz = 0.0 ; h=1e30 do i=1,n_points particle_id(i) = i @@ -1631,12 +1636,18 @@ subroutine read_density_file() enddo mask = 0 - is_ghost = 0 + !call test_duplicate_particles(n_points, particle_id, x,y,z, massgas,massdust,rho,rhodust, is_ghost) + ! We do not trum particles for now limits(1) = minval(x) ; limits(2) = maxval(x) - limits(3) = minval(y) ; limits(2) = maxval(y) - limits(5) = minval(z) ; limits(2) = maxval(z) + limits(3) = minval(y) ; limits(4) = maxval(y) + limits(5) = minval(z) ; limits(6) = maxval(z) + + write(*,*) "# Model limits :" + write(*,*) "x =", limits(1), limits(2) + write(*,*) "y =", limits(3), limits(4) + write(*,*) "z =", limits(5), limits(6) check_previous_tesselation = .false. From 21ca586253cc13d8165427c20635243053784557 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Thu, 23 May 2024 07:31:12 +1000 Subject: [PATCH 06/13] Increasing number of neighbours --- src/Voronoi.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Voronoi.f90 b/src/Voronoi.f90 index 3a7c5197c..10b283158 100644 --- a/src/Voronoi.f90 +++ b/src/Voronoi.f90 @@ -210,7 +210,7 @@ subroutine Voronoi_tesselation(n_points, particle_id, x,y,z,h, vx,vy,vz, is_ghos logical, intent(in) :: check_previous_tesselation - integer, parameter :: max_neighbours = 30 ! maximum number of neighbours per cell (to build temporary neighbours list) + integer, parameter :: max_neighbours = 100 ! maximum number of neighbours per cell (to build temporary neighbours list) real(kind=dp), dimension(:), allocatable :: x_tmp, y_tmp, z_tmp, h_tmp integer, dimension(:), allocatable :: SPH_id, SPH_original_id From 98d75a038bec7550dd44f278eb351078193efbf5 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 14:36:04 +1000 Subject: [PATCH 07/13] Merging Voronoi fits and SPH interfaces --- src/SPH2mcfost.f90 | 13 +- src/density.f90 | 405 +++++++++++++++++++++--------------------- src/dust_transfer.f90 | 7 +- src/init_mcfost.f90 | 25 ++- src/mhd2mcfost.f90 | 6 +- src/parameters.f90 | 2 +- 6 files changed, 233 insertions(+), 225 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 233adb630..0265d8d40 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -4,7 +4,7 @@ module SPH2mcfost use constantes use utils use sort, only : find_kth_smallest_inplace - use density, only : normalize_dust_density, reduce_density + use density, only : normalize_dust_density, reduce_density, read_Voronoi_fits_file use read_phantom, only : read_phantom_bin_files, read_phantom_hdf_files use sort, only : index_quicksort use stars, only : compute_stellar_parameters @@ -89,7 +89,10 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) else if (lascii_SPH_file) then write(*,*) "Performing SPH2mcfost setup" write(*,*) "Reading SPH density file: "//trim(SPH_file) - call read_ascii_SPH_file(iunit,SPH_file, x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,particle_id, ndusttypes,n_SPH,ierr) + call read_ascii_SPH_file(iunit,SPH_file, x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,& + particle_id, ndusttypes,n_SPH,ierr) + else if (ldensity_file) then + call read_Voronoi_fits_file(SPH_file, x,y,z,h,vx,vy,vz,particle_id,massgas,n_SPH) else call error("Unknown SPH structure.") endif @@ -101,7 +104,7 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) if (allocated(massdust)) deallocate(massdust) endif - if ((.not.lfix_star).and.(lphantom_file .or. lgadget2_file)) call compute_stellar_parameters() + if (.not.lfix_star) call compute_stellar_parameters() ! Model limits call read_SPH_limits_file(SPH_limits_file, SPH_limits) @@ -109,8 +112,8 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) ! Voronoi tesselation check_previous_tesselation = (.not. lrandomize_Voronoi) call SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, & - T_gas, massgas,massdust,rho,rhodust,SPH_grainsizes, SPH_limits, check_previous_tesselation, is_ghost, & - ldust_moments, dust_moments, mass_per_H, mask=mask) + T_gas, massgas,massdust,rho,rhodust,SPH_grainsizes, SPH_limits, check_previous_tesselation, & + is_ghost, ldust_moments, dust_moments, mass_per_H, mask=mask) deallocate(x,y,z,h) if (allocated(vx)) deallocate(vx,vy,vz) diff --git a/src/density.f90 b/src/density.f90 index 2a1aa62ee..d5b79bf57 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -19,7 +19,7 @@ module density public :: densite_gaz, masse_gaz, surface_density, densite_gaz_midplane, densite_pouss, masse, icell_not_empty, & define_density, define_density_wall3d, define_dust_density, read_density_file, is_density_file_Voronoi, & densite_seb_charnoz2, densite_seb_charnoz, remove_species, read_sigma_file, normalize_dust_density, & - reduce_density + reduce_density, read_Voronoi_fits_file private @@ -950,8 +950,8 @@ subroutine is_density_file_Voronoi() ! Opening file readwrite=0 - write(*,*) "Checking "//trim(density_file) - call ftopen(unit,density_file,readwrite,blocksize,status) + write(*,*) "Checking "//trim(density_files(1)) + call ftopen(unit,density_files(1),readwrite,blocksize,status) if (status /= 0) call error("density file needed") !nfound=0 @@ -1012,12 +1012,10 @@ subroutine read_density_file() real, dimension(:,:,:), allocatable :: sph_gas_dens ! (n_rad,nz,n_az) real, dimension(:,:,:,:), allocatable :: sph_V ! (n_rad,nz,n_az,3) real, dimension(:), allocatable :: a_sph, n_a_sph, log_a_sph, log_n_a_sph ! n_a - real, dimension(:,:), allocatable :: xyz_sph real(kind=dp), dimension(:,:,:,:), allocatable :: sph_dens_dp, sph_V_dp real(kind=dp), dimension(:,:,:), allocatable :: sph_gas_dens_dp real(kind=dp), dimension(:), allocatable :: a_sph_dp - real(kind=dp), dimension(:,:), allocatable :: xyz_sph_dp real(kind=dp) :: f logical :: lread_gas_density, lread_gas_velocity @@ -1025,12 +1023,6 @@ subroutine read_density_file() type(star_type), dimension(:), allocatable :: etoile_old integer :: i_etoile, n_etoiles_old - real(dp), dimension(6) :: limits - real(dp), dimension(:), allocatable :: x,y,z, vx,vy,vz, h - integer, dimension(:), allocatable :: is_ghost, mask, particle_id - logical :: check_previous_tesselation - integer :: n_points - if (lplanet_az) then ! Planet is along x-axis by default in wakeflow RT_n_az = 1 RT_az_min = planet_az @@ -1044,10 +1036,10 @@ subroutine read_density_file() ! Get an unused Logical Unit Number to use to open the FITS file. call ftgiou(unit,status) - write(*,*) "Reading density file : "//trim(density_file) + write(*,*) "Reading density file : "//trim(density_files(1)) readwrite=0 - call ftopen(unit,density_file,readwrite,blocksize,status) + call ftopen(unit,density_files(1),readwrite,blocksize,status) if (status /= 0) call error("density file needed") ! Do we read the gas density ? @@ -1165,65 +1157,55 @@ subroutine read_density_file() write(*,*) "Reading dust density ..." call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) - if ((nfound /= 1) .and. (nfound /= 3) .and. (nfound /= 4)) then - write(*,*) "I found", nfound, "axis instead of 1, 3 or 4" - call error('failed to read the NAXIS keyword in HDU 1 of '//trim(density_file)//' file') + if ((nfound /= 3) .and. (nfound /= 4)) then + write(*,*) "I found", nfound, "axis instead of 3 or 4" + call error('failed to read the NAXIS keyword in HDU 1 of '//trim(density_files(1))//' file') endif - if (nfound == 1) then ! Voronoi mesh - lVoronoi = .true. - l3D = .true. - n_cells = naxes(1) - npixels = naxes(1) - n_a = 1 ! for now - - allocate(sph_dens(n_cells,1,1,n_a), a_sph(n_a), n_a_sph(n_a)) - else ! Structured grid - if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then - write(*,*) "# fits_file vs mcfost_grid" - write(*,*) naxes(1), n_rad - write(*,*) naxes(2), nz - write(*,*) naxes(3), n_az - !write(*,*) naxes(4), n_a - call error(trim(density_file)//" does not have the right dimensions in HDU 1.") - endif + if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then + write(*,*) "# fits_file vs mcfost_grid" + write(*,*) naxes(1), n_rad + write(*,*) naxes(2), nz + write(*,*) naxes(3), n_az + !write(*,*) naxes(4), n_a + call error(trim(density_files(1))//" does not have the right dimensions in HDU 1.") + endif - if (nfound == 3) then - n_a = 1 - write(*,*) "No grain size found" - npixels=naxes(1)*naxes(2)*naxes(3) - else - if (.not.lvariable_dust) then - call warning("Forcing variable dust") - lvariable_dust = .true. - endif - n_a = naxes(4) - write(*,*) n_a, "grain sizes found" - npixels=naxes(1)*naxes(2)*naxes(3)*naxes(4) + if (nfound == 3) then + n_a = 1 + write(*,*) "No grain size found" + npixels=naxes(1)*naxes(2)*naxes(3) + else + if (.not.lvariable_dust) then + call warning("Forcing variable dust") + lvariable_dust = .true. endif + n_a = naxes(4) + write(*,*) n_a, "grain sizes found" + npixels=naxes(1)*naxes(2)*naxes(3)*naxes(4) + endif - if (lvariable_dust) then - p_n_cells = n_cells - p_n_rad=n_rad ; p_nz = nz - else - p_n_cells = 1 - p_n_rad=1 ; p_nz=1 - endif + if (lvariable_dust) then + p_n_cells = n_cells + p_n_rad=n_rad ; p_nz = nz + else + p_n_cells = 1 + p_n_rad=1 ; p_nz=1 + endif - if (naxes(2) == 2*nz) then - l3D_file = .true. - else - write(*,*) "The density file only has positive z, making it symmetric" - l3D_file = .false. - endif + if (naxes(2) == 2*nz) then + l3D_file = .true. + else + write(*,*) "The density file only has positive z, making it symmetric" + l3D_file = .false. + endif - if (l3D_file) then - allocate(sph_dens(n_rad,2*nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) - else - allocate(sph_dens(n_rad,nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) - endif - sph_dens = 0.0 ; a_sph = 0.0 ; n_a_sph = 0.0 + if (l3D_file) then + allocate(sph_dens(n_rad,2*nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) + else + allocate(sph_dens(n_rad,nz,n_az,n_a), a_sph(n_a), n_a_sph(n_a)) endif + sph_dens = 0.0 ; a_sph = 0.0 ; n_a_sph = 0.0 bitpix = 0 call ftgkyj(unit,"bitpix",bitpix,comment,status) @@ -1232,15 +1214,12 @@ subroutine read_density_file() if (bitpix==-32) then call ftgpve(unit,group,firstpix,npixels,nullval,sph_dens,anynull,status) else if (bitpix==-64) then - if (lVoronoi) then - allocate(sph_dens_dp(n_cells,1,1,n_a)) + if (l3D_file) then + allocate(sph_dens_dp(n_rad,2*nz,n_az,n_a)) else - if (l3D_file) then - allocate(sph_dens_dp(n_rad,2*nz,n_az,n_a)) - else - allocate(sph_dens_dp(n_rad,nz,n_az,n_a)) - endif + allocate(sph_dens_dp(n_rad,nz,n_az,n_a)) endif + sph_dens_dp = 0.0_dp call ftgpvd(unit,group,firstpix,npixels,nullval,sph_dens_dp,anynull,status) sph_dens = real(sph_dens_dp,kind=sp) @@ -1255,49 +1234,6 @@ subroutine read_density_file() sph_dens = sph_dens/maxval(sph_dens) ! normalization avant d'ajouter une constante sph_dens = max(sph_dens,1e10*tiny_real) - - if (lVoronoi) then - !--------------------------------------------------------- - ! hdu2 is particle positions - !--------------------------------------------------------- - ! move to next hdu - call ftmrhd(unit,1,hdutype,status) - nfound=1 - ! Check dimensions - naxes(:) = 0 - call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) - if (nfound /= 2) then - write(*,*) 'HDU 2 has', nfound, 'dimensions.' - call error('did not find 2 dimension in HDU 2') - endif - if ((naxes(1) /= 3)) then - write(*,*) "HDU2 dim 1 is ", naxes(1), "instead of 3" - call error("HDU 2 does not have the right dimension") - endif - if ((naxes(2) /= n_cells)) then - write(*,*) "HDU2 dim 2 is ", naxes(3), "instead of ",n_cells - call error("HDU 2 does not have the right dimension") - endif - - npixels=naxes(1)*naxes(2) - - bitpix=0 - call ftgkyj(unit,"bitpix",bitpix,comment,status) - - ! read_image - if (bitpix==-32) then - allocate(xyz_sph(3,n_cells),xyz_sph_dp(3,n_cells)) - call ftgpve(unit,group,firstpix,npixels,nullval,xyz_sph,anynull,status) - xyz_sph_dp = real(xyz_sph,kind=dp) - deallocate(xyz_sph) - else if (bitpix==-64) then - allocate(xyz_sph_dp(3,n_cells)) - call ftgpvd(unit,group,firstpix,npixels,nullval,xyz_sph_dp,anynull,status) - else - call error("cannot read bitpix in fits file") - endif - endif - if (n_a > 1) then write(*,*) "Reading grain sizes ..." read_n_a = 0 @@ -1458,7 +1394,7 @@ subroutine read_density_file() if (nfound /= 3) then write(*,*) 'Gas density:' write(*,*) "I found", nfound, "axis instead of 3, status=", status - call error('failed to read the NAXISn keywords of '//trim(density_file)//' file.') + call error('failed to read the NAXISn keywords of '//trim(density_files(1))) endif if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az) ) then @@ -1466,7 +1402,7 @@ subroutine read_density_file() write(*,*) naxes(1), n_rad write(*,*) naxes(2), nz write(*,*) naxes(3), n_az - call error(trim(density_file)//" does not have the right dimensions.") + call error(trim(density_files(1))//" does not have the right dimensions.") endif npixels=naxes(1)*naxes(2)*naxes(3) @@ -1531,7 +1467,7 @@ subroutine read_density_file() if (nfound /= 4) then write(*,*) 'Gas velocity:' write(*,*) "I found", nfound, "axis instead of 4" - call error('failed to read the NAXISn keywords of '//trim(density_file)//' file.') + call error('failed to read the NAXISn keywords of '//trim(density_files(1))//' file.') endif if ((naxes(1) /= n_rad).or.((naxes(2) /= nz).and.(naxes(2) /= 2*nz)).or.(naxes(3) /= n_az).or.(naxes(4) /= 3) ) then @@ -1540,7 +1476,7 @@ subroutine read_density_file() write(*,*) naxes(2), nz write(*,*) naxes(3), n_az write(*,*) naxes(4), 3 - call error(trim(density_file)//" does not have the right dimensions.") + call error(trim(density_files(1))//" does not have the right dimensions.") endif npixels=naxes(1)*naxes(2)*naxes(3)*naxes(4) @@ -1574,89 +1510,38 @@ subroutine read_density_file() call ftfiou(unit, status) - write(*,*) "lVoronoi", lVoronoi - ! Passing the densities and velocities to mcfost arrays - if (lVoronoi) then - if (lread_gas_density) then - do icell=1,n_cells - densite_gaz(icell) = sph_gas_dens(icell,1,1) - enddo - else - do icell=1,n_cells - densite_gaz(icell) = sph_dens(icell,1,1,1) - enddo - endif - else - do k=1, n_az - do j=j_start,nz - if (l3D_file) then - jj = j - if (jj > 0) then - jj = nz + jj - else - jj = nz+1 + jj - endif + do k=1, n_az + do j=j_start,nz + if (l3D_file) then + jj = j + if (jj > 0) then + jj = nz + jj else - jj = abs(j) + jj = nz+1 + jj endif - if (j==0) then - !densite_gaz(cell_map(i,j,k)) =0.0 - else - do i=1, n_rad - icell = cell_map(i,j,k) - if (lread_gas_density) then - densite_gaz(icell) = sph_gas_dens(i,jj,k) - else - densite_gaz(icell) = sph_dens(i,jj,k,1) ! gaz = plus petites particules - endif - - if (lread_gas_velocity) then - vfield3d(icell,:) = sph_V(i,jj,k,:) - if ((.not.l3D_file).and.(j<0)) vfield3d(icell,3) = - vfield3d(icell,3) - endif - enddo ! i - endif ! j==0 - enddo !j - enddo ! k - endif - - ! Performing tesselation - if (lVoronoi) then - n_points = n_cells - allocate(x(n_points),y(n_points),z(n_points),vx(n_points),vy(n_points),vz(n_points), & - h(n_points), is_ghost(n_points), mask(n_points), particle_id(n_points)) - vx = 0.0 ; vy = 0.0 ; vz = 0.0 ; h=1e30 - - do i=1,n_points - particle_id(i) = i - x(i) = xyz_sph_dp(1,i) - y(i) = xyz_sph_dp(2,i) - z(i) = xyz_sph_dp(3,i) - enddo - - mask = 0 - !call test_duplicate_particles(n_points, particle_id, x,y,z, massgas,massdust,rho,rhodust, is_ghost) - - - ! We do not trum particles for now - limits(1) = minval(x) ; limits(2) = maxval(x) - limits(3) = minval(y) ; limits(4) = maxval(y) - limits(5) = minval(z) ; limits(6) = maxval(z) - - write(*,*) "# Model limits :" - write(*,*) "x =", limits(1), limits(2) - write(*,*) "y =", limits(3), limits(4) - write(*,*) "z =", limits(5), limits(6) - - check_previous_tesselation = .false. - - call Voronoi_tesselation(n_points, particle_id, x,y,z,h,vx,vy,vz, is_ghost, mask, limits, check_previous_tesselation) - !deallocate(x,y,z) - write(*,*) "Using n_cells =", n_cells - endif - + else + jj = abs(j) + endif + if (j==0) then + !densite_gaz(cell_map(i,j,k)) =0.0 + else + do i=1, n_rad + icell = cell_map(i,j,k) + if (lread_gas_density) then + densite_gaz(icell) = sph_gas_dens(i,jj,k) + else + densite_gaz(icell) = sph_dens(i,jj,k,1) ! gaz = plus petites particules + endif + if (lread_gas_velocity) then + vfield3d(icell,:) = sph_V(i,jj,k,:) + if ((.not.l3D_file).and.(j<0)) vfield3d(icell,3) = - vfield3d(icell,3) + endif + enddo ! i + endif ! j==0 + enddo !j + enddo ! k ! Calcul de la masse de gaz de la zone mass = 0. @@ -1780,6 +1665,124 @@ end subroutine read_density_file !********************************************************** +subroutine read_Voronoi_fits_file(filename, x,y,z,h,vx,vy,vz,particle_id, massgas,n_points) + + + character(len=*), intent(in) :: filename + + real(dp), intent(out), dimension(:), allocatable :: x,y,z,h,vx,vy,vz,massgas + integer, intent(out), dimension(:), allocatable :: particle_id + integer, intent(out) :: n_points + + real, dimension(:,:), allocatable :: xyz_sp + real(kind=dp), dimension(:,:), allocatable :: xyz_dp + real, dimension(:), allocatable :: mass_sp + + integer :: status, readwrite, unit, blocksize,nfound,group,firstpix,npixels,hdutype, bitpix, nullval, i + logical :: anynull + character(len=80) :: comment + integer, dimension(4) :: naxes + + + write(*,*) "Reading Voronoi fits file: "//trim(filename) + + ! Lecture donnees + status=0 + ! Get an unused Logical Unit Number to use to open the FITS file. + call ftgiou(unit,status) + + readwrite=0 + call ftopen(unit,filename,readwrite,blocksize,status) + if (status /= 0) call error("density file needed") + + status = 0 + group=1 + firstpix=1 + nullval=-999 + call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) + + if ((nfound /= 1)) then + write(*,*) "I found", nfound, "axis instead of 1" + call error('failed to read the NAXIS keyword in HDU 1 of '//trim(filename)) + endif + + n_points = naxes(1) + npixels = naxes(1) + + allocate(massgas(n_points), x(n_points),y(n_points),z(n_points), h(n_points), particle_id(n_points)) + h = 0.0_dp + + bitpix = 0 + call ftgkyj(unit,"bitpix",bitpix,comment,status) + + if (bitpix==-32) then + allocate(mass_sp(n_points)) + call ftgpve(unit,group,firstpix,npixels,nullval,mass_sp,anynull,status) + massgas = mass_sp + deallocate(mass_sp) + else if (bitpix==-64) then + call ftgpvd(unit,group,firstpix,npixels,nullval,massgas,anynull,status) + else + call error("cannot read bitpix in fits file") + endif + + + !--------------------------------------------------------- + ! hdu2 is particle positions + !--------------------------------------------------------- + ! move to next hdu + call ftmrhd(unit,1,hdutype,status) + nfound=1 + ! Check dimensions + naxes(:) = 0 + call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) + if (nfound /= 2) then + write(*,*) 'HDU 2 has', nfound, 'dimensions.' + call error('did not find 2 dimension in HDU 2') + endif + if ((naxes(1) /= 3)) then + write(*,*) "HDU2 dim 1 is ", naxes(1), "instead of 3" + call error("HDU 2 does not have the right dimension") + endif + if ((naxes(2) /= n_points)) then + write(*,*) "HDU2 dim 2 is ", naxes(3), "instead of ", n_points + call error("HDU 2 does not have the right dimension") + endif + + npixels=naxes(1)*naxes(2) + + bitpix=0 + call ftgkyj(unit,"bitpix",bitpix,comment,status) + + ! read_image + if (bitpix==-32) then + allocate(xyz_sp(3,n_cells)) + call ftgpve(unit,group,firstpix,npixels,nullval,xyz_sp,anynull,status) + x = xyz_sp(1,:) + y = xyz_sp(2,:) + z = xyz_sp(3,:) + deallocate(xyz_sp) + else if (bitpix==-64) then + allocate(xyz_dp(3,n_cells)) + call ftgpvd(unit,group,firstpix,npixels,nullval,xyz_dp,anynull,status) + x = xyz_dp(1,:) + y = xyz_dp(2,:) + z = xyz_dp(3,:) + deallocate(xyz_dp) + else + call error("cannot read bitpix in fits file") + endif + + do i=1,n_points + particle_id(i) = i + enddo + + return + +end subroutine read_Voronoi_fits_file + +!********************************************************** + subroutine normalize_dust_density(disk_dust_mass) ! Normalize the dust density to follow the slope in nbre_grains ! and the mass in each dust population @@ -2097,10 +2100,10 @@ subroutine densite_Seb_Charnoz2() call ftgiou(unit,status) write(*,*) "Density structure from Seb. Charnoz" ; - write(*,*) "Reading density file : "//trim(density_file) + write(*,*) "Reading density file : "//trim(density_files(1)) readwrite=0 - call ftopen(unit,density_file,readwrite,blocksize,status) + call ftopen(unit,density_files(1),readwrite,blocksize,status) if (status /= 0) call error("density file needed") group=1 @@ -2109,13 +2112,13 @@ subroutine densite_Seb_Charnoz2() ! determine the size of density file call ftgknj(unit,'NAXIS',1,10,naxes,nfound,status) - if (nfound /= 2) call error('failed to read the NAXISn keywords of '//trim(density_file)//' file.') + if (nfound /= 2) call error('failed to read the NAXISn keywords of '//trim(density_files(1))//' file.') if ((naxes(1) /= n_rad).or.(naxes(2) /= nz) ) then write(*,*) "# fits_file vs mcfost_grid" write(*,*) naxes(1), n_rad write(*,*) naxes(2), nz - call error(trim(density_file)//" does not have the right dimensions.") + call error(trim(density_files(1))//" does not have the right dimensions.") endif npixels=naxes(1)*naxes(2) diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index 41e98a21a..608fcfaf3 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -114,11 +114,8 @@ subroutine transfert_poussiere() if (lVoronoi) then if (lmhd_voronoi) then call setup_mhd_to_mcfost() !uses sph_to_voronoi - else if (ldensity_file) then - call allocate_densities(n_cells + n_etoiles) - call read_density_file() else - call setup_SPH2mcfost(density_file, limits_file, n_SPH, extra_heating) + call setup_SPH2mcfost(density_files(1), limits_file, n_SPH, extra_heating) endif call setup_grid() else ! Setting up a regular grid @@ -139,7 +136,7 @@ subroutine transfert_poussiere() call read_athena_model() else if (lsphere_model) then !on a structured spherical grid - call read_spherical_model(density_file) + call read_spherical_model(density_files(1)) else if (lmodel_1d) then !1d spherically symmetric "stellar atmosphere" models call setup_model1d_to_mcfost() else if (lidefix) then diff --git a/src/init_mcfost.f90 b/src/init_mcfost.f90 index edca940d4..316df9c7c 100644 --- a/src/init_mcfost.f90 +++ b/src/init_mcfost.f90 @@ -820,7 +820,8 @@ subroutine initialisation_mcfost() i_arg = i_arg + 1 ldensity_file=.true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 case("-start_step") i_arg = i_arg + 1 @@ -841,19 +842,22 @@ subroutine initialisation_mcfost() lVoronoi = .true. l3D = .true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 case ("-model_1d") i_arg = i_arg + 1 lmodel_1d = .true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 case("-sphere_mesh") i_arg = i_arg + 1 lsphere_model = .true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 case("-zeeman_polarisation") call error("Zeeman polarisation not yet!") @@ -961,7 +965,8 @@ subroutine initialisation_mcfost() lVoronoi = .true. l3D = .true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 if (.not.llimits_file) limits_file = "phantom.limits" case("-SPH_amin") @@ -985,7 +990,8 @@ subroutine initialisation_mcfost() lVoronoi = .true. l3D = .true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 if (.not.llimits_file) limits_file = "gadget2.limits" case("-limits_file","-limits") @@ -1118,7 +1124,8 @@ subroutine initialisation_mcfost() i_arg = i_arg + 1 lread_Seb_Charnoz2=.true. call get_command_argument(i_arg,s) - density_file = s + allocate(density_files(1)) + density_files(1) = s i_arg = i_arg + 1 case("-only_top") i_arg = i_arg+1 @@ -1503,13 +1510,13 @@ subroutine initialisation_mcfost() write(*,*) "------------------------------------------------" call warning(" THERE ARE PROBABLY SOME CHECKS TO DO MORE ") write(*,*) "------------------------------------------------" - call read_model_1d(density_file) + call read_model_1d(density_files(1)) endif if (lsphere_model) then !could be 3d or 2d (2.5d). Depends on flag l3D or N_az>1 n_zones = 1 disk_zone(1)%geometry = 2 - call read_spherical_grid_parameters(density_file) + call read_spherical_grid_parameters(density_files(1)) endif if (lidefix) then diff --git a/src/mhd2mcfost.f90 b/src/mhd2mcfost.f90 index 5bd63969e..2f8c3cf09 100644 --- a/src/mhd2mcfost.f90 +++ b/src/mhd2mcfost.f90 @@ -61,8 +61,6 @@ subroutine setup_mhd_to_mcfost() n_points = 0 ! to avoid compiler warning !needed for Voronoi - if (allocated(density_files)) deallocate(density_files) - allocate(density_files(1)); density_files(1) = density_file if (lignore_dust) then ndusttypes = 0 @@ -71,7 +69,7 @@ subroutine setup_mhd_to_mcfost() call error("Dust not handled yet for pluto models!") endif - cmd = "wc -l "//trim(density_file)//" > ntest.txt" + cmd = "wc -l "//trim(density_files(1))//" > ntest.txt" call appel_syst(cmd,syst_status) open(unit=1,file="ntest.txt",status="old") read(1,*) N_points @@ -81,7 +79,7 @@ subroutine setup_mhd_to_mcfost() N_points = N_points + n_etoiles - open(unit=1,file=density_file, status="old") + open(unit=1,file=density_files(1), status="old") call read_line(1, FormatLine, inputline, Nread) lvelocity_file = .false. diff --git a/src/parameters.f90 b/src/parameters.f90 index bc9e397a0..59d128ee6 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -222,7 +222,7 @@ module parametres ! Vertical scaling of the envelope real :: z_scaling_env - character(len=512) :: density_file, sigma_file, grain_size_file, limits_file + character(len=512) :: sigma_file, grain_size_file, limits_file character(len=512), dimension(:), allocatable :: density_files integer :: n_phantom_files From 4d2454f7720d8b3d7db482702d0fa27d3f9b6418 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 17:28:02 +1000 Subject: [PATCH 08/13] Cleaning code --- src/SPH2mcfost.f90 | 25 ++++++++++++------------- src/Voronoi.f90 | 2 +- src/density.f90 | 3 ++- src/dust_transfer.f90 | 2 +- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 0265d8d40..88c02a19c 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -15,14 +15,12 @@ module SPH2mcfost contains - subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) + subroutine setup_SPH2mcfost(extra_heating) use read_gadget2, only : read_gadget2_file use io_phantom_utils, only : get_error_text use utils, only : read_comments - character(len=512), intent(in) :: SPH_file, SPH_limits_file - integer, intent(out) :: n_SPH integer, parameter :: iunit = 1 @@ -36,7 +34,7 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) integer, dimension(:), allocatable :: is_ghost real(dp), dimension(6) :: SPH_limits real :: factor - integer :: ndusttypes, ierr, i, ilen + integer :: ndusttypes, ierr, i, ilen, n_SPH logical :: check_previous_tesselation, ldust_moments real(dp) :: mass_per_H @@ -84,15 +82,15 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) endif else if (lgadget2_file) then write(*,*) "Performing gadget2mcfost setup" - write(*,*) "Reading Gadget-2 density file: "//trim(SPH_file) - call read_gadget2_file(iunit,SPH_file, x,y,z,h,massgas,rho,rhodust,ndusttypes,n_SPH,ierr) + write(*,*) "Reading Gadget-2 density file: "//trim(density_files(1)) + call read_gadget2_file(iunit,density_files(1), x,y,z,h,massgas,rho,rhodust,ndusttypes,n_SPH,ierr) else if (lascii_SPH_file) then write(*,*) "Performing SPH2mcfost setup" - write(*,*) "Reading SPH density file: "//trim(SPH_file) - call read_ascii_SPH_file(iunit,SPH_file, x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,& + write(*,*) "Reading SPH density file: "//trim(density_files(1)) + call read_ascii_SPH_file(iunit,density_files(1), x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,& particle_id, ndusttypes,n_SPH,ierr) else if (ldensity_file) then - call read_Voronoi_fits_file(SPH_file, x,y,z,h,vx,vy,vz,particle_id,massgas,n_SPH) + call read_Voronoi_fits_file(density_files(1), x,y,z,h,vx,vy,vz,particle_id,massgas,n_SPH) else call error("Unknown SPH structure.") endif @@ -107,7 +105,7 @@ subroutine setup_SPH2mcfost(SPH_file,SPH_limits_file, n_SPH, extra_heating) if (.not.lfix_star) call compute_stellar_parameters() ! Model limits - call read_SPH_limits_file(SPH_limits_file, SPH_limits) + call read_SPH_limits_file(limits_file, SPH_limits) ! Voronoi tesselation check_previous_tesselation = (.not. lrandomize_Voronoi) @@ -205,7 +203,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g use, intrinsic :: ieee_arithmetic integer, intent(in) :: n_SPH, ndusttypes - real(dp), dimension(n_SPH), intent(inout) :: x,y,z,h,massgas!,rho, !move rho to allocatable, assuming not always allocated + real(dp), dimension(:), allocatable, intent(inout) :: x,y,z,h,massgas!,rho, !move rho to allocatable, assuming not always allocated real(dp), dimension(:), allocatable, intent(inout) :: rho real(dp), dimension(:), allocatable, intent(inout) :: vx,vy,vz ! dimension n_SPH or 0 real(dp), dimension(:), allocatable, intent(in) :: T_gas @@ -380,6 +378,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g mass = mass * g_to_Msun if (lforce_Mgas) then ! Todo : testing for now, use a routine to avoid repeting code + write(*,*) "Forcing gas mass to value in parameter file rather" ! Normalisation if (mass > 0.0) then ! pour le cas ou gas_to_dust = 0. facteur = disk_zone(1)%diskmass * disk_zone(1)%gas_to_dust / mass @@ -387,10 +386,10 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g ! Somme sur les zones pour densite finale do icell=1,n_cells densite_gaz(icell) = densite_gaz(icell) * facteur - masse_gaz(icell) = densite_gaz(icell) * masse_mol_gaz * volume(icell) * AU3_to_m3 + masse_gaz(icell) = masse_gaz(icell) * facteur enddo ! icell else - call error('Gas mass is 0') + call error('Gas mass is 0 in hydero file') endif endif diff --git a/src/Voronoi.f90 b/src/Voronoi.f90 index 10b283158..3a7c5197c 100644 --- a/src/Voronoi.f90 +++ b/src/Voronoi.f90 @@ -210,7 +210,7 @@ subroutine Voronoi_tesselation(n_points, particle_id, x,y,z,h, vx,vy,vz, is_ghos logical, intent(in) :: check_previous_tesselation - integer, parameter :: max_neighbours = 100 ! maximum number of neighbours per cell (to build temporary neighbours list) + integer, parameter :: max_neighbours = 30 ! maximum number of neighbours per cell (to build temporary neighbours list) real(kind=dp), dimension(:), allocatable :: x_tmp, y_tmp, z_tmp, h_tmp integer, dimension(:), allocatable :: SPH_id, SPH_original_id diff --git a/src/density.f90 b/src/density.f90 index d5b79bf57..5c605c75d 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -1710,7 +1710,7 @@ subroutine read_Voronoi_fits_file(filename, x,y,z,h,vx,vy,vz,particle_id, massga npixels = naxes(1) allocate(massgas(n_points), x(n_points),y(n_points),z(n_points), h(n_points), particle_id(n_points)) - h = 0.0_dp + h = 1e-6 * huge_dp bitpix = 0 call ftgkyj(unit,"bitpix",bitpix,comment,status) @@ -1776,6 +1776,7 @@ subroutine read_Voronoi_fits_file(filename, x,y,z,h,vx,vy,vz,particle_id, massga do i=1,n_points particle_id(i) = i enddo + SPH_keep_particles = 1.0 return diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index 608fcfaf3..7f430a8f8 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -115,7 +115,7 @@ subroutine transfert_poussiere() if (lmhd_voronoi) then call setup_mhd_to_mcfost() !uses sph_to_voronoi else - call setup_SPH2mcfost(density_files(1), limits_file, n_SPH, extra_heating) + call setup_SPH2mcfost(extra_heating) endif call setup_grid() else ! Setting up a regular grid From ba3c9cfbfde9182a9cb12e4e0e1bb729f84311b7 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 17:30:55 +1000 Subject: [PATCH 09/13] Removing unused code --- src/SPH2mcfost.f90 | 128 +-------------------------------------------- 1 file changed, 2 insertions(+), 126 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 88c02a19c..0eaf96c4a 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -218,8 +218,6 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g integer, dimension(:), allocatable, intent(out) :: is_ghost - - logical :: lwrite_ASCII = .false. ! produce an ASCII file for yorick logical :: use_single_grain real, allocatable, dimension(:) :: a_SPH, log_a_SPH, rho_dust @@ -259,54 +257,6 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g write(*,*) "Found", n_SPH, " hydro sites." endif - if (lwrite_ASCII) then - ! Write the file for the grid version of mcfost - !- N_part: total number of particles - ! - r_in: disk inner edge in AU - ! - r_out: disk outer edge in AU - ! - p: surface density exponent, Sigma=Sigma_0*(r/r_0)^(-p), p>0 - ! - q: temperature exponent, T=T_0*(r/r_0)^(-q), q>0 - ! - m_star: star mass in solar masses - ! - m_disk: disk mass in solar masses (99% gas + 1% dust) - ! - H_0: disk scale height at 100 AU, in AU - ! - rho_d: dust density in g.cm^-3 - ! - flag_ggrowth: T with grain growth, F without - ! - ! - ! N_part lines containing: - ! - x,y,z: coordinates of each particle in AU - ! - h: smoothing length of each particle in AU - ! - s: grain size of each particle in �m - ! - ! Without grain growth: 2 lines containing: - ! - n_sizes: number of grain sizes - ! - (s(i),i=1,n_sizes): grain sizes in �m - ! OR - ! With grain growth: 1 line containing: - ! - s_min,s_max: smallest and largest grain size in �m - - open(unit=1,file="SPH_phantom.txt",status="replace") - write(1,*) size(x) - write(1,*) minval(sqrt(x**2 + y**2)) - write(1,*) maxval(sqrt(x**2 + y**2)) - write(1,*) 1 ! p - write(1,*) 0.5 ! q - write(1,*) 1.0 ! mstar - write(1,*) 1.e-3 !mdisk - write(1,*) 10 ! h0 - write(1,*) 3.5 ! rhod - write(1,*) .false. - !rhoi = massoftype(itypei)*(hfact/hi)**3 * udens ! g/cm**3 - - do icell=1,size(x) - write(1,*) x(icell), y(icell), z(icell), 1.0, 1.0 - enddo - - write(1,*) 1 - write(1,*) 1.0 - close(unit=1) - endif - if (abs(maxval(SPH_limits)) < tiny_real) then write(*,*) "Selecting spatial range which contains" write(*,*) SPH_keep_particles*100, "% of particles in each dimension" @@ -403,7 +353,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g iSPH = Voronoi(icell)%id mass = mass + massgas(iSPH) * dust_moments(3,iSPH) * 12.*amu/mass_per_H enddo - write(*,*) "Dust mass in phantom dump is ", real(mass), "Msun" + write(*,*) "Dust mass in hydro model is ", real(mass), "Msun" lvariable_dust = .true. allocate(grainsize_f(n_grains_tot),gsize(n_grains_tot),dN_ds(n_grains_tot),& @@ -490,7 +440,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g ! this is required if ndusttypes == 1, but I do it in any case allocate(a_SPH(ndusttypes+1),log_a_SPH(ndusttypes+1),rho_dust(ndusttypes+1)) - write(*,*) "Found the following grain sizes in SPH calculation:" + write(*,*) "Found the following grain sizes in hydro calculations:" do l=1, ndusttypes if (dust_pop(1)%porosity > tiny_real) then if (l==1) call warning("Grain sizes are adjusted for porosity") @@ -996,80 +946,6 @@ subroutine read_ascii_SPH_file(iunit,filename,x,y,z,h,vx,vy,vz,T_gas,massgas,rho end subroutine read_ascii_SPH_file - -subroutine test_voro_star(x,y,z,h,vx,vy,vz,T_gas,massgas,rhogas,rhodust,particle_id,ndusttypes,n_sph) - use naleat, only : seed, stream, gtype -#include "sprng_f.h" - - real :: rand, rand2, rand3 - real(dp), intent(out), dimension(:), allocatable :: x,y,z,h,rhogas,massgas,vx,vy,vz,T_gas - real(dp), intent(out), dimension(:,:), allocatable :: rhodust - integer, intent(out), dimension(:), allocatable :: particle_id - integer, intent(out) :: ndusttypes, n_SPH - real, parameter :: rmi = 2.2, rmo = 3.0 !unit of etoile(1)%r - integer :: alloc_status, i, id - real(kind=dp) :: rr, tt, pp, sintt, costt, vol - - - !force nb_proc to 1 here, but we don't want to set nb_proc = 1 for the rest of the calc. - if (allocated(stream)) deallocate(stream) - allocate(stream(1)) - stream = 0.0 - do i=1, 1 - !write(*,*) gtype, i-1,nb_proc,seed,SPRNG_DEFAULT - !init_sprng(gtype, i-1,nb_proc,seed,SPRNG_DEFAULT) - stream(i) = init_sprng(gtype, i-1,nb_proc,seed,SPRNG_DEFAULT) - enddo - - lignore_dust = .true. - ndusttypes = 0 - llimits_file = .false. - limits_file = "" - lascii_sph_file = .false. - lphantom_file = .false. - lgadget2_file = .false. - - if (allocated(x)) then - deallocate(x,y,z,h,vx,vy,vz,T_gas,massgas,rhogas,rhodust,particle_id) - endif - - n_sph = 100000 - - alloc_status = 0 - allocate(x(n_SPH),y(n_SPH),z(n_SPH),h(n_SPH),massgas(n_SPH),rhogas(n_SPH),particle_id(n_sph), & - vx(n_sph), vy(n_sph), vz(n_sph), T_gas(n_sph), stat=alloc_status) - if (alloc_status /=0) call error("Allocation error in phanton_2_mcfost") - - id = 1 - do i=1, n_sph - particle_id(i) = i - rand = sprng(stream(id)) - rand2 = sprng(stream(id)) - rand3 = sprng(stream(id)) - costt = (2*rand2-1) - sintt = sqrt(1.0 - costt*costt) - pp = pi * (2*rand3-1) - rr = ( rmi + (rmo - rmi) * rand )! * sintt*sintt - x(i) = rr * sintt * cos(pp) * etoile(1)%r - y(i) = rr * sintt * sin(pp) * etoile(1)%r - z(i) = rr * costt * etoile(1)%r - vol = (rand + 1) * 0.025 * etoile(1)%r**3 * (4.0/3.0) * pi - vx(i) = 0.0 - vy(i) = 0.0 - vz(i) = 0.0 - T_gas(i) = 1000.0 - massgas(i) = 1d-10 * (vol * AU_to_m**3) * kg_to_Msun!Msun - rhogas(i) = massgas(i) - h(i) = 3.0 * rmo * etoile(1)%r - enddo - - - deallocate(stream) - - -return -end subroutine - !********************************************************* ! subroutine read_Judith_file(iunit,filename,x,y,z,h,massgas,rhogas,rhodust,ndusttypes,n_SPH,ierr) From 7feba79feb7fd5e0a46652ecbc6fd77bc3e7a64d Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 17:35:15 +1000 Subject: [PATCH 10/13] Removing asciii SPH interface --- src/SPH2mcfost.f90 | 53 ++------------------------------------------- src/init_mcfost.f90 | 15 ------------- src/parameters.f90 | 2 +- 3 files changed, 3 insertions(+), 67 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 0eaf96c4a..78bdf7ba9 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -84,11 +84,6 @@ subroutine setup_SPH2mcfost(extra_heating) write(*,*) "Performing gadget2mcfost setup" write(*,*) "Reading Gadget-2 density file: "//trim(density_files(1)) call read_gadget2_file(iunit,density_files(1), x,y,z,h,massgas,rho,rhodust,ndusttypes,n_SPH,ierr) - else if (lascii_SPH_file) then - write(*,*) "Performing SPH2mcfost setup" - write(*,*) "Reading SPH density file: "//trim(density_files(1)) - call read_ascii_SPH_file(iunit,density_files(1), x,y,z,h,vx,vy,vz,T_gas,massgas,rho,rhodust,& - particle_id, ndusttypes,n_SPH,ierr) else if (ldensity_file) then call read_Voronoi_fits_file(density_files(1), x,y,z,h,vx,vy,vz,particle_id,massgas,n_SPH) else @@ -763,6 +758,7 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga end subroutine hydro_to_Voronoi_atomic + !********************************************************* subroutine test_duplicate_particles(n_SPH, particle_id, x,y,z, massgas,massdust,rho,rhodust, is_ghost) ! Filtering particles at the same position @@ -818,7 +814,7 @@ subroutine test_duplicate_particles(n_SPH, particle_id, x,y,z, massgas,massdust, if (nkill > 0) then write(*,*) - write(*,*) nkill, "SPH particles were flagged as ghosts and merged" + write(*,*) nkill, "hydro sites were flagged as ghosts and merged" write(*,*) endif @@ -903,51 +899,6 @@ end subroutine delete_Hill_sphere !********************************************************* - subroutine read_ascii_SPH_file(iunit,filename,x,y,z,h,vx,vy,vz,T_gas,massgas,rhogas,rhodust,particle_id, ndusttypes,n_SPH,ierr) - - integer, intent(in) :: iunit - character(len=*), intent(in) :: filename - real(dp), intent(out), dimension(:), allocatable :: x,y,z,h,rhogas,massgas,vx,vy,vz,T_gas - real(dp), intent(out), dimension(:,:), allocatable :: rhodust - integer, intent(out), dimension(:), allocatable :: particle_id - integer, intent(out) :: ndusttypes, n_SPH,ierr - - integer :: syst_status, alloc_status, ios, i - character(len=512) :: cmd - - ierr = 0 - - cmd = "wc -l "//trim(filename)//" > ntest.txt" - call appel_syst(cmd,syst_status) - open(unit=1,file="ntest.txt",status="old") - read(1,*) n_SPH - close(unit=1) - ndusttypes = 1 - - write(*,*) "n_SPH read_test_ascii_file = ", n_SPH - - alloc_status = 0 - allocate(x(n_SPH),y(n_SPH),z(n_SPH),h(n_SPH),massgas(n_SPH),rhogas(n_SPH),particle_id(n_sph), rhodust(ndusttypes,n_SPH), & - vx(n_sph), vy(n_sph), vz(n_sph), T_gas(n_sph), stat=alloc_status) - if (alloc_status /=0) call error("Allocation error in phanton_2_mcfost") - - open(unit=1, file=filename, status='old', iostat=ios) - do i=1, n_SPH - read(1,*) x(i), y(i), z(i), h(i), T_gas(i), massgas(i), vx(i), vy(i), vz(i) - rhogas(i) = massgas(i) - particle_id(i) = i - enddo - - write(*,*) "MinMax=", minval(massgas), maxval(massgas) - - write(*,*) "Using stars from mcfost parameter file" - - return - - end subroutine read_ascii_SPH_file - - !********************************************************* - ! subroutine read_Judith_file(iunit,filename,x,y,z,h,massgas,rhogas,rhodust,ndusttypes,n_SPH,ierr) ! ! !The 3D cube is 20 cells in colatitude direction, 215 radially, and 680 azimuthally. diff --git a/src/init_mcfost.f90 b/src/init_mcfost.f90 index 316df9c7c..a13191924 100644 --- a/src/init_mcfost.f90 +++ b/src/init_mcfost.f90 @@ -133,7 +133,6 @@ subroutine set_default_variables() lforce_Mgas = .false. lforce_SPH_amin = .false. lforce_SPH_amax = .false. - lascii_SPH_file = .false. lgadget2_file=.false. llimits_file = .false. lsigma_file = .false. @@ -959,16 +958,6 @@ subroutine initialisation_mcfost() density_files(i) = s enddo if (.not.llimits_file) limits_file = "phantom.limits" - case("-ascii_SPH") - i_arg = i_arg + 1 - lascii_SPH_file = .true. - lVoronoi = .true. - l3D = .true. - call get_command_argument(i_arg,s) - allocate(density_files(1)) - density_files(1) = s - i_arg = i_arg + 1 - if (.not.llimits_file) limits_file = "phantom.limits" case("-SPH_amin") lforce_SPH_amin = .true. i_arg = i_arg + 1 @@ -1586,10 +1575,6 @@ subroutine initialisation_mcfost() write(*,*) 'Input file read successfully' -! if ((lsphere_model.or.lmodel_1d).and.(lascii_sph_file.or.lphantom_file)) then -! call error("Cannot use Phantom and MHD files at the same time presently.") -! end if - ! Correction sur les valeurs du .para if (lProDiMo) then if ((lsed_complete).or.(tab_wavelength(1:7) /= "ProDiMo")) then diff --git a/src/parameters.f90 b/src/parameters.f90 index 59d128ee6..fa3295a52 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -133,7 +133,7 @@ module parametres logical :: lapprox_diffusion, lcylindrical, lspherical, llinear_rgrid, lVoronoi, is_there_disk, lno_backup logical :: laverage_grain_size, lisotropic, lno_scattering, lqsca_equal_qabs, lonly_diff_approx, lforce_diff_approx logical :: ldensity_file, lsigma_file, lvelocity_file, lphantom_file, lphantom_multi, lphantom_avg - logical :: lgadget2_file, lascii_SPH_file, llimits_file, lforce_SPH_amin, lforce_SPH_amax, lmcfost_lib + logical :: lgadget2_file, llimits_file, lforce_SPH_amin, lforce_SPH_amax, lmcfost_lib logical :: lweight_emission, lcorrect_density, lProDiMo2mcfost, lProDiMo2mcfost_test, lastrochem, lML logical :: lspot, lforce_PAH_equilibrium, lforce_PAH_out_equilibrium, lchange_Tmax_PAH, lISM_heating, lcasa, lJy, lforce_Mgas integer :: ISR_model ! 0 : no ISM radiation field, 1 : ProDiMo, 2 : Bate & Keto From 07c40b5042c282ffd599246dc33608df5d22aa0d Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 17:40:48 +1000 Subject: [PATCH 11/13] A bit more cleaning --- src/SPH2mcfost.f90 | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 78bdf7ba9..ba250996f 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -638,12 +638,10 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g end subroutine SPH_to_Voronoi - - !**************************************************************************** - ! copy additional parameters from hydro code needed for atomic line transfer !**************************************************************************** - subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massgas,mask) + subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massgas,mask) + ! copy additional parameters from hydro code needed for atomic line transfer ! ************************************************************************************ ! ! n_sph : number of points in the input model ! particle_id : index of the particle / cell centre in the input model @@ -659,8 +657,6 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga use disk_physics, only : compute_othin_sublimation_radius use mem use elements_type, only : wght_per_H, read_abundance - ! use mhd2mcfost, only : alloc_atomrt_grid, nHtot, ne, & - ! v_char, lmagnetized, vturb, T, icompute_atomRT, lcalc_ne use grid, only : alloc_atomrt_grid, nHtot, ne, v_char, lmagnetized, vturb, T, icompute_atomRT, lcalc_ne, & check_for_zero_electronic_density @@ -729,12 +725,7 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga write(*,*) "Read ", size(pack(icompute_atomRT,mask=icompute_atomRT==0)), " transparent zones" write(*,*) "Read ", size(pack(icompute_atomRT,mask=icompute_atomRT<0)), " dark zones" - ! if (lmagnetized) then - ! - ! endif - - - call check_for_zero_electronic_density + call check_for_zero_electronic_density() write(*,*) "Maximum/minimum velocities in the model (km/s):" write(*,*) "|Vx|", vxmax*1d-3, vxmin*1d-3 @@ -748,14 +739,16 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga write(*,*) maxval(vturb)/1d3, minval(vturb, mask=icompute_atomRT>0)/1d3 write(*,*) "Maximum/minimum Temperature in the model (K):" - write(*,*) MAXVAL(T), MINVAL(T,mask=icompute_atomRT>0) + write(*,*) maxval(T), minval(T,mask=icompute_atomRT>0) write(*,*) "Maximum/minimum Hydrogen total density in the model (m^-3):" - write(*,*) MAXVAL(nHtot), MINVAL(nHtot,mask=icompute_atomRT>0) + write(*,*) maxval(nHtot), minval(nHtot,mask=icompute_atomRT>0) if (.not.lcalc_ne) then write(*,*) "Maximum/minimum ne density in the model (m^-3):" - write(*,*) MAXVAL(ne), MINVAL(ne,mask=icompute_atomRT>0) + write(*,*) maxval(ne), minval(ne,mask=icompute_atomRT>0) endif + return + end subroutine hydro_to_Voronoi_atomic !********************************************************* From a3022092b7752384f167a64dead7906173d943f0 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Tue, 28 May 2024 18:18:35 +1000 Subject: [PATCH 12/13] Adding an artificial smoothling length when reading Voronoi fits file --- src/density.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/density.f90 b/src/density.f90 index 5c605c75d..d7db86d5d 100644 --- a/src/density.f90 +++ b/src/density.f90 @@ -1773,6 +1773,9 @@ subroutine read_Voronoi_fits_file(filename, x,y,z,h,vx,vy,vz,particle_id, massga call error("cannot read bitpix in fits file") endif + ! Defining a smoothing length + h = 0.02 * sqrt(x*x+y*y+z*z) + do i=1,n_points particle_id(i) = i enddo From c100edc4760b51bbdf56a3335dbddc49fe7b238a Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Thu, 30 May 2024 10:07:48 +1000 Subject: [PATCH 13/13] Bug fix: non-allocated arrays --- src/SPH2mcfost.f90 | 5 ++++- src/dust_ray_tracing.f90 | 6 ------ 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index ba250996f..d9dc565ac 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -86,6 +86,8 @@ subroutine setup_SPH2mcfost(extra_heating) call read_gadget2_file(iunit,density_files(1), x,y,z,h,massgas,rho,rhodust,ndusttypes,n_SPH,ierr) else if (ldensity_file) then call read_Voronoi_fits_file(density_files(1), x,y,z,h,vx,vy,vz,particle_id,massgas,n_SPH) + ldust_moments = .false. + ndusttypes=0 else call error("Unknown SPH structure.") endif @@ -116,7 +118,8 @@ subroutine setup_SPH2mcfost(extra_heating) if (lemission_atom) then call hydro_to_Voronoi_atomic(n_SPH,T_gas,vturb,massgas,mass_ne_on_massgas,atomic_mask) endif - deallocate(massgas,rho) + deallocate(massgas) + if (allocated(rho)) deallocate(rho) if (allocated(vturb)) deallocate(vturb) if (allocated(mass_ne_on_massgas)) deallocate (mass_ne_on_massgas) if (allocated(atomic_mask)) deallocate(atomic_mask) diff --git a/src/dust_ray_tracing.f90 b/src/dust_ray_tracing.f90 index f55b28adc..0489af6e8 100644 --- a/src/dust_ray_tracing.f90 +++ b/src/dust_ray_tracing.f90 @@ -290,12 +290,6 @@ subroutine init_directions_ray_tracing() enddo enddo - - do ibin=1, RT_n_incl - write(*,*) ibin, tab_RT_incl(ibin), -tab_uv_rt(ibin), tab_w_rt(ibin) - enddo - - if (.not. allocated(star_position)) then allocate(star_position(n_etoiles,RT_n_incl,RT_n_az,2), star_vr(n_etoiles,RT_n_incl,RT_n_az), stat=alloc_status) if (alloc_status > 0) call error('Allocation error star_position')