From 9837796bd228fbd2c34b81ec9571ddcab5d950b5 Mon Sep 17 00:00:00 2001 From: Thomas Ceulemans Date: Tue, 9 Jul 2024 12:09:25 +0200 Subject: [PATCH] Implemented other workaround for masering values; Now bounds the minimum line opacity instead. --- .../4_extra_options/NLTE_in_low_density.rst | 21 ++++++++++-- src/bindings/pybindings.cpp | 2 ++ .../lineProducingSpecies.tpp | 20 +++++++---- src/model/parameters/parameters.hpp | 16 +++++---- src/solver/solver.tpp | 34 ++++++++----------- .../analytic/all_zero_single_ray.py | 16 ++++----- .../analytic/density_distribution_1D.py | 16 ++++----- 7 files changed, 72 insertions(+), 53 deletions(-) diff --git a/docs/src/4_extra_options/NLTE_in_low_density.rst b/docs/src/4_extra_options/NLTE_in_low_density.rst index 8c963671..70b44b67 100644 --- a/docs/src/4_extra_options/NLTE_in_low_density.rst +++ b/docs/src/4_extra_options/NLTE_in_low_density.rst @@ -5,12 +5,27 @@ In low density regimes, the level populations might start to maser. This means that the resulting line opacity (and optical depth) can become close to zero or even negative. Unfortunately, our solvers cannot handle this situation too well, and might produce unphysical results. -In more recent versions of the code (starting from version 0.5.3), a workaround has been implemented to avoid this situation by resetting specific levels at specific places to LTE. +This is why we implement a workaround by setting the minimum allowed value for the line opacity (starting from version 0.7.0). +This value is set to 1e-13 [W sr−1 m−3] by default, but can be changed by the user. Higher values might overestimate the impact of the less dense model regions on the intensity, while lower value can lead to numerical artifacts in the NLTE intensities. + +.. code-block:: python + + parameters.min_line_opacity = 1e-13#default + +In versions 0.5.3 until 0.6.0, the code implemented another workaround to avoid this situation by resetting specific levels at specific places to LTE. This might work in some cases (with the caveat of some artifacts being produced in the resulting spectrum), but might as well produce new artifacts in when doing NLTE line radiative transfer in extremely low density regions. -To restore the previous behavior (without the new workaround), one can call +To restore that behavior (without the new workaround), one can call + +.. code-block:: python + + parameters.population_inversion_fraction = 1.01#previous default + parameters.min_line_opacity = -1 + +Finally, to restore the behavior of the code before version 0.5.3 (not recommended when using NLTE), one can call .. code-block:: python - parameters.population_inversion_fraction = -1 \ No newline at end of file + parameters.population_inversion_fraction = -1 + parameters.min_line_opacity = -1 \ No newline at end of file diff --git a/src/bindings/pybindings.cpp b/src/bindings/pybindings.cpp index 5befbbdc..b5e50d4c 100644 --- a/src/bindings/pybindings.cpp +++ b/src/bindings/pybindings.cpp @@ -288,6 +288,8 @@ PYBIND11_MODULE(core, module) { .def_readwrite("pop_prec", &Parameters::pop_prec, "Required precision for ALI.") .def_readwrite("min_opacity", &Parameters::min_opacity, "Minimum opacity that will be assumed in the solver.") + .def_readwrite("min_line_opacity", &Parameters::min_line_opacity, + "Minimum line opacity that will be assumed in the solver.") .def_readwrite("min_dtau", &Parameters::min_dtau, "Minimum optical depth increment that will be assumed in the solver.") .def_readwrite("population_inversion_fraction", &Parameters::population_inversion_fraction, diff --git a/src/model/lines/lineProducingSpecies/lineProducingSpecies.tpp b/src/model/lines/lineProducingSpecies/lineProducingSpecies.tpp index 3be04e6d..00b61cf3 100644 --- a/src/model/lines/lineProducingSpecies/lineProducingSpecies.tpp +++ b/src/model/lines/lineProducingSpecies/lineProducingSpecies.tpp @@ -21,9 +21,10 @@ inline Size LineProducingSpecies ::index(const Size p, const Size i) const { /// @return line emissivity for cell p and transition k ////////////////////////////////////////////////////////// inline Real LineProducingSpecies ::get_emissivity(const Size p, const Size k) const { - const Size i = index(p, linedata.irad[k]); + const Size i = index(p, linedata.irad[k]); + const Real freq = linedata.frequency[k]; - return HH_OVER_FOUR_PI * linedata.A[k] * population(i); + return freq * HH_OVER_FOUR_PI * linedata.A[k] * population(i); } /// Getter for the line opacity @@ -32,10 +33,13 @@ inline Real LineProducingSpecies ::get_emissivity(const Size p, const Size k) co /// @return line opacity for cell p and transition k /////////////////////////////////////////////////////// inline Real LineProducingSpecies ::get_opacity(const Size p, const Size k) const { - const Size i = index(p, linedata.irad[k]); - const Size j = index(p, linedata.jrad[k]); + const Size i = index(p, linedata.irad[k]); + const Size j = index(p, linedata.jrad[k]); + const Real freq = linedata.frequency[k]; - return HH_OVER_FOUR_PI * (population(j) * linedata.Ba[k] - population(i) * linedata.Bs[k]); + return std::max( + freq * HH_OVER_FOUR_PI * (population(j) * linedata.Ba[k] - population(i) * linedata.Bs[k]), + parameters->min_line_opacity); } /// set_LTE_level_populations @@ -435,8 +439,9 @@ inline void LineProducingSpecies::update_using_statistical_equilibrium( for (Size k = 0; k < linedata.nrad; k++) { for (Size m = 0; m < lambda.get_size(p, k); m++) { + const Real freq = linedata.frequency[k]; // line frequency const Size nr = lambda.get_nr(p, k, m); - const long double v_IJ = -lambda.get_Ls(p, k, m) * get_opacity(p, k); + const long double v_IJ = -lambda.get_Ls(p, k, m) * get_opacity(p, k) / freq; // Note: we define our transition matrix as the transpose of R in the // paper. @@ -628,8 +633,9 @@ inline void LineProducingSpecies::update_using_statistical_equilibrium_sparse( // Approximated Lambda operator for (Size k = 0; k < linedata.nrad; k++) { for (Size m = 0; m < lambda.get_size(p, k); m++) { + const Real freq = linedata.frequency[k]; // line frequency const Size nr = lambda.get_nr(p, k, m); - const long double v_IJ = -lambda.get_Ls(p, k, m) * get_opacity(p, k); + const long double v_IJ = -lambda.get_Ls(p, k, m) * get_opacity(p, k) / freq; if (nr != p) { throw std::runtime_error("ERROR: non-local Approximated Lambda operator."); diff --git a/src/model/parameters/parameters.hpp b/src/model/parameters/parameters.hpp index fa453fcb..2462964d 100644 --- a/src/model/parameters/parameters.hpp +++ b/src/model/parameters/parameters.hpp @@ -73,13 +73,15 @@ struct Parameters { Real max_width_fraction = 0.35; // max diff is +-2.5% // Real max_width_fraction = 0.3;//max diff is // +-2% - Real convergence_fraction = 0.995; - Real min_rel_pop_for_convergence = 1.0e-10; - Real pop_prec = 1.0e-6; - Real min_opacity = 1.0e-26; - Real min_dtau = 1.0e-15; - Real population_inversion_fraction = 1.01; // threshold factor for population inversion required - // for LTE to be used; set this higher than 1 + Real convergence_fraction = 0.995; + Real min_rel_pop_for_convergence = 1.0e-10; + Real pop_prec = 1.0e-6; + Real min_opacity = 1.0e-26; + long double min_line_opacity = 1.0e-13; + Real min_dtau = 1.0e-15; + Real population_inversion_fraction = + -1.0; // previously 1.01; // threshold factor for population inversion required + // for LTE to be used; set this higher than 1 bool store_intensities = false; // bool use_Ng_acceleration = true;//Not used, // so may safely be removed diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp index 85cbeb2c..5f96ed30 100644 --- a/src/solver/solver.tpp +++ b/src/solver/solver.tpp @@ -992,7 +992,7 @@ accel inline void Solver ::get_eta_and_chi(const Model& model, const Size // Set line emissivity and opacity for (Size l = 0; l < model.parameters->nlines(); l++) { const Real diff = freq - model.lines.line[l]; - const Real prof = freq * gaussian(model.lines.inverse_width(p, l), diff); + const Real prof = gaussian(model.lines.inverse_width(p, l), diff); eta += prof * model.lines.emissivity(p, l); chi += prof * model.lines.opacity(p, l); @@ -1038,7 +1038,7 @@ accel inline void Solver ::get_eta_and_chi(const Model& model, const const Real diff = freq - *freq_sort_l; // should be equal to the // previous line of code const Real inv_width = model.lines.inverse_width(p, l); - const Real prof = freq * gaussian(model.lines.inverse_width(p, l), diff); + const Real prof = gaussian(model.lines.inverse_width(p, l), diff); eta += prof * model.lines.emissivity(p, l); chi += prof * model.lines.opacity(p, l); } @@ -1058,7 +1058,7 @@ template <> accel inline void Solver ::get_eta_and_chi( const Model& model, const Size p, const Size l, const Real freq, Real& eta, Real& chi) const { const Real diff = freq - model.lines.line[l]; - const Real prof = freq * gaussian(model.lines.inverse_width(p, l), diff); + const Real prof = gaussian(model.lines.inverse_width(p, l), diff); eta = prof * model.lines.emissivity(p, l); chi = prof * model.lines.opacity(p, l) + model.parameters->min_opacity; @@ -1558,11 +1558,9 @@ inline void Solver ::compute_S_dtau_line_integrated(Model& model, Size Real& Snext) { dtau = compute_dtau_single_line(model, currpoint, nextpoint, lineidx, currfreq, nextfreq, dZ); Scurr = model.lines.emissivity(currpoint, lineidx) - / (model.lines.opacity(currpoint, lineidx) - + model.parameters->min_opacity); // current source - Snext = - model.lines.emissivity(nextpoint, lineidx) - / (model.lines.opacity(nextpoint, lineidx) + model.parameters->min_opacity); // next source + / model.lines.opacity(currpoint, lineidx); // current source + Snext = model.lines.emissivity(nextpoint, lineidx) + / model.lines.opacity(nextpoint, lineidx); // next source // note: due to interaction with dtau when computing all // sources individually, we do need to recompute Scurr and // Snext for all position increments @@ -1600,12 +1598,10 @@ inline void Solver ::compute_S_dtau_line_integrated(Model& model, Size cur for (Size l = 0; l < model.parameters->nlines(); l++) { Real line_dtau = compute_dtau_single_line(model, currpoint, nextpoint, l, currfreq, nextfreq, dZ); - Real line_Scurr = - model.lines.emissivity(currpoint, l) - / (model.lines.opacity(currpoint, l) + model.parameters->min_opacity); // current source + Real line_Scurr = model.lines.emissivity(currpoint, l) + / model.lines.opacity(currpoint, l); // current source Real line_Snext = - model.lines.emissivity(nextpoint, l) - / (model.lines.opacity(nextpoint, l) + model.parameters->min_opacity); // next source + model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l); // next source sum_dtau += line_dtau; sum_dtau_times_Scurr += line_dtau * line_Scurr; sum_dtau_times_Snext += line_dtau * line_Snext; @@ -1688,12 +1684,10 @@ inline void Solver ::compute_S_dtau_line_integrated(Model& model, Si Real line_dtau = compute_dtau_single_line(model, currpoint, nextpoint, l, currfreq, nextfreq, dZ); - Real line_Scurr = - model.lines.emissivity(currpoint, l) - / (model.lines.opacity(currpoint, l) + model.parameters->min_opacity); // current source + Real line_Scurr = model.lines.emissivity(currpoint, l) + / model.lines.opacity(currpoint, l); // current source Real line_Snext = - model.lines.emissivity(nextpoint, l) - / (model.lines.opacity(nextpoint, l) + model.parameters->min_opacity); // next source + model.lines.emissivity(nextpoint, l) / model.lines.opacity(nextpoint, l); // next source sum_dtau += line_dtau; sum_dtau_times_Scurr += line_dtau * line_Scurr; sum_dtau_times_Snext += line_dtau * line_Snext; @@ -1972,8 +1966,8 @@ inline Real Solver ::compute_dtau_single_line(Model& model, Size curridx, Size n // opacity is stored divided by the linefreq, so multiply // by it - const Real curr_line_opacity = linefreq * model.lines.opacity(curridx, lineidx); - const Real next_line_opacity = linefreq * model.lines.opacity(nextidx, lineidx); + const Real curr_line_opacity = model.lines.opacity(curridx, lineidx); + const Real next_line_opacity = model.lines.opacity(nextidx, lineidx); // if frequencies are equal, division by zero (due to the // optical depth formula) happens if we were not to use diff --git a/tests/benchmarks/analytic/all_zero_single_ray.py b/tests/benchmarks/analytic/all_zero_single_ray.py index fc56381b..0a0805d9 100644 --- a/tests/benchmarks/analytic/all_zero_single_ray.py +++ b/tests/benchmarks/analytic/all_zero_single_ray.py @@ -66,12 +66,12 @@ def velocity(i): model.thermodynamics.temperature.gas .set( temp * np.ones(npoints)) model.thermodynamics.turbulence.vturb2.set((turb/magritte.CC)**2 * np.ones(npoints)) - model = setup.set_Delaunay_neighbor_lists (model) - model = setup.set_Delaunay_boundary (model) - model = setup.set_boundary_condition_CMB (model) - model = setup.set_uniform_rays (model) - model = setup.set_linedata_from_LAMDA_file(model, lamdaFile) - model = setup.set_quadrature (model) + setup.set_Delaunay_neighbor_lists (model) + setup.set_Delaunay_boundary (model) + setup.set_boundary_condition_CMB (model) + setup.set_uniform_rays (model) + setup.set_linedata_from_LAMDA_file(model, lamdaFile) + setup.set_quadrature (model) model.write() @@ -176,8 +176,8 @@ def u_ (x): plt.show() #error bounds are chosen somewhat arbitrarily, based on previously obtained results; this should prevent serious regressions. - FIRSTORDER_AS_EXPECTED=(np.max(error_u_0s)<2.0e-6) - FEAUTRIER_AS_EXPECTED=(np.max(error_u_2f)<2.0e-6) + FIRSTORDER_AS_EXPECTED=(np.max(error_u_0s)<6.42e-6) + FEAUTRIER_AS_EXPECTED=(np.max(error_u_2f)<6.41e-6) if not FIRSTORDER_AS_EXPECTED: print("First order solver max error too large: ", np.max(error_u_0s)) diff --git a/tests/benchmarks/analytic/density_distribution_1D.py b/tests/benchmarks/analytic/density_distribution_1D.py index fff3b982..4a9887ab 100644 --- a/tests/benchmarks/analytic/density_distribution_1D.py +++ b/tests/benchmarks/analytic/density_distribution_1D.py @@ -71,12 +71,12 @@ def nTT (r): model.thermodynamics.temperature.gas .set( temp * np.ones(npoints)) model.thermodynamics.turbulence.vturb2.set((turb/magritte.CC)**2 * np.ones(npoints)) - model = setup.set_Delaunay_neighbor_lists (model) - model = setup.set_Delaunay_boundary (model) - model = setup.set_boundary_condition_CMB (model) - model = setup.set_rays_spherical_symmetry (model) - model = setup.set_linedata_from_LAMDA_file(model, lamdaFile) - model = setup.set_quadrature (model) + setup.set_Delaunay_neighbor_lists (model) + setup.set_Delaunay_boundary (model) + setup.set_boundary_condition_CMB (model) + setup.set_rays_spherical_symmetry (model) + setup.set_linedata_from_LAMDA_file(model, lamdaFile) + setup.set_quadrature (model) model.write() @@ -238,8 +238,8 @@ def plot (r, p): #setting 'b' not yet used for testing if a_or_b == 'a': #error bounds are chosen somewhat arbitrarily, based on previously obtained results; this should prevent serious regressions. - FEAUTRIER_AS_EXPECTED=(np.mean(error_u_2f)<3e-4) - FIRSTORDER_AS_EXPECTED=(np.mean(error_u_0s)<2.7e-4) + FEAUTRIER_AS_EXPECTED=(np.mean(error_u_2f)<4.4e-4) + FIRSTORDER_AS_EXPECTED=(np.mean(error_u_0s)<4.1e-4) if not FIRSTORDER_AS_EXPECTED: print("First order solver mean error too large: ", np.mean(error_u_0s))