24 #include "sll_assert.h"
25 #include "sll_memory.h"
26 #include "sll_working_precision.h"
59 logical :: feet_inside
63 procedure, pass(charac) :: compute_characteristics => &
74 process_outside_point, &
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
88 sll_int32,
intent(in),
optional :: maxiter
89 sll_real64,
intent(in),
optional :: tol
90 logical,
intent(in),
optional :: feet_inside
93 sll_allocate(charac, ierr)
101 process_outside_point, &
114 process_outside_point, &
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
125 process_outside_point
127 sll_int32,
intent(in),
optional :: maxiter
128 sll_real64,
intent(in),
optional :: tol
129 logical,
intent(in),
optional :: feet_inside
133 if (
present(eta_min))
then
134 charac%eta_min = eta_min
136 charac%eta_min = 0._f64
138 if (
present(eta_max))
then
139 charac%eta_max = eta_max
141 charac%eta_max = 1._f64
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'
152 select case (bc_type)
153 case (sll_p_periodic)
155 case (sll_p_set_to_limit)
158 print *,
'#bad value of boundary condition'
159 print *,
'#in initialize_trapezoid_1d_charac'
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'
171 charac%A_interp => a_interp
173 if (
present(maxiter))
then
174 charac%maxiter = maxiter
176 charac%maxiter = 1000
179 if (
present(tol))
then
182 charac%tol = 1.e-12_f64
185 if (
present(feet_inside))
then
186 charac%feet_inside = feet_inside
188 charac%feet_inside = .true.
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
211 sll_real64 :: eta_min
212 sll_real64 :: eta_max
215 eta_min = charac%eta_min
216 eta_max = charac%eta_max
218 sll_assert(
size(a) >= npts)
219 sll_assert(
size(input) >= npts)
220 sll_assert(
size(output) >= npts)
227 call charac%A_interp%compute_interpolants(a)
237 x2 = input(j) - dt*a(j)
240 do while (iter < charac%maxiter .and. abs(x2_old - x2) > charac%tol)
243 if ((x2 <= eta_min) .or. (x2 >= eta_max))
then
244 x2_i = charac%process_outside_point(x2, eta_min, eta_max)
248 x2 = input(j) - 0.5_f64*dt*(charac%A_interp%interpolate_from_interpolant_value(x2_i) + a(j))
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)
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)
Describe different boundary conditions.
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.