Add Linear Compton Scattering

Author: TomGoffrey

epoch1d/src/constants.F90

@@ -229,6 +229,13 @@ 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_thomson = 8.0_num * pi / 3.0_num * classical_re**2
+  REAL(num), PARAMETER :: inv_mc0 = 1.0_num / mc0
+  REAL(num), PARAMETER :: inv_m0c2 = 1.0_num / m0c2
+  REAL(num), PARAMETER :: pire2 = pi * classical_re**2
+  REAL(num), PARAMETER :: quarter_pire2 = 0.25_num * pire2
 #endif
 
   ! Constants used for bremsstrahlung with plasma screening

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_LCS = .FALSE.
+      use_LCS_diff = .TRUE.
     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,22 @@ SUBROUTINE qed_deck_finalise
     END IF
 
     IF (use_qed) need_random_state = .TRUE.
+
+    use_binary_collisions = use_LCS
+
+    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
 #else
     IF (use_qed) THEN
       IF (rank == 0) THEN
@@ -87,7 +108,6 @@ SUBROUTINE qed_deck_finalise
       CALL abort_code(c_err_pp_options_missing)
     END IF
 #endif
-
   END SUBROUTINE qed_deck_finalise
 
 
@@ -172,6 +192,16 @@ FUNCTION qed_block_handle_element(element, value) RESULT(errcode)
       RETURN
     END IF
 
+    IF(str_cmp(element, 'linear_compton_scattering')) THEN
+      use_LCS = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF(str_cmp(element, 'LCS_differential_cross')) THEN
+      use_LCS_diff = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
     errcode = c_err_unknown_element
 #endif
 

epoch1d/src/epoch1d.F90

@@ -214,7 +214,8 @@ PROGRAM pic
       ! .FALSE. this time to use load balancing threshold
       IF (use_balance) CALL balance_workload(.FALSE.)
       CALL push_particles
-      IF (use_particle_lists) THEN
+
+      IF (use_particle_lists .OR. use_binary_collisions) THEN
         ! Check whether this is a step with collisions or collisional ionisation
         collision_step = (MODULO(step, coll_n_step) == coll_n_step - 1) &
           .AND. use_collisions
@@ -223,7 +224,8 @@ PROGRAM pic
 
         ! After this line, the particles can be accessed on a cell by cell basis
         ! Using the particle_species%secondary_list property
-        IF (use_split .OR. collision_step .OR. coll_ion_step) THEN
+        IF (use_split .OR. collision_step .OR. coll_ion_step &
+          .OR. use_binary_collisions) THEN
           CALL reorder_particles_to_grid
         END IF
 
@@ -239,7 +241,15 @@ PROGRAM pic
         ! Early beta version of particle splitting operator
         IF (use_split) CALL split_particles
 
-        IF (use_split .OR. collision_step .OR. coll_ion_step) THEN
+#ifdef PHOTONS
+        IF (use_binary_collisions) THEN
+          CALL do_binary_collisions
+        END IF
+#endif
+
+        IF (use_split .OR. collision_step .OR. coll_ion_step &
+          .OR. use_binary_collisions) THEN
+
           CALL reattach_particles_to_mainlist
         END IF
       END IF