!> Matching-plane outer-sheath response table の immutable snapshot と補間。 module bem_matching_plane_response use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use bem_kinds, only: dp, i32, i64 use bem_mpi, only: mpi_context use bem_matching_plane_contract, only: & matching_plane_response_input_count, & matching_plane_response_output_count, & matching_plane_input_displacement, & matching_plane_input_photoelectron_outward_flux, & matching_plane_input_photoelectron_mean_normal_energy, & matching_plane_input_electron_outward_flux, & matching_plane_input_ion_outward_flux, & matching_plane_output_matching_potential, & matching_plane_output_electron_inward_flux, & matching_plane_output_ion_inward_flux, & matching_plane_output_electron_access_potential, & matching_plane_output_ion_access_potential, & matching_plane_output_photoelectron_barrier_potential implicit none private ! The table and analytic model share the same response coordinates. public :: matching_plane_response_input_count public :: matching_plane_response_output_count public :: matching_plane_input_displacement public :: matching_plane_input_photoelectron_outward_flux public :: matching_plane_input_photoelectron_mean_normal_energy public :: matching_plane_input_electron_outward_flux public :: matching_plane_input_ion_outward_flux public :: matching_plane_output_matching_potential public :: matching_plane_output_electron_inward_flux public :: matching_plane_output_ion_inward_flux public :: matching_plane_output_electron_access_potential public :: matching_plane_output_ion_access_potential public :: matching_plane_output_photoelectron_barrier_potential integer(i32), parameter, public :: matching_plane_response_ok = 0_i32 integer(i32), parameter, public :: matching_plane_response_invalid_argument = 1_i32 integer(i32), parameter, public :: matching_plane_response_io_error = 2_i32 integer(i32), parameter, public :: matching_plane_response_invalid_header = 3_i32 integer(i32), parameter, public :: matching_plane_response_invalid_metadata = 4_i32 integer(i32), parameter, public :: matching_plane_response_invalid_row = 5_i32 integer(i32), parameter, public :: matching_plane_response_invalid_grid = 6_i32 integer(i32), parameter, public :: matching_plane_response_out_of_range = 7_i32 character(len=*), parameter, public :: matching_plane_response_query_csv_header = & 'displacement_c_m2,'// & 'photoelectron_outward_number_flux_m2_s,'// & 'photoelectron_outward_mean_normal_energy_ev,'// & 'electron_outward_number_flux_m2_s,'// & 'ion_outward_number_flux_m2_s' character(len=*), parameter, public :: matching_plane_response_csv_header = & matching_plane_response_query_csv_header//','// & 'matching_potential_v,'// & 'electron_inward_number_flux_m2_s,'// & 'ion_inward_number_flux_m2_s,'// & 'electron_access_potential_v,'// & 'ion_access_potential_v,'// & 'photoelectron_barrier_potential_v' character(len=52), parameter :: input_names(matching_plane_response_input_count) = & [character(len=52) :: & 'displacement_c_m2', & 'photoelectron_outward_number_flux_m2_s', & 'photoelectron_outward_mean_normal_energy_ev', & 'electron_outward_number_flux_m2_s', & 'ion_outward_number_flux_m2_s'] !> 検証済み response table。components は外部から変更できない。 type, public :: matching_plane_response_table_type private logical :: loaded = .false. character(len=:), allocatable :: source_path real(dp) :: matching_plane_z_m = 0.0_dp integer(i32) :: axis_sizes(matching_plane_response_input_count) = 0_i32 real(dp), allocatable :: axes(:, :) real(dp), allocatable :: response_values(:, :) contains procedure, public :: evaluate => evaluate_matching_plane_response procedure, public :: get_axis_data => get_matching_plane_axis_data procedure, public :: get_matching_plane_z => get_matching_plane_z procedure, public :: get_source_path => get_matching_plane_source_path procedure, public :: is_loaded => matching_plane_response_is_loaded end type matching_plane_response_table_type type :: matching_plane_response_snapshot_entry character(len=:), allocatable :: path type(matching_plane_response_table_type) :: table end type matching_plane_response_snapshot_entry type(matching_plane_response_snapshot_entry), allocatable, save :: response_snapshots(:) public :: get_matching_plane_response_snapshot public :: load_matching_plane_response_mpi public :: reset_matching_plane_response_snapshot_cache interface !> root が読み込んだテーブルを全 rank に配信する。 module subroutine load_matching_plane_response_mpi(path, mpi, table, status, message) character(len=*), intent(in) :: path type(mpi_context), intent(in) :: mpi type(matching_plane_response_table_type), intent(out) :: table integer(i32), intent(out) :: status character(len=*), intent(out) :: message end subroutine load_matching_plane_response_mpi module subroutine read_matching_plane_response_csv(path, table, status, message) character(len=*), intent(in) :: path type(matching_plane_response_table_type), intent(out) :: table integer(i32), intent(out) :: status character(len=*), intent(out) :: message end subroutine read_matching_plane_response_csv pure module subroutine accept(status, message) integer(i32), intent(out) :: status character(len=*), intent(out) :: message end subroutine accept pure module subroutine reject(code, text, status, message) integer(i32), intent(in) :: code character(len=*), intent(in) :: text integer(i32), intent(out) :: status character(len=*), intent(out) :: message end subroutine reject end interface contains !> path ごとに最初の検証済み内容を固定し、同一 process 内で共有する。 subroutine get_matching_plane_response_snapshot(path, table, status, message) character(len=*), intent(in) :: path type(matching_plane_response_table_type), intent(out) :: table integer(i32), intent(out) :: status character(len=*), intent(out) :: message integer :: entry table = matching_plane_response_table_type() call ensure_matching_plane_response_snapshot(path, entry, status, message) if (status /= matching_plane_response_ok) return table = response_snapshots(entry)%table end subroutine get_matching_plane_response_snapshot subroutine ensure_matching_plane_response_snapshot(path, entry_index, status, message) character(len=*), intent(in) :: path integer, intent(out) :: entry_index integer(i32), intent(out) :: status character(len=*), intent(out) :: message type(matching_plane_response_snapshot_entry), allocatable :: grown(:) type(matching_plane_response_table_type) :: loaded_table character(len=:), allocatable :: normalized_path integer :: entry, old_size entry_index = 0 call accept(status, message) normalized_path = trim(path) if (len(normalized_path) == 0) then call reject( & matching_plane_response_invalid_argument, 'matching-plane response path must not be empty.', status, message & ) return end if if (allocated(response_snapshots)) then do entry = 1, size(response_snapshots) if (response_snapshots(entry)%path == normalized_path) then entry_index = entry return end if end do old_size = size(response_snapshots) else old_size = 0 end if call read_matching_plane_response_csv(normalized_path, loaded_table, status, message) if (status /= matching_plane_response_ok) return allocate (grown(old_size + 1)) do entry = 1, old_size call move_alloc(response_snapshots(entry)%path, grown(entry)%path) call move_matching_plane_response_table(response_snapshots(entry)%table, grown(entry)%table) end do grown(old_size + 1)%path = normalized_path call move_matching_plane_response_table(loaded_table, grown(old_size + 1)%table) call move_alloc(grown, response_snapshots) entry_index = old_size + 1 end subroutine ensure_matching_plane_response_snapshot !> allocatable componentsをcopyせずcache entryへ移す。 subroutine move_matching_plane_response_table(source, destination) type(matching_plane_response_table_type), intent(inout) :: source type(matching_plane_response_table_type), intent(out) :: destination destination = matching_plane_response_table_type() destination%loaded = source%loaded destination%matching_plane_z_m = source%matching_plane_z_m destination%axis_sizes = source%axis_sizes if (allocated(source%source_path)) call move_alloc(source%source_path, destination%source_path) if (allocated(source%axes)) call move_alloc(source%axes, destination%axes) if (allocated(source%response_values)) call move_alloc(source%response_values, destination%response_values) source = matching_plane_response_table_type() end subroutine move_matching_plane_response_table !> 独立 run や unit test の開始時に immutable snapshot cache を解放する。 subroutine reset_matching_plane_response_snapshot_cache() if (allocated(response_snapshots)) deallocate (response_snapshots) end subroutine reset_matching_plane_response_snapshot_cache !> 5入力から6出力を多重線形補間する。singleton axis は inactive dimension として扱う。 subroutine evaluate_matching_plane_response(self, input, output, status, message) class(matching_plane_response_table_type), intent(in) :: self real(dp), intent(in) :: input(matching_plane_response_input_count) real(dp), intent(out) :: output(matching_plane_response_output_count) integer(i32), intent(out) :: status character(len=*), intent(out) :: message integer(i32) :: lower(matching_plane_response_input_count) integer(i32) :: upper(matching_plane_response_input_count) integer(i32) :: index_at_corner(matching_plane_response_input_count) real(dp) :: upper_weight(matching_plane_response_input_count) real(dp) :: query_value, range_tolerance real(dp) :: corner_weight integer(i64) :: linear_index, stride integer :: axis, corner output = 0.0_dp call accept(status, message) if (.not. self%loaded) then call reject( & matching_plane_response_invalid_argument, 'matching-plane response table is not loaded.', status, message & ) return end if if (any(.not. ieee_is_finite(input))) then call reject( & matching_plane_response_invalid_argument, 'matching-plane response query must be finite.', status, message & ) return end if do axis = 1, matching_plane_response_input_count if (self%axis_sizes(axis) == 1_i32) then ! A singleton explicitly removes this dimension from the response. Its ! coordinate documents the generation point but does not restrict queries. lower(axis) = 1_i32 upper(axis) = 1_i32 upper_weight(axis) = 0.0_dp cycle end if query_value = input(axis) range_tolerance = 64.0_dp*epsilon(1.0_dp)*max( & 1.0_dp, abs(self%axes(1, axis)), & abs(self%axes(self%axis_sizes(axis), axis)) & ) if (query_value < self%axes(1, axis) - range_tolerance .or. & query_value > self%axes(self%axis_sizes(axis), axis) + range_tolerance) then call reject_out_of_range( & axis, query_value, self%axes(1, axis), self%axes(self%axis_sizes(axis), axis), status, message & ) return end if ! 丸め幅だけ端点を越えた値は端点そのものとして補間する。この幅を ! 越える外挿は上で拒否し、物理的な範囲外入力をclampしない。 if (query_value < self%axes(1, axis)) query_value = self%axes(1, axis) if (query_value > self%axes(self%axis_sizes(axis), axis)) then query_value = self%axes(self%axis_sizes(axis), axis) end if call locate_axis_interval( & self%axes(:self%axis_sizes(axis), axis), query_value, lower(axis), upper(axis), upper_weight(axis) & ) end do do corner = 0, 2**matching_plane_response_input_count - 1 corner_weight = 1.0_dp do axis = 1, matching_plane_response_input_count if (btest(corner, axis - 1)) then if (self%axis_sizes(axis) == 1_i32) then corner_weight = 0.0_dp exit end if index_at_corner(axis) = upper(axis) corner_weight = corner_weight*upper_weight(axis) else index_at_corner(axis) = lower(axis) corner_weight = corner_weight*(1.0_dp - upper_weight(axis)) end if end do if (corner_weight == 0.0_dp) cycle linear_index = 1_i64 stride = 1_i64 do axis = 1, matching_plane_response_input_count linear_index = linear_index + int(index_at_corner(axis) - 1_i32, i64)*stride stride = stride*int(self%axis_sizes(axis), i64) end do output = output + corner_weight*self%response_values(:, int(linear_index)) end do if (any(.not. ieee_is_finite(output))) then output = 0.0_dp call reject( & matching_plane_response_invalid_grid, & 'matching-plane response interpolation produced a non-finite output.', status, message & ) end if end subroutine evaluate_matching_plane_response !> フィードバックの範囲と尺度を決めるために補間軸を返す。 subroutine get_matching_plane_axis_data( & self, axis_sizes, axis_values, matching_plane_z_m, status, message & ) class(matching_plane_response_table_type), intent(in) :: self integer(i32), allocatable, intent(out) :: axis_sizes(:) real(dp), allocatable, intent(out) :: axis_values(:) real(dp), intent(out) :: matching_plane_z_m integer(i32), intent(out), optional :: status character(len=*), intent(out), optional :: message integer :: axis, first, last, total_axis_values matching_plane_z_m = 0.0_dp if (.not. self%loaded) then allocate (axis_sizes(0), axis_values(0)) call assign_optional_status( & matching_plane_response_invalid_argument, 'matching-plane response table is not loaded.', status, message & ) return end if allocate (axis_sizes(matching_plane_response_input_count), source=self%axis_sizes) total_axis_values = sum(self%axis_sizes) allocate (axis_values(total_axis_values)) first = 1 do axis = 1, matching_plane_response_input_count last = first + self%axis_sizes(axis) - 1 axis_values(first:last) = self%axes(:self%axis_sizes(axis), axis) first = last + 1 end do matching_plane_z_m = self%matching_plane_z_m call assign_optional_status(matching_plane_response_ok, '', status, message) end subroutine get_matching_plane_axis_data !> CSV metadata のmatching-plane高さを返す。 subroutine get_matching_plane_z(self, matching_plane_z_m, status, message) class(matching_plane_response_table_type), intent(in) :: self real(dp), intent(out) :: matching_plane_z_m integer(i32), intent(out), optional :: status character(len=*), intent(out), optional :: message matching_plane_z_m = 0.0_dp if (.not. self%loaded) then call assign_optional_status( & matching_plane_response_invalid_argument, 'matching-plane response table is not loaded.', status, message & ) return end if matching_plane_z_m = self%matching_plane_z_m call assign_optional_status(matching_plane_response_ok, '', status, message) end subroutine get_matching_plane_z !> snapshot のcache keyとなったpathを返す。 subroutine get_matching_plane_source_path(self, path, status, message) class(matching_plane_response_table_type), intent(in) :: self character(len=:), allocatable, intent(out) :: path integer(i32), intent(out), optional :: status character(len=*), intent(out), optional :: message if (.not. self%loaded) then path = '' call assign_optional_status( & matching_plane_response_invalid_argument, 'matching-plane response table is not loaded.', status, message & ) return end if path = self%source_path call assign_optional_status(matching_plane_response_ok, '', status, message) end subroutine get_matching_plane_source_path pure logical function matching_plane_response_is_loaded(self) result(loaded) class(matching_plane_response_table_type), intent(in) :: self loaded = self%loaded end function matching_plane_response_is_loaded pure subroutine locate_axis_interval(axis_values, value, lower, upper, upper_weight) real(dp), intent(in) :: axis_values(:), value integer(i32), intent(out) :: lower, upper real(dp), intent(out) :: upper_weight integer :: low_index, high_index, middle if (value == axis_values(size(axis_values))) then lower = int(size(axis_values) - 1, i32) upper = int(size(axis_values), i32) upper_weight = 1.0_dp return end if low_index = 1 high_index = size(axis_values) - 1 do while (low_index <= high_index) middle = (low_index + high_index)/2 if (value < axis_values(middle)) then high_index = middle - 1 else if (value >= axis_values(middle + 1)) then low_index = middle + 1 else lower = int(middle, i32) upper = int(middle + 1, i32) upper_weight = (value - axis_values(middle))/(axis_values(middle + 1) - axis_values(middle)) return end if end do ! evaluate側のclosed-range検証後なので、ここへ到達するのは丸め境界だけ。 lower = 1_i32 upper = 2_i32 upper_weight = 0.0_dp end subroutine locate_axis_interval subroutine reject_out_of_range(axis, value, lower, upper, status, message) integer, intent(in) :: axis real(dp), intent(in) :: value, lower, upper integer(i32), intent(out) :: status character(len=*), intent(out) :: message character(len=512) :: text write (text, '(a,a,a,es16.8,a,es16.8,a,es16.8,a)') & 'matching-plane query ', trim(input_names(axis)), '=', value, ' is outside [', lower, ',', upper, '].' call reject(matching_plane_response_out_of_range, trim(text), status, message) end subroutine reject_out_of_range pure subroutine assign_optional_status(code, text, status, message) integer(i32), intent(in) :: code character(len=*), intent(in) :: text integer(i32), intent(out), optional :: status character(len=*), intent(out), optional :: message if (present(status)) status = code if (present(message)) message = text end subroutine assign_optional_status end module bem_matching_plane_response