Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_characteristics_1d_trapezoid.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
18 !http://en.wikipedia.org/wiki/Trapezoidal_rule_%28differential_equations%29
19 
23 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24 #include "sll_assert.h"
25 #include "sll_memory.h"
26 #include "sll_working_precision.h"
27 
29  sll_p_periodic, &
30  sll_p_set_to_limit
31 
37 
38  use sll_m_interpolators_1d_base, only: &
40 
41  implicit none
42 
43  public :: &
46 
47  private
48 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
49 
51  sll_int32 :: npts
52  sll_real64 :: eta_min
53  sll_real64 :: eta_max
54  procedure(sll_i_signature_process_outside_point_1d), pointer, nopass :: &
55  process_outside_point
56  class(sll_c_interpolator_1d), pointer :: a_interp
57  sll_int32 :: maxiter
58  sll_real64 :: tol
59  logical :: feet_inside
60 
61  contains
62  procedure, pass(charac) :: init => initialize_trapezoid_1d_charac
63  procedure, pass(charac) :: compute_characteristics => &
66 
67 contains
69  Npts, &
70  A_interp, &
71  bc_type, &
72  eta_min, &
73  eta_max, &
74  process_outside_point, &
75  maxiter, &
76  tol, &
77  feet_inside) &
78  result(charac)
79 
80  type(sll_t_trapezoid_1d_charac), pointer :: charac
81  sll_int32, intent(in) :: npts
82  sll_int32, intent(in), optional :: bc_type
83  sll_real64, intent(in), optional :: eta_min
84  sll_real64, intent(in), optional :: eta_max
85  procedure(sll_i_signature_process_outside_point_1d), optional :: &
86  process_outside_point
87  class(sll_c_interpolator_1d), target :: a_interp
88  sll_int32, intent(in), optional :: maxiter
89  sll_real64, intent(in), optional :: tol
90  logical, intent(in), optional :: feet_inside
91  sll_int32 :: ierr
92 
93  sll_allocate(charac, ierr)
95  charac, &
96  npts, &
97  a_interp, &
98  bc_type, &
99  eta_min, &
100  eta_max, &
101  process_outside_point, &
102  maxiter, &
103  tol, &
104  feet_inside)
105 
106  end function sll_f_new_trapezoid_1d_charac
108  charac, &
109  Npts, &
110  A_interp, &
111  bc_type, &
112  eta_min, &
113  eta_max, &
114  process_outside_point, &
115  maxiter, &
116  tol, &
117  feet_inside)
118 
119  class(sll_t_trapezoid_1d_charac) :: charac
120  sll_int32, intent(in) :: npts
121  sll_int32, intent(in), optional :: bc_type
122  sll_real64, intent(in), optional :: eta_min
123  sll_real64, intent(in), optional :: eta_max
124  procedure(sll_i_signature_process_outside_point_1d), optional :: &
125  process_outside_point
126  class(sll_c_interpolator_1d), target :: A_interp
127  sll_int32, intent(in), optional :: maxiter
128  sll_real64, intent(in), optional :: tol
129  logical, intent(in), optional :: feet_inside
130 
131  charac%Npts = npts
132 
133  if (present(eta_min)) then
134  charac%eta_min = eta_min
135  else
136  charac%eta_min = 0._f64
137  end if
138  if (present(eta_max)) then
139  charac%eta_max = eta_max
140  else
141  charac%eta_max = 1._f64
142  end if
143 
144  if (present(process_outside_point)) then
145  charac%process_outside_point => process_outside_point
146  else if (.not. (present(bc_type))) then
147  print *, '#provide boundary condition'
148  print *, '#bc_type or process_outside_point function'
149  print *, '#in initialize_trapezoid_1d_charac'
150  stop
151  else
152  select case (bc_type)
153  case (sll_p_periodic)
154  charac%process_outside_point => sll_f_process_outside_point_periodic
155  case (sll_p_set_to_limit)
156  charac%process_outside_point => sll_f_process_outside_point_set_to_limit
157  case default
158  print *, '#bad value of boundary condition'
159  print *, '#in initialize_trapezoid_1d_charac'
160  stop
161  end select
162  end if
163 
164  if ((present(process_outside_point)) .and. (present(bc_type))) then
165  print *, '#provide either process_outside_point or bc_type'
166  print *, '#and not both'
167  print *, '#in initialize_trapezoid_1d_charac'
168  stop
169  end if
170 
171  charac%A_interp => a_interp
172 
173  if (present(maxiter)) then
174  charac%maxiter = maxiter
175  else
176  charac%maxiter = 1000
177  end if
178 
179  if (present(tol)) then
180  charac%tol = tol
181  else
182  charac%tol = 1.e-12_f64
183  end if
184 
185  if (present(feet_inside)) then
186  charac%feet_inside = feet_inside
187  else
188  charac%feet_inside = .true.
189  end if
190 
191  end subroutine initialize_trapezoid_1d_charac
192 
194  charac, &
195  A, &
196  dt, &
197  input, &
198  output)
199 
200  class(sll_t_trapezoid_1d_charac) :: charac
201  sll_real64, dimension(:), intent(in) :: a
202  sll_real64, intent(in) :: dt
203  sll_real64, dimension(:), intent(in) :: input
204  sll_real64, dimension(:), intent(out) :: output
205  sll_int32 :: j
206  sll_real64 :: x2
207  sll_real64 :: x2_old
208  sll_real64 :: x2_i !i for inside, so that interpolation is possible
209  sll_int32 :: iter
210  sll_int32 :: npts
211  sll_real64 :: eta_min
212  sll_real64 :: eta_max
213 
214  npts = charac%Npts
215  eta_min = charac%eta_min
216  eta_max = charac%eta_max
217 
218  sll_assert(size(a) >= npts)
219  sll_assert(size(input) >= npts)
220  sll_assert(size(output) >= npts)
221 
222  ! [YG: 8 Aug 2016]
223  ! Optional arguments 'eta_coords=input' and 'size_eta_coords=Npts'
224  ! not passed because:
225  ! 1. most interpolators do not implement such an option
226  ! 2. 'input' is (usually) the same mesh used to initialize the interpolator
227  call charac%A_interp%compute_interpolants(a)
228 ! call charac%A_interp%compute_interpolants( &
229 ! A, &
230 ! input, &
231 ! Npts)
232 
233  do j = 1, npts
234  !We start from Y(t_{n+1})=y_j
235  !and look for Y(t_n) = Yn
236  !Yn = y_j-(A(y_j)+A(Yn))*dt/2
237  x2 = input(j) - dt*a(j)
238  x2_old = 0._f64
239  iter = 0
240  do while (iter < charac%maxiter .and. abs(x2_old - x2) > charac%tol)
241  x2_old = x2
242  x2_i = x2
243  if ((x2 <= eta_min) .or. (x2 >= eta_max)) then
244  x2_i = charac%process_outside_point(x2, eta_min, eta_max)
245  else
246  x2_i = x2
247  end if
248  x2 = input(j) - 0.5_f64*dt*(charac%A_interp%interpolate_from_interpolant_value(x2_i) + a(j))
249  iter = iter + 1
250  !if(j==1)then
251  ! print *,'#x2=',x2
252  !endif
253  end do
254  if (iter == charac%maxiter .and. abs(x2_old - x2) > charac%tol) then
255  print *, '#not enough iterations for compute_trapezoid_1d_charac', &
256  iter, abs(x2_old - x2)
257  !stop
258  end if
259  if (((x2 <= eta_min) .or. (x2 >= eta_max)) .and. (charac%feet_inside)) then
260  x2 = charac%process_outside_point(x2, eta_min, eta_max)
261  end if
262  output(j) = x2
263  end do
264  !print *,'input=',input(1),input(Npts)
265  !print *,'output=',output(1),output(Npts)
266  end subroutine compute_trapezoid_1d_charac
267 
Abstract class for characteristic derived type.
function, public sll_f_process_outside_point_set_to_limit(eta, eta_min, eta_max)
function, public sll_f_process_outside_point_periodic(eta, eta_min, eta_max)
computes the characteristic with trapezoidal rule
subroutine initialize_trapezoid_1d_charac(charac, Npts, A_interp, bc_type, eta_min, eta_max, process_outside_point, maxiter, tol, feet_inside)
subroutine compute_trapezoid_1d_charac(charac, A, dt, input, output)
type(sll_t_trapezoid_1d_charac) function, pointer, public sll_f_new_trapezoid_1d_charac(Npts, A_interp, bc_type, eta_min, eta_max, process_outside_point, maxiter, tol, feet_inside)
Module for 1D interpolation and reconstruction.
Abstract class for 1D interpolation and reconstruction.
    Report Typos and Errors