internal_force_without_visc Subroutine

public pure subroutine internal_force_without_visc(particle, pairs, acc, n)

Internal force function without viscosity
内力求解

Arguments

Type IntentOptional Attributes Name
type(particle_type), intent(inout) :: particle

区域

type(nnps_pairs_type), intent(in) :: pairs

粒子对

real(kind=rk), intent(inout) :: acc(:,:)

加速度

integer, intent(in) :: n

实弹性粒子数, numbers(1)


Contents


Source Code

    pure subroutine internal_force_without_visc(particle, pairs, acc, n)
        type(particle_type), intent(inout) :: particle  !! 区域
        type(nnps_pairs_type), intent(in) :: pairs   !! 粒子对
        real(rk), intent(inout) :: acc(:, :)    !! 加速度
        integer, intent(in) :: n                !! 实弹性粒子数, numbers(1)
        integer :: i
        real(rk) :: h(2)

        ! ----------------------------- 根据 EoS 计算压强信息 -------------------------------- !
        do concurrent(i=1:n)

            select case (particle%itype(i))  ! TOCHECK: 虚粒子是否需要基于 EoS 计算压强
            case (50:59); call eow_obj%p_water(particle%rho(i), particle%p(i), particle%c(i))             ! 淡水
            case (110:119); call p_ideal_gas(particle%rho(i), particle%u(i), particle%p(i), particle%c(i))  ! 理想气体
            end select

        end do

        ! ------------------------------------ 计算压力 ------------------------------------ !
        do i = 1, pairs%len
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i))
                if (pairs%contact_type(i) /= 2) cycle ! 仅弹性实粒子参与内力计算

                ! 用到了虚粒子的压强对实粒子的影响,防止边界处实粒子有较大的、指向边界的初始压强合力
                h(:) = -(particle%p(ip) + particle%p(jp))*pairs%wdwdx(2:3, i)/ &
                       &(particle%rho(ip)*particle%rho(jp))
                acc(:, ip) = acc(:, ip) + h(:)*particle%mass(jp)
                acc(:, jp) = acc(:, jp) - h(:)*particle%mass(ip)
            end associate
        end do

    end subroutine internal_force_without_visc