Merge branch '4.20-devel'

Author: TomGoffrey

epoch1d/Makefile

@@ -234,6 +234,9 @@ DEFINES := $(DEFINE)
 # Don't generate any output at all. Useful for benchmarking.
 #DEFINES += $(D)NO_IO
 
+# Track ionisation/recombination rates on each particle
+#DEFINES += $(D)TRANSITION_RATES
+
 # Specify charge and mass per particle rather than per species.
 # This is a legacy option which will be removed soon.
 #DEFINES += $(D)PER_PARTICLE_CHARGE_MASS
@@ -304,7 +307,7 @@ FC_INFO := $(shell ${FC} --version 2>/dev/null \
 
 
 SRCFILES = background_collisions.F90 balance.F90 bethe_heitler.F90 \
-  boundary.F90 bremsstrahlung.F90 calc_df.F90 collisions.F90 \
+	boundary.F90 bremsstrahlung.F90 calc_df.F90 collisions.F90 \
 	collision_ionise.F90 constants.F90 current_smooth.F90 custom_deck.f90 \
 	custom_laser.f90 custom_parser.f90 deck.F90 deck_boundaries_block.f90 \
 	deck_bremsstrahlung_block.F90 deck_collision_block.F90 \
@@ -319,10 +322,11 @@ SRCFILES = background_collisions.F90 balance.F90 bethe_heitler.F90 \
 	mpi_routines.F90 mpi_subtype_control.f90 numerics.f90 particle_id_hash.F90 \
 	particle_migration.F90 particle_pointer_advance.f90 particle_temperature.F90 \
 	particles.F90 partlist.F90 photons.F90 prefetch.F90 probes.F90 \
-	random_generator.f90 redblack_module.f90 setup.F90 shape_functions.F90 \
-	shared_data.F90 shunt.F90 simple_io.F90 split_particle.F90 stack.f90 \
-	strings.f90 strings_advanced.f90 terminal_controls.F90 timer.f90 \
-	tokenizer_blocks.f90 utilities.f90 version_data.F90 welcome.F90 window.F90
+	random_generator.f90 recombination.F90 redblack_module.f90 \
+	secondary_list.F90 setup.F90 shape_functions.F90 shared_data.F90 shunt.F90 \
+	simple_io.F90 stack.f90 strings.f90 strings_advanced.f90 \
+	terminal_controls.F90 timer.f90 tokenizer_blocks.f90 utilities.f90 \
+	version_data.F90 welcome.F90 window.F90
 
 OBJFILES := $(SRCFILES:.f90=.o)
 OBJFILES := $(OBJFILES:.F90=.o)
@@ -489,7 +493,7 @@ deck_particle_probe_block.o: deck_particle_probe_block.F90 probes.o \
   strings_advanced.o utilities.o
 deck_qed_block.o: deck_qed_block.F90 strings_advanced.o utilities.o
 deck_species_block.o: deck_species_block.F90 partlist.o setup.o simple_io.o \
-  strings_advanced.o utilities.o
+  shared_data.o strings_advanced.o utilities.o
 deck_stencil_block.o: deck_stencil_block.f90 fields.o strings_advanced.o
 deck_subset_block.o: deck_subset_block.F90 particle_id_hash.o \
   strings_advanced.o utilities.o
@@ -502,8 +506,8 @@ dist_fn.o: dist_fn.F90 mpi_subtype_control.o particles.o $(SDFMOD)
 epoch1d.o: epoch1d.F90 background_collisions.o balance.o bremsstrahlung.o \
 	calc_df.o collisions.o collision_ionise.o current_smooth.o deck.o \
 	diagnostics.o fields.o finish.o helper.o ic_module.o injectors.o ionise.o \
-	mpi_routines.o particle_migration.o particles.o photons.o setup.o \
-	split_particle.o welcome.o window.o
+	mpi_routines.o particle_migration.o particles.o photons.o recombination.o \
+	secondary_list.o setup.o welcome.o window.o
 	$(FC) -c $(FFLAGS) $(FLTCONS) -o $(OBJDIR)/$@ $(PREPROFLAGS) $<
 evaluate.o: evaluate.F90 shunt.o stack.o
 evaluator_blocks.o: evaluator_blocks.F90 custom_parser.o stack.o strings.o
@@ -516,8 +520,8 @@ helper.o: helper.F90 balance.o boundary.o deltaf_loader.o partlist.o simple_io.o
 ic_module.o: ic_module.f90 helper.o shared_data.o
 injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.o \
 	partlist.o random_generator.o utilities.o
-ionise.o: ionise.F90 boundary.o numerics.o partlist.o random_generator.o \
-  utilities.o
+ionise.o: ionise.F90 boundary.o calc_df.o numerics.o partlist.o \
+  random_generator.o utilities.o
 iterators.o: iterators.F90 particle_id_hash.o particle_pointer_advance.o \
   partlist.o
 laser.o: laser.f90 custom_laser.o evaluate.o
@@ -532,19 +536,20 @@ particle_temperature.o: particle_temperature.F90 constants.o evaluate.o \
   random_generator.o
 particles.o: particles.F90 boundary.o partlist.o prefetch.o
 partlist.o: partlist.F90 particle_id_hash.o random_generator.o shared_data.o
-photons.o: photons.F90 partlist.o
+photons.o: photons.F90 collisions.o partlist.o utilities.o
 prefetch.o: prefetch.F90 shared_data.o
 probes.o: probes.F90 partlist.o $(SDFMOD)
 random_generator.o: random_generator.f90
+recombination.o: recombination.F90 calc_df.o collisions.o collision_ionise.o
 redblack_module.o: redblack_module.f90 partlist.o
+secondary_list.o: secondary_list.F90 boundary.o
 setup.o: setup.F90 balance.o boundary.o constants.o fields.o helper.o \
   injectors.o laser.o mpi_routines.o mpi_subtype_control.o particle_id_hash.o \
-  shunt.o split_particle.o timer.o version_data.o welcome.o window.o $(SDFMOD)
+  shunt.o timer.o version_data.o welcome.o window.o $(SDFMOD)
 shape_functions.o: shape_functions.F90 constants.o
 shared_data.o: shared_data.F90 constants.o $(SDFMOD)
 shunt.o: shunt.F90 evaluator_blocks.o tokenizer_blocks.o utilities.o
 simple_io.o: simple_io.F90 boundary.o mpi_subtype_control.o
-split_particle.o: split_particle.F90 boundary.o
 stack.o: stack.f90 constants.o
 strings.o: strings.f90 shared_data.o
 strings_advanced.o: strings_advanced.f90 evaluate.o shunt.o

epoch1d/example_decks/ionisation.deck

@@ -24,6 +24,7 @@ end:control
 #
 #   use_collisions = T       # Collisions must be on for collisional ionisation
 #   collisional_ionisation = T             # Switches on collisional ionisation
+#   use_nanbu = T
 #
 # end:collisions
 
@@ -34,6 +35,7 @@ begin:constant
   den_max = 10.0 * critical(omega)
   den_preplasma = 1.0e18
   len_scale = 5.0 * micron
+  ion_charge_state = 22 # Initial ionisation level of ions
 end:constant
 
 
@@ -47,13 +49,15 @@ begin:species
   name = Electron
   charge = -1.0
   mass = 1.0
+  number_density = den_max * ion_charge_state
   nparticles = nx * 5
+  identify:electron
 end:species
 
 
 begin:species
   name = Palladium
-  charge = 22.0
+  charge = ion_charge_state
   atomic_no = 46
   mass = 1836.2 * 105.0
 
@@ -67,13 +71,6 @@ begin:species
 end:species
 
 
-begin:species
-  name = Electron
-  number_density = number_density(Palladium) * 22.0
-  temperature = 1e6
-end:species
-
-
 begin:output
 
   # Simulated time between output dumps

epoch1d/src/constants.F90

@@ -125,7 +125,8 @@ MODULE constants
   INTEGER, PARAMETER :: c_err_io_error = 2**15
   INTEGER, PARAMETER :: c_err_bad_setup = 2**16
   INTEGER, PARAMETER :: c_err_window = 2**17
-  INTEGER, PARAMETER :: c_err_max = 17
+  INTEGER, PARAMETER :: c_err_qed = 2**18
+  INTEGER, PARAMETER :: c_err_max = 18
 
   CHARACTER(LEN=*), PARAMETER :: c_err_char(0:c_err_max) = (/ &
       'unknown_block           ', &
@@ -145,7 +146,8 @@ MODULE constants
       'pp_options_wrong        ', &
       'io_error                ', &
       'bad_setup               ', &
-      'window                  '/)
+      'window                  ', &
+      'qed                     '/)
 
   INTEGER, PARAMETER :: c_ds_first = 1
   INTEGER, PARAMETER :: c_ds_last = 2
@@ -229,6 +231,14 @@ MODULE constants
   REAL(num), PARAMETER :: alpha_f = 7.297352575523020256850802729527158e-3_num
   ! tau_c = h_bar / (m0 * c**2)
   REAL(num), PARAMETER :: tau_c = 1.288088667367242662108649212042082e-21_num
+  REAL(num), PARAMETER :: classical_re = 0.25_num / pi / epsilon0 / m0 &
+      * (q0 / c)**2
+  REAL(num), PARAMETER :: sigma_lBW_max = pi * classical_re**2 * &
+          0.6817055055017870382984600045421994441648264608312_num
+  REAL(num), PARAMETER :: inv_c = 1.0_num / c
+  REAL(num), PARAMETER :: inv_mc0_sq = 1.0_num / mc0 / mc0
+  REAL(num), PARAMETER :: half_pire2 = 0.5_num * pi * classical_re**2
+  REAL(num), PARAMETER :: quarter_pire2 = 0.25_num * pi * classical_re**2
 #endif
 
   ! Constants used for bremsstrahlung with plasma screening
@@ -294,6 +304,7 @@ MODULE constants
   INTEGER(i8), PARAMETER :: c_def_use_mpi3 = 2**25
   INTEGER(i8), PARAMETER :: c_def_bremsstrahlung = 2**26
   INTEGER(i8), PARAMETER :: c_def_probe_time = 2**27
+  INTEGER(i8), PARAMETER :: c_def_transition_rates = 2**28
 
   ! Stagger types
   INTEGER, PARAMETER :: c_stagger_ex = c_stagger_face_x
@@ -530,6 +541,8 @@ MODULE constants
   INTEGER, PARAMETER :: c_func_drifty = 45
   INTEGER, PARAMETER :: c_func_driftz = 46
   INTEGER, PARAMETER :: c_func_arctan2 = 47
+  INTEGER, PARAMETER :: c_func_min = 48
+  INTEGER, PARAMETER :: c_func_max = 49
 
   INTEGER, PARAMETER :: c_func_custom_lowbound = 4096
 
@@ -633,7 +646,12 @@ MODULE constants
   INTEGER, PARAMETER :: c_dump_part_opdepth_brem = 71
   INTEGER, PARAMETER :: c_dump_probe_time        = 72
   INTEGER, PARAMETER :: c_dump_cou_log           = 73
-  INTEGER, PARAMETER :: num_vars_to_dump         = 73
+  INTEGER, PARAMETER :: c_dump_part_rate_fi      = 74
+  INTEGER, PARAMETER :: c_dump_part_rate_ci      = 75
+  INTEGER, PARAMETER :: c_dump_part_rate_dr      = 76
+  INTEGER, PARAMETER :: c_dump_part_rate_rr      = 77
+  INTEGER, PARAMETER :: c_dump_part_rate_3br     = 78
+  INTEGER, PARAMETER :: num_vars_to_dump         = 78
 
   INTEGER, PARAMETER :: c_subset_random     = 1
   INTEGER, PARAMETER :: c_subset_gamma_min  = 2

epoch1d/src/deck/deck_collision_block.F90

@@ -139,7 +139,7 @@ SUBROUTINE collision_deck_finalise
 
     END IF
 
-    IF (use_collisional_ionisation) use_particle_lists = .TRUE.
+    use_particle_lists = use_particle_lists .OR. use_collisional_ionisation
 
   END SUBROUTINE collision_deck_finalise
 

epoch1d/src/deck/deck_control_block.F90

@@ -331,6 +331,15 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
         .OR. str_cmp(element, 'bsi')) THEN
       use_bsi = as_logical_print(value, element, errcode)
 
+    ELSE IF (str_cmp(element, 'use_radiative_recombination')) THEN
+      use_radiative_recombination = as_logical_print(value, element, errcode)
+
+    ELSE IF (str_cmp(element, 'use_dielectronic_recombination')) THEN
+      use_dielectronic_recombination = as_logical_print(value, element, errcode)
+
+    ELSE IF (str_cmp(element, 'use_three_body_recombination')) THEN
+      use_three_body_recombination = as_logical_print(value, element, errcode)
+
     ELSE IF (str_cmp(element, 'particle_tstart')) THEN
       particle_push_start_time = as_real_print(value, element, errcode)
 
@@ -427,6 +436,9 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
     ELSE IF (str_cmp(element, 'smooth_iterations')) THEN
       smooth_its = as_integer_print(value, element, errcode)
 
+    ELSE IF (str_cmp(element, 'recombine_n_step')) THEN
+      recombine_n_step = as_integer_print(value, element, errcode)
+
     ELSE IF (str_cmp(element, 'smooth_compensation')) THEN
       IF (as_logical_print(value, element, errcode)) smooth_comp_its = 1
 

epoch1d/src/deck/deck_io_block.F90

@@ -611,6 +611,23 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode)
       elementselected = c_dump_part_work_z_total
 #endif
 
+#ifdef TRANSITION_RATES
+    ELSE IF (str_cmp(element, 'rate_fi')) THEN
+      elementselected = c_dump_part_rate_fi
+
+    ELSE IF (str_cmp(element, 'rate_ci')) THEN
+      elementselected = c_dump_part_rate_ci
+
+    ELSE IF (str_cmp(element, 'rate_dr')) THEN
+      elementselected = c_dump_part_rate_dr
+
+    ELSE IF (str_cmp(element, 'rate_rr')) THEN
+      elementselected = c_dump_part_rate_rr
+
+    ELSE IF (str_cmp(element, 'rate_3br')) THEN
+      elementselected = c_dump_part_rate_3br
+#endif
+
     ELSE IF (str_cmp(element, 'ex')) THEN
       elementselected = c_dump_ex
 

epoch1d/src/deck/deck_qed_block.F90

@@ -45,6 +45,10 @@ SUBROUTINE qed_deck_initialise
       use_radiation_reaction = .TRUE.
       produce_photons = .FALSE.
       photon_dynamics = .FALSE.
+      use_binary_collisions = .FALSE.
+      use_LBW = .FALSE.
+      use_LBW_diff = .TRUE.
+      LBW_amp_factor = 1.0_num
     END IF
 #endif
 
@@ -57,6 +61,7 @@ SUBROUTINE qed_deck_finalise
     INTEGER :: io, iu
 #ifdef PHOTONS
     LOGICAL :: exists
+    INTEGER :: j
 
     IF (deck_state == c_ds_first) RETURN
 
@@ -74,6 +79,23 @@ SUBROUTINE qed_deck_finalise
     END IF
 
     IF (use_qed) need_random_state = .TRUE.
+
+    use_binary_collisions = use_LBW
+    IF (use_binary_collisions) THEN
+      DO j = 1, n_species
+        IF (species_list(j)%species_type == c_species_id_photon) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+        IF (species_list(j)%species_type == c_species_id_electron) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+        IF (species_list(j)%species_type == c_species_id_positron) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+    END DO
+  END IF
+
+  lbw_amp_factor = MAX(lbw_amp_factor, 1.0_num)
 #else
     IF (use_qed) THEN
       IF (rank == 0) THEN
@@ -87,7 +109,6 @@ SUBROUTINE qed_deck_finalise
       CALL abort_code(c_err_pp_options_missing)
     END IF
 #endif
-
   END SUBROUTINE qed_deck_finalise
 
 
@@ -172,6 +193,21 @@ FUNCTION qed_block_handle_element(element, value) RESULT(errcode)
       RETURN
     END IF
 
+    IF(str_cmp(element, 'linear_breit_wheeler')) THEN
+      use_LBW = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF(str_cmp(element, 'LBW_differential_cross')) THEN
+      use_LBW_diff = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF (str_cmp(element, 'amplify_LBW_factor')) THEN
+      LBW_amp_factor = as_real_print(value, element, errcode)
+      RETURN
+    END IF
+
     errcode = c_err_unknown_element
 #endif