module params_equi_dommaschk_m !! Module for parameter reading for the Dommaschk equilibrium type ! [*] W. Dommaschk, "Representations for vacuum potentials in ! stellarators", Computer Physics Communications 40, pg. 203 (1986) use, intrinsic :: iso_fortran_env, only: IOSTAT_END, IOSTAT_EOR use screen_io_m, only: get_stdout use precision_m, only: FP use parallax_build_info_m, only: git_repo_path use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_ERR_PARAMETERS implicit none private integer, parameter :: str_buffer_len = 256 !! Maximum length of character parameters integer, parameter :: DEFAULT_L_POL = 4 !! Default poloidal periodicity upper bound integer, parameter :: DEFAULT_M_TOR = 1 !! Default toroidal periodicity upper bound. These are used to define !! explicit-shape default arrays. ! Common equilibrium parameters real(FP) :: x0 = 1.0_FP !! Radial position of a point on the magnetic axis, normalized real(FP) :: y0 = 0.0_FP !! Vertical position of a point on the magnetic axis, normalized real(FP) :: phi0 = 0.0_FP !! Toroidal position of a point on the magnetic axis. In combination !! the parameters x0 and y0, this specifies the position where the !! normalized toroidal magnetic field should equal 1.0. real(FP) :: xmin = 0.80_FP !! Lower x box limit real(FP) :: xmax = 1.20_FP !! Upper x box limit real(FP) :: ymin = -0.20_FP !! Lower y box limit real(FP) :: ymax = 0.20_FP !! Upper y box limit real(FP) :: rhomin = 0.0_FP !! Lower rho (normalized psi) limit (unused for Dommaschk) real(FP) :: rhomax = 1.1_FP !! Upper rho (normalized psi) limit (unused for Dommaschk) ! Specific Dommaschk parameters ! NOTE: The default specifies a simple rotated ellipse stellarator ! with island chain as seen in Fig. 2a of [*]. integer :: l_pol = DEFAULT_L_POL !! Poloidal periodicity upper bound integer :: m_tor_consecutive = DEFAULT_M_TOR !! Toroidal periodicity upper bound integer :: num_field_periods = 5 !! Number of field periods real(FP), dimension(:,:,:), allocatable, target :: fitting_coef_local !! Coefficients determining amplitudes of Dommaschk fourier modes ! NOTE: Since the shape of the fitting coefficient array is determined by ! the parameters l_pol and m_tor_consecutive, its initialization is ! handled differently. If the read routine is never called, and thus ! the fitting_coef_local is left unallocated, then the getter routine ! returns the default array below. Otherwise the allocated ! fitting_coef_local array is returned. real(FP), dimension(1:4, 0:DEFAULT_M_TOR, 0:DEFAULT_L_POL), target :: & fitting_coef_default = & reshape([0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, & 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, & 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, & 0.00_FP, 0.00_FP, 0.00_FP, 1.40_FP, 1.40_FP, 0.00_FP, & 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, & 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, 0.00_FP, & 0.00_FP, 19.25_FP, 0.00_FP, 0.00_FP], & shape=[4, DEFAULT_M_TOR + 1, DEFAULT_L_POL + 1]) !! Fitting coefficients for the default configuration for unit testing logical :: bndry_from_file = .true. !! Switch whether to use naive box limits or polygons from file for !! defining the boundaries of the equilibrium district character(len=str_buffer_len) :: inner_bndry_file = "" !! Path to NetCDF file which includes the desired inner boundary polygon. !! By default no inner boundary is used. character(len=str_buffer_len) :: inner_bndry_var = "" !! Variable name of inner polygon from the NetCDF file above character(len=str_buffer_len) :: outer_bndry_file = "default" !! Path to NetCDF file which includes the desired outer boundary polygon. !! The default is set for CI/CD. character(len=str_buffer_len) :: outer_bndry_var = "polygon_015" !! Variable name of outer polygon from the NetCDF file above logical :: exclusion_from_file = .true. !! Switch whether to use box limits or polygons from file for !! defining exclusion region character(len=str_buffer_len) :: exclusion_file = "default" !! Path to NetCDF file which includes the desired exclusion polygon. !! The default is set for CI/CD. character(len=str_buffer_len) :: exclusion_var = "exclusion_001" !! Variable name of exclusion polygon from the NetCDF file above ! NOTE: The boundary polygon NetCDF files should minimally have the ! following characteristics: ! - Dimensions: n_vertices, n_planes ! - Variable phi_array(n_planes) ! * Toroidal angles in the first field period of the given ! equilibrium. These should be evenly spaced and sorted in ! increasing order. ! - Variable for polygon (n_planes, 2, n_vertices) ! * Ordered polygon vertices in RZ space for every plane. Polygon ! should be closed (first point equal to last). logical :: rho_from_file = .false. !! Whether to calculate an approximate rho from flux surface polygons !! read from file character(len=str_buffer_len) :: rho_file = "" !! Path to NetCDF file which includes the desired flux surface polygons. ! NOTE: The flux surface polygon NetCDF file should minimally have the ! following characteristics ~in addition to~ those listed for the ! boundary polygon file above: ! - Dimensions: n_surf, n_surf_excluded ! - Variable surf_excluded(n_surf_excluded) ! * Poincare surface numbers which were not used to make polygons ! - Variables polygon_{000..n_surf}(n_planes, 2, n_vertices) ! * Ordered polygon vertices in RZ space for every plane. Polygon ! should be closed (first point equal to last). ! * The polygon numbers should be in I3.3 form, spanning 0 to ! n_surf, excluding the numbers in surf_excluded public :: read_params_dommaschk, write_params_dommaschk public :: get_dommaschk_x0, get_dommaschk_y0, get_dommaschk_phi0, & get_dommaschk_xmin, get_dommaschk_xmax, & get_dommaschk_ymin, get_dommaschk_ymax, & get_dommaschk_rhomin, get_dommaschk_rhomax, & get_dommaschk_l_pol, get_dommaschk_m_tor_consecutive, & get_dommaschk_num_field_periods, get_dommaschk_fitting_coef, & get_dommaschk_bndry_from_file, get_dommaschk_inner_bndry_file, & get_dommaschk_inner_bndry_var, get_dommaschk_outer_bndry_file, & get_dommaschk_outer_bndry_var, & get_dommaschk_exclusion_from_file, get_dommaschk_exclusion_file, & get_dommaschk_exclusion_var, get_dommaschk_rho_from_file, & get_dommaschk_rho_file namelist / params_equi_dommaschk / x0, y0, phi0, xmin, xmax, ymin, ymax, & rhomin, rhomax, & l_pol, m_tor_consecutive, & num_field_periods, & bndry_from_file, & inner_bndry_file, inner_bndry_var, & outer_bndry_file, outer_bndry_var, & exclusion_from_file, & exclusion_file, exclusion_var, & rho_from_file, rho_file real(FP), dimension(:,:,:), pointer :: fitting_coef !! Pointer to fitting coefficients, required for the parameter namelist namelist / params_equi_dommaschk_fitting_coef / fitting_coef contains pure real(FP) function get_dommaschk_x0() get_dommaschk_x0 = x0 end function pure real(FP) function get_dommaschk_y0() get_dommaschk_y0 = y0 end function pure real(FP) function get_dommaschk_phi0() get_dommaschk_phi0 = phi0 end function pure real(FP) function get_dommaschk_xmin() get_dommaschk_xmin = xmin end function pure real(FP) function get_dommaschk_xmax() get_dommaschk_xmax = xmax end function pure real(FP) function get_dommaschk_ymin() get_dommaschk_ymin = ymin end function pure real(FP) function get_dommaschk_ymax() get_dommaschk_ymax = ymax end function pure real(FP) function get_dommaschk_rhomin() get_dommaschk_rhomin = rhomin end function pure real(FP) function get_dommaschk_rhomax() get_dommaschk_rhomax = rhomax end function pure integer function get_dommaschk_l_pol() get_dommaschk_l_pol = l_pol end function pure integer function get_dommaschk_m_tor_consecutive() get_dommaschk_m_tor_consecutive = m_tor_consecutive end function pure integer function get_dommaschk_num_field_periods() get_dommaschk_num_field_periods = num_field_periods end function function get_dommaschk_fitting_coef() result(ptr) real(FP), dimension(:,:,:), pointer :: ptr if (allocated(fitting_coef_local)) then ptr => fitting_coef_local else ptr => fitting_coef_default end if end function pure logical function get_dommaschk_bndry_from_file() get_dommaschk_bndry_from_file = bndry_from_file end function pure function get_dommaschk_inner_bndry_file() result(res) character(len=:), allocatable :: res if (inner_bndry_file == "default") then res = git_repo_path // & "/res/dommaschk_equilibria/flux_polygons_dommaschk_default.nc" else res = trim(inner_bndry_file) endif end function pure function get_dommaschk_inner_bndry_var() result(res) character(len=:), allocatable :: res res = trim(inner_bndry_var) end function pure function get_dommaschk_outer_bndry_file() result(res) character(len=:), allocatable :: res if (outer_bndry_file == "default") then res = git_repo_path // & "/res/dommaschk_equilibria/flux_polygons_dommaschk_default.nc" else res = trim(outer_bndry_file) endif end function pure function get_dommaschk_outer_bndry_var() result(res) character(len=:), allocatable :: res res = trim(outer_bndry_var) end function pure logical function get_dommaschk_exclusion_from_file() get_dommaschk_exclusion_from_file = exclusion_from_file end function pure function get_dommaschk_exclusion_file() result(res) character(len=:), allocatable :: res if (exclusion_file == "default") then res = git_repo_path // & "/res/dommaschk_equilibria/exclusion_polygons_dommaschk_default.nc" else res = trim(exclusion_file) endif end function pure function get_dommaschk_exclusion_var() result(res) character(len=:), allocatable :: res res = trim(exclusion_var) end function pure logical function get_dommaschk_rho_from_file() get_dommaschk_rho_from_file = rho_from_file end function pure function get_dommaschk_rho_file() result(res) character(len=:), allocatable :: res if (rho_file == "default") then res = git_repo_path // & "/res/dommaschk_equilibria/flux_polygons_dommaschk_default.nc" else res = trim(rho_file) endif end function subroutine read_params_dommaschk(filename) character(len=*), intent(in) :: filename integer :: funit, io_error open(newunit=funit, file=filename, status='old', action='read', & iostat=io_error) if (io_error /= 0) then call handle_error("Opening parameter file "//filename//" failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if read(funit, nml=params_equi_dommaschk, iostat=io_error) ! END and EOR are returned if namelist can not be found. We treat ! these cases as non fatal and write a notification. if (io_error == IOSTAT_END .or. io_error == IOSTAT_EOR) then write(get_stdout(), *) & "Info: Namelist params_equi_dommaschk not present." elseif (io_error /= 0) then call handle_error("Reading Dommaschk equi parameters failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if if (allocated(fitting_coef_local)) deallocate(fitting_coef_local) allocate(fitting_coef_local(1:4, 0:m_tor_consecutive, 0:l_pol), & source=0.0_FP) fitting_coef => fitting_coef_local read(funit, nml=params_equi_dommaschk_fitting_coef, & iostat=io_error) ! END and EOR are returned if namelist can not be found. We treat ! these cases as non fatal and write a notification. if (io_error == IOSTAT_END .or. io_error == IOSTAT_EOR) then write(get_stdout(), *) & "Info: Namelist params_equi_dommaschk_fitting_coef not present." elseif (io_error /= 0) then call handle_error("Reading Dommaschk equi fitting coefficients & &failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if close(funit) ! If the periodicity of the equilibrium is the same as the default, ! and no non-zero fitting coefficients were provided, the default ! fitting coefficients are returned if(m_tor_consecutive == DEFAULT_M_TOR .and. & l_pol == DEFAULT_L_POL .and. & all(fitting_coef_local == 0.0_FP)) then fitting_coef_local = fitting_coef_default endif end subroutine subroutine write_params_dommaschk(filename) character(len=*), intent(in) :: filename integer :: funit, io_error open(newunit=funit, file=filename, action='write', position='append', & iostat=io_error) if (io_error /= 0) then call handle_error("Opening parameter file "//filename//" failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if write(funit, nml=params_equi_dommaschk, iostat=io_error) if (io_error /= 0) then call handle_error("Writing Dommaschk equi parameters failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if fitting_coef => get_dommaschk_fitting_coef() write(funit, nml=params_equi_dommaschk_fitting_coef, & iostat=io_error) if (io_error /= 0) then call handle_error("Writing Dommaschk equi fitting coefficients & &failed!", & PARALLAX_ERR_PARAMETERS, __LINE__, __FILE__, & additional_info=& error_info_t("Iostat:", [io_error])) end if close(funit) end subroutine end module