!> author: 左志华 !> date: 2022-07-04 !> version: alpha !> !> Smoothed kernel function <br> !> 光滑核函数 module spx_smoothed_kernel_function use spx_kinds, only: rk use seakeeping_module, only: pi use spx_math, only: r23 use spx_configuration, only: spc_obj use fffc_module, only: terminal_obj use spx_env, only: rgn_obj implicit none !> 光滑核函数 type smoothed_kernel_function_type real(rk) :: selfden !! 粒子的自身物理量加权系数, 一般用于计算粒子密度 real(rk) :: hsml !! 光滑长度 integer :: scale = 2 !! 光滑长度缩尺比例 real(rk) :: scale_hsml !! 缩尺后的光滑长度 character(:), allocatable :: type !! 光滑函数类型 procedure(kernel_fcn), pointer :: kernel => cubic_spline_kernel !! 光滑核函数句柄 real(rk), private :: factor !! 核函数常量 contains procedure :: init end type smoothed_kernel_function_type type(smoothed_kernel_function_type) :: skf_obj !! 光滑核函数对象 abstract interface !> Kernel function <br> !> 光滑核函数 pure subroutine kernel_fcn(skf, r, dx, w, dwdx) import :: rk, smoothed_kernel_function_type class(smoothed_kernel_function_type), intent(in) :: skf !! Smoothed kernel function <br> !! 光滑核函数 real(rk), intent(in) :: r !! Euclidean distance between two points <br> !! 欧氏距离 real(rk), intent(in) :: dx(2) !! dimensional distance between two points <br> !! 对应坐标轴距离 real(rk), intent(out) :: w !! kernel value <br> !! 核函数 real(rk), intent(out) :: dwdx(2) !! derivative of kernel value <br> !! 核函数导数 end subroutine kernel_fcn end interface contains !> 建立光滑核函数句柄 subroutine init(self) class(smoothed_kernel_function_type), intent(inout) :: self real(rk) :: hv(2) self%hsml = rgn_obj%hsml select case (spc_obj%smoothed_kernel_function) case ('cubic-spline') self%scale = 2 self%kernel => cubic_spline_kernel self%factor = 15.0_rk/(7.0_rk*Pi*self%hsml*self%hsml) case ('quintic-spline') self%scale = 3 self%kernel => quintic_spline_kernel self%factor = 7.0_rk/(478.0_rk*Pi*self%hsml*self%hsml) case default self%scale = 2 self%kernel => cubic_spline_kernel self%factor = 15.0_rk/(7.0_rk*Pi*self%hsml*self%hsml) call terminal_obj%warning('Unknown kernel type, use cubic-spline kernel by default') end select call self%kernel(0.0_rk, [0.0_rk, 0.0_rk], self%selfden, hv) self%scale_hsml = self%hsml*self%scale end subroutine init !> Smoothed kernel function: Cubic spline <br> !> 三次样条曲线 (monaghan 1985) <br> !> 计算光滑函数 \( W_{ij} \) 及其导数 \( \frac{dW}{dx_{ij}} \) 的子例程 pure subroutine cubic_spline_kernel(skf, r, dx, w, dwdx) class(smoothed_kernel_function_type), intent(in) :: skf real(rk), intent(in) :: r !! Euclidean distance between two points <br> !! 欧氏距离 real(rk), intent(in) :: dx(2) !! Dimensional distance between two points <br> !! 对应坐标轴距离 real(rk), intent(out) :: w !! kernel value <br> !! 核函数值 real(rk), intent(out) :: dwdx(2) !! Derivative of kernel value <br> !! 核函数导数 real(rk) :: q q = r/skf%hsml if (q <= 1.0_rk) then w = skf%factor*(r23 - q*q + q**3/2) dwdx = skf%factor*(-2.0_rk + 1.5_rk*q)*dx/(skf%hsml*skf%hsml) ! 计算机原则，先乘后除 else if (q > 1.0_rk) then w = skf%factor*(2.0_rk - q)**3/6.0_rk dwdx = -skf%factor*(2.0_rk - q)**2*dx/(2.0_rk*skf%hsml*r) end if end subroutine cubic_spline_kernel !> Smoothed kernel function: Quintic spline <br> !> 四次样条曲线 (Liu 2010) pure subroutine quintic_spline_kernel(skf, r, dx, w, dwdx) class(smoothed_kernel_function_type), intent(in) :: skf real(rk), intent(in) :: r !! Euclidean distance between two points <br> !! 欧氏距离 real(rk), intent(in) :: dx(2) !! Dimensional distance between two points <br> !! 对应坐标轴距离 real(rk), intent(out) :: w !! kernel value <br> !! 核函数值 real(rk), intent(out) :: dwdx(2) !! Derivative of kernel value <br> !! 核函数导数 real(rk) :: q q = r/skf%hsml if (q <= 1.0_rk) then w = skf%factor*((3 - q)**5 - 6*(2 - q)**5 + 15*(1 - q)**5) dwdx = skf%factor*(-120.0_rk + 120.0_rk*q - 50.0_rk*q**2)*dx/skf%hsml**2 else if (q > 1.0_rk .and. q <= 2.0_rk) then w = skf%factor*((3 - q)**5 - 6*(2 - q)**5) dwdx = skf%factor*(-5.0_rk*(3 - q)**4 + 30.0_rk*(2 - q)**4)*dx/(skf%hsml*r) else if (q > 2.0_rk .and. q <= 3.0_rk) then w = skf%factor*((3 - q)**5) dwdx = skf%factor*(-5.0_rk*(3 - q)**4)*dx/(skf%hsml*r) else w = 0.0_rk dwdx = 0.0_rk end if end subroutine quintic_spline_kernel end module spx_smoothed_kernel_function