sample_reservoir_face_particles Subroutine

public subroutine sample_reservoir_face_particles(box_min, box_max, inject_face, pos_low, pos_high, drift_velocity, m_particle, temperature_k, batch_duration, x, v, barrier_normal_energy, vmin_normal, position_jitter_dt, apply_barrier_energy_shift)

上流リザーバ境界から流入する粒子群を面注入としてサンプルする。

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: box_min(3)
real(kind=dp), intent(in) :: box_max(3)
character(len=*), intent(in) :: inject_face

注入面識別子（x_low/x_high/y_low/y_high/z_low/z_high）。
real(kind=dp), intent(in) :: pos_low(3)
real(kind=dp), intent(in) :: pos_high(3)
real(kind=dp), intent(in) :: drift_velocity(3)
real(kind=dp), intent(in) :: m_particle

粒子1個あたりの質量 [kg]。
real(kind=dp), intent(in) :: temperature_k

粒子1個あたりの質量 [kg]。 温度 [K]。
real(kind=dp), intent(in) :: batch_duration

粒子1個あたりの質量 [kg]。 温度 [K]。 1バッチの物理時間長 [s]（現在は妥当性チェックのみ）。
real(kind=dp), intent(out) :: x(:,:)
real(kind=dp), intent(out) :: v(:,:)
real(kind=dp), intent(in), optional :: barrier_normal_energy

法線方向のエネルギー障壁 2 q Δφ / m [m^2/s^2]（省略時 0）。
real(kind=dp), intent(in), optional :: vmin_normal

法線速度の下限 [m/s]（省略時は barrier_normal_energy から自動導出）。
real(kind=dp), intent(in), optional :: position_jitter_dt

初期位置に速度方向で与えるランダムジッタ時間幅[s]（省略時は 0）。
logical, intent(in), optional :: apply_barrier_energy_shift

.true. のとき法線速度へ障壁エネルギー変換を適用する。

Calls

proc~~sample_reservoir_face_particles~~CallsGraph proc~sample_reservoir_face_particles sample_reservoir_face_particles proc~sample_shifted_maxwell_velocities sample_shifted_maxwell_velocities proc~sample_reservoir_face_particles->proc~sample_shifted_maxwell_velocities

Called by

proc~~sample_reservoir_face_particles~~CalledByGraph proc~sample_reservoir_face_particles sample_reservoir_face_particles proc~sample_species_state sample_species_state proc~sample_species_state->proc~sample_reservoir_face_particles proc~init_particle_batch_from_config init_particle_batch_from_config proc~init_particle_batch_from_config->proc~sample_species_state proc~init_particles_from_config init_particles_from_config proc~init_particles_from_config->proc~sample_species_state

Source Code

  subroutine sample_reservoir_face_particles( &
    box_min, box_max, inject_face, pos_low, pos_high, drift_velocity, m_particle, temperature_k, batch_duration, x, v, &
    barrier_normal_energy, vmin_normal, position_jitter_dt, apply_barrier_energy_shift &
    )
    real(dp), intent(in) :: box_min(3), box_max(3)
    character(len=*), intent(in) :: inject_face
    real(dp), intent(in) :: pos_low(3), pos_high(3), drift_velocity(3)
    real(dp), intent(in) :: m_particle, temperature_k, batch_duration
    real(dp), intent(out) :: x(:, :)
    real(dp), intent(out) :: v(:, :)
    real(dp), intent(in), optional :: barrier_normal_energy
    real(dp), intent(in), optional :: vmin_normal
    real(dp), intent(in), optional :: position_jitter_dt
    logical, intent(in), optional :: apply_barrier_energy_shift

    integer :: i, axis_n, axis_t1, axis_t2
    real(dp) :: boundary_value, inward_normal(3), sigma, u_n, vn_floor, barrier, vn_inf, jitter_dt
    real(dp), allocatable :: u(:, :), tau(:)
    logical :: apply_energy_shift

    if (size(x, 1) /= 3 .or. size(v, 1) /= 3) error stop "reservoir particle arrays must have first dimension 3"
    if (size(x, 2) /= size(v, 2)) error stop "reservoir x/v size mismatch"
    if (batch_duration < 0.0_dp) error stop "batch_duration must be >= 0"
    if (present(position_jitter_dt)) then
      if (position_jitter_dt < 0.0_dp) error stop "position_jitter_dt must be >= 0"
      jitter_dt = position_jitter_dt
    else
      jitter_dt = 0.0_dp
    end if
    apply_energy_shift = .true.
    if (present(apply_barrier_energy_shift)) apply_energy_shift = apply_barrier_energy_shift
    if (size(x, 2) == 0) return

    call resolve_face_geometry(box_min, box_max, inject_face, axis_n, boundary_value, inward_normal)
    call resolve_face_axes(inject_face, axis_t1, axis_t2)

    sigma = sqrt(k_boltzmann*temperature_k/m_particle)
    u_n = dot_product(drift_velocity, inward_normal)
    barrier = 0.0_dp
    if (present(barrier_normal_energy)) barrier = barrier_normal_energy
    vn_floor = 0.0_dp
    if (barrier > 0.0_dp) vn_floor = sqrt(barrier)
    if (present(vmin_normal)) vn_floor = max(vn_floor, max(0.0_dp, vmin_normal))

    call sample_shifted_maxwell_velocities(drift_velocity, m_particle, v, temperature_k=temperature_k)
    call sample_flux_weighted_normal_component(u_n, sigma, v(axis_n, :), vmin_normal=vn_floor)
    do i = 1, size(v, 2)
      if (apply_energy_shift) then
        vn_inf = v(axis_n, i)
        vn_inf = sqrt(max(0.0_dp, vn_inf*vn_inf - barrier))
        v(axis_n, i) = inward_normal(axis_n)*vn_inf
      else
        v(axis_n, i) = inward_normal(axis_n)*v(axis_n, i)
      end if
    end do

    allocate (u(2, size(x, 2)))
    call random_number(u)
    if (jitter_dt > 0.0_dp) then
      allocate (tau(size(x, 2)))
      call random_number(tau)
    end if

    do i = 1, size(x, 2)
      x(:, i) = 0.0_dp
      x(axis_n, i) = boundary_value
      x(axis_t1, i) = pos_low(axis_t1) + (pos_high(axis_t1) - pos_low(axis_t1))*u(1, i)
      x(axis_t2, i) = pos_low(axis_t2) + (pos_high(axis_t2) - pos_low(axis_t2))*u(2, i)
      if (jitter_dt > 0.0_dp) x(:, i) = x(:, i) + v(:, i)*(tau(i)*jitter_dt)
      x(:, i) = x(:, i) + inward_normal*1.0d-12
    end do
  end subroutine sample_reservoir_face_particles