!> Online Zhao evaluator から既存形式の matching-plane 応答表を生成する。 module bem_matching_plane_response_generator use bem_kinds, only: dp, i32, i64 use bem_app_config_types, only: app_config use bem_matching_plane_response, only: matching_plane_response_input_count, & matching_plane_response_output_count, & matching_plane_response_query_csv_header, & matching_plane_response_csv_header use bem_matching_plane_response_provider, only: matching_plane_response_provider_type, & matching_plane_provider_ok use bem_filesystem, only: atomic_rename, filesystem_success use bem_mpi, only: mpi_context use bem_string_utils, only: lower_ascii use bem_matching_plane_query_io, only: read_matching_plane_query_csv, & matching_plane_query_ok, matching_plane_query_io_error, & matching_plane_query_invalid_grid implicit none private integer(i32), parameter, public :: matching_plane_generator_ok = matching_plane_query_ok integer(i32), parameter, public :: matching_plane_generator_invalid_argument = 1_i32 integer(i32), parameter, public :: matching_plane_generator_io_error = matching_plane_query_io_error integer(i32), parameter, public :: matching_plane_generator_invalid_grid = matching_plane_query_invalid_grid integer(i32), parameter, public :: matching_plane_generator_evaluation_failure = 4_i32 public :: generate_matching_plane_zhao_response_table contains !> 5列query gridを評価し、既存table backend互換の11列CSVを書く。 subroutine generate_matching_plane_zhao_response_table(cfg, query_path, output_path, status, message) type(app_config), intent(in) :: cfg character(len=*), intent(in) :: query_path, output_path integer(i32), intent(out) :: status character(len=*), intent(out) :: message type(matching_plane_response_provider_type) :: provider type(mpi_context) :: serial_mpi real(dp), allocatable :: queries(:, :), responses(:, :) integer(i32) :: provider_status integer :: row, output_unit, ios, close_ios, rename_status character(len=512) :: provider_message character(len=:), allocatable :: temporary_output_path call accept_generator(status, message) if (trim(lower_ascii(cfg%surface_current%model)) /= 'matching_plane_quasistatic' .or. & trim(lower_ascii(cfg%surface_current%response_backend)) /= 'zhao_online') then call reject_generator( & matching_plane_generator_invalid_argument, & 'response generation requires matching_plane_quasistatic with response_backend="zhao_online".', & status, message & ) return end if if (trim(lower_ascii(cfg%surface_current%zhao_root_selection)) == 'continuation') then call reject_generator( & matching_plane_generator_invalid_argument, & 'response-table generation does not define an accepted endpoint for Zhao continuation.', & status, message & ) return end if if (len_trim(query_path) == 0 .or. len_trim(output_path) == 0) then call reject_generator( & matching_plane_generator_invalid_argument, & 'query and output paths must not be empty.', status, message & ) return end if if (trim(query_path) == trim(output_path)) then call reject_generator( & matching_plane_generator_invalid_argument, & 'query and output paths must be different.', status, message & ) return end if temporary_output_path = trim(output_path)//'.beach-zhao-response.tmp' if (trim(query_path) == temporary_output_path) then call reject_generator( & matching_plane_generator_invalid_argument, & 'query path conflicts with the response generator temporary path.', status, message & ) return end if serial_mpi = mpi_context() call provider%initialize(cfg, serial_mpi, provider_status, provider_message) if (provider_status /= matching_plane_provider_ok) then call reject_generator( & matching_plane_generator_evaluation_failure, & 'Zhao response initialization failed: '//trim(provider_message), status, message & ) return end if call read_matching_plane_query_csv( & query_path, matching_plane_response_query_csv_header, queries, status, message & ) if (status /= matching_plane_generator_ok) return call validate_cartesian_query_grid(queries, status, message) if (status /= matching_plane_generator_ok) return allocate (responses(matching_plane_response_output_count, size(queries, 2))) do row = 1, size(queries, 2) call provider%evaluate( & queries(:, row), serial_mpi, responses(:, row), provider_status, provider_message & ) if (provider_status /= matching_plane_provider_ok) then call reject_generator( & matching_plane_generator_evaluation_failure, & 'Zhao response evaluation failed at query row '//trim(integer_text(row))//': '// & trim(provider_message), status, message & ) return end if end do open (newunit=output_unit, file=temporary_output_path, status='replace', action='write', iostat=ios) if (ios /= 0) then call reject_generator( & matching_plane_generator_io_error, & 'could not open matching-plane response output: '//trim(output_path), status, message & ) return end if write (output_unit, '(a,es24.16)', iostat=ios) '# matching_plane_z_m=', cfg%sim%box_max(3) if (ios == 0) write (output_unit, '(a)', iostat=ios) matching_plane_response_csv_header do row = 1, size(queries, 2) if (ios /= 0) exit write (output_unit, '(11(es24.16,:,","))', iostat=ios) queries(:, row), responses(:, row) end do close (output_unit, iostat=close_ios) if (ios == 0) ios = close_ios if (ios /= 0) then call reject_generator( & matching_plane_generator_io_error, & 'failed while writing matching-plane response output.', status, message & ) return end if call atomic_rename(temporary_output_path, trim(output_path), rename_status) if (rename_status /= filesystem_success) then call reject_generator( & matching_plane_generator_io_error, & 'failed to atomically publish matching-plane response output.', status, message & ) end if end subroutine generate_matching_plane_zhao_response_table subroutine validate_cartesian_query_grid(queries, status, message) real(dp), intent(in) :: queries(:, :) integer(i32), intent(out) :: status character(len=*), intent(out) :: message real(dp), allocatable :: axes(:, :) logical, allocatable :: filled(:) integer(i32) :: axis_sizes(matching_plane_response_input_count) integer :: axis, row, candidate, coordinate_index integer(i64) :: point_count, linear_index, stride logical :: found call accept_generator(status, message) if (size(queries, 1) /= matching_plane_response_input_count .or. size(queries, 2) <= 0) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane query grid has invalid dimensions.', status, message & ) return end if if (any(queries(2:5, :) < 0.0_dp)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane query grid contains negative flux or energy values.', status, message & ) return end if allocate (axes(size(queries, 2), matching_plane_response_input_count)) axes = 0.0_dp axis_sizes = 0_i32 do axis = 1, matching_plane_response_input_count do row = 1, size(queries, 2) found = .false. do candidate = 1, axis_sizes(axis) if (axes(candidate, axis) == queries(axis, row)) then found = .true. exit end if end do if (found) cycle axis_sizes(axis) = axis_sizes(axis) + 1_i32 axes(axis_sizes(axis), axis) = queries(axis, row) end do if (axis == 2 .and. .not. any(axes(:axis_sizes(axis), axis) == 0.0_dp)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'the matching-plane photoelectron-flux axis must include zero.', status, message & ) return end if if (axis == 3 .and. axis_sizes(axis) /= 1_i32) then call reject_generator( & matching_plane_generator_invalid_grid, & 'online Zhao table generation requires a singleton photoelectron-energy axis.', & status, message & ) return end if if (axis == 3 .and. any(axes(:axis_sizes(2), 2) > 0.0_dp) .and. & axes(1, axis) <= 0.0_dp) then call reject_generator( & matching_plane_generator_invalid_grid, & 'positive photoelectron flux requires a positive singleton energy node.', & status, message & ) return end if if (axis >= 4 .and. & (axis_sizes(axis) /= 1_i32 .or. axes(1, axis) /= 0.0_dp)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'online Zhao ambient outward axes must be singleton zero nodes.', status, message & ) return end if end do point_count = 1_i64 do axis = 1, matching_plane_response_input_count if (point_count > huge(point_count)/int(axis_sizes(axis), i64)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane Cartesian query size overflowed.', status, message & ) return end if point_count = point_count*int(axis_sizes(axis), i64) end do if (point_count /= int(size(queries, 2), i64)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane queries must form one complete Cartesian product.', status, message & ) return end if allocate (filled(size(queries, 2)), source=.false.) do row = 1, size(queries, 2) linear_index = 1_i64 stride = 1_i64 do axis = 1, matching_plane_response_input_count coordinate_index = 0 do candidate = 1, axis_sizes(axis) if (axes(candidate, axis) == queries(axis, row)) then coordinate_index = candidate exit end if end do if (coordinate_index == 0) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane query could not be mapped to a canonical axis.', status, message & ) return end if linear_index = linear_index + int(coordinate_index - 1, i64)*stride stride = stride*int(axis_sizes(axis), i64) end do if (filled(int(linear_index))) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane query grid contains a duplicate coordinate.', status, message & ) return end if filled(int(linear_index)) = .true. end do if (.not. all(filled)) then call reject_generator( & matching_plane_generator_invalid_grid, & 'matching-plane queries must form one complete Cartesian product.', status, message & ) end if end subroutine validate_cartesian_query_grid function integer_text(value) result(text) integer, intent(in) :: value character(len=32) :: text write (text, '(i0)') value end function integer_text subroutine accept_generator(status, message) integer(i32), intent(out) :: status character(len=*), intent(out) :: message status = matching_plane_generator_ok message = '' end subroutine accept_generator subroutine reject_generator(code, text, status, message) integer(i32), intent(in) :: code character(len=*), intent(in) :: text integer(i32), intent(out) :: status character(len=*), intent(out) :: message status = code message = '' message = trim(text) end subroutine reject_generator end module bem_matching_plane_response_generator