compute_laplace_derivatives Subroutine

public subroutine compute_laplace_derivatives(plan, r, deriv)

ラプラス核の多重微分係数を計算する。 1/sqrt(|r|^2 + s^2) の全次数の偏微分 ∂^α f を漸化式で求める。 f = a^{-1/2}（a = |r|^2 + s^2）は a ∂_i f = -x_i f を満たし、Leibniz 則で ∂^β を当てると D(β+e_i) = -[ x_i D(β) + β_i D(β-e_i) + 2Σ_j β_j x_j D(β+e_i-e_j) + Σ_j β_j(β_j-1) D(β+e_i-2e_j) ] / a となる。multi-index は次数昇順に並んでいるので、1 回の走査で全次数が決まる。 旧来の多項式乗算版 compute_laplace_derivatives_reference と同じ値を返す（契約テストで検査）。

Arguments

Type IntentOptional Attributes Name
type(fmm_plan_type), intent(in) :: plan

多重指数テーブルを持つ FMM plan。
real(kind=dp), intent(in) :: r(3)
real(kind=dp), intent(out) :: deriv(:)

Called by

proc~~compute_laplace_derivatives~~CalledByGraph proc~compute_laplace_derivatives compute_laplace_derivatives proc~core_complete_plan_impl core_complete_plan_impl proc~core_complete_plan_impl->proc~compute_laplace_derivatives proc~core_build_panel_plan_impl core_build_panel_plan_impl proc~core_build_panel_plan_impl->proc~core_complete_plan_impl proc~core_build_plan_impl core_build_plan_impl proc~core_build_plan_impl->proc~core_complete_plan_impl

Source Code

  subroutine compute_laplace_derivatives(plan, r, deriv)
    type(fmm_plan_type), intent(in) :: plan
    real(dp), intent(in) :: r(3)
    real(dp), intent(out) :: deriv(:)
    real(dp) :: a, soft2, value
    integer(i32) :: idx, axis, j, alpha(3), beta(3), lower(3)

    if (size(deriv) /= plan%nderiv) error stop 'Derivative buffer size mismatch.'
    soft2 = plan%options%softening*plan%options%softening
    a = dot_product(r, r) + soft2
    if (a <= tiny(1.0d0)) error stop 'Cannot expand Laplace kernel at r=0.'

    deriv(plan%deriv_map(0_i32, 0_i32, 0_i32)) = 1.0d0/sqrt(a)
    do idx = 1_i32, plan%nderiv
      alpha = plan%deriv_alpha(:, idx)
      if (sum(alpha) == 0_i32) cycle
      axis = 1_i32
      do while (alpha(axis) == 0_i32)
        axis = axis + 1_i32
      end do
      beta = alpha
      beta(axis) = beta(axis) - 1_i32
      value = r(axis)*deriv(plan%deriv_map(beta(1), beta(2), beta(3)))
      if (beta(axis) > 0_i32) then
        lower = beta
        lower(axis) = lower(axis) - 1_i32
        value = value + real(beta(axis), dp)*deriv(plan%deriv_map(lower(1), lower(2), lower(3)))
      end if
      do j = 1_i32, 3_i32
        if (beta(j) == 0_i32) cycle
        lower = alpha
        lower(j) = lower(j) - 1_i32
        value = value + 2.0d0*real(beta(j), dp)*r(j)*deriv(plan%deriv_map(lower(1), lower(2), lower(3)))
        if (beta(j) > 1_i32) then
          lower(j) = lower(j) - 1_i32
          value = value + real(beta(j)*(beta(j) - 1_i32), dp)*deriv(plan%deriv_map(lower(1), lower(2), lower(3)))
        end if
      end do
      deriv(idx) = -value/a
    end do
  end subroutine compute_laplace_derivatives