diff --git a/epoch1d/Makefile b/epoch1d/Makefile index dc633fe65..790450480 100644 --- a/epoch1d/Makefile +++ b/epoch1d/Makefile @@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.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 + partlist.o photons.o laser.o: laser.f90 custom_laser.o evaluate.o mpi_routines.o: mpi_routines.F90 helper.o mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o diff --git a/epoch1d/src/constants.F90 b/epoch1d/src/constants.F90 index 46bfba58a..56ed134e4 100644 --- a/epoch1d/src/constants.F90 +++ b/epoch1d/src/constants.F90 @@ -651,7 +651,8 @@ MODULE constants 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_dump_part_qed_chi = 79 + INTEGER, PARAMETER :: num_vars_to_dump = 79 INTEGER, PARAMETER :: c_subset_random = 1 INTEGER, PARAMETER :: c_subset_gamma_min = 2 diff --git a/epoch1d/src/deck/deck_io_block.F90 b/epoch1d/src/deck/deck_io_block.F90 index d0952b3dd..6d6a43f0c 100644 --- a/epoch1d/src/deck/deck_io_block.F90 +++ b/epoch1d/src/deck/deck_io_block.F90 @@ -574,6 +574,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode) #ifdef PHOTONS ELSE IF (str_cmp(element, 'optical_depth')) THEN elementselected = c_dump_part_opdepth + + ELSE IF (str_cmp(element, 'qed_chi')) THEN + elementselected = c_dump_part_qed_chi #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch1d/src/io/diagnostics.F90 b/epoch1d/src/io/diagnostics.F90 index 85a3e1d94..53577bedb 100644 --- a/epoch1d/src/io/diagnostics.F90 +++ b/epoch1d/src/io/diagnostics.F90 @@ -635,6 +635,9 @@ SUBROUTINE output_routines(step, force_write) ! step = step index #ifdef PHOTONS CALL write_particle_variable(c_dump_part_opdepth, code, & 'Optical depth', '', it_output_real) + + CALL write_particle_variable(c_dump_part_qed_chi, code, & + 'Chi', '', it_output_real) #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) CALL write_particle_variable(c_dump_part_qed_energy, code, & diff --git a/epoch1d/src/io/iterators.F90 b/epoch1d/src/io/iterators.F90 index 52654ad7a..91cd776c9 100644 --- a/epoch1d/src/io/iterators.F90 +++ b/epoch1d/src/io/iterators.F90 @@ -18,6 +18,7 @@ MODULE iterators USE particle_pointer_advance USE partlist USE particle_id_hash_mod + USE photons IMPLICIT NONE @@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param) TYPE(particle_list), POINTER, SAVE :: current_list INTEGER :: part_count, ndim REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge - +#ifdef PHOTONS + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x +#endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) END IF @@ -346,6 +349,27 @@ FUNCTION it_output_real(array, npoint_it, start, param) array(part_count) = cur%optical_depth cur => cur%next END DO + CASE (c_dump_part_qed_chi) + IF (current_species%species_type == c_species_id_photon) THEN + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_x = cur%part_pos - x_grid_min_local + norm = c / cur%particle_energy + dir_x = cur%part_p(1) * norm + dir_y = cur%part_p(2) * norm + dir_z = cur%part_p(3) * norm + part_e = cur%particle_energy / m0 / c**2 + part_count = part_count + 1 + array(part_count) = calculate_chi(part_x, dir_x, dir_y, & + dir_z, part_e) + cur => cur%next + END DO + ELSE + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_count = part_count + 1 + array(part_count) = 0.0_num + cur => cur%next + END DO + END IF #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90 index e4cfbf5ca..bb8e83680 100644 --- a/epoch1d/src/physics_packages/photons.F90 +++ b/epoch1d/src/physics_packages/photons.F90 @@ -606,7 +606,8 @@ SUBROUTINE qed_update_optical_depth current%optical_depth = current%optical_depth & - delta_optical_depth_photon(chi_val, part_e) - ! If optical depth dropped below zero generate pair... + ! If optical depth dropped below zero generate pair + ! Photon will be deleted IF (current%optical_depth < 0.0_num) THEN CALL generate_pair(current, chi_val, photon_species, & breit_wheeler_electron_species, breit_wheeler_positron_species) @@ -679,7 +680,10 @@ FUNCTION calculate_eta(part_x, part_ux, part_uy, part_uz, & REAL(num) :: e_at_part(3), b_at_part(3) REAL(num) :: beta_x, beta_y, beta_z REAL(num) :: flperp(3), i_e, tau0, roland_eta - REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e + REAL(num) :: moduclip, moduclip2, u_dot_e + REAL(num), PARAMETER :: lambdac = h_bar / mc0 + REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac & + / (2.0_num * alpha_f * m0 * c**3)) CALL field_at_particle(part_x, e_at_part, b_at_part) @@ -690,10 +694,6 @@ FUNCTION calculate_eta(part_x, part_ux, part_uy, part_uz, & beta_y = part_uy / gamma_rel beta_z = part_uz / gamma_rel - lambdac = h_bar / mc0 - - coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3)) - u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) & + part_uz * e_at_part(3)) / moduclip2 diff --git a/epoch2d/Makefile b/epoch2d/Makefile index 84212ce4d..e84ad7dff 100644 --- a/epoch2d/Makefile +++ b/epoch2d/Makefile @@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.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 + partlist.o photons.o laser.o: laser.f90 custom_laser.o evaluate.o mpi_routines.o: mpi_routines.F90 helper.o mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o diff --git a/epoch2d/src/constants.F90 b/epoch2d/src/constants.F90 index 810d9b48b..1bcdcf0af 100644 --- a/epoch2d/src/constants.F90 +++ b/epoch2d/src/constants.F90 @@ -653,7 +653,8 @@ MODULE constants 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_dump_part_qed_chi = 79 + INTEGER, PARAMETER :: num_vars_to_dump = 79 INTEGER, PARAMETER :: c_subset_random = 1 INTEGER, PARAMETER :: c_subset_gamma_min = 2 diff --git a/epoch2d/src/deck/deck_io_block.F90 b/epoch2d/src/deck/deck_io_block.F90 index c9c392d1b..0ba447bd3 100644 --- a/epoch2d/src/deck/deck_io_block.F90 +++ b/epoch2d/src/deck/deck_io_block.F90 @@ -575,6 +575,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode) #ifdef PHOTONS ELSE IF (str_cmp(element, 'optical_depth')) THEN elementselected = c_dump_part_opdepth + + ELSE IF (str_cmp(element, 'qed_chi')) THEN + elementselected = c_dump_part_qed_chi #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch2d/src/io/diagnostics.F90 b/epoch2d/src/io/diagnostics.F90 index 0637579f2..1c367e6ca 100644 --- a/epoch2d/src/io/diagnostics.F90 +++ b/epoch2d/src/io/diagnostics.F90 @@ -644,6 +644,9 @@ SUBROUTINE output_routines(step, force_write) ! step = step index #ifdef PHOTONS CALL write_particle_variable(c_dump_part_opdepth, code, & 'Optical depth', '', it_output_real) + + CALL write_particle_variable(c_dump_part_qed_chi, code, & + 'Chi', '', it_output_real) #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) CALL write_particle_variable(c_dump_part_qed_energy, code, & diff --git a/epoch2d/src/io/iterators.F90 b/epoch2d/src/io/iterators.F90 index e9d6941ac..d9c290b2a 100644 --- a/epoch2d/src/io/iterators.F90 +++ b/epoch2d/src/io/iterators.F90 @@ -18,6 +18,7 @@ MODULE iterators USE particle_pointer_advance USE partlist USE particle_id_hash_mod + USE photons IMPLICIT NONE @@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param) TYPE(particle_list), POINTER, SAVE :: current_list INTEGER :: part_count, ndim REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge - +#ifdef PHOTONS + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y +#endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) END IF @@ -346,6 +349,28 @@ FUNCTION it_output_real(array, npoint_it, start, param) array(part_count) = cur%optical_depth cur => cur%next END DO + CASE (c_dump_part_qed_chi) + IF (current_species%species_type == c_species_id_photon) THEN + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_x = cur%part_pos(1) - x_grid_min_local + part_y = cur%part_pos(2) - y_grid_min_local + norm = c / cur%particle_energy + dir_x = cur%part_p(1) * norm + dir_y = cur%part_p(2) * norm + dir_z = cur%part_p(3) * norm + part_e = cur%particle_energy / m0 / c**2 + part_count = part_count + 1 + array(part_count) = calculate_chi(part_x, part_y, & + dir_x, dir_y, dir_z, part_e) + cur => cur%next + END DO + ELSE + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_count = part_count + 1 + array(part_count) = 0.0_num + cur => cur%next + END DO + END IF #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch2d/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90 index 6e7f40d6b..3ff3cc741 100644 --- a/epoch2d/src/physics_packages/photons.F90 +++ b/epoch2d/src/physics_packages/photons.F90 @@ -610,7 +610,8 @@ SUBROUTINE qed_update_optical_depth current%optical_depth = current%optical_depth & - delta_optical_depth_photon(chi_val, part_e) - ! If optical depth dropped below zero generate pair... + ! If optical depth dropped below zero generate pair + ! Photon will be deleted IF (current%optical_depth < 0.0_num) THEN CALL generate_pair(current, chi_val, photon_species, & breit_wheeler_electron_species, breit_wheeler_positron_species) @@ -683,7 +684,10 @@ FUNCTION calculate_eta(part_x, part_y, part_ux, part_uy, part_uz, & REAL(num) :: e_at_part(3), b_at_part(3) REAL(num) :: beta_x, beta_y, beta_z REAL(num) :: flperp(3), i_e, tau0, roland_eta - REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e + REAL(num) :: moduclip, moduclip2, u_dot_e + REAL(num), PARAMETER :: lambdac = h_bar / mc0 + REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac & + / (2.0_num * alpha_f * m0 * c**3)) CALL field_at_particle(part_x, part_y, e_at_part, b_at_part) @@ -694,10 +698,6 @@ FUNCTION calculate_eta(part_x, part_y, part_ux, part_uy, part_uz, & beta_y = part_uy / gamma_rel beta_z = part_uz / gamma_rel - lambdac = h_bar / mc0 - - coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3)) - u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) & + part_uz * e_at_part(3)) / moduclip2 diff --git a/epoch3d/Makefile b/epoch3d/Makefile index 29f71489f..70cdca79b 100644 --- a/epoch3d/Makefile +++ b/epoch3d/Makefile @@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.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 + partlist.o photons.o laser.o: laser.f90 custom_laser.o evaluate.o mpi_routines.o: mpi_routines.F90 helper.o mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o diff --git a/epoch3d/src/constants.F90 b/epoch3d/src/constants.F90 index 8fe305fad..d28c46d4e 100644 --- a/epoch3d/src/constants.F90 +++ b/epoch3d/src/constants.F90 @@ -654,7 +654,8 @@ MODULE constants 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_dump_part_qed_chi = 79 + INTEGER, PARAMETER :: num_vars_to_dump = 79 INTEGER, PARAMETER :: c_subset_random = 1 INTEGER, PARAMETER :: c_subset_gamma_min = 2 diff --git a/epoch3d/src/deck/deck_io_block.F90 b/epoch3d/src/deck/deck_io_block.F90 index 0b8e85069..d545aca7a 100644 --- a/epoch3d/src/deck/deck_io_block.F90 +++ b/epoch3d/src/deck/deck_io_block.F90 @@ -576,6 +576,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode) #ifdef PHOTONS ELSE IF (str_cmp(element, 'optical_depth')) THEN elementselected = c_dump_part_opdepth + + ELSE IF (str_cmp(element, 'qed_chi')) THEN + elementselected = c_dump_part_qed_chi #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch3d/src/io/diagnostics.F90 b/epoch3d/src/io/diagnostics.F90 index 640ca14fc..0fe270567 100644 --- a/epoch3d/src/io/diagnostics.F90 +++ b/epoch3d/src/io/diagnostics.F90 @@ -653,6 +653,9 @@ SUBROUTINE output_routines(step, force_write) ! step = step index #ifdef PHOTONS CALL write_particle_variable(c_dump_part_opdepth, code, & 'Optical depth', '', it_output_real) + + CALL write_particle_variable(c_dump_part_qed_chi, code, & + 'Chi', '', it_output_real) #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) CALL write_particle_variable(c_dump_part_qed_energy, code, & diff --git a/epoch3d/src/io/iterators.F90 b/epoch3d/src/io/iterators.F90 index e9d6941ac..56b4ccdd0 100644 --- a/epoch3d/src/io/iterators.F90 +++ b/epoch3d/src/io/iterators.F90 @@ -18,6 +18,7 @@ MODULE iterators USE particle_pointer_advance USE partlist USE particle_id_hash_mod + USE photons IMPLICIT NONE @@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param) TYPE(particle_list), POINTER, SAVE :: current_list INTEGER :: part_count, ndim REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge - +#ifdef PHOTONS + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y, part_z +#endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) END IF @@ -346,6 +349,29 @@ FUNCTION it_output_real(array, npoint_it, start, param) array(part_count) = cur%optical_depth cur => cur%next END DO + CASE (c_dump_part_qed_chi) + IF (current_species%species_type == c_species_id_photon) THEN + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_x = cur%part_pos(1) - x_grid_min_local + part_y = cur%part_pos(2) - y_grid_min_local + part_z = cur%part_pos(3) - z_grid_min_local + norm = c / cur%particle_energy + dir_x = cur%part_p(1) * norm + dir_y = cur%part_p(2) * norm + dir_z = cur%part_p(3) * norm + part_e = cur%particle_energy / m0 / c**2 + part_count = part_count + 1 + array(part_count) = calculate_chi(part_x, part_y, part_z, & + dir_x, dir_y, dir_z, part_e) + cur => cur%next + END DO + ELSE + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_count = part_count + 1 + array(part_count) = 0.0_num + cur => cur%next + END DO + END IF #endif #if defined(PHOTONS) || defined(BREMSSTRAHLUNG) diff --git a/epoch3d/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90 index d45420cf3..ae0ab0d62 100644 --- a/epoch3d/src/physics_packages/photons.F90 +++ b/epoch3d/src/physics_packages/photons.F90 @@ -612,7 +612,8 @@ SUBROUTINE qed_update_optical_depth current%optical_depth = current%optical_depth & - delta_optical_depth_photon(chi_val, part_e) - ! If optical depth dropped below zero generate pair... + ! If optical depth dropped below zero generate pair + ! Photon will be deleted IF (current%optical_depth < 0.0_num) THEN CALL generate_pair(current, chi_val, photon_species, & breit_wheeler_electron_species, breit_wheeler_positron_species) @@ -685,7 +686,10 @@ FUNCTION calculate_eta(part_x, part_y, part_z, part_ux, part_uy, part_uz, & REAL(num) :: e_at_part(3), b_at_part(3) REAL(num) :: beta_x, beta_y, beta_z REAL(num) :: flperp(3), i_e, tau0, roland_eta - REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e + REAL(num) :: moduclip, moduclip2, u_dot_e + REAL(num), PARAMETER :: lambdac = h_bar / mc0 + REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac & + / (2.0_num * alpha_f * m0 * c**3)) CALL field_at_particle(part_x, part_y, part_z, e_at_part, b_at_part) @@ -696,10 +700,6 @@ FUNCTION calculate_eta(part_x, part_y, part_z, part_ux, part_uy, part_uz, & beta_y = part_uy / gamma_rel beta_z = part_uz / gamma_rel - lambdac = h_bar / mc0 - - coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3)) - u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) & + part_uz * e_at_part(3)) / moduclip2