Main function of KORC.
The main program contains the calls to the main functions and subroutines. Also, it contains the variables that control the behavior of the core of KORC and all other external/optional modules.
Subroutine initialize_communications in korc_hpc that initializes MPI and OpenMP communications.
Subroutine timing_KORC in korc_hpc that times the execution of any parallel sections of KORC.
Subroutine initialize_HDF5 in korc_HDF5 that initializes HDF5 library.
Subroutine initialize_korc_parameters in korc_initialize that initializes paths and KORC parameters through load_korc_params on MPI processes.
Subroutine initialize_fields in korc_fields that initializes parameters of the EM fields, either analytically or from an external HDF5 file. Reads in &analytical_fields_params and &externalPlasmaModel namelists from input file.
Subroutine initialize_profiles in korc_profiles that initializes parameters of the plasma profiles, either analytically or from an external HDF5 file. Reads in &plasmaProfiles namelist from input file. Only initialized if collisions (params%collisions==T) are present.
Subroutine initialize_particles in korc_initialize that initializes particle parameters from &lasma_species namelist, allocates arrays for individual particles, including location, velocity, local EM fields and plasma profiles, etc., and calls initial_energy_pitch_dist to assign particles' energy and pitch angle according to the chosen distribution.
Subroutine compute_charcs_plasma_params in korc_units calculates the characteristic plasma parameters params%cpp that are used for normalizations. Also finds the maximum non-relativistic and relativistic cyclotron frequencies to be used for setting the timstep for the time-evolution algorithms.
Subroutine initialize_collision_params in korc_collisions that initializes collision parameters for the SS (single-species) and MS (multiple-species) data types, reading in namefiles from the KORC input file. MS reads in namelist &CollisionParamsMultipleSpecies while SS reads in namelist &CollisionParamsSingleSpecies.
Subroutine define_time_step in korc_initialize either loads time-stepping parameters for a restart, or defines new parameters based on a maximum timestep set by the inverse of the relativistic cyclotron frequency.
Subroutine normalize_variables in korc_units normalizes variables consistent with characteristic plasma parameters calculated in compute_charcs_plasma_params.
Subroutine normalize_collisions_params in korc_collisions that normalizes collision parameters for the SS (single-species) and MS (multiple-species) data types.
Subroutine define_collisions_time_step in korc_collisions that sets subcycling iteration number for collisions based off of the collision frequency model used.
Subroutine initialize_fields_interpolant in korc_interp calls EZspline subroutines EZspline_init for memory allocation and boundary condition setup and EZspline_setup to compute the necessary cubic coefficients needed for subsequent field interpolations. The magnetic field can be defined in terms of an axisymmetric scalar flux function, axisymmetric field, or 3D field, while the electric field can be defined as an axisymmetric or 3D field.
Subroutine initialize_profiles_interpolant in korc_interp calls EZspline subroutines EZlinear_init for axisymmetric (flux-surface quantities) or EZspline_init for 3D profiles for memory allocation and boundary condition setup and EZspline_setup to compute the necessary cubic coefficients needed for subsequent field interpolations. Only initialized if collisions (params%collisions==T) are present for ne, Te, Zeff
Subroutine set_up_particles_ic in korc_initialize calls subroutines to prescribe initial conditions or load them from file for a restart. Initial spatial values are prescribed with intitial_spatial_distribution in korc_spatial_distribution and initial velocity values are prescribed with initial_gyro_distribution in korc_velocity_distribution.
Subroutines save_simulation_parameters in korc_HDF5 and save_collision_params in korc_collisions call subroutines to save simulation and collision parameters.
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS) | :: | params | Contains the parameters that control the core of KORC: time steping, output list, etc. |
|||
| type(SPECIES), | DIMENSION(:), ALLOCATABLE | :: | spp | Contains the initial parameters of each species, which can be different electrons with different distribution functions. |
||
| type(FIELDS) | :: | F | and electric fields, or in the case of using external fields it contains the data used in the interpolations. See korc_fields for details. |
|||
| type(PROFILES) | :: | P | or in the case of using external fields it contains the data used in the interpolations. See korc_profiles for details. |
|||
| integer(kind=ip) | :: | it | Time iteration |
|||
| integer | :: | mpierr |
program main
!! @note Main function of KORC. @endnote
!! The main program contains the calls to the main functions and subroutines.
!! Also, it contains the variables that control
!! the behavior of the core of KORC and all other external/optional modules.
use korc_types
use korc_units
use korc_hpc
use korc_HDF5
use korc_fields
use korc_ppusher
use korc_interp
use korc_collisions
use korc_initialize
use korc_finalize
use korc_profiles
use korc_synthetic_camera
use korc_binning_diagnostic
implicit none
TYPE(KORC_PARAMS) :: params
!! Contains the parameters that control the core of KORC:
!! time steping, output list, etc.
TYPE(SPECIES), DIMENSION(:), ALLOCATABLE :: spp
!! Contains the initial parameters of each species, which
!! can be different electrons with different
!! distribution functions.
TYPE(FIELDS) :: F
!! F: Contains the parameters of the analytical magnetic
!! and electric fields, or in the case of using
!! external fields it contains the data used in the interpolations.
!!See [[korc_fields(module)]] for details.
TYPE(PROFILES) :: P
!! P: Contains the parameters of the analytical plasma profiles,
!! or in the case of using external
!! fields it contains the data used in the interpolations.
!! See [[korc_profiles(module)]] for details.
INTEGER(ip) :: it
!! Time iteration
INTEGER :: mpierr
call initialize_communications(params)
!!<h2>Order of KORC operations</h2>
!!
!!<h3>Communication and Timing</h3>
!! <h4>1\. Parallel Communications</h4>
!!
!! Subroutine [[initialize_communications]] in [[korc_hpc]] that
!! initializes MPI and OpenMP communications.
call timing_KORC(params)
!! <h4>2\. Timers</h4>
!!
!! Subroutine [[timing_KORC]] in [[korc_hpc]] that times the
!! execution of any parallel sections of KORC.
! * * * INITIALIZATION STAGE * * *!
call initialize_HDF5()
!!<h3>Initialization</h3>
!!
!! <h4>1\. HDF5</h4>
!!
!! Subroutine [[initialize_HDF5]] in [[korc_HDF5]] that initializes
!! HDF5 library.
call initialize_korc_parameters(params)
!! <h4>2\. Initialize korc parameters</h4>
!!
!! Subroutine [[initialize_korc_parameters]] in [[korc_initialize]] that
!! initializes paths and KORC parameters through [[load_korc_params]]
!! on MPI processes.
call initialize_fields(params,F)
!! <h4>3\. Initialize fields</h4>
!!
!! Subroutine [[initialize_fields]] in [[korc_fields]] that initializes
!! parameters of the EM fields, either analytically or from an external HDF5
!! file. Reads in &analytical_fields_params and
!! &externalPlasmaModel namelists from input file.
call initialize_profiles(params,P,F)
!! <h4>4\. Initialize Profiles</h4>
!!
!! Subroutine [[initialize_profiles]] in [[korc_profiles]] that initializes
!! parameters of the plasma profiles, either analytically or from an
!! external HDF5
!! file. Reads in &plasmaProfiles namelist from input file.
!! Only initialized if collisions (params%collisions==T) are present.
call initialize_particles(params,F,P,spp) ! Initialize particles
!! <h4>5\. Initialize Particle Velocity Phase Space</h4>
!!
!! Subroutine [[initialize_particles]] in [[korc_initialize]] that
!! initializes particle parameters from &lasma_species namelist,
!! allocates arrays for individual particles, including location, velocity,
!! local EM fields and plasma profiles, etc., and
!! calls [[initial_energy_pitch_dist]] to assign particles' energy and pitch
!! angle according to the chosen distribution.
! write(6,'("init eta: ",E17.10)') spp(1)%vars%eta
call initialize_synthetic_camera(params,F)
!! <h4>7\. Initialize Synthetic Cameras</h4>
call initialize_binning_diagnostic(params)
!! <h4>8\. Initialize Binning Diagnostic</h4>
call compute_charcs_plasma_params(params,spp,F)
!! <h4>9\. Compute Characteristic Plasma Parameters</h4>
!!
!! Subroutine [[compute_charcs_plasma_params]] in [[korc_units]] calculates
!! the characteristic plasma parameters params%cpp that are used for normalizations.
!! Also finds the maximum non-relativistic and relativistic cyclotron frequencies
!! to be used for setting the timstep for the time-evolution algorithms.
call initialize_collision_params(params)
!! <h4>6\. Initialize Collision Parameters</h4>
!!
!! Subroutine [[initialize_collision_params]] in [[korc_collisions]] that
!! initializes collision parameters for the SS (single-species) and MS
!! (multiple-species) data types, reading in namefiles from the KORC input file.
!! MS reads in namelist &CollisionParamsMultipleSpecies while SS reads in
!! namelist &CollisionParamsSingleSpecies.
call define_time_step(params)
!! <h4>10\. Define Time Step</h4>
!!
!! Subroutine [[define_time_step]] in [[korc_initialize]] either loads
!! time-stepping parameters for a restart, or defines new parameters based
!! on a maximum timestep
!! set by the inverse of the relativistic cyclotron frequency.
call initialize_particle_pusher(params)
!! <h4>11\. Initialize Particle Pusher</h4>
if (params%SC_E) then
call define_SC_time_step(params,F)
end if
call normalize_variables(params,spp,F,P)
!! <h4>12\. Normalize Variables</h4>
!!
!! Subroutine [[normalize_variables]] in [[korc_units]] normalizes
!! variables consistent with characteristic plasma parameters
!! calculated in [[compute_charcs_plasma_params]].
call normalize_collisions_params(params)
!! <h4>13\. Normalize Collision Parameters </h4>
!!
!! Subroutine [[normalize_collisions_params]] in [[korc_collisions]] that
!! normalizes collision parameters for the SS (single-species) and MS
!! (multiple-species) data types.
call define_collisions_time_step(params)
!! <h4>14\. Define Collision Time Step</h4>
!!
!! Subroutine [[define_collisions_time_step]] in [[korc_collisions]] that
!! sets subcycling iteration number for collisions based off of the collision
!! frequency model used.
! *** *** *** *** *** *** *** *** *** *** *** *** ***
! *** BEYOND THIS POINT VARIABLES ARE DIMENSIONLESS ***
! *** *** *** *** *** *** *** *** *** *** *** *** ***
call initialize_fields_interpolant(params,F)
!! <h4>15\. Initialize Fields Interpolant</h4>
!!
!! Subroutine [[initialize_fields_interpolant]] in [[korc_interp]] calls
!! EZspline
!! subroutines EZspline_init for memory allocation and boundary condition
!! setup
!! and EZspline_setup to compute the necessary cubic coefficients needed
!! for subsequent
!! field interpolations. The magnetic field can be defined in terms of an
!! axisymmetric
!! scalar flux function, axisymmetric field, or 3D field, while the
!! electric field
!! can be defined as an axisymmetric or 3D field.
call initialize_profiles_interpolant(params,P)
!! <h4>16\. Initialize Profiles Interpolant</h4>
!!
!! Subroutine [[initialize_profiles_interpolant]] in [[korc_interp]]
!! calls EZspline
!! subroutines EZlinear_init for axisymmetric (flux-surface quantities) or
!! EZspline_init for 3D profiles for memory allocation and boundary
!! condition setup
!! and EZspline_setup to compute the necessary cubic coefficients needed
!! for subsequent
!! field interpolations.
!! Only initialized if collisions (params%collisions==T) are present for
!! ne, Te, Zeff
if (params%mpi_params%rank .EQ. 0) then
write(6,'("* * * * INITIALIZING INITIAL CONDITIONS * * * *")')
end if
call set_up_particles_ic(params,F,spp,P)
if (params%mpi_params%rank .EQ. 0) then
write(6,'("* * * * * * * * * * * * * * * * * * * * * * * *",/)')
end if
! write(6,'("post ic eta: ",E17.10)') spp(1)%vars%eta
!! <h4>17\. Set Particle Initial Conditions</h4>
!!
!! Subroutine [[set_up_particles_ic]] in [[korc_initialize]] calls
!! subroutines to prescribe initial conditions or load them
!! from file for a restart. Initial spatial values are prescribed with
!! [[intitial_spatial_distribution]] in [[korc_spatial_distribution]] and
!! initial velocity values are prescribed with [[initial_gyro_distribution]]
!! in [[korc_velocity_distribution]].
! if (minval(spp(1)%vars%Y(:,1)).lt.1._rp/params%cpp%length) stop 'error with init'
! * * * INITIALIZATION STAGE * * *
! * * * SAVING INITIAL CONDITION AND VARIOUS SIMULATION PARAMETERS * * * !
call save_simulation_parameters(params,spp,F,P)
call save_collision_params(params)
!! <h4>18\. Save Simulation and Collision Parameters</h4>
!!
!! Subroutines [[save_simulation_parameters]] in [[korc_HDF5]] and
!! [[save_collision_params]] in [[korc_collisions]] call
!! subroutines to save simulation and collision parameters.
! write(6,'("GC init eta: ",E17.10)') spp(1)%vars%eta
if (.NOT.(params%restart.OR.params%proceed.or.params%reinit)) then
if (params%orbit_model(1:2).eq.'FO') then
call FO_init(params,F,spp,.true.,.false.)
else if (params%orbit_model(1:2).eq.'GC') then
call GC_init(params,F,spp)
end if
if (params%SC_E) then
if (params%field_model(1:1).eq.'A') then
call init_SC_E1D(params,F,spp(1))
else if (params%field_model(1:1).eq.'E') then
call init_SC_E1D_FS(params,F,spp(1))
end if
end if
else
call get_fields(params,spp(1)%vars,F)
if (params%SC_E) then
if (params%field_model(1:1).eq.'A') then
call reinit_SC_E1D(params,F)
else if (params%field_model(1:1).eq.'E') then
call reinit_SC_E1D_FS(params,F)
end if
end if
end if
if (.NOT.(params%restart.OR.params%proceed)) then
call save_simulation_outputs(params,spp,F) ! Save initial condition
call synthetic_camera(params,spp)
call binning_diagnostic(params,spp)
end if
! * * * SAVING INITIAL CONDITION AND VARIOUS SIMULATION PARAMETERS * * * !
! write(6,'("pre ppusher loop eta: ",E17.10)') spp(1)%vars%eta
call timing_KORC(params)
if (params%orbit_model(1:2).eq.'FO'.and.params%field_eval.eq.'eqn') then
do it=params%ito,params%t_steps,params%t_skip
call adv_FOeqn_top(params,F,P,spp)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'FO'.and.params%field_eval.eq.'interp') then
call FO_init(params,F,spp,.false.,.true.)
! Initial half-time particle push
do it=params%ito,params%t_steps,params%t_skip
call adv_FOinterp_top(params,F,P,spp)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
! if (params%orbit_model.eq.'FO'.and.params%field_eval.eq.'interp') then
! do it=params%ito,params%t_steps
! params%time = params%init_time+REAL(it,rp)*params%dt
! params%it = it
! if ( modulo(it,params%output_cadence) .EQ. 0_ip ) then
! call advance_FOinterp_vars(params,spp,params%dt, &
! .TRUE.,.FALSE.)
! call save_simulation_outputs(params,spp)
! call synthetic_camera(params,spp) ! Synthetic camera
! call binning_diagnostic(params,spp) ! Binning diagnostic
! call save_restart_variables(params,spp)
! else
! call advance_FOinterp_vars(params,spp,params%dt, &
! .FALSE.,.FALSE.)
! end if
! end do
! end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'eqn') then
do it=params%ito,params%t_steps,params%t_skip*params%t_it_SC
call adv_GCeqn_top(params,F,P,spp)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip*params%t_it_SC,rp)*params%dt
params%it = it-1_ip+params%t_skip*params%t_it_SC
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.params%field_model(10:12).eq.'PSI'.and. &
params%SC_E) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_psi_top_FS(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.(params%field_model(10:12).eq.'PSI'.OR. &
params%field_model(12:14).eq.'PSI').and. &
(.not.params%SC_E).and.(.not.F%Dim2x1t)) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_psi_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.(params%field_model(10:12).eq.'PSI'.OR. &
params%field_model(12:14).eq.'PSI').and. &
(.not.params%SC_E).and.F%Dim2x1t.and.(.not.F%ReInterp_2x1t)) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_psi2x1t_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.(params%field_model(10:12).eq.'PSI'.OR. &
params%field_model(12:14).eq.'PSI').and. &
(.not.params%SC_E).and.F%Dim2x1t.and.F%ReInterp_2x1t) then
if (params%mpi_params%rank .EQ. 0) then
write(6,*) 'time',F%X%PHI(F%ind_2x1t)*params%cpp%time
end if
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_psiwE_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
F%ind_2x1t=F%ind_2x1t+1_ip
if (params%mpi_params%rank .EQ. 0) then
write(6,*) 'time',F%X%PHI(F%ind_2x1t)*params%cpp%time
end if
call initialize_fields_interpolant(params,F)
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.F%dBfield) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_2DBdB_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
F%axisymmetric_fields.and.(params%field_model(10:12).eq.'2DB'.or. &
params%field_model(12:13).eq.'2D').and..not.(F%dBfield)) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_B2D_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
.not.(F%axisymmetric_fields).and.(F%dBfield).and. &
(params%field_model(10:14).eq.'3DBdB')) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_3DBdB_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
.not.(F%axisymmetric_fields).and.(F%dBfield).and. &
.not.(params%field_model(10:14).eq.'3DBdB')) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_3DBdB1_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
if (params%orbit_model(1:2).eq.'GC'.and.params%field_eval.eq.'interp'.and. &
.not.(F%axisymmetric_fields).and..not.(F%dBfield)) then
do it=params%ito,params%t_steps,params%t_skip
call adv_GCinterp_B_top(params,spp,P,F)
params%time = params%init_time &
+REAL(it-1_ip+params%t_skip,rp)*params%dt
params%it = it-1_ip+params%t_skip
call save_simulation_outputs(params,spp,F)
call synthetic_camera(params,spp) ! Synthetic camera
call binning_diagnostic(params,spp) ! Binning diagnostic
call save_restart_variables(params,spp,F)
end do
end if
call timing_KORC(params)
! * * * FINALIZING SIMULATION * * *
call finalize_HDF5()
call finalize_interpolants(params)
! DEALLOCATION OF VARIABLES
call deallocate_variables(params,F,spp)
call deallocate_collisions_params(params)
call finalize_communications(params)
! * * * FINALIZING SIMULATION * * *
if (params%mpi_params%rank .EQ. 0) then
write(6,'("KORC ran successfully!")')
end if
end program main