From 72059366fbcb142f9ba9fc656520dd55244ee743 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 13 Mar 2025 15:34:23 +0100 Subject: [PATCH 1/9] modularize and lint FLR2 --- src/mephit_flr2.f90 | 1049 ++++++++++++++++++++++--------------------- src/mephit_iter.F90 | 2 + 2 files changed, 529 insertions(+), 522 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 4c6ce37..0fd31d9 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -1,535 +1,540 @@ -! +module mephit_flr2 + +contains + 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 -! - 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) & - * (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) & - * (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) -! + 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 + + 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) & + * (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) & + * (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 + !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 - symbI=Imn + ! 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 + 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 -!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 - -! Date: Thu, 13 Mar 2025 15:39:03 +0100 Subject: [PATCH 2/9] put implicit none and public/private at top of FLR2 module --- src/mephit_flr2.f90 | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 0fd31d9..d8f7aed 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -1,4 +1,9 @@ module mephit_flr2 + implicit none + + private + + public :: response_current contains @@ -34,8 +39,6 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & ! 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 - implicit none - integer, parameter :: mnmax=3 double precision, parameter :: pi=3.14159265358979d0 double precision, parameter :: c=2.9979d10 @@ -251,8 +254,6 @@ end subroutine getIfunc subroutine W2_arr (x1_in,x2_in,Imn) - implicit none; - 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') @@ -426,8 +427,6 @@ subroutine progonka(isw_f,npoi,x,q, & ! Output: ! phi(npoi) - (double complex) array of the solution $\Phi$ - implicit none - integer :: isw_f,npoi,i double precision :: dxp,dxm,dxt @@ -517,8 +516,6 @@ subroutine first_deriv(npoi,x,f,df_dx) ! df_dx(npoi) - (double precision) array of derivative values - implicit none - integer :: npoi,i double precision :: dxp,dxm,dxt From 018862d0a46d1e7c2c1b89eb0efa4ef727ab11ef Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 13 Mar 2025 15:42:38 +0100 Subject: [PATCH 3/9] use standard types in FLR2 module --- src/mephit_flr2.f90 | 122 ++++++++++++++++++++++---------------------- 1 file changed, 62 insertions(+), 60 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index d8f7aed..862ba64 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -1,4 +1,6 @@ module mephit_flr2 + use iso_fortran_env, only: dp => real64 + implicit none private @@ -20,53 +22,53 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & ! 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 + ! 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 + ! 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) - (double complex) parallel current density divided by uperturbed magnetic field module - ! Phi_m(npoi) - (double complex) perturbation of the electrostatic potential + ! parcur_over_b0(npoi) - (complex) parallel current density divided by uperturbed magnetic field module + ! Phi_m(npoi) - (complex) perturbation of the electrostatic potential 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) + real(dp), parameter :: pi=3.14159265358979d0 + real(dp), parameter :: c=2.9979d10 + real(dp), parameter :: e_charge=4.8032d-10 + real(dp), parameter :: e_mass=9.1094d-28 + real(dp), parameter :: p_mass=1.6726d-24 + real(dp), parameter :: ev=1.6022d-12 + complex(dp), 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, & + real(dp) :: 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(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 + complex(dp), 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 + complex(dp), dimension(0:mnmax,0:mnmax) :: symbI + real(dp), dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i + complex(dp), dimension(:), allocatable :: a2_in,a2_out,a0 + complex(dp), dimension(:), allocatable :: b2_in,b2_out,b0 + complex(dp), dimension(:), allocatable :: c2_in,c2_out,c0 + complex(dp), 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)) @@ -201,9 +203,9 @@ end subroutine response_current 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 + real(dp) :: x1,x2,z + complex(dp) :: denom + complex(dp), dimension(0:mnmax,0:mnmax) :: symbI,Imn ! if(.true.) then if(.false.) then @@ -277,7 +279,7 @@ end subroutine hypergeometric1f1_kummer_modified_0_ada complex(dpc), dimension(0:3,0:3) :: Imn complex(dpc) :: t1, t2, F11m, x1, x2 !complex(8), parameter :: I = cmplx(0.0d0, 1.0d0, 8), one = cmplx(1.0d0, 0.0d0, 8); - double complex, parameter :: I = (0.0d0, 1.0d0), one = (1.0d0, 0.0d0) + complex(dp), parameter :: I = (0.0d0, 1.0d0), one = (1.0d0, 0.0d0) integer :: l,nmax,m,n; @@ -419,34 +421,34 @@ subroutine progonka(isw_f,npoi,x,q, & ! 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) - (double precision) array of independent variable, need not to be equidistant - ! c_1(npoi) - (double complex) array of coefficient $C_1$ - ! c_2(npoi) - (double complex) array of coefficient $C_2$ - ! b(npoi) - (double complex) array of source $B$ + ! 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) - (double complex) array of the solution $\Phi$ + ! phi(npoi) - (complex) array of the solution $\Phi$ integer :: isw_f,npoi,i - double precision :: dxp,dxm,dxt - double complex :: denom + real(dp) :: dxp,dxm,dxt + complex(dp) :: denom - double precision, dimension(npoi) :: x - double complex, dimension(npoi) :: a2_in,a2_out,a0, & + real(dp), dimension(npoi) :: x + complex(dp), dimension(npoi) :: a2_in,a2_out,a0, & b2_in,b2_out,b0, & c2_in,c2_out,c0, & d2_in,d2_out,d0, & q,f,g - double complex, dimension(:), allocatable :: alp,bet,quelle - double complex, dimension(:,:), allocatable :: wsecder,eqmat + complex(dp), dimension(:), allocatable :: alp,bet,quelle + complex(dp), dimension(:,:), allocatable :: wsecder,eqmat 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 + ! 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 @@ -509,18 +511,18 @@ subroutine first_deriv(npoi,x,f,df_dx) ! 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 + ! x(npoi) - (real) array of independent variable, need not to be equidistant + ! f(npoi) - (real) array of function values ! Output: - ! df_dx(npoi) - (double precision) array of derivative values + ! df_dx(npoi) - (real) array of derivative values integer :: npoi,i - double precision :: dxp,dxm,dxt + real(dp) :: dxp,dxm,dxt - double precision, dimension(npoi) :: x,f,df_dx + real(dp), dimension(npoi) :: x,f,df_dx do i=2,npoi-1 dxp=x(i+1)-x(i) From 4fe952b05cbddd35617d247ed6c65559f4d62fe5 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 13 Mar 2025 16:02:20 +0100 Subject: [PATCH 4/9] clean parameters & separate interfaces in FLR2 --- src/mephit_flr2.f90 | 64 ++++++++++++++++++--------------------------- src/mephit_util.f90 | 15 ++++++----- 2 files changed, 34 insertions(+), 45 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 862ba64..8435207 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -7,6 +7,24 @@ module mephit_flr2 public :: response_current + integer, parameter :: mnmax = 3 + + 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 + contains subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & @@ -41,14 +59,7 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & ! parcur_over_b0(npoi) - (complex) parallel current density divided by uperturbed magnetic field module ! Phi_m(npoi) - (complex) perturbation of the electrostatic potential - integer, parameter :: mnmax=3 - real(dp), parameter :: pi=3.14159265358979d0 - real(dp), parameter :: c=2.9979d10 - real(dp), parameter :: e_charge=4.8032d-10 - real(dp), parameter :: e_mass=9.1094d-28 - real(dp), parameter :: p_mass=1.6726d-24 - real(dp), parameter :: ev=1.6022d-12 - complex(dp), parameter :: imun=(0.d0,1.d0) + use mephit_util, only: imun, pi, c => clight, e_charge => elem_charge, e_mass => m_e, p_mass => m_p integer :: isw_Phi_m,mpol,ntor,npoi,i real(dp) :: am_i,z_i,Rtor @@ -201,7 +212,6 @@ end subroutine response_current !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine getIfunc(x1,x2,symbI) - integer, parameter :: mnmax=3 integer :: m,n real(dp) :: x1,x2,z complex(dp) :: denom @@ -255,47 +265,23 @@ end subroutine getIfunc !------------------------------------------------------------------------------ subroutine W2_arr (x1_in,x2_in,Imn) - - 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 :: dp = 8 - integer, parameter :: dpc = 8 - real(dp), parameter :: pi = 3.141592653589793238462643383279502884197_dp real(dp) :: x1_in,x2_in; - complex(dpc), dimension(0:3,0:3) :: Imn - complex(dpc) :: t1, t2, F11m, x1, x2 - !complex(8), parameter :: I = cmplx(0.0d0, 1.0d0, 8), one = cmplx(1.0d0, 0.0d0, 8); + complex(dp), dimension(0:3,0:3) :: Imn + complex(dp) :: t1, t2, F11m, x1, x2 complex(dp), parameter :: I = (0.0d0, 1.0d0), one = (1.0d0, 0.0d0) - integer :: l,nmax,m,n; + integer :: l,m,n; real(dp) :: F_im, F_re; - complex(dpc), allocatable, dimension(:,:,:) :: W2; - - nmax=3 + complex(dp), allocatable, dimension(:,:,:) :: W2; - !allocate (W2(0:1, 0:3, -Nmax:Nmax)) - allocate (W2(0:3, 0:3, 0:nmax)) + allocate (W2(0:3, 0:3, 0:mnmax)) W2 = (0.0d0, 0.0d0) x1 = dcmplx(x1_in,0.d0) t1 = x1**2 - do l = 0, nmax ! does not work for 3, 2 and 1! + do l = 0, mnmax ! does not work for 3, 2 and 1! x2 = x2_in + I*l 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 From acdbd9ebc39644d2795757c7bea127869345061c Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 13 Mar 2025 16:15:26 +0100 Subject: [PATCH 5/9] minor cleanup in FLR2 module --- src/mephit_flr2.f90 | 58 +++++++++++++++++++-------------------------- 1 file changed, 25 insertions(+), 33 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 8435207..db84306 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -190,7 +190,7 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & c0(i) = d0(i)*factor_of_Phi c2_in(i) = d2_in(i)*factor_of_Phi c2_out(i) = c2_in(i) - enddo + end do ! Solve Laplace equation: @@ -210,7 +210,7 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & end subroutine response_current - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + subroutine getIfunc(x1,x2,symbI) integer :: m,n real(dp) :: x1,x2,z @@ -247,38 +247,33 @@ subroutine getIfunc(x1,x2,symbI) 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 + end do + end do 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)) + end if - endif + end if end subroutine getIfunc - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - ! Calculates array of W2 special functions. + subroutine W2_arr(x1_in,x2_in,Imn) + real(dp) :: x1_in,x2_in complex(dp), dimension(0:3,0:3) :: 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; + integer :: l,m,n - real(dp) :: F_im, F_re; - complex(dp), allocatable, dimension(:,:,:) :: W2; + real(dp) :: F_im, F_re + complex(dp), allocatable, dimension(:,:,:) :: W2 allocate (W2(0:3, 0:3, 0:mnmax)) W2 = (0.0d0, 0.0d0) - x1 = dcmplx(x1_in,0.d0) + x1 = cmplx(x1_in, 0.d0, dp) t1 = x1**2 do l = 0, mnmax ! does not work for 3, 2 and 1! @@ -286,15 +281,15 @@ subroutine W2_arr (x1_in,x2_in,Imn) x2 = x2_in + I*l - t2 = - I*x2 + t1; + t2 = - I*x2 + t1 - ! call Hypergeometric1F1_kummer_modified_0_ada (real(one+t2), aimag(one+t2), real(t1), aimag(t1), F_re, F_im); + ! 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; + ! 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); + 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; + F11m = F_re + I*F_im !asymptotic form of expressions: has no fake singularity at kp=0 @@ -362,13 +357,11 @@ subroutine W2_arr (x1_in,x2_in,Imn) do m=0,3 do n=0,m-1 Imn(m,n)=Imn(n,m) - enddo - enddo + end do + end do end subroutine W2_arr - !------------------------------------------------------------------------------ - subroutine progonka(isw_f,npoi,x,q, & a2_in,a2_out,a0, & @@ -452,7 +445,7 @@ subroutine progonka(isw_f,npoi,x,q, & eqmat(:,i) = a2_out(i)*wsecder(:,i) & + a2_in(i-1:i+1)*wsecder(:,i) eqmat(0,i) = eqmat(0,i) + a0(i) - enddo + end do eqmat(0,1)=(1.d0,0.d0) eqmat(0,npoi)=(1.d0,0.d0) @@ -464,13 +457,13 @@ subroutine progonka(isw_f,npoi,x,q, & 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 - enddo + 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) - enddo + end do f=f*dble(isw_f) @@ -483,13 +476,12 @@ subroutine progonka(isw_f,npoi,x,q, & + sum(d2_in(i-1:i+1)*q(i-1:i+1)*wsecder(:,i)) & + d0(i)*q(i) - enddo + end do deallocate(wsecder,eqmat,alp,bet,quelle) end subroutine progonka - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine first_deriv(npoi,x,f,df_dx) @@ -515,7 +507,7 @@ subroutine first_deriv(npoi,x,f,df_dx) 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 + 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)) From 1192d91f56bb69fae54b294b7cb72f080abf6941 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 13 Mar 2025 16:36:04 +0100 Subject: [PATCH 6/9] declare intent in FLR2 module --- src/mephit_flr2.f90 | 72 +++++++++++++++++++++++++-------------------- 1 file changed, 40 insertions(+), 32 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index db84306..8c641cb 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -27,10 +27,10 @@ end subroutine hypergeometric1f1_kummer_modified_0_ada contains - 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) + 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 @@ -61,8 +61,9 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & use mephit_util, only: imun, pi, c => clight, e_charge => elem_charge, e_mass => m_e, p_mass => m_p - integer :: isw_Phi_m,mpol,ntor,npoi,i - real(dp) :: am_i,z_i,Rtor + integer, intent(in) :: isw_Phi_m, mpol, ntor, npoi + integer :: i + 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 @@ -70,9 +71,10 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & complex(dp) :: F_m,Gtor_m,Htor complex(dp) :: F_me,Gtor_me,Htore complex(dp) :: F_mi,Gtor_mi,Htori - real(dp), dimension(npoi) :: psi,qsaf,bcovar_phi,Phi_0,avR2nabpsi2, & - dens_e,temp_e,temp_i,anu_e,anu_i - complex(dp), dimension(npoi) :: bpsi_over_bphi,parcur_over_b0,Phi_m + real(dp), dimension(npoi), intent(in) :: psi, qsaf, bcovar_phi, Phi_0, avR2nabpsi2, & + dens_e, temp_e, temp_i, anu_e, anu_i + complex(dp), dimension(npoi), intent(in) :: bpsi_over_bphi + complex(dp), dimension(npoi), intent(out) :: parcur_over_b0, Phi_m complex(dp), dimension(0:mnmax,0:mnmax) :: symbI real(dp), dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i @@ -211,11 +213,13 @@ subroutine response_current(isw_Phi_m,mpol,ntor,npoi,am_i,z_i,Rtor, & end subroutine response_current - subroutine getIfunc(x1,x2,symbI) - integer :: m,n - real(dp) :: x1,x2,z + 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) :: symbI,Imn + complex(dp), dimension(0:mnmax, 0:mnmax), intent(out) :: symbI + complex(dp), dimension(0:mnmax, 0:mnmax) :: Imn ! if(.true.) then if(.false.) then @@ -259,9 +263,9 @@ end subroutine getIfunc !> Calculates array of W2 special functions. - subroutine W2_arr(x1_in,x2_in,Imn) - real(dp) :: x1_in,x2_in - complex(dp), dimension(0:3,0:3) :: Imn + 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) @@ -363,12 +367,12 @@ subroutine W2_arr(x1_in,x2_in,Imn) 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) + 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$ @@ -408,17 +412,19 @@ subroutine progonka(isw_f,npoi,x,q, & ! Output: ! phi(npoi) - (complex) array of the solution $\Phi$ - integer :: isw_f,npoi,i + integer, intent(in) :: isw_f, npoi + integer :: i real(dp) :: dxp,dxm,dxt complex(dp) :: denom - real(dp), dimension(npoi) :: x - complex(dp), dimension(npoi) :: a2_in,a2_out,a0, & - b2_in,b2_out,b0, & - c2_in,c2_out,c0, & - d2_in,d2_out,d0, & - q,f,g + real(dp), dimension(npoi), intent(in) :: x + complex(dp), dimension(npoi), 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(npoi), intent(out) :: f, g complex(dp), dimension(:), allocatable :: alp,bet,quelle complex(dp), dimension(:,:), allocatable :: wsecder,eqmat @@ -483,7 +489,7 @@ subroutine progonka(isw_f,npoi,x,q, & end subroutine progonka - subroutine first_deriv(npoi,x,f,df_dx) + subroutine first_deriv(npoi, x, f, df_dx) ! Computes first derivative df_dx of function f over variable x @@ -496,11 +502,13 @@ subroutine first_deriv(npoi,x,f,df_dx) ! df_dx(npoi) - (real) array of derivative values - integer :: npoi,i + integer, intent(in) :: npoi + integer :: i real(dp) :: dxp,dxm,dxt - real(dp), dimension(npoi) :: x,f,df_dx + real(dp), dimension(npoi), intent(in) :: x, f + real(dp), dimension(npoi), intent(out) :: df_dx do i=2,npoi-1 dxp=x(i+1)-x(i) From cda21de2ac9d8552f5cc3dda4af56ee84915e40d Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Fri, 14 Mar 2025 15:28:58 +0100 Subject: [PATCH 7/9] check most array sizes in FLR2 module --- src/mephit_flr2.f90 | 77 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 67 insertions(+), 10 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 8c641cb..c84bdf7 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -59,6 +59,7 @@ subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & ! parcur_over_b0(npoi) - (complex) parallel current density divided by uperturbed magnetic field module ! Phi_m(npoi) - (complex) perturbation of the electrostatic potential + 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 integer, intent(in) :: isw_Phi_m, mpol, ntor, npoi @@ -71,10 +72,10 @@ subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & complex(dp) :: F_m,Gtor_m,Htor complex(dp) :: F_me,Gtor_me,Htore complex(dp) :: F_mi,Gtor_mi,Htori - real(dp), dimension(npoi), intent(in) :: psi, qsaf, bcovar_phi, Phi_0, avR2nabpsi2, & + 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(npoi), intent(in) :: bpsi_over_bphi - complex(dp), dimension(npoi), intent(out) :: parcur_over_b0, Phi_m + complex(dp), dimension(:), intent(in) :: bpsi_over_bphi + complex(dp), dimension(:), intent(out) :: parcur_over_b0, Phi_m complex(dp), dimension(0:mnmax,0:mnmax) :: symbI real(dp), dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i @@ -83,6 +84,24 @@ subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & complex(dp), dimension(:), allocatable :: c2_in,c2_out,c0 complex(dp), dimension(:), allocatable :: d2_in,d2_out,d0 + character(len=14), dimension(13), parameter :: array_names = [character(len=14) :: & + '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'] + integer, dimension(13) :: 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), & + size(bpsi_over_bphi), size(parcur_over_b0), size(Phi_m)] + do i = 1, 13 + if (npoi /= array_sizes(i)) then + call logger%msg_arg_size('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 + 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)) @@ -272,9 +291,8 @@ subroutine W2_arr(x1_in, x2_in, Imn) integer :: l,m,n real(dp) :: F_im, F_re - complex(dp), allocatable, dimension(:,:,:) :: W2 + complex(dp), dimension(0:3, 0:3, 0:mnmax) :: W2 - allocate (W2(0:3, 0:3, 0:mnmax)) W2 = (0.0d0, 0.0d0) x1 = cmplx(x1_in, 0.d0, dp) @@ -412,22 +430,44 @@ subroutine progonka(isw_f, npoi, x, q, & ! 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(npoi), intent(in) :: x - complex(dp), dimension(npoi), intent(in) :: q, & + 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(npoi), intent(out) :: f, g + 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 @@ -501,14 +541,31 @@ subroutine first_deriv(npoi, x, f, df_dx) ! 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(npoi), intent(in) :: x, f - real(dp), dimension(npoi), intent(out) :: df_dx + real(dp), dimension(:), intent(in) :: x, f + real(dp), dimension(:), intent(out) :: df_dx + + if (npoi /= size(x)) then + call logger%msg_arg_size('first_deriv', 'npoi', 'size(x)', npoi, size(x)) + if (logger%err) call logger%write_msg + error stop + end if + if (npoi /= size(x)) then + call logger%msg_arg_size('first_deriv', 'npoi', 'size(f)', npoi, size(f)) + if (logger%err) call logger%write_msg + error stop + end if + if (npoi /= size(x)) then + call logger%msg_arg_size('first_deriv', 'npoi', 'size(df_dx)', npoi, size(df_dx)) + if (logger%err) call logger%write_msg + error stop + end if do i=2,npoi-1 dxp=x(i+1)-x(i) From 0e533637990d8492207bd933dd9aaeba67095f5a Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Wed, 19 Mar 2025 17:36:47 +0100 Subject: [PATCH 8/9] compute FLR2 coeffs only once --- src/mephit_flr2.f90 | 427 +++++++++++++++++++++++++++++++------------- src/mephit_iter.F90 | 96 ++++++---- 2 files changed, 366 insertions(+), 157 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index c84bdf7..0f6f8f0 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -5,9 +5,20 @@ module mephit_flr2 private - public :: response_current + ! types and associated procedures + public :: flr2_t, flr2_init, flr2_deinit, flr2_write, flr2_read, & + flr2_coeff, flr2_response_current - integer, parameter :: mnmax = 3 + 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) & @@ -25,45 +36,177 @@ subroutine hypergeometric1f1_kummer_modified_0_ada(b_re, b_im, z_re, z_im, f_re, end subroutine hypergeometric1f1_kummer_modified_0_ada end interface + integer, parameter :: mnmax = 3 + contains - subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & + 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, & - 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 - (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 - ! 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 + 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 - integer, intent(in) :: isw_Phi_m, mpol, ntor, npoi - integer :: i + 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 @@ -74,40 +217,34 @@ subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & 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(:), intent(in) :: bpsi_over_bphi - complex(dp), dimension(:), intent(out) :: parcur_over_b0, Phi_m complex(dp), dimension(0:mnmax,0:mnmax) :: symbI real(dp), dimension(:), allocatable :: dens_i,dPhi_0_dpsi,derpar,A1e,A2e,A1i,A2i - complex(dp), dimension(:), allocatable :: a2_in,a2_out,a0 - complex(dp), dimension(:), allocatable :: b2_in,b2_out,b0 - complex(dp), dimension(:), allocatable :: c2_in,c2_out,c0 - complex(dp), dimension(:), allocatable :: d2_in,d2_out,d0 - character(len=14), dimension(13), parameter :: array_names = [character(len=14) :: & + 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', & - 'bpsi_over_bphi', 'parcur_over_b0', 'Phi_m'] - integer, dimension(13) :: array_sizes + '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), & - size(bpsi_over_bphi), size(parcur_over_b0), size(Phi_m)] - do i = 1, 13 + 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('response_current', & + 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)) - 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): @@ -151,85 +288,126 @@ subroutine response_current(isw_Phi_m, mpol, ntor, npoi, am_i, z_i, Rtor, & A1i=A1i-1.5d0*A2i - do i=1,npoi + 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) + 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: + ! 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 + ! 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: + ! susceptibility functions for electrons: - call getIfunc(x_1,x_2,symbI) + 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)) + 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 + ! 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: + ! susceptibility functions for ions: - call getIfunc(x_1,x_2,symbI) + 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_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 + 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) + 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) - 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) + 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) - 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) + 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) - c0(i) = d0(i)*factor_of_Phi - c2_in(i) = d2_in(i)*factor_of_Phi - c2_out(i) = c2_in(i) + 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 - ! Solve Laplace equation: + 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) - 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) + use mephit_conf, only: logger - 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) + 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 - end subroutine response_current + ! 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) @@ -551,21 +729,18 @@ subroutine first_deriv(npoi, x, f, df_dx) real(dp), dimension(:), intent(in) :: x, f real(dp), dimension(:), intent(out) :: df_dx - if (npoi /= size(x)) then - call logger%msg_arg_size('first_deriv', 'npoi', 'size(x)', npoi, size(x)) - if (logger%err) call logger%write_msg - error stop - end if - if (npoi /= size(x)) then - call logger%msg_arg_size('first_deriv', 'npoi', 'size(f)', npoi, size(f)) - if (logger%err) call logger%write_msg - error stop - end if - if (npoi /= size(x)) then - call logger%msg_arg_size('first_deriv', 'npoi', 'size(df_dx)', npoi, size(df_dx)) - if (logger%err) call logger%write_msg - error stop - end if + 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) diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index 8a977a9..e0e1cf2 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,16 +1255,32 @@ 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: response_current + use mephit_flr2, only: flr2_t, flr2_response_current type(RT0_t), intent(in) :: Bn type(polmodes_t), intent(inout) :: jmnpar_over_Bmod + type(flr2_t), intent(in) :: flr2 complex(dp), intent(out) :: Phi_mn(:, mesh%m_res_min:) complex(dp), intent(out) :: Phi_aligned_mn(:, mesh%m_res_min:) integer :: m, m_res, kf_min @@ -1270,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 @@ -1412,17 +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: response_current + 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 @@ -1434,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) @@ -1454,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 // '/') From 1a1c8124f27af637a2c3dddb3ad78060fbe8dc2e Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 20 Mar 2025 15:04:33 +0100 Subject: [PATCH 9/9] correct indexing for pert. elec. potential --- src/mephit_iter.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index e0e1cf2..255b280 100644 --- a/src/mephit_iter.F90 +++ b/src/mephit_iter.F90 @@ -1281,8 +1281,8 @@ subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_ type(RT0_t), intent(in) :: Bn type(polmodes_t), intent(inout) :: jmnpar_over_Bmod type(flr2_t), intent(in) :: flr2 - complex(dp), intent(out) :: Phi_mn(:, mesh%m_res_min:) - complex(dp), intent(out) :: Phi_aligned_mn(:, mesh%m_res_min:) + 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