sll_m_characteristics_2d_explicit_euler.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_characteristics_2d_explicit_euler
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_boundary_condition_descriptors, only: &
      sll_p_periodic, &
      sll_p_set_to_limit
 
   use sll_m_characteristics_2d_base, only: &
      sll_f_process_outside_point_periodic, &
      sll_f_process_outside_point_set_to_limit, &
      sll_i_signature_process_outside_point, &
      sll_c_characteristics_2d_base
 
   implicit none
 
   public :: &
      sll_f_new_explicit_euler_2d_charac, &
      sll_t_charac_2d_explicit_euler
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type, extends(sll_c_characteristics_2d_base) :: sll_t_charac_2d_explicit_euler
      sll_int32                               :: npts1
      sll_int32                               :: npts2
      sll_real64                              :: eta1_min
      sll_real64                              :: eta1_max
      sll_real64                              :: eta2_min
      sll_real64                              :: eta2_max
      procedure(sll_i_signature_process_outside_point), pointer, nopass    :: &
         process_outside_point1
      procedure(sll_i_signature_process_outside_point), pointer, nopass    :: &
         process_outside_point2
 
   contains
      procedure, pass(charac) :: init => initialize_explicit_euler_2d_charac
      procedure, pass(charac) :: compute_characteristics => &
         compute_explicit_euler_2d_charac
   end type sll_t_charac_2d_explicit_euler
 
contains
   function sll_f_new_explicit_euler_2d_charac( &
      Npts1, &
      Npts2, &
      bc_type_1, &
      bc_type_2, &
      eta1_min, &
      eta1_max, &
      eta2_min, &
      eta2_max, &
      process_outside_point1, &
      process_outside_point2) &
      result(charac)
 
      type(sll_t_charac_2d_explicit_euler), pointer :: charac
      sll_int32, intent(in) :: npts1
      sll_int32, intent(in) :: npts2
      sll_int32, intent(in), optional :: bc_type_1
      sll_int32, intent(in), optional :: bc_type_2
      sll_real64, intent(in), optional  :: eta1_min
      sll_real64, intent(in), optional  :: eta1_max
      sll_real64, intent(in), optional  :: eta2_min
      sll_real64, intent(in), optional  :: eta2_max
      procedure(sll_i_signature_process_outside_point), optional    :: &
         process_outside_point1
      procedure(sll_i_signature_process_outside_point), optional    :: &
         process_outside_point2
      sll_int32 :: ierr
 
      sll_allocate(charac, ierr)
      call initialize_explicit_euler_2d_charac( &
         charac, &
         npts1, &
         npts2, &
         bc_type_1, &
         bc_type_2, &
         eta1_min, &
         eta1_max, &
         eta2_min, &
         eta2_max, &
         process_outside_point1, &
         process_outside_point2)
 
   end function sll_f_new_explicit_euler_2d_charac
 
   subroutine initialize_explicit_euler_2d_charac( &
      charac, &
      Npts1, &
      Npts2, &
      bc_type_1, &
      bc_type_2, &
      eta1_min, &
      eta1_max, &
      eta2_min, &
      eta2_max, &
      process_outside_point1, &
      process_outside_point2)
 
      class(sll_t_charac_2d_explicit_euler) :: charac
      sll_int32, intent(in) :: npts1
      sll_int32, intent(in) :: npts2
      sll_int32, intent(in), optional :: bc_type_1
      sll_int32, intent(in), optional :: bc_type_2
      sll_real64, intent(in), optional  :: eta1_min
      sll_real64, intent(in), optional  :: eta1_max
      sll_real64, intent(in), optional  :: eta2_min
      sll_real64, intent(in), optional  :: eta2_max
      procedure(sll_i_signature_process_outside_point), optional    :: &
         process_outside_point1
      procedure(sll_i_signature_process_outside_point), optional    :: &
         process_outside_point2
 
      charac%Npts1 = npts1
      charac%Npts2 = npts2
 
      if (present(eta1_min)) then
         charac%eta1_min = eta1_min
      else
         charac%eta1_min = 0._f64
      end if
      if (present(eta1_max)) then
         charac%eta1_max = eta1_max
      else
         charac%eta1_max = 1._f64
      end if
      if (present(eta2_min)) then
         charac%eta2_min = eta2_min
      else
         charac%eta2_min = 0._f64
      end if
 
      if (present(eta2_max)) then
         charac%eta2_max = eta2_max
      else
         charac%eta2_max = 1._f64
      end if
 
      !charac%process_outside_point1 => process_outside_point1
      !charac%process_outside_point2 => process_outside_point2
 
      if (present(process_outside_point1)) then
         charac%process_outside_point1 => process_outside_point1
      else if (.not. (present(bc_type_1))) then
         print *, '#provide boundary condition'
         print *, '#bc_type_1 or process_outside_point1 function'
         print *, '#in initialize_explicit_euler_2d_charac'
         stop
      else
         select case (bc_type_1)
         case (sll_p_periodic)
            charac%process_outside_point1 => sll_f_process_outside_point_periodic
         case (sll_p_set_to_limit)
            charac%process_outside_point1 => sll_f_process_outside_point_set_to_limit
         case default
            print *, '#bad value of boundary condition'
            print *, '#in initialize_sll_t_characteristics_2d_explicit_euler_2d'
            stop
         end select
      end if
 
      if ((present(process_outside_point1)) .and. (present(bc_type_1))) then
         print *, '#provide either process_outside_point1 or bc_type_1'
         print *, '#and not both'
         print *, '#in initialize_sll_t_characteristics_2d_explicit_euler_2d'
         stop
      end if
 
      if (present(process_outside_point2)) then
         charac%process_outside_point2 => process_outside_point2
      else if (.not. (present(bc_type_2))) then
         print *, '#provide boundary condition'
         print *, '#bc_type_2 or process_outside_point1 function'
         stop
      else
         select case (bc_type_2)
         case (sll_p_periodic)
            charac%process_outside_point2 => sll_f_process_outside_point_periodic
         case (sll_p_set_to_limit)
            charac%process_outside_point2 => sll_f_process_outside_point_set_to_limit
         case default
            print *, '#bad value of boundary condition'
            print *, '#in initialize_sll_t_characteristics_2d_explicit_euler_2d'
            stop
         end select
      end if
 
      if ((present(process_outside_point2)) .and. (present(bc_type_2))) then
         print *, '#provide either process_outside_point2 or bc_type_2'
         print *, '#and not both'
         print *, '#in initialize_sll_t_characteristics_2d_explicit_euler_2d'
         stop
      end if
 
   end subroutine initialize_explicit_euler_2d_charac
 
   subroutine compute_explicit_euler_2d_charac( &
      charac, &
      A1, &
      A2, &
      dt, &
      input1, &
      input2, &
      output1, &
      output2)
 
      class(sll_t_charac_2d_explicit_euler) :: charac
      sll_real64, dimension(:, :), intent(in) :: a1
      sll_real64, dimension(:, :), intent(in) :: a2
      sll_real64, intent(in) :: dt
      sll_real64, dimension(:), intent(in) ::  input1
      sll_real64, dimension(:), intent(in) ::  input2
      sll_real64, dimension(:, :), intent(out) :: output1
      sll_real64, dimension(:, :), intent(out) :: output2
      sll_int32 :: i
      sll_int32 :: j
      sll_int32 :: npts1
      sll_int32 :: npts2
      sll_real64 :: eta1_min
      sll_real64 :: eta1_max
      sll_real64 :: eta2_min
      sll_real64 :: eta2_max
 
      npts1 = charac%Npts1
      npts2 = charac%Npts2
      eta1_min = charac%eta1_min
      eta1_max = charac%eta1_max
      eta2_min = charac%eta2_min
      eta2_max = charac%eta2_max
 
      sll_assert(size(a1, 1) >= charac%Npts1)
      sll_assert(size(a1, 2) >= charac%Npts2)
      sll_assert(size(a2, 1) >= charac%Npts1)
      sll_assert(size(a2, 2) >= charac%Npts2)
      sll_assert(size(input1) >= charac%Npts1)
      sll_assert(size(input2) >= charac%Npts2)
      sll_assert(size(output1, 1) >= charac%Npts1)
      sll_assert(size(output1, 2) >= charac%Npts2)
      sll_assert(size(output2, 1) >= charac%Npts1)
      sll_assert(size(output2, 2) >= charac%Npts2)
 
      do j = 1, npts2
         do i = 1, npts1
            output1(i, j) = input1(i) - dt*a1(i, j)
            if ((output1(i, j) <= eta1_min) .or. (output1(i, j) >= eta1_max)) then
               output1(i, j) = charac%process_outside_point1(output1(i, j), eta1_min, eta1_max)
            end if
            output2(i, j) = input2(j) - dt*a2(i, j)
            if ((output2(i, j) <= eta2_min) .or. (output2(i, j) >= eta2_max)) then
               output2(i, j) = charac%process_outside_point2(output2(i, j), eta2_min, eta2_max)
            end if
         end do
      end do
 
   end subroutine compute_explicit_euler_2d_charac
 
end module sll_m_characteristics_2d_explicit_euler