bem_boundary.f90 Source File

!> シミュレーションボックス境界（流出/反射/面内再配置反射/周期）を適用するモジュール。
module bem_boundary
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  use bem_kinds, only: dp, i32
  use bem_types, only: particle_domain_config, bc_open, bc_reflect, bc_periodic, bc_redistributed_reflect
  use bem_status_codes, only: boundary_event_ok, boundary_event_invalid_geometry
  implicit none
  private

  integer(i32), parameter :: boundary_face_all = 63_i32

  type, public :: boundary_event_type
    logical :: has_event = .false.
    real(dp) :: fraction = 0.0_dp
    integer(i32) :: face_mask = 0_i32
    integer(i32) :: face_bc(6) = -1_i32
  end type boundary_event_type

  public :: apply_ray_box_boundary
  public :: find_first_boundary_event
  public :: apply_escape_reflect_periodic_event
contains

  !> 閉じたbox内の始点から候補終点へ向かう線分について、最初のbox面交差を返す。
  subroutine find_first_boundary_event(cfg, x0, x1, event, status)
    class(particle_domain_config), intent(in) :: cfg
    real(dp), intent(in) :: x0(3), x1(3)
    type(boundary_event_type), intent(out) :: event
    integer(i32), intent(out) :: status

    real(dp) :: delta(3), candidate_fraction(3), first_fraction, tie_tolerance
    integer(i32) :: axis, candidate_face(3), n_candidate

    event = boundary_event_type()
    status = boundary_event_ok
    if (.not. all(ieee_is_finite(x0)) .or. .not. all(ieee_is_finite(x1))) then
      status = boundary_event_invalid_geometry
      return
    end if
    if (.not. cfg%use_box) return
    if (.not. valid_event_box(cfg)) then
      status = boundary_event_invalid_geometry
      return
    end if
    if (any(x0 < cfg%box_min) .or. any(x0 > cfg%box_max)) then
      status = boundary_event_invalid_geometry
      return
    end if

    delta = x1 - x0
    if (.not. all(ieee_is_finite(delta))) then
      status = boundary_event_invalid_geometry
      return
    end if

    event%face_bc = boundary_face_conditions(cfg)
    candidate_fraction = 0.0_dp
    candidate_face = 0_i32
    n_candidate = 0_i32
    do axis = 1_i32, 3_i32
      if (delta(axis) > 0.0_dp .and. x1(axis) >= cfg%box_max(axis)) then
        n_candidate = n_candidate + 1_i32
        candidate_fraction(n_candidate) = (cfg%box_max(axis) - x0(axis))/delta(axis)
        candidate_face(n_candidate) = 2_i32*axis
      else if (delta(axis) < 0.0_dp .and. x1(axis) <= cfg%box_min(axis)) then
        n_candidate = n_candidate + 1_i32
        candidate_fraction(n_candidate) = (cfg%box_min(axis) - x0(axis))/delta(axis)
        candidate_face(n_candidate) = 2_i32*axis - 1_i32
      end if
    end do
    if (n_candidate == 0_i32) return
    if (.not. all(ieee_is_finite(candidate_fraction(1:n_candidate))) .or. &
        any(candidate_fraction(1:n_candidate) < 0.0_dp) .or. &
        any(candidate_fraction(1:n_candidate) > 1.0_dp)) then
      event = boundary_event_type()
      status = boundary_event_invalid_geometry
      return
    end if

    first_fraction = minval(candidate_fraction(1:n_candidate))
    if (first_fraction == 0.0_dp) first_fraction = 0.0_dp
    tie_tolerance = 64.0_dp*epsilon(1.0_dp)*max(1.0_dp, abs(first_fraction))
    event%has_event = .true.
    event%fraction = first_fraction
    do axis = 1_i32, n_candidate
      if (abs(candidate_fraction(axis) - first_fraction) <= tie_tolerance) then
        event%face_mask = ior(event%face_mask, shiftl(1_i32, candidate_face(axis) - 1_i32))
      end if
    end do
  end subroutine find_first_boundary_event

  !> 最初のbox eventへescape/reflect/redistributed-reflect/periodicを軸順序に依存せず適用する。
  !! face_bc を渡すと、box 設定の代わりにその face 別境界条件を event と照合する（外部障壁で一面を反射に
  !! 置き換える場合など）。省略時は cfg の bc_low/bc_high。
  subroutine apply_escape_reflect_periodic_event( &
    cfg, event, x, v, alive, escaped, status, redistribution_uniform, face_bc &
    )
    class(particle_domain_config), intent(in) :: cfg
    type(boundary_event_type), intent(in) :: event
    real(dp), intent(inout) :: x(3), v(3)
    logical, intent(inout) :: alive
    logical, intent(out) :: escaped
    integer(i32), intent(out) :: status
    real(dp), intent(in), optional :: redistribution_uniform(3)
    integer(i32), intent(in), optional :: face_bc(6)

    real(dp) :: x_work(3), v_work(3)
    integer(i32) :: axis, low_face, high_face, low_bc, high_bc, expected_face_bc(6)
    logical :: alive_work, has_open, has_redistributed_reflect

    escaped = .false.
    status = boundary_event_ok
    if (.not. cfg%use_box .or. .not. valid_event_box(cfg)) then
      status = boundary_event_invalid_geometry
      return
    end if
    if (.not. all(ieee_is_finite(x)) .or. .not. all(ieee_is_finite(v))) then
      status = boundary_event_invalid_geometry
      return
    end if
    if (.not. event%has_event) then
      if (event%face_mask /= 0_i32 .or. event%fraction /= 0.0_dp) status = boundary_event_invalid_geometry
      return
    end if
    if (.not. ieee_is_finite(event%fraction) .or. event%fraction < 0.0_dp .or. event%fraction > 1.0_dp) then
      status = boundary_event_invalid_geometry
      return
    end if
    if (event%face_mask <= 0_i32 .or. iand(event%face_mask, boundary_face_all) /= event%face_mask) then
      status = boundary_event_invalid_geometry
      return
    end if
    expected_face_bc = boundary_face_conditions(cfg)
    if (present(face_bc)) then
      if (.not. all(valid_boundary_condition(face_bc))) then
        status = boundary_event_invalid_geometry
        return
      end if
      expected_face_bc = face_bc
    end if
    if (any(event%face_bc /= expected_face_bc)) then
      status = boundary_event_invalid_geometry
      return
    end if
    do axis = 1_i32, 3_i32
      low_face = 2_i32*axis - 1_i32
      high_face = 2_i32*axis
      if (btest(event%face_mask, low_face - 1_i32) .and. btest(event%face_mask, high_face - 1_i32)) then
        status = boundary_event_invalid_geometry
        return
      end if
    end do

    x_work = x
    v_work = v
    alive_work = alive
    has_open = .false.
    has_redistributed_reflect = .false.
    do axis = 1_i32, 3_i32
      low_face = 2_i32*axis - 1_i32
      high_face = 2_i32*axis
      low_bc = event%face_bc(low_face)
      high_bc = event%face_bc(high_face)
      if (btest(event%face_mask, low_face - 1_i32)) has_open = has_open .or. low_bc == bc_open
      if (btest(event%face_mask, high_face - 1_i32)) has_open = has_open .or. high_bc == bc_open
      if (btest(event%face_mask, low_face - 1_i32)) then
        has_redistributed_reflect = has_redistributed_reflect .or. low_bc == bc_redistributed_reflect
      end if
      if (btest(event%face_mask, high_face - 1_i32)) then
        has_redistributed_reflect = has_redistributed_reflect .or. high_bc == bc_redistributed_reflect
      end if
    end do
    if (has_redistributed_reflect .and. .not. has_open) then
      if (.not. present(redistribution_uniform)) then
        status = boundary_event_invalid_geometry
        return
      end if
      if (.not. valid_redistribution_uniform(redistribution_uniform)) then
        status = boundary_event_invalid_geometry
        return
      end if
    end if

    if (has_open) then
      alive_work = .false.
      escaped = .true.
    else
      do axis = 1_i32, 3_i32
        low_face = 2_i32*axis - 1_i32
        high_face = 2_i32*axis
        low_bc = event%face_bc(low_face)
        high_bc = event%face_bc(high_face)
        if (btest(event%face_mask, low_face - 1_i32)) then
          call apply_surviving_event_axis( &
            axis, low_bc, .false., cfg%box_min(axis), cfg%box_max(axis), x_work, v_work, status &
            )
        else if (btest(event%face_mask, high_face - 1_i32)) then
          call apply_surviving_event_axis( &
            axis, high_bc, .true., cfg%box_min(axis), cfg%box_max(axis), x_work, v_work, status &
            )
        end if
        if (status /= boundary_event_ok) return
      end do
      if (has_redistributed_reflect) then
        call redistribute_event_position(cfg, event%face_mask, redistribution_uniform, x_work)
      end if
    end if
    if (.not. all(ieee_is_finite(x_work)) .or. .not. all(ieee_is_finite(v_work))) then
      escaped = .false.
      status = boundary_event_invalid_geometry
      return
    end if

    x = x_work
    v = v_work
    alive = alive_work
  end subroutine apply_escape_reflect_periodic_event

  !> 反射または周期境界の単一軸作用をwork stateへ適用する。
  subroutine apply_surviving_event_axis(axis, bc, high_side, box_min, box_max, x, v, status)
    integer(i32), intent(in) :: axis, bc
    logical, intent(in) :: high_side
    real(dp), intent(in) :: box_min, box_max
    real(dp), intent(inout) :: x(3), v(3)
    integer(i32), intent(inout) :: status
    real(dp) :: inset

    inset = boundary_inset(box_min, box_max)

    select case (bc)
    case (bc_reflect, bc_redistributed_reflect)
      if (high_side) then
        x(axis) = box_max - inset
      else
        x(axis) = box_min + inset
      end if
      v(axis) = -v(axis)
    case (bc_periodic)
      if (high_side) then
        x(axis) = box_min + inset
      else
        x(axis) = box_max - inset
      end if
    case default
      status = boundary_event_invalid_geometry
    end select
  end subroutine apply_surviving_event_axis

  !> Event APIで扱える有限かつ正幅のbox設定かを返す。
  logical function valid_event_box(cfg) result(valid)
    class(particle_domain_config), intent(in) :: cfg
    real(dp) :: inset
    integer(i32) :: axis

    valid = .false.
    if (.not. all(ieee_is_finite(cfg%box_min)) .or. .not. all(ieee_is_finite(cfg%box_max))) return
    if (any(cfg%box_max <= cfg%box_min)) return
    if (.not. all(valid_boundary_condition(cfg%bc_low)) .or. &
        .not. all(valid_boundary_condition(cfg%bc_high))) return
    do axis = 1_i32, 3_i32
      inset = boundary_inset(cfg%box_min(axis), cfg%box_max(axis))
      if (.not. ieee_is_finite(inset) .or. inset <= 0.0_dp) return
      if (cfg%box_min(axis) + inset >= cfg%box_max(axis) - inset) return
    end do
    valid = .true.
  end function valid_event_box

  !> 境界作用後の座標をnormal範囲かつevent toleranceより十分内側へ置く。
  pure real(dp) function boundary_inset(box_min, box_max) result(inset)
    real(dp), intent(in) :: box_min, box_max
    real(dp) :: span, scale

    span = box_max - box_min
    scale = max(abs(box_min), abs(box_max), span, tiny(1.0_dp))
    inset = max(64.0_dp*epsilon(1.0_dp)*scale, spacing(scale))
    inset = min(0.25_dp*span, inset)
  end function boundary_inset

  !> 境界条件enum配列の各値が既知かを返す。
  elemental logical function valid_boundary_condition(bc) result(valid)
    integer(i32), intent(in) :: bc

    valid = bc == bc_open .or. bc == bc_reflect .or. bc == bc_periodic .or. bc == bc_redistributed_reflect
  end function valid_boundary_condition

  !> 再配置反射に供給された各軸の一様乱数が閉区間 `[0,1]` にあるかを返す。
  pure logical function valid_redistribution_uniform(uniform) result(valid)
    real(dp), intent(in) :: uniform(3)

    valid = all(ieee_is_finite(uniform)) .and. all(uniform >= 0.0_dp) .and. all(uniform <= 1.0_dp)
  end function valid_redistribution_uniform

  !> 同時eventの法線軸を固定し、それ以外の座標をbox面内で一様再配置する。
  pure subroutine redistribute_event_position(cfg, face_mask, uniform, x)
    class(particle_domain_config), intent(in) :: cfg
    integer(i32), intent(in) :: face_mask
    real(dp), intent(in) :: uniform(3)
    real(dp), intent(inout) :: x(3)
    integer(i32) :: axis
    real(dp) :: inset

    do axis = 1_i32, 3_i32
      if (btest(face_mask, 2_i32*axis - 2_i32) .or. btest(face_mask, 2_i32*axis - 1_i32)) cycle
      inset = boundary_inset(cfg%box_min(axis), cfg%box_max(axis))
      x(axis) = cfg%box_min(axis) + inset + &
                uniform(axis)*(cfg%box_max(axis) - cfg%box_min(axis) - 2.0_dp*inset)
    end do
  end subroutine redistribute_event_position

  !> face bit順にbox設定の境界条件を並べる。
  pure function boundary_face_conditions(cfg) result(face_bc)
    class(particle_domain_config), intent(in) :: cfg
    integer(i32) :: face_bc(6)
    integer(i32) :: axis

    do axis = 1_i32, 3_i32
      face_bc(2_i32*axis - 1_i32) = cfg%bc_low(axis)
      face_bc(2_i32*axis) = cfg%bc_high(axis)
    end do
  end function boundary_face_conditions

  !> photo raycast の ray が box 面を越えた点に反射・周期境界を適用する。粒子の境界処理は event API
  !! （find_first_boundary_event / apply_escape_reflect_periodic_event）が担い、これは ray 専用。
  !! ray は面内再配置の乱数を持たないので、流出面と再配置反射の面では ray を止める（alive=.false.）。
  !! box が無効なら何もせず、幅 0 の軸があれば ray を止める。
  subroutine apply_ray_box_boundary(cfg, x, dir, alive)
    class(particle_domain_config), intent(in) :: cfg
    real(dp), intent(inout) :: x(3), dir(3)
    logical, intent(inout) :: alive

    real(dp), parameter :: eps = 1.0d-12
    integer(i32) :: axis
    real(dp) :: span

    if (.not. cfg%use_box) return
    if (.not. alive) return

    do axis = 1, 3
      span = cfg%box_max(axis) - cfg%box_min(axis)
      if (span <= 0.0_dp) then
        alive = .false.
        return
      end if

      do while (x(axis) < cfg%box_min(axis))
        call apply_ray_face(cfg%bc_low(axis))
        if ((.not. alive) .or. x(axis) >= cfg%box_min(axis)) exit
      end do
      if (.not. alive) return

      do while (x(axis) > cfg%box_max(axis))
        call apply_ray_face(cfg%bc_high(axis))
        if ((.not. alive) .or. x(axis) <= cfg%box_max(axis)) exit
      end do
      if (.not. alive) return
    end do

  contains

    subroutine apply_ray_face(bc)
      integer(i32), intent(in) :: bc

      select case (bc)
      case (bc_reflect)
        call reflect_one_side(axis, cfg%box_min(axis), cfg%box_max(axis), eps, x, dir)
      case (bc_periodic)
        x(axis) = modulo(x(axis) - cfg%box_min(axis), span) + cfg%box_min(axis)
        x(axis) = min(cfg%box_max(axis) - eps, max(cfg%box_min(axis) + eps, x(axis)))
      case default
        alive = .false.
      end select
    end subroutine apply_ray_face
  end subroutine apply_ray_box_boundary

  !> 反射境界と同じ鏡像位置・法線速度反転を適用する。
  subroutine reflect_one_side(axis, box_min, box_max, eps, x, v)
    integer(i32), intent(in) :: axis
    real(dp), intent(in) :: box_min, box_max, eps
    real(dp), intent(inout) :: x(3), v(3)

    if (x(axis) < box_min) then
      x(axis) = box_min + (box_min - x(axis))
    else if (x(axis) > box_max) then
      x(axis) = box_max - (x(axis) - box_max)
    end if
    x(axis) = min(box_max - eps, max(box_min + eps, x(axis)))
    v(axis) = -v(axis)
  end subroutine reflect_one_side

end module bem_boundary