sample_reservoir_face_particle_counter Subroutine

public subroutine sample_reservoir_face_particle_counter(box_min, box_max, inject_face, pos_low, pos_high, drift_velocity, m_particle, temperature_k, barrier_normal_energy, vmin_normal, rng_seed, rng_counter, x, v)

上流リザーバ境界から流入する粒子を 1 つ、counter-based 乱数で面注入としてサンプルする。 同じ (seed, counter) は常に同じ粒子を返すので、どの rank・thread が生成しても結果は同じ。 速度の法線成分は flux-weighted に draw 1、面内位置は draw 2, 3、熱速度は draw 100 以降を使う。

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
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
real(kind=dp), intent(in) :: temperature_k
real(kind=dp), intent(in) :: barrier_normal_energy

法線エネルギー障壁 2 q Δφ / m [m^2/s^2]（0 なら無し）。
real(kind=dp), intent(in) :: vmin_normal

法線エネルギー障壁 2 q Δφ / m [m^2/s^2]（0 なら無し）。 無限遠法線速度の下限 [m/s]。
integer(kind=i32), intent(in) :: rng_seed
integer(kind=i64), intent(in) :: rng_counter(:)
real(kind=dp), intent(out) :: x(3)
real(kind=dp), intent(out) :: v(3)

Calls

proc~~sample_reservoir_face_particle_counter~~CallsGraph proc~sample_reservoir_face_particle_counter sample_reservoir_face_particle_counter proc~counter_uniforms counter_uniforms proc~sample_reservoir_face_particle_counter->proc~counter_uniforms proc~resolve_face_axes resolve_face_axes proc~sample_reservoir_face_particle_counter->proc~resolve_face_axes proc~resolve_face_geometry resolve_face_geometry proc~sample_reservoir_face_particle_counter->proc~resolve_face_geometry proc~sample_flux_weighted_normal_component sample_flux_weighted_normal_component proc~sample_reservoir_face_particle_counter->proc~sample_flux_weighted_normal_component proc~sample_standard_normal sample_standard_normal proc~sample_reservoir_face_particle_counter->proc~sample_standard_normal none~next_uniform counter_stream_type%next_uniform proc~counter_uniforms->none~next_uniform interface~random_uniform random_uniform proc~sample_standard_normal->interface~random_uniform proc~sample_standard_normal->none~next_uniform pcg64_random_double53 pcg64_random_double53 none~next_uniform->pcg64_random_double53 proc~open_uniform_from_double53 open_uniform_from_double53 none~next_uniform->proc~open_uniform_from_double53

Source Code

  subroutine sample_reservoir_face_particle_counter( &
    box_min, box_max, inject_face, pos_low, pos_high, drift_velocity, m_particle, temperature_k, &
    barrier_normal_energy, vmin_normal, rng_seed, rng_counter, x, v &
    )
    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, barrier_normal_energy, vmin_normal
    integer(i32), intent(in) :: rng_seed
    integer(i64), intent(in) :: rng_counter(:)
    real(dp), intent(out) :: x(3), v(3)
    integer :: axis_n, axis_t1, axis_t2
    real(dp) :: boundary_value, inward_normal(3), sigma, u_n, vn_floor, vn_inf, vn(1), z(3, 1), uniform(3)

    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)
    if (m_particle <= 0.0_dp) error stop "particle mass [kg] must be > 0"
    if (temperature_k < 0.0_dp) error stop "temperature [K] must be >= 0"
    sigma = sqrt(k_boltzmann*temperature_k/m_particle)
    u_n = dot_product(drift_velocity, inward_normal)
    vn_floor = 0.0_dp
    if (barrier_normal_energy > 0.0_dp) vn_floor = sqrt(barrier_normal_energy)
    vn_floor = max(vn_floor, max(0.0_dp, vmin_normal))

    call sample_standard_normal(z, sigma_cutoff=default_velocity_sigma_cutoff, rng_seed=rng_seed, rng_counter=rng_counter)
    v = sigma*z(:, 1) + drift_velocity
    ! Draw 1 picks the normal speed and draws 2 and 3 the position on the face.
    call counter_uniforms(rng_seed, rng_counter, 1_i64, uniform)
    call sample_flux_weighted_normal_component( &
      u_n, sigma, vn, vmin_normal=vn_floor, sigma_cutoff=default_velocity_sigma_cutoff, uniform=uniform(1:1) &
      )
    ! A negative barrier (accelerating sheath) raises the normal speed.
    vn_inf = sqrt(max(0.0_dp, vn(1)*vn(1) - barrier_normal_energy))
    v(axis_n) = inward_normal(axis_n)*vn_inf

    x = 0.0_dp
    ! One representable step inward avoids a zero-time boundary event without an untracked flight segment.
    x(axis_n) = nearest(boundary_value, inward_normal(axis_n))
    x(axis_t1) = pos_low(axis_t1) + (pos_high(axis_t1) - pos_low(axis_t1))*uniform(2)
    x(axis_t2) = pos_low(axis_t2) + (pos_high(axis_t2) - pos_low(axis_t2))*uniform(3)
  end subroutine sample_reservoir_face_particle_counter