external_densities.F90

Go to the documentation of this file.

!! Copyright (C) 2019 R. Jestaedt, H. Appel
!! Copyright (C) 2020 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 external_densities_oct_m
  use debug_oct_m
  use global_oct_m
  use math_oct_m
  use states_mxll_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use tdfunction_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                               &
    get_rs_density_ext,                   &
    external_current_init,                &
    external_current_calculation
 
  integer, parameter, public ::           &
    EXTERNAL_CURRENT_PARSER      = 0,     &
    external_current_td_function = 1
 
contains
 
  !----------------------------------------------------------
  subroutine get_rs_density_ext(st, space, mesh, time, rs_current_density_ext)
    type(states_mxll_t), intent(inout) :: st
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: time
    complex(real64),     optional, intent(inout) :: rs_current_density_ext(:,:)
 
    real(real64), allocatable :: current(:,:,:)
 
    push_sub(get_rs_density_ext)
 
    safe_allocate(current(1:mesh%np, 1:space%dim, 1))  
 
    call external_current_calculation(st, space, mesh, time, current(:, :, 1))
    call build_rs_current_state(current(:, :, 1), mesh, rs_current_density_ext, st%ep, mesh%np)
 
    safe_deallocate_a(current)
 
    pop_sub(get_rs_density_ext)
  end subroutine get_rs_density_ext
 
 
  !----------------------------------------------------------
  subroutine external_current_init(st, space, namespace, mesh)
    type(states_mxll_t), intent(inout) :: st
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(namespace_t),   intent(in)    :: namespace
 
    type(block_t)        :: blk
    integer              :: ip, il, nlines, ncols, idir, ierr
    real(real64)         :: j_vector(space%dim), dummy(space%dim), xx(1:space%dim), rr, omega
    character(len=1024)  :: tdf_expression, phase_expression
 
 
    push_sub(external_current_init)
 
    call profiling_in('EXTERNAL_CURRENT_INIT')
 
    !%Variable UserDefinedMaxwellExternalCurrent
    !%Type block
    !%Section Maxwell
    !%Description
    !%
    !% Example:
    !%
    !% <tt>%UserDefinedMaxwellExternalCurrent
    !% <br>&nbsp;&nbsp; current_parser      | "expression_x_dir1" | "expression_y_dir1" | "expression_z_dir1"
    !% <br>&nbsp;&nbsp; current_parser      | "expression_x_dir2" | "expression_y_dir2" | "expression_z_dir2"
    !% <br>&nbsp;&nbsp; current_td_function | "amplitude_j0_x" | "amplitude_j0_y" | "amplitude_j0_z" | omega |
    !%                  envelope_td_function_name | phase
    !% <br>%</tt>
    !%
    !% Description about UserDefinedMaxwellExternalCurrent follows
    !%
    !%Option current_parser 0
    !% description follows
    !%Option current_td_function 1
    !% description follows
    !%End
 
    if (parse_block(namespace, 'UserDefinedMaxwellExternalCurrent', blk) == 0) then
 
      ! find out how many lines (i.e. states) the block has
      nlines = parse_block_n(blk)
 
      st%external_current_number = nlines
      safe_allocate(st%external_current_modus(1:nlines))
      safe_allocate(st%external_current_string(1:space%dim, 1:nlines))
      safe_allocate(st%external_current_amplitude(1:mesh%np, 1:space%dim, 1:nlines))
      safe_allocate(st%external_current_td_function(1:nlines))
      safe_allocate(st%external_current_omega(1:nlines))
      safe_allocate(st%external_current_td_phase(1:nlines))
 
      ! read all lines
      do il = 1, nlines
        ! Check that number of columns is four, five, six or seven.
        ncols = parse_block_cols(blk, il - 1)
        if ((ncols  /=  4) .and. (ncols /= 5) .and. (ncols /= 6) .and. (ncols /= 7)) then
          message(1) = 'Each line in the MaxwellExternalCurrent block must have'
          message(2) = 'four, five, six or or seven columns.'
          call messages_fatal(2, namespace=namespace)
        end if
 
        call parse_block_integer(blk, il - 1, 0, st%external_current_modus(il))
 
        if (st%external_current_modus(il) == external_current_parser) then
          ! parse formula string
          do idir = 1, space%dim
            call parse_block_string(blk, il - 1, idir, st%external_current_string(idir, il), convert_to_c=.true.)
          end do
        else if (st%external_current_modus(il) == external_current_td_function) then
          do idir = 1, space%dim
            call parse_block_string(blk, il - 1, idir, st%external_current_string(idir, il), convert_to_c=.true.)
            do ip = 1, mesh%np
              call mesh_r(mesh, ip, rr, coords = xx)
              call parse_expression(j_vector(idir), dummy(idir), st%dim, xx, rr, m_zero, &
                st%external_current_string(idir, il))
              st%external_current_amplitude(ip, idir, il) = j_vector(idir)
            end do
          end do
          call parse_block_float(blk, il-1, 4, omega, unit_one/units_inp%time)
          st%external_current_omega(il) = omega
          call parse_block_string(blk, il-1, 5, tdf_expression)
          call tdf_read(st%external_current_td_function(il), namespace, trim(tdf_expression), ierr)
          if (parse_block_cols(blk, il-1) > 6) then
            call parse_block_string(blk, il-1, 6, phase_expression)
            call tdf_read(st%external_current_td_phase(il), namespace, trim(phase_expression), ierr)
            if (ierr /= 0) then
              write(message(1),'(3A)') 'Error in the "', trim(tdf_expression), '" field defined in the TDExternalFields block:'
              write(message(2),'(3A)') 'Time-dependent phase function "', trim(phase_expression), '" not found.'
              call messages_warning(2, namespace=namespace)
            end if
          else
            call tdf_init(st%external_current_td_phase(il))
          end if
        end if
      end do
      call parse_block_end(blk)
    end if
 
    call profiling_out('EXTERNAL_CURRENT_INIT')
 
    pop_sub(external_current_init)
  end subroutine external_current_init
 
  !----------------------------------------------------------
  subroutine external_current_calculation(st, space, mesh, time, current)
    type(states_mxll_t), intent(in)    :: st
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    real(real64),        intent(in)    :: time
    real(real64),        intent(inout) :: current(:,:)
 
    integer :: ip, jn, idir
    real(real64)   :: xx(space%dim), rr, tt, j_vector(space%dim), dummy(space%dim), amp(space%dim)
    complex(real64)   :: exp_arg
    real(real64)   :: tmp_amp, phase
 
    push_sub(external_current_calculation)
 
    call profiling_in("EXTERNAL_CURRENT_CALC")
 
    current(:,:) = m_zero
    do jn = 1, st%external_current_number
      if (st%external_current_modus(jn) == external_current_parser) then
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, coords = xx)
          do idir = 1, space%dim
            tt = time
            call parse_expression(j_vector(idir), dummy(idir), space%dim, xx, rr, tt, &
              trim(st%external_current_string(idir,jn)))
          end do
          current(ip, :) = current(ip, :) + j_vector
        end do
 
      else if (st%external_current_modus(jn) == external_current_td_function) then
        exp_arg = st%external_current_omega(jn) * time + tdf(st%external_current_td_phase(jn),time)
        phase = real(exp(-m_zi*exp_arg), real64)
        tmp_amp = tdf(st%external_current_td_function(jn), time)
        do ip = 1, mesh%np
          amp = st%external_current_amplitude(ip, :, jn) * tmp_amp
          current(ip, :) = current(ip, :) + amp(:) * phase
        end do
      end if
    end do
 
    call profiling_out("EXTERNAL_CURRENT_CALC")
 
    pop_sub(external_current_calculation)
  end subroutine external_current_calculation
 
end module external_densities_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: