spx_external_force.f90 Source File

!> author: 左志华
!> date: 2022-07-02
!> license: BSD 3-Clause
!>
!> External forces <br>
!> 外力：1. 重力；2. 边界排斥力。
module spx_external_force

    use spx_kinds, only: rk
    use seakeeping_module, only: g
    use spx_particle, only: particle_type
    use spx_pif_h5part, only: pif_obj
    use spx_nnps_pairs, only: nnps_pairs_type
    use spx_configuration, only: spc_obj
    use spx_env, only: rgn_obj
    implicit none

    private
    public :: exf_obj

    !> 外力求解
    type external_force_type
        procedure(external_force_fcn), pointer :: run
        real(rk) :: r0  !! 排斥力半径
        real(rk) :: dd  !! 速度最大值的平方量级
    contains
        procedure :: init
    end type external_force_type

    type(external_force_type) :: exf_obj

    integer, parameter :: p1 = 12, p2 = 4   !! p1, p2 are used to calculate the boundary repulsion <br>
                                            !! p1, p2 用来计算边界排斥力
    abstract interface

        !> Solve the external force <br>
        !> 外力求解，涉及重力和边界力
        subroutine external_force_fcn(self, particle, pairs, acc, n)
            import :: particle_type, nnps_pairs_type, rk, external_force_type
            class(external_force_type), intent(inout) :: self
            type(particle_type), intent(in) :: particle     !! particle <br>
                                                            !! 区域
            type(nnps_pairs_type), intent(in) :: pairs           !! Pairs <br>
                                                            !! 粒子对
            real(rk), intent(inout) :: acc(:, :)            !! Acceleration <br>
                                                            !! 加速度
            integer, intent(in) :: n                        !! 实粒子数, nreal
        end subroutine external_force_fcn

    end interface

contains

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

        if (spc_obj%has_gravity) then
            self%run => external_force_with_gravity
        else
            self%run => external_force_without_gravity
        end if

        self%r0 = 0.9_rk*rgn_obj%hsml
        self%dd = (spc_obj%c0/2)**2

    end subroutine init

    !> 含重力、排斥力的外力求解
    subroutine external_force_with_gravity(self, particle, pairs, acc, n)
        class(external_force_type), intent(inout) :: self
        type(particle_type), intent(in) :: particle     !! 区域
        type(nnps_pairs_type), intent(in) :: pairs      !! 粒子对
        real(rk), intent(inout) :: acc(:, :)            !! 加速度 \( F_{external} = -g + F_{boundary} \)
        integer, intent(in) :: n                        !! 实粒子数 \( N_{real} \)
        integer :: i

        ! -------------------------------- 计算重力 --------------------------------- !
        acc(2, :n) = acc(2, :n) - g   ! 仅对实体粒子起作用，虚粒子不需要求解加速度而运动

        ! ----------------------------- 计算边界排斥力 ------------------------------- !
        do i = 1, pairs%len
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i), r => pairs%rdx(3*i - 2))
                select case (pairs%contact_type(i))
                case (3:4)
                    if (r < self%r0) then
                        associate (f => ((self%r0/r)**p1 - (self%r0/r)**p2)/(r*r), &
                                   dx => pairs%rdx(3*i - 1:3*i))
                            acc(:, ip) = acc(:, ip) + self%dd*dx(:)*f
                            acc(:, jp) = acc(:, jp) - self%dd*dx(:)*f
                        end associate
                    end if
                end select
            end associate
        end do

    end subroutine external_force_with_gravity

    !> 排斥力的外力求解 (无重力)
    pure subroutine external_force_without_gravity(self, particle, pairs, acc, n)
        class(external_force_type), intent(inout) :: self
        type(particle_type), intent(in) :: particle      !! 区域
        type(nnps_pairs_type), intent(in) :: pairs       !! 粒子对
        real(rk), intent(inout) :: acc(:, :)             !! 加速度 \( F_{external} = F_{boundary} \)
        integer, intent(in) :: n                         !! 实粒子数 \( N_{real} \)
        integer :: i

        ! ----------------------------- 计算边界排斥力 ------------------------------- !
        do i = 1, pairs%len
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i), r => pairs%rdx(3*i - 2))
                select case (pairs%contact_type(i))
                case (3:4)
                    if (r < self%r0) then
                        associate (f => ((self%r0/r)**p1 - (self%r0/r)**p2)/(r*r), &
                                   dx => pairs%rdx(3*i - 1:3*i))
                            acc(:, ip) = acc(:, ip) + self%dd*dx(:)*f
                            acc(:, jp) = acc(:, jp) - self%dd*dx(:)*f
                        end associate
                    end if
                end select
            end associate
        end do

    end subroutine external_force_without_gravity

end module spx_external_force