module diagnostics_scalar_m !! Module for computing scalar diagnostics for the BRAGINSKII model use netcdf use MPI use comm_handler_m, only : comm_handler_t use variable_m, only : variable_t use precision_grillix_m, only : GP use perf_m, only : perf_start, perf_stop use error_handling_grillix_m, only: handle_error_netcdf, handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use mesh_cart_m, only: mesh_cart_t use parallel_map_m, only : parallel_map_t use inplane_operators_m, only : inplane_operators_t use equilibrium_m, only: equilibrium_t use penalisation_m, only : penalisation_t use polars_m, only : polars_t use equilibrium_storage_m, only : equilibrium_storage_t use parallel_target_flux_m, only : parallel_target_flux_t use diagnostics_helpers_m, only: inner_product, total_volume use diagnostics_group_m, only : diagnostics_group_t use diagnostics_packet_m, only : diagnostics_packet_t use params_brag_model_m, only : beta implicit none type, public, extends(diagnostics_group_t) :: & diagnostics_scalar_t contains procedure, public :: process_to_file procedure, public :: project_scalar end type contains subroutine process_to_file(self, diag_packet, comm_handler, equi, equi_on_cano, equi_on_stag, mesh_cano, & mesh_stag, map, penalisation_cano, penalisation_stag, polars_cano, polars_stag, & tau, ne, te, ti, pot, vort, upar, jpar, apar, apar_fluct, & neutrals_dens, neutrals_parmom, neutrals_pressure, & isnaps, idiag, start_new_file) !! Main diagnostics computation routine class(diagnostics_scalar_t), intent(inout) :: self !! Instance of class type(diagnostics_packet_t), intent(inout) :: diag_packet !! Diagnostic packet type(comm_handler_t), intent(in) :: comm_handler !! MPI communication handler class(equilibrium_t), intent(in) :: equi !! Magnetic equilibrium class(equilibrium_storage_t), intent(in) :: equi_on_cano !! Equilibrium storage on canonical plane enabling faster performance at certain locations class(equilibrium_storage_t), intent(in) :: equi_on_stag !! Equilibrium storage on staggered plane enabling faster performance at certain locations type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_cano !! Penalisation (canonical) type(penalisation_t), intent(in) :: penalisation_stag !! Penalisation (staggered) type(polars_t), intent(in) :: polars_cano !! Polar grid and operators (canonical) type(polars_t), intent(in) :: polars_stag !! Polar grid and operators (staggered) real(GP), intent(in) :: tau !! Time type(variable_t), intent(in) :: ne !! Electron density type(variable_t), intent(in) :: te !! Electron temperature type(variable_t), intent(in) :: ti !! Ion temperature type(variable_t), intent(in) :: pot !! Electrostatic potential type(variable_t), intent(in) :: vort !! Generalised vorticity type(variable_t), intent(in) :: upar !! Parallel ion velocity type(variable_t), intent(in) :: jpar !! Parallel current type(variable_t), intent(in) :: apar !! Parallel electromagnetic potential type(variable_t), intent(in) :: apar_fluct !! Parallel electromagnetic potential fluctation type(variable_t), intent(in) :: neutrals_dens !! Neutrals density type(variable_t), intent(in) :: neutrals_parmom !! Neutrals parallel momentum type(variable_t), intent(in) :: neutrals_pressure !! Neutrals pressure integer, intent(in) :: isnaps !! Braginskii model snapshot number integer, intent(in) :: idiag !! Diagnostic snapshot number logical, intent(in) :: start_new_file !! Flag whether to write to new file integer :: ind !! Reference index of last initialized diagnostic integer :: l real(GP), dimension(mesh_cano%get_n_points()) :: aux_cano real(GP), dimension(mesh_stag%get_n_points()) :: aux_stag call perf_start('diagnostics_braginskii_scalar') call self%allocate_diags(n_diags=52) ! << ADJUST WHEN ADDING DIAGNOSTICS ! Time call self%init_next_diagnostic(1, 'dim_scalar', 'tau', ind) self%diags(ind)%vals(1) = tau ! Volume (calculated on cano and stag) call self%init_next_diagnostic(1, 'dim_scalar', 'volume_cano', ind) self%diags(ind)%vals = total_volume(comm_handler, mesh_cano, map, penalisation_cano, & is_stag=.false., use_full_domain=.false.) call self%init_next_diagnostic(1, 'dim_scalar', 'volume_stag', ind) self%diags(ind)%vals = total_volume(comm_handler, mesh_stag, map, penalisation_stag, & is_stag=.true., use_full_domain=.false.) ! Volume averages of base variables call self%init_next_diagnostic(1, 'dim_scalar', 'ne_avg', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & ne%vals, is_stag=.false., mode='VOLAVG') call self%init_next_diagnostic(1, 'dim_scalar', 'te_avg', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & te%vals, is_stag=.false., mode='VOLAVG') call self%init_next_diagnostic(1, 'dim_scalar', 'ti_avg', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & ti%vals, is_stag=.false., mode='VOLAVG') call self%init_next_diagnostic(1, 'dim_scalar', 'pot_avg', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & pot%vals, is_stag=.false., mode='VOLAVG') call self%init_next_diagnostic(1, 'dim_scalar', 'upar_avg', ind) call self%project_scalar(comm_handler, mesh_stag, map, penalisation_stag, ind, & upar%vals, is_stag=.true., mode='VOLAVG') call self%init_next_diagnostic(1, 'dim_scalar', 'neutrals_dens_avg', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & neutrals_dens%vals, is_stag=.false., mode='VOLAVG') ! Volume integrals of external sources call self%init_next_diagnostic(1, 'dim_scalar', 'src_ne_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_ext_ne, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'src_te_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_ext_te, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'src_ti_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_ext_ti, is_stag=.false., mode='VOLINT') ! Volume integrals of heating power and stored energy change rate call self%init_next_diagnostic(1, 'dim_scalar', 'P_heat_e_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_heat_e, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_heat_i_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_heat_i, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'ddt_W_MHD', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ddt_eth, is_stag=.false., mode='VOLINT') ! Volume integrals of radiated power due to impurities and main plasma recombination call self%init_next_diagnostic(1, 'dim_scalar', 'P_rad_rec_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_rad_rec, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_rad_imp_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_rad_imp, is_stag=.false., mode='VOLINT') ! Volume integral of thermal and electron kinetic energy densities call self%init_next_diagnostic(1, 'dim_scalar', 'W_MHD', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%eth, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'ekin_perp_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ekin_perp, is_stag=.false., mode='VOLINT') ! Volume integral of energy sources due to neutrals interactions call self%init_next_diagnostic(1, 'dim_scalar', 'P_rad_iz_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_rad_iz, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_neut_e_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_neut_e, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_neut_i_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_neut_i, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_ext_perp_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_ext_perp, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_neut_perp_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_neut_perp, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'P_ext_par_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_ext_par, is_stag=.false., mode='VOLINT') ! NOTE: P_neut_par_int is computed from a part living on the canonical mesh and ! one part living on the staggered. Therefore, we have to add those two parts. call self%init_next_diagnostic(1, 'dim_scalar', 'P_neut_par_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%p_neut_par_cano, is_stag=.false., mode='VOLINT') call self%project_scalar(comm_handler, mesh_stag, map, penalisation_stag, ind, & diag_packet%p_neut_par_stag, is_stag=.true., mode='VOLINT', & accumulate=.true.) ! Volume integrals of changerates of base variables call self%init_next_diagnostic(1, 'dim_scalar', 'ddt_ne_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ddt_ne, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'ddt_te_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ddt_te, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'ddt_ti_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ddt_ti, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'ddt_neutrals_dens_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%ddt_neutrals_dens, is_stag=.false., mode='VOLINT') ! Plasma-neutrals sources call self%init_next_diagnostic(1, 'dim_scalar', 'src_iz_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_iz_ne, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'src_rec_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_rec_ne, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'src_rcy_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_neut_rcy, is_stag=.false., mode='VOLINT') ! Numerical sources call self%init_next_diagnostic(1, 'dim_scalar', 'src_floor_ne_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_floor_ne, is_stag=.false., mode='VOLINT') call self%init_next_diagnostic(1, 'dim_scalar', 'src_floor_neutrals_dens_int', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%src_floor_neutrals_dens, is_stag=.false., mode='VOLINT') ! Boundary fluxes ! ... parallel call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_ne_fwd', ind) self%diags(ind)%vals = diag_packet%parflux_ne_fwd call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_ne_bwd', ind) self%diags(ind)%vals = diag_packet%parflux_ne_bwd call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_ne_fwd_weighted', ind) self%diags(ind)%vals = diag_packet%parflux_ne_fwd_weighted call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_ne_bwd_weighted', ind) self%diags(ind)%vals = diag_packet%parflux_ne_bwd_weighted call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_neutrals_fwd', ind) self%diags(ind)%vals = diag_packet%parflux_neutrals_fwd call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_neutrals_bwd', ind) self%diags(ind)%vals = diag_packet%parflux_neutrals_bwd call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_neutralsdiss_fwd', ind) self%diags(ind)%vals = diag_packet%parflux_neutralsdiss_fwd call self%init_next_diagnostic(1, 'dim_scalar', 'parflux_neutralsdiss_bwd', ind) self%diags(ind)%vals = diag_packet%parflux_neutralsdiss_bwd ! ...perpendicular call self%init_next_diagnostic(1, 'dim_scalar', 'perpflux_ne_exb', ind) self%diags(ind)%vals = diag_packet%perpflux_ne_exb call self%init_next_diagnostic(1, 'dim_scalar', 'perpflux_ne_dia', ind) self%diags(ind)%vals = diag_packet%perpflux_ne_dia call self%init_next_diagnostic(1, 'dim_scalar', 'perpflux_neutrals_outer', ind) self%diags(ind)%vals = diag_packet%perpflux_neutrals_outer call self%init_next_diagnostic(1, 'dim_scalar', 'perpflux_neutrals_volint', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & diag_packet%perpflux_neutrals_volint, is_stag=.false., mode='VOLINT') ! ------------------------ ! Energy dissipation terms ! -- Thermal energy, electrons !$omp parallel do default(none) private(l) shared(mesh_cano, diag_packet, ne, te, aux_cano) do l = 1, mesh_cano%get_n_points() aux_cano(l) = 1.5_GP * (diag_packet%diss_ne(l) * te%vals(l) + & diag_packet%diss_te(l) * ne%vals(l)) enddo !$omp end parallel do call self%init_next_diagnostic(1, 'dim_scalar', 'ethe_diss', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & aux_cano, is_stag=.false., mode='VOLINT') ! -- Thermal energy, ions !$omp parallel do default(none) private(l) shared(mesh_cano, diag_packet, ne, ti, aux_cano) do l = 1, mesh_cano%get_n_points() aux_cano(l) = 1.5_GP * (diag_packet%diss_ne(l) * ti%vals(l) + & diag_packet%diss_ti(l) * ne%vals(l)) enddo !$omp end parallel do call self%init_next_diagnostic(1, 'dim_scalar', 'ethi_diss', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & aux_cano, is_stag=.false., mode='VOLINT') ! -- Perpendicular kinetic energy !$omp parallel do default(none) private(l) shared(mesh_cano, diag_packet, ne, pot, aux_cano) do l = 1, mesh_cano%get_n_points() aux_cano(l) = diag_packet%ekin_perp(l) / ne%vals(l) * diag_packet%diss_ne(l) & - pot%vals(l) * diag_packet%diss_vort(l) enddo !$omp end parallel do call self%init_next_diagnostic(1, 'dim_scalar', 'ekinperp_diss', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & aux_cano, is_stag=.false., mode='VOLINT') ! -- Parallel kinetic energy ! Note: Here we add quantities on staggered and canonical mesh_cano ! The calculation is done separatly for both meshes and then added !$omp parallel default(none) private(l) & !$omp shared(equi, mesh_cano, mesh_stag, diag_packet, upar, aux_cano, aux_stag) !$omp do do l = 1, mesh_cano%get_n_points() aux_cano(l) = 0.5_GP * diag_packet%upar_cano(l)**2 * diag_packet%diss_ne(l) enddo !$omp end do !$omp do do l = 1, mesh_stag%get_n_points() aux_stag(l) = diag_packet%diss_upar(l) * diag_packet%ne_stag(l) * upar%vals(l) enddo !$omp end do !$omp end parallel call self%init_next_diagnostic(1, 'dim_scalar', 'ekinpar_diss', ind) call self%project_scalar(comm_handler, mesh_cano, map, penalisation_cano, ind, & aux_cano, is_stag=.false., mode='VOLINT') call self%project_scalar(comm_handler, mesh_stag, map, penalisation_stag, ind, & aux_stag, is_stag=.true., mode='VOLINT', accumulate=.true.) ! -- Electromagnetic energy !$omp parallel do default(none) private(l) shared(mesh_stag, diag_packet, jpar, aux_stag, beta) do l = 1, mesh_stag%get_n_points() aux_stag(l) = beta * jpar%vals(l) * diag_packet%diss_apar(l) enddo !$omp end parallel do call self%init_next_diagnostic(1, 'dim_scalar', 'eem_diss', ind) call self%project_scalar(comm_handler, mesh_stag, map, penalisation_stag, ind, & aux_stag, is_stag=.true., mode='VOLINT') call self%write_diagnostics(tau, isnaps, idiag, start_new_file) deallocate(self%diags) call perf_stop('diagnostics_braginskii_scalar') end subroutine subroutine project_scalar(self, comm_handler, mesh, map, penalisation, ind, u, is_stag, mode, accumulate) !! Project input onto scalar diagnostic dimension nrho class(diagnostics_scalar_t), intent(inout) :: self !! Instance of class type(comm_handler_t), intent(in) :: comm_handler !! MPI communication handler type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation integer, intent(in) :: ind !! Work index in diagnostics array real(GP), dimension(:), intent(in) :: u !! Input array logical, intent(in) :: is_stag !! Logical is true if located on the staggered mesh character(len=*), intent(in) :: mode !! Projection mode, select from : !! 'VOLAVG' - volume average !! 'VOLINT' - volume integral logical, intent(in), optional :: accumulate !! Logical is true if the result is added to the variable logical :: accumulation_mode = .false. real(GP) :: aux if (present(accumulate)) then accumulation_mode = accumulate endif select case(mode) case('VOLAVG') aux = inner_product(comm_handler, mesh, & map, penalisation, & is_stag, u, & use_vol_avg=.true., & use_abs=.false., & use_full_domain=.false.) case('VOLINT') aux = inner_product(comm_handler, mesh, & map, penalisation, & is_stag, u, & use_vol_avg=.false., & use_abs=.false., & use_full_domain=.false.) case default call handle_error('Invalid mode for scalar projection', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select if (accumulation_mode) then self%diags(ind)%vals(1) = self%diags(ind)%vals(1) + aux else self%diags(ind)%vals(1) = aux endif end subroutine end module