stencil_stargeneral.F90

Go to the documentation of this file.

!! Copyright (C) 2016 U. De Giovannini, H Huebener; 2019 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 stencil_stargeneral_oct_m
  use affine_coordinates_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use messages_oct_m
  use profiling_oct_m
  use stencil_oct_m
  use sort_oct_m
 
  implicit none
 
  private
  public ::                     &
    stencil_stargeneral_size_lapl, &
    stencil_stargeneral_extent,    &
    stencil_stargeneral_get_lapl,  &
    stencil_stargeneral_pol_lapl,  &
    stencil_stargeneral_get_arms
 
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_get_arms(this, dim, coord_system)
    type(stencil_t),            intent(inout) :: this
    integer,                    intent(in)    :: dim
    class(coordinate_system_t), intent(in)    :: coord_system
 
    real(real64)   :: theta(dim), basis_vectors(dim, dim)
    integer :: arm(dim)
    integer :: ii, jj, kk
    real(real64)   :: norm
 
    integer, allocatable :: cand_vecs(:,:), idx(:)
    real(real64), allocatable :: cand_norms(:)
    integer :: n_cand, max_cand
    logical :: active(3)
    integer :: n_found
    integer :: det
    integer :: v1_xy, v1_xz, v1_yz
    integer :: v2_xy, v2_xz, v2_yz
    integer :: v3_xy, v3_xz, v3_yz
    integer :: i1, i2, i3, i, j, k
 
    push_sub(stencil_stargeneral_get_arms)
 
    assert(dim <= 3)
 
    this%stargeneral%narms = 0
    this%stargeneral%arms = 0
 
    if (dim == 1) then
      !we are done
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
    select type (coord_system)
    class is (affine_coordinates_t)
      basis_vectors = coord_system%basis%vectors
    class default
      message(1) = "Stencil stargeneral can currently only be used with affine coordinate systems."
      call messages_fatal(1)
    end select
 
    if (dim == 2) then
      !get the angle between the primitive vectors
      theta(1) = acos(dot_product(basis_vectors(:, 1), basis_vectors(:, 2)))
 
      if (theta(1) < m_pi*m_half) then
        this%stargeneral%narms = this%stargeneral%narms + 1
        arm = (/1,-1/)
        this%stargeneral%arms(this%stargeneral%narms, 1:dim) = arm
      else if (theta(1) > m_pi*m_half) then
        this%stargeneral%narms = this%stargeneral%narms + 1
        arm = (/1,+1/)
        this%stargeneral%arms(this%stargeneral%narms, 1:dim) = arm
      end if
      !if theta == pi/2 we do not need additional arms
 
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
    ! dim>2
 
    !We first count how many arms we need
    theta(1) = acos(dot_product(basis_vectors(:, 1), basis_vectors(:, 2)))
    if (abs(theta(1) - m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
    theta(2) = acos(dot_product(basis_vectors(:, 2), basis_vectors(:, 3)))
    if (abs(theta(2)-m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
    theta(3) = acos(dot_product(basis_vectors(:, 3), basis_vectors(:, 1)))
    if (abs(theta(3) - m_pi*m_half) > 1e-8_real64) this%stargeneral%narms = this%stargeneral%narms + 1
 
 
    !Cubic cell, no extra arms
    if (this%stargeneral%narms == 0) then
      pop_sub(stencil_stargeneral_get_arms)
      return
    end if
 
    ! Determine active off-diagonal components
    ! Mapping:
    ! active(1): xy (theta 12) -> theta(1)
    ! active(2): xz (theta 13) -> theta(3)
    ! active(3): yz (theta 23) -> theta(2)
    active = .false.
    if (abs(theta(1) - m_pi*m_half) > 1e-8_real64) active(1) = .true.
    if (abs(theta(3) - m_pi*m_half) > 1e-8_real64) active(2) = .true.
    if (abs(theta(2) - m_pi*m_half) > 1e-8_real64) active(3) = .true.
 
    ! Generate candidates in range [-2, 2]
    max_cand = 5**3
    safe_allocate(cand_vecs(3, max_cand))
    safe_allocate(cand_norms(max_cand))
    safe_allocate(idx(max_cand))
    n_cand = 0
 
    do ii= -2,2
      do jj = -2,2
        do kk = -2,2
          ! Skip origin
          if (ii == 0 .and. jj == 0 .and. kk == 0) cycle
 
          ! Skip axes
          if (ii /= 0 .and. jj == 0 .and. kk == 0) cycle
          if (ii == 0 .and. jj /= 0 .and. kk == 0) cycle
          if (ii == 0 .and. jj == 0 .and. kk /= 0) cycle
 
          ! Skip orthogonal planes if angle is 90
          if (.not. active(1) .and. kk == 0) cycle ! xy plane, but x-y orthogonal
          if (.not. active(3) .and. ii == 0) cycle ! yz plane, but y-z orthogonal
          if (.not. active(2) .and. jj == 0) cycle ! xz plane, but x-z orthogonal
 
          ! We only store one vector from the pair (v, -v) as they define the same line.
          ! To ensure uniqueness, we enforce that the first non-zero component is negative.
          if (ii > 0) cycle
          if (ii == 0 .and. jj > 0) cycle
          if (ii == 0 .and. jj == 0 .and. kk > 0) cycle
 
          norm = sum((ii*basis_vectors(:, 1) + jj*basis_vectors(:, 2) + kk*basis_vectors(:, 3))**2)
 
          n_cand = n_cand + 1
          cand_vecs(1, n_cand) = ii
          cand_vecs(2, n_cand) = jj
          cand_vecs(3, n_cand) = kk
          cand_norms(n_cand) = norm
        end do
      end do
    end do
 
    ! Sort candidates by norm
    call robust_sort_by_abs(cand_norms(1:n_cand), cand_vecs(:, 1:n_cand), idx(1:n_cand))
 
    ! Select narms vectors
    n_found = 0
 
    ! We iterate through sorted candidates
    ! Case 1: narms = 1
    ! This is a 3D system with two orthogonal directions.
    ! We only need to have a non-zero component along the non-orthogonal direction.
    if (this%stargeneral%narms == 1) then
      do i = 1, n_cand
        i1 = idx(i)
        ! Check if this vector has non-zero component for the active direction
        if (active(1) .and. cand_vecs(1,i1)*cand_vecs(2,i1) /= 0) n_found = 1
        if (active(2) .and. cand_vecs(1,i1)*cand_vecs(3,i1) /= 0) n_found = 1
        if (active(3) .and. cand_vecs(2,i1)*cand_vecs(3,i1) /= 0) n_found = 1
 
        if (n_found == 1) then
          this%stargeneral%arms(1, 1:3) = cand_vecs(:, i1)
          exit
        end if
      end do
    end if
 
    ! Case 2: narms = 2
    ! This corresponds to a 3D system, with one orthogonal direction.
    ! We just need to pick two linearly independent vectors.
    if (this%stargeneral%narms == 2) then
      do i = 1, n_cand
        i1 = idx(i)
        do j = i + 1, n_cand
          i2 = idx(j)
 
          ! Check linear independence (cross product non-zero)
          if ( (cand_vecs(2, i1) * cand_vecs(3, i2) - cand_vecs(3, i1) * cand_vecs(2, i2)) /= 0 .or. &
            (cand_vecs(3, i1) * cand_vecs(1, i2) - cand_vecs(1, i1) * cand_vecs(3, i2)) /= 0 .or. &
            (cand_vecs(1, i1) * cand_vecs(2, i2) - cand_vecs(2, i1) * cand_vecs(1, i2)) /= 0 ) then
 
            this%stargeneral%arms(1, 1:3) = cand_vecs(:, i1)
            this%stargeneral%arms(2, 1:3) = cand_vecs(:, i2)
            n_found = 2
            exit
          end if
        end do
        if (n_found == 2) exit
      end do
    end if
 
    ! Case 3: narms = 3 (Triclinic)
    if (this%stargeneral%narms == 3) then
      do i = 1, n_cand
        i1 = idx(i)
        do j = i + 1, n_cand
          i2 = idx(j)
          do k = j + 1, n_cand
            i3 = idx(k)
 
            ! We check if the off-diagonal terms of the outer product of the vectors
            ! form a linearly independent set. This is required for the matrix of the
            ! linear system to determine the weights to be non-singular.
            ! See Natan et al., PRB 78, 075109 (2008), eq. 19.
            ! The matrix is M_ij = (n_l)_i * (n_l)_j, where n_l are the arm vectors.
            ! The off-diagonal terms are xy, xz, yz.
            ! We form the 3x3 matrix where each row corresponds to one arm vector,
            ! and the columns correspond to the xy, xz, yz components.
 
            v1_xy = cand_vecs(1, i1) * cand_vecs(2, i1)
            v1_xz = cand_vecs(1, i1) * cand_vecs(3, i1)
            v1_yz = cand_vecs(2, i1) * cand_vecs(3, i1)
 
            v2_xy = cand_vecs(1, i2) * cand_vecs(2, i2)
            v2_xz = cand_vecs(1, i2) * cand_vecs(3, i2)
            v2_yz = cand_vecs(2, i2) * cand_vecs(3, i2)
 
            v3_xy = cand_vecs(1, i3) * cand_vecs(2, i3)
            v3_xz = cand_vecs(1, i3) * cand_vecs(3, i3)
            v3_yz = cand_vecs(2, i3) * cand_vecs(3, i3)
 
            det = v1_xy * ( v2_xz * v3_yz - v3_xz * v2_yz ) &
              - v1_xz * ( v2_xy * v3_yz - v3_xy * v2_yz ) &
              + v1_yz * ( v2_xy * v3_xz - v3_xy * v2_xz )
 
            if (det /= 0) then
              this%stargeneral%arms(1, 1:3) = cand_vecs(:, i1)
              this%stargeneral%arms(2, 1:3) = cand_vecs(:, i2)
              this%stargeneral%arms(3, 1:3) = cand_vecs(:, i3)
              n_found = 3
              exit
            end if
          end do
          if (n_found == 3) exit
        end do
        if (n_found == 3) exit
      end do
    end if
 
    if (n_found < this%stargeneral%narms) then
      write(message(1), '(a, i0, a, i0)') 'Stencil stargeneral could not find enough independent arms. Found ', &
        n_found, ' needed ', this%stargeneral%narms
      call messages_fatal(1)
    end if
 
    safe_deallocate_a(cand_vecs)
    safe_deallocate_a(cand_norms)
    safe_deallocate_a(idx)
 
    pop_sub(stencil_stargeneral_get_arms)
  end subroutine stencil_stargeneral_get_arms
 
 
  ! ---------------------------------------------------------
  integer function stencil_stargeneral_size_lapl(this, dim, order) result(n)
    type(stencil_t),     intent(inout) :: this
    integer, intent(in) :: dim
    integer, intent(in) :: order
 
    push_sub(stencil_stargeneral_size_lapl)
 
    !normal star
    n = 2*dim*order + 1
 
    ! star general
    n = n + 2 * order * this%stargeneral%narms
 
 
    pop_sub(stencil_stargeneral_size_lapl)
  end function stencil_stargeneral_size_lapl
 
 
  ! ---------------------------------------------------------
  integer function stencil_stargeneral_extent(dir, order)
    integer, intent(in) :: dir
    integer, intent(in) :: order
 
    integer :: extent
 
    push_sub(stencil_stargeneral_extent)
 
    extent = 0
    if (dir >= 1 .and. dir <= 3) then
      if (order <= 2) then
        extent = 2
      else
        extent = order
      end if
    end if
    stencil_stargeneral_extent = extent
 
    pop_sub(stencil_stargeneral_extent)
  end function stencil_stargeneral_extent
 
 
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_get_lapl(this, dim, order)
    type(stencil_t), intent(inout) :: this
    integer,         intent(in)    :: dim
    integer,         intent(in)    :: order
 
    integer :: i, j, n
    logical :: got_center
 
    push_sub(stencil_stargeneral_get_lapl)
 
    call stencil_allocate(this, dim, stencil_stargeneral_size_lapl(this, dim, order))
 
    n = 1
    select case (dim)
    case (1)
      n = 1
      do i = 1, dim
        do j = -order, order
          if (j == 0) cycle
          n = n + 1
          this%points(i, n) = j
        end do
      end do
    case (2)
      n = 1
      do i = 1, dim
        do j = -order, order
          if (j == 0) cycle
          n = n + 1
          this%points(i, n) = j
        end do
      end do
 
      do j = -order, order
        if (j == 0) cycle
        do i = 1, this%stargeneral%narms
          n = n + 1
          this%points(1:2, n) = this%stargeneral%arms(i, 1:2)*j
        end do
      end do
 
    case (3)
      got_center = .false.
 
      n = 0
      do i = 1, dim
        do j = -order, order
 
          ! count center only once
          if (j == 0) then
            if (got_center) then
              cycle
            else
              got_center = .true.
            end if
 
          end if
          n = n + 1
          this%points(i, n) = j
        end do
      end do
 
      do j = -order, order
        if (j == 0) cycle
        do i = 1, this%stargeneral%narms
          n = n + 1
          this%points(1:3, n) = this%stargeneral%arms(i, 1:3)*j
        end do
      end do
 
    end select
 
    call stencil_init_center(this)
 
    pop_sub(stencil_stargeneral_get_lapl)
  end subroutine stencil_stargeneral_get_lapl
 
 
 
 
  ! ---------------------------------------------------------
  subroutine stencil_stargeneral_pol_lapl(this, dim, order, pol)
    type(stencil_t), intent(in) :: this
    integer, intent(in)          :: dim
    integer, intent(in)          :: order
    integer, intent(out)         :: pol(:,:)
 
    integer :: i, j, n, j1
    logical :: zeros_treated(dim)
 
    push_sub(stencil_stargeneral_pol_lapl)
 
    n = 1
    select case (dim)
    case (1)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
    case (2)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
 
      do j = 1, 2*order
        do i = 1, this%stargeneral%narms
          n = n + 1
          if (sum(this%stargeneral%arms(i,1:dim)) == 0) then
            pol(1:2, n) = (/j,1/)
          else
            pol(1:2, n) = (/1,j/)
          end if
        end do
      end do
 
    case (3)
      n = 1
      pol(:,:) = 0
      do i = 1, dim
        do j = 1, 2*order
          n = n + 1
          pol(i, n) = j
        end do
      end do
 
      ! first, treat arms with zeros: they should have polynomials with zeros in the
      ! corresponding dimension
      zeros_treated = .false.
      do i = 1, this%stargeneral%narms
        if (this%stargeneral%arms(i, 1) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/0, 1, j1/)
          end do
          zeros_treated(1) = .true.
        else if (this%stargeneral%arms(i, 2) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/j1, 0, 1/)
          end do
          zeros_treated(2) = .true.
        else if (this%stargeneral%arms(i, 3) == 0) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/1, j1, 0/)
          end do
          zeros_treated(3) = .true.
        end if
      end do ! end loop on number of arms
      ! then treat the remaining arms, use the polynomials that have not been used yet
      do i = 1, this%stargeneral%narms - count(zeros_treated)
        if (.not. zeros_treated(1)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/0, 1, j1/)
          end do
          zeros_treated(1) = .true.
        else if (.not. zeros_treated(2)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/j1, 0, 1/)
          end do
          zeros_treated(2) = .true.
        else if (.not. zeros_treated(3)) then
          do j1 = 1, 2*order
            n = n + 1
            pol(1:3, n) = (/1, j1, 0/)
          end do
          zeros_treated(3) = .true.
        end if
      end do ! end loop on number of arms
    end select
 
    pop_sub(stencil_stargeneral_pol_lapl)
  end subroutine stencil_stargeneral_pol_lapl
 
 
end module stencil_stargeneral_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: