module create_buffer_m !! Contains routines to create buffer function use netcdf use precision_grillix_m, only : GP 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 parallelisation_setup_m, only : is_rank_info_writer ! from PARALLAX use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use elementary_functions_m, only : step_hermite use descriptors_m, only : DISTRICT_CORE, DISTRICT_CLOSED, DISTRICT_SOL, DISTRICT_WALL, & DISTRICT_PRIVFLUX, DISTRICT_DOME, DISTRICT_OUT implicit none public :: create_buffer public :: buffer_write_netcdf private :: create_buffer_zeronull private :: create_buffer_singlenull private :: create_buffer_sample_mms private :: create_buffer_none contains subroutine create_buffer(equi, mesh, buffer_select, buffer_paramfile, cf_buffer) !! Creates buffer function class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane character(len = *), intent(in) :: buffer_select !! Type of buffer function character(len = *), intent(in) :: buffer_paramfile !! Filename where to read buffer function parameters from real(GP), dimension(mesh%get_n_points()), intent(out) :: cf_buffer !! Buffer function if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'Creating buffer function' endif select case(buffer_select) case('ZERONULL') if (is_rank_info_writer) then write(get_stdout(),*)'ZERONULL' endif call create_buffer_zeronull(equi, mesh, buffer_paramfile, cf_buffer) case('SINGLENULL') if (is_rank_info_writer) then write(get_stdout(),*)'SINGLENULL' endif call create_buffer_singlenull(equi, mesh, buffer_paramfile, cf_buffer) case('SAMPLE_MMS') if (is_rank_info_writer) then write(get_stdout(),*)'SAMPLE_MMS' endif call create_buffer_sample_mms(equi, mesh, cf_buffer) case default if (is_rank_info_writer) then write(get_stdout(),*)'none' endif ! cf_buffer = 0 call create_buffer_none(equi, mesh, cf_buffer) end select if (is_rank_info_writer) then write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'' endif end subroutine subroutine create_buffer_zeronull(equi, mesh, filename, cf_buffer) !! Creates buffer function for zeronull equilibria, e.g. circular without private flux region !! Step function based on hermite functions class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane character(len = *), intent(in) :: filename !! Filename where to read buffer function parameters from real(GP), dimension(mesh%get_n_points()), intent(out) :: cf_buffer !! Buffer function real(GP) :: width_core !! Width of core buffer region in units of (rhomax-rhomin) real(GP) :: width_edge !! Width of edge buffer region real(GP) :: step_width_core !! Decay length of core buffer region real(GP) :: step_width_edge !! Decay length of edge buffer region integer :: io_error character(len=256) :: io_errmsg integer :: l real(GP) :: phi, x, y, rhon, rhowidth namelist / buffer_zeronull / width_core, width_edge, step_width_core, step_width_edge ! Read input from file ------------------------------------------------ open(unit = 25, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif ! default values width_core = 0.0_GP width_edge = 0.0_GP step_width_core = 0.0_GP step_width_edge = 0.0_GP read(25, nml = buffer_zeronull, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif rhowidth = equi%rhomax - equi%rhomin phi = mesh%get_phi() !$OMP PARALLEL PRIVATE(l, x, y, rhon) !$OMP DO do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rhon = (equi%rho(x, y, phi) - equi%rhomin) / rhowidth cf_buffer(l) = step_hermite(1.0_GP-width_edge, rhon, step_width_edge) cf_buffer(l) = cf_buffer(l) + 1.0_GP - step_hermite(width_core, rhon, step_width_core) enddo !$OMP END DO !$OMP END PARALLEL close(25) end subroutine subroutine create_buffer_singlenull(equi, mesh, filename, cf_buffer) !! Creates buffer function for single-null equilibria, i.e. with a single private flux region !! Step function based on hermite functions class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane character(len = *), intent(in) :: filename !! Filename where to read buffer function parameters from real(GP), dimension(mesh%get_n_points()), intent(out) :: cf_buffer !! Buffer function real(GP) :: rho_buff_private !! Private flux rho, where buffer regions starts real(GP) :: width_core !! Width of core buffer region in units of (rhomax-rhomin) real(GP) :: width_sol !! Width of sol buffer region in units of (rhomax-rhomin) real(GP) :: width_private !! Width of private flux buffer region in units of (1-rhomin_privflux) real(GP) :: step_width_core !! Decay length of core buffer region real(GP) :: step_width_sol !! Decay length of edge buffer region real(GP) :: step_width_private !! Decay length of private flux buffer region in units of (1-rhomin_privflux) integer :: io_error character(len=256) :: io_errmsg integer :: l, district real(GP) :: phi, x, y, rhow_can, rhow_pf, rho, rhon namelist / buffer_singlenull / rho_buff_private, width_core, width_sol, width_private, & step_width_core, step_width_sol, step_width_private ! Read input from file ------------------------------------------------ open(unit = 25, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif ! default values rho_buff_private = 0.0_GP width_core = 0.0_GP width_sol = 0.0_GP width_private = 0.0_GP step_width_core = 0.0_GP step_width_sol = 0.0_GP step_width_private = 0.0_GP read(25, nml = buffer_singlenull, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif rhow_can = equi%rhomax - equi%rhomin rhow_pf = 1.0_GP - rho_buff_private phi = mesh%get_phi() !$OMP PARALLEL PRIVATE(l, x, y, rho, district, rhon) !$OMP DO do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) district = mesh%get_district(l) select case (district) case (DISTRICT_CORE, DISTRICT_CLOSED, DISTRICT_SOL, DISTRICT_WALL) rhon = (rho - equi%rhomin) / rhow_can cf_buffer(l) = step_hermite(1.0_GP-width_sol, rhon, step_width_sol) cf_buffer(l) = cf_buffer(l) + 1.0_GP - step_hermite(width_core, rhon, step_width_core) case (DISTRICT_PRIVFLUX, DISTRICT_DOME) rhon = (rho - rho_buff_private) / rhow_pf cf_buffer(l) = 1.0_GP - step_hermite(width_private, rhon, step_width_private) case (DISTRICT_OUT) cf_buffer(l) = 0.0_GP case default !$omp critical call handle_error('District not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('District: ',[district])) !$omp end critical end select enddo !$OMP END DO !$OMP END PARALLEL close(25) end subroutine subroutine create_buffer_sample_mms(equi, mesh, cf_buffer) !! Returns a smooth and analytic buffer function used for MMS verification class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane real(GP), dimension(mesh%get_n_points()), intent(out) :: cf_buffer !! Buffer function integer :: l real(GP) :: rhocenter, rhowidth, phi, x, y, rho rhocenter = (equi%rhomin+equi%rhomax)*0.5_GP rhowidth = equi%rhomax - equi%rhomin phi = mesh%get_phi() !$OMP PARALLEL PRIVATE(l, x, y, rho) !$OMP DO do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) cf_buffer(l) = ( rho - rhocenter )**2 / (0.5_GP * rhowidth)**2 enddo !$OMP END DO !$OMP END PARALLEL end subroutine subroutine create_buffer_none(equi, mesh, cf_buffer) !! Returns buffer function with zeros class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane real(GP), dimension(mesh%get_n_points()), intent(out) :: cf_buffer !! Buffer function integer l !$OMP PARALLEL !$OMP DO do l = 1, mesh%get_n_points() cf_buffer(l) = 0.0_GP enddo !$OMP END DO !$OMP END PARALLEL end subroutine subroutine buffer_write_netcdf(mesh, cf_buffer, fgid) !! Writes cf_buffer to netcdf file type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane real(GP), dimension(mesh%get_n_points()), intent(in) :: cf_buffer !! Buffer function integer, intent(in) :: fgid !! Group or file id integer :: nf90_stat integer :: n_points integer :: iddim_n_points, id_cf_buffer n_points = mesh%get_n_points() ! define dimensions nf90_stat = nf90_def_dim(fgid, 'n_points', n_points, iddim_n_points) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! define variables nf90_stat = nf90_def_var(fgid, 'cf_buffer', NF90_DOUBLE, iddim_n_points, id_cf_buffer ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! end of definition nf90_stat = nf90_enddef(fgid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) ! put variables nf90_stat = nf90_put_var(fgid, id_cf_buffer, cf_buffer ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) end subroutine end module