module init_zhdanov_m !! Sets initial state for ZHDANOV model use comm_handler_m, only : comm_handler_t use parallelisation_setup_m, only : is_rank_info_writer use precision_grillix_m, only : GP use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t use params_zhdanov_general_m, only : params_zhdanov_general_t use params_zhdanov_species_m, only : params_zhdanov_species_t use variable_m, only: variable_t use multistep_m, only : karniadakis_t use compute_electron_density_m, only : compute_electron_density use mms_zhdanov_m, only : MmsSolDens ! From PARALLAX use screen_io_m, only : get_stdout use constants_m, only : PI use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t implicit none public :: build_zhdanov_initial_state interface module subroutine init_zhdanov_state_profile(comm_handler, & equi, mesh, map, & penalisation, & params_species, & paramfile_zhdanov_species, & dens_species) !! Builds profile state type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation type(params_zhdanov_species_t), intent(in) :: params_species !! Species paramters for the Zhdanov model character(len=*), intent(in) :: paramfile_zhdanov_species !! Filename, where parameters are read from type(variable_t), intent(inout) :: dens_species !! Species density end subroutine end interface contains subroutine build_zhdanov_initial_state(comm_handler, equi, mesh, map, & penalisation, & params_general, & params_species, & paramfile_zhdanov_species, & tau, isnaps, & dens_species) !! 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 !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation type(params_zhdanov_general_t), intent(in) :: params_general !! General parameters related with ZHDANOV model type(params_zhdanov_species_t), intent(in) :: params_species !! Species paramters for the Zhdanov model character(len=*), intent(in) :: paramfile_zhdanov_species !! Filename, where parameters are read from real(GP), intent(out) :: tau !! Time of initial state integer, intent(out) :: isnaps !! Snapshot frame number type(variable_t), intent(inout) :: dens_species !! Species density integer :: ierr,io_error character(len=256) :: io_errmsg character(LEN=32) :: init_select integer :: l real(GP) :: x, y, phi real(GP), dimension(mesh%get_n_points()) :: charge_dens ! Check for MMS initialisation if (params_general%mms_on) then init_select = 'MMS' else init_select = params_species%init_select endif select case(init_select) case('PROFILE') isnaps = 1 tau = 0.0_GP if (is_rank_info_writer) then write(get_stdout(),*)'Setting up profile as initial state' endif ! Initialise variables call dens_species%init(comm_handler, mesh%get_n_points(), & staggered = .false.) if (is_rank_info_writer) then write(get_stdout(),*)' ...setting up profiles for dens_species' endif if (comm_handler%get_species() /= 0) then ! Initialisation of ion species call init_zhdanov_state_profile(comm_handler, & equi, mesh, map, & penalisation, params_species, & paramfile_zhdanov_species, & dens_species) else ! Initialisation of electron species with zero ! -> will be set at the end with MPI_reduce !$omp parallel do private(l) do l = 1, mesh%get_n_points() dens_species%vals(l) = 0.0_GP enddo !$omp end parallel do endif case('MMS') isnaps = 1 tau = 0.0_GP !Initialise variables call dens_species%init(comm_handler, mesh%get_n_points(), & staggered = .false.) phi = mesh%get_phi() !$omp parallel do private(l, x, y) do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) dens_species%vals(l) = MmsSolDens(equi,& comm_handler%get_species(), x, y, phi, tau) enddo !$omp end parallel do return case default write(get_stdout(),*)'error(build_initial_state): init_select not & valid',params_species%init_select error stop end select call compute_electron_density(comm_handler, mesh, & params_species, dens_species%vals) end subroutine subroutine initialise_timesteppers(comm_handler, equi, mesh,& params_general, tstep_continuity) !! Initialises the timesteppers type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(params_zhdanov_general_t), intent(in) :: params_general !! General parameters related with ZHDANOV model type(karniadakis_t), intent(inout) :: tstep_continuity !! Time integrator for continuity equation ! Initialise timesteppers. The order of the timestepping scheme ! will automatically be reduced for the first time-steps integer :: iorder, l real(GP) :: tau, x, y, phi, dtau_eps real(GP), dimension(mesh%get_n_points(), params_general%tstep_order) ::& log_dens_species_first, d_log_dens_species_first if (.not.params_general%mms_on) then call tstep_continuity%init(mesh%get_n_points(),& params_general%tstep_order, params_general%dtau) else phi = mesh%get_phi() dtau_eps = params_general%dtau * 1.0E-4_GP !$omp parallel private(l, x, y) do iorder = 1, params_general%tstep_order tau = -1.0_GP * params_general%dtau * iorder !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) log_dens_species_first(l, iorder) = & log( MmsSolDens(equi, comm_handler%get_species(),& x, y, phi, tau) ) d_log_dens_species_first(l, iorder) = ( MmsSolDens(equi,& comm_handler%get_species(), x, y, phi, tau+dtau_eps) -& MmsSolDens(equi, comm_handler%get_species(), x, y, phi,& tau-dtau_eps) ) /(2*dtau_eps) / & MmsSolDens(equi, comm_handler%get_species(),& x, y, phi, tau) enddo !$omp end do enddo !$omp end parallel call tstep_continuity%init(mesh%get_n_points(),& params_general%tstep_order, params_general%dtau, & log_dens_species_first, d_log_dens_species_first) endif end subroutine end module