-
Notifications
You must be signed in to change notification settings - Fork 8
Add optional frequency dependence to PowerLaw #70
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.
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
Merged
Changes from all commits
Commits
Show all changes
5 commits
Select commit
Hold shift + click to select a range
9fe9174
Add optional nu dependence to PowerLaw
pdmullen 237229b
Add scattering version of powerlaw opacity.
RyanWollaeger c72fec9
Add additional optional parameters to photon power-law opacities.
RyanWollaeger 68c72ce
Re-order optional arguments.
RyanWollaeger a934fe9
+ Add test of scattering powerlaw (nearly same as absorption).
RyanWollaeger File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,25 @@ | ||
| ## Goal | ||
|
|
||
| Extend photon `PowerLaw` opacities to support an optional frequency power law | ||
| while preserving the legacy gray behavior when no frequency exponent is | ||
| provided. | ||
|
|
||
| ## Functional Plan | ||
|
|
||
| - Add optional `nu_exp` and `nu_ref` parameters to photon `PowerLaw`. | ||
| - Apply the frequency power law only to the monochromatic opacity and | ||
| monochromatic emissivity APIs. | ||
| - Preserve the existing legacy constructor calls by keeping the default | ||
| frequency behavior equivalent to `nu_exp = 0`. | ||
| - Keep the `NonCGSUnits` wrapper path working for the frequency-dependent | ||
| `PowerLaw`. | ||
| - Add focused regression coverage that checks: | ||
| - the monochromatic absorption coefficient matches the expected | ||
| `(\nu / \nu_ref)^{\nu_exp}` scaling | ||
| - the monochromatic emissivity remains `j_\nu = \alpha_\nu B_\nu` | ||
| - the same frequency dependence survives the `NonCGSUnits` wrapper | ||
|
|
||
| ## Scope Boundaries | ||
|
|
||
| - Do not add integrated `Emissivity` or `NumberEmissivity` support for | ||
| `nu_exp != 0`. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
125 changes: 125 additions & 0 deletions
125
singularity-opac/photons/powerlaw_s_opacity_photons.hpp
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,125 @@ | ||
| // ====================================================================== | ||
| // © 2026. Triad National Security, LLC. All rights reserved. This | ||
| // program was produced under U.S. Government contract | ||
| // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which | ||
| // is operated by Triad National Security, LLC for the U.S. | ||
| // Department of Energy/National Nuclear Security Administration. All | ||
| // rights in the program are reserved by Triad National Security, LLC, | ||
| // and the U.S. Department of Energy/National Nuclear Security | ||
| // Administration. The Government is granted for itself and others | ||
| // acting on its behalf a nonexclusive, paid-up, irrevocable worldwide | ||
| // license in this material to reproduce, prepare derivative works, | ||
| // distribute copies to the public, perform publicly and display | ||
| // publicly, and to permit others to do so. | ||
| // ====================================================================== | ||
|
|
||
| #ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ | ||
| #define SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ | ||
|
|
||
| // This file was partly copied from a file made in part with generative AI. | ||
|
|
||
| #include <cassert> | ||
| #include <cmath> | ||
| #include <cstdio> | ||
|
|
||
| #include <ports-of-call/portability.hpp> | ||
| #include <singularity-opac/base/opac_error.hpp> | ||
|
|
||
| namespace singularity { | ||
| namespace photons { | ||
|
|
||
| template <typename pc = PhysicalConstantsCGS> | ||
| class PowerLawSOpacity { | ||
| public: | ||
| PowerLawSOpacity() = default; | ||
| PowerLawSOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp, | ||
| const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., | ||
| const Real rho_ref = 1., const Real rho_off = 0., | ||
| const Real temp_ref = 1., const Real temp_off = 0., | ||
| const Real avg_particle_mass = 1.) | ||
| : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp), | ||
| rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), temp_off_(temp_off), | ||
| nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), apm_(avg_particle_mass) { | ||
| if (!(nu_ref_ > 0.)) { | ||
| OPAC_ERROR("PowerLawSOpacity: nu_ref must be positive"); | ||
| } | ||
| if (!(temp_ref_ > 0.)) { | ||
| OPAC_ERROR("PowerLawSOpacity: temp_ref must be positive"); | ||
| } | ||
| if (!(nu_ref_ > 0.)) { | ||
| OPAC_ERROR("PowerLawSOpacity: nu_ref must be positive"); | ||
| } | ||
| if (rho_off_ < 0.) { | ||
| OPAC_ERROR("PowerLawSOpacity: rho_off must be nonnegative"); | ||
| } | ||
| if (temp_off_ < 0.) { | ||
| OPAC_ERROR("PowerLawSOpacity: temp_off must be nonnegative"); | ||
| } | ||
| if (nu_off_ < 0.) { | ||
| OPAC_ERROR("PowerLawSOpacity: nu_off must be nonnegative"); | ||
| } | ||
| } | ||
|
|
||
| PowerLawSOpacity GetOnDevice() { return *this; } | ||
| PORTABLE_INLINE_FUNCTION | ||
| int nlambda() const noexcept { return 0; } | ||
| PORTABLE_INLINE_FUNCTION | ||
| void PrintParams() const noexcept { | ||
| printf("Power law scattering opacity. kappa0 = %g rho_exp = %g temp_exp = %g " | ||
| "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" | ||
| "temp_ref = %g temp_off = %g avg particle mass = %g\n", | ||
| kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, | ||
| nu_off_, rho_ref_, rho_off_, temp_ref_, temp_off_, | ||
| apm_); | ||
| } | ||
| inline void Finalize() noexcept {} | ||
|
|
||
| PORTABLE_INLINE_FUNCTION | ||
| Real TotalCrossSection(const Real rho, const Real temp, const Real nu, | ||
| Real *lambda = nullptr) const { | ||
| return OpacityScale_(rho, temp, nu); | ||
| } | ||
|
|
||
| PORTABLE_INLINE_FUNCTION | ||
| Real DifferentialCrossSection(const Real rho, const Real temp, const Real nu, | ||
| const Real mu, Real *lambda = nullptr) const { | ||
| // assumed isotropic, elastic | ||
| return OpacityScale_(rho, temp, nu) / (4. * M_PI); | ||
| } | ||
|
|
||
| PORTABLE_INLINE_FUNCTION | ||
| Real TotalScatteringCoefficient(const Real rho, const Real temp, | ||
| const Real nu, Real *lambda = nullptr) const { | ||
| return (rho / apm_) * OpacityScale_(rho, temp, nu); | ||
| } | ||
|
|
||
| private: | ||
| PORTABLE_INLINE_FUNCTION | ||
| Real OpacityScale_(const Real rho, const Real temp, const Real nu) const { | ||
| return OpacityPrefactor_(rho, temp) * std::pow((nu + nu_off_) / nu_ref_, nu_exp_); | ||
| } | ||
|
|
||
| PORTABLE_INLINE_FUNCTION | ||
| Real OpacityPrefactor_(const Real rho, const Real temp) const { | ||
| const Real rhom = (rho + rho_off_) / rho_ref_; | ||
| const Real tempm = (temp + temp_off_) / temp_ref_; | ||
| return kappa0_ * std::pow(rhom, rho_exp_) * std::pow(tempm, temp_exp_); | ||
| } | ||
|
|
||
| Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. | ||
| Real rho_exp_; // Power law index of density | ||
| Real temp_exp_; // Power law index of temperature | ||
| Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 | ||
| Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 | ||
| Real temp_ref_; // Temperature normalization for temp_exp. Units of K | ||
| Real temp_off_; // Temperature offset (same units as temp_ref). Units of K | ||
| Real nu_exp_; // Power law index of frequency | ||
| Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s | ||
| Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s | ||
| Real apm_; // Mean molecular weight. Units of g | ||
| }; | ||
|
|
||
| } // namespace photons | ||
| } // namespace singularity | ||
|
|
||
| #endif // SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this brings in units... can you double check that all the conversions are threaded through correctly to get the correctly unit-informed Planck and Boltzmann constant?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good catch - will do...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since this is the CGS side of the class, it seems
nuandtempshould be in CGS by the time they reach this function, so that leaves the question of whetherpc::handpc::kbare properly non-dimensionalizing the exponential argument (or equivalently, what the actual units ofnuis - e.g. energy or 1/time). If thenuhere is assumed consistent with thenutaken by theEmissivityPerNuPerOmegafunction, it seems this should be correct:EmissivityPerNuPerOmegacallsThermalDistributionOfNuwhich takesnudirectly and computes the nondimensional parameterx = pc::h * nu / (pc::kb * temp)as in the stimulated emission formula here. Let me know if this reasoning misses something.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Following on this, noticed that the parameters
rho_refandtemp_refmust be in CGS when the constructor is invoked. That is maybe a bit cumbersome - noted this in the comments next to these class members in each of the power law classes (also made a note thatnu_refhas to have units of 1/s, which I think is true for the same reason unlesskappa0absorbs the unit).