!> 1バッチ内で固定する静電場 snapshot。 module bem_electrostatic_snapshot use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use bem_kinds, only: dp, i32 use bem_constants, only: eps0, pi use bem_types, only: mesh_type, sim_config, particle_domain_config, bc_periodic use bem_mpi, only: mpi_context use bem_field_solver, only: field_solver_type use bem_physics_config_types, only: field_physics_config, periodic2_physics_config, panel_kernel_config, & validate_upper_fourier_height, upper_fourier_auto_height_range, & physics_config_ok use bem_string_utils, only: lower_ascii use bem_runtime_modes, only: lower_boundary_model_id, lower_boundary_e_bottom_zero, & lower_boundary_symmetric_vacuum, lower_boundary_e_top_from_sheath, lower_boundary_unknown use bem_periodic_zero_mode_plan, only: periodic_zero_mode_plan_type, periodic_zero_mode_state_type, & build_periodic_zero_mode_plan, refresh_periodic_zero_mode_state, & symmetric_vacuum_bottom_field, periodic_zero_mode_ok use bem_periodic_zero_mode_eval, only: eval_periodic_zero_mode, zero_mode_trace_plus use bem_coulomb_fmm_periodic_nonzero_reference, only: eval_periodic_nonzero_panel_reference use bem_coulomb_fmm_periodic_nonzero_upper_vacuum, only: & periodic_nonzero_upper_vacuum_plan_type use bem_panel_geometry, only: panel_geometry_type, init_panel_geometry, panel_geometry_ok use bem_panel_kernel, only: panel_potential_field, panel_side_principal_value implicit none private type, public :: electrostatic_diagnostics_type logical :: upper_fourier_switch_active = .false. real(dp) :: upper_fourier_switch_height_m = 0.0_dp integer(i32) :: upper_fourier_switch_mode_layers = 0_i32 real(dp) :: upper_fourier_switch_selected_mismatch_v = 0.0_dp real(dp) :: lower_boundary_field_v_m = 0.0_dp logical :: split_periodic_active = .false. logical :: top_reference_available = .false. integer(i32) :: top_reference_last_batch = -1_i32 real(dp) :: top_reference_simulated_time = 0.0_dp real(dp) :: top_reference_z_high = 0.0_dp integer(i32) :: top_reference_sample_n = 0_i32 real(dp) :: top_reference_potential_mean = 0.0_dp real(dp) :: top_reference_potential_std = 0.0_dp real(dp) :: top_reference_potential_min = 0.0_dp real(dp) :: top_reference_potential_max = 0.0_dp real(dp) :: gauss_residual = 0.0_dp character(len=64) :: status = 'not_applicable' logical :: periodic_cache_hit = .false. integer(i32) :: periodic_operator_build_count = 0_i32 character(len=128) :: periodic_cache_fingerprint = '' character(len=512) :: periodic_cache_path = '' end type electrostatic_diagnostics_type type, public :: electrostatic_snapshot_type type(field_solver_type) :: nonzero_solver real(dp) :: prescribed_e(3) = [0.0_dp, 0.0_dp, 0.0_dp] real(dp) :: prescribed_phi_origin(3) = [0.0_dp, 0.0_dp, 0.0_dp] logical :: use_zero_mode = .false. logical :: use_panel_spectral_reference = .false. logical :: use_cached_kneq0 = .false. logical :: use_upper_panel_fourier_retry = .false. real(dp) :: periodic_length(2) = [0.0_dp, 0.0_dp] real(dp) :: periodic_origin(2) = [0.0_dp, 0.0_dp] integer(i32) :: reference_mode_layers = 0_i32 integer(i32) :: panel_quadrature_order = 0_i32 type(periodic_zero_mode_plan_type) :: zero_plan type(periodic_zero_mode_state_type) :: zero_state type(periodic_nonzero_upper_vacuum_plan_type) :: upper_fourier_plan type(periodic2_physics_config) :: periodic_options !> periodic_options%lower_boundary_model を init で解決した値。zero mode の更新ごとに文字列を比べない。 integer(i32) :: lower_boundary_id = lower_boundary_unknown type(electrostatic_diagnostics_type) :: diagnostics real(dp) :: gauss_residual = 0.0_dp logical :: matching_plane_gauge_active = .false. real(dp) :: matching_plane_z = 0.0_dp real(dp) :: matching_plane_phi = 0.0_dp real(dp) :: matching_plane_top_field = 0.0_dp contains procedure :: init => init_electrostatic_snapshot procedure :: refresh => refresh_electrostatic_snapshot procedure :: refresh_zero_mode => refresh_snapshot_zero_mode procedure :: eval_local_e => eval_snapshot_local_e procedure :: eval_upper_panel_fourier_e => eval_snapshot_upper_panel_fourier_e procedure :: eval_upper_panel_fourier_phi => eval_snapshot_upper_panel_fourier_phi procedure :: eval_local_phi => eval_snapshot_local_phi procedure :: eval_local_phi_without_primary_self => eval_snapshot_local_phi_without_primary_self procedure :: compute_mesh_potential => compute_snapshot_mesh_potential procedure :: measure_kneq0_potential_step => measure_snapshot_kneq0_potential_step procedure :: get_diagnostics => get_snapshot_diagnostics procedure :: set_matching_plane_gauge => set_snapshot_matching_plane_gauge procedure :: set_matching_plane_top_field => set_snapshot_matching_plane_top_field procedure :: measure_switch_mismatch => measure_snapshot_switch_mismatch end type electrostatic_snapshot_type interface module subroutine init_snapshot_upper_fourier_switch(self, mesh, sim) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim end subroutine init_snapshot_upper_fourier_switch module subroutine measure_snapshot_switch_mismatch(self, mesh, sim, phi_rms, ez_rms) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim real(dp), intent(out) :: phi_rms, ez_rms end subroutine measure_snapshot_switch_mismatch module subroutine eval_snapshot_local_e(self, mesh, position, electric_field) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: position(3) real(dp), intent(out) :: electric_field(3) end subroutine eval_snapshot_local_e module subroutine eval_snapshot_upper_panel_fourier_e(self, mesh, position, electric_field, available) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: position(3) real(dp), intent(out) :: electric_field(3) logical, intent(out) :: available end subroutine eval_snapshot_upper_panel_fourier_e module subroutine eval_snapshot_upper_panel_fourier_phi(self, mesh, position, potential, available) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: position(3) real(dp), intent(out) :: potential logical, intent(out) :: available end subroutine eval_snapshot_upper_panel_fourier_phi module subroutine eval_snapshot_local_phi(self, mesh, sim, position, potential) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh class(particle_domain_config), intent(in) :: sim real(dp), intent(in) :: position(3) real(dp), intent(out) :: potential end subroutine eval_snapshot_local_phi module subroutine eval_snapshot_local_phi_without_primary_self(self, mesh, sim, element, potential) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim integer(i32), intent(in) :: element real(dp), intent(out) :: potential end subroutine eval_snapshot_local_phi_without_primary_self module subroutine compute_snapshot_mesh_potential(self, mesh, sim, potential) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim real(dp), intent(out) :: potential(:) end subroutine compute_snapshot_mesh_potential module subroutine measure_snapshot_kneq0_potential_step( & self, mesh, candidate_charge, max_abs_delta_phi_v, delta_phi_v & ) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: candidate_charge(:) real(dp), intent(out) :: max_abs_delta_phi_v real(dp), intent(out), optional :: delta_phi_v(:) end subroutine measure_snapshot_kneq0_potential_step end interface contains !> mpi は cached_kneq0 の root operator を rank で分担して作るときの文脈。全 rank が同じ引数で呼ぶ。 subroutine init_electrostatic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config, mpi) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim type(field_physics_config), intent(in), optional :: field_config type(periodic2_physics_config), intent(in), optional :: periodic_config type(panel_kernel_config), intent(in), optional :: panel_config type(mpi_context), intent(in), optional :: mpi self%prescribed_e = sim%e0 self%prescribed_phi_origin = 0.0_dp self%use_zero_mode = .false. self%use_panel_spectral_reference = .false. self%use_cached_kneq0 = .false. self%use_upper_panel_fourier_retry = .false. self%gauss_residual = 0.0_dp self%matching_plane_gauge_active = .false. self%matching_plane_z = 0.0_dp self%matching_plane_phi = 0.0_dp self%matching_plane_top_field = 0.0_dp self%diagnostics = electrostatic_diagnostics_type() if (present(periodic_config)) then if (trim(lower_ascii(periodic_config%nonzero_mode_backend)) == 'panel_spectral_reference') then call init_panel_spectral_snapshot(self, mesh, sim, periodic_config) return end if if (trim(lower_ascii(periodic_config%nonzero_mode_backend)) == 'cached_kneq0') then if (.not. present(field_config) .or. .not. present(panel_config)) then error stop 'cached_kneq0 requires typed field and panel configuration.' end if call init_cached_periodic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config, mpi) return end if end if if (present(field_config) .and. present(periodic_config) .and. present(panel_config)) then call self%nonzero_solver%init(mesh, sim, field_config, periodic_config, panel_config, mpi=mpi) else call self%nonzero_solver%init(mesh, sim) end if end subroutine init_electrostatic_snapshot !> 電荷から場の状態を更新する。`mpi` を与えると FMM の M2L を rank で分担する（結果は同じ）。 subroutine refresh_electrostatic_snapshot(self, mesh, mpi) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(mpi_context), intent(in), optional :: mpi if (.not. self%use_panel_spectral_reference) call self%nonzero_solver%refresh(mesh, mpi) if (self%upper_fourier_plan%nmode > 0_i32) call self%upper_fourier_plan%refresh(mesh%q_elem) call self%refresh_zero_mode(mesh) end subroutine refresh_electrostatic_snapshot !> Refresh only the physical mean field after a charge commit; no FMM update is needed for receipts. subroutine refresh_snapshot_zero_mode(self, mesh) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh if (self%use_zero_mode) then if (self%matching_plane_gauge_active) then call refresh_periodic_zero_mode_state( & self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), & self%matching_plane_z, self%matching_plane_phi, self%zero_state & ) else call refresh_periodic_zero_mode_state( & self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), & self%zero_plan%break_z(1), 0.0_dp, self%zero_state & ) end if end if end subroutine refresh_snapshot_zero_mode !> 周期k=0成分の電位ゲージをmatching planeへ移す。 !! 非ゼロモードのsolver/cacheは変更しない。次回refreshでもこのゲージを保持する。 subroutine set_snapshot_matching_plane_gauge(self, mesh, z_high, phi_high) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: z_high, phi_high if (.not. self%use_zero_mode) then error stop 'matching-plane gauge requires an active periodic zero mode.' end if if (.not. ieee_is_finite(z_high) .or. .not. ieee_is_finite(phi_high)) then error stop 'matching-plane gauge values must be finite.' end if if (z_high <= self%zero_plan%break_z(self%zero_plan%nbreak)) then error stop 'matching-plane gauge must lie strictly above the charged geometry.' end if self%matching_plane_gauge_active = .true. self%matching_plane_z = z_high self%matching_plane_phi = phi_high call refresh_periodic_zero_mode_state( & self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), & z_high, phi_high, self%zero_state & ) end subroutine set_snapshot_matching_plane_gauge !> Set E_z(H+) from the outer root. The selected zero-mode boundary model consumes it on refresh. subroutine set_snapshot_matching_plane_top_field(self, mesh, field_v_m) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh real(dp), intent(in) :: field_v_m if (.not. ieee_is_finite(field_v_m)) error stop 'matching-plane top field must be finite.' self%matching_plane_top_field = field_v_m call self%refresh_zero_mode(mesh) end subroutine set_snapshot_matching_plane_top_field subroutine get_snapshot_diagnostics(self, diagnostics) class(electrostatic_snapshot_type), intent(in) :: self type(electrostatic_diagnostics_type), intent(out) :: diagnostics diagnostics = self%diagnostics diagnostics%gauss_residual = self%gauss_residual if (self%use_zero_mode) diagnostics%lower_boundary_field_v_m = self%zero_state%e_bottom end subroutine get_snapshot_diagnostics subroutine init_panel_spectral_snapshot(self, mesh, sim, periodic_config) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim type(periodic2_physics_config), intent(in) :: periodic_config integer(i32) :: status character(len=256) :: message if (trim(lower_ascii(periodic_config%zero_mode_policy)) /= 'exclude_k0' .or. & .not. supported_lower_boundary(periodic_config%lower_boundary_model)) then error stop 'panel_spectral_reference requires exclude_k0 and a supported lower boundary model.' end if if (.not. sim%use_box .or. any(sim%bc_low(1:2) /= bc_periodic) .or. & any(sim%bc_high(1:2) /= bc_periodic) .or. sim%bc_low(3) == bc_periodic .or. & sim%bc_high(3) == bc_periodic) then error stop 'panel_spectral_reference requires x/y periodic and z open.' end if self%periodic_length = sim%box_max(1:2) - sim%box_min(1:2) self%periodic_origin = sim%box_min(1:2) if (any(self%periodic_length <= 0.0_dp)) then error stop 'panel_spectral_reference requires positive periods.' end if call build_periodic_zero_mode_plan( & mesh, product(self%periodic_length), self%zero_plan, status, message & ) if (status /= periodic_zero_mode_ok) then error stop 'panel_spectral_reference zero-mode plan: '//trim(message) end if self%periodic_options = periodic_config self%lower_boundary_id = lower_boundary_model_id(periodic_config%lower_boundary_model) call refresh_periodic_zero_mode_state( & self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), & self%zero_plan%break_z(1), 0.0_dp, self%zero_state & ) self%reference_mode_layers = periodic_config%reference_mode_layers self%panel_quadrature_order = periodic_config%panel_quadrature_order self%use_panel_spectral_reference = .true. self%use_zero_mode = .true. self%diagnostics%split_periodic_active = .true. self%diagnostics%status = 'panel_spectral_reference' end subroutine init_panel_spectral_snapshot subroutine init_cached_periodic_snapshot(self, mesh, sim, field_config, periodic_config, panel_config, mpi) class(electrostatic_snapshot_type), intent(inout) :: self type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim type(field_physics_config), intent(in) :: field_config type(periodic2_physics_config), intent(in) :: periodic_config type(panel_kernel_config), intent(in) :: panel_config type(mpi_context), intent(in), optional :: mpi integer(i32) :: status character(len=256) :: message if (trim(lower_ascii(periodic_config%zero_mode_policy)) /= 'exclude_k0' .or. & .not. supported_lower_boundary(periodic_config%lower_boundary_model)) then error stop 'cached_kneq0 requires exclude_k0 and a supported lower boundary model.' end if if (.not. sim%use_box .or. any(sim%bc_low(1:2) /= bc_periodic) .or. & any(sim%bc_high(1:2) /= bc_periodic) .or. sim%bc_low(3) == bc_periodic .or. & sim%bc_high(3) == bc_periodic) then error stop 'cached_kneq0 snapshot currently requires x/y periodic and z open.' end if call validate_upper_fourier_height(sim, periodic_config, maxval(mesh%bb_max(3, :)), status, message) if (status /= physics_config_ok) error stop trim(message) call build_periodic_zero_mode_plan( & mesh, product(sim%box_max(1:2) - sim%box_min(1:2)), self%zero_plan, status, message & ) if (status /= periodic_zero_mode_ok) error stop 'cached_kneq0 zero-mode plan: '//trim(message) self%periodic_options = periodic_config self%lower_boundary_id = lower_boundary_model_id(periodic_config%lower_boundary_model) call refresh_periodic_zero_mode_state( & self%zero_plan, mesh%q_elem, zero_mode_bottom_field(self, mesh%q_elem), & self%zero_plan%break_z(1), 0.0_dp, self%zero_state & ) call self%nonzero_solver%init(mesh, sim, field_config, periodic_config, panel_config, mpi=mpi) self%use_cached_kneq0 = .true. self%use_zero_mode = .true. self%periodic_length = sim%box_max(1:2) - sim%box_min(1:2) self%periodic_origin = sim%box_min(1:2) self%reference_mode_layers = periodic_config%reference_mode_layers self%panel_quadrature_order = periodic_config%panel_quadrature_order self%diagnostics%split_periodic_active = .true. self%diagnostics%status = 'cached_kneq0' self%diagnostics%periodic_cache_hit = self%nonzero_solver%fmm_core_plan%periodic_cache_hit self%diagnostics%periodic_operator_build_count = & self%nonzero_solver%fmm_core_plan%periodic_operator_build_count self%diagnostics%periodic_cache_fingerprint = & self%nonzero_solver%fmm_core_plan%periodic_cache_fingerprint self%diagnostics%periodic_cache_path = self%nonzero_solver%fmm_core_plan%periodic_cache_path if (periodic_config%upper_fourier_switch /= 'off') then call init_snapshot_upper_fourier_switch(self, mesh, sim) else if (trim(lower_ascii(sim%multiple_box_events_retry_backend)) == 'upper_panel_fourier') then call self%upper_fourier_plan%build( & mesh, self%periodic_length(1), self%periodic_length(2), self%reference_mode_layers & ) call self%upper_fourier_plan%refresh(mesh%q_elem) end if self%use_upper_panel_fourier_retry = trim(lower_ascii(sim%multiple_box_events_retry_backend)) == 'upper_panel_fourier' end subroutine init_cached_periodic_snapshot pure real(dp) function zero_mode_bottom_field(self, charge) result(field) class(electrostatic_snapshot_type), intent(in) :: self real(dp), intent(in) :: charge(:) select case (self%lower_boundary_id) case (lower_boundary_e_bottom_zero) field = 0.0_dp case (lower_boundary_e_top_from_sheath) field = self%matching_plane_top_field - sum(charge)/(eps0*self%zero_plan%area_xy) case (lower_boundary_symmetric_vacuum) field = symmetric_vacuum_bottom_field(self%zero_plan, charge) case default error stop 'unsupported periodic2 lower boundary model.' end select end function zero_mode_bottom_field pure logical function supported_lower_boundary(model) result(supported) character(len=*), intent(in) :: model select case (trim(lower_ascii(model))) case ('e_bottom_zero', 'symmetric_vacuum', 'e_top_from_sheath') supported = .true. case default supported = .false. end select end function supported_lower_boundary end module bem_electrostatic_snapshot