electric_field_at Subroutine

public subroutine electric_field_at(mesh, r, softening, e)

全要素電荷を点電荷近似で総和し、softening付きで観測点 r の電場ベクトルを返す。

Arguments

Type IntentOptional Attributes Name
type(mesh_type), intent(in) :: mesh

要素重心 centers と要素電荷 q_elem を保持したメッシュ情報。
real(kind=dp), intent(in) :: r(3)
real(kind=dp), intent(in) :: softening

特異点回避のために距離2乗へ加える softening 長さ [m]。
real(kind=dp), intent(out) :: e(3)

Source Code

  subroutine electric_field_at(mesh, r, softening, e)
    type(mesh_type), intent(in) :: mesh
    real(dp), intent(in) :: r(3)
    real(dp), intent(in) :: softening
    real(dp), intent(out) :: e(3)

    integer(i32) :: i
    real(dp) :: soft2, r2, inv_r3, ex, ey, ez
    real(dp) :: rx, ry, rz, dx, dy, dz, qi

    ex = 0.0d0
    ey = 0.0d0
    ez = 0.0d0
    soft2 = softening*softening
    rx = r(1)
    ry = r(2)
    rz = r(3)

    !$omp simd reduction(+:ex,ey,ez) private(dx,dy,dz,r2,inv_r3,qi)
    do i = 1, mesh%nelem
      dx = rx - mesh%center_x(i)
      dy = ry - mesh%center_y(i)
      dz = rz - mesh%center_z(i)
      r2 = dx*dx + dy*dy + dz*dz + soft2
      inv_r3 = 1.0d0/(sqrt(r2)*r2)
      qi = mesh%q_elem(i)
      ex = ex + qi*inv_r3*dx
      ey = ey + qi*inv_r3*dy
      ez = ez + qi*inv_r3*dz
    end do

    e(1) = k_coulomb*ex
    e(2) = k_coulomb*ey
    e(3) = k_coulomb*ez
  end subroutine electric_field_at