gyrokinetics / gs2 / 2021218321

!> This module provides the external interface to gs2. It contains functions

!> to initialize gs2, functions to run gs2, functions to finalize gs2 and

!> functions to override/tweak gs2 parameters.

!>

!> The main entry point is [[run_gs2]], which takes care of initialising all the

!> subsystems, solving the equations (either the time advance or eigensolver),

!> and then finalising all the subsystems. This looks something like:

!>

!>

!>     ! An object holding the state of the simulationg

!>     type(gs2_program_state_type) :: state

!>     ! Initialisation of subsystems

!>     call initialize_gs2(state)

!>     call initialize_equations(state)

!>     ! Solve the equations

!>     if (state%do_eigsolve) then

!>       call run_eigensolver(state)

!>     else

!>       call evolve_equations(state, state%nstep)

!>     end if

!>     ! Clean up

!>     call finalize_equations(state)

!>     call finalize_gs2(state)

!>

!> You can manipulate the [[gs2_program_state]] before running gs2.

!>

!> It is very important to note that just because this interface is presented

!> in an object-oriented way, it does not, unfortunately, mean that the entire

!> program state, i.e. data arrays etc, is encapsulated in the gs2_program_state

!> object. In fact most data is stored globally as module public variables.

!> The gs2_program_state object is provided for convenience, as a way of

!> monitoring the execution state of gs2, and because it is hoped that in the

!> future gs2 will, in fact, be thread safe and have proper data encapsulation.

!> A particular consequence of this is that only one instance of the

!> gs2_program_state object should exist on each process, i.e. in a given

!> memory space.

module gs2_main

  use gs2_init, only: init_type, init_level_list

  use optimisation_configuration, only: optimisation_type

  use constants, only: run_name_size

  use exit_codes, only: exit_code, EXIT_NOT_REQUESTED

  use mp, only: mp_comm_null

  implicit none

  public :: run_gs2, finish_gs2, parse_command_line, gs2_program_state_type, reset_time_step

  public :: initialize_gs2, finalize_equations, finalize_gs2

  public :: initialize_equations, evolve_equations, run_eigensolver

  public :: prepare_optimisations_overrides, finalize_overrides, initialize_wall_clock_timer

  public :: prepare_initial_values_overrides, set_initval_overrides_to_current_vals

  !> Holds data pertaining to some of the main algorithm timers.

  type gs2_timers_type

    real :: init(2) = 0.

    real :: advance(2) = 0.

    real :: timestep(2) = 0.

    real :: finish(2) = 0.

    real :: total(2) = 0.

    real :: diagnostics(2)=0.

    real :: eigval(2)=0.

    real :: main_loop(2) = 0.

  end type gs2_timers_type

  !> The object which specifies and records the gs2 program state.

  !> Some attributes of the object, like mp_comm_external, are

  !> designed to be directly manipulated by the user, and some are

  !> designed to store program state information and be manipulated

  !> internally, like `init%level`.

  type gs2_program_state_type

     !> A type for keeping track of the current initialization level

     !> of gs2, as well as storing all the overrides. See

     !> gs2_init::init_type for more information.

     type(init_type) :: init

     !> Do not set manually. If fewer than the available number of

     !> processors are being used, this is true for the processors that

     !> are active and false for those that lie idle.

     logical :: included = .true.

     !> Derived type [[gs2_timers_type]] holding some of the

     !> timer data for this run.

     type(gs2_timers_type) :: timers

     !> The exit flag is set to true by any part of the main timestep

     !> loop that wants to cause the loop to exit

     logical :: exit = .false.

     !> Whether the run has converged to a stationary state

     logical :: converged = .false.

     !> Whether to run the eigenvalue solver or not. Set equal to the

     !> input value in intialize_equations. Typically only important

     !> for the standalone gs2 program

     logical :: do_eigsolve = .false.

     !> This parameter is set equal to run_parameters::nstep in

     !> initialize_equations and is the maximum number of timesteps

     !> that can be run without reinitalising the diagnostics.  Do not

     !> modify. We should make this private but it is used in testing

     !> currently so we don't.

     integer :: nstep

     !> Gets set to the final value of istep in evolve equations. Any

     !> future calls to evolve_equations will increment this

     !> further. A call to finalize_equations will set it back to

     !> 0. Note that setting this manually is NOT advised.

     integer :: istep_end = 0

     !> Verbosity at which we print out debug messages

     integer :: verb = 3

     !> Parameters pertaining to cases when gs2 is being used as a

     !> library. external_job_id is not to be confused with the parameter

     !> job in mp, which identifies the subjob if running in list mode

     !> or with nensembles > 1

     integer :: external_job_id = -1

     !> is_external_job should be set to true when GS2 is being used

     !> as a library. Perhaps this could just check external_job_id instead?

     logical :: is_external_job = .false.

     !> If true, print full timing breakdown.

     logical :: print_full_timers = .true.

     !> If true, print run time or full timing breakdown, depending on

     !> the value of print_full_timers

     logical :: print_times = .true.

     !> MPI communicator. This MUST be set.

     integer :: mp_comm = mp_comm_null

     !> This is set in initialize_gs2 to the number of procs actually used

     integer :: nproc_actual = -1

     !> If true then we end up passing [[gs2_program_state_type:run_name]]

     !> as the base file name.

     logical :: run_name_external = .false.

     !> The run name to pass to [[init_file_utils]] if `run_name_external` is

     !> set to true. GS2 then uses this in place of any input file name passed

     !> at the command line.

     character(run_name_size) :: run_name = 'gs'

     !> If true then skip the call to diagnostics calculations in

     !> evolve_equations (main loop and after main loop)

     logical :: skip_diagnostics = .false.

     !> If true then ignore any requests to change the time step in

     !> evolve_equations as a part of the time advance algorithm

     logical :: dont_change_timestep = .false.

     !> Whether this is a list mode run

     logical :: list = .false.

     !> The number of identical runs happening simultaneously (used

     !> for ensemble averaging). Cannot be used in conjunction with

     !> list mode.

     integer :: nensembles = 1

     !> Optimisation configuration

     type(optimisation_type) :: optim

     !> Reason for end of simulation

     type(exit_code) :: exit_reason = EXIT_NOT_REQUESTED

   contains

     procedure :: write_used_inputs_file => write_used_inputs_file_bound

  end type gs2_program_state_type

  private

contains

  !> Parse some basic command line arguments. Currently just 'version' and 'help'.

  !>

  !> This should be called before anything else, but especially before initialising MPI.

    use git_version_mod, only: get_git_version

    use build_config, only : formatted_build_config

    use ingen_mod, only : run_ingen

    use mp, only: proc0, init_mp, finish_mp

    use file_utils, only: init_file_utils, finish_file_utils

    integer :: arg_count, arg_n, arg_length

    logical :: list

    character(len=*), parameter :: nl = new_line('a')

    character(len=*), parameter :: usage = &

         "gs2 [--version|-v] [--help|-h] [--build-config] [--check-input] [input file]" // nl // nl // &

         "GS2: A fast, flexible physicist's toolkit for gyrokinetics" // nl // &

         "For more help, see the documentation at https://gyrokinetics.gitlab.io/gs2/" // nl // &

         "or create an issue https://bitbucket.org/gyrokinetics/gs2/issues?status=open" // nl // &

         nl // &

         "  -h, --help           Print this message" // nl // &

         "  -v, --version        Print the GS2 version" // nl // &

         "  --build-config       Print the current build configuration" // nl // &

         "  --check-input FILE   Check input FILE for errors" // nl // &

         "  --parse-input FILE   Check we can parse the full input FILE"

        ! Next argument should be the input file

        end if

        ! Next argument should be the input file

        end if

      end if

    end do

  !> Run whatever checks we can early in the initialisation process

  !> which might lead us to abort. Intended to help save compute resource

  !> by halting as early as we can.

    use config_collection, only: gs2_config_type

    use gs2_diagnostics, only: check_restart_file_writeable

    use gs2_save, only: set_restart_file

    use init_g, only: add_restart_dir_to_file, set_restart_dir

    use file_utils, only: run_name

    logical :: logic_dummy

    integer :: ind_slash

    ! Let us try populating a config collection from file so that

    ! we get an early abort if there is a problem with parsing the

    ! provided input file. It is important to note that the values

    ! stored here may not be correct due to inaccurate smart defaults

    ! at this stage.

    ! Check we can write restart files etc. if needed

    associate(init_g => configs%init_g_config, diag => configs%diagnostics_config)

      ! Use logic_dummy for extra_files as we can't change save_distfn now and it won't

      ! lead to an abort anyway so suppress warning message.

      call check_restart_file_writeable(diag%file_safety_check, diag%save_for_restart, &

    end associate

  !> Starts the global wall clock timer used by check time. This is

  !> useful if you have multiple copies of gs2 running but you don't

  !> want to start them all at the same time, but you want them all to

  !> respect avail_cpu_time

    use job_manage, only: init_checktime

  !> Initialise message passing, open the input file, split the

  !> communicator if running in list mode or if nensembles > 1,

  !> initialize the timers. After calling this function, gs2 reaches

  !> init\%level = basic. If it is desired to provide an external

  !> commuicator or set the input file name externally, these should

  !> be set before calling this function.

    use file_utils, only: init_file_utils, run_name, run_name_target

    use gs2_init, only: init_gs2_init

    use job_manage, only: time_message, job_fork

    use job_manage, only: init_checktime

    use mp, only: init_mp, broadcast, job, nproc, proc0, iproc, get_processor_name

    use mp, only: use_nproc, included, mp_comm, mp_comm_null

    use redistribute, only: using_measure_scatter

    use unit_tests, only: debug_message, set_job_id

    use runtime_tests, only: verbosity

    use standard_header, only: standard_header_type, get_default_header, set_default_header

    use optionals, only: get_option_with_default

#ifdef OPENMP

    use omp_lib, only: omp_get_num_threads

#endif

    implicit none

    type(gs2_program_state_type), intent(inout) :: state

    !> Header for files with build and run information

    type(standard_header_type), intent(in), optional :: header_in

    !> If true, don't print start-up messages/header

    logical, intent(in), optional :: quiet

    ! Actual value for optional `header` input

    ! Internal value for optional `quiet` input

    logical :: quiet_value

#ifdef OPENMP

    integer :: nthread

#endif

       error stop "ERROR: Called initialize_gs2 twice &

    end if

    !Ensure the job id reported in debug messages is set to something before we call debug message. Note

    !that this may get set to a different value later in the Initialisation.

    ! If optimisations say to use a subset of available processors

    ! create the communicator for this here.

       end if

    else

    end if

    end if

      ! Report number of processors being used

         end if

         ! Pluralise processor

#ifdef OPENMP

         !$OMP PARALLEL

         nthread = omp_get_num_threads()

         !$OMP END PARALLEL

         write(*, '(A, I0, A)',advance="no") " with ", nthread, " thread"

         ! Pluralise thread

         if (nthread > 1) write(*, '(a)', advance="no") "s"

#endif

         ! Make sure we still get a blank line

      end if

      ! Call init_file_utils, ie. initialize the run_name, checking

      !  whether we are doing a list of runs.

      ! Until the logic of init_file_utils is fixed we set trin_run

      ! when an external file name has been provided... this prevents

      ! it overriding the name from the command line.

      call init_file_utils (state%list, trin_run = state%run_name_external, &

    end if

       call debug_message(state%verb, &

    end if

    ! Pass run name to other procs

    ! Check for early aborts

    !Set using_measure_scatter to indicate we want to use in "gather/scatter" timings

    ! Initialize the gs2 initialization system

  !> Initialize all the modules which are used to evolve the

  !> equations. After calling this function, gs2 reaches `init%level = full`.

    use job_manage, only: time_message

    use run_parameters, only: nstep, do_eigsolve

    use unit_tests, only: debug_message

    implicit none

    type(gs2_program_state_type), intent(inout) :: state

    ! This triggers initializing of all the grids, all the physics parameters

    ! and all the modules which solve the equations and the diagnostics

    ! Set defaults. These are typically only important

    ! for the standalone gs2 program

  end subroutine initialize_equations

  !> Run the initial value solver. nstep_run must

  !> be less than or equal to `state%nstep`, which is

  !> set from the input file. The cumulative number of

  !> steps that can be run cannot exceed `state%nstep`,

  !> Examples:

  !>

  !>      call evolve_equations(state, state%nstep) ! Basic run

  !>

  !> Note that these two calls:

  !>

  !>      call evolve_equations(state, state%nstep/2)

  !>      call evolve_equations(state, state%nstep/2)

  !>

  !> will in general have the identical effect to the single call above.

  !> There are circumstances when this is not true,

  !> for example, if the first of the two calls to evolve_equations

  !> exits without running nstep/2 steps (perhaps because the

  !> growth rate has converged).

  !>

  !> This example will cause an error because the total number

  !> of steps exceeds `state%nstep`:

  !>

  !>      call evolve_equations(state, state%nstep/2)

  !>      call evolve_equations(state, state%nstep/2)

  !>      call evolve_equations(state, state%nstep/2)

    use collisions, only: vnmult

    use dist_fn_arrays, only: gnew

    use fields, only: advance

    use gs2_diagnostics, only: nsave, loop_diagnostics

    use gs2_diagnostics_new, only: run_diagnostics, gnostics

    use gs2_reinit, only: check_time_step

    use gs2_save, only: gs2_save_for_restart

    use gs2_time, only: user_time, update_time, write_dt, user_dt

    use gs2_time, only: code_dt, code_dt_prev1, code_dt_prev2, code_dt_max

    use job_manage, only: time_message, checkstop, checktime

    use mp, only: proc0, scope, subprocs, broadcast

    use run_parameters, only: reset, has_phi, has_apar, has_bpar, nstep, max_sim_time

    use run_parameters, only: avail_cpu_time, margin_cpu_time, use_old_diagnostics

    use run_parameters, only: progress_frequency, ncheck_stop, save_timer_statistics

    use unit_tests, only: debug_message

    use standard_header, only: standard_header_type, get_default_header

    use exit_codes, only: EXIT_MAX_SIM_TIME

    use file_utils, only: open_output_file, close_output_file

    type(gs2_program_state_type), intent(inout) :: state

    type(standard_header_type), optional, intent(in) :: header

    integer :: istep, progress_unit

    integer, intent(in) :: nstep_run

    logical :: temp_initval_override_store, should_save_restart, show_progress

    integer :: istep_loop_start, istep_loop_max, istep_offset, nsteps_left

    real, dimension(2) :: progress_time

    real :: time_per_step

    character(len = 256) :: progress_message

    ! Actual value for optional `header` input

    call debug_message(state%verb, &

    ! Make sure exit is false before entering timestep loop

    ! We run for nstep_run iterations, starting from whatever istep we got

    ! to in previous calls to this function. Note that calling

    ! finalize_equations resets state%istep_end

         ! Should we set an exit_reason here?

       end if

       ! Try to take a time step

       !If the advance call has triggered a reset then change the timestep.

       !Currently this is triggered if immediate_reset is true and the nonlinear

       !cfl limit is exceeded.

         ! Now we've reset set the flag back to .false.

         ! If we've changed the timestep and have not decided to exit then

         ! we need to call advance to actually complete this timestep.

       end if

       call debug_message(state%verb+1, &

       else

       end if

            has_phi, has_apar, has_bpar, code_dt, code_dt_prev1, &

         else 

         end if

       end if

         call debug_message(state%verb+1, &

         !Note this should only trigger a reset for timesteps too small

         !as timesteps too large are already handled when immediate_reset is

         !.true. If immediate_reset is .false. then this is where we handle

         !all checking of the time step.

         call debug_message(state%verb+1, &

         !If something has triggered a reset then reset here

            ! Switch off the reset flag now, otherwise it can still be active

            ! on the next iteration.

         end if

       end if

       call debug_message(state%verb+1, &

       ! Check if the next step would exceed the maximum requested simulation time,

       ! and if so flag that we wish to exit. Only check this if we're not already exiting.

       end if

          ! Have to split this from previous logic as we shouldn't call mod if

          ! progress_frequency is zero and we can't rely on short circuit logic

             write(progress_message, &

                  '("Done step ",I0," of ",I0," : Time/step ",G10.3," s ETA ",G10.3," s")') &

             ! Reset progress timer

          end if

       end if

       end if

    end do

    end if

    ! Try to ensure that we write the time dependent diagnostics at the end of the simulation

    ! even in cases where nstep does not have nwrite as a factor.

       ! Note when we have completed all our iterations we pass istep - 1 as we assume istep has been 

       ! incremented once from the last step actually taken. In cases where we explicitly exit the loop

       ! we need to pass istep instead.

       else

       end if

    end if

    ! We also restore the value of the override switch to what it

    ! was when this function was called.

    use gs2_reinit, only: prepare_reset_time_step, do_reset_time_step, time_reinit

    use gs2_init, only: init_type, init, init_level_list

    use run_parameters, only: reset

    use job_manage, only: time_message

    implicit none

    type(init_type), intent(inout) :: current_init

    integer, intent(in) :: istep

    logical, intent(in out) :: my_exit

    integer, intent (in), optional :: job_id

    logical :: reset_in

    !First disable the reset flag so we can call

    !routines needed in reinit

    end if

  end subroutine reset_time_step

  !> Employ the eigenvalue solver to seek eigenmodes using SLEPc

#ifdef WITH_EIG

#ifdef FIELDS_EIG

    use fields_eigval, only: init_eigval, finish_eigval, BasicSolve

#else

    use eigval, only: init_eigval, finish_eigval, BasicSolve

#endif

#endif 

    use job_manage, only: time_message

    use mp, only: mp_abort

    type(gs2_program_state_type), intent(inout) :: state

#ifdef WITH_EIG

   call time_message(.false.,state%timers%total,' Total')

   !Start timer

   call time_message(.false.,state%timers%eigval,' Eigensolver')

   !Initialise slepc and the default/input eigensolver parameters

   call init_eigval

   !Create a solver based on input paramters, use it to solve and

   !then destroy it.

   call BasicSolve

   !Tidy up

   call finish_eigval

   !Stop timer

   call time_message(.false.,state%timers%eigval,' Eigensolver')

#else

#endif

  end subroutine run_eigensolver

  !> FIXME : Add documentation  

    use job_manage, only: time_message

    use unit_tests, only: debug_message

    implicit none

    type(gs2_program_state_type), intent(inout) :: state

  end subroutine finalize_equations

  !> FIXME : Add documentation  

    use file_utils, only: finish_file_utils

    use gs2_init, only: finish_gs2_init

    use job_manage, only: time_message

    use mp, only: finish_mp, proc0, barrier, unsplit_all, included, mp_comm_null

    use unit_tests, only: debug_message

    use standard_header, only: date_iso8601

    use optionals, only: get_option_with_default

    use run_parameters, only: save_timer_statistics

    implicit none

    type(gs2_program_state_type), intent(inout) :: state

    !> If true, don't print start-up messages/header

    logical, intent(in), optional :: quiet

    ! Internal value for optional `quiet` input

    logical :: quiet_value

        write  (*,*) "ERROR: Called finalize_gs2 at the &

          & wrong init_level (perhaps you have called finalize_gs2 &

          & without calling initialize_gs2, or without calling &

      end if

    end if

            state%init%opt_ov%old_comm /= mp_comm_null) &

    end if

  !> FIXME : Add documentation  

    type(gs2_program_state_type), intent(inout) :: state

  end subroutine prepare_optimisations_overrides

  !> FIXME : Add documentation  

    use fields, only: force_maxwell_reinit

    use gs2_init, only: init_level_list

    use gs2_reinit, only: in_memory

    use gs2_layouts, only: g_lo

    use kt_grids, only: naky, ntheta0

    use theta_grid, only: ntgrid

    use dist_fn_arrays, only: gexp_3

    implicit none

    type(gs2_program_state_type), intent(inout) :: state

    ! Initialize to the level below so that overrides are triggered

    call state%init%initval_ov%initialise(ntgrid, ntheta0, naky, g_lo%llim_proc, &

         g_lo%ulim_alloc, force_maxwell_reinit=force_maxwell_reinit, &

  end subroutine prepare_initial_values_overrides

  !> FIXME : Add documentation  

    use dist_fn_arrays, only: gnew, gexp_1, gexp_2, gexp_3

    use gs2_save, only: gs2_save_for_restart

    use mp, only: proc0, broadcast, mp_abort

    use collisions, only: vnmult

    use file_utils, only: error_unit

    use run_parameters, only: has_phi, has_apar, has_bpar

    use fields, only:  force_maxwell_reinit

    use fields_arrays, only: phinew, aparnew, bparnew

    use gs2_time, only: user_time, code_dt, code_dt_prev1, code_dt_prev2, code_dt_max

    use overrides, only: initial_values_overrides_type

    use antenna, only: dump_ant_amp

    use array_utils, only: copy

    implicit none

    type(initial_values_overrides_type), intent(inout) :: initval_ov

      call mp_abort("Trying to set initial value overrides &

      & before they are initialized... have you called &

      ! Do not use the value from [old_iface_]state%init%initval_ov%force_maxwell_reinit = initval_ov%force_maxwell_reinit as

      ! follows:

      !     if (.not. initval_ov%force_maxwell_reinit) then

      ! as was done previously. This is because this causes problems in gs2_init so we use the value from

      ! fields::force_maxwell_reinit directly (see comments in gs2_init for details). To use here the value from

      ! initval_ov%force_maxwell_reinit could result in inconsistent values being used in different parts of the code. Therefore, we

      ! now use the value from fields::force_maxwell_reinit directly here too.  Now that the reference to

      ! initval_ov%force_maxwell_reinit here has also been removed, this component is no longer referenced anywhere in the codebase.

      ! Therefore, the force_maxwell_reinit component has been completely removed from the initial_values_overrides_type derived

      ! type to avoid confusion in the future.

      end if

    else ! if(.not.in_memory)then

      !Should really do this with in_memory=.true. as well but

      !not sure that we really need to as we never read in the dumped data.

      call gs2_save_for_restart (gnew, user_time, vnmult, &

    endif

  !> FIXME : Add documentation  

    type(gs2_program_state_type), intent(inout) :: state

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !! Private subroutines

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !> Reset timers data to 0

    use fft_work, only: time_fft

    use fields_arrays, only: time_field, time_field_mpi