build_periodic_shell_mean_plan Subroutine

public subroutine build_periodic_shell_mean_plan(periodic_axes, free_axis, cell_center, cell_len, image_layers, softening, z_min, z_max, spacing, vertices, plan)

線源ごとの s を高さの格子に表として持つ。

Arguments

Type IntentOptional Attributes Name
integer(kind=i32), intent(in) :: periodic_axes(2)
integer(kind=i32), intent(in) :: free_axis
real(kind=dp), intent(in) :: cell_center(2)
real(kind=dp), intent(in) :: cell_len(2)
integer(kind=i32), intent(in) :: image_layers
real(kind=dp), intent(in) :: softening
real(kind=dp), intent(in) :: z_min
real(kind=dp), intent(in) :: z_max
real(kind=dp), intent(in) :: spacing
real(kind=dp), intent(in) :: vertices(:,:,:)
type(periodic_shell_mean_plan_type), intent(inout) :: plan

Calls

proc~~build_periodic_shell_mean_plan~~CallsGraph proc~build_periodic_shell_mean_plan build_periodic_shell_mean_plan proc~periodic_shell_mean_point_excess periodic_shell_mean_point_excess proc~build_periodic_shell_mean_plan->proc~periodic_shell_mean_point_excess proc~reset_periodic_shell_mean_plan reset_periodic_shell_mean_plan proc~build_periodic_shell_mean_plan->proc~reset_periodic_shell_mean_plan proc~periodic_shell_mean_excess periodic_shell_mean_excess proc~periodic_shell_mean_point_excess->proc~periodic_shell_mean_excess

Called by

proc~~build_periodic_shell_mean_plan~~CalledByGraph proc~build_periodic_shell_mean_plan build_periodic_shell_mean_plan proc~core_complete_plan_impl core_complete_plan_impl proc~core_complete_plan_impl->proc~build_periodic_shell_mean_plan 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 build_periodic_shell_mean_plan( &
    periodic_axes, free_axis, cell_center, cell_len, image_layers, softening, z_min, z_max, spacing, vertices, plan &
    )
    integer(i32), intent(in) :: periodic_axes(2), free_axis, image_layers
    real(dp), intent(in) :: cell_center(2), cell_len(2), softening, z_min, z_max, spacing
    real(dp), intent(in) :: vertices(:, :, :)
    type(periodic_shell_mean_plan_type), intent(inout) :: plan
    integer(i32) :: nsrc, k, i, q, nq
    real(dp) :: point(3), weight, s, dsdh, z
    logical :: is_point

    call reset_periodic_shell_mean_plan(plan)
    if (z_max < z_min .or. spacing <= 0.0_dp) error stop 'periodic shell-mean plan has an invalid height range.'
    nsrc = int(size(vertices, 3), i32)
    plan%periodic_axes = periodic_axes
    plan%free_axis = free_axis
    plan%cell_lo = cell_center - 0.5_dp*cell_len
    plan%cell_len = cell_len
    plan%image_layers = image_layers
    plan%soft2 = softening*softening
    plan%nz = max(2_i32, int(ceiling((z_max - z_min)/spacing), i32) + 1_i32)
    plan%z0 = z_min
    plan%dz = max(z_max - z_min, tiny(1.0_dp))/real(plan%nz - 1_i32, dp)
    allocate (plan%value(plan%nz, nsrc), plan%slope(plan%nz, nsrc))
    plan%value = 0.0_dp
    plan%slope = 0.0_dp
    !$omp parallel do default(shared) private(i, k, q, nq, point, weight, s, dsdh, z, is_point) schedule(static)
    do i = 1_i32, nsrc
      is_point = all(vertices(:, 2, i) == vertices(:, 1, i)) .and. all(vertices(:, 3, i) == vertices(:, 1, i))
      nq = merge(1_i32, 7_i32, is_point)
      do q = 1_i32, nq
        if (is_point) then
          point = vertices(:, 1, i)
          weight = 1.0_dp
        else
          point = dunavant_a(q)*vertices(:, 1, i) + dunavant_b(q)*vertices(:, 2, i) + &
                  (1.0_dp - dunavant_a(q) - dunavant_b(q))*vertices(:, 3, i)
          weight = dunavant_w(q)
        end if
        do k = 1_i32, plan%nz
          z = plan%z0 + real(k - 1_i32, dp)*plan%dz
          call periodic_shell_mean_point_excess(plan, point, z, s, dsdh)
          plan%value(k, i) = plan%value(k, i) + weight*s
          plan%slope(k, i) = plan%slope(k, i) + weight*dsdh
        end do
      end do
    end do
    !$omp end parallel do
    plan%ready = .true.
  end subroutine build_periodic_shell_mean_plan