sll_m_utilities.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!  Authors :
!     CALVI project team
!
!  This code SeLaLib (for Semi-Lagrangian-Library)
!  is a parallel library for simulating the plasma turbulence
!  in a tokamak.
!
!  This software is governed by the CeCILL-B license
!  under French law and abiding by the rules of distribution
!  of free software.  You can  use, modify and redistribute
!  the software under the terms of the CeCILL-B license as
!  circulated by CEA, CNRS and INRIA at the following URL
!  "http://www.cecill.info".
!**************************************************************
 
module sll_m_utilities
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_errors.h"
#include "sll_working_precision.h"
 
   implicit none
 
   public :: &
      sll_p_byte_size, &
      sll_s_compute_bloc, &
      sll_s_compute_mesh_from_bloc, &
      sll_s_display_matrix_2d_integer, &
      sll_s_int2string, &
      sll_f_is_even, &
      sll_f_is_power_of_two, &
      sll_s_mpe_decomp1d, &
      sll_s_pfenvelope, &
      sll_o_display, &
      sll_o_factorial, &
      sll_s_new_file_id, &
      sll_f_query_environment, &
      sll_s_new_array_linspace
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  intrinsic :: selected_int_kind ! this line gives an error, why?
 
   ! Tentative implementation of a standard-compliant way to get the
   ! memory footprint of a variable. This is our yardstick...
   !
   ! selected_int_kind(r) returns the default integer scalar that is the kind
   ! type parameter value for an integer data type able to represent all
   ! integer values n in the range -10^(-r) < n < 10^r, where r is a scalar
   ! integer. If more than one is available, a kind with least decimal exponent
   ! range is chosen (and least kind value if several have least decimal
   ! exponent range). If no corresponding kind is availalble, the result is -1.
   ! (Metcalf & Reid. F90/95 2nd ed. p. 176).
   !
   ! We are, maybe dangerously, relying on the common practice of many compilers
   ! of using the kind values to indicate the number of bytes of storage
   ! occupied by a value. But this is not mandated by the standard. For
   ! instance a compiler that only has available 4-byte integers may still
   ! support kind values of 1, 2 and 4 to 'ease portability' from other
   ! platforms. The size of k1 will become our basic 'yardstick' to measure
   ! the size of a memory footprint of a variable. When we ask for the size
   ! of 'var', the answer will be given in terms of how many 'yardsticks'
   ! are needed to represent 'var'. The implications of this assumption
   ! need to be checked further.
 
   integer, parameter :: sll_p_byte_size = selected_int_kind(0)
 
   interface sll_o_display
      module procedure sll_s_display_matrix_2d_integer
      module procedure display_matrix_2d_real
      module procedure display_vector_integer
      module procedure display_vector_real
   end interface sll_o_display
 
   logical :: flag = .true.
 
   interface sll_o_factorial
      module procedure factorial_int32, factorial_int64
   end interface sll_o_factorial
 
contains
 
   function sll_f_is_power_of_two(n)
      sll_int64, intent(in) :: n
      logical               :: sll_f_is_power_of_two
 
      intrinsic :: not, iand
 
      if ((n > 0) .and. (0 .eq. (iand(n, (n - 1))))) then
         sll_f_is_power_of_two = .true.
      else
         sll_f_is_power_of_two = .false.
      end if
   end function sll_f_is_power_of_two
 
   function sll_f_is_even(n)
      sll_int32, intent(in) :: n
      logical               :: sll_f_is_even
 
      intrinsic :: modulo
 
      if (modulo(n, 2) .eq. 0) then
         sll_f_is_even = .true.
      else
         sll_f_is_even = .false.
      end if
   end function sll_f_is_even
 
   function factorial_int32(n) result(fac)
      sll_int32, intent(in) :: n
      sll_int64             :: fac
 
      sll_int64 :: acc
      sll_int64 :: i
 
      if (n < 0) then
         print *, 'ERROR, factorial_int32(): n < 0 or n > 20, if latter, ', &
            'function will overflow a 64-bit integer. n =  ', n
      end if
 
      acc = 1
      if (n >= 1) then
         do i = int(n, kind=i64), 1, -1
            acc = acc*i
         end do
      end if
      ! case n == 0 is already taken care of. No protection for negative input.
      fac = acc
   end function factorial_int32
 
   function factorial_int64(n) result(fac)
      sll_int64, intent(in) :: n
      sll_int64             :: fac
 
      sll_int64 :: acc
      sll_int64 :: i
 
      if ((n < 0) .or. (n > 20)) then
         print *, 'ERROR, factorial_int64(): either a negative n was passed: ', &
            'or n > 20, which will overflow a 64-bit integer. n = ', n
      end if
 
      acc = 1
      if (n >= 1) then
         do i = n, 1, -1
            acc = acc*i
         end do
      end if
      ! case n == 0 is already taken care of.
      fac = acc
   end function factorial_int64
 
   subroutine sll_s_int2string(istep, cstep)
      integer, intent(in) :: istep
      character(len=*), intent(out) :: cstep
 
      integer          :: l
      character(len=8) :: str_fmt
 
      l = len(cstep)
 
      if (istep >= 0 .and. istep < 10**l) then
         str_fmt = "(I0"//char(l + 48)//"."//char(l + 48)//")"
         write (cstep, str_fmt) istep
      else
         sll_warning('sll_s_int2string', 'index is negative or too big')
         print *, 'index =', istep, ' cstep length = ', l
         cstep = 'xxxx'
      end if
 
   end subroutine sll_s_int2string
 
   subroutine sll_s_new_file_id(file_id, error)
      sll_int32, intent(out) :: file_id   
      sll_int32, intent(out) :: error     
 
      logical :: lopen
 
      error = 1
 
      do 100 file_id = 20, 99
 
         inquire (unit=file_id, opened=lopen)
         if (lopen) then
            cycle
         else
            open (file_id, status='SCRATCH', err=100)
            close (file_id, status='DELETE', err=100)
            error = 0
            exit
         end if
 
100      continue
 
         !SLL_ASSERT(error == 0)
 
         end subroutine sll_s_new_file_id
 
         subroutine display_vector_real(array, real_format)
            sll_real64, dimension(:), intent(in) :: array
            character(len=*), intent(in) :: real_format
 
            character(len=20) :: display_format
            sll_int32         :: n, i
 
            n = size(array, 1)
 
            write (display_format, "('(''|''',i4,a,''' |'')')") n, real_format
 
            write (*, *)
            write (*, display_format) (array(i), i=1, n)
            write (*, *)
 
         end subroutine display_vector_real
 
         subroutine display_vector_integer(array, integer_format)
            sll_int32, dimension(:), intent(in) :: array
            character(len=*), intent(in) :: integer_format
 
            character(len=20) :: display_format
            sll_int32         :: n, i
 
            n = size(array)
 
            write (display_format, "('(''|''',i4,a,''' |'')')") n, integer_format
 
            write (*, *)
            write (*, display_format) (array(i), i=1, n)
            write (*, *)
 
         end subroutine display_vector_integer
 
         subroutine display_matrix_2d_real(array, real_format)
            sll_real64, dimension(:, :), intent(in) :: array
            character(len=*), intent(in) :: real_format
 
            character(len=20) :: display_format
            sll_int32         :: n1, n2, i, j
 
            n1 = size(array, 1)
            n2 = size(array, 2)
 
            write (display_format, "('(''|''',i4,a,''' |'')')") n2, real_format
 
            write (*, *)
            do i = 1, n1
               write (*, display_format) (array(i, j), j=1, n2)
            end do
            write (*, *)
 
         end subroutine display_matrix_2d_real
 
         subroutine sll_s_display_matrix_2d_integer(array, integer_format)
            sll_int32, dimension(:, :), intent(in) :: array
            character(len=*), intent(in) :: integer_format
 
            character(len=20) :: display_format
            sll_int32         :: n1, n2, i, j
 
            n1 = size(array, 1)
            n2 = size(array, 2)
 
            write (display_format, "('(''|''',i4,a,''' |'')')") n2, integer_format
 
            write (*, *)
            do i = 1, n1
               write (*, display_format) (array(i, j), j=1, n2)
            end do
            write (*, *)
 
         end subroutine sll_s_display_matrix_2d_integer
 
         subroutine initialize_file(data_file_id, thf_file_id)
            sll_int32, intent(out) :: data_file_id 
            sll_int32, intent(out) :: thf_file_id  
 
            character(len=*), parameter :: this_sub_name = "initialize_file"
            character(len=72)           :: filename
            integer                     :: IO_stat
            sll_int32                   :: error
 
            call get_command_argument(1, filename)
 
            call sll_s_new_file_id(data_file_id, error)
            open (data_file_id, file=trim(filename), iostat=io_stat)
            if (io_stat /= 0) sll_error(this_sub_name, "Miss argument file.nml")
 
            call sll_s_new_file_id(thf_file_id, error)
            open (thf_file_id, file="thf.dat", iostat=io_stat, position='append')
            if (io_stat /= 0) sll_error(this_sub_name, "Cannot open file thf.dat")
 
            rewind(thf_file_id)
            close (thf_file_id)
 
         end subroutine initialize_file
 
         subroutine time_history(file_id, desc, fformat, array)
            sll_int32, intent(in) :: file_id 
            character(3), intent(in) :: desc
            character(14), intent(in) :: fformat
            sll_real64, dimension(:), intent(in) :: array   
 
            if (desc(1:3) == "thf") then
               open (file_id, file="thf.dat", position='append')
               if (flag) then
                  rewind(file_id)
                  flag = .false.
               end if
               write (file_id, fformat) array
               close (file_id)
            else
               write (*, *) desc, " not recognized"
            end if
 
         end subroutine time_history
 
!------------------------------------------------------------------------
!------------------------------------------------------------------------
         subroutine sll_s_mpe_decomp1d(n, numprocs, myid, s, e)
            sll_int32, intent(in)  :: n
            sll_int32, intent(in)  :: numprocs
            sll_int32, intent(in)  :: myid
            sll_int32, intent(out) :: s
            sll_int32, intent(out) :: e
 
            sll_int32 :: nlocal
            sll_int32 :: deficit
 
            nlocal = n/numprocs
            s = myid*nlocal + 1
            deficit = mod(n, numprocs)
            s = s + min(myid, deficit)
            if (myid < deficit) then
               nlocal = nlocal + 1
            end if
            e = s + nlocal - 1
            if (e > n .or. myid == numprocs - 1) e = n
 
         end subroutine sll_s_mpe_decomp1d
 
         subroutine sll_s_pfenvelope(S, &
                                     t, &
                                     tflat, &
                                     tL, &
                                     tR, &
                                     twL, &
                                     twR, &
                                     t0, &
                                     turn_drive_off)
 
            sll_real64, intent(out) :: s
            sll_real64, intent(in)  :: t
            sll_real64, intent(in)  :: tflat
            sll_real64, intent(in)  :: tl
            sll_real64, intent(in)  :: tr
            sll_real64, intent(in)  :: twl
            sll_real64, intent(in)  :: twr
            sll_real64, intent(in)  :: t0
            logical, intent(in)  :: turn_drive_off
 
            sll_real64 :: epsilon
 
            ! The envelope function is defined such that it is zero at t0,
            ! rises to 1 smoothly, stay constant for tflat, and returns
            ! smoothly to zero.
            if (turn_drive_off) then
               epsilon = 0.5_f64*(tanh((t0 - tl)/twl) - tanh((t0 - tr)/twr))
               s = 0.5_f64*(tanh((t - tl)/twl) - tanh((t - tr)/twr)) - epsilon
               s = s/(1 - epsilon)
            else
               epsilon = 0.5_f64*(tanh((t0 - tl)/twl) + 1.0_f64)
               s = 0.5_f64*(tanh((t - tl)/twl) + 1.0_f64) - epsilon
               s = s/(1.0_f64 - epsilon)
            end if
            if (s < 0) then
               s = 0.0_f64
            end if
            s = s + 0.0_f64*tflat ! for use of unused
            return
 
         end subroutine sll_s_pfenvelope
 
         subroutine sll_s_compute_bloc(bloc_coord, bloc_index, N)
 
            sll_real64, intent(inout)  :: bloc_coord(2)
            sll_int32, intent(inout)  :: bloc_index(3)
            sll_int32, intent(in)     :: n
 
            sll_real64 :: a, b
            sll_int32  :: i1, i2, n_coarse, n_local, n_fine
 
            a = bloc_coord(1)
            b = bloc_coord(2)
 
            !case of uniform mesh with refined zone
            !we have a coarse mesh with N_coarse
            !N=i1+N_local*(i2-i1)+N_coarse-i2
            !N_fine=N_local*(i2-i1)
            !x(i1)=i1/N_coarse x(i1+N_fine)=i2/N_coarse
 
            if ((bloc_index(1) == 1) .and. (bloc_index(3) == 1)) then
 
               n_local = bloc_index(2)
               n_coarse = floor(real(n, f64)/(1._f64 + (b - a)*(real(n_local, f64) - 1._f64)))
               if (n_local /= 1) then
                  i2 = (n - n_coarse)/(n_local - 1)
               else
                  i2 = floor((b - a)*real(n_coarse, f64))
               end if
               n_coarse = n - i2*(n_local - 1)
               i1 = floor(a*real(n_coarse, f64))
               i2 = i2 + i1
               bloc_index(1) = i1
               n_fine = n_local*(i2 - i1)
               bloc_index(2) = n_fine
               bloc_index(3) = n - i1 - n_fine
               bloc_coord(1) = real(i1, f64)/real(n_coarse, f64)
               bloc_coord(2) = real(i2, f64)/real(n_coarse, f64)
 
               print *, '#uniform fine mesh would be:', n_coarse*n_local
               print *, '#N_coarse=', n_coarse
               print *, '#saving:', real(n, f64)/real(n_coarse*n_local, f64)
               print *, '#new x(i1),x(i1+N_fine)=', bloc_coord(1), bloc_coord(2)
               print *, '#error for x(i1),x(i1+N_fine)=', bloc_coord(1) - a, bloc_coord(2) - b
               print *, '#i1,i1+N_fine,N_fine,N=', i1, i1 + n_fine, n_fine, n
 
            else
 
               print *, 'case in sll_s_compute_bloc not implemented yet'
 
            end if
 
         end subroutine sll_s_compute_bloc
 
         subroutine sll_s_compute_mesh_from_bloc(bloc_coord, bloc_index, node_positions)
 
            sll_real64, intent(in)  :: bloc_coord(2)
            sll_int32, intent(in)  :: bloc_index(3)
            sll_real64, dimension(:), intent(out) :: node_positions
 
            sll_int32  :: i, i1, i2, n
            sll_real64 :: dx
 
            n = bloc_index(1) + bloc_index(2) + bloc_index(3)
            i1 = bloc_index(1)
            i2 = i1 + bloc_index(2)
            node_positions(1:n + 1) = -1._f64
            node_positions(1) = 0._f64
            node_positions(i1 + 1) = bloc_coord(1)
            node_positions(i2 + 1) = bloc_coord(2)
            node_positions(n + 1) = 1._f64
 
            !piecewise linear mapping (maybe enhanced like in complete mesh)
            if (bloc_index(1) .ne. 0) then
               dx = bloc_coord(1)/real(bloc_index(1), f64)
               do i = 2, bloc_index(1)
                  node_positions(i) = (real(i, f64) - 1._f64)*dx
               end do
            end if
            if (bloc_index(2) .ne. 0) then
               dx = (bloc_coord(2) - bloc_coord(1))/real(bloc_index(2), f64)
               do i = 2, bloc_index(2)
                  node_positions(i + i1) = bloc_coord(1) + (real(i, f64) - 1._f64)*dx
               end do
            end if
            if (bloc_index(3) .ne. 0) then
               dx = (1._f64 - bloc_coord(2))/real(bloc_index(3), f64)
               do i = 2, bloc_index(3)
                  node_positions(i + i2) = bloc_coord(2) + (real(i, f64) - 1._f64)*dx
               end do
            end if
         end subroutine sll_s_compute_mesh_from_bloc
 
         logical function sll_f_query_environment(env_variable, default_value)
 
            character(len=*), intent(in) :: env_variable
            logical, intent(in) :: default_value
 
            character(len=255) :: env_str
 
            sll_f_query_environment = default_value
 
            call get_environment_variable(env_variable, value=env_str)
 
            if (len_trim(env_str) > 0) then
               select case (trim(env_str))
               case ("1", "ON", "TRUE", "on", "true")
                  sll_f_query_environment = .true.
               case ("0", "OFF", "FALSE", "off", "false")
                  sll_f_query_environment = .false.
               end select
            end if
 
         end function sll_f_query_environment
 
         !-----------------------------------------------------------------------------
         !-----------------------------------------------------------------------------
         pure subroutine sll_s_new_array_linspace(array, vmin, vmax, endpoint, step)
 
            integer, parameter :: wp = f64
 
            real(wp), intent(out)           :: array(:)
            real(wp), intent(in)           :: vmin
            real(wp), intent(in)           :: vmax
            logical, intent(in), optional :: endpoint
            real(wp), intent(out), optional :: step
 
            integer  :: i
            integer  :: np
            integer  :: nc
            real(wp) :: a
            real(wp) :: b
 
            ! Read number of points from given array
            np = size(array)
 
            ! Calculate number of intervals in domain (cells) based on 'endpoint' flag
            ! By default, we assume that endpoint = .true.
            if (present(endpoint)) then
               nc = merge(np - 1, np, endpoint)
            else
               nc = np - 1
            end if
 
            ! Set first value to vmin in all cases
            array(1) = vmin
 
            ! Calculate internal array values using linear blending of vmin and vmax,
            ! in order to minimize round-off
            a = vmin/real(nc, wp)
            b = vmax/real(nc, wp)
            do i = 2, nc
               array(i) = a*real(nc + 1 - i, f64) + b*real(i - 1, f64)
            end do
 
            ! If endpoint=.true., set last value to vmax
            if (np == nc + 1) array(np) = vmax
 
            ! If required, return step size
            if (present(step)) then
               step = b - a
            end if
 
         end subroutine sll_s_new_array_linspace
         !-----------------------------------------------------------------------------
 
         end module sll_m_utilities