module model_braginskii_m !! Implementation of BRAGINSKII model use MPI use netcdf use comm_handler_m, only : comm_handler_t use parallelisation_setup_m, only : is_rank_info_writer use precision_grillix_m, only : GP, MPI_GP, GP_EPS use error_handling_grillix_m, only: handle_error_netcdf, handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_CMD, GRILLIX_WRN_TSTEP_NO_SUCCESS use static_data_m, only : equi, & equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & hsolver_cano => field_solver_cano, & hsolver_stag => field_solver_stag, & map, & penalisation_cano => parbnd_immersed_cano, & penalisation_stag => parbnd_immersed_stag, & polars_cano, polars_stag, & parflux_utils_cano, parflux_utils_stag, & perp_bnd_flux_cano, perp_bnd_flux_stag use variable_m, only : variable_t use gyroviscosity_m, only : gyroviscosity_t use multistep_m, only : karniadakis_t use inplane_operators_m, only : inplane_operators_t use init_braginskii_m, only : build_braginskii_initial_state, & initialise_timesteppers, initialise_timesteppers_with_mms, & init_braginskii_halos use boundaries_braginskii_m, only : boundaries_braginskii_t use create_buffer_m, only : create_buffer, buffer_write_netcdf use timestep_braginskii_m, only : timestep_braginskii use checkpoint_monitor_m, only : checkpoint_monitor_t use snapshot_m, only : snapshot_t, char_arr_t use sources_external_m, only : sources_external_t use mms_braginskii_m, only : mms_diagnostics_braginskii, & mms_sol_braginskii_ne, & mms_sol_braginskii_pot, & mms_sol_braginskii_vort, & mms_sol_braginskii_upar, & mms_sol_braginskii_jpar, & mms_sol_braginskii_apar, & mms_sol_braginskii_te, & mms_sol_braginskii_ti use checksum_braginskii_m, only : checksum_braginskii ! From NEUTRALS module use init_neutrals_m, only : build_neutrals_initial_state, & initialise_neutrals_timesteppers, & init_neutrals_halos use boundaries_neutrals_m, only : boundaries_neutrals_t use neutrals_module_m, only : neutrals_module_t use mms_neutrals_m, only : mms_diagnostics_neutrals, & mms_sol_neutrals_dens, & mms_sol_neutrals_parmom, & mms_sol_neutrals_pressure ! From DIAGNOSTICS module use diagnostics_braginskii_m, only : diagnostics_braginskii_t ! From MAGNETIC_SHIFT module use model_apar_shift_m, only : apar_shift_t ! From PARALLAX use perf_m, only : perf_print, perf_reset use screen_io_m, only : get_stdout ! IOL module #ifdef ENABLE_IOL use iol_source_m, only : iol_source_t #endif ! Parameters use facade_params_m, only : & tinfo_size, & path_diag_lineout, & read_paths, write_paths, & read_all_params_braginskii, write_all_params_braginskii, & read_all_params_neutrals, write_all_params_neutrals use params_feature_selection_m, only : & stop_after_init, neutrals_on, mms_on, checksum_on, iol_on, & analyse_snapshots_on, analyse_snapshots_start, analyse_snapshots_end, evol_apar_shift_on use params_brag_switches_m, only : swb_flutter_main use params_tstep_m, only : & dtau, tau_snaps, tau_fin, tstep_order, tstep_dbgout use params_brag_model_m, only : & rhos, delta, beta use params_brag_boundaries_perp_m, only : bnddescr_te_core, bnddescr_ti_core use descriptors_braginskii_m, only : BND_BRAGTYPE_ZONAL_NEUMANN use params_brag_buffer_m, only : & buffer_select, buffer_path, buffer_dbgout use params_brag_numdiss_m, only : & perpdiss_exclude_buffer use params_brag_init_select_m, only : & brag_iselect, brag_init_path use params_neut_init_select_m, only : & neut_iselect use params_apar_shift_m, only : & apar_shift_iselect implicit none public :: model_braginskii private :: timestep_debug_output contains subroutine model_braginskii(comm_handler) type(comm_handler_t), intent(in) :: comm_handler !! Communication handler character(len=3) :: plane_c type(neutrals_module_t) :: neutrals !! Neutrals module object type(inplane_operators_t), target :: opsinplane_cano !! In-plane operators (canonical) type(inplane_operators_t), pointer :: opsinplane_stag !! In-plane operators (staggered) !! Points to opsinplane_cano for axisymmetric case !! and to opsinplane_stag_mem for non-axisymmetric case type(inplane_operators_t), target, allocatable :: opsinplane_stag_mem !! In-plane operators (staggered) !! only allocated for non-axisymetric case type(boundaries_braginskii_t) :: boundaries !! Boundaries related with BRAGINSKII model type(boundaries_neutrals_t) :: boundaries_neutrals !! Boundaries related with neutrals module type(diagnostics_braginskii_t) :: diagnostics !! Diagnostics type(checkpoint_monitor_t) :: clock !! Main checkpoint monitor type(apar_shift_t) :: ashift !! Apar shift object real(GP) :: tau !! Time type(variable_t) :: ne !! Electron density type(variable_t) :: pot !! Electrostatic potential type(variable_t) :: vort !! Generalised vorticity type(variable_t) :: upar !! Parallel ion velocity type(variable_t) :: jpar !! Parallel current type(variable_t) :: apar !! Parallel electromagnetic potential type(variable_t) :: apar_fluct !! Fluctation of apar used for flutter operators type(variable_t) :: te !! Electron temperature type(variable_t) :: ti !! Ion temperature type(gyroviscosity_t) :: gstress !! Gyroviscosity type(sources_external_t) :: sources_external !! External sources real(GP), dimension(mesh_cano%get_n_points_inner()) :: src_ne !! Particle source values on inner mesh points real(GP), dimension(mesh_cano%get_n_points_inner()) :: src_te !! Electron temperature source values on inner mesh points real(GP), dimension(mesh_cano%get_n_points_inner()) :: src_ti !! Ion temperature source values on inner mesh points real(GP), dimension(mesh_stag%get_n_points_inner()) :: src_upar !! Parallel momentum source values on inner mesh points real(GP), dimension(mesh_cano%get_n_points_inner()) :: src_vort !! Vorticity source values on inner mesh points type(variable_t) :: neutrals_dens !! Neutrals density type(variable_t) :: neutrals_parmom !! Neutrals parallel momentum type(variable_t) :: neutrals_pressure !! Neutrals pressure type(karniadakis_t) :: tstep_continuity !! Time integrator for continuity equation type(karniadakis_t) :: tstep_vort !! Time-step integrator for vorticity equation type(karniadakis_t) :: tstep_parmom !! Time-step integrator for parallel momentum equation type(karniadakis_t) :: tstep_ohm !! Time-step integrator for Ohm's law type(karniadakis_t) :: tstep_etemp !! Time-step integrator for electron temperature equation type(karniadakis_t) :: tstep_itemp !! Time-step integrator for ion temperature equation type(karniadakis_t) :: tstep_neutrals_dens !! Time integrator for neutrals continuity equation type(karniadakis_t) :: tstep_neutrals_parmom !! Time integrator for neutrals parallel momentum equation type(karniadakis_t) :: tstep_neutrals_pressure !! Time integrator for neutrals pressure equation #ifdef ENABLE_IOL type(iol_source_t) :: iol_source !! Ion orbit loss source (sink) logical :: iol_source_updated real(GP) :: tau_iol character(len=3), parameter :: snapsdir_iol = 'iol' character(len=:), allocatable :: snapsfile_iol integer :: cmd_err, ierr #endif real(GP), dimension(mesh_cano%get_n_points()), target :: cf_buffer_cano, cf_diss_cano !! Buffer and dissipation envleope function (canonical) real(GP), dimension(:), pointer :: cf_buffer_stag, cf_diss_stag !! Buffer amd dissipation envelope function (staggered) !! Points to cf_buffer_cano for axisymmetric case !! and to cf_buffer_cano_mem for non-axisymmetric case real(GP), allocatable, dimension(:), target :: cf_buffer_stag_mem, cf_diss_stag_mem !! Buffer and dissipation envelope function (staggered) !! only allocated for non-axisymetric case integer :: isnaps, isnaps_apar_shift, isnaps_neutrals, tstep_count integer, dimension(tinfo_size) :: tinfo !! Info status (integers) real(GP), dimension(tinfo_size,2) :: rinfo !! Info status (float) logical :: success_plasma, success_neutrals !! Success flags for timesteps integer :: nvars type(char_arr_t), allocatable, dimension(:) :: varnames type(snapshot_t) :: snaps_braginskii type(snapshot_t) :: snaps_neutrals integer :: nf90_id, nf90_grpid, nf90_stat integer :: l, ki if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'START RUNNING BRAGINSKII MODEL' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'' endif write(plane_c,'(I3.3)')comm_handler%get_plane() ! Read parameters call read_paths('params_braginskii.in') if (is_rank_info_writer) then call write_paths() endif call read_all_params_braginskii() if (neutrals_on) then call read_all_params_neutrals() endif if (is_rank_info_writer) then call write_all_params_braginskii() if (neutrals_on) then call write_all_params_neutrals() endif endif ! Initialise in-plane operators ----------------------------------------------------------- call opsinplane_cano%init(equi, mesh_cano, equi_on_cano, rhos, & norm_flutter=delta*beta*swb_flutter_main) if (equi%is_axisymmetric()) then opsinplane_stag => opsinplane_cano else allocate(opsinplane_stag_mem) call opsinplane_stag_mem%init(equi, mesh_stag, equi_on_stag, rhos, & norm_flutter=delta*beta*swb_flutter_main) opsinplane_stag => opsinplane_stag_mem endif ! Initialise sources ---------------------------------------------------------------------- call sources_external%init() ! Create buffer zones and dissipation envelope functions ---------------------------------- call create_buffer(equi, mesh_cano, buffer_select, buffer_path, cf_buffer_cano) !$omp parallel default(none) private(l) & !$omp shared(mesh_cano, perpdiss_exclude_buffer, cf_buffer_cano, cf_diss_cano) if (perpdiss_exclude_buffer) then !$omp do do l = 1, mesh_cano%get_n_points() cf_diss_cano(l) = 1.0_GP - cf_buffer_cano(l) enddo !$omp end do else !$omp do do l = 1, mesh_cano%get_n_points() cf_diss_cano(l) = 1.0_GP enddo !$omp end do endif !$omp end parallel if (equi%is_axisymmetric()) then cf_buffer_stag => cf_buffer_cano cf_diss_stag => cf_diss_cano if ( (buffer_dbgout) .and. (is_rank_info_writer) ) then nf90_stat = nf90_create('buffer.nc', NF90_NETCDF4+NF90_CLOBBER, nf90_id ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call buffer_write_netcdf(mesh_cano, cf_buffer_cano, nf90_id) nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif else allocate(cf_buffer_stag_mem(mesh_stag%get_n_points())) call create_buffer(equi, mesh_stag, buffer_select, buffer_path, cf_buffer_stag_mem) cf_buffer_stag => cf_buffer_stag_mem allocate(cf_diss_stag_mem(mesh_stag%get_n_points())) !$omp parallel default(none) private(l) & !$omp shared(mesh_stag, perpdiss_exclude_buffer, cf_buffer_stag, cf_diss_stag_mem) if (perpdiss_exclude_buffer) then !$omp do do l = 1, mesh_stag%get_n_points() cf_diss_stag_mem(l) = 1.0_GP - cf_buffer_stag(l) enddo !$omp end do else !$omp do do l = 1, mesh_stag%get_n_points() cf_diss_stag_mem(l) = 1.0_GP enddo !$omp end do endif !$omp end parallel cf_diss_stag => cf_diss_stag_mem if (buffer_dbgout) then nf90_stat = nf90_create('buffer_plane'//plane_c//'.nc', NF90_NETCDF4, nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) nf90_stat = nf90_def_grp(nf90_id, 'canonical', nf90_grpid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call buffer_write_netcdf(mesh_cano, cf_buffer_cano, nf90_grpid) nf90_stat = nf90_def_grp(nf90_id, 'staggered', nf90_grpid) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call buffer_write_netcdf(mesh_stag, cf_buffer_stag, nf90_grpid) nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif endif ! Set snapsfiles --------------------------------------------------------------------------- if (mms_on) then nvars = 16 allocate(varnames(nvars)) varnames(9)%str = 'mms_ne' varnames(10)%str = 'mms_te' varnames(11)%str = 'mms_ti' varnames(12)%str = 'mms_pot' varnames(13)%str = 'mms_vort' varnames(14)%str = 'mms_upar' varnames(15)%str = 'mms_jpar' varnames(16)%str = 'mms_apar' else nvars = 8 allocate(varnames(nvars)) endif varnames(1)%str = 'ne' varnames(2)%str = 'te' varnames(3)%str = 'ti' varnames(4)%str = 'pot' varnames(5)%str = 'vort' varnames(6)%str = 'upar' varnames(7)%str = 'jpar' varnames(8)%str = 'apar' call snaps_braginskii%init(comm_handler, mesh_cano, mesh_stag, 'snapsdir_braginskii', nvars, & varnames) deallocate(varnames) if ( neutrals_on ) then if (mms_on) then nvars = 6 allocate(varnames(nvars)) varnames(4)%str = 'mms_neutrals_dens' varnames(5)%str = 'mms_neutrals_parmom' varnames(6)%str = 'mms_neutrals_pressure' else nvars = 3 allocate(varnames(nvars)) endif varnames(1)%str = 'neutrals_dens' varnames(2)%str = 'neutrals_parmom' varnames(3)%str = 'neutrals_pressure' call snaps_neutrals%init(comm_handler, mesh_cano, mesh_stag, 'snapsdir_neutrals', nvars, & varnames) deallocate(varnames) endif ! Setup boundary conditions --------------------------------------------------------------- !! Boundary conditions for the braginskii model call boundaries%init(mesh_cano, mesh_stag) !! Boundary conditions for the neutrals model if ( neutrals_on ) then call boundaries_neutrals%init(mesh_cano, mesh_stag) endif ! Build initial state --------------------------------------------------------------------- call build_braginskii_initial_state(comm_handler, & equi, equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & hsolver_cano, hsolver_stag, & map, & penalisation_cano, penalisation_stag, & polars_cano, polars_stag, & opsinplane_cano, opsinplane_stag, & brag_iselect, brag_init_path, & snaps_braginskii, & boundaries, & isnaps, tau, & ne, pot, vort, & upar, jpar, apar, apar_fluct, & te, ti) call init_braginskii_halos(comm_handler, ne, pot, vort, upar, jpar, apar, apar_fluct, te, ti) if ( neutrals_on ) then call build_neutrals_initial_state(comm_handler, equi, & mesh_cano, mesh_stag, & 'params_neutrals.in', & snaps_neutrals, & isnaps_neutrals, & ti, & neutrals_dens, & neutrals_parmom, & neutrals_pressure) call init_neutrals_halos(comm_handler, neutrals_dens, neutrals_parmom, neutrals_pressure) endif ! Initialise gyroviscosity ---------------------------------------------------------------- call gstress%init(comm_handler, mesh_cano, mesh_stag) ! Initialise checkpoint monitors ---------------------------------------------------------- call clock%init(major_interval=tau_snaps, & major_tally_offset=isnaps) ! Create netcdf snapsfiles ---------------------------------------------------------------- if (brag_iselect /= 'CONTINUE') then call write_braginskii() endif if ( (neutrals_on) .and. (neut_iselect /= 'CONTINUE') ) then call write_neutrals() endif if ( evol_apar_shift_on ) then call ashift%init(comm_handler, equi, mesh_cano, mesh_stag, apar, tau, isnaps_apar_shift) call ashift%eval_fluct(comm_handler, equi, mesh_stag, apar, apar_fluct) ! Write initial state if ( apar_shift_iselect /= 'CONTINUE' ) then call ashift%write_snaps(comm_handler, isnaps, tau) endif endif ! Initialise diagnostics modules ---------------------------------------------------------- call diagnostics%init(comm_handler, path_diag_lineout, equi, mesh_cano, mesh_stag, polars_cano, & polars_stag, tau, isnaps+1, ne, te, ti, pot, upar, jpar, apar, apar_fluct, & neutrals_dens, neutrals_parmom, neutrals_pressure) if (stop_after_init) then if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'INITIAL STATE SET UP SUCCESSFULLY' write(get_stdout(),*)'code returns as stop_after_init=true, you may want to check now the initial state' write(get_stdout(),*)'To continue simulation set init_type=CONTINUE and stop_after_init=false' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'' endif return endif ! Initialise plasma time-step integrators --------------------------------------------------- if (.not.mms_on) then ! Call tstep%init on each timestepper call initialise_timesteppers(mesh_cano, mesh_stag, tstep_order, dtau, & tstep_continuity, tstep_vort, tstep_parmom, & tstep_ohm, tstep_etemp, tstep_itemp) else ! Fill the storage with the MMS solution and derivative for t < 0, then initialise call initialise_timesteppers_with_mms(equi, mesh_cano, mesh_stag, tstep_order, dtau, & tstep_continuity, tstep_vort, tstep_parmom, & tstep_ohm, tstep_etemp, tstep_itemp) endif ! Initialise neutrals and related time-step integrators ----------------------------------- if ( neutrals_on ) then call neutrals%init(comm_handler, equi, mesh_cano, mesh_stag, map, penalisation_cano, ne, & neutrals_dens, neutrals_parmom, neutrals_pressure, isnaps_neutrals, tau) call initialise_neutrals_timesteppers(mesh_cano, mesh_stag, tstep_order, dtau, & tstep_neutrals_dens, tstep_neutrals_parmom, & tstep_neutrals_pressure) endif #ifdef ENABLE_IOL ! Initialise Ion-Orbit-Loss model---------------------------------------------------------- if (iol_on) then call iol_source%init('params_iol.in', equi, mesh_cano, polars_cano) ! Create folder for iol data and snapsfile if (is_rank_info_writer) then call EXECUTE_COMMAND_LINE('mkdir '//trim(snapsdir_iol), cmdstat = cmd_err) if (cmd_err /= 0) then call handle_error('Mkdir failed', & GRILLIX_ERR_CMD, __LINE__, __FILE__, & error_info_t('Cmd_err: ',[cmd_err])) endif endif call MPI_BARRIER(comm_handler%get_comm_cart(), ierr) snapsfile_iol = snapsdir_iol//'/source_iol.nc' ! Same file for all ranks if (iol_source%get_init_continue()) then if (comm_handler%get_plane() == 0) then nf90_stat = nf90_open(snapsfile_iol, NF90_NETCDF4 + NF90_WRITE, nf90_id ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call iol_source%read_netcdf(nf90_id, tau = tau_iol) nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif if (is_rank_info_writer) then write(get_stdout(),*) write(get_stdout(),*)'Continuing with IOL source from tau = ', tau_iol write(get_stdout(),*) endif call MPI_BARRIER(comm_handler%get_comm_cart(), ierr) endif call iol_source%display() endif #endif ! Snapshots-analysing loop ---------------------------------------------------------------------- if ( analyse_snapshots_on ) then if (is_rank_info_writer) then write(get_stdout(),*) write(get_stdout(),*)'start analysing snapshots .....' write(get_stdout(),*) endif ! Try to read and analyse every snapshot in the snaps folder do ki = analyse_snapshots_start, analyse_snapshots_end call analyse_snapshots(ki) enddo if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'SNAPSHOTS ANALYSIS RUN SUCCESFULLY' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'' endif return ! exit model_braginskii routine endif ! Time-stepping loop ---------------------------------------------------------------------- if (is_rank_info_writer) then write(get_stdout(),*) write(get_stdout(),*)'start timestepping .....' write(get_stdout(),*) endif tstep_count = 0 do while (tau <= tau_fin + dtau) tstep_count = tstep_count + 1 ! Set sources to zero !$omp parallel default(none), private(ki) & !$omp shared(mesh_cano, mesh_stag, & !$omp src_ne, src_te, src_ti, src_upar, src_vort) !$omp do do ki = 1, mesh_cano%get_n_points_inner() src_ne(ki) = 0.0_GP src_te(ki) = 0.0_GP src_ti(ki) = 0.0_GP src_vort(ki) = 0.0_GP enddo !$omp end do !$omp do do ki = 1, mesh_stag%get_n_points_inner() src_upar(ki) = 0.0_GP enddo !$omp end do !$omp end parallel ! Reset solver information tinfo = 0 rinfo = 0.0_GP success_neutrals = .true. success_plasma = .true. call sources_external%eval(comm_handler, equi_on_cano, mesh_cano, polars_cano, tau, & ne, pot, vort, upar, jpar, apar, te, ti, & src_ne, src_te, src_ti) ! Timestepping if ( neutrals_on ) then call neutrals%eval_sources(comm_handler, equi, equi_on_cano,mesh_cano, mesh_stag, map, & ne, pot, upar, te, ti, & neutrals_dens, neutrals_parmom, neutrals_pressure, & src_ne, src_upar, src_te, src_ti, src_vort) endif #ifdef ENABLE_IOL if (iol_on) then if ( (tstep_count == 1) .and. iol_source%get_init_continue() ) then ! For continuation, do not perform update of IOL source at first timestep call iol_source%driver(comm_handler, equi, mesh_cano, polars_cano, tau, dtau, ne, pot, ti, & src_vort, iol_source_updated, force_no_update = .true.) else call iol_source%driver(comm_handler, equi, mesh_cano, polars_cano, tau, dtau, ne, pot, ti, & src_vort, iol_source_updated) endif ! Write out IOL source if ( (iol_source_updated) .and. (comm_handler%get_plane() == 0) ) then if ( (tstep_count == 1) .and. (.not.iol_source%get_init_continue()) ) then nf90_stat = nf90_create(snapsfile_iol, NF90_NETCDF4 + NF90_NOCLOBBER, nf90_id ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call iol_source%write_netcdf(tau, fexist = .false., fgid = nf90_id) else nf90_stat = nf90_open(snapsfile_iol, NF90_NETCDF4 + NF90_WRITE, nf90_id ) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) call iol_source%write_netcdf(tau, fexist = .true., fgid = nf90_id) endif nf90_stat = nf90_close(nf90_id) call handle_error_netcdf(nf90_stat, __LINE__, __FILE__) endif endif #endif call timestep_braginskii(comm_handler, & equi, equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & hsolver_cano, hsolver_stag, & map, & penalisation_cano, penalisation_stag, & polars_cano, polars_stag, & opsinplane_cano, opsinplane_stag, & boundaries, & cf_buffer_cano, cf_buffer_stag, & cf_diss_cano, cf_diss_stag, & tau, & ne, src_ne, tstep_continuity, & pot, vort, src_vort, tstep_vort, & upar, src_upar, tstep_parmom, & jpar, apar, apar_fluct, tstep_ohm, & te, src_te, tstep_etemp, & ti, src_ti, tstep_itemp, & gstress, & tinfo, rinfo, success_plasma) if ( neutrals_on ) then call neutrals%timestep(comm_handler, & equi, equi_on_cano, equi_on_stag, & mesh_cano, mesh_stag, & hsolver_cano, hsolver_stag, & map, & penalisation_cano, penalisation_stag, & parflux_utils_cano, parflux_utils_stag, & perp_bnd_flux_cano, perp_bnd_flux_stag, & opsinplane_cano, opsinplane_stag, & boundaries_neutrals, tau, & ne, pot, upar, jpar, apar_fluct, te, ti, & neutrals_dens, neutrals_parmom, neutrals_pressure, & tstep_neutrals_dens, tstep_neutrals_parmom, & tstep_neutrals_pressure, & tinfo(8:11), rinfo(8:11,1), success_neutrals) endif if ( evol_apar_shift_on ) then call ashift%timestep(comm_handler, mesh_stag, tau, apar) call ashift%eval_fluct(comm_handler, equi, mesh_stag, apar, apar_fluct) endif ! Advance time tau = tau + dtau ! Diagnostics call diagnostics%drive( & comm_handler, equi, equi_on_cano, equi_on_stag, mesh_cano, mesh_stag, map, & penalisation_cano, penalisation_stag, polars_cano, polars_stag, parflux_utils_cano, parflux_utils_stag, & perp_bnd_flux_cano, perp_bnd_flux_stag, opsinplane_cano, opsinplane_stag, boundaries, boundaries_neutrals, & cf_buffer_cano, cf_buffer_stag, cf_diss_cano, cf_diss_stag, isnaps+1, tau, & ne, te, ti, pot, vort, upar, jpar, apar, apar_fluct, gstress, & neutrals_dens, neutrals_parmom, neutrals_pressure, & src_ne, src_te, src_ti, src_upar, src_vort, & sources_external) ! Check for success of timestep, write output in case of failure if ( (.not.success_neutrals).or.(.not.success_plasma) ) then call timestep_debug_output(comm_handler, tstep_count, tau, tinfo, rinfo) isnaps = 99999 call handle_error('Timestep returned without success', & GRILLIX_WRN_TSTEP_NO_SUCCESS, __LINE__, __FILE__, & error_info_t('tstep_count, isnaps, tau: ', & [tstep_count, isnaps],[tau])) call write_braginskii() if ( neutrals_on ) then isnaps_neutrals = 99999 call write_neutrals() endif if ( evol_apar_shift_on ) then isnaps_apar_shift = 99999 call ashift%write_snaps(comm_handler, isnaps_apar_shift, tau) endif return ! exit model_braginskii routine endif ! Advance checkpoint monitors call clock%drive(tau, dtau) ! Update isnaps isnaps = clock%get_major_tally() if ( neutrals_on ) then isnaps_neutrals = clock%get_major_tally() endif if ( evol_apar_shift_on ) then isnaps_apar_shift = clock%get_major_tally() endif ! Main snapshot/mms checkpoint if (clock%on_major_checkpoint()) then ! Write information to screen if (is_rank_info_writer) then write(get_stdout(),*)'Writing snapshot at -----------------------------------' write(get_stdout(),733)tau, tstep_count, isnaps 733 FORMAT( 3X,'tau = ',ES14.6E3 /, & 3X,'tstep_count = ',I8 /, & 3X,'isnaps = ',I8) if ( neutrals_on ) then write(get_stdout(),734)isnaps_neutrals 734 FORMAT(3X,'isnaps_neutrals = ',I8) endif write(get_stdout(),*)'-------------------------------------------------------' endif ! Write snapshots call write_braginskii() if (neutrals_on) then call write_neutrals() endif if ( evol_apar_shift_on ) then call ashift%write_snaps(comm_handler, isnaps_apar_shift, tau) endif ! Write MMS diagnostics if (mms_on) then call mms_diagnostics_braginskii(comm_handler, equi, mesh_cano, mesh_stag, map, tau, & ne, pot, vort, upar, jpar, apar, te, ti) if ( neutrals_on ) then call mms_diagnostics_neutrals(comm_handler, equi, mesh_cano, mesh_stag, map, tau, & neutrals_dens, neutrals_parmom, neutrals_pressure) endif endif ! Perform checksum (for CI/CD) if (checksum_on) then call checksum_braginskii(comm_handler, mesh_cano, tau, ne, pot, vort, upar, jpar, apar, te, ti) endif ! Write performance statistics call perf_print() call perf_reset() endif if (tstep_dbgout) then call timestep_debug_output(comm_handler, tstep_count, tau, tinfo, rinfo) endif enddo if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'BRAGINSKII MODEL RUN SUCCESFULLY' write(get_stdout(),*)'-------------------------------------------------------' write(get_stdout(),*)'' endif contains subroutine write_braginskii() !! Writes Braginskii fields to snapshots integer :: l real(GP) :: phi_cano, phi_stag, x, y type(variable_t) :: mms_ne, mms_pot, mms_vort, mms_upar, mms_jpar, mms_apar, mms_te, mms_ti if (mms_on) then ! Compute analytic MMS solution phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() call mms_ne%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call mms_pot%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call mms_vort%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call mms_upar%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call mms_jpar%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call mms_apar%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call mms_te%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call mms_ti%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) !$omp parallel default(none) private(l, x, y) & !$omp shared(equi, tau, mesh_cano, phi_cano, mesh_stag, phi_stag, & !$omp mms_ne, mms_te, mms_ti, mms_pot, mms_vort, & !$omp mms_upar, mms_jpar, mms_apar) !$omp do do l = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) mms_ne%vals(l) = mms_sol_braginskii_ne(equi, x, y, phi_cano, tau) mms_te%vals(l) = mms_sol_braginskii_te(equi, x, y, phi_cano, tau) mms_ti%vals(l) = mms_sol_braginskii_ti(equi, x, y, phi_cano, tau) mms_pot%vals(l) = mms_sol_braginskii_pot(equi, x, y, phi_cano, tau) mms_vort%vals(l) = mms_sol_braginskii_vort(equi, x, y, phi_cano, tau) 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) mms_upar%vals(l) = mms_sol_braginskii_upar(equi, x, y, phi_stag, tau) mms_jpar%vals(l) = mms_sol_braginskii_jpar(equi, x, y, phi_stag, tau) mms_apar%vals(l) = mms_sol_braginskii_apar(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp end parallel call snaps_braginskii%write_snaps(comm_handler, isnaps, tau, & [ne , te, ti, pot, vort, upar, jpar, apar, & mms_ne , mms_te, mms_ti, mms_pot, mms_vort, mms_upar, mms_jpar, mms_apar]) else call snaps_braginskii%write_snaps(comm_handler, isnaps, tau, & [ne , te, ti, pot, vort, upar, jpar, apar]) endif end subroutine subroutine write_neutrals() !! Writes Neutrals fields to snapshots integer :: l real(GP) :: phi_cano, phi_stag, x, y type(variable_t) :: mms_neutrals_dens, mms_neutrals_parmom, mms_neutrals_pressure if (mms_on) then ! Compute analytic MMS solution phi_cano = mesh_cano%get_phi() phi_stag = mesh_stag%get_phi() call mms_neutrals_dens%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) call mms_neutrals_parmom%init(comm_handler, mesh_stag%get_n_points(), staggered = .true.) call mms_neutrals_pressure%init(comm_handler, mesh_cano%get_n_points(), staggered = .false.) !$omp parallel default(none) private(l, x, y) & !$omp shared(equi, tau, mesh_cano, phi_cano, mesh_stag, phi_stag, & !$omp mms_neutrals_dens, mms_neutrals_parmom, mms_neutrals_pressure) !$omp do do l = 1, mesh_cano%get_n_points() x = mesh_cano%get_x(l) y = mesh_cano%get_y(l) mms_neutrals_dens%vals(l) = mms_sol_neutrals_dens(equi, x, y, phi_cano, tau) mms_neutrals_pressure%vals(l) = mms_sol_neutrals_pressure(equi, x, y, phi_cano, tau) 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) mms_neutrals_parmom%vals(l) = mms_sol_neutrals_parmom(equi, x, y, phi_stag, tau) enddo !$omp end do !$omp end parallel call snaps_neutrals%write_snaps(comm_handler, isnaps_neutrals, & tau, & [neutrals_dens, neutrals_parmom, neutrals_pressure, & mms_neutrals_dens, mms_neutrals_parmom, mms_neutrals_pressure]) else call snaps_neutrals%write_snaps(comm_handler, isnaps_neutrals, & tau, & [neutrals_dens, neutrals_parmom, neutrals_pressure]) endif end subroutine subroutine analyse_snapshots(i_analyse_snaps) !! Read and analyse snapshots, then write out diagnostics results integer, intent(in) :: i_analyse_snaps !! Snapshot number for analysising real(GP) :: tau_neutrals, tau_rms type(variable_t), allocatable, dimension(:) :: tmp_var logical :: fexist character(len=5) :: isnaps_c character(len=:), allocatable :: snapsfile integer :: l ! Local variables imported from snapshots type(variable_t) :: ne, te, ti, pot, vort, & upar, jpar, apar, apar_fluct, & neutrals_dens, neutrals_parmom, & neutrals_pressure ! Check whether this snapshot exists write(isnaps_c,'(I5.5)')i_analyse_snaps snapsfile = 'snapsdir_braginskii' // '/snaps_t' // isnaps_c // '.nc' inquire(file=snapsfile, exist=fexist) if (.not.fexist) then return ! exit subroutine analyse_snapshots endif ! Read snaps_braginskii folder allocate(tmp_var(8)) call snaps_braginskii%read_snaps(comm_handler, i_analyse_snaps, tau, tmp_var) ne = tmp_var(1) te = tmp_var(2) ti = tmp_var(3) pot = tmp_var(4) vort = tmp_var(5) upar = tmp_var(6) jpar = tmp_var(7) apar = tmp_var(8) deallocate(tmp_var) call apar_fluct%init(comm_handler, mesh_stag%get_n_points(), staggered=.true., init_vals=apar%vals) call init_braginskii_halos(comm_handler, ne, pot, vort, upar, jpar, apar, apar_fluct, te, ti) ! Read snaps_neutrals folder if ( neutrals_on ) then allocate(tmp_var(3)) call snaps_neutrals%read_snaps(comm_handler, i_analyse_snaps, tau_neutrals, tmp_var) neutrals_dens = tmp_var(1) neutrals_parmom = tmp_var(2) neutrals_pressure = tmp_var(3) deallocate(tmp_var) if ( abs(tau_neutrals-tau) > GP_EPS ) then call handle_error('conflicting tau between snaps_braginskii and snaps_neutrals', & GRILLIX_ERR_CMD, __LINE__, __FILE__, & error_info_t('tau, tau_neutrals: ',[tau, tau_neutrals])) endif call init_neutrals_halos(comm_handler, neutrals_dens, neutrals_parmom, neutrals_pressure) endif ! Read snaps_apar_shift folder if ( evol_apar_shift_on ) then allocate(tmp_var(1)) call ashift%snaps_apar_shift%read_snaps(comm_handler, i_analyse_snaps, tau_rms, tmp_var) !$omp parallel default(none) private(l) shared(mesh_stag, ashift, tmp_var) !$omp do do l = 1, mesh_stag%get_n_points() ashift%apar_shift%vals(l) = tmp_var(1)%vals(l) enddo !$omp end do !$omp end parallel deallocate(tmp_var) call ashift%eval_fluct(comm_handler, equi, mesh_stag, apar, apar_fluct) if ( abs(tau_rms-tau) > GP_EPS ) then call handle_error('conflicting tau between snaps_braginskii and snaps_apar_shift', & GRILLIX_ERR_CMD, __LINE__, __FILE__, & error_info_t('tau, tau_rms: ',[tau, tau_rms])) endif endif if (is_rank_info_writer) then write(get_stdout(),*)'Analysing snapshot at--------------------' write(get_stdout(),735)i_analyse_snaps, tau 735 FORMAT(1X,'isnaps = ',I8 /, & 1X,'tau = ',ES14.6E3) write(get_stdout(),*)'-----------------------------------------' endif ! Do diagnostics call diagnostics%drive(comm_handler, equi, equi_on_cano, equi_on_stag, mesh_cano, mesh_stag, map, & penalisation_cano, penalisation_stag, polars_cano, polars_stag, & parflux_utils_cano, parflux_utils_stag, perp_bnd_flux_cano, perp_bnd_flux_stag, & opsinplane_cano, opsinplane_stag, boundaries, boundaries_neutrals, & cf_buffer_cano, cf_buffer_stag, cf_diss_cano, cf_buffer_stag, i_analyse_snaps, tau, & ne, te, ti, pot, vort, upar, jpar, apar, apar_fluct, gstress, & neutrals_dens, neutrals_parmom, neutrals_pressure, & src_ne, src_te, src_ti, src_upar, src_vort, & sources_external) call perf_print() call perf_reset() end subroutine end subroutine subroutine timestep_debug_output(comm_handler, tstep_count, tau, tinfo, rinfo) !! Writes debug output from timestep (e.g. convergence results of solvers) to file use params_brag_pardiss_model_m, only : heatflux_model, landau_numlorentzians_e, landau_numlorentzians_i type(comm_handler_t), intent(in) :: comm_handler !! Communicators integer, intent(in) :: tstep_count !! Number of timestep real(GP), intent(in) :: tau !! Time integer, intent(in), dimension(tinfo_size) :: tinfo !! Info status (integers) real(GP), intent(in), dimension(tinfo_size,2) :: rinfo !! Info status (float) integer :: iplane, nplanes, ierr, recvcount, il, ilr integer, allocatable, dimension(:,:) :: tinfo_gather real(GP), allocatable, dimension(:,:,:) :: rinfo_gather character(len=:),allocatable :: wrt_adv ! Gather tinfo and rinfo on plane 0 nplanes = comm_handler%get_nplanes() recvcount = nplanes*tinfo_size allocate(tinfo_gather(tinfo_size, nplanes)) allocate(rinfo_gather(tinfo_size, 2, nplanes)) call MPI_Gather(tinfo, tinfo_size, MPI_Integer, tinfo_gather, tinfo_size, MPI_Integer, 0, & comm_handler%get_comm_planes(), ierr) call MPI_Gather(rinfo, 2*tinfo_size, MPI_GP, rinfo_gather, 2*tinfo_size, MPI_GP, 0, & comm_handler%get_comm_planes(), ierr) ! Write debug output open(unit = 55, file = 'tstep_debug.out', status = 'unknown', position = 'append') if (is_rank_info_writer) then write(55,*)repeat('-',100) write(55,501)tstep_count, tau 501 FORMAT('tstep_count = ', I10,3X,'tau = ', ES15.7E3,X) write(55,*)'Statistics of 2D solves' write(55,502)'iplane', 'pot', 'pot_zon', 'ohm', 'ohm_pen', 'n_dens', 'n_parmom', 'n_press', 'n_smooth' 502 FORMAT(3X,A6,8('|',X,A14)) do iplane = 1, nplanes write(55,503)iplane-1, tinfo_gather(1,iplane), rinfo_gather(1,1,iplane), & tinfo_gather(2,iplane), rinfo_gather(2,1,iplane), & tinfo_gather(3,iplane), rinfo_gather(3,1,iplane), & tinfo_gather(4,iplane), rinfo_gather(4,1,iplane), & tinfo_gather(8,iplane), rinfo_gather(8,1,iplane), & tinfo_gather(9,iplane), rinfo_gather(9,1,iplane), & tinfo_gather(10,iplane), rinfo_gather(10,1,iplane), & tinfo_gather(11,iplane), rinfo_gather(11,1,iplane) 503 FORMAT(3X,I6,8('|',X,I4,X,ES9.2E2)) enddo write(55,*)'Statistics of 3D Braginskii solves' write(55,504)'niter', 'pseudores', 'true res' 504 FORMAT(10X,A5,X,2(A9,X)) if (heatflux_model /= 'LANDAU') then write(55,505)'te:', tinfo_gather(5,1), rinfo_gather(5,1,1), rinfo_gather(5,2:2,1) if (bnddescr_te_core == BND_BRAGTYPE_ZONAL_NEUMANN) then ! Preceded by 7 braginskii + 4 neutrals + 12 lorentzian_e + 12 lorentzian_i (=35) indices write(55,505)'te2:', tinfo_gather(36,1), rinfo_gather(36,1,1), rinfo_gather(36,2:2,1) endif write(55,505)'ti:', tinfo_gather(6,1), rinfo_gather(6,1,1), rinfo_gather(6,2:2,1) if (bnddescr_ti_core == BND_BRAGTYPE_ZONAL_NEUMANN) then ! Preceded by 7 braginskii + 4 neutrals + 12 lorentzian_e + 12 lorentzian_i + te2 (=36) indices write(55,505)'ti2:', tinfo_gather(37,1), rinfo_gather(37,1,1), rinfo_gather(37,2:2,1) endif endif write(55,505)'upar:', tinfo_gather(7,1), rinfo_gather(7,1,1), rinfo_gather(7,2:2,1) 505 FORMAT(3X,A6,X,I5,2(X,ES9.2E2)) if (heatflux_model == 'LANDAU') then write(55,*)'Statistics of Landau solves' write(55,'(3X,A6)', advance='no')'qe:' do il = 1, landau_numlorentzians_e if (il == landau_numlorentzians_e) then wrt_adv = 'yes' else wrt_adv = 'no' endif ilr = il + 11 ! Shift by previous 7 braginskii + 4 neutrals (=11) indices write(55,507, advance=wrt_adv)tinfo(ilr), rinfo(ilr,1), rinfo(ilr,2) enddo 507 FORMAT(X,I4,2(X,ES9.2E2)) write(55,'(3X,A6)', advance='no')'qi:' do il = 1, landau_numlorentzians_i if (il == landau_numlorentzians_i) then wrt_adv = 'yes' else wrt_adv = 'no' endif ilr = il + 23 ! Shift by previous 7 braginskii + 4 neutrals + 12 lorentzian_e (=23) indices write(55,507, advance=wrt_adv)tinfo(ilr), rinfo(ilr,1), rinfo(ilr,2) enddo endif write(55,*)repeat('-',100) endif close(55) end subroutine end module