spx_dynamics.f90 Source File

!> 动力学例程，求解单步加速度
!> @note 粒子加速度受多个粒子影响，不可简单并行加速
module spx_dynamics

    use spx_kinds, only: rk
    use spx_env, only: rgn_obj, backup, float_obj, dt
    use spx_pif_h5part, only: pif_obj
    use spx_configuration, only: spc_obj
    use spx_internal_force, only: internal_force
    use spx_external_force, only: exf_obj
    use spx_artifical_force, only: artificial_viscosity
    use spx_smoothed_kernel_function, only: skf_obj
    use spx_density_summation, only: dsm_obj
    use spx_nnps_pairs, only: pairs, nnps_grid
    implicit none

    private
    public :: solver
    logical :: is_first = .true. !! flag for first time step

contains

    !> 求解实弹性、刚性粒子单步的加速度（受力）
    !> @note 修改 rgn_obj%acc
    subroutine solver(loc, vel, acc, t)
        real(rk), intent(in) :: loc(:), vel(:), t
        real(rk), intent(out) :: acc(:)
        integer :: i

        rgn_obj%acc(:, :pif_obj%nreal) = 0.0_rk   ! 实粒子的加速度归零

        call nnps_grid%query(skf_obj%scale_hsml, pairs%items, pairs%rdx, pif_obj%numbers(4))
        call pairs%calc(itype=pif_obj%ptc%itype)

        call dsm_obj%density(rgn_obj, pairs, n=pif_obj%numbers(1))              ! 密度求和
        call internal_force(rgn_obj, pairs, rgn_obj%acc, n=pif_obj%numbers(1))  ! 内力计算(压力、粘性)
        call exf_obj%run(rgn_obj, pairs, rgn_obj%acc, n=pif_obj%nreal)          ! 外力计算

        if (spc_obj%has_artificial_viscosity) call artificial_viscosity(rgn_obj, pairs, rgn_obj%acc, &
                                                                        n=pif_obj%numbers(1)) ! 人工粘性
        call rgn_obj%restore_backup(backup, istart=pif_obj%nreal + 1, iend=pif_obj%numbers(4))

        ! @todo 手动更新刚体粒子的速度、位移，无法直接调用 leapfrog 积分器，尝试优化代码逻辑
        if (is_first) then
            !$omp parallel do private(i) schedule(dynamic)
            do i = 1, pif_obj%numbers(5)
                call float_obj(i)%get_acc(pif_obj%ptc)
                call float_obj(i)%update_particle_init(pif_obj%ptc, dt)
            end do
            is_first = .false.
        else
            !$omp parallel do private(i) schedule(dynamic)
            do i = 1, pif_obj%numbers(5)
                call float_obj(i)%get_acc(pif_obj%ptc)
                call float_obj(i)%update_particle_step(pif_obj%ptc, dt)
            end do
        end if

    end subroutine solver

end module spx_dynamics