absorbing_boundaries.F90

Go to the documentation of this file.

!! Copyright (C) 2015 U. De Giovannini
!!
!! 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 absorbing_boundaries_oct_m
  use box_oct_m
  use box_cylinder_oct_m
  use box_sphere_oct_m
  use box_parallelepiped_oct_m
  use cube_function_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use io_function_oct_m
  use, intrinsic :: iso_fortran_env
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::               &
    absorbing_boundaries_init,              &
    absorbing_boundaries_end,               &
    absorbing_boundaries_t
 
  type absorbing_boundaries_t
    private
    integer,            public :: abtype
    real(real64), allocatable, public :: mf(:)
    type(cube_function_t)      :: cf
    logical                    :: ab_user_def
    real(real64), allocatable         :: ab_ufn(:)
  end type absorbing_boundaries_t
 
  integer, public, parameter :: &
    NOT_ABSORBING       = 0, &
    mask_absorbing      = 1, &
    imaginary_absorbing = 2
 
contains
 
  ! ---------------------------------------------------------
  subroutine absorbing_boundaries_init(this, namespace, space, gr)
    type(absorbing_boundaries_t), intent(out) :: this
    type(namespace_t),            intent(in)  :: namespace
    class(space_t),               intent(in)  :: space
    type(grid_t),                 intent(in)  :: gr
 
    integer             :: ip
    real(real64)        :: bounds(space%dim, 2)
    integer             :: cols_abshape_block, imdim, maxdim
 
    real(real64)        :: xx(space%dim), rr
    real(real64)        :: ufn_re, ufn_im
    character(len=1024) :: user_def_expr
 
    real(real64), allocatable  :: mf(:)
    real(real64)        :: abheight, abwidth, abwidth_def
    type(block_t)       :: blk
 
    push_sub(absorbing_boundaries_init)
 
    bounds = m_zero
 
    this%ab_user_def = .false.
 
    !%Variable AbsorbingBoundaries
    !%Type integer
    !%Default not_absorbing
    !%Section Time-Dependent::Absorbing Boundaries
    !%Description
    !% To improve the quality of the spectra by avoiding the formation of
    !% standing density waves, one can make the boundaries of the simulation
    !% box absorbing.
    !%Option not_absorbing 0
    !% Reflecting boundaries.
    !%Option mask 1
    !% Absorbing boundaries with a mask function.
    !%Option cap 2
    !% Absorbing boundaries with a complex absorbing potential.
    !%End
    call parse_variable(namespace, 'AbsorbingBoundaries', not_absorbing, this%abtype)
    if (.not. varinfo_valid_option('AbsorbingBoundaries', this%abtype)) then
      call messages_input_error(namespace, 'AbsorbingBoundaries')
    end if
 
    if (this%abtype /= not_absorbing) then
      call messages_print_with_emphasis(msg='Absorbing Boundaries', namespace=namespace)
 
      !%Variable ABCapHeight
      !%Type float
      !%Default -0.2 a.u.
      !%Section Time-Dependent::Absorbing Boundaries
      !%Description
      !% When <tt>AbsorbingBoundaries = cap</tt>, this is the height of the imaginary potential.
      !%End
      if (this%abtype == imaginary_absorbing) then
        call parse_variable(namespace, 'ABCapHeight', -0.2_real64, abheight, units_inp%energy)
      end if
 
      !%Variable ABShape
      !%Type block
      !%Section Time-Dependent::Absorbing Boundaries
      !%Description
      !% Set the shape of the absorbing boundaries. Here you can set the inner
      !% and outer bounds by setting the block as follows:
      !%
      !% <tt>%ABShape
      !% <br>&nbsp;&nbsp; inner | outer | "user-defined"
      !% <br>%</tt>
      !%
      !% The optional 3rd column is a user-defined expression for the absorbing
      !% boundaries. For example, <math>r</math> creates a spherical absorbing zone for
      !% coordinates with <math>{\tt inner} < r < {\tt outer}</math>, and <math>z</math> creates an
      !% absorbing plane.
      !% Note, values <tt>outer</tt> larger than the box size may lead in these cases to
      !% unexpected reflection behaviours.
      !% If no expression is given, the absorbing zone follows the edges of the
      !% box (not valid for user-defined box).
      !%End
 
      cols_abshape_block = 0
      if (parse_block(namespace, 'ABShape', blk) < 0) then
        message(1) = "Input: ABShape not specified. Using default values for absorbing boundaries."
        call messages_info(1, namespace=namespace)
 
        select type (sb=>gr%box)
        type is (box_sphere_t)
          bounds(1,1) = sb%radius/m_two
          bounds(1,2) = sb%radius
        type is (box_parallelepiped_t)
          bounds(:,1) = sb%half_length/m_two
          bounds(:,2) = sb%half_length
        type is (box_cylinder_t)
          bounds(1,2) = sb%half_length
          bounds(2,2) = sb%radius
          bounds(1,1) = bounds(1,1)/m_two
          bounds(2,1) = bounds(2,2)/m_two
        class default
          bounds(:,1) = m_zero
          bounds(:,2) = 0.4_real64
        end select
      else
        cols_abshape_block = parse_block_cols(blk, 0)
 
        select case (cols_abshape_block)
        case (2)
          call parse_block_float(blk, 0, 0, bounds(1,1), units_inp%length)
          call parse_block_float(blk, 0, 1, bounds(1,2), units_inp%length)
          ! Apply single inner/outer pair to all dimensions.
          bounds(:,1) = bounds(1,1)
          bounds(:,2) = bounds(1,2)
 
          select type (sb=>gr%box)
          type is (box_sphere_t)
            if (bounds(1,2) > sb%radius) then
              bounds(1,2) = sb%radius
            end if
            message(1) = "Info: using spherical absorbing boundaries."
 
          type is (box_parallelepiped_t)
            do imdim = 1, space%dim
              if (bounds(imdim,2) > sb%half_length(imdim)) then
                bounds(imdim,2) = sb%half_length(imdim)
              end if
            end do
            message(1) = "Info: using cubic absorbing boundaries."
 
          type is (box_cylinder_t)
            if (bounds(1,2) > sb%half_length) then
              bounds(1,2) = sb%half_length
            end if
            if (bounds(2,2) > sb%radius) then
              bounds(2,2) = sb%radius
            end if
            message(1) = "Info: using cylindrical absorbing boundaries."
          end select
          call messages_info(1, namespace=namespace)
 
        case (3)
          this%ab_user_def = .true.
          safe_allocate(this%ab_ufn(1:gr%np))
          this%ab_ufn = m_zero
          call parse_block_float( blk, 0, 0, bounds(1,1), units_inp%length)
          call parse_block_float( blk, 0, 1, bounds(1,2), units_inp%length)
          call parse_block_string(blk, 0, 2, user_def_expr)
          do ip = 1, gr%np
            xx = units_from_atomic(units_inp%length, gr%x(:, ip))
            rr = norm2(xx)
            call parse_expression(ufn_re, ufn_im, space%dim, xx, rr, m_zero, user_def_expr)
            this%ab_ufn(ip) = ufn_re
          end do
          message(1) = "Input: using user-defined function from expression:"
          write(message(2),'(a,a)') '   F(x,y,z) = ', trim(user_def_expr)
          call messages_info(2, namespace=namespace)
        case default
          message(1) = "Input: ABShape block must have at least 2 columns."
          call messages_fatal(1, namespace=namespace)
        end select
 
        call parse_block_end(blk)
      end if
 
      !%Variable ABWidth
      !%Type float
      !%Section Time-Dependent::Absorbing Boundaries
      !%Description
      !% Specifies the boundary width. For a finer control over the absorbing boundary
      !% shape use ABShape.
      !%End
      abwidth_def = bounds(1,2) - bounds(1,1)
      call parse_variable(namespace, 'ABWidth', abwidth_def, abwidth, units_inp%length)
      bounds(:, 1) = bounds(:, 2) - abwidth
 
      select type (sb=>gr%box)
      type is (box_sphere_t)
        maxdim = 1
      type is (box_cylinder_t)
        maxdim = 2
      class default
        if (this%ab_user_def) then
          maxdim = 1
        else
          maxdim = space%dim
        end if
      end select
      write(message(1),'(a,2a,9f12.6)') &
        " Absorbing boundary spans from [", trim(units_abbrev(units_inp%length)), '] =', &
        (units_from_atomic(units_inp%length, bounds(imdim,1)), imdim=1,maxdim)
      write(message(2),'(a,2a,9f12.6)') &
        " to [", trim(units_abbrev(units_inp%length)), '] =', &
        (units_from_atomic(units_inp%length, bounds(imdim,2)), imdim=1,maxdim)
      call messages_info(2, namespace=namespace)
 
      ! generate boundary function
      safe_allocate(mf(1:gr%np))
      call absorbing_boundaries_generate_mf(this, space, gr, bounds, mf)
 
      ! mask or cap
      safe_allocate(this%mf(1:gr%np))
 
      select case (this%abtype)
      case (mask_absorbing)
        this%mf = m_one - mf
      case (imaginary_absorbing)
        this%mf(:) = abheight * mf(:)
      end select
 
      if (debug%info) call absorbing_boundaries_write_info(this, gr, namespace, space)
 
      call messages_print_with_emphasis(namespace=namespace)
 
      safe_deallocate_a(mf)
    end if
 
    pop_sub(absorbing_boundaries_init)
  end subroutine absorbing_boundaries_init
 
  ! ---------------------------------------------------------
  subroutine absorbing_boundaries_end(this)
    type(absorbing_boundaries_t),   intent(inout) :: this
 
    push_sub(absorbing_boundaries_end)
 
    if (this%abtype /= not_absorbing) then
      safe_deallocate_a(this%mf)
    end if
 
    pop_sub(absorbing_boundaries_end)
  end subroutine absorbing_boundaries_end
 
  ! ---------------------------------------------------------
  subroutine absorbing_boundaries_write_info(this, mesh, namespace, space)
    type(absorbing_boundaries_t), intent(in) :: this
    class(mesh_t),                intent(in) :: mesh
    type(namespace_t),            intent(in) :: namespace
    class(space_t),               intent(in) :: space
 
    integer :: err
 
    if (space%dim < 3) return
 
    push_sub(absorbing_boundaries_write_info)
 
    if (this%abtype == mask_absorbing) then
      call dio_function_output(io_function_fill_how("VTK"), "./td.general", "mask", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneZ"), "./td.general", "mask", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
    else if (this%abtype == imaginary_absorbing) then
      call dio_function_output(io_function_fill_how("VTK"), "./td.general", "cap", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneZ"), "./td.general", "cap", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneX"), "./td.general", "cap", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
      call dio_function_output(io_function_fill_how("PlaneY"), "./td.general", "cap", namespace, space, mesh, &
        this%mf(1:mesh%np), unit_one, err)
    end if
 
    pop_sub(absorbing_boundaries_write_info)
  end subroutine absorbing_boundaries_write_info
 
  ! ---------------------------------------------------------
  subroutine absorbing_boundaries_generate_mf(this, space, gr, bounds, mf)
    type(absorbing_boundaries_t), intent(inout) :: this
    class(space_t), intent(in) :: space
    type(grid_t), intent(in) :: gr
    real(real64), intent(in) :: bounds(1:space%dim, 1:2)
    real(real64), contiguous, intent(inout) :: mf(:)
 
    integer :: ip, dir
    real(real64)   :: width(space%dim)
    real(real64)   :: xx(space%dim), rr, dd, ddv(space%dim), tmp(space%dim)
 
    push_sub(absorbing_boundaries_generate_mf)
 
    mf = m_zero
 
    ! generate the boundaries on the mesh
    width = bounds(:, 2) - bounds(:,1)
 
    do ip = 1, gr%np
      xx = gr%x(:, ip)
 
      if (this%ab_user_def) then
        dd = this%ab_ufn(ip) - bounds(1,1)
        if (dd > m_zero) then
          if (this%ab_ufn(ip) < bounds(1,2)) then
            mf(ip) = sin(dd * m_pi / (m_two * (width(1))))**2
          else
            mf(ip) = m_one
          end if
        end if
 
      else ! this%ab_user_def == .false.
 
        select type (sb => gr%box)
        type is (box_sphere_t)
          rr = norm2(xx - sb%center)
          dd = rr -  bounds(1,1)
          if (dd > m_zero) then
            if (dd  <  width(1)) then
              mf(ip) = sin(dd * m_pi / (m_two * (width(1))))**2
            else
              mf(ip) = m_one
            end if
          end if
 
        type is (box_parallelepiped_t)
          ! We are filling from the center opposite to the spherical case
          tmp = m_one
          mf(ip) = m_one
          ddv = abs(xx - sb%center) -  bounds(:, 1)
          do dir = 1, space%dim
            if (ddv(dir) > m_zero) then
              if (ddv(dir)  <  width(dir)) then
                tmp(dir) = m_one - sin(ddv(dir) * m_pi / (m_two * (width(dir))))**2
              else
                tmp(dir) = m_zero
              end if
            end if
            mf(ip) = mf(ip) * tmp(dir)
          end do
          mf(ip) = m_one - mf(ip)
 
        type is (box_cylinder_t)
          rr = norm2(xx(2:space%dim) - sb%center(2:space%dim))
          tmp = m_one
          mf(ip) = m_one
          ddv(1) = abs(xx(1) - sb%center(1)) - bounds(1,1)
          ddv(2) = rr         - bounds(2,1)
          do dir = 1, 2
            if (ddv(dir) > m_zero) then
              if (ddv(dir)  <  width(dir)) then
                tmp(dir) = m_one - sin(ddv(dir) * m_pi / (m_two * (width(dir))))**2
              else
                tmp(dir) = m_zero
              end if
            end if
            mf(ip) = mf(ip) * tmp(dir)
          end do
          mf(ip) = m_one - mf(ip)
 
        class default
          ! Other box shapes are not implemented
          assert(.false.)
        end select
      end if
    end do
 
    if (this%ab_user_def) then
      safe_deallocate_a(this%ab_ufn)
    end if
 
    pop_sub(absorbing_boundaries_generate_mf)
  end subroutine absorbing_boundaries_generate_mf
 
end module absorbing_boundaries_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: