Skip to content

Minor changes Patrick #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pazathoth merged 6 commits into from
Apr 29, 2025
Merged

Minor changes Patrick #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
95f0879
remove .gitkeep in conflict with symlink
pazathoth Apr 17, 2025
e616949
use explicit temporary arrays
pazathoth Apr 17, 2025
07f492e
remove obsolete shielding projection auxiliaries
pazathoth Apr 24, 2025
8a734bf
specify FLR2 coefficient units
pazathoth Apr 24, 2025
769df29
default shielding factor to 1, apply only if refined
pazathoth Apr 24, 2025
64951ea
inline finite element projections
pazathoth Apr 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 21 additions & 6 deletions data/mephit_g000001.0001_TCFP_hip.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,29 @@
config%refinement_scheme = 0
! 0: iMHD, 1: FLRE 2nd order
config%currn_model = 0
! 0: Viktor Nemov's code, 1: GPEC, 2: pre-computed Fourier transform
config%vac_src = 2
! prefactor in Biot-Savart law, e.g. 5 windings and conversion from A to statA / c_0
config%Biot_Savart_prefactor = 5.0d-1
!> file containing coil field data
config%coil_file = '/proj/plasma/CODE/lainer_p/MEPHIT/run/AUG_B_coils.h5'
!> file containing coil currents
config%currents_file = '/proj/plasma/DATA/BALANCE/COIL/33353/33353.2900_coil_markl.dat'
!> file containing equilibrium density profile
config%dens_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_ne_PED_MMARKL_rho_pol.dat'
!> file containing equilibrium electron temperature profile
config%temp_e_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Te_PED_MMARKL_rho_pol.dat'
!> file containing equilibrium ion temperature profile
config%temp_i_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Ti_PED_MMARKL_rho_pol.dat'
!> file containing equilibrium radial electric field profile
config%E_r_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353_2900_Er.dat'
! mass number of ions
config%m_i = 2.0
! toroidal mode number of perturbation
config%n = 2
! maximum number of iteration steps
! maximum number of (preconditioned) fixed-point iteration steps
config%niter = 50
! maximum relative error in fixed-point iteration
! maximum relative error in (preconditioned) fixed-point iteration
config%iter_rel_err = 1d-12
! maximum number of iterations in Arnoldi method
config%nkrylov = 100
Expand Down Expand Up @@ -49,8 +67,5 @@
&arrays
! relative refinement around resonances
refinement(4) = 1.414213562373095d0
! deleted flux surfaces around resonances
deletions(4) = 1
! coefficients for sheet currents
sheet_current_factor(4) = 1d4
Delta_rad_res(4) = 0.1125d0
/
9 changes: 0 additions & 9 deletions data/mephit_g33353_2900_EQH.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,13 +90,4 @@
Delta_rad_res(9) = 5.347772320733589e-02
Delta_rad_res(10) = 6.066328245075899e-02
Delta_rad_res(11) = 3.464590203985631e-02
sheet_current_factor(3) = 1d5
sheet_current_factor(4) = 1d5
sheet_current_factor(5) = 1d5
sheet_current_factor(6) = 1d5
sheet_current_factor(7) = 1d5
sheet_current_factor(8) = 1d5
sheet_current_factor(9) = 1d5
sheet_current_factor(10) = 0d0
sheet_current_factor(11) = 0d0
/
50 changes: 24 additions & 26 deletions data/mephit_illustration.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,28 +3,24 @@
config%log_level = 4
! 0: single iteration, 1: direct iterations, 2: preconditioned iterations
config%runmode = 2
! use non-resonant test case
config%nonres = .false.
! average over quadrilaterals for non-resonant test case
config%quad_avg = .true.
! 0: EPS paper profile, 1: parabolic profile, 2: EQDSK profile
! 0: EPS 2018 paper profile, 1: parabolic profile, 2: EQDSK profile
config%pres_prof = 2
! 0: only PS current, 1: rot B_0, 2: EQDSK profile
config%curr_prof = 2
! 0: surface integral, 1: field line integral, 2: EQDSK profile
config%q_prof = 1
! 0: iMHD, 1: KiLCA interface
! 0: geometric series around resonance, 1: Gaussian profile
config%refinement_scheme = 0
! 0: iMHD, 1: FLRE 2nd order
config%currn_model = 1
! 0: Viktor Nemov's code, 1: GPEC, 2: pre-computed Fourier transform
config%vac_src = 2
! mass number of ions
config%m_i = 2.0
! prefactor in Biot-Savart law, e.g. 5 windings and conversion from A to statA / c_0
config%Biot_Savart_prefactor = 5.0d-1
!> file containing coil field data
config%coil_file = '/proj/plasma/CODE/lainer_p/MEPHIT/run/AUG_B_coils.h5'
!> file containing coil currents
config%currents_file = '/proj/plasma/DATA/BALANCE/COIL/33353/33353.2900_coil_markl.dat'
!> file containing coil geometries
config%coil_file = '/temp/lainer_p/MEPHIT-run/AUG_B_coils.h5'
!> file containing equilibrium density profile
config%dens_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_ne_PED_MMARKL_rho_pol.dat'
!> file containing equilibrium electron temperature profile
Expand All @@ -33,11 +29,13 @@
config%temp_i_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Ti_PED_MMARKL_rho_pol.dat'
!> file containing equilibrium radial electric field profile
config%E_r_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353_2900_Er.dat'
! mass number of ions
config%m_i = 2.0
! toroidal mode number of perturbation
config%n = 2
! maximum number of iteration steps
! maximum number of (preconditioned) fixed-point iteration steps
config%niter = 75
! maximum relative error in fixed-point iteration
! maximum relative error in (preconditioned) fixed-point iteration
config%iter_rel_err = 1d-9
! maximum number of iterations in Arnoldi method
config%nkrylov = 100
Expand All @@ -47,13 +45,22 @@
config%ritz_rel_err = 1d-9
! maximum distance between flux surfaces along theta = 0
config%max_Delta_rad = 1.75d0
! enable damping of Pfirsch-Schlueter current
! maximum number of points per flux surface (0: no maximum)
config%pol_max = 0
! minimum number of points per flux surface (0: no minimum)
config%pol_min = 0
! ignore resonance position where q = 1, which is usually spurious
config%ignore_q1_res = .true.
! enable damping of Pfirsch-Schlueter current near resonances
config%damp = .true.
! number of points in resonance sweeep
! maximal number of eigenvectors to be saved
config%max_eig_out = 10
! maximum poloidal mode number for Fourier transform of results
config%m_max = 24
! number of points in sweep over electron fluid resonance
config%resonance_sweep = 0

config%debug_mfem = .false.
config%pol_max = 0
! offset added to radial electric field
config%offset_E_r = 0.0d0
/
&arrays
refinement(3) = 1.618d0
Expand Down Expand Up @@ -83,13 +90,4 @@
Delta_rad_res(9) = 5.347772320733589e-02
Delta_rad_res(10) = 6.066328245075899e-02
Delta_rad_res(11) = 3.464590203985631e-02
sheet_current_factor(3) = 1d5
sheet_current_factor(4) = 1d5
sheet_current_factor(5) = 1d5
sheet_current_factor(6) = 1d5
sheet_current_factor(7) = 1d5
sheet_current_factor(8) = 1d5
sheet_current_factor(9) = 1d5
sheet_current_factor(10) = 0d0
sheet_current_factor(11) = 0d0
/
Empty file removed run/.gitkeep
View file
Open in desktop
Empty file.
6 changes: 3 additions & 3 deletions src/mephit_conf.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ module mephit_conf

integer :: m_min, m_max

!> Width of refined flux surfaces around resonances in cm.
!> Width of refined flux surfaces around resonances in cm. Defaults to 0.
real(dp), dimension(:), allocatable :: Delta_rad_res

!> Number of additional fine flux surfaces. Defaults to 0.
Expand All @@ -237,7 +237,7 @@ module mephit_conf
!> Width ratio of neighbouring refined flux surfaces. Defaults to 0 (no refinement).
real(dp), dimension(:), allocatable :: refinement

!> Free parameters setting the magnitudes of sheet currents. Defaults to 0.
!> Free parameters setting the magnitudes of sheet currents. Only applied with refinement. Defaults to 1.
real(dp), dimension(:), allocatable :: sheet_current_factor

contains
Expand Down Expand Up @@ -372,7 +372,7 @@ subroutine config_delayed_read(config, filename, m_min, m_max)
Delta_rad_res = 0d0
add_fine = 0
refinement = 0d0
sheet_current_factor = 0d0
sheet_current_factor = 1d0
open(newunit = fid, file = filename)
read(fid, nml = arrays)
close(fid)
Expand Down
28 changes: 13 additions & 15 deletions src/mephit_flr2.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,44 +102,42 @@ subroutine flr2_write(flr2, file, dataset)
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 = '?', &
flr2%a2_in, lbound(flr2%a2_in), ubound(flr2%a2_in), unit = '1', &
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 = '?', &
flr2%a2_out, lbound(flr2%a2_out), ubound(flr2%a2_out), unit = '1', &
comment = 'Poisson equation coefficient')
call h5_add(h5id_root, trim(adjustl(dataset)) // '/a0', &
flr2%a0, lbound(flr2%a0), ubound(flr2%a0), unit = '?', &
flr2%a0, lbound(flr2%a0), ubound(flr2%a0), unit = 'Mx^-2', &
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 = '?', &
flr2%b2_in, lbound(flr2%b2_in), ubound(flr2%b2_in), unit = 'statV Mx^-1', &
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 = '?', &
flr2%b2_out, lbound(flr2%b2_out), ubound(flr2%b2_out), unit = 'statV Mx^-1', &
comment = 'Poisson equation coefficient')
call h5_add(h5id_root, trim(adjustl(dataset)) // '/b0', &
flr2%b0, lbound(flr2%b0), ubound(flr2%b0), unit = '?', &
flr2%b0, lbound(flr2%b0), ubound(flr2%b0), unit = 'statV Mx^-3', &
comment = 'Poisson equation coefficient')
! note: statA with c_0 = 1
call h5_add(h5id_root, trim(adjustl(dataset)) // '/c2_in', &
flr2%c2_in, lbound(flr2%c2_in), ubound(flr2%c2_in), unit = '?', &
flr2%c2_in, lbound(flr2%c2_in), ubound(flr2%c2_in), unit = 'statA statV^-1 Mx', &
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 = '?', &
flr2%c2_out, lbound(flr2%c2_out), ubound(flr2%c2_out), unit = 'statA statV^-1 Mx', &
comment = 'Poisson equation coefficient')
call h5_add(h5id_root, trim(adjustl(dataset)) // '/c0', &
flr2%c0, lbound(flr2%c0), ubound(flr2%c0), unit = '?', &
flr2%c0, lbound(flr2%c0), ubound(flr2%c0), unit = 'statA statV^-1 M^-1', &
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 = '?', &
flr2%d2_in, lbound(flr2%d2_in), ubound(flr2%d2_in), unit = 'statA', &
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 = '?', &
flr2%d2_out, lbound(flr2%d2_out), ubound(flr2%d2_out), unit = 'statA', &
comment = 'Poisson equation coefficient')
call h5_add(h5id_root, trim(adjustl(dataset)) // '/d0', &
flr2%d0, lbound(flr2%d0), ubound(flr2%d0), unit = '?', &
flr2%d0, lbound(flr2%d0), ubound(flr2%d0), unit = 'statA Mx^-2', &
comment = 'Poisson equation coefficient')
call h5_close(h5id_root)
end subroutine flr2_write
Expand Down
25 changes: 15 additions & 10 deletions src/mephit_iter.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1063,8 +1063,7 @@ 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
L1_t, L1_init, L1_deinit, RT0_t, RT0_init, RT0_deinit
use mephit_flr2, only: flr2_t
type(perteq_t), intent(inout) :: perteq
type(fdm_t), intent(in) :: fdm
Expand Down Expand Up @@ -1285,7 +1284,8 @@ subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_
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)
real(dp) :: Bmnrho_over_B0theta(1:mesh%nflux), rtemp(1:mesh%nflux)
complex(dp) :: Bmnpsi_over_B0phi(1:mesh%nflux)
type(vec_polmodes_t) :: Bmn

kf_min = 0 ! max(1, mesh%res_ind(mesh%m_res_min) / 4)
Expand All @@ -1297,13 +1297,18 @@ subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_
call RT0_poloidal_modes(Bn, Bmn)
do m = mesh%m_res_min, mesh%m_res_max
m_res = -equil%cocos%sgn_q * m
call resample1d(fs_half%psi, Bmn%coeff_rad(m_res, :)%Re / fs_half%q * &
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)
Bmnrho_over_B0theta(:) = Bmn%coeff_rad(m_res, :)%Re
call resample1d(fs_half%psi, Bmnrho_over_B0theta / fs_half%q * &
equil%cocos%sgn_dpsi, fs%psi(1:), rtemp, 3)
Bmnpsi_over_B0phi%Re = rtemp
Bmnrho_over_B0theta(:) = Bmn%coeff_rad(m_res, :)%Im
call resample1d(fs_half%psi, Bmnrho_over_B0theta / fs_half%q * &
equil%cocos%sgn_dpsi, fs%psi(1:), rtemp, 3)
Bmnpsi_over_B0phi%Im = rtemp
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)
Bmnpsi_over_B0phi, jmnpar_over_Bmod%coeff(m_res, 1:), Phi_mn(1:, m))
Phi_aligned_mn(1:, m) = imun * Bmnpsi_over_B0phi * &
fs%q(1:) * dPhi0_dpsi%y(1:) / (m_res + mesh%n * fs%q(1:))
jmnpar_over_Bmod%coeff(m_res, :kf_min) = (0d0, 0d0) ! suppress spurious current near axis
end do
call vec_polmodes_deinit(Bmn)
Expand All @@ -1323,7 +1328,7 @@ subroutine compute_shielding_current(pn, jmnpar_over_Bmod)
call L1_poloidal_modes(pn, pmn)
do m = mesh%m_res_min, mesh%m_res_max
m_res = -equil%cocos%sgn_q * m
if (abs(conf_arr%sheet_current_factor(m)) > 0d0) then
if (abs(conf_arr%refinement(m)) > 1d0) then
kf = mesh%res_ind(m)
jmnpar_over_Bmod%coeff(m_res, kf-1:kf) = &
(pmn%coeff(m_res, kf) - pmn%coeff(m_res, kf-1)) * &
Expand Down
Loading
Loading