module model_apar_shift_m !! Implementation of apar_shift model use MPI use comm_handler_m, only : comm_handler_t use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use precision_grillix_m, only : GP, MPI_GP, GP_EPS use variable_m, only : variable_t use snapshot_m, only : snapshot_t, char_arr_t use multistep_m, only : karniadakis_t use snapshots_helpers_m, only : snaps_average use parallelisation_setup_m, only : is_rank_info_writer ! From PARALLAX use perf_m, only : perf_start, perf_stop use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use screen_io_m, only : get_stdout ! From parameters use params_tstep_m, only : tstep_order, dtau use params_apar_shift_m, only : apar_shift_iselect, tau_filter_width, & nsnaps_start, nsnaps_end implicit none real(GP), dimension(:), allocatable, private, save :: dapar_shift !! Explicit changerate of apar_shift, always equal to zero type, public :: apar_shift_t !! Type containing necessary objects and subroutines for apar_shift type(variable_t), public :: apar_shift !! Background shift in apar type(snapshot_t), public :: snaps_apar_shift !! Snapshots for apar_shift type(karniadakis_t), private :: tstep_apar_shift !! Time-step integrator for apar_shift contains procedure, public :: init => initialise_apar_shift procedure, public :: eval_fluct => evaluate_apar_fluct procedure, public :: write_snaps procedure, public :: timestep => timestep_apar_shift end type contains subroutine initialise_apar_shift(self, comm_handler, equi, mesh_cano, mesh_stag, & apar, tau, isnaps_apar_shift, apar_shift_firsts, dapar_shift_firsts) !! Initialise apar shift variable, snapshot, and timestepper class(apar_shift_t), intent(inout) :: self !! Instance of class 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) within poloidal plane type(mesh_cart_t), intent(inout) :: mesh_stag !! Mesh (staggered) within poloidal plane type(variable_t), intent(in) :: apar !! Parallel electromagnetic potential real(GP), intent(in) :: tau !! Time integer, intent(out) :: isnaps_apar_shift !! Snapshot of apar shift model real(GP), dimension(mesh_stag%get_n_points(), tstep_order-1), intent(in), optional :: apar_shift_firsts !! Initial storage of apar_shift at previous timesteps to achieve high order accuracy real(GP), dimension(mesh_stag%get_n_points(), tstep_order-1), intent(in), optional :: dapar_shift_firsts !! Initial storage of dapar_shift at previous timesteps to achieve high order accuracy integer :: nvars, l type(char_arr_t), allocatable, dimension(:) :: varnames real(GP) :: tau_apar_shift type(variable_t), allocatable, dimension(:) :: tmp_var type(variable_t) :: apar_avg ! Set snapsfiles nvars = 1 allocate(varnames(nvars)) varnames(1)%str = 'apar_shift' call self%snaps_apar_shift%init(comm_handler, mesh_cano, mesh_stag, 'snapsdir_apar_shift', nvars, varnames) deallocate(varnames) ! Set variables call self%apar_shift%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) ! Build initial state select case(apar_shift_iselect) case('CONTINUE') ! Continue from the last apar_shift snapshot allocate(tmp_var(1)) call self%snaps_apar_shift%read_snaps_last(comm_handler, isnaps_apar_shift, tau_apar_shift, tmp_var) self%apar_shift = tmp_var(1) deallocate(tmp_var) if (abs(tau_apar_shift-tau) >= GP_EPS) then call handle_error('Cannot continue due to conflict tau', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('tau, tau_apar_shift: ',[tau, tau_apar_shift])) endif if (is_rank_info_writer) then write(get_stdout(),*)'Continue to evolve apar_shift at--------------------' write(get_stdout(),201)isnaps_apar_shift, tau_apar_shift 201 FORMAT(1X,'isnaps_apar_shift = ',I8 /, & 1X,'tau_apar_shift = ',ES14.6E3) write(get_stdout(),*)'----------------------------------------------------' endif case('APAR') ! Set apar_shift as the latest apar !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, self, apar) !$omp do do l = 1, mesh_stag%get_n_points() self%apar_shift%vals(l) = apar%vals(l) enddo !$omp end do !$omp end parallel if (is_rank_info_writer) then write(get_stdout(),*)'Build the initial profile apar_shift = apar----------' endif case('AVERAGE') ! Set apar_shift as the averaged apar call snaps_average(comm_handler, mesh_stag, 'snapsdir_braginskii', & nsnaps_start, nsnaps_end, 'apar', self%apar_shift%vals) if (is_rank_info_writer) then write(get_stdout(),*)'Average apar over the existing snapshots --------------' write(get_stdout(),202)nsnaps_start, nsnaps_end 202 FORMAT(1X,'from snapshot ',I8 /, & 1X,'until snapshot ',I8) write(get_stdout(),*)'Build the initial profile apar_shift as the average----' endif case default call handle_error('apar_shift_iselect not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(apar_shift_iselect)) end select ! Initialise time-step integrators if(.not. allocated(dapar_shift) ) allocate(dapar_shift(mesh_stag%get_n_points())) !$omp parallel default(none) private(l) shared(mesh_stag, dapar_shift) !$omp do do l = 1, mesh_stag%get_n_points() dapar_shift(l) = 0.0_GP enddo !$omp end do !$omp end parallel if ( (present(apar_shift_firsts)) .and. (present(dapar_shift_firsts)) ) then call self%tstep_apar_shift%init(mesh_stag%get_n_points(), tstep_order, dtau, & apar_shift_firsts, dapar_shift_firsts) else call self%tstep_apar_shift%init(mesh_stag%get_n_points(), tstep_order, dtau) endif end subroutine subroutine evaluate_apar_fluct(self, comm_handler, equi, mesh_stag, apar, apar_fluct) !! Evaluate the fluctuating part of parallel electromagnetic potential class(apar_shift_t), intent(inout) :: self !! Instance of class type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh_stag !! Mesh (staggered) within poloidal plane type(variable_t), intent(in) :: apar !! Full parallel electromagnetic potential returned from Ohm's law type(variable_t), intent(inout) :: apar_fluct !! Fluctuating part of apar to be used by flutter terms integer :: l, ierr real(GP), dimension(mesh_stag%get_n_points()) :: axisymmetric_shift ! Store the toroidally averaged apar_shift if (equi%is_axisymmetric()) then call MPI_allreduce(self%apar_shift%vals, axisymmetric_shift, mesh_stag%get_n_points(), & MPI_GP, MPI_SUM, comm_handler%get_comm_planes(), ierr) !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, comm_handler, axisymmetric_shift, & !$omp apar, apar_fluct) !$omp do do l = 1, mesh_stag%get_n_points() axisymmetric_shift(l) = axisymmetric_shift(l) / comm_handler%get_nplanes() apar_fluct%vals(l) = apar%vals(l) - axisymmetric_shift(l) enddo !$omp end do !$omp end parallel else !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, apar, self, apar_fluct) !$omp do do l = 1, mesh_stag%get_n_points() apar_fluct%vals(l) = apar%vals(l) - self%apar_shift%vals(l) enddo !$omp end do !$omp end parallel endif end subroutine subroutine write_snaps(self, comm_handler, isnaps_apar_shift, tau) !! Writes apar_shift to snapshots class(apar_shift_t), intent(inout) :: self !! Instance of class type(comm_handler_t), intent(in) :: comm_handler !! Communicators integer, intent(in) :: isnaps_apar_shift !! Snapshot of apar shift model real(GP), intent(in) :: tau !! Time if (is_rank_info_writer) then write(get_stdout(),*)'Writing apar_shift snapshot at -----------------------------' write(get_stdout(),203)tau, isnaps_apar_shift 203 FORMAT( 3X,'tau = ',ES14.6E3 /, & 3X,'isnaps_apar_shift = ',I8) endif call self%snaps_apar_shift%write_snaps(comm_handler, isnaps_apar_shift, & tau, [self%apar_shift]) end subroutine subroutine timestep_apar_shift(self, comm_handler, mesh_stag, tau, apar) !! Advances apar_shift from t to t+1 class(apar_shift_t), intent(inout) :: self !! Instance of class type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) within poloidal plane real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: apar !! Parallel electromagnetic potential at timestep t real(GP) :: tau_adv, dtau_bdf real(GP), dimension(mesh_stag%get_n_points()) :: apar_shift_adv ! apar_shift explicit changerate and advanced in time at t+1 integer :: l call perf_start('timestep_apar_shift') tau_adv = tau + dtau ! Force apar_shift = apar without timestepping if tau_filter_width<=0 if (tau_filter_width < GP_EPS) then !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, self, apar) !$omp do do l = 1, mesh_stag%get_n_points() self%apar_shift%vals(l) = apar%vals(l) enddo !$omp end do !$omp end parallel call perf_stop('timestep_apar_shift') return endif ! Implicit timestepping ----------------------------- call self%tstep_apar_shift%advance(dapar_shift, self%apar_shift%vals, apar_shift_adv) dtau_bdf = self%tstep_apar_shift%get_dtau_bdf() !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, self, apar_shift_adv, apar, tau_filter_width, dtau_bdf) !$omp do do l = 1, mesh_stag%get_n_points() apar_shift_adv(l) = (apar_shift_adv(l) + dtau_bdf * apar%vals(l) / tau_filter_width ) & / (1.0_GP + dtau_bdf / tau_filter_width) enddo !$omp end do !$omp end parallel ! Formal advancement of time ------------------------------------------- call self%tstep_apar_shift%shift_storage( [self%apar_shift%vals, dapar_shift] ) !$omp parallel default(none) private(l) shared(mesh_stag, self, apar_shift_adv) !$omp do do l = 1, mesh_stag%get_n_points() self%apar_shift%vals(l)=apar_shift_adv(l) enddo !$omp end do !$omp end parallel call perf_stop('timestep_apar_shift') end subroutine end module