@@ -30,14 +30,32 @@ module m_thinc
3030 real (wp) :: gq3_pts(3 ) = [- 5e-1_wp * 0.7745966692414834_wp , 0._wp , 5e-1_wp * 0.7745966692414834_wp ]
3131 !> Weights: 5 / 18 , 8 / 18 , 5 / 18
3232 real (wp) :: gq3_wts(3 ) = [5._wp / 18._wp , 8._wp / 18._wp , 5._wp / 18._wp ]
33- $:GPU_DECLARE(create= ' [gq3_pts, gq3_wts]' )
33+ !> ln(2 )
34+ real (wp) :: ln2 = 0.6931471805599453_wp
35+ $:GPU_DECLARE(create= ' [gq3_pts, gq3_wts, ln2]' )
3436
3537 real (wp), allocatable, dimension (:,:,:,:) :: mthinc_nhat !> Unit normal vector
3638 real (wp), allocatable, dimension (:,:,:) :: mthinc_d !> Interface position parameter
3739 $:GPU_DECLARE(create= ' [mthinc_nhat, mthinc_d]' )
3840
3941contains
4042
43+ !> @brief Stable computation of ln(cosh (x))
44+ function f_log_cosh (x ) result(res)
45+
46+ $:GPU_ROUTINE(parallelism= ' [seq]' )
47+ real (wp), intent (in ) :: x
48+ real (wp) :: res, ax
49+
50+ ax = abs (x)
51+ if (ax > 20._wp ) then
52+ res = ax - ln2
53+ else
54+ res = ax + log (1._wp + exp (- 2._wp * ax)) - ln2
55+ end if
56+
57+ end function f_log_cosh
58+
4159 !> @brief Stable difference: log_cosh(a+ h) - log_cosh(a- h) = 2 * atanh(tanh (a)* tanh (h)). Avoids catastrophic cancellation when h
4260 !! is small relative to a.
4361 function f_log_cosh_diff (a , h ) result(res)
0 commit comments