!> シース数値モデルと app_config / 注入ランタイムの橋渡しを行うモジュール。 module bem_sheath_runtime use bem_kinds, only: dp, i32 use bem_constants, only: k_boltzmann, qe use bem_types, only: sim_config use bem_app_config_types, only: app_config, particle_species_spec use bem_string_utils, only: lower_ascii use bem_config_helpers, only: & resolve_inject_face, resolve_inward_normal, species_number_density_m3, species_temperature_k use bem_sheath_model_core, only: & sheath_model_species, zhao_params_type, zhao_local_state_type, & solve_no_photo_floating_potential, build_zhao_params, solve_zhao_unknowns, sample_zhao_state_at_z, & resolve_species_drift_speed, zhao_electron_vmin_normal, zhao_photo_vmin_normal, zhao_photo_emit_current_density implicit none type :: sheath_injection_context logical :: enabled = .false. logical :: has_photo_species = .false. logical :: has_local_reservoir_profile = .false. character(len=32) :: model = 'none' character(len=1) :: branch = ' ' integer(i32) :: electron_species = 0_i32 integer(i32) :: ion_species = 0_i32 integer(i32) :: photo_species = 0_i32 character(len=16) :: reference_face = '' integer(i32) :: reference_axis = 0_i32 real(dp) :: reference_coordinate = 0.0d0 real(dp) :: reference_inward_normal(3) = 0.0d0 real(dp) :: reservoir_plane_distance_m = 0.0d0 real(dp) :: reservoir_phi_v = 0.0d0 real(dp) :: phi0_v = 0.0d0 real(dp) :: phi_m_v = 0.0d0 real(dp) :: n_swe_inf_m3 = 0.0d0 real(dp) :: electron_number_density_m3 = 0.0d0 real(dp) :: electron_vmin_normal = 0.0d0 real(dp) :: ion_number_density_m3 = 0.0d0 real(dp) :: ion_normal_speed_mps = 0.0d0 real(dp) :: photo_emit_current_density_a_m2 = 0.0d0 real(dp) :: photo_vmin_normal = 0.0d0 end type sheath_injection_context public :: sheath_injection_context public :: resolve_sheath_injection_context contains subroutine resolve_sheath_injection_context(cfg, ctx) type(app_config), intent(in) :: cfg type(sheath_injection_context), intent(out) :: ctx integer(i32) :: electron_idx, ion_idx, photo_idx type(particle_species_spec) :: spec_e_cfg, spec_i_cfg, spec_p_cfg type(sheath_model_species) :: spec_e, spec_i, spec_p type(zhao_params_type) :: params type(zhao_local_state_type) :: local_state real(dp) :: phi0_v, phi_m_v, n_swe_inf_m3 real(dp) :: n_e_inf_m3, n_i_inf_m3, n_phe_ref_m3 real(dp) :: t_swe_ev, t_phe_ev, v_d_electron_mps, v_d_ion_mps real(dp) :: reference_coordinate, reference_inward_normal(3) integer :: reference_axis character(len=32) :: model character(len=1) :: branch ctx = sheath_injection_context() model = trim(lower_ascii(cfg%sim%sheath_injection_model)) ctx%model = model if (model == 'none') return call detect_sheath_species(cfg, electron_idx, ion_idx, photo_idx) spec_e_cfg = cfg%particle_species(electron_idx) spec_i_cfg = cfg%particle_species(ion_idx) spec_e = to_sheath_model_species(spec_e_cfg) spec_i = to_sheath_model_species(spec_i_cfg) ctx%enabled = .true. ctx%electron_species = electron_idx ctx%ion_species = ion_idx ctx%reference_face = trim(lower_ascii(spec_e_cfg%inject_face)) call resolve_sheath_reference_plane( & cfg%sim, ctx%reference_face, reference_axis, reference_coordinate, reference_inward_normal & ) ctx%reference_axis = int(reference_axis, i32) ctx%reference_coordinate = reference_coordinate ctx%reference_inward_normal = reference_inward_normal n_e_inf_m3 = spec_e%number_density_m3 n_i_inf_m3 = spec_i%number_density_m3 t_swe_ev = spec_e%temperature_k*k_boltzmann/qe v_d_electron_mps = resolve_species_drift_speed( & spec_e, trim(lower_ascii(cfg%sim%sheath_electron_drift_mode)), & reference_inward_normal) v_d_ion_mps = resolve_species_drift_speed( & spec_i, trim(lower_ascii(cfg%sim%sheath_ion_drift_mode)), & reference_inward_normal) if (v_d_electron_mps <= 0.0d0) error stop 'sheath electron drift must point inward.' if (v_d_ion_mps <= 0.0d0) error stop 'sheath ion drift must point inward.' select case (model) case ('floating_no_photo') call solve_no_photo_floating_potential(spec_e, spec_i, reference_inward_normal, phi0_v) ctx%branch = 'N' ctx%phi0_v = phi0_v ctx%phi_m_v = phi0_v ctx%n_swe_inf_m3 = n_e_inf_m3 ctx%electron_number_density_m3 = n_e_inf_m3 ctx%electron_vmin_normal = sqrt(max(0.0d0, 2.0d0*qe*max(0.0d0, -phi0_v))/abs(spec_e%m_particle)) return case ('zhao_auto', 'zhao_a', 'zhao_b', 'zhao_c') if (photo_idx <= 0_i32) then error stop 'Zhao sheath injection requires one enabled negative-q photo_raycast species.' end if spec_p_cfg = cfg%particle_species(photo_idx) spec_p = to_sheath_model_species(spec_p_cfg) if (abs(spec_p_cfg%q_particle) <= 0.0d0) error stop 'Zhao photoelectron species must have non-zero q_particle.' ctx%photo_species = photo_idx ctx%has_photo_species = .true. t_phe_ev = spec_p%temperature_k*k_boltzmann/qe n_phe_ref_m3 = cfg%sim%sheath_photoelectron_ref_density_cm3*1.0d6 call build_zhao_params( & cfg%sim%sheath_alpha_deg, n_i_inf_m3, n_phe_ref_m3, t_swe_ev, t_phe_ev, v_d_electron_mps, v_d_ion_mps, & spec_i%m_particle, abs(spec_e%m_particle), params & ) call solve_zhao_unknowns(model, params, phi0_v, phi_m_v, n_swe_inf_m3, branch) ctx%branch = branch ctx%phi0_v = phi0_v ctx%phi_m_v = phi_m_v ctx%n_swe_inf_m3 = n_swe_inf_m3 ctx%electron_number_density_m3 = n_swe_inf_m3 ctx%electron_vmin_normal = zhao_electron_vmin_normal(branch, phi0_v, phi_m_v, abs(spec_e%m_particle)) ctx%photo_vmin_normal = zhao_photo_vmin_normal(branch, phi0_v, phi_m_v, abs(spec_p%m_particle)) ctx%photo_emit_current_density_a_m2 = zhao_photo_emit_current_density(branch, phi0_v, phi_m_v, spec_p_cfg%q_particle, params) if (cfg%sim%has_sheath_reference_coordinate) then call sample_zhao_reservoir_state(cfg%sim, spec_e_cfg, ctx, params, branch, phi0_v, phi_m_v, n_swe_inf_m3, local_state) ctx%has_local_reservoir_profile = .true. ctx%reservoir_phi_v = local_state%phi_v ctx%electron_number_density_m3 = local_state%electron_source_density_m3 if (local_state%swe_reflected_active) then ctx%electron_vmin_normal = 0.0d0 else ctx%electron_vmin_normal = local_state%vcut_swe_mps end if ctx%ion_number_density_m3 = local_state%n_swi_m3 ctx%ion_normal_speed_mps = local_state%v_i_mps end if case default error stop 'Unknown sim.sheath_injection_model in runtime.' end select end subroutine resolve_sheath_injection_context subroutine detect_sheath_species(cfg, electron_idx, ion_idx, photo_idx) type(app_config), intent(in) :: cfg integer(i32), intent(out) :: electron_idx, ion_idx, photo_idx integer(i32) :: s character(len=16) :: mode electron_idx = 0_i32 ion_idx = 0_i32 photo_idx = 0_i32 do s = 1_i32, cfg%n_particle_species if (.not. cfg%particle_species(s)%enabled) cycle mode = trim(lower_ascii(cfg%particle_species(s)%source_mode)) select case (mode) case ('reservoir_face') if (cfg%particle_species(s)%q_particle < 0.0d0) then if (electron_idx == 0_i32) electron_idx = s else if (cfg%particle_species(s)%q_particle > 0.0d0) then if (ion_idx == 0_i32) ion_idx = s end if case ('photo_raycast') if (cfg%particle_species(s)%q_particle < 0.0d0 .and. photo_idx == 0_i32) photo_idx = s end select end do if (electron_idx <= 0_i32) error stop 'sheath injection requires one enabled negative-q reservoir_face species.' if (ion_idx <= 0_i32) error stop 'sheath injection requires one enabled positive-q reservoir_face species.' if (trim(lower_ascii(cfg%particle_species(electron_idx)%inject_face)) & /= trim(lower_ascii(cfg%particle_species(ion_idx)%inject_face))) then error stop 'sheath electron/ion reservoir species must share the same inject_face.' end if end subroutine detect_sheath_species subroutine sample_zhao_reservoir_state(sim, spec_e, ctx, p, branch, phi0_v, phi_m_v, n_swe_inf_m3, state) type(sim_config), intent(in) :: sim type(particle_species_spec), intent(in) :: spec_e type(sheath_injection_context), intent(inout) :: ctx type(zhao_params_type), intent(in) :: p character(len=1), intent(in) :: branch real(dp), intent(in) :: phi0_v, phi_m_v, n_swe_inf_m3 type(zhao_local_state_type), intent(out) :: state integer :: axis real(dp) :: boundary_value, distance_m, axis_sign call resolve_inject_face(sim%box_min, sim%box_max, spec_e%inject_face, axis, boundary_value) axis_sign = -ctx%reference_inward_normal(axis) distance_m = (boundary_value - ctx%reference_coordinate)*axis_sign if (distance_m < -1.0d-12) then error stop 'sheath_reference_coordinate must lie on or inside the shared reservoir_face boundary.' end if ctx%reservoir_plane_distance_m = max(0.0d0, distance_m) call sample_zhao_state_at_z(p, branch, phi0_v, phi_m_v, n_swe_inf_m3, ctx%reservoir_plane_distance_m, state) end subroutine sample_zhao_reservoir_state subroutine resolve_sheath_reference_plane(sim, inject_face, axis, reference_coordinate, inward_normal) type(sim_config), intent(in) :: sim character(len=*), intent(in) :: inject_face integer, intent(out) :: axis real(dp), intent(out) :: reference_coordinate real(dp), intent(out) :: inward_normal(3) real(dp) :: boundary_value call resolve_inward_normal(inject_face, inward_normal) call resolve_inject_face(sim%box_min, sim%box_max, inject_face, axis, boundary_value) if (sim%has_sheath_reference_coordinate) then reference_coordinate = sim%sheath_reference_coordinate else reference_coordinate = boundary_value end if end subroutine resolve_sheath_reference_plane type(sheath_model_species) function to_sheath_model_species(spec_cfg) result(spec) type(particle_species_spec), intent(in) :: spec_cfg spec%q_particle = spec_cfg%q_particle spec%m_particle = spec_cfg%m_particle if (spec_cfg%has_number_density_m3 .or. spec_cfg%has_number_density_cm3) then spec%number_density_m3 = species_number_density_m3(spec_cfg) end if spec%temperature_k = species_temperature_k(spec_cfg) spec%drift_velocity = spec_cfg%drift_velocity end function to_sheath_model_species end module bem_sheath_runtime