bem_field_solver.f90 Source File

!> 粒子位置での電場評価を direct / treecode / fmm で切り替える場ソルバ。
module bem_field_solver
!$ use omp_lib
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  use bem_kinds, only: dp, i32
  use bem_constants, only: k_coulomb
  use bem_types, only: mesh_type, sim_config, bc_periodic
  use bem_coulomb_fmm_core, only: fmm_options_type, fmm_plan_type, fmm_state_type
  use bem_string_utils, only: lower_ascii
  use bem_physics_config_types, only: field_physics_config, periodic2_physics_config, panel_kernel_config
  implicit none
  private

  integer(i32), parameter, public :: field_solver_fmm_expansion_order = 4_i32

  type :: field_solver_type
    character(len=16) :: mode = 'direct'
    character(len=16) :: field_bc_mode = 'free'
    character(len=16) :: field_normalization = 'si'
    real(dp) :: field_length_scale = 1.0d0
    real(dp) :: field_origin(3) = 0.0d0
    real(dp) :: field_inv_length_scale = 1.0d0
    real(dp) :: field_output_scale = k_coulomb
    real(dp) :: potential_output_scale = k_coulomb
    real(dp) :: theta = 0.5d0
    integer(i32) :: leaf_max = 16_i32
    integer(i32) :: min_nelem = 256_i32
    logical :: use_periodic2 = .false.
    integer(i32) :: periodic_axes(2) = 0_i32
    real(dp) :: periodic_len(2) = 0.0d0
    integer(i32) :: periodic_image_layers = 1_i32
    character(len=16) :: periodic_far_correction = 'none'
    real(dp) :: periodic_ewald_alpha = 0.0d0
    integer(i32) :: periodic_ewald_layers = 4_i32
    character(len=256) :: periodic_cache_dir = '.beach_cache/periodic2'
    real(dp) :: periodic_generation_tolerance = 1.0d-8
    real(dp) :: target_box_min(3) = 0.0d0
    real(dp) :: target_box_max(3) = 0.0d0
    logical :: tree_ready = .false.
    integer(i32) :: nelem = 0_i32
    integer(i32) :: max_node = 0_i32
    integer(i32) :: nnode = 0_i32
    integer(i32), allocatable :: elem_order(:)
    integer(i32), allocatable :: node_start(:), node_count(:)
    integer(i32), allocatable :: child_count(:), child_idx(:, :), child_octant(:, :)
    integer(i32), allocatable :: node_depth(:)
    integer(i32) :: node_max_depth = 0_i32
    integer(i32), allocatable :: node_level_start(:), node_level_nodes(:)
    real(dp), allocatable :: node_center(:, :)
    real(dp), allocatable :: node_half_size(:, :)
    real(dp), allocatable :: node_radius(:)
    real(dp), allocatable :: node_q(:), node_abs_q(:)
    real(dp), allocatable :: node_qx(:), node_qy(:), node_qz(:)
    real(dp), allocatable :: node_charge_center(:, :)
    logical :: fmm_core_ready = .false.
    type(fmm_options_type) :: fmm_core_options = fmm_options_type()
    type(fmm_plan_type) :: fmm_core_plan
    type(fmm_state_type) :: fmm_core_state = fmm_state_type()
  contains
    procedure :: init => init_field_solver
    procedure :: refresh => refresh_field_solver
    procedure :: eval_e => eval_e_field_solver
    procedure :: eval_potential => eval_potential_field_solver
    procedure :: compute_mesh_potential => compute_mesh_potential_field_solver
    procedure :: compute_cached_kneq0_mesh_potential_step => compute_cached_kneq0_mesh_potential_step_field_solver
  end type field_solver_type

  public :: field_solver_type
  public :: resolve_field_solver_mode
  public :: resolve_field_solver_tree_params

  interface
    !> 設定とメッシュから電場ソルバを初期化する。
    module subroutine init_field_solver(self, mesh, sim, field_config, periodic_config, panel_config)
      class(field_solver_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
    end subroutine init_field_solver

    !> 要素数に応じて treecode の代表パラメータを推定する。
    pure module subroutine estimate_auto_tree_params(nelem, theta, leaf_max)
      integer(i32), intent(in) :: nelem
      real(dp), intent(out) :: theta
      integer(i32), intent(out) :: leaf_max
    end subroutine estimate_auto_tree_params

    !> solver mode、要素数、明示overrideから実際に使うtree/FMMパラメータを解決する。
    pure module subroutine resolve_field_solver_tree_params(nelem, sim, theta, leaf_max)
      integer(i32), intent(in) :: nelem
      type(sim_config), intent(in) :: sim
      real(dp), intent(out) :: theta
      integer(i32), intent(out) :: leaf_max
    end subroutine resolve_field_solver_tree_params

    !> solver指定と要素数から実際に使用する direct/treecode/fmm mode を解決する。
    module function resolve_field_solver_mode(nelem, sim) result(mode)
      integer(i32), intent(in) :: nelem
      type(sim_config), intent(in) :: sim
      character(len=16) :: mode
    end function resolve_field_solver_mode

    !> 現在の要素電荷から treecode/FMM モーメントを再計算する。
    module subroutine refresh_field_solver(self, mesh)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
    end subroutine refresh_field_solver

    !> 観測点 `r` の電場を設定されたソルバで評価する。
    module subroutine eval_e_field_solver(self, mesh, r, e)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: r(3)
      real(dp), intent(out) :: e(3)
    end subroutine eval_e_field_solver

    !> 観測点 `r` の電位を設定されたソルバで評価する。
    module subroutine eval_potential_field_solver(self, mesh, sim, r, phi)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      type(sim_config), intent(in) :: sim
      real(dp), intent(in) :: r(3)
      real(dp), intent(out) :: phi
    end subroutine eval_potential_field_solver

    !> メッシュ重心を octree 分割して treecode トポロジを構築する。
    module subroutine build_tree_topology(self, mesh)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
    end subroutine build_tree_topology

    !> 要素添字区間を1ノードとして登録し、必要なら8分割で子ノードを作る。
    recursive module subroutine build_node(self, mesh, node_idx, start_idx, end_idx, depth)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      integer(i32), intent(in) :: node_idx, start_idx, end_idx, depth
    end subroutine build_node

    !> treecode 判定に使う8分木の octant 添字を返す。
    pure module function octant_index(x, y, z, center) result(oct)
      real(dp), intent(in) :: x, y, z
      real(dp), intent(in) :: center(3)
      integer(i32) :: oct
    end function octant_index

    !> ノードを再帰走査し、葉では direct 総和、遠方は monopole 近似を適用する。
    recursive module subroutine traverse_node(self, mesh, node_idx, rx, ry, rz, ex, ey, ez)
      class(field_solver_type), intent(in) :: self
      type(mesh_type), intent(in) :: mesh
      integer(i32), intent(in) :: node_idx
      real(dp), intent(in) :: rx, ry, rz
      real(dp), intent(inout) :: ex, ey, ez
    end subroutine traverse_node

    !> ノードを再帰走査し、葉では direct 総和、遠方は monopole 近似で電位を加算する。
    recursive module subroutine traverse_potential_node(self, mesh, node_idx, rx, ry, rz, phi_sum)
      class(field_solver_type), intent(in) :: self
      type(mesh_type), intent(in) :: mesh
      integer(i32), intent(in) :: node_idx
      real(dp), intent(in) :: rx, ry, rz
      real(dp), intent(inout) :: phi_sum
    end subroutine traverse_potential_node

    !> ノード半径・距離と電荷符号の一貫性から遠方近似を判定する。
    module function accept_node(self, node_idx, rx, ry, rz) result(accept_it)
      class(field_solver_type), intent(in) :: self
      integer(i32), intent(in) :: node_idx
      real(dp), intent(in) :: rx, ry, rz
      logical :: accept_it
    end function accept_node

    !> ノード配列を要求サイズで確保し、未使用要素をゼロ初期化する。
    module subroutine ensure_tree_capacity(self, max_node_needed)
      class(field_solver_type), intent(inout) :: self
      integer(i32), intent(in) :: max_node_needed
    end subroutine ensure_tree_capacity

    !> source tree ノードを深さごとの連続バケットへ並べ替える。
    module subroutine rebuild_source_level_cache(self)
      class(field_solver_type), intent(inout) :: self
    end subroutine rebuild_source_level_cache

    !> treecode 作業配列を解放する。
    module subroutine reset_tree_storage(self)
      class(field_solver_type), intent(inout) :: self
    end subroutine reset_tree_storage

    !> FMM のパネル幾何と電荷状態を構築・更新する。
    module subroutine refresh_fmm_solver(self, mesh)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
    end subroutine refresh_fmm_solver

    !> メッシュ重心での電位を計算する（FMM/direct 自動切替）。
    module subroutine compute_mesh_potential_field_solver(self, mesh, sim, potential_v)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      type(sim_config), intent(in) :: sim
      real(dp), intent(out) :: potential_v(:)
    end subroutine compute_mesh_potential_field_solver

    !> cached k/=0 演算子を電荷増分へ作用させ、要素重心での電位増分を返す。
    !!
    !! 評価中だけ FMM state を `charge_step` へ更新し、返る前に必ず
    !! 現在の `mesh%q_elem` へ戻す。cached_kneq0 FMM 以外は受理しない。
    module subroutine compute_cached_kneq0_mesh_potential_step_field_solver(self, mesh, charge_step, potential_step_v)
      class(field_solver_type), intent(inout) :: self
      type(mesh_type), intent(in) :: mesh
      real(dp), intent(in) :: charge_step(:)
      real(dp), intent(out) :: potential_step_v(:)
    end subroutine compute_cached_kneq0_mesh_potential_step_field_solver

  end interface

end module bem_field_solver