sll_m_advection_1d_CSL.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".
!**************************************************************
 
! in development; should be at least cubic splines
! attached with computation of characteristics
 
module sll_m_advection_1d_csl
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_advection_1d_base, only: &
      sll_c_advector_1d
 
   use sll_m_boundary_condition_descriptors, only: &
      sll_p_periodic, &
      sll_p_set_to_limit, &
      sll_p_user_defined
 
   use sll_m_characteristics_1d_base, only: &
      sll_f_process_outside_point_periodic, &
      sll_f_process_outside_point_set_to_limit, &
      sll_i_signature_process_outside_point_1d, &
      sll_c_characteristics_1d_base
 
   use sll_m_interpolators_1d_base, only: &
      sll_c_interpolator_1d
 
   implicit none
 
   public :: &
      sll_f_new_csl_1d_advector
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type, extends(sll_c_advector_1d) :: csl_1d_advector
 
      class(sll_c_interpolator_1d), pointer  :: interp
      class(sll_c_characteristics_1d_base), pointer  :: charac
      sll_real64, dimension(:), pointer :: eta_coords
      sll_real64, dimension(:), pointer :: eta_coords_unit
      sll_real64, dimension(:), pointer :: eta_coords_unit_back
      sll_real64, dimension(:), pointer :: charac_feet
      sll_real64, dimension(:), pointer :: charac_feet_i
      sll_real64, dimension(:), pointer :: buf1d
      sll_real64, dimension(:), pointer :: buf1d_out
      sll_int32 :: npts
      sll_int32 :: bc_type
      procedure(sll_i_signature_process_outside_point_1d), pointer, nopass    :: &
         process_outside_point
 
   contains
      procedure, pass(adv) :: initialize => &
         initialize_csl_1d_advector
      procedure, pass(adv) :: advect_1d => &
         csl_advect_1d
      procedure, pass(adv) :: advect_1d_constant => &
         csl_advect_1d_constant
      procedure, pass(adv) :: delete => delete_csl_1d_advector
   end type csl_1d_advector
 
contains
   function sll_f_new_csl_1d_advector( &
      interp, &
      charac, &
      Npts, &
      eta_min, &
      eta_max, &
      eta_coords, &
      bc_type, &
      process_outside_point) &
      result(adv)
      type(csl_1d_advector), pointer :: adv
      class(sll_c_interpolator_1d), pointer :: interp
      class(sll_c_characteristics_1d_base), pointer  :: charac
      sll_int32, intent(in) :: npts
      sll_real64, intent(in), optional :: eta_min
      sll_real64, intent(in), optional :: eta_max
      sll_real64, dimension(:), pointer, optional :: eta_coords
      sll_int32, intent(in), optional :: bc_type
      procedure(sll_i_signature_process_outside_point_1d), optional    :: &
         process_outside_point
      sll_int32 :: ierr
 
      sll_allocate(adv, ierr)
 
      call initialize_csl_1d_advector( &
         adv, &
         interp, &
         charac, &
         npts, &
         eta_min, &
         eta_max, &
         eta_coords, &
         bc_type, &
         process_outside_point)
 
   end function sll_f_new_csl_1d_advector
 
   subroutine initialize_csl_1d_advector( &
      adv, &
      interp, &
      charac, &
      Npts, &
      eta_min, &
      eta_max, &
      eta_coords, &
      bc_type, &
      process_outside_point)
      class(csl_1d_advector), intent(inout) :: adv
      class(sll_c_interpolator_1d), pointer :: interp
      class(sll_c_characteristics_1d_base), pointer  :: charac
      sll_int32, intent(in) :: npts
      sll_real64, intent(in), optional :: eta_min
      sll_real64, intent(in), optional :: eta_max
      sll_real64, dimension(:), pointer, optional :: eta_coords
      sll_int32, intent(in), optional :: bc_type
      procedure(sll_i_signature_process_outside_point_1d), optional    :: &
         process_outside_point
      sll_int32 :: ierr
      sll_int32 :: i
      sll_real64 :: delta_eta
 
      adv%Npts = npts
      adv%interp => interp
      adv%charac => charac
      sll_allocate(adv%eta_coords(npts), ierr)
      sll_allocate(adv%eta_coords_unit(npts), ierr)
      sll_allocate(adv%eta_coords_unit_back(npts), ierr)
      sll_allocate(adv%buf1d(npts), ierr)
      sll_allocate(adv%buf1d_out(npts), ierr)
 
      sll_allocate(adv%charac_feet(npts), ierr)
      sll_allocate(adv%charac_feet_i(npts), ierr)
 
      if (present(eta_min) .and. present(eta_max)) then
         if (present(eta_coords)) then
            print *, '#provide either eta_coords or eta_min and eta_max'
            print *, '#and not both in subroutine initialize_BSL_1d_advector'
            stop
         else
            delta_eta = (eta_max - eta_min)/real(npts - 1, f64)
            do i = 1, npts
               adv%eta_coords(i) = eta_min + real(i - 1, f64)*delta_eta
            end do
         end if
      else if (present(eta_coords)) then
         if (size(eta_coords) < npts) then
            print *, '#bad size for eta_coords in initialize_BSL_1d_advector'
            stop
         else
            adv%eta_coords(1:npts) = eta_coords(1:npts)
         end if
      else
         print *, '#Warning, we assume eta_min = 0._f64 eta_max = 1._f64'
         delta_eta = 1._f64/real(npts - 1, f64)
         do i = 1, npts
            adv%eta_coords(i) = real(i - 1, f64)*delta_eta
         end do
      end if
 
      if (present(process_outside_point)) then
         adv%process_outside_point => process_outside_point
         adv%bc_type = sll_p_user_defined
      else if (.not. (present(bc_type))) then
         print *, '#provide boundary condition'
         print *, '#bc_type or process_outside_point function'
         print *, '#in initialize_CSL_1d_advector'
         stop
      else
         adv%bc_type = bc_type
         select case (bc_type)
         case (sll_p_periodic)
            adv%process_outside_point => sll_f_process_outside_point_periodic
         case (sll_p_set_to_limit)
            adv%process_outside_point => sll_f_process_outside_point_set_to_limit
         case default
            print *, '#bad value of boundary condition'
            print *, '#in initialize_CSL_1d_advector'
            stop
         end select
      end if
 
      if ((present(process_outside_point)) .and. (present(bc_type))) then
         print *, '#provide either process_outside_point or bc_type'
         print *, '#and not both'
         print *, '#in initialize_CSL_1d_advector'
         stop
      end if
 
      adv%eta_coords_unit(1:npts) = &
         (adv%eta_coords(1:npts) - adv%eta_coords(1)) &
         /(adv%eta_coords(npts) - adv%eta_coords(1))
 
   end subroutine initialize_csl_1d_advector
 
   subroutine csl_advect_1d( &
      adv, &
      A, &
      dt, &
      input, &
      output)
      class(csl_1d_advector) :: adv
      sll_real64, dimension(:), intent(in) :: a
      sll_real64, intent(in) :: dt
      sll_real64, dimension(:), intent(in) :: input
      sll_real64, dimension(:), intent(out) :: output
      sll_int32 :: npts
      sll_real64 :: mean
      sll_int32 :: i
      sll_real64 :: eta_min
      sll_real64 :: eta_max
 
      npts = adv%Npts
 
      eta_min = adv%eta_coords(1)
      eta_max = adv%eta_coords(npts)
 
      call adv%charac%compute_characteristics( &
         a, &
         dt, &
         adv%eta_coords, &
         adv%charac_feet)
 
      !print *,'#eta_min=',eta_min,eta_max
      !do i=1,Npts
      !  print *,i,adv%charac_feet(i)
      !enddo
      !stop
      do i = 1, npts
         if ((adv%charac_feet(i) <= eta_min) .or. (adv%charac_feet(i) >= eta_max)) then
            adv%charac_feet_i(i) = &
               adv%process_outside_point(adv%charac_feet(i), eta_min, eta_max)
         else
            adv%charac_feet_i(i) = adv%charac_feet(i)
         end if
      end do
 
      adv%buf1d(1:npts - 1) = input(1:npts - 1)
 
!    do i=1,Npts
!      !x1 = adv%eta_coords(1)+real(i-1,f64)*(adv%eta_coords(Npts)-adv%eta_coords(1))/real(Npts-1,f64)
!      x1 = adv%eta_coords_unit(i)
!      adv%buf1d(i) = cos(2._f64*sll_p_pi*x1)
!    enddo
 
      call function_to_primitive_csl(adv%buf1d, adv%eta_coords_unit, npts - 1, mean)
 
      !adv%buf1d(1:Npts) = adv%buf1d(1:Npts)*(adv%eta_coords(Npts)-adv%eta_coords(1))
 
      print *, mean
      if (mean < 0.95) then
         print *, '#mean value=', mean
         !print *,'#charac=',adv%charac_feet(Npts)-adv%charac_feet(1)
         !do i=1,Npts-1
         !  print *,i,(adv%charac_feet(i+1)-adv%charac_feet(i))
         !enddo
 
         stop
      end if
      !print *,'#mean=',mean
!
!    do i=1,Npts
!      !print *,i,input(i),adv%buf1d(i),adv%eta_coords_unit(i),adv%eta_coords(i)
!      print *,adv%eta_coords_unit(i),cos(2._f64*sll_p_pi*adv%eta_coords_unit(i)),adv%buf1d(i), &
!      sin(2._f64*sll_p_pi*adv%eta_coords_unit(i))/(2._f64*sll_p_pi)
!      !-sin(adv%eta_coords(i)-adv%eta_coords(1))/(2._f64*sll_p_pi)
!    enddo
 
      !do i=1,Npts
      !  print *,i,adv%eta_coords(i),adv%charac_feet(i)
      !enddo
 
!    stop
 
!    call adv%interp%compute_interpolants( &
!      input, &
!      adv%eta1_coords, &
!      adv%Npts1, &
!      adv%eta2_coords, &
!      adv%Npts2 )
 
      call adv%interp%interpolate_array( &
         npts, &
         adv%buf1d, &
         adv%charac_feet_i, &
         adv%buf1d_out)
 
      !adv%buf1d_out(1:Npts) = adv%buf1d_out(1:Npts)/(adv%eta_coords(Npts)-adv%eta_coords(1))
 
      adv%eta_coords_unit_back(1:npts) = &
         (adv%charac_feet(1:npts) - adv%charac_feet(1)) &
         /(adv%charac_feet(npts) - adv%charac_feet(1))
 
      call primitive_to_function_csl( &
         adv%buf1d_out, &
         adv%eta_coords_unit, &
         !adv%eta_coords_unit, &
         adv%eta_coords_unit_back, &
         npts - 1, &
         mean)
 
      output(1:npts - 1) = adv%buf1d_out(1:npts - 1)
 
      output(npts) = output(1) !warning only valid in periodic case
 
   end subroutine csl_advect_1d
 
   subroutine csl_advect_1d_constant( &
      adv, &
      A, &
      dt, &
      input, &
      output)
      class(csl_1d_advector) :: adv
      sll_real64, intent(in) :: a
      sll_real64, intent(in) :: dt
      sll_real64, dimension(:), intent(in) :: input
      sll_real64, dimension(:), intent(out) :: output
      sll_real64, dimension(:), allocatable :: a1
      sll_int32 :: ierr
 
      !this version is not optimized
 
      sll_allocate(a1(adv%Npts), ierr)
 
      a1 = a
 
      call adv%charac%compute_characteristics( &
         a1, &
         dt, &
         adv%eta_coords, &
         adv%charac_feet)
 
!    call adv%interp%compute_interpolants( &
!      input, &
!      adv%eta1_coords, &
!      adv%Npts1, &
!      adv%eta2_coords, &
!      adv%Npts2 )
 
      call adv%interp%interpolate_array( &
         adv%Npts, &
         input, &
         adv%charac_feet, &
         output)
 
      sll_deallocate_array(a1, ierr)
 
   end subroutine csl_advect_1d_constant
 
   subroutine function_to_primitive_csl(f, node_positions, N, M)
      sll_real64, dimension(:), intent(inout) :: f
      sll_real64, dimension(:), intent(in) :: node_positions
      sll_int32, intent(in):: n
      sll_real64, intent(out)::m
      sll_int32::i
      sll_real64::tmp, tmp2
 
      !from f compute the mean
      m = 0._f64
      do i = 1, n
         m = m + f(i)*(node_positions(i + 1) - node_positions(i))
      end do
 
      f(1) = (f(1) - m)*(node_positions(2) - node_positions(1))
      tmp = f(1)
      f(1) = 0._f64
      do i = 2, n!+1
         f(i) = (f(i) - m)*(node_positions(i + 1) - node_positions(i))
         tmp2 = f(i)
         f(i) = f(i - 1) + tmp
         tmp = tmp2
      end do
      f(n + 1) = f(n) + tmp
 
      !print *,M,f(1),f(N+1)
 
   end subroutine function_to_primitive_csl
 
   subroutine primitive_to_function_csl( &
      f, &
      node_positions, &
      node_positions_back, &
      N, &
      M)
      sll_real64, dimension(:), intent(inout) :: f
      sll_real64, dimension(:), intent(in) :: node_positions
      sll_real64, dimension(:), intent(in) :: node_positions_back
      sll_int32, intent(in):: n
      sll_real64, intent(in)::m
      sll_int32::i
      sll_real64::tmp!,tmp2
 
      tmp = f(1)
      do i = 1, n - 1
         f(i) = f(i + 1) - f(i) + m*(node_positions_back(i + 1) - node_positions_back(i))
      end do
      f(n) = tmp - f(n) + m*(node_positions_back(n + 1) - node_positions_back(n))
 
      !from mean compute f
      do i = 1, n
         f(i) = f(i)/(node_positions(i + 1) - node_positions(i))
      end do
      !f(N+1) = f(1)
   end subroutine primitive_to_function_csl
 
   subroutine delete_csl_1d_advector(adv)
      class(csl_1d_advector), intent(inout) :: adv
      sll_int32 :: ierr
      sll_deallocate(adv%eta_coords, ierr)
      sll_deallocate(adv%eta_coords_unit, ierr)
      sll_deallocate(adv%eta_coords_unit_back, ierr)
      sll_deallocate(adv%charac_feet, ierr)
      sll_deallocate(adv%charac_feet_i, ierr)
      sll_deallocate(adv%buf1d, ierr)
      sll_deallocate(adv%buf1d_out, ierr)
   end subroutine delete_csl_1d_advector
 
end module sll_m_advection_1d_csl