From 3ee0240c2a76119a6dab429466f6f44dda1d8485 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Mon, 10 Aug 2026 17:02:07 +0200 Subject: [PATCH 01/16] Refactored particle initialisation in kinetic_main, rejection functions have their own module and bundling type, this commit compiles --- particles/initialisers/initialisers_RE.f90 | 85 ++---- particles/initialisers/initialisers_base.f90 | 63 +++-- .../mod_import_experimental_dist.f90 | 6 +- .../initialisers/mod_initialise_particles.f90 | 254 ++++++++++-------- particles/initialisers/mod_rej_f.f90 | 217 +++++++++++++++ 5 files changed, 419 insertions(+), 206 deletions(-) rename particles/{ => initialisers}/mod_import_experimental_dist.f90 (99%) create mode 100644 particles/initialisers/mod_rej_f.f90 diff --git a/particles/initialisers/initialisers_RE.f90 b/particles/initialisers/initialisers_RE.f90 index a69e15da3e..e93ad59cc7 100644 --- a/particles/initialisers/initialisers_RE.f90 +++ b/particles/initialisers/initialisers_RE.f90 @@ -5,66 +5,17 @@ module initialisers_RE use mod_particle_sim use mod_rng use initialisers_base - use constants, only: EL_CHG, ATOMIC_MASS_UNIT, SPEED_OF_LIGHT, MASS_ELECTRON, TWOPI + use constants, only: EL_CHG, ATOMIC_MASS_UNIT, SPEED_OF_LIGHT, MASS_ELECTRON, TWOPI, MU_ZERO use mod_pusher_tools, only: get_orthonormals use mod_coordinate_transforms, only: vector_cylindrical_to_cartesian use mod_model_settings - use phys_module, only: CENTRAL_DENSITY, CENTRAL_MASS, ATOMIC_MASS_UNIT, MU_ZERO + use phys_module, only: CENTRAL_DENSITY, CENTRAL_MASS use equil_info + use mod_rej_f implicit none contains - -! Quick and rough function to sample markers based on RZ-coordinates -pure function RZ_pdf(var) result(p) - real*8, intent(in) :: var(2) ! var(1)=j - real*8 :: p - real*8 :: minor_r - real*8 :: R_ax, Z_ax - - R_ax = ES%R_axis - Z_ax = ES%Z_axis - - minor_r = sqrt((var(1)-Z_ax)**2 + (var(2)-R_ax)**2) - - p = 1 / (1 + minor_r)**2 -end function RZ_pdf - -pure function analytical_pdf(var) result(p) - real*8, intent(in) :: var(2) ! var(1)=j - real*8 :: p - real*8 :: minor_r - real*8 :: nu - real*8 :: R_ax, Z_ax - real*8 :: LCFS_a - - R_ax = ES%R_axis - Z_ax = ES%Z_axis - LCFS_a = ES%LCFS_a - nu = 2.d0 - - minor_r = sqrt((var(1)-Z_ax)**2 + (var(2)-R_ax)**2) - - p = (1.d0 - (minor_r/LCFS_a)**2)**nu - -end function analytical_pdf - -! Quick and rough function to sample markers proportionally to toroidal current density -pure function current_pdf(var) result(p) - real*8, intent(in) :: var(2) ! var(2)=j, var(1) = R - !real*8, intent(in) :: var(1) ! var(1)=j - real*8 :: p - real*8 :: jzmin, jzmax - - !> temporarily hard coded, but should be able to be obtained from fluid restart file - jzmax = 3.0 / 10.0 !1.173 / 10 - jzmin = 0.0001239 / 11.0 !0.0003166 / 11 - - p = (var(2)/var(1)-jzmin)/(jzmax-jzmin) - -end function current_pdf - -subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, std_energy) +subroutine initialise_gaussian_re(sim, group_num, rng, space_pdf, energy, pitch, std_energy) use phys_module, only: tstep_particles use mod_kinetic_relativistic use mod_sampling, only: boxmueller_transform @@ -72,7 +23,7 @@ subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, st type(particle_sim), intent(inout) :: sim integer, intent(in) :: group_num class(type_rng), intent(in) :: rng - character(len=50), intent(in) :: init_pdf + type(spatial_pdf), optional, intent(in) :: space_pdf real*8, intent(in) :: energy, pitch ! Kinetic energy in units of eV and pitch real*8, optional, intent(in) :: std_energy real*8, allocatable :: p_tot(:), p_par(:), p_perp(:) @@ -83,20 +34,14 @@ subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, st integer :: num_part real*8, allocatable :: ran_uniform(:), ran_gaussian(:) - select case (trim(init_pdf)) - case ("RZ") - call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-2,-1], transform=RZ_pdf) - case ("current") - call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-1,var_zj], transform=current_pdf) - case ("analytical") - call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-2,-1], transform=analytical_pdf) - case default - if (sim%my_id == 0) then - write(*,*) "ERROR: ", trim(init_pdf), " is not a valid pdf/transform function for " - write(*,*) " for group '", sim%groups(group_num)%id, "' when using the 'basic_initialization' function" - endif - stop 1 - end select + if (present(space_pdf)) then + call initialise_particles(sim%groups(group_num)%particles, & + sim%fields%node_list, sim%fields%element_list, rng, & + variables=space_pdf%vars, transform=space_pdf%f) + else + call initialise_particles(sim%groups(group_num)%particles, & + sim%fields%node_list, sim%fields%element_list, rng) + endif num_part = size(sim%groups(group_num)%particles,1) @@ -156,6 +101,6 @@ subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, st deallocate(p_par) deallocate(p_perp) -end subroutine basic_initialization +end subroutine initialise_gaussian_re -end module initialisers_RE \ No newline at end of file +end module initialisers_RE diff --git a/particles/initialisers/initialisers_base.f90 b/particles/initialisers/initialisers_base.f90 index 21e667d17b..10fcc5a42c 100644 --- a/particles/initialisers/initialisers_base.f90 +++ b/particles/initialisers/initialisers_base.f90 @@ -6,6 +6,7 @@ module initialisers_base use mod_particle_types use constants use mod_interp + use mod_rej_f implicit none private public initialise_particles, no_transform, adjust_particle_weights @@ -26,13 +27,6 @@ subroutine find_RZ(node_list,element_list,R_find,Z_find,R_out,Z_out,ielm_out,s_o integer, intent(inout) :: ielm_out integer, intent(out) :: ifail end subroutine find_RZ - function rej_f(n, P, gradP) - implicit none - integer, intent(in) :: n - real*8, dimension(n), intent(in) :: P - real*8, dimension(3,n), intent(in) :: gradP - real*4 :: rej_f - end function rej_f function real_f(n_x,x,st,time,i_elm,fields,x_min,x_max,& n_real_param,real_param,n_int_param,int_param) use mod_fields, only: fields_base @@ -101,18 +95,19 @@ subroutine initialise_particles(particles, node_list, element_list, & class(particle_base), dimension(:), intent(inout) :: particles type(type_node_list), intent(in) :: node_list type(type_element_list), intent(in) :: element_list - class(type_rng), intent(in) :: rng !< What type of random number generator to use. Is re-seeded in the subroutine. + class(type_rng), intent(in) :: rng !< What type of random number generator to use. Is re-seeded in the subroutine. integer, dimension(:), intent(in), optional :: variables !< Which variables from JOREK to use. If absent, sample uniformly. - real*8, external, optional :: transform !< Merge variables into a single criterium between 0 and 1 for rej. sampling - !< Special values: 0 = 1, -1 = R, -2 = Z, -3 = Phi. Must be in ascending order! - real*8, intent(in), optional :: f !< Weighting factor: f=0 indicates uniform weights, f=1 indicates uniform distribution - !< (particle weight proportional to transform(P) at that point.) If omitted take f=0. - real*8, dimension(2), intent(in), optional :: Rbound, Zbound, Phibound !< Between which coordinates to sample (RZPhi). - !< if omitted, determine automatically from node_list + procedure(rej_f), optional :: transform !< Merge variables into a single criterium between 0 and 1 for rej. sampling + !< Special values: 0 = 1, -1 = R, -2 = Z, -3 = Phi. Must be in ascending order! + real*8, intent(in), optional :: f !< Weighting factor: f=0 indicates uniform weights, f=1 indicates uniform distribution + !< (particle weight proportional to transform(P) at that point.) If omitted take f=0. + real*8, dimension(2), intent(in), optional :: Rbound, Zbound, Phibound !< Between which coordinates to sample (RZPhi). + !< if omitted, determine automatically from node_list logical, intent(in), optional :: rng_n_streams_round_off_in !< round-off the rng n_streams at 2**ceil ! Internal variables real*8 :: R, Z, phi, s, t, DUMMY_REAL + real*8 :: R_i, R_s, R_t, Z_i, Z_s, Z_t, xjac real*8 :: Rbox(2), Zbox(2), Phibox(2) integer :: i, j, k, ifail real*8 :: ran(7) @@ -123,10 +118,11 @@ subroutine initialise_particles(particles, node_list, element_list, & integer :: my_id, n_mpi integer :: seed logical :: rng_n_streams_round_off - real*8, dimension(:), allocatable :: P + real*8, dimension(:), allocatable :: P + real*8, dimension(:,:), allocatable :: gradP class(type_rng), allocatable, dimension(:) :: rngs ! The RNGs for all the threads - integer, dimension(:), allocatable :: i_to_find - logical, dimension(:), allocatable :: not_found + integer, dimension(:), allocatable :: i_to_find + logical, dimension(:), allocatable :: not_found ostart = 0.d0 oend = 0.d0 @@ -142,6 +138,8 @@ subroutine initialise_particles(particles, node_list, element_list, & end if ! Get the number of mhd variables to use allocate(P(size(variables,1))) + allocate(gradP(3, size(variables,1))) + gradP = 0.d0 n_mhd = count(variables .gt. 0) n_geom = size(variables, 1) - n_mhd else @@ -217,7 +215,8 @@ subroutine initialise_particles(particles, node_list, element_list, & #endif !$omp shared(particles, node_list, element_list, Rbox, Zbox, PhiBox, variables, & !$omp rngs, n_threads, n_streams, seed, my_id, n_mhd, n_geom, i_to_find, not_found) & - !$omp private(j, i, R, Z, phi, i_elm, s, t, ifail, seq, ran, i_thread, P, DUMMY_REAL) + !$omp private(j, i, R, Z, phi, i_elm, s, t, ifail, seq, ran, i_thread, P, DUMMY_REAL, gradP, & + !$omp R_i, R_s, R_t, Z_i, Z_s, Z_t, xjac) i_thread = 0 !$ i_thread=omp_get_thread_num() !$omp do schedule(static) @@ -232,19 +231,30 @@ subroutine initialise_particles(particles, node_list, element_list, & if (present(variables)) then ! Select the mhd variables requested if (n_mhd .ge. 1) then - call interp_0(node_list,element_list,i_elm,variables(n_geom+1:n_geom+n_mhd),n_mhd,s,t,phi,P(n_geom+1:n_geom+n_mhd)) + call interp_PRZ(node_list, element_list, i_elm, variables(n_geom+1:n_geom+n_mhd), & + n_mhd, s, t, phi, P(n_geom+1:n_geom+n_mhd), gradP(1, n_geom+1:n_geom+n_mhd), & + gradP(2, n_geom+1:n_geom+n_mhd), gradP(3, n_geom+1:n_geom+n_mhd), R_i, R_s, R_t, & + Z_i, Z_s, Z_t) + + !> Calculate R/Z derivs for gradP if needed by the rej_f (interp_PRZ gives s, t derivs) + xjac = R_s*Z_t - R_t*Z_s + do k = 1, n_mhd + gradP(1:2,n_geom+k) = [Z_t*gradP(1,n_geom+k) - Z_s*gradP(2,n_geom+k), & + R_s*gradP(2,n_geom+k) - R_t*gradP(1,n_geom+k)]/xjac + enddo end if + do k=1,n_geom select case (variables(k)) - case (0); P(k) = 1.d0 - case (-1); P(k) = R - case (-2); P(k) = Z - case (-3); P(k) = phi + case (0); P(k) = 1.d0 ; gradP(:,k) = 0.d0 ! 0 + case (-1); P(k) = R ; gradP(:,k) = [1.d0, 0.d0, 0.d0] ! R + case (-2); P(k) = Z ; gradP(:,k) = [0.d0, 1.d0, 0.d0] ! Z + case (-3); P(k) = phi ; gradP(:,k) = [0.d0, 0.d0, 1.d0] ! phi end select end do if (present(transform)) then - if (ran(4) .lt. transform(p)) then + if (ran(4) .lt. transform(size(p,1), p, gradP)) then particles(j)%x = [R, Z, phi] particles(j)%i_elm = i_elm particles(j)%st = [s, t] @@ -275,6 +285,11 @@ subroutine initialise_particles(particles, node_list, element_list, & not_found = .true. end do + if (present(variables)) then + deallocate(P) + deallocate(gradP) + endif + call cpu_time(t1) !$ oend = omp_get_wtime() write(*,'(i5,A,2f12.4)') my_id, ' Time particle initialize cpu/wall :',t1-t0, oend-ostart diff --git a/particles/mod_import_experimental_dist.f90 b/particles/initialisers/mod_import_experimental_dist.f90 similarity index 99% rename from particles/mod_import_experimental_dist.f90 rename to particles/initialisers/mod_import_experimental_dist.f90 index dcfc42decd..a5fd2b36e4 100644 --- a/particles/mod_import_experimental_dist.f90 +++ b/particles/initialisers/mod_import_experimental_dist.f90 @@ -1,6 +1,10 @@ module mod_import_experimental_dist use mod_rng - use particle_tracer + use mod_particle_sim + use mod_fields + use mod_boris + use mod_pcg32_rng + use initialisers_base use data_structure use mod_particle_types use constants diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index a2d5759f90..9d93f91041 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -3,16 +3,41 @@ !> This module mainly handles the interface between the initialisers and the simulation !> The function of the specific initialiser subroutines are contained in the !> "initialisers_*.f90" files +!> +!> init_function -> 1st dispatch, choses global initialiser. +!> Could be only phase-space (e.g. maxwell) or both real and phase-space (e.g. experimental) +!> init_pdf -> 2nd dispatch. Ignored if init_function includes real-space init (e.g. experimental) +!> Optional real-space distribution function. Implemented via rejection sampling (see mod_rej_f). +!> particle_type -> Determines valid init_functions +!> +!> Valid init_functions: +!> 'maxwell' -> Maxwellian Energy distribution, uniform pitch. Spatial dist set by init_pdf (or uniform if no init_pdf specified) +!> Temperature specified by T_maxwell in input file. Only valid for particle_kinetic_leapfrog +!> 'gaussian_re' -> Monoenergetic or Gaussian energy distribution, delta function in pitch. +!> Spatial dist set by init_pdf (or uniform if no init_pdf specified) +!> re_energy, re_std_energy, re_pitch specified in input file. Only valid for particle_kinetic_relativistic +!> 'experimental' -> Samples from external F(R, Z, energy, pitch) distribution, expects 'experimental_dist.h5' hdf5 file with the correct format +!> real-space dist is implicit in file, init_pdf is ignored. Only valid for particle_kinetic_leapfrog +!> +!> Any other init_function produces an explicit error at initialisation time. module mod_initialise_particles + use mod_rej_f use initialisers_RE use initialisers_base + use mod_import_experimental_dist use phys_module, only: part_group_configs, type_part_group_config, n_part_groups use mod_particle_group_id, only: matching_part_config_indices - use equil_info - + use mpi + + implicit none + integer :: ierr !> mpi error code + contains + !> Top-level initialisation loop over all particle groups + !> atm skips ics/ncs groups (assume particles born via other methods) + !> delegates everything else to initialise_group subroutine initialise_particles_for_sim(sim) implicit none class(particle_sim), intent(inout) :: sim @@ -20,139 +45,146 @@ subroutine initialise_particles_for_sim(sim) do i=1, n_part_groups ! loop over part_groups_in_use j = matching_part_config_indices(i) ! get the matching part_group_config index - - if (sim%my_id == 0) write(*,*) "----- Initialising particles for group '", part_group_configs(j)%id, "' with coupling scheme '", part_group_configs(j)%coupling_scheme, "' -----" select case(trim(part_group_configs(j)%coupling_scheme)) - - !> initialisation of runaway electrons - case('rep') - call initialise_group_RE(sim, i) - !> initialisation of energetic particles - case ('epc', 'epp', 'epf') - call initialise_group_EP(sim, i) - !> ics and ncs schemes don't need to initialise particles - case ('ics', 'ncs') - if (sim%my_id == 0) write(*,*) " Initialisation skipped for ncs / ics, particles are not all initialised at once" - !> default case - give error - case default - if (sim%my_id == 0) write(*,*) "ERROR : No particle coupling scheme selected for group '", part_group_configs(j)%id, "'" - stop 1 - end select - enddo - - end subroutine initialise_particles_for_sim - - subroutine initialise_group_RE(sim, group_num) - use mod_pcg32_rng - - implicit none - class(particle_sim), intent(inout) :: sim - integer, intent(in) :: group_num - character(len=50) :: init_function_name, init_pdf_name - type(type_part_group_config) :: config - - config = part_group_configs(matching_part_config_indices(group_num)) - init_function_name = config%init_function - init_pdf_name = config%init_pdf - - select case (trim(init_function_name)) - case ("basic") - !> Call initialiser subroutine - if (sim%my_id == 0) write(*,*) " Using the 'basic_initialization' function, with PDF: ", trim(init_pdf_name) - call basic_initialization(sim, group_num, pcg32_rng(), init_pdf_name, config%re_energy, config%re_pitch, config%re_std_energy) - - !> Set particle charge and weight - select type (particles => sim%groups(group_num)%particles) - type is (particle_kinetic_relativistic) - particles(:)%q = -1 !> default electron charge - particles(:)%weight = config%num_re / config%n_particles - end select - - if (sim%my_id == 0) then - write(*,*) "----- Finished initialisation for group '", config%id, "' with coupling scheme '", config%coupling_scheme, "' -----" - write(*,*) "" + case('ics', 'ncs') + if (sim%my_id == 0) then + write(*,*) "Particle initialisation skipped for ics/ncs, particles are born elsewhere" endif case default if (sim%my_id == 0) then - write(*,*) "ERROR: ", trim(init_function_name), " is not a valid initialisation function " - write(*,*) " for group '", config%id, "' with coupling scheme: '", config%coupling_scheme, "'" + write(*,*) "----- Initialising group '", part_group_configs(j)%id, "' -----" + write(*,*) " init_function : '", trim(part_group_configs(j)%init_function), "'" + write(*,*) " init_pdf : '", trim(part_group_configs(j)%init_pdf), "'" endif - stop 1 - end select - - end subroutine initialise_group_RE + call initialise_group(sim, i) + endselect + enddo + end subroutine initialise_particles_for_sim - subroutine initialise_group_EP(sim, group_num) + !> Initialise a single particle group using namelist defined init_function and init_pdf + subroutine initialise_group(sim, group_num) use mod_pcg32_rng implicit none + + !> i/o vars class(particle_sim), intent(inout) :: sim integer, intent(in) :: group_num - character(len=50) :: init_function_name - type(type_part_group_config) :: config - real*8 :: T_maxwell - integer :: n_phi_planes_in - real*8 :: n_particles_total + + !> internal vars + type(type_part_group_config) :: config + type(spatial_pdf) :: space_pdf + logical :: exists config = part_group_configs(matching_part_config_indices(group_num)) - init_function_name = config%init_function - T_maxwell = config%T_maxwell - n_phi_planes_in = config%n_phi_planes - n_particles_total = config%n_particles_total - - select case (trim(init_function_name)) - case ("maxwell") - - !> Initialise particles via a maxwellian distribution - if (sim%my_id == 0) write(*,*) " Using the initialise_particles_h_mu_psi_phiplanes subroutine" - call initialise_particles_H_mu_psi_phiplanes(sim%groups(group_num)%particles, sim%fields, pcg32_rng(),sim%groups(group_num)%mass, & - uniform_space=.true., uniform_space_rej_f = f_toroidal_flux, & - uniform_space_rej_vars=[1], charge=1, T_maxwell=T_maxwell, n_phi_planes_in=n_phi_planes_in) - - !> Adjust particle weights so sum(weights) = n_particles_total - call adjust_particle_weights(sim%groups(group_num)%particles, n_particles_total) - if (sim%my_id == 0) then - write(*,*) "----- Finished initialisation for group '", config%id, "' with coupling scheme '", config%coupling_scheme, "' -----" - write(*,'(A,ES12.3)') " Number of super particles : ", config%n_particles - write(*,'(A,ES12.3)') " Total number of particles : ", n_particles_total - write(*,'(A,ES12.3)') " Particle weights adjusted to : ", sim%groups(group_num)%particles(1)%weight - write(*,'(A,ES12.3,A)') " Temperature of Maxwellian : ", T_maxwell, "[eV]" - write(*,*) "" + + !> Compatability check : init_function x particle_type + select type(particles => sim%groups(group_num)%particles) + type is (particle_kinetic_relativistic) + if (trim(config%init_function) /= 'gaussian_re') then + write(*,*) "ERROR : particle_kinetic_relativistic requred init_function='gaussian_re'" + write(*,*) " got init_function='",trim(config%init_function),"' for group '", trim(config%id), "'" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif - case default - if (sim%my_id == 0) then - write(*,*) "ERROR: ", trim(init_function_name), " is not a valid initialisation function " - write(*,*) " for group '", config%id, "' with coupling scheme: '", config%coupling_scheme, "'" + type is (particle_kinetic_leapfrog) + if (trim(config%init_function) == 'gaussian_re') then + write(*,*) "ERROR (initialise_group): particle_kinetic_leapfrog incompatible with init_function='gaussian_re'" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif - stop 1 end select - end subroutine initialise_group_EP + !> Set real-space pdf + space_pdf = spatial_pdf_from_name(trim(config%init_pdf)) - !> rejection function to produce EP spatial - !> distribution from ITPA TAE benchmark - !> A. Könies et al 2018 Nucl. Fusion 58 126027 - !> https://doi.org/10.1088/1741-4326/aae4e6 - pure function f_toroidal_flux(n, P, grad_P) result(f) - integer, intent(in) :: n - real*8, intent(in) :: P(n), grad_P(3,n) - real*8 :: s, psi_norm, coeff(0:3) - real*4 :: f + !> Select phase-space initialiser and sample in both real- and phase-space + select case(trim(config%init_function)) - ! central densiy should be 1.44131x10^17 + !> EPs : Maxwellian distribution + case ('maxwell') + if (sim%my_id == 0) then + write(*,*) " Sampler: initialise_particles_H_mu_psi_phiplanes" + write(*, '(A,ES12.3,A)') " T_maxwell : ", config%T_maxwell, " [eV]" + write(*, '(A,I0)') " n_phi_planes : ", config%n_phi_planes + endif - coeff(0)=0.49123 - coeff(1)=0.298228 - coeff(2)=0.198739 - coeff(3)=0.521298 + if (associated(space_pdf%f)) then + call initialise_particles_H_mu_psi_phiplanes( & + sim%groups(group_num)%particles, sim%fields, pcg32_rng(), & + sim%groups(group_num)%mass, uniform_space=.true., & + uniform_space_rej_f=space_pdf%f, & + uniform_space_rej_vars=space_pdf%vars, charge=1, & + T_maxwell=config%T_maxwell, n_phi_planes_in=config%n_phi_planes) + else + call initialise_particles_H_mu_psi_phiplanes( & + sim%groups(group_num)%particles, sim%fields, pcg32_rng(), & + sim%groups(group_num)%mass, uniform_space=.true., & + charge=1, T_maxwell=config%T_maxwell, & + n_phi_planes_in=config%n_phi_planes) + endif - psi_norm = max((P(1) - ES%Psi_axis) / ( ES%Psi_bnd - ES%Psi_axis),0.d0) + !> EPs : experimental distribution - f(R,Z,energy,pitch), expects a .h5 file - see mod_import_experimental_dist + case ('experimental') + if (sim%my_id == 0) write(*,*) " Sampler: import_particles - experimental distribution from 'experimental_dist.h5' file" - s = 0.957 * psi_norm + 0.043 * psi_norm**2 + !> check for existence of experimental_dist.h5 file + inquire(file="experimental_dist.h5", exist=exists) + if (.not. exists) then + if (sim%my_id == 0) write(*,*) "ERROR (initialise_group): init_function='experimental' requires 'experimental_dist.h5' file in working directory" + stop 1 + endif - f = coeff(3)*exp(-coeff(2)/coeff(1)*(tanh((sqrt(s)-coeff(0))/coeff(2)))) + call import_particles(sim%groups(group_num)%particles, sim%fields, & + "experimental_dist.h5", pcg32_rng(), sim%groups(group_num)%mass, & + n_phi_planes_in=config%n_phi_planes, fraction_phi_planes=1.d0) - end function f_toroidal_flux + !> REs : Monoenergetic or Gaussian energy distribution + case('gaussian_re') + if (sim%my_id == 0) then + write(*,*) " Sampler: gaussian_re_initialisation" + write(*,'(A,ES12.3,A)') " Energy : ", config%re_energy, " [eV]" + write(*,'(A,ES12.3,A)') " std : ", config%re_std_energy, " [eV]" + write(*,'(A,ES12.3)') " Pitch : ", config%re_pitch + endif + + if (associated(space_pdf%f)) then + call initialise_gaussian_re(sim, group_num, pcg32_rng(), & + space_pdf, config%re_energy, config%re_pitch, config%re_std_energy) + else + call initialise_gaussian_re(sim, group_num, pcg32_rng(), & + energy=config%re_energy, pitch=config%re_pitch, std_energy=config%re_std_energy) + endif + + !> Unknown init functions + case default + if (sim%my_id == 0) then + write(*,*) "ERROR (initialise_group): unknown init_function '", & + trim(config%init_function), "' for group '", trim(config%id), "'" + write(*,*) " Valid init_functions: 'maxwell', 'gaussian_re', 'experimental'" + endif + stop 1 + + end select + + !> Finalize weights and charge + select type(particles => sim%groups(group_num)%particles) + type is (particle_kinetic_relativistic) + particles(:)%q = -1 + end select + + !> Finalize weights + sim%groups(group_num)%particles(:)%weight = 1.d0 + call adjust_particle_weights(sim%groups(group_num)%particles, config%n_particles_total) + + !> Summary print + if (sim%my_id == 0) then + write(*,'(A,ES12.3)') " Number of super particles : ", config%n_particles + write(*,'(A,ES12.3)') " Number of actual particles : ", config%n_particles_total + write(*,'(A,ES12.3)') " Particle weight : ", sim%groups(group_num)%particles(1)%weight + write(*,*) "----- Finished initialisation for group '", trim(config%id), "'-----" + write(*,*) "" + endif + + end subroutine initialise_group end module mod_initialise_particles \ No newline at end of file diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 new file mode 100644 index 0000000000..6010b9836a --- /dev/null +++ b/particles/initialisers/mod_rej_f.f90 @@ -0,0 +1,217 @@ +!> Module defining abstract interfaces and implemtation for +!> real space rejection functions used during particle initialisation. +!> +!> A rejection funciton (rej_f) maps field values at a sample point to +!> an acceptance probability in [0,1]. The spatial_pdf derived type bundles +!> a procedure pointer with the variable list it expects, so both pieces of +!> information always travel together +!> +!> ## Variable index convention (vars(:) entries) +!> +!> Positive index -> JOREK variable number (model-dependent) +!> -1 -> R (major radius) +!> -2 -> Z (height) +!> -3 -> Phi (toroidal angle) +!> +!> The P(:) array passed to a rej_f function is ordered to match vars(:), +!> so a function must be paired with the vars array it was written for. +!> Use spatial_pdf_from_name() to get a correctly-paired spatial_pdf. + +module mod_rej_f + use equil_info ! (R_axis,Z_axis), etc.. + use mod_model_settings + + implicit none + private + + public :: spatial_pdf + public :: rej_f + public :: itpa_tae_pdf, RZ_pdf, analytical_pdf, current_pdf + public :: spatial_pdf_from_name + + ! ============================================================================= + ! Interface + ! ============================================================================= + + !> Primary interface for spatial rejection functions used in particle init + !> + !> @param n Number of field values in P and gradP + !> @param P Field values at the sample point, ordered as vars(:) + !> @param gradP Gradients (3,n); may be unused + !> @return Acceptance probability in [0,1] + abstract interface + function rej_f(n, P, gradP) + implicit none + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: gradP + real*4 :: rej_f + end function rej_f + end interface + + ! ============================================================================= + ! spatial_pdf -> Type bundling rejection f and vars it expects + ! ============================================================================= + + !> Bundle type pairing a spatial rejection function with + !> the JOREK / (R,Z,Phi) variable list it expects + !> + !> Always construct via spatial_pdf_from_name() or by setting + !> both components together - the procedure f is written to + !> expect exactly size(vars) values in P(:), ordered to match vars(:). + !> + !> Example (custom profile) + !> + !> type(spatial_pdf) :: pdf + !> pdf%f => my_rej_function + !> pdf%vars = [var_psi, -1] ! psi and R + type :: spatial_pdf + procedure(rej_f), nopass, pointer :: f => null() + integer, allocatable :: vars(:) + end type spatial_pdf + + + contains + ! ============================================================================= + ! Constructor + ! ============================================================================= + + !> Returns a correctly-paired spatial_pdf type + !> + !> Recongnised names: + !> 'itpa_tae' - reproduce EP dist in ITPA TAE benchmark + !> 'RZ' - weight by 1/(1+r_minor)^2 + !> 'analytical' - weight by (1-(r/a)^2)^nu + !> 'current' - weight by normalised j_tor + !> 'none' - no rejection (f => null, vars unallocated) + function spatial_pdf_from_name(name) result(pdf) + character(len=*), intent(in) :: name + type(spatial_pdf) :: pdf + + select case(trim(name)) + + case('itpa_tae') + pdf%f => itpa_tae_pdf + pdf%vars = [var_psi] + + case('RZ') + pdf%f => RZ_pdf + pdf%vars = [-2, -1] + + case('analytical') + pdf%f => analytical_pdf + pdf%vars = [-2, -1] + + case('current') + !> NOTE : var_zj = 0 in fullMHD (j_tor is not a stored variable) + !> This rej_f is only valid for models where var_zj > 0 (e.g. model600) + if (var_zj == 0) then + write(*,*) "ERROR (mod_rej_f): 'current' spatial PDF requires model" + write(*,*) " with var_zj > 0. j_tor is not a stored variable in this model" + stop 1 + endif + pdf%f => current_pdf + pdf%vars = [-1, var_zj] + + case('none') + !> Leave f => null() and vars unallocated + !> caller should handle this case by skipping use of rej_f + !> no rej_f should be passed to initialiser, sampling will be uniform in space + + case default + write(*,*) "ERROR (mod_rej_f): Unknown spatial PDF name '", trim(name), "'" + write(*,*) " Valid names: 'itpa_tae', 'RZ', 'analytical', " + write(*,*) " 'current', 'none'" + stop 1 + end select + end function spatial_pdf_from_name + + + ! ============================================================================= + ! Rejection functions + ! ============================================================================= + + + !> rejection function to produce EP spatial, distribution from ITPA TAE benchmark + !> A. Könies et al 2018 Nucl. Fusion 58 126027, https://doi.org/10.1088/1741-4326/aae4e6 + !> + !> Expected vars : [var_psi] + !> P(1) = psi + pure function itpa_tae_pdf(n, P, gradP) result(f) + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: gradP + real*4 :: f + + !> internal vars + real*8 :: s, psi_norm, coeff(0:3) + + coeff(0)=0.49123 + coeff(1)=0.298228 + coeff(2)=0.198739 + coeff(3)=0.521298 + + psi_norm = max((P(1) - ES%Psi_axis) / ( ES%Psi_bnd - ES%Psi_axis),0.d0) + s = 0.957 * psi_norm + 0.043 * psi_norm**2 + + f = real(coeff(3)*exp(-coeff(2)/coeff(1)*(tanh((sqrt(s)-coeff(0))/coeff(2)))), 4) + end function itpa_tae_pdf + + !> Weight as 1/(1 + r_minor)^2 - broad profile peaked on axis + !> + !> Expected vars : [-2, -1] + !> P(1) = Z + !> P(2) = R + pure function RZ_pdf(n, P, gradP) result(f) + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: gradP + real*4 :: f + + real*8 :: minor_r + + minor_r = sqrt((P(1) - ES%Z_axis)**2 + (P(2) - ES%R_axis)**2) + f = real(1.d0 / (1.d0 + minor_r)**2, 4) + end function RZ_pdf + + !> Analytically prescribed profile (1 - (r_minor/a)^2)^nu, nu=2 + !> + !> Expected vars : [-2, -1] + !> P(1) = Z + !> P(2) = R + pure function analytical_pdf(n, P, gradP) result(f) + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: gradP + real*4 :: f + + real*8, parameter :: nu = 2.d0 + real*8 :: minor_r + + minor_r = sqrt((P(1) - ES%Z_axis)**2 + (P(2) - ES%R_axis)**2) + f = real((1.d0 - (minor_r/ES%LCFS_a)**2)**nu, 4) + f = max(f, 0.0e0) !> clamp outside LCFS + end function analytical_pdf + + !> Weight proportional to normalised toroidal current density j_tor + !> + !> Expected vars : [-1, var_zj] + !> P(1) = R + !> P(2) = j_tor + !> + !> NOTE : Won't work for fMHD models, see note in constructor above + !> NOTE : jzmin and jzmax are currently hardcoded, needs better handling + pure function current_pdf(n, P, gradP) result(f) + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: gradP + real*4 :: f + + !> TODO: make jzmin.jzmax namelist params + real*8, parameter :: jzmax = 3.0d0 / 10.0d0 + real*8, parameter :: jzmin = 1.239d-4 / 11.0d0 + + f = real((P(2)/P(1) - jzmin) / (jzmax - jzmin), 4) + f = max(f, 0.0e0) + end function current_pdf +end module mod_rej_f From 22c299b62a22f7c0430a9892feebedff053c7b1a Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Mon, 10 Aug 2026 17:09:43 +0200 Subject: [PATCH 02/16] Updated stop clauses to be MPI_ABORT in particle initialiser routines --- particles/initialisers/mod_initialise_particles.f90 | 6 +++--- particles/initialisers/mod_rej_f.f90 | 6 ++++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 9d93f91041..6b1dbdffb8 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -32,7 +32,7 @@ module mod_initialise_particles implicit none integer :: ierr !> mpi error code - + contains !> Top-level initialisation loop over all particle groups @@ -131,7 +131,7 @@ subroutine initialise_group(sim, group_num) inquire(file="experimental_dist.h5", exist=exists) if (.not. exists) then if (sim%my_id == 0) write(*,*) "ERROR (initialise_group): init_function='experimental' requires 'experimental_dist.h5' file in working directory" - stop 1 + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif call import_particles(sim%groups(group_num)%particles, sim%fields, & @@ -162,7 +162,7 @@ subroutine initialise_group(sim, group_num) trim(config%init_function), "' for group '", trim(config%id), "'" write(*,*) " Valid init_functions: 'maxwell', 'gaussian_re', 'experimental'" endif - stop 1 + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) end select diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index 6010b9836a..333739b2e1 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -20,8 +20,10 @@ module mod_rej_f use equil_info ! (R_axis,Z_axis), etc.. use mod_model_settings + use mpi implicit none + integer :: ierr ! mpi error code private public :: spatial_pdf @@ -108,7 +110,7 @@ function spatial_pdf_from_name(name) result(pdf) if (var_zj == 0) then write(*,*) "ERROR (mod_rej_f): 'current' spatial PDF requires model" write(*,*) " with var_zj > 0. j_tor is not a stored variable in this model" - stop 1 + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif pdf%f => current_pdf pdf%vars = [-1, var_zj] @@ -122,7 +124,7 @@ function spatial_pdf_from_name(name) result(pdf) write(*,*) "ERROR (mod_rej_f): Unknown spatial PDF name '", trim(name), "'" write(*,*) " Valid names: 'itpa_tae', 'RZ', 'analytical', " write(*,*) " 'current', 'none'" - stop 1 + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) end select end function spatial_pdf_from_name From b0ae17aebe5b0d92d041d41caf52234fdb8c81d7 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 13:51:12 +0200 Subject: [PATCH 03/16] Removed redundant num_re --- algexpr2fort | Bin 0 -> 21240 bytes communication/broadcast_phys.f90 | 2 -- generate_code | 0 models/mod_log_params.f90 | 2 +- models/phys_module.f90 | 3 +- models/preset_parameters.f90 | 5 +-- particles/initialisers/initialisers_RE.f90 | 4 +-- .../initialisers/mod_initialise_particles.f90 | 20 +++++------ particles/mod_coupling_settings.f90 | 31 ++++++++++++++++++ reg_tests/testcases/particle_rep_600/input | 4 +-- 10 files changed, 51 insertions(+), 20 deletions(-) create mode 100755 algexpr2fort create mode 100644 generate_code diff --git a/algexpr2fort b/algexpr2fort new file mode 100755 index 0000000000000000000000000000000000000000..ba29d1d0ac04a1f5246a219c9d67dfcbdc3382a6 GIT binary patch literal 21240 zcmeHPdyE^$d7mZkbSIy6x`%9^Y&W7DKjhlv9Z%Hxbhc!ko;*j^TdIoIUN6ZdxfZ#k zc9*Bq6-^ZQQK>J;eZeWxMnLm-(_A9+8AlZ!(B4N4(-YzL27;*^-l3g>U6{;X- zz;MVHlU*;NiWaJV3^%HBGGx6*elapeilrrBSXKpN7__6?kfjY~ss&SmoF0IviB>W~iwQ#>-ZADm7lt$IF&o zUx+VEPsOKF3AdVlu!x>q0LUs*dLv#L-VXmS^hus%-}NBy z$bzEb6rrV~UF1|LDu!i?O0HVBJps1s<#YGlM@?jk$U?z8y;-l+EaGf#a2iz2a4M!} zxnPhmh`TwXY}Ab0NhPb+OgmdPNHsA*%B6f^i(V@0DhASLGAIdi+A9XM770(vg48U% z?wTlgVPcxNDd`xrV(V7Xt~#dfopf5oWKwEjH3`04#MvUVAxs9USgq8EkhvFgh;D;A zRX6e#!;?;BDRREm-hFdBcIuOf%_+g?B;4-ZbKvmN!?XMKX|y8hF%kKT;>_cSMmR*& z{w8qz;d8%%uGb?56|a7)g>D%Y>&YBR!|(IxB&1KN^dG;1E*lYdDf#t>(~?j9SFt4h zUVzVIEeXcbX(^|EPkoP5EWqb}PqIw`KG%Wx69K;bUn3^cbbucmQ_=xG<=su)2y`RR zjX*a7-3W9e@c%ag-&^~U_cKp@yFc>_z26uYBJ<3W7g>HS^VFa9Uy(ss-ttYbm+yKD z=d~Nsz=$r9?An{lNO%1$(KMC0_L`*sn&=48OOpOF(KMyGc1hBIN;FMnu3eP$OGMKY z=Gq0&$1>;NzWO&nu6_=um!F8G#TRJ4a`jh<_W`iv(qmT&iM1P^1}+4ZKlv|n=fC!N z=KMD^PrY^h;L+X5rR0|~ufBR=2*tf=6z>16RX<>U;qnPgj56Ek?2TR~j}(H{Um(Wv zn}rKk3Qb$8|I1SUli!ix{PoN$-+4Im%6rj__QlNCKJbPc9{QArusxOc?6$v4Afx)I ztn4Bx_Ql@61V?-9%8L9*|DS)pSq3u6eDd2v5rpK@2g}R2ZZ8g^(7*Kg>;DDn)o0!V z=!vDptIr6C^Hibp&yeiuEUEi?FVO_~-8S7c0On`pEGXm?XL0^JC7 zBhZaNHv-)VbR*D>KsN&22y`RRjllou2+;G8N;Pkm-EpH_G#6^lBt0ui6t>dS5h1p3 zkL`pwRy1wXG3ZfQ%yMH@yL=`Fy%?SoS%z(t&v+J|S;c%$iP%FAaZq;Q;m7h?pQ4)t z{LEv6QaB4p&ndLqq7TwDU%Cm0j=%#K2 zx)JC`pc{d11iBIEMxYykZUnj!_#Z~#;&578)D-0P#6{(Y>E930uY~BYhUmW!(f<^p z-w)9PDkIMCt`Pk|h~6Kf9}CgN5dHBG{aYdWQi!H^)F>$+{whSj9-`k4(ccTv^o|^* z8-0KURa#MuP6B(xmy&!*5~UpKdor4 z2dyh9(L24AIDb4*mG!y4DlPeV*ehwy7t2XWzpms`>3{LW@0EwMoKn|&kAlY)EGRgy z;Ij&b?f357ITPE2H=pl~C6kFQiHX?c#AIsXf#lZMro(1FmN7gfh$oZx3Sa#Qz4ICa z4kNtWiFaLRkI$rP$;qu#;_&eqr*7MpU5w!!(Xdmo2x7{-n!a9OeOBpQN@{pm-M&Z6sMYJg6rDxnS+aZ#v!4$;Ly z|K2-$r?%NhHuiRWOI^EmI`1kO!*)76CfxiXet4X?U5yv*{ER9Y!*)8H$88^y7w552 z`Jfm;d)JB8jjOT#UwijP-vDzga%f#70?Xbp zh`DJ+-xk^`ZPkTUCn95zH8(?a-1N}c1hMI)W#_likh}j_Y>&T(Vsc_)T8+4|c*SsX zrEQD$yy>JsRVoWO*Q&Lg?E!IJ(?4&Pz%_knrV>e?pGr&uG3Ln%XCh#eG63yryihU7 zDX@iV4G3tQp|>SDSpb%;SpHc$F5dqYRt#c2!n3FTklx|Qv`^(K6d>f6Bj6U@tRb&yVZl|gDwX9mHcAvN7?n@cr7AD_Eud1J&ypN4 zZ#h#b%E`y6c9qkH$w_S6 zquqK=$t-V;FPe609Fs-NFcTP`d{oXWpc75ERn~tNPxpE!P#^QyVW5yM=ZM^O}D3cU5bjE2HQ+ysB zy6VJvGgc|M!>rAwvzCo*2zxH5-+t>F?^C1kqi8z*V4u2I$gOI@2Ydql=0!86rE0}Q z$xX*xz{rCqdgD*#3-TE(z)g-SQTlH@c>ArN3CGn^j z5WG$e3ZpdLRI|cX{u&W$%O7k^Vq5js3ZBom@-e05`DZ{76~lt}^#<73T{j|@+Uk!A zzK_tVzfSP|gI4}cqTTv6DsC3MFRxV}Grx2?$dY&ztMyD<{(3?0d52R}@Li2oNmTIN zk1+ps@x&A9Rwk&ZpjN&zMa4%}GQo)oernz-j*3_-Ls3z&u>=2(4t!b}tV~gH=L)6< zPLKHRa#(OP7#d2-{V1js2)5Htf*%osE7oSy;A>ba@w$fZ8|~`g=YaAfhSR{BNzSz4 zehsfm^Sh}L|fP5-?I{1MbQJP_Xx{xB9J;r^0U z`b%n@;emTv@+ptjM2pXn$jWi%r@p!O!ImeuSKH;m#Z;Uij`p z&he6oDz@9>=A4dcp9s!f`~{k&{?W*#2eYUr5gS=j#X7-w|MhnzE_fay5C zH&Cxs&cGx<>2&w5MYUQel&eOb%cbjk4$nTiTi?BZ7k$KF@4f>&X7}j__Ut*b`>1|& zcE`Tm;O_d#{j-nG?SyNr(e>RKWs});*q_Jg`q4*svJ2i=-}nwfI5{dj78~%LgIlcsUX%D3>Fvv0vT92*rYhg_@hg;PRi;5$S65()GwS|o~ux+^d12|m>(Ter&j zgyCVCm#ur~n2EfZtry|ku!~L7vJ2IQD4TW6d8TOb$T;H3D6w(9nshZ+#u(7>@7P73 zZW-}@G=gS&7+Ue)@29b7S|a>CBc87)vD9UpAkDPsSrXe%D|`NaPYkkP@>fDl&J@s~ zJFC)W_AV${iksL3o8GGsINiuFG{#eEEDI?erUe(?|62__Ph?dr0lu= zVSD9H><2>6>e!ywORp(=n%_{OI#QfD|1n}B9mhecY|rbe7gPcHc_8UhA|Ju@I0$+c z$M(GbdR-Nm|LeH?tjF*q?CIGa)4Wc*t_%*U%=j8{`MC@<&!nQVJ+J4|WPmi`gB9`T z`acah+3!$#yzZ+h`wIbk%=%Wu@0Y=&m~j1snrF~oIi>JQMupWE;(Q)Q&_1T@W6HkP z(DV_-X7}RPLiRkLyr2}H3SBU^XScr@vgdW_I(6arIW_Bt%YRYXbNqS!B4WP9IUh96 z;rw5Kgo?=ayv{ur^Gz`Lx1@0Xe++5Rp4Y=mG2iU>TCAJbkJ(E=263qF(pWe1cdo+i zjcJBof$iFWJ+Esn-w!4d>-kiZW#YCD`iufzGfWfMxEnq$Ks$MGAmrTViyuh*}S zgNej?KD9FcZ$S{WZ>;Mle2WYlD#+z#_)f^4*Z&!1-zk4~OY14}#pUO9zzJo~`LI5f zHJIK9K|H125V8H<2FFL(o?#S7JNwTlhuc|E`C*db2oSzDq{8;;Qqk>U-X^~_nn~qxm`Z>^^o&${=*&Y OA5$IY`H(?~6#otCeaAWg literal 0 HcmV?d00001 diff --git a/communication/broadcast_phys.f90 b/communication/broadcast_phys.f90 index 52debf0e7c..cf04247cfd 100644 --- a/communication/broadcast_phys.f90 +++ b/communication/broadcast_phys.f90 @@ -924,7 +924,6 @@ subroutine broadcast_phys(my_id) call MPI_PACK(part_group_configs%homma2020_alpha, n_part_groups_max, MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr) call MPI_PACK(part_group_configs%ics_group_idx, n_part_groups_max, MPI_INTEGER,buffer,bufsize,position,MPI_COMM_WORLD,ierr) - call MPI_PACK(part_group_configs%num_re, n_part_groups_max, MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr) call MPI_PACK(part_group_configs%re_energy, n_part_groups_max, MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr) call MPI_PACK(part_group_configs%re_std_energy, n_part_groups_max, MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr) call MPI_PACK(part_group_configs%re_pitch, n_part_groups_max, MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr) @@ -1917,7 +1916,6 @@ subroutine broadcast_phys(my_id) call MPI_UNPACK(buffer,bufsize,position,part_group_configs%homma2020_alpha, n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,bufsize,position,part_group_configs%ics_group_idx, n_part_groups_max,MPI_INTEGER,MPI_COMM_WORLD,ierr) - call MPI_UNPACK(buffer,bufsize,position,part_group_configs%num_re, n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_energy, n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_std_energy, n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr) call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_pitch, n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr) diff --git a/generate_code b/generate_code new file mode 100644 index 0000000000..e69de29bb2 diff --git a/models/mod_log_params.f90 b/models/mod_log_params.f90 index 50c4aabe16..1d96a7bf35 100644 --- a/models/mod_log_params.f90 +++ b/models/mod_log_params.f90 @@ -1142,7 +1142,7 @@ subroutine log_parameters(my_id, short) ! rep (runaway electrons, only pressure coupling for now) ----- if (sim%groups(group_num)%coupling_scheme .eq. 'rep') then - write(*,REAL_FMT) 'num_re, ',part_group_configs(group_num)%num_re + write(*,REAL_FMT) 'n_particles_total, ',part_group_configs(group_num)%n_particles_total write(*,REAL_FMT) 're_energy, ',part_group_configs(group_num)%re_energy write(*,REAL_FMT) 're_std_energy, ',part_group_configs(group_num)%re_std_energy write(*,REAL_FMT) 're_pitch, ',part_group_configs(group_num)%re_pitch diff --git a/models/phys_module.f90 b/models/phys_module.f90 index 89a8f2a5c7..8754e9519f 100644 --- a/models/phys_module.f90 +++ b/models/phys_module.f90 @@ -1112,7 +1112,6 @@ module phys_module ! ================ for runaway electrons ('rep' coupling scheme) particles =============== - real*8 :: num_re !< number of runaway electrons in the group real*8 :: re_energy !< energy [eV] of the runaway electrons in the group real*8 :: re_std_energy !< standard deviation of the energy [eV] of the runaway electrons in the group real*8 :: re_pitch !< pitch between RE momentum and magnetic field line (i.e. p_re_par/p_re_tot) @@ -1123,6 +1122,8 @@ module phys_module !< for example n_phi_planes=4, n_particles=1e4 then only 250 particles are initialised !< each particle is then copied multiple (3) times around the torus with angle 2pi/n_phi_planes (= pi/2) !< if n_phi_planes=int*n_period then projected particle quantities are initialised as 0 for n_tor>1 + + ! =============== for rep and epf =================== real*8 :: n_particles_total !< Total number of particles to simulate (ie sum(weights)) !!NOT n_particles - total number of super/numeric-particles diff --git a/models/preset_parameters.f90 b/models/preset_parameters.f90 index 8732e1a6e2..f600d39569 100644 --- a/models/preset_parameters.f90 +++ b/models/preset_parameters.f90 @@ -924,8 +924,7 @@ subroutine preset_parameters part_group_configs(:)%homma2020_alpha = 1.5d0 part_group_configs(:)%ics_group_idx = -1 -!----- specific to rep -part_group_configs(:)%num_re = 0.d0 +!----- specific to rep part_group_configs(:)%re_energy = 0.d0 part_group_configs(:)%re_std_energy = 0.d0 part_group_configs(:)%re_pitch = 0.d0 @@ -933,6 +932,8 @@ subroutine preset_parameters !----- specific to epf part_group_configs(:)%T_maxwell = 0.d0 part_group_configs(:)%n_phi_planes = 0 + +!----- specific to rep and epf part_group_configs(:)%n_particles_total = 0.d0 diff --git a/particles/initialisers/initialisers_RE.f90 b/particles/initialisers/initialisers_RE.f90 index e93ad59cc7..eea5bedf0d 100644 --- a/particles/initialisers/initialisers_RE.f90 +++ b/particles/initialisers/initialisers_RE.f90 @@ -15,7 +15,7 @@ module initialisers_RE implicit none contains -subroutine initialise_gaussian_re(sim, group_num, rng, space_pdf, energy, pitch, std_energy) +subroutine initialise_re_gaussian(sim, group_num, rng, space_pdf, energy, pitch, std_energy) use phys_module, only: tstep_particles use mod_kinetic_relativistic use mod_sampling, only: boxmueller_transform @@ -101,6 +101,6 @@ subroutine initialise_gaussian_re(sim, group_num, rng, space_pdf, energy, pitch, deallocate(p_par) deallocate(p_perp) -end subroutine initialise_gaussian_re +end subroutine initialise_re_gaussian end module initialisers_RE diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 6b1dbdffb8..3848a9b470 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -13,7 +13,7 @@ !> Valid init_functions: !> 'maxwell' -> Maxwellian Energy distribution, uniform pitch. Spatial dist set by init_pdf (or uniform if no init_pdf specified) !> Temperature specified by T_maxwell in input file. Only valid for particle_kinetic_leapfrog -!> 'gaussian_re' -> Monoenergetic or Gaussian energy distribution, delta function in pitch. +!> 're_gaussian' -> Monoenergetic or Gaussian energy distribution, delta function in pitch. !> Spatial dist set by init_pdf (or uniform if no init_pdf specified) !> re_energy, re_std_energy, re_pitch specified in input file. Only valid for particle_kinetic_relativistic !> 'experimental' -> Samples from external F(R, Z, energy, pitch) distribution, expects 'experimental_dist.h5' hdf5 file with the correct format @@ -82,14 +82,14 @@ subroutine initialise_group(sim, group_num) !> Compatability check : init_function x particle_type select type(particles => sim%groups(group_num)%particles) type is (particle_kinetic_relativistic) - if (trim(config%init_function) /= 'gaussian_re') then - write(*,*) "ERROR : particle_kinetic_relativistic requred init_function='gaussian_re'" + if (trim(config%init_function) /= 're_gaussian') then + write(*,*) "ERROR : particle_kinetic_relativistic requred init_function='re_gaussian'" write(*,*) " got init_function='",trim(config%init_function),"' for group '", trim(config%id), "'" call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif type is (particle_kinetic_leapfrog) - if (trim(config%init_function) == 'gaussian_re') then - write(*,*) "ERROR (initialise_group): particle_kinetic_leapfrog incompatible with init_function='gaussian_re'" + if (trim(config%init_function) == 're_gaussian') then + write(*,*) "ERROR (initialise_group): particle_kinetic_leapfrog incompatible with init_function='re_gaussian'" call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif end select @@ -139,19 +139,19 @@ subroutine initialise_group(sim, group_num) n_phi_planes_in=config%n_phi_planes, fraction_phi_planes=1.d0) !> REs : Monoenergetic or Gaussian energy distribution - case('gaussian_re') + case('re_gaussian') if (sim%my_id == 0) then - write(*,*) " Sampler: gaussian_re_initialisation" + write(*,*) " Sampler: re_gaussian_initialisation" write(*,'(A,ES12.3,A)') " Energy : ", config%re_energy, " [eV]" write(*,'(A,ES12.3,A)') " std : ", config%re_std_energy, " [eV]" write(*,'(A,ES12.3)') " Pitch : ", config%re_pitch endif if (associated(space_pdf%f)) then - call initialise_gaussian_re(sim, group_num, pcg32_rng(), & + call initialise_re_gaussian(sim, group_num, pcg32_rng(), & space_pdf, config%re_energy, config%re_pitch, config%re_std_energy) else - call initialise_gaussian_re(sim, group_num, pcg32_rng(), & + call initialise_re_gaussian(sim, group_num, pcg32_rng(), & energy=config%re_energy, pitch=config%re_pitch, std_energy=config%re_std_energy) endif @@ -160,7 +160,7 @@ subroutine initialise_group(sim, group_num) if (sim%my_id == 0) then write(*,*) "ERROR (initialise_group): unknown init_function '", & trim(config%init_function), "' for group '", trim(config%id), "'" - write(*,*) " Valid init_functions: 'maxwell', 'gaussian_re', 'experimental'" + write(*,*) " Valid init_functions: 'maxwell', 're_gaussian', 'experimental'" endif call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) diff --git a/particles/mod_coupling_settings.f90 b/particles/mod_coupling_settings.f90 index 9d06adeb67..81b7d2cce8 100644 --- a/particles/mod_coupling_settings.f90 +++ b/particles/mod_coupling_settings.f90 @@ -50,6 +50,7 @@ subroutine check_compatibility_and_determine_coupling_schemes() case ('rep') call check_no_epf_params(group_num) call check_no_ics_ncs_params(group_num) + call check_compatibility_rep(group_num) use_rep = .true. case ('epc') write(*,*) "ERROR: coupling scheme 'epc' is not yet implemented" @@ -124,6 +125,36 @@ subroutine check_compatibility_ics(group_num) end subroutine check_compatibility_ics +subroutine check_compatibility_rep(group_num) + implicit none + integer :: group_num + + !> Check initialisation parameters + if (trim(part_group_configs(group_num)%init_function) .eq. 're_gaussian') then + if (part_group_configs(group_num)%re_energy .eq. 0.d0) then + write(*,*) "ERROR: re_gaussian initialisation chosen, but no energy supplied" + write(*,*) " please set part_group_configs()%re_energy" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + if (part_group_configs(group_num)%re_std_energy .eq. 0) then + write(*,*) "ERROR: re_gaussian initialisation chosen, but re_std_energy = 0" + write(*,*) " please set part_group_configs()%re_std_energy" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + if (part_group_configs(group_num)%re_pitch .eq. 0) then + write(*,*) "ERROR: re_gaussian initialisation chosen, but re_pitch = 0" + write(*,*) " please set part_group_configs()%re_pitch" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + endif + if (part_group_configs(group_num)%n_particles_total .eq. 0.d0) then + write(*,*) "ERROR: n_particles_total = 0, this is how weights are set" + write(*,*) " please set part_group_configs()%n_particles_total" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + +end subroutine check_compatibility_rep + subroutine check_compatibility_epf(group_num) implicit none integer :: group_num diff --git a/reg_tests/testcases/particle_rep_600/input b/reg_tests/testcases/particle_rep_600/input index ade62c6d69..fdb952a802 100644 --- a/reg_tests/testcases/particle_rep_600/input +++ b/reg_tests/testcases/particle_rep_600/input @@ -113,9 +113,9 @@ part_group_configs(1)%type = 'particle_kinetic_relativistic' part_group_configs(1)%n_particles = 1e2 part_group_configs(1)%mass = 0.000548579870184805 - part_group_configs(1)%init_function = 'basic' + part_group_configs(1)%init_function = 're_gaussian' part_group_configs(1)%init_pdf = 'RZ' - part_group_configs(1)%num_re = 1.175403e17 + part_group_configs(1)%n_particles_total = 1.175403e17 part_group_configs(1)%re_energy = 9.8e4 part_group_configs(1)%re_std_energy = 5e2 part_group_configs(1)%re_pitch = -1 From 75c9d1af78b7e7e1af4c24dc3ec18066f17e484b Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 14:09:04 +0200 Subject: [PATCH 04/16] Safer mpi ierr handling in mod_initialise_particles --- particles/diagnostics/mod_neutral_density.f90 | 1 + particles/initialisers/mod_initialise_particles.f90 | 7 +++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/particles/diagnostics/mod_neutral_density.f90 b/particles/diagnostics/mod_neutral_density.f90 index 5eb1a56ee6..7d84f80685 100644 --- a/particles/diagnostics/mod_neutral_density.f90 +++ b/particles/diagnostics/mod_neutral_density.f90 @@ -4,6 +4,7 @@ module mod_neutral_density use mod_interp, only: mode_moivre use mod_basisfunctions use particle_tracer + use phys_module, only: n_part_groups !$ use omp_lib implicit none diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 3848a9b470..0f86a4852d 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -27,11 +27,13 @@ module mod_initialise_particles use mod_import_experimental_dist use phys_module, only: part_group_configs, type_part_group_config, n_part_groups use mod_particle_group_id, only: matching_part_config_indices - use mpi + use mpi, only: MPI_ABORT, MPI_COMM_WORLD implicit none - integer :: ierr !> mpi error code + private + + public :: initialise_particles_for_sim contains @@ -76,6 +78,7 @@ subroutine initialise_group(sim, group_num) type(type_part_group_config) :: config type(spatial_pdf) :: space_pdf logical :: exists + integer :: ierr config = part_group_configs(matching_part_config_indices(group_num)) From c760231cf223d89b51cc2c0ebcade92c093a76c9 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 14:41:11 +0200 Subject: [PATCH 05/16] Added extra compatibility checks before initialising particles --- particles/initialisers/mod_initialise_particles.f90 | 2 +- particles/mod_coupling_settings.f90 | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 0f86a4852d..d57124b024 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -49,7 +49,7 @@ subroutine initialise_particles_for_sim(sim) j = matching_part_config_indices(i) ! get the matching part_group_config index select case(trim(part_group_configs(j)%coupling_scheme)) - case('ics', 'ncs') + case('ics', 'ncs', 'non') if (sim%my_id == 0) then write(*,*) "Particle initialisation skipped for ics/ncs, particles are born elsewhere" endif diff --git a/particles/mod_coupling_settings.f90 b/particles/mod_coupling_settings.f90 index 81b7d2cce8..7f621d9750 100644 --- a/particles/mod_coupling_settings.f90 +++ b/particles/mod_coupling_settings.f90 @@ -184,6 +184,14 @@ subroutine check_compatibility_epf(group_num) endif !> Check initialisation parameters + if (trim(part_group_configs(group_num)%init_function) .eq. 'experimental') then + if (part_group_configs(group_num)%n_phi_planes .eq. 0) then + write(*,*) "ERROR: Maxwell initialisation chosen, but n_phi_planes = 0" + write(*,*) " needs to be at least 1, please set part_group_configs()%n_phi_planes" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + endif + if (trim(part_group_configs(group_num)%init_function) .eq. 'maxwell') then if (part_group_configs(group_num)%T_maxwell .eq. 0.d0) then write(*,*) "ERROR: Maxwell initialisation chosen, but no temperature supplied" @@ -196,6 +204,7 @@ subroutine check_compatibility_epf(group_num) call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif endif + if (part_group_configs(group_num)%n_particles_total .eq. 0.d0) then write(*,*) "ERROR: n_particles_total = 0, this is how weights are set" write(*,*) " please set part_group_configs()%n_particles_total" From 6af9ba70baffd527ac5d66613c4611023e882d25 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 14:41:38 +0200 Subject: [PATCH 06/16] Fixed handling outside LCFS for analytical space pdf in particle init --- particles/initialisers/mod_rej_f.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index 333739b2e1..a030f34e4a 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -191,8 +191,7 @@ pure function analytical_pdf(n, P, gradP) result(f) real*8 :: minor_r minor_r = sqrt((P(1) - ES%Z_axis)**2 + (P(2) - ES%R_axis)**2) - f = real((1.d0 - (minor_r/ES%LCFS_a)**2)**nu, 4) - f = max(f, 0.0e0) !> clamp outside LCFS + f = real(max(1.d0 - (minor_r/ES%LCFS_a)**2, 0.d0)**nu, 4) end function analytical_pdf !> Weight proportional to normalised toroidal current density j_tor From 88e37870200639ee3a6968e4336c05e815e9f561 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 17:36:42 +0200 Subject: [PATCH 07/16] Updated initialiser subroutines to exepct spatial_pdf type rather than rej function and variables seperately, evaluating the rej_f has its own subroutine now too --- particles/initialisers/initialisers_RE.f90 | 10 +- particles/initialisers/initialisers_base.f90 | 362 +++++++----------- .../initialisers/mod_initialise_particles.f90 | 28 +- particles/initialisers/mod_rej_f.f90 | 52 ++- 4 files changed, 172 insertions(+), 280 deletions(-) diff --git a/particles/initialisers/initialisers_RE.f90 b/particles/initialisers/initialisers_RE.f90 index eea5bedf0d..ea538a22e3 100644 --- a/particles/initialisers/initialisers_RE.f90 +++ b/particles/initialisers/initialisers_RE.f90 @@ -34,14 +34,8 @@ subroutine initialise_re_gaussian(sim, group_num, rng, space_pdf, energy, pitch, integer :: num_part real*8, allocatable :: ran_uniform(:), ran_gaussian(:) - if (present(space_pdf)) then - call initialise_particles(sim%groups(group_num)%particles, & - sim%fields%node_list, sim%fields%element_list, rng, & - variables=space_pdf%vars, transform=space_pdf%f) - else - call initialise_particles(sim%groups(group_num)%particles, & - sim%fields%node_list, sim%fields%element_list, rng) - endif + call initialise_particles(sim%groups(group_num)%particles, & + sim%fields%node_list, sim%fields%element_list, rng, space_pdf) num_part = size(sim%groups(group_num)%particles,1) diff --git a/particles/initialisers/initialisers_base.f90 b/particles/initialisers/initialisers_base.f90 index 10fcc5a42c..ad83e5919c 100644 --- a/particles/initialisers/initialisers_base.f90 +++ b/particles/initialisers/initialisers_base.f90 @@ -9,7 +9,7 @@ module initialisers_base use mod_rej_f implicit none private - public initialise_particles, no_transform, adjust_particle_weights + public initialise_particles, adjust_particle_weights, eval_rej_f public set_velocity_from_T, domain_bounding_box, initialise_particles_H_mu_psi public initialise_particles_H_mu_psi_phiplanes public set_particle_weights_canonical_maxwellian, normalize_with_projection @@ -79,11 +79,82 @@ subroutine part_inout_s(p_inout,n_x,x,time,fields,n_real_param,& end interface contains +!> Evaluate the variables a spatial rejection function needs at one sample point. +!> +!> Fills P(:), and gradP(:,:) when needs_grad, in the order given by variables(:), +!> following the index convenction documented in mod_rej_f: +!> >0: JOREK variable number +!> 0: constant 1, -1: R, -2: Z, -3: Phi +!> variables(:) must be sorted ascending, so the n_geom non-positive entries +!> come first and the n_mhd JOREK variables last. +!> +!> gradP is (d/dR, d/dZ, d/dphi) and is set to zero when needs_grad is .false. +pure function eval_rej_vars(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf) result(f) + + !> i/o vars + type(type_node_list), intent(in) :: node_list + type(tyep_element_list), intent(in) :: element_list + integer, intent(in) :: i_elm + real*8, intent(in) :: s, t !> local element coordinates + real*8, intent(in) :: R, Z, phi !> global cylindrical coordinates + type(spatial_pdf), intent(in) :: space_pdf + real*4 :: f !> acceptance probability + + !> internal vars + real*8, parameter :: EPS_JAC = 1.d-12 !> below this element is degenerate + integer :: n_geom, n_mhd, k + real*8 :: P(size(space_pdf%vars)), gradP(3,size(space_pdf%vars)) + real*8 :: P_s(count(space_pdf%vars > 0)), P_t(count(space_pdf%vars > 0)) + real*8 :: P_phi(count(space_pdf%vars > 0)) + real*8 :: R_i, R_S, R_t, Z_i, Z_s, Z_t, xjac, inv_xjac + + n_mhd = count(space_pdf%vars > 0) + n_geom = size(space_pdf%vars) - n_mhd + gradP = 0.d0 + + !> Geomteric vars + do k = 1, n_geom + select case (space_pdf%vars(k)) + case(0); P(k) = 1.d0 + case(-1); P(k) = R ; if (needs_grad) gradP(1,k) = 1.d0 + case(-2); P(k) = Z ; if (needs_grad) gradP(2,k) = 1.d0 + case(-1); P(k) = phi ; if (needs_grad) gradP(3,k) = 1.d0 + endselect + enddo + + !> MHD Variables : values only unless the rej_f asked for gradients + !> interpolated within the element + if (n_mhd .ge. 1) then + if (space_pdf%needs_grad) then + call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & + s, t, phi, P(n_geom+1:), P_s, P_t, P_phi, R_i, R_s, R_t, Z_i, Z_s, Z_t) + + !> interp_PRZ gives s,t derivatives - transform to R,Z + !> degenerate element (e.g axis point) -> leave gradient zero + xjac = R_s*Z_t - R_t*Z_s + inv_xjac = 0.d0 + if (abs(xjac) > EPS_JAC) inv_xjac = 1.d0/xjac + + do k = 1, n_mhd + gradP(1,n_geom+k) = ( Z_t*P_s(k) - Z_s*P_t(k)) * inv_xjac + gradP(2,n_geom+k) = (-R_t*P_s(k) + R_s*P_t(k)) * inv_xjac + gradP(3,n_geom+k) = P_phi(k) + enddo + else + call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & + s, t, phi, P(n_geom+1), R_i, Z_i) + endif + endif + + f = space_pdf%f(size(space_pdf%vars), P, gradP) +end function eval_rej_vars + + !> Set positions for particles by rejection sampling from geometric and mhd !> variables after collecting with transform, within Rbound, Zbound and Phibound !> if present. See [[test_rejection_sampling]] for examples. subroutine initialise_particles(particles, node_list, element_list, & - rng, variables, transform, f, Rbound, Zbound, Phibound, & + rng, space_pdf, Rbound, Zbound, Phibound, & rng_n_streams_round_off_in) use mpi use mod_sampling @@ -96,30 +167,22 @@ subroutine initialise_particles(particles, node_list, element_list, & type(type_node_list), intent(in) :: node_list type(type_element_list), intent(in) :: element_list class(type_rng), intent(in) :: rng !< What type of random number generator to use. Is re-seeded in the subroutine. - integer, dimension(:), intent(in), optional :: variables !< Which variables from JOREK to use. If absent, sample uniformly. - procedure(rej_f), optional :: transform !< Merge variables into a single criterium between 0 and 1 for rej. sampling - !< Special values: 0 = 1, -1 = R, -2 = Z, -3 = Phi. Must be in ascending order! - real*8, intent(in), optional :: f !< Weighting factor: f=0 indicates uniform weights, f=1 indicates uniform distribution - !< (particle weight proportional to transform(P) at that point.) If omitted take f=0. + type(spatial_pdf), intent(in), optional :: space_pdf real*8, dimension(2), intent(in), optional :: Rbound, Zbound, Phibound !< Between which coordinates to sample (RZPhi). !< if omitted, determine automatically from node_list logical, intent(in), optional :: rng_n_streams_round_off_in !< round-off the rng n_streams at 2**ceil ! Internal variables real*8 :: R, Z, phi, s, t, DUMMY_REAL - real*8 :: R_i, R_s, R_t, Z_i, Z_s, Z_t, xjac real*8 :: Rbox(2), Zbox(2), Phibox(2) - integer :: i, j, k, ifail + integer :: i, j, ifail real*8 :: ran(7) integer :: i_elm real*8 :: t0, t1, ostart, oend integer :: seq, n_streams, n_threads, i_thread - integer :: n_geom, n_mhd integer :: my_id, n_mpi integer :: seed - logical :: rng_n_streams_round_off - real*8, dimension(:), allocatable :: P - real*8, dimension(:,:), allocatable :: gradP + logical :: rng_n_streams_round_off, use_rej class(type_rng), allocatable, dimension(:) :: rngs ! The RNGs for all the threads integer, dimension(:), allocatable :: i_to_find logical, dimension(:), allocatable :: not_found @@ -131,21 +194,12 @@ subroutine initialise_particles(particles, node_list, element_list, & call MPI_COMM_SIZE(MPI_COMM_WORLD, n_mpi, ifail) rng_n_streams_round_off = .false. if(present(rng_n_streams_round_off_in)) rng_n_streams_round_off = rng_n_streams_round_off_in - if (present(variables)) then - if (.not. present(transform)) then - write(*,*) "ERROR: if variables are present in set_particle_position_rejection_sampling transform must also be present" - call MPI_ABORT(MPI_COMM_WORLD, 10, ifail) - end if - ! Get the number of mhd variables to use - allocate(P(size(variables,1))) - allocate(gradP(3, size(variables,1))) - gradP = 0.d0 - n_mhd = count(variables .gt. 0) - n_geom = size(variables, 1) - n_mhd - else - n_mhd = 0 - n_geom = 0 - end if + + use_rej = .false. + if (present(space_pdf)) then + call check_spatial_pdf(space_pdf, "initialise_particles") + use_rej = associated(space_pdf%f) + endif ! Setup bounding boxes call domain_bounding_box(node_list, element_list, Rbox(1), Rbox(2), Zbox(1), Zbox(2)) @@ -213,10 +267,9 @@ subroutine initialise_particles(particles, node_list, element_list, & #else !$omp parallel default(none) & #endif - !$omp shared(particles, node_list, element_list, Rbox, Zbox, PhiBox, variables, & - !$omp rngs, n_threads, n_streams, seed, my_id, n_mhd, n_geom, i_to_find, not_found) & - !$omp private(j, i, R, Z, phi, i_elm, s, t, ifail, seq, ran, i_thread, P, DUMMY_REAL, gradP, & - !$omp R_i, R_s, R_t, Z_i, Z_s, Z_t, xjac) + !$omp shared(particles, node_list, element_list, Rbox, Zbox, PhiBox, space_pdf, use_rej, & + !$omp rngs, n_threads, n_streams, seed, my_id, i_to_find, not_found) & + !$omp private(j, i, R, Z, phi, i_elm, s, t, ifail, seq, ran, i_thread, DUMMY_REAL) i_thread = 0 !$ i_thread=omp_get_thread_num() !$omp do schedule(static) @@ -228,53 +281,18 @@ subroutine initialise_particles(particles, node_list, element_list, & call find_RZ(node_list,element_list,R,Z,DUMMY_REAL,DUMMY_REAL,i_elm,s,t,ifail) if (ifail .eq. 0) then - if (present(variables)) then - ! Select the mhd variables requested - if (n_mhd .ge. 1) then - call interp_PRZ(node_list, element_list, i_elm, variables(n_geom+1:n_geom+n_mhd), & - n_mhd, s, t, phi, P(n_geom+1:n_geom+n_mhd), gradP(1, n_geom+1:n_geom+n_mhd), & - gradP(2, n_geom+1:n_geom+n_mhd), gradP(3, n_geom+1:n_geom+n_mhd), R_i, R_s, R_t, & - Z_i, Z_s, Z_t) - - !> Calculate R/Z derivs for gradP if needed by the rej_f (interp_PRZ gives s, t derivs) - xjac = R_s*Z_t - R_t*Z_s - do k = 1, n_mhd - gradP(1:2,n_geom+k) = [Z_t*gradP(1,n_geom+k) - Z_s*gradP(2,n_geom+k), & - R_s*gradP(2,n_geom+k) - R_t*gradP(1,n_geom+k)]/xjac - enddo - end if - - do k=1,n_geom - select case (variables(k)) - case (0); P(k) = 1.d0 ; gradP(:,k) = 0.d0 ! 0 - case (-1); P(k) = R ; gradP(:,k) = [1.d0, 0.d0, 0.d0] ! R - case (-2); P(k) = Z ; gradP(:,k) = [0.d0, 1.d0, 0.d0] ! Z - case (-3); P(k) = phi ; gradP(:,k) = [0.d0, 0.d0, 1.d0] ! phi - end select - end do - - if (present(transform)) then - if (ran(4) .lt. transform(size(p,1), p, gradP)) then - particles(j)%x = [R, Z, phi] - particles(j)%i_elm = i_elm - particles(j)%st = [s, t] - select type (pa => particles(j)) - type is (particle_kinetic_leapfrog) - pa%v = ran(5:7) ! save other components of this point for velocity init in a later routine - end select - not_found(i) = .false. - end if - end if - else - particles(j)%x = [r, z, phi] - particles(j)%i_elm = i_elm - particles(j)%st = [s, t] - select type (pa => particles(j)) - type is (particle_kinetic_leapfrog) - pa%v = ran(5:7) ! save other components of this point for velocity init in a later routine - end select - not_found(i) = .false. + !> Reject this position if it fails the spatial pdf + if (use_rej) then + if (ran(4) .ge. eval_rej_f(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf)) cycle end if + particles(j)%x = [R, z, phi] + particles(j)%i_elm = i_elm + particles(j)%st = [s, t] + select type(pa => particles(j)) + type is (particle_kinetic_leapfrog) + pa%v = ran(5:7) !> Save other components of this point for velocity init in a later routine + endselect + not_found(i) = .false. end if enddo !$omp end do @@ -285,11 +303,6 @@ subroutine initialise_particles(particles, node_list, element_list, & not_found = .true. end do - if (present(variables)) then - deallocate(P) - deallocate(gradP) - endif - call cpu_time(t1) !$ oend = omp_get_wtime() write(*,'(i5,A,2f12.4)') my_id, ' Time particle initialize cpu/wall :',t1-t0, oend-ostart @@ -615,7 +628,7 @@ end function rejection_funct_gpdf !> Set Psi_transform to transform from [0,1] to your desired range subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_maxwell, & Theta_transform, Psi_transform, alpha, E_max, include_vpar, uniform_space, & - uniform_space_rej_f, uniform_space_rej_vars, cor, charge, rng_n_streams_round_off_in) + space_pdf, cor, charge, rng_n_streams_round_off_in) use mod_rng use mod_fields use mod_random_seed @@ -638,10 +651,7 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma real*8, intent(in), optional :: alpha !< Make more fast (>0) or slow (<0) particles and weigh them appropriately real*8, intent(in), optional :: E_max !< If alpha=1 we select particles from a block-distribution, up to E_max logical, intent(in), optional :: include_vpar !< Initialize particles with local parallel velocity - logical, intent(in), optional :: uniform_space !< Do not - !< use {psi,theta}_transform if present but use rejection sampling in RZ - procedure(rej_f), optional :: uniform_space_rej_f !< Merge variables into a single criterium between 0 and 1 for rej. sampling - !< Special values: 0 = 1, -1 = R, -2 = Z, -3 = Phi. Must be in ascending order! + type(spatial_pdf), intent(in), optional :: space_pdf !< spatial pdf to reject against when uniform_space is set. If absent the RZ sampling is uniform integer, dimension(:), intent(in), optional :: uniform_space_rej_vars !< Variables to use for uniform_space_rej_f type(coronal), intent(in), optional :: cor !< Coronal equilibrium datatype for this particle. If unset, do not alter q integer, intent(in), optional :: charge !< Use this if cor is not present @@ -663,12 +673,10 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma real*8, dimension(1) :: P, P_s, P_t, P_phi #endif - real*8, dimension(:), allocatable :: P2 - real*8, dimension(:,:), allocatable :: grad_P2 - real*8 :: R_s, R_t, Z_s, Z_t, R_i, Z_i, xjac + real*8 :: R_s, R_t, Z_s, Z_t real*8 :: s, t, u_init_max, temp, u real*8 :: psi_axis, R_axis, Z_axis, s_axis, t_axis - integer :: i_elm, i, j, k, ifail, my_id, n_mpi, ierr, n_mhd, n_geom + integer :: i_elm, i, j, ifail, my_id, n_mpi, ierr real*8, dimension(fields%element_list%n_elements,2) :: psi_minmax_list real*8, allocatable, dimension(:,:) :: rans class(particle_base), dimension(:), allocatable :: particles_tmp @@ -677,7 +685,7 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma real*8 :: Rbox(2), Zbox(2), DUMMY_R, DUMMY_Z integer :: blocksize, prev_blocksize, particles_to_do_local, particles_done_local integer :: to_find, n_tries_now, n_found - logical :: all_done, init_uniform_space, my_include_vpar, rng_n_streams_round_off + logical :: all_done, init_uniform_space, my_include_vpar, rng_n_streams_round_off, use_rej real*8 :: my_alpha rng_n_streams_round_off = .false. @@ -697,25 +705,11 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma my_include_vpar = .false. end if - if (present(uniform_space_rej_f)) then - - if (.not. present(uniform_space_rej_vars)) then - write(*,*) "ERROR: if sampling function f is present variables must be given" - call MPI_ABORT(MPI_COMM_WORLD, 10, ifail) - end if - - ! Get the number of mhd variables to use - allocate(P2(size(uniform_space_rej_vars,1))) - - n_mhd = count(uniform_space_rej_vars .gt. 0) - n_geom = size(uniform_space_rej_vars, 1) - n_mhd - - allocate(grad_P2(3,size(uniform_space_rej_vars,1))) - - else - n_mhd = 0 - n_geom = 0 - end if + use_rej = .false. + if (present(space_pdf)) then + call check_spatial_pdf(space_pdf, "initialise_particles_H_mu_psi") + use_rej = associated(space_pdf%f) + endif if (present(alpha)) then write(*,*) "alpha not implemented yet" @@ -821,10 +815,10 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma !$omp parallel do default(none) & !$omp shared(particles_tmp, psimax, psimin, found, F0, cor, mass, charge, T_Maxwell, & !$omp fields, psi_minmax_list, rans, R_axis, Z_axis, blocksize, & - !$omp my_include_vpar, central_density, init_uniform_space, Rbox, Zbox, uniform_space_rej_vars, n_geom, n_mhd) & + !$omp my_include_vpar, central_density, init_uniform_space, Rbox, Zbox, space_pdf, use_rej) & #endif - !$omp private(i, psi, theta, phi, i_elm, s, t, R, Z, R_s, R_t, Z_s, Z_t, P2, & - !$omp R_i, Z_i, xjac, grad_P2, u, particle_kinetic_tmp, v2, v_par, & + !$omp private(i, psi, theta, phi, i_elm, s, t, R, Z, R_s, R_t, Z_s, Z_t, & + !$omp u, particle_kinetic_tmp, v2, v_par, & #ifdef fullmhd !$omp A3, AR, AZ, A3_R, A3_Z, AR_Z, AR_p, AZ_R, AZ_P, Fprof, & #endif @@ -838,44 +832,9 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma call transform_uniform_cylindrical([ran(3),ran(4),ran(5)], Rbox, Zbox, [0.d0,TWOPI], R, Z, phi) call find_RZ(fields%node_list, fields%element_list,R,Z,DUMMY_R,DUMMY_Z,i_elm,s,t,ifail) - if (present(uniform_space_rej_f) .and. i_elm .gt. 0) then - - do k=1,n_geom - select case (uniform_space_rej_vars(k)) - case (0); P2(k) = 1.d0; - case (-1); P2(k) = R; - case (-2); P2(k) = Z; - case (-3); P2(k) = phi; - end select - end do - - do k=1,n_geom - select case (uniform_space_rej_vars(k)) - case (0); grad_P2(:,k) = 0.d0; ! 0 - case (-1); grad_P2(:,k) = [1.d0,0.d0,0.d0]; ! R - case (-2); grad_P2(:,k) = [0.d0,1.d0,0.d0]; ! Z - case (-3); grad_P2(:,k) = [0.d0,0.d0,1.d0]; ! phi - end select - end do - - if (n_mhd .ge. 1) then - - call interp_PRZ(fields%node_list, fields%element_list,i_elm, & - uniform_space_rej_vars(n_geom+1:n_geom+n_mhd),n_mhd,s,t,phi, & - P2(n_geom+1:n_geom+n_mhd), grad_P2(1,n_geom+1:n_geom+n_mhd), & - grad_P2(2,n_geom+1:n_geom+n_mhd), grad_P2(3,n_geom+1:n_geom+n_mhd), & - R_i, R_s, R_t, Z_i, Z_s, Z_t) - - xjac = R_s*Z_t - R_t*Z_s - - do k=1,n_mhd - grad_P2(1:2,n_geom+k) = [Z_t * grad_P2(1,n_geom+k) - Z_s * grad_P2(2,n_geom+k), & - -R_t * grad_P2(1,n_geom+k) + R_s * grad_P2(2,n_geom+k)]/xjac - end do - - end if - - if (uniform_space_rej_f(size(uniform_space_rej_vars), P2, grad_P2) .lt. ran(7)) i_elm = 0 + !> Reject this position if it fails the spatial pdf + if (use_rej .and. i_elm .gt. 0) then + if (eval_rej_f(fields%node_list, fields%element_list, i_elm, s, t, phi, R, Z, space_pdf) .lt. ran(7)) i_elm = 0 end if else @@ -1053,7 +1012,7 @@ end subroutine initialise_particles_H_mu_psi !! multiple marker particles per guiding centre particle. subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, mass, T_maxwell, & Theta_transform, Psi_transform, alpha, E_max, include_vpar, uniform_space, & - uniform_space_rej_f, uniform_space_rej_vars, cor, charge, n_phi_planes_in, & + space_pdf, cor, charge, n_phi_planes_in, & n_gyro_orbit_in, rng_n_streams_round_off_in) use mod_rng use mod_fields @@ -1079,9 +1038,7 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, logical, intent(in), optional :: include_vpar !< Initialize particles with local parallel velocity logical, intent(in), optional :: uniform_space !< Do not !< use {psi,theta}_transform if present bute rejection sampling in RZ - procedure(rej_f), optional :: uniform_space_rej_f !< Merge variables into a single criterium between 0 and 1 for rej. sampling - !< Special values: 0 = 1, -1 = R, -2 = Z, - Phi. Must be in ascending order! - integer, dimension(:), intent(in), optional :: uniform_space_rej_vars !< Variables to use for uniform_space_rej_f + type(spatial_pdf), intent(in), optional :: space_pdf !< spatial pdf to reject against when uniform_space is set. if absent the RZ sampling is uniform type(coronal), intent(in), optional :: cor !< Coronal equilibrium datatype for this particle. If unset, do not alter q integer, intent(in), optional :: charge !< Use this if cor is not present real*8, intent(in), optional :: T_Maxwell !< constant Maxwellian temperature [eV] @@ -1103,12 +1060,10 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, #else real*8, dimension(1) :: P, P_s, P_t, P_phi #endif - real*8, dimension(:), allocatable :: P2 - real*8, dimension(:,:), allocatable :: grad_P2 - real*8 :: R_s, R_t, Z_s, Z_t, R_i, Z_i, xjac + real*8 :: R_s, R_t, Z_s, Z_t real*8 :: s, t, u_init_max, temp, u, v2, v_par real*8 :: psi_axis, R_axis, Z_axis, s_axis, t_axis - integer :: i_elm, i, j, k, ifail, my_id, n_mpi, ierr, n_mhd, n_geom + integer :: i_elm, i, j, ifail, my_id, n_mpi, ierr real*8, dimension(fields%element_list%n_elements,2) :: psi_minmax_list real*8, allocatable, dimension(:,:) :: rans class(particle_base), dimension(:), allocatable :: particles_tmp @@ -1118,7 +1073,7 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, integer :: blocksize, prev_blocksize, particles_to_do_local, particles_done_local, blocksize_tmp integer :: to_find, n_tries_now, n_found logical :: all_done, init_uniform_space, my_include_vpar - logical :: init_phiplanes, init_gyro_orbit, rng_n_streams_round_off + logical :: init_phiplanes, init_gyro_orbit, rng_n_streams_round_off, use_rej real*8 :: my_alpha integer :: n_phi_planes,i_phi_planes, n_gyro_orbit, i_gyro_temp, i_gyro_orbit @@ -1158,24 +1113,11 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, n_gyro_orbit=1 !In this case, the loops aren't there endif - if (present(uniform_space_rej_f)) then - if (.not. present(uniform_space_rej_vars)) then - write(*,*) "ERROR: if sampling function f is present variables must be given" - call MPI_ABORT(MPI_COMM_WORLD, 10, ifail) - end if - - ! Get the number of mhd variables to use - allocate(P2(size(uniform_space_rej_vars,1))) - - n_mhd = count(uniform_space_rej_vars .gt. 0) - n_geom = size(uniform_space_rej_vars, 1) - n_mhd - - allocate(grad_P2(3,size(uniform_space_rej_vars,1))) - - else - n_mhd = 0 - n_geom = 0 - end if + use_rej = .false. + if (present(space_pdf)) then + call check_spatial_pdf(space_pdf, "initialise_particles_H_mu_psi_phiplanes") + use_rej = associated(space_pdf%f) + endif if (present(alpha)) then write(*,*) "alpha not implemented yet" @@ -1290,10 +1232,10 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, !$omp parallel do default(none) & !$omp shared(particles_tmp, psimax, psimin, found, F0, cor, mass, charge, T_Maxwell, & !$omp fields, psi_minmax_list, rans, R_axis, Z_axis, blocksize, init_phiplanes,init_gyro_orbit, n_gyro_orbit,blocksize_tmp,& - !$omp my_include_vpar, central_density, init_uniform_space, Rbox, Zbox, uniform_space_rej_vars, n_geom, n_mhd,n_phi_planes,my_id) & + !$omp my_include_vpar, central_density, init_uniform_space, Rbox, Zbox, space_pdf, use_rej, n_phi_planes,my_id) & #endif - !$omp private(i, psi, theta, phi, i_elm, s, t, R, Z, R_s, R_t, Z_s, Z_t, P2, & - !$omp R_i, Z_i, xjac, grad_P2, u, particle_kinetic_tmp, v2, v_par, & + !$omp private(i, psi, theta, phi, i_elm, s, t, R, Z, R_s, R_t, Z_s, Z_t, & + !$omp u, particle_kinetic_tmp, v2, v_par, & #ifdef fullmhd !$omp A3, AR, AZ, A3_R, A3_Z, AR_Z, AR_p, AZ_R, AZ_P, Fprof, & #endif @@ -1315,44 +1257,9 @@ subroutine initialise_particles_H_mu_psi_phiplanes(particles, fields, rng_base, call find_RZ(fields%node_list, fields%element_list,R,Z,DUMMY_R,DUMMY_Z,i_elm,s,t,ifail) - if (present(uniform_space_rej_f) .and. i_elm .gt. 0) then - - do k=1,n_geom - select case (uniform_space_rej_vars(k)) - case (0); P2(k) = 1.d0; - case (-1); P2(k) = R; - case (-2); P2(k) = Z; - case (-3); P2(k) = phi; - end select - end do - - do k=1,n_geom - select case (uniform_space_rej_vars(k)) - case (0); grad_P2(:,k) = 0.d0; ! 0 - case (-1); grad_P2(:,k) = [1.d0,0.d0,0.d0]; ! R - case (-2); grad_P2(:,k) = [0.d0,1.d0,0.d0]; ! Z - case (-3); grad_P2(:,k) = [0.d0,0.d0,1.d0]; ! phi - end select - end do - - if (n_mhd .ge. 1) then - - call interp_PRZ(fields%node_list, fields%element_list,i_elm, & - uniform_space_rej_vars(n_geom+1:n_geom+n_mhd),n_mhd,s,t,phi, & - P2(n_geom+1:n_geom+n_mhd), grad_P2(1,n_geom+1:n_geom+n_mhd), & - grad_P2(2,n_geom+1:n_geom+n_mhd), grad_P2(3,n_geom+1:n_geom+n_mhd), & - R_i, R_s, R_t, Z_i, Z_s, Z_t) - - xjac = R_s*Z_t - R_t*Z_s - - do k=1,n_mhd - grad_P2(1:2,n_geom+k) = [Z_t * grad_P2(1,n_geom+k) - Z_s * grad_P2(2,n_geom+k), & - -R_t * grad_P2(1,n_geom+k) + R_s * grad_P2(2,n_geom+k)]/xjac - end do - - end if - - if (uniform_space_rej_f(size(uniform_space_rej_vars), P2, grad_P2) .lt. ran(7)) i_elm = 0 + !> reject this position if it fails the spatial pdf + if (use_rej .and. i_elm .gt. 0) then + if (eval_rej_f(fields%node_list, fields%element_list, i_elm, s, t, phi, R, Z, space_pdf) .lt. ran(7)) i_elm = 0 end if else @@ -1778,19 +1685,6 @@ subroutine domain_bounding_box(node_list, element_list, Rmin, Rmax, Zmin, Zmax) end subroutine domain_bounding_box -!> Dummy function to use when no transform is desired. Copies the first parameter -!> into the output or sets out to 1. -pure function no_transform(in) result(out) - real*8, dimension(:), intent(in) :: in - real*8 :: out - if (size(in,1) .gt. 0) then - out = in(1) - else - out = 1.d0 - end if -end function no_transform - - !> Adjust weights on all particles to have the correct number of atoms in total subroutine adjust_particle_weights(particles, num_atoms_total) use mpi diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index d57124b024..6e5e9d1128 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -111,20 +111,11 @@ subroutine initialise_group(sim, group_num) write(*, '(A,I0)') " n_phi_planes : ", config%n_phi_planes endif - if (associated(space_pdf%f)) then - call initialise_particles_H_mu_psi_phiplanes( & - sim%groups(group_num)%particles, sim%fields, pcg32_rng(), & - sim%groups(group_num)%mass, uniform_space=.true., & - uniform_space_rej_f=space_pdf%f, & - uniform_space_rej_vars=space_pdf%vars, charge=1, & - T_maxwell=config%T_maxwell, n_phi_planes_in=config%n_phi_planes) - else - call initialise_particles_H_mu_psi_phiplanes( & - sim%groups(group_num)%particles, sim%fields, pcg32_rng(), & - sim%groups(group_num)%mass, uniform_space=.true., & - charge=1, T_maxwell=config%T_maxwell, & - n_phi_planes_in=config%n_phi_planes) - endif + call initialise_particles_H_mu_psi_phiplanes( & + sim%groups(group_num)%particles, sim%fields, pcg32_rng(), & + sim%groups(group_num)%mass, uniform_space=.true., & + space_pdf=space_pdf, charge=1, & + T_maxwell=config%T_maxwell, n_phi_planes_in=config%n_phi_planes) !> EPs : experimental distribution - f(R,Z,energy,pitch), expects a .h5 file - see mod_import_experimental_dist case ('experimental') @@ -150,13 +141,8 @@ subroutine initialise_group(sim, group_num) write(*,'(A,ES12.3)') " Pitch : ", config%re_pitch endif - if (associated(space_pdf%f)) then - call initialise_re_gaussian(sim, group_num, pcg32_rng(), & - space_pdf, config%re_energy, config%re_pitch, config%re_std_energy) - else - call initialise_re_gaussian(sim, group_num, pcg32_rng(), & - energy=config%re_energy, pitch=config%re_pitch, std_energy=config%re_std_energy) - endif + call initialise_re_gaussian(sim, group_num, pcg32_rng(), & + space_pdf, config%re_energy, config%re_pitch, config%re_std_energy) !> Unknown init functions case default diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index a030f34e4a..3dff352894 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -23,13 +23,12 @@ module mod_rej_f use mpi implicit none - integer :: ierr ! mpi error code private public :: spatial_pdf public :: rej_f public :: itpa_tae_pdf, RZ_pdf, analytical_pdf, current_pdf - public :: spatial_pdf_from_name + public :: spatial_pdf_from_name, check_spatial_pdf ! ============================================================================= ! Interface @@ -39,10 +38,13 @@ module mod_rej_f !> !> @param n Number of field values in P and gradP !> @param P Field values at the sample point, ordered as vars(:) - !> @param gradP Gradients (3,n); may be unused + !> @param gradP (3,n) = (d/dR, d/dZ, d/dphi) of each entry of P. + !> Only computed when needs_grad = .true.; zero otherwise. !> @return Acceptance probability in [0,1] + !> + !> Must be pure : it is called from inside the OpenMP sampling loops. abstract interface - function rej_f(n, P, gradP) + pure function rej_f(n, P, gradP) implicit none integer, intent(in) :: n real*8, dimension(n), intent(in) :: P @@ -65,11 +67,13 @@ end function rej_f !> Example (custom profile) !> !> type(spatial_pdf) :: pdf - !> pdf%f => my_rej_function - !> pdf%vars = [var_psi, -1] ! psi and R + !> pdf%f => my_rej_function + !> pdf%vars = [var_psi, -1] ! psi and R + !> pdf%needs_grad = .true. ! only if my_rej_function requires gradP type :: spatial_pdf procedure(rej_f), nopass, pointer :: f => null() integer, allocatable :: vars(:) + logical :: needs_grad = .false. !> Does f require gradients of any vars? end type spatial_pdf @@ -86,23 +90,35 @@ end function rej_f !> 'analytical' - weight by (1-(r/a)^2)^nu !> 'current' - weight by normalised j_tor !> 'none' - no rejection (f => null, vars unallocated) - function spatial_pdf_from_name(name) result(pdf) - character(len=*), intent(in) :: name - type(spatial_pdf) :: pdf + !> + !> node_list/element_list are the equilibrium the pdf is calibrated against. + !> only 'current' uses them at the moment, but any profile normalised to + !> equilibrium fields will need them + function spatial_pdf_from_name(name, node_list, element_list) result(pdf) + character(len=*), intent(in) :: name + type(type_node_list), intent(in) :: node_list + type(type_element_list), intent(in) :: element_list + type(spatial_pdf) :: pdf + + integer :: ierr select case(trim(name)) case('itpa_tae') - pdf%f => itpa_tae_pdf - pdf%vars = [var_psi] + pdf%f => itpa_tae_pdf + pdf%vars = [var_psi] + pdf%needs_grad = .false. case('RZ') - pdf%f => RZ_pdf - pdf%vars = [-2, -1] + pdf%f => RZ_pdf + pdf%vars = [-2, -1] + pdf%needs_grad = .false. case('analytical') - pdf%f => analytical_pdf - pdf%vars = [-2, -1] + pdf%f => analytical_pdf + pdf%vars = [-2, -1] + pdf%needs_grad = .false. + case('current') !> NOTE : var_zj = 0 in fullMHD (j_tor is not a stored variable) @@ -112,8 +128,10 @@ function spatial_pdf_from_name(name) result(pdf) write(*,*) " with var_zj > 0. j_tor is not a stored variable in this model" call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif - pdf%f => current_pdf - pdf%vars = [-1, var_zj] + pdf%f => current_pdf + pdf%vars = [-1, var_zj] + pdf%needs_grad = .false. + case('none') !> Leave f => null() and vars unallocated From 94f538e60803b13c5ba2c6671fef1583bb8e724b Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Thu, 20 Aug 2026 17:37:09 +0200 Subject: [PATCH 08/16] Added calculation of j_tor/R bounds for current_pdf --- particles/initialisers/mod_rej_f.f90 | 82 ++++++++++++++++++++++++++-- 1 file changed, 77 insertions(+), 5 deletions(-) diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index 3dff352894..b0195f51d7 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -23,6 +23,12 @@ module mod_rej_f use mpi implicit none + + !> Normalisation range for current_pdf, in units of j_tor/R + !> Set from the equilibrium by calibrate_current_pdf() when the 'current' + !> pdf is constructed - never use current_pdf without that call. + real*8 :: jzmin = 0.d0, jzmax = 0.d0 + private public :: spatial_pdf @@ -128,6 +134,7 @@ function spatial_pdf_from_name(name, node_list, element_list) result(pdf) write(*,*) " with var_zj > 0. j_tor is not a stored variable in this model" call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif + call calibrate_current_pdf(node_list, element_list) pdf%f => current_pdf pdf%vars = [-1, var_zj] pdf%needs_grad = .false. @@ -146,6 +153,75 @@ function spatial_pdf_from_name(name, node_list, element_list) result(pdf) end select end function spatial_pdf_from_name + !> Set the j_tor/R normalisation range used by current_pdf, once per run. + !> + !> find_variable_mimax returns the extrema of one variable over an element's edges + !> j_tor and R are bounded independently and combined afterwardds, which can only widen the range. + subroutine calibrate_current_pdf(node_list, element_list) + use mod_newton_methods, only: find_variable_minmax + + !> i/o vars + type(type_node_list), intent(in) :: node_list + type(type_element_list), intent(in) :: element_list + + !> internal vars + real*8 :: zj_lo, zj_hi, R_lo, R_hi, e_lo, e_hi + integer :: i_elm, my_id, ierr + + zj_lo = 1.d10; zj_hi = -1.d10 + R_lo = 1.d10; R_hi = -1.d10 + + do i_elm = 1, element_list%n_elements + call find_variable_minmax(node_list, element_list, i_elm, var_zj, e_lo, e_hi) + zj_lo = min(zj_lo, e_lo); zj_hi = max(zj_hi, e_hi) + call find_variable_minmax(node_list, element_list, i_elm, -1, e_lo, e_hi) + R_lo = min(R_lo, e_lo); R_hi = max(R_hi, e_hi) + enddo + + !> R>0 on any grid, so j_tor/R is extremal at an R endpoint; + jzmin = min(zj_lo/R_lo, zj_lo/R_hi) + jzmax = max(zj_hi/R_lo, zj_hi/R_hi) + + call MPI_COMM_RANK(MPI_COMM_WORLD, my_id, ierr) + + if (jzmax .le. jzmin) then + if (my_id == 0) then + write(*,*) "ERROR (mod_rej_f): 'current' spatial PDF found no spread in j_tor/R" + write(*,*) " j_tor range ", zj_lo, zj_hi + write(*,*) " R range ", R_lo, R_hi + endif + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + + if (my_id == 0) then + write(*,'(A,2ES12.3)') " j_tor range : ", zj_lo, zj_hi + write(*,'(A,2Es12.3)') " j_tor/R range : ", jzmin, jzmax + endif + end subroutine calibrate_current_pdf + + !> Check that a spatial_pdf satisfies the invariants the samplers rely on. + !> Does nothing for an inavtive pdf (f => null(), ie init_pdf = 'none') + subroutine check_spatial_pdf(pdf, caller) + type(spatial_pdf), intent(in) :: pdf + character(len=*), intent(in) :: caller !> name to quote in error message + + integer :: ierr + + if (.not. associated(pdf%f)) return + + if (.not. allocated(pdf%vars)) then + write(*,*) "ERROR (", caller, "): spatial_pdf has a rejection function but no vars(:)" + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + + !> Ascendinf order puts the geometric entries (<= 0) first, which is how + !> the samplers split vars(:) into geometric and JOREK variables + if (any(pdf%vars(2:) < pdf%vars(:size(pdf%vars)-1))) then + write(*,*) "ERROR (", caller, "): spatial_pdf vars(:) must be in ascending order, got ", pdf%vars + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) + endif + end subroutine check_spatial_pdf + ! ============================================================================= ! Rejection functions @@ -226,11 +302,7 @@ pure function current_pdf(n, P, gradP) result(f) real*8, dimension(3,n), intent(in) :: gradP real*4 :: f - !> TODO: make jzmin.jzmax namelist params - real*8, parameter :: jzmax = 3.0d0 / 10.0d0 - real*8, parameter :: jzmin = 1.239d-4 / 11.0d0 - f = real((P(2)/P(1) - jzmin) / (jzmax - jzmin), 4) - f = max(f, 0.0e0) + f = min(max(f, 0.0e0), 1.0e0) end function current_pdf end module mod_rej_f From 98d87400deef6d18054763a41a7d028fc850e8d1 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Fri, 21 Aug 2026 12:06:38 +0200 Subject: [PATCH 09/16] Fixed a few typos/bugs --- particles/initialisers/initialisers_base.f90 | 17 +++++++++-------- .../initialisers/mod_initialise_particles.f90 | 2 +- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/particles/initialisers/initialisers_base.f90 b/particles/initialisers/initialisers_base.f90 index ad83e5919c..91aab1895c 100644 --- a/particles/initialisers/initialisers_base.f90 +++ b/particles/initialisers/initialisers_base.f90 @@ -89,11 +89,11 @@ subroutine part_inout_s(p_inout,n_x,x,time,fields,n_real_param,& !> come first and the n_mhd JOREK variables last. !> !> gradP is (d/dR, d/dZ, d/dphi) and is set to zero when needs_grad is .false. -pure function eval_rej_vars(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf) result(f) +pure function eval_rej_f(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf) result(f) !> i/o vars type(type_node_list), intent(in) :: node_list - type(tyep_element_list), intent(in) :: element_list + type(type_element_list), intent(in) :: element_list integer, intent(in) :: i_elm real*8, intent(in) :: s, t !> local element coordinates real*8, intent(in) :: R, Z, phi !> global cylindrical coordinates @@ -116,9 +116,9 @@ pure function eval_rej_vars(node_list, element_list, i_elm, s, t, phi, R, Z, spa do k = 1, n_geom select case (space_pdf%vars(k)) case(0); P(k) = 1.d0 - case(-1); P(k) = R ; if (needs_grad) gradP(1,k) = 1.d0 - case(-2); P(k) = Z ; if (needs_grad) gradP(2,k) = 1.d0 - case(-1); P(k) = phi ; if (needs_grad) gradP(3,k) = 1.d0 + case(-1); P(k) = R ; if (space_pdf%needs_grad) gradP(1,k) = 1.d0 + case(-2); P(k) = Z ; if (space_pdf%needs_grad) gradP(2,k) = 1.d0 + case(-3); P(k) = phi ; if (space_pdf%needs_grad) gradP(3,k) = 1.d0 endselect enddo @@ -142,12 +142,12 @@ pure function eval_rej_vars(node_list, element_list, i_elm, s, t, phi, R, Z, spa enddo else call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & - s, t, phi, P(n_geom+1), R_i, Z_i) + s, t, phi, P(n_geom+1:), R_i, Z_i) endif endif f = space_pdf%f(size(space_pdf%vars), P, gradP) -end function eval_rej_vars +end function eval_rej_f !> Set positions for particles by rejection sampling from geometric and mhd @@ -651,8 +651,9 @@ subroutine initialise_particles_H_mu_psi(particles, fields, rng_base, mass, T_ma real*8, intent(in), optional :: alpha !< Make more fast (>0) or slow (<0) particles and weigh them appropriately real*8, intent(in), optional :: E_max !< If alpha=1 we select particles from a block-distribution, up to E_max logical, intent(in), optional :: include_vpar !< Initialize particles with local parallel velocity + logical, intent(in), optional :: uniform_space !< Do not + !< use {psi,theta}_transform if present bute rejection sampling in RZ type(spatial_pdf), intent(in), optional :: space_pdf !< spatial pdf to reject against when uniform_space is set. If absent the RZ sampling is uniform - integer, dimension(:), intent(in), optional :: uniform_space_rej_vars !< Variables to use for uniform_space_rej_f type(coronal), intent(in), optional :: cor !< Coronal equilibrium datatype for this particle. If unset, do not alter q integer, intent(in), optional :: charge !< Use this if cor is not present real*8, intent(in), optional :: T_Maxwell !< constant Maxwellian temperature [eV] diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 6e5e9d1128..31efcbf771 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -98,7 +98,7 @@ subroutine initialise_group(sim, group_num) end select !> Set real-space pdf - space_pdf = spatial_pdf_from_name(trim(config%init_pdf)) + space_pdf = spatial_pdf_from_name(trim(config%init_pdf), sim%fields%node_list, sim%fields%element_list) !> Select phase-space initialiser and sample in both real- and phase-space select case(trim(config%init_function)) From ede436e0c1fb04484cf9a0820b155133482a6403 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Fri, 21 Aug 2026 12:44:33 +0200 Subject: [PATCH 10/16] Moved eval_rej_f into the rejection function module --- particles/initialisers/initialisers_base.f90 | 73 +------------------ particles/initialisers/mod_rej_f.f90 | 76 +++++++++++++++++++- 2 files changed, 76 insertions(+), 73 deletions(-) diff --git a/particles/initialisers/initialisers_base.f90 b/particles/initialisers/initialisers_base.f90 index 91aab1895c..b958762fc6 100644 --- a/particles/initialisers/initialisers_base.f90 +++ b/particles/initialisers/initialisers_base.f90 @@ -9,7 +9,7 @@ module initialisers_base use mod_rej_f implicit none private - public initialise_particles, adjust_particle_weights, eval_rej_f + public initialise_particles, adjust_particle_weights public set_velocity_from_T, domain_bounding_box, initialise_particles_H_mu_psi public initialise_particles_H_mu_psi_phiplanes public set_particle_weights_canonical_maxwellian, normalize_with_projection @@ -79,77 +79,6 @@ subroutine part_inout_s(p_inout,n_x,x,time,fields,n_real_param,& end interface contains -!> Evaluate the variables a spatial rejection function needs at one sample point. -!> -!> Fills P(:), and gradP(:,:) when needs_grad, in the order given by variables(:), -!> following the index convenction documented in mod_rej_f: -!> >0: JOREK variable number -!> 0: constant 1, -1: R, -2: Z, -3: Phi -!> variables(:) must be sorted ascending, so the n_geom non-positive entries -!> come first and the n_mhd JOREK variables last. -!> -!> gradP is (d/dR, d/dZ, d/dphi) and is set to zero when needs_grad is .false. -pure function eval_rej_f(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf) result(f) - - !> i/o vars - type(type_node_list), intent(in) :: node_list - type(type_element_list), intent(in) :: element_list - integer, intent(in) :: i_elm - real*8, intent(in) :: s, t !> local element coordinates - real*8, intent(in) :: R, Z, phi !> global cylindrical coordinates - type(spatial_pdf), intent(in) :: space_pdf - real*4 :: f !> acceptance probability - - !> internal vars - real*8, parameter :: EPS_JAC = 1.d-12 !> below this element is degenerate - integer :: n_geom, n_mhd, k - real*8 :: P(size(space_pdf%vars)), gradP(3,size(space_pdf%vars)) - real*8 :: P_s(count(space_pdf%vars > 0)), P_t(count(space_pdf%vars > 0)) - real*8 :: P_phi(count(space_pdf%vars > 0)) - real*8 :: R_i, R_S, R_t, Z_i, Z_s, Z_t, xjac, inv_xjac - - n_mhd = count(space_pdf%vars > 0) - n_geom = size(space_pdf%vars) - n_mhd - gradP = 0.d0 - - !> Geomteric vars - do k = 1, n_geom - select case (space_pdf%vars(k)) - case(0); P(k) = 1.d0 - case(-1); P(k) = R ; if (space_pdf%needs_grad) gradP(1,k) = 1.d0 - case(-2); P(k) = Z ; if (space_pdf%needs_grad) gradP(2,k) = 1.d0 - case(-3); P(k) = phi ; if (space_pdf%needs_grad) gradP(3,k) = 1.d0 - endselect - enddo - - !> MHD Variables : values only unless the rej_f asked for gradients - !> interpolated within the element - if (n_mhd .ge. 1) then - if (space_pdf%needs_grad) then - call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & - s, t, phi, P(n_geom+1:), P_s, P_t, P_phi, R_i, R_s, R_t, Z_i, Z_s, Z_t) - - !> interp_PRZ gives s,t derivatives - transform to R,Z - !> degenerate element (e.g axis point) -> leave gradient zero - xjac = R_s*Z_t - R_t*Z_s - inv_xjac = 0.d0 - if (abs(xjac) > EPS_JAC) inv_xjac = 1.d0/xjac - - do k = 1, n_mhd - gradP(1,n_geom+k) = ( Z_t*P_s(k) - Z_s*P_t(k)) * inv_xjac - gradP(2,n_geom+k) = (-R_t*P_s(k) + R_s*P_t(k)) * inv_xjac - gradP(3,n_geom+k) = P_phi(k) - enddo - else - call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & - s, t, phi, P(n_geom+1:), R_i, Z_i) - endif - endif - - f = space_pdf%f(size(space_pdf%vars), P, gradP) -end function eval_rej_f - - !> Set positions for particles by rejection sampling from geometric and mhd !> variables after collecting with transform, within Rbound, Zbound and Phibound !> if present. See [[test_rejection_sampling]] for examples. diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index b0195f51d7..d9be42c021 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -21,6 +21,8 @@ module mod_rej_f use equil_info ! (R_axis,Z_axis), etc.. use mod_model_settings use mpi + use mod_interp + use data_structure implicit none @@ -33,8 +35,10 @@ module mod_rej_f public :: spatial_pdf public :: rej_f + public :: spatial_pdf_from_name + public :: check_spatial_pdf + public :: eval_rej_f public :: itpa_tae_pdf, RZ_pdf, analytical_pdf, current_pdf - public :: spatial_pdf_from_name, check_spatial_pdf ! ============================================================================= ! Interface @@ -222,6 +226,76 @@ subroutine check_spatial_pdf(pdf, caller) endif end subroutine check_spatial_pdf + !> Evaluate the variables a spatial rejection function needs at one sample point. + !> + !> Fills P(:), and gradP(:,:) when needs_grad, in the order given by variables(:), + !> following the index convenction documented in mod_rej_f: + !> >0: JOREK variable number + !> 0: constant 1, -1: R, -2: Z, -3: Phi + !> variables(:) must be sorted ascending, so the n_geom non-positive entries + !> come first and the n_mhd JOREK variables last. + !> + !> gradP is (d/dR, d/dZ, d/dphi) and is set to zero when needs_grad is .false. + pure function eval_rej_f(node_list, element_list, i_elm, s, t, phi, R, Z, space_pdf) result(f) + + !> i/o vars + type(type_node_list), intent(in) :: node_list + type(type_element_list), intent(in) :: element_list + integer, intent(in) :: i_elm + real*8, intent(in) :: s, t !> local element coordinates + real*8, intent(in) :: R, Z, phi !> global cylindrical coordinates + type(spatial_pdf), intent(in) :: space_pdf + real*4 :: f !> acceptance probability + + !> internal vars + real*8, parameter :: EPS_JAC = 1.d-12 !> below this element is degenerate + integer :: n_geom, n_mhd, k + real*8 :: P(size(space_pdf%vars)), gradP(3,size(space_pdf%vars)) + real*8 :: P_s(count(space_pdf%vars > 0)), P_t(count(space_pdf%vars > 0)) + real*8 :: P_phi(count(space_pdf%vars > 0)) + real*8 :: R_i, R_S, R_t, Z_i, Z_s, Z_t, xjac, inv_xjac + + n_mhd = count(space_pdf%vars > 0) + n_geom = size(space_pdf%vars) - n_mhd + gradP = 0.d0 + + !> Geomteric vars + do k = 1, n_geom + select case (space_pdf%vars(k)) + case(0); P(k) = 1.d0 + case(-1); P(k) = R ; if (space_pdf%needs_grad) gradP(1,k) = 1.d0 + case(-2); P(k) = Z ; if (space_pdf%needs_grad) gradP(2,k) = 1.d0 + case(-3); P(k) = phi ; if (space_pdf%needs_grad) gradP(3,k) = 1.d0 + endselect + enddo + + !> MHD Variables : values only unless the rej_f asked for gradients + !> interpolated within the element + if (n_mhd .ge. 1) then + if (space_pdf%needs_grad) then + call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & + s, t, phi, P(n_geom+1:), P_s, P_t, P_phi, R_i, R_s, R_t, Z_i, Z_s, Z_t) + + !> interp_PRZ gives s,t derivatives - transform to R,Z + !> degenerate element (e.g axis point) -> leave gradient zero + xjac = R_s*Z_t - R_t*Z_s + inv_xjac = 0.d0 + if (abs(xjac) > EPS_JAC) inv_xjac = 1.d0/xjac + + do k = 1, n_mhd + gradP(1,n_geom+k) = ( Z_t*P_s(k) - Z_s*P_t(k)) * inv_xjac + gradP(2,n_geom+k) = (-R_t*P_s(k) + R_s*P_t(k)) * inv_xjac + gradP(3,n_geom+k) = P_phi(k) + enddo + else + call interp_PRZ(node_list, element_list, i_elm, space_pdf%vars(n_geom+1:), n_mhd, & + s, t, phi, P(n_geom+1:), R_i, Z_i) + endif + endif + + f = space_pdf%f(size(space_pdf%vars), P, gradP) + end function eval_rej_f + ! ============================================================================= ! Rejection functions From 2e311bb27bd49fae5769747718b6fc71e279e87d Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Fri, 21 Aug 2026 13:48:23 +0200 Subject: [PATCH 11/16] Removed accidentally commited files --- algexpr2fort | Bin 21240 -> 0 bytes generate_code | 0 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100755 algexpr2fort delete mode 100644 generate_code diff --git a/algexpr2fort b/algexpr2fort deleted file mode 100755 index ba29d1d0ac04a1f5246a219c9d67dfcbdc3382a6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 21240 zcmeHPdyE^$d7mZkbSIy6x`%9^Y&W7DKjhlv9Z%Hxbhc!ko;*j^TdIoIUN6ZdxfZ#k zc9*Bq6-^ZQQK>J;eZeWxMnLm-(_A9+8AlZ!(B4N4(-YzL27;*^-l3g>U6{;X- zz;MVHlU*;NiWaJV3^%HBGGx6*elapeilrrBSXKpN7__6?kfjY~ss&SmoF0IviB>W~iwQ#>-ZADm7lt$IF&o zUx+VEPsOKF3AdVlu!x>q0LUs*dLv#L-VXmS^hus%-}NBy z$bzEb6rrV~UF1|LDu!i?O0HVBJps1s<#YGlM@?jk$U?z8y;-l+EaGf#a2iz2a4M!} zxnPhmh`TwXY}Ab0NhPb+OgmdPNHsA*%B6f^i(V@0DhASLGAIdi+A9XM770(vg48U% z?wTlgVPcxNDd`xrV(V7Xt~#dfopf5oWKwEjH3`04#MvUVAxs9USgq8EkhvFgh;D;A zRX6e#!;?;BDRREm-hFdBcIuOf%_+g?B;4-ZbKvmN!?XMKX|y8hF%kKT;>_cSMmR*& z{w8qz;d8%%uGb?56|a7)g>D%Y>&YBR!|(IxB&1KN^dG;1E*lYdDf#t>(~?j9SFt4h zUVzVIEeXcbX(^|EPkoP5EWqb}PqIw`KG%Wx69K;bUn3^cbbucmQ_=xG<=su)2y`RR zjX*a7-3W9e@c%ag-&^~U_cKp@yFc>_z26uYBJ<3W7g>HS^VFa9Uy(ss-ttYbm+yKD z=d~Nsz=$r9?An{lNO%1$(KMC0_L`*sn&=48OOpOF(KMyGc1hBIN;FMnu3eP$OGMKY z=Gq0&$1>;NzWO&nu6_=um!F8G#TRJ4a`jh<_W`iv(qmT&iM1P^1}+4ZKlv|n=fC!N z=KMD^PrY^h;L+X5rR0|~ufBR=2*tf=6z>16RX<>U;qnPgj56Ek?2TR~j}(H{Um(Wv zn}rKk3Qb$8|I1SUli!ix{PoN$-+4Im%6rj__QlNCKJbPc9{QArusxOc?6$v4Afx)I ztn4Bx_Ql@61V?-9%8L9*|DS)pSq3u6eDd2v5rpK@2g}R2ZZ8g^(7*Kg>;DDn)o0!V z=!vDptIr6C^Hibp&yeiuEUEi?FVO_~-8S7c0On`pEGXm?XL0^JC7 zBhZaNHv-)VbR*D>KsN&22y`RRjllou2+;G8N;Pkm-EpH_G#6^lBt0ui6t>dS5h1p3 zkL`pwRy1wXG3ZfQ%yMH@yL=`Fy%?SoS%z(t&v+J|S;c%$iP%FAaZq;Q;m7h?pQ4)t z{LEv6QaB4p&ndLqq7TwDU%Cm0j=%#K2 zx)JC`pc{d11iBIEMxYykZUnj!_#Z~#;&578)D-0P#6{(Y>E930uY~BYhUmW!(f<^p z-w)9PDkIMCt`Pk|h~6Kf9}CgN5dHBG{aYdWQi!H^)F>$+{whSj9-`k4(ccTv^o|^* z8-0KURa#MuP6B(xmy&!*5~UpKdor4 z2dyh9(L24AIDb4*mG!y4DlPeV*ehwy7t2XWzpms`>3{LW@0EwMoKn|&kAlY)EGRgy z;Ij&b?f357ITPE2H=pl~C6kFQiHX?c#AIsXf#lZMro(1FmN7gfh$oZx3Sa#Qz4ICa z4kNtWiFaLRkI$rP$;qu#;_&eqr*7MpU5w!!(Xdmo2x7{-n!a9OeOBpQN@{pm-M&Z6sMYJg6rDxnS+aZ#v!4$;Ly z|K2-$r?%NhHuiRWOI^EmI`1kO!*)76CfxiXet4X?U5yv*{ER9Y!*)8H$88^y7w552 z`Jfm;d)JB8jjOT#UwijP-vDzga%f#70?Xbp zh`DJ+-xk^`ZPkTUCn95zH8(?a-1N}c1hMI)W#_likh}j_Y>&T(Vsc_)T8+4|c*SsX zrEQD$yy>JsRVoWO*Q&Lg?E!IJ(?4&Pz%_knrV>e?pGr&uG3Ln%XCh#eG63yryihU7 zDX@iV4G3tQp|>SDSpb%;SpHc$F5dqYRt#c2!n3FTklx|Qv`^(K6d>f6Bj6U@tRb&yVZl|gDwX9mHcAvN7?n@cr7AD_Eud1J&ypN4 zZ#h#b%E`y6c9qkH$w_S6 zquqK=$t-V;FPe609Fs-NFcTP`d{oXWpc75ERn~tNPxpE!P#^QyVW5yM=ZM^O}D3cU5bjE2HQ+ysB zy6VJvGgc|M!>rAwvzCo*2zxH5-+t>F?^C1kqi8z*V4u2I$gOI@2Ydql=0!86rE0}Q z$xX*xz{rCqdgD*#3-TE(z)g-SQTlH@c>ArN3CGn^j z5WG$e3ZpdLRI|cX{u&W$%O7k^Vq5js3ZBom@-e05`DZ{76~lt}^#<73T{j|@+Uk!A zzK_tVzfSP|gI4}cqTTv6DsC3MFRxV}Grx2?$dY&ztMyD<{(3?0d52R}@Li2oNmTIN zk1+ps@x&A9Rwk&ZpjN&zMa4%}GQo)oernz-j*3_-Ls3z&u>=2(4t!b}tV~gH=L)6< zPLKHRa#(OP7#d2-{V1js2)5Htf*%osE7oSy;A>ba@w$fZ8|~`g=YaAfhSR{BNzSz4 zehsfm^Sh}L|fP5-?I{1MbQJP_Xx{xB9J;r^0U z`b%n@;emTv@+ptjM2pXn$jWi%r@p!O!ImeuSKH;m#Z;Uij`p z&he6oDz@9>=A4dcp9s!f`~{k&{?W*#2eYUr5gS=j#X7-w|MhnzE_fay5C zH&Cxs&cGx<>2&w5MYUQel&eOb%cbjk4$nTiTi?BZ7k$KF@4f>&X7}j__Ut*b`>1|& zcE`Tm;O_d#{j-nG?SyNr(e>RKWs});*q_Jg`q4*svJ2i=-}nwfI5{dj78~%LgIlcsUX%D3>Fvv0vT92*rYhg_@hg;PRi;5$S65()GwS|o~ux+^d12|m>(Ter&j zgyCVCm#ur~n2EfZtry|ku!~L7vJ2IQD4TW6d8TOb$T;H3D6w(9nshZ+#u(7>@7P73 zZW-}@G=gS&7+Ue)@29b7S|a>CBc87)vD9UpAkDPsSrXe%D|`NaPYkkP@>fDl&J@s~ zJFC)W_AV${iksL3o8GGsINiuFG{#eEEDI?erUe(?|62__Ph?dr0lu= zVSD9H><2>6>e!ywORp(=n%_{OI#QfD|1n}B9mhecY|rbe7gPcHc_8UhA|Ju@I0$+c z$M(GbdR-Nm|LeH?tjF*q?CIGa)4Wc*t_%*U%=j8{`MC@<&!nQVJ+J4|WPmi`gB9`T z`acah+3!$#yzZ+h`wIbk%=%Wu@0Y=&m~j1snrF~oIi>JQMupWE;(Q)Q&_1T@W6HkP z(DV_-X7}RPLiRkLyr2}H3SBU^XScr@vgdW_I(6arIW_Bt%YRYXbNqS!B4WP9IUh96 z;rw5Kgo?=ayv{ur^Gz`Lx1@0Xe++5Rp4Y=mG2iU>TCAJbkJ(E=263qF(pWe1cdo+i zjcJBof$iFWJ+Esn-w!4d>-kiZW#YCD`iufzGfWfMxEnq$Ks$MGAmrTViyuh*}S zgNej?KD9FcZ$S{WZ>;Mle2WYlD#+z#_)f^4*Z&!1-zk4~OY14}#pUO9zzJo~`LI5f zHJIK9K|H125V8H<2FFL(o?#S7JNwTlhuc|E`C*db2oSzDq{8;;Qqk>U-X^~_nn~qxm`Z>^^o&${=*&Y OA5$IY`H(?~6#otCeaAWg diff --git a/generate_code b/generate_code deleted file mode 100644 index e69de29bb2..0000000000 From e654e6e061865d29623d79ca71f4e1549fe1d104 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Fri, 21 Aug 2026 13:59:10 +0200 Subject: [PATCH 12/16] Removed redundant comment --- particles/initialisers/mod_rej_f.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/particles/initialisers/mod_rej_f.f90 b/particles/initialisers/mod_rej_f.f90 index d9be42c021..3e005c3f9d 100644 --- a/particles/initialisers/mod_rej_f.f90 +++ b/particles/initialisers/mod_rej_f.f90 @@ -369,7 +369,6 @@ end function analytical_pdf !> P(2) = j_tor !> !> NOTE : Won't work for fMHD models, see note in constructor above - !> NOTE : jzmin and jzmax are currently hardcoded, needs better handling pure function current_pdf(n, P, gradP) result(f) integer, intent(in) :: n real*8, dimension(n), intent(in) :: P From 2bc631a56257076cc3f38378279aa293baa0cca3 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Fri, 21 Aug 2026 13:59:42 +0200 Subject: [PATCH 13/16] Fixed seed_particles example --- particles/examples/seed_particles.f90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/particles/examples/seed_particles.f90 b/particles/examples/seed_particles.f90 index cad88c5828..4078c53828 100644 --- a/particles/examples/seed_particles.f90 +++ b/particles/examples/seed_particles.f90 @@ -4,9 +4,12 @@ program seed_particles use particle_tracer use mod_particle_io +use mod_rej_f +use initialisers_base implicit none type(event) :: fieldreader +type(spatial_pdf) :: space_pdf integer :: i, j ! Start up MPI, jorek @@ -16,6 +19,10 @@ program seed_particles fieldreader = event(read_jorek_fields_interp_linear(basename='jorek', i=-1)) call with(sim, fieldreader) +! Set up spatial pdf type +space_pdf%f => f_psi_inside +space_pdf%vars = [1] + ! Set up particles sim%groups(:)%Z = -2 sim%groups(:)%mass = 2.d0 !< atomic mass units @@ -26,16 +33,16 @@ program seed_particles call initialise_particles_H_mu_psi(sim%groups(i)%particles, & sim%fields, pcg32_rng(), & sim%groups(i)%mass, & - uniform_space=.true., uniform_space_rej_f=f_psi_inside, & - uniform_space_rej_vars=[1], charge=1) + uniform_space=.true., space_pdf=space_pdf, charge=1) sim%groups(i)%particles(:)%weight = 1.0 end do call write_simulation_hdf5(sim, 'part_restart.h5') contains pure function f_psi_inside(n, P, grad_P) result(f) - integer, intent(in) :: n - real*8, intent(in) :: P(n), grad_P(3,n) + integer, intent(in) :: n + real*8, dimension(n), intent(in) :: P + real*8, dimension(3,n), intent(in) :: grad_P real*4 :: f f = 5e-3 if (P(1) .lt. -0.26 .and. P(1) .gt. -0.35) f = 1e0 From 1d87554613070dd18f70d38bc738feeaeef509c5 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Tue, 25 Aug 2026 12:30:43 +0100 Subject: [PATCH 14/16] Fixed broken use statements, mostly mod_initialise_particles -> initialisers_base --- particles/examples/plasma_volume.f90 | 2 +- particles/examples/v_ExB_hist.f90 | 2 +- particles/particle_tracer.f90 | 1 + particles/tests/mod_fieldline_spec_mpi_test.f90 | 1 + particles/tests/mod_particle_projection_spec_mpi_test.f90 | 4 ++-- 5 files changed, 6 insertions(+), 4 deletions(-) diff --git a/particles/examples/plasma_volume.f90 b/particles/examples/plasma_volume.f90 index ebe55ce931..0e6f01cbd7 100644 --- a/particles/examples/plasma_volume.f90 +++ b/particles/examples/plasma_volume.f90 @@ -12,7 +12,7 @@ program plasma_volume use mod_random_seed, only: random_seed use mod_sobseq_rng, only: sobseq_rng use constants, only: TWOPI, PI -use mod_initialise_particles, only: domain_bounding_box +use initialisers_base, only: domain_bounding_box use domains use hdf5_io_module use mod_sampling, only: transform_uniform_cylindrical diff --git a/particles/examples/v_ExB_hist.f90 b/particles/examples/v_ExB_hist.f90 index 6041469db8..b7ef790360 100644 --- a/particles/examples/v_ExB_hist.f90 +++ b/particles/examples/v_ExB_hist.f90 @@ -15,7 +15,7 @@ program v_ExB_hist use mod_sobseq_rng, only: sobseq_rng use mod_random_seed, only: random_seed use constants, only: TWOPI -use mod_initialise_particles, only: domain_bounding_box +use initialisers_base, only: domain_bounding_box use mod_math_operators, only: cross_product use domains use equil_info, only:find_xpoint diff --git a/particles/particle_tracer.f90 b/particles/particle_tracer.f90 index 1f7b668b9a..e385c84dbb 100644 --- a/particles/particle_tracer.f90 +++ b/particles/particle_tracer.f90 @@ -5,6 +5,7 @@ module particle_tracer use mod_particle_types use mod_event use mod_initialise_particles +use initialisers_base ! IO use mod_io_actions diff --git a/particles/tests/mod_fieldline_spec_mpi_test.f90 b/particles/tests/mod_fieldline_spec_mpi_test.f90 index 2e43677eb1..d75db44d7f 100644 --- a/particles/tests/mod_fieldline_spec_mpi_test.f90 +++ b/particles/tests/mod_fieldline_spec_mpi_test.f90 @@ -6,6 +6,7 @@ module mod_fieldline_spec_mpi_test use mod_particle_types use mod_sobseq_rng use mod_initialise_particles +use initialisers_base use mod_fields_linear use mod_fieldline_euler use mod_neighbours diff --git a/particles/tests/mod_particle_projection_spec_mpi_test.f90 b/particles/tests/mod_particle_projection_spec_mpi_test.f90 index f69fb468f0..a62eacfb5a 100644 --- a/particles/tests/mod_particle_projection_spec_mpi_test.f90 +++ b/particles/tests/mod_particle_projection_spec_mpi_test.f90 @@ -345,7 +345,7 @@ subroutine project_n(rank,master,n_tasks,node_list,element_list,proj_f_proj,& use mpi_mod use data_structure use mod_rng, only: type_rng - use mod_initialise_particles, only: initialise_particles + use initialisers_base, only: initialise_particles use mod_particle_sim, only: particle_sim use mod_project_particles, only: projection,new_projection, write_particle_distribution_to_h5 use mod_rhs_projections, only: proj_f @@ -438,7 +438,7 @@ subroutine rhs_convergence(rank,n_tasks,node_list,element_list,& use constants, only: TWOPI use data_structure, only: type_node_list,type_element_list use mod_rng, only: type_rng - use mod_initialise_particles, only: initialise_particles + use initialisers_base, only: initialise_particles use mod_particle_sim, only: particle_sim use mod_project_particles, only: projection, new_projection, sample_rhs use mod_rhs_projections, only: proj_f From 65df06ee8765ea8aeeb8e02fa18faefd57858978 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Tue, 22 Sep 2026 11:40:23 +0100 Subject: [PATCH 15/16] Added use mod_interp to particle example files - should fix compile error related to calling interp_PRZ --- particles/examples/E_diagnostic.f90 | 1 + .../examples/test_generalised_initialisation_gc_H_mu_psi.f90 | 1 + 2 files changed, 2 insertions(+) diff --git a/particles/examples/E_diagnostic.f90 b/particles/examples/E_diagnostic.f90 index 9de5bacd28..84786df5b8 100644 --- a/particles/examples/E_diagnostic.f90 +++ b/particles/examples/E_diagnostic.f90 @@ -12,6 +12,7 @@ module write_E_time use mod_event, only: event, action use equil_info, only:find_xpoint + use mod_interp implicit none type, extends(action) :: save_E diff --git a/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90 b/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90 index d06c0b3aec..6059c055ef 100644 --- a/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90 +++ b/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90 @@ -4,6 +4,7 @@ program test_generalised_initialisation_gc_H_mu_psi !> initialize_particles_H_mu_psi implemented in mod_initialise_particles use constants, only: TWOPI,PI,ATOMIC_MASS_UNIT,EL_CHG use particle_tracer +use mod_interp implicit none type(pcg32_rng) :: rng_pcg32 type(event) :: field_reader,particle_writer From 1aae8905ca8c0614164d1cb0d240ff497002ead1 Mon Sep 17 00:00:00 2001 From: James Carpenter Date: Tue, 6 Oct 2026 15:30:37 +0200 Subject: [PATCH 16/16] Addressing Edo's comments. Namely adding default branch to init_function x particle_type check --- .../runaway_electrons/rep_tutorial.md | 2 +- particles/diagnostics/mod_neutral_density.f90 | 2 +- .../initialisers/mod_initialise_particles.f90 | 16 +++++++++++++--- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md b/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md index 858dd80a72..e46f67033a 100644 --- a/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md +++ b/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md @@ -29,7 +29,7 @@ Apart from the jorek_restart.h5, an part_restart.h5 file is required by any simu part_group_configs(1)%type = 'particle_kinetic_relativistic' part_group_configs(1)%n_particles = 1e6 ! number of super particles part_group_configs(1)%mass = 0.000548579870184805 - part_group_configs(1)%num_re = 1.175403e16 ! number of physical REs, doesn't work if initialized according to current + part_group_configs(1)%n_particles_total = 1.175403e16 ! number of physical REs, doesn't work if initialized according to current part_group_configs(1)%use_kin_recombination = .false. ! Projecting (particle field -> MHD) smoothing diff --git a/particles/diagnostics/mod_neutral_density.f90 b/particles/diagnostics/mod_neutral_density.f90 index 7d84f80685..b1c23cca2b 100644 --- a/particles/diagnostics/mod_neutral_density.f90 +++ b/particles/diagnostics/mod_neutral_density.f90 @@ -4,7 +4,6 @@ module mod_neutral_density use mod_interp, only: mode_moivre use mod_basisfunctions use particle_tracer - use phys_module, only: n_part_groups !$ use omp_lib implicit none @@ -17,6 +16,7 @@ module mod_neutral_density !> first calculates the neutral density from all particles with ncs, then projects it subroutine get_neutral_density(sim, neutral_density_proj) + use phys_module, only: n_part_groups implicit none class(particle_sim), target, intent(inout) :: sim diff --git a/particles/initialisers/mod_initialise_particles.f90 b/particles/initialisers/mod_initialise_particles.f90 index 31efcbf771..e29eff5e33 100644 --- a/particles/initialisers/mod_initialise_particles.f90 +++ b/particles/initialisers/mod_initialise_particles.f90 @@ -86,15 +86,25 @@ subroutine initialise_group(sim, group_num) select type(particles => sim%groups(group_num)%particles) type is (particle_kinetic_relativistic) if (trim(config%init_function) /= 're_gaussian') then - write(*,*) "ERROR : particle_kinetic_relativistic requred init_function='re_gaussian'" - write(*,*) " got init_function='",trim(config%init_function),"' for group '", trim(config%id), "'" + if (sim%my_id .eq. 0) then + write(*,*) "ERROR : particle_kinetic_relativistic requred init_function='re_gaussian'" + write(*,*) " got init_function='",trim(config%init_function),"' for group '", trim(config%id), "'" + endif call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif type is (particle_kinetic_leapfrog) if (trim(config%init_function) == 're_gaussian') then - write(*,*) "ERROR (initialise_group): particle_kinetic_leapfrog incompatible with init_function='re_gaussian'" + if (sim%my_id .eq. 0) then + write(*,*) "ERROR : particle_kinetic_leapfrog incompatible with init_function='re_gaussian'" + endif call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) endif + class default + if (sim%my_id .eq. 0) then + write(*,*) "ERROR : particle type may not be compatible with chosen init_function" + write(*,*) " see subroutine initialise_group in mod_initialise_particles to change behaviour" + endif + call MPI_ABORT(MPI_COMM_WORLD, 1, ierr) end select !> Set real-space pdf