poisson_nfft.F90

Go to the documentation of this file.

!! Copyright (C) 2026 Octopus developers
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module poisson_nfft_oct_m
  use debug_oct_m
  use finufft_oct_m
  use global_oct_m
  use lalg_basic_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use poisson_cutoff_oct_m
  use profiling_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::               &
    poisson_nfft_t,      &
    poisson_nfft_init,   &
    dpoisson_nfft_solve, &
    zpoisson_nfft_solve, &
    poisson_nfft_end
 
  type poisson_nfft_t
    private
    type(finufft_t)           :: finufft
    integer                   :: nn(3)
    real(real64), allocatable :: kernel(:,:,:)
  end type poisson_nfft_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine poisson_nfft_init(this, namespace, space, mesh, threshold)
    type(poisson_nfft_t), intent(out) :: this
    type(namespace_t),     intent(in)  :: namespace
    class(space_t),        intent(in)  :: space
    class(mesh_t),         intent(in)  :: mesh
    real(real64),          intent(in)  :: threshold
 
    type(finufft_t) :: options
    real(real64)    :: bounds(2, space%dim), aa(3), spacing, spacing_min, r_c, ll(3)
    integer         :: idir
 
    push_sub(poisson_nfft_init)
 
#ifndef HAVE_NFFT
    message(1) = 'PoissonSolver = nfft requires FINUFFT but that library was not linked.'
    call messages_fatal(1, namespace=namespace)
#endif
 
    assert(space%dim == 3)
 
    ! Every domain would have to reduce the whole mode grid to complete the quadrature,
    ! which does not fit for a large mesh
    if (mesh%parallel_in_domains) then
      call messages_not_implemented('nfft Poisson solver with domain parallelization', namespace=namespace)
    end if
 
    ! The density lives in the box, so these half widths about the origin bound it,
    ! whatever the shape of the box and however the mesh samples it
    bounds = mesh%box%bounds()
    aa = m_zero
    do idir = 1, space%dim
      aa(idir) = maxval(abs(bounds(:, idir)))
    end do
 
    ! The coarsest sampling of the mesh sets how far the modes reach, and only the
    ! mesh knows it
    if (mesh%use_curvilinear) then
      spacing = maxval(mesh%vol_pp(1:mesh%np))**(m_one/m_three)
      spacing_min = minval(mesh%vol_pp(1:mesh%np))**(m_one/m_three)
    else
      spacing = mesh%vol_pp(1)**(m_one/m_three)
      spacing_min = spacing
    end if
 
    ! The kernel truncated at the diameter of the box reproduces free boundary
    ! conditions, provided the mode grid keeps the periodic images beyond the cutoff.
    ! The modes reach as far as the coarsest sampling of the mesh, which is as far as a
    ! quadrature over the mesh points goes: beyond it the sum aliases and the potential
    ! degrades.
    r_c = m_two*norm2(aa)
    ll = m_two*aa + r_c
    this%nn = 2*ceiling(m_half*ll/spacing)  ! even, so that the modes straddle G = 0
    ll = this%nn*spacing                    ! only grows, so the images stay out
 
    call poisson_nfft_build_3d_0d(this, ll, r_c)
 
    call finufft_set_tolerance(options, threshold)
    call finufft_init(this%finufft, options, this%nn, space%dim, [mesh%np, 1, 1])
    call finufft_set_nodes(this%finufft, mesh%x_t(1:mesh%np, 1), mesh%x_t(1:mesh%np, 2), &
      mesh%x_t(1:mesh%np, 3), ll)
 
    call messages_write("Info: NFFT Poisson solver")
    call messages_new_line()
    call messages_write("      Fourier modes             =")
    do idir = 1, 3
      call messages_write(this%nn(idir))
      if (idir < 3) call messages_write(" x ")
    end do
    call messages_new_line()
    call messages_write("      Mesh spacing, min .. max  =")
    call messages_write(units_from_atomic(units_out%length, spacing_min))
    call messages_write(" ..")
    call messages_write(units_from_atomic(units_out%length, spacing))
    call messages_write(' [' // trim(units_abbrev(units_out%length)) // ']')
    call messages_new_line()
    call messages_write("      Cutoff radius             =")
    call messages_write(units_from_atomic(units_out%length, r_c))
    call messages_write(' [' // trim(units_abbrev(units_out%length)) // ']')
    call messages_info(namespace=namespace)
 
    pop_sub(poisson_nfft_init)
  end subroutine poisson_nfft_init
 
  ! ---------------------------------------------------------
  subroutine poisson_nfft_build_3d_0d(this, ll, r_c)
    type(poisson_nfft_t), intent(inout) :: this
    real(real64),         intent(in)    :: ll(3)
    real(real64),         intent(in)    :: r_c
 
    real(real64) :: dg(3), gg(3), gmod2
    integer      :: i1, i2, i3
 
    push_sub(poisson_nfft_build_3d_0d)
 
    dg = m_two*m_pi/ll
 
    safe_allocate(this%kernel(1:this%nn(1), 1:this%nn(2), 1:this%nn(3)))
 
    ! modes -N/2 .. N/2-1 per direction, the ordering FINUFFT uses
    !$omp parallel do private(i1, i2, gg, gmod2)
    do i3 = 1, this%nn(3)
      gg(3) = (i3 - this%nn(3)/2 - 1)*dg(3)
      do i2 = 1, this%nn(2)
        gg(2) = (i2 - this%nn(2)/2 - 1)*dg(2)
        do i1 = 1, this%nn(1)
          gg(1) = (i1 - this%nn(1)/2 - 1)*dg(1)
 
          gmod2 = sum(gg**2)
          if (gmod2 > m_epsilon) then
            this%kernel(i1, i2, i3) = poisson_cutoff_3d_0d(sqrt(gmod2), r_c)/gmod2
          else
            this%kernel(i1, i2, i3) = m_half*r_c**2
          end if
        end do
      end do
    end do
    !$omp end parallel do
 
    ! mode -N/2 has no partner at +N/2, and that asymmetry leaves the convolution of a
    ! real density with a spurious imaginary part, so drop it
    this%kernel(1, :, :) = m_zero
    this%kernel(:, 1, :) = m_zero
    this%kernel(:, :, 1) = m_zero
 
    call lalg_scal(this%nn(1), this%nn(2), this%nn(3), m_four*m_pi/product(ll), this%kernel)
 
    pop_sub(poisson_nfft_build_3d_0d)
  end subroutine poisson_nfft_build_3d_0d
 
  ! ---------------------------------------------------------
  subroutine dpoisson_nfft_solve(this, mesh, pot, rho)
    type(poisson_nfft_t),    intent(in)  :: this
    class(mesh_t),            intent(in)  :: mesh
    real(real64), contiguous, intent(out) :: pot(:)
    real(real64), contiguous, intent(in)  :: rho(:)
 
    complex(real64), allocatable :: values(:)
    integer :: ip
 
    push_sub(dpoisson_nfft_solve)
 
    safe_allocate(values(1:mesh%np))
 
    if (mesh%use_curvilinear) then
      !$omp parallel do
      do ip = 1, mesh%np
        values(ip) = mesh%vol_pp(ip)*rho(ip)
      end do
      !$omp end parallel do
    else
      !$omp parallel do
      do ip = 1, mesh%np
        values(ip) = rho(ip)
      end do
      !$omp end parallel do
      call lalg_scal(mesh%np, mesh%vol_pp(1), values)
    end if
 
    call poisson_nfft_convolve(this, values)
 
    !$omp parallel do
    do ip = 1, mesh%np
      pot(ip) = real(values(ip), real64)
    end do
    !$omp end parallel do
 
    safe_deallocate_a(values)
 
    pop_sub(dpoisson_nfft_solve)
  end subroutine dpoisson_nfft_solve
 
  ! ---------------------------------------------------------
  subroutine zpoisson_nfft_solve(this, mesh, pot, rho)
    type(poisson_nfft_t),       intent(in)  :: this
    class(mesh_t),               intent(in)  :: mesh
    complex(real64), contiguous, intent(out) :: pot(:)
    complex(real64), contiguous, intent(in)  :: rho(:)
 
    complex(real64), allocatable :: values(:)
    integer :: ip
 
    push_sub(zpoisson_nfft_solve)
 
    safe_allocate(values(1:mesh%np))
 
    if (mesh%use_curvilinear) then
      !$omp parallel do
      do ip = 1, mesh%np
        values(ip) = mesh%vol_pp(ip)*rho(ip)
      end do
      !$omp end parallel do
    else
      call lalg_copy(mesh%np, rho, values)
      call lalg_scal(mesh%np, mesh%vol_pp(1), values)
    end if
 
    call poisson_nfft_convolve(this, values)
 
    call lalg_copy(mesh%np, values, pot)
 
    safe_deallocate_a(values)
 
    pop_sub(zpoisson_nfft_solve)
  end subroutine zpoisson_nfft_solve
 
  ! ---------------------------------------------------------
  subroutine poisson_nfft_convolve(this, values)
    type(poisson_nfft_t),       intent(in)    :: this
    complex(real64), contiguous, intent(inout) :: values(:)
 
    complex(real64), allocatable :: modes(:,:,:)
    integer :: i1, i2, i3
 
    push_sub(poisson_nfft_convolve)
 
    call profiling_in("POISSON_NFFT_SOLVE")
 
    safe_allocate(modes(1:this%nn(1), 1:this%nn(2), 1:this%nn(3)))
 
    call finufft_nodes_to_modes(this%finufft, values, modes)
 
    !$omp parallel do private(i1, i2)
    do i3 = 1, this%nn(3)
      do i2 = 1, this%nn(2)
        do i1 = 1, this%nn(1)
          modes(i1, i2, i3) = this%kernel(i1, i2, i3)*modes(i1, i2, i3)
        end do
      end do
    end do
    !$omp end parallel do
 
    call finufft_modes_to_nodes(this%finufft, modes, values)
 
    safe_deallocate_a(modes)
 
    call profiling_out("POISSON_NFFT_SOLVE")
 
    pop_sub(poisson_nfft_convolve)
  end subroutine poisson_nfft_convolve
 
  ! ---------------------------------------------------------
  subroutine poisson_nfft_end(this)
    type(poisson_nfft_t), intent(inout) :: this
 
    push_sub(poisson_nfft_end)
 
    call finufft_end(this%finufft)
    safe_deallocate_a(this%kernel)
 
    pop_sub(poisson_nfft_end)
  end subroutine poisson_nfft_end
 
end module poisson_nfft_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: