spx_density_summation.f90 Source File

!> author: 左志华
!> date: 2022-06-17
!>
!> 针对实粒子进行密度求和：<br>
!> 1. 连续密度法；<br>
!> 2. 正则密度求和法；<br>
!> 3. 密度求和法。
module spx_density_summation

    use spx_kinds, only: rk
    use spx_particle, only: particle_type
    use spx_nnps_pairs, only: nnps_pairs_type
    use spx_smoothed_kernel_function, only: skf_obj
    use spx_configuration, only: spc_obj
    implicit none

    private
    public :: dsm_obj

    !> 密度求和
    type density_summation
        procedure(density_summation_fcn), pointer, nopass :: density  !! 密度求和函数指针
    contains
        procedure :: init
    end type density_summation

    type(density_summation) :: dsm_obj  !! 密度求和对象

    abstract interface
        !> Density summation <br>
        !> 密度求和
        subroutine density_summation_fcn(particle, pairs, n)
            import :: particle_type, nnps_pairs_type
            type(particle_type), intent(inout) :: particle      !! particle <br>
                                                                !! 区域
            type(nnps_pairs_type), intent(in) :: pairs               !! particle pairs <br>
                                                                !! 粒子对
            integer, intent(in) :: n                            !! 实弹性粒子数, numbers(1)
        end subroutine density_summation_fcn

    end interface

contains

    !> 初始化
    subroutine init(self)
        class(density_summation), intent(inout) :: self

        select case (spc_obj%density_summation_method)
        case ('continuous')
            self%density => continuous_density_summation
        case ('normalized')
            self%density => normalized_density_summation
        case ('denormalized')
            self%density => denormalized_density_summation
        case default
            self%density => denormalized_density_summation
        end select

    end subroutine init

    !> 对所有实粒子粒子进行非正则化密度求和
    subroutine denormalized_density_summation(particle, pairs, n)
        type(particle_type), intent(inout) :: particle      !! 区域
        type(nnps_pairs_type), intent(in) :: pairs          !! 粒子对
        integer, intent(in) :: n                            !! 实弹性粒子数 \( N_{real} \)
        integer :: i

        particle%rho(:n) = skf_obj%selfden*particle%mass(:n)              ! 计算自身对密度的影响，并初始化当前步长密度值

        do i = 1, pairs%len                        ! 计算密度求和
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i))

                if (ip <= n .and. jp <= n) then
                    if (particle%itype(ip) /= particle%itype(jp)) cycle
                    particle%rho(ip) = particle%rho(ip) + particle%mass(jp)*pairs%wdwdx(1, i)
                    particle%rho(jp) = particle%rho(jp) + particle%mass(ip)*pairs%wdwdx(1, i)
                end if

            end associate
        end do

    end subroutine denormalized_density_summation

    !> Normalized density summation <br>
    !> 对所有粒子进行正则化密度求和法
    !> @note 计算量较密度求和法大
    pure subroutine normalized_density_summation(particle, pairs, n)
        type(particle_type), intent(inout) :: particle  !! particle <br>
                                                        !! 区域
        type(nnps_pairs_type), intent(in) :: pairs           !! particle pairs <br>
                                                        !! 粒子对
        integer, intent(in) :: n                        !! 实弹性粒子数, numbers(1)
                                                        !! @note 用以定义本地数组长度
        integer :: i
        real(rk) :: rho_tmp(n), wi(n)                   !! 分母部分

        rho_tmp(:n) = particle%rho(:n)                  ! 临时存储密度值

        particle%rho(:n) = skf_obj%selfden*particle%mass(:n)
        wi(:n) = particle%rho(:n)/rho_tmp(:n)

        do concurrent(i=1:pairs%len)            ! 计算密度求和
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i))

                if (ip <= n .and. jp <= n) then
                    if (particle%itype(ip) /= particle%itype(jp)) cycle
                    particle%rho(ip) = particle%rho(ip) + particle%mass(jp)*pairs%wdwdx(1, i)
                    particle%rho(jp) = particle%rho(jp) + particle%mass(ip)*pairs%wdwdx(1, i)
                    wi(ip) = wi(ip) + particle%mass(jp)*pairs%wdwdx(1, i)/rho_tmp(jp)
                    wi(jp) = wi(jp) + particle%mass(ip)*pairs%wdwdx(1, i)/rho_tmp(ip)
                end if

            end associate
        end do

        particle%rho(:n) = particle%rho(:n)/wi(:n)  ! 正则化密度

    end subroutine normalized_density_summation

    !> Continuous density summation <br>
    !> 基于连续性方程的连续密度法，计算密度的增量以在时间积分算法中更新密度值
    !> @todo itype 相同才进行密度求和？
    pure subroutine continuous_density_summation(particle, pairs, n)
        type(particle_type), intent(inout) :: particle  !! particle <br>
                                                        !! 区域
        type(nnps_pairs_type), intent(in) :: pairs           !! particle pairs <br>
                                                        !! 粒子对
        integer, intent(in) :: n                        !! 实弹性粒子数, numbers(1)
        integer :: i
        real(rk) :: dvel(2)

        particle%drho(:n) = 0.0_rk
        do concurrent(i=1:pairs%len)
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i))

                if (ip <= n .and. jp <= n) then
                    if (particle%itype(ip) /= particle%itype(jp)) cycle
                    dvel(:) = particle%vel(:, ip) - particle%vel(:, jp)

                    associate (vcc => sum(dvel(:)*pairs%wdwdx(2:3, i)))
                        particle%drho(ip) = particle%drho(ip) + particle%mass(jp)*vcc
                        particle%drho(jp) = particle%drho(jp) + particle%mass(ip)*vcc
                    end associate
                end if

            end associate
        end do

    end subroutine continuous_density_summation

end module spx_density_summation