79 integer,
public,
parameter :: &
93 integer :: mode = tdf_empty
94 real(real64) :: t0 =
m_zero
95 real(real64) :: tau0 =
m_zero
96 real(real64) :: tau1 =
m_zero
97 real(real64) :: a0 =
m_zero
98 real(real64) :: omega0 =
m_zero
99 real(real64) :: dt =
m_zero
100 real(real64) :: init_time =
m_zero
101 real(real64) :: final_time =
m_zero
103 integer :: nfreqs = 0
105 type(spline_t) :: amplitude
106 character(len=1024) :: expression
107 real(real64),
allocatable :: val(:)
108 real(real64),
allocatable :: valww(:)
109 type(fft_t) :: fft_handler
124 subroutine tdf_read(f, namespace, function_name, ierr)
125 type(tdf_t),
intent(inout) :: f
126 type(namespace_t),
intent(in) :: namespace
127 character(len=*),
intent(in) :: function_name
128 integer,
intent(out) :: ierr
131 integer :: nrows, i, function_type
132 character(len=1024) :: row_name, filename, function_expression
133 real(real64) :: a0, tau0, t0, tau1
208 if (
parse_block(namespace,
'TDFunctions', blk) /= 0)
then
215 row_loop:
do i = 1, nrows
217 if (trim(row_name) == trim(function_name))
then
225 select case (function_type)
271 type(
tdf_t),
intent(in) ::
f
279 type(
tdf_t),
intent(in) ::
f
286 real(real64)
pure function
tdf_dt(
f)
287 type(
tdf_t),
intent(in) ::
f
295 type(
tdf_t),
intent(inout) ::
f
310 type(
tdf_t),
intent(in) ::
f
322 type(
tdf_t),
intent(inout) ::
f
323 real(real64),
intent(in) :: a0, omega0
338 type(
tdf_t),
intent(inout) ::
f
339 real(real64),
intent(in) :: a0, omega0, t0, tau0
356 type(
tdf_t),
intent(inout) ::
f
357 real(real64),
intent(in) :: a0, omega0, t0, tau0
374 type(
tdf_t),
intent(inout) ::
f
375 real(real64),
intent(in) :: a0, omega0, t0, tau0, tau1
393 type(
tdf_t),
intent(inout) ::
f
394 character(len=*),
intent(in) :: expression
399 f%expression = trim(expression)
408 type(
tdf_t),
intent(inout) ::
f
409 character(len=*),
intent(in) :: filename
410 type(namespace_t),
intent(in) :: namespace
411 integer,
intent(out) :: ierr
413 integer :: iunit, lines, j
414 real(real64) :: dummy
415 real(real64),
allocatable :: t(:), am(:)
422 iunit = io_open(trim(filename), namespace, action=
'read', status=
'old')
425 call io_skip_header(iunit)
428 read(iunit, *, err=100,
end=100) dummy, dummy
433 call io_skip_header(iunit)
436 write(message(1),
'(a,a)')
"No data found in time-dependent function file ", trim(filename)
437 message(2) =
"Please check that the header properly starts by '#'"
438 call messages_fatal(2, namespace=namespace)
442 safe_allocate( t(1:lines))
443 safe_allocate(am(1:lines))
445 read(iunit, *) t(j), am(j)
450 f%final_time = t(lines)
452 call spline_init(
f%amplitude)
453 call spline_fit(lines, t, am,
f%amplitude, m_zero)
456 safe_deallocate_a(am)
465 type(
tdf_t),
intent(inout) ::
f
466 integer,
intent(in) :: niter
467 real(real64),
intent(in) :: dt
468 real(real64),
intent(in) :: omegamax
469 real(real64),
optional,
intent(in) :: initval
470 integer,
optional,
intent(in) :: rep
472 integer :: n(3), optimize_parity(3)
480 safe_allocate(
f%val(1:niter+1))
481 if (
present(initval))
then
489 f%final_time =
f%dt *
f%niter
491 if (omegamax > m_zero)
then
492 bigt =
f%final_time -
f%init_time
493 f%nfreqs = int(bigt * omegamax / (m_two * m_pi)) + 1
495 f%nfreqs =
f%niter/2+1
498 n(1:3) = (/
f%niter, 1, 1 /)
500 optimize_parity(1:3) = -1
501 call fft_init(
f%fft_handler, n, 1, fft_real, fftlib_fftw,
optimize, optimize_parity)
503 if (
present(rep))
then
506 safe_allocate(
f%valww(1:2*(
f%nfreqs-1)+1))
508 safe_deallocate_a(
f%val)
520 type(
tdf_t),
intent(in) ::
f
521 real(real64),
intent(inout) :: wgrid(:)
531 df = m_two * m_pi / (
f%final_time-
f%init_time)
536 message(1) =
"Illegal mode in tdf_fourier_grid."
537 call messages_fatal(1)
547 type(
tdf_t),
intent(inout) ::
f
550 complex(real64),
allocatable :: tmp(:)
554 safe_allocate(tmp(1:
f%niter/2+1))
558 f%val(1) = m_half*(
f%val(1)+
f%val(
f%niter+1))
559 f%val(
f%niter+1) =
f%val(1)
560 call dfft_forward(
f%fft_handler,
f%val(1:
f%niter), tmp)
561 tmp = tmp *
f%dt *
sqrt(m_one/(
f%final_time-
f%init_time))
563 safe_allocate(
f%valww(1:2*(
f%nfreqs-1)+1))
564 f%valww(1) = real(tmp(1), real64)
566 f%valww(j) = (
sqrt(m_two)) * real(tmp(j), real64)
568 do j =
f%nfreqs+1, 2*
f%nfreqs-1
569 f%valww(j) = - (
sqrt(m_two)) * aimag(tmp(j-
f%nfreqs+1))
572 safe_deallocate_a(tmp)
573 safe_deallocate_a(
f%val)
582 type(
tdf_t),
intent(inout) ::
f
585 complex(real64),
allocatable :: tmp(:)
589 safe_allocate(tmp(1:
f%niter/2+1))
593 tmp(j) = cmplx((
sqrt(m_two)/m_two)*
f%valww(j), -(
sqrt(m_two)/m_two)*
f%valww(j+
f%nfreqs-1), real64)
595 safe_allocate(
f%val(1:
f%niter+1))
596 call dfft_backward(
f%fft_handler, tmp,
f%val(1:
f%niter))
597 f%val(
f%niter+1) =
f%val(1)
598 f%val =
f%val *
f%niter *
sqrt(m_one/(
f%final_time-
f%init_time))
601 safe_deallocate_a(tmp)
602 safe_deallocate_a(
f%valww)
611 type(
tdf_t),
intent(inout) ::
f
618 s = sum(
f%valww(2:
f%nfreqs))
619 f%valww(2:
f%nfreqs) =
f%valww(2:
f%nfreqs) - s/
f%nfreqs
629 type(
tdf_t),
intent(inout) ::
f
633 assert(abs(
f%valww(1)) <= m_epsilon)
643 type(
tdf_t),
intent(inout) ::
f
644 real(real64),
intent(in) ::
values(:)
652 f%valww(1:2*
f%nfreqs-1) =
values(1:2*
f%nfreqs-1)
655 f%valww(2:2*
f%nfreqs-1) =
values(1:2*
f%nfreqs-2)
664 type(
tdf_t),
intent(inout) ::
f
665 integer,
intent(in) :: index
666 real(real64),
intent(in) :: val
676 f%valww(index+1) = val
685 type(
tdf_t),
intent(inout) ::
f
686 real(real64),
optional,
intent(in) :: fdotf
688 type(c_ptr) :: random_gen_pointer
690 real(real64) :: fdotf_, nrm
691 real(real64),
allocatable :: e(:)
701 message(1) =
"Illegal value for f%mode in tdf_set_random."
702 call messages_fatal(1)
704 safe_allocate(e(1:n))
706 if (mpi_world%is_root())
then
707 if (
present(fdotf))
then
713 call loct_ran_init(random_gen_pointer)
716 e(i) = loct_ran_gaussian(random_gen_pointer, m_one)
719 e =
sqrt(fdotf_) * e/ nrm
722 e(1:
f%nfreqs-1) = e(1:
f%nfreqs-1) - sum(e(1:
f%nfreqs-1))/(
f%nfreqs-1)
724 e =
sqrt(fdotf_) * e/ nrm
727 call loct_ran_end(random_gen_pointer)
730 call mpi_world%bcast(e(1), n, mpi_double_precision, 0)
741 type(
tdf_t),
intent(inout) :: f
742 integer,
optional,
intent(in) :: niter
743 real(real64),
optional,
intent(in) :: dt
744 real(real64),
optional,
intent(in) :: omegamax
748 real(real64),
allocatable :: val(:)
759 safe_allocate(val(1:niter+1))
765 assert(
present(niter))
767 assert(
present(omegamax))
770 safe_deallocate_a(val)
779 real(real64)
pure function
tdfi(
f, i) result(y)
780 type(
tdf_t),
intent(in) ::
f
781 integer,
intent(in) :: i
800 real(real64) function
tdft(
f, t) result(y)
801 type(
tdf_t),
intent(in) ::
f
802 real(real64),
intent(in) :: t
804 real(real64),
allocatable :: timearray(:), valarray(:)
805 real(real64) :: r, fre, fim, tcu
814 y =
f%a0 *
cos(
f%omega0*t)
818 r =
exp(-(t -
f%t0)**2 / (m_two*
f%tau0**2))
819 y =
f%a0 * r *
cos(
f%omega0 * t)
824 if (abs(t -
f%t0) <=
f%tau0)
then
825 r =
cos((m_pi / 2) * ((t - 2 *
f%tau0 -
f%t0) /
f%tau0))
827 y =
f%a0 * r *
cos(
f%omega0 * t)
830 if (t >
f%t0-
f%tau0/m_two-
f%tau1 .and. t <=
f%t0-
f%tau0 / m_two)
then
831 r = (t - (
f%t0 -
f%tau0/m_two -
f%tau1)) /
f%tau1
832 else if (t>
f%t0-
f%tau0/m_two .and. t <=
f%t0+
f%tau0 / m_two)
then
834 else if (t>
f%t0+
f%tau0/m_two .and. t <=
f%t0+
f%tau0 / m_two+
f%tau1)
then
835 r = (
f%t0 +
f%tau0/m_two +
f%tau1 - t) /
f%tau1
839 y =
f%a0 * r *
cos(
f%omega0 * t)
842 if (t >=
f%init_time .and. t <=
f%final_time)
then
843 y = spline_eval(
f%amplitude, t)
850 il = int(t/
f%dt)+1; iu = il+1
852 safe_allocate(timearray(1:4))
853 safe_allocate(valarray(1:4))
856 timearray = (/ m_zero,
f%dt, m_two*
f%dt, m_three*
f%dt /)
857 valarray = (/
f%val(1),
f%val(2),
f%val(3),
f%val(4) /)
858 elseif (il >=
f%niter)
then
859 timearray = (/ (
f%niter-3)*
f%dt, (
f%niter-2)*
f%dt, (
f%niter-1)*
f%dt,
f%niter*
f%dt /)
860 valarray = (/
f%val(
f%niter-2),
f%val(
f%niter-1),
f%val(
f%niter),
f%val(
f%niter+1) /)
862 timearray = (/ (il-2)*
f%dt, (il-1)*
f%dt, il*
f%dt, (il+1)*
f%dt /)
863 valarray(:) = (/
f%val(il-1),
f%val(il),
f%val(il+1),
f%val(il+2) /)
866 call interpolate(timearray, valarray, t, y)
868 safe_deallocate_a(valarray)
869 safe_deallocate_a(timearray)
872 tcu = units_from_atomic(units_inp%time, t)
873 call parse_expression(fre, fim,
't', tcu,
f%expression)
887 type(
tdf_t),
intent(inout) ::
f
893 call spline_end(
f%amplitude)
895 call fft_end(
f%fft_handler)
897 safe_deallocate_a(
f%valww)
898 call fft_end(
f%fft_handler)
901 safe_deallocate_a(
f%val)
910 type(
tdf_t),
intent(inout) :: fout
911 type(
tdf_t),
intent(in) :: fin
925 fout%omega0 = fin%omega0
926 fout%niter = fin%niter
927 fout%final_time = fin%final_time
928 fout%init_time = fin%init_time
929 fout%expression = fin%expression
930 fout%nfreqs = fin%nfreqs
932 fout%amplitude = fin%amplitude
935 safe_allocate(fout%val(1:fout%niter+1))
937 call fft_copy(fin%fft_handler, fout%fft_handler)
940 safe_allocate(fout%valww(1:2*fout%nfreqs-1))
941 fout%valww = fin%valww
942 call fft_copy(fin%fft_handler, fout%fft_handler)
953 real(real64),
intent(in) :: alpha
954 type(
tdf_t),
intent(inout) ::
f
964 f%valww = alpha*
f%valww
966 call spline_times(alpha,
f%amplitude, m_zero)
976 real(real64),
intent(in) :: omega
977 type(
tdf_t),
intent(inout) ::
f
987 do j = 1,
f%niter + 1
988 t =
f%init_time + (j-1)*
f%dt
989 f%val(j) =
f%val(j) *
cos(omega*t)
998 subroutine tdf_write(f, iunit, namespace)
999 type(
tdf_t),
intent(in) ::
f
1000 integer,
optional,
intent(in) :: iunit
1001 type(namespace_t),
optional,
intent(in) :: namespace
1007 select case (
f%mode)
1009 write(message(1),
'(6x,a)')
'Mode: continuous wave.'
1010 write(message(2),
'(6x,a,f10.4,a)')
'Amplitude: ',
f%a0,
' [a.u]'
1013 write(message(1),
'(6x,a)')
'Mode: Gaussian envelope.'
1014 write(message(2),
'(6x,a,f10.4,a)')
'Amplitude: ',
f%a0,
' [a.u]'
1015 write(message(3),
'(6x,a,f10.4,3a)')
'Width: ', units_from_atomic(units_out%time,
f%tau0), &
1016 ' [', trim(units_abbrev(units_out%time)),
']'
1017 write(message(4),
'(6x,a,f10.4,3a)')
'Middle t: ', units_from_atomic(units_out%time,
f%t0), &
1018 ' [', trim(units_abbrev(units_out%time)),
']'
1021 write(message(1),
'(6x,a)')
'Mode: cosinoidal envelope.'
1022 write(message(2),
'(6x,a,f10.4,a)')
'Amplitude: ',
f%a0,
' [a.u]'
1023 write(message(3),
'(6x,a,f10.4,3a)')
'Width: ', units_from_atomic(units_out%time,
f%tau0), &
1024 ' [', trim(units_abbrev(units_out%time)),
']'
1025 write(message(4),
'(6x,a,f10.4,3a)')
'Middle t: ', units_from_atomic(units_out%time,
f%t0), &
1026 ' [', trim(units_abbrev(units_out%time)),
']'
1029 write(message(1),
'(6x,a)')
'Mode: trapezoidal envelope.'
1030 write(message(2),
'(6x,a,f10.4,a)')
'Amplitude: ',
f%a0,
' [a.u]'
1031 write(message(3),
'(6x,a,f10.4,3a)')
'Width: ', units_from_atomic(units_out%time,
f%tau0), &
1032 ' [', trim(units_abbrev(units_out%time)),
']'
1033 write(message(4),
'(6x,a,f10.4,3a)')
'Middle t: ', units_from_atomic(units_out%time,
f%t0), &
1034 ' [', trim(units_abbrev(units_out%time)),
']'
1035 write(message(5),
'(6x,a,f10.4,3a)')
'Ramp time: ', units_from_atomic(units_out%time,
f%tau1), &
1036 ' [', trim(units_abbrev(units_out%time)),
']'
1039 write(message(1),
'(6x,a)')
'Mode: time-dependent function read from file.'
1042 write(message(1),
'(6x,a)')
'Mode: time-dependent function stored in a numerical array.'
1045 write(message(1),
'(6x,a)')
'Mode: time-dependent function parsed from the expression:'
1046 write(message(2),
'(6x,a)')
' f(t) = '//trim(
f%expression)
1050 if (abs(
f%omega0) > m_epsilon)
then
1052 write(message(n_msg),
'(6x,a,f10.4,3a)')
'Frequency: ', units_from_atomic(units_out%energy,
f%omega0), &
1053 ' [', trim(units_abbrev(units_out%energy)),
']'
1056 call messages_info(n_msg, iunit=iunit, namespace=namespace)
1080 assert(
f%mode == g%mode)
1082 select case (
f%mode)
1090 fg = fg +
f%val(2*i-2+1)*g%val(2*i-2+1) + m_four*
f%val(2*i-1+1)*g%val(2*i-1+1) +
f%val(2*i+1)*g%val(2*i+1)
1092 fg = fg *
f%dt / m_three
1094 if (mod(
f%niter, 2).eq.1)
then
1095 fg = fg + m_half * (
f%val(
f%niter)*g%val(
f%niter) +
f%val(
f%niter+1) * g%val(
f%niter+1)) *
f%dt
1099 fg = dot_product(
f%valww, g%valww)
1101 assert(abs(
f%valww(1)) <= m_epsilon)
1102 fg = dot_product(
f%valww, g%valww)
1104 do i = 1,
f%niter + 1
1106 fg = fg +
tdf(
f, i) *
tdf(g, i)
1122 real(real64) function
tdf_diff(
f, g) result (fg)
1123 type(
tdf_t),
intent(in) ::
f, g
1126 type(
tdf_t) :: fminusg
1131 assert(
f%mode == g%mode)
1135 select case (
f%mode)
1138 fminusg%val(i) = fminusg%val(i) - g%val(i)
1141 do i = 1, 2*(
f%nfreqs-1)+1
1142 fminusg%valww(i) = fminusg%valww(i) - g%valww(i)
1145 message(1) =
"Illegal value for f%mode in tdf_diff"
1146 call messages_fatal(1)
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
Fast Fourier Transform module. This module provides a single interface that works with different FFT ...
real(real64), parameter, public m_zero
This module is intended to contain "only mathematical" functions and procedures.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public tdf_numerical_to_fourier(f)
integer, parameter, public tdf_from_file
integer, parameter, public tdf_cosinoidal
subroutine, public tdf_init_cw(f, a0, omega0)
integer pure function, public tdf_nfreqs(f)
logical function, public tdf_is_empty(f)
integer, parameter, public tdf_numerical
subroutine tdf_set_numericalr1(f, index, val)
subroutine tdf_set_numericalr(f, values)
subroutine, public tdf_init_fromfile(f, filename, namespace, ierr)
integer pure function, public tdf_niter(f)
integer, parameter, public tdf_empty
subroutine, public tdf_init_cosinoidal(f, a0, omega0, t0, tau0)
subroutine, public tdf_init_gaussian(f, a0, omega0, t0, tau0)
real(real64) pure function tdfi(f, i)
subroutine, public tdf_numerical_to_zerofourier(f)
subroutine, public tdf_end(f)
integer, parameter, public tdf_zero_fourier
subroutine, public tdf_write(f, iunit, namespace)
subroutine, public tdf_copy(fout, fin)
real(real64) function, public tdf_diff(f, g)
subroutine, public tdf_init_trapezoidal(f, a0, omega0, t0, tau0, tau1)
integer, parameter, public tdf_cw
real(real64) function tdft(f, t)
subroutine, public tdf_cosine_multiply(omega, f)
subroutine, public tdf_set_random(f, fdotf)
subroutine, public tdf_fourier_grid(f, wgrid)
integer, parameter, public tdf_trapezoidal
subroutine, public tdf_init(f)
subroutine, public tdf_init_numerical(f, niter, dt, omegamax, initval, rep)
subroutine, public tdf_to_numerical(f, niter, dt, omegamax)
subroutine, public tdf_fourier_to_numerical(f)
subroutine, public tdf_zerofourier_to_numerical(f)
real(real64) function, public tdf_dot_product(f, g)
integer, parameter, public tdf_gaussian
integer, parameter, public tdf_from_expr
subroutine, public tdf_read(f, namespace, function_name, ierr)
This function initializes "f" from the TDFunctions block.
integer, parameter, public tdf_fourier_series
subroutine, public tdf_scalar_multiply(alpha, f)
real(real64) pure function, public tdf_dt(f)
subroutine, public tdf_init_fromexpr(f, expression)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
static double f(double w, void *p)
real(real64) function values(xx)