Skip to content

fast density2 #159

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cpinte merged 27 commits into from
Mar 24, 2026
Merged

fast density2 #159

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
63da8e9
First try at splitting dust density
cpinte Sep 20, 2025
fbe4e72
Bug fix: incorrect index when sublimating dust
cpinte Sep 20, 2025
de9c0ab
Merge branch 'main' into fast_density
cpinte Dec 15, 2025
e62f98a
Trying to fix unit tests (incorrect dust grain index)
cpinte Jan 5, 2026
5fdcffa
Fixing dust density normalisation
cpinte Jan 6, 2026
96a62c8
Optimise loops by hoisting dust_density access and moving lvariable_d…
cpinte Jan 6, 2026
96b4b46
Bug fix: incorrect grain index for SPH models
cpinte Jan 7, 2026
61ad4d4
Cleaning Benjamin's model input
cpinte Jan 8, 2026
983e48d
Simplying output with new density arrays
cpinte Jan 8, 2026
3c7888e
Simplifying dust_ray_tracing
cpinte Jan 10, 2026
0522883
Simplifying dust_prop
cpinte Jan 10, 2026
a26507d
Trying in debug mode instead of release
cpinte Jan 10, 2026
5b94005
Merge branch 'fast_density' of github.com:cpinte/mcfost into fast_den…
cpinte Jan 10, 2026
f9dbed1
Use integer index instead of pointer to avoid ifort OpenMP issues
cpinte Jan 10, 2026
094c812
Revert "Trying in debug mode instead of release"
cpinte Jan 11, 2026
46ff4dc
Merge branch 'main' into fast_density
cpinte Jan 11, 2026
78d5f23
Merge branch 'main' into fast_density
cpinte Jan 11, 2026
652ffcd
Simplifying thermal_emission
cpinte Jan 11, 2026
7f3f08a
Simplifying density
cpinte Jan 11, 2026
eaf243e
Cleaning code, removing references to densite_pouss
cpinte Jan 12, 2026
9d75a72
Fixing bug with multiple populations
cpinte Jan 13, 2026
18f5e3c
Merge branch 'main' into fast_density2
cpinte Jan 20, 2026
609adf2
Merge branch 'main' into fast_density2
cpinte Jan 21, 2026
d46b1c9
Swapping loops in normalize_dust_density + dust density per zone
cpinte Mar 21, 2026
6da3739
Merge branch 'fast_density2' of github.com:cpinte/mcfost into fast_de…
cpinte Mar 21, 2026
86b2611
Bug fix: zone-based dust_density allocation and physical scaling fixes
cpinte Mar 22, 2026
ae9549e
Merge branch 'main' into fast_density2
cpinte Mar 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions src/ML_prodimo.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,9 @@ end subroutine save_J_ML

subroutine xgb_compute_features()

use grains, only : grain, r_grain, n_grains_tot
use grains, only : grain, r_grain, n_grains_tot, nbre_grains
use optical_depth, only : compute_column
use density, only : densite_pouss
use density, only : dust_density
use cylindrical_grid, only : r_grid, z_grid
use molecular_emission, only : densite_gaz
use temperature, only : Tdust
Expand All @@ -166,12 +166,12 @@ subroutine xgb_compute_features()
!--- Moments de la distribution de grain
mask_not_PAH(:) = .not.grain(:)%is_PAH
do icell=1, n_cells
N = sum(densite_pouss(:,icell),mask=mask_not_PAH)
N = sum(dust_density(:,icell) * nbre_grains(:),mask=mask_not_PAH)
Copy link
Copy Markdown

@gemini-code-assist gemini-code-assist Bot Mar 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

high: This line calculates N by summing the product of dust_density and nbre_grains. It's crucial to ensure that these arrays are properly aligned and that this calculation accurately represents the intended physical quantity (number density of grains). Consider adding a comment explaining the physical units of the resulting N and how it relates to the input arrays.

N_grains(0,icell) = N
if (N > 0) then
N_grains(1,icell) = sum(densite_pouss(:,icell) * r_grain(:),mask=mask_not_PAH) / N
N_grains(2,icell) = sum(densite_pouss(:,icell) * r_grain(:)**2,mask=mask_not_PAH) / N
N_grains(3,icell) = sum(densite_pouss(:,icell) * r_grain(:)**3,mask=mask_not_PAH) / N
N_grains(1,icell) = sum(dust_density(:,icell) * nbre_grains(:) * r_grain(:),mask=mask_not_PAH) / N
N_grains(2,icell) = sum(dust_density(:,icell) * nbre_grains(:) * r_grain(:)**2,mask=mask_not_PAH) / N
N_grains(3,icell) = sum(dust_density(:,icell) * nbre_grains(:) * r_grain(:)**3,mask=mask_not_PAH) / N
else
N_grains(1,icell) = 0.0
N_grains(2,icell) = 0.0
Expand Down
2 changes: 1 addition & 1 deletion src/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@ SOURCES = mcfost_env.f90 parameters.f90 constants.f90 healpix_mod.f90 sha.f90 me
read_fargo3d.f90 hdf5_utils.f90 utils_hdf5.f90 read_athena++.f90 \
readVTK.f90 read_idefix.f90 read_pluto.f90 voigts.f90 \
coated_sphere.f90 dust_prop.f90 molecular_emission.f90 PAH.f90 \
read1d_models.f90 read_spherical_grid.f90 \
read_1d_models.f90 read_spherical_grid.f90 \
input.f90 benchmarks.f90 atom_type.f90 \
wavelengths_gas.f90 broad.f90 \
read_param.f90 dust_ray_tracing.f90 uplow.f90 abo.f90 \
Expand Down
38 changes: 19 additions & 19 deletions src/SPH2mcfost.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
! mask : integer array, 1 if a SPH particle will be made transparent, 2 if deleted before tesselation
! ************************************************************************************ !
use Voronoi_grid
use density, only : densite_gaz, masse_gaz, densite_pouss, masse
use density, only : densite_gaz, masse_gaz, dust_density, masse
use grains, only : n_grains_tot, M_grain
use disk_physics, only : compute_othin_sublimation_radius
use mem
Expand Down Expand Up @@ -367,7 +367,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
!$omp default(none) &
!$omp private(icell,iSPH,use_single_grain,rhoi,ki,err,ierr,lambsol,rho_monomers,a,i,norm,mdust) &
!$omp shared(n_cells,Voronoi,densite_gaz,n_grains_tot,r_grain,mass_factor,dN_ds,N_monomers) &
!$omp shared(densite_pouss,dust_moments,masse_gaz,volume) &
!$omp shared(dust_density,nbre_grains,dust_moments,masse_gaz,volume) &
!$omp reduction(+:N_pb)
!$omp do schedule(dynamic,1)
do icell=1,n_cells
Expand All @@ -389,19 +389,19 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
endif
enddo

densite_pouss(:,icell) = rho_monomers(:) * dN_ds(:)
dust_density(:,icell) = rho_monomers(:) * dN_ds(:)

! Simple approximation : we assume 1 single grain size
if (use_single_grain) then
densite_pouss(:,icell) = 0._dp
dust_density(:,icell) = 0._dp
if (dust_moments(1,iSPH) > tiny_dp) then
a = a0 * dust_moments(2,iSPH)/dust_moments(1,iSPH)
i = locate(1.0_dp*r_grain(:),a)
densite_pouss(i,icell) = 1.0
dust_density(i,icell) = 1.0
endif
endif
else ! iSPH == 0, star
densite_pouss(:,icell) = 0._dp
dust_density(:,icell) = 0._dp
endif
enddo ! icell
!$omp end do
Expand All @@ -417,7 +417,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g

mass = 0.0_dp
do l=1,n_grains_tot
mass=mass + (densite_pouss(l,icell) *1.0_dp) * M_grain(l)
mass=mass + (dust_density(l,icell) * nbre_grains(l) *1.0_dp) * M_grain(l)
enddo !l
mass = mass * volume(icell)

Expand All @@ -426,7 +426,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g

if (mass > tiny_dp) then
factor = mdust/ mass
densite_pouss(:,icell) = densite_pouss(:,icell) * factor
dust_density(:,icell) = dust_density(:,icell) * factor
endif
enddo !icell
call normalize_dust_density(mdust_tot) ! this should only calculates the array masse
Expand Down Expand Up @@ -532,17 +532,17 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
l=1
do k=1,n_grains_tot
if (r_grain(k) < a_SPH(1)) then ! small grains
densite_pouss(k,icell) = rho_dust(1)
dust_density(k,icell) = rho_dust(1) / nbre_grains(k)
Copy link
Copy Markdown

@gemini-code-assist gemini-code-assist Bot Mar 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

medium: The division by nbre_grains(k) here and in the following lines seems intended to normalize the dust density. It's important to verify that this normalization is consistent with the overall dust mass budget and that it doesn't introduce any unintended side effects, especially when dealing with multiple dust populations. Add a comment explaining the purpose of this division and how it affects the final dust density.

else if (r_grain(k) > a_SPH(ndusttypes+1)) then ! large grains
densite_pouss(k,icell) = rho_dust(ndusttypes+1)
dust_density(k,icell) = rho_dust(ndusttypes+1) / nbre_grains(k)
else ! interpolation
if (r_grain(k) > a_sph(l+1)) l = l+1
f = (log(r_grain(k))-log_a_sph(l))/(log_a_sph(l+1)-log_a_sph(l))
densite_pouss(k,icell) = rho_dust(l) + f * (rho_dust(l+1) - rho_dust(l))
dust_density(k,icell) = (rho_dust(l) + f * (rho_dust(l+1) - rho_dust(l))) / nbre_grains(k)
endif
enddo !k
else ! iSPH == 0, star
densite_pouss(:,icell) = 0.
dust_density(:,icell) = 0.
endif
enddo ! icell

Expand All @@ -556,14 +556,14 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g

do icell=1,n_cells
masse(icell) = 0.
dust_density(1,icell) = densite_gaz(icell)
Copy link
Copy Markdown

@gemini-code-assist gemini-code-assist Bot Mar 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

medium: It seems that dust_density(1,icell) is assigned the value of densite_gaz(icell). This might be valid only for specific cases where dust and gas are assumed to be well-mixed. Add a comment to clarify the assumption behind this assignment and whether it holds true for all simulation scenarios.

do k=1,n_grains_tot
densite_pouss(k,icell) = densite_gaz(icell) * nbre_grains(k)
masse(icell) = masse(icell) + densite_pouss(k,icell) * M_grain(k) * volume(icell)
masse(icell) = masse(icell) + dust_density(1,icell) * nbre_grains(k) * M_grain(k) * volume(icell)
enddo
enddo
masse(:) = masse(:) * AU3_to_cm3
f = 1./disk_zone(1)%gas_to_dust * sum(masse_gaz)/sum(masse)
densite_pouss(:,:) = densite_pouss(:,:) * f
dust_density(:,:) = dust_density(:,:) * f
masse(:) = masse(:) * f
endif ! ndusttypes == 0

Expand Down Expand Up @@ -816,7 +816,7 @@ end subroutine test_duplicate_particles
subroutine delete_masked_particles()

use Voronoi_grid
use density, only : densite_gaz, masse_gaz, densite_pouss, masse
use density, only : densite_gaz, masse_gaz, dust_density, masse

integer :: icell, k

Expand All @@ -827,7 +827,7 @@ subroutine delete_masked_particles()
masse_gaz(icell) = 0.
densite_gaz(icell) = 0.
masse(icell) = 0.
densite_pouss(:,icell) = 0.
dust_density(:,icell) = 0.
endif
enddo

Expand All @@ -842,7 +842,7 @@ end subroutine delete_masked_particles
subroutine delete_Hill_sphere()

use Voronoi_grid
use density, only : densite_gaz, masse_gaz, densite_pouss, masse
use density, only : densite_gaz, masse_gaz, dust_density, masse

integer :: istar, icell, n_delete
real(kind=dp) :: d2, r_Hill2, r_hill, dx, dy, dz
Expand Down Expand Up @@ -874,7 +874,7 @@ subroutine delete_Hill_sphere()
masse_gaz(icell) = 0.
densite_gaz(icell) = 0.
masse(icell) = 0.
densite_pouss(:,icell) = 0.
dust_density(:,icell) = 0.
n_delete = n_delete + 1
endif
enddo cell_loop
Expand Down
Loading
Loading