diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 4c6ce37..0f6f8f0 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -1,535 +1,757 @@ -! - subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & - psi,qsaf,bcovar_phi,Phi_0,avR2nabpsi2, & - dens_e,temp_e,temp_i,anu_e,anu_i, & - bpsi_over_bphi,parcur_over_b0,Phi_m) -! -! All units are Gaussian if not otherwise mentioned -! -! Input :: -! isw_Phi_m - (integer) switch for the potential perturbation in the parallel current: -! 0 - set it to zero in the current, 1 - use it as is in the current -! mpol - (integer) poloidal mode number -! ntor - (integer) toroidal mode number -! npoi - (integer) number of radial grid points -! am_i - (double precision) ion mass number (mass in units of proton mass) -! Z_i - (double precision) ion charge number (charge in units of elementary charge) -! Rtor - (double precision) reference major radius -! psi(npoi) - (double precision) unperturbed poloidal flux $\psi$ -! qsaf(npoi) - (double precision) safety factor -! bcovar_phi(npoi) - (double precision) co-variant unperturbed toroidal magnetic field component -! Phi_0(npoi) - (double precision) equilibrium electrostatic potential -! avR2nabpsi2(npoi) - (double precision) average over the polodal angle of $R^2 |\nabla \psi|^2$ -! dens_e(npoi) - (double precision) electron density -! temp_e(npoi) - (double precision) electron temperature -! temp_i(npoi) - (double precision) ion temperature -! anu_e(npoi) - (double precision) electron collision frequency -! anu_i(npoi) - (double precision) ion collision frequency -! bpsi_over_bphi(npoi) - (double complex) psi-component of the perturbation magnetic field divided by -! contra-variant component of the unperturbed magnetic field -! Output: -! parcur_over_b0(npoi) - (double complex) parallel current density divided by uperturbed magnetic field module -! Phi_m(npoi) - (double complex) perturbation of the electrostatic potential -! +module mephit_flr2 + use iso_fortran_env, only: dp => real64 + implicit none -! - integer, parameter :: mnmax=3 - double precision, parameter :: pi=3.14159265358979d0 - double precision, parameter :: c=2.9979d10 - double precision, parameter :: e_charge=4.8032d-10 - double precision, parameter :: e_mass=9.1094d-28 - double precision, parameter :: p_mass=1.6726d-24 - double precision, parameter :: ev=1.6022d-12 - double complex, parameter :: imun=(0.d0,1.d0) -! - integer :: isw_Phi_m,mpol,ntor,npoi,i - double precision :: am_i,z_i,Rtor - double precision :: e_e,e_i,omega_E,x_1,x_2,v_T,rho2_factor - double precision :: switch_flr_e,switch_flr_i - double precision :: switch_cur_e,switch_cur_i - double complex :: factor_of_Phi - double complex :: F_m,Gtor_m,Htor - double complex :: F_me,Gtor_me,Htore - double complex :: F_mi,Gtor_mi,Htori - double precision, dimension(npoi) :: psi,qsaf,bcovar_phi,Phi_0,avR2nabpsi2, & - dens_e,temp_e,temp_i,anu_e,anu_i - double complex, dimension(npoi) :: bpsi_over_bphi,parcur_over_b0,Phi_m -! - double complex, dimension(0:mnmax,0:mnmax) :: symbI - double precision, dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i - double complex, dimension(:), allocatable :: a2_in,a2_out,a0 - double complex, dimension(:), allocatable :: b2_in,b2_out,b0 - double complex, dimension(:), allocatable :: c2_in,c2_out,c0 - double complex, dimension(:), allocatable :: d2_in,d2_out,d0 -! - allocate(dens_i(npoi),dPhi_0_dpsi(npoi),derpar(npoi)) - allocate(A1e(npoi),A2e(npoi),A1i(npoi),A2i(npoi)) - allocate(a2_in(npoi),a2_out(npoi),a0(npoi)) - allocate(b2_in(npoi),b2_out(npoi),b0(npoi)) - allocate(c2_in(npoi),c2_out(npoi),c0(npoi)) - allocate(d2_in(npoi),d2_out(npoi),d0(npoi)) -! -! For playing around: -! switch for FLR effects in electrons (1 - on, 0 - off): - switch_flr_e=1.d0 -! switch for FLR effects in ions (1 - on, 0 - off): - switch_flr_i=1.d0 -! switch for the electron contribution to the parallel current (1 - on, 0 - off): - switch_cur_e=1.d0 -! switch for the ion contribution to the parallel current (1 - on, 0 - off): - switch_cur_i=1.d0 -! -! electron and ion charges and ion density: - e_e=-e_charge - e_i=e_charge*z_i - dens_i=dens_e/z_i -! - call first_deriv(npoi,psi,Phi_0,dPhi_0_dpsi) -! -! -! thermodynamic forces for electrons: -! - call first_deriv(npoi,psi,dens_e,derpar) -! - A1e=derpar/dens_e+e_e*dPhi_0_dpsi/temp_e -! - call first_deriv(npoi,psi,temp_e,derpar) -! - A2e=derpar/temp_e - A1e=A1e-1.5d0*A2e -! -! -! thermodynamic forces for ions: -! - call first_deriv(npoi,psi,dens_i,derpar) -! - A1i=derpar/dens_i+e_i*dPhi_0_dpsi/temp_i -! - call first_deriv(npoi,psi,temp_i,derpar) -! - A2i=derpar/temp_i - A1i=A1i-1.5d0*A2i -! -! - do i=1,npoi -! - factor_of_Phi=imun*(mpol+qsaf(i)*ntor)/(qsaf(i)*dPhi_0_dpsi(i)) - omega_E=ntor*c*dPhi_0_dpsi(i) -! -! functions F_m and Gtor_m: -! -! electron contribution: - v_T=sqrt(temp_e(i)/e_mass) - x_1=(mpol+qsaf(i)*ntor)*v_T/(qsaf(i)*Rtor*anu_e(i)) - x_2=-omega_E/anu_e(i) - rho2_factor=(e_mass*c*v_T/(e_e*bcovar_phi(i)))**2*switch_flr_e -! -! susceptibility functions for electrons: -! - call getIfunc(x_1,x_2,symbI) -! - F_me = e_e*dens_e(i)*v_T**2/anu_e(i) & - * (symbI(1,1)*(A1e(i)+A2e(i))+0.5d0*symbI(1,3)*A2e(i)) - Gtor_me = 0.5d0*rho2_factor*e_e*dens_e(i)*v_T**2/anu_e(i) & - * (symbI(1,1)*(A1e(i)+2.d0*A2e(i))+0.5d0*symbI(1,3)*A2e(i)) - Htore = 0.5d0*rho2_factor*e_e*dens_e(i)*Rtor**2/dPhi_0_dpsi(i) & + + private + + ! types and associated procedures + public :: flr2_t, flr2_init, flr2_deinit, flr2_write, flr2_read, & + flr2_coeff, flr2_response_current + + type :: flr2_t + private + + integer :: npoi = 0, m_min = 0, m_max = 0 + real(dp), dimension(:), allocatable, public :: psi + complex(dp), dimension(:, :), allocatable, public :: a2_in, a2_out, a0 + complex(dp), dimension(:, :), allocatable, public :: b2_in, b2_out, b0 + complex(dp), dimension(:, :), allocatable, public :: c2_in, c2_out, c0 + complex(dp), dimension(:, :), allocatable, public :: d2_in, d2_out, d0 + end type flr2_t + + interface + subroutine hypergeometric1f1_cont_fract_1_modified_0_ada(b_re, b_im, z_re, z_im, f_re, f_im) & + bind(C, name = 'hypergeometric1f1_cont_fract_1_modified_0_ada') + use iso_c_binding, only: c_double + real(c_double), intent(in) :: b_re, b_im, z_re, z_im + real(c_double), intent(out) :: f_re, f_im + end subroutine hypergeometric1f1_cont_fract_1_modified_0_ada + + subroutine hypergeometric1f1_kummer_modified_0_ada(b_re, b_im, z_re, z_im, f_re, f_im) & + bind(C, name = 'hypergeometric1f1_kummer_modified_0_ada') + use iso_c_binding, only: c_double + real(c_double), intent(in) :: b_re, b_im, z_re, z_im + real(c_double), intent(out) :: f_re, f_im + end subroutine hypergeometric1f1_kummer_modified_0_ada + end interface + + integer, parameter :: mnmax = 3 + +contains + + subroutine flr2_init(flr2, npoi, m_min, m_max) + type(flr2_t), intent(inout) :: flr2 + integer, intent(in) :: npoi, m_min, m_max + + call flr2_deinit(flr2) + flr2%npoi = npoi + flr2%m_min = m_min + flr2%m_max = m_max + allocate(flr2%psi(npoi)) + allocate(flr2%a2_in(npoi, m_min:m_max)) + allocate(flr2%a2_out(npoi, m_min:m_max)) + allocate(flr2%a0(npoi, m_min:m_max)) + allocate(flr2%b2_in(npoi, m_min:m_max)) + allocate(flr2%b2_out(npoi, m_min:m_max)) + allocate(flr2%b0(npoi, m_min:m_max)) + allocate(flr2%c2_in(npoi, m_min:m_max)) + allocate(flr2%c2_out(npoi, m_min:m_max)) + allocate(flr2%c0(npoi, m_min:m_max)) + allocate(flr2%d2_in(npoi, m_min:m_max)) + allocate(flr2%d2_out(npoi, m_min:m_max)) + allocate(flr2%d0(npoi, m_min:m_max)) + end subroutine flr2_init + + subroutine flr2_deinit(flr2) + type(flr2_t), intent(inout) :: flr2 + + flr2%npoi = 0 + flr2%m_min = 0 + flr2%m_max = 0 + if (allocated(flr2%psi)) deallocate(flr2%psi) + if (allocated(flr2%a2_in)) deallocate(flr2%a2_in) + if (allocated(flr2%a2_out)) deallocate(flr2%a2_out) + if (allocated(flr2%a0)) deallocate(flr2%a0) + if (allocated(flr2%b2_in)) deallocate(flr2%b2_in) + if (allocated(flr2%b2_out)) deallocate(flr2%b2_out) + if (allocated(flr2%b0)) deallocate(flr2%b0) + if (allocated(flr2%c2_in)) deallocate(flr2%c2_in) + if (allocated(flr2%c2_out)) deallocate(flr2%c2_out) + if (allocated(flr2%c0)) deallocate(flr2%c0) + if (allocated(flr2%d2_in)) deallocate(flr2%d2_in) + if (allocated(flr2%d2_out)) deallocate(flr2%d2_out) + if (allocated(flr2%d0)) deallocate(flr2%d0) + end subroutine flr2_deinit + + subroutine flr2_write(flr2, file, dataset) + use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close + type(flr2_t), intent(in) :: flr2 + character(len = *), intent(in) :: file + character(len = *), intent(in) :: dataset + integer(HID_T) :: h5id_root + + call h5_open_rw(file, h5id_root) + call h5_create_parent_groups(h5id_root, trim(adjustl(dataset)) // '/') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/npoi', flr2%npoi, & + comment = 'number of flux surfaces') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/m_min', flr2%m_min, & + comment = 'minimum poloidal mode number including sign') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/m_max', flr2%m_max, & + comment = 'maximum poloidal mode number including sign') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', & + flr2%psi, lbound(flr2%psi), ubound(flr2%psi), unit = 'Mx', & + comment = 'poloidal flux') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', & + flr2%psi, lbound(flr2%psi), ubound(flr2%psi), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/a2_in', & + flr2%a2_in, lbound(flr2%a2_in), ubound(flr2%a2_in), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/a2_out', & + flr2%a2_out, lbound(flr2%a2_out), ubound(flr2%a2_out), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/a0', & + flr2%a0, lbound(flr2%a0), ubound(flr2%a0), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/b2_in', & + flr2%b2_in, lbound(flr2%b2_in), ubound(flr2%b2_in), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/b2_out', & + flr2%b2_out, lbound(flr2%b2_out), ubound(flr2%b2_out), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/b0', & + flr2%b0, lbound(flr2%b0), ubound(flr2%b0), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/c2_in', & + flr2%c2_in, lbound(flr2%c2_in), ubound(flr2%c2_in), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/c2_out', & + flr2%c2_out, lbound(flr2%c2_out), ubound(flr2%c2_out), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/c0', & + flr2%c0, lbound(flr2%c0), ubound(flr2%c0), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/d2_in', & + flr2%d2_in, lbound(flr2%d2_in), ubound(flr2%d2_in), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/d2_out', & + flr2%d2_out, lbound(flr2%d2_out), ubound(flr2%d2_out), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_add(h5id_root, trim(adjustl(dataset)) // '/d0', & + flr2%d0, lbound(flr2%d0), ubound(flr2%d0), unit = '?', & + comment = 'Poisson equation coefficient') + call h5_close(h5id_root) + end subroutine flr2_write + + subroutine flr2_read(flr2, file, dataset) + use hdf5_tools, only: HID_T, h5_open, h5_get, h5_close + type(flr2_t), intent(inout) :: flr2 + character(len = *), intent(in) :: file + character(len = *), intent(in) :: dataset + integer(HID_T) :: h5id_root + integer :: npoi, m_min, m_max + + call h5_open(file, h5id_root) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/npoi', npoi) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/m_min', m_min) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/m_max', m_max) + call flr2_init(flr2, npoi, m_min, m_max) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/psi', flr2%psi) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/a2_in', flr2%a2_in) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/a2_out', flr2%a2_out) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/a0', flr2%a0) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/b2_in', flr2%b2_in) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/b2_out', flr2%b2_out) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/b0', flr2%b0) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/c2_in', flr2%c2_in) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/c2_out', flr2%c2_out) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/c0', flr2%c0) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/d2_in', flr2%d2_in) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/d2_out', flr2%d2_out) + call h5_get(h5id_root, trim(adjustl(dataset)) // '/d0', flr2%d0) + call h5_close(h5id_root) + end subroutine flr2_read + + + ! All units are Gaussian if not otherwise mentioned + + ! Input/Output :: + ! flr2 - (flr2_t) Poisson equation coefficients + ! + ! Input :: + ! mpol_min - (integer) minimum poloidal mode number + ! mpol_max - (integer) maximum poloidal mode number + ! ntor - (integer) toroidal mode number + ! npoi - (integer) number of radial grid points + ! am_i - (real) ion mass number (mass in units of proton mass) + ! Z_i - (real) ion charge number (charge in units of elementary charge) + ! Rtor - (real) reference major radius + ! psi(npoi) - (real) unperturbed poloidal flux $\psi$ + ! qsaf(npoi) - (real) safety factor + ! bcovar_phi(npoi) - (real) co-variant unperturbed toroidal magnetic field component + ! Phi_0(npoi) - (real) equilibrium electrostatic potential + ! avR2nabpsi2(npoi) - (real) average over the polodal angle of $R^2 |\nabla \psi|^2$ + ! dens_e(npoi) - (real) electron density + ! temp_e(npoi) - (real) electron temperature + ! temp_i(npoi) - (real) ion temperature + ! anu_e(npoi) - (real) electron collision frequency + ! anu_i(npoi) - (real) ion collision frequency + subroutine flr2_coeff(flr2, mpol_min, mpol_max, ntor, npoi, am_i, z_i, Rtor, & + psi, qsaf, bcovar_phi, Phi_0, avR2nabpsi2, & + dens_e, temp_e, temp_i, anu_e, anu_i) + + use mephit_conf, only: logger + use mephit_util, only: imun, pi, c => clight, e_charge => elem_charge, e_mass => m_e, p_mass => m_p + + type(flr2_t), intent(inout) :: flr2 + integer, intent(in) :: mpol_min, mpol_max, ntor, npoi + integer :: i, mpol + real(dp), intent(in) :: am_i, Z_i, Rtor + real(dp) :: e_e,e_i,omega_E,x_1,x_2,v_T,rho2_factor + real(dp) :: switch_flr_e,switch_flr_i + real(dp) :: switch_cur_e,switch_cur_i + complex(dp) :: factor_of_Phi + complex(dp) :: F_m,Gtor_m,Htor + complex(dp) :: F_me,Gtor_me,Htore + complex(dp) :: F_mi,Gtor_mi,Htori + real(dp), dimension(:), intent(in) :: psi, qsaf, bcovar_phi, Phi_0, avR2nabpsi2, & + dens_e, temp_e, temp_i, anu_e, anu_i + + complex(dp), dimension(0:mnmax,0:mnmax) :: symbI + real(dp), dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i + + character(len=11), dimension(10), parameter :: array_names = [character(len=11) :: & + 'psi', 'qsaf', 'bcovar_phi', 'Phi_0', 'avR2nabpsi2', & + 'dens_e', 'temp_e', 'temp_i', 'anu_e', 'anu_i'] + integer, dimension(10) :: array_sizes + array_sizes = [ & + size(psi), size(qsaf), size(bcovar_phi), size(Phi_0), size(avR2nabpsi2), & + size(dens_e), size(temp_e), size(temp_i), size(anu_e), size(anu_i)] + do i = 1, 10 + if (npoi /= array_sizes(i)) then + call logger%msg_arg_size('flr2_coeff', & + 'npoi', 'size(' // trim(array_names(i)) // ')', npoi, array_sizes(i)) + if (logger%err) call logger%write_msg + error stop + end if + end do + if (npoi /= flr2%npoi) then + call logger%msg_arg_size('flr2_coeff', 'npoi', 'flr2%npoi', npoi, flr2%npoi) + if (logger%err) call logger%write_msg + error stop + end if + + flr2%psi(:) = psi + allocate(dens_i(npoi),dPhi_0_dpsi(npoi),derpar(npoi)) + allocate(A1e(npoi),A2e(npoi),A1i(npoi),A2i(npoi)) + + ! For playing around: + ! switch for FLR effects in electrons (1 - on, 0 - off): + switch_flr_e=1.d0 + ! switch for FLR effects in ions (1 - on, 0 - off): + switch_flr_i=1.d0 + ! switch for the electron contribution to the parallel current (1 - on, 0 - off): + switch_cur_e=1.d0 + ! switch for the ion contribution to the parallel current (1 - on, 0 - off): + switch_cur_i=1.d0 + + ! electron and ion charges and ion density: + e_e=-e_charge + e_i=e_charge*z_i + dens_i=dens_e/z_i + + call first_deriv(npoi,psi,Phi_0,dPhi_0_dpsi) + + + ! thermodynamic forces for electrons: + + call first_deriv(npoi,psi,dens_e,derpar) + + A1e=derpar/dens_e+e_e*dPhi_0_dpsi/temp_e + + call first_deriv(npoi,psi,temp_e,derpar) + + A2e=derpar/temp_e + A1e=A1e-1.5d0*A2e + + + ! thermodynamic forces for ions: + + call first_deriv(npoi,psi,dens_i,derpar) + + A1i=derpar/dens_i+e_i*dPhi_0_dpsi/temp_i + + call first_deriv(npoi,psi,temp_i,derpar) + + A2i=derpar/temp_i + A1i=A1i-1.5d0*A2i + + + do mpol = mpol_min, mpol_max + do i = 1, npoi + + factor_of_Phi=imun*(mpol+qsaf(i)*ntor)/(qsaf(i)*dPhi_0_dpsi(i)) + omega_E=ntor*c*dPhi_0_dpsi(i) + + ! functions F_m and Gtor_m: + + ! electron contribution: + v_T=sqrt(temp_e(i)/e_mass) + x_1=(mpol+qsaf(i)*ntor)*v_T/(qsaf(i)*Rtor*anu_e(i)) + x_2=-omega_E/anu_e(i) + rho2_factor=(e_mass*c*v_T/(e_e*bcovar_phi(i)))**2*switch_flr_e + + ! susceptibility functions for electrons: + + call getIfunc(x_1,x_2,symbI) + + F_me = e_e*dens_e(i)*v_T**2/anu_e(i) & + * (symbI(1,1)*(A1e(i)+A2e(i))+0.5d0*symbI(1,3)*A2e(i)) + Gtor_me = 0.5d0*rho2_factor*e_e*dens_e(i)*v_T**2/anu_e(i) & + * (symbI(1,1)*(A1e(i)+2.d0*A2e(i))+0.5d0*symbI(1,3)*A2e(i)) + Htore = 0.5d0*rho2_factor*e_e*dens_e(i)*Rtor**2/dPhi_0_dpsi(i) & * (A1e(i)+2.5d0*A2e(i)) -! -! ion contribution: - v_T=sqrt(temp_i(i)/(am_i*p_mass)) - x_1=(mpol+qsaf(i)*ntor)*v_T/(qsaf(i)*Rtor*anu_i(i)) - x_2=-omega_E/anu_i(i) - rho2_factor=(am_i*p_mass*c*v_T/(e_i*bcovar_phi(i)))**2*switch_flr_i -! -! susceptibility functions for ions: -! - call getIfunc(x_1,x_2,symbI) -! - F_mi = e_i*dens_i(i)*v_T**2/anu_i(i) & - * (symbI(1,1)*(A1i(i)+A2i(i))+0.5d0*symbI(1,3)*A2i(i)) - Gtor_mi = 0.5d0*rho2_factor*e_i*dens_i(i)*v_T**2/anu_i(i) & - * (symbI(1,1)*(A1i(i)+2.d0*A2i(i))+0.5d0*symbI(1,3)*A2i(i)) - Htori = 0.5d0*rho2_factor*e_i*dens_i(i)*Rtor**2/dPhi_0_dpsi(i) & + + ! ion contribution: + v_T=sqrt(temp_i(i)/(am_i*p_mass)) + x_1=(mpol+qsaf(i)*ntor)*v_T/(qsaf(i)*Rtor*anu_i(i)) + x_2=-omega_E/anu_i(i) + rho2_factor=(am_i*p_mass*c*v_T/(e_i*bcovar_phi(i)))**2*switch_flr_i + + ! susceptibility functions for ions: + + call getIfunc(x_1,x_2,symbI) + + F_mi = e_i*dens_i(i)*v_T**2/anu_i(i) & + * (symbI(1,1)*(A1i(i)+A2i(i))+0.5d0*symbI(1,3)*A2i(i)) + Gtor_mi = 0.5d0*rho2_factor*e_i*dens_i(i)*v_T**2/anu_i(i) & + * (symbI(1,1)*(A1i(i)+2.d0*A2i(i))+0.5d0*symbI(1,3)*A2i(i)) + Htori = 0.5d0*rho2_factor*e_i*dens_i(i)*Rtor**2/dPhi_0_dpsi(i) & * (A1i(i)+2.5d0*A2i(i)) -! - F_m = F_me + F_mi - Gtor_m = Gtor_me + Gtor_mi - Htor = Htore + Htori -! - b0(i) = -4.d0*pi*(mpol+qsaf(i)*ntor)*F_m/(qsaf(i)*omega_E*avR2nabpsi2(i)) - b2_in(i) = -4.d0*pi*(mpol+qsaf(i)*ntor)*Gtor_m/(qsaf(i)*omega_E) - b2_out(i) = b2_in(i) -! - a0(i) = -b0(i)*factor_of_Phi - a2_in(i) = -b2_in(i)*factor_of_Phi+4.d0*pi*Htor - a2_out(i) = (1.d0,0.d0)+a2_in(i) -! - d0(i) = -(F_me*switch_cur_e+F_mi*switch_cur_i)/bcovar_phi(i) - d2_in(i) = -(Gtor_me*switch_cur_e+Gtor_mi*switch_cur_i)*avR2nabpsi2(i)/bcovar_phi(i) - d2_out(i) = d2_in(i) -! - c0(i) = d0(i)*factor_of_Phi - c2_in(i) = d2_in(i)*factor_of_Phi - c2_out(i) = c2_in(i) - enddo -! -! Solve Laplace equation: -! - call progonka(isw_Phi_m,npoi, & - psi,bpsi_over_bphi, & - a2_in,a2_out,a0, & - b2_in,b2_out,b0, & - c2_in,c2_out,c0, & - d2_in,d2_out,d0, & - Phi_m,parcur_over_b0) -! - deallocate(dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i) - deallocate(a2_in,a2_out,a0) - deallocate(b2_in,b2_out,b0) - deallocate(c2_in,c2_out,c0) - deallocate(d2_in,d2_out,d0) -! - end subroutine response_current - -!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -subroutine getIfunc(x1,x2,symbI) - integer, parameter :: mnmax=3 - integer :: m,n - double precision :: x1,x2,z - double complex :: denom - double complex, dimension(0:mnmax,0:mnmax) :: symbI,Imn -! -! if(.true.) then - if(.false.) then -! collisionless case: - symbI=(0.d0,0.d0) - z=x2/(sqrt(2.d0)*x1) -! - symbI(0,0)=sqrt(2.d0)*exp(-z**2)/abs(x1) - symbI(1,0)=symbI(0,0)*x2/x1 - symbI(1,1)=symbI(1,0)*x2/x1 -! - symbI(2,0)=symbI(0,0)*(x2/x1)**2 - symbI(2,1)=symbI(1,0)*(x2/x1)**2 - symbI(3,0)=symbI(2,1) - symbI(3,1)=symbI(1,1)*(x2/x1)**2 -! - symbI(2,2)=symbI(2,0)*(x2/x1)**2 - symbI(3,2)=symbI(2,1)*(x2/x1)**2 - symbI(3,3)=symbI(3,1)*(x2/x1)**2 - else -! collisional case: -! - call W2_arr(x1,x2,Imn) -! - if(.true.) then -! if(.false.) then -! energy conservation: - denom=(1.d0,0.d0)-Imn(0,0)+(2.d0,0.d0)*Imn(2,0)-Imn(2,2) - do m=0,3 - do n=0,3 - symbI(m,n)=Imn(m,n)+(Imn(m,0)-Imn(m,2))*(Imn(n,0)-Imn(n,2))/denom - enddo - enddo + + F_m = F_me + F_mi + Gtor_m = Gtor_me + Gtor_mi + Htor = Htore + Htori + + flr2%b0(i, mpol) = -4.d0 * pi * (mpol + qsaf(i) * ntor) * F_m / (qsaf(i) * omega_E * avR2nabpsi2(i)) + flr2%b2_in(i, mpol) = -4.d0 * pi * (mpol + qsaf(i) * ntor) * Gtor_m / (qsaf(i) * omega_E) + flr2%b2_out(i, mpol) = flr2%b2_in(i, mpol) + + flr2%a0(i, mpol) = -flr2%b0(i, mpol) * factor_of_Phi + flr2%a2_in(i, mpol) = -flr2%b2_in(i, mpol) * factor_of_Phi + 4.d0 * pi * Htor + flr2%a2_out(i, mpol) = (1.d0, 0.d0) + flr2%a2_in(i, mpol) + + flr2%d0(i, mpol) = -(F_me * switch_cur_e + F_mi * switch_cur_i) / bcovar_phi(i) + flr2%d2_in(i, mpol) = -(Gtor_me * switch_cur_e + Gtor_mi * switch_cur_i) * avR2nabpsi2(i) / bcovar_phi(i) + flr2%d2_out(i, mpol) = flr2%d2_in(i, mpol) + + flr2%c0(i, mpol) = flr2%d0(i, mpol) * factor_of_Phi + flr2%c2_in(i, mpol) = flr2%d2_in(i, mpol) * factor_of_Phi + flr2%c2_out(i, mpol) = flr2%c2_in(i, mpol) + end do + end do + + deallocate(dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i) + end subroutine flr2_coeff + + + ! Input/Output :: + ! flr2 - (flr2_t) Poisson equation coefficients + ! + ! Input :: + ! isw_Phi_m - (integer) switch for the potential perturbation in the parallel current: + ! 0 - set it to zero in the current, 1 - use it as is in the current + ! mpol - (integer) poloidal mode number + ! npoi - (integer) number of radial grid points + ! psi(npoi) - (real) unperturbed poloidal flux $\psi$ + ! bpsi_over_bphi(npoi) - (complex) psi-component of the perturbation magnetic field divided by + ! contra-variant component of the unperturbed magnetic field + ! Output :: + ! parcur_over_b0(npoi) - (complex) parallel current density divided by uperturbed magnetic field module + ! Phi_m(npoi) - (complex) perturbation of the electrostatic potential + subroutine flr2_response_current(flr2, isw_Phi_m, mpol, npoi, bpsi_over_bphi, parcur_over_b0, Phi_m) + + use mephit_conf, only: logger + + type(flr2_t), intent(in) :: flr2 + integer, intent(in) :: isw_Phi_m, mpol, npoi + complex(dp), dimension(:), intent(in) :: bpsi_over_bphi + complex(dp), dimension(:), intent(out) :: parcur_over_b0, Phi_m + + integer :: i + character(len=14), dimension(3), parameter :: array_names = & + [character(len=14) :: 'bpsi_over_bphi', 'parcur_over_b0', 'Phi_m'] + integer, dimension(3) :: array_sizes + array_sizes = [size(bpsi_over_bphi), size(parcur_over_b0), size(Phi_m)] + do i = 1, 3 + if (npoi /= array_sizes(i)) then + call logger%msg_arg_size('flr2_response_current', & + 'npoi', 'size(' // trim(array_names(i)) // ')', npoi, array_sizes(i)) + if (logger%err) call logger%write_msg + error stop + end if + end do + if (npoi /= flr2%npoi) then + call logger%msg_arg_size('flr2_response_current', 'npoi', 'flr2%npoi', npoi, flr2%npoi) + if (logger%err) call logger%write_msg + error stop + end if + + ! Solve Poisson equation: + + call progonka(isw_Phi_m, npoi, & + flr2%psi, bpsi_over_bphi, & + flr2%a2_in(:, mpol), flr2%a2_out(:, mpol), flr2%a0(:, mpol), & + flr2%b2_in(:, mpol), flr2%b2_out(:, mpol), flr2%b0(:, mpol), & + flr2%c2_in(:, mpol), flr2%c2_out(:, mpol), flr2%c0(:, mpol), & + flr2%d2_in(:, mpol), flr2%d2_out(:, mpol), flr2%d0(:, mpol), & + Phi_m, parcur_over_b0) + end subroutine flr2_response_current + + + subroutine getIfunc(x1, x2, symbI) + integer :: m, n + real(dp), intent(in) :: x1, x2 + real(dp) :: z + complex(dp) :: denom + complex(dp), dimension(0:mnmax, 0:mnmax), intent(out) :: symbI + complex(dp), dimension(0:mnmax, 0:mnmax) :: Imn + + ! if(.true.) then + if(.false.) then + ! collisionless case: + symbI=(0.d0,0.d0) + z=x2/(sqrt(2.d0)*x1) + + symbI(0,0)=sqrt(2.d0)*exp(-z**2)/abs(x1) + symbI(1,0)=symbI(0,0)*x2/x1 + symbI(1,1)=symbI(1,0)*x2/x1 + + symbI(2,0)=symbI(0,0)*(x2/x1)**2 + symbI(2,1)=symbI(1,0)*(x2/x1)**2 + symbI(3,0)=symbI(2,1) + symbI(3,1)=symbI(1,1)*(x2/x1)**2 + + symbI(2,2)=symbI(2,0)*(x2/x1)**2 + symbI(3,2)=symbI(2,1)*(x2/x1)**2 + symbI(3,3)=symbI(3,1)*(x2/x1)**2 else - symbI=Imn - endif -!write (1234,*) x1,real(Imn(0,0)),dimag(Imn(0,0)),real(Imn(2,0)),dimag(Imn(2,0)), & -!real(Imn(2,2)),dimag(Imn(2,2)),real(Imn(1,0)-Imn(1,2)),dimag(Imn(1,0)-Imn(1,2)) -! - endif -! -end subroutine getIfunc -!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -! Calculates array of W2 special functions. + subroutine W2_arr(x1_in, x2_in, Imn) + real(dp), intent(in) :: x1_in, x2_in + complex(dp), dimension(0:3, 0:3), intent(out) :: Imn + complex(dp) :: t1, t2, F11m, x1, x2 + complex(dp), parameter :: I = (0.0d0, 1.0d0), one = (1.0d0, 0.0d0) + + integer :: l,m,n + + real(dp) :: F_im, F_re + complex(dp), dimension(0:3, 0:3, 0:mnmax) :: W2 + + W2 = (0.0d0, 0.0d0) + + x1 = cmplx(x1_in, 0.d0, dp) + t1 = x1**2 + + do l = 0, mnmax ! does not work for 3, 2 and 1! + + + x2 = x2_in + I*l + + t2 = - I*x2 + t1 + + ! call Hypergeometric1F1_kummer_modified_0_ada (real(one+t2), aimag(one+t2), real(t1), aimag(t1), F_re, F_im) + + ! F11m = F_re + I*F_im + + call hypergeometric1f1_cont_fract_1_modified_0_ada (real(one+t2), aimag(one+t2), real(t1), aimag(t1), F_re, F_im) + + F11m = F_re + I*F_im + + !asymptotic form of expressions: has no fake singularity at kp=0 + + W2(0,0,l) = & + (-I)*(-x2**2 + 2 + 5*x1**2 - & + (3*I)*x2*(1 + x1**2) + (3 + F11m)*x1**4)/ & + ((x2 + I*x1**2)*(x2 + I*(1 + x1**2))* & + (x2 + I*(2 + x1**2))) + + W2(0,1,l) = & + (-I)*x1*(x2*(-1 + F11m*x1**2) - & + I*(2 + 3*x1**2+ x1**4))/ & + ((x2 + I*x1**2)*(x2 + I*(1 + x1**2))* & + (x2 + I*(2 + x1**2))) + + W2(0,2,l) = & + -((x2 + I)*(2 + 3*x1**2 - & + I*x2*(1 - F11m*x1**2) + x1**4))/ & + ((x2 + I*x1**2)*(x2 + I*(1 + x1**2))* & + (x2 + I*(2 + x1**2))) + + W2(0,3,l) = & + (-I)*x1*(F11m*x2**3- (3 + 2*F11m)*x2 + & + I*x2**2*(3*F11m - x1**2) - & + I*(6 + (7 + 2*F11m)*x1**2 + x1**4))/ & + ((x2+ I*x1**2)*(x2+ I*(1 + x1**2))* & + (x2+ I*(2 + x1**2))) + + W2(1,0,l) = & + (-I)*x1*(x2*(-1 + F11m*x1**2) - & + I*(2 + 3*x1**2 + x1**4))/ & + ((x2+ I*x1**2)*(x2+ I*(1 + x1**2))* & + (x2+ I*(2 + x1**2))) + + W2(1,1,l) = & + (-I)*x2*(x2*(-1 + F11m*x1**2) - & + I*(2 + 3*x1**2 + x1**4))/ & + ((x2+ I*x1**2)*(x2+ I*(1 + x1**2))* & + (x2+ I*(2 + x1**2))) + + W2(1,2,l) = & + (-I)*x1*(F11m*x2**3+ I*x2**2*(F11m - x1**2) - & + x2*(3 + 2*x1**2) - & + I*(2 + 3*x1**2 + x1**4))/ & + ((x2+ I*x1**2)*(x2+ I*(1 + x1**2))* & + (x2+ I*(2 + x1**2))) + + W2(1,3,l) = & + x2*((-I)*F11m*x2**3+ I*(3 + 2*F11m)*x2 - & + 6 - (7 + 2*F11m)*x1**2 + & + x2**2*(3*F11m - x1**2) - x1**4)/ & + ((x2+ I*x1**2)*(x2+ I*(1 + x1**2))* & + (x2+ I*(2 + x1**2))) + + + end do + + W2(2,2,0) = I*x1*W2(1,2,1) + 2.d0*W2(1,1,1) - I*x1*W2(1,2,0) + W2(0,2,0) + W2(2,2,1) = I*x1*W2(1,2,2) + 2.d0*W2(1,1,2) - I*x1*W2(1,2,1) + W2(0,2,1) + W2(2,3,0) = I*x1*W2(2,2,1) + 2.d0*W2(1,2,1) - I*x1*W2(2,2,0) + 2.0d0*W2(1,2,0) + W2(2,3,1) = I*x1*W2(1,3,2) + 3.d0*W2(1,2,2) - I*x1*W2(1,3,1) + W2(0,3,1) + W2(3,3,0) = I*x1*W2(2,3,1) + 3.d0*W2(2,2,1) - I*x1*W2(2,3,0) + 2.0d0*W2(1,3,0) + + Imn = W2(:,:,0) + do m=0,3 + do n=0,m-1 + Imn(m,n)=Imn(n,m) + end do + end do + + end subroutine W2_arr + + + subroutine progonka(isw_f, npoi, x, q, & + a2_in, a2_out, a0, & + b2_in, b2_out, b0, & + c2_in, c2_out, c0, & + d2_in, d2_out, d0, & + f, g) + + ! Solves 2-nd order ODE for the unknown $f$ + + !$$ + ! a_2^{out}(x) \difp{^2}{x^2} f(x) + \difp{^2}{x^2} a_2^{in}(x) f(x) + a_0(x) f(x) = S(x) + !$$ + + ! where the source term is given by derivatives of a given function $q(x)$ + + !$$ + ! S(x) = b_2^{out}(x) \difp{^2}{x^2} q(x) + \difp{^2}{x^2} b_2^{in}(x) q(x) + b_0(x) q(x) + !$$ + + ! Boundary conditions at both interval ends are + + !$$ + ! a_0 f = b_0 q + !$$ + + ! Using the result, computes the following function + + !$$ + ! g(x) = c_2^{out}(x) \difp{^2}{x^2} f(x) + \difp{^2}{x^2} c_2^{in}(x) f(x) + c_0(x) f(x) + ! + d_2^{out}(x) \difp{^2}{x^2} q(x) + \difp{^2}{x^2} d_2^{in}(x) q(x) + d_0(x) q(x) + !$$ + + + ! Input: + ! isw_f - (integer) switch for account of f(x) in g(x): + ! 0 - skip the first line in the expression for g(x), 1 - compute g(x) as is + ! npoi - (integer) number of grid points + ! x(npoi) - (real) array of independent variable, need not to be equidistant + ! c_1(npoi) - (complex) array of coefficient $C_1$ + ! c_2(npoi) - (complex) array of coefficient $C_2$ + ! b(npoi) - (complex) array of source $B$ + + ! Output: + ! phi(npoi) - (complex) array of the solution $\Phi$ + + use mephit_conf, only: logger + + integer, intent(in) :: isw_f, npoi + integer :: i + + real(dp) :: dxp,dxm,dxt + complex(dp) :: denom + + real(dp), dimension(:), intent(in) :: x + complex(dp), dimension(:), intent(in) :: q, & + a2_in, a2_out, a0, & + b2_in, b2_out, b0, & + c2_in, c2_out, c0, & + d2_in, d2_out, d0 + complex(dp), dimension(:), intent(out) :: f, g + complex(dp), dimension(:), allocatable :: alp,bet,quelle + complex(dp), dimension(:,:), allocatable :: wsecder,eqmat + + character(len=6), dimension(16), parameter :: array_names = [character(len=6) :: & + 'x', 'q', & + 'a2_in', 'a2_out', 'a0', 'b2_in', 'b2_out', 'b0', & + 'c2_in', 'c2_out', 'c0', 'd2_in', 'd2_out', 'd0', & + 'f', 'g'] + integer, dimension(16) :: array_sizes + array_sizes = [ & + size(x), size(q), & + size(a2_in), size(a2_out), size(a0), size(b2_in), size(b2_out), size(b0), & + size(c2_in), size(c2_out), size(c0), size(d2_in), size(d2_out), size(d0), & + size(f), size(g)] + do i = 1, 16 + if (npoi /= array_sizes(i)) then + call logger%msg_arg_size('progonka', & + 'npoi', 'size(' // trim(array_names(i)) // ')', npoi, array_sizes(i)) + if (logger%err) call logger%write_msg + error stop + end if + end do + + allocate(wsecder(-1:1,npoi),eqmat(-1:1,npoi),alp(npoi),bet(npoi),quelle(npoi)) + + wsecder(:,1)=0.d0 + wsecder(:,npoi)=0.d0 + + ! weights wsecder are defined so that the second derivative of function F(x) at the point x(i) is + ! ddF(i)=sum(F(i-1:i+1)*wsecder(:,i)) + + do i=2,npoi-1 + dxp = x(i+1)-x(i) + dxm = x(i)-x(i-1) + dxt = dxp+dxm + wsecder(-1,i) = 2.d0/(dxm*dxt) + wsecder( 0,i) = -2.d0/(dxm*dxp) + wsecder( 1,i) = 2.d0/(dxp*dxt) + + quelle(i) = b2_out(i)*sum(q(i-1:i+1)*wsecder(:,i)) & + + sum(b2_in(i-1:i+1)*q(i-1:i+1)*wsecder(:,i)) & + + b0(i)*q(i) + + eqmat(:,i) = a2_out(i)*wsecder(:,i) & + + a2_in(i-1:i+1)*wsecder(:,i) + eqmat(0,i) = eqmat(0,i) + a0(i) + end do + + eqmat(0,1)=(1.d0,0.d0) + eqmat(0,npoi)=(1.d0,0.d0) + + alp(2)=(0.d0,0.d0) + bet(2)=b0(1)*q(1)/a0(1) + + do i=2,npoi-1 + denom=eqmat(0,i)+alp(i)*eqmat(-1,i) + alp(i+1)=-eqmat(1,i)/denom + bet(i+1)=(quelle(i)-eqmat(-1,i)*bet(i))/denom + end do + + f(npoi)=b0(npoi)*q(npoi)/a0(npoi) + + do i=npoi-1,1,-1 + f(i)=alp(i+1)*f(i+1)+bet(i+1) + end do + + f=f*dble(isw_f) + + do i=2,npoi-1 + + g(i) = c2_out(i)*sum(f(i-1:i+1)*wsecder(:,i)) & + + sum(c2_in(i-1:i+1)*f(i-1:i+1)*wsecder(:,i)) & + + c0(i)*f(i) & + + d2_out(i)*sum(q(i-1:i+1)*wsecder(:,i)) & + + sum(d2_in(i-1:i+1)*q(i-1:i+1)*wsecder(:,i)) & + + d0(i)*q(i) + + end do + + deallocate(wsecder,eqmat,alp,bet,quelle) + end subroutine progonka -! -!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -! - subroutine first_deriv(npoi,x,f,df_dx) -! -! Computes first derivative df_dx of function f over variable x -! -! Input: -! npoi - (integer) number of grid points -! x(npoi) - (double precision) array of independent variable, need not to be equidistant -! f(npoi) - (double precision) array of function values -! -! Output: -! df_dx(npoi) - (double precision) array of derivative values -! -! - implicit none -! - integer :: npoi,i -! - double precision :: dxp,dxm,dxt -! - double precision, dimension(npoi) :: x,f,df_dx -! - do i=2,npoi-1 - dxp=x(i+1)-x(i) - dxm=x(i)-x(i-1) - dxt=dxp+dxm - df_dx(i)=((f(i+1)-f(i))*dxm**2-(f(i-1)-f(i))*dxp**2)/(dxp*dxm*dxt) - enddo -! - df_dx(1)=(df_dx(2)*(x(3)-x(1))+df_dx(3)*(x(1)-x(2)))/(x(3)-x(2)) - df_dx(npoi)=(df_dx(npoi-1)*(x(npoi-2)-x(npoi))+df_dx(npoi-2)*(x(npoi)-x(npoi-1)))/(x(npoi-2)-x(npoi-1)) -! + + + subroutine first_deriv(npoi, x, f, df_dx) + + ! Computes first derivative df_dx of function f over variable x + + ! Input: + ! npoi - (integer) number of grid points + ! x(npoi) - (real) array of independent variable, need not to be equidistant + ! f(npoi) - (real) array of function values + + ! Output: + ! df_dx(npoi) - (real) array of derivative values + + use mephit_conf, only: logger + + integer, intent(in) :: npoi + integer :: i + + real(dp) :: dxp,dxm,dxt + + real(dp), dimension(:), intent(in) :: x, f + real(dp), dimension(:), intent(out) :: df_dx + + character(len=5), dimension(3), parameter :: array_names = & + [character(len=5) :: 'x', 'f', 'df_dx'] + integer, dimension(3) :: array_sizes + array_sizes = [size(x), size(f), size(df_dx)] + do i = 1, 3 + if (npoi /= array_sizes(i)) then + call logger%msg_arg_size('first_deriv', & + 'npoi', 'size(' // trim(array_names(i)) // ')', npoi, array_sizes(i)) + if (logger%err) call logger%write_msg + error stop + end if + end do + + do i=2,npoi-1 + dxp=x(i+1)-x(i) + dxm=x(i)-x(i-1) + dxt=dxp+dxm + df_dx(i)=((f(i+1)-f(i))*dxm**2-(f(i-1)-f(i))*dxp**2)/(dxp*dxm*dxt) + end do + + df_dx(1)=(df_dx(2)*(x(3)-x(1))+df_dx(3)*(x(1)-x(2)))/(x(3)-x(2)) + df_dx(npoi)=(df_dx(npoi-1)*(x(npoi-2)-x(npoi))+df_dx(npoi-2)*(x(npoi)-x(npoi-1)))/(x(npoi-2)-x(npoi-1)) + end subroutine first_deriv + +end module mephit_flr2 diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index 0e0657f..255b280 100644 --- a/src/mephit_iter.F90 +++ b/src/mephit_iter.F90 @@ -109,7 +109,9 @@ end function FEM_test subroutine mephit_run(runmode, config, suffix) bind(C, name = 'mephit_run') use iso_c_binding, only: c_int, c_ptr use input_files, only: gfile + use field_sub, only : read_field_input use geqdsk_tools, only: geqdsk_read, geqdsk_classify, geqdsk_standardise + use hdf5_tools, only: h5_init, h5overwrite use mephit_util, only: C_F_string, init_field, geqdsk_scale, geqdsk_export_hdf5, geqdsk_import_hdf5, & save_symfluxcoord, load_symfluxcoord use mephit_conf, only: conf, config_read, config_export_hdf5, conf_arr, logger, & @@ -117,8 +119,7 @@ subroutine mephit_run(runmode, config, suffix) bind(C, name = 'mephit_run') use mephit_mesh, only: equil, mesh, generate_mesh, mesh_write, mesh_read, write_cache, read_cache, & resample_profiles, write_profiles_hdf5, read_profiles_hdf5 use mephit_pert, only: generate_vacfield, vac, vac_init, vac_write, vac_read - use hdf5_tools, only: h5_init, h5overwrite - use field_sub, only : read_field_input + use mephit_flr2, only: flr2_t, flr2_write, flr2_read, flr2_deinit integer(c_int), intent(in), value :: runmode type(c_ptr), intent(in), value :: config type(c_ptr), intent(in), value :: suffix @@ -126,6 +127,7 @@ subroutine mephit_run(runmode, config, suffix) bind(C, name = 'mephit_run') integer(c_int) :: runmode_flags logical :: meshing, preconditioner, iterations type(fdm_t) :: fdm + type(flr2_t) :: flr2 type(perteq_t) :: perteq type(precond_t) :: precond @@ -197,17 +199,21 @@ subroutine mephit_run(runmode, config, suffix) bind(C, name = 'mephit_run') if (preconditioner) then call FDM_compute_matrix(fdm) call FDM_write(fdm, datafile, 'iter/FDM') - call debug_initial_iteration(perteq, fdm) - call precond_compute(precond, perteq, fdm) + call compute_FLR2_coeff(flr2) + call FLR2_write(flr2, datafile, 'iter/FLR2') + call debug_initial_iteration(perteq, fdm, flr2) + call precond_compute(precond, perteq, fdm, flr2) call precond_write(precond, datafile, 'iter') else call FDM_read(fdm, datafile, 'iter/FDM') + call FLR2_read(flr2, datafile, 'iter/FLR2') call precond_read(precond, datafile, 'iter') end if if (iterations) then - call perteq_iterate(perteq, precond, fdm) + call perteq_iterate(perteq, precond, fdm, flr2) end if call FDM_deinit(fdm) + call FLR2_deinit(flr2) call precond_deinit(precond) call perteq_deinit(perteq) end if @@ -377,15 +383,17 @@ pure function precond_apply(precond, vec) matmul(transpose(conjg(precond%eigvecs(:, 1:precond%nritz))), vec))) end function precond_apply - subroutine precond_compute(precond, perteq, fdm) + subroutine precond_compute(precond, perteq, fdm, flr2) use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close use mephit_conf, only: conf, logger, datafile, cmplx_fmt, runmode_precon use mephit_util, only: arnoldi_break use mephit_mesh, only: mesh use mephit_pert, only: vac, RT0_init, RT0_deinit, RT0_write, RT0_tor_comp_from_zero_div, RT0_L2int + use mephit_flr2, only: flr2_t type(precond_t), intent(inout) :: precond type(perteq_t), intent(inout) :: perteq type(fdm_t), intent(in) :: fdm + type(flr2_t), intent(in) :: flr2 integer :: i, j, info integer(HID_T) :: h5id_root complex(dp), allocatable :: start_vector(:), eigvals(:), Yr(:, :) @@ -398,7 +406,7 @@ subroutine precond_compute(precond, perteq, fdm) perteq%Bn%DOF(:) = vac%Bn%DOF perteq%Bn%comp_phi(:) = vac%Bn%comp_phi call compute_presn(perteq, fdm, conf%damp) - call compute_currn(perteq, fdm, conf%damp, .false.) + call compute_currn(perteq, fdm, flr2, conf%damp, .false.) call compute_magfn(perteq) allocate(start_vector, source = perteq%Bn%DOF) ! calculate eigenvectors -- system dimension: number of non-redundant edges in core plasma @@ -467,7 +475,7 @@ subroutine precond_compute(precond, perteq, fdm) call RT0_init(Bn_prev, mesh%nedge, mesh%ntri) Bn_prev%DOF(:) = perteq%Bn%DOF call compute_presn(perteq, fdm, conf%damp) - call compute_currn(perteq, fdm, conf%damp, .false.) + call compute_currn(perteq, fdm, flr2, conf%damp, .false.) call compute_magfn(perteq) perteq%Bn%DOF(:) = perteq%Bn%DOF - precond_apply(precond, perteq%Bn%DOF - Bn_prev%DOF) call RT0_tor_comp_from_zero_div(perteq%Bn) @@ -487,21 +495,23 @@ subroutine next_iteration_arnoldi(old_val, new_val) perteq%Bn%DOF(:) = old_val call RT0_tor_comp_from_zero_div(perteq%Bn) call compute_presn(perteq, fdm, conf%damp) - call compute_currn(perteq, fdm, conf%damp, .false.) + call compute_currn(perteq, fdm, flr2, conf%damp, .false.) call compute_magfn(perteq) new_val(:) = perteq%Bn%DOF end subroutine next_iteration_arnoldi end subroutine precond_compute - subroutine perteq_iterate(perteq, precond, fdm) + subroutine perteq_iterate(perteq, precond, fdm, flr2) use ieee_arithmetic, only: ieee_value, ieee_quiet_nan use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close use mephit_conf, only: conf, logger, datafile, runmode_single, currn_model_kilca use mephit_mesh, only: mesh use mephit_pert, only: vac, L1_write, RT0_init, RT0_deinit, RT0_tor_comp_from_zero_div, RT0_L2int + use mephit_flr2, only: flr2_t type(perteq_t), intent(inout) :: perteq type(precond_t), intent(inout) :: precond type(fdm_t), intent(in) :: fdm + type(flr2_t), intent(in) :: flr2 integer :: kiter, niter, maxiter integer(HID_T) :: h5id_root real(dp), allocatable :: L2int_Bn_diff(:) @@ -554,7 +564,7 @@ subroutine perteq_iterate(perteq, precond, fdm) call perteq_write('("iter/", a, "' // postfix // '")', & ' (after iteration)', presn = perteq%pn, presmn = perteq%pn) end if - call compute_currn(perteq, fdm, conf%damp, .false.) + call compute_currn(perteq, fdm, flr2, conf%damp, .false.) call compute_magfn(perteq) perteq%Bn%DOF(:) = perteq%Bn%DOF + vac%Bn%DOF if (precond%nritz > 0) then @@ -627,11 +637,13 @@ subroutine perteq_iterate(perteq, precond, fdm) call RT0_deinit(Bn_diff) end subroutine perteq_iterate - subroutine debug_initial_iteration(perteq, fdm) + subroutine debug_initial_iteration(perteq, fdm, flr2) use mephit_conf, only: conf use mephit_pert, only: vac + use mephit_flr2, only: flr2_t type(perteq_t), intent(inout) :: perteq type(fdm_t), intent(in) :: fdm + type(flr2_t), intent(in) :: flr2 if (conf%debug_initial) then perteq%Bn%DOF(:) = vac%Bn%DOF @@ -646,11 +658,11 @@ subroutine debug_initial_iteration(perteq, fdm) call perteq_write('("debug_MFEM_initial/", a)', & ' (initial iteration)', presn = perteq%pn, presmn = perteq%pn) #endif - call compute_currn(perteq, fdm, .false., .true.) + call compute_currn(perteq, fdm, flr2, .false., .true.) perteq%Bn%DOF(:) = vac%Bn%DOF perteq%Bn%comp_phi(:) = vac%Bn%comp_phi call compute_presn(perteq, fdm, .true.) - call compute_currn(perteq, fdm, .true., .true.) + call compute_currn(perteq, fdm, flr2, .true., .true.) end if end subroutine debug_initial_iteration @@ -1046,14 +1058,16 @@ subroutine debug_MDE(group, presn, magfn, currn_perp, currn_par) deallocate(lorentz, pn, grad_pn, Bn_psi_contravar, jnpar_Bmod, grad_jnpar_Bmod, div_jnperp, div_jnperp_RT0) end subroutine debug_MDE - subroutine compute_currn(perteq, fdm, apply_damping, debug_initial) + subroutine compute_currn(perteq, fdm, flr2, apply_damping, debug_initial) use mephit_conf, only: conf, currn_model_kilca, currn_model_mhd, logger, datafile use mephit_mesh, only: mesh use mephit_pert, only: polmodes_t, polmodes_init, polmodes_write, polmodes_deinit, & L1_t, L1_init, L1_interp, L1_deinit, RT0_t, RT0_init, RT0_deinit, & RT0_interp, RT0_project_pol_comp, RT0_project_tor_comp + use mephit_flr2, only: flr2_t type(perteq_t), intent(inout) :: perteq type(fdm_t), intent(in) :: fdm + type(flr2_t), intent(in) :: flr2 logical, intent(in) :: apply_damping logical, intent(in) :: debug_initial type(polmodes_t) :: resonant_jmnpar_over_Bmod @@ -1079,8 +1093,8 @@ subroutine compute_currn(perteq, fdm, apply_damping, debug_initial) case (currn_model_mhd) call compute_shielding_current(perteq%pn, resonant_jmnpar_over_Bmod) case (currn_model_kilca) - call compute_kilca_current(perteq%Bn, resonant_jmnpar_over_Bmod, & - perteq%Phi_mn, perteq%Phi_aligned_mn) + call compute_flr2_current(perteq%Bn, resonant_jmnpar_over_Bmod, & + flr2, perteq%Phi_mn, perteq%Phi_aligned_mn) case default write (logger%msg, '("unknown response current model selection", i0)') conf%currn_model if (logger%err) call logger%write_msg @@ -1241,17 +1255,34 @@ subroutine helical_current_from_parallel_current(jmnpar_over_Bmod, fdm, jnpar_ov call L1_deinit(coeff_jnperp) end subroutine helical_current_from_parallel_current - subroutine compute_kilca_current(Bn, jmnpar_over_Bmod, Phi_mn, Phi_aligned_mn) + subroutine compute_flr2_coeff(flr2) + use mephit_conf, only: conf + use mephit_util, only: ev2erg + use mephit_mesh, only: equil, mesh, fs, dens_e, temp_e, temp_i, Phi0, dPhi0_dpsi, nu_i, nu_e + use mephit_flr2, only: flr2_t, flr2_init, flr2_coeff + type(flr2_t), intent(inout) :: flr2 + integer :: m_min, m_max + + m_min = min(-equil%cocos%sgn_q * mesh%m_res_min, -equil%cocos%sgn_q * mesh%m_res_max) + m_max = max(-equil%cocos%sgn_q * mesh%m_res_min, -equil%cocos%sgn_q * mesh%m_res_max) + call flr2_init(flr2, mesh%nflux, m_min, m_max) + call flr2_coeff(flr2, m_min, m_max, mesh%n, mesh%nflux, conf%m_i, conf%Z_i, & + mesh%R_O, fs%psi(1:), fs%q(1:), fs%F(1:), Phi0%y(1:), mesh%avg_R2gradpsi2(1:), & + dens_e%y(1:), temp_e%y(1:) * ev2erg, temp_i%y(1:) * ev2erg, nu_e%y(1:), nu_i%y(1:)) + end subroutine compute_flr2_coeff + + subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_mn) use mephit_conf, only: conf - use mephit_util, only: imun, ev2erg, resample1d - use mephit_mesh, only: equil, mesh, fs, fs_half, & - dens_e, temp_e, temp_i, Phi0, dPhi0_dpsi, nu_i, nu_e + use mephit_util, only: imun, resample1d + use mephit_mesh, only: equil, mesh, fs, fs_half, dPhi0_dpsi use mephit_pert, only: polmodes_t, & vec_polmodes_t, vec_polmodes_init, vec_polmodes_deinit, RT0_poloidal_modes + use mephit_flr2, only: flr2_t, flr2_response_current type(RT0_t), intent(in) :: Bn type(polmodes_t), intent(inout) :: jmnpar_over_Bmod - complex(dp), intent(out) :: Phi_mn(:, mesh%m_res_min:) - complex(dp), intent(out) :: Phi_aligned_mn(:, mesh%m_res_min:) + type(flr2_t), intent(in) :: flr2 + complex(dp), intent(out) :: Phi_mn(0:, mesh%m_res_min:) + complex(dp), intent(out) :: Phi_aligned_mn(0:, mesh%m_res_min:) integer :: m, m_res, kf_min complex(dp) :: Bmnpsi_over_B0phi(0:mesh%nflux) type(vec_polmodes_t) :: Bmn @@ -1269,15 +1300,13 @@ subroutine compute_kilca_current(Bn, jmnpar_over_Bmod, Phi_mn, Phi_aligned_mn) equil%cocos%sgn_dpsi, fs%psi(1:), Bmnpsi_over_B0phi(1:)%Re, 3) call resample1d(fs_half%psi, Bmn%coeff_rad(m_res, :)%Im / fs_half%q * & equil%cocos%sgn_dpsi, fs%psi(1:), Bmnpsi_over_B0phi(1:)%Im, 3) - call response_current(1, m_res, mesh%n, mesh%nflux, conf%m_i, conf%Z_i, & - mesh%R_O, fs%psi(1:), fs%q(1:), fs%F(1:), Phi0%y(1:), mesh%avg_R2gradpsi2(1:), & - dens_e%y(1:), temp_e%y(1:) * ev2erg, temp_i%y(1:) * ev2erg, nu_e%y(1:), nu_i%y(1:), & + call flr2_response_current(flr2, 1, m_res, mesh%nflux, & Bmnpsi_over_B0phi(1:), jmnpar_over_Bmod%coeff(m_res, 1:), Phi_mn(1:, m)) Phi_aligned_mn(:, m) = imun * Bmnpsi_over_B0phi * fs%q * dPhi0_dpsi%y / (m_res + mesh%n * fs%q) jmnpar_over_Bmod%coeff(m_res, :kf_min) = (0d0, 0d0) ! suppress spurious current near axis end do call vec_polmodes_deinit(Bmn) - end subroutine compute_kilca_current + end subroutine compute_flr2_current subroutine compute_shielding_current(pn, jmnpar_over_Bmod) use mephit_conf, only: conf, conf_arr @@ -1411,16 +1440,18 @@ subroutine resonance_sweep(file, group) use mephit_mesh, only: equil, mesh, cache, fs, dens_e, temp_e, temp_i, Phi0, nu_i, nu_e, E_r use mephit_pert, only: L1_t, L1_init, L1_deinit, polmodes_t, polmodes_init, polmodes_deinit, & L1_sum_poloidal_modes, compute_Ires + use mephit_flr2, only: flr2_t, flr2_init, flr2_deinit, flr2_coeff, flr2_response_current character(len = *), intent(in) :: file character(len = *), intent(in) :: group character(len = len_trim(group)) :: grp type(polmodes_t) :: jmnpar_over_Bmod type(L1_t) :: jnpar_over_Bmod + type(flr2_t) :: flr2 real(dp) :: Delta_E_r(conf%resonance_sweep), E_r_zero(conf%resonance_sweep) complex(dp) :: Bmnpsi_over_B0phi(0:mesh%nflux), & Phi_mn(0:mesh%nflux, mesh%m_res_min:mesh%m_res_max), & Imn_res(mesh%m_res_min:mesh%m_res_max, conf%resonance_sweep) - integer :: k, m, m_res, kf_min, kf_max + integer :: k, m, m_res, m_min, m_max, kf_min, kf_max integer(HID_T) :: h5id_root if (logger%info) then @@ -1432,16 +1463,20 @@ subroutine resonance_sweep(file, group) call polmodes_init(jmnpar_over_Bmod, conf%m_max, mesh%nflux) jmnpar_over_Bmod%coeff(:, :) = (0d0, 0d0) call L1_init(jnpar_over_Bmod, mesh%npoint) + m_min = min(-equil%cocos%sgn_q * mesh%m_res_min, -equil%cocos%sgn_q * mesh%m_res_max) + m_max = max(-equil%cocos%sgn_q * mesh%m_res_min, -equil%cocos%sgn_q * mesh%m_res_max) + call flr2_init(flr2, mesh%nflux, m_min, m_max) kf_min = mesh%res_ind(mesh%m_res_min) - 1 kf_max = mesh%res_ind(mesh%m_res_max) Delta_E_r = linspace(-E_r%y(kf_min), -E_r%y(kf_max), conf%resonance_sweep, 0, 0) do k = 1, conf%resonance_sweep E_r_zero(k) = zeroin(fs%psi(kf_min), fs%psi(kf_max), E_r_shifted, 1d-9) + call flr2_coeff(flr2, m_min, m_max, mesh%n, mesh%nflux, conf%m_i, conf%Z_i, & + mesh%R_O, fs%psi(1:), fs%q(1:), fs%F(1:), Phi0%y(1:) - fs%rsmall(1:) * Delta_E_r(k), mesh%avg_R2gradpsi2(1:), & + dens_e%y(1:), temp_e%y(1:) * ev2erg, temp_i%y(1:) * ev2erg, nu_e%y(1:), nu_i%y(1:)) do m = mesh%m_res_min, mesh%m_res_max m_res = -equil%cocos%sgn_q * m - call response_current(1, m_res, mesh%n, mesh%nflux, conf%m_i, conf%Z_i, & - mesh%R_O, fs%psi(1:), fs%q(1:), fs%F(1:), Phi0%y(1:) - fs%rsmall(1:) * Delta_E_r(k), mesh%avg_R2gradpsi2(1:), & - dens_e%y(1:), temp_e%y(1:) * ev2erg, temp_i%y(1:) * ev2erg, nu_e%y(1:), nu_i%y(1:), & + call flr2_response_current(flr2, 1, m_res, mesh%nflux, & Bmnpsi_over_B0phi(1:), jmnpar_over_Bmod%coeff(m_res, 1:), Phi_mn(1:, m)) end do call L1_sum_poloidal_modes(jmnpar_over_Bmod, jnpar_over_Bmod) @@ -1452,6 +1487,7 @@ subroutine resonance_sweep(file, group) end do call polmodes_deinit(jmnpar_over_Bmod) call L1_deinit(jnpar_over_Bmod) + call flr2_deinit(flr2) grp = trim(group) call h5_open_rw(file, h5id_root) call h5_create_parent_groups(h5id_root, grp // '/') diff --git a/src/mephit_util.f90 b/src/mephit_util.f90 index 3c6bb23..85b967b 100644 --- a/src/mephit_util.f90 +++ b/src/mephit_util.f90 @@ -22,12 +22,15 @@ module mephit_util C_F_string ! module variables - public :: pi, clight, ev2erg, imun - - real(dp), parameter :: pi = 4d0 * atan(1d0) - real(dp), parameter :: clight = 2.99792458d10 !< Speed of light in cm sec^-1. - real(dp), parameter :: ev2erg = 1.6021766d-12 !< convert eV to erg - complex(dp), parameter :: imun = (0.0_dp, 1.0_dp) !< Imaginary unit in double precision. + public :: imun, pi, clight, ev2erg, m_e, m_p, elem_charge + + complex(dp), parameter :: imun = (0.0_dp, 1.0_dp) !< imaginary unit in double precision + real(dp), parameter :: pi = 4d0 * atan(1d0) !< pi in double precision + real(dp), parameter :: clight = 2.99792458d10 !< speed of light in cm sec^-1 + real(dp), parameter :: ev2erg = 1.6021766d-12 !< convert eV to erg + real(dp), parameter :: m_e = 9.1093837139d-28 !< electron mass in g + real(dp), parameter :: m_p = 1.67262192595d-24 !< proton mass in g + real(dp), parameter :: elem_charge = 4.80320425d-10 !< elementary charge in statC ! types and interfaces