find_first_boundary_event Subroutine

public subroutine find_first_boundary_event(cfg, x0, x1, event, status)

閉じたbox内の始点から候補終点へ向かう線分について、最初のbox面交差を返す。

Arguments

Type IntentOptional Attributes Name
type(sim_config), intent(in) :: cfg
real(kind=dp), intent(in) :: x0(3)
real(kind=dp), intent(in) :: x1(3)
type(boundary_event_type), intent(out) :: event
integer(kind=i32), intent(out) :: status

Called by

proc~~find_first_boundary_event~~CalledByGraph proc~find_first_boundary_event find_first_boundary_event proc~advance_particle_step advance_particle_step proc~advance_particle_step->proc~find_first_boundary_event proc~advance_particle_step_upper_panel_fourier advance_particle_step_upper_panel_fourier proc~advance_particle_step_upper_panel_fourier->proc~find_first_boundary_event proc~resolve_particle_boundary_candidate resolve_particle_boundary_candidate proc~resolve_particle_boundary_candidate->proc~find_first_boundary_event

Source Code

  subroutine find_first_boundary_event(cfg, x0, x1, event, status)
    type(sim_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