module iol_source_m !! Compute ion orbit loss source/sink coupling to ion-orbit-loss-standalone !! see https://gitlab.mpcdf.mpg.de/phoenix/ion-orbit-loss-standalone !! [R.W. Brzozowski III, Phys. Plasmas 26:042511 (2019), !! https://doi.org/10.1063/1.5075613] ! Note: ! This module does not completely follow the code and naming conventions, ! but rather uses names and conventions of the ion-orbit-loss-standalone. ! For usage of the library and the detailed meaning of variables, we refer ! to that library use MPI use netcdf use constants_m, only : two_pi use precision_grillix_m, only : GP, GP_NAN, GP_LARGEST, MPI_FP use screen_io_m, only : get_stdout use error_handling_grillix_m, only: handle_error_netcdf, handle_error, & error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER use physical_constants_m, only : mass_hydrogen_SI, elementary_charge_SI, & speed_of_light_SI use comm_handler_m, only : comm_handler_t use polars_m, only : polars_t use variable_m, only: variable_t use perf_m, only : perf_start, perf_stop use equilibrium_m, only : equilibrium_t use numerical_equilibrium_m, only : numerical_t use mesh_cart_m, only : mesh_cart_t use coords_polar_m, only : polar_to_cart, cart_to_polar use polar_operators_m, only : map_cart_to_polar, map_polar_to_cart use safety_factor_m, only : safety_factor ! From IOL use iol_m, only : iol_eval implicit none real(GP), parameter, private :: phi_axsym = 0.0_GP !! This is a placeholder for calling equilibria routines with 3D interfaces !! IOL module only works for axisymmetric configurations private :: flip_array_direction type, public :: iol_source_t !! Contains information/routines to compute ion-orbit-loss effects !! See interface description of ion-orbit-loss standalone routines integer, private:: imax !! Number of radial grid points integer, private :: jmax !! Number of poloidal grid points integer, private :: signBt !! Sign of toroidal field following AUG conventions integer, private :: signIp !! Sign of toroidal current following AUG conventions real(GP), private :: RX !! R-coordinate of X-point [m] real(GP), private :: zX !! Z-coordinate of X-point [m] real(GP), private :: BphiX !! Signed value of toroidal field at X-point [T], & !! following AUG convention real(GP), private :: psirange_conf_red !! Difference of poloidal flux between magnetic axis and separatrix & !! [Weber=T*m^2] integer, private :: idx_imp !! Poloidal index of inner midplane (HFS) integer, private :: idx_X !! Poloidal index of X-point real(GP), private :: invs_aspc !! Approximate inverse aspect ratio (a/R_0) real(GP), private :: R_magaxis !! Radial position of magnetiox axis [m] real(GP), private :: Rmin_m !! Slope of linear fit for Rmin(psi) w/ Rmin in [m] and psi !! in [Wb/2pi] [2pi m/Wb] !! TODO: WHAT ARE THE UNITS HERE, HOWW TO COMPUTE THIS ??? @Robert ???? real(GP), private :: Rmin_b !! Intercept for linear fit for Rmin(psi) w/ Rmin n [m] !! and psi in [Wb/2pi] [m] !! TODO: WHAT ARE THE UNITS HERE, HOWW TO COMPUTE THIS ??? @Robert ???? real(GP), private :: mass !! Ion mass [kg] real(GP), private :: charge !! Ion charge [C] integer, private :: kmax !! Number of points in pitch angle space integer, private :: lmax !! Number of points in energy space integer, private :: mmax !! Number of points to consider as IOL for R_loss > R_X real(GP), private :: er_mult !! Scaling factor for radial electric field real(GP), private :: ti_mult !! Scaling factor for temperature real(GP), private, allocatable, dimension(:,:) :: Rcoord !! R on grid [m] real(GP), private, allocatable, dimension(:,:) :: zcoord !! Z on grid [m] real(GP), private, allocatable, dimension(:,:) :: Babs !! Absolute value of magnetic field on grid [T] real(GP), private, allocatable, dimension(:,:) :: Bphi !! Signed toroidal magnetic field on grid [T], following AUG convention real(GP), private, allocatable, dimension(:,:) :: psi !! Polodial flux on grid [Wb =T*m^2] real(GP), private, allocatable, dimension(:) :: sep_Rcoord !! Radial positions along separatrix [m] real(GP), private, allocatable, dimension(:) :: sep_cw !! Clockwise poloidal distance from point (j,imax) to X-point [m] real(GP), private, allocatable, dimension(:) :: sep_ccw !! Counterclockwise poloidal distance from point (j,imax) to X-point [m] real(GP), private, allocatable, dimension(:) :: sep_cw_mod !! Corrected (for poloidal stagnation) clockwise poloidal distance !! from point (j,i) to X-point [m] real(GP), private, allocatable, dimension(:) :: sep_ccw_mod !! Corrected (for poloidal stagnation) counterclockwise !! poloidal distance from point (j,i) to integer, private :: Xplane_n !! Number of points for low-field side plane at X-point real(GP), private, allocatable, dimension(:) :: xplane_R !! R coordinate for low-field side plane at X-point [m] real(GP), private, allocatable, dimension(:) :: xplane_psi !! Poloidal flux for low-field side plane at X-point [Weber = T*m^2] real(GP), private, allocatable, dimension(:) :: xplane_Babs !! Absolute value for low-field side plane at X-point [T] real(GP), private, allocatable, dimension(:) :: xplane_Bphi !! Signed toroidal magnetic field component !! for low-field side plane at X-point [T] real(GP), private, allocatable, dimension(:) :: safetyf !! Safety factor integer, private :: nrho_sp !! Number of collocation points for spline interpolation !! (for radial electric field, suggested >~ imax) real(GP), private, allocatable, dimension(:) :: rhopts_sp !! Radial grid positions rho_pol, normalised from [0..1] real(GP), private, allocatable, dimension(:) :: Rmincontour_sp !! Position of minimum R on flux surface [m] ! Further auxiliary variables ---------------------------------- real(GP), private :: psix !! Poloidal flux at separatrix [T*m^2] integer, private, allocatable, dimension(:) :: indmesh_sp !! Indices on mesh corresponding to spline interpolation grid real(GP), private, allocatable, dimension(:) :: sep_Zcoord !! Vertical positions along separatrix [m] real(GP), private :: bnrm !! Normalisation of magnetic field (usally B on axis) [T] real(GP), private :: ne_nrm !! Normalisation factor for density [m^{-3}] real(GP), private :: ti_nrm !! Normalisation factor for ion temperature !! and electrostatic potential [eV] real(GP), private :: tau_nrm !! Time normalisation R_0/c_s [s] integer, private :: nsep_corr_cw !! Number of points that are corrected around the X-point !! in clockwise direction !! at finding separatrix integer, private :: nsep_corr_ccw !! Number of points that are corrected around the X-point !! in counter clockwise direction !! at finding separatrix ! Variables for time averaging integer, private :: nt_av = 0 !! Number of timepoints that the fields are currently averaged over real(GP), private :: tau_av_frame !! Time that fields are averaged over logical, private :: init_continue !! Switch if simulation continued from existing state, i.e. !! IOL source initially read from existing file ! Internally stored fields, averaged over time real(GP), private, allocatable, dimension(:,:) :: dens !! Ion density on iol grid [m^-3] real(GP), private, allocatable, dimension(:,:) :: temp_ion !! Ion temperature on iol grid [J] real(GP), private, allocatable, dimension(:,:) :: pot !! Electrostatic potential on iol grid [V] !! (naming convention here different from IOL-standalone !! to avoid confusion with toroidal angle) real(GP), private, allocatable, dimension(:) :: pot_sp !! Electrostatic potential on spline grid [V] !! (naming convention here different from IOL-standalone !! to avoid confusion with toroidal angle) real(GP), private, allocatable, dimension(:) :: src_iol !! Ion Orbit loss source/sink on inner mesh points of GRILLIX grid contains ! Initialisation routines procedure, public :: init => init_iol procedure, private :: read_params_iol procedure, private :: grid_init procedure, private :: spline_grid_init procedure, private :: xplane_init procedure, private :: flip ! Getter procedure, public :: get_init_continue ! Source comoputation routines procedure, public :: driver => driver_iol procedure, private :: add_fields_to_averaged_internal_fields procedure, private :: update_iol_source ! I/O routines procedure, public :: display => display_iol_source procedure, public :: write_netcdf => write_netcdf_iol_source procedure, public :: read_netcdf => read_netcdf_iol_source procedure, public :: write_ascii => write_ascii_iol_source ! Destructor final :: destructor_iol end type interface module subroutine init_iol(self, filename, equi, mesh, polars) !! Initialises IOL source character(len=*), intent(in) :: filename !! Filename, to read parameters from class(iol_source_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh type(polars_t), intent(in) :: polars !! Polar grid and operators end subroutine module subroutine read_params_iol(self, filename, equi) !! Reads IOL parameters character(len=*), intent(in) :: filename !! Filename, to read parameters from class(iol_source_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium end subroutine module subroutine grid_init(self, equi, polars) !! Initialises IOL grid ( = polar grid + extended information) class(iol_source_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(polars_t), intent(in) :: polars !! Polar grid and operators end subroutine module subroutine spline_grid_init(self, equi, mesh, polars) !! Initialises spline grid class(iol_source_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh type(polars_t), intent(in) :: polars !! Polar grid and operators end subroutine module subroutine xplane_init(self, equi, mesh) !! Initialises xplane, i.e. horizontal plane intersecting X-point class(iol_source_t), intent(inout) :: self !! Instance of the type class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane end subroutine module subroutine flip(self) !! Flips grids vertically to account for flipped equilibria class(iol_source_t), intent(inout) :: self !! Instance of the type end subroutine module subroutine driver_iol(self, comm_handler, equi, mesh, polars, & tau, dtau, ne, pot, ti, src_vort, update_performed, force_no_update) !! Driver routine for computing IOL sources. !! Takes care of time averaging procedure, i.e. !! (Re-)computes IOL sources at time frames tau_average !! based on internal average fields. !! Adds contribution of fields to current internal average fields class(iol_source_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Communication handler class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(polars_t), intent(in) :: polars !! Polar grid and operators real(GP), intent(in) :: tau !! Current time of simulation real(GP), intent(in) :: dtau !! Timestep of simulation type(variable_t), intent(in) :: ne !! Electron density type(variable_t), intent(in) :: pot !! Electrostatic potential type(variable_t), intent(in) :: ti !! Ion temperature real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: & src_vort !! Vorticity source, on output with IOL contribution added logical, intent(out) :: update_performed !! Indicates if update of IOL source was performed logical, intent(in), optional :: force_no_update !! Switch that enforces that IOL source is not updated end subroutine module subroutine add_fields_to_averaged_internal_fields(self, comm_handler, & equi, mesh, polars, ne, pot, ti) !! Adds current fields ne, pot and ti, to current internal !! time averaged fields used for IOL source evaluation class(iol_source_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Communication handler class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(polars_t), intent(in) :: polars !! Polar grid and operators type(variable_t), intent(in) :: ne !! Electron density type(variable_t), intent(in) :: pot !! Electrostatic potential type(variable_t), intent(in) :: ti !! Ion temperature end subroutine module subroutine update_iol_source(self, comm_handler, equi, mesh, polars) !! Internally updates the iol source with the current averaged fields !! Here is the interface to the IOL-standalone library !! see https://gitlab.mpcdf.mpg.de/phoenix/ion-orbit-loss-standalone !! for its interface description class(iol_source_t), intent(inout) :: self !! Instance of the type type(comm_handler_t), intent(in) :: comm_handler !! Communication handler class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(polars_t), intent(in) :: polars !! Polar grid and operators end subroutine module subroutine display_iol_source(self) !! Displays parameters for iol_source class(iol_source_t), intent(inout) :: self !! Instance of the type end subroutine module subroutine write_netcdf_iol_source(self, tau, fexist, fgid) !! Writes current IOL source to NETCDF file class(iol_source_t), intent(in) :: self !! Instance of the type real(GP), intent(in) :: tau !! Timestamp of current IOL source logical, intent(in) :: fexist !! Indicates if file integer, intent(in) :: fgid !! File or group id number of existing NETCDF file end subroutine module subroutine read_netcdf_iol_source(self, fgid, nsnaps_req, nsnaps, tau) !! Reads IOL source from NETCDF file class(iol_source_t), intent(inout) :: self !! Instance of the type integer, intent(in) :: fgid !! File or group id number of existing NETCDF file integer, intent(in), optional :: nsnaps_req !! Number of requested snapshot integer, intent(out), optional :: nsnaps !! Number of snapshot returned real(GP), intent(out), optional :: tau !! Tme stamp of snapshot end subroutine module subroutine write_ascii_iol_source(self, dirpath) !! Writes information of IOL to ASCII format !! Useful for running test case with IOL-standalone !! see https://gitlab.mpcdf.mpg.de/phoenix/ion-orbit-loss-standalone class(iol_source_t), intent(inout) :: self !! Instance of the type character(len=*), intent(in) :: dirpath !! Directory, where files are written end subroutine end interface contains pure logical function get_init_continue(self) !! Getter routine class(iol_source_t), intent(in) :: self !! Instance of the type get_init_continue = self%init_continue end function subroutine destructor_iol(self) !! Deallocates stuff related with iol source type(iol_source_t), intent(inout) :: self !! Instance of the type if (allocated(self%Rcoord)) deallocate(self%Rcoord) if (allocated(self%zcoord)) deallocate(self%zcoord) if (allocated(self%Babs)) deallocate(self%Babs) if (allocated(self%Bphi)) deallocate(self%Bphi) if (allocated(self%psi)) deallocate(self%psi) if (allocated(self%sep_Rcoord)) deallocate(self%sep_Rcoord) if (allocated(self%sep_cw)) deallocate(self%sep_cw) if (allocated(self%sep_ccw)) deallocate(self%sep_ccw) if (allocated(self%sep_cw_mod)) deallocate(self%sep_cw_mod) if (allocated(self%sep_ccw_mod)) deallocate(self%sep_ccw_mod) if (allocated(self%xplane_R)) deallocate(self%xplane_R) if (allocated(self%xplane_psi)) deallocate(self%xplane_psi) if (allocated(self%xplane_Babs)) deallocate(self%xplane_Babs) if (allocated(self%xplane_Bphi)) deallocate(self%xplane_Bphi) if (allocated(self%safetyf)) deallocate(self%safetyf) if (allocated(self%rhopts_sp)) deallocate(self%rhopts_sp) if (allocated(self%Rmincontour_sp)) deallocate(self%Rmincontour_sp) if (allocated(self%indmesh_sp)) deallocate(self%indmesh_sp) if (allocated(self%sep_Zcoord)) deallocate(self%sep_Zcoord) if (allocated(self%dens)) deallocate(self%dens) if (allocated(self%temp_ion)) deallocate(self%temp_ion) if (allocated(self%pot)) deallocate(self%pot) if (allocated(self%pot_sp)) deallocate(self%pot_sp) if (allocated(self%src_iol)) deallocate(self%src_iol) end subroutine subroutine flip_array_direction(nd, arr) !! Flips direction of array entries !! Needed to handle flipped equilibria integer, intent(in) :: nd !! Dimension of array real(GP), dimension(nd), intent(inout) :: arr !! Array, on output revertede real(GP), dimension(nd) :: wrk integer :: i !$OMP PARALLEL PRIVATE(i) !$OMP DO do i = 1, nd wrk(i) = arr(i) enddo !$OMP END DO !$OMP DO do i=2, nd arr(i) = wrk(nd-i+2) enddo !$OMP END DO !$OMP END PARALLEL end subroutine end module