fast density2 #159
fast density2 #159
Conversation
- Correcting dust_density access: Replaced direct indexing with conditional logic or the p_k/p_l pointers in dust_prop.f90, thermal_emission.f90, output.f90, io_prodimo.f90, and read1d_models.f90. - Updating OpenMP directives: Ensured thread safety for the new pointers in thermal_emission.f90. - Fixing Normalization Logic: Corrected density.f90 to use volume integrals and properly handle grain size distributions when lvariable_dust is true.
…ust checks outside loops
This reverts commit a26507d.
Summary of ChangesHello @cpinte, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed! This pull request focuses on a significant refactoring and optimization of dust density calculations within the codebase. The primary goal is to enhance performance and maintainability by standardizing variable names, optimizing data structures for dust properties, and streamlining density-related computations. These changes aim to improve the efficiency of various modules that interact with dust density, such as those involved in radiative transfer, SPH model conversion, and thermal emission. Highlights
Using Gemini Code AssistThe full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips. Invoking Gemini You can request assistance from Gemini at any point by creating a comment using either
Customization To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a Limitations & Feedback Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here. You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension. Footnotes
|
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request primarily refactors the dust density representation across the codebase, renaming the variable densite_pouss to dust_density and introducing a nbre_grains factor to explicitly handle particle counts. The dust_density array's first dimension now dynamically adapts based on whether lvariable_dust is true (n_grains_tot) or false (n_pop), optimizing memory usage and calculations when dust properties are uniform. This change impacts numerous subroutines in files like ML_prodimo.f90, SPH2mcfost.f90, density.f90, disk_physics.f90, dust_prop.f90, dust_ray_tracing.f90, atom_transfer.f90, input.f90, io_prodimo.f90, mem.f90, optical_depth.f90, output.f90, read_athena++.f90, read_fargo3d.f90, read_idefix.f90, read_pluto.f90, read_spherical_grid.f90, and thermal_emission.f90. Additionally, several instances of read1d_models.f90 were renamed to read_1d_models.f90 for consistency, and a debugging write statement was removed from density.f90. Review comments highlighted the need to remove debugging write statements and a commented-out stop in src/density.f90, and to remove unused variable declarations (p_k) in src/thermal_emission.f90 and src/dust_prop.f90 to improve code clarity. An uninitialized variable k in src/dust_prop.f90 was also identified as a potential source of undefined behavior.
|
|
||
| integer :: p_k_val, k | ||
|
|
||
| p_k_val = merge(k, grain(k)%pop, lvariable_dust) |
There was a problem hiding this comment.
The variable k is used here before it has been initialized. This will lead to undefined behavior or a crash. This line seems to be a leftover from a copy-paste and should be removed, as p_k_val is correctly assigned inside the loops later on.
| ! 1 grain -> 100 grains | ||
|
|
||
| ! 0.37s -> 1.9s | ||
| write(*,*) "A" | ||
|
|
||
| call define_gas_density() | ||
| ! 0.8s after OpenMP -> 2.4s | ||
| write(*,*) "B" | ||
|
|
||
| call define_dust_density() | ||
| ! 1.4s -> 40.7s | ||
| write(*,*) "C" | ||
| !stop | ||
|
|
||
|
|
||
|
|
There was a problem hiding this comment.
These lines appear to be for debugging purposes, including performance timings, write statements, and a commented-out stop. They should be removed before merging.
call define_gas_density()
call define_dust_density()
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method | ||
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used within this subroutine. It should be removed to improve code clarity.
integer :: l, l1, l2, Ti, T1, T2, k, heating_method
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method | ||
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used within this subroutine. It should be removed to improve code clarity.
integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda | ||
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used within this subroutine. It should be removed to improve code clarity.
integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda
There was a problem hiding this comment.
Pull request overview
This PR refactors dust density handling in mcfost (a radiative transfer code), primarily renaming densite_pouss to dust_density and optimizing performance by conditionally allocating arrays based on whether dust properties vary spatially (lvariable_dust flag). The PR includes bug fixes for grain indexing, module name corrections (read1d_models → read_1d_models), and comment/spelling updates.
Changes:
- Refactored dust density array from
densite_pouss(n_grains_tot, n_cells)todust_density(p_n_grains, n_cells)where p_n_grains is either n_grains_tot or n_pop depending on lvariable_dust flag - Updated all references throughout the codebase to use the new array with proper indexing and multiplication by
nbre_grains(k) - Fixed module name from read1d_models to read_1d_models for consistency
Reviewed changes
Copilot reviewed 23 out of 23 changed files in this pull request and generated 10 comments.
Show a summary per file
| File | Description |
|---|---|
| src/density.f90 | Core refactoring of dust density calculations and array structure; includes debugging code that needs removal |
| src/thermal_emission.f90 | Updated grain selection and emission calculations; contains critical pointer logic bug |
| src/dust_prop.f90 | Updated opacity calculations with new dust_density array; contains initialization bug |
| src/dust_ray_tracing.f90 | Updated thermal emission calculations with proper grain indexing |
| src/mem.f90 | Modified memory allocation to conditionally size dust_density based on lvariable_dust |
| src/output.f90 | Updated output routines; contains array bounds bug in taille_moyenne_grains |
| src/read_spherical_grid.f90 | Sets lvariable_dust=false for uniform dust models |
| src/read_*.f90 | Updated all model readers to use dust_density instead of densite_pouss |
| src/io_prodimo.f90 | Updated ProDiMo output with proper lvariable_dust handling |
| src/SPH2mcfost.f90 | Updated SPH reader with correct dust density normalization |
| src/Makefile | Updated module name reference |
Comments suppressed due to low confidence (1)
src/read_1d_models.f90:257
- Inconsistent indentation: this line uses spaces while the surrounding code in this function uses tabs. The indentation style should be consistent throughout the file.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| write(*,*) "A" | ||
|
|
||
| call define_gas_density() | ||
| ! 0.8s after OpenMP -> 2.4s | ||
| write(*,*) "B" | ||
|
|
||
| call define_dust_density() | ||
| ! 1.4s -> 40.7s | ||
| write(*,*) "C" | ||
| !stop |
There was a problem hiding this comment.
This appears to be debugging code that should be removed before merging. The write statements and commented-out stop statement are not appropriate for production code.
| write(*,*) "A" | |
| call define_gas_density() | |
| ! 0.8s after OpenMP -> 2.4s | |
| write(*,*) "B" | |
| call define_dust_density() | |
| ! 1.4s -> 40.7s | |
| write(*,*) "C" | |
| !stop | |
| call define_gas_density() | |
| ! 0.8s after OpenMP -> 2.4s | |
| call define_dust_density() | |
| ! 1.4s -> 40.7s |
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method | ||
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used in this function. It should be removed to avoid confusion.
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method, p_k | |
| integer :: l, l1, l2, Ti, T1, T2, k, heating_method |
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda | ||
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used in this function. It should be removed to avoid confusion.
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda, p_k | |
| integer :: l, l1, l2, T_int, T1, T2, k, lambda0, ilambda |
| CDF = CDF + C_abs(k,lambda) * dust_density(p_k,icell) * nbre_grains(k) | ||
| if (CDF > prob) exit | ||
| enddo | ||
| else ! We start from the end of the grain size distribution | ||
| prob = (1.0-aleat) * norm | ||
| CDF = 0.0 | ||
| do k=kend, kstart, -1 | ||
| CDF = CDF + C_abs(k,lambda) * densite_pouss(k,icell) | ||
| CDF = CDF + C_abs(k,lambda) * dust_density(p_k,icell) * nbre_grains(k) | ||
| if (CDF > prob) exit |
There was a problem hiding this comment.
The pointer p_k is being dereferenced inside the loop, but when lvariable_dust is false, p_k points to k1 which is always 1. This means the code accesses dust_density(1, icell) but uses k as the index for C_abs(k, lambda) and nbre_grains(k). This mismatch will cause incorrect results when lvariable_dust is false. The intended behavior appears to be accessing dust_density(grain(k)%pop, icell) when lvariable_dust is false, which should be computed inside the loop rather than using a pointer.
| real(kind=dp), allocatable :: T_tmp(:,:,:), ne_tmp(:,:,:), vt_tmp(:,:,:) | ||
| real(kind=dp) :: mass | ||
|
|
||
| lvariable_dust = .false. ! the model has the same dust everywhere |
There was a problem hiding this comment.
The lvariable_dust flag is set to false but there's no corresponding initialization to true in other reader files. This could lead to inconsistent behavior depending on which reader is used. Consider adding a comment explaining when this should be true versus false, or ensure all reader files explicitly set this flag.
| lvariable_dust = .false. ! the model has the same dust everywhere | |
| ! Spherical models assume uniform dust properties throughout the grid. | |
| ! Set lvariable_dust = .false. here; reader routines that support | |
| ! spatially varying dust must explicitly set lvariable_dust = .true. | |
| lvariable_dust = .false. |
| do l=1, n_grains_tot | ||
| a_moyen(icell) = a_moyen(icell) + densite_pouss(l,icell) * r_grain(l)**2 | ||
| somme = somme + densite_pouss(l,icell) | ||
| a_moyen(icell) = a_moyen(icell) + dust_density(l,icell) * nbre_grains(l) * r_grain(l)**2 | ||
| somme = somme + dust_density(l,icell) * nbre_grains(l) | ||
| enddo | ||
| a_moyen(icell) = a_moyen(icell) / somme | ||
| a_moyen(icell) = sqrt(a_moyen(icell) / somme) |
There was a problem hiding this comment.
When lvariable_dust is false, dust_density is allocated with first dimension n_pop (not n_grains_tot), but this loop iterates over all grains (1 to n_grains_tot). This will cause an array bounds error when l > n_pop. The code should use the grain population index when lvariable_dust is false.
| integer, intent(inout) :: lambda | ||
|
|
||
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method | ||
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method, p_k |
There was a problem hiding this comment.
The variable p_k is declared but never used in this function. It should be removed to avoid confusion.
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method, p_k | |
| integer :: l, l1, l2, T_int, T1, T2, k, kmin, kmax, lambda0, ilambda, heating_method |
| if (lTemp) then | ||
| spectre_emission_cumul(0) = 0.0 | ||
| ! Fonction de r�partition �mssion | ||
| ! Fonction de répartition émssion |
There was a problem hiding this comment.
Typo in the French comment: "émssion" should be "émission" (missing 'i').
| ! Fonction de répartition émssion | |
| ! Fonction de répartition émission |
|
|
||
| integer :: p_k_val, k | ||
|
|
||
| p_k_val = merge(k, grain(k)%pop, lvariable_dust) |
There was a problem hiding this comment.
The variable p_k_val is initialized before the loop using an uninitialized value of k on line 1297. This line should be removed since p_k_val is correctly computed inside each loop (lines 1307 and 1316).
| p_k_val = merge(k, grain(k)%pop, lvariable_dust) |
| if (l_RE(k,icell)) CDF = CDF + C_abs(k,lambda) * dust_density(p_k,icell) * nbre_grains(k) | ||
| if (CDF > prob) exit | ||
| enddo | ||
| else ! We start from the end of the grain size distribution | ||
| prob = (1.0-aleat) * norm | ||
| CDF = 0.0 | ||
| do k=kend, kstart, -1 | ||
| if (l_RE(k,icell)) CDF = CDF + C_abs(k,lambda) * densite_pouss(k,icell) | ||
| if (l_RE(k,icell)) CDF = CDF + C_abs(k,lambda) * dust_density(p_k,icell) * nbre_grains(k) | ||
| if (CDF > prob) exit |
There was a problem hiding this comment.
Same issue as lines 2005-2013: when lvariable_dust is false, p_k points to a constant k1=1, which doesn't correctly access the dust density. The code should compute the proper index inside the loop based on grain(k)%pop when lvariable_dust is false.
1 participant
Type of PR:
Bug fix / new physics / modification to existing code / other
Description:
Give a short description of the change.
Testing:
Describe how you have tested the change
Did you run the botscheck that the code and comments follow the code of conduct? yes/no
Did you update relevant documentation in the docs directory? yes/no