!> `bem_field_solver` の初期化・設定補助手続きを実装する submodule。 submodule(bem_field_solver) bem_field_solver_config use bem_coulomb_fmm_core, only: build_plan, update_state, destroy_plan, destroy_state, fmm_options_type, fmm_state_type implicit none contains !> 設定値から direct/treecode/fmm の実行モードを確定し、必要なら木構造を初期化する。 module procedure init_field_solver character(len=16) :: requested_mode, field_bc_mode integer(i32) :: axis, n_periodic real(dp) :: span, fmm_softening real(dp), allocatable :: src_pos(:, :) self%fmm_core_plan = fmm_plan_type() self%fmm_core_state = fmm_state_type() call destroy_plan(self%fmm_core_plan) call destroy_state(self%fmm_core_state) self%fmm_use_core = .false. self%fmm_core_ready = .false. call reset_tree_storage(self) self%softening = sim%softening self%field_normalization = lower_ascii(trim(sim%field_normalization)) self%field_length_scale = resolve_field_length_scale(mesh, sim) self%field_origin = resolve_field_origin(mesh, sim, self%field_normalization) self%field_inv_length_scale = 1.0d0/self%field_length_scale self%field_output_scale = k_coulomb*self%field_inv_length_scale*self%field_inv_length_scale self%potential_output_scale = k_coulomb*self%field_inv_length_scale self%min_nelem = sim%tree_min_nelem self%nelem = mesh%nelem self%use_periodic2 = .false. self%periodic_axes = 0_i32 self%periodic_len = 0.0d0 self%periodic_image_layers = max(0_i32, sim%field_periodic_image_layers) self%periodic_far_correction = lower_ascii(trim(sim%field_periodic_far_correction)) self%periodic_ewald_alpha = sim%field_periodic_ewald_alpha self%periodic_ewald_layers = max(0_i32, sim%field_periodic_ewald_layers) self%target_box_min = 0.0d0 self%target_box_max = 0.0d0 requested_mode = lower_ascii(trim(sim%field_solver)) field_bc_mode = lower_ascii(trim(sim%field_bc_mode)) self%field_bc_mode = field_bc_mode select case (trim(field_bc_mode)) case ('free') self%periodic_far_correction = 'none' continue case ('periodic2') if (trim(requested_mode) /= 'fmm') then error stop 'sim.field_bc_mode must be "free" unless sim.field_solver="fmm".' end if if (.not. sim%use_box) then error stop 'sim.field_bc_mode="periodic2" requires sim.use_box=true.' end if n_periodic = 0_i32 do axis = 1_i32, 3_i32 if ((sim%bc_low(axis) == bc_periodic) .neqv. (sim%bc_high(axis) == bc_periodic)) then error stop 'periodic2 requires bc_low(axis)=bc_high(axis)=periodic for periodic axes.' end if if (sim%bc_low(axis) == bc_periodic) then n_periodic = n_periodic + 1_i32 if (n_periodic <= 2_i32) self%periodic_axes(n_periodic) = axis end if end do if (n_periodic /= 2_i32) then error stop 'sim.field_bc_mode="periodic2" requires exactly two periodic axes.' end if do axis = 1_i32, 2_i32 span = sim%box_max(self%periodic_axes(axis)) - sim%box_min(self%periodic_axes(axis)) if (span <= 0.0d0) error stop 'periodic2 requires positive box length on periodic axes.' self%periodic_len(axis) = span end do self%target_box_min = sim%box_min self%target_box_max = sim%box_max select case (trim(self%periodic_far_correction)) case ('auto') self%periodic_far_correction = 'none' case ('none') continue case ('m2l_root_oracle') continue case default error stop 'periodic2 far correction supports "auto", "none", '// & 'or "m2l_root_oracle" only.' end select self%use_periodic2 = .true. case default error stop 'Unknown sim.field_bc_mode in field solver init.' end select select case (trim(requested_mode)) case ('direct') self%mode = 'direct' case ('treecode') self%mode = 'treecode' case ('fmm') self%mode = 'fmm' case ('auto') if (mesh%nelem >= self%min_nelem) then self%mode = 'treecode' else self%mode = 'direct' end if case default error stop 'Unknown sim.field_solver in field solver init.' end select if (trim(self%mode) == 'fmm' .and. sim%use_box) then if (any(sim%box_max <= sim%box_min)) then error stop 'sim.use_box=true requires positive box extents for fmm dual-target.' end if self%target_box_min = sim%box_min self%target_box_max = sim%box_max end if if (trim(self%mode) == 'treecode' .or. trim(self%mode) == 'fmm') then call estimate_auto_tree_params(mesh%nelem, self%theta, self%leaf_max) if (sim%has_tree_theta) self%theta = sim%tree_theta if (sim%has_tree_leaf_max) self%leaf_max = sim%tree_leaf_max else self%theta = sim%tree_theta self%leaf_max = sim%tree_leaf_max end if self%fmm_core_options = fmm_options_type() fmm_softening = self%softening*self%field_inv_length_scale self%fmm_core_options%theta = self%theta self%fmm_core_options%leaf_max = self%leaf_max self%fmm_core_options%order = 4_i32 self%fmm_core_options%softening = fmm_softening self%fmm_core_options%use_periodic2 = self%use_periodic2 self%fmm_core_options%periodic_far_correction = self%periodic_far_correction self%fmm_core_options%periodic_axes = self%periodic_axes self%fmm_core_options%periodic_len = self%periodic_len*self%field_inv_length_scale self%fmm_core_options%periodic_image_layers = self%periodic_image_layers self%fmm_core_options%periodic_ewald_alpha = self%periodic_ewald_alpha*self%field_length_scale self%fmm_core_options%periodic_ewald_layers = self%periodic_ewald_layers self%fmm_core_options%target_box_min = (self%target_box_min - self%field_origin)*self%field_inv_length_scale self%fmm_core_options%target_box_max = (self%target_box_max - self%field_origin)*self%field_inv_length_scale if (trim(self%mode) == 'fmm') then select case (trim(self%periodic_far_correction)) case ('auto', 'none', 'm2l_root_oracle') continue case default error stop 'FMM core supports periodic far correction "auto", "none", or "m2l_root_oracle" only.' end select self%fmm_use_core = .true. if (mesh%nelem > 0_i32) then call build_core_source_positions(mesh, src_pos, self%field_inv_length_scale, self%field_origin) call build_plan(self%fmm_core_plan, src_pos, self%fmm_core_options) call update_state(self%fmm_core_plan, self%fmm_core_state, mesh%q_elem) deallocate (src_pos) self%fmm_core_ready = self%fmm_core_plan%built .and. self%fmm_core_state%ready call sync_core_plan_view(self) end if return end if if (trim(self%mode) == 'treecode' .and. mesh%nelem > 0_i32) then call build_tree_topology(self, mesh) call refresh_field_solver(self, mesh) else call reset_tree_storage(self) self%tree_ready = .false. end if end procedure init_field_solver !> FMM core 内部で使う長さ正規化スケール `L0` [m] を解決する。 real(dp) function resolve_field_length_scale(mesh, sim) result(length_scale) type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim real(dp) :: bbox_min(3), bbox_max(3), span(3) character(len=16) :: mode mode = lower_ascii(trim(sim%field_normalization)) select case (trim(mode)) case ('si') length_scale = 1.0d0 case ('length') length_scale = sim%field_length_scale case ('box') if (.not. sim%use_box) error stop 'sim.field_normalization="box" requires sim.use_box=true.' span = sim%box_max - sim%box_min if (any(span <= 0.0d0)) error stop 'sim.field_normalization="box" requires positive sim.box extents.' length_scale = maxval(span) case ('mesh') if (mesh%nelem <= 0_i32) then length_scale = sim%field_length_scale else bbox_min = min(minval(mesh%v0, dim=2), min(minval(mesh%v1, dim=2), minval(mesh%v2, dim=2))) bbox_max = max(maxval(mesh%v0, dim=2), max(maxval(mesh%v1, dim=2), maxval(mesh%v2, dim=2))) span = bbox_max - bbox_min length_scale = maxval(span) end if case default error stop 'Unknown sim.field_normalization in field solver init.' end select if (length_scale <= 0.0d0) error stop 'Resolved field normalization length must be > 0.' end function resolve_field_length_scale !> FMM core 内部で使う座標原点 [m] を解決する。 function resolve_field_origin(mesh, sim, normalization) result(origin) type(mesh_type), intent(in) :: mesh type(sim_config), intent(in) :: sim character(len=*), intent(in) :: normalization real(dp) :: origin(3) character(len=16) :: mode mode = lower_ascii(trim(normalization)) select case (trim(mode)) case ('box') origin = sim%box_min case ('mesh') if (mesh%nelem <= 0_i32) then origin = 0.0d0 else origin = min(minval(mesh%v0, dim=2), min(minval(mesh%v1, dim=2), minval(mesh%v2, dim=2))) end if case default origin = 0.0d0 end select end function resolve_field_origin !> 要素数レンジに応じた treecode/FMM の `theta` と `leaf_max` 推奨値を返す。 module procedure estimate_auto_tree_params if (nelem < 1500_i32) then theta = 0.40d0 leaf_max = 12_i32 else if (nelem < 10000_i32) then theta = 0.50d0 leaf_max = 16_i32 else if (nelem < 50000_i32) then theta = 0.58d0 leaf_max = 20_i32 else theta = 0.65d0 leaf_max = 24_i32 end if end procedure estimate_auto_tree_params end submodule bem_field_solver_config