module init_template_m !! Creates an initial state use precision_grillix_m, only : GP use parallelisation_setup_m, only : is_rank_info_writer use screen_io_m, only : get_stdout use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER use elementary_functions_m, only : step_hermite, gaussian use comm_handler_m, only : comm_handler_t use static_data_m, only : mesh_cano, equi_on_cano, masks_cano, mesh_stag, masks_stag, map use state_template_m, only : np_max, dens, parmom, pion, tau, tstep, isnaps, read_template_state_netcdf, & update_derived_fields use params_template_m, only : path_snapsdir, path_init_select use mms_template_m, only: mms_sol_dens_template, mms_sol_parmom_template, mms_sol_pion_template implicit none public :: init_template private :: init_template_profile private :: init_template_blob contains subroutine init_template(comm_handler) !! Creates an initial state depending on init_type type(comm_handler_t), intent(in) :: comm_handler !! Communication handler real(GP) :: x, y, phi_cano, phi_stag integer :: l, io_unit, io_error character(len=256) :: io_errmsg character(len=32) :: init_type = 'CONTINUE' !! Type of initial state, possible values: !! - 'CONTINUE' : Continues from last available snapshot file (default) !! - 'PROFILE': Creates initial state based on density and pressure profiles namelist / init_template_type / init_type if (is_rank_info_writer) then write(get_stdout(),*)new_line('a') // repeat('-',80) write(get_stdout(),*)'Creating initial state' endif ! Read init_type open(newunit=io_unit, file=path_init_select, status='old', action='read', iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) read(io_unit, nml=init_template_type, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) close(io_unit) ! Display init_type if(is_rank_info_writer) write(get_stdout(),nml=init_template_type) select case(init_type) case('CONTINUE') call read_template_state_netcdf(comm_handler, dirname=trim(path_snapsdir)) if (is_rank_info_writer) then write(get_stdout(),*)'Continuing from snapshot file' write(get_stdout(),501)isnaps, tau, tstep 501 FORMAT(3X,'isnaps = ',I5,'; tau = 'ES20.12E3,'; tstep = ',I10) endif case('PROFILE') call init_template_profile() tau = 0.0_GP tstep = 0 if (is_rank_info_writer) then write(get_stdout(),*)'Using PROFILE as initial state' endif case('MMS') ! Set initial state to MMS solution tau = 0.0_GP tstep = 0 phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() !$omp parallel default(none) & !$omp private(l, x, y) & !$omp shared(np_max, mesh_cano, masks_cano, phi_cano, mesh_stag, masks_stag, phi_stag, & !$omp tau, dens, mms_sol_dens_template, parmom, mms_sol_parmom_template, pion, mms_sol_pion_template) !$omp do do l = 1, np_max if (masks_cano%is_inner(l) .or. masks_cano%is_ghost(l) .or. masks_cano%is_boundary(l)) then x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) dens(l) = mms_sol_dens_template(x, y, phi_cano, tau) pion(l) = mms_sol_pion_template(x, y, phi_cano, tau) endif if (masks_stag%is_inner(l) .or. masks_stag%is_ghost(l) .or. masks_stag%is_boundary(l)) then x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) parmom(l) = mms_sol_parmom_template(x, y, phi_stag, tau) endif enddo !$omp end do !$omp end parallel if (is_rank_info_writer) then write(get_stdout(),*)'Using MMS as initial state' endif case('BLOB') call init_template_blob() tau = 0.0_GP tstep = 0 if (is_rank_info_writer) then write(get_stdout(),*)'Using BLOB as initial state' endif case default call handle_error('init_type not valid:', GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(init_type)) end select ! Compute derived fields call update_derived_fields(comm_handler) if (is_rank_info_writer) then write(get_stdout(),*)repeat('-',80) // new_line('a') endif end subroutine subroutine init_template_profile() !! Builds initial state, based on profile for density and pressures integer :: io_unit, io_error character(len=256) :: io_errmsg character(len=32) :: profile_type = 'STEP_HERMITE' ! Parameters for step_hermite init state real(GP) :: dens_rho_top = 0.9_GP real(GP) :: dens_rho_btm = 1.0_GP integer :: l real(GP) :: rho, rho0_dens, wrho_dens namelist / init_template_profile_type / profile_type namelist / init_template_step_hermite / dens_rho_top, dens_rho_btm ! Read profile_type open(newunit=io_unit, file=path_init_select, status='old', action='read', iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) read(io_unit, nml=init_template_profile_type, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) close(io_unit) ! Display profile type if(is_rank_info_writer) write(get_stdout(),nml=init_template_profile_type) ! Init density and pressure profiles select case(profile_type) case('STEP_HERMITE') ! Sigmoid function based on step hermite open(newunit=io_unit, file=path_init_select, status='old', action='read', iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) read(io_unit, nml=init_template_step_hermite, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) close(io_unit) if(is_rank_info_writer) write(get_stdout(),nml=init_template_step_hermite) rho0_dens = 0.5_GP * (dens_rho_btm + dens_rho_top) wrho_dens = 2.0_GP * (dens_rho_btm - dens_rho_top) !$omp parallel default(none) private(l, rho) & !$omp shared(np_max, equi_on_cano, masks_cano, dens, rho0_dens, wrho_dens) !$omp do do l = 1, np_max if (masks_cano%is_inner(l) .or. masks_cano%is_boundary(l) .or. masks_cano%is_ghost(l)) then rho = equi_on_cano%rho(l) dens(l) = 1.0_GP - step_hermite(rho0_dens, rho, wrho_dens) endif enddo !$omp enddo !$omp end parallel case default call handle_error('profile_type not valid:', GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(profile_type)) end select end subroutine subroutine init_template_blob() !! Builds initial state, based on blob for density integer :: io_unit, io_error character(len=256) :: io_errmsg integer :: l real(GP) :: x, y, phi ! Parameters for blob init state real(GP) :: xc = 1.18_GP real(GP) :: yc = 0.0_GP real(GP) :: phic = 0.0_GP real(GP) :: wx = 0.075_GP real(GP) :: wy = 0.075_GP real(GP) :: wphi = 1.0E-10_GP real(GP) :: amp_bck = 1.0_GP real(GP) :: amp_blob = 1.0_GP namelist / init_blob / xc, yc, wx, wy, phic, wphi, amp_bck, amp_blob phi = mesh_cano%get_phi() ! Read profile_type open(newunit=io_unit, file=path_init_select, status='old', action='read', iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) read(io_unit, nml=init_blob, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) close(io_unit) ! Display blob parameters if(is_rank_info_writer) write(get_stdout(),nml=init_blob) !$omp parallel default(none) private(l, x, y) & !$omp shared(np_max, mesh_cano, masks_cano, dens, phi, xc, yc, wx, wy, phic, wphi, amp_bck, amp_blob) !$omp do do l = 1, np_max if (masks_cano%is_inner(l) .or. masks_cano%is_boundary(l) .or. masks_cano%is_ghost(l)) then x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) dens(l) = amp_bck + amp_blob * gaussian(xc, wx, x) * gaussian(yc, wy, y) * gaussian(phic, wphi, phi) endif enddo !$omp enddo !$omp end parallel end subroutine end module