Regular low-speed collision substep in stost (#615)

## Cause
In `collis_alp::stost` the speed p is advanced by Euler-Maruyama with
the substep bounded by `fraction*p/|fpeff|` and
`fraction**2*p**2/(6*dpp)`. Near p = 0 the speed process is a 3-D
Bessel-type process (`fpeff ~ 2*dpp/p`). In a cold plasma and from a
subthermal start, about one call in 1e6 to 1e7 lets the discrete walk
descend geometrically to p ~ 1e-10 within one call. The substep then
falls below ulp of the remaining time and `stost` returns `ierr = 4`, so
the caller rejects the step. Below p ~ 1e-8 the `plim` floor in
`coleff_local` also weakens the repulsive drift.

## Fix
For `iswmode` 1 and 2, below `p_switch = 1/max(velrat)` over the field
species with `efcolf > 0` (thermal speed of the slowest species over
v0), the substep is taken in Cartesian velocity space. This is the same
generator, but regular at p = 0:

`v_new = v + A*vhat*h + sqrt(2*Dpar*h)*xi1*vhat +
sqrt(2*Dperp*h)*(xi2*e2 + xi3*e3)`

with `Dpar = dpp`, `Dperp = dhh*p**2` and `A = fpeff - 2*Dperp/p`,
evaluated without the 1/p cancellation in the new helper
`coleff_cartesian` (next to `coleff_local`/`onseff`). Ito check: `dp =
(A + 2*Dperp/p) dt + sqrt(2*Dpar) dW = fpeff dt + ...` and the pitch
generator is `Dperp/p**2 d/dlambda (1-lambda**2) d/dlambda`. Above
`p_switch`, and whenever `velrat` is zero (synthetic tests), the
existing scheme is unchanged.

## Tests
New `test_collision_low_speed` (labels `unit;collisions;physics`), cold
H plasma at 20 eV with 2.5% C6+, n_e = 2e19 m^-3:
- Start at p = 1e-9 and 1e-12 must return `ierr = 0` (the old code
returns `ierr = 4` immediately, checked once by hand).
- First moments of one substep on both sides of `p_switch` match
`fpeff*h` and `-2*dhh*lambda*h` within 5 sigma.
- 4096 markers from p = 1e-6 thermalize to the Maxwellian `<p**2>`,
`<lambda> = 0`, `<lambda**2> = 1/3` within 5 sigma, and subthermal
markers never error over 2000 seeds.

Verified on macOS arm64 (gfortran 16, Release): `ctest -L collisions`
passes (12/12, including `test_collision_substeps` and the new test).
`golden_record_collision` could not run there, because
`golden_record.sh` uses bash 4 `exec {fd}` redirection.

Chris&AI

Author: krystophny

src/collis_alphas.f90

@@ -118,6 +118,38 @@ subroutine onseff(v, dp, dh, dpd)
         end if
     end subroutine onseff
 
+    !> Cartesian-velocity coefficients of the same collision operator, regular
+    !> at p = 0: parallel diffusion dpar, perpendicular diffusion dperp and the
+    !> drift a along the velocity, with fpeff = a + 2*dperp/p.
+    subroutine coleff_cartesian(p, efcolf_loc, velrat_loc, enrat_loc, &
+                                dpar, dperp, a)
+        real(wp), intent(in) :: p
+        real(wp), intent(in) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts)
+        real(wp), intent(out) :: dpar, dperp, a
+
+        real(wp), parameter :: cons = 0.75225278d0
+        integer :: i
+        real(wp) :: x, dp, dh, dpd, g
+
+        dpar = 0.0_wp
+        dperp = 0.0_wp
+        a = 0.0_wp
+        do i = 1, nsorts
+            x = p*velrat_loc(i)
+            call onseff(x, dp, dh, dpd)
+            if (x < 0.01d0) then
+                g = -2.0d0*cons*x
+            else if (x > 6.0d0) then
+                g = -2.0d0/x**2
+            else
+                g = (dpd - 2.0d0*dh)/x
+            end if
+            dpar = dpar + dp*efcolf_loc(i)
+            dperp = dperp + dh*efcolf_loc(i)
+            a = a + efcolf_loc(i)*(velrat_loc(i)*g - 2.0d0*dp*p*enrat_loc(i))
+        end do
+    end subroutine coleff_cartesian
+
     subroutine loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, &
                             tempe, ealpha, v0, dchichi, slowrate, &
                             dchichi_norm, slowrate_norm, test_mass, test_charge)
@@ -288,10 +320,11 @@ subroutine stost(z, dtauc, iswmode, ierr, max_substep_fraction, s_profile)
         integer, parameter :: max_substeps = 100000
         real(wp) :: initial(5), fraction, remaining, step, p, dpp, dhh, fpeff
         real(wp) :: alam, angle, azimuth
+        real(wp) :: p_switch, dpar, dperp, adrift, vhat(2), vnew(3), xi(3)
         real(wp) :: efcolf_loc(nsorts), velrat_loc(nsorts), enrat_loc(nsorts)
         real :: ur
-        integer :: substep
-        logical :: pitch, energy, diffusion
+        integer :: substep, i
+        logical :: pitch, energy, diffusion, cartesian_ok
 
         initial = z
         ierr = 0
@@ -324,9 +357,56 @@ subroutine stost(z, dtauc, iswmode, ierr, max_substep_fraction, s_profile)
         energy = iswmode /= 4
         diffusion = iswmode < 3
         remaining = dtauc
+        ! Below the thermal speed of the slowest field species the radial
+        ! process is Bessel-like and the Euler walk in p can collapse towards
+        ! p = 0. There the substep is taken in Cartesian velocity space.
+        cartesian_ok = iswmode <= 2 .and. any(efcolf_loc > 0.0_wp .and. &
+                                              velrat_loc > 0.0_wp)
+        p_switch = 0.0_wp
+        if (cartesian_ok) p_switch = 1.0_wp/maxval(velrat_loc, &
+                                                   mask=efcolf_loc > 0.0_wp)
 
         do substep = 1, max_substeps
             p = z(4)
+            if (cartesian_ok .and. p < p_switch) then
+                call coleff_cartesian(p, efcolf_loc, velrat_loc, enrat_loc, &
+                                      dpar, dperp, adrift)
+                if (.not. all(is_finite([dpar, dperp, adrift])) .or. &
+                    dpar < 0.0_wp .or. dperp < 0.0_wp) then
+                    ierr = 1
+                    exit
+                end if
+                step = remaining
+                if (max(dpar, dperp) > 0.0_wp) step = min(step, &
+                    fraction**2*p_switch**2/(6.0_wp*max(dpar, dperp)))
+                if (abs(adrift) > 0.0_wp) &
+                    step = min(step, fraction*p_switch/abs(adrift))
+                if (step <= 0.0_wp .or. remaining - step == remaining) then
+                    ierr = 4 ! Cannot represent further time progress.
+                    exit
+                end if
+                alam = z(5)
+                vhat = [alam, sqrt(1.0_wp - alam**2)]
+                do i = 1, 3
+                    call getran(0, ur)
+                    xi(i) = real(ur, wp)
+                end do
+                vnew(1:2) = p*vhat + (adrift*step + &
+                            sqrt(2.0_wp*dpar*step)*xi(1))*vhat + &
+                            sqrt(2.0_wp*dperp*step)*xi(2)*[-vhat(2), vhat(1)]
+                vnew(3) = sqrt(2.0_wp*dperp*step)*xi(3)
+                z(4) = sqrt(sum(vnew**2))
+                if (pitch .and. is_finite(z(4)) .and. z(4) > 0.0_wp) &
+                    z(5) = max(-1.0_wp, min(1.0_wp, vnew(1)/z(4)))
+                if (.not. all(is_finite(z)) .or. z(4) <= 0.0_wp .or. &
+                    abs(z(5)) > 1.0_wp) then
+                    ierr = 1
+                    exit
+                end if
+                remaining = remaining - step
+                if (remaining == 0.0_wp) return
+                cycle
+            end if
             call coleff_local(p, efcolf_loc, velrat_loc, enrat_loc, dpp, dhh, fpeff)
             if (.not. all(is_finite([dpp, dhh, fpeff])) .or. &
                 dpp < 0.0_wp .or. dhh < 0.0_wp) then

test/tests/CMakeLists.txt

@@ -535,6 +535,12 @@ add_test(NAME test_collision_substeps COMMAND test_collision_substeps.x)
 set_tests_properties(test_collision_substeps PROPERTIES
     LABELS "unit;collisions;physics" TIMEOUT 180)
 
+add_executable(test_collision_low_speed.x test_collision_low_speed.f90)
+target_link_libraries(test_collision_low_speed.x simple)
+add_test(NAME test_collision_low_speed COMMAND test_collision_low_speed.x)
+set_tests_properties(test_collision_low_speed PROPERTIES
+    LABELS "unit;collisions;physics" TIMEOUT 180)
+
 add_executable (test_profiles.x test_profiles.f90)
 target_link_libraries(test_profiles.x simple)
 add_test(NAME test_profiles COMMAND test_profiles.x)

test/tests/test_collision_low_speed.f90

@@ -0,0 +1,143 @@
+!> Low-speed behaviour of collis_alp::stost in a cold plasma. Below the thermal
+!> speed of the slowest field species the radial speed process is Bessel-like,
+!> and the Euler walk in p can collapse to p ~ 1e-10, where it cannot represent
+!> further time progress (ierr = 4). The Cartesian substep is regular there.
+program test_collision_low_speed
+    use, intrinsic :: iso_fortran_env, only: dp => real64
+    use collis_alp, only: stost, coleff, efcolf, velrat, enrat, efcolf_grid, &
+                          velrat_grid, enrat_grid, N_S_GRID, &
+                          compute_collision_coeffs
+    implicit none
+
+    real(dp), parameter :: ealpha = 55.0e3_dp, temp = 20.0_dp
+    real(dp) :: p_switch, rate, v0
+    integer :: nseed
+    integer, allocatable :: seed(:)
+
+    call random_seed(size=nseed)
+    allocate (seed(nseed))
+    call setup_cold_plasma()
+    p_switch = 1.0_dp/maxval(velrat)
+    rate = sum(efcolf*enrat)*2.0_dp*0.75225278_dp
+    print '(a,es10.3,a,es10.3)', 'p_switch = ', p_switch, ', rate = ', rate
+
+    call check_tiny_speed_start()
+    call check_generator_moments()
+    call check_maxwellian()
+    print *, 'PASS: low-speed collision substep'
+contains
+
+    subroutine setup_cold_plasma()
+        ! H test ion in a 20 eV H plasma with 2.5% C6+, n_e = 2e19 m^-3.
+        real(dp) :: ni1, ni2
+        integer :: node
+        ni1 = 2.0e13_dp/(1.0_dp + 6.0_dp*0.025_dp)
+        ni2 = 0.025_dp*ni1
+        call compute_collision_coeffs(1.0_dp, 12.0_dp, 1.0_dp, 6.0_dp, ni1, ni2, &
+                                      temp, temp, temp, ealpha, v0, &
+                                      efcolf, velrat, enrat)
+        do node = 1, N_S_GRID
+            efcolf_grid(:, node) = efcolf
+            velrat_grid(:, node) = velrat
+            enrat_grid(:, node) = enrat
+        end do
+    end subroutine setup_cold_plasma
+
+    subroutine check_tiny_speed_start()
+        ! Before the Cartesian substep, a start at p = 1e-9 or 1e-12 ended
+        ! immediately with ierr = 4 (substep ~ p**2 below ulp of the interval).
+        real(dp), parameter :: starts(2) = [1.0e-9_dp, 1.0e-12_dp]
+        integer :: k, mode, ierr
+        real(dp) :: z(5)
+        do mode = 1, 2
+            do k = 1, 2
+                z = [0.5_dp, 0.2_dp, 0.3_dp, starts(k), 0.4_dp]
+                call stost(z, 6.76_dp, mode, ierr)
+                if (ierr /= 0) error stop 'Tiny-speed start returned an error'
+                if (.not. all(z == z) .or. z(4) <= 0.0_dp .or. &
+                    abs(z(5)) > 1.0_dp) error stop 'Tiny-speed start left bad state'
+            end do
+        end do
+    end subroutine check_tiny_speed_start
+
+    subroutine check_generator_moments()
+        ! One substep from either side of p_switch: first moments of p and
+        ! lambda must match fpeff*h and -2*dhh*lambda*h.
+        integer, parameter :: n = 400000
+        real(dp), parameter :: lam = 0.5_dp
+        real(dp) :: fac(3) = [0.5_dp, 0.98_dp, 1.02_dp]
+        real(dp) :: p, h, dpp, dhh, fpeff, z(5), dp1, dl1, m1, m2, s1, s2
+        integer :: k, i, ierr
+        seed = 77
+        call random_seed(put=seed)
+        do k = 1, size(fac)
+            p = fac(k)*p_switch
+            if (k == 3) p = 1.02_dp*p_switch
+            call coleff(p, dpp, dhh, fpeff)
+            h = 0.5_dp*min(0.01_dp/dhh, 0.1_dp*p/abs(fpeff), &
+                           0.01_dp*p**2/(6.0_dp*dpp))
+            m1 = 0.0_dp
+            m2 = 0.0_dp
+            s1 = 0.0_dp
+            s2 = 0.0_dp
+            do i = 1, n
+                z = [0.5_dp, 0.2_dp, 0.3_dp, p, lam]
+                call stost(z, h, 1, ierr)
+                if (ierr /= 0) error stop 'Generator substep failed'
+                dp1 = z(4) - p
+                dl1 = z(5) - lam
+                m1 = m1 + dp1
+                s1 = s1 + dp1**2
+                m2 = m2 + dl1
+                s2 = s2 + dl1**2
+            end do
+            m1 = m1/n
+            m2 = m2/n
+            s1 = sqrt((s1/n - m1**2)/n)
+            s2 = sqrt((s2/n - m2**2)/n)
+            print '(a,f6.3,4es12.4)', 'moments p/p_switch, dp: ', p/p_switch, &
+                m1, fpeff*h, m2, -2.0_dp*dhh*lam*h
+            if (abs(m1 - fpeff*h) > 5.0_dp*s1) &
+                error stop 'Speed drift disagrees with fpeff'
+            if (abs(m2 + 2.0_dp*dhh*lam*h) > 5.0_dp*s2) &
+                error stop 'Pitch drift disagrees with -2*dhh*lambda'
+        end do
+    end subroutine check_generator_moments
+
+    subroutine check_maxwellian()
+        integer, parameter :: n = 4096
+        real(dp) :: z(5), sp2, sl, sl2, target, sig, duration
+        integer :: i, ierr
+        duration = 6.0_dp/rate
+        seed = 20261023
+        call random_seed(put=seed)
+        sp2 = 0.0_dp
+        sl = 0.0_dp
+        sl2 = 0.0_dp
+        do i = 1, n
+            z = [0.5_dp, 0.2_dp, 0.3_dp, 1.0e-6_dp, 0.0_dp]
+            call stost(z, duration, 1, ierr)
+            if (ierr /= 0) error stop 'Thermalization ended in an error'
+            sp2 = sp2 + z(4)**2/n
+            sl = sl + z(5)/n
+            sl2 = sl2 + z(5)**2/n
+        end do
+        target = 1.5_dp*temp/ealpha
+        ! Maxwellian: var(E) = 1.5*T**2, var(lambda) = 1/3, var(lambda**2) = 4/45.
+        sig = sqrt(1.5_dp)*temp/ealpha/sqrt(real(n, dp))
+        print '(a,3es12.4)', '<p^2>, target, sigma: ', sp2, target, sig
+        print '(a,2es12.4)', '<lambda>, <lambda^2>: ', sl, sl2
+        if (abs(sp2 - target) > 5.0_dp*sig) error stop '<p^2> off Maxwellian'
+        if (abs(sl) > 5.0_dp*sqrt(1.0_dp/3.0_dp/n)) error stop '<lambda> not 0'
+        if (abs(sl2 - 1.0_dp/3.0_dp) > 5.0_dp*sqrt(4.0_dp/45.0_dp/n)) &
+            error stop '<lambda^2> not 1/3'
+        ! Subthermal starts never error, from many seeds.
+        do i = 1, 2000
+            seed = 1000 + 7*i
+            call random_seed(put=seed)
+            z = [0.5_dp, 0.2_dp, 0.3_dp, 0.5_dp*p_switch, 0.3_dp]
+            call stost(z, 6.76_dp, 1, ierr)
+            if (ierr /= 0) error stop 'Subthermal marker errored'
+        end do
+    end subroutine check_maxwellian
+end program test_collision_low_speed