double_grid.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 X. Andrade
!! Copyright (C) 2026 N. Tancogne-Dejean
!!
!! 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 double_grid_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
!$ use omp_lib
  use parser_oct_m
  use profiling_oct_m
  use splines_oct_m
  use submesh_oct_m
 
  implicit none
 
  private
 
  public ::                  &
    double_grid_t,           &
    double_grid_init,        &
    double_grid_end,         &
    double_grid_margin,      &
    double_grid_transform
 
  type double_grid_t
    ! Components are public by default
    logical :: enabled = .false.
    integer :: order   = 0
    integer :: divisor = 0
    integer :: nn = 0
    integer :: interp_min = 0
    integer :: interp_max = 0
    real(real64), allocatable :: co(:, :)
  end type double_grid_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine double_grid_init(this, namespace, dim)
    type(double_grid_t), intent(out) :: this
    type(namespace_t),   intent(in)  :: namespace
    integer,             intent(in)  :: dim
 
    integer :: aa, ll
    real(real64), allocatable :: points(:)
 
    push_sub(double_grid_init)
 
    !%Variable DoubleGrid
    !%Type logical
    !%Default no
    !%Section Mesh
    !%Description
    !% Enables the Ono-Hirose double-grid technique [T. Ono and K. Hirose,
    !% <i>Phys. Rev. Lett.</i> <b>82</b>, 5016 (1999)]: localized functions of the
    !% pseudopotentials (Kleinman-Bylander projectors, ultrasoft augmentation
    !% charges) are evaluated on a fine grid and folded back onto the simulation
    !% grid with the transpose of a high-order interpolation, at construction time.
    !% This increases the accuracy of ultrasoft pseudopotentials on coarse grids
    !% and reduces the egg-box effect. Only available in 3D.
    !%End
    if (dim == 3) then
      call parse_variable(namespace, 'DoubleGrid', .false., this%enabled)
    else
      this%enabled = .false.
    end if
 
    if (.not. this%enabled) then
      pop_sub(double_grid_init)
      return
    end if
 
    call messages_experimental('DoubleGrid', namespace=namespace)
 
    !%Variable DoubleGridOrder
    !%Type integer
    !%Default 9
    !%Section Mesh
    !%Description
    !% Order of the Lagrange interpolation whose transpose is used by the
    !% double-grid technique. Must be an odd number; the integration error
    !% decreases as O(h^(order+1)). Low-order interpolation is not recommended.
    !%End
    call parse_variable(namespace, 'DoubleGridOrder', 9, this%order)
    if (mod(this%order, 2) == 0 .or. this%order < 3) then
      call messages_input_error(namespace, 'DoubleGridOrder', 'must be an odd number >= 3')
    end if
 
    !%Variable DoubleGridDivisor
    !%Type integer
    !%Default 3
    !%Section Mesh
    !%Description
    !% Number of fine-grid points per coarse-grid spacing and axis used by the
    !% double-grid technique. Must be an odd number.
    !%End
    call parse_variable(namespace, 'DoubleGridDivisor', 3, this%divisor)
    if (mod(this%divisor, 2) == 0 .or. this%divisor < 3) then
      call messages_input_error(namespace, 'DoubleGridDivisor', 'must be an odd number >= 3')
    end if
 
    this%nn = (this%divisor - 1)/2
    this%interp_min = -this%order/2
    this%interp_max =  this%order/2 + 1
 
    ! Lagrange weights of the coarse points at l*h, l = interp_min:interp_max,
    ! for interpolation at the fine offset a*h/divisor, a = 1:nn.
    ! In units of h/divisor the coarse points sit at divisor*l - a relative to
    ! the fine point, so we evaluate the interpolant at 0.
    safe_allocate(this%co(this%interp_min:this%interp_max, 1:this%nn))
    safe_allocate(points(1:this%order + 1))
    do aa = 1, this%nn
      do ll = this%interp_min, this%interp_max
        points(ll - this%interp_min + 1) = real(this%divisor*ll - aa, real64)
      end do
      call interpolation_coefficients(this%order + 1, points, m_zero, this%co(:, aa))
    end do
    safe_deallocate_a(points)
 
    pop_sub(double_grid_init)
  end subroutine double_grid_init
 
  ! ---------------------------------------------------------
  subroutine double_grid_end(this)
    type(double_grid_t), intent(inout) :: this
 
    push_sub(double_grid_end)
 
    safe_deallocate_a(this%co)
    this%enabled = .false.
 
    pop_sub(double_grid_end)
  end subroutine double_grid_end
 
  ! ---------------------------------------------------------
  real(real64) function double_grid_margin(this, mesh) result(margin)
    type(double_grid_t), intent(in) :: this
    class(mesh_t),       intent(in) :: mesh
 
    integer :: i1, i2, i3
 
    margin = m_zero
    if (.not. this%enabled) return
 
    do i1 = -1, 1, 2
      do i2 = -1, 1, 2
        do i3 = -1, 1, 2
          margin = max(margin, norm2(mesh%coord_system%to_cartesian(mesh%spacing(1:3)*[i1, i2, i3])))
        end do
      end do
    end do
    margin = this%interp_max*margin
 
  end function double_grid_margin
 
  ! ---------------------------------------------------------
  subroutine double_grid_transform(this, sm, nf, spl, ll, mm, rcut, feff, cart_dir)
    type(double_grid_t), intent(in)  :: this
    type(submesh_t),     intent(in)  :: sm
    integer,             intent(in)  :: nf
    type(spline_t),      intent(in)  :: spl(:)
    integer,             intent(in)  :: ll(:)
    integer,             intent(in)  :: mm(:)
    real(real64),        intent(in)  :: rcut(:)
    real(real64),        intent(out) :: feff(:, :)
    integer, optional,   intent(in)  :: cart_dir(:)
 
    integer :: ip, jf
    real(real64) :: ylm
    real(real64), allocatable :: rcut_eff(:)
    integer, allocatable :: cdir(:)
 
    push_sub(double_grid_transform)
 
    if (sm%np == 0) then
      pop_sub(double_grid_transform)
      return
    end if
 
    assert(sm%np > 0)
    assert(size(feff, 1) >= sm%np)
    assert(size(spl) >= nf .and. size(ll) >= nf .and. size(mm) >= nf .and. size(rcut) >= nf)
 
    ! Never evaluate a spline beyond its radial range.
    safe_allocate(rcut_eff(1:nf))
    do jf = 1, nf
      rcut_eff(jf) = min(rcut(jf), spline_range_max(spl(jf)))
    end do
 
    safe_allocate(cdir(1:nf))
    cdir = 0
    if (present(cart_dir)) cdir(1:nf) = cart_dir(1:nf)
 
    if (.not. this%enabled) then
 
      do jf = 1, nf
        do ip = 1, sm%np
          if (sm%r(ip) <= rcut_eff(jf)) then
            feff(ip, jf) = spline_eval(spl(jf), sm%r(ip))
            if (ll(jf) >= 0) then
              call ylmr_real(sm%rel_x(1:3, ip), ll(jf), mm(jf), ylm)
              feff(ip, jf) = feff(ip, jf)*ylm
            end if
            if (cdir(jf) > 0) feff(ip, jf) = feff(ip, jf)*sm%rel_x(cdir(jf), ip)
          else
            feff(ip, jf) = m_zero
          end if
        end do
      end do
 
    else
 
      call double_grid_transform_kernel(this, sm, nf, spl, ll, mm, rcut_eff, cdir, feff)
 
    end if
 
    safe_deallocate_a(rcut_eff)
    safe_deallocate_a(cdir)
 
    pop_sub(double_grid_transform)
  end subroutine double_grid_transform
 
  ! ---------------------------------------------------------
  subroutine double_grid_transform_kernel(this, sm, nf, spl, ll, mm, rcut, cdir, feff)
    type(double_grid_t), intent(in)  :: this
    type(submesh_t),     intent(in)  :: sm
    integer,             intent(in)  :: nf
    type(spline_t),      intent(in)  :: spl(:)
    integer,             intent(in)  :: ll(:)
    integer,             intent(in)  :: mm(:)
    real(real64),        intent(in)  :: rcut(:)
    integer,             intent(in)  :: cdir(:)
    real(real64),        intent(out) :: feff(:, :)
 
    integer :: ip, jf, idir, io, io1, io2, io3, il1
    integer :: ia1, ia2, ia3, noff, ngrp, igrp, nthr, ithr
    integer :: kmin(1:3), kmax(1:3), key(1:3)
    integer, allocatable :: isub(:, :, :), offs(:, :), gfirst(:)
    real(real64) :: hh(1:3), chi(1:3), rmax, anchor_reach, rel_ref(1:3)
    real(real64), allocatable :: xstep(:, :), xoff(:, :), wgt(:, :), fbuf(:, :, :)
    integer, allocatable :: tgt(:, :), ntgt(:)
    type(c_ptr) :: acc1
    type(c_ptr), allocatable :: acc(:)
 
    call profiling_in("DOUBLE_GRID")
    push_sub(double_grid_transform_kernel)
 
    assert(.not. sm%mesh%use_curvilinear)
    assert(sm%mesh%box%dim == 3)
 
    hh(1:3) = sm%mesh%spacing(1:3)
    rel_ref(1:3) = sm%rel_x(1:3, 1)
    rmax = maxval(rcut(1:nf))
 
    ! Integer keys of the local submesh points, relative to point 1, and their
    ! bounding box. Coordinate differences are mapped to the primitive space,
    ! where mesh points (and their periodic replicas) are exact multiples of
    ! the spacing.
    kmin = 0
    kmax = 0
    do ip = 1, sm%np
      chi = sm%mesh%coord_system%from_cartesian(sm%rel_x(1:3, ip) - rel_ref(1:3))
      key = nint(chi/hh)
      assert(all(abs(chi/hh - key) < 1e-6_real64))
      kmin = min(kmin, key)
      kmax = max(kmax, key)
    end do
 
    safe_allocate(isub(kmin(1):kmax(1), kmin(2):kmax(2), kmin(3):kmax(3)))
    isub = 0
    do ip = 1, sm%np
      chi = sm%mesh%coord_system%from_cartesian(sm%rel_x(1:3, ip) - rel_ref(1:3))
      key = nint(chi/hh)
      assert(isub(key(1), key(2), key(3)) == 0)
      isub(key(1), key(2), key(3)) = ip
    end do
 
    ! Cartesian steps of one grid cell per axis, and of the fine offsets.
    safe_allocate(xstep(1:3, 1:3))
    do idir = 1, 3
      chi = m_zero
      chi(idir) = hh(idir)
      xstep(1:3, idir) = sm%mesh%coord_system%to_cartesian(chi)
    end do
 
    noff = this%divisor**3
    safe_allocate(offs(1:3, 1:noff))
    safe_allocate(xoff(1:3, 1:noff))
    io = 0
    do io3 = -this%nn, this%nn
      do io2 = -this%nn, this%nn
        do io1 = -this%nn, this%nn
          io = io + 1
          offs(1:3, io) = [io1, io2, io3]
          xoff(1:3, io) = (xstep(1:3, 1)*io1 + xstep(1:3, 2)*io2 + xstep(1:3, 3)*io3)/this%divisor
        end do
      end do
    end do
 
    ! A whole anchor site can be skipped when even its nearest fine point lies
    ! beyond every support radius.
    anchor_reach = rmax + 1e-12_real64
    do io = 1, noff
      anchor_reach = max(anchor_reach, rmax + norm2(xoff(1:3, io)) + 1e-12_real64)
    end do
 
    ! Per-axis scatter targets and weights for each fine offset o: the sites
    ! l*sign(o) of the mirrored stencil with weights co(l, |o|); for o = 0 the
    ! interpolation is exact at the site itself (single target, weight 1).
    safe_allocate(tgt(1:this%order + 1, -this%nn:this%nn))
    safe_allocate(wgt(1:this%order + 1, -this%nn:this%nn))
    safe_allocate(ntgt(-this%nn:this%nn))
    do io1 = -this%nn, this%nn
      if (io1 == 0) then
        ntgt(io1) = 1
        tgt(1, io1) = 0
        wgt(1, io1) = m_one
      else
        ntgt(io1) = this%order + 1
        do il1 = this%interp_min, this%interp_max
          tgt(il1 - this%interp_min + 1, io1) = il1*sign(1, io1)
          wgt(il1 - this%interp_min + 1, io1) = this%co(il1, abs(io1))
        end do
      end if
    end do
 
    ! Consecutive functions sharing a radial spline and cutoff form a group, whose
    ! spline is evaluated once per fine point.
    safe_allocate(gfirst(1:nf + 1))
    ngrp = 1
    gfirst(1) = 1
    do jf = 2, nf
      if (spline_same(spl(jf), spl(jf - 1)) .and. abs(rcut(jf) - rcut(jf - 1)) < m_epsilon) cycle
      ngrp = ngrp + 1
      gfirst(ngrp) = jf
    end do
    gfirst(ngrp + 1) = nf + 1
 
    nthr = 1
!$  nthr = omp_get_max_threads()
 
    ! Each thread owns a GSL accelerator per spline, as the shared ones are
    ! modified on evaluation. The work split keeps the result deterministic
    if (ngrp >= nthr) then
 
      feff(1:sm%np, 1:nf) = m_zero
      !$omp parallel do schedule(dynamic) private(ia1, ia2, ia3, acc1)
      do igrp = 1, ngrp
        call spline_accel_init(acc1)
        do ia3 = kmin(3) - this%interp_max, kmax(3) + this%interp_max
          do ia2 = kmin(2) - this%interp_max, kmax(2) + this%interp_max
            do ia1 = kmin(1) - this%interp_max, kmax(1) + this%interp_max
              call scatter_anchor(ia1, ia2, ia3, igrp, igrp, [acc1], feff)
            end do
          end do
        end do
        call spline_accel_end(acc1)
      end do
      !$omp end parallel do
 
    else
 
      safe_allocate(fbuf(1:sm%np, 1:nf, 0:nthr - 1))
      !$omp parallel private(ithr, igrp, ia1, ia2, ia3, acc)
      ithr = 0
!$    ithr = omp_get_thread_num()
      fbuf(1:sm%np, 1:nf, ithr) = m_zero
      safe_allocate(acc(1:ngrp))
      do igrp = 1, ngrp
        call spline_accel_init(acc(igrp))
      end do
      !$omp do schedule(static) collapse(2)
      do ia3 = kmin(3) - this%interp_max, kmax(3) + this%interp_max
        do ia2 = kmin(2) - this%interp_max, kmax(2) + this%interp_max
          do ia1 = kmin(1) - this%interp_max, kmax(1) + this%interp_max
            call scatter_anchor(ia1, ia2, ia3, 1, ngrp, acc, fbuf(:, :, ithr))
          end do
        end do
      end do
      !$omp end do
      do igrp = 1, ngrp
        call spline_accel_end(acc(igrp))
      end do
      safe_deallocate_a(acc)
      !$omp end parallel
 
      feff(1:sm%np, 1:nf) = fbuf(1:sm%np, 1:nf, 0)
      do ithr = 1, nthr - 1
        feff(1:sm%np, 1:nf) = feff(1:sm%np, 1:nf) + fbuf(1:sm%np, 1:nf, ithr)
      end do
      safe_deallocate_a(fbuf)
 
    end if
 
    feff(1:sm%np, 1:nf) = feff(1:sm%np, 1:nf)/real(noff, real64)
 
    safe_deallocate_a(isub)
    safe_deallocate_a(offs)
    safe_deallocate_a(xoff)
    safe_deallocate_a(xstep)
    safe_deallocate_a(tgt)
    safe_deallocate_a(wgt)
    safe_deallocate_a(ntgt)
    safe_deallocate_a(gfirst)
 
    pop_sub(double_grid_transform_kernel)
    call profiling_out("DOUBLE_GRID")
 
  contains
 
    subroutine scatter_anchor(ia1, ia2, ia3, ig1, ig2, accg, fout)
      integer,      intent(in)    :: ia1, ia2, ia3, ig1, ig2
      type(c_ptr),  intent(in)    :: accg(:)
      real(real64), intent(inout) :: fout(:, :)
 
      integer :: io, ig, jf, il1, il2, il3, it1, it2, it3
      real(real64) :: xanchor(1:3), xfine(1:3), rr, rad, ylm, val, ww3, ww23
 
      xanchor(1:3) = rel_ref(1:3) + xstep(1:3, 1)*ia1 + xstep(1:3, 2)*ia2 + xstep(1:3, 3)*ia3
      if (norm2(xanchor) > anchor_reach) return
 
      do io = 1, noff
        xfine(1:3) = xanchor(1:3) + xoff(1:3, io)
        rr = norm2(xfine)
        if (rr > rmax) cycle
 
        do ig = ig1, ig2
          if (rr > rcut(gfirst(ig))) cycle
          rad = spline_eval_accel(spl(gfirst(ig)), rr, accg(ig - ig1 + 1))
 
          do jf = gfirst(ig), gfirst(ig + 1) - 1
            val = rad
            if (ll(jf) >= 0) then
              call ylmr_real(xfine, ll(jf), mm(jf), ylm)
              val = val*ylm
            end if
            if (cdir(jf) > 0) val = val*xfine(cdir(jf))
            if (abs(val) <= m_epsilon) cycle
 
            ! Scatter to the tensor-product stencil around the anchor.
            do il3 = 1, ntgt(offs(3, io))
              it3 = ia3 + tgt(il3, offs(3, io))
              if (it3 < kmin(3) .or. it3 > kmax(3)) cycle
              ww3 = wgt(il3, offs(3, io))*val
              do il2 = 1, ntgt(offs(2, io))
                it2 = ia2 + tgt(il2, offs(2, io))
                if (it2 < kmin(2) .or. it2 > kmax(2)) cycle
                ww23 = wgt(il2, offs(2, io))*ww3
                do il1 = 1, ntgt(offs(1, io))
                  it1 = ia1 + tgt(il1, offs(1, io))
                  if (it1 < kmin(1) .or. it1 > kmax(1)) cycle
                  if (isub(it1, it2, it3) > 0) then
                    fout(isub(it1, it2, it3), jf) = fout(isub(it1, it2, it3), jf) &
                      + wgt(il1, offs(1, io))*ww23
                  end if
                end do
              end do
            end do
          end do
        end do
 
      end do
    end subroutine scatter_anchor
  end subroutine double_grid_transform_kernel
 
end module double_grid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: