Accuracy target for higher precision and tweaks for check_accuracy stability

camb/_config.py

@@ -12,6 +12,9 @@ class _config:
     # print feedback if > 0 (note in Jupyter notebook this will appear in the terminal, not the notebook)
     FeedbackLevel = import_property(c_int, "config", "FeedbackLevel")
 
+    # enable targeted accuracy improvements if > 0, AccuracyTarget being SO-like.
+    AccuracyTarget = import_property(c_int, "config", "AccuracyTarget")
+
     # print additional timing and progress (when FeedbackLevel>0)
     DebugMsgs = import_property(c_bool, "config", "DebugMsgs")
 

camb/tests/camb_test.py

@@ -557,25 +557,21 @@ def testPowers(self):
         pk_interp2 = PKnonlin2.P(z, kh)
         self.assertTrue(np.sum((pk_interp / pk_interp2 - 1) ** 2) < 0.005)
 
-        pars.NonLinearModel.set_params(halofit_version="mead")
-        _, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
-        self.assertAlmostEqual(pk[0][160], 814.9, delta=0.5)
+        # The nonlinear matter grid can move when CMB source sampling changes; compare at a fixed physical k/h.
+        def pk_at_fixed_kh(halofit_version, kh_value=0.44223076105117803):
+            pars.NonLinearModel.set_params(halofit_version=halofit_version)
+            kh, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
+            return np.exp(np.interp(np.log(kh_value), np.log(kh), np.log(pk[0])))
 
-        pars.NonLinearModel.set_params(halofit_version="mead2016")
-        _, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
-        self.assertAlmostEqual(pk[0][160], 814.9, delta=0.5)
+        self.assertAlmostEqual(pk_at_fixed_kh("mead"), 814.9, delta=0.5)
 
-        pars.NonLinearModel.set_params(halofit_version="mead2015")
-        _, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
-        self.assertAlmostEqual(pk[0][160], 791.3, delta=0.5)
+        self.assertAlmostEqual(pk_at_fixed_kh("mead2016"), 814.9, delta=0.5)
 
-        pars.NonLinearModel.set_params(halofit_version="mead2020")
-        _, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
-        self.assertAlmostEqual(pk[0][160], 815.8, delta=0.5)
+        self.assertAlmostEqual(pk_at_fixed_kh("mead2015"), 791.3, delta=0.5)
 
-        pars.NonLinearModel.set_params(halofit_version="mead2020_feedback")
-        _, _, pk = results.get_nonlinear_matter_power_spectrum(params=pars)
-        self.assertAlmostEqual(pk[0][160], 799.0, delta=0.5)
+        self.assertAlmostEqual(pk_at_fixed_kh("mead2020"), 815.8, delta=0.5)
+
+        self.assertAlmostEqual(pk_at_fixed_kh("mead2020_feedback"), 799.0, delta=0.5)
 
         lmax = 4000
         pars.set_for_lmax(lmax)

docs/changelog/HighL_accuracy_stability.md

@@ -0,0 +1,96 @@
+# Planck 2018 Accuracy Stability
+
+## Summary
+
+This investigation targeted
+
+```bash
+camb check_accuracy inifiles/planck_2018.ini --set-for-lmax 4000 --lens-potential-accuracy 4
+```
+
+with separate treatment of CMB power spectra and matter power. The final changes
+avoid broad accuracy-boost floors and instead adjust the specific source-grid,
+lensing-correlation, and late photon-hierarchy ranges that control the observed
+instability.
+
+## Main Conclusions
+
+- The low-ell EE instability is controlled by low-`k` CMB source sampling, not by
+  high-`L` lensing physics. A direct cap on the low-`k` logarithmic source step is
+  preferable to increasing a global integration or source boost.
+- The high-`L` lensed CMB instability at `lmax=4000` is controlled by
+  correlation-function angular sampling. A separate `ThetaSampleBoost` floor is
+  cleaner than changing `LensAccuracyBoost`, which also affects unrelated
+  lensing work.
+- The matter-power failure is not fixed by increasing output
+  `transfer_k_per_logint`; sweeps up to very dense output sampling left the same
+  interpolation-stable mismatch. The sensitive calculation is the shared CMB
+  source `q` grid used by high-precision transfer output, plus the late photon
+  hierarchy truncation time.
+- The high-precision transfer source-grid spacing should follow an explicit
+  dense `transfer_k_per_logint` request. For automatic sampling
+  (`transfer_k_per_logint = 0`), 20 points per log interval is sufficient for the
+  low-`q` source grid in this case.
+- Extra refinement of the smooth high-`q` source tail via `dksmooth` was tested
+  and found unnecessary after the low/intermediate source-grid fixes.
+- At `lmax=6000` and `8000`, larger theta sampling did not resolve the remaining
+  CMB differences, even with higher lensing accuracy. Those failures appear to be
+  a separate high-`lmax` issue and are not addressed by this change set.
+
+## Code Changes
+
+- `fortran/cmbmain.f90`
+  - Caps the low-`k` CMB source logarithmic step for accurate reionization
+    polarization.
+  - Uses a combined high-precision transfer/CMB-source condition for the source
+    grid reused by transfer output.
+  - Sets the low-`q` source log density to
+    `max(20, transfer_k_per_logint)` for high-precision transfer output.
+  - Refines the two linear source intervals, `dkn1` and `dkn2`, by a factor of
+    `1.5` for this shared high-precision transfer source grid.
+  - Caps `qmax_log` at the horizon-entry switch to avoid an unnecessary
+    intermediate interval when the logarithmic range already reaches that scale.
+  - Increases the low-`k` scalar integration grid minimum from 10 to 11 log
+    samples for accurate reionization polarization.
+
+- `fortran/equations.f90`
+  - Isolates the high-precision transfer switch change to late photon multipole
+    truncation, `tau_switch_no_phot_multpoles`.
+  - Leaves late neutrino multipole truncation and massive-neutrino switch timing
+    on the original `TimeSwitchBoost` behavior.
+
+- `fortran/lensing.f90`
+  - Adds `ThetaSampleBoost` as a separate angular-sampling factor for lensed CMB
+    correlation functions.
+  - Applies a floor of `1.8` only for `CP%Max_l > 3500`, avoiding a broad
+    `LensAccuracyBoost` increase.
+
+- `camb/tests/camb_test.py`
+  - Updates the nonlinear matter-power regression check to compare at a fixed
+    physical `k/h` by interpolation, rather than checking the grid-fragile
+    `pk[0][160]` index.
+
+## Tests
+
+- Original target command passes with automatic transfer sampling:
+  - matter power max error: `0.0008865`
+  - standard run timing: about `8.33s` CPU, `1.08s` wall in the final check run
+- Explicit dense transfer sampling also passes:
+  - command uses `transfer_k_per_logint = 50`
+  - matter power max error: `0.0007711`
+  - standard run timing: about `11.59s` CPU, `1.50s` wall
+- `python -m unittest camb.tests.camb_test` passes: 20 tests.
+- `git diff --check` passes for the touched source and test files.
+
+## Negative Tests
+
+- Removing the `dkn2` refinement fails with matter-power max error `0.001343`
+  near `k/h = 0.104`.
+- Removing the `dkn1` refinement fails with matter-power max error `0.001279`
+  near `k/h = 0.076`.
+- Refining both `dkn1` and `dkn2` by only `1.25` passes but with too little
+  margin: matter-power max error `0.0009885`.
+- Refining `dksmooth` by `4` passes but is not needed and is slower; it was
+  removed from the final patch.
+- Neutrino-only late radiation truncation was not sufficient for the matter
+  power stability target; photon-only late truncation was sufficient.

fortran/camb.f90

@@ -887,6 +887,7 @@ subroutine CAMB_CommandLineRun(InputFile)
     highL_unlensed_cl_template = Ini%Read_String_Default( &
         'highL_unlensed_cl_template', highL_unlensed_cl_template)
     call Ini%Read('number_of_threads', ThreadNum)
+    call Ini%Read('AccuracyTarget', AccuracyTarget)
     call Ini%Read('DebugParam', DebugParam)
     call Ini%Read('feedback_level', FeedbackLevel)
     if (Ini%HasKey('DebugMsgs')) call Ini%Read('DebugMsgs', DebugMsgs)

fortran/cmbmain.f90

@@ -850,12 +850,15 @@ subroutine SetkValuesForSources
     implicit none
     real(dl) dlnk0, dkn1, dkn2, q_switch, q_cmb, dksmooth
     real(dl) qmax_log
-    real(dl) SourceAccuracyBoost
+    real(dl) SourceAccuracyBoost, transferLogDensity
+    logical highPrecisionTransferSources
     !     set k values for which the sources for the anisotropy and
     !     polarization will be calculated. For low values of k we
     !     use a logarithmic spacing. closed case dealt with by SetClosedkValues
 
     SourceAccuracyBoost = CP%Accuracy%AccuracyBoost * CP%Accuracy%SourcekAccuracyBoost
+    highPrecisionTransferSources = AccuracyTarget > 0 .and. CP%Want_CMB .and. CP%WantScalars .and. CP%WantTransfer &
+        .and. CP%Transfer%high_precision
     if (CP%WantScalars .and. CP%Reion%Reionization .and. CP%Accuracy%AccuratePolarization) then
         dlnk0=2._dl/10/SourceAccuracyBoost
         !Need this to get accurate low l polarization
@@ -865,17 +868,31 @@ subroutine SetkValuesForSources
     end if
 
     if (CP%Accuracy%AccurateReionization) dlnk0 = dlnk0/2
+    !Low-l reionization polarization benefits from a modest cap on the low-k logarithmic source grid.
+    if (AccuracyTarget > 0 .and. CP%Want_CMB .and. CP%WantScalars .and. CP%Reion%Reionization .and. &
+        CP%Accuracy%AccuratePolarization .and. CP%Accuracy%AccurateReionization) dlnk0 = min(dlnk0, 0.075_dl)
+    if (highPrecisionTransferSources) then
+        !High-precision transfer output reuses this CMB source q grid up to q_cmb.
+        !Use 20/log interval for automatic transfer sampling, or the requested denser transfer grid.
+        transferLogDensity = max(20._dl, real(CP%Transfer%k_per_logint, dl))
+        dlnk0 = min(dlnk0, 1._dl/transferLogDensity)
+    end if
 
     dkn1=0.6_dl/State%taurst/SourceAccuracyBoost
     dkn2=0.9_dl/State%taurst/SourceAccuracyBoost/1.2
+    if (highPrecisionTransferSources) then
+        dkn1 = dkn1/1.5_dl
+        dkn2 = dkn2/1.5_dl
+    end if
     if (CP%WantTensors .or. CP%WantVectors) then
         dkn1=dkn1  *0.8_dl
         dlnk0=dlnk0/2 !*0.3_dl
         dkn2=dkn2*0.85_dl
     end if
 
-    qmax_log = dkn1/dlnk0
     q_switch = 2*6.3/State%taurst
+    !If the low-k logarithmic grid reaches the horizon-entry switch, skip the intermediate linear interval.
+    qmax_log = min(dkn1/dlnk0, q_switch)
     !Want linear spacing for wavenumbers which come inside horizon
     !Could use sound horizon, but for tensors that is not relevant
 
@@ -889,7 +906,7 @@ subroutine SetkValuesForSources
     associate(Evolve_q => ThisSources%Evolve_q)
         call Evolve_q%Init()
         call Evolve_q%Add_delta(qmin, qmax_log, dlnk0, IsLog = .true.)
-        if (qmax > qmax_log) &
+        if (qmax > qmax_log .and. q_switch > qmax_log) &
             call Evolve_q%Add_delta(qmax_log, min(qmax,q_switch), dkn1)
         if (qmax > q_switch) then
             call Evolve_q%Add_delta(q_switch, min(q_cmb,qmax), dkn2)
@@ -1271,6 +1288,8 @@ subroutine SetkValuesForInt(ThisCT)
         !then no-lognum*dk0 linearly spaced at dk0 up to no*dk0
         !then at dk up to max_k_dk, then dk2 up to qmax_int
         lognum=nint(10*IntSampleBoost)
+        if (AccuracyTarget > 0 .and. CP%Want_CMB .and. CP%Accuracy%AccuratePolarization &
+            .and. CP%Accuracy%AccurateReionization) lognum = max(lognum, 11)
         dlnk1=1._dl/lognum
         no=nint(600*IntSampleBoost)
         dk0=1.8_dl/State%curvature_radius/State%chi0/IntSampleBoost

fortran/config.f90

@@ -7,6 +7,8 @@ module config
 
     integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
 
+    integer :: AccuracyTarget = 1 !if >0 enable targeted accuracy improvements
+
     logical :: output_file_headers = .true.
 
     logical :: DebugMsgs =.false. !Set to true to view progress and timing

fortran/equations.f90

@@ -259,12 +259,13 @@ recursive subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,rk_setti
     real(dl) cs2, opacity, dopacity
     real(dl) tau_switch_ktau, tau_switch_nu_massless, tau_switch_nu_massive, next_switch
     real(dl) tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles,tau_switch_nu_nonrel
-    real(dl) noSwitch, smallTime
+    real(dl) noSwitch, smallTime, switchBoost, latePhotSwitchBoost
     !Sources
     real(dl) tau_switch_saha, Delta_TM, xe,a,tau_switch_evolve_TM
 
     noSwitch= State%tau0+1
     smallTime =  min(tau, 1/EV%k_buf)/100
+    switchBoost = CP%Accuracy%AccuracyBoost*CP%Accuracy%TimeSwitchBoost
 
     tau_switch_ktau = noSwitch
     tau_switch_no_nu_multpoles= noSwitch
@@ -302,13 +303,20 @@ recursive subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,rk_setti
     if (CP%DoLateRadTruncation) then
         if (.not. EV%no_nu_multpoles) & !!.and. .not. EV%has_nu_relativistic .and. tau_switch_nu_massless ==noSwitch)  &
             tau_switch_no_nu_multpoles= &
-            max(15/EV%k_buf*CP%Accuracy%AccuracyBoost*CP%Accuracy%TimeSwitchBoost,&
+            max(15/EV%k_buf*switchBoost,&
             min(State%taurend,EV%ThermoData%matter_verydom_tau))
 
-        if (.not. EV%no_phot_multpoles .and. (.not.CP%WantCls .or. &
-            EV%k_buf>0.03*CP%Accuracy%AccuracyBoost*CP%Accuracy%TimeSwitchBoost)) &
-            tau_switch_no_phot_multpoles =max(15/EV%k_buf,State%taurend)*CP%Accuracy%AccuracyBoost &
-            *CP%Accuracy%TimeSwitchBoost
+        if (.not. EV%no_phot_multpoles) then
+
+            !High-precision transfer matter power is sensitive to late photon hierarchy truncation.
+            if (AccuracyTarget > 0 .and. CP%WantTransfer .and. CP%Transfer%high_precision) then
+                latePhotSwitchBoost = max(switchBoost, 2._dl)
+            else
+                latePhotSwitchBoost = switchBoost
+            end if
+            if (.not.CP%WantCls .or. EV%k_buf>0.03*latePhotSwitchBoost) &
+                tau_switch_no_phot_multpoles =max(15/EV%k_buf,State%taurend)*latePhotSwitchBoost
+        end if
     end if
 
     next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &

fortran/lensing.f90

@@ -168,7 +168,7 @@ subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
     logical :: short_integral_range
     real(dl) range_fac
     logical, parameter :: approx = .false.
-    real(dl) theta_cut(lmax), LensAccuracyBoost
+    real(dl) theta_cut(lmax), LensAccuracyBoost, ThetaSampleBoost
     Type(TTimer) :: Timer
 
     !$ integer  OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS
@@ -178,6 +178,12 @@ subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
     associate(lSamp => State%CLData%CTransScal%ls, CP=>State%CP)
 
         LensAccuracyBoost = CP%Accuracy%AccuracyBoost*CP%Accuracy%LensingBoost
+        ThetaSampleBoost = LensAccuracyBoost
+        !High-l lensed spectra need denser correlation-function angular sampling.
+        if (CP%Max_l > 3500) then
+            ThetaSampleBoost = ThetaSampleBoost*1.3_dl
+            if (AccuracyTarget > 0) ThetaSampleBoost = max(ThetaSampleBoost, 1.8_dl)
+        end if
         max_lensed_ix = lSamp%nl-1
         do while(lSamp%l(max_lensed_ix) > CP%Max_l - lensed_convolution_margin)
             max_lensed_ix = max_lensed_ix -1
@@ -188,10 +194,9 @@ subroutine CorrFuncFullSkyImpl(State,CL,CLout,CPP,lmin,lmax)
         if (allocated(CLout%Cl_lensed)) deallocate(CLout%Cl_lensed)
         allocate(CLout%Cl_lensed(lmin:CLout%lmax_lensed,1:4), source = 0._dl)
 
-        npoints = CP%Max_l  * 2 *LensAccuracyBoost
+        npoints = CP%Max_l  * 2 * ThetaSampleBoost
         short_integral_range = .not. CP%Accuracy%AccurateBB
         dtheta = const_pi / npoints
-        if (CP%Max_l > 3500) dtheta=dtheta/1.3
         apodize_point_width = nint(0.003 / dtheta)
         npoints = int(const_pi/dtheta)
         if (short_integral_range) then