diff --git a/epoch1d/Makefile b/epoch1d/Makefile index 413225167..a80b9b8f6 100644 --- a/epoch1d/Makefile +++ b/epoch1d/Makefile @@ -536,7 +536,7 @@ 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 prefetch.o: prefetch.F90 shared_data.o probes.o: probes.F90 partlist.o $(SDFMOD) random_generator.o: random_generator.f90 diff --git a/epoch1d/src/constants.F90 b/epoch1d/src/constants.F90 index 66819ea64..c598ab714 100644 --- a/epoch1d/src/constants.F90 +++ b/epoch1d/src/constants.F90 @@ -229,6 +229,15 @@ 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 diff --git a/epoch1d/src/deck/deck_qed_block.F90 b/epoch1d/src/deck/deck_qed_block.F90 index 0dc13d7d5..71aa989df 100644 --- a/epoch1d/src/deck/deck_qed_block.F90 +++ b/epoch1d/src/deck/deck_qed_block.F90 @@ -45,6 +45,9 @@ 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 +60,7 @@ SUBROUTINE qed_deck_finalise INTEGER :: io, iu #ifdef PHOTONS LOGICAL :: exists + INTEGER :: j IF (deck_state == c_ds_first) RETURN @@ -74,6 +78,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 +107,6 @@ SUBROUTINE qed_deck_finalise CALL abort_code(c_err_pp_options_missing) END IF #endif - END SUBROUTINE qed_deck_finalise @@ -172,6 +191,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 diff --git a/epoch1d/src/epoch1d.F90 b/epoch1d/src/epoch1d.F90 index 4f5ae5bdd..793b4db9a 100644 --- a/epoch1d/src/epoch1d.F90 +++ b/epoch1d/src/epoch1d.F90 @@ -216,7 +216,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 @@ -227,7 +228,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 (collision_step .OR. coll_ion_step .OR. recombine_step) THEN + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reorder_particles_to_grid END IF @@ -242,7 +244,13 @@ PROGRAM pic IF (recombine_step) CALL run_recombination - IF (collision_step .OR. coll_ion_step .OR. recombine_step) THEN +#ifdef PHOTONS + IF (use_binary_collisions) THEN + CALL do_binary_collisions + END IF +#endif + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reattach_particles_to_mainlist END IF END IF diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90 index 1197e903b..d19c011d1 100644 --- a/epoch1d/src/physics_packages/photons.F90 +++ b/epoch1d/src/physics_packages/photons.F90 @@ -16,11 +16,14 @@ MODULE photons #ifdef PHOTONS + USE collisions USE partlist + IMPLICIT NONE - SAVE + REAL(num) :: sig2cdt_dV_lcs + REAL(num) :: cdt_dV REAL(num) :: part_pos_global, gamma_global, eta_global CONTAINS @@ -67,6 +70,10 @@ SUBROUTINE setup_qed_module END DO END IF + cdt_dV = c * dt / dx + sig2cdt_dV_lcs = 2.0_num * sigma_thomson * cdt_dV + + END SUBROUTINE setup_qed_module @@ -212,7 +219,6 @@ FUNCTION check_qed_variables() END FUNCTION check_qed_variables - SUBROUTINE setup_tables_qed ! Reads files epsilon.table, log_chi.table, energy_split.table @@ -557,7 +563,7 @@ SUBROUTINE qed_update_optical_depth current%optical_depth = & current%optical_depth - delta_optical_depth(eta, gamma_rel) ENDIF - + #ifdef TRIDENT_PHOTONS current%optical_depth_tri = current%optical_depth_tri & - delta_optical_depth_tri(eta, gamma_rel) @@ -920,7 +926,7 @@ SUBROUTINE generate_photon(generating_electron, iphoton, eta) samp_photon = .TRUE. IF (photon_sample_fraction < 1.0_num) THEN - samp_photon = random() < photon_sample_fraction + samp_photon = random() < photon_sample_fraction END IF ! This will only create photons that have energies above a user specified @@ -1407,5 +1413,432 @@ FUNCTION find_value_from_table(x_in, p_value, nx, ny, x, y, p_table) END FUNCTION find_value_from_table + + + SUBROUTINE do_binary_collisions + + + INTEGER :: is, js + INTEGER(i8) :: ix + TYPE(particle_list), POINTER :: p_list1, p_list2 + TYPE(particle_list) :: splitted_lcs_photons, splitted_lcs_leptons + + IF (use_LCS) THEN + DO is = 1, n_species + IF (species_list(is)%species_type /= c_species_id_photon) CYCLE + DO js = 1, n_species + IF (species_list(js)%species_type == c_species_id_electron & + .OR. species_list(js)%species_type == c_species_id_positron) THEN + DO ix = 1, nx + + CALL create_empty_partlist(splitted_lcs_photons) + CALL create_empty_partlist(splitted_lcs_leptons) + + p_list1 => species_list(is)%secondary_list(ix) + p_list2 => species_list(js)%secondary_list(ix) + + CALL linear_Compton_scattering( & + p_list1, p_list2, is, js, ix, & + splitted_lcs_photons, splitted_lcs_leptons) + + IF (splitted_lcs_photons%count > 0) THEN + CALL append_partlist(species_list(is & + )%secondary_list(ix), splitted_lcs_photons) + END IF + + IF (splitted_lcs_leptons%count > 0) THEN + CALL append_partlist(species_list(js & + )%secondary_list(ix), splitted_lcs_leptons) + END IF + END DO ! do ix = 1, nx + END IF ! js being lepton + END DO ! js + END DO ! is + END IF ! if use_LCS + + END SUBROUTINE do_binary_collisions + + + + SUBROUTINE linear_Compton_scattering(p_list_i, p_list_j, & + ispe, jspe, ixx, & + splitted_phot_list, splitted_lept_list) + + !!! i is always photon, j is always lepton + + TYPE(particle_list), INTENT(INOUT) :: p_list_i, p_list_j + INTEGER, INTENT(IN) :: ispe, jspe + INTEGER(i8), INTENT(IN) :: ixx + TYPE(particle_list), INTENT(INOUT) :: splitted_phot_list + TYPE(particle_list), INTENT(INOUT) :: splitted_lept_list + + INTEGER(I8) :: icount, jcount + REAL(num) :: q_i, q_j, P_max, N_max + INTEGER :: N_coll + + TYPE(particle), POINTER :: current_i, current_j + REAL(num) :: i_Pmax + INTEGER :: N_scattered + + REAL(num) :: weight_i, weight_j + + REAL(num), DIMENSION(3) :: p_phot_lab_si, p_lept_lab_si + REAL(num), DIMENSION(3) :: p_phot_lab , p_lept_lab + REAL(num) :: gamma_lept_lab, energy_phot_lab, energy_phot_0 + REAL(num) :: sigma_rest, sigma_lab + REAL(num) :: P_coll + + REAL(num), DIMENSION(3) :: n_v, p_phot_0_si + REAL(num) :: rand_phi, rand_mu + REAL(num) :: energy_phot_0_sca + REAL(num), DIMENSION(3) :: e1, e2, e3 + REAL(num), DIMENSION(3) :: p_phot_0_sca + + REAL(num), DIMENSION(3) :: p_phot_lab_sca, p_phot_lab_sca_si + REAL(num) :: energy_phot_lab_sca, energy_phot_lab_sca_si + + REAL(num), DIMENSION(3) :: p_lept_lab_sca_si + + TYPE(particle), POINTER :: splitted_particle + + + ! If there aren't enough particles to collide, then don't bother + icount = p_list_i%count + jcount = p_list_j%count + IF (icount < 1 .OR. jcount < 1) RETURN + + !!! Determine how many macro-macro pairings are up to be checked + + ! maximal possible collisional probability P_max + q_i = max_weight(p_list_i) + q_j = max_weight(p_list_j) + P_max = sig2cdt_dV_lcs * MAX(q_i, q_j) + + ! determine how many macro-particle pairs are up to collide (N_coll) + N_max = P_max * icount * jcount + + If (random()< (N_max-FLOOR(N_max))) THEN + N_coll = CEILING(N_max) + ELSE + N_coll = FLOOR(N_max) + END IF + + IF (N_coll > MIN(icount, jcount)) THEN + PRINT*, "Too many LCS collisions." + STOP + END IF + + !!! Check the collision of these N_coll pairings of macro-photon + + IF (N_coll > 0) THEN + ! shuffle particle list only if there are pairs to check + CALL shuffle_particle_list_random(p_list_i) + CALL shuffle_particle_list_random(p_list_j) + + current_i => p_list_i%head + current_j => p_list_j%head + + i_Pmax = 1.0_num/P_max + N_scattered = 0 + + ELSE + RETURN ! N_coll = 0, no pair to check + + END IF + + DO WHILE (N_scattered < N_coll) + + !!! calculate joint collisional probability P_coll + + weight_i = current_i%weight + weight_j = current_j%weight + + ! particle momentum in lab frame in S.I. + p_phot_lab_si = current_i%part_p + p_lept_lab_si = current_j%part_p + ! (notice in this case, next_i and next_j are not defined) + + + ! particle momentum in lab frame norm. by mc + p_phot_lab = p_phot_lab_si * inv_mc0 + p_lept_lab = p_lept_lab_si * inv_mc0 + + ! lorentz factor of lepton + gamma_lept_lab = SQRT(1.0_num + DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon energy in lab frame norm. by mc2 + energy_phot_lab = current_i%particle_energy * inv_m0c2 + + ! photon energy in lepton rest frame norm. by mc2 + energy_phot_0 = gamma_lept_lab * energy_phot_lab & + - DOT_PRODUCT(p_phot_lab, p_lept_lab) + + ! compton cross section in lepton rest frame + sigma_rest = lcs_cross_sec(energy_phot_0) + + ! compton cross section in lab frame + sigma_lab = sigma_rest * energy_phot_0 & + / energy_phot_lab / gamma_lept_lab + + ! collisional probability after modification by P_max + P_coll = sigma_lab * cdt_dV * MAX(weight_i, weight_j) * i_Pmax + + If (random() < P_coll) THEN + + !!! Now, scatter these two macro-particles. + + ! unit vector of lepton velocity in lab frame + n_v = p_lept_lab / SQRT(DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon momentum in rest frame in SI + p_phot_0_si = p_phot_lab_si + (gamma_lept_lab-1.0_num) * & + DOT_PRODUCT(p_phot_lab_si,n_v)*n_v & + - p_lept_lab_si*energy_phot_lab + + IF (use_LCS_diff) THEN + + ! random azimuthal angle in c.o.m. + rand_phi = 2.0_num * pi * random() + + ! random polar angle in c.o.m. (cosine of this random angle) + rand_mu = random_polar_lcs(energy_phot_0, sigma_rest) + + ELSE + ! uniform distribution on sphere surface + rand_phi = 2.0_num * pi * random() + rand_mu = 2.0_num * random() - 1.0_num + + ! Notice unlike LBW and EPA, for LCS, there's one-to-one correspondance + ! between polar angle and scattered photon energy. + ! Therefore, the calculated scattered photon energy is only valid + ! for that particular polar angle w.r.t. collisional axis. + ! So, we still need to calculate the collisional axis (even though + ! the distribution is uniform). We cannot just assign a random + ! angle uniform in S^2, while using energy_phot_0_sca as its magnitude. + + END IF + + ! scattered photon energy in rest frame norm. by mc2 + energy_phot_0_sca = energy_phot_0 & + / (1.0_num+energy_phot_0*(1.0_num-rand_mu)) + + ! orthonormal basis (e1, e2, e3) s.t. e1//photon momentum in rest frame + CALL get_orthonormal(p_phot_0_si, e1, e2, e3) + + ! scattered photon momentum in rest frame norm. by mc + p_phot_0_sca = energy_phot_0_sca * & + ( e1*rand_mu & + + e2*SQRT(1.0_num-rand_mu**2)*COS(rand_phi) & + + e3*SQRT(1.0_num-rand_mu**2)*SIN(rand_phi)) + + !!! scattered photon momentum and energy in lab frame + ! Here the plus sign is not a typo: transform w.r.t. -n_v + p_phot_lab_sca = p_phot_0_sca + (gamma_lept_lab-1.0_num) & + * DOT_PRODUCT(p_phot_0_sca,n_v)*n_v & + + p_lept_lab*energy_phot_0_sca + + energy_phot_lab_sca = SQRT(DOT_PRODUCT(p_phot_lab_sca, & + p_phot_lab_sca)) + + p_phot_lab_sca_si = p_phot_lab_sca * mc0 + energy_phot_lab_sca_si = energy_phot_lab_sca * m0c2 + + ! scattered lepton momentum in lab frame + + p_lept_lab_sca_si = p_lept_lab_si + p_phot_lab_si & + - p_phot_lab_sca_si + + + !!! Now, Split particle if needed + + IF ((weight_i/weight_j) >1.000001_num) THEN ! photon is larger + + current_i%weight = weight_j + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_i - weight_j + splitted_particle%part_pos = current_i%part_pos + splitted_particle%part_p = current_i%part_p + splitted_particle%particle_energy = current_i%particle_energy + splitted_particle%optical_depth = current_i%optical_depth + CALL add_particle_to_partlist(splitted_phot_list, splitted_particle) + + ELSE IF ((weight_j/weight_i) >1.000001_num) THEN ! lepton is larger + + current_j%weight = weight_i + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_j - weight_i + splitted_particle%part_pos = current_j%part_pos + splitted_particle%part_p = current_j%part_p + splitted_particle%optical_depth = current_j%optical_depth + CALL add_particle_to_partlist(splitted_lept_list, splitted_particle) + + END IF + + + !!! Update particle momentum + + current_i%part_p = p_phot_lab_sca_si + current_i%particle_energy = energy_phot_lab_sca_si + + current_j%part_p = p_lept_lab_sca_si + + END IF ! random() < P_coll + + ! Scattered finished, move pointer to next particle, increment counter + current_i => current_i%next + current_j => current_j%next + N_scattered = N_scattered + 1 + + END DO ! While more to scatter + + + END SUBROUTINE linear_Compton_scattering + + + + FUNCTION max_weight(p_list) + + TYPE(particle_list), INTENT(IN) :: p_list + TYPE(particle), POINTER :: current + REAL(num) :: max_weight + + max_weight = 0.0_num + + current => p_list%head + + DO WHILE (ASSOCIATED(current)) + + max_weight = MAX(max_weight, current%weight) + current => current%next + + END DO + + END FUNCTION max_weight + + + SUBROUTINE get_orthonormal(vec, e1, e2, e3) + ! output a set of 3D orthonormal basis (e1,e2,e3) s.t. e1 // vec + + REAL(num), DIMENSION(3), INTENT(IN) :: vec + REAL(num), DIMENSION(3), INTENT(OUT) :: e1, e2, e3 + + + e1 = vec/ SQRT(DOT_PRODUCT(vec, vec)) + + IF (ABS(e1(2))<=c_tiny .AND. ABS(e1(3))<=c_tiny) THEN + ! if e1 // (1,0,0), cross product e1 with (0,1,0) to get e2 + e2(1) = -e1(3) + e2(2) = 0 + e2(3) = e1(1) + ELSE + ! cross product e1 with (1,0,0) to get e2 + e2(1) = 0 + e2(2) = e1(3) + e2(3) = -e1(2) + END IF + e2 = e2 / SQRT(DOT_PRODUCT(e2, e2)) + + e3(1) = e1(2)*e2(3) - e1(3)*e2(2) + e3(2) = e1(3)*e2(1) - e1(1)*e2(3) + e3(3) = e1(1)*e2(2) - e1(2)*e2(1) + + END SUBROUTINE get_orthonormal + + + + FUNCTION lcs_cross_sec(e) + + REAL(num), INTENT(IN) :: e + REAL(num) :: lcs_cross_sec + + lcs_cross_sec = pire2 * ( (1.0_num-2.0_num/e & + -2.0_num/e**2)*LOG(1.0_num+2.0_num*e) + 0.5_num & + + 4.0_num/e - 0.5_num/(1.0_num+2.0_num*e)**2 ) / e + + END FUNCTION lcs_cross_sec + + + + FUNCTION random_polar_lcs(en, sigma) + + ! This function generates (the cosine of) a random polar angle + ! following the distribution of the differential cross section + ! of the lcs process in the c.o.m. frame, by inverse transform + ! sampling method using bisection method. + + REAL(num), INTENT(IN) :: en, sigma + + REAL(num) :: rnd_cdf + REAL(num) :: lb, ub, mid_point + REAL(num) :: cdf_err_lb, cdf_err_mp + REAL(num) :: random_polar_lcs + + rnd_cdf = random() + + ! upper bound, lower bound, and mid point of bisection method. + lb = -1.0_num + ub = 1.0_num + mid_point = 0.0_num + + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + DO WHILE (ABS(cdf_err_mp) > tolerance_cdf & + .AND. (ABS(ub-lb)>tolerance_cos_angle)) + + IF (ABS(SIGN(1.0_num, cdf_err_lb) - SIGN(1.0_num, cdf_err_mp)) & + < 0.5_num) THEN + ! lb and mp have same sign + ! meaning correct value is between mp and ub + lb = mid_point + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + ELSE + ! lb and mp have different sign + ! meaning correct value is between lb and mp + ub = mid_point + END IF + + mid_point = (lb+ub)*0.5_num + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + END DO + + random_polar_lcs = mid_point + + END FUNCTION random_polar_lcs + + + FUNCTION lcs_polar_cdf_err(mu, rnd, e, sig) + + REAL(num), INTENT(IN) :: mu, rnd, e, sig + REAL(num) :: big_bracket, k + REAL(num) :: lcs_polar_cdf_err + + k = 1.0_num+e*(1.0_num-mu) + + big_bracket = (1.0_num-2.0_num/e-2.0_num/e**2) & + * LOG(k) - 0.5_num/k**2 -((1.0_num+2.0_num*e) & + /e**2)/k + (1.0_num-mu)/e +((1.0_num+e)/e)**2 & + - 0.5_num + + lcs_polar_cdf_err = pire2 * big_bracket / e / sig - rnd + + END FUNCTION lcs_polar_cdf_err + + + FUNCTION epa_cross_sec(v) + + REAL(num), INTENT(IN) :: v + REAL(num) :: epa_cross_sec + + ! input v is dimensionless lepton velocity in c.o.m. + + epa_cross_sec = quarter_pire2 * (1.0_num - v**2) / v * & + ( (3.0_num-v**4) / v * LOG((1.0_num+v)/(1.0_num-v)) + & + 2.0_num * (v**2-2.0_num) ) + + END FUNCTION epa_cross_sec #endif END MODULE photons diff --git a/epoch1d/src/shared_data.F90 b/epoch1d/src/shared_data.F90 index 012a6572b..9f9189441 100644 --- a/epoch1d/src/shared_data.F90 +++ b/epoch1d/src/shared_data.F90 @@ -614,8 +614,16 @@ MODULE shared_data CHARACTER(LEN=string_length) :: qed_table_location LOGICAL :: use_continuous_emission = .FALSE., use_classical_emission=.FALSE. REAL(num) :: photon_sample_fraction = 1.0_num + + LOGICAL :: use_LCS = .FALSE. + LOGICAL :: use_LCS_diff = .TRUE. + + REAL(num), PARAMETER :: tolerance_cdf = 1.0e-6_num + REAL(num), PARAMETER :: tolerance_cos_angle = 1.0e-6_num + #endif LOGICAL :: use_qed = .FALSE. + LOGICAL :: use_binary_collisions = .FALSE. #ifdef BREMSSTRAHLUNG !---------------------------------------------------------------------------- diff --git a/epoch2d/Makefile b/epoch2d/Makefile index e8ddef662..338641f0c 100644 --- a/epoch2d/Makefile +++ b/epoch2d/Makefile @@ -536,7 +536,7 @@ 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 prefetch.o: prefetch.F90 shared_data.o probes.o: probes.F90 partlist.o $(SDFMOD) random_generator.o: random_generator.f90 diff --git a/epoch2d/src/constants.F90 b/epoch2d/src/constants.F90 index 8606df132..26e7db103 100644 --- a/epoch2d/src/constants.F90 +++ b/epoch2d/src/constants.F90 @@ -231,6 +231,15 @@ 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 diff --git a/epoch2d/src/deck/deck_qed_block.F90 b/epoch2d/src/deck/deck_qed_block.F90 index 09b76ea34..25a4acd04 100644 --- a/epoch2d/src/deck/deck_qed_block.F90 +++ b/epoch2d/src/deck/deck_qed_block.F90 @@ -45,6 +45,9 @@ 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 +60,7 @@ SUBROUTINE qed_deck_finalise INTEGER :: io, iu #ifdef PHOTONS LOGICAL :: exists + INTEGER :: j IF (deck_state == c_ds_first) RETURN @@ -74,6 +78,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 +107,6 @@ SUBROUTINE qed_deck_finalise CALL abort_code(c_err_pp_options_missing) END IF #endif - END SUBROUTINE qed_deck_finalise @@ -172,6 +191,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 diff --git a/epoch2d/src/epoch2d.F90 b/epoch2d/src/epoch2d.F90 index cca400047..ce0ee0de2 100644 --- a/epoch2d/src/epoch2d.F90 +++ b/epoch2d/src/epoch2d.F90 @@ -216,7 +216,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 @@ -227,7 +228,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 (collision_step .OR. coll_ion_step .OR. recombine_step) THEN + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reorder_particles_to_grid END IF @@ -242,7 +244,14 @@ PROGRAM pic IF (use_recombination) CALL run_recombination - IF (collision_step .OR. coll_ion_step .OR. recombine_step) THEN +#ifdef PHOTONS + IF (use_binary_collisions) THEN + CALL do_binary_collisions + END IF +#endif + + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reattach_particles_to_mainlist END IF END IF diff --git a/epoch2d/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90 index 04cf89b85..16380af79 100644 --- a/epoch2d/src/physics_packages/photons.F90 +++ b/epoch2d/src/physics_packages/photons.F90 @@ -16,11 +16,14 @@ MODULE photons #ifdef PHOTONS + + USE collisions USE partlist IMPLICIT NONE - SAVE + REAL(num) :: sig2cdt_dV_lcs + REAL(num) :: cdt_dV REAL(num) :: part_pos_global, gamma_global, eta_global CONTAINS @@ -67,6 +70,9 @@ SUBROUTINE setup_qed_module END DO END IF + cdt_dV = c * dt / dx / dy + sig2cdt_dV_lcs = 2.0_num * sigma_thomson * cdt_dV + END SUBROUTINE setup_qed_module @@ -933,7 +939,7 @@ SUBROUTINE generate_photon(generating_electron, iphoton, eta) samp_photon = .TRUE. IF (photon_sample_fraction < 1.0_num) THEN - samp_photon = random() < photon_sample_fraction + samp_photon = random() < photon_sample_fraction END IF ! This will only create photons that have energies above a user specified @@ -974,7 +980,7 @@ FUNCTION calculate_photon_energy(rand_seed, eta, generating_gamma) REAL(num) :: calculate_photon_energy REAL(num), INTENT(IN) :: rand_seed, eta, generating_gamma - REAL(num) :: eta_min, chi_tmp, chi_final + REAL(num) :: eta_min, chi_tmp, chi_final eta_min = 10.0_num**MINVAL(log_eta) ! In the classical case, always use the spectrum at minimum eta @@ -1420,5 +1426,421 @@ FUNCTION find_value_from_table(x_in, p_value, nx, ny, x, y, p_table) END FUNCTION find_value_from_table + + + SUBROUTINE do_binary_collisions + + INTEGER :: is, js + INTEGER(i8) :: ix, iy + TYPE(particle_list), POINTER :: p_list1, p_list2 + TYPE(particle_list) :: splitted_lcs_photons, splitted_lcs_leptons + + IF (use_LCS) THEN + DO is = 1, n_species + IF (species_list(is)%species_type /= c_species_id_photon) CYCLE + DO js = 1, n_species + IF (species_list(js)%species_type == c_species_id_electron & + .OR. species_list(js)%species_type == c_species_id_positron) THEN + DO ix = 1, nx + DO iy = 1, ny + + CALL create_empty_partlist(splitted_lcs_photons) + CALL create_empty_partlist(splitted_lcs_leptons) + + p_list1 => species_list(is)%secondary_list(ix,iy) + p_list2 => species_list(js)%secondary_list(ix,iy) + + CALL linear_Compton_scattering( & + p_list1, p_list2, is, js, ix, iy,& + splitted_lcs_photons, splitted_lcs_leptons) + + IF (splitted_lcs_photons%count > 0) THEN + CALL append_partlist(species_list(is & + )%secondary_list(ix,iy), splitted_lcs_photons) + END IF + + IF (splitted_lcs_leptons%count > 0) THEN + CALL append_partlist(species_list(js & + )%secondary_list(ix,iy), splitted_lcs_leptons) + END IF + END DO ! do iy = 1, ny + END DO ! do ix = 1, nx + END IF ! js being lepton + END DO ! js + END DO ! is + END IF ! if use_LCS + + END SUBROUTINE do_binary_collisions + + + + SUBROUTINE linear_Compton_scattering(p_list_i, p_list_j, & + ispe, jspe, ixx, iyy, & + splitted_phot_list, splitted_lept_list) + + !!! i is always photon, j is always lepton + + TYPE(particle_list), INTENT(INOUT) :: p_list_i, p_list_j + INTEGER, INTENT(IN) :: ispe, jspe + INTEGER(i8), INTENT(IN) :: ixx, iyy + TYPE(particle_list), INTENT(INOUT) :: splitted_phot_list + TYPE(particle_list), INTENT(INOUT) :: splitted_lept_list + + INTEGER(I8) :: icount, jcount + REAL(num) :: q_i, q_j, P_max, N_max + INTEGER :: N_coll + + TYPE(particle), POINTER :: current_i, current_j + REAL(num) :: i_Pmax + INTEGER :: N_scattered + + REAL(num) :: weight_i, weight_j + + REAL(num), DIMENSION(3) :: p_phot_lab_si, p_lept_lab_si + REAL(num), DIMENSION(3) :: p_phot_lab , p_lept_lab + REAL(num) :: gamma_lept_lab, energy_phot_lab, energy_phot_0 + REAL(num) :: sigma_rest, sigma_lab + REAL(num) :: P_coll + + REAL(num), DIMENSION(3) :: n_v, p_phot_0_si + REAL(num) :: rand_phi, rand_mu + REAL(num) :: energy_phot_0_sca + REAL(num), DIMENSION(3) :: e1, e2, e3 + REAL(num), DIMENSION(3) :: p_phot_0_sca + + REAL(num), DIMENSION(3) :: p_phot_lab_sca, p_phot_lab_sca_si + REAL(num) :: energy_phot_lab_sca, energy_phot_lab_sca_si + + REAL(num), DIMENSION(3) :: p_lept_lab_sca_si + + TYPE(particle), POINTER :: splitted_particle + + + ! If there aren't enough particles to collide, then don't bother + icount = p_list_i%count + jcount = p_list_j%count + IF (icount < 1 .OR. jcount < 1) RETURN + + !!! Determine how many macro-macro pairings are up to be checked + + ! maximal possible collisional probability P_max + q_i = max_weight(p_list_i) + q_j = max_weight(p_list_j) + P_max = sig2cdt_dV_lcs * MAX(q_i, q_j) + + ! determine how many macro-particle pairs are up to collide (N_coll) + N_max = P_max * icount * jcount + + If (random()< (N_max-FLOOR(N_max))) THEN + N_coll = CEILING(N_max) + ELSE + N_coll = FLOOR(N_max) + END IF + + IF (N_coll > MIN(icount, jcount)) THEN + PRINT*, "Too many LCS collisions." + STOP + END IF + + !!! Check the collision of these N_coll pairings of macro-photon + + IF (N_coll > 0) THEN + ! shuffle particle list only if there are pairs to check + CALL shuffle_particle_list_random(p_list_i) + CALL shuffle_particle_list_random(p_list_j) + + current_i => p_list_i%head + current_j => p_list_j%head + + i_Pmax = 1.0_num/P_max + N_scattered = 0 + + ELSE + RETURN ! N_coll = 0, no pair to check + + END IF + + DO WHILE (N_scattered < N_coll) + + !!! calculate joint collisional probability P_coll + + weight_i = current_i%weight + weight_j = current_j%weight + + ! particle momentum in lab frame in S.I. + p_phot_lab_si = current_i%part_p + p_lept_lab_si = current_j%part_p + ! (notice in this case, next_i and next_j are not defined) + + + ! particle momentum in lab frame norm. by mc + p_phot_lab = p_phot_lab_si * inv_mc0 + p_lept_lab = p_lept_lab_si * inv_mc0 + + ! lorentz factor of lepton + gamma_lept_lab = SQRT(1.0_num + DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon energy in lab frame norm. by mc2 + energy_phot_lab = current_i%particle_energy * inv_m0c2 + + ! photon energy in lepton rest frame norm. by mc2 + energy_phot_0 = gamma_lept_lab * energy_phot_lab & + - DOT_PRODUCT(p_phot_lab, p_lept_lab) + + ! compton cross section in lepton rest frame + sigma_rest = lcs_cross_sec(energy_phot_0) + + ! compton cross section in lab frame + sigma_lab = sigma_rest * energy_phot_0 & + / energy_phot_lab / gamma_lept_lab + + ! collisional probability after modification by P_max + P_coll = sigma_lab * cdt_dV * MAX(weight_i, weight_j) * i_Pmax + + If (random() < P_coll) THEN + + !!! Now, scatter these two macro-particles. + + ! unit vector of lepton velocity in lab frame + n_v = p_lept_lab / SQRT(DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon momentum in rest frame in SI + p_phot_0_si = p_phot_lab_si + (gamma_lept_lab-1.0_num) * & + DOT_PRODUCT(p_phot_lab_si,n_v)*n_v & + - p_lept_lab_si*energy_phot_lab + + IF (use_LCS_diff) THEN + + ! random azimuthal angle in c.o.m. + rand_phi = 2.0_num * pi * random() + + ! random polar angle in c.o.m. (cosine of this random angle) + rand_mu = random_polar_lcs(energy_phot_0, sigma_rest) + + ELSE + ! uniform distribution on sphere surface + rand_phi = 2.0_num * pi * random() + rand_mu = 2.0_num * random() - 1.0_num + + ! Notice unlike LBW and EPA, for LCS, there's one-to-one correspondance + ! between polar angle and scattered photon energy. + ! Therefore, the calculated scattered photon energy is only valid + ! for that particular polar angle w.r.t. collisional axis. + ! So, we still need to calculate the collisional axis (even though + ! the distribution is uniform). We cannot just assign a random + ! angle uniform in S^2, while using energy_phot_0_sca as its magnitude. + + END IF + + ! scattered photon energy in rest frame norm. by mc2 + energy_phot_0_sca = energy_phot_0 & + / (1.0_num+energy_phot_0*(1.0_num-rand_mu)) + + ! orthonormal basis (e1, e2, e3) s.t. e1//photon momentum in rest frame + CALL get_orthonormal(p_phot_0_si, e1, e2, e3) + + ! scattered photon momentum in rest frame norm. by mc + p_phot_0_sca = energy_phot_0_sca * & + ( e1*rand_mu & + + e2*SQRT(1.0_num-rand_mu**2)*COS(rand_phi) & + + e3*SQRT(1.0_num-rand_mu**2)*SIN(rand_phi)) + + !!! scattered photon momentum and energy in lab frame + ! Here the plus sign is not a typo: transform w.r.t. -n_v + p_phot_lab_sca = p_phot_0_sca + (gamma_lept_lab-1.0_num) & + * DOT_PRODUCT(p_phot_0_sca,n_v)*n_v & + + p_lept_lab*energy_phot_0_sca + + energy_phot_lab_sca = SQRT(DOT_PRODUCT(p_phot_lab_sca, & + p_phot_lab_sca)) + + p_phot_lab_sca_si = p_phot_lab_sca * mc0 + energy_phot_lab_sca_si = energy_phot_lab_sca * m0c2 + + ! scattered lepton momentum in lab frame + + p_lept_lab_sca_si = p_lept_lab_si + p_phot_lab_si & + - p_phot_lab_sca_si + + + !!! Now, Split particle if needed + + IF ((weight_i/weight_j) >1.000001_num) THEN ! photon is larger + + current_i%weight = weight_j + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_i - weight_j + splitted_particle%part_pos = current_i%part_pos + splitted_particle%part_p = current_i%part_p + splitted_particle%particle_energy = current_i%particle_energy + splitted_particle%optical_depth = current_i%optical_depth + CALL add_particle_to_partlist(splitted_phot_list, splitted_particle) + + ELSE IF ((weight_j/weight_i) >1.000001_num) THEN ! lepton is larger + + current_j%weight = weight_i + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_j - weight_i + splitted_particle%part_pos = current_j%part_pos + splitted_particle%part_p = current_j%part_p + splitted_particle%optical_depth = current_j%optical_depth + CALL add_particle_to_partlist(splitted_lept_list, splitted_particle) + + END IF + + + !!! Update particle momentum + + current_i%part_p = p_phot_lab_sca_si + current_i%particle_energy = energy_phot_lab_sca_si + + current_j%part_p = p_lept_lab_sca_si + + END IF ! random() < P_coll + + ! Scattered finished, move pointer to next particle, increment counter + current_i => current_i%next + current_j => current_j%next + N_scattered = N_scattered + 1 + + END DO ! While more to scatter + + + END SUBROUTINE linear_Compton_scattering + + + + FUNCTION max_weight(p_list) + + TYPE(particle_list), INTENT(IN) :: p_list + TYPE(particle), POINTER :: current + REAL(num) :: max_weight + + max_weight = 0.0_num + + current => p_list%head + + DO WHILE (ASSOCIATED(current)) + + max_weight = MAX(max_weight, current%weight) + current => current%next + + END DO + + END FUNCTION max_weight + + + + SUBROUTINE get_orthonormal(vec, e1, e2, e3) + ! output a set of 3D orthonormal basis (e1,e2,e3) s.t. e1 // vec + + REAL(num), DIMENSION(3), INTENT(IN) :: vec + REAL(num), DIMENSION(3), INTENT(OUT) :: e1, e2, e3 + + + e1 = vec/ SQRT(DOT_PRODUCT(vec, vec)) + + IF (ABS(e1(2))<=c_tiny .AND. ABS(e1(3))<=c_tiny) THEN + ! if e1 // (1,0,0), cross product e1 with (0,1,0) to get e2 + e2(1) = -e1(3) + e2(2) = 0 + e2(3) = e1(1) + ELSE + ! cross product e1 with (1,0,0) to get e2 + e2(1) = 0 + e2(2) = e1(3) + e2(3) = -e1(2) + END IF + e2 = e2 / SQRT(DOT_PRODUCT(e2, e2)) + + e3(1) = e1(2)*e2(3) - e1(3)*e2(2) + e3(2) = e1(3)*e2(1) - e1(1)*e2(3) + e3(3) = e1(1)*e2(2) - e1(2)*e2(1) + + END SUBROUTINE get_orthonormal + + + + FUNCTION lcs_cross_sec(e) + + REAL(num), INTENT(IN) :: e + REAL(num) :: lcs_cross_sec + + lcs_cross_sec = pire2 * ( (1.0_num-2.0_num/e & + -2.0_num/e**2)*LOG(1.0_num+2.0_num*e) + 0.5_num & + + 4.0_num/e - 0.5_num/(1.0_num+2.0_num*e)**2 ) / e + + END FUNCTION lcs_cross_sec + + + + FUNCTION random_polar_lcs(en, sigma) + + ! This function generates (the cosine of) a random polar angle + ! following the distribution of the differential cross section + ! of the lcs process in the c.o.m. frame, by inverse transform + ! sampling method using bisection method. + + REAL(num), INTENT(IN) :: en, sigma + + REAL(num) :: rnd_cdf + REAL(num) :: lb, ub, mid_point + REAL(num) :: cdf_err_lb, cdf_err_mp + REAL(num) :: random_polar_lcs + + rnd_cdf = random() + + ! upper bound, lower bound, and mid point of bisection method. + lb = -1.0_num + ub = 1.0_num + mid_point = 0.0_num + + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + DO WHILE (ABS(cdf_err_mp) > tolerance_cdf & + .AND. (ABS(ub-lb)>tolerance_cos_angle)) + + IF (ABS(SIGN(1.0_num, cdf_err_lb) - SIGN(1.0_num, cdf_err_mp)) & + < 0.5_num) THEN + ! lb and mp have same sign + ! meaning correct value is between mp and ub + lb = mid_point + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + ELSE + ! lb and mp have different sign + ! meaning correct value is between lb and mp + ub = mid_point + END IF + + mid_point = (lb+ub)*0.5_num + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + END DO + + random_polar_lcs = mid_point + + END FUNCTION random_polar_lcs + + + + FUNCTION lcs_polar_cdf_err(mu, rnd, e, sig) + + REAL(num), INTENT(IN) :: mu, rnd, e, sig + REAL(num) :: big_bracket, k + REAL(num) :: lcs_polar_cdf_err + + k = 1.0_num+e*(1.0_num-mu) + + big_bracket = (1.0_num-2.0_num/e-2.0_num/e**2) & + * LOG(k) - 0.5_num/k**2 -((1.0_num+2.0_num*e) & + /e**2)/k + (1.0_num-mu)/e +((1.0_num+e)/e)**2 & + - 0.5_num + + lcs_polar_cdf_err = pire2 * big_bracket / e / sig - rnd + + END FUNCTION lcs_polar_cdf_err #endif END MODULE photons diff --git a/epoch2d/src/shared_data.F90 b/epoch2d/src/shared_data.F90 index 5c1678e5a..8571376c3 100644 --- a/epoch2d/src/shared_data.F90 +++ b/epoch2d/src/shared_data.F90 @@ -637,8 +637,15 @@ MODULE shared_data CHARACTER(LEN=string_length) :: qed_table_location LOGICAL :: use_continuous_emission = .FALSE., use_classical_emission=.FALSE. REAL(num) :: photon_sample_fraction = 1.0_num + + LOGICAL :: use_LCS = .FALSE. + LOGICAL :: use_LCS_diff = .TRUE. + REAL(num), PARAMETER :: tolerance_cdf = 1.0e-6_num + REAL(num), PARAMETER :: tolerance_cos_angle = 1.0e-6_num + #endif LOGICAL :: use_qed = .FALSE. + LOGICAL :: use_binary_collisions = .FALSE. #ifdef BREMSSTRAHLUNG !---------------------------------------------------------------------------- diff --git a/epoch3d/Makefile b/epoch3d/Makefile index 837cdde5f..4b9c5bbee 100644 --- a/epoch3d/Makefile +++ b/epoch3d/Makefile @@ -536,7 +536,7 @@ 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 prefetch.o: prefetch.F90 shared_data.o probes.o: probes.F90 partlist.o $(SDFMOD) random_generator.o: random_generator.f90 diff --git a/epoch3d/src/constants.F90 b/epoch3d/src/constants.F90 index 18ea9537c..2a15ccd37 100644 --- a/epoch3d/src/constants.F90 +++ b/epoch3d/src/constants.F90 @@ -232,6 +232,15 @@ 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 diff --git a/epoch3d/src/deck/deck_qed_block.F90 b/epoch3d/src/deck/deck_qed_block.F90 index 09b76ea34..5b0991d6a 100644 --- a/epoch3d/src/deck/deck_qed_block.F90 +++ b/epoch3d/src/deck/deck_qed_block.F90 @@ -45,6 +45,9 @@ 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 +60,7 @@ SUBROUTINE qed_deck_finalise INTEGER :: io, iu #ifdef PHOTONS LOGICAL :: exists + INTEGER :: j IF (deck_state == c_ds_first) RETURN @@ -74,6 +78,23 @@ 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 diff --git a/epoch3d/src/epoch3d.F90 b/epoch3d/src/epoch3d.F90 index 67754908e..b745287c2 100644 --- a/epoch3d/src/epoch3d.F90 +++ b/epoch3d/src/epoch3d.F90 @@ -216,7 +216,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 @@ -227,7 +228,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 (collision_step .OR. coll_ion_step .OR. recombine_step) THEN + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reorder_particles_to_grid END IF @@ -242,7 +244,14 @@ PROGRAM pic IF (recombine_step) CALL run_recombination - IF (collision_step .OR. coll_ion_step .OR. recombine_step) THEN +#ifdef PHOTONS + IF (use_binary_collisions) THEN + CALL do_binary_collisions + END IF +#endif + + IF (collision_step .OR. coll_ion_step .OR. recombine_step & + .OR. use_binary_collisions) THEN CALL reattach_particles_to_mainlist END IF END IF diff --git a/epoch3d/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90 index db5804098..fc44d85ea 100644 --- a/epoch3d/src/physics_packages/photons.F90 +++ b/epoch3d/src/physics_packages/photons.F90 @@ -16,11 +16,14 @@ MODULE photons #ifdef PHOTONS + + USE collisions USE partlist IMPLICIT NONE - SAVE + REAL(num) :: sig2cdt_dV_lcs + REAL(num) :: cdt_dV REAL(num) :: part_pos_global, gamma_global, eta_global CONTAINS @@ -67,6 +70,9 @@ SUBROUTINE setup_qed_module END DO END IF + cdt_dV = c * dt / dx/dy/dz + sig2cdt_dV_lcs = 2.0_num * sigma_thomson * cdt_dV + END SUBROUTINE setup_qed_module @@ -946,7 +952,7 @@ SUBROUTINE generate_photon(generating_electron, iphoton, eta) samp_photon = .TRUE. IF (photon_sample_fraction < 1.0_num) THEN - samp_photon = random() < photon_sample_fraction + samp_photon = random() < photon_sample_fraction END IF ! This will only create photons that have energies above a user specified @@ -969,7 +975,7 @@ SUBROUTINE generate_photon(generating_electron, iphoton, eta) IF (use_continuous_emission) THEN ! Calculate photon weight from synchrotron emission rate sync_rate = delta_optical_depth(eta, generating_gamma) ! This already has *dt included - new_photon%weight = generating_electron%weight * sync_rate & + new_photon%weight = generating_electron%weight * sync_rate & / photon_sample_fraction ELSE new_photon%weight = generating_electron%weight / photon_sample_fraction @@ -987,7 +993,7 @@ FUNCTION calculate_photon_energy(rand_seed, eta, generating_gamma) REAL(num) :: calculate_photon_energy REAL(num), INTENT(IN) :: rand_seed, eta, generating_gamma - REAL(num) :: eta_min, chi_tmp, chi_final + REAL(num) :: eta_min, chi_tmp, chi_final eta_min = 10.0_num**MINVAL(log_eta) ! In the classical case, always use the spectrum at minimum eta @@ -1433,5 +1439,422 @@ FUNCTION find_value_from_table(x_in, p_value, nx, ny, x, y, p_table) END FUNCTION find_value_from_table + + + SUBROUTINE do_binary_collisions + + INTEGER :: is, js + INTEGER(i8) :: ix, iy, iz + TYPE(particle_list), POINTER :: p_list1, p_list2 + TYPE(particle_list) :: splitted_lcs_photons, splitted_lcs_leptons + + IF (use_LCS) THEN + DO is = 1, n_species + IF (species_list(is)%species_type /= c_species_id_photon) CYCLE + DO js = 1, n_species + IF (species_list(js)%species_type == c_species_id_electron & + .OR. species_list(js)%species_type == c_species_id_positron) THEN + DO ix = 1, nx + DO iy = 1, ny + DO iz = 1, nz + + CALL create_empty_partlist(splitted_lcs_photons) + CALL create_empty_partlist(splitted_lcs_leptons) + + p_list1 => species_list(is)%secondary_list(ix,iy,iz) + p_list2 => species_list(js)%secondary_list(ix,iy,iz) + + CALL linear_Compton_scattering( & + p_list1, p_list2, is, js, ix, iy, iz, & + splitted_lcs_photons, splitted_lcs_leptons) + + IF (splitted_lcs_photons%count > 0) THEN + CALL append_partlist(species_list(is & + )%secondary_list(ix,iy,iz), splitted_lcs_photons) + END IF + + IF (splitted_lcs_leptons%count > 0) THEN + CALL append_partlist(species_list(js & + )%secondary_list(ix,iy,iz), splitted_lcs_leptons) + END IF + END DO ! do iz = 1, nz + END DO ! do iy = 1, ny + END DO ! do ix = 1, nx + END IF ! js being lepton + END DO ! js + END DO ! is + END IF ! if use_LCS + + END SUBROUTINE do_binary_collisions + + + + SUBROUTINE linear_Compton_scattering(p_list_i, p_list_j, & + ispe, jspe, ixx, iyy, izz, & + splitted_phot_list, splitted_lept_list) + + !!! i is always photon, j is always lepton + + TYPE(particle_list), INTENT(INOUT) :: p_list_i, p_list_j + INTEGER, INTENT(IN) :: ispe, jspe + INTEGER(i8), INTENT(IN) :: ixx, iyy, izz + TYPE(particle_list), INTENT(INOUT) :: splitted_phot_list + TYPE(particle_list), INTENT(INOUT) :: splitted_lept_list + + INTEGER(I8) :: icount, jcount + REAL(num) :: q_i, q_j, P_max, N_max + INTEGER :: N_coll + + TYPE(particle), POINTER :: current_i, current_j + REAL(num) :: i_Pmax + INTEGER :: N_scattered + + REAL(num) :: weight_i, weight_j + + REAL(num), DIMENSION(3) :: p_phot_lab_si, p_lept_lab_si + REAL(num), DIMENSION(3) :: p_phot_lab , p_lept_lab + REAL(num) :: gamma_lept_lab, energy_phot_lab, energy_phot_0 + REAL(num) :: sigma_rest, sigma_lab + REAL(num) :: P_coll + + REAL(num), DIMENSION(3) :: n_v, p_phot_0_si + REAL(num) :: rand_phi, rand_mu + REAL(num) :: energy_phot_0_sca + REAL(num), DIMENSION(3) :: e1, e2, e3 + REAL(num), DIMENSION(3) :: p_phot_0_sca + + REAL(num), DIMENSION(3) :: p_phot_lab_sca, p_phot_lab_sca_si + REAL(num) :: energy_phot_lab_sca, energy_phot_lab_sca_si + + REAL(num), DIMENSION(3) :: p_lept_lab_sca_si + + TYPE(particle), POINTER :: splitted_particle + + + ! If there aren't enough particles to collide, then don't bother + icount = p_list_i%count + jcount = p_list_j%count + IF (icount < 1 .OR. jcount < 1) RETURN + + !!! Determine how many macro-macro pairings are up to be checked + + ! maximal possible collisional probability P_max + q_i = max_weight(p_list_i) + q_j = max_weight(p_list_j) + P_max = sig2cdt_dV_lcs * MAX(q_i, q_j) + + ! determine how many macro-particle pairs are up to collide (N_coll) + N_max = P_max * icount * jcount + + If (random()< (N_max-FLOOR(N_max))) THEN + N_coll = CEILING(N_max) + ELSE + N_coll = FLOOR(N_max) + END IF + + IF (N_coll > MIN(icount, jcount)) THEN + PRINT*, "Too many LCS collisions." + STOP + END IF + + !!! Check the collision of these N_coll pairings of macro-photon + + IF (N_coll > 0) THEN + ! shuffle particle list only if there are pairs to check + CALL shuffle_particle_list_random(p_list_i) + CALL shuffle_particle_list_random(p_list_j) + + current_i => p_list_i%head + current_j => p_list_j%head + + i_Pmax = 1.0_num/P_max + N_scattered = 0 + + ELSE + RETURN ! N_coll = 0, no pair to check + + END IF + + DO WHILE (N_scattered < N_coll) + + !!! calculate joint collisional probability P_coll + + weight_i = current_i%weight + weight_j = current_j%weight + + ! particle momentum in lab frame in S.I. + p_phot_lab_si = current_i%part_p + p_lept_lab_si = current_j%part_p + ! (notice in this case, next_i and next_j are not defined) + + + ! particle momentum in lab frame norm. by mc + p_phot_lab = p_phot_lab_si * inv_mc0 + p_lept_lab = p_lept_lab_si * inv_mc0 + + ! lorentz factor of lepton + gamma_lept_lab = SQRT(1.0_num + DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon energy in lab frame norm. by mc2 + energy_phot_lab = current_i%particle_energy * inv_m0c2 + + ! photon energy in lepton rest frame norm. by mc2 + energy_phot_0 = gamma_lept_lab * energy_phot_lab & + - DOT_PRODUCT(p_phot_lab, p_lept_lab) + + ! compton cross section in lepton rest frame + sigma_rest = lcs_cross_sec(energy_phot_0) + + ! compton cross section in lab frame + sigma_lab = sigma_rest * energy_phot_0 & + / energy_phot_lab / gamma_lept_lab + + ! collisional probability after modification by P_max + P_coll = sigma_lab * cdt_dV * MAX(weight_i, weight_j) * i_Pmax + + If (random() < P_coll) THEN + + !!! Now, scatter these two macro-particles. + + ! unit vector of lepton velocity in lab frame + n_v = p_lept_lab / SQRT(DOT_PRODUCT(p_lept_lab, p_lept_lab)) + + ! photon momentum in rest frame in SI + p_phot_0_si = p_phot_lab_si + (gamma_lept_lab-1.0_num) * & + DOT_PRODUCT(p_phot_lab_si,n_v)*n_v & + - p_lept_lab_si*energy_phot_lab + + IF (use_LCS_diff) THEN + + ! random azimuthal angle in c.o.m. + rand_phi = 2.0_num * pi * random() + + ! random polar angle in c.o.m. (cosine of this random angle) + rand_mu = random_polar_lcs(energy_phot_0, sigma_rest) + + ELSE + ! uniform distribution on sphere surface + rand_phi = 2.0_num * pi * random() + rand_mu = 2.0_num * random() - 1.0_num + + ! Notice unlike LBW and EPA, for LCS, there's one-to-one correspondance + ! between polar angle and scattered photon energy. + ! Therefore, the calculated scattered photon energy is only valid + ! for that particular polar angle w.r.t. collisional axis. + ! So, we still need to calculate the collisional axis (even though + ! the distribution is uniform). We cannot just assign a random + ! angle uniform in S^2, while using energy_phot_0_sca as its magnitude. + + END IF + + ! scattered photon energy in rest frame norm. by mc2 + energy_phot_0_sca = energy_phot_0 & + / (1.0_num+energy_phot_0*(1.0_num-rand_mu)) + + ! orthonormal basis (e1, e2, e3) s.t. e1//photon momentum in rest frame + CALL get_orthonormal(p_phot_0_si, e1, e2, e3) + + ! scattered photon momentum in rest frame norm. by mc + p_phot_0_sca = energy_phot_0_sca * & + ( e1*rand_mu & + + e2*SQRT(1.0_num-rand_mu**2)*COS(rand_phi) & + + e3*SQRT(1.0_num-rand_mu**2)*SIN(rand_phi)) + + !!! scattered photon momentum and energy in lab frame + ! Here the plus sign is not a typo: transform w.r.t. -n_v + p_phot_lab_sca = p_phot_0_sca + (gamma_lept_lab-1.0_num) & + * DOT_PRODUCT(p_phot_0_sca,n_v)*n_v & + + p_lept_lab*energy_phot_0_sca + + energy_phot_lab_sca = SQRT(DOT_PRODUCT(p_phot_lab_sca, & + p_phot_lab_sca)) + + p_phot_lab_sca_si = p_phot_lab_sca * mc0 + energy_phot_lab_sca_si = energy_phot_lab_sca * m0c2 + + ! scattered lepton momentum in lab frame + + p_lept_lab_sca_si = p_lept_lab_si + p_phot_lab_si & + - p_phot_lab_sca_si + + + !!! Now, Split particle if needed + + IF ((weight_i/weight_j) >1.000001_num) THEN ! photon is larger + + current_i%weight = weight_j + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_i - weight_j + splitted_particle%part_pos = current_i%part_pos + splitted_particle%part_p = current_i%part_p + splitted_particle%particle_energy = current_i%particle_energy + splitted_particle%optical_depth = current_i%optical_depth + CALL add_particle_to_partlist(splitted_phot_list, splitted_particle) + + ELSE IF ((weight_j/weight_i) >1.000001_num) THEN ! lepton is larger + + current_j%weight = weight_i + + CALL create_particle(splitted_particle) + splitted_particle%weight = weight_j - weight_i + splitted_particle%part_pos = current_j%part_pos + splitted_particle%part_p = current_j%part_p + splitted_particle%optical_depth = current_j%optical_depth + CALL add_particle_to_partlist(splitted_lept_list, splitted_particle) + + END IF + + + !!! Update particle momentum + + current_i%part_p = p_phot_lab_sca_si + current_i%particle_energy = energy_phot_lab_sca_si + + current_j%part_p = p_lept_lab_sca_si + + END IF ! random() < P_coll + + ! Scattered finished, move pointer to next particle, increment counter + current_i => current_i%next + current_j => current_j%next + N_scattered = N_scattered + 1 + + END DO ! While more to scatter + + + END SUBROUTINE linear_Compton_scattering + + + + FUNCTION max_weight(p_list) + + TYPE(particle_list), INTENT(IN) :: p_list + TYPE(particle), POINTER :: current + REAL(num) :: max_weight + + max_weight = 0.0_num + + current => p_list%head + + DO WHILE (ASSOCIATED(current)) + + max_weight = MAX(max_weight, current%weight) + current => current%next + + END DO + + END FUNCTION max_weight + + + + SUBROUTINE get_orthonormal(vec, e1, e2, e3) + ! output a set of 3D orthonormal basis (e1,e2,e3) s.t. e1 // vec + + REAL(num), DIMENSION(3), INTENT(IN) :: vec + REAL(num), DIMENSION(3), INTENT(OUT) :: e1, e2, e3 + + + e1 = vec/ SQRT(DOT_PRODUCT(vec, vec)) + + IF (ABS(e1(2))<=c_tiny .AND. ABS(e1(3))<=c_tiny) THEN + ! if e1 // (1,0,0), cross product e1 with (0,1,0) to get e2 + e2(1) = -e1(3) + e2(2) = 0 + e2(3) = e1(1) + ELSE + ! cross product e1 with (1,0,0) to get e2 + e2(1) = 0 + e2(2) = e1(3) + e2(3) = -e1(2) + END IF + e2 = e2 / SQRT(DOT_PRODUCT(e2, e2)) + + e3(1) = e1(2)*e2(3) - e1(3)*e2(2) + e3(2) = e1(3)*e2(1) - e1(1)*e2(3) + e3(3) = e1(1)*e2(2) - e1(2)*e2(1) + + END SUBROUTINE get_orthonormal + + + + FUNCTION lcs_cross_sec(e) + + REAL(num), INTENT(IN) :: e + REAL(num) :: lcs_cross_sec + + lcs_cross_sec = pire2 * ( (1.0_num-2.0_num/e & + -2.0_num/e**2)*LOG(1.0_num+2.0_num*e) + 0.5_num & + + 4.0_num/e - 0.5_num/(1.0_num+2.0_num*e)**2 ) / e + + END FUNCTION lcs_cross_sec + + + + FUNCTION random_polar_lcs(en, sigma) + + ! This function generates (the cosine of) a random polar angle + ! following the distribution of the differential cross section + ! of the lcs process in the c.o.m. frame, by inverse transform + ! sampling method using bisection method. + + REAL(num), INTENT(IN) :: en, sigma + + REAL(num) :: rnd_cdf + REAL(num) :: lb, ub, mid_point + REAL(num) :: cdf_err_lb, cdf_err_mp + REAL(num) :: random_polar_lcs + + rnd_cdf = random() + + ! upper bound, lower bound, and mid point of bisection method. + lb = -1.0_num + ub = 1.0_num + mid_point = 0.0_num + + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + DO WHILE (ABS(cdf_err_mp) > tolerance_cdf & + .AND. (ABS(ub-lb)>tolerance_cos_angle)) + + IF (ABS(SIGN(1.0_num, cdf_err_lb) - SIGN(1.0_num, cdf_err_mp)) & + < 0.5_num) THEN + ! lb and mp have same sign + ! meaning correct value is between mp and ub + lb = mid_point + cdf_err_lb = lcs_polar_cdf_err(lb, rnd_cdf, en, sigma) + ELSE + ! lb and mp have different sign + ! meaning correct value is between lb and mp + ub = mid_point + END IF + + mid_point = (lb+ub)*0.5_num + cdf_err_mp = lcs_polar_cdf_err(mid_point, rnd_cdf, en, sigma) + + END DO + + random_polar_lcs = mid_point + + END FUNCTION random_polar_lcs + + + FUNCTION lcs_polar_cdf_err(mu, rnd, e, sig) + + REAL(num), INTENT(IN) :: mu, rnd, e, sig + REAL(num) :: big_bracket, k + REAL(num) :: lcs_polar_cdf_err + + k = 1.0_num+e*(1.0_num-mu) + + big_bracket = (1.0_num-2.0_num/e-2.0_num/e**2) & + * LOG(k) - 0.5_num/k**2 -((1.0_num+2.0_num*e) & + /e**2)/k + (1.0_num-mu)/e +((1.0_num+e)/e)**2 & + - 0.5_num + + lcs_polar_cdf_err = pire2 * big_bracket / e / sig - rnd + + END FUNCTION lcs_polar_cdf_err #endif END MODULE photons diff --git a/epoch3d/src/shared_data.F90 b/epoch3d/src/shared_data.F90 index 190663e33..0eb897b2e 100644 --- a/epoch3d/src/shared_data.F90 +++ b/epoch3d/src/shared_data.F90 @@ -659,8 +659,15 @@ MODULE shared_data CHARACTER(LEN=string_length) :: qed_table_location LOGICAL :: use_continuous_emission = .FALSE., use_classical_emission=.FALSE. REAL(num) :: photon_sample_fraction = 1.0_num + + LOGICAL :: use_LCS = .FALSE. + LOGICAL :: use_LCS_diff = .TRUE. + REAL(num), PARAMETER :: tolerance_cdf = 1.0e-6_num + REAL(num), PARAMETER :: tolerance_cos_angle = 1.0e-6_num #endif + LOGICAL :: use_qed = .FALSE. + LOGICAL :: use_binary_collisions = .FALSE. #ifdef BREMSSTRAHLUNG !----------------------------------------------------------------------------