From 7eeb34c7dc8f3ead08f0f4cb58da8a027ad60f51 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Fri, 5 Apr 2024 14:04:04 +0200 Subject: [PATCH] src/gray_core.f90: refine realtime termination check MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In some very rare cases the residual power check can stop the integration too early. For example, this extremely unlikely chain of events was observed in a few cases while performing a parameter sweep (~3000 simulations): 1. n⊥ converges on wrong branch ⇒ α jumps ⇒ dP/ds jumps (common) 2. dP/ds jumps close to the peak (unlikely) 3. the peak of dP/ds is close to 0.6P₀ (unlikely) 4. dP/ds jumps to dP/ds_max + 10⁻¹² ⇒ wrong peak found (very unlikely) 5. peak is at exactly the current point ⇒ three-point interpolation fails (unlikely) --- src/gray_core.f90 | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index ff3cba1..3ec99b0 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -529,24 +529,36 @@ contains end do end if - if (params%raytracing%realtime .and. pow <= 0.4_wp_*p0ray(1)) then - ! Compute the approximate position of the absorption peak + if (params%raytracing%realtime) then + ! Check whether we are past the absorption peak block - use utils, only : vmax, parabola_vertex - real(wp_) :: power, peak(2) - integer :: index + logical :: past_peak + integer :: tail + ! We assume so if dP/ds has been decreasing in the last 8 steps + tail = max(2, i-8) + past_peak = i > 8 .and. all(dpds(tail:i) - dpds(tail-1:i-1) < 0) - ! Find maximum in dP/ds - call vmax(dpds, i, power, index) + if (past_peak .or. pow <= 0.1_wp_*p0ray(1)) then + ! Compute the approximate position of the absorption peak + block + use utils, only : vmax, parabola_vertex + real(wp_) :: power, peak(2) + integer :: index - ! Interpolate the peak with a parabola - peak = parabola_vertex(psjki(1, index - 1:index + 1), & - dpds(index - 1:index + 1)) + ! Find maximum in dP/ds + call vmax(dpds, i, power, index) + + ! Interpolate the peak with a parabola + peak = parabola_vertex(psjki(1, index - 1:index + 1), & + dpds(index - 1:index + 1)) + + results%rho_peak = sqrt(peak(1)) + end block + ! Stop propagation loop + exit + end if - results%rho_peak = sqrt(peak(1)) end block - ! Stop propagation loop - exit end if ! print ray positions for j=nrayr in local reference system