module init_neutrals_m !! Sets initial state for neutrals module use MPI use netcdf use comm_handler_m, only : comm_handler_t use parallelisation_setup_m, only : is_rank_info_writer use netcdf use precision_grillix_m, only : GP, GP_EPS use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER use parallel_map_m, only : parallel_map_t use variable_m, only: variable_t use multistep_m, only : karniadakis_t use snapshot_m, only : snapshot_t use mms_neutrals_m, only : mms_sol_neutrals_dens, mms_sol_neutrals_parmom, & mms_sol_neutrals_pressure ! From PARALLAX use constants_m, only : PI, TWO_PI use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t ! Parameters use params_neut_init_select_m, only : & neut_iselect implicit none public :: build_neutrals_initial_state public :: initialise_neutrals_timesteppers public :: init_neutrals_halos contains subroutine build_neutrals_initial_state(comm_handler, equi, & mesh_cano, mesh_stag, & filename, & snaps_neutrals, & isnaps_neutrals, & ti, & neutrals_dens, & neutrals_parmom, & neutrals_pressure) !! Builds initial state type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(inout) :: mesh_stag !! Mesh (staggered) character(len=*), intent(in) :: filename !! Filename, where parameters are read from type(snapshot_t), intent(in) :: snaps_neutrals !! Snapshot, where initial state is possibly read from integer, intent(out) :: isnaps_neutrals !! Snapshot of neutral model type(variable_t), intent(in) :: ti !! Ion temperature to use as initial neutrals temperature, !! if no valid initial state parameters for neutral pressure are provided type(variable_t), intent(inout) :: neutrals_dens !! Neutrals density type(variable_t), intent(inout) :: neutrals_parmom !! Neutrals parallel momentum type(variable_t), intent(inout) :: neutrals_pressure !! Neutrals pressure namelist / init_homogeneous / & neutrals_dens_init_val, & neutrals_parmom_init_val, & neutrals_pressure_init_val integer :: io_error character(len=256) :: io_errmsg integer :: l real(GP) :: tau_neutrals, x, y, phi_cano, phi_stag real(GP) :: neutrals_dens_init_val, neutrals_parmom_init_val, neutrals_pressure_init_val type(variable_t), allocatable, dimension(:) :: tmp_var select case(neut_iselect) case('MMS') ! Set initial state to MMS solution ----------------------------------------------- isnaps_neutrals = 1 tau_neutrals = 0.0_GP phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() call neutrals_dens%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call neutrals_parmom%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call neutrals_pressure%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) !$omp parallel default(none) private(l, x, y) & !$omp shared(equi, mesh_cano, mesh_stag, & !$omp neutrals_dens, neutrals_parmom, neutrals_pressure, & !$omp phi_cano, phi_stag, tau_neutrals) !$omp do do l = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) neutrals_dens%vals(l) = mms_sol_neutrals_dens(equi, x, y, phi_cano, tau_neutrals) neutrals_pressure%vals(l) = mms_sol_neutrals_pressure(equi, x, y, phi_cano, tau_neutrals) enddo !$omp end do !$omp do do l = 1, mesh_stag%get_n_points() x = mesh_stag%get_x(l) y = mesh_stag%get_y(l) neutrals_parmom%vals(l) = mms_sol_neutrals_parmom(equi, x, y, phi_stag, tau_neutrals) enddo !$omp end do !$omp end parallel case('HOMOGENEOUS') ! Read parameters of homogeneous initial state ------------------------------------ isnaps_neutrals = 1 open(unit = 20, 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 neutrals_dens_init_val = 0.0_GP neutrals_parmom_init_val = 0.0_GP neutrals_pressure_init_val = 0.0_GP read(20, nml = init_homogeneous, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif if (is_rank_info_writer) then write(get_stdout(),*)'Creating homegenous as initial state for neutrals -----' write(get_stdout(),202) neutrals_dens_init_val, & neutrals_parmom_init_val, neutrals_pressure_init_val 202 FORMAT(1X,'neutrals_dens_init_val = ',ES14.6E3 /, & 1X,'neutrals_parmom_init_val = ',ES14.6E3 /, & 1X,'neutrals_pressure_init_val = ',ES14.6E3 ) write(get_stdout(),*)'-------------------------------------------------------' endif close(20) call neutrals_dens%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call neutrals_parmom%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call neutrals_pressure%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) !$omp parallel default(none) private(l) & !$omp shared(mesh_cano, mesh_stag, ti, neutrals_dens, neutrals_parmom, neutrals_pressure, & !$omp neutrals_dens_init_val, neutrals_parmom_init_val, neutrals_pressure_init_val) !$omp do do l = 1, mesh_cano%get_n_points() neutrals_dens%vals(l) = neutrals_dens_init_val enddo !$omp end do !$omp do do l = 1, mesh_stag%get_n_points() neutrals_parmom%vals(l) = neutrals_parmom_init_val enddo !$omp end do if (neutrals_pressure_init_val > GP_EPS) then !$omp do do l = 1, mesh_cano%get_n_points() neutrals_pressure%vals(l) = neutrals_pressure_init_val enddo !$omp end do else !$omp do do l = 1, mesh_cano%get_n_points() neutrals_pressure%vals(l) = ti%vals(l) * neutrals_dens_init_val enddo !$omp end do endif !$omp end parallel case('CONTINUE') allocate(tmp_var(3)) call snaps_neutrals%read_snaps_last(comm_handler, isnaps_neutrals, tau_neutrals, tmp_var) neutrals_dens = tmp_var(1) neutrals_parmom = tmp_var(2) neutrals_pressure = tmp_var(3) deallocate(tmp_var) if (is_rank_info_writer) then write(get_stdout(),*)'Continuing simulation (neutrals) at--------------------' write(get_stdout(),201)isnaps_neutrals, tau_neutrals 201 FORMAT(1X,'isnaps_neutrals = ',I8 /, & 1X,'tau_neutrals = ',ES14.6E3) write(get_stdout(),*)'-------------------------------------------------------' endif case default call handle_error('Init_select not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(neut_iselect)) end select end subroutine subroutine initialise_neutrals_timesteppers(mesh_cano, mesh_stag, tstep_order, dtau, & tstep_neutrals_dens, tstep_neutrals_parmom, & tstep_neutrals_pressure) !! Initialises the timesteppers. This automatically will use a lower-order method for the first few timesteps !! to fill the storage type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) integer, intent(in) :: tstep_order !! Order of time-stepping scheme real(GP), intent(in) :: dtau !! Time-step size (base, not bdf) type(karniadakis_t), intent(inout) :: tstep_neutrals_dens !! Time integrator for neutrals density equation type(karniadakis_t), intent(inout) :: tstep_neutrals_parmom !! Time integrator for neutrals parallel momentum equation type(karniadakis_t), intent(inout) :: tstep_neutrals_pressure !! Time integrator for neutrals pressure equation integer :: np_cano, np_stag np_cano = mesh_cano%get_n_points() np_stag = mesh_stag%get_n_points() call tstep_neutrals_dens%init(np_cano, tstep_order, dtau) call tstep_neutrals_parmom%init(np_stag, tstep_order, dtau) call tstep_neutrals_pressure%init(np_cano, tstep_order, dtau) end subroutine subroutine init_neutrals_halos(comm_handler, neutrals_dens, neutrals_parmom, neutrals_pressure) !! Allocates halo points of Neutrals fields type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(variable_t), intent(inout) :: neutrals_dens !! Neutrals density type(variable_t), intent(inout) :: neutrals_parmom !! Neutrals parallel momentum type(variable_t), intent(inout) :: neutrals_pressure !! Neutrals pressure call neutrals_dens%create_halos(comm_handler, 1, 1) call neutrals_parmom%create_halos(comm_handler, 0, 1) call neutrals_pressure%create_halos(comm_handler, 1, 1) end subroutine end module